Wormholes and masses for Goldstone bosons

There exist non-trivial stationary points of the Euclidean action for an axion particle minimally coupled to Einstein gravity, dubbed wormholes. They explicitly break the continuos global shift symmetry of the axion in a non-perturbative way, and generate an effective potential that may compete with QCD depending on the value of the axion decay constant. In this paper, we explore both theoretical and phenomenological aspects of this issue. On the theory side, we address the problem of stability of the wormhole solutions, and we show that the spectrum of the quadratic action features only positive eigenvalues. On the phenomenological side, we discuss, beside the obvious application to the QCD axion, relevant consequences for models with ultralight dark matter, black hole superradiance, and the relaxation of the electroweak scale. We conclude discussing wormhole solutions for a generic coset and the potential they generate.


Introduction
The explicit breaking of global symmetries due to gravitational effects is a topic of great relevance in physics, both for theoretical and phenomenological reasons. Its study might shed light on the reason why gravity is so different from the other fundamental interactions, and what rôle did these special features play in the history of the Universe.
Let us focus for simplicity and definiteness on the case of a U (1) Peccei-Quinn (PQ) global symmetry that is spontaneously broken at some scale f a by the vacuum expectation value (VEV) of a complex scalar field Φ. By introducing the polar parametrization Φ = ρe ıφ/fa , the VEV of the complex field is |Φ| = f a , and the axion can be identified with the angular direction φ -the Nambu-Goldstone boson arising from the spontaneous symmetry breaking. Below the scale f a , the axion enjoys the non-linearly realized global U (1) for arbitrary values of the parameter α describing the shift in field space. The conserved Noether current implied by the symmetry in eq. (1) is J µ = f a ∂ µ φ, and the associated charge Q generates the symmetry transformation. For the purpose of the present discussion, it is important to notice that the subgroup of eq. (1) given by the discrete shift represents a gauge symmetry of the system. This symmetry is inherent to the mere definition of the axion field in terms of an angular dynamical variable. Referring for definiteness to the U (1) case discussed above, the discrete shift in eq. (2) maps the complex field Φ into itself, and thus corresponds to a redundancy in the description of the physical degrees of freedom. At the QCD scale Λ QCD , Yang-Mills instanton effects break the continuos shift in eq. (1) into its discrete subgroup φ → φ + 2kπf a generating the potential minimized at φ = 0 in agreement with the axion solution of the strong CP problem. What does gravity do with global symmetries, like the axion shift symmetry in eq. (1)? A number of 'folk theorems' state the impossibility to have global symmetries in a consistent theory of quantum gravity. To support this hypothesis, there exists a general semi-classical argument based on black hole physics. The reasoning is that if some amount of global charge is thrown into a black hole, then the subsequent thermal decay of the black hole into photons and gravitons via Hawking radiation destroys the global charge without any possibility to reconstruct it in the final state, thus defining a process by means of which the global charge is violated. To this respect, one also understands the crucial difference between global and local symmetries. Local symmetries obey the Gauss's law, and any observer outside the black hole can determine its charge. The electric charge of a Kerr-Newman black hole represents an explicit example of this fact. On the contrary, there is no Gauss's law associated with global symmetries, and when charged particles are thrown into the black hole, there is no way to track them from the outside: The global charge appears to be deleted, in obvious contradiction to its conservation.
In addition to the aforementioned theorem related to black hole physics, various arguments in perturbative string theory [1] -where global symmetries on the world-sheet become gauge symmetries in the target space -and AdS/CFT [2] -where global symmetries on the boundary correspond to gauge symmetries in the bulk -seem to corroborate this common belief.
Motivated by these arguments, many authors studied the consequences of the explicit breaking induced by higher-dimensional operators suppressed by the Planck scale M Pl = 1/G N of the form [3,4] V gravity (Φ) = g 2m+n with mass-dimension 2m + n, and generic coupling g 2m+n that in principle is of order unity. Needless to say, the impact of these operators is catastrophic. To give an idea, in the QCD axion case including dimension-5 symmetry breaking operator induced by Planck scale physics would require a coupling of order |g 5 | 10 −55 in order to avoid dangerous CP violation effects for f a ∼ 10 12 GeV. This is of course unacceptable, since it introduces a fine-tuning by far more severe than the one the axion claims to solve. The higher-dimensional operators in eq. (4) therefore represent, if present, a threat for the axion solution of the strong CP problem. To solve the issue, it is possible to tailor suitable extensions of the simple U (1) model discussed at the beginning of this section. For instance, if the axion global symmetry is an accidental symmetry descending from an exact (discrete) gauge symmetry, then the problematic higher-dimensional operators can be forbidden to very high order [5][6][7][8][9][10][11][12] (a similar conclusion remains true considering the axion global symmetry as an accidental symmetry descending from exact discrete global symmetries [13,14]).
However, before engineering possible solutions, it is important to pose the following question: Does gravity really generate power-suppressed symmetry breaking operators like those in eq. (4)?
Let us tackle the problem from the simplest perspective, and consider the action describing an axion field minimally coupled to Einstein gravity. In this setup, it is possible to show that the global symmetry in eq. (1) remains intact at any finite order in a perturbative expansion in G N , and no power-suppressed operators are generated. Albeit surprising, this is not in contrast with the arguments presented at the beginning of this section about the breaking of global symmetries due to gravity. What is the origin of this conundrum? First of all, the reader should keep in mind that the 'folk theorems' mentioned above have to do with black holes, that is with non-perturbative objects. This fact suggests that the breaking of global symmetries induced by gravity is to be found at the non-perturbative level, and the lack of perturbative effects is not, after all, surprising. At the same time, this may sound discouraging -in particular from a phenomenological perspective -since referring to non-perturbative gravitational effects seems to be nothing but a vague and unclear indication. Is there a computable non-perturbative effect showing indisputably that gravity breaks explicitly global symmetries?
A careful analysis confirms the expectation sketched above, and gives a positive answer to the aforementioned question. Gravity breaks the global symmetry in eq. (1) down to its discrete version in eq. (2) in a non-perturbative way [15]. In more detail, Euclidean wormhole solutions [16] swallow axionic charge, break the shift symmetry it generates, and create a potential for the axion that may compete, depending on the value of the axion decay constant, with the QCD effect in eq. (3). These effects are, as promised, controllably calculable -up to some extent that we shall discuss in detail -in the context of Einstein gravity, thus relying surprisingly little on our ignorance of its ultraviolet (UV) completion.
This paper is structured as follows. In section 2 we review the Euclidean wormhole solutions [16] in the context of Einstein gravity. Section 3 contains the computation for the spectrum of the quadratic action, and in section 4 the effective potential for a generic U (1) Goldstone boson generated by Euclidean wormholes can be found. This concludes the theoretical part of this work. In section 5 we illustrate several phenomenological implications, including the QCD axion, the case of ultralight scalar dark matter, and the relaxation of the electroweak scale. In section 6 possible UV completions are discussed, and their impact on the existence of wormhole solutions assessed. In section 7, we derive wormhole solutions for a generic coset and look into O(n + 1)/O(n) as an specific example. Conclusions can be found in 8. Since part of this letter reviews existing material, for quick reference main results can be found in green boxes.

The Euclidean wormhole solution
In this section we review the wormhole instanton solution first discussed in [16]. The computation we are about to discuss is not novel and can be found in the literature (see, e.g., [15,[17][18][19][20]). However, we consider useful and important to re-elaborate its derivation both for completeness and to highlight the points that will matter most for the rest of this work.

The bulk action
The starting point is the Euclidean action of a three-form H coupled to Einstein gravity 1 where M Pl = 1/ √ G N 1.22 × 10 19 GeV is the Planck mass, and F ≡ 1/(3! f 2 a ), with f a the PQ scale. The (dimensionless) pseudo-scalar axion field θ is related to the three-form H via the relation H µνρ = f 2 a µνρσ (∂ σ θ).
The variation of S E w.r.t. the metric is and the corresponding Einstein equation reads By taking the trace of eq. (7) we find −R = 8πG N FH 2 . The equation of motion takes the form where we write the three form as the field strength of an antisymmetric tensor B µν , i.e. H µνρ = ∂ µ B νρ + ∂ ρ B µν + ∂ ν B ρµ . In terms of the angular axion field θ the Euclidean action in eq. (5) reads The canonically normalized axion field is φ ≡ f a θ. Notice that the kinetic term in eq. (8) has the 'wrong' sign if compared to the Euclidean action of an ordinary scalar. At the technical level, this is related to the Levi-Civita contraction µνρσ µνρλ = (−1) t 3! δ λ σ , where t is the number of negative eigenvalues of the metric: Continuing from Minkowski (t = 1) to Euclidean (t = 0) space using the dual description, the axion kinetic term flips its overall sign. 2 We refer to appendix A for a more detailed explanation and careful derivation.
The same property appears at the level of the energy-momentum tensor. From eq. (7), we find where we used the contraction αβµν αβλη = (−1) t 2 δ µ λ δ ν η − δ µ η δ ν λ with again t = 0 in the Euclidean case. Note that this energy-momentum tensor has the opposite overall sign if compared with the one of an ordinary real scalar. The flipped sign in eqs. (8,11) plays a crucial rôle in what follows. The mere existence of wormhole solutions -as we shall see -is a peculiar property of massless free pseudo-scalar (scalar) fields, which admit a dual description in terms of a two-index antisymmetric tensor (pseudo-tensor) whose field strength is the three-form H in eq. (5). Indeed, even if we focus on the case of a pseudo-scalar axion field, it is important to remark that the construction we put forward in this section remains valid also for a genuine scalar field.
We now look for a spherically symmetric solution with the following ansatz ds 2 = α 2 (r)dr 2 + β 2 (r)dΩ 2 D−1,K , where r is the time coordinate in Euclidean space, and the D − 1 space has curvature K. Without loss of generality, we write dΩ 2 D−1,K =ĝ ij dϕ i dϕ j , with metric determinant √ g = αβ 3 √ĝ . K = 1 corresponds 2 Notice that is a genuine tensor under general coordinate transformations, and we have the relation where the tensor density ε of weight −1 is the usual total antisymmetric Levi-Civita symbol, defined to be ±1 and 0 in all frames. In eq. (9), |g| = √ g if t = 0 (or more generally if t is even), and |g| = √ −g if t = 1 (or more generally if t is odd). Finally, in full generality, we have  (17). The wormhole throat is characterized by the axion charge n (see eq. (18)). This characteristic allows for the possibility to interpret the wormhole as the combination of an instanton with axion charge +n and an anti-instanton with axion charge −n (right panel).
On the contrary, the Bianchi identity µνρσ (∂ ρ H σµν ) = 0 provides a non-trivial constraint. Using the field ansatz, and bearing in mind that ijk = β 3 √ĝ ε ijk (see footnote 2), we find ∂ r (H ijk ) = ∂ r (Hβ 3 √ĝ ε ijk ) = 0. As a direct consequence, we have the r dependence H(r) = H 0 /β 3 (r), with H 0 dimensionless constant. To make contact with the original solution proposed in [16], we now restrict the analysis to the three-sphere S 3 , with K = 1 and β 2 (r) = r 2 . The constant H 0 can be fixed by normalizing the integral of the field strength over S 3 . We have where -as we shall see later -n is an integer number. We can now use eq. (14). Solving for α, we find The metric in eqs. (12,17) describes the Giddings-Strominger wormhole configuration. The wormhole connects two asymptotically flat Euclidean region, as shown in the left panel of fig. 1. At any fixed Euclidean time r, the cross-section of the wormhole is a three-sphere. The length scale L in eq. (17) defines the minimal radius of the wormhole throat. Let us note that there is no singularity at L = 0; this is indeed nothing but a coordinate singularity, as we shall discuss later.
Using the relation H µνρ = f 2 a µνρσ (∂ σ θ) = µνρσ J σ and the field ansatz in eq. (12) we can compute the axion charge stored in a cross-section of the wormhole throat. We find The wormhole throat is therefore characterized by the axion charge n. This characteristic offers the possibility to consider the wormhole as the combination of an instanton with axion charge +n and an anti-instanton with axion charge −n ( fig. 1, right panel). Furthermore, as we shall discuss at the end of this section, n is quantized, and it only takes integer values. As discussed in [16], separating the wormhole in two halves has an important topological interpretation. On the one hand, the wormhole interpolates between two R 3 geometries widely separated in time. On the other one, the instanton (anti-instanton) represents a topology change from Instead of considering the wormhole as connecting two separate asymptotic Euclidean spaces -a configuration that has an unclear physical interpretation -it is possible to consider wormholes that connect back to the same space, as shown in fig. 2. By slicing the wormhole in half, the instanton/anti-instanton represents, as before, a topology change from R 3 to R 3 ⊕ S 3 and viceversa.
From eq. (5), we can compute the action of an instanton configuration with charge n. Using R = −8πG N FH 2 , we find 3 This is the right moment to do some useful dimensional analysis, and it is convenient to restore the appropriate powers of (while keeping c = 1 First, we can exploit the analogy with the electroweak theory. The Fermi constant G F (or equivalently the electroweak VEV) is related to the ratio between the electroweak mass (defined here by the mass of the W boson, and corresponding to the physical threshold of energy E ∼ M W at which the new degrees of freedom in the electroweak sector become dynamical) and the electroweak coupling constant g W by means of Similarly, in gravity we expect the relation where we introduced a string mass scale M S and a string coupling g S . If the UV completion of general relativity is weakly coupled, we expect new degrees of freedom to occur even below the Planck scale. If, on the contrary, the UV completion of general relativity is strongly coupled, g S 1, we expect the new dynamics of quantum gravity to occur above M Pl . We shall return on this point later, discussing the validity of the wormhole solutions (see section 6.2). Second, from eq. (19) we recover the correct dimensionality of an instanton action since we have that On general ground, this also explains why instantons are non-perturbative. If g → 0, exp(−S inst / ) goes to zero faster then any power of g, and instanton effects cannot be captured at any finite order in perturbation theory.
Finally, notice that the wormholes with action given in eq. (19) realize the Weak Gravity Conjecture (WGC) proposed in [21]. The generalization of the WGC to the case of axions states that for an axion with decay constant f a there exists an instanton with action S inst M Pl /f a . In the context of Einstein gravity, this is precisely the wormhole with n = 1 and S inst 0.38 M Pl /f a .
Let us now explore in more detail our solution. The wormhole metric is conformally equivalent to a flat metric on S 3 × R. The coordinate transformation brings the metric in eq. (21) into the form Notice that r = L is mapped to τ = 0 and in this coordinate system is plain to see that there is no singularity and the wormhole connects the two asymptotically flat regions at τ = ±∞. With respect to this metric, the Einstein equations and the equation of motion for the axion field, which can be obtained from taking α = β = a in eqs. (14,15), are where we defined κ ≡ 8π/M 2 Pl . From the first Einstein equation, we find . Integrating over τ , we find In the left panel of fig. 3 we show the instanton solution in eq. (26). Asymptotically, the axion field configuration approaches the constant values f 2 a κθ(τ )/ √ 6 τ →±∞ → ±π/4, and features a sharp transition in between. Notice that the asymptotic limit is exponentially fast, tanh(τ ) τ →±∞ → ± 1 − e −2τ . Since instantons are welllocalized objects, there are also approximate solutions consisting of strings of widely separated instantons and anti-instantons centered at τ 1 , . . . , τ m . Notice that the possibility to center the instanton-anti-instanton pair at an arbitrary position in τ is due to the fact that the wormhole solution in eqs. (23,25) possesses a 'time'translation symmetry. Indeed, the most general solution of the field equations is , and it depends on the arbitrary constant τ 0 . As we shall see in section 3.1, the presence of this 'time'-translation symmetry will imply the existence of a zero eigenvalue in the spectrum of the quadratic action describing fluctuations around the wormhole solution. For illustrative purposes, we show a multi-instanton solution in the right panel of fig. 3. Since the m instantons are widely separated, the classical action is just mS inst .
Finally, we close this section with the following remark. As already anticipated, the axion charge in eq. (18) turns out to be quantized. This property can be understood as follows. We consider the action for the axion field in the background metric of eq. (23), and, after introducing the new variable dτ /dσ ≡ a 2 and considering radially symmetric field configurations θ(τ (σ)), we find With respect to the new variable, the equation of motion of the axion field is simply given by d 2 θ/dσ 2 = 0. Equivalently, using eqs. (23,25), we find dθ/dσ = const = n/2π 2 f 2 a . We can now expand eq. (27) around the instanton solution, θ → θ + δθ. We find, using the equation of motion, the following change in the action In the presence of wormholes, the action is no longer invariant under the generic shift θ → θ + δθ. We now restrict to the discrete transformations δθ = 2kπ, with k ∈ Z. Since θ is an angular variable, we know that θ → θ + 2kπ corresponds to the same physical point, as discussed in section 1. This is a gauge redundancy, and, as such, it cannot be broken. The simple observation is that quantum mechanical transition amplitudes involve e ıS θ , and the phase exponential must be gauge invariant. This means that the only change δS θ that can be tolerated corresponds to an integer multiple of 2π. We therefore find the quantization condition n ∈ Z.
The quantization of the wormhole charge Q -and the subsequent invariance of the action under the transformation θ → θ + 2kπ, with k ∈ Z -is a crucial point, since it is related to the fact that in the presence of wormholes the continuos shift symmetry in eq. (1) is broken down to its discrete subgroup in eq. (2). We shall elaborate more on this point in section 4.
The quantization of the wormhole charge Q is tightly related to the periodicity of the axion field, and it has a profound mathematical explanation. It is indeed possible to show that n can only be either an integer or zero. The point is that if the flux of the three-form H through S 3 is non-zero, then in general a unique two-form potential B -an antisymmetric tensor whose filed strength is H according to H µνρ = ∂ µ B νρ + ∂ ν B ρµ + ∂ ρ B µν , see appendix A -covering the whole sphere S 3 does not exist. Only if n is an integer the ambiguity can be bypassed by defining two potentials B 1,2 , covering the upper and lower halves of S 3 , and related to each other by a gauge transformation in an overlap region which is topologically S 2 [24][25][26]. The argument is very similar to the Dirac quantization of the electric charge in the presence of a magnetic monopole. We discuss in more detail this point in appendix A.3. The analogy is indeed deeper, and the Dirac quantization condition corresponds to the quantization of the axion charge over a spatial slice S 3 [27,28].

The Gibbons-Hawking-York boundary term
In the presence of a manifold with boundary, the Einstein-Hilbert action in eq. (5) must be supplemented by a boundary term so that the variational principle is well-defined [29][30][31]. To make this point clear, let us only focus on the term involving the Ricci scalar for the metric ansatz in eq. (23) with R = 6(a − a )/a 3 The presence of the second derivative prevents the correct formulation of the variational principle. 4 It is therefore necessary to add the Gibbons-Hawking-York boundary term, whose general form is given by 4 As a prototype, consider the action S = t 2 In order to extract the equation of motion by means of the variational principle, it is necessary to specify boundary data prescribing the values of the variation at the boundaries δq(t 2 ), δq(t 1 ). Let us now suppose to add a total derivative term The equation of motion remains unchanged, since a total derivative does not affect the dynamics. However, in order to apply the variational principle, four boundary data are needed, namely δq(t 2 ), δq(t 1 ), δq(t 2 ), δq(t 1 ). This is inconsistent since fixing position and velocity violates the uncertainty principle. It is therefore necessary to add a boundary term which cancels the extra total derivative inL. This example mimics the typical situation in general relativity -and, in particular, the case of eq. (29) -with the rôle ofL played by the Einstein-Hilbert action. In fact, the Ricci scalar has second derivative of the metric, and its variation generates a total derivative that is precisely canceled by the addition of the Gibbons-Hawking-York boundary term.
where h is the determinant of the induced metric on the boundary manifold ∂M , K is the trace of the extrinsic curvature tensor, and K 0 is the trace of the extrinsic curvature tensor of the same boundary embedded in flat spacetime. The explicit computation of S GHY is straightforward. The trace of the extrinsic curvature tensor is K = ∇ µ n µ = ∂ µ n µ + Γ µ µν n ν , where n µ is the unit normal to ∂M . The boundary we are interested in are hypersurfaces of constant τ , and the corresponding unit normal vector is n τ = 1/a. By direct computation, we find K = 3a /a 2 . As far as K 0 is concerned, we use a = e τ for the flat metric, and we find K 0 = 3/a. All in all, the Gibbons-Hawking-York boundary term reduces to and it cancels the total derivative in eq. (29). The Gibbons-Hawking-York boundary term must be included for a correct formulation of the gravitational action in the presence of a manifold with boundary, and this is the case for the half-wormhole configuration. It is therefore important to evaluate the topological contribution of eq. (33) to the bulk action in eq. (19). Using the explicit solution in eq. (23), we find where only the boundary at the wormhole throat gives a non-vanishing contribution. Taking into account the boundary term, the action of an instanton configuration with charge n is Half-wormhole action 3 Quadratic action and the absence of negative modes Instantons are responsible for one of the most important effect in Quantum Mechanics, that is the tunneling through a potential barrier. The importance of such effect lies in the fact that it changes the vacuum of the ������ ��������� system, and there exist two paradigmatic situations in which this effect is particularly manifest [22,23].
In the left panel of fig. 4 we show a one-dimensional quantum system featuring a potential with a double well, of the form V (x) = λ(x 2 − k 2 ) 2 . In perturbation theory, the ground state has two degenerate minima at x = ±k, which are equivalent because of the reflection symmetry of the system (dashed black lines). Instantons lift the degeneracy, and the energy difference between the true ground state and the first excited state (solid black lines -respectively, the anti-symmetric and symmetric combination of the two vacua at x = ±k) is controlled by the instanton action, according to the qualitative scaling E 1 (g) − E 0 (g) ∼ e −const/g 2 .
In the right panel of fig. 4 we show a potential with a metastable vacuum, of the form V (x) = mω 2 x 2 +λx 3 . In perturbation theory -assuming small λ -one can compute the energy of the ground state of the system at, in principle, any order (solid black line). However, perturbation theory can not describe the tunneling through the barrier, and the instability of the vacuum. Again, the tunneling is captured by instanton effects, that in this case introduce a small imaginary part in the ground state energy, E 0 (g) = Re E 0 (g) + ıIm E 0 (g), with Im E 0 (g) ∼ e −const/g 2 , which in turn becomes a decay probability in the time evolution e ıE0(g)t . In this case the instanton is called bounce.
In order to understand the true nature of the instanton, it is necessary to inspect the action at the quadratic order. Let us clarify this point with an explicit example. In Euclidean time, the quantum probability amplitude from the initial state |i at time −T /2 to the final state |f at time T /2 -that corresponds, for instance, to the transition amplitude between the two vacua in the quantum systems discussed in fig. 4 -is given by where N is a normalization constant, S E [x cl (τ )] is the Euclidean action computed for the classical instanton solution x cl (τ ), and the prefactor encodes the quantum fluctuation around the classical trajectory. The latter, as a result of a Gaussian integration, takes the shape of a determinant with a prime denoting the subtraction of zero eigenvalues. If we take T to be very large, the lowest energy state dominates the sum, and we can extract the ground state energy from the computation of the path integral. The spectrum of the differential operator δ 2 S E /δx 2 is therefore crucial. In the presence of one (or an odd number of) negative eigenvalue(s), there is an imaginary contribution to the vacuum energy, and the instanton plays the rôle of a bounce. If the spectrum of the quadratic action is positive (or if the number of negative eigenvalues is even) the vacuum energy is real, and the instanton describes the transition between degenerate vacua, lifting their degeneracy as discussed in the case of the double well. The analogy with one-dimensional systems in Quantum Mechanics suits particularly well our case since the wormhole solutions discussed in section 2 are characterized by one relevant dimension, that is the time coordinate in the Euclidean 4D space.

Bulk action and boundary terms
In this section we study, motivated by the previous discussion, the spectrum of the quadratic action describing perturbations around the Euclidean wormhole solution. To this end, we shall exploit the formalism developed in the context of cosmological perturbation theory for the analysis of fluctuations of a scalar field ϕ 0 coupled to gravity in a Friedmann-Robertson-Walker metric [32]. This study requires the canonical formalism in order to correctly get rid of the non-dynamical degrees of freedom, as in any other gauge theory. The canonical formulation, although non-explicitly, preserves gauge invariance.
In the following, we start considering metric perturbations in Minkowski space around the background line element ds 2 = a 2 (η)(−dη 2 + γ ij dx i dx j ). Including fluctuations, in full generality we write ds 2 = a 2 (η) −(1 + 2A)dη 2 + 2B i dx i dη + (γ ij + h ij )dx i dx j . This perturbed line element can be further analyzed by means of the scalar-vector-tensor decomposition, which separates the fluctuations into components according to their transformations under spatial rotations. In particular, the three-vector B i splits into the gradient of a scalar, and a divergence-free vector B i = ∂ i B +B i while the symmetric tensor h ij can be decomposed into a scalar part, a divergence-free vector, and a trace-free transverse tensor h ij = −2ψγ ij +2∂ i ∂ j E +2∂ (iÊj) +t ij . 5 The advantage of this decomposition is that scalar, vector and tensor Einstein equations are decoupled at the linear order and can therefore be treated separately. Vector perturbations are pure gauge modes, and they carry no dynamics [32]. As far as tensor perturbations are concerned, they do not couple to scalar modes, and we can therefore rely on standard computations of gravitational waves in instanton background [33]. We found a positive spectrum, and we summarize our computation in appendix C.
We focus here on scalar perturbations. The perturbed line element is therefore where we simultaneously expand the scalar field, assumed to be function of η, as ϕ 0 + δϕ. The second order Lagrangian in terms of the perturbations in eq. (37) can be found in e.g. [32], but it has, at this level, non-dynamical degrees of freedom and gauge redundancy. First A and B are not dynamical, that is they do not have a 'kinetic term'; they behave like A 0 in a 'conventional' gauge theory. In addition, there is a two-parameter gauge transformation, corresponding to diffeomorphism acting on scalar perturbations, given by δx µ = (λ 0 , ∂ i λ), and we have where we indicate with˙the derivative w.r.t. the time variable η. Out of the 5 perturbations A, B, E, ψ, δϕ, 2 are non dynamical, and 2 can be removed via a gauge transformation which leaves us with 1 d.o.f. All in all, we have the following quadratic action in Hamiltonian form pq − H (see [32] for details of the derivation) where the dynamical variable Ψ and its conjugate momentum Π Ψ in terms of the metric perturbations in eq. (37) are Ψ ≡ ψ + (ȧ/a)δϕ/φ 0 and Π Ψ ≡ Π ψ − 2a 2 √ γ(∆ + 3K)δϕ/κφ 0 , where Π ψ is the conjugate momentum of ψ and ∆ is the Laplacian on the three-sphere associated with γ ij . The background solutions satisfy ȧ a In the following we focus on the spatially homogeneous O(4) invariant modes, and we formally set ∆ = 0. 6 As discussed in [34,35], the homogeneous modes are those that are potentially responsible for the presence of negative eigenvalues in the spectrum of the quadratic action, and they need to be carefully analyzed. Indeed, in [36] the authors claimed the existence of a negative mode among the homogenous fluctuations around the Giddings-Strominger wormholes. If true, this result would imply that wormholes are bounce solutions mediating unstable vacuum decay via tunneling transition. Motivated by the necessity to further explore and check this result, we now turn to discuss our own analysis. For completeness, in appendix C we discuss the spectrum of inhomogeneous scalar perturbations.
Restricting to the homogeneous modes, eq. (39) becomes 6 In a closed Universe, K = 1, the line element on S 3 is dl 2 = γ ij dx i dx j = dψ 2 + sin 2 ψ(dφ 2 + sin 2 φdϕ 2 ). The scalar spherical harmonics Q are the usual spherical harmonics on S 2 and Π n l (ψ) are the Fock harmonics. The eigenfunction corresponding to n = 1 with zero eigenvalue defines the spatially homogeneous modes. In a open Universe, K = −1, the Laplacian is again zero for the spatially homogeneous modes, and takes the values −p 2 − 1, with p 2 > 0, for the continuum of square integrable modes.
The first difficulty one faces is the conformal factor problem [37]. 7 The conformal factor problem can be spotted in the Hamiltonian of eq. (41) by using the equation of motion 6ȧ 2 − κa 2φ2 0 = −6Ka 2 . As a consequence, we find a 'wrong sign' kinetic term since we will be looking at the closed case K = 1. The literature often treats this problem by just turning Ψ into the imaginary axis but this, other than pretty arbitrary, will lead to augmented confusion when convoluted with the rotation to Euclidean. Instead the most sensible way in our humble opinion is to perform a canonical (symplectic) transformation and study the resulting system as done in [34]. This will correspond roughly to swapping Ψ and Π Ψ and then performing another transformation, explicitly and in a single step which is symplectic since and leads to the action where we discarded boundary terms since the equation above is the starting point for a canonical formulation, as discussed in sec. 2.2 and in accordance with [35]. The canonical transformation above might appear to have no other virtue than to simplify the Hamiltonian, so it is worth pausing and looking at its connection with the original variables and physical interpretation. Indeed in ref. [32] a similar canonical transformation is used such that the resulting coordinate is proportional to Bardeen's potential. In our case we have, after substituting in the canonical transformation: which, under a gauge transformation as in eq. (38) stays invariant. The gauge invariance of the variable speaks to the success of the proccess of removing redundancies and is a valuable check on the procedure here employed. 7 The natural definition of path integral in pure gravity involves the Euclidean partition function where the sum runs over all metrics on a four-dimensional manifold M with boundary ∂M. The Euclidean path integral is ill-defined due to the fact that the scalar curvature can be made arbitrarily negative. The physics behind is related to the fact that the gravitational potential energy is negative because gravity is attractive. A concrete way to look at the problem is to perform the conformal transformation gµν →gµν = Ω 2 gµν which brings the Euclidean action into the form The Euclidean action can be made arbitrarily negative by choosing a rapidly varying conformal factor. In turn, this implies that the path integral diverges since one has to integrate over all possible Ω. An alternative formulation of the problem, closer in spirit to the content of this section, consists in performing an expansion of the Euclidean action around a fixed on-shell metric g obeying the classical field equations Rµν [g] = 0. At the quadratic order in the fluctuation hµν we have [38] δS where we decomposed the fluctuation into its trace and a trace-free part, hµν = 1 4 hgµν + h T µν . The quadratic action in eq. (44) reveals that the trace part of the perturbation is unbounded from below, since the eigenvalue equation for the Laplace-Beltrami operator, −∇ 2 u = λu, admits only positive eigenvalues, λ 0, in a compact Riemannian manifold (a result known as the Lichnerowicz-Obata theorem).
Notice that we can write this action, since there is no 'potential' in the Hamiltonian for q, as with p = ∂H /∂q. It would seem that the discussion of the sign of eigenvalues ends here: they are ≥ 0 from construction. The only non-trivial step to corroborate this conclusion is the rotation to Euclidean, which as it turns out changes nothing. However we pursue here further to find the explicit eigenvalues of the spectrum. After usingq = ∂H /∂p to write p in terms of q,q, repeated use of equations of motion and integration by parts leads to We note that the bulk term can be generalized to the Lorentz invariant form d 4 x √ −g(∂ µ q∂ µ q + . . . ). As the last step in Minkowski we substitute q → q/a, and we find Finally, going to the Euclidean, τ = iη and with () = d/dτ (52) At this point, an important comment is in order. In addition to the usual analytical continuation, in eq. (52) we performed the formal substitution ϕ 0 → ıf a θ. This is because the starting point of our analysis, that is the quadratic action in eq. (39), was derived in [32] considering a generic real scalar field ϕ 0 . In our case, on the contrary, we are interested in the case of an axion field φ (or, more generically, in the case of a Goldstone boson which admits a three-form description). Going from Minkowski to Euclidean, in the axion case one gets an extra minus sign whenever a term quadratic in θ appears (as discussed in section 2.1 and appendix A.2), and the imaginary factor in ϕ 0 → ıf a θ precisely accounts for this issue. For instance, this replacement is indeed crucial to match eq. (24) starting from the corresponding Minkowski version in eq. (40).
Substituting the wormhole solution for the background fields a and κf 2 a θ 2 as given in eq. (25), and considering the case of a closed Universe K = 1, we have the following quadratic action and, consequently, the following operator in the bulk Homogeneous scalar fluctuations around the wormhole background This action can be written as the square of a generalized momentum, as the analysis in Minkowski indicated, explicitly Remarkably, the spectral problem for the quadratic action describing fluctuations around the wormhole background reduced to the eigenvalue problem for a Schrödinger-type operator. In particular, the differential operator in eq. (54) describes the one-dimensional motion of a particle subject to the Pöschl-Teller potential (see appendix B).
The possible presence of negative eigenvalues is therefore related to the existence of bound states. Indeed, in general O has a discreet spectrum of bound states and for higher energies a continuum. The eigenvalues have definite 'partiy' under τ → −τ and this will be determining. In figure 5 we show the odd (left panel) and even (right panel) discrete eigenvalues. Considering the eigenvalue equation Oq λ (τ ) = λq λ (τ ), they Figure 5: Spectrum of the differential operator in eq. (54) that coincides with the Pöschl-Teller potential well in Quantum Mechanics. The green band corresponds to the continuum spectrum. In the left (right) panel we show the parity odd (parity even) eigenfunction q − λ=0 (τ ) (q + λ=−8 (τ )). We do not show the parity even eigenfunction q + λ=0 (τ ) since it is not square-integrable.

correspond to
Odd eigenfunction : Even eigenfunctions : At this stage, one would naïvely conclude that the differential operator in eq. (54) has one negative eigenvalue; this is, the reader might have noticed, in contradiction with eq. (55). The resolution of this conflict requires closer inspection, in particular, one needs to check the behaviour of the eigenfunctions at the boundaries in eq. (52) and the transformation property under parity. Let us start from the latter. Parity under τ → −τ (∼ η → −η) can be assigned to the perturbations in eq. (37), and carried on to the expression for q in eq. (48) we see that it has odd parity, a reminder of the fact that the canonical transformation has 'swapped' momentum and coordinate. In eqs. (56,57), the only eigenfunction that is compatible with this requirement is q − λ=0 (τ ), with eigenvalue λ = 0. Based on this parity argument, we therefore discard the eigenfunctions q + λ=−8 (τ ) and q + λ=0 (τ ). Second, we have to check the behaviour of the eigenfunctions at the boundaries in eq. (53). This can be done in two ways: i) direct computation in eq. (53) where one finds that the boundary term cancel at τ ± ∞ for q − λ=0 (it does cancel at τ = 0 as well if one instead considers a half-wormhole), ii) alternatively substitution of q − λ=0 in eq. (55) yields 0 which comprises both bulk and boundary contributions in eq (53). In this regard, let us note that the even solutions put in eq. (55) yield ∞ when integrating around the throat of the wormhole. In this sense only q − λ=0 is an eigenfunction of both eq. (53) and eq. (55) which adds to the evidence in favour of the odd eigenfunction.
Notice that otherwise the presence of a zero eigenvalue can be inferred from the 'time'-translation symmetry for the wormhole solution, as discussed in section 2.1.
All in all, we found that homogeneous scalar perturbations give a positive contribution to the fluctuation determinant in eq. (36). This result supports the interpretation of the wormhole as an instanton mediating tunneling transitions between degenerate vacua. Equipped with this result, we can now move to compute the effective potential generated by gravity.
Before proceeding, let us comment about the discrepancy with the result presented in [36] where a negative eigenvalue was found. The metric studied in [36] has the form ds 2 = N 2 (ρ)dρ 2 + R 2 (ρ)dΩ 2 , with the wormhole solution corresponding to N (ρ) = 1 and R(ρ) = R wh (ρ), where R wh (ρ) satisfies the equation The authors of [36] focused the analysis only on homogeneous perturbations; they defined the perturbed metric element as ds 2 = [1 + n(ρ)] 2 dρ 2 + [R wh (ρ) + r(ρ)] 2 , and -contrary to the explicit gauge-invariant formulation in this work -they fixed the gauge with the choice n(ρ) = 0. This procedure may lead to incorrect results since perturbations are truncated before the gauge fixing, thus preventing from the possibility to fix all possible gauge degrees of freedom. Furthermore, in [36] the conformal factor problem was circumvented by using the Gibbons-Hawking-Perry rotation [37]. In synthesis, the 'wrong' negative kinetic term in the quadratic action for the perturbation r(ρ) changes its sign as a consequence of the replacement r → ır, with the new r real. As already mentioned above, this is an ad hoc prescription without a clear physical interpretation, and its naïve application may lead to misleading results [38].

Effective potential from Euclidean wormholes
The non-linearly realized global U (1) symmetry φ → φ + αf a is generated by the axion charge Q. The ground state -corresponding to the bottom of the mexican hat potential -is degenerate along the angular direction, and vacua are described by continuously connected field configurations with minimum energy.
Half-wormholes induce quantum tunneling between classical vacuum states with different axion charge Q. The situation is schematically represented in fig. 6. As Euclidean time passes by, in the presence of an instanton (anti-instanton) an observer on R 3 experiences a change ∆Q = −n (∆Q = +n) since there is a net flux of axion charge equal to +n (−n) through the wormhole throat [16,17,39]. Eventually, we will consider Figure 6: Half-wormholes as instantons in Euclidean space. In the left (right) panel we show an instanton (anti-instanton) mediating the topology change R 3 → R 3 ⊕ S 3 , and carrying away an axion charge +n (−n). The axion charge is overall conserved but an observer on R 3 experiences, as Euclidean time passes by, a change ∆Q = −n (∆Q = +n). The non-conservation of the axion charge Q on R 3 implies an explicit breaking of the symmetry it generates.
the case with n = 1, since the instanton action in eq. (19) is minimized. Consequently, from the point of view of the observer on R 3 , the axion charge is not conserved, and the associated symmetry explicitly broken down to the discrete gauge symmetry φ → φ + 2kπf a , with k ∈ Z, which remains intact in the presence of wormhole instantons, as discussed in section 2.1.
In the presence of wormhole instantons the situation is similar to that of a particle moving in a onedimensional periodic potential. There are infinite classical minima at x = j, and each minimum corresponds to a degenerate ground state |j , as shown in fig. 7. Instantons can begin at any initial position, x = j, and go to the next one, x = j + 1. As a result of these tunneling transitions, the true vacuum of the system is a superposition of the degenerate ground states |j . It is instructive to work out this analogy in more detail. In particular, we can compute the probability for the tunneling process |j → |k summing over all the possible instanton and anti-instanton configurations. We find [22] where ω ≡ V (0), K is the determinant factor describing quantum fluctuations around the classical instanton trajectory. Notice that for a single instanton/anti-instanton path with action S 0 the contribution to the transition amplitude is [22] ω π e −ωT /2 Ke −S0 .
Eq. (58) include multi-instanton/anti-instanton solutions using the dilute-gas approximation. According to this approximation, strings of widely separated instantons and anti-instantons centered at τ 1 , . . . , τ n , and satisfying the condition −T /2 < τ 1 < · · · < τ n < T /2, are distributed arbitrarily along the time direction. Multi-instantons are not exact classical solutions, but they represent the leading term for the tunneling amplitude between distant wells. The contribution of a multi-instanton solution consisting in n well-separated objects takes the same form as in eq. (59) but with K → K n , S 0 → nS 0 . A similar result, with n substituted byn, is valid for an anti-instanton string. In addition, the integration over the freely-distributed position gives the factor The Kronecker delta in eq. (58) takes into account the fact that the total number of instantons minus the total number of anti-instantons must equal the change in x between the initial and final position eigenstates.
Using the Fourier series representation we find From eq. (62), we read the energy eigenstates (Bloch waves, using the language of solid state systems) and eigenvalues The periodic potential contains a translation symmetry -analogue in rôle to the gauge symmetry θ → θ + 2kπ in eq. (2), left unbroken by wormhole instantons -that forces the eigenstates to be shift-invariant. It is indeed possible to introduce an operator T that generates an elementary translation T |j = |j + 1 . T commutes with the Hamiltonian, and both operators can be diagonalized simultaneously. Each energy eigenstate in eq. (63) is characterized by an angle ζ, eigenvalue of T . It is indeed immediate to check that the states |ζ are eigenstates of T , with T |ζ = e ıζ |ζ . A rigorous treatment in the context of wormhole physics was proposed in [39][40][41] (see also [42][43][44]). The Figure 8: The computation of the transition amplitude between a state |n i with n i instantons to a state |n f with n f instantons requires the sum over all possible four-geometries. We show here the most general one [40,41] which includes m non-interacting wormholes, n wormholes that are emitted and reabsorbed, n i −m half-wormholes in the initial and n f −m half-wormholes in the final state. From a local perspective, an effective potential on the space-time volume M is generated (see text for details).
goal of these papers was to compute the influence of wormholes on physics at low energy, or at distances large compared to the thickness of the wormhole throat, L. The key idea is that at energy scales below L −1 wormholes can be integrated out. 8 In the resulting low-energy effective theory, the Lagrangian density takes the form The second term captures the effect of topological fluctuations due to wormhole physics. In full generality, the sum over i represents the possible presence of instantons of different type (for instance, instantons with 8 Notice that the validity of this approach relies on the fact that large wormholes, with size L > L, can be safely neglected. Naïvely, by looking at the action in eq. (19), one would conclude that this is indeed the case: The action of a wormhole of size L scales as L 2 M 2 Pl , and larger wormholes will be strongly suppressed by e −S . It is important to remark that this conclusion might not be so obvious. Wormholes received in the past a lot of attention in connection with a possible solution of the cosmological constant problem [45,46]. The major objection to this proposal, known precisely as 'large wormhole problem', refers to the fact that large wormholes are also characterized by large densities in the Euclidean space, and, despite the suppression provided by their action, they may lead to strong non-local interactions over arbitrarily large distances, a prediction of fatal consequences [47]. A possible solution to this obstruction was proposed in [48,49]. The point is that if the wormhole solutions carry a non-zero value of some conserved charge -as in the case we are considering in this work -large wormholes are destabilized by small ones as a consequence of the charge non-conserving interactions introduced by the latter (see eq. (68) below), and, as a result, only small wormholes survive. This argument was refined in [50] in the axion case.
. . ] are generic functions of fields and their derivatives, and A i are combinations of creation and annihilation operators describing emission and absorption of half-wormhole geometries of type i. Notice that A i do not depend on space-time position since wormholes do not carry any momentum. Furthermore, wormholes do not carry away any quantity coupled to gauge fields, and O i must be a Lorentz scalar, singlet under charge and color. However, as we learned in section 2, wormholes carry off axion charges, and as a consequence we expect the operator O i (x) to transform non-trivially under the global U (1) symmetry. We shall return on this point later.
Before proceeding, let us give a more quantitative understanding. The fact that integration over wormhole geometries gives rise to the structure in eq. (64) can be understood with an explicit computation [40,41]. In practice, one can compute the probability for the transition between a state |n i with n i instantons to a state |n f with n f instantons, n f |e −HT |n i . For simplicity we do not distinguish here between instantons and anti-instantons, and we only consider instantons of unit charge. The computation should include a sum over all possible four-geometries, and in fig. 8 we show the most general of such configurations [40,41]. This configuration includes m disconnected wormholes that do not interact with the Euclidean space-time, and n wormholes that are emitted and reabsorbed (see also fig. 2). The amplitude for this geometry is weighted by the factor e −Sinst(2n+n f +ni−2m) , where S inst is the instanton action (with unit charge) computed in eq. (35). Furthermore, we indicate with K √ gd 4 x the amplitude for inserting a single wormhole end in an infinitesimal volume of M. For simplicity, we treat K as a constant while, in general, it contains combinations of fields defined on M. The final result -in analogy with eq. (58) -is where the volume factor V T = M d 4 x √ g comes from an integration over the location of each instanton in the dilute gas approximation. The crucial observation in [40,41] is that the same result presented in eq. (65) can be directly obtained from the left-hand side of eq. (65) using the Hamiltonian H = Ke −Sinst V (a + a † ), after introducing creation and annihilation operators a † and a subject to the commutation relation [a, a † ] = 1, and defined by |n = 1/n!(a † ) n |0 . This result shows the validity of the assumption made in eq. (64). 9 The effective wormhole action -writing explicitly the presence of field operators O i (x) in the spirit of eq. (64) -is therefore [39] where a † q + a −q (a † −q + a q ) describes the creation of a half-wormhole with charge q (−q) or, equivalently, the annihilation of a half-wormhole with charge −q (q), as illustrated in fig. 9. The operators a q and a † q obey the usual commutation relations [a q , a q ] = [a † q , a † q ] = 0 and [a q , a † q ] = δ qq . As anticipated before, the operator O q (x) transforms non-trivially under the global U (1) symmetry, and it is possible to . As shown in [39], this transformation property ensures conservation of the global U (1) charge when the sum over all possible topologies is considered (that is, for instance, R 3 ⊕ S 3 for the creation of a single instanton, see fig. 6).
Finally, in order to derive a viable effective action, it is crucial to understand the rôle of the creation and annihilation operators in eq. (66). To this end, we can apply the lesson learned from Quantum Mechanics. At the beginning of section 3, we discussed the double-well potential. In perturbation theory, the system has two degenerate vacua | ± k , and, after including instanton solutions, the Hamiltonian is diagonalized by | + k ± | − k , with the true vacuum corresponding to the antisymmetric combination. This example shows that in the presence of instanton with no negative modes the quantum vacuum is constructed as a coherent superposition of classical vacua. The same conclusion remains valid for the periodic potential studied at the beginning of this section, with the true vacuum defined by the superposition in eq. (63). A similar situation arises in Quantum Field Theory. The QCD Lagrangian has a discrete set of degenerate classical minima, labelled by integers n -dubbed winding number -and indicated with |n . Indeed, considering the Chern-Simons current K CS µ and the corresponding charge K CS it is possible to show that integer values K CS = n correspond to non-equivalent pure gauge field configurations with zero energy. Yang-Mills instantons tunnel between two such classical vacua, resulting in a non-zero matrix element n + 1|e −HT |n . It means that in the presence of instantons the |n vacua do not diagonalize the Hamiltonian, in analogy with eq. (62) and eq. (65). On the contrary, the true vacuum is the coherent superposition |θ = n e ınθ |n . We expect the same situation in the presence of gravitational instantons.
To this end, we introduce the operators C q ≡ a † q + a −q and C † q ≡ a † −q + a q with commutation relations It is possible to simultaneously diagonalize the operators C q and C † q , and we define the eigenvalue equations C q |α = α q |α , C † q |α = α * q |α . Tunneling transitions bring the system in the coherent state (a † −q + a q )|α = α q e ıδq |α , defined in analogy with the harmonic oscillator as the eigenstate of the annihilation operator with complex eigenvalueα q ≡ α q e ıδq . We can therefore replace creation and annihilation operators in eq. (66) with their eigenvalues, and obtain Effective potential from Euclidean wormholes where in the last line we explicitly wrote the most general combination of operators singlet under the U (1) global symmetry. In eq. (69), we used as mass suppression scale the cutoff 1/L, and a, b and α q are expected to be order one numbers. In eq. (68) the crucial point is that breaking effects are always proportional to the factor e −Sinst , a trademark representing their non-perturbative origin. We remark that this point seems to be quite often overlooked in the literature, especially as far as phenomenological applications are concerned.
Very often breaking effects generated by gravity are treated in a naïve way, along the line of the discussion outlined in the introduction, see eq. (4). On the contrary, the presence of the suppression factor e −Sinst plays a crucial rôle. In the second part of this paper we shall discuss in detail the most relevant phenomenological applications.

Phenomenological implications
In this section we discuss several phenomenological situations characterized by the presence of light axions with a decay constant such that non-perturbative gravity effects become relevant.

The QCD axion
We start discussing the implication of eq. (68) for the QCD axion. We consider only wormholes with unit charge, since they give the dominant contribution. The effective axion potential, including both QCD and non-perturbative gravitational contributions, is schematically given by where we estimated the prefactor of the gravitational contribution to be of order (1/L) 4 . This is nothing but an order-of-magnitude estimate based on dimensional arguments but it does not affect the final result since the exponential term e −Sinst dominates. In this regard, it is worth emphasizing the size of this suppression.
Considering the benchmark value f a = 10 10 GeV for a QCD axion in the classical window, we have S inst 1.7 × 10 8 . The explicit breaking generated by gravity is therefore negligible for all purposes. Let us explore the implications of eq. (70) in more detail. First, notice that we defined the axion field so that the low energy QCD contribution to the axion potential is minimized at φ = 0, and we redefined accordingly the phase δ 1 in eq. (68). The phase shift δ 1 arises from a mismatch between the gravitational and the low energy QCD term. In the absence of CP violation in the gravitational sector, the phase δ 1 is just proportional to Arg Det[Y u Y d ], with Y u,d the Yukawa matrices for up-and down-type quarks, and it arises from the chiral rotation needed to move the phase of the fermion mass matrix into the θ-term. 10 In the 10 The possibility to have CP violation from gravitational effects is a subtle question, and, to the best of our knowledge, the situation is the following. Gravity is CP-conserving at the perturbative level, and the only source of CP violation may arise in connection with non-perturbative physics. In the Standard Model plus gravity, CP-violating effects are related to the three terms whereRµνρσ = µναβ R αβ ρσ /2. Notice that, in addition to the standard Ricci term of the Einstein-Hilbert action, one needs to add the topological contribution proportional to the CP-odd combination RµνρσR µνρσ . As well known, these three terms are total derivatives, and they do not contribute to the equation of motion. However, in QCD θ QCD becomes a fundamental physical parameter at the non-perturbative level. The reason is ultimately related to the fact that the third homotopy group of SU (3) C is non-trivial, since π 3 [SU (N )] = Z, and non-equivalent pure gauge filed configurations with zero energy fall into topologically distinct classes. QED in Minkowski space does not possess this property since π 3 [U (1)] = 0. However, in the presence of gravity this is not true anymore. Eq. (71) admits non-perturbative gravitational instanton solutions, dubbed Eguchi-Hanson instantons [51,52], that provide a non-trivial background metric in which θ EM becomes physical. The Eguchi-Hanson gravitational instantons are not related to wormhole physics, and in general their contribution is suppressed by the large action S EH = P 2 π/α EM , where P is the electric charge of the instanton [53]. In this paper we do not investigate such non-perturbative solutions. Recently, in [54][55][56][57][58] CP-violating effects related to RR were explored using the language of differential forms. In the three-form formulation the strong CP problem is equivalent to the dynamical generation of a mass gap for the Chern-Simons three-form of QCD, C. This is exactly what the axion solution does, since it provides a pseudo-scalar degrees of freedom that is eaten up by C which in turn becomes massive. In this picture, the only way to re-introduce the strong CP problem in the presence of the axion is to re-establish a massless pole in the propagator of C. This can be accomplished by introducing an additional massless three-form, C G , since in this case the axion will be able to give a mass only to a linear combination of C and C G . Gravity naturally provides a candidate for C G , that can be identified with the gravitational Chern-Simons three-form whose field strength equals RR [54].
following, we interpret δ 1 as an arbitrary order one phase. In the regions shaded in red non-perturbative gravity dominates. In the right panel, the dashed blue line represents the contribution to the axion mass generated by QCD, computed according to [59].
The effect of the gravitational correction is twofold: It shifts the axion mass, and produces a non-zero θ eff . We find where we remind the reader that the action scales with M Pl /f a . The gravitational correction to the axion mass features a very interesting property if compared with the contribution generated by QCD. As f a increases, the QCD axion mass Λ 4 QCD /f 2 a decreases but the gravitational term becomes more and more important, and eventually it overcomes the former. This means that we expect a lower bound on the axion mass. This is shown in the right panel of fig. 10 where we shaded in red the region in which non-perturbative gravitational corrections dominate. We find the lower bound on the axion mass (or, equivalently, the upper bound on the axion decay constant) Lower bound on the QCD axion mass m a 4.8 × 10 −10 eV , f a 10 16 GeV .
In the left panel of fig. 10 we show the impact of the non-zero θ eff in eq. (73). For reference, we consider the experimental bound θ eff < 10 −10 (horizontal gray line). We find that if f a 10 16 GeV then the contribution induced by gravity becomes too large. Notice that, by accident, the bounds on f a extracted from m a and θ eff are of the same order, and assuming δ 1 = 0 does not change quantitatively our conclusion.
In fig. 11 we show the impact of the bound in eq. (74) on the QCD axion parameter space. As customary, the condition f a ≶ H I /2π, with H I the Hubble expansion rate at the end of inflation distinguishes between the cases in which the PQ symmetry is broken before (>) or after (<) the end of inflation. The first possibility defines the so-called anthropic axion window, in which the possibility to reproduce the observed dark matter relic density Ω a h 2 0.1 relies on a fine-tuned choice of the initial misalignment angle θ in (black dot-dashed lines in fig. 11 for different θ in ). Notice that a substantial part of the anthropic axion window is excluded by the presence of axion isocurvature perturbations in the early Universe (blue region with vertical meshes) that we Figure 11: QCD axion parameter space (see text for details). The region shaded in red with horizontal meshes is excluded by non-perturbative gravitational effects, eq. (74).
compute following [60,61] (see also [62,63] for recent progress for non-standard cosmologies and generalization to axion-like particles). The region f a < H I /2π defines the classical axion window. In the region (i) shaded in brown with diagonal meshes Ω a h 2 > 0.1, and the Universe is over-closed. The dashed brown lines reproduce 1% and 10% of the observed dark matter abundance. In this region a sizable contribution from decay of topological defects is expected [64]. For completeness, we also show the region excluded by white dwarf cooling time [65] and upper limit on the tensor-to-scalar ratio [66]. The bound in eq. (74) applies in the anthropic axion window, and disfavors a region of the parameter space (shaded in red with horizontal meshes) previously allowed. We remark that the bound in eq. (74) was derived without specific assumption but only considering minimal coupling of the axion field with Einstein gravity.
We now move to briefly discuss few scenarios of phenomenological relevance in which non-perturbative gravitational corrections to the mass of the QCD axion may play an important rôle.

QCD axions and black hole superradiance
Spin and mass measurements of stellar-size black holes exclude the QCD axion mass window [67] 6 × 10 −13 m a [eV] 2 × 10 −11 , corresponding to 3 × 10 17 f a [GeV] 10 19 . This is because superradiance effects become efficient [68,69] when the Compton wavelength of the axion is comparable with the horizon size of the black hole. The axionic field forms a quasi-stationary configuration around the black hole at the expense of its rotational energy, giving birth to a quasi-bound system that shares remarkable similarities -such as energy orbitals and level transitions -with the hydrogen atom. From λ Compton = h/m a c ∼ R, we have where R is the typical radius of a stellar-mass black hole, thus justifying the mass range excluded in eq. (75). At the two sides of this interval, and within the mass range m a = [10 −14 , 10 −10 ] eV, stellar black hole superradiance in the presence of the QCD axion may produce in the next few years spectacular signatures -both direct and indirect -in gravitational wave detectors such as Advanced LIGO [70,71]. Indirect signatures refer to the observation of gaps in the spin-mass distribution of final state black holes produced by binary black hole mergers. Direct signatures refer to monochromatic gravitational wave signals produced during the dissipation of the scalar condensate after the superradiant condition is saturated. As far as direct signatures are concerned, a careful assess of the detection prospects in Advanced LIGO and LISA was recently proposed in [72,73]. The outcome of the analysis is that, considering optimistic astrophysical models for black hole populations, the gravitational wave signal produced by superradiant clouds of scalar bosons with mass in the range m a = [2 × 10 −13 , 10 −12 ] eV is observable -i.e. it is characterized by a signalto-noise ratio larger than the experimental threshold -by Advanced LIGO. Notice that this region seems to be ruled out if one considers at face value the bound in eq. (75). However, it is worth emphasizing that the bound in eq. (75) is most likely only indicative since it is based on black hole spin measurements that are extracted indirectly from X-ray observations of accretion disks in X-ray binaries. We only have very few of such measurements at our disposal, and it is difficult to extract a bound with significant statistical confidence. The discussion of this matter however escapes the scope of the present work.
As clear from the right panel of fig. 10 and from the parameter space in fig. 11, the QCD axion in the mass range m a = [10 −14 , 10 −10 ] eV violates the bound in eq. (74), and fits into a region where nonperturbative gravitational effects dominate over QCD. We therefore conclude that -working under the very same hypothesis, that is an axion minimally coupled to Einstein gravity -the phenomenologically interesting mass range m a = [10 −14 , 10 −10 ] eV motivated by black hole superradiance is theoretically forbidden for the QCD axion. This is not, of course, a lapidary conclusion. To be more optimistic, observing the QCD axion in connection with black hole mergers at the Advanced LIGO could imply an evidence for modifications of Einstein gravity. Alternatively, a scalar boson different from the QCD axion could still leave its imprint in the texture of gravitational wave signatures. In this last case, the mass term generated by gravity can be used  fig. 12. Clearly, if f a 10 16 GeV it is possible to span, thanks to the exponential factor e −Sinst a large range of allowed values, including the regions favored by the analysis in [72,73].

Axion stars as black hole seeds
The idea is to consider axion stars as black hole seeds, which are supermassive in the case of ultralight axions [74]. By studying numerically the collapse of axion stars in the context of general relativity, the authors of [74] identified the critical value of the axion decay constant f TP ∼ 0.06 M Pl above which black hole formation occurs. The mass of the typical black hole formed from axion star collapse is M BH ∼ 3.4 (f a /0.12 M Pl ) 1.2 M .
As before, the condition f a f TP ∼ 0.06 M Pl is not compatible with the bound in eq. (74). We stress that this bound was derived considering a QCD axion minimally coupled to Einstein gravity, therefore without any additional assumptions compared to the ones in [74]. Said differently, it would be interesting to include in the analysis of [74] non-perturbative gravitational corrections to the axion potential since they play a relevant rôle in the region of parameter space in which black hole formation may occur.

Ultralight scalars as cosmological dark matter
Dark matter must behave sufficiently classically to be confined on galaxy scales. If we suppose dark matter to be a boson with mass m a and velocity v, we can require to behave classically down to the typical size of Milky Way satellite galaxies, and obtain the condition For a typical halo with size R ∼ kpc, and mass M ∼ 10 9 M , we expect a virial velocity v ∼ G N M/R ∼ 70 km/s. From eq. (77), we extract a lower limit for the mass of bosonic dark matter of about 10 −22 eV. The dark matter saturating this value is known as fuzzy dark matter (FDM) [75]. This kind of ultralight dark matter could form a Bose-Einstein condensate on galactic scales, providing a possible solution to the tensions that arise when the standard cold dark matter paradigm is probed into the deep non-linear regime at redshift z ∼ 0 [76][77][78]. Apart from this motivation, an ultralight dark matter condensate features peculiar observational astrophysical properties, and it catalyzed increasing attention in the dark matter community (see, e.g., [79][80][81][82][83][84][85][86]). From a particle physics perspective, the extreme lightness of the FDM is well suited by an axion-like particle. Furthermore, reproducing the observed value of relic abundance via the misalignment mechanism gives a clue about the value of its decay constant f a .
For completeness, let us quickly sketch the computation. We consider here an axion-like field a(t, x) = f a θ(t, x) with potential V (a) = µ 4 (1 − cos θ). In the Friedmann-Robertson-Walker cosmology ds 2 = −dt 2 + R 2 (t)d x 2 the evolution of the axion field is given byä(t, x) + 3Hȧ(t, x) − a(t, x)/R 2 + dV (a)/da = 0, where H ≡Ṙ/R, and the dot indicates derivative w.r.t. time. Neglecting higher order in the potential, we havë a(t, x) + 3Hȧ(t, x) − a(t, x)/R 2 + m 2 a a(t, x) = 0, with m 2 a ≡ µ 4 /f 2 a . As customary, we can define the time t 1 at which the condition m a = 3H(t 1 ) is satisfied. Notice that for simplicity we are considering an axion mass that is, to a first approximation, temperature-independent. The time t 1 separates two regimes.
For t < t 1 , we can neglect the mass term, H m a . Introducing the Fourier decomposition a(t, k) = d 3 xe i k· x a(t, x) we obtainä(t, k) + 3Hȧ(t, k) − (k 2 /R 2 )a(t, k) = 0. The Fourier modes a(t, k) separate into modes outside (that is for k/R H) and inside (that is for k/R H) the horizon. Inside the horizon, we can not neglect the k 2 /R 2 term. Solving the equation of motion for these modes, it is possible to show that the corresponding solutions oscillate with frequency k/R, and the amplitude decreases with time as 1/R. The modes that are confined inside the horizons until t t 1 can therefore be neglected. As far as the modes outside the horizon are concerned, we can neglect the k 2 /R 2 term, and the equation of motion is solved by a(t, k) = c 1 ( k) + c 2 ( k)t −1/2 . The modes that are confined outside the horizon until t t 1 are therefore frozen at some constant initial value, and they are collectively called zero modes a 0 .
Let us now move to discuss the second regime, t t 1 , focusing on the zero modes. When the axion mass becomes non-negligible, the evolution of the zero modes is described by the equationä 0 + 3Hȧ 0 + m 2 a a 0 = 0. The zero modes, previously frozen, at t t 1 start oscillate with frequency set by m a . From the definition of axion energy density 2ρ a0 =ȧ 2 0 +m 2 a a 2 0 we find, using the equation of motion,ρ a0 = −3Hȧ 2 0 . Since a 0 oscillates with period m a , we can average over an oscillation, and estimateȧ 0 m a a 0 . Consequently, ρ a0 m 2 a a 2 0 anḋ ρ a0 −3Hρ a0 . After integration, the latter equation gives the scaling ρ a0 (t) ∼ m a /R 3 (t). This is the key argument defining the physics of the misalignment mechanism: During the cosmological evolution of the axion field, for t t 1 the number of axions in a comoving volume is conserved. We can therefore write the energy density stored in the axion zero modes today as where the dependence from the so-called initial misalignment angle θ in is manifest. The rest of the computation makes use of some basic thermodynamic concepts to rephrase eq. (78) in terms of observable quantities. From entropy conservation, g s * (T )R 3 T 3 = const, we have where h ≡ H 0 /(100 km/s/Mpc), with H 0 = 67.8 km/s/Mpc [87], is the reduced Hubble constant, T 1 is the temperature at time t 1 , T 0 = 2.726 • K is the present temperature of the Universe, ρ crit = 3H 2 0 /8πG N is the critical energy density, g * (T ) (g s * (T )) is the effective number of degrees of freedom (effective number of degrees of freedom in entropy) at temperature T , and ρ R (T 0 ) is the present energy density in radiation. In eq. (79) s(T 0 ) is the entropy density at present time, and the relation between s(T 0 ) and ρ R (T 0 ) follows from The energy density in radiation can be estimated from the two Planck measurements [87] of the current matter density Ω m and the redshift of radiation-matter equality z eq . From ρ m (z eq ) = ρ R (z eq ), which is equivalent to Ω m (z eq + 1) 3 = Ω R (z eq + 1) 4 , and using Ω m = 0.3, z eq = 3365, it follows that Ω R h 2 4.3 × 10 −5 .
The value of T 1 follows from the condition m a = 3H(t 1 ) taking into account that we can write the Hubble parameter as a function of the temperature as H(T ) = π 2 g * (T )/90M 2 Pl T 2 . Parametrically, we have the relation M Pl m a T 2 1 . For a FDM candidate with m a 10 −22 eV, it follows that T 1 10 3 eV. At such small value of the temperature, we can approximate g * (T 1 ) ≈ 1. Furthermore, we also have g * (T 0 ) g s * (T 0 ). All in all, defining Ω a ≡ ρ a0 (t 0 )/ρ crit , and assuming θ in 1, we find the parametric order-of-magnitude scaling Eq. (81) points towards a very specific interplay between f a and m a , and it is interesting to further investigate its origin. In [85] this relation was justified referring to string theory models in which one expects an explicit breaking of the axionic shift symmetry due to the existence of worldsheet or membrane instantons [88]. Even without invoking string theory constructions, we point out that the mass term generated by gravity in eq. (68) provides a nice explanation for the relic abundance in eq. (81). Indeed, using m 2 a (1/L) 4 e −Sinst /f 2 a , the observed abundance can be reproduced with f a 8 × 10 15 GeV and m a 2.5 × 10 −18 eV. Notice that this value falls within the mass range that will be explored by the LISA gravitational wave interferometer, see fig. 12. This is an interesting observation, since it shows that a dark matter candidate with nothing more than gravitational interactions may explain in a natural way the lightness of its mass.
We finally note that, as fig. 12 shows, a value of the axion decay constant above f a 10 16 GeV would generate an unacceptably large contribution in eq. (81).

Comments on the rôle of possible UV completions
In this section we discuss possible UV completions of the axion Goldstone mode, and their consequences on wormhole physics.

Dynamical radial mode
We start our investigation considering the simplest UV completion of the axion theory. This is the case of a U (1) global symmetry that is spontaneously broken by the VEV |Φ| = f a / √ 2 of a complex scalar field Φ(x) = [f (x)/ √ 2]e ıφ(x)/fa . The mexican hat potential responsible for the spontaneous symmetry breaking is The radial field f (x) gets a mass M 2 f = λf 2 a . The physics discussed in the previous sections corresponds to the case in which the radial mode does not participate in the dynamics, and its value remains frozen at f = f a . We expect this situation to be strictly true if M f > 1/L. For the critical value f a 10 16 GeV that we identified in section 5.1, we have 1/L 10 18 GeV. It is therefore hard to imagine that the radial field f plays no dynamical rôle, since we expect its mass -for reasonable small coupling λ -to lie below the scale set by the wormhole throat. In order to validate our conclusions, it is necessary to generalize the wormhole solutions, and compute their action, to the case in which f is a dynamical field [15,95]. In order to make contact with the existing literature, we look for spherically symmetric solutions, and we use the Euclidean metric ansatz ds 2 = dr 2 + R 2 (r)dΩ 2 3,1 . The Ricci scalar is R = −6[−1 + (R ) 2 + RR ]/R 2 , and the Einstein field equations, together with the equation of motion for the radial field f , are [17,95] f where we already substituted f 2 (θ ) 2 = n 2 /4π 4 f 2 R 6 , with n the quantized axion charge along the wormhole throat. Instead of eqs. (82,83), in our numerical investigation we found more useful to use the combination We introduce the dimensionless variables [15,95] ρ ≡ rM Pl 3λ 8π and the differential equations become with Q ≡ n 2 λ 2 /8π 4 . Finally, in terms of the dimensionless variables defined in eq. (85), the action -including the Gibbons-Hawking-York boundary term -reads We solve eq. (86) by means of a shooting method that we validate against the results of [15]. The procedure goes as follows. First, for a given initial guess for F (ρ = 0) we compute A(0) from the equation Second, we solve the system in eq. (86) with the four boundary data F (0), A(0), A (0) = 0, F (0) = 0. Finally, we check the asymptotic condition F (∞) = F a , and we tune the initial guess F (0) until we find it satisfied. We show our result in fig. 13, where we focused on the critical value f a = 10 16 GeV and wormhole with unit charge n = 1. In the left (right) panel we show the dimensionless metric function A(ρ) (the dimensionless radial field F (ρ)) as a function of the rescaled distance ρ for three different values of λ. As far as the geometry is concerned, we see that, asymptotically, the solution recovers the flat space A(ρ) = ρ while for small ρ deviations describing the wormhole geometry emerge. Going through the wormhole throat, the radial field F (ρ) substantially deviates w.r.t. the its asymptotic value F a . As expected, including the radial mode with mass M f < 1/L as dynamical degree of freedom alters the wormhole solution found for non-propagating f . It is therefore important to quantify such deviation by computing the value of  Figure  the action in eq. (87). We show our result in the left panel of fig. 14, in which we plot the on-shell action S as a function of the mass of the radial mode M f . From eq. (35), we have that the instanton action with the radial mode frozen at f a = 10 16 GeV is S inst 170. Including the dynamics of the radial mode results in a net decrease of the action, in agreement with the result of [15]. Let us give a qualitative understanding of this effect. As derived in eq. (19), the wormhole action is proportional to M 2 Pl L 2 , and it decreases if the size of its throat gets smaller. The inclusion of a light dynamical radial mode has precisely this effect, as can be seen from the left panel in fig. 14: Going towards smaller values of λ the size of the wormhole shrinks, and, consequently, the wormhole action decreases. For illustrative purposes, in fig. 14 we extend the plot to large values of λ in order to show that, as M f moves towards the threshold 1/L above which it can be safely integrated out, the action increases.
The inclusion of the radial mode in the U (1) model shows that the presence of dynamical degrees of freedom belonging to the UV completion of the axion Lagrangian does not weaken the strength of the nonperturbative global symmetry breaking induced by gravity. On the contrary, the total action in fig. 14 is significantly smaller if compared to the one in which the radial mode does not fluctuate around its VEV. We therefore conclude that the bound derived in section 5.1 presumably represents a conservative estimate, since extra dynamics related to the UV completion of the axion Lagrangian may even strengthen the global symmetry breaking.

Higher-curvature operators and string theory
In this section we explore the impact of higher-curvature operators on the wormhole solution. To be more concrete, we include the Gauss-Bonnet term that represents the 4 th order in the derivative expansion of the gravitational action. We limit our analysis to the Gauss-Bonnet term since it is the unique ghost-free quadratic curvature invariant [96]. Furthermore, as we shall discuss later, the presence of the effective higher-curvature correction in eq. (88) is a peculiar prediction of many string theories -bosonic [96], heterotic [97], and type-I string [98] (while it vanishes in type-II superstring theory [99]). In D = 4 dimensions, the Gauss-Bonnet term reduces to a total derivative. As a consequence, it does not enter in the equation of motion but it gives a non-vanishing topological contribution to the wormhole action. Before proceeding, we notice that the validity of the derivative expansion in the gravitational action imposes the cutoff Λ ∼ M Pl / γ/4π. It is therefore natural to require the condition Λ > 1/L, in order to ensure the validity of the derivative expansion all the way down the wormhole throat. We find that this condition corresponds to To fix ideas, if f a = 10 16 GeV we have γ 1600. We now compute the value of the Gauss-Bonnet term on the wormhole solution. Using the metric in eq. (21), a direct computation gives and the action in eq. (88) gives Let us focus on the case of the QCD axion, discussed in section 5.1. At the critical value f a = 10 16 GeV, the instanton action is S inst 170. As a rule of thumb, we can say that if γ 170 the higher-curvature topological correction due to the Gauss-Bonnet term in eq. (91) significantly alters the instanton action, increasing its value, and making its impact on axion physics harmless. Although qualitatively correct, the real point here is to understand the physical implications of this numerical condition. On a pure gravitational ground, γ is nothing but a dimensionless number, and it seems difficult to justify the relation γ 170. Reinterpreting the effective action in eq. (88) having in mind the UV physics responsible for its generation would be much more useful. The lack of a quantum theory of gravity, however, makes any argument highly speculative. For lack of a better choice, we can ask string theory to help us with this task. In string theory the Gauss-Bonnet term often appears with the identification [100,101] where α = l 2 S is related to the string scale l S . We can gain more insight considering the case of heterotic string theory in which the string gauge coupling g S , the string mass scale M S , and the parameter α are related by [102] in close analogy with the relation discussed in eq. (20). From eq. (92), we find γ = 8π 2 /g 2 S . Notice that, presented in this form, the topological contribution to the action precisely matches the typical instanton action in the Yang-Mills case. We are now in the position to compare γ = 8π 2 /g 2 S and S inst . We show our result in the right panel of fig. 14, in which we compare the half-wormhole action S inst computed at the critical value f a = 10 16 GeV that we identified in section 5.1, with the action γ = 8π 2 /g 2 S as a function of the string coupling. On the top y-axis we put the values of the string mass scale M S , computed according to eq. (93). The region shaded in gray corresponds to the condition M S < 1/L.
Clearly, in the presence of a weak string coupling our computation can be completely invalidated. The string mass scale becomes lower than L −1 -thus breaking the effective description -and the topological action generated by string theory dominates over the gravitational instanton contribution.
However, we also stress that in the presence of a moderately strong string coupling the UV completion of GR cannot fix the problem since the gravitational instanton contribution dominates over the topological string term.
The lack of a quantum theory of gravity makes any speculation quite far-fetched, and the only intent of the plot in the right panel of fig. 14 is to provide a fair example in which the gravitational action based on Einstein gravity captures the relevant non-perturbative gravitational corrections.

Higher-dimensional operators
As discussed in the introduction, many authors proposed and tailored suitable extensions of the PQ symmetry in order to protect it from the presence of power-suppressed higher dimensional operators like those in eq. (4). Clearly, these attempts were motivated by a naïve understanding of the breaking effects generated by gravity. As revisited in this paper, gravity breaks global symmetries at the non-perturbative level, and the effective operators originated from this physics are always suppressed by the exponential of the wormhole action, as exemplified in eqs. (68,69).
However, it is conceivable that the U (1) PQ symmetry arises as an accidental global symmetry from a more fundamental theory. In this realization, because of its accidental nature, higher-dimensional operators may source an explicit breaking. In order to preserve the solution of the strong CP problem without introducing any degree of fine-tuning, it is necessary to protect the accidental symmetry from breaking effects induced by dangerous irrelevant operators. This can be achieved in a natural way by supporting the accidental PQ symmetry with a gauge symmetry. In this way, as we shall see, one can prevent the presence of symmetry breaking operators up to (in principle arbitrarily) high dimensions, thus obtaining a global PQ symmetry of very good quality, though accidental.
In this context, it is interesting to see what happens to non-perturbative breaking effects. In order to answer this question, let us discuss the specific construction put forward in [8]. The model features the presence of two sectors with two anomalous PQ symmetries, U (1) PQ and U (1) PQ , realized through the complex fields √ 2φ = f a e ıã/fa and √ 2φ = f b e ıb/f b . Notice that for simplicity we do not consider explicitly the radial components. The domain of the two phase fields areã/f a = [0, 2π) andb/f b = [0, 2π). The crucial observations are the following.
• It is possible to define a linear combination of U (1) PQ and U (1) PQ that is free from QCD anomaly.
Consequently, it can be promoted to a gauge symmetry, U (1) gauge PQ .
• In the limit in which the only interactions between the two sectors are those dictated by the aforementioned gauge symmetry, the theory possesses an accidental U (1) symmetry, and delivers a massless Goldstone field.
The gauge symmetry U (1) gauge PQ can be characterized as follows. The phases of φ and φ , with gauge charges q and q , transform under U (1) gauge PQ according to the shifts and from the kinetic Lagrangian we get where the mass of the gauge field is m 2 brings the Lagrangian into the neat form The field b is the would-be Goldstone boson eaten by the massive vector field A µ . The theory in eq. (97) enjoys an accidental (non-linearly realized) U (1) global symmetry under which the massless Goldstone field a remains shift-invariant. The Goldstone field a plays the rôle of the axion, and its continuos global shift symmetry is analogue to the one discussed in eq. (1). In addition, and again in close analogy with the discussion outlined in the introduction, a subgroup of this continuos shift symmetry is gauged, meaning that there is a phase redundancy in the definition of a. In order to visualize this property, we consider, following [8], the specific values q = 2 and q = 3. Notice that these two integers are relatively prime, meaning that their greatest common divisor is 1. This particular choice does not change the conclusion of the present discussion but it simplifies the formulas. In the right panel of fig. 15 we show in solid red the gauge orbits described in the field space (ã/f a ,b/f b ) by the gauge transformations in eq. (94). For instance, starting from the pointã/f a = 0,b/f b = 0 in field space ('pure gauge' configuration), by changing the value of the gauge parameter α one moves along the corresponding red line in the direction of the red arrow. At the pointã/f a = 4π/3,b/f b = 2π, because of the 2π-periodicity inb/f b , the gauge orbits reappears atã/f a = 4π/3,b/f b = 0. Similarly, because of the 2π-periodicity inã/f a , the gauge orbit connects the points in field spaceã/f a = 2π,b/f b = π andã/f a = 0, b/f b = π.
The axion field a = (q 2 f 2 a + q 2 f 2 b ) −1/2 (q f bã − qf ab ) is gauge invariant, and describes the direction orthogonal to the gauge orbits. The crucial points is that configurations in field space connected by gauge orbits are physically equivalent, since related by gauge transformations. As a consequence, by crossing the gauge orbits the axion field experiences a periodicity, and it is possible to show that the domain of a is given by [8,104] In parallel with eq. (2), we therefore conclude that the discrete symmetry represents a gauge symmetry of the axion field a. Finally, since the U (1) gauge PQ gauge symmetry must be anomaly free, the anomalous QCD coupling should involve the combination The advantage of this construction is the following. The accidental U (1) symmetry remains unbroken in the limit in which one only considers gauge interactions between the two sectors. More generally, additional interactions between the two sectors have no reason to respect the accidental U (1) symmetry but they must be invariant under the gauge symmetry U (1) gauge PQ . This means that operators of the form may source a breaking of the accidental U (1) symmetry. The suppression scale Λ may or may not coincide with the Planck scale. O 1 and O 2 are generic operators with mass-dimension d O1 and d O2 made out of fields in the two sectors, and they have equal and opposite U (1) gauge PQ charge in order to respect the gauge symmetry. The lowest dimensional U (1) symmetry breaking operator depends on the charge assignment under U (1) gauge PQ , and, in the context of explicit models [8], it is possible to show that judicious choices suppress U (1) breaking effects down to an acceptable level.
Having set the framework, we can now go back to the original motivation of this section, and explore the rôle and impact of non-perturbative wormhole solutions.
As noticed before, eq. (97) contains an accidental (non-linearly realized) U (1) global symmetry under which the massless Goldstone field a remains shift-invariant. The subgroup a → a + 2kπF a , k ∈ Z is a gauge symmetry. The conditions for the existence of Euclidean wormhole solutions discussed in section 2 are therefore fulfilled, and non-perturbative gravitational effects lead to an explicit breaking of the axionic shift symmetry generating a potential as in eq. (68). The wormhole action in eq. (35) is controlled by the decay constant F a that sets the periodicity of a. In conclusion, the outcome of this qualitative discussion is that the protection due to the gauge symmetry U (1) gauge PQ may suppress irrelevant operators like those in eq. (101) but cannot remove non-perturbative breaking terms. The idea behind this construction is to use gauge symmetry to protect the axion by means of a local discrete symmetry; however, since wormhole solutions already respect this symmetry, their presence is unaffected. We expect this argument to be fairly general.

Relaxation of the electroweak scale and the clockwork mechanism
The relaxation of the electroweak scale addresses the issue of naturalness from a cosmological perspective [89]. The idea that lies at the heart of the model is that the mass parameter in the Higgs potential V (|H|) = µ 2 H † H + λ(H † H) 2 changes during inflation as a function of the classical value of a slow-rolling axion-like scalar field, the relaxion φ. To achieve a large separation of scales between the electroweak scale v and the cutoff M of the theory one needs a compact field space of size 2πF , with F > M/g. The condition M v implies g 1, and, for reasonably large values of M , one comes up against the problem of having trans-Planckian field excursions, typically F > M/g > M Pl [90].
An elegant solution to this problem, dubbed clockwork mechanism, was proposed in [91,92]. In a nutshell, the clockwork is a theory where a U (1) N +1 global symmetry is explicitly broken in such a way to preserve a single U (1) φ symmetry under which the original U (1) i=0,...,N fields φ i rotate synchronously, δφ i = qδφ i+1 , where q is the clockwork factor. This construction delivers a single Goldstone boson with large field range f eff ≈ q N f f , where f is the symmetry breaking scale of the U (1) N +1 symmetry. As already noticed in [90], the effective operators generated by gravity in eq. (68) with F > M Pl would spoil the simplest realization of the relaxation mechanism by introducing large corrections to the slow-roll potential of φ. It is therefore interesting to address the following question: Is the clockwork mechanism protected against non-perturbative gravitational corrections? To answer this question we can start from a simple realization of the mechanism with only two axions [91]. The Lagrangian density of the model is where the two shift symmetries φ i=1,2 /f → φ i=1,2 /f i=1,2 + c i=1,2 are explicitly broken by the potential V (φ 1 , φ 2 ) down to the residual U (1) φ symmetry under which the two axion fields transform according to with c arbitrary continuos parameter. In full generality, we are considering different (but comparable , as one can see by diagonalizing the mass matrix The canonically normalized massive direction is The crucial point is to identify the periodicity of the flat direction φ. The length of a periodic flat direction is determined by the minimal shift ∆φ under which the field configuration φ comes back to its original value. We can visualize the periodicity of the flat direction φ by means of a simple plot, as shown in the left panel Figure 15: Left panel. The red lines follow the U (1) φ orbits described by the symmetry transformation in eq. (103) with q = 5. For illustrative purposes, we also show the direction in field space corresponding to the axion φ and the massive mode φ H . Right panel. The red lines follow the gauge orbits described by the symmetry transformation in eq. (94). The axion direction is orthogonal to the gauge orbits (black arrow, see section 6.3 for details).
of fig. 15. Let us start from the point in field space (φ 1 = 0, φ 2 = 0), with φ = 0, and follow the U (1) φ orbits (solid red lines) obtained from the symmetry transformation in eq. (103) by continuously changing the parameter c. For illustration, we take q = 5. Because of the periodicities φ i=1,2 ≡ φ i=1,2 + 2πf i=1,2 , the U (1) φ orbits are wrapped in the elementary domain φ 1 /f 1 = [0, 2π], −φ 2 /f 2 = [0, 2π]. Notice that φ H is invariant under the symmetry transformation in eq. (103), and it describes the direction orthogonal to the U (1) φ orbits. The axion direction, on the contrary, is aligned with the U (1) φ orbits (black arrows in fig. 15). The axion field comes back to the origin after a distance in field space which is enhanced by the winding number q if compared with the original periodicities 2πf 1 ∼ 2πf 2 .
In terms of mass eigenstates, the Lagrangian density takes the simple form with m 2 H defined by the non-zero eigenvalue in eq. (104), and where . . . represents higher order in the potential for φ H . The massless scalar field φ corresponding to the unbroken U (1) φ symmetry admits, in full analogy with the construction we put forward in section 2, wormhole solutions whose quantization condition is set by the periodicity 2πf eff . We remark that φ represents the only direction in field space that can accommodate wormhole solutions as presented here. To understand this point one might consider wormhole solutions for the fields φ 1,2 with decay constants f 1,2 f ef f defined by the Lagrangian in eq. (102) and see if the potential V (φ 1 , φ 2 ) can be treated as a small perturbation. However, it is possible to see 11 that gravitational instantons computed for the free field φ 1,2 do not minimize the potential V (φ 1 , φ 2 ). As a consequence of the latter, instead of being a small correction, the potential gives a large contribution to the action when integrating over large scales away from the throat, and this makes these solutions, on balance, irrelevant.
This conclusion can be generalized to the case of N + 1 axions. The Lagrangian density is The periodicity of the flat direction is given by [91,92] where in the equality we assumed for simplicity equal decay constant f i=0,...,N ≡ f . As before, it is possible to construct wormhole solutions in connection with the Goldstone boson φ associated with the unbroken non-linearly realized U (1) φ symmetry. The periodicity of φ, ∆φ = 2πf eff ≈ 2πq N f , is completely fixed by the clockwork construction, and it sets the quantization condition discussed in section 2.1. From this discussion, it is clear that the instanton action is controlled by the ratio M Pl /f eff , and that nonperturbative breaking effects prevent from the possibility to have trans-Planckian field excursions: Gravity breaks explicitly the U (1) direction left unbroken by the clockwork construction.
However, in order to provide a rock-solid conclusion, it is important to keep in mind -in close analogy with what discussed in section 6.1 -the possible rôle played by UV dynamics. For instance, the Lagrangian density in eq. (107) could be generated by the dynamics describing N +1 complex scalar fields Φ i = [(ρ i +f i )/ √ 2]e ıφi/fi subject to the potential with m i ∼ f i , λ i ∼ 1, after integrating out the radial modes. In this setup, one gets Λ i ≡ ( i f i f 3 i+1 /2) 1/4 . The situation in similar in spirit to the setup explored in section 6.1 with the additional complication given by the presence of more fields and explicit breaking terms controlled by the order parameters i . The results of section 6.1 suggest that the presence of extra fields does modify the geometry of the wormhole solution. However, this could result -as explicitly proved in section 6.1 considering the presence of a dynamical radial mode -in a decrease of the wormhole action, with a consequent tightening of the amount of symmetry breaking.
Within the present study, no firm conclusion can be established. Nevertheless, we argue that the possible presence of Euclidean wormhole solutions represents the most concrete threat against the possibility to engineer trans-Planckian field excursions via the clockwork mechanism. We think this is an interesting open question, and we leave it for future investigation.

Wormhole solutions for a generic Goldstone coset
Let us start with a lightning review of the treatment of Goldstone bosons from the coset G/H in the formalism of CCWZ [93,94]. The broken part of the Lie algebra will be denoted by the generators X A whereas the unbroken symmetry generators are T a with the following structure where the last equation follows for compact groups and we take the generators to be hermitian. The Goldstones parametrize the broken part of the group; one can write η = exp(iXθ) with η ∈ G/H in terms of n = dim(G)− dim(H) θ coordinates although different parametrizations are useful and often used. The transformation properties of η are with G ∈ G and h ∈ H. This expression follows from the result that any element of the group can be factorized into a broken times an unbroken element and the composition law Gη = η h(η). The basic building block is the Lie-dragged derivative which is a one-form ω living in the broken Lie algebra, or tangent space of the broken group manifold. Explicitly where ω A = ω A µ dx µ . In terms of the θ coordinates one has the vierbein and metric (we use the convention of summed repeated indices unless otherwise stated) In order to extract the symmetry transformations in terms of the θ coordinates we make eq. (111) infinitesimal and we project out the non-linear piece F (η) by restricting to broken generators only. We find where by X we denote projection onto the broken generators. The above result yields the Killing vectors in θ-space (δ θ i = A ξ i A ); there are n Killing vectors ξ A associated to the coset and dim(H) Killing vectors ξ a for the unbroken symmetries. The transformation of the one-form ω A under the unbroken part of the group is specially simple since the symmetry is linearly realized where we have used the second relation in eq. (110). In coordinate space the Killing vectors by construction satisfy The Noether currents for the broken generators read where we introduced the Goldstones decay constant f and equivalent expressions apply for the unbroken currents. We note that, if a coordinate system can be found for which ξ i A e B i = δ A B , the Noether currents read √ gω µ A and their conservation can be written as the equation of motion for a one-form although the Bianchi identity is not satisfied in general.
We now turn to the action. By using the variational principle, one should be able to find non-trivial solutions of the equation of motion, and subsequently obtain the corresponding non-perturbative action by substituting back the solution in the original action -in parallel with what was done in section 2.1. It is often remarked how the axion gets an extra sign in its kinetic term when turning the action from Minkowski to Euclidean; here we observe that regardless of the sign, starting from the variational principle in Euclidean space yields a vanishing wormhole action, as we shall show below. This is remedied either using a three-form description -as done in section 2.1 -or introducing Lagrange multipliers [17,103]. Since a general coset does not seemingly accept a three-form description we opt for the second option and write an action where κ = 8πG N , we have suppressed scalar indices which otherwise are summed univocally to build invariants and the Lagrange multiplier's (λ A (x)) coefficient is the divergence of the Noether currents (J µ = √ gf 2 ξG∂ µ θ) which vanishes on-shell. Here we have used θ coordinates yet the following derivation does not require explicit expressions for the metric G or the Killing vectors; in terms of the one form ω the kinetic term of the Goldstones reads simply ω µ A ω A µ /2.

The variational principle yields the equations of motion δS
Taking the trace in Einstein equations, if we omit λ terms, yields −R + κf 2 ∂ µ θG∂ µ θ = 0 and a vanishing action (regardless of the sign that one can change with κf 2 ). The rôle of the extra λ term, as we shall see in the following, is exactly to turn this result into a properly defined variational principle for wormhole solutions with non-vanishing action. In this regard, our result is the generalization of the construction put forward in [17] to the case of a generic coset. Most of the computation now involves manipulations of the equations of motion in order to put them in a clean and readable form, and we proceed as follows. In the equation for the Goldstones we perform the substitution δθ i = ∂ ν λ B ξ i B , and, after repeated use of the definition of Killing vectors with respect to the metric in eq. (117), we find where the last term is antisymmetric in the summed indices A, B and proportional to the structure constants. This means that, for spherically symmetric solutions which depend only on the radius r as the ones we are interested in, such a term cancels, since ∂ r λ A ∂ r λ B is symmetric in A, B. From now on, we focus on this type of solutions, θ(r) i , λ A (r), and we adopt the spherically symmetric metric in eq. (12). As remarked above eq. (123) reduces in this case to solved by ξ A ∂ µ λ A = ∂ r θ. On passing, we note that this substitution in Einstein equations yields R + κf 2 ∂ r θG∂ r θ = 0, and a non-zero on-shell action. This relation for λ back into the Goldstone equation of motion eq. (121) yields, substituting δθ = ∂ r θ We now use the eq. (122) for the variation of the Lagrange multiplier to substitute ∂ r ξ. Eq. (125) above becomes and, after introducing a total derivative which can be integrated to find where for consistency the constant C ≥ 0. Finally, we substitute the Goldstone boson kinetic term in eq. (128) back in the Einsteins equations which, remarkably, have the same structure and solution for the metric as found in the axion case, see eqs. (14)(15). The positive-definite constant C can be given in terms of the n Noether currents which are constants of motion √ gξ A G∂ r θ = C A ; this in certain cases takes a simple form as the sum of the squares but a general expression for all cases is not known to us. Let us now move to discuss the most important aspect of this computation, that is the potential that these wormhole solutions generate as a consequence of the explicit symmetry breaking that they induce. As emphasized in the discussion we put forward in the first part of this letter, the crucial point in this respect is to understand the correct quantization condition. In the simple case of the U (1) symmetry, the quantization condition is related to the periodicity of the axion field. In the generalization to the coset space G/H, the fact that the Goldstone manifold is compact discretizes the possible charges and, as a consequence, the wormhole action. However, the discussion of the geometry a generic coset can be forbiddingly complicated. To shed light on the issue, we shall consider in the rest of this section an explicit example. In order to make an educated choice, as ever, symmetry is helpful, in particular the unbroken symmetry. In this respect the maximal unbroken symmetry allowed is O(n) for n Goldstones, which is the largest linearly realised symmetry that the associated tangent space, ω A µ , admits. A breaking pattern that yields this unbroken symmetry is O(n + 1)/O(n), which we will adopt as an explicit case to exemplify our discussion. One has n broken generators X A and n(n − 1)/2 unbroken ones T a that read -we will use antisymmetric generators and therefore there is no need for factors of ı in the following: The Goldstone fields can be taken as θ i with i = 1, .., n but, motivated by the symmetries of the unbroken group, we shall use in the following spherical coordinates θ = ρ u(ϕ), with u · u = 1. We therefore have one radius ρ and n − 1 angles ϕ i . The Goldstone matrix reads where we note that ρ + 2π gives the same group element as ρ. This yields the one-form which, in turn, gives a kinetic term for the Goldstones where (∂u/∂ϕ) 2 is the metric in a S n−1 sphere. Using eq. (118), the Noether currents can be found to be and in this case we can find the explicit connection between kinetic term and Noether currents by taking the sum over the unbroken currents We now proceed in parallel to section 2.1. We denote an instanton solution withρ ,φ, and we consider a variation of the fields δρ, δϕ i around it. The first variation of the action reads, after using the equations of motion for the background: If we take the variations to be independent from the solid angle Ω 3 , we can rewrite δS in terms of charges Q A(a) = dΩ 3,1 J r A(a) . We find where there is a sum over broken Q A and unbroken Q a charges. The explicit form of this expression is not relevant, what is important to note is that if gravity is to respect the unbroken symmetry, which certainly would be the case if it were gauged, the solutions would have Q a = 0 and therefore ϕ coordinates would not get a potential. On the other hand the radius ρ gets a contribution proportional to broken charges Q A . Furthermore, the periodicity of this field, in full analogy with the axion case, discretizes the possible charges to beū A Q A = k, with k ∈ Z. We can expect therefore a potential generated for the radial coordinate as follows which nevertheless produces a mass term for all Goldstone bosons. To visualize this fact, notice that the minimum of eq. (140) sits at ρ = 0 were the spherical coordinate system is singular. It is therefore better to go back to θ coordinates. We find This means gravity produces a mass for all Goldstone bosons. Although obtained in the specific case of the O(n + 1)/O(n) symmetry breaking pattern, it is tempting to generalize this result to a general coset. We postpone the verification of this conjecture to future investigation.

Conclusions and outlook
The explicit breaking of global symmetries induced by gravitational effects is a mesmerising phenomenon. Furthermore, far from academic, this issue could have direct implication in phenomenology.
In this paper we considered the specific case of the global shift symmetry of a Goldstone boson, and we explored the non-perturbative breaking due to Euclidean gravitational instantons in the context of Einstein gravity (aka wormholes).
This field has been studied in a decades-long endeavour, and in this regard the main novelties of our analysis are the following.
• Theoretical analysis. On the theory side, we addressed the problem of stability of wormhole solutions by computing the spectrum of the quadratic action, and we found a positive spectrum. This result is of fundamental importance since only in the absence of negative eigenvalues the gravitational instantons mediate tunneling transitions between degenerate vacua, in parallel with other known situations both in Quantum Mechanics and Quantum Field Theory. In turn, this property allows to consistently compute the effective potential generated by gravitational instantons. The latter, as a consequence of its nonperturbative nature, is characterized by the suppression factor e −Sinst , where S inst is the wormhole action which scales with the Goldstone decay constant and Planck mass as S inst ∝ M Pl /f a .
• Phenomenological analysis. On the phenomenological side, we focused our analysis on the compelling case of the QCD axion. By computing non-perturbative gravitational corrections to the QCD axion potential, we found the following lower bound on the mass of the QCD axion m a 4.8 × 10 −10 eV , f a 10 16 GeV .
As an application, we discussed important consequences related to black hole superradiance, and we showed that the mass range m a = [10 −14 , 10 −10 ] eV motivated by experimental searches and phenomenological considerations is theoretically disfavored.
In addition to the QCD axion, we discussed the case of ultralight scalars as cosmological dark matter. We showed that non-perturbative gravitational effects due to Einstein gravity generate a mass term in agreement with the estimate of the relic abundance based on the misalignment mechanism. In numbers, we found a dark matter candidate with mass m a 2.5 × 10 −18 eV corresponding to the decay constant f a 8 × 10 15 GeV.
Finally, we discussed to rôle of non-perturbative gravitational corrections in connection with the relaxation of the electroweak scale and the clockwork mechanism.
We concluded with a derivation of the of wormhole solution for a generic Goldstone coset, and we inspected the breaking pattern O(n + 1)/O(n) to find that all Goldstone bosons in the coset acquire a mass via gravitational effects.

A Dual Formulation of broken symmetry phases A.1 Axions and forms
Let us start considering the theory of a massless Goldstone field θ described by the Lagrangian density For simplicity, we restrict here to flat space. It is straightforward to check that L θ can be obtained from after eliminating the non-dynamical current J µ by means of its algebraic equation of motion J µ = −if 2 a (∂ µ θ). Equivalently, after integration by parts, we have where now the Goldstone boson field is non-dynamical, implying the constraint ∂ µ J µ = 0. The Lagrangian density in eq. (143) is therefore nothing but L J = J µ J µ /2f 2 a , subject to the current conservation constraint ∂ µ J µ = 0. The latter can be formally solved by introducing the antisymmetric tensor field B, with J µ ≡ − i 2 µνρσ (∂ ν B ρσ ). From these simple manipulations follows that the Goldstone theory admits a description in terms of an antisymmetric tensor field. Using the condition J µ = −if 2 a (∂ µ θ), we have where H is the dual strength of the antisymmetric tensor field B, H µνρ = ∂ µ B νρ + ∂ ν B ρµ + ∂ ρ B µν . Notice that this definition implies the Bianchi identity µνρσ (∂ ρ H σµν ) = 0. Eq. (146), f 2 a µαβλ (∂ µ θ) = H αβλ , incarnates the dual description of the field θ in terms of a three-form. Notice that at this level of the analysis we did not specify the parity transformation of the Goldstone field. On the contrary, the construction in eqs. (143-146) remains valid both for scalar and pseudo-scalar bosons. From the dual relation 2f 2 a (∂ µ θ) = µνρσ (∂ ν B ρσ ), it follows that if θ is a pseudo-scalar (scalar) field, then the two-form B ρσ is forced to transform as a tensor (pseudo-tensor). Without additional input one can not tell one option from the other, both perfectly consistent (for instance, in string-inspired situation the two-form B ρσ can be identified with the Kalb-Ramond field that, after compactification in 4-dimensions, is an even parity tensor; having this specific dual picture in mind, the field θ transforms as a pseudo-scalar).
In the dual picture, the free Lagrangian density in eq. (145) takes the form The antisymmetric tensor B ρσ has 6 independent components, but the gauge symmetry with the extra redundancy Λ σ = Λ σ + ∂ σ λ correctly reduces to one the number of dynamical degrees of freedom.

A.2 From Minkowski to Euclidean space
Let us consider the action in curved space-time, with flat metric η µν = diag(+1, −1, −1, −1) describing a massless Goldstone boson field θ minimally coupled to Einstein gravity. As discussed in section 2.1, the field θ admits a dual description in terms of an antisymmetric two-form with field strength H µνρ given by H µνρ = f 2 a µνρσ (∂ σ θ). In Lorentzian space-time, we have the following identity with F defined right below eq. (5). The minus sign follows from the Levi-Civita contraction µνρσ µνρλ = (−1) t 3! δ λ σ , with t = 1 in Minkowski space-time. The action in eq. (149) takes the form We are now ready to go from Minkowski -with coordinates x ≡ (t, x), and metric g µν (x) -to Euclidean space-time -with coordinatesx ≡ (t E , x), and metricg µν (x) -by means of the Wick rotation t E ≡ ıt. The Euclidean action is related to eq. (149) via S E ≡ −ı S| t=−ıtE . The analytical continuation can be seen as the coordinate transformation x ≡ (t, x) →x ≡ (t E , x) = (λt, x), with λ = ı. We can therefore write the general coordinate transformation from which it follows that g = λ 2g = −g, and d 4 x √ −g = d 4x √g (−ı). The Euclidean version of eq. (151) is which coincides, after renamingg, with eq. (5). In Euclidean space, the duality relation in eq. (150) reads where we now used µνρσ µνρλ = (−1) t 3! δ λ σ , with t = 0 in Euclidean space-time. Consequently, we reconstruct the Euclidean action in eq. (8).
Notice that it is crucial to define the Euclidean action by analytical continuation of eq. (151) rather then eq. (149). Indeed, the analytical continuation of eq. (149) would generate the Euclidean action that features the 'wrong' minus sign in front of the Goldstone kinetic term. Finally, we stress again that the construction we put forward in this appendix, as well as the consequent wormhole solutions, does not rely on the transformation properties under parity of the Goldstone boson. This suggests that wormhole solutions can be defined for both free massless scalar and pseudo-scalar fields.

A.3 Axion charge quantization
In this appendix we provide a more quantitative picture of the quantization condition discussed at the end of section 2.1. To this end, we exploit the analogy with the quantization of the electric charge in the presence of a magnetic monopole.
The magnetic field of a putative magnetic monopole placed at the origin of R 3 is where in the second equation the magnetic charge g M corresponds to the magnetic flux through the sphere S 2 . Bearing in mind the relation with the electromagnetic field strength, F ij = − ijk B k , one can already identify the analogy with the three-form field H ijk = ijk H 0 = ijk n/2π 2 r 3 and the computation of the flux in eq. (18). Of course, in the electromagnetic case we are dealing with a two-form field strength F µν = ∂ µ A ν − ∂ ν A µ while in the axion case the field strength is the three-form H µνρ = ∂ µ B νρ + ∂ ν B ρµ + ∂ ρ B µν but the comparison is evident. The rôle of the electromagnetic potential A µ is played by the antisymmetric tensor B µν . The obstacle is that it is not possible to find a vector potential such that ∇ × A = B since the monopole field is not divergenceless. The best that one can do is to define (in ordinary spherical coordinates) the vector potentials defined, respectively, in R 1 : θ ∈ [0, π/2 + δ) and R 2 : θ ∈ (π/2 − δ, π] with an overlap region π/2 − δ < θ < π/2 + δ. In their domains, both potentials satisfy ∇ × A 1,2 = B but they are singular in the complementary region. This is the same ambiguity one encounters for the two-form potential B in the axion case. The only way to obtain a consistent picture is to show that the two vector potentials describe the same physics in the overlap region. Said differently, A 1 and A 2 must be related by a gauge transformation A → A = A + ∇χ. From A 1 − A 2 = ∇χ one immediately finds χ ≡ g M φ/2π. Notice that the function χ(φ) is not continuous since the azimuthal angle is defined modulo 2π. This is actually crucial for the existence of a non-zero flux in eq. (156) since Dirac quantization proceeds as follows. We consider a test particle of mass m and charge q in the field of the magnetic monopole. Its wavefunction satisfies the Schrödinger equation −/ h 2 /2m( ∇+ıe A) 2 ψ = ı/ h∂ψ/∂t, with e = q// h, and the gauge transformation A → A = A+ ∇χ transforms ψ into ψ = e −ıeχ ψ = e −ıegMφ/2π ψ. Since φ = 0 and φ = 2π are the same physical point, consistency of Quantum Mechanics requires the wavefunction ψ to be single valued. The quantization condition eg M = 2πk, with k ∈ Z, then follows. Dirac quantization involves the product eg M of electric and magnetic charge. The electric charge of a particle is given by the surface integral at spatial infinity over a S 2 sphere of the Hodge dual of the electromagnetic field strength two-form F µν . Conversely, the magnetic charge of a particle is given by the surface integral at spatial infinity over a sphere S 2 of the electromagnetic field strength two-form F µν , as already discussed in eq. (156).
The quantization in the axion case is based on the fact that potentials of higher rank admit the following generalization. Formally, in D dimensions, the 'electric' charge Q E is the integral over the sphere S D−(p+2) of the Hodge dual of the (p + 2)-form field strength. For the axion three-form H, p = 1, and in D = 4 dimensions Q E = S1 * H. The 'magnetic' charge, as already noticed, is Q M = S3 H = n (see eq. (18)), that in our case corresponds to the axion charge flowing through the wormhole throat. Dirac quantization therefore generalizes to Q E n = 2πk. If Q E = 0, the quantization of n trivially follows. In terms of the Goldstone field, we have Q E = S1 * H = S1 dθ. The key observation is that this integral encompasses a closed path. This means that Q E can be different from zero only if θ is multi-valued, and this is exactly the case for the Goldstone field since it is a phase, and possesses the periodicity θ = θ + 2π. This concludes the argument, and shows the fundamental importance of the periodicity of the Goldstone boson field for the quantization of the axion charge along the wormhole throat.

B The Pöschl-Teller potential
We want to solve the spectral problem for the differential operator  Eq. (160) is describes the one-dimensional motion of a particle in the Pöschl-Teller potential [106]. In full generality, the dynamics in the Pöschl-Teller potential is described by the differential equation −u (x) − α 2 ξ(ξ−1) cosh 2 (αx) u(x) = 2Eu(x). Eq. (160) corresponds to α = 1, ξ = 5/2. By changing variable y ≡ cosh 2 x we find y(1 − y)u (y) + 1 2 − y u (y) − 15 16y After the rescaling u(y) = y 5/4 v(y) the problem reduces to the hypergeometric differential equation with ξ = 5/2. We have the general solution where 2 F 1 (a, b; c; z) is the Gaussian or ordinary hypergeometric function, and we introduced the short-hand notation We can separate the two independent solutions. If B = 0, A = 1, we have the even solution u even (x) = cosh ξ x 2 F 1 a, b; If A = 0, B = ı we have the odd solution u odd (x) = cosh ξ x sinh x 2 F 1 a + 1 2 , b + 1 2 ; 3 2 ; − sinh 2 x .
In order to identify the bound states with negative energy solutions, we need to study the asymptotic behavior. Using the trigonometric definitions sinh x = (e x − e −x )/2, cosh x = (e x + e −x )/2 we find In the presence of negative energy states, we have The coefficient of e −ı √ 2E must therefore vanish in order to ensure the correct asymptotic behavior. This condition leads to quantization of energy. The Gamma function diverges for negative integer numbers, and we can impose this condition to cancel the coefficient of e −ı √ 2E in both u even (x) and u odd (x). For the even functions we find whereas for the even functions we find with m = 0, 1, 2, . . . . With ξ = 5/2 we only have one negative (even) state with energy 2E (even) 0 = −9/4 or, equivalently, λ = −8. Next, we have two states with λ = 0 corresponding to 2E Even eigenfunctions : u λ=−8 even (τ ) = 1 cosh 3/2 (2τ ) , u λ=0 even (τ ) = 3 − cosh(4τ ) 2 cosh 3/2 (2τ ) .

C Spectrum of tensor and inhomogeneous scalar perturbations
In this appendix we compute the spectrum of tensor and inhomogeneous scalar perturbations of the quadratic action expanded around the wormhole background.

C.1 Tensor perturbations
We exploit the results of [33] in which tensor fluctuations around instanton background were studied. The perturbed line element, in Euclidean space, is ds 2 = a 2 (τ ) dτ 2 + (γ ij + t ij )dx i dx j , where γ ij is the background metric on the three-sphere and a(τ ) the background conformal factor. As discussed in section 3.1, t ij is a transverse trace-free symmetric tensor. The Euclidean quadratic action takes the form [33] δS (2) E tensor where covariant derivatives and contractions are defined w.r.t. the metric γ. By introducing the rescaling t ij ≡ at ij , and integrating by parts, one finds [33] δS (2) E tensor with the Schrödinger-type operator where in the last step we substituted the wormhole background solution. Notice that the quadratic action in eq. (175) does not present the conformal factor problem. This is indeed correct since, as reviewed in section 3.1, the quadratic action for the trace-free part of the tensor fluctuations is not unbounded from below. The differential operatorK acts on the τ variable while the Laplacian operator ∆ acts on the threedimensional space. We can therefore diagonalize them independently and add their spectra. The differential operatorK does not possess negative eigenvalues since the potential is a positive function and there are no bound states. The tensor harmonics (G ij ) n lm (ψ, φ, ϕ) are tensor eigenfucntions of the Laplacian operator on lm (G ij ) n lm (ψ, φ, ϕ), with C (n) lm arbitrary constants. From this general discussion it follows that the tensor Laplacian −∆ in eq. (175) has positive eigenvalues, starting from λ = 6. We conclude that tensor fluctuations around the wormhole background have a positive spectrum.

C.2 Inhomogeneous scalar perturbations
We now move to discuss inhomogeneous scalar perturbations. In section 3.1 we focused on scalar homogeneous perturbations. A separated analysis is required since, as well known [34,35], inhomogeneous perturbations have positive Euclidean action and do not suffer from the conformal factor problem. We exploit the results of [35]. In the Euclidean space, the quadratic action takes in form with the Schrödinger-type operator As before, the differential operatorÔ acts on the τ variable while the Laplacian operator ∆ acts on the threedimensional space. We can therefore diagonalize them independently and add their spectra. Furthermore, we already know that −∆ has positive eigenvalues starting from λ = 3 (ignoring the homogeneous mode). In eq. (177) q is the only physical dynamical variable describing scalar perturbations after gauge degrees of freedom are removed [35]. By defining 2τ ≡ x, the differential equationÔu(τ ) = λu(τ ) in eq. (178) takes the form This is again a Pöschl-Teller potential with α = 1, ξ = 3/2. We can repeat the analysis presented in appendix B. With ξ = 3/2, we have a double-degenerate negative eigenvalue E even/odd 0 = −1/4 but it corresponds to the positive value λ = 3, that is the starting point of the spectrum of the differential operatorÔ. We therefore conclude that the spectrum of the bulk operator in eq. (177) is fully positive.