Hydrogen Axion Star: Metallic Hydrogen Bound to a QCD Axion BEC

As a cold dark matter candidate, the QCD axion may form Bose-Einstein condensates, called axion stars, with masses around $10^{-11}\,M_{\odot}$. In this paper, we point out that a brand new astrophysical object, a Hydrogen Axion Star (HAS), may well be formed by ordinary baryonic matter becoming gravitationally bound to an axion star. We study the properties of the HAS and find that the hydrogen cloud has a high pressure and temperature in the center and is likely in the liquid metallic hydrogen state. Because of the high particle number densities for both the axion star and the hydrogen cloud, the feeble interaction between axion and hydrogen can still generate enough internal power, around $10^{13}~\mbox{W}\times(m_a/5~\mbox{meV})^4$, to make these objects luminous point sources. High resolution ultraviolet, optical and infrared telescopes can discover HAS via black-body radiation.


Introduction
in experiments like CASPEr [28] (see Ref. [29] for a recent review). ADMX in particular has recently claimed the first bounds of the QCD preferred region of axion parameter space [22], though the application of their results to the QCD axion relies on assumptions about the abundance and local distribution of dark matter that we examine critically below.
Another potential avenue to probe axion dark matter exists that does not necessarily rely on axion coupling to the SM exists. Assuming that the end of inflation scale is above PQ-breaking scale and before the QCD phase transition, the axion field has O(1) fluctuations on a length scale of H −1 10 km [30]. This large variation on small distance scales, along with the high occupancy number of axions, allows for axions to relax into a high density, compact Bose-Einstein condensate (BEC) clump conventionally known as an axion star or minicluster [31][32][33]. The exact properties of such a clump [30,[34][35][36][37][38][39][40][41], as well its stability [42][43][44][45] has come under scrutiny recently. Our determination is that the most likely final state for primordial QCD axion stars is a dense axion star with meter scale radius and order 10 19 kg (10 −11 M ) mass. Such a conclusion remains to be fully tested via numerical simulations that are beyond the scope of our work.
Rather, we focus on the observational manifestations that can result from the likely structure of primordial axion stars. Some potential signals specific to axion stars have been studied. Axion stars are sufficiently dense that multi-axion interactions are possible, yielding boosted signals from mutliaxion processes n a → 2 a and photons from n a → 2 γ, where n > 2 is the number of incoming axions in the process. Collisions of axion stars with a dense baryonic object such as a neutron star could also generate photons and could potentially explain fast radio bursts [46,47]. Collapsing of axion stars could undergo maser-like axion production.
In this article, we study signals generated by a primordial hydrogen cloud that forms around the axion star, forming a HAS. The axion star, being a massive, dense object, seeds small scale baryonic structure formation at a very early time. Hydrogen and helium gas gravitationally collapse on it and reach a cloud in hydrostatic quasi-equilibrium, balanced by thermal pressure. As the universe expands, several changes happen. The gas surrounding the cloud gets colder and collapses further, leading to growth of the cloud. The cloud properties change from insulating gas to conducting fluid, allowing the conduction band electrons to convert axions and heat up the cloud. The hot cloud cools back down as it expands, as the generated heat is spread wider within the cloud. The hot cloud radiates blackbody photons from its surface, which yields the observable signal with which we are concerned.
The properties of such clouds at low redshifts below z = 10 are dependent on the axion star surroundings and are no longer uniform due to reionization, collapse into galaxies, and heating from astrophysical sources. Nevertheless, so long as the cloud isn't destroyed by a collision with an unrelated molecular cloud, the power emitted should remain unchanged at roughly 10 13 W (for 10 −3 eV m a 10 −2 eV) and be detectable by a telescope with an appropriate frequency band. We study the potential sensitivity of such a telescope to a nearby baryon-clouded axion star, axion stars in a dwarf spheroidal satellite and diffuse emissions from baryon-clouded axion stars. The best sensitivity would come from a nearby axion star point source, which could likely be within a few hundred astronomical units of the Earth in distance. The signal should be much larger than the diffuse background if the axion star is located in a direction away from the galactic plane and should be readily visible to a high resolution The remainder of this article is structured as follows. In Section 2, we review the properties of QCD axion dark matter and determine the most relevant parameter space for the axion of QCD or any new strong dynamics sector. Then, in Section 3, we review the formation of axion stars and determine the most likely properties of primordial axion stars. We turn to the structure of the baryonic cloud around axion stars in Section 4. We study the potential for observing HAS in Section 5. The consequences of an axion star inside a solar system planet are studied in Sectin 6. Finally, in Sections 7 and 8, we discuss uncertainties in our analysis and conclude. Appendix A presents some details of the conversion of axions to heat in a metal.

Axion as Cold Dark Matter
In this section, we present relevant aspects of early universe axion cosmology, as well as the large scale features of axions as dark matter. In particular, we discuss the various mechanisms for axion production in the early universe and the parameter space for axions as dark matter. We briefly present a model that generalizes the axion parameter space beyond the QCD axion, allowing for a wider range of axion phenomenology. Helpful reviews of axion cosmology can be found in Refs. [8,9].

Axion Couplings and Lifetime
In order for the axion to be a viable cold dark matter candidate, it must first and foremost be stable.
Axion couplings to SM particles are expected to be extremely weak, though model dependent. The relevant interactions for our discussion are Here, ψ represents SM fermions. The coefficient c γ is related to the triangle anomaly of U (1) PQ with two electromagnetic currents. For certain models, the anomaly could be highly suppressed such that the axion lifetime can be increased. The axion decay width is with the decay lifetime of (3)

Energy Density from the Misalignment Mechanism
In order for an axion star to form, the PQ symmetry breaking scale must be below the inflation scale.
In parts of parameter space for which the PQ breaking happens before or during inflation, the small scale density perturbations required for primordial axion stars are washed out and the axion density is nearly uniform on small length scales of order the Hubble length. On the other hand, if PQ breaking happens after inflation, then initial angles for the axion field are only uniform on a scale of order Hubble, with different Hubble patches having different angles sampled uniformly. As a result, the non-linear perturbations of axion field can easily generate localized clumps or miniclusters. This sets an upper bound on the decay constant of f a H I /(2π) ≈ 10 13 GeV [52]. There is also a lower bound of f a 10 8 GeV from stellar energy-loss limits [10,53]. So, for a light axion below O(1 MeV), the allowed f a range is For the axion misalignment mechanism and in the standard cosmology (see Ref. [54] for other non-standard cosmology models), the axion field equation of motion is for θ(x) ≡ a(x)/f a and with the Hubble rate given by in the radiation dominated era. Here, g * (T ) is the temperature-dependent number of relativistic degrees of freedom and the Planck scale M pl = 1.22 × 10 19 GeV. At earlier times with higher temperature, one can ignore the axion mass and have a time-independent solution for θ(x) (the other solution proportional to 1/ √ t decays away by later times and can be ignored). The axion mass term becomes important at a temperature of T 1 with 3 H(T 1 ) = m a (T 1 ), which requires us to know the temperature-dependent axion mass.
Based on a dilute instanton gas calculation [8,55,56], the axion mass increases as the temperature decreases with a power-law behavior. An updated study based on the interacting instanton liquid model shows a slightly higher mass for temperatures above the confinement scale [21]. More recently, a calculation based on lattice simulations has shown that the axion mass could increase much faster than the conventional power-like behavior, although with a large uncertainty [57,58]. In this study, we will simply follow the traditional diluted instanton gas calculation in Ref. [56] to estimate the oscillation temperature.
At zero-temperature, the axion mass is where Λ QCD is the QCD strong coupling scale and N 0 f is the number of light quarks with m i < Λ QCD . At higher temperature, following the calculation in Ref. [56], the axion mass as a function of T is Here, N f as the number of flavor below T and the function I(T, Λ QCD , N f ) can be found in Eq. (6.15) of Ref. [55]. After numerical integration, the axion mass can be fit by the following formula The numerical values of α and β for different numbers of light flavors are listed in Table 1.  The misalignment mechanism works as follows. Given the initial constant value of the axion field in a given Hubble patch and the negligible axion potential well before the QCD phase transition, the energy density of the axion field is initially negligible. At a temperature T 1 , the axion mass becomes larger than Hubble and the axion field begins to roll toward its minimum, while the energy density in the axion field grows for this time-dependent Hamiltonian. Around the QCD phase transition temperature, the energy density in the axion field per comoving volume freezes out and the axion field has relaxed to its minimum. The resulting energy density of axions at that time is where f (θ i ) takes into account of the anharmonic terms in the axion potential and other O(1) effects.
The energy density scales with the entropy density to its present day value of [9] ρ where the entropy density is s(T ) = 2π 2 g * S (T )T 3 /45 and the Hubble patch averaged initial angle dependence is given by θ 2 i f (θ i ) ≈ 8.77 [59]. The dark matter energy fraction parameter Ω a = ρ a /ρ crit with ρ crit = 3H 2 0 M 2 pl /8π, to be matched to the experimental value of Ω CDM h 2 = 0.1199 ± 0.0022 and h = 0.678 ± 0.009 [60] if the axion is the only component of dark matter.

Energy Density from Axion String Decay
As alluded to in the previous section, the value of the axion field is only constant on roughly the Hubble scale. At the PQ phase transition, topological defects may form on that length scale. We focus on the contributions to the axion abundance due to axion string defects. 1 Axion strings can be produced in early universe by the Kibble mechanism [63,64] and contribute non-negligibly to the energy density for the U (1) PQ breaking after inflation. Two types of decaying strings produce axions: loop strings and long strings. Detailed axion production results rely on the typical loop size and the back-reaction rate [9,21]. Depending the axion energy spectrum from string radiating axions, the ratio of axion production from string decays to that from mis-alignment can vary up to two orders of magnitude. The abundance of axions generated nevertheless has the same axion mass and decay constant parametric dependence [21], so that we can simply parametrize the theoretical uncertainty in string decay. We therefore define the axion parameter-independent ratio r string = ρ string−decay ρ mis−align , and take the summation of the two contributions as the total axion energy density. In Fig. 1, we show the parameters m a and f a required for axions to account for all the cold dark matter in the universe.
Generically, if the contribution from string decays is increased, a lower value of f a is needed for a fixed value of m a . We also show the SM QCD axion relation of m a f a ∼ (0.1 GeV) 2 in the blue straight line of Fig. 1.   The decay constant f a is bounded from above by the Hubble scale at the end of inflation [52] and from below by the stellar energy loss [53]. The parameter space on the right side of EBL & X-ray curves are excluded for c γ = 1, where c γ is defined in Eq. (1). The region on the right side of the gray and straight line labeled by τ a > 1.9 × 10 25 s is excluded by reionization effects on CMB observables. The gray and dashed lines are for c γ = 0.01 and a smaller coupling of axion to two photons. The two benchmark points have: (m a , f a ) = (5 × 10 −5 eV, 2 × 10 11 GeV) (star) and (5 × 10 −3 eV, 2 × 10 9 GeV) (cloud).
Also shown in Fig. 1 are the constraints on the decay constant from stellar energy-loss limits [53] and requiring it below the Hubble scale at the end of inflation [52]. For a heavier axion mass between 1-100 eV, even though the axion lifetime is greater than the age of our Universe, the axion still decays fast enough to have non-trivial contributions to the extragalactic background light (EBL) and X-ray.
Using the summarized results in Refs. [65,66], we show the constraints on our axion model parameters in gray and curved lines. Furthermore, the decays of axions in the early Universe add additional contribution to the radiation energy and reionization effects. The axion particles decay 100% into electromagnetic components. The cosmic microwave background (CMB) observations impose a nontrivial bound on its decay lifetime to be τ 0 a > 1.9 × 10 25 s. We translate this constraint as gray and straight lines for two different choices of c γ = 1 and 0.01 in a later Eq. (3). All together, the axion can account for all of cold dark matter and satisfy experimental constraints if its mass falls in the range O(10 −5 ) − O(100) eV. Generally, we consider models in which SM QCD generates the dominant potential for the axion. It may be phenomenologically interesting, however, to consider a more general parameter space with more massive axions. In order to move away from the QCD axion line in the plot, additional contributions to the axion potential are required, which we briefly discuss below.

Heavy Axion Model
In a typical axion model, the U (1) PQ symmetry is spontaneously broken such that the corresponding pseudo-Nambu-Goldstone boson φ has the shift symmetry of φ → φ + f a under the infinitesimal U (1) PQ transformation. The U (1) PQ is also anomalous under the SM QCD interaction with a non-zero triangular anomaly, A c PQ ≡ Tr(t a t b Q PQ ). The QCD instanton effects generate an effective potential for the axion field V (φ) ∼ m 2 π f 2 π cos(φ/f a − θ). Minimizing this potential, one has φ = f a θ, which makes the physical strong-CP parameter θ eff = θ − φ/f a zero and solve the strong-CP problem.
For the "invisible axion" models with f a f π , the mixings of axion with η , η and π 0 are suppressed by f π /f a and are small. The axion mass is model-dependent and related to the quark masses as [67] So, the axion mass m a and its decay constant f a are not independent parameters and are related by SM QCD scales, up to order-one model parameters.
Beyond the minimum SM QCD axion models, additional models aim to increase the axion mass.
For instance, in Ref. [68], a hidden QCD sector is considered to have additional instanton contributions to the axion mass (see also Ref. [69] for a recent study). The effective potential for the pseudo-Nambu- where c and c are model-dependent number of order unity and θ is the θ-angle of the hidden QCD sector. To solve the strong-CP problem, one also need to require θ = θ . This can be achieved by assuming a Z 2 symmetry in the UV theory such that the θ-angles in both sectors are identical. After imposing the vacuum expectation value (VEV) of φ field, we have the potential for the axion field as For the parameter space with c m 2 π f 2 π > c m 2 π f 2 π , the axion mass is dominated by the hidden sector contributions with m a f a ≈ √ c m π f π . Absent a picture of hidden QCD dynamics, we should treat m a and f a as independent parameters. As seen above, requiring axions to account for the entire dark matter density fixes one relation between m a and f a , up to theoretical uncertainties.
To calculate the axion relic energy density, it is also important to know the time-dependent mass, which cares about the hidden quark masses in the hidden QCD sector. One nature choice based on the Z 2 symmetry is to have the hidden QCD with N c = 3 with the total number of flavor N f = 6.
The Yukawa couplings in the hidden sector are identical to the visible sector: y i f = y i f . On the other hand, the QCD confinement and electroweak scales in both sectors could be different. Generically, we will assume v > v and Λ QCD > Λ QCD to have the relation c m 2 π f 2 π > c m 2 π f 2 π in Eq. (15).

Axion Stars
For the two main non-thermal axion production mechanisms, the axions are extremely cold after the QCD phase transition 2 and may form a Bose-Einstein condensate. Furthermore, because of the large inhomogeneity of axion energy density at or below the oscillation temperature T 1 , the group of axions that are causally connected within the Hubble volume at T 1 can be gravitationally bound to form axion stars, sometimes also called axion miniclusters [31,70,71]. Numerical simulations are in general needed to answer questions like the fraction of axions belong an axion star and the mass and size of the axion stars. Some initial attempts to answer those questions can be found in Refs. [70,72]. In this paper, we will mainly concentrate on the potential phenomenological consequence of axion stars and will only have estimations of some axion star properties at the order-of-magnitude level. The general picture of axion stars arises from the following timeline: • Reheating after inflation places the axion field at a temperature above the PQ phase transition that occurs at T PQ ∼ f PQ . There are small inhomogeneities due to quantum fluctuations in the (now decayed) inflaton.
• At T PQ , the axion field (the radial field) rolls down to it's minimum in a random U (1) PQ phase direction that is homogeneous on a scale of roughly H(T PQ ) −1 . Excitations along the flat direction at the bottom of the Mexican-hat potential are the axions. At this stage, they can be treated as a massless classical scalar field in an expanding background. Differences in the vacuum phase throughout space lead to topological defects, namely strings.
• As the Universe cools below the phase transition, the axion strings can split off loops that decay and can also radiate axions. Production is dominated at low temperatures and so low energies.
The thermal abundance of axions is subdominant.
• Around the oscillation temperature T 1 defined by 3H(T 1 ) = m a (T 1 ), the potential for the axions becomes sufficiently large that they roll down toward their minimum, damped by the Hubble friction. After a few oscillations around the minimum, the axion field relaxes to its minimum value. This rolling causes copious axion production. In addition, as the mass becomes large, the string network annihilates and radiates to now massive axions and disappears.
• At the time of the QCD phase transition, where T QCD is slightly below T 1 , the axion string network has disappeared, the axion mass is frozen at its zero-temperature value, and the abun- The axion quartic self-interaction, and later the gravitational self-interaction, are sufficient to cause the axion to locally relax into a BEC [30]. The resulting BEC axion star contains at most the axions within roughly a Hubble volume at the QCD phase transition and is coherent at most on a Hubble length at the QCD phase transition, though we will show that the typical size of the axion star is far below that. The ground state of the system is set by the gravitational and self-interaction potential of the clump. Almost all the axions in the collapsing region end up in the ground state, as it is a BEC state. We now estimate the properties of the resulting axion clump. After the number of axions per comoving volume has frozen, the number density of axions is fixed by requiring that they account for the total present day energy density of dark matter Ω CDM h 2 ≈ 0.12. For this to be the case, the number density must satisfy for the temperature from around T 1 to the present day temperature T 0 = 2.7 K, since the axion number per comoving volume is approximately conserved. The number of axions in a Hubble volume at the time of the QCD phase transition is thus Here, we have used g * S (T 0 ) = 3.91, g * S (T 1 ) ≈ g * (T 1 ) and T 1 ≈ 1 GeV from solving the relation m a (T 1 ) = 3H(T 1 ) in the dilute instanton gas approximation. This number of axions is roughly the amount available for gravitational collapse. We assume that the final number of axions in an axion star is around N a , even through additional numerical simulations are required to know the exact relation.
That is, we assume that the fraction of axions collapsing into axion stars is 100%. 3 In addition, the above estimate assumes an average Hubble patch. Clearly, some patches will be more or less dense than others; without such an inhomogeneity, collapse is not even possible. On the other hand, the variations compared to the mean are at most O(1) since the height of the axion potential is finite. The axion-star total mass is simply which scales as (m a f a ) −3/2 ∼ Λ −3 QCD . In Fig. 2 we show the distributions of N a and M a for different axion masses after requiring the axions to account for all dark matter. In Fig. 3, we show the population of axion stars using both the averaged dark matter energy density in the Universe and in the local solar system.
Relaxation can begin once the axion acquires a significant mass. The relaxation time to reach a BEC state can be estimated by the time scale for evolution of the occupancy number in mode k. As shown in [30], this evolution time is given by where n ave is the average axion number density and G N = 6.7×10 −39 GeV −2 is Newton's gravitational constant. The number density is at least the average number density of axions in the universe, n a (T ), which at T 1 is given by n a (T 1 ) ∼ 5 × 10 49 m −3 for an axion mass of m a = 5 × 10 −3 eV. At a reference wavenumber determined by the Hubble scale at the T 1 , relaxation is dominated by self-interactions; the relaxation time determined by gravitational interactions is more than three orders of magnitude larger. We then find that the relaxation time is given by for a QCD axion. Note that so the axion star relaxation is rapid around the QCD phase transition. Even though the relaxation due to self-interactions becomes slower, the gravitational relaxation of modes on the Hubble scale becomes more important and helps ensure rapid relaxation at all times. The states we consider thus fall into the "dense" regime studied in Ref. [38]. In this regime, the gravitational interactions are balanced not by quantum pressure, but rather by the self-interaction potential. The mass of the axion star, even in the dense regime, is dominated by the mass of the axions in the star, such that M a ≈ m a N a . The radius and mass of the axion star are roughly related which is independent of m a for fixed m a f a . Following Ref. [38], we solve for the wavefunction numerically under the Thomas-Fermi approximation, which is a good approximation in the region of interest.
The total energy is given by where Φ is the gravitational potential and V eff is the effective self-interaction potential for the axion in the non-relativistic limit. The potential can be derived in the non-relativistic limit by expanding the axion field a as a(x, t) = [ψ(x)e −imat + ψ * (x)e imat ]/ √ 2m a and collecting the terms that are time independent, i.e. slowly varying. The resulting potential, as in [43,45], is given by where J 0 (x) is a Bessel function. The equations of motion are the Schröedinger equation combined with the Poisson equation for gravitational potential Here, the space-independent chemical potential, µ, does not enter our approximation for the wavefunction. In the Thomas-Fermi approximation, the kinetic energy term is dropped, so that a single differential equation for ψ is obtained, namely The wavefunction depends only on the radial coordinate r under the assumptions we have made. The boundary conditions are taken to be |ψ(0)| = ψ 0 = n 1/2 0 with n 0 as the central number density and ψ (0) = 0. As clarified in Ref. [45], if V eff (ψ 2 0 ) < 0, then the density increases away from the origin and there is no solution that goes to 0 at infinity, that is no physical solution exists, so there are only viable solutions for V eff > 0. The second derivative of the potential is given by a regularized hypergeometric withn ≡ 2|ψ| 2 /m a f 2 a ; this has an infinite number of zeroes, leading to an infinite number of solution branches (one for each zero). If the axion star collapse is sufficiently slow, it should collapse down to the least dense stable state, which is the described by a radius in Eq. (22). If the collapse process is more violent, the axion star could shed mass or skip down to an even more dense state. A full simulation of the dynamics of axion star collapse is required to determine the ultimate axion star state, but for the purposes of this work we assume that the collapse is sufficiently slow that the least dense of the dense axion star states is reached. The axion star profile depends only on the axion star total mass, which we fix as above for the QCD axion and observed dark matter abundance to where M a ,enc is the enclosed mass of axions, which can be determined from its density profile in Fig. 4. We therefore solve for the equilibrium state of the baryons. We are also interested in this state after recombination and so consider a fluid of hydrogen. It turns out that the hydrogen will become sufficiently dense as to reach a state of high-conductivity metallic hydrogen. In the presence of this metallic hydrogen, the axions will convert off conduction band electrons, generating a non-trivial amount of heat. We therefore must solve for thermal, hydrostatic, and heat flow equilibrium.
Thermal equilibrium enforces that there is a locally well-defined temperature T , so that we can make use thermodynamic properties of hydrogen such as its equation of state for pressure P (ρ, T ), where ρ is the local mass density. Assuming that there is no net flow in the hydrogen fluid, the condition for hydrostatic equilibrium is then where Φ is the gravitational potential. The self-gravitation of the hydrogen gas will turn out to be negligible, so that the gravitational force per unit mass is given by We assume an isotropic distribution for axions inside the axion star as above. At the high pressures reached for baryons in the core of the axion star, there is a series of phase transitions during which the hydrogen gas goes from an atomic gas to a molecular gas (second order), then a gas to a liquid (second order) and then from a molecular liquid to an atomic metallic liquid (generally first order) [75]. For a relatively cold axion star, it is possible that some of the "cloud" will actually be solid hydrogen.
Nevertheless, the core will generally be quite hot and the solid layer very thin. We neglect any solid phases in our calculations below. The gas states are well described by the ideal gas equation of state where M gas is the mass of the gas particles, either m p for atomic gas or 2m p for molecular gas. Here, m p is the proton mass. For the gas far away from the axion star, it has an ideal gas equation of state with the temperature given by T m . At early times, the interactions of free electrons with CMB photons via Compton scattering maintains a uniform gas temperature that follows the CMB photon temperature. At later times, after redshift z t ≈ 137 [76], corresponding to a CMB photon temperature T γ t ≈ 0.032 eV, the gas temperature declines adiabatically as T m = [(1 + z)/(1 + z t )] 2 T γ t and is colder than CMB photons until photoionization heating when the first stars start to burn.
For the liquid states surrounding the axion star, we assume a common-power law equation of state We fit the parameters to a combination of calculations and data [77,78], finding ρ 0 = 51.33 kg/m 3 , T 0 = 1000 K, γ ρ = 1.36, γ T = −0.89 .
We approximate the gas to liquid phase transition as a first order transition at (ρ/ρ 0 ) γρ (T /T 0 ) γ T = 1.
This transition essentially defines the boundary of the baryon component of the axion star, as the density of the baryons grows very rapidly near the boundary, generally surpassing the required density to reach a liquid metallic state. We therefore define R bound to be the radius for which where T surf is the surface temperature of the axion star, which we revisit below.
The metallic liquid phase has a large DC conductivity of order σ 0 ∼ 2500 S/cm [79] (in natural units, S/cm = 0.0075 eV). The conductivity can be related to the absorption cross-section for axions.
The details of this calculation are presented in the Appendix A and have been studied by Refs. [80][81][82]. If the axions are sufficiently light, then frequency dependent effects in the conductivity can be neglected. The relaxation rate in metallic hydrogen is quite large, such that the onset of mass dependence should be well above our benchmark masses [83]. We make this approximation below.
For the two axion benchmarks, this leads to power generated in the core of order (see Appendix A for the relation of conductivity σ 0 and the axion absorption rate Γ a abs. ) The generated power is emitted as blackbody radiation near the surface of the HAS. Heat transportation within the HAS is rather complex and depends on the exact phase of the matter inside the HAS as well as the heat transportation properties. It is, however, likely that radiative or convective transportation will dominate. We do not need to solve this exactly, but rather can approximate the HAS as a black body, solving simply for the boundary condition that requires that the total generated power is equal to the total power emitted at the surface where T surf is the surface temperature of the HAS; σ SB is the Stefan-Boltzmann constant; the emissivity is chosen to be = 1. In order to do so, we must assume that the equations of state inside the HAS remain unchanged. If the temperature inside the HAS were to get so high that the fluid inside undergoes an additional phase transition, say to a plasma, then the total power generated could be modified.
We leave the detailed properties of the interior of the HAS, such as the exact temperature profile and the opacity, undetermined, assuming only that we know the density profile and conductivity.
The conductivity is assumed to be that of fluid metallic hydrogen as determined by Ref. [79]. We also assume the fluid hydrogen equation of state Eq. (32) holds throughout the interior of the HAS.
Aside from ensuring that our approximations regarding the density and conductivity are justified, the structure of the interior of the HAS has little bearing on the main observables we consider.   Table 2: Properties of the baryonic cloud boundary surrounding an axion star at different stopping redshifts in the evolution of the universe. effects like collapse into galaxies and reionization could have a significant effect on the clouds surrounding axion stars. We therefore do not attempt to evolve the HAS beyond this point and ascribe a substantial uncertainty to the final surface temperature of the cloud. It is also worth noting that if the HAS were to collide with a molecular cloud, it could be significantly disturbed in a way reminiscent of the collision of bullet clusters. The HAS cloud would likely get absorbed into the molecular cloud, while the axion core would move through unperturbed. To see this, we perform a rough estimate of the scales involved here. For a HAS near the Earth, the velocity relative to any nearby objects is of order 100 km/s, so that the kinetic energy transfer in a collision with a gas molecule is of order ∆E ∼ m H v 2 ∼ 2 × 10 −17 J. Meanwhile, the binding energy for a gas molecule to the axion star is of which is much less. It remains to determine the probability that a collision between a HAS and an unrelated gas cloud leads to most of the hydrogen in the HAS interacting. The number of interactions for a single hydrogen atom of HAS to interact with the gas cloud is of order π a 2 0 R GC n GC , where a 0 is the hydrogen Bohr radius, R GC is the radius of the gas cloud, n GC is the density of the gas cloud. To get an idea of the probability of destruction in a given event, we consider two benchmark cases. First, consider a collision with a giant molecular cloud, which has a density of 10 4 -10 6 cm −3 and a radius of 2-100 parsecs [84]. The resulting number of interactions per hydrogen atom is then at least roughly 5 × 10 6 , so that all the hydrogen in the HAS is likely to interact and disrupt the hydrogen cloud of the axion star. For a smaller system, like a solar system which has a diffuse density of order 5 cm −3 and a radius of order 100 AU [85], the expected number of interactions is order 1, meaning that an appreciable fraction of the baryonic cloud may survive such an interaction. Therefore, a fly-by HAS is still possible to behave as a point-like source.

HAS Detection
Radiation from the hot hydrogen fluid surrounding axion stars could be seen by a suitably high resolution telescope. Provided that the HAS fluid reaches a sufficient density to become metallic and that the metallic fluid reaches the edge of the axion core, the power generated by the HAS is independent of the radius of the baryonic cloud. In particular this means that the total power emitted is insensitive to the ultimate surface temperature of the HAS. On the other hand, the hotter the HAS, the lower wavelength (higher frequency) the emitted radiation will be. Depending on the final temperature of the HAS, the ideal wavelength for a sensitive telescope is in the UV through IR range.
Given an HAS at a distance D from the detector, surface temperature T , and baryonic cloud radius R bound , the rate of photons passing through a detector is given by There are several different strategies one could use to search for HAS photon emissions. One could look for a nearby HAS, which should appear as a fairly bright point source in the sky. The apparent magnitude (defined with respect to the brightness of Vega) of such a source would be m = 18.9 − 2.5 log 10 P 2 × 10 13 W 300 AU D Given the local density of dark matter, we expect there to be axion stars and HAS within roughly 300 AU of the Earth, which is near enough to give a detectable flux. If the telescope points away from the galactic plane, the estimated diffuse background can be relatively small [86]. The flux for a nearby HAS is shown in Fig. 6. Depending on the telescope angular resolution and observation time, the HAS located at a distance of the solar-system size can be observed.  [86]. We present the signal flux for two different choices of axion mass, m a = 5 meV (left panel) and m a = 0.5 meV (right panel). The total flux scales as m 4 a for a fixed z s , while the temperature for a lower axion mass is reduced due to decreased power generation. As a result, the peak wavelength is larger for a lower mass.
One could also look at an area of the sky rich in dark matter and low in background. The best such option given the frequencies emitted by the HAS is a baryon-poor dwarf spheroidal galaxy. The dwarf spheroidal with the smallest luminosity to mass ratio is the Ursa Major II dwarf, which has a luminosity of ∼ 4000 L , a distance of 30 kpc and a mass of around 5 × 10 6 M [87,88]. Given that mass, there are a total of roughly 2.8 × 10 17 axion stars in the galaxy, leading to a total luminosity of For our benchmark model, this luminosity is comparable to that of the dwarf galaxy, so a study of the spectrum of dwarf spheroidal could be sensitive to HAS.
Finally, one could simply look for diffuse light coming from HAS. The constraints from searches for diffuse photons should be fairly weak. To see that this is the case, we perform a simplified conversion of HAS radiation to effective dark matter decay into a pair of photons. While a dark matter particle decaying to a pair of photons would yield monochromatic photons of energy m DM /2, our axion decays to photons with energy peaked at a temperature T ∼ 1000 K -50, 000 K corresponding to photons of energy roughly k B T ∼ 0.09 eV -4.3 eV, such that we can take the effective mass to be m eff c 2 ∼ 0.18 eV -8.6 eV. The total rate of photon emission per axion is given by roughly the total power emitted divided by the typical energy of emitted photons and the total number of axions in the axion star, Γ eff P a /(k B T N a ), giving an effective lifetime in the range τ eff ∼ 2.9×10 24 s -1.5×10 26 s.
Constraints in the m a -τ a plane for monochromatic photons have been presented in Ref. [66]. No constraint is found for the effective mass and lifetime ranges above.

Axion Stars in Planets
Axion stars are have fairly low masses and, as such, would be small gravitational sources in a mature solar system. A remote axion star, isolated from baryons would be a small gravitational source and would only be potentially observable by decay into photons [43,45,89]. In this work, we have focused on signals from axion stars surrounded by baryonic matter. We have already discussed the potential for observing the hot baryonic cloud around a nearby axion star in the previous section, but another possibility is that axion stars end up inside planets or stars themselves. This scenario is less likely, as the processes for axion star capture are fairly weak and a planetary collision with other objects could cause the axion star to be pushed away from the center of the planet. In this section, we nevertheless briefly discuss the prospects and signals for such a scenario.
An axion star that ended up in a proto-stellar cloud would likely get stripped of its baryonic cloud.
It therefore would have few means to lose energy and end up in the center of a planet. There are processes with three body gravitational interactions that could cause a (proto-)solar system to capture an axion star [90][91][92][93][94][95][96]. With some degree of uncertainty, the conclusion of studies of such processes indicate a fairly small abundance of bound dark matter. If the axion star did end up captured by the proto-stellar cloud, it could end up seeding the formation of a planet, or the star itself, as it would be a large density perturbation within the cloud [97]. As seen above and in Appendix A, axions can convert off conduction band electrons into phonons and generate heat. A planet containing an axion star would have an additional heat source due to axion conversion in the metallic planetary core, which is likely to be a conductor. The heating occurs much as in a HAS, except with a smaller conductivity for which the frequency dependence becomes important.
If an axion star was at the core of the Earth, the heat generated by axion conversions could be comparable to the internal heat, 47±2 TW, generated by radioactive decays and primordial heat in the Earth [98]. This is particularly true for high mass axions, which have an enhanced conversion rate in iron. The relaxation rate in iron is significantly lower than in hydrogen such that frequency dependence is important for even the QCD axion. We therefore use the frequency-dependent conductivity for iron [99] to estimate the power generated by an axion star in the core of the Earth, which is shown in Fig. 7. For a heavier axion mass of 5 meV, the additional heat from axion absorption is around 1 TW and just one order of magnitude below the Earth total internal heat. For a lighter axion mass of 0.05 meV, both the suppression factor of m 2 a /f 2 a and the frequency-dependent Iron conductivity make the generated heat negligible. If the axion stars exist in other smaller size planets, exoplanets or moons, the axion-star-generated heat could be easier to notice.
In addition, the presence of an axion star at the center of a planet would alter the density-radius relation probed by seismology. By measuring seismic wave propagation at various points along the surface of the Earth, a reasonably accurate picture of the interior of the Earth has been inferred [100,101]. With sufficient resolution of the core, a large enough discrepancy between observation and reasonable models could be detected. Such a discrepancy could be explained by an axion star core, which would tend to increase the density of baryonic matter in the core of the axion star.

Discussion
We have found a very interesting and new class of astrophysical objects based on the well-motivated QCD axion model, but there remain several largely uncertain components of our analysis. The preferred axion parameter space depends both on a better understanding of the thermal axion mass and the generation of axions from the decay of topological defects. The choices of axion mass and decay constant have large, fourth power effect on the power generated by axion star interactions with a surrounding metal, so the viability of the signal presented in this work relies on a large axion mass implied by copious production in axion string decays. In addition, a complete understanding of the formation of primordial axion stars requires a full numerical simulation of axions around the time of the QCD phase transition. We have assumed that all the axions collapse into axion stars for our calculations, but axion star formation efficiency will be less than one. Some simulations have shown a formation efficiency of 30% to 50% [72], but it would be interesting to verify this finding with modern computing power and techniques. Furthermore, the process of collapse could cause the axion star to shed axions or to collapse to an even more dense state than that assumed in this work. Other possibilities due to axion self-interactions, such as a bosenova [43], also cannot be ruled out at this time. Finally, though the early evolution of the HAS is relatively easy to understand using quasi-static evolution of a system in hydrostatic and radiative equilibrium, the evolution of HAS beyond the time of galaxy and star formation remains challenging. Reionization may heat up the gas cloud surrounding the axion star.
Presence of gas with heavier atoms and molecules surrounding the axion star could cause it to accrete more gas or other massive particles. The reverse is also possible: a surrounding hot gas could cause the axion star to partially evaporate. In addition, the axion star is fairly fragile to collisions with an unrelated dense cloud of gas, which could strip the baryonic axion star cloud as the axion star core passes through undisturbed. A complete study of such astrophysical effects is beyond the scope of this work.
Despite the large number of open questions, we believe that signals due to hot matter surrounding axion stars deserve exploration using telescopes. We have found that axion stars could be quite nearby, within a few hundred astronomic units of the Earth. This leads to a point source signal that should be visible with high resolution telescopes. Due to the large uncertainty in the temperature of the axion stars, a wide frequency coverage is in order. We have found that the blackbody peak frequency may be in the range of ultraviolet to infrared if the baryonic cloud properties are frozen in at some early time. A "smoking gun" for ascertaining our signal would be HAS emissions from a low-baryon region of the sky, such as a dwarf spheroidal galaxy. Under certain assumptions, the total expected axion star luminosity in a dwarf galaxy could be comparable to or larger than the total observed luminosity, leading to constraints on the axion star properties.

Conclusions
The QCD axion remains a compelling solution to both the strong CP and dark matter puzzles. We have found that on small distance scales, the cosmology of the axion is interesting and leads to new potential signals. In particular, the following picture of QCD axion cosmology emerges: • Before the QCD phase transition, the axion field is uniform on roughly the Hubble scale, with a network of axion strings that are decaying into axions.
• Just before the QCD phase transition, the thermal axion mass becomes important and causes both the axion field to roll down to its minimum, producing axions, and the remaining long axion strings to annihilate away.
• After the QCD phase transition, the comoving density of axions is frozen, but there are large density contrasts on the scale of Hubble, around 10 km.
• The axions can rapidly relax down to their local ground state, forming a compact, dense, meter scale BEC, called axion stars. The typical mass of the axion star is roughly 10 −11 M ≈ 10 19 kg.
• Around the time of recombination, the primordial hydrogen gas begins to collapse into the gravitational potential well of the axion stars, reaching a dense, metallic fluid state and forming combined objects called Hydrogen Axion Stars (HAS).
• In this epoch, axions can interact with conduction band electrons, generating heat in the hydrogen fluid and causing it to blackbody radiate with a peak wavelength in the UV.
• As the universe expands, the gas surrounding the axion stars cools and further accretes onto the HAS, causing the hydrogen fluid to expand and cool. The peak wavelength of the radiation increases toward the IR. The mass of the hydrogen cloud is subdominant to that of the axions in the HAS, but can grow to be of order 10 13 kg, comparable to the mass of an asteroid, but distinguishable by its composition. The radius can grow to be several kilometers.
• At late times, when reionization and star formation occur, the properties of the HAS depend on its local surroundings and the fate of the HAS becomes more uncertain. The final state of the HAS depends on the redshift at which the HAS evolution stops.
The HAS is likely to be fairly hot and to have sufficient power output to be observed by telescopes sensitive to photons in the UV to IR range, yielding observational signals. In our local vicinity, axion stars making up a predominant component of the dark matter would have of a number density of order (300 AU) −3 , such that there can be nearby HAS visible as point sources to high resolution telescopes.
In addition, the combined HAS luminosity from faint dwarf spheroidal galaxies such as Ursa Major II (with order 10 17 HAS) can be comparable to the observed luminosity of such objects. A rich new phenomenology is open for exploration and discovery.

Acknowledgments
We would like to thank Eric Braaten

A Absorption of Axions in Metals
In this appendix, we determine the absorption cross-section for photons in a material and relate this to the absorption cross-section for an axion in the same material.
The absorption cross-section of photons in the material is related to the macroscopic electromagnetic properties of the material, which can be parametrized in terms of the complex dielectric constant, index of refraction, or conductivity (see Ref. [102][103][104] for detail). We choose to work with the latter for reasons that will become clear soon.
We need to obtain a relation between the conductivityσ = σ 1 + iσ 2 and the absorption crosssection σ abs . The approach we use is the optical theorem. For a particle P with four-momentum q (not necessarily on-shell), the optical theorem says that Im M[P (q) → P (q)] = q 2 Γ(q) , where Γ is the decay rate of P . The amplitude M[P (q) → P (q)] equals to the product of the polarization tensor Π µν (q) = e 2 d 4 x e iq·x 0|T [J µ † (x)J ν (0)]|0 and the corresponding polarization vectors for the state. 4 We can therefore write the polarization tensor as where L,T are longitudinal and transverse polarization four vectors respectively, satisfying q · L,T = 0.
Thus, if we determine the polarization tensor for the photon in the medium, we can determine the rate Γ at which photons get absorbed by the medium.
To determine the polarization tensor, we begin by determining a similar quantity, R µν (x, y), the kernel that transforms A into J − e J µ (x) = − d 4 yR µν (x, y) A ν (y) .
In momentum space, this relation can be written as since R µν (x, y) can be written as a function of x − y. The kernel can be determined by explicitly calculating the expectation value of the current operator J as where |Ω(t) is the time-evolved vacuum state, |0 is the vacuum state at t → ±∞, and U is the usual interaction picture evolution operator Since U is unitary and the vacuum state at t → ±∞ is the same, we can write U (∞, t) = U † (t, −∞), so that at leading order in perturbation theory we find The interaction Lagrangian is H int = −e d 3 xJ ν (x)A ν (x) as usual and we can write t −∞ dt = dt θ(t − t ), so that R µν = −i e 2 0|[J µ (x), J ν (y)]|0 θ(x 0 − y 0 ) .
We can see that this is closely related to the polarization tensor, differing just by the operator ordering.
In momentum space, this corresponds to a difference in the prescription for dealing with the poles, such that Re R µν (q) = Re Π µν (q) , Im R µν (q) = sgn(q 0 ) Im Π µν (q) .
Since we are only interested in q 0 > 0 here, we can just determine R µν .
To do so, we first rewrite J µ and A µ making use of charge conservation and gauge invariance respectively. First, in momentum space, the continuity equation allows us to write The current density J, by the definition of conductivity, is given by −eJ =σE, whereσ = σ 1 + iσ 2 .
Furthermore, the electric field E = −i(q A 0 − q 0 A), so that the current four-vector is Since we are calculating a physical quantity ultimately, our final result should be gauge invariant.
Nevertheless, it is helpful to pick a gauge for our calculation. One convenient gauge is the Coulomb gauge with q · A = 0, so that −eJ 0 = −iσq 2 A 0 /q 0 , L · A = 0 and L · A = 0 L A 0 = |q|A 0 / q 2 . We can then go ahead and solve the various components of Eq. (43). We begin with the 0 component.
leading to For the spatial components, with our choice of gauge, Since we have λ λ λ = L L + as usual for a vector polarization and q · A = 0, we find leading to Assuming that q is real, we find the decay rates from Eq. (40) for longitudinal and transverse modes are given by Further, if we take the approximation that |q| q 0 , which is certainly true for the kinematics of axion absorption, we find that Γ γ abs. = Γ L = Γ T = σ 1 .
In terms of the absorption cross-section in the medium, which we assume is coming entirely from absorption by free electrons in the conductor, we find that Γ γ abs. = n e σ abs v rel , where n e is the number density of free electrons, σ abs is the absorption cross-section, and v rel is the relative velocity of the incoming particles.
Given the photon absorption rate in the material, the axion absorption rate is determined by a simple rescaling [80] Γ a abs. = 3 m 2 a 4παf 2 a Γ γ abs. , where α is the fine-structure constant. Eq. (58) thus suffices to determine the axion absorption rate in a metal once the conductivity is known.