Looking at the NANOGrav Signal Through the Anthropic Window of Axion-Like Particles

We explore the inflationary dynamics leading to formation of closed domain walls in course of evolution of an axion like particle (ALP) field whose Peccei-Quinn-like phase transition occurred well before inflationary epoch. Evolving after inflation, the domain walls may leave their imprint in stochastic gravitational waves background, in the frequency range accessible for the pulsar timing array measurements. We derive the characteristic strain power spectrum produced by the distribution of the closed domain walls and relate it with the recently reported NANOGrav signal excess. We found that the slope of the frequency dependence of the characteristic strain spectrum generated by the domain walls is very well centered inside the range of the slopes in the signal reported by the NANOGrav. Analyzing the inflationary dynamics of the ALP field, in consistency with the isocurvature constraint, we revealed those combinations of the parameters where the signal from the inflationary induced ALPs domain walls could saturate the amplitude of the NANOGrav excess. The evolution of big enough closed domain walls may incur in formation of wormholes with the walls escaping into baby universes. We studied the conditions, when closed walls escaped into baby universes could leave a detectable imprint in the stochastic gravitational waves background.

The solution [1] based on the assumption of the existence of a global U(1) PQ Peccei-Quinn (PQ) symmetry that is broken, at scale F , to a discrete subgroup Z N , which generates a pseudo-Nambu-Goldstone (PNG) boson, known as axion [6]. Due to the QCD anomaly the axion obtains a periodic potential with N distinct minima and hence acquire a mass specifically connected with the PQ scale as ∝ Λ 2 QCD /F . Such axion has been baptized the QCD axion [6][7][8][9]. It is suitable to set the PQ symmetry breaking scale much above the electroweak scale, so that the axion should have a very reduced coupling to the Standard Model particles, which makes the axion "invisible" rendering it a very promising candidate to act as dark matter (DM) [6][7][8][9][10].
The axion may naturally realized in string theory [11,16,17] where, moreover, the compactification always generates PQ-like symmetries which generically produce the multitude of axions and ALPs. The number and properties of the ALPs are defined by the specific model being considered. However, typically those string models which incorporate the QCD axion produce ALPs with PQ-like scale of about or higher than the GUT scale and masses similar to that of the QCD axion [11].
In analogy to the QCD axion, a high PQ-like symmetry breaking scale assumes that this symmetry is broken before inflation. Thus, during inflation, the ALP behaves as a massless spectator field, which is made uniform trough the scale of the observable universe. The value of the ALP field, over this biggest scale, is defined by some random displacement θ 0 from the minimum of its potential, which is called the misalignment angle. The misalignment angle θ 0 determines the amplitude of the ALP field, once its mass exceeds the Hubble rate and the friction term does non prevent it anymore from oscillations about its minima. These coherent oscillations of the Bose-Einstein condensate are pressureless, representing a dust like substance acting as a cold DM (CDM) in the Universe [3][4][5].
Either the QCD axion [7] or ALPs [7,12,13] with above GUT scale of the PQ symmetry breaking might be accomodated only if we inflated from a rare patch of the space with very small misalignment angle θ 0 << 1. The unlikely value can be explained by involving of anthropic selection based on the arguments that a high axion or ALP density should be hostile to the development of life [18][19][20][21][22]. In this case, it is said that the axion (or ALP) parameters are located in the anthropic window.
A massless inflationary ALP spectator acquires the quantum fluctuations imprinted by the de Sitter expansion. When the ALP obtains its mass, well after the end of inflation, these fluctuations, being converted into isocurvature fluctuations, become dynamically relevant.
The isocurvature fluctuations are uncorrelated with the adiabatic fluctuations inherited by all other matter and radiation from the fluctuations of the inflaton field [21][22][23][24][25][26][27]. The main observational effect of the ALP isocurvature fluctuations consists in modification of the Sachs-Wolfe plateau at small cosmic microwave background (CMB) multipole orders corresponding to the scales comparable to the size of the last scattering surface. So far, there is no trace of the isocurvature fluctuations in existing CMB data [28]. Therefore, it would be interesting to study opportunities, if any, that ALPs from anthropic window could show up in other types of observations.
The inflationary quantum fluctuations [29,30] of the ALP field should lead to substantial deviations from the initial misalignment angle [31], at the late stages of the inflation, when scales much smaller than that ones relevant for CMB observations were exiting the inflating Hubble patches. The deviations can be so large that, in some domains, the ALP field is already localized in a position to start oscillations about its different minimum, once it is liberated from the Hubble friction. In this case the neighboring domain with ALP oscillating about different minima must be interpolated by a closed domain wall, with stress energy density defined by the scale of PQ-like symmetry breaking and the mass of the ALP [32,33]. The abundance of such walls of a given size is defined by the initial misalignment angle θ 0 and the ratio of the PQ-like scale to the inflationary Hubble rate. Since, here, we are talking about scales much smaller than those ones relevant for the CMB observations, the inhomogeneities, caused by dramatic difference in values of the initial amplitudes of the ALP oscillations, cannot be observed as the isocurvature fluctuations, discussed above. Depending on its initial size, the evolution of a closed domain wall in the FRW Universe consists in either its collapse into a primordial black hole (PBH) [34][35][36][37][38][39][40][41] or formation of a wormhole with the wall escaping into a baby universe. In both cases of wall's evolution, a certain mass is set in a motion, which may lead to a production of gravitational waves (GWs). Indeed, a collapsing domain should release its asphericity in course of its entering under the Hubble radius, so that the biggest wall's fragments make one or few oscillations radiating their energy in form of GWs with characteristic frequency of about the Hubble rate. A closed wall forming a wormhole, due to its negative pressure, may induce sound waves in the interior bulk fluid leading also to emission of GWs in similar frequency band.
data of 45 pulsars timing array (PTA) measurements [42]. The signal can be interpreted as a stochastic GWs background 1 at frequency around 31 nHz. In the mean time, it is still not clear whether the detected signal is of truly origin from the stochastic GWs background because of the tension with previous PTA measurements in similar frequency range and also there has been not detected the quadrupole correlations 2 as yet.
Nevertheless, it was quickly realized, that apart the common astrophysical sources, such as super massive black holes binaries (SMBHBs), the stochastic GWs interpreted signal of the NANOGrav [42], might be an aftersound of the vacuum rebuilding processes which could take place in the early Universe. These are the formation of cosmic strings [46][47][48][49][50], first order phase transitions [51,52], the formation of domain walls [53,54], turbulent motions occurring at QCD phase transitions [55,56] and cosmological inflation [57].
In this paper we estimate the spectra of GWs produced by closed domain walls, the formation of which is induced in course of inlationary dynamics of ALPs with parameters, particularly but not exclusively, related to the anthropic window. We found that the walls collapsing into PBHs, at the stage of occurring of sphericity, can generate the stochastic GWs background with the characteristic strain power spectrum slope remarkably well centered within the range of slopes reported in the NANOGrav signal. Analyzing the inflationary dynamics of the ALP field, in consistency with the isocurvature constraint, we reviled such combinations of the parameters where the signal from the inflationary induced ALPs walls could saturate the amplitude of the NANOGrav excess. As a by product of the analyzes, we have elucidated the conditions, when closed domain walls escaped into baby universes could leave a detectable imprint in the stochastic gravitational waves background.
The rest of the paper is organized as follows. In Section 2, we describe the mechanism of origin of closed domain walls in presence of ALP field during inflation and their subsequent evolution. In Section 3 and Section 4 we derive the spectra of stochastic GWs background generated by collapsing and escaping domain walls, respectively. In Section 5, the derived spectra are related to the spectral slope and the amplitude of the signal reported by the NANOGrav. In Section 6, we analyze the parameters of ALP in the light its inflationary dynamics and saturation of the NANOGrav signal amplitude. In Section 7, we define the combinations of parameters in scenarios where ALPs could manifest themselves in the NANOGrav signal, being compatible with DM density and isocurvature constraints. In Section 8, we discuss possible signatures of domain walls escaping into baby universes in stochastic GWs background signals.
In Section 9, we check the consistency of the PBHs production by closed domain walls versus existing PBHs constraints. We conclude in Section 10.
1 PTA measurements for detection of very long gravitational waves have been proposed in [43,44]. 2 The quadrupol correlations is a smoking gun of stochastic GWs background in PTA measurements [45].

Inducing domain walls formation in presence of ALP while inflating
We consider potential of a complex scalar field φ = r φ e iθ with self-interaction constant λ whose U(1) PQ symmetry is broken at the scale obeying the condition which implies that the radial mass m r = √ λF exceeds the inflationary Hubble rate H inf , so that the field φ resides in the ground state during inflation. The ALP is then the pseudo Nambu-Goldtone (PNG) boson of the broken U(1) PQ symmetry Unlike in case of the QCD axion, the mass of ALP field a should not have specific connection to the scale of the U(1) PQ symmetry breaking F . In case of realization of axion scenario in string theory, compactifications always generate PQ like symmetries [11,16,17] and generically produce a multitude of ALPs.
We assume that ALP potential, being generated by some exotic, strongly interacting, sector, may be approximated as where Λ is a scale set by the ALP coupling to instantons and may span a wide range from QCD scale to the string scale, with each ALP associated with its own gauge group [16]. The ALP's mass reads so that the ALP is described by two out of the three parameters F , m θ and Λ. If string theory can produce a QCD axion with large scale of PQ symmetry breaking and hence with small mass, one can expect that the spectrum of ALPs should include masses within few orders of magnitude of that one of the QCD axion. String theory models, as for example are described in [11,16,17], can accommodate many ALPs with F values about the GUT scale.
The masses of such ALPs should be homogeneously distributed on a log scale [16], with several ALPs per energy decade.
During inflation and some period of Friedman-Robertson-Walker (FRW) epoch, the energy density (4) is negligible starting playing a significant role in the dynamics of phase θ(x) at the moment when the mass of the ALP (5) overcomes the Hubble rate. The potential (4) possess a discrete set of degenerate minima corresponding to the phase values θ min = 0, ±2π, ±4π . . . .
After the inflation ended, as soon as the m θ started to exceed the Hubble rate, the field θ(x) becomes oscillating about the minima, so that the energy stored in the potential (4) gets converted via misalignment angle mechanism [3][4][5] into the ALPs in form of Bose-Einstein condensate behaving as a cosmological cold DM (CDM) 3 (see discussion in Section 7).
At the condition of vanishing m θ , the inflationary dynamics of the phase is driven by the quantum fluctuations of magnitude [18,24,25,31] δθ ≈ 1 2π taking place each Hubble time H −1 inf . The amplitude of these fluctuation freezes out due to a large friction term in the equation of motion of the massless PNG spectating the de Sitter background, whereas their wavelength grows exponentially. Such process resembles one dimensional Brownian motion of variable θ, inducing, at each inflationary e-fold, a classical increment/decrement of the phase by factor δθ given in (6). Therefore, once the inflation began in a single, causally connected domain of the horizon size H −1 inf , being filled with initial phase value θ 0 < π, it progressed producing exponentially growing domains with phase values being more and more separated from the initial phase value θ 0 . It is unavoidable, that in some domains the phase values has grown above π, so that one expects to see a Swiss cheese like picture where the domains with phase value θ > π are inserted into a space remaining filled in with θ < π. Therefore, in the domains where the inflationary dynamics has induced the phase value θ > π the oscillations of the ALP field will occur about the minimum in θ min = 2π, while in the surrounding space with θ < π the oscillations will chose the minimum in θ min = 0. It is well known [58], that in such a setup, in case of periodic potential (4), two above vaccua should be interpolated by the sine-Gordon kink solution a(z) = 4F tan −1 exp(m θ z) interpreted as a domain wall of width m −1 θ posed perpendicularly to z axis, with the stress Therefore, one expects that every boundary separating a domain with phase value θ > π from the space of θ < π traces the surface of a closed domain wall to be formed in some time after the end of inflation, when θ(x) starts to oscillate. Since after inflation, during FRW epoch, the wall of size R(t end ) is simply conformally stretched by the expansion of the universe where a(t) is the scale factor and t end corresponds to the end of the inflationary epoch, the size distribution of such a sort of inflationary induced walls' seeding contours can be directly mapped into the size distribution of the domain walls. In what follows, we derive this size distribution.
A domain wall's seeding contour of radius R ≈ H −1 inf emerging during inflation of total duration t inf at the time moment t s , when the universe is still have to inflate over ∆N s = H inf (t inf − t s ) = N inf − N s e-folds, is getting stretched in course of the expansion as The number of contours created in a comoving volume dV within e-fold interval dN s is given by where Γ s is the contours' formation rate per Hubble time-space volume H −4 inf , which, in general, should depend on formation instant N s and is defined by the inflationary dynamics of the ALP spectator field (see the detailed discussion in Section 6). Expressing N s from (10), one can write down the number distribution of seeding contours with respect to their physical radius R dN = Γ s e 3Ns dV R 4 dR.
Thus, the number density in physical inflationary volume dV inf = e 3Ns dV reads At the end of inflation the distribution (13) spans the range of scales from R min where N π is the number of e-folds when the angular degree of freedom θ overcame π, since the beginning of inflation lasting enough N inf e-folds needed to solve horizon and flatness problem 4 .
Once the phase oscillations begin and a domain wall is materialized, replacing the surface of its seeding contour, it will stay essentially at rest relative to the Hubble flow until its size becomes comparable with the Hubble radius R(t H ) ≈ H −1 defined by the Hubble rate H =ȧ/a. Evidently, at the entering into the horizon the surface tension forces developed in the wall tend to minimize its surface. The dynamics of spherical domain wall has been studied in [64] in asymptotically flat space. It was shown, that in this case the internal domain wall metric is Minkowski while the external one is Schwartzshild, so that the wall being initially bigger than its Schwartzshild radius always collapses into a BH singularity. Thus, the ALP field induced domain wall should obtain a spherical geometry and contract toward the center. The regime (2) assumes that the ALP's coupling to the Standard Model particles is strongly suppressed by the large scale F , so that the ALP wall interacts very weakly with matter and hence nothing can prevent the wall finally to be localized within its gravitational radius and hereby deposit its energy into a PBH [32,33].
A more general case of a false vacuum bubble surrounded by true vacuum, and the study the motion of the domain wall at the boundary of these two regions are presented in [65][66][67]. In particularly, it was shown that if the energy density inside the bubble is vanished and the unique source of gravitation field is the domain wall, as for example occurs when the domain wall emerges from a "white hole" singularity, it will expand unbounded to escape in a baby universe. The baby universe is initially connected to the asymptotically flat patch of space by a wormhole, which eventually is getting pinched off leading to changes in topology.
A similar picture of evolution may take place in case of a wall exceeding by its mass the energy of fluid inside the Hubble radius, at the entering. Indeed, as it is described in [61,62], due to its repealing nature, the wall, with dominating gravitational effects, pushes away the radiation/matter fluid, creating rarefied layers in the vicinity of its both interior and exterior surfaces. Since the exterior FRW Universe continues to follow the Hubble expansion, while the wall extends exponentially in proper time, it has to create a wormhole, through which it escapes into a baby universe. Finally, the wormhole is getting pinched off, within about its light crossing time, so that observers on either side of the throat of the wormhole is seeing a BH forming, probably having unequal masses on both sides.
Here we notice that, in general, along with the phase fluctuations (6), one has to expect the radial fluctuations which might kick the radial component r φ out of its vacuum state r φ = F , over the top of the Mexican hat potential (1). This would lead to the formation of strings along the lines in space where r φ = 0, so that one arrived to a realization of the scenario of the system of walls bounded by strings [58,59]. The radial fluctuations, in the context of string formation during inflation have been investigated both, analytically and numerically in [31]. In particular, it was elucidated that the strings would be formed if the Hubble rate were exceeding the radius of the Mexican hat potential. More precisely, if the condition (2) is violated, during inflation 5 . In some sense, there is a "Hawking like" temperature of order of H inf which takes over to drive a thermal like phase transition with ALP string formation. In contrary, r φ sits in its vacuum r φ = F whenever (2) stays in tact, for a reasonable value of λ 10 −2 and simple inflationary model [31]. Therefore, in the regime Q FH inf > 1 relevant for the analysis presented bellow, only closed domain walls can be formed, since the formation of strings is prevented.
Actually, the closed walls can decay due to quantum nucleation of holes bounded by strings on their surfaces [58,60]. The process of the hole nucleation can be described semi-classically by an instanton so that the decay probability is expressed as [60] where µ πf 2 is the string tension and A is a non infinitesimally small factor which can be calculated from the analysis of small perturbations around the instanton solution. Thus, applying (8), one can see that the probability (14) is suppressed by factor exp( which is very small for a large hierarchy of the scales, F/Λ >> 1, which is the case for the parameters of ALPs, discussed bellow. Thus, one can affirm that the closed ALP's domain walls considered here are stable relative to the quantum nucleation of holes bounded by strings.
The process of occurring of the spherical shape of a collapsing closed domain wall should be accompanied by a production of gravitational waves [68][69][70]. Also, one might assume that sound waves induced in the interior radiation fluid by walls escaping into wormholes may also contribute as a source of stochastic gravitational waves background.
Bellow, we study the domain walls formed as a result of the dynamics of ALP field found itself in the regime of the inflationary spectator, as a source of primordial stochastic GWs background, in both collapsing and expanding regimes outlined above. We derive the spectra of the GWs generated by domain walls of sizes distributed with respect to (13). The spectra, being related to the NANOGrav signal [42], elucidate the detectability of the evidence of existing of ALPs, including those characterized by parameters relevant for the anthropic window.

Stochastic gravitational waves background generated by the collapsing ALP's walls
The ceeding contours of closed ALP's domain walls, formed in the way described above, have complex structure which is characterized by inflections (or folds) in all possible scales including those comparable with the size of the walls as whole objects (see a detailed discussion in [41]).
After the materialization, as soon as a folded fragment of a closed domain wall becomes causally connected, namely it finds itself within a Hubble horizon comparable to the scale of the inflection's curvature, this fragment is getting stretched due to the wall's tension forces, so that its shape is smoothed out within the Horizon scale. Therefore, once a close domain wall finds itself as whole causally connected unit, being localized within a Hubble horizon which encompasses it as complete object, the wall will mostly contains fragments inflected (folded) in scales which are as big as the scale of the horizon. The stretching of these biggest fragments leads to oscillations in scales comparable with the size of the whole closed wall. Since the ALP's walls, for the PQ-like scale of symmetry breaking considered here, are coupled extremely weekly with the SM particles media the oscillations can be dumped releasing their energy in form of GWs contributing into the stochastic GWs background, which might be observed today. Finally the closed domain wall obtains a spherical shape in course of its one or, may be, few oscillations within the scale of the Hubble radius 6 .
Here, we are interested in the presently existing stochastic background of GWs, created by domain walls described above, as the dimensionless fraction of the critical density expressed by the energy of GWs in units of logarithmic interval of frequency where ρ gw is the energy density of GWs in unit frequency. The energy density in comoving volume is just the total energy deposited in GWs where P gw (t, f ) is the total GWs energy power emitted at time t into unit range of frequencies.
The frequency f today had been emitted as frequency f = (1 + z)f , so that we can write This background can be expressed as its power spectral density [71] S The pulsar timing arrays like NANOGrav [42] use so called characteristic strain The power of the emission of the GWs by an individual domain wall can be estimated using the formula of quadrupole radiation power by a massive object [72] where I is the quadrupole moment of the object and C 2 = 45. We assume that just in the course of a complete entering of a closed wall of size R into the matchable Hubble radius, its biggest fragment(s) pass through one or few oscillations within the horizon scale. For such a fragment of a domain wall the quadrupole moment is estimated as where f ≈ R −1 is the frequency of the GW emitted within a Hubble horizon of size R.
Therefore, in this case, the power (20) can be expressed as The GWs power P gw (t, f ) is contributed by domain walls that collapse having radii between R and R + dR at Hubble crossing. For such domain walls, collapsing before mater-radiation equality t H < t eq , the Hubble crossing is t H = R/2, so that their number density evolves as where t > t H . Thus, one can write Substituting (24) in (17) we can change the integration variable using to get For the ΛCDM cosmology we have where the function G(z) accounts for changes of number of relativistic degrees of freedom at early time, H 0 = 100h km/s/Mpc. As the standard parameters we will use Ω Λ = 0.69, Ω m = 0.31, We argue in Section 6 that the scales, relevant for the dynamics of domain walls in the current study, emerge well before equality epoch. Thus, one can integrate (26) in the radiation dominated area ignoring the change of the number of degrees of freedom. The ratio of the current scale factor a 0 to that one corresponded to the equality of matter and radiation a eq occurred at red shift z eq is given by In the radiation dominated regime one can write so that where the contribution of radiation to the Hubble rate is given by Therefore, the integral (26) can be reduced to Finally, one can express the fraction (15) of the stochastic GWs background generated by the collapsing domain walls as follows 4 Gravitational wave signal induced by walls escaping into baby universes Due to its repulsive nature the domain wall, which dominates by its mass over the energy of the interior fluid at the Hubble radius entering instant, pushes away from its interior surface the radiation bulk fluid and then inflates out forming a wormhole [61,62]. We assume that the radiation fluid should respond on this act of repulsion by sound waves of characteristic wavelength H −1 which set in a motion a fraction of the mass contained inside the Hubble radius. This motion should generate GWs of frequency f ≈ H, which will give a contribution into stochastic GWs background 7 .
In such a setup the quadrupole moment in (20) can be estimated as where M b (H) is the mass of the radiation fluid that would be contained inside a dominating wall at the instant of its entering into the Hubble radius and κ is the fraction of this mass which 7 A similar kind of contribution is usually considered in connection with GW signal generated by a first order phase transitions, see for example [73].
is set in a motion by the sound waves. Therefore, the power (20) generated by an individual wall can be expressed as In case of size independent κ, which is used bellow, the power (35) is frequency independent as well. For the domain walls, escaping before mater-radiation equality t H < t eq , and hence pushing away the interior radiation fluid inside the Hubble radius t H = R/2, the differential power can be estimated in the analogy to (24), and reads Further, proceeding in the way it is done to derive (32), we arrive to Finally, one can express the fraction (15) of the stochastic GWs background left over by the domain walls escaping into wormholes as follows

Connecting to the NANOGrav signal
The NANOGrav [42] is sensitive to the characteristic strain (19) of GWs background presented in terms of power low spectrum where f yr = 1yr −1 = 31 nHz, A yr B is the amplitude at f yr . The spectrum (39) is obtained in direct measurements of the timing-residual cross-power density, whose slope is parametrized as [42]. The fit (39) is performed to thirty bins within frequency range from 2.5 nHz to 90 nHz. However, the excess is reported in first five signal dominated bins spanning the range from 2.5 nHz to 12 nHz, while the higher frequency bins are assumed to be white noise dominated. The signal excess of the strain (39) is reported for the parameters range and 4.5 ≤ γ ≤ 6.5 at 68% confidence level (C.L.).
Using the formulas (18), (19), (28) and (33) we calculate the signal strain of GWs generated by the ALP field induced collapsing domain walls where Therefore, comparing with (39) we find that slope parameter of the spectrum of stochastic GWs background generated in course of the evolution of collapsing domain walls, induced by ALP inflationary dynamics, corresponds to γ = 5.5, which exhibits a remarkable agreement with the range of values (41) reported in NANOGrav signal [42]. The stress energy density (8) gives rise the estimate so that (40) implies 8 where we put A yr B = C A × 10 −15 so that C A = 0.16 ÷ 1. Thus, the axion mass value needed to saturates the amplitude (40) of the NANOGrav signal lies at the level of where M Pl is taken for the Plank mass. Formula (47) can be presented as where, to evaluate the numerical pre-factor, we applied the upper bound as inferred from Plank results [28], while the ratio Q FH inf is still kept as a variable which controls Γ s during the inflationary dynamics of the ALP field, as we will see in Section 6.
Using (18), (19), (28) along with (38) one can also estimate the characteristic strain of the signal generated by the acoustic response of bulk radiation fluid on the formation of wormholes by domain wall with dominated gravitational dynamics Therefore, with respect to (39), the slope parameter of the acoustic response spectrum corresponds to γ = 3.5, which is more than 1σ off from the central value of the range (41). Obviously, the total spectrum, around f yr , can be comparably contributed by both h cw (f ) (44) and h cb (f ) (50), only in the case of valid relation which is not the case for the ALPs parameters localized in Section 7.
6 Inflationary dynamics of the ALP field so that, as follows from (29), The condition above simply demands that the ALP mass m θ should exceed the blue shifted characteristic scale of the NANOGrav, expressed as (1 + z osc )f yr . Therefore, the lower bound on the mass of the ALP reads We assume that the inflation begins in a single Hubble volume containing a phase value θ 0 < π which gets a random kick of magnitude δθ, given by (6), from vacuum quantum fluctuation of Fourier modes leaving the horizon. In other words, the process of the phase variation can be interpreted as a one dimension Brownian motion (random walk) of step length δθ (6) per Hubble time [18,29,31]. Let us consider, in such a setup, the inflationary dynamics of the phase to define the formation rate Γ s . For every value of initial separation ∆θ π = |π − θ 0 |, one can define the probability density p(∆θ π , N ) that a Brownian path of the fluctuating phase will cross π first time at a given e-fold N , in a given Hubble volume H −3 inf . In general formulation, the first passage probability density is calculated as one that a Brownian path, starting in x 0 at t = 0, would cross the origin for the first time at time t. This first-passage probability derived using random walk [74] is given by where D is interpreted as diffusion coefficient expressing the randomness.
In the language of the inflationary dynamics of the phase θ with respect to e-folds, we may use the simplest one dimensional Pearson walk, which implies D=1/2, and re-define the special degree of freedom in (55) as the phase difference measured in steps δθ, given by (6), so that Since the quantum fluctuations increment (decrement) the phase value by δθ, after each e-fold, one has to use the e-fold as the time variable, in (55), so that t := N . In this setup the expression (55) is converted into which will be used bellow to account for the probability of the domain walls' seeding contours formation in the number density distribution (13).
As the first-passage happened, after N e-fold from the beginning of inflation, the Hubble volume, where it took place, became to be filled with the phase value in the vicinity of π, so that the separation of the phase value from π should obey ∆θ first ≤ δθ. In this position of the phase, the probability density of π crossing, in a Hubble volume, at each e-fold, should be calculated from (57) as Therefore, the probability density Γ s , corresponding to the horizon exit scale e-fold N s , can be factorized as Γ s ≈ p(∆θ π , N s − 1)p(δθ, 1).
Obviously, that for all smaller scales, which correspond to later horizon exit e-folds N l > N s , Γ s can be taken as a constant, provided that p(∆θ π , N s − 1) << 1.
e-fold, after the beginning of inflation. Using t eq = 51.1 kyr [76] 9 one arrives to the estimate N yr ≈ 18.8.
Based on the above consideration, one can express the ALP's mass value (48), saturating the amplitude of NANOGrav signal, as a function of the initial phase separation ∆θ π and the scale of the PQ symmetry breaking Q FH inf measured in units of the inflationary Hubble rate The expression (61) is visualized in Fig.1 as the function of Q FH inf , for C A = 1 and three distinct values of the phase separation, namely ∆θ π = π/2 (left panel, dashed line), ∆θ π = π (left panel, solid line) and ∆θ π = 0.1 (right panel).

Which ALPs may show up in the NANOGrav signal
Let us look at the ALP's inflationary dynamics related to the NANOGrav signal, as described above, from the point of view of its simultaneous compatibility with constraints on the DM density and isocurvature fluctuations.
In the biggest scale, emerged at the beginning of the inflation, from a single causally connected domain of size H −1 inf , the ALP field is uniform and frozen at the initial phase value θ 0 , called misalignment angle. The ALP remains frozen until the equality (52) gets satisfied, which takes place when the temperature of the Universe drops bellow certain value T osc . In general, m θ should depend on the temperature, which, for example, in the case of the QCD axion, is defined by the relation between T osc and Λ at which the instanton effects become important (see the most advance lattice calculations in [75]). For ALPs, such kind of dependence should be driven by the details of its gauge couplings and cannot be obtained though a straightforward generalization of the axion case. Bellow, when we quantify the cosmological abundance of ALPs, we treat the mass as being independent of temperature. Our analysis is not sensitive to this assumption.
The total density of ALPs produced in the regime of the inflationary spectator ALP field, namely under the condition (2), is observationally constrained by the CDM density. In particular, in the regime (2), it is required that the mass of ALP and its misalignment angle, defined by θ 0 , would be mutually fine-tuned. Such kind of fine-tinning can be substantiated with invoking of the anthropic reasoning, which initially has been applied to the QCD axion [18,19].
The anthropic selection, applied to the QCD axion [18,19], implies that the misalignment angle should be small (θ 0 << 1) for PQ symmetry breaking scale relevant for the inflationary QCD axion spectator. In the version of the anthropic selection demanding that any selected variables should have values near the maximal consistent with existence of the human being, the CDM should completely consist of the QCD axions, so that the anthropic window of parameter space is well defined for the QCD axion [18,19]. However, in this anthropic window a large isocurvature mode is strongly preferred [21,22,25,26], which is at odds with observations, as we discuss bellow. The above constraints apply equaly to QCD axion and ALPs. In particular, a multitude of ALPs is expected to be produced in string theory framework attempting to realize the axion scenario [11,16,17], where all ALPs can contribute into CDM content and produce isocurvature perturbations. In this context, it is instructive to investigate, which ALPs from the multitude, possibly existing in string theory, could be seen as the source of stochastic GWs background being able to satisfy the results of ALPs related NANOGrav signal compilation presented in Section 5 and Section 6.
In Section 5 we established that the slope parameter of the spectrum of stochastic GWs generated by the domain walls, created in the ALP inflationary dynamics, exhibits a remarkably good agreement with that one reported by the NANOGrav. Further, we specify the values of the ALP's scale F , its mass m θ and the misalignment angle θ 0 which would lead to the saturation of the NANOGrav signal amplitude (40) simultaneously satisfying the anthropic and isocurvature constrains. In our study, we follow the QCD axion cosmology to describe the origins of the constraints and quote the relevant formulas. Except for the differences due to the treatment of the temperature-dependent mass, the formulas apply equally well for the QCD axion and for ALPs.
When the ALP field starts to oscillate, the ALPs are produced via vacuum misalignment mechanism [3][4][5] in form of cold Bose-Einstein condensate with local number density defined by the initial phase value θ 0 where f (θ 2 ) 10 is a correction for anharmonic effects of the axion-like potential. The number density n a red shifts in the same manner as the entropy density, so that where s = (2π 2 /45)g * s T 3 osc is the thermal entropy density when the oscillations begin, while s 0 = (2π 2 /45)g * s0 T 3 0 is the nowadays entropy density g * s0 = 3.91 and T 0 = 2.73K. The factor R γ indicates a possible entropy production after the axion begins to oscillate. Eventually, the above outlined misalignment ALPs abundance can be expressed as (see for example [14,15]) This implies, that in the regime when the PQ-like symmetry remains broken during the inflation and afterwards, ALPs in neV mass range will contribute 100% into CDM, provided that the at k = 0.05 Mpc −1 , as inferred in [28]. The details of estimation technique of the rate (65), for QCD axion and ALP, are given in [25]. In particular, the relative contribution (65) can be parametrized as follows where R ALP = Ω ALP /Ω CDM , is the ALPs' fraction of the total CDM content of the Universe, and R 2 (k) and S 2 ALP (k) are the adiabatic and the entropy power spectra, respectively. The curvature power spectrum R 2 (k) , for the adiabatic mode of fluctuations of inflaton reads where is the first inflationary slow-roll parameter [30] and the subscript indicates that the quantities are evaluated at k = aH inf . Provided, that the quantum fluctuations of the phase (6), are imprinted with a nearly scale invariant spectrum one obtains [25] Therefore, the isocurvature contribution (66) can be related to the fundamental inflationary and ALP's parameters by where it is assumed that α 1, as follows from (65).
Taking into account (64), one can express the isocurvature fraction (70) as follows During inflation, the primordial tensor perturbations are generated, setting the initial amplitude for GWs which oscillate after horizon entry. The spectrum of these tensor perturbations is conveniently specified as by the tensor fraction r = R 2 T (k) / R 2 (k) . In the slow-roll approximation [30] r = 16 (72) is experimentally bounded to r < 0.058 [76] (at k = 0.002 Mpc −1 ), which we use to estimate the upper limit on the first slow-roll parameter Evaluating (71), under the condition of saturation of the NANOGrav amplitude (40), which implies that m θ = m A θ (61), and for the upper bounds (49) and (73), one can frame three scenario in which the ALP might contribute, by its inflationary dynamical induced formation of the domain walls, into the stochastic GWs indicated in the NANOGrav observations.
(i) In the first scenario, the initial value of the phase, during inflation, might be chosen in the vicinity of the origin θ 0 ≈ 10 −2 (∆θ π ≈ π), which would correspond to the anthropic 10 0 Figure 2: The fraction of ALP-type isocurvature perturbations as the function of Q FH inf (solid line) evaluated under the condition that the ALP mass saturates the NANOGrav signal, as shown in Fig.1. The upper panel corresponds to the initial phase value θ 0 ≈ 10 −2 (δθ π ≈ π). The middle panel shows the functional dependence for θ 0 = π/2 (δθ π = π/2). The lower panel corresponds to the initial phase value in the vicinity of π, namely θ 0 = π − 0.1 (δθ π = 0.1). The horizontal, dashed line indicates the bound (65), deduced from the Plank measurements.
window usually discussed in the context of the QCD axion. The bound (65), as it is displayed in Fig.2 (upper panel), allows to accommodate ALPs with their inflationary dynamics driven by the ratio Q FH inf ≥ 6.1 (74) which corresponds 11 , at bound (49), to In this scenario, according to (61), as it is shown in left panel of Fig.1, ALP with detectable domain walls GWs signature should have mass and contribute, as it follows from (64), of the total CDM content of the Universe. Eventually, applying (8), one obtains the value of the stress energy density of the domain walls (ii) The second scenario assumes that the initial phase value θ 0 = π/2 (that is O(1)), and, as it is shown in Fig.2 (middle panel), accommodates the ALP with ratio which corresponds to F ≥ 7.0 × 10 14 GeV.
In this scenario, the ALP with detectable domain walls implies as it follows from (61) (see left panel of Fig.1). Thus, according to (64) which implies that the ALPs constitutes total CDM budget of the Universe. The value of the stress energy of the domain walls reads (iii) In the third scenarios, we put the initial phase value to θ 0 = π − 0.1 (∆θ π = 0.1), which implies the initial phase position θ 0 in a vicinity of the domain wall formation crossing point. Thus, the isocurvature bound (65), as shown in lower panel of Fig.2, can accommodate the ALP with ratio which would correspond to the ALP mass, saturating the NANOGrav signal amplitude (see Fig.1, right panel), However, according to (64), the ALPs with such mass and initial phase, in the vicinity of π, are definitely overabundant. If we increase the mass limit up to which implies that (see Fig.1, right panel) and the ALPs will constitute the total CDM content of the Universe (R ALP 1) and give a negligible contribution into the ALP-like isocurvature fraction, which reads In this case, one expects the highest possible value of the stress energy density of the domain walls, namely σ 5.7 × 10 18 GeV 3 .
Thus, in the context of detectability of ALPs from their multitude, possibly produced in string theory [11,16,17] along with the QCD axion, one infers the following. If one of the ALPs has the misalignment angle value set to the same one, θ 0 ≈ 10 −2 , imposed by the anthropic window of the QCD axion, it may saturate the NANOGrav excess amplitude even having a negligible contribution into the CDM density (77), at the same time being consistent with the isocurvature constraint. If one of the ALPs has large value of the misalignment angle, θ 0 ≈ 1, it should dominate the CDM density (82) and the isocurvature contributions to be interpreted as the source of the saturation of the NANOGrav amplitude. In case of one of the ALPs has θ 0 value very close to π, which is, however, as probable as the QCD axion anthropically selected θ 0 , it may manifest itself in the NANOGrav, in the condition of the dominant contribution into the CDM density and negligible contribution into the isocurvature perturbations (89).
However, in this case, instead of collapsing the domain walls of NANOGrav scale range should form wormholes, as we elucidate in Section 8.
Let us notice that, according to (113), at N ≈ 60, the domain wall formation crossing point is reached within almost 100% of the horizon exiting scales, for each of three above considered scenario. However, these scales, corresponding to the very end of inflation, are definitely much smaller than the width of domain walls typical for the scenario described above, so that the dynamics of the resulting vacuum-like objects should be much different than that one of the domain walls discussed above.

Domain walls escaping into baby universes
Once the wall becomes encompassed by a Hubble radius it collapses, within about a Hubble time, into a PBH of mass M PBH = ξM w , where ξ 1 indicates the fraction of the wall energy deposited into the BH. The wall stress energy tension is a source of repulsive gravity [58,64,65,68] and hence should maintain a negative pressure in case of a domination of the wall's material inside of encompassing it Hubble radius. The collapse can take place unless the the mass of the wall inside a given Hubble radius exceed the mass of its matter content which happens when the Hubble rate reaches the value Thus, for stress energy density (90) relevant for the ALP parameters localized in scenario (iii), discussed in Section 7, the dominating wall "enters" 12 the horizon at Hubble rate H w(iii) ≈ 10 −9 eV. Comparing this rate with that one at which f yr enters into Hubble radius (54) one can see that the walls of size corresponding to the frequencies f w(iii) 300 nHz are able to collapse into PBHs. Bigger walls entered the Hubble radius starts expanding faster than the background reaching eventually the inflationary vacuum and develop wormholes to baby universes [61,62].
Such wormholes are seen as BHs in the FRW Universe. Since this frequency is much higher than the frequencies of the signal bins of the NANOGrav, the ALP domain walls, described in (iii), should mostly form wormholes. Provided that the condition (51) fails in each scenario 12 It is more likely to say, that the size of the wall coincides with the Hubble radius.
of Section 7, the acoustic response signal (50) induced by these escaping walls cannot saturate the amplitude (28). Therefore, most likely the ALP field with parameters specified in scenario (iii) cannot be a source of the signal indicated by the NANOGrav. In scenarios (i) and (ii), the dominating wall enters the horizon at Hubble rate H w(i) ≈ 10 −11 eV, which corresponds to the GW signal frequency of f w(i) 3 nHz. Therefore, the domain walls induced by inflationary dynamics of the ALP fields specified in scenarios (i) and (ii) should collapse into PBHs within almost the whole range of frequency support of the NANOGrav. The slope of the strain spectrum (44), which might be generated by these walls, γ = 5.5, exhibits a remarkable agreement with the central value of the slope range (41) reported by NANOGrav [42], in its excess frequency bins. Generically, it would be reasonable to expect a sort of spectral feature expressed as the slope change from γ = 3.5, for frequencies bellow f w(i) 3 nHz to γ = 5.5 for higher frequencies. Such kind of feature might indicate the transition from the regime when bigger walls, "entering" their Hubble radii, were escaping into baby universes and those smaller ones which were collapsing into PBHs. However, because of the amplitude inconsistency between collapsing (44) and escaping (50) strains, caused by the violation of condition (51), 13 one might expect to see only a spectrum of slope γ = 5.5 with amplitude almost abruptly attenuating at frequencies bellow f w(i) 3 nHz.
Currently, NANOGrav has about T ≈ 15 years of high precision timing observations from many pulsars, where every pulsar is observed each ∆t 1÷3 weeks with integration time about 20 minutes. Thus, the sensitivity band to GW frequency can be defined as 1/T < f < ∆t/2, which spans the range from 2 nHz to 1 µHz. This means, that in scenarios (i) and (ii), the frequency range indicating the domain walls forming wormholes, may be, just barely reaches the lowest frequency bin accessible for the NANOGrav. Therewith, it would be hard to expect to obtain a real testification on wormholes formation by ALP field induced domain walls.
The boundary between a PBH collapsing and wormhole escaping domain wall can be characterized by the mass of wall material contained within the size H w , which reads their masses cannot exceed the energy of the radiation fluid contained in a respective Hubble radius [61], provided that horizon entering took place before the equality epoch.
Closing the section, we notice, that the observational access to the signature of domain walls escaping into baby universes may be more encouraging in the context of the spontaneous nucleation of domain walls on the inflationary stage, considered in [61][62][63]. In this setup, the nucleation rate Γ s ∝ exp(−S E ) is defined by the action S E of the semi-classical tunneling pass, To instance, for higher reference frequency f 0 > f yr , one may represent (50) as Therefore, the signature of domain wall escaping into baby universes might be seen with a detector of GWs having sensitive to amplitude in the characteristic strain power spectrum like (39), measuring the slope value to be close to that one of the acoustic response, γ = 3.5. Presumably, Gaia and THEIA can have some potential to provide sensitive measurements in much higher frequency band [77].

Consistency of collapsing walls evolution with PBHs constraints
There is no domain wall problem related to their production mechanism considered in this paper. The walls which create wormholes export the domain wall problem into baby universes.
The collapsing walls form PBHs before their contribution into the energy density becomes large enough to contradict the observational constraints.
The dynamics of spherical domain wall has been studied in [64]. The result of this study implies that a closed domain wall being initially bigger than its Schwartzshild radius always collapses into a BH. Thus, the ALP field induced domain wall after obtaining of a spherical geometry contracts toward the center. Since the ALP's coupling to the Standard Model particles is strongly suppressed by the large scale F , so that the ALP wall interacts very weakly with matter and hence nothing can prevent the wall finally to be localized within a volume of size comparable to its width [32,33]. The Schwarzschild radius of a wall of total mass M w can be expressed trough the radius of the wall R w and the ALP parameters as follows Provided that the wall contracts down to about the size comparable to its width m −1 θ , one naturally assumes that a PBH will be formed under the condition R S > m −1 θ , which reads Since, the NANOGrav frequencies would correspond to the size of the emitting object one can impose the following constraint on the ALP's mass where we used (53) to evolve (98) and (99). The numerical value of the constraint (100) reads m θ(i) ≈ 8.5 × 10 −7 eV, m θ(ii) ≈ 2.5 × 10 −7 eV and m θ(ii) ≈ 8.5 × 10 −10 eV for scenario (i), (ii) and (iii) respectively. Thus, one can affirm that in all three scenarios described in Section 7, the condition (100) is well satisfied, which implies that the domain walls producing GWs signals of frequencies about f yr , are capable to deposit their energy under their Schwarzschild radius and hence convert it into PBHs. It is obvious, that the above statement is still in tact also for frequencies one order of magnitude bellow and above f yr .
In what follows, we verify the consistency of PBHs abundance of different masses produced by collapsing closed domain walls, with existing constraints [40].
The mass distribution of the BHs is defined by the size distribution of the walls' seeding contours (13) scaled with respect to the expansion of the Universe, which is given in (23) and reads as at the equality time. A useful characteristic of this distribution, which allows to compare the PBHs yield with the constraints on their abundance in different mass range [40] 14 , is the mass density of PBHs per logarithmic mass interval in units of the total density of the Universe where ρ eq = M 2 Pl /(6πt 2 eq ) is the matter density at the time of equality. Using (101), we can obtain 14 In particular, see Fig. 18 in [40]. so that (102) can be converted into to be compared with the model used in [40] to quote the constraints on the density fraction β deposited in PBHs at the moment of their formation. In particular, independently of the NANOGrav signal, one can provide an upper bound, shown in the left panel of Fig. 3, on the combination σΓ 4/3 s for PBHs of different masses formed in a domain walls collapse process, provided that each wall converts its total mass into a PBH. One can see that the restrictions (see left panel in Fig. 3) become gap-like tougher, namely σΓ 4/3 s ≤ 10 −7 , in a range of masses exceeding the boundary mass M w (H w ) = 17M quoted above.
In general, the relation (91) implies that different combinations of R and σ may be composed into one and the same mass of the PBH. This degeneracy is removed if we fix the reference frequency f yr of the NANOGrav signal. Thus, using the chain of relations (104), (91), R = ct, and the connection 1 + z = (1 + z eq )t 1/2 eq /t 1/2 (107) one arrives to and The expressions (108) and (109) together with constraints presented in [40] impose the bound on Γ s , which is shown by dogleg lines in the right panel of Fig. 3.
Using (44) and (40) one can infer the mutual correlation which NANOGrav signal imposes on σ and Γ s , in case of saturation of the signal amplitude This equation is shown by the straight line in the right panel of Fig. 3.
Thus, one can see that the PBHs contribution from the collapsing ALP domain walls saturating with their GWs emission the NANOGrav signal is well bellow the allowed limit on the PBHs abundance.

Conclusions
In this paper, we have explored the impact of the evolution of closed domain walls, created in course of dynamics of ALP field spectating the inflation, on the stochastic GWs background in the frequency range accessible for the PTA measurements.
As long as the Universe is inflating, an axion or ALP field owing, by definition, the discrete set of vacuum states, with PQ (-like) symmetry broken before inflation, experiences the quantum fluctuations. In some patches of the space, an accumulation of the fluctuations can move the field, from its initial value (misalignment angle), so much that after inflation it will be located in a vacuum state differing from that one corresponding to its value occurred at the beginning of the inflation. Thus, due to inflationary expansion and provided that the probability for the field to acquire the value corresponding to another vacuum state is small enough, one expects to see a Swiss cheese like picture where the domains with underrepresented vacuum will be inserted into a space remaining filled with the initial one. Two different neighboring vacuum states should be interpolated by a domain wall, represented by the sine-Gordon kink solution (7) in the case of axions or ALPs. Obviously, the Swiss cheese structure guaranties that the produced domain walls have closed, although irregular, shapes 15 .
Evolving after inflation, the domain walls, depending on their sizes, either collapse forming PBHs or escape into baby universes leaving wormholes anyway observed as BHs in our Universe.
Therethrough, there is no domain wall problem in either scenario. The collapsing walls tend to decrease their entropy leading to a smoothing out of their surfaces, and hence radiating our the energy in form of GWs, with characteristic frequency of about the Hubble rate established during their collapsing instants. We have estimated the characteristic strain power spectrum produced by the size distribution of the collapsing closed domain walls and relate it with the recently reported NANOGrav signal excess obtained in PTA measurements. It is remarkable, that the slope of the frequency dependence of the strain spectrum γ = 5.5, generated by such domain walls, is very well centered inside the range of the slopes in the signal obtained by the NANOGrav.
Analyzing the inflationary dynamics of the ALP field, in consistency with the isocurvature constraint, we defined those combinations of its parameters where the signal from the inflationary induced ALPs domain walls could saturate the amplitude of the NANOGrav excess.
In particular, if an ALP has the misalignment angle value set to the same one, imposed by the anthropic window of the QCD axion, it may saturate the NANOGrav excess amplitude even having a negligible contribution into the CDM density, at the same time being consistent with the isocurvature constraint. Provided that string theory models that include the QCD axion, generically at high scale of PQ symmetry breaking, might also incorporate other ALPs which can leave their imprint in the stochastic GWs background, with the preferred by the NANOGrav excess parameters of the characteristic strain power spectrum and without overproduction of CDM and/or isocuvature modes.
Those walls exceeding by their mass the energy of the bulk fluid inside the Hubble radius will dominate by their repulsive nature and hence push the fluid away, leaving two nearly empty layers in their vicinity. Such walls, keeping exponentially extending, escape into baby universes forming wormholes in our Universe, which are observed as BHs of masses ranging from planetary to some of LIGO/Virgo scales [80,81], depending on the parameters of ALPs. We assumed that being pushed away the bulk radiation fluid responds by acoustics waves which may set in a motion a fraction of the mass in Hubble radius. This mass, in a motion, could generate the GWs signal with the strain spectrum slope and the amplitude different from that ones expected from the collapsing walls. It was found that in the context of domain walls produces by ALPs, it is quite unlikely, while in principle possible, to trace baby universes in PTA measurements. However, a detectable stochastic GWs imprint of the phenomena can be specified for other types of phase transitions which might occur during inflation.
At some final inflationary e-folds, corresponding to much smaller than the NANOGrav scales, the probability of the ALP field to move its value into another vacuum may become large, which may correspond to percolated domain walls system considered in [82] and induce small scale density inhomogeneities in the distribution of the coherent oscillations of the ALPs Bose-Eisenstein condensate. Possible scenarios of evolution and observational signatures of small scale inhomogeneities in axionic Bose-Eisenstein condensate have been studied in great details in the cosmology of QCD axion (see for example [83][84][85][86][87][88]).

A First passage e-fold
Providing the condition p(∆θ π , N ) ≈ 1, in (57), one can define how many e-folds would it take the phase to cross π, solving the following equation One can put the term 3 2 ln N in the parenthesis to a constant C 1 ≈ 4.5 ÷ 6.2, provided that it changes quite slow with increasing of N in the range N ≈ 20 ÷ 60. Under this assumption, the solution of (111) can be approximated by the root of quadratic equation, given by N π ≈ (24π 2 ∆θ 2 π Q 2 FH inf − (0.9 ln ∆θ π + ln Q FH inf − C 1 )) 1/2 6 − (0.9 ln ∆θ π + ln Q FH inf − C 1 ) 2 6 .
(112) For reasonable choose of parameter, 1 < Q FH inf 10 and ∆θ π π/2 and provided that typically accepted number of e-folds needed for inflation N inf 60, one can neglect relevant terms in (112), so that the solution is reduced to