Axionic domain walls at Pulsar Timing Arrays: QCD bias and particle friction

The recent results from the Pulsar Timing Array (PTA) collaborations show the first evidence for the detection of a stochastic background of gravitational waves at the nHz frequencies. This discovery has profound implications for the physics of both the late and the early Universe. In fact, together with the interpretation in terms of supermassive black hole binaries, many sources in the early Universe can provide viable explanations as well. In this paper, we study the gravitational wave background sourced by a network of axion-like-particle (ALP) domain walls at temperatures around the QCD crossover, where the QCD-induced potential provides the necessary bias to annihilate the network. Remarkably, this implies a peak amplitude at frequencies around the sensitivity range of PTAs. We extend previous analysis by taking into account the unavoidable friction on the network stemming from the topological coupling of the ALP to QCD in terms of gluon and pion reflection off the domain walls at high and low temperatures, respectively. We identify the regions of parameter space where the network annihilates in the scaling regime ensuring compatibility with the PTA results, as well as those where friction can be important and a more detailed study around the QCD crossover is required.

A DW network in the early Universe can be a powerful source of gravitational waves, as shown quantitatively in numerical simulations of the corresponding field theory [70][71][72][73].The NANOGrav collaboration has in fact interpreted this signal as coming from a DW network in the scaling regime [1], identifying the ranges of temperature and energy density of the network compatible with the data.We will employ their results in our analysis.
DWs arise naturally in ALP models where the discrete subgroup of the U (1) which is preserved by the anomaly undergoes spontaneous symmetry breaking.These DWs can be topologically stable and eventually dominate the energy density of the Universe [49,96], in conflict with standard cosmology.This problem can be avoided by introducing a (small) explicit breaking of the discrete symmetry (the so-called bias), leading to the annihilation of the network [96,97].The size of the bias is generally a new independent parameter, which determines the phenomenology of the network and in particular the emitted SGWB.
As we shall discuss, a natural scenario to motivate the signal from ALP DWs at the PTA frequencies is the one where the ALP couples to QCD, as already suggested e.g. in [44,47,68,72,98], so that the QCD-induced potential itself acts as a bias.In this case, however, friction from the QCD sector in the thermal plasma is inevitable.Determining the impact of this friction force for the ALP domain wall interpretation of the PTA data is the main goal of our study (see [68,69,[99][100][101] for previous studies of friction effects on DW networks).In fact, when friction dominates the network departs from scaling and the corresponding SGWB can significantly change with respect to the one observed in numerical simulations [70][71][72][73] where plasma effects are not included.
Our main result is that friction acting on the DWs from the QCD sector can in general be relevant for a significant part of the model parameter space capable of reproducing the NANOGrav signal, and that a departure from the scaling regime (on which the PTA interpretation is based) is possible, even though a more detailed analysis around the QCD crossover is required to completely settle the issue.We also find parameter space compatible with the SGWB at PTAs where friction can be safely neglected.
Note that our study focuses on a minimal realization of ALP DWs with QCD-induced bias, where the heavy ALP cannot be the QCD axion.Nevertheless our analysis on friction effects can apply equally to non-minimal scenarios where the heavy axion is actually the one responsible to solve the strong CP problem, see e.g.[50,102].
Additionally, while the study presented in this paper focuses on the unavoidable coupling with QCD, friction could also impact the DW network if other ALP-SM (modeldependent) couplings are present.For instance, in Ref. [68] it was shown that if the ALP couples to muons1 the resulting friction can actually play a role at PTA frequencies provided that this interaction has the right strength.
The structure of the paper is as follows.In the next section we review the basics of ALP domain walls and SGWBs.Then, in Section 3 we show that ALP DWs whose bias is induced by QCD are a natural candidate to explain the PTA signal, by comparison with the model-independent analysis in [1].In Section 4 we study the impact of friction from the QCD sector, analyzing the contribution from gluons and pions at high and low temperatures.In section 5 we present our results showing the relevance of friction for the ALP parameter space.

ALP domain walls
We consider domain walls that generically arise in ALP models.We introduce our notation and setup.We start off by considering a dark QCD non abelian group SU (N ) and a U (1) Peccei-Quinn symmetry anomalous under SU (N ).The corresponding pseudo-Nambu-Goldstone boson, a, plays the role of the ALP in our analysis.The Lagrangian for the ALP contains the following term where α d is the dark gauge coupling constant, v is the U (1) breaking vev, and G d is the dark gauge boson field strength, which is contracted with its dual Gd .In terms of the quantities above it is useful to define the ALP decay constant f a as In terms of this quantity the Lagrangian is defined as usual as For N d = 1/2 the domain wall number is one and the vacuum manifold for the ALP potential induced by the dark gauge theory is trivial, namely it contains only one minimum as a = 0 and a = 2πf a = 2πv are to be identified.For N DW > 1 however the vacuum consists of disconnected points, and stable domain wall solutions exist interpolating between neighboring minima.
The most important features of the ALP potential are captured by the following structure which encodes the discrete Z N DW of the ALP Lagrangian with respect to the dark sector: where m a is the ALP mass.The ALP potential is defined in the range a ∈ [0, 2πv) = [0, 2πN DW f a ) and it then supports N DW − 1 degenerate and inequivalent minima.The simple cosine potential allows to obtain analytical domain wall solutions, with energy per unit surface (tension) given by In general, the ALP potential could differ significantly from the cosine shape if additional (pseudo) Nambu-Goldstone bosons exist below the dark QCD confinement scale, as it is the case for the QCD axion, see e.g.[103].In any case, the profiles above can still capture most of the relevant physics of the ALP domain wall solution2 .
Cosmology evolution of DW and SGWB spectrum We assume that the PQ breaking scale v is smaller than the reheating temperature and we focus on the case in which the domain wall number is larger than one, so that the DW network is stable [96,104]. 3At the scale of the PQ phase transition global strings associated to U (1) are formed according to the Kibble mechanism [48], see [109][110][111][112] for recent work.Subsequently, when the ALP potential becomes cosmologically relevant, that is when H(T f ) ∼ m a , domain walls can be considered formed.The precise relation to determine T f should take into account the temperature dependence of the axion potential generated by the dark-QCD sector.For an ALP decay constant below the Planck scale, one can show that T f √ m a f a , so we consider √ m a f a as an estimate of T f that will anyway play no important role in our study.
Soon after DW formation, the energy density of the resulting string-wall hybrid network is soon dominated by the walls [69,72].The DWs reach then the so-called scaling regime where the energy density of the DW network redshifts as ρ DW ∼ σH, corresponding to O(1) DWs per Hubble patch and mildly relativistic average velocity [113][114][115][116][117][118].
A scaling network of DWs will eventually dominate the energy density of the Universe, in contrast to standard cosmology [49,96,119].This occurs when ρ DW ∼ 3H 2 M 2 Pl , or in terms of temperature where we have assumed radiation domination, and we have used ρ DW = 2σHA with A = 0.8 from numerical simulations [70][71][72][73].
In order to collapse the DW network before domination, one may add a bias potential ∆V that breaks explicitly the Z N DW symmetry.We shall then define T * as the the annihilation temperature, where T * > T dom for consistency.
At the time of annihilation, the DW energy density normalized to radiation is given by This definition of α * is inspired by analogous studies of first order phase transitions, see e.g.[120], and as we shall see it is directly related to the strength of the GW emission.
The annihilation temperature may be estimated by balancing the curvature pressure with the energy difference induced by the bias, namely ∆V ∼ σ DW /R ∼ ρ DW , where R is the correlation length of the network.One finds where we used the condition ∆V = C ann ρ DW with C ann 2 from numerical simulations [70][71][72][73].The bias is in principle a free parameter that should be added to the model, and the phenomenology of the DW network can change drastically depending on its size, In general, one can expect this to come from quantum gravity effects that make the starting U (1) global symmetry only approximate [121][122][123][124][125][126][127][128].On the other hand, the size of the bias can be predicted if it is dynamically generated.In the next section we will explore the possibility that such bias is in fact generated by QCD.
The DW network in the scaling regime has been proven by numerical simulations [70,72,73] to generate a large SGWB Ω gw (f ) with broken power law in frequency.The signal is dominated by the last moment of emission, so it depends explicitly on T * .The signal redshifted today has the form (2.11) where we have normalized the numerical coefficient to the values obtained in numerical simulations [70][71][72][73].The SGWB spectrum of DW is determined by two parameters, the tension (see Eq. (2.7)) and T *4 .Eq. (2.11) shows that the later the DW network annihilates, the larger the GW signal is.As we can see, the best T * for PTA frequencies is in the ballpark of the QCD scale.
3 The QCD potential as the natural bias for DWs at PTAs Let us now consider the effect of the QCD-induced potential on the ALP model illustrated above.This comes from the anomalous coupling between the ALP and the gluons, where N c is the color anomaly from fermions charged under QCD.In general, N c and N d are two independent numbers.Whenever these numbers are not coprime, the degeneracy in the vacuum manifold is lifted and domain walls become metastable.
The contribution to the ALP potential from QCD at low energy can be captured within to chiral perturbation theory, see e.g.[103,130].One finds the following potential for the ALP-pion system: Notice that since f a = f a , the periodicity of the QCD potential is generally misaligned with respect to the one of the dark-QCD potential.The interactions in (3.2) follow from an ALP-dependent rephasing of the light up and down quarks that removes the ALP from the topological term, q → q exp(iγ 5 a 2fa Q a ) and Q a proportional to the identity (Tr Q a = 1).Assuming that the QCD contribution is very small compared to the dark QCD one (m π f π m a f a ), the size of the bias is generically given by |∆V k | ∼ m 2 π f 2 π .However, when the two sectors are almost aligned, for instance when N c /N d = 1 + with 1, the bias can become parameterically small: and the life time of the network may be parameterically enhanced.In our analysis we will keep as a free parameter, keeping in mind that a scenario with 1 requires somewhat large or fine-tuned values of N c and N d .
Let us now turn to discuss how temperature corrections modify the size of the QCDinduced ALP potential.At very high temperatures above QCD confinement, the ALP potential is expected to behave as [103,[130][131][132][133] with n 7 and χ 1/4 0 75.6 MeV, even though some uncertainty on these parameter still remains (see e.g.[133] and references therein).Similarly to the low-temperature case, some approximate alignment between QCD and the dark QCD can lead to a parametric suppression of the natural bias From this simple estimate we can already draw some conclusions in the light of the recent PTA results.In Ref. [1], the collaboration has performed a bayesian analysis on the NANOGrav data for the DW interpretation.The results were displayed in a two dimensional plane of T * vs α * as 1 and 2 sigma contours, as shown in Figure 1, for the case of DWs as the only source contributing to the GW signal.
In order to compare the NANOGrav contours with the scenario we are discussing, we can first use equations (2.8) and (2.9) to relate directly the fraction of energy density to the annihilation temperature, for a given bias potential, Plugging in ∆V ∼ 2 m 2 π f 2 π in the equation above we obtain a line in the T * vs α * plane.Notice that each point on this line corresponds to a specific domain wall tension.
We display in Figure 1 the lines of two representative cases for , showing how the QCD-induced bias can accommodate the NANOGrav data.Notice that since ≤ 1 in this minimal realization, we cannot access the entire region favoured by the NANOGrav analysis in this model.
The results from the NANOGrav collaboration are obtained in a model-independent way assuming that the DW network is in the scaling regime at annihilation.In the following section, we shall investigate whether this assumption is compatible with the natural QCD bias.
4 The impact of friction from QCD In this section we study the friction acting on the domain wall as a consequence of the reflection of particles in the plasma.In particular, we are interested in friction effects close to the annihilation temperature of the network in the range relevant for Fig. 1.If friction dominates before annihilation, the DW network is not in scaling and the predictions for the GW spectrum can change significantly [68,99], possibly jeopardizing the PTA interpretation.The irreducible friction on ALP domain walls in scenarios which the QCD bias comes from gluons and in general from hadrons at low temperatures.
The effect of friction is usually parameterized by defining a friction length f which feeds in the total damping scale of the network d as where H is the Hubble parameter and Here v w is the average velocity of the network, and ∆P is the pressure on the wall from interactions with the plasma.The definition above takes into account that one generally expects ∆P ∝ v w , at least for moderate velocities.The pressure can be computed from an integral over the particle thermal distribution involving the reflection coefficient [68,134] where a = ±1 for FD or BE statistics respectively, g counts the number of d.o.f, and f (v) is the thermal distribution in the wall rest frame The reflection coefficient R(p z ) can be computed by solving the quantum mechanical reflection for a particle scattering off the ALP wall.For temperatures such that 1/ f 3H the DW network deviates from scaling and enters a friction dominated regime, where one expects suppressed GW signals [68,99].
The temperatures of interest are again in the ballpark of the QCD crossover.We face the problem following a very simple strategy 5 and we study two regimes at temperatures above and below the QCD scale: • For T 2 GeV we compute friction by considering the scattering of gluons off the ALP DW through the interaction (3.1).
• For T 60 MeV we employ chiral perturbation theory and the main source of friction is induced by scattering of pions off the ALP DW.
Even if we cannot compute the friction in the intermediate temperature regime, we will be able to draw interesting conclusions concerning the DW dynamics around the QCD crossover.

Friction from gluons
At high temperatures, the contribution to ∆P comes from gluons reflecting off the ALP domain wall.We stress that this effect is unavoidable in the scenario in which QCD provides the bias collapsing the network.
In the simplified picture of friction as coming from one-to-one particle reflection off the ALP wall, we can neglect as a first approximation the gluon self interactions and work at the linear order in the field fluctuations.In this limit, one recovers independent abelian equations of motion for each gluon degree of freedom.Additionally, when the ALP wall comes from the simplest cosine potential, a reasonable approximation for the reflection coefficient can be obtained analytically [137] (see also [100]).The reflection probability for a negative helicity gluon is given by where ρ = p z /m a , with p z the gluon momentum in the direction orthogonal to the wall, and The reflection for positive helicity gives a quantitatively very similar result as in (4.5) in the case of interest.At momenta much below the inverse width of the domain wall ∼ m a (which also sets the height of the potential/well seen by the gluon depending on the helicity) particles have a finite probability of being reflected ∝ β 2 .The pressure induced by this type of interaction may be computed according to (4.3) and grows with the temperature as where g = 2 • 8 for gluons.For temperatures T m a an even larger fraction of particles is simply transmitted and the ∝ T 4 behavior is tamed to a much slower increase ∝ T , ∆P = 2 Our calculation of the friction at low temperatures is of course limited by QCD becoming non-perturbative.In our analysis we will push this description down to T 2 GeV, keeping in mind that corrections should be expected at the low end of this region.Notice also that there is a source of model dependence given by the ratio of the QCD and dark QCD anomaly, N c /N d .Clearly, when N c → 0 the ALP is decoupled from QCD and indeed (4.7) and (4.8) yield a vanishing contribution.In the following we shall take N c /N d = O(1), keeping in mind that decoupling the ALP from QCD is anyway not compatible with the generation of the required bias term.
Let us also notice that the pressure from gluon reflection is generally much bigger than the high-temperature bias induced by the QCD instantons in (3.5).This crucially implies that the network can consistently reach a friction-dominated regime well before annihilation begins.

Friction from pions
In order to evaluate the friction from pions we refer to chiral perturbation theory and consider the potential in (3.2).Let us first stress that in the scenario of interest where the QCD contribution is not aligned with the potential induced by the dark QCD, the pion mass will change in the different vacua.This simply signals that the degeneracy has been indeed removed.Taking this into account, pressure from the pions is expected to be of the same order of the potential bias.
However, this pressure is not what we are interested in, as this would only determine the terminal velocity during the domain wall network collapse.Instead, the question we wish to address is whether the ALP interaction with pions could in principle turn a scaling network (where the bias is by definition irrelevant) into a friction-dominated evolution.
To determine this, we shall evaluate the pion pressure in a system in which the QCD potential and the dark-QCD potential are aligned.The pressure comes then from pion reflection off the domain wall, with the pion mass being the same on both sides.This has to be interpreted as a lower bound on the pion pressure.
We then set f a = f a only for this specific calculation, and we additionally assume a large hierarchy between the ALP and the pion mass, m a m π .This allows us to neglect the backreaction of the pion on the ALP domain wall solution following its own dark-QCD potential.We then set the ALP to its z-dependent background given in (2.5), and study the motion of the π 0 around it.
First of all, we have to take into account that the ALP background induces a zdependent background on the π 0 as well, given as the solution of The structure of the vacua implies that the profile for π b is such that π b (−∞) = 0 (where a(−∞) = 0), and π b (+∞)/f π = π (where a(−∞)/2f a = π).
In the assumed hierarchy m a m π , the ALP background varies on scales much shorter than the inverse pion mass.Thus, the equation above can be split on the two sides of the ALP wall as and where we have used the asymptotics for a(z).Once the solution for π b (z − ) is obtained for z < 0, the solution for z > 0 is given by As we can see, the equations of motion for π b are similar to the ones for the sine-Gordon model, where the potential would instead be ∝ cos(φ) for both positive and negative z.
Our case is however distinct to the sine-Gordon as both the equation of motion and the boundary conditions for π b are different.Contrary to the sine-Gordon model where particle excitations are exactly (self) reflectionless, we then generically expect a non-zero reflection coefficient6 .The shape of π b can be obtained by solving numerically the equation of motion.To this end we employ a relaxation algorithm to the extreme approximation of the ALP as a step function as in (4.10) and (4.11), as well as to the actual ALP profile with m π /m a = 0.1.We find that as long as the ALP mass is hierarchically larger than m π , the solution for π b is practically independent of m a .Qualitatively, one has even though corrections are clearly visible in Fig. 2 (left panel) as (4.13) does not in fact solve the equation of motion.Notice that the resulting ALP-π 0 domain wall has structure at two different scales, m −1 a and m −1 π , similarly to the case of η -π 0 domain walls in pure QCD [138].
Given our background solution, we can study the reflection probability for small oscillations around it.Writing π 0 = π b (z) + δπ 0 (x), and π 0 (x) = f (z)e iEt−ikxx−ikyy , we have The reflection coefficient is evaluated by solving (4.14) numerically.The result is again independent of m a as long as m a m π , and it is well approximated by where c 8 (see right panel of Fig. 2).For p z m π we are able to identify an exponential drop as expected when the momentum of the scattering particle is of the same order as the inverse wall width ∼ m −1 π .However, this kinematic region is irrelevant for our analysis, as we apply the pion Lagrangian only at temperatures T 60 MeV where high-momentum excitations are Boltzmann suppressed.
Using our result for the reflection coefficient and (4.3), we can straightforwardly evaluate the pressure from pion reflection in the alignment limit: where we have included g π = 3 expecting a similar contribution from the charged pions.

Implications of friction for PTAs
In this final section we summarize our results by indicating in the ALP parameter space where deviations from the scaling regime of the DW network are to be expected at temperatures around annihilation, possibly affecting the SGWB signal.
Notice that for temperatures 60 MeV < T < 2 GeV the pressure from the hadronic sector is not calculable within our simple approach.However, to extract information on this intermediate temperature range, we can look at the pressure at higher temperature (from gluons) and at lower temperature (from pions).Notice that for T < 60 MeV we consider the pressure from pions that would act on the domain walls as if they were to be still around (see the previous section for the details of this argument).In practice, the NANOGrav data suggests annihilation temperatures T * > 60 MeV for our model with the QCD bias, so that the would-be pion pressure is only useful for this extrapolation.For instance, if both the gluon and would-be pion pressure were to dominate in their temperature range of validity, we would conclude that annihilation in the intermediate temperature range is very likely to occur during friction domination.Otherwise, if friction dominates in only one of the two calculable regimes, a more detailed analysis around the QCD crossover is needed.
We now illustrate this strategy by presenting in Fig. 3 two benchmark points characterized by representative choices of the model parameters.The right panel shows a benchmark for which the signal from scaling domain walls and QCD-induced bias can explain the GW signal at PTAs (α * = α obs ).As we can see, the network does actually enter friction domination around T ∼ 100 GeV driven by the gluon scattering.Friction remains dominant also at T ∼ 2 GeV, which we take as the edge for the validity of the gluon calculation.However, at temperatures T ∼ 60 MeV the would-be pressure from the pions is insufficient to drive the DW network away from scaling.Therefore, it is possible that the period of friction domination ends in the gray region where neither of our calculation is applicable and the DW network goes back to scaling just before annihilation, which in this benchmark is predicted around T ann ∼ 124 MeV7 .Points of this type are shown in Fig. 4 in the purple region (labelled by gluon friction), where we suggest that a more refined analysis is needed in order to establish the viability to explain the PTA signal.
In the right panel of Fig. 3 we show instead a benchmark point for which the gluon and would-be pion are both satisfying the friction domination condition.In this case, it is very likely that the domain walls collapse without ever going back to the scaling regime, Figure 3. Benchmark points illustrating our friction-domination analysis.Left: Friction domination occurs for temperatures 100 GeV > T > 2 GeV driven by gluon pressure, as the properly normalized pressure (blue line) overcomes Hubble (green line).At lower temperature, the would-be pion pressure is unable to drive friction domination.Whether the network will have enough time to go back to scaling at T ann remains uncertain.If this is the case, this benchmark point is able to explain the PTA signal (α * = α obs ).Right: Gluon and would-be pion pressure are both capable of inducing friction domination, and thus it is very likely that the network never goes back to scaling above the annihilation temperature, identified as the crossing between the properly normalized bias (orange line) and Hubble.Points of this kind, however, require a relatively small domain wall tension and would not be able to explain the GWs observation even if they were to annihilate in the scaling regime (α * < α obs ).The anomaly coefficients have been chosen as N c /N d = 1.5.
with strong implications for the GW signal.However, points of this kind where friction dominates in both our calculable regions require a relatively small tension, and therefore cannot explain the GW signal observed at PTAs even if the network were to annihilate in the scaling regime (emphasized in the right panel as α * < α obs ).Points of this kind are found in the pion and gluon friction region in Fig. 4.
Let us now comment on the overall results shown in Fig. 4 as a scan over the (m a , f a ) parameter space.Additionally to the regions mentioned above, we see that there exists some parameter space for m a < 3 GeV where only the would-be pion friction is able to induce fricition domination, while gluons do not.This is understood by noticing that the gluon reflection becomes more and more suppressed as m a is lowered, see e.g.Eq. (4.7).On the other hand, as long as m π m a the would-be pion pressure is independent of m a .This, combined with the fact that gluons need to face a faster Hubble expansion at higher temperatures leads to the only-pion region in Fig. 4. Notice also that our scan does not extend to points with m a < 1 GeV as the approximation m a m π used in Sec.4.2 would break down.
Together with the colored regions indicating the impact of friction, we also show in (dark) gray the parameter space where domain walls come to dominate the energy energy of the Universe before annihilation for the choice = 0.26 ( = 1) for the bias in Eq. (3.4).These points are excluded from our analysis.
The parameter space that can fit the NANOGrav data if the network annihilates in the scaling regime is shown by the light blue band for = 1 and by the narrower orange Figure 4. Scan over the (m a , f a ) parameter space summarizing the results of our analysis.The light blue (orange) band indicates the parameter space that is compatible with the NANOGrav data in the ALP model with a QCD bias considered here with = 1 ( = 0.26).As the observed GW background is rather large, both our signal bands are not too far from domain wall domination, shown in the upper right corner by the dark gray (gray) region for = 1 ( = 0.26).The other colored regions highlight the relevance of friction.The purple region corresponds to the parameter space where gluon friction dominates over Hubble at T = 2 GeV, where we take α s = 0.2 and N c /N d = 1.5.This is the lowest temperature where the gluon computation can be trusted, see also Fig. 3. On the other hand, the would-be pion pressure is evaluated at T = 60 MeV and provides information about the friction in the confined phase, see text for details.The region where the would-be pion pressure can induce friction domination is shown in yellow, and its intersection with the gluon friction region is shown by the pink color.The implication for the ALP domain wall interpretation of the PTA data is as follows: for relatively light ALPs with m a < 10 GeV it is fair to assume that the network annihilates in the scaling regime, so that the signal bands shown here can indeed explain the NANOGrav data.On the other hand, for m a > 10 GeV friction is shown to be important at least to the right of the QCD crossover, and a more detailed analysis is required to assess the viability of this interpretation.band for = 0.26.These signal bands follow straightforwardly from the results shown in Fig. 1.As we can see, both these regions are not too far from domain wall domination, as expected given that the preferred values for the network energy density at annihilation are rather large, α ∼ 0.1.
The intersection of these signal bands and our friction regions provides the main result of our analysis, which we now summarize.Most of the parameter space compatible with the NANOGrav data implies friction domination from gluons at temperatures T > 2 GeV.However, the would-be pion pressure at low temperatures is not big enough to conclude that the network will be friction dominated at annihilation as well.We nevertheless suggest that a more detailed analysis is needed to ensure viability of these points.On the other hand, for a relatively light ALP with m a < 10 GeV we find no evidence for friction domination around the QCD crossover, and thus these points can be viable candidates to explain the PTA data.Even though our analysis cannot extend for m a < 1 GeV, we expect this conclusion to apply also for lighter ALPs.
Before concluding this section, let us mention again that our results only take into account the inevitable friction on the DW network in scenarios with a QCD bias, and that additional, although model-dependent, interactions with the other SM particles can provide important sources of friction as well (see [68]).

Conclusions
The results from PTAs have opened a new era of exploration of the Universe by providing the first evidence of a stochastic background of gravitational waves.One possible cosmological explanation for this signal is a network of DWs annihilating around the temperature of the QCD crossover.An important class of DWs predicted in BSM theories is the one arising in ALP models.These axionic DWs can accommodate the PTA signal if their annihilation is determined by the potential for the ALP that is dynamically induced by QCD as a consequence of the ALP-gluon coupling.However it is important to remark that the prediction of SGWB generated by the DWs is based on numerical simulations that neglect the interaction of the DW with the cosmic plasma, so that the network reaches the scaling regime.
Our main observation is that in scenarios where the DW annihilation is induced by QCD effects, there is an unavoidable source of friction exerted by QCD states scattering off the DWs.Our results are summarized in Fig 4 .We have identified the portion of the ALP parameter space where friction can be important, even though for domain wall tensions capable of explaining the NANOGrav data we cannot firmly conclude whether friction will be dominant at the annihilation temperature.This is because of lack of calculability around the QCD crossover in our simplified approach, and a more refined analysis would be then required.On the other hand, we were able to identify the region of ALP parameter space, namely m a 10 GeV, where friction is negligible and the ALP DW interpretation of the the NANOGrav signal is unaffected.

Figure 1 .
Figure 1.One and two sigma contours for the DW interpretation of the signal as provided by the [1] collaboration (blue and yellow dots).The prediction of a DW network with QCD induced bias (∆V ∼ 2 m 2 π f 2 π ) are displayed as lines with varying DW tension.

Figure 2 .
Figure 2. Left: Results for the pion profile π b (z) (blue) in the presence of a realistic ALP background (orange) with m π /m a = 0.1 obtained numerically with a relaxation algorithm.The dashed red line shows the qualitative behavior discussed in the text.The pion profile varies on a scale ∼ m −1 π much larger than the ALP background ∼ m −1 a .Right Reflection coefficient evaluated numerically (blue) and the approximation (4.15) with c 8 (orange).