Cosmological Aspects of the Clockwork Axion

The clockwork axion refers to a family of aligned multi-axion models that lead to an exponential hierarchy between the scale of Peccei-Quinn symmetry breaking and the scale of the axion decay constant. The clockworking can bring the Peccei-Quinn-scale particles to within reach of collider experiments.In this work we are interested in whether cosmological observations impose any new constraints on the clockwork axion. If the universe reheats above the scale of Peccei-Quinn breaking, then the ensuing cosmological phase transition produces a network of topological defects, which have a qualitatively different behavior from the string-wall network in the usual axion models. We estimate the relic abundances of axion dark matter and dark radiation that arise from the emission of axions by the defect network, and we infer a constraint on the scale of Peccei-Quinn breaking and the mass spectrum. We find that the defect contribution to the axion dark matter relic abundance is generally negligible. However, the defect production of relativistic axion dark radiation becomes significant if the scale of Peccei-Quinn symmetry breaking is larger than $100 \, {\rm TeV}$, and measurements of $\Delta N_{\rm eff}$ provide a new probe of this class of models.


Introduction
The axion [1][2][3][4][5][6][7][8] is an elegant solution to the strong CP problem that also naturally provides a candidate for the dark matter. (See Ref. [9] for a general review.) The axion solution introduces new particles and interactions at a scale, f pq , set by the spontaneous breaking of the global U(1) Peccei-Quinn (PQ) symmetry, and the axion is the corresponding pseudo-Goldstone boson. Generally f pq is well above the weak scale, and upon integrating out the PQ-scale particles, the axion acquires dimension-five interactions with the gluons (and possibly other Standard Model particles as well) with a coupling g agg ≡ α s /(2πf g ), which defines the parameter f g . Due to instanton effects in quantum chromodynamics (QCD), this interaction induces a potential for the axion, which lifts its mass to [10,11] m a 5.6 µeV f g 10 12 GeV (1) In standard axion models one usually finds f g = f pq /N DW , where the QCD domain wall number, N DW , is typically an O(1) integer that grow linearly with the number of PQ-scale colored fermions. Consequently, the most well-studied axion models predict the scale of the axion-gluon coupling and the scale of PQ breaking to be comparable, f g ∼ f pq .
In this work, we study the cosmological implications of the following relation: f pq = zf g with z 1 and N DW = 1 .
In general aligned axion models with multiple pseudoscalar fields [12] allow one to break the canonical relation between f pq and f g . The clockwork mechanism [13][14][15] provides a concrete benchmark model in which Eq. (2) can be realized with an exponential hierarchy. The mechanism, which will be reviewed in Sec. 2, introduces N + 1 complex scalar fields with global U (1) symmetries that are spontaneously broken at a scale f pq generating N + 1 Goldstone bosons. Additionally, a specific scalar potential explicitly breaks N of the symmetries lifting all of the flat directions save one, which corresponds to an axion with axion-gluon coupling strength set by f g ∼ q N f pq for integer q. The model with q = 3 is renormalizable and well-defined [15], and we will use it for all the studies in this paper.
The regime f pq f g is particularly interesting for both particle phenomenology and cosmology. Whereas f g must remain larger than roughly 10 8 GeV to satisfy astrophysical constraints on axion-gluon interactions (for a review see Ref. [16]), the scale of PQ breaking could be as low as f pq ∼ 1 TeV. This puts the new PQ-scale particles within reach of collider experiments [17]. If the primordial plasma reaches temperatures above f pq , which is quite reasonable for compelling models of inflation and reheating, then PQ-symmetry breaking is accomplished through a cosmological phase transition. It is well known from standard QCD axion models that the PQ-breaking phase transition produces a network of topological defects [18], namely global strings and domain walls, which efficiently radiate the Goldstone axion [19][20][21][22][23] and contribute to the axion relic abundance. Previous studies of these defect networks have generally assumed f pq ∼ f g , but this relation is badly broken for the clockwork axion, and moreover, the presence of N + 1 complex scalar fields participating in the phase transition gives the defect network a much richer structure than what is found in the usual QCD axion cosmology.
In this article we are interested in the structure of the defect network and its contribution to the axion dark matter and dark radiation relic abundances. For the usual QCD axion, the collapse of the defect network at the QCD epoch (t = t QCD ) is one of the dominant contributions to the axion dark matter relic abundance [16], which corresponds to an energy density ρ a (t QCD ) ∼ ε σH QCD where ε measures how efficiently axions are emitted from the domain walls, σ ∼ m a f 2 g is the surface tension of the axion-field domain walls, and H QCD = H(t QCD ) is the Hubble parameter. Since the domain walls are composed of the axion field, the emission of axion particles is very efficient and typically ε ∼ 1 [19][20][21][22][23]. In the clockwork axion model we will see that the domain walls are built from the additional pseudoscalar (gear) fields with mass m G , and the wall tension is σ ∼ m G f 2 pq . Therefore it is reasonable to guess ρ a (t QCD ) ∼ ε m G f 2 pq H QCD , which corresponds to a relic abundance of axion dark matter given by This estimate suggests that the observed dark matter relic abundance (Ω DM h 2 0.12) can be explained if the clockworking lowers the PQ scale to within reach of collider experiments. Moreover, Eq. (3) suggests a strong upper bound on m G and f pq in order to avoid producing too much axion dark matter. However, it is not immediately obvious how large is the factor ε , which represents the efficiency with which axions are emitted from the domain walls, which are now composed of the gear fields rather than the axion field itself.
To confirm or refute the naive estimate, we have performed a more careful calculation of the axion emission from the clockwork defect network. Since the pseudoscalar (gear) fields that form the domain walls only have a weak coupling to the axion, it is found that axion emission in the clockwork model is generally much less efficient than for the usual QCD axion where the axion field itself composes the walls; that is, Eq. (40) reveals that ε ε ∼ 1 for Hubble-scale domain walls at the QCD epoch. The result is a suppression of the axion dark matter relic abundance arising from the defect network. (The relic abundance from the misalignment mechanism depends on the axion decay constant f g directly, and it is unaffected by the clockworking.) On the other hand, we also find that axion emission from small-scale structure on the domain walls is still efficient, and this can give rise to a population of relativistic axions that contribute to the dark radiation.
The remainder of this article is organized in the following way. We briefly review the clockwork axion in Sec. 2 with an emphasis on understanding the topology of the field space. In Sec. 3 we focus on the cosmology of the topological defect network and its contribution to the axion dark matter and dark radiation relic abundances. We summarize our results in Sec. 4 and briefly discuss some more general cosmological implications.

The Clockwork Axion
The clockwork axion is a class of multi-axion models [12] with an exponentially large hierarchy between the PQ-breaking scale, f pq , and the axion decay constant, f a . Different implementations were explored in the original literature [13][14][15] and subsequent work [24][25][26][27][28][29]. In this section we briefly review the model that was proposed in Ref. [15].
The model consists in a family of N + 1 complex scalar fields that are denoted by φ n (x) with n ∈ {0, 1, · · · , N }. The properties and interactions of these fields are determined by the The four model parameters are the dimensionless couplings λ and , the energy scale f pq , and the integer N . The first term in the Lagrangian (4) is invariant under the symmetry with θ n ∈ [0, 2π) for n ∈ {0, 1, · · · , N }. The second term in Eq. (4) explicitly breaks this symmetry down to the subgroup with θ ∈ [0, 2π) and q ≡ 3. We see that field φ n has charge q N −n under the U(1) pq symmetry. Notice also that the U(1) pq symmetry has a family of discrete subgroups, which correspond to taking θ = 2πq k−N with k ∈ {0, 1, · · · , N − 1} in Eq. (6). The scalar potential in Eq. (4) leads to spontaneous symmetry breaking. The N + 1 complex scalar fields acquire nonzero vacuum expectation values φ n = (f pq / √ 2) 1 + O( ) where the -dependent shift is not independent of n. Consequently the U(1) N +1 symmetry is spontaneously broken without leaving any unbroken subgroup. Since the -dependent term in Eq. (4) explicitly breaks U(1) N +1 down to U(1) pq , the spectrum contains N massive pseudo-Goldstone bosons and 1 massless Goldstone boson.
In order to understand the vacuum structure of this theory, it is useful to parametrize the complex scalar fields φ n in terms of real scalar fields ρ n and the pseudoscalar fields π n by writing 1 Eq. (4) has an approximate (flavor) permutation symmetry acting on the φ n . This symmetry is a consequence of not allowing for generic λ n , n , f n , and in general need not to be present. However, it simplifies the model and the description of the defects dynamics in the early universe. The symmetry is broken by the -term and by the coupling to PQ matter fields. When the dynamics is perturbative in and y, we will exploit the flavor symmetry to understand the scaling of some observables with the parameter of the models. The discussion of the effects of the breaking of this symmetry can be found in [27].    p X H h a a R L 9 u 9 E z Z T 3 c 1 U E p 2 J 4 5 u 9 r C / J f 2 q j C c m 9 c C 2 0 r B M 1 v F 5 W V p G j o o g 0 6 F Q 4 4 y n k A j D s R b q X 8 j D n G M X R 2 Z 0 u c + U D F T E 8 d / J B G z + I j Z 2 x h L u P P w v N Y z V u j b / 4 v F h a p k O m G x t P 7 / T 4 E 2 X D w c Z B + 2 e k f D N v q V 8 l b 8 p 5 8 I C n Z J Q f k E z k m G e H k i v w k 1 + R X 5 y b a i t 5 E 7 2 6 t U a f N v C J 3 J t r + A 2 a s v M w = < / l a t e x i t > Figure 1: For the model with N + 1 = 2 complex scalar fields the pseudoscalar field space (π 0 , π 1 ) is topologically equivalent to a 2-torus. For = 0 and QCD instanton effects neglected, there is no potential for (π 0 , π 1 ); taking = 0 induces a potential, leaving a single flat direction, which is represented by the black line; and at the QCD epoch, instanton effects lift also this flat direction. Along the flat direction (black line) the axion field covers a ∈ [0, 2πf a ) (13) while the field π n completes q N −n cycles (10) corresponding to its PQ charge (6). The three blue lines indicate three π 0 -strings, each connected to one domain wall; the red line indicates a single π 1 -string connected to three domain walls; and the black line indicates an a-string. See discussion in Sec. 3.2.
The full configuration space is covered when ρ n (x) ∈ [−f pq , ∞) and π n (x) ∈ [0, 2πf pq ). Note that the Lagrangian is left invariant under the discrete gauge symmetries that take π n → π n + 2πf pq , which correspond to the trivial transformations φ n → e 2πi φ n = φ n . After modding out by this discrete gauge symmetry, the configuration space for the π n becomes compact. For instance, if N = 0 then the configuration space of π 0 is topologically equivalent to a circle; if N = 1 then the configuration space of (π 0 , π 1 ) is a torus; and so on. The N = 1 case is illustrated in Fig. 1.
In the vacuum where φ n = f pq / √ 2 we identify the two mass parameters, which set the scales of the real scalar and pseudoscalar mass spectra. The N + 1 real scalar fields ρ n have masses on the order of m ρ with a splitting on the order of m G . We consider more carefully the spectrum of the N + 1 pesudoscalar fields π n , since the massless mode will be identified with the QCD axion. Let us denote the mass eigenstate pseudoscalar fields by A i ; we identify the massless axion as a ≡ A 0 and the N massive gears as G i ≡ A i for i ∈ {1, 2, · · · , N }. A linear transformation relates the original field basis π n and the mass-diagonalized field basis where one should remember that π n are evaluated modulo the discrete gauge symmetry π n → π n + 2πf pq . The elements of the orthogonal matrix O are given by [24] O n0 = C q n and where q = 3 and C ≈ 8/9 0.94. The pseudoscalar mass spectrum is given by where G i for i ∈ {1, 2, · · · , N } are the N pseudoscalar gear fields, and a is the massless axion. The specific wave-function for the massless mode ensures a very large field range for the axion. Indeed, since we must allow π n (x) ∈ [0, 2πf pq ) in order to span the full configuration space, the property O n0 1 for n = N 1 requires the axion field range be exponentially large. This implies that the axion field has a range where In other words, the transformation a → a + 2πf a is a discrete gauge symmetry, e.g. a field configuration with a ≥ 2πf a is redundant with a configuration that has 0 ≤ a < 2πf a . (Similar considerations give the allowed field range for the gears, G i .) For the case N = 1 this calculation can be visualized concretely as in Fig. 1. Due to the clockworking (13), we can have f a ∼ 10 9 − 10 11 GeV in the classic axion mass window with f pq ∼ 10 TeV and N ∼ 10 − 15, which potentially puts the new PQ-scale particles within reach of collider experiments. In order to interpret the Goldstone boson a as the QCD axion, we must let the φ n couple to colored fermions. Let ψ r and ψ r c represent a pair of vector-like colored fermions in representations r and r c of SU(3) c . We introduce a Yukawa coupling between these fermions and the scalar field at site-l, which is written as yφ l ψ r ψ r c . As observed in [17], a phenomenologically "safe" choice corresponds to localizing the QCD fermions on the last site of the chain, which corresponds to l = N . Any other choice would lead to an exponentially large domain wall number N DW (see below). Therefore, for l = N , integrating out the heavy colored fermions, ψ r and ψ r c , induces the color anomaly, which is where T (ψ) is the index of the color representation (T (3) = 1/2), and d ψ the dimension of the SU(2) L representation (d ψ = 1 for singlet). Writing π N in terms of the axion a using Eqs. (10), (11), and (13), and setting the gear fields to zero, the color anomaly becomes The mismatch between the period of the discrete gauge symmetry (2πf a ) and the period of the effective QCD theta parameter (2πf g ) is the domain wall number, N DW . A reasonable choice is to consider ψ r in the same representation as the u R or d R quarks, which gives N DW = 1 and avoids color stable relics with mass m ψ = yf pq , which are a potential disaster if PQ breaks after inflation. With this setup we can compute the axion mass; it is given by Eq. (1) where f g = f a .

Comparison with other multi-axion models
In this work, we focus on the model described by Eq. (4), since it is perhaps the simplest implementation of the clockwork axion. However, there exist compelling variations on the minimal model that do not require the introduction of N ≈ 10 fundamental scalar fields. The low energy particle phenomenology is not appreciably changed from the minimal model to these variations, but the cosmology can differ. We briefly discuss generalizations and their implications in this section. Several alternative UV extensions of the clockwork axion were recently suggested by the authors of Ref. [30]. The common feature of these variations is that the explicit symmetrybreaking term is allowed to arise dynamically. In effect, the parameter in Eq. (4) becomes time dependent. A first implementation supposes that the symmetry-breaking ( ) term arises from strong dynamics. In this scenario, is effectively zero at high energy scales, but it takes on a fixed, nonzero value below a certain energy scale corresponding to the confinement of the new strong force. A second implementation supposes that the symmetry-breaking term arises from a Froggat-Nielsen-type mechanism. The term is replaced by a non-renormalizable operator, such as n φ * n (Σ n,n+1 /Λ)φ 3 n+1 where Λ is the cutoff of the effective theory, and Σ n,n+1 transforms as a bi-fundamental. After Σ gets a vacuum expectation value, the symmetry-breaking term in Eq. (4) is generated with ∼ Σ /Λ. This hypothesis naturally accommodates small values, 1, since it is related to the ratio of hierarchical energy scales. Provided that the new particles and interactions are not within reach of collider experiments, then the predictions for particle physics observables are not significantly changed from the predictions of the effective theory in Eq. (4). However, the cosmology can be somewhat different, because is effectively time dependent. If PQ-symmetry breaking occurs while = 0 then the network of topological defects is consists of cosmic strings, without any domain walls. This is simpler than the scenario we describe next in Sec. 3 where the strings are bound to walls during the prior from PQ breaking until the QCD epoch . Provided that becomes nonzero before the QCD epoch, and assuming that m G > H at this time, the strings become bound by domain walls. Then the predictions for cosmological observables, like the axion relic abundance, are not significantly changed from the scenario we describe in Sec. 3. Thus for the remainder of this article, we focus on the specific implementation of the clockwork axion that is defined by Eq. (4).

Cosmological effects of the clockwork axion
Clockworking allows the PQ-breaking scale to be much smaller than the axion decay constant (13), f pq f a . Since most studies of axion cosmology assume f pq ∼ f a , it is interesting to investigate whether new phenomena can arise when this relation is badly broken. (We will see shortly that the topological defect network is qualitatively different from the standard axion cosmology.) However, more to the point, the PQ-breaking scale is a threshold for new physics: the masses of the scalar fields and PQ fermions are all set by f pq times a dimensionless coupling. Whether or not these states are accessible to laboratory experiments depends sensitively the value of f pq . In this section we explore how cosmological evolution constrains the PQ-breaking scale.
The scenario we are going to discuss is illustrated in Fig. 2. If f pq is smaller than both the reheating temperature and the Hubble scale during inflation, and if the fields of Eq. (4) are in thermal equilibrium with the SM plasma, then a cosmological phase transition occurs when the plasma temperature decreases below T ∼ f pq . During the phase transition the symmetries of the clockwork model become broken, and a network of topological defects forms [31]. In this section we discuss the rich structure and dynamics of the defect network, and we calculate the axion relic abundance that is emitted from the defect network.

Formation of the defect network
Here we are interested in the formation of the defect network at t = t pq , and therefore we can neglect the potential induced by QCD instanton effects, which is negligible at temperatures T ∼ f pq Λ QCD . It is illustrative to first consider the model with = 0 in Eq. (4): in this case the theory admits N + 1 "flavors" of global cosmic strings [18]. We will refer to these strings collectively as π-strings and individually as π n -strings, since we have used π n /f pq to denote the phase of φ n in Eq. (8). For instance, along a trajectory in spacetime that encircles a π n -string with winding number w ∈ Z the phase of φ n evolves smoothly from 0 to 2πw. For = 0 the additional interactions in Eq. (4) break the U(1) N +1 global symmetry to its U(1) pq subgroup, which is then spontaneously broken through the cosmological phase transition. Naively this suggests that we should only consider one string solution, but in the regime where the term can be treated as a perturbation, it is still necessary to consider the N + 1 π-strings, since their formation during the cosmological phase transition is unaffected by a perturbatively small explicit symmetry breaking term. What differs now is that the π-strings become bound to one another by domain walls [18]. Along a trajectory through spacetime that crosses a domain wall, one of the gear fields G i has a nontrivial step-like profile, and therefore we denote this domain wall as a G i -wall.
In addition to the N + 1 flavors of π-strings, the model also admits a global cosmic string, which we call the a-string, associated with the spontaneous breaking of U(1) pq . Since U(1) pq is not broken explicitly the a-string is topologically stable, 2 unlike the π-strings, and also the a-string is not connected to any domain walls. As one traces a trajectory around the a-string in real space, the axion field (13) evolves from a = 0 to a = 2πf a in its configuration space. Simultaneously the field π n for n ∈ {0, 1, · · · N } passes through q N −n complete cycles of π n ∈ [0, 2πf pq ), which can be seen from Eqs. (10) and (11). In other words, the a-string carries the same topological charge as a particular composition of π-strings that combines an n-string with winding number q N −n for each of the n ∈ {0, 1, · · · N }. As noted in Ref. [32], this observation lets us view the a-string as a "bundle" of π-strings that are localized in a region of space that is comparable to the string thickness and much smaller than the inter-string separation. See Fig. 1 to visualize concretely how the string solutions interpolate through the vacuum manifold.
At the QCD epoch the axion mass is lifted by instanton effects, and the residual U(1) pq symmetry is broken. We assume that the breaking of U(1) pq leaves no unbroken subgroup (N DW = 1), and the potential has a single global minimum. Consequently the domain wall no longer interpolates between degenerate vacua, but rather the vacuum with lower energy pushes the wall away. In this way, the entire string-wall network collapses at the QCD epoch, and any energy it carried is converted into particle emission.

Network of cosmic strings
In this section we neglect the presence of the domain walls, and we treat the defect network as N + 1 independent and non-interacting cosmic string networks, which are distinguished by the strength of their coupling to the axion. We are especially interested in the emission of axions from the string network and their contribution to the axion relic abundance [21][22][23].
The reader may question whether it is reasonable to neglect the presence of the domain walls, which can pull on the cosmic strings and thereby affect their motions. This approach is certainly justified if the gear mass ( term in Eq. (4)) arises dynamically at a scale that is much lower than the scale of PQ breaking; for a model that concretely implements this idea, see Ref. [30]. However, if the gear mass is already present at the PQ-breaking phase transition, the situation becomes less clear. We return to discuss the domain walls in Sec. 3.3. Nevertheless, in the end we will see that even with this favorable assumption, the axion relic abundance arising from cosmic strings is subdominant to other contributions.

Cosmic string solution
A global string [18,19] arises in the theory of a complex scalar field φ(x) with a symmetry breaking potential V = λ(|φ| 2 − f 2 /2). The string solution with winding number w ∈ {±1, ±2, ±3, · · · } can be parametrized in cylindrical coordinates as 2λf is the mass of φ, r is the distance from the center of the string, θ is the azimuthal angle, F (x) is the radial profile function, and δ s ∼ 1/m is the thickness of the string core. The profile function satisfies F (x) = Ax |w| for x 1 and F (x) → 1 for x 1 independent of w. The string tension, µ, measures the string's rest energy per unit length, and it is given by The second term in the integrand leads to the well-known logarithmic divergence in the tension of global strings. It is customary to cut off the radial integral at the typical inter-string separation length scale, r ∼ R. Additionally, it is illustrative to split the integral into two pieces at x = 1 corresponding approximately to the edge of the string core at r = 1/m ∼ δ s . Specifying to the clockwork axion model, we replace f → f pq and m → m ρ . Thus the tension of a π-string with winding number w is estimated as where µ core ≈ πf 2 pq arises from x ∈ [0, 1), and the second, larger term arises from x ∈ [1, ∞). The logarithmic factor can be quite large since the typical inter-string separation is given by the Hubble length scale; for m ρ = 1 TeV and R = 1/H QCD we have log(m ρ R) ∼ 50.
As we discussed in Sec. 3.1, the a-string can be interpreted as a composition of the various π nstrings each with winding number q N −n ; see also Ref. [32]. Therefore, Eq. (18) lets us estimate the tension of the a-string (with winding number ±1) to be where we have used Eq. (13) in the last equality. Here the approximation requires that the integral of Eq. (17) is done in a common region of space for all the string of the "bundle." Notice that this estimate agrees with the standard expectation for an axion string, µ a ∼ f 2 a [18,19]. We can also understand this result by noting that the axion field has to pass through a large field excursion, a ∈ [0, 2πf a ), as we circumnavigate the string, and this large field gradient translates into a large energy density and a large tension.

Formation of the a-string is prevented
Although the formation of defects depends crucially on initial conditions and dynamics, we argue that -in the parameter space of interest -the PQ-breaking phase transition produces a network of π-strings connected by G-walls, and that the a-string does not form, neither at the PQ phase transition nor later during the evolution of the network. First we argue that the initial conditions at the time of the PQ-breaking phase transition inhibit the formation of the a-string. During the phase transition the scalar fields φ n develop a tachyonic mass of order m ρ , and they acquire vacuum expectation values according to Eq. (4). In the regime m G m ρ where U(1) N +1 becomes a good symmetry, the fields φ n evolve almost independently, each one "seeing" an identical tachyonic mass of order m ρ , and the defect network that forms consists of π-strings connected by G-walls. In the other regime m G ∼ m ρ where U(1) N +1 is badly broken, the symmetry-breaking potential induces an "anisotropic" contribution to the tachyonic mass of O(m G ), which biases the path of the φ n toward the vacuum manifold, i.e. the flat direction associated with the axion field, and the resulting defect network consists of a-strings without any walls. We are interested in the parameter regime where λ so that m G m ρ , and we expect that the PQ-breaking phase transition produces a string-wall network rather than a network of a-strings.
Next we argue that the a-strings are not formed during the subsequent evolution of the string-wall network. Recall that the a-string has the same topological charge as a composition of q N −n strings of type π n , and for the parameters of interest, q = 3 and N ∼ 15 so that q N ∼ 10 7 . In order to construct an a-string from the string-wall network, it is necessary to combine an enormous number of π-strings. Typically there are only O(1 − 10) long strings in a Hubble volume at any time. Therefore, for N 1 it is very improbable that the a-string will form during the evolution of the string-wall network. In fact a numerical simulation of the evolution of the string-wall network appears in Ref. [32]. For the model with two scalar fields (N + 1 = 2) they observe that the a-string is able to form from the collapse of the string-wall network, but for the model with three scalar fields (N + 1 = 3) the a-string does not form, and it is reasonable to expect that formation of the a-string is also prevented for models with more than three scalar fields.
Therefore, in the remainder of this article, we assume that the PQ-breaking phase transition creates a network of π-strings connected by G-walls. If a network of a-strings were formed instead, then the cosmology would be unchanged from the usual QCD model, and specifically the PQ-scale gear particles would play no role. (We assume the gears decay quickly to gluons through the interaction in Eq. (14).)

Evolution of the string network
If the a-string does not form, as we have argued above, then we can focus on the π-strings. In this section we discuss briefly the dynamical evolution of the string network, and in the next section we estimate the efficiency of axion emission.
A network of π-strings forms at the PQ-breaking phase transition. Initially particles in the plasma scatter frequently on the strings leading to an effective friction that damps their motion.
Since the strings interact most strongly with radial modes and gears, and since these particles go out of equilibrium soon after the phase transition, the friction force quickly becomes negligible, and the string network begins to evolve freely under the pull of its tension. At this time the network contains at least a few long strings that cross the Hubble volume. When long strings intersect one another (or self-intersect), they create large Hubble-scale loops, and subsequently smaller loops are created from the intersections of larger ones. After a few Hubble times of free evolution, these dynamics bring the string network into the scaling regime where new loops are continuously formed at a roughly fixed fraction of the Hubble scale. A loop of size L oscillates with a period of approximately T = L/2 under the pull of its tension. These oscillations cause string segments to be accelerated, which leads to particle emission and gravitational wave radiation. The associated energy loss causes the loop to lose energy and shrink. Eventually the loop size becomes comparable to its thickness, δ s ∼ m −1 ρ , and it finally decays. Then the energy density of the π n -string network at time t is given by [18] ρ (n) where the first term arises from the Hubble-scale long strings, and the second term arises from the string loops. Here N (n) long is the number of long strings, µ is the string tension (18), H = H(t) is the Hubble parameter at time t, and P (n) is the average power emitted by the loop during each oscillation period. The derivation of Eq. (20) assumes that P (n) is independent of t and the loop length L. In general N (n) long and P (n) can be different for the various flavors of π n -strings, whereas µ is universal. We expect that N (n) long is either zero or O(1), because the "chopping" of long strings prevents N (n) long 1. There is a model-independent contribution to the power from gravitational wave emission, which is parametrically P (n) ∼ G N µ 2 with G N Newton's constant, but the emission of Goldstone bosons typically dominates as we discuss in the next section. One can understand the enhancement factor, (P (n) /µ) −1/2 ≥ 1, because loops that do not emit efficiently will live longer and contribute to ρ str .

Axion emission from strings
Let us now consider the emission of axions as the string network evolves from the PQ-breaking phase transition to the QCD epoch. It is useful to first recall the axion emission calculation for standard QCD axion strings.
For the usual QCD axion, one estimates the energy density in axions at time t as ρ a (t) ∼ ε × µ a H 2 × log t/t pq , for example by requiring self-consistency with the scaling solution [18,23,33,34]. The first factor, ε ∼ P/µ ∼ O(1), accounts for the high efficiency with which axion strings emit the pseudo-Goldstone axions. The second factor, µ a H 2 , is the total energy density in the network of strings with tension µ a that are in the scaling regime (20). The logarithmic factor accounts for integrating the axion emission from time t pq to time t. In particular, note that ρ a ∝ µ a ∼ f 2 a .
By translating these results to the QCD clockwork axion model, we can immediately infer the axion abundance. For models with very small explicit breaking, 1 such that L m −1 G , each π n -string will very efficiently emit their corresponding pseudo-Goldstone boson, π n , and we expect an energy density ρ πn ∝ µH 2 where µ ∼ f 2 pq is the tension of the π n -string. Since the axion a is a linear combination of the π n 's, we also expect that the axion energy density is Since µ ∼ f 2 pq ∼ q −2N f 2 a from the clockworking, the predicted axion abundance is exponentially suppressed with respect to the usual QCD axion case. (We will see in Sec. 3.5 that the corresponding axion dark matter is also negligible compared to the misalignment contribution.) The preceding analysis assumes that the gear mass is negligible compared to the loop size, which implies that string loops efficiently radiate all the π n . The more realistic scenario corresponds to the opposite regime, m −1 G L, and here the calculation of ρ a is more subtle, because the emission of massive gears, G i , is exponentially suppressed, and the emission of the massless axion, a, depends on the strength of its coupling to each different flavor of π n -string. Since the loops do not radiate so efficiently, they also live longer, and they carry more energy at late times [P (n) /µ 1 in Eq. (20)]. These effects parametrically enhance the axion relic abundance, but nevertheless the enhancement does not compensate the exponential suppression, f 2 pq /f 2 a ∼ q −2N 1, and the axion emission from the string network is still much smaller than the standard QCD axion result.

Network of gear-field domain walls
In this section we discuss the network of domain walls that connects the π-strings, which are built from the gear fields, and we calculate the relic abundance of axions emitted by the walls.

Domain wall solution
In general domain wall solutions arise in theories for which the scalar potential has two or more minima separated by barriers [18]. If the potential can be approximated by a quartic polynomial V = (λ/4)(ϕ 2 − v 2 ) 2 , then the domain wall solution is parametrized as ϕ = ±v tanh(mz/2) where z is the spatial coordinate normal to the wall (assumed planar), ± is for kink and antikink solutions, and δ w ∼ 1/m is the wall thickness, which is related to the mass m = √ 2λ v. The wall's surface tension σ measures its energy per unit area and for this model the tension is σ = 8/9 mv 2 . Alternatively, if the potential is written as V = m 2 f 2 cos(ϕ/f ) then the domain wall solutions take the form ϕ = ± 4f arctan e mz − πf where δ w ∼ 1/m and σ = 8mf 2 .
For the clockwork axion model, we set f → f pq and m → m G to estimate the surface tension of the G-walls as and the wall thickness is δ w ∼ 1/m G . (Actually, since the gear spectrum is split (12) each wall will have a different tension and width, but this feature does not play a significant role in the dynamics, and we simply use Eq. (22) to estimate the tension for all of the walls.) Figure 3: For the model with three pseudoscalar fields (π 0 , π 1 , and π 2 ) we show the scalar potential over slices of the field space, and we overlay the trajectories corresponding to the π 0 -, π 1 -, and π 2strings. The potential is As we encircle a π 0 -string in spacetime, the fields pass along the horizontal white line (left panel) from the local minimum (green dot) to the saddle point (red cross), and back to the local minimum; therefore, the π 0 -string is connected to 1 domain wall. Similarly the π 2 -string, which corresponds to the vertical white line (left panel), connects to 3 domain walls, and the π 1 -string, which is shown on the right panel, connects to 4 domain walls. This discussion readily generalizes to N + 1 > 3. (The "local minima" are all connected by a flat direction, which is the axion field.) A given G i -wall in the network can either form a closed bubble, or it can terminate at one of the π n -strings. Moreover, the π n -strings may connect to multiple G-walls. Let N (n) DW denote the number of G-walls connected to a π n -string. From the structure of the interactions in Eq. (4), we see that This situation is illustrated in Fig. 3, and the caption provides further explanation; see also Ref. [32].

Evolution of the wall network
Given the rich structure of the clockwork axion string-wall defect network, it is challenging to study its evolution, either analytically or numerically. However, we gain some understanding from simple physical arguments. Since strings are pulled by multiple walls with different tensions (22), we expect that the string-wall network will begin to collapse, even before the QCD epoch. For example, consider a particular string that is connected to several walls. It is energetically preferable for the string to move in a direction that shortens the wall with the largest tension. In other words, the high-tension walls start to pull together the pair of strings that they are connecting. This attractive force accelerates the string segment, which causes it to radiate as we discussed in Sec. 3.2. The energy loss into radiation causes the string-wall system to shrink. Eventually the strings meet, and they either annihilate (if they had opposite winding number) or they merge into a new hybrid string with a composition of the winding numbers of the constituent π-strings. In this way, we expect that the string-wall network will partially collapse. This dynamical behavior has recently been observed in numerical simulations of aligned (clockwork) axion models with 2 and 3 complex scalar fields. For the model with two scalar fields (N + 1 = 2), the simulations of Ref. [32] show that the string-wall network completely collapses to form a network of hybrid strings, which are precisely the a-string that we discussed in Sec. 3.1. For model with three scalar fields (N + 1 = 3), the simulations suggest that a partial collapse will occur in which the π 0 -strings, which connect to only a single wall, are merged with the π 1 -strings to form hybrid strings, but these hybrid strings do not merge with the π 2 -strings.
It is difficult to assess how far does the partial collapse proceed before the system reaches a stationary scaling behavior. Surely the network cannot completely collapse for models with many fields, N 1, because this would require a network of a-strings to form, but we have already argued in Sec. 3.2.2 that this does not occur. Therefore, we expect that the partial collapse will suppress the energy carried by π n -strings and their connected walls with low n few, and we assume that the network evolution is unaffected by the partial collapse for n 1. Since our primary interested is to calculate the axion relic abundance, we do not expect more than a factor of O(N ) ∼ 10 uncertainty.
After the partial collapse has occurred, we expect that the subsequent evolution of the stringwall network is similar to the string-wall network evolution in a QCD axion model with N DW > 1 [35]. Specifically, we assume that the network reaches the scaling regime where new walls are created from the "chopping up" of larger walls. In the scaling regime, we estimate the energy density of the wall network at time t as [18] where σ i ∼ m G i f 2 pq is the tension of the G i -wall (22), and H(t) is the Hubble parameter. Since the energy of a wall is proportional to its area, the energy density of the network is dominated by the largest Hubble-scale walls, which makes this estimate more robust.

Axion emission from walls
In this section we estimate the energy density of axions that are emitted from the network of G-walls. A formalism for calculating particle emission from domain walls has been developed in Ref. [20], and we loosely follow that approach here. 3 (The reader may also recall the calculation of electromagnetic radiation from an oscillating sheet of charge, which follows a similar approach.) We first identify the interaction between the wall-forming field and the radiation, which correspond to the gear G i and the axion a. Next we assume a profile for the G i and solve the field equation for a to determine the radiated energy.
Let us inspect the Lagrangian (4) to identify the axion-gear coupling. Since the axion has an exact shift symmetry before the QCD epoch, interactions with the other pseudoscalar particles can arise in general only through higher derivative operators. By integrating out the n th radial mode, ρ n , one generates the operator (∂π n ) 4 . We write this operator in the mass basis using Eq. (10), and we focus on the single-axion term, which is The coefficient 4 is defined by b n,ijk ≡ 4O n0 O ni O nj O nk , and using Eq. (11) we estimate b n,ijk ∼ q −n N −3/2 . In the approximation where the single axion term is the dominant one, the axion field equation becomes In the vicinity of a G-wall the right side of Eq. (26) becomes a source for the axion field, and the associated axion radiation is given by the solution evaluated far away from the wall. In order to solve for the axion emission we must first know how the gear fields are evolving. In principle a numerical simulation of the string-wall network evolution can measure G i (x) and solve Eq. (26) directly to determine the amount of axion emission. Since this technique is not available to us, we take a semi-analytical approach instead, and we comment on the main sources of uncertainties at the end.
Without detailed information on the network evolution, we need to make an ansatz for the trajectory of the domain wall; this is an essential ingredient, because a static domain wall does not radiate. In the cosmological setting, the wall will be accelerated by the pull of its own tension and the tension of the strings to which it connects. Motivated by the stationary domain 3 The emission of axions from axion-walls has been considered previously in various references, including Refs. [36][37][38][39], where it is generally assumed that an O(1) fraction of the energy density carried by the domain walls goes into non-relativistic axions. We cannot apply this assumption to calculate axion emission from the gear-walls, because we expect that most of the energy goes into gears rather than axions. 4 The term with two axion and two gear fields has the coefficient For n 1 this term is suppressed by an additional factor of q −n compared to the single axion term in Eq. (25). Even for n ∼ 1 we expect that the two axion term will give a subdominant contribution to the axion emission calculation [40].
wall solution that was discussed in Sec. 3.3.1, we parametrize an accelerated wall by writing 5 which represents a wall located at z = z 0 (t, x) with width δ w ∼ 1/m G . The wall is not planar, but rather it has oscillatory (standing wave) perturbations that can be written as where A ω (t) is the displacement amplitude at time t. In the vicinity of this domain wall, the axion field evolves subject to Parametrically, the effective coupling strength is λ G ∼ (N 3/2 f 2 pq m 2 ρ ) −1 since the n = 0 term dominates the sum over n.
Formally the solution of Eq. (29) is written as where ∆ R (x, x ) is the retarded Green's function. It is now a straightforward (but lengthy!) exercise to evaluate this integral. To determine the amount of axion radiation, we focus on the regime z → ∞ where we find 6 up to an O(1) numerical coefficient. (The field amplitude approaches a constant far from the wall, as in the case of electromagnetic radiation from a sheet of charge, and unlike the case of radiation from a point charge, which gives instead a familiar 1/r behavior.) As the wall oscillates it radiates axions with a energy flux T 03 , which has the units of power per area. The energy flux is estimated as where the angled brackets indicate time averaging over one oscillation period. We see that axion emission is more efficient for walls with high frequency oscillations, which follows from the derivative interaction between a and the G i .
The oscillating domain wall emits axions at the expense of reducing its own kinetic energy and decreasing the amplitude of its oscillations. The domain wall's kinetic energy per unit area is (1/2)σA 2 ω ω 2 , since A ω ω is the effective oscillation velocity and the tension σ ∼ m G f 2 pq is the mass per unit area. Axion emission leads to d dt which gives the evolution of A ω (t). Using the expression for T 03 from Eq. (32), the solution is where A ω (t 1 ) is the initial amplitude of the mode with frequency ω at an arbitrary time t pq ≤ t 1 < t. We have also defined the dimensionless parameter v ω ≡ A ω (t 1 )ω, which represents the initial speed of the oscillation with frequency ω. From Eq. (34) we see that the mode with frequency ω has a "lifetime" given by For t−t 1 < τ ω the mode amplitude is effectively constant, but for t−t 1 > τ ω the mode amplitude decreases like A ω ∼ 1/ √ t as axion emission damps the wall's oscillations. The initial oscillation spectrum, parametrized by A ω (t 1 ) or v ω = A ω (t 1 )ω, represents the largest source of uncertainty in our calculation. A calculation of this spectrum is not analytically tractable, given the chaotic and dynamic nature of the string-wall network. The two essential questions are: how does A ω (t 1 ) depend on ω and what it the magnitude of A ω (t 1 )? To address these questions, we simply assume that as modes on the domain wall enter the horizon, they have a displacement amplitude that is some fraction α < 1 of the Hubble scale at that time. In other words, the mode with frequency ω enters the horizon at time t ω such that In effect this assumption lets us treat the product v ω = A ω (t ω ) ω ∼ α ≤ 1 as a constant that is static and independent of ω. We do not attempt to calculate α from first principles, but we have no reason to expect that it should not be O(1). From Eq. (34) we identify a characteristic frequency, This characteristic frequency sets a smoothing length scale l * (t) = 1/ω * (t) that is important for understanding the structure of the domain wall and its axion emission. If we consider the spectrum of oscillations on the domain wall at time t, the high-frequency (small-scale) modes with ω > ω * (t) are absent, because they were damped away at earlier times due to axion emission, and the low-frequency (large-scale) modes with ω < ω * (t) are still present, because they do not yet emit axions efficiently. Therefore the modes with ω = ω * (t) will give the dominant contribution to the axion emission at time t. At the QCD epoch we have ω * (t QCD ) ∼ MeV/v ω for N = 15, m ρ = 10 TeV, and m G = 100 GeV. These modes are much smaller than the horizon scale (H QCD /ω * (t QCD ) ∼ 10 −16 v ω ), but also much larger than the wall thickness (m G /ω * (t QCD ) ∼ 10 3 v ω ). Collectively the domain walls form a network that radiates axions with a power density P a (t), which has the units of power per volume. As we saw in Sec. 3.3.2, the horizon-scale domain walls carry most of the energy, and we can estimate the power density as P a (t) ∼ T 03 H. The energy density of radiated axions satisfieṡ ρ a + 4Hρ a = P a (38) where the term 4Hρ a appears because the axions are massless prior to the QCD epoch. We integrate the evolution equation for ρ a to find the energy density of emitted axions at the QCD epoch. During the radiation-dominated era we have H(t) = 1/(2t), and the solution is where ω * is evaluated at t QCD , and σ ∼ m G f 2 pq is the wall's surface tension. Eq. (39) is a main result of this work.
Inspecting Eq. (39), we see that the quantity in square brackets peaks at ω/ω * 1.5 where its value is 0.22. If the wall is oscillating with a low frequency, ω ω * (t QCD ), then there has not yet been enough time for axion emission to occur, and ρ a is suppressed like (ω/ω * ) 2 . On the other hand, if ω ω * (t QCD ) then the axion emission has occurred long before the QCD epoch, and redshifting has diluted the axion energy density leading to a suppression of (ω/ω * ) −2 .
In the Introduction we asserted that ρ a (t QCD ) ∼ ε m G f 2 pq H QCD (3), and now by comparing with Eq. (39) we identify the efficiency factor as where we have used σ ∼ m G f 2 pq . For the usual QCD axion, the dark matter relic abundance is dominated by axion emission from Hubble-scale domain wall oscillations at the QCD epoch. To evaluate this contribution in the clockwork axion model, we evaluate the efficiency factor with ω ∼ H QCD to find that ε is suppressed by (H QCD /m G ) 1. Therefore we can anticipate the result that appears in Sec. 3.5; namely, the axion dark matter emission from oscillating gearfield domain walls is negligible. On the other hand, in the formula for energy flux (32) we have seen that high frequency oscillations on the domain wall lead to a larger axion emission. For the modes with ω ∼ ω * (t QCD ) we find ρ a (t QCD ) ∼ 0.1v 2 ω σH QCD , which corresponds to extracting an O(0.1) fraction of the wall's kinetic energy as axion emission. However these axions are highly boosted since E a (t QCD ) ∼ ω * (t QCD ) H QCD , and they remain relativistic during the epoch of recombination. Therefore they should be treated as a dark radiation, rather than a dark matter.

Collapse of the string-wall network
The axion mass arises at the QCD epoch when instanton effects induce a potential for π N ; see Eq. (14). Provided that N DW = 1, the flat directions are lifted leaving only a single, unique vacuum. As a result the string-wall network collapses, and the energy is liberated as mostly gears and axions. Here we estimate the amount of energy that goes into the axion.
The instanton-induced axion potential causes a mixing between the gear fields and the axion. The mixing with the i th gear is estimated as for i ∈ {1, 2, · · · , N }. Despite the enhancement from q N 1, this factor is typically very small due to m 2 a m 2 G . The mixing (41) controls the efficiency with which the defect network can emit axions. The total energy density in the defect network at the QCD epoch is dominated by the Hubble-scale domain walls, which we estimate with Eq. (24) to be ∼ σH QCD where σ ∼ m G f 2 pq is the surface tension. Therefore the axion energy density that arises from the collapsing string-wall network at the QCD epoch is estimated to be Here we have used Eq. (13) to write q N f pq ∼ f a . Since the usual QCD axion gives ρ a ∼ m a f 2 a H QCD , we see that the axion emission here is suppressed by (m a /m G ) 3 . The axions emitted from Hubble-scale walls are expected to have an energy E a ∼ H QCD so that they become nonrelativistic soon after the QCD epoch [41]. The remainder of the energy is emitted mostly into gears G i that decay to Standard Model radiation through their interaction with gluons (14). The associated entropy injection heats the SM plasma momentarily, but since the energy of the defect network is negligible compared to the energy of the plasma, the heating is not significant.

Axion relic abundance
In this section we calculate the relic abundance of dark matter and dark radiation in the clockwork axion model. In the subsections, we enumerate the various contributions to ρ a (t QCD ).
The population of nonrelativistic axions contributes to the axion dark matter relic abundance, which is parametrized by Ω a h 2 = ρ a (t 0 )/3M 2 plH 2 0 whereH 0 ≡ 100 km/sec/Mpc and ρ a (t 0 ) is the energy density of nonrelativistic axions today. Between the QCD epoch and today the energy density of nonrelativistic axions redshifts as ρ a (t 0 ) = ρ a (t QCD ) a(t 0 )/a(t QCD ) −3 . Assuming adiabatic expansion of the cosmological plasma from the QCD epoch until today, the axion relic abundance is written as where g * S (t QCD ) 20, T QCD 0.2 GeV, g * S (t 0 ) 3.91, and T 0 0.234 meV. The measured dark matter relic abundance is Ω DM h 2 0.12 [42]. The population of relativistic axions contributes to the axion dark radiation relic abundance, which is parametrized by ∆N eff (t) = ρ a (t)/[2(7/8)(π 2 /30)(4/11) 4/3 T (t) 4 ]. The presence of a dark radiation during the epoch of recombination (t = t rec ) is strongly constrained by observations of the cosmic microwave background, ∆N eff 0.1 [42]. Assuming adiabatic expansion from the QCD epoch until recombination, we have where g * S (t QCD ) 20, T QCD 0.2 GeV, g * S (t rec ) 3.91, and T rec 0.3 eV.

Axion relics from cosmic strings
Eq. (21) gives the energy density of axions emitted by the network of π-strings from the PQ phase transition until the QCD epoch. As we have discussed in Sec. 3.2, the axions emitted closer to the PQ epoch remain relativistic at late times, whereas the axions emitted later become nonrelativistic soon after the QCD epoch. Both population are roughly equally abundant, which is the source of the logarithmic enhancement in Eq. (21), and therefore we use ρ a (t QCD ) ∼ µH 2 QCD log t QCD /t pq in Eqs. (43) and (44) to estimate As we discussed in Sec. 3.2.4, the tension of the π-strings is smaller than the usual axion string tension by the clockworking factor, (f pq /f a ) 2 ∼ q −2N 1. Since the tension is smaller, there is less energy available for axion emission, and the relic abundance from strings is suppressed. From these formulas we infer that the clockwork axion model predicts too much dark matter (Ω a h 2 > Ω DM h 2 0.12) for f pq > 9 × 10 13 GeV, and it predicts too much dark radiation (∆N eff > 0.1) for f pq > 2 × 10 17 GeV. However, the calculation breaks down for f pq > f a , since there is no clockworking, and we will see that f a < 10 11 GeV is required to avoid overproducing axion dark matter from misalignment in Sec. 3.5.4.

Axion dark radiation from domain walls
Eq. (39) gives the energy density of axions emitted by the network of G-walls from the PQ phase transition until the QCD epoch. We have discussed in Sec. 3.3 that domain wall oscillations with ω ∼ H QCD produce nonrelativistic axion dark matter that contributes to Ω a h 2 whereas oscillations with ω ∼ ω * (t QCD ) H QCD produce relativistic axion dark radiation that contributes to ∆N eff . Using Eqs. (39) and (40) in Eqs. (43) and (44) we estimate (48) As we anticipated in Sec. 3.3, the axion dark matter relic abundance arising from oscillation of Hubble-scale gear-field domain walls is very suppressed as compared with the usual QCD axion cosmology, and we will see that it is negligible when compared with the dark matter produced through misalignment. On the other hand, the predicted dark radiation in Eq. (48) becomes comparable to the observational upper limit, ∆N eff < 0.1, for sufficiently large m G and f pq . This translates into a constraint on the (m G , f pq ) parameter space appearing in Fig. 4. Recall that v ω ≤ 1 parametrizes the spectrum of perturbations on the domain wall on small scales. Due to the complicated structure and dynamics of the string-wall network, which we have discussed in Sec. 3.3, it is challenging to calculate the axion emission precisely. The axion relic abundance may be suppressed compared to these estimates if more energy is lost in other ways, for instance by the emission of gears.

Axion dark matter from string-wall collapse
Eq. (42) gives the energy density of axions emitted during the collapse of the string-wall network at the QCD epoch. Since this axions are predominantly produced from Hubble-scale defects, they are produced with very low momentum and contribute to the dark matter. The predicted dark matter relic abundance is given by Here we see that the axion production during the collapse of the string-wall network is entirely negligible, whereas for the usual QCD axion model this is one of the dominant contributions to the axion dark matter relic abundance.

Axion dark matter from misalignment
Since the axion mass is negligible at the PQ phase transition, the axion field is randomly and uniformly distributed across the vacuum manifold, a ∈ [0, 2πf a ) or θ a = a/f a ∈ [0, 2π). (Here we assume N DW = 1 so that f g = f a .) In general the local value of θ a at some point in space will be misaligned with the value of θ a that minimizes the axion potential, which arises at the QCD In the orange shaded region there is too much axion dark matter (Ω a h 2 > Ω DM h 2 0.12) produced from misalignment with f pq = f a . In this work we have calculated the axion emission from oscillating gear-field domain walls (48), and in the blue region we find that this produces too much dark radiation (∆N eff > 0.1). Uncertainties in this calculation are large and difficult to quantify. Assuming efficient axion emission from small-scale oscillating structure on the gear-field domain walls (v ω ∼ 1), everything above the blue-dotted line is excluded.
phase transition. This misalignment [43][44][45] corresponds to a local potential energy density of approximately ρ a = m a (t) 2 f 2 a (1 − cos θ a ) where m a (t) is the effective axion mass at time t, which grows rapidly during the QCD phase transition. Subsequently the axion field begins to oscillate and behave like pressureless dust. The corresponding relic abundance of axion dark matter is given by [46] Ω a h 2 0.2 f a 10 11 GeV where the value of θ a is averaged. Note that it is f a and not f pq that controls the misalignment contribution to the axion dark matter relic abundance. To avoid producing too much dark matter (Ω a h 2 > Ω DM h 2 0.12), it is necessary that f a 1 × 10 11 GeV, which translates into a constraint on the parameter space appearing in Fig. 4.  Figure 5: The canonical parameter space of QCD axion cosmology (left), and how this parameter space is modified due to the clockworking (right). In the right panel we fix the clockworking factor, q N , and vary f pq ∼ f a /q N . The diagonal boundary corresponds to f pq ∼ H I /2π. Clockworking widens the classic axion window, f a ∈ (10 9 , 10 11 ) GeV, allowing lower H I . The various shaded regions are excluded: the yellow region is excluded by CMB limits on isocurvature [11]; the blue region is excluded by limits on white dwarf cooling time [9] (observations constrain the axion-photon and axion-electron couplings, and the mapping to axion-gluon coupling is model-dependent); the gray region is excluded by CMB limits on inflationary gravitational waves (r < 0.1) [42]; and the red region is excluded by measurements of the dark matter relic abundance (assuming the axion is cosmologically stable).

Conclusion
In this work we have explored some of the cosmological implications of the clockwork axion with a focus on the network of topological defects and their contribution to the axion relic abundance. The primary cosmological consequence of clockworking, which lowers f pq compared to f a , is to make it easier for the PQ symmetry to be restored when the universe reheats after inflation. As shown in Fig. 5, the clockworking widens the "classic axion window" in which the PQ symmetry becomes broken during a cosmological phase transition, and a network of topological defects is formed. This behavior is well known from the usual QCD axion cosmology where a network of cosmic strings forms at the time of the PQ phase transition, and the strings become bounded by axion-field domain walls at the QCD epoch. In the clockwork axion model, on the other hand, the PQ phase transition gives rise to both O(N ) "flavors" of cosmic strings as well as a network of domain walls composed of the new pseudoscalar gear fields; see the discussion in Sec. 3.1.
The string-wall defect network emits axions as it evolves from the PQ phase transition until it eventually collapses at the QCD epoch. Given the rich structure of the network, it is challenging to make robust predictions for its dynamics. Assuming that the defect network reaches the socalled scaling regime, we calculate the relic abundance of relativistic axions as dark radiation and nonrelativistic axions as dark matter. Our main results are given by the formulas in Sec. 3.5 as well as Fig. 4, and we summarize the key points here.
• Clockworking lowers the tension of the cosmic strings from µ ∼ f 2 a to µ ∼ f 2 pq . Conse-quently the cosmic strings carry less energy, and the axion relic abundance arising from string emission is suppressed by (f pq /f a ) 2 ∼ q −2N . (This suppression is partially offset by a competing enhancement factor, because strings loops that emit less efficiently will live longer, but the net effect still leaves a negligible axion relic abundance.) • In the clockwork axion model, domain walls composed of the gear fields form already at the time of PQ symmetry breaking. The surface tension of these walls is σ ∼ m G f 2 pq , which can be comparable to the tension of axion-field walls in the usual axion cosmology, σ ∼ m a f 2 a , because m G m a even though f pq f a . This comparison suggests that the relic abundance of axion dark matter can be comparable to the usual estimates (3). However, the axion interacts only very weakly with the domain walls through the dimension-8 operator, (∂a)(∂G) 3 /f 2 pq m 2 ρ . Consequently the relic abundance of axion dark matter arising from Hubble-scale domain walls at the QCD epoch, given by Eq. (47), is suppressed by (f pq /f a ) 2 (m G /m ρ ) 4 (H QCD /m a ) compared to the usual axion cosmology, which makes it totally negligible. However, the presence of small-scale oscillating features on the domain walls produces a population of relativistic axions that contribute to dark radiation; see Eq. (48) and the blue region in Fig. 4. This calculation carries a large uncertainty associated with the spectrum of oscillations on the domain wall, which would be interesting to measure with a lattice simulation of the defect network evolution.
• The defect network collapses at the QCD epoch when the axion mass arises from instanton effects. At this time the axion potential induces a mixing between the axion and the gear fields, which provides a new avenue for axion emission. However, the mixing is related to the ratio of the axion mass and the gear mass, and the resultant axion relic abundance is suppressed by (m a /m G ) 3 , which makes it negligibly small. Most of the energy from the defect network collapse goes into emission of gear-field particles, which decay to Standard Model radiation.
In this article we have focused on the "simplest" example of a clockwork axion model [15] that is reviewed in Sec. 2. Our results are readily generalized to other aligned axion models with a large hierarchy between f pq and f a . Our key assumptions have been that (1) f pq is low enough that the PQ symmetry is restored during reheating after inflation and subsequently broken during a cosmological phase transition, (2) the complex structure of the defect network prohibits the formation of strings formed from the axion field, (3) the axion interacts weakly with the (gear) fields that form the domain walls, and (4) these interactions can be described by low energy effective field theory to calculate axion emission from the defects. More generally, in Sec. 2.1 we discussed various UV extensions of the clockwork axion following Ref. [30]. In general there can be a delay between the formation of the string network at the PQ phase transition and the later formation of the domain walls. Our results for the axion relic abundance are insensitive to the UV physics, because most of the axion production is occurring at late times, at the QCD epoch.