Axion dark matter and the Lattice

First I will review the QCD theta problem and the Peccei-Quinn solution, with its new particle, the axion. I will review the possibility of the axion as dark matter. If PQ symmetry was restored at some point in the hot early Universe, it should be possible to make a definite prediction for the axion mass if it constitutes the Dark Matter. I will describe progress on one issue needed to make this prediction -- the dynamics of axionic string-wall networks and how they produce axions. Then I will discuss the sensitivity of the calculation to the high temperature QCD topological susceptibility. My emphasis is on what temperature range is important, and what level of precision is needed.


Overview
The axion [1][2][3][4] is a proposed particle, the angular excitation of a new "Peccei-Quinn" (PQ) field ϕ that would solve the strong CP problem [5][6][7] and which is also a very interesting dark matter candidate [8][9][10], thereby solving two puzzles with one mechanism. That's why I think it well motivated to study the axion as a dark matter candidate. The axion model has one undetermined parameter, the vacuum value of ϕ, f a ; the axion mass m a scales as f −1 a . The value of f a also plays a role in determining the amount of axion dark matter produced in the early Universe. So do some nontrivial dynamics which we will explain in detail below. If we can understand the nontrivial dynamics of the axion field during cosmology, that lets us find a fixed relation between f a and the (measured) dark matter abundance. It therefore allows a clean determination of the axion mass, under some simple assumptions. This is valuable in the experimental search for the axion and it motivates us to solve the cosmological axion dynamics.
Supposing we make the following assumptions: • The axion exists.
• The axion field starts out "random" (in a sense we will define precisely below) either during or shortly after inflation (or whatever physics featured in the Universe at very high energy density).
• Gravity and the Universe's energy budget followed the "standard" picture (General Relativity and the known Standard Model species dominating the energy density) at temperatures T < 2 GeV.
1. We need to know the temperature dependence of the QCD topological susceptibility at temperatures between 540 and 1150 MeV.
2. We need to control the axion field dynamics during this temperature range, to solve for the efficiency of axion production.
The session following this talk addresses item 1 and the remainder of this talk will lay out the groundwork more completely and will then address item 2.

Strong CP problem and the axion
Let me refresh your memories on the strong CP problem. There are two gauge-invariant dimension-4 scalars which can enter the gauge-field part of the QCD Lagrangian: The latter term is P and T odd because it contains the antisymmetric tensor µναβ . The operator which Θ multiplies is the topological density and it integrates to the instanton number. Its zero-momentum two-point function defines the topological susceptibility, . . . T means the expectation value in the thermal ensemble at temperature T . Such a term is strongly constrained by the absence of a measured neutron electric dipole moment. The experimental limit of [11] |d n,meas | < 2.9 × 10 −26 e cm (3) contradicts the lattice results for the dipole moment from Θ [12] d n = −3.8 × 10 −16 e cm × Θ unless |Θ| < 10 −10 . At this conference we saw new results for the lattice Θ-dependent dipole moment which show that the above result may be too high and its error bar was certainly underestimated (see [13,14] and these proceedings). However it is clear that the absence of a neutron electric dipole moment places an extremely tight constraint on Θ. This is hard to understand because we know P and T are not fundamental symmetries; and any physics at a high scale which violates them generically gives rise to a Θ which does not decrease as we move to lower scales. For instance, consider a very heavy Dirac quark species Q, with Lagrangian Note that m can be complex, since (QP L Q) † =QP R Q. But an imaginary part is T and P odd, since the role of left and right projector, P L and P R , switch under parity and because T is antiunitary. We can remove this mass through a chiral rotation of Q, at the cost of reintroducing it, via the Fujikawa mechanism [15,16], as a shift in the value of Θ, 2 EPJ Web of Conferences 175, 01009 (2018) https://doi.org/10.1051/epjconf/201817501009 Lattice 2017 Therefore even a very heavy quark can influence the P and T symmetry properties of low energy QCD.
But what if there is a symmetry forbidding the mass term for this quark? For instance, suppose P L Q is charge 1 and P R Q is charge-0 under some global U(1) symmetry? Then the mass term breaks this symmetry, but a complex scalar ϕ with charge 1 under the symmetry could induce a mass via a Yukawa interaction and a vacuum value. The possible Lagrangian terms for such a scalar are The combination yϕ plays the role of m * in the previous case. But now Arg ϕ is a dynamical quantity. We will be interested in temperatures around 1 GeV and ϕ varying on scales of the Hubble scale at that time -tens of meters! Therefore from the point of view of QCD we can take ϕ to be spaceindependent, and perform an Arg ϕ dependent rotation on Q, making the theta term Here we have absorbed the phase in y into a phase redefinition of ϕ. We can also absorb Θ in the same way, so that Arg ϕ alone plays the role of Θ-angle. For notational compactness we will henceforth write ϕ = ve iθ a , with θ a = Arg ϕ. This specific way of coupling QCD topology to a complex scalar is called the KSVZ axion [17,18]. There are other mechanisms but the low-energy phenomenology is essentially identical and this mechanism is particularly clear to understand. From the point of view of QCD, the Θ-angle is replaced by a possibly spacetime-varying dynamical field θ a . What about from the point of view of the field ϕ? Since we want physics on the meter length scale, we can integrate out QCD, leading to an effective potential: with Ω the volume of space included in the path integration. In the second line we have made a dilute instanton approximation, which is that the integration exponentiates over the two-point function of the topological density, controlled by the topological susceptibility χ(T ) introduced already in Eq. (2). This is not a good approximation for large θ a and low temperatures [19], but it works well when instantons are dilute, which is true for T > 500 MeV, and for small values of θ a , which will be all we encounter below this temperature. So we can actually use this approximation all the time. Independent of this approximation, it is easy to see that the effective potential is smallest (the θ a choice is most energetically favored) for θ a = 0 and therefore when P and T symmetry are restored. Note that θ a only enters as the coefficient in a complex phase, in an otherwise real and positive integral. The integral is maximized, and the free energy minimized, if the phase is always unity. Any nonzero value of θ a gives rise to phase cancellations and therefore suppresses the partition function, raising the free energy. Although we derived it in Euclidean space, we can also use this effective potential in Minkowski space to study the spacetime evolution of the ϕ field. In summary, the Minkowski effective Lagrangian for the ϕ field is We will use this to determine the dynamics of the field in the next sections.

Axion in Cosmology
Let us see what happens to the axion field during cosmological evolution.  Figure 1. Cartoon of what we know about the topological susceptibility as a function of temperature (log-log plot). At low temperatures, chiral perturbation theory gives us the zero-temperature limit and small-temperature behavior. At large temperatures, perturbation theory predicts a power-law falloff with a power near 8.
The form of Eq. (11) makes it clear that, in order to study the axion's role in cosmology, we are going to need to know the temperature dependence of the topological susceptibility χ(T ). It does not yet tell us what temperature range will be interesting. Figure 1 shows our knowledge at the cartoon level. At low temperature or vacuum, chiral perturbation theory works and [19] At high temperatures we have conventional perturbation theory, which forecasts [20] that χ(T ) ∝ T −7−N f /3 . However the exact coefficient is sensitive to the physics of electric screening and is not known accurately. This is why we need lattice results for this quantity! Recently there have been several [21][22][23][24][25][26][27], which give generally compatible results but generally at temperatures below 600 MeV (or in the quenched approximation). It takes new techniques to reach higher temperatures, and only one recent paper [27] achieves this, reaching temperatures of 1500 MeV. Also, one group [28,29] finds results which are discrepant with the others, indicating that the matter is not yet settled. Here we will assume that the results of [27] are correct. This may well be the case, but we leave the discussion of the relative merits of these approaches and results to the panel, who have more expertise. Needless to say it would be valuable to know definitively that χ(T ) is well determined.

Space-uniform axion field
So let's assume for now that we know χ(T ). For simplicity let us also assume that the axion takes the same value everywhere in space, θ a (x, t) = θ a (t). It is simplest to work in terms of conformal time, so the metric is g µν = a 2 (t)η µν with a the scale factor. (Later we will use a to represent the lattice spacing. This is actually the same thing, since we will work in comoving coordinates; the lattice spacing is proportional to the scale factor and we may as well use a proportionality of 1.) In 4 EPJ Web of Conferences 175, 01009 (2018) https://doi.org/10.1051/epjconf/201817501009 Lattice 2017 the radiation era a(t) ∝ t and T ∝ t −1 . The radial component of ϕ = ve iθ a is inactive, v = f a , and the angular part obeys whereχ(t −1 ) is a rescaled form of the susceptibility. This leads to damped, anharmonic oscillations. The oscillations start roughly at the time t * when t 2 * χ (t −1 * ) = t −2 * , or in physical units, when m a ≡ χ(T )/ f 2 a obeys m a t * = 1 or equivalently m a /H = 1 with H the Hubble scale. After this time the oscillations accelerate as t 2χ increases, and they damp away. The damping arises both from the 2∂ t θ a /t term (Hubble drag) and from the time variation of the susceptibility. After several oscillations the axion particle number becomes an approximate adiabatic invariant, with number density parametrically of form (t * /t) 2 ( f 2 a /t * ). We see that the number density is quadratic in f a , while the axion mass is m a ∝ f −1 a . Because χ(T ) is a very strong function of temperature, t * depends only weakly on f a , and so the generated axion energy density is almost linear in f a .
Therefore, the larger the value of f a , the larger the produced axion abundance. However the axion abundance also depends on the unknown initial angle θ a (t = 0). Therefore the dark matter density depends on two variables and it is impossible to make a clean prediction for the value of f a . We can make a baseline prediction, however, by averaging over the value of the starting angle θ a (t = 0). Doing so, one finds the axion mass should be 32 µeV, and t * corresponds to a temperature of T * = 1.6 GeV.

Space-random axion field
It is far more likely that the Universe started out with a spatially random value for θ a , with no correlations on scales longer than the Hubble scale. Arguments for this picture are presented in [30] and are summarized as follows: • It is likely that inflation occurs with a high scale, H 2 > f 2 a /60. In this case, over 60 efoldings of inflation, quantum fluctuations stretched (squeezed) by inflation into classical fluctuations would randomize the value of the axion field over the course of inflation. The observation of cosmological tensor modes would more-or-less settle this issue.
• After inflation, the Universe reheats to a temperature which can be as high as T rh ∼ 0.1 Hm pl .
Even if H f a , if the reheat temperature is T rh > f a , there would be thermal symmetry restoration for ϕ. Then when the temperature falls below this scale, ϕ would independently take on a vacuum value at different points in space, which would be uncorrelated.
• The case where inflation and reheating are both low-scale is actually tightly constrained by the absence of observed isocurvature fluctuations (different fluctuations in θ a than in the radiation temperature), which require roughly H < 10 −5 f a . Most inflation model-builders would consider this rather unlikely.
I emphasize that we do not know that θ a was randomized in the early Universe (assuming the axion exists). But it appears likely, and it motivates studying the consequences. I will also assume that the axion makes up all of the dark matter in the universe, so we may equate the final axion matter density with the dark matter density, which is known to obey ρ dm /s = 0.39 eV with s the entropy density [31].
The space-inhomogeneous case is much more complicated than the space-homogeneous case. Nevertheless, in the remainder of the presentation I will show how to solve it.

Axion string/wall network
The ϕ field varies with amplitude of order f a ∼ 10 11 GeV over a length scale controlled by H ∼ T 2 /m pl ∼ 10 −18 GeV. This huge hierarchy makes the dynamics those of a classical field to extremely high accuracy. The Lagrangian Eq. (11) (times t 2 to account for Hubble expansion) and resulting classical equations of motion are easy to put on the lattice and solve as a function of time, from random initial conditions. In broad brushstrokes, our approach is to do just this, evolving the system until only small fluctuations in θ a remain and their evolution has become adiabatic. Then we integrate the associated axion number, and compare it to the result of the angle-averaged misalignment baseline.
In fact such a simulation is not sufficient, because of the large hierarchy in Eq. (11) between the mass scale m ∼ 10 11 GeV of radial excitations and the mass scale m a = χ(T )/ f a ∼ H ∼ 10 −18 GeV of angular fluctuations. The simulations have to take place at the m a scale, which means that the radial-mass scale cannot be resolved. Naively this should not matter, as radial excitations should decouple. But it does matter, because the theory contains topological string defects which play a role in the dynamics, and the string tension depends logarithmically on the ratio m/m a . Let's explain this in a little more detail.

String defects
First note that θ a is only defined modulo 2π. Therefore in traversing a circle, θ a might return to its starting value, but it might only return modulo 2π, that is, ∂ i θ a dx i = 2πN. The integer N is a winding number which counts a "flux" of string defects through the circle. If we deform a loop, N can only change when the loop passes through a singularity in the θ a field. The locus of these singularities defines the axionic cosmic string. Figure 2. A 2D slice of a simulation, with the complex field ϕ represented by a field of arrows with a length and direction. Going around the blue circle, the arrow direction winds by 2π. The associated string defect is at the center of the circle; another string defect is farther down and to the right.
We illustrate the idea with Figure 2, which shows a 2D slice out of a simulation, representing the complex field as a field of arrows with length and direction. The field direction has singularities where the arrows have zero length; going around the singularity, the direction of ϕ revolves by ±2π. The singular point extends in 3D into a line where the field has zero value; any loop circling this line will have the direction of ϕ revolve by ±2π as the one circles around the string.
Such a defect -essentially a vortex in the ϕ field -is called an axionic cosmic string, and it is topologically stable; no local changes to the value of ϕ can cause it to disappear. If PQ symmetry is restored in the early Universe, then θ a starts out uncorrelated at widely separated points and will generically begin with a dense network of these strings (the Kibble mechanism for string production [32]). The strings evolve, straightening out, chopping off loops, and otherwise reducing their density, arriving at a scaling solution [33] where the length of string per unit volume scales with time t as t −2 . They may play a dominant role in establishing axion production in the scenario under discussion [34].
Let us analyze the structure of a string in a little more detail. Consider a straight string along the z axis; in polar (z, r, φ) coordinates the string equations of motion are solved by √ 2ϕ = v(r) f a e iφ , with v(r) 1 for all r 1/m; so θ a = φ (up to a constant which we can remove by our choice of x-axis). The string's energy is dominated by the gradient energy due to the space variation of θ a : where the integral over r is cut off at small r by the scale where v(r) 1 (the string core), and at large distances by the scale where the string is not alone in the Universe but its field is modified by other strings or effects; this should be the larger of H and m a . We define κ = ln(m/H) as the log of this scale ratio. Now m is at most f a ∼ 10 11 GeV, and to ensure that the radial particles decay by the scale of 1 GeV we need m > 10 3 GeV. Therefore κ ∈ [48, 67]. This logarithm, κ, controls several aspects of the strings' dynamics. It controls the string tension, as we just saw. More relevant, while the string tension is πκ f 2 a , the string's interactions with the long-range ϕ field scale as f 2 a without the κ factor. Therefore the string's long-range interactions become less important, relative to the string evolution under tension, as κ gets larger. The long-range interactions are responsible for energy radiation from the strings, as well as for long-range, often attractive, interactions between strings. Since these effects tend to deplete and straighten out the string network, the large-κ theory will have denser, kinkier strings. Indeed, in the large κ limit the string behavior should go over to that of local (Nambu-Goto) strings [35]. Unfortunately, a numerical implementation must resolve the length scale m, ma ≤ 1, and cannot exceed m/H ∼ 1000; numerical studies of the scalar field system have κ < 7, nearly an order of magnitude too small.

Wall defects
Besides the strings, there are also wall defects. These occur late in the simulation when m a H. The potential term χ(T )(1 − cos θ a ) then forces θ a 0 nearly everywhere -modulo 2π. But suppose some region has θ a 0 and another has θ a 2π. There must be some 2D surface between them with θ a π. This is a wall defect. The region near the defect where θ a differs significantly away from its minimum has thickness ∼ 1/m a , which is easily resolved on the lattice. The surface tension of such a surface turns out to be 8m a f 2 a . A 2D slice of a configuration, illustrating such a domain wall attached to a string, is shown in Figure 3.
These wall defects are not a problem to simulate. But they play an important role in the dynamics. Every string has θ a take every value [0, 2π] as one goes around the string. That includes θ a = π. Therefore every string is attached to a domain wall. When m a becomes m a κH, the force from the domain wall tension becomes large enough to pull around the strings, leading to the collapse of the string network and the annihilation of all strings. It is only after this network collapse that one can speak about axion number. Because of the factor κ in the needed tension, a large-κ simulation will feature a more persistent string network. One final problem for scalar-only simulations, pointed out in [36], is that the domain walls actually lose even their metastability as soon as m 2 a /m 2 > 1/39. This drives up the required size of simulations so that large m can be achieved.

Simulating high-tension strings
We see from the previous section that simulations of the ϕ field alone are not reliable. Although one can make the scale m very heavy compared to H, m a , the string tension depends logarithmically on this scale, and is nearly a factor of 10 too small. This profoundly affects the dynamics of the string network, and therefore renders the results unreliable. We need a method to simulate high-tension strings coupled to θ a . We found such a method in [37] and present it here.

Effective theory
We are interested in the large-scale structure of string networks and the infrared behavior of any (pseudo)Goldstone modes they radiate. For these purposes it is not necessary to keep track of all physics down to the scale of the string core. Rather, it is sufficient to describe the desired IR behavior with an effective theory of the strings and the Goldstone modes around them. This consists of replacing the physics very close to the string core with an equivalent set of physics. It has long been known how to do this [35]. The string cores are described by the Nambu-Goto action [38][39][40], which describes the physics generated by the string tension arising close to the string core. The physics of the Goldstone mode is described by a Lagrangian containing the scalar field's phase. And they are  [41,42]: Here σ is an affine parameter describing the string's location y µ (σ, t), v µ = (1,ẏ) = dy µ /dt is the string velocity, and H µνα and A µν are the Kalb-Ramond field strength and tensor potential, which are a dual representation of θ a . Effectively L NG tracks the effects of the string tension, which we nameκπ f 2 a , stored locally along its length. Next, L GS says that the axion angle propagates under a free wave equation, as expected for a Goldstone boson, and its decay constant is f a . And L KR incorporates the interaction between strings and axions, also controlled by f a . The interaction can be summarized by saying that the string forces θ a to wind by 2π in going around the string (in the same sense that the eJ µ A µ interaction in electrodynamics can be summarized by saying that it enforces that the electric flux emerging from a charge is e).
It should be emphasized that in writing these equations, we are implicitly assuming a separation scale r min ; at larger distances from a string r > r min we consider ∇ϕ energy to be associated with θ a ; for r < r min the gradient energy is considered as part of the string tension [35], meaning thatκ incorporates all tension contributions from scales shorter than r min .
Any other set of UV physics which reduces to the effective description of Eq. (18) would present an equally valid way to study this string network. Our plan is to find a model without a large scale hierarchy, such that the IR behavior is also described by Eq. (18) with a large value for the string tension. Optimally, we want a model which is easy to simulate on the lattice with a spacing not much smaller than r min . Reading Eq. (19) through Eq. (23) in order, the model must have Goldstone bosons with a decay constant f a and strings with a large and tunable tension T str =κπ f 2 a , withκ 1. There can be other degrees of freedom, but only if they are very heavy (with mass m ∼ r −1 min ), and we will be interested in the limit that their mass goes to infinity. Finally, the string must have the correct Kalb-Ramond charge. Provided everything is derived from an action, this will be true if the Goldstone boson mode always winds by 2π around a loop which circles a string.

The model
We do this by writing down a model of two scalar fields ϕ 1 , ϕ 2 , each with a U(1) phase symmetry. A linear combination of the phases is gauged; specifically, the fields are given electrical charges q 1 ∈ Z and q 2 = q 1 −1 under a single U(1) gauge field. The orthogonal phase combination represents a global U(1) symmetry which will give rise to our Goldstone bosons. The role of the gauge symmetry will be to attach an abelian-Higgs string onto every global string, which will enhance the string tension. The added degrees of freedom are all massive off the string, achieving our intended effective description. The model falls under the general rubric of "frustrated cosmic strings" [43], but our motivation and some specifics (particularly our initial conditions) are different.

Specifically, the Lagrangian is
For simplicity we will specialize to the case The model has 6 degrees of freedom; two from each scalar and two from the gauge boson. Symmetry breaking, ϕ 1 = e iθ 1 v 1 √ 2 and ϕ 2 = e iθ 2 v 2 √ 2, spontaneously breaks both U(1) symmetries and leaves five massive and one massless degrees of freedom. Specifically, expanding about a vacuum configuration, the fluctuations and their masses are We see that the choices in Eq. (25) have made all heavy masses equal. 1 To clarify, note that a gauge transformation A µ → A µ + ∂ µ ω changes θ 1 → θ 1 + q 1 ω and θ 2 → θ 2 + q 2 ω. Therefore the linear combination of θ 1 , θ 2 fluctuations with δθ 1 ∝ q 1 and δθ 2 ∝ q 2 is precisely the combination which can be shifted into A µ by a gauge change, and is therefore the combination which is "eaten" by the A-field to become the third massive degree of freedom. The remaining phase difference q 2 θ 1 − q 1 θ 2 is gauge invariant, q 2 θ 1 − q 1 θ 2 → ω q 2 (θ 1 + q 1 ω) − q 1 (θ 2 + q 2 ω) = q 2 θ 1 − q 1 θ 2 + 0ω (31) and represents a global, Goldstone-boson mode.

The strings
We initialize ϕ 1 and ϕ 2 with the same space-random initial phase, which ensures that all strings will have each scalar wind by 2π, and the strings will have global charge q 1 − q 2 = 1. To find the tension of such a string, we write the Ansatz The value ofκ is thereforē   Figure 6 shows how the density of the string network depends on the string tension. The network density is expressed in terms of the dimensionless scaling variable ξ, We see that it increases by over a factor of 3 as one goes from scalar-only (κ 7) to high-tension (κ 57) simulations. It is not clear that our lattices are large enough to see the onset of scaling behavior for the highest-tension networks we studied. We certainly don't see scaling behavior for the abelian Higgs simulations, but this is another story.

Numerical implementation
I will not insult this audience by explaining how to implement a bosonic U(1) theory on the lattice. I use the noncompact formulation of U(1) and a next-nearest neighbor (a 2 ) improved action. There are no issues of renormalization of parameters because we are studying classical field theory. A subtlety in implementing electric fields with improvement is handled as in [44]. We pause only to mention one subtlety in how we implement χ(1 − cos θ a ). The function cos θ a , with θ a = q 2 θ 1 − q 1 θ 2 , is highly singular near the string core. Such singular behavior creates problems under space discretization, so we have to "round off" the behavior inside the string core, with the substitution The smoothing function F(r) is chosen such that F(1) = 1, F (1) = 0, F(0) = 0, and F (r) is continuous. This modification only affects the dynamics inside string cores, but the χ(t)(1 − cos θ a ) term is much smaller than the other potential terms there. Indeed, this term is everywhere small and it is only important because it operates over much of the lattice volume, while the radial potential has an effect only in the tiny fraction of the lattice corresponding to string cores.

Results
We studied this model [45] on 2048 3 lattices using a single compute node containing two Xeon Phi (KNL) processors. After systematic studies of lattice spacing and continuum limits, we also studied how the network evolution and the axion production depend on the string tension. The most notable result is that the axion production is actually smaller for the case with strings than the angle-average of the "misalignment" mechanism of Eq. (14). This contradicts the "conventional wisdom" [46] that axion production should be the sum of a misalignment contribution, a string contribution, and a wall contribution. We claim that this view double counts; the energy in domain walls is the energy of field misalignment, from values θ a ∼ π. This energy represents most of the potential axion production from misalignment. But the energy in the wall network is mostly absorbed by the strings when the walls pull on the strings, giving their energy to string velocity. Then the strings chop up into loops, which appear to produce axions quite inefficiently.
Combining our numbers with expressions for the energy budget g * and topological susceptibility from [27], and the observed dark matter density, we find m a = 26.2 ± 3.4 µeV.

Topology: temperature range
We have emphasized that an input value for χ(T ) from the lattice is essential to deriving these results. But we don't need χ(T ) at all temperatures; some are more important than others. We will extend the results of [45] by investigating over what temperature range the topological susceptibility is really needed.
At sufficiently high temperature χ(T ) is small. Therefore χ(T )(1 − cos θ a ) plays little role in the field dynamics. Its importance is controlled by the combination m a t which rises with time as approximately t 5.8 . Therefore at high temperature, large errors in χ(T ), or no value at all, is not a problem. To study this, we replace χ(T ) with the following "chopped" form: Technically it is t 2 χ(T ) we chop in this way, because of the t 2 factor in, eg, Eq. (14). We then study the axion production as a function of the temperature T chop . When n ax ceases to depend on T chop , we know that χ(T chop ) is not relevant. But so long as the result with T chop is different than the unchopped limit, we need χ(T ) at that temperature. We see from Fig.8 that above about 1150 MeV, it makes little difference if we change χ(T ). But below this value, we are quite sensitive. So it is necessary to determine χ(T ) all the way up to 1150 MeV.
On the IR side, perhaps surprisingly, there is also a range of temperatures where the susceptibility is not important. That is because the strings and walls are gone and the θ a value is small and oscillating rapidly, with m a H. Then the axion number is an approximate adiabatic invariant, which reacts smoothly to Hubble expansion and mass shifts and does not feel the anharmonicity of the potential, 1 − cos θ a θ 2 a /2. All that matters is the final value of χ(T → 0), which is well known. To find out where this temperature range starts, we make a similar modification: That is, we freeze χ(T ) (really, t 2 χ(T )) from rising after some point, and see if that changes the axion number produced. The results are shown in Fig. 9. They indicate that the susceptibility is irrelevant below about 550 MeV; the behavior above that temperature is important. Therefore the axion dynamics is sensitive to χ(T ) in the range 550 to 1150 MeV; outside of that range it is not.

Conclusions
The axion is a well motivated hypothetical particle, because it might explain two mysteries -the P and T symmetry of QCD, and the nature of the Dark Matter -with a single mechanism and particle. The physics of the axion in cosmology is rich, governed by a network of string defects which are swept together when domain walls form. Its explication requires two new pieces of physics. We need to understand this network evolution and axion production better. And we need to know the QCD topological susceptibility as a function of temperature, χ(T ), because it sets the tension of the domain walls and controls the physics which destroys the string network. I have presented the latest details on the network evolution, introducing a new technique which allows simulation of high-tension strings without excessive numerical resources. The results indicate rather inefficient axion production and therefore a rather light axion compared to previous studies, with m a = 26.2 ± 3.4 µeV.
It remains to form a consensus in the lattice community that the topological susceptibility is well measured. My work indicates that the susceptibility is needed in a temperature range from 550 to 1150 MeV. Below this range the evolution is adiabatic. Above this range the susceptibility does not yet play a role in axion dynamics. Clearly it is challenging to determine the susceptibility at such high temperatures. But I am confident this is a challenge which the lattice community will accept with gusto.