Gravitational Lensing Signatures of Axion Dark Matter Minihalos in Highly Magnified Stars

TL;DR

This paper proposes a novel method to detect axion dark matter minihalos through their subtle gravitational lensing effects on highly magnified stars near caustics in galaxy clusters, revealing potential signatures of axion dark matter.

Contribution

It introduces a new observational approach to identify axion minihalos via photometric monitoring of caustic transiting stars, focusing on their collective lensing signatures in galaxy clusters.

Findings

Axion minihalos produce surface density fluctuations detectable through microlensing light curve irregularities.

Minihalos with masses around 10^{-15} to 10^{-8} solar masses can cause measurable lensing effects.

Most minihalos are likely stable over cosmic timescales, avoiding disruption.

Abstract

Axions are a viable candidate for Cold Dark Matter (CDM) which should generically form minihalos of sub-planetary masses from white-noise isocurvature density fluctuations if the Peccei-Quinn phase transition occurs after inflation. Despite being denser than the larger halos formed out of adiabatic fluctuations from inflation, axion minihalos have surface densities much smaller than the critical value required for gravitational lensing to produce multiple images or high magnification, and hence are practically undetectable as lenses in isolation. However, their lensing effect can be enhanced when superposed near critical curves of other lenses. We propose a method to detect them through photometric monitoring of recently discovered caustic transiting stars behind cluster lenses, under extreme magnification factors $\mu \gtrsim 10^3$--$10^4$ as the lensed stars cross microlensing caustics induced by intracluster stars. For masses of the first gravitationally collapsed minihalos in the range $\sim 10^{-15}$--$10^{-8}\,h^{-1}\,M_\odot$, we show that axion minihalos in galaxy clusters should collectively produce subtle surface density fluctuations of amplitude $\sim 10^{-4}$--$10^{-3}$ on projected length scales of $\sim 10$--$10^4\,$AU, which imprint irregularities in the microlensing light curves of caustic transiting stars. We estimate that, inside a cluster halo and over the age of the Universe, most of these minihalos are likely to avoid dynamic disruption by encounters with stars or other minihalos.

Equations

$$M_0=(4\pi /3){\left(\pi /k_0\right)}^3{\bar{\rho }}_{a0},$$

$$P_m(k,z)$$

$$k_J=(16\pi Ga{\bar{\rho }}_{a0}{)}^{1/4}{\left(m_a/\hslash \right)}^{1/2}=6.7\times 10^{10}{\mathrm{Mpc}}^{-1}(1+z{)}^{-1/4}{\left(\frac{{\Omega }_m{h}^2}{0.12}\right)}^{1/4}{\left(\frac{m_a}{10^{-4}\mathrm{eV}}\right)}^{1/2}.$$

$$1+z_{\mathrm{coll}}(M)={\left(18/5{\pi }^2\right)}^{1/2}\left({\nu }_{\mathrm{med}}/{\delta }_c\right)\left(1+z_{\mathrm{eq}}\right){\left(M_0/M\right)}^{1/2}=0.24\left(1+z_{\mathrm{eq}}\right){\left(M_0/M\right)}^{1/2}.$$

$${\rho }_s(M)\approx 0.24M_{\odot }{h}^2/{\mathrm{pc}}^3{\left(\frac{\nu }{{\nu }_{\mathrm{med}}}\right)}^3\frac{{\Omega }_a}{0.3}{\left(\frac{1+z_{\mathrm{eq}}}{3400}\right)}^3{\left(\frac{M_0}{10^{-10}{M}_{\odot }/h}\right)}^{3/2}{\left(\frac{M}{10^{-6}{M}_{\odot }/h}\right)}^{-3/2},$$

$$r_s(M)\approx 1530h^{-1}\mathrm{AU}{\left(\frac{1+z_{\mathrm{eq}}}{3400}\frac{\nu }{{\nu }_{\mathrm{med}}}\right)}^{-1}{\left(\frac{{\Omega }_a}{0.3}\right)}^{-1/3}{\left(\frac{M_0}{10^{-10}{M}_{\odot }/h}\right)}^{-1/2}{\left(\frac{M}{10^{-6}{M}_{\odot }/h}\right)}^{5/6}.$$

$${\kappa }_s(M)\approx 7\times 10^{-7}\frac{D_{\mathrm{eff}}}{1\mathrm{Gpc}}{\left(\frac{\nu }{{\nu }_{\mathrm{med}}}\right)}^2{\left(\frac{{\Omega }_a}{0.3}\right)}^{2/3}{\left(\frac{1+z_{\mathrm{eq}}}{3400}\right)}^2\left(\frac{M_0}{10^{-10}{M}_{\odot }/h}\right){\left(\frac{M}{10^{-6}{M}_{\odot }/h}\right)}^{-2/3}.$$

$$m_a>{\left[\frac{(2\pi \hslash {)}^2}{G{\rho }_s{r}_s^4}\right]}^{1/2}\approx 8\times 10^{-12}\mathrm{eV}h{\left(\frac{{\rho }_s}{10^4{M}_{\odot }{h}^2{\mathrm{pc}}^{-3}}\right)}^{-1/2}{\left(\frac{r_s}{10^2{h}^{-1}\mathrm{AU}}\right)}^{-2}.$$

$$\frac{\mathrm{d}f(M,t)}{\mathrm{d}\ln M}=\sqrt{\frac{2}{\pi }}\nu (M,t)\exp \left(-\frac{{\nu }^2(M,t)}{2}\right)\left|\frac{\partial \ln \sigma (M,t)}{\partial \ln M}\right|.$$

$$P_{\rho }(q)=\int \mathrm{d}\ln Mn(M){\left|{\tilde{\rho }}^{(h)}\left(q;M\right)\right|}^2=\overline{\rho }\int \frac{\mathrm{d}M}{M^2}\frac{\mathrm{d}f(M)}{\mathrm{d}\ln M}{\left|{\tilde{\rho }}^{(h)}\left(q;M\right)\right|}^2.$$

$${\tilde{\rho }}^{(h)}(q;M):=4\pi {\int }_0^{\infty }{r}^2\mathrm{d}r\frac{\sin \left(qr\right)}{qr}{\rho }^{(h)}\left(r;M\right).$$

$$P_{\Sigma }(q_{\perp })$$

$$P_{{\delta }_{2d}}(q_{\perp })=\frac{1}{\overline{\Sigma }}\left(\int \frac{\mathrm{d}M}{M^2}\frac{\mathrm{d}f(M)}{\mathrm{d}\ln M}{\left|{\tilde{\rho }}^{(h)}\left(q_{\perp };M\right)\right|}^2\right).$$

$$P_{\kappa }(q_{\perp })=\frac{{\overline{\Sigma }}_{\mathrm{cl}}}{{\Sigma }_c^2}\int \frac{\mathrm{d}M}{M^2}\frac{\mathrm{d}f(M)}{\mathrm{d}\ln M}{\left|{\tilde{\rho }}^{(h)}\left(q_{\perp };M\right)\right|}^2.$$

$$\frac{\mathrm{d}\tau (M)}{\mathrm{d}\ln M}=\frac{\pi r_s^2(M){\overline{\Sigma }}_{\mathrm{cl}}}{M}\frac{\mathrm{d}f(M)}{\mathrm{d}\ln M}.$$

$$t_{\mathrm{df}}\approx t_{\mathrm{orb}}\left(\frac{M}{m_1}\right){\left(\ln \frac{M}{m_1}\right)}^{-1},$$

$$\frac{\Delta E_1}{E_1}\sim \frac{{\sigma }_1^2}{{\sigma }^2}.$$

$${\Gamma }_{\mathrm{enc}}(m_1;M)\sim \frac{1}{m_1}\frac{\mathrm{d}{f}_{\mathrm{sub}}(m_1;M)}{\mathrm{d}\ln m_1}\frac{R}{t_{\mathrm{orb}}(M)}\frac{r_1^{2}M}{R^3},$$

$$t_{\mathrm{dispt}}(m_1;M)\sim t_{\mathrm{orb}}(M)\frac{R}{r_1}{\left(\frac{\mathrm{d}{f}_{\mathrm{sub}}(m_1;M)}{\mathrm{d}\ln m_1}\right)}^{-1},$$

$$t_{\mathrm{dispt}}(m_1;M)\sim t_{\mathrm{df}}(m_1;M){\left(\frac{{\sigma }_1}{\sigma }\right)}^2{\left(\frac{\mathrm{d}{f}_{\mathrm{sub}}(m_1;M)}{\mathrm{d}\ln m_1}\right)}^{-1}\ln \frac{M}{m_1}.$$

$$t_{\mathrm{dispt}}$$

$$\frac{\Delta v}{{\sigma }_1}\simeq \frac{G{m}_{\star }{r}_1}{b^2\sigma {\sigma }_1}.$$

$$\frac{b_{\min }}{r_1}$$

$$t_{\mathrm{dispt},\star }$$

$${\mu }_{\mathrm{pk}}$$

$$|\tilde{\mathbf{d}}|\sim {\kappa }_{\star }^{3/2}/{\theta }_{\star }\sim 350{\mathrm{arcsec}}^{-1}{\left(\frac{{\kappa }_{\star }}{0.01}\right)}^{3/2}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{1/2}{\left(\frac{D_{\mathrm{eff}}}{1\mathrm{Gpc}}\right)}^{-1/2}.$$

$$a(t_0)=\left(\frac{T_{\gamma 0}}{T_0}\right){\left[\left(2+6\frac{7}{8}{\left(\frac{T_{\nu 0}}{T_{\gamma 0}}\right)}^3\right)\frac{1}{g_{\ast s}}\right]}^{1/3},$$

$$M_0=10^{-9}{M}_{\odot }{h}^2\left(\frac{{\Omega }_a}{0.27}\right){\left(\frac{T_0}{1.3\mathrm{GeV}}\right)}^{-3}\left(\frac{g_{\ast s}}{78}\right){\left(\frac{g_{\ast ϵ}}{78}\right)}^{-3/2}.$$

$$h^{\prime \prime }+\frac{2+3a}{2a(1+a)}{h}^{\prime }$$

$${\delta }_m^{\prime }$$

$${\delta }_r^{\prime }$$

Full text

cosmology, gravitational lensing, dark matter

1 Introduction

Gravitational lensing can reveal self-gravitating objects containing only dark matter (DM), which are difficult to detect by other means because they have expelled or have never accreted baryons. Most DM theories predict the presence of these objects on mass scales much smaller than that of galaxies. Direct detection of these structures would provide crucial clues to the physical origin of DM. A common problem is that these structures are predicted to have surface densities well below the critical value, so they by themselves cannot act as gravitational lenses capable of creating multiple images or inducing substantial image magnifications and distortions, and their effects are often unnoticeable.

However, when these small-scale objects are superposed on a large-scale lens near lensing critical curves, very small surface density perturbations may result in large observable changes (Minor et al., 2017; Dai et al., 2018b). Recently, the first highly magnified stars (Miralda-Escude, 1991) have been discovered in the field of massive clusters of galaxies (Kelly et al., 2018; Chen et al., 2019; Kaurov et al., 2019). Near the critical curve of the cluster, individual superluminous stars in the background become visible to our most powerful telescopes when magnified by a factor $\sim 10^{2}$–$10^{3}$. Microlensing by intracluster stars in the lensing cluster introduces fast variability in the magnification (Venumadhav et al., 2017; Diego et al., 2018; Oguri et al., 2018), turning highly magnified stars into a probe not only of the abundance and mass function of intracluster stellar objects, but also of any small-scale density inhomogeneity in the DM.

Axions are a promising DM candidate produced by non-thermal mechanisms in the early Universe. For a famous example, the QCD axion is motivated by the Peccei-Quinn mechanism (Peccei & Quinn, 1977), which explains the absence of a CP violating term in the Quantum Chromodynamics (QCD) lagrangian (Weinberg, 1978; Wilczek, 1978; Kim, 1979; Shifman et al., 1980; Dine et al., 1981) as implied by the null value of the neutron electric dipole moment (Pendlebury et al., 2015). Generic axions can arise from other spontaneous broken U(1) symmetries (Chikashige et al., 1981; Gelmini & Roncadelli, 1981; Wilczek, 1982; Berezhiani & Khlopov, 1990; Jaeckel, 2014) or from low-energy effective theories emerging from string theory (Witten, 1984; Conlon, 2006; Cicoli et al., 2012; Georgi et al., 1981; Choi et al., 2009; Dias et al., 2014), and may also account for the DM.

Cosmological models involving axions as the DM predict the existence of very dense, low-mass dark matter halos through gravitational instability (Hogan & Rees, 1988; Kolb & Tkachev, 1994) if the Peccei-Quinn symmetry breaking occurs after inflation. However, the surface density of these minihalos would be too low to make them generally detectable through gravitational lensing. Despite a number of proposed lensing methods to detect small-scale axion structures (Kolb & Tkachev, 1996; Fairbairn et al., 2017, 2018), so far the range of axion minihalo densities and masses do not seem accessible with any foreseeable observational techniques.

In this paper we propose the first realistic astrophysical method to detect axion minihalos through gravitational lensing. Highly magnified stars lensed by clusters of galaxies can have apparent fluxes affected by the presence of axion minihalos when they are crossing one of the micro-caustics produced by intracluster stars. While the intracluster stars disrupt the cluster critical curve into a network of micro-critical curves on angular scales of $\sim 10\,\mu{\rm as}$, axion minihalos having sub-planetary masses and solar system sizes would produce finer magnification variations on much smaller scales of $\sim 10^{2}$–$10^{3}$ nano-arcseconds. When crossing a micro-caustic, the observed light curve of a highly magnified star should be altered by surface density fluctuations caused by axion minihalos inside the cluster lens. This paper aims to predict the lensing impacts of axion minihalos on highly magnified stars, which can be monitored through dedicated observations with our most powerful space-borne or ground-based optical and infrared telescopes (Windhorst et al., 2018; Diego, 2019).

We outline this work as the following. In Section 2, we start with a review on primordial isocurvature density fluctuations expected to arise from a Peccei-Quinn phase transition in axion DM scenarios. Then in Section 3, we derive the physical properties of axion minihalos that collapse from primordial isocurvature density fluctuations, and show that they are significantly denser and more compact than halos of comparable masses expected in the standard cosmology. In Section 4, we discuss the mass function of axion minihalos. In Section 5, we calculate the size of surface mass density fluctuations due to the cumulative effect of many substructure minihalos along the line of sight through the cluster halo. In Section 6, we estimate possible dynamic disruption of axion minihalos by mutual encounters or encounters with stars. In Section 7, we study the lensing signatures of axion minihalos imprinted on highly magnified stars, and show that sizable irregularities are expected in their light curves during micro-caustic crossings. Additional discussion is presented in Section 8, and we conclude in Section 9. We discuss some technical details in a few Appendices. We use the conventional notation in cosmology that many physical quantities are expressed in a way that scales with the dimensionless Hubble parameter $h$.

2 Isocurvature density fluctuations in axion cosmology

Axions are a viable particle candidate for non-thermal DM, namely DM that never interacted substantially with the baryonic matter in the Universe since it was created. The axion DM hypothesis is that a global U(1) symmetry of a complex scalar field $\varphi$ was spontaneously broken in the very early Universe (Peccei & Quinn, 1977) and gave rise to the axion particle $\phi$ as the Goldstone mode associated with the angle of $\varphi$. If the global U(1) symmetry breaking occurs after inflation during the era of radiation domination, causally disconnected Hubble patches settle down to different vacuum expectation values of $\phi$.

The axion can acquire a non-zero mass $m_{a}$ through non-perturbative effects such as instantons (Gross et al., 1981) which can tilt the potential so that it is no longer flat with respect to $\phi$. The axion mass is in general temperature dependent, which induces fluctuations of order unity in the axion field energy density as the field configuration oscillates around the potential minimum $\phi=0$. This, however, does not occur immediately after the axion acquires its mass because the Hubble timescale is shorter than $\hbar/(m_{a}\,c^{2})$ in the beginning and the axion is dynamically “frozen”. During this time, the field configuration is smoothed over the instantaneous horizon scale by the Kibble mechanism (Kibble, 1976).

The situation changes when the cosmic expansion eventually slows down to be comparable to the axion field oscillation time scale, $3\,H(t_{0})\simeq m_{a}(t_{0})\,c^{2}/\hbar$. The axion field, which is coherent on the horizon scale at this epoch, $\lambda_{H}(t_{0})\sim c/H(t_{0})$, starts to oscillate around $\phi=0$ and to contribute to the cosmic matter density budget. This generates isocurvature perturbations in addition to the inflationary adiabatic perturbations. The QCD axion starts to oscillate when the Universe cools to a temperature $T_{0}=T(t_{0})\simeq 1.3\,{\rm GeV}$ (Buschmann et al., 2019), and the axion mass evolves with the temperature from $m_{a}(T_{0})\sim 6\times 10^{-9}\,{\rm eV}$ to $m_{a}=2.5\times 10^{-5}\,{\rm eV}$ at zero temperature (di Cortona et al., 2016; Buschmann et al., 2019), although the exact numbers are subject to theoretical uncertainties (Borsányi et al., 2016). A generic axion particle can have a different mass for which the onset of field oscillation is at a different epoch.

Assuming Gaussian statistics, the total power spectrum for the matter overdensity is the sum of the standard inflationary adiabatic power spectrum $P_{\rm cdm}$ and the axion isocurvature power spectrum $P_{\rm iso}$ 111Throughout the paper we adopt the simplifying terminology “isocurvature fluctuation” to mean the perturbation in the ratio of the axion density to total density, with a constant baryon-to-photon ratio, which should be distinguished from other forms of isocurvature modes.. At the initial time $t_{0}$, the latter can be modeled as a white-noise power spectrum (Hogan & Rees, 1988; Fairbairn et al., 2018) $P_{\rm iso}(k)=\Theta(k_{0}-k)\,24\,\pi^{2}/(5\,k^{3}_{0})$, where $k_{0}=a(t_{0})\,H(t_{0})/c$ is a cutoff comoving wave number associated with the horizon scale at $t_{0}$. Requiring a DM overdensity variance of order unity on the scale $k_{0}$ at $t_{0}$ fixes the normalization of $P_{\rm iso}$. The total amount of axion mass enclosed within a spherical volume of a comoving radius $\pi/k_{0}$ defines a characteristic mass scale

[TABLE]

where $\bar{\rho}_{a0}$ is the present day mean density in axions. This sets the mass scale of the first collapsed objects from the isocurvature perturbations. For the QCD axion, we obtain $M_{0}\simeq 5\times 10^{-10}\,M_{\odot}$ based on the axion mass parameterization used in Buschmann et al. (2019) and the the effective number of relativistic degrees of freedom presented by Husdal (2016) (see Appendix A). After accounting for a factor of $4\,\pi^{4}/3\approx 130$ larger in our definition, this number is compatible with the characteristic masses for the first gravitationally collapsed axion minihalos quoted in other references (Davidson & Schwetz, 2016; Hardy, 2017; Fairbairn et al., 2018). We therefore estimate that the typical mass for the first collapsed axion minihalos is $\sim 0.01\,M_{0}$, and is $\sim 5\times 10^{-12}\,M_{\odot}$ for QCD axions. For generic axion DM, a wide range of values for $M_{0}$ are possible.

The white noise model is not likely to be a good approximation around the cutoff scale $k_{0}$, where the axion field fluctuations are subject to complicated nonlinear dynamics and can be significantly non-Gaussian (Kolb & Tkachev, 1996). Numerical simulations find intricate axion clumps on scales shorter than $k_{0}^{-1}$ (Vaquero et al., 2018; Buschmann et al., 2019). Unrelated to gravitational instability, these clumps emerge soon after the onset of field oscillation at $t_{0}$ as a result of complicated Klein-Gordon dynamics, with mass scales as much as four orders of magnitude smaller than the analytic estimate $\sim 0.01\,M_{0}$ for the smallest gravitationally collapsed minihalos (Vaquero et al., 2018; Buschmann et al., 2019). These clumps are likely to be contained within the first gravitationally collapsed minihalos, and we distinguish them from axion minihalos, which are the focus of study throughout this paper. The Gaussian white-noise model for minihalos should work well in the shot-noise regime $k\ll k_{0}$, regardless of the formation of smaller clumps.

With favorable parameter values for the axion cosmology, $k_{0}$ typically translates into a length scale many orders of magnitude smaller than the range of scales probed through the Cosmic Microwave Background anisotropies and the large-scale clustering of galaxies. Therefore, the isocurvature modes that dominate over the adiabatic modes enter the horizon well before the epoch of radiation-matter equality at $z_{\rm eq}\simeq 3400$, so they have wave numbers $k\gg k_{\rm eq}=a(z_{\rm eq})\,H(z_{\rm eq})/c$. Analytic and numerical calculations of linear perturbations have shown that sub-horizon isocurvature matter modes start to grow substantially only when $z>z_{\rm eq}$ (Efstathiou & Bond, 1986). Adiabatic modes behave differently during the era of radiation domination: upon entering the horizon, the oscillation of the photon-baryon fluid as an acoustic wave generates an initial peculiar velocity in the DM, leading to a comoving displacement of the DM that grows logarithmically with the scale factor and is responsible for the shape of the linear CDM power spectrum at $k>k_{\rm eq}$. In comparison, isocurvature modes induce no such initial peculiar velocities, and grow only through the self gravitation of the DM after horizon entry and the photon-baryon acoustic oscillation (see Appendix B for more details).

Since the white noise isocurvature power spectrum is only meant to be a crude model, we adopt the approximation that all the isocurvature modes grow by the same factor, equal to the usual $\Lambda$CDM growth factor $D_{+}(z)$. We neglect the growth of isocurvature modes at $z<z_{\rm eq}$, noting that this small growth is in any case a nearly constant factor over the relevant range of scales that can be reabsorbed into the definition of $k_{0}$. We can then write the linear power spectrum for the total matter overdensity (including both isocurvature and adiabatic modes, but neglecting the gravitational influence of baryons on the isocurvature component) at redshift $z$ as

[TABLE]

This is valid for $z<z_{\rm eq}$, where $z_{i}$ is an arbitrarily chosen reference redshift during the era of matter domination satisfying $1\ll z_{i}\ll z_{\rm eq}$.

We note that any primordial adiabatic fluctuations of large amplitude and on small scales would be subject to stringent constraints on entropy production from acoustic damping in the primordial plasma (Jeong et al., 2014; Inomata et al., 2016). Since modes of $k\sim 10^{2}$–$10^{8}\,h\,{\rm Mpc}^{-1}$ enter the horizon well before the epoch of radiation-matter equality, they would have rapidly dissipated. By contrast, large CDM isocurvature fluctuations are allowed because density fluctuations during the radiation era $z\gg z_{\rm eq}$ are highly suppressed. Neither are these axion isocurvature modes subject to Big Bang Nucleosynthesis constraints, which apply to baryonic isocurvature modes of $k\lesssim 10^{9}\,h\,{\rm Mpc}^{-1}$ (Inomata et al., 2018).

In Figure 1, we plot the total matter overdensity power spectrum linearly extrapolated to $z=0$ for $M_{0}$ in the range $10^{-10}$–$10^{-6}\,M_{\odot}/h$. On comoving scales shorter than the inverse of $k\sim 10^{4}$–$10^{8}\,h/{\rm Mpc}$, the isocurvature contribution can dominate the total power. This means that isocurvature perturbations lead to earlier gravitational collapse on those small scales. In fact, on scales $k\sim k_{0}$, fluctuations are initially of order unity and they reach non-linear collapse near the epoch of radiation-matter equality. Note that the sharp reduction of power at scale $k_{0}$ in Figure 1 is expected to be smoothed once axion dynamics on scales shorter than $k_{0}$ are accurately calculated.

On sufficiently small scales, gravitational instability is prohibited by the quantum degeneracy pressure of the axion field. The comoving wave number of the effective Jeans scale is given by (Fairbairn et al., 2018)

[TABLE]

Structure formation becomes efficient at $z<z_{\rm eq}\simeq 3400$, and the expected QCD axion has a mass in the range $m_{a}\sim 10^{-6}$–$10^{-4}\,$eV at zero temperature, so Eq. (3) suggests that quantum degeneracy pressure is negligible for the gravitational instability of the isocurvature modes shown in Figure 1.

3 Size and density of axion minihalos

During the epoch of matter domination, collapsed objects form after the matter overdensity amplitude reaches a value near unity. Around the epoch of radiation-matter equality, the first collapsed objects form on a comoving scale $\sim 1/k_{0}$ with masses $\sim 0.01\,M_{0}$. Similar to that in the standard cosmology, hierarchical assembly of collapsed objects proceeds in a bottom-up fashion toward larger length scales (and smaller comoving wave numbers). Increasingly larger halos are built up from mergers of smaller ones and accretion of smooth matter, and they develop an internal structure of orbiting satellite halos, with at least a fraction of them surviving dynamical disruption. This structure may have multiple levels of orbiting satellites within orbiting satellites. The hierarchical process takes place first during $z_{\rm eq}>z\gtrsim 20$ on very small scales, where the white-noise isocurvature power spectrum dominates. Then, at $z\lesssim 20$, the inflationary adiabatic power spectrum is dominant, and structure formation proceeds in a way identical to that of the standard CDM cosmology, except that a fraction of the DM is clumped in the form of vast numbers of very small and dense orbiting minihalos.

Consider the collapse of a spherical region at redshift $z_{\rm coll}$, which has a comoving Lagrangian radius $R$ and an enclosed mass $M=(4\pi/3)\,\bar{\rho}_{a0}\,R^{3}$. Exact spherical symmetry in collapse is unrealistic for predicting the correct density profile of the collapsed halos. Numerical studies of collapse from Gaussian random overdensities using N-body simulations suggest that halos formed through hierarichical assembly in general follow the phenomenological Navarro-Frenk-White (NFW) density profile (Navarro et al., 1996, 1997) $\rho(r)=\rho_{s}/[(r/r_{s})(1+r/r_{s})^{2}]$, which is defined by a scale density $\rho_{s}$ and a scale radius $r_{s}$.

The literature on CDM halos has commonly used three quantities to characterize an NFW halo: the virial mass $M_{200}$, the virial radius $r_{200}$, and the concentration parameter $c_{200}=r_{200}/r_{s}$. Instead, we treat $\rho_{s}$ and $r_{s}$ as fundamental parameters which are closely tied to the time of halo collapse. Numerical studies suggest a “two step” picture for the formation and evolution of halos (especially at low peak heights): upon collapse, the halo first undergoes a phase of rapid growth to reach a typical concentration factor $c_{\star}\simeq 4$ (Zhao et al., 2003, 2009); after that, halo growth slows down while its inner profile stabilizes (Bullock et al., 2001; Ludlow et al., 2013). During the second phase, which is often referred to as pseudo evolution, $\rho_{s}$ and $r_{s}$ hardly change if measured in proper units, while the virial mass and the virial radius still increase due to mass accretion onto the outer parts. The picture of pseudo evolution applies particularly well to isolated halos that grow slowly from the low-density surrounding material, which are typically of low mass compared to the characteristic collapsing mass at each epoch. Using $\rho_{s}$ and $r_{s}$ is convenient for quantifying the gravitational lensing effect; they characterize the halo’s inner region, which is much denser and more compact than the entire virialized region at high halo concentration. We therefore describe NFW halos using $\rho_{s}$ and $r_{s}$, assumed to be stationary from the onset of the pseudo-evolution phase.

The linear matter overdensity smoothed over a scale $R$ has a variance $\sigma^{2}(M,z)=\int\,\mathrm{d}\ln k\,$ $[k^{3}\,P_{m}(k,z)]/(2\pi^{2})\,|W(k,R)|^{2}$, where $W(k,R)=3\,(kR)^{-3}\,[\sin(kR)-(kR)\,\cos(kR)]$ is the spherical top-hat window function. At each redshift, one can derive the mass scale $M$ and the comoving radius scale $R=R(M)$ that collapse from any given peak height $\nu(M,z)=\delta_{c}/\sigma(M,z)$, where $\delta_{c}=1.686$ is the threshold overdensity during the era of matter domination. Halos of mass $M$ collapse over a broad range of redshifts, but we can define a universal median collapse redshift corresponding to a fluctuation with $\nu=\nu_{\rm med}=0.67$. For any $M$, a characteristic mass-dependent collapse redshift $z_{\rm coll}(M)$ follows by solving $\nu_{\rm med}=\delta_{c}/\sigma(M,z_{\rm coll})$. When isocurvature fluctuations dominate the overdensity power spectrum in Eq. (2), the median collapse redshift is derived to be

[TABLE]

To estimate the NFW fundamental parameters, we assume that the collapsed NFW halo has a universal concentration parameter $c_{\star}=r_{200}/r_{s}=4$, its enclosed mass within $r_{200}$ equals $M$, and the mean density within $r_{200}$ is 200 times the cosmic mean $\bar{\rho}_{a}(z_{\rm coll})=\bar{\rho}_{a0}\,(1+z_{\rm coll})^{3}$. For $c_{\star}=4$, we find $\rho_{s}\approx 5271\,\bar{\rho}_{a}(z_{\rm coll})$, $r_{s}\approx 0.0265\,[M/\bar{\rho}_{a}(z_{\rm coll})]^{1/3}$, and $M\approx 10.2\,\rho_{s}\,r^{3}_{s}$. When the white-noise isocurvature power spectrum dominates, the resultant halos, which we refer to as minihalos, obey the following scaling relations with the mass scale $M$,

[TABLE]

and

[TABLE]

It is useful at this point to compare the characteristic surface density of axion minihalos, $\rho_{s}(M)\,r_{s}(M)$, with the critical surface density in gravitational lensing, required for an isolated lens to produce multiple images and high magnification. The critical surface density is given by $\Sigma_{\rm crit}=c^{2}/(4\pi\,G\,D_{\rm eff})$, where $D_{\rm eff}=D_{L}\,D_{LS}/D_{S}$ for angular diameter distances $D_{L}$ to the lens, $D_{S}$ to the source, and $D_{LS}$ from the lens to the source. The typical lensing convergence of an axion minihalo near the scale radius, $\kappa_{s}(M)\approx\rho_{s}(M)\,r_{s}(M)/\Sigma_{\rm crit}$, is tiny,

[TABLE]

Even for dense minihalos with $M=10^{-10}\,M_{\odot}/h$, acting as individual lenses at cosmological distances, $\kappa_{s}(M)$ is no greater than $\sim 10^{-3}$. Axion minihalos can therefore not produce substantial lensing effects by themselves, but we will show in this paper that they can be detected when acting in conjunction with a galaxy cluster lens and intracluster stars near a lensing critical curve.

In Figure 2, we check whether the above simple peak height prescription gives reasonably good estimates for the fundamental parameters $\rho_{s}$ and $r_{s}$ of halos at low masses. We first generate random NFW halos in the standard $\Lambda$CDM cosmology within a range of redshifts $0<z<15$ and a range of virial masses $-16<\log(M_{200}\,h/M_{\odot})<10$ sampled from N-body simulations (Diemer & Joyce, 2018). Next, we calculate their $r_{s}$ and $\rho_{s}$ using a number of median concentration-mass relations from the literature, which are empirically calibrated to simulations. When we scatter plot the halos on the $(r_{s},\,\rho_{s})$ plane, they all lie within a squeezed band, which is particularly narrow at very small halo masses. This indicates that the notion of pseudo evolution is reasonably valid for low mass halos. For a comparison, we then analytically estimate $\rho_{s}$ and $r_{s}$ of median halos using the peak height prescription: we identify peaks with a median height $\nu=\nu_{\rm med}$ at various redshifts, and then calculate using the relations $\rho_{s}\approx 5271\,\bar{\rho}_{a}(z_{\rm coll})$ and $r_{s}\approx 0.0265\,[M/\bar{\rho}_{a}(z_{\rm coll})]^{1/3}$, where at any given collapse redshift $z_{\rm coll}$ we solve the top-hat mass $M$ from $\nu(M,z_{\rm coll})=\delta_{c}/\sigma(M,z_{\rm coll})=\nu_{\rm med}$. On the $(r_{s},\,\rho_{s})$ plane, this produces a single curve parameterized by the collapse redshift $z_{\rm coll}$ or alternatively the top hat mass $M$. At small halo masses, this curve aligns reasonably well with the narrow band of halos identified from numerical simulations.

With this justification, hereafter we rely on the above peak height prescription to estimate halo size and density in axion cosmology. We show some sample curves in Figure 2. We see that in axion cosmology, sufficiently small halos can be many orders of magnitude denser than their $\Lambda$CDM counterparts (Diemand et al., 2005; Berezinsky et al., 2006; Diemand et al., 2006).

The inner region of the collapse axion halo has a typical de Broglie wavelength $\lambda_{\rm dB}=(2\,\pi\,\hbar)/(m_{a}\,\sigma)$, where we estimate velocity dispersion as $\sigma\sim(G\,\rho_{s}\,r_{s})^{1/2}$. The quantum degeneracy pressure is dynamically unimportant if $\lambda_{\rm dB}\ll r_{s}$, which requires

[TABLE]

Whenever this is valid, the axion halo can be modeled as a self-gravitating system made of classical collisionless particles (if axion self-interaction is negligible). In Figure 2, we overplot contours for this minimum axion mass, such that for any given value of $m_{a}$ halos located to the upper right side of the contour on the $(r_{s},\,\rho_{s})$ plane are not dynamically affected by the degeneracy pressure. For the mass range $m_{a}\sim 10^{-6}$–$10^{-4}\,$eV (at zero temperature) which brackets the theoretically favorable mass for the QCD axion, axion minihalos formed from isocurvature fluctuations are not affected by the degeneracy pressure.

4 Mass function of axion minihalos

The Press-Schechter (PS) formalism (Press & Schechter, 1974) predicts that at a given time $t$ the collapsed halos have a mass function $n(M,t)=(\bar{\rho}_{m}(t)/M^{3})\,\mathrm{d}f(M,t)/\mathrm{d}\ln M$. We define $n(M,t)$ to be the proper volume number density of halos per logarithmic interval of mass $M$. The mass fraction per logarithmic interval of mass is

[TABLE]

Here $\bar{\rho}_{m}(t)$ is the mean matter density, and $\nu(M,t)=\delta_{c}/\sigma(M,t)$ is the peak height.

The PS mass function excludes small halos that have been incorporated into larger ones as subhalos. Some subhalos are disrupted, while others may survive for long. Both isolated halos and subhalos act as perturber lenses, so as far as gravitational lensing is concerned, Eq. (9) sets a lower limit to the fraction of mass locked up in collapsed, dense objects.

In Figure 3, we adopt the PS formalism to compute $\mathrm{d}f/\mathrm{d}\log M$ at various redshifts. We consider different values of $M_{0}$ for the axion cosmology and compare the results to those of the standard cosmology. As an example, for $M_{0}=10^{-10}\,h^{-1}\,M_{\odot}$, by $z=100$ virtually all DM is locked up in collapsed axion minihalos of masses $M<10^{-6}\,h^{-1}\,M_{\odot}$, well before any galactic scale halos form. As hierarchical assembly proceeds, small minihalos gradually assemble into larger halos. By $z=10$, about $50\%$ of DM resides in isolated minihalos of masses $M<10^{-4}\,h^{-1}\,M_{\odot}$, and a substantial fraction has already assembled into the first dwarf-galaxy halos from the nonlinear collapse of adiabatic fluctuations. By $z=1$, around the time when the most massive galaxy clusters are forming, isolated minihalos of masses $M<10^{-4}\,h^{-1}\,M_{\odot}$ still contribute about $10\%$ of all DM. This provides a lower limit on the mass budget of minihalos of $M<10^{-4}\,h^{-1}\,M_{\odot}$ in the Universe, because we have excluded subhalos of similar masses whose contribution depends on the substructure survival fraction and is not accounted for by the PS formalism.

Figure 3 shows that structure formation in axion cosmology differs substantially from that in the standard cosmology because the collapse times of the very low-mass halos shift to earlier epochs. In the standard cosmology, practically no DM structures have collapsed on any scales by $z=100$. The lowest mass halos form only at $z\sim 20$, at an epoch not well separated from the formation times of the first halos of dwarf galaxy scale. For the axion cosmology with $M_{0}=10^{-10}\,h^{-1}\,M_{\odot}$, isolated halos with $M\lesssim 10^{-4}\,h^{-1}\,M_{\odot}$ account for $\sim 30\%$ of all DM. By $z=1$, this fraction drops to $\sim 10\%$ as more of them merge into larger halos. On the other hand, the DM mass fraction (including substructure) residing in large halos converges between the standard cosmology and the axion cosmology. The PS calculation indicates that in axion cosmology, the mass fraction in tiny isolated minihalos with $M<10^{-4}\,h^{-1}\,M_{\odot}$ is much higher than that in standard cosmology. Moreover, the low-mass axion minihalos are much denser, so they are much more likely to survive as substructure when they merge into larger halos. When the value of $M_{0}$ is increased, for example, to $M_{0}=10^{-6}\,h^{-1}\,M_{\odot}$, the characteristic minihalo masses shift to larger masses; the aforementioned conclusions remain qualitatively true after numbers are appropriately rescaled.

The PS calculation includes only isolated minihalos, whereas observations of caustic transiting stars behind lensing clusters are predominantly probing DM subhalos inside galaxy clusters rather than intervening field halos in intergalactic space (c.f. discussion in Section 5). Nevertheless, we make the following argument to translate the above results to the subhalo mass fraction inside a galaxy cluster, up to caveats regarding minihalos surviving dynamic disruption in the intracluster environment. Most of the mass of a galaxy cluster accretes at $z\sim 1$ from the intergalactic medium, where the abundance of small DM halos shortly before the accretion to the cluster is well described by the PS mass function. Being much denser than the cluster halo, the axion minihalos are unlikely to be disrupted once they merge into the cluster. The fractional mass function of minihalos inside the cluster halo can therefore be approximately described by Eq. (9). This approximation is still a lower limit to the number of minihalos, because the PS formalism does not include minihalos that collapse and survive as subhalos in intermediate mass halos, and are incorporated into the cluster when the intermediate mass object merges into the cluster.

Using this justification, we shall use Eq. (9) to estimate the power spectrum of the surface density projected along the line of sight through a cluster halo in Section 5. The various theoretical uncertainties regarding dynamical disruption of minihalos and the multiple levels of subhalos inside subhalos will affect our result only by a moderate factor reflecting the fraction of DM in axion minihalos.

5 Cumulative effect of minihalos along a line of sight

The main goal of this paper is to analyse the gravitational lensing impact of axion minihalos in the line of sight to a source observed at high magnification. As mentioned, axion minihalos are difficult to detect with gravitational lensing because of their low surface density. However, small surface density fluctuations $\Delta\Sigma$ produce large effects near critical curves, when the lensing magnification is as high as $\sim\Sigma_{\rm crit}/\Delta\Sigma$.

Near the critical curve of a galaxy cluster, the total cluster surface density is close to $\Sigma_{\rm crit}$, whereas the surface density of an individual minihalo is much lower. This implies that if, as argued above, axion minihalos contain a substantial fraction of all the DM, their area covering factor is much larger than unity. Therefore, axion minihalos must move through each other on many high-speed encounters as they orbit inside the cluster halo. At the same time, surface density fluctuations near a lensing critical curve result from the superposition of many axion minihalos, and can therefore be treated as a Gaussian random field, with all the information being included in the power spectrum. We evaluate the surface density power spectrum in this section.

5.1 Surface density power spectrum

A useful starting point of the calculation is to consider that all the relevant density fluctuations are contained in a statistically homogeneous slab of thickness $L$, which is much larger than any transverse scales of interest. As shown in Appendix C, the power spectrum of the surface density field $\Sigma$ is in this case related to that of the volumetric density field $\rho$ by $P_{\Sigma}(q_{\perp})=P_{\rho}(q_{\perp})\,L$, where $q_{\perp}$ is the two-dimensional Fourier wave number conjugate to the transverse length scale in proper units. This also assumes that the power spectrum is not too blue, obeying $\mathrm{d}\ln P_{\rho}(q)/\mathrm{d}\ln q<1$, where $q$ denotes the three-dimensional Fourier wave number.

If all mass is locked into subhalos of some mass $M$ and number density per unit $\ln M$ given by $n(M)=(\bar{\rho}/M)\,\mathrm{d}f(M)/\mathrm{d}\ln M$, the volumetric density power spectrum is (see Cooray & Sheth, 2002),

[TABLE]

where $\overline{\rho}$ is the mean density, and

[TABLE]

is the Fourier transform of the halo density profile (assumed to be spherical symmetric) at halo mass $M$. In principle, the subhalo mass distribution does not have to be the same as the field distribution in Eq. (9) considered in the previous section, but by the argument put forward in Section 4 we will use the approximation here that they are the same.

We now consider a long slab along the line of sight through the cluster halo. The total length of the slab is of the order of the cluster virialization scale $\sim$Mpc. The scale of the relevant transverse modes is comparable to the radii of the axion minihalos, which is orders of magnitude smaller. Since the coarse-grained density profile varies substantially from the inner region to the outer region of the cluster halo, it is invalid to assume that the volumetric density field is statistically homogeneous when deriving the power spectrum of the surface density field across the slab. Therefore, we divide the slab along the line of sight into many sub-slabs. Each sub-slab is sufficiently thin for the density field to have approximately homogeneous statistics, but is thick enough for the derivation of Appendix C to be valid. The surface density field through the entire slab has a power spectrum equal to the uncorrelated sum of the contributions from all sub-slabs.

We make another simplifying premise that all sub-slabs share the same fractional mass distribution as a function of halo mass, $\mathrm{d}f/\mathrm{d}\ln M$. Then the halo mass function is linearly proportional to the local mean density in each sub-slab. According to Eq. (10), the power spectrum for the surface density field is given by a summation over all sub-slabs, or an integral along the line of sight in the limit of a large number of sub-slabs,

[TABLE]

where $\int\,\mathrm{d}L\,\left(\cdots\right)$ is the line-of-sight integral. The fractional surface overdensity, defined as $\delta_{\rm 2d}:=\Sigma/\overline{\Sigma}$, has a power spectrum

[TABLE]

Under our assumptions, only the mean surface density $\overline{\Sigma}$ is needed for computing $P_{\delta_{\rm 2d}}(q_{\perp})$, provided that the (uniform) fractional mass distribution $\mathrm{d}f/\mathrm{d}\ln M$ and the halo density profiles are known.

Apart from minihalos residing in the cluster halo, intervening minihalos free floating in the intergalactic space are also expected. Indeed, in substructure tests utilizing galaxy scale lensing, intervening field halos integrated over a cosmological distance can be equally important or even dominant compared to subhalos (Despali et al., 2017). By comparison, the amount of substructure through a rich galaxy cluster is substantially larger. Assuming similar mass fractions of minihalos inside the cluster halo and in the intergalactic space, the relative importance of the two contributions can be inferred from comparing the surface density through the cluster $\overline{\Sigma}_{\rm cl}$ and the integrated DM surface density everywhere else along the line of sight $\overline{\Sigma}_{\rm los}$. For the two known highly magnified stars behind MACS J1149 (Kelly et al., 2018) and behind MACS J0416 (Chen et al., 2019; Kaurov et al., 2019), we find $\overline{\Sigma}_{\rm los}\approx 0.1$–$0.2\,\overline{\Sigma}_{\rm cl}$, which indicates that intergalactic minihalos should be subdominant.

It is convenient to measure the surface mass density in units of $\Sigma_{\rm crit}$, which gives the lensing convergence $\kappa=\Sigma/\Sigma_{\rm crit}$. Neglecting intergalactic minihalos, the lensing convergence has a power spectrum,

[TABLE]

We emphasize that this calculation is conservative because it only counts minihalos formed in isolation prior to cluster formation and excludes minihalos who themselves orbit intermediate mass subhalos (i.e. subhalos of subhalos).

In Figure 4, we present a calculation of the convergence power spectrum for the specific case of the line of sight toward the highly magnified star LS1 behind MACS J1149. In the calculation, we assume that axion minihalos inside the cluster halo obey the PS mass function Eq. (9), but neglect additional contributions to the clumpy surface density due to those minihalos having their own satellite minihalos.

In the left panel of Figure 4, we compute $\Delta_{\kappa}(q_{\perp}):=\left[q^{2}_{\perp}\,P_{\kappa}(q_{\perp})/(2\pi)\right]^{1/2}$, which quantifies the characteristic fluctuation in the lensing convergence $\kappa$ as a function of the transverse length scale. For $M_{0}=10^{-10}\,M_{\odot}/h$, $\kappa$ can fluctuate at the level of $\Delta_{\kappa}\sim 10^{-4}$ on scales $r_{\perp}=(2\pi)/q_{\perp}=10$–$10^{4}\,{\rm AU}/h$, which, as we will explain in Section 7, are probed by the light curves of the highly magnified stars during events of microlensing peaks. The right panel of Figure 4 shows that at a fixed wave number $q_{\perp}$, the convergence fluctuation originates mainly from minihalos with a typical scale radius $r_{s}\simeq r_{\perp}=(2\pi)/q_{\perp}$.

For an increased value for the parameter $M_{0}=10^{-6}\,M_{\odot}/h$, axion minihalos of $M\lesssim 10^{-2}\,M_{\odot}/h$ form earlier and thus are more compact. As a result, the typical convergence fluctuation on those scales can increase to $\Delta_{\kappa}\sim 10^{-3}$. If instead $M_{0}$ is several orders of magnitudes smaller, axion minihalos in the mass range of our interest become less dense, and the resultant convergence fluctuation is reduced.

The left panel of Figure 4 indicates convergence fluctuations at the level of $\Delta_{\kappa}\sim 10^{-4}$–$10^{-3}$ for a broad range $M_{0}=10^{-13}$–$10^{-6}\,M_{\odot}/h$, which translates into the mass scale of the earliest gravitational collapse $\sim 0.01\,M_{0}=10^{-15}$–$10^{-8}\,M_{\odot}/h$. The results are not very sensitive to the choice of the cluster formation redshift $z_{f}$ at which the PS mass function is evaluated.

5.2 Area covering factor

Assuming a uniform $\mathrm{d}f/\mathrm{d}\ln M$ throughout the cluster halo, we can calculate the optical depth to intersecting a minihalo of mass $M$ within the scale radius $r_{s}(M)$. Taking a geometric cross section $\pi\,r^{2}_{s}$, the differential optical depth is

[TABLE]

This can be interpreted as the area covering factor.

Figure 5 plots this optical depth as a function of the halo mass, for $M_{0}=10^{-10}$ and $10^{-6}\,M_{\odot}/h$, respectively. Taking $M_{0}=10^{-10}\,M_{\odot}/h$ as the example, we find that minihalos formed from isocurvature density fluctuations $M\sim 10^{-10}$–$10^{-6}\,M_{\odot}/h$ generally have $\mathrm{d}\tau/\mathrm{d}\log M\gg 1$. Within this mass range, as increasingly massive (and hence physically bigger) minihalos are considered, a single line of sight traversing the entire cluster halo intersects an increasingly larger number of minihalos. Since these minihalos are the major contributors to convergence fluctuations on scales of $\sim 10$–$10^{4}\,{\rm AU}/h$, the convergence field, by the Central Limit Theorem, is locally well described by a Gaussian random field with an isotropic power spectrum Eq. (14). The super-critical area covering fraction for the axion minihalos is in sharp contrast to their volume occupation fraction, which is tiny because of the characteristic density of minihalos being many orders of magnitude higher than that of the cluster halo.

As for the more massive halos $M\gtrsim 10^{-6}\,M_{\odot}/h$ formed from the adiabatic density fluctuations, the area covering factor starts to decrease. That part of the curve is shown as dotted, because the area covering factor is certainly overestimated toward larger $M$. Applying the PS mass function $\mathrm{d}f/\mathrm{d}\ln M$ to intracluster substructures severely overestimates the abundance of those halos inside the cluster halo as simulations suggest that many of them do not survive disruption. This is because the density hierarchy between those halos and the cluster halo is substantially narrower, and because across several decades of masses those halos collapse at similar redshifts and therefore all have similar densities (Diemand et al., 2006). It should be clear to readers that those halos are not the focus of this study.

In the case of $M_{0}=10^{-6}\,M_{\odot}/h$, we reach a similar conclusion for minihalos of masses $M\sim 10^{-4}$–$10^{-2}\,M_{\odot}/h$, although the overlap in projection for the least massive minihalos ($M\lesssim 10^{-4}\,M_{\odot}$) are not high enough to render the convergence field completely Gaussian. Finally, we note that our evaluation of Eq. (15) has left out minihalos inside subhalos or minihalos into higher levels of substructure hierarchy. Including those will further Gaussianize the convergence field and increase its power spectrum.

6 Dynamical disruption of minihalos inside clusters

We now discuss the dynamical processes that determine the survival rate of axion minihalos orbiting within larger CDM halos. We consider minihalos of mass $m_{1}$ and radius $r_{1}$ orbiting within a parent halo of mass $M$ and radius $R$. Those account for a fraction of the total host halo mass $\mathrm{d}f_{\rm sub}(m_{1};M)/\mathrm{d}\ln m_{1}$. We distinguish this from the PS mass fraction $\mathrm{d}f(m_{1})/\mathrm{d}\ln m_{1}$, which neither accounts for dynamic destruction of subhalos, nor for the fact that minihalos of mass $m_{1}$ exist at multiple levels of substructure hierarchy (subhalos of subhalos, and so on), all of which participate in dynamical processes.

We discuss three important dynamical processes that determine the destruction rate of minihalos:

Spiraling to the host halo center by dynamical friction, and tidal disruption at the high central halo mass density (which may be enhanced by baryonic processes in galactic centers such as disk and bulge formation). 2. 2.

Dynamical heating during high-speed encounters with other minihalos. 3. 3.

Dynamical heating during high-speed encounters with baryonic structures, including stars and molecular clouds.

We will show that among these destruction mechanisms, the most important is the second one if only DM is present, and the third one after stars have formed. Our estimates suggest that, over the age of the Universe, minihalos of our concern are unlikely to be disrupted by a significant fraction in a typical galaxy cluster when stars in member galaxies are taken into account. Some useful order-of-magnitude formulae are collected in Appendix D.

6.1 Dynamical friction and tidal disruption near the halo center

For a subhalo of mass $m_{1}$ orbiting inside a parent halo of mass $M$, the timescale for orbital decay due to dynamical friction is (Binney & Tremaine, 2011)

[TABLE]

After a few $t_{\rm df}$’s, the subhalo can spiral into the very inner part of the parent halo where they are destroyed by the tidal force of the dense central mass concentration.

Eq. (16) suggests that this mechanism is effective on minihalos if the parent halo is a larger minihalo so that $M/m_{1}$ is not too large. The orbital timescale is evaluated in Eq. (D). For example, a minihalo with $m_{1}=10^{-8}\,M_{\odot}/h$ can sink to the center of a parent minihalo with $M=10^{-6}\,M_{\odot}/h$ within $\sim 3\,$Gyr. Precise understanding of this hierarchical assembly and destruction of minihalos from white noise initial density fluctuations will require a dedicated study that employs numerical simulations.

However, many minihalos should evade dynamic friction constraints. First, the PS analysis in Section 4 has shown that a significant fraction of minihalos can be directly bound into either the cluster halo or its galactic-scale subhalos without first being assembled into a moderately larger parent halo or a series of slightly larger parent halos. For these minihalos, the humongous $M/m_{1}$ ratio ensures immunity to dynamic friction.

Second, even if a minihalo is first bound into a moderately larger parent halo and starts to spiral in, the parent may be subsequently disrupted by an even larger halo. In this case, the increased $M/m_{1}$ ratio allows the minihalo to survive for a longer time. This is especially relevant if the first parent halo is one on the portion of the colored curve in Figure 2 that levels off. In this case, this halo is very susceptible to disruption because subsequent hierarchical mergers take place quickly and involve small density differences. If the original minihalo is quickly “liberated” into orbiting a very large and low density parent halo, then $t_{\rm df}$ becomes certainly too long compared to the age of the Universe.

A subhalo can be destroyed even without spiraling in, simply by being randomly placed into a highly radial orbit toward the very inner part of a larger halo. This can lead to a high-speed encounter with either the central DM density cusp, or possibly with other highly condensed form of matter and compact objects. We will see in the following subsections that these fine-tuned encounters are not the most important processes compared to encounters with other minihalos and ordinary stars.

6.2 High-speed encounters with other subhalos

The large optical depth calculated for Eq. (15) of Section 5 and the fact that over the cluster’s age minihalos should have completed many orbits inside the cluster halo imply that direct encounters between minihalos are frequent.

We now consider the dynamical heating experienced by a subhalo of mass $m_{1}$ and internal velocity dispersion $\sigma_{1}$, when passing close to or moving through another subhalo of mass $m_{2}$ due to a random encounter. We focus on high-speed encounters for which the typical encounter relativity velocity222Our chosen notation $\sigma$ is motivated because the typical encounter velocity is set by the internal velocity dispersion of the host halo. This should not be confused with the peak height $\sigma(M,t)$ of the PS formalism. $\sigma$ is much greater than other velocity scales (Aguilar & White, 1985; Binney & Tremaine, 2011).

We first consider the case $m_{2}\simeq m_{1}$. Then, when the two subhalos fly through each other at an impact parameter $b$ smaller than the subhalo scale radii $r_{1}$, the velocity perturbation induced on any particle of the subhalos is $\Delta v\simeq\sigma^{2}_{1}/\sigma$, which is a small fraction of $\sigma_{1}$. This is because the tidal perturbation during the encounter is comparable to the orbital acceleration within one halo, but the duration of the encounter is shorter than the orbital time by a factor $\sigma_{1}/\sigma$. The fractional amount by which the internal subhalo energy varies owing to dynamical heating at each encounter with $b<r_{1}$ is therefore

[TABLE]

Dynamical heating scales as the second power of $\Delta v/\sigma_{1}$ if $\Delta v/\sigma_{1}\ll 1$. At the level of individual bound particles of halo $m_{1}$, the leading order fractional change in the particle’s energy is linearly proportional to $\Delta v/\sigma_{1}$, but can be either positive or negative depending on the orbital phase during the encounter. After the halo relaxes back to equilibrium, the average energy injection scales as $(\Delta v/\sigma_{1})^{2}$.

For unequal encounters, the minimum impact parameter is the larger $r_{s}$ of the two halos. Larger impact parameters are less important, because the induced velocity perturbations cale as $b^{-2}$, and the internal dynamical heating as $b^{-4}$, whereas the rate of encounters within $b$ increases only as $b^{2}$.

In the case $m_{1}<m_{2}$, an impact parameter $b$ smaller than the scale radius $r_{2}$ of the perturber halo $m_{2}$ leads to an encounter with the inner density cusp. For an NFW slope of the halo inner density profile, the cusp mass scales as $b^{2}$. The velocity perturbation $\Delta v$ is therefore comparable to that if $b\simeq r_{2}$, but the probability for a smaller impact parameter is suppressed as $\propto b^{2}$. This implies that inner density cusps are not the dominant sources of dynamic heating, and it is justified to take the larger $r_{s}$ of the two halos to be the minimum impact parameter.

At the minimum impact parameter, the induced velocity perturbation is maximized for $m_{2}\simeq m_{1}$. Encounters with a less massive perturber, with $m_{2}<m_{1}$, has $\Delta v\propto m_{2}$. For encounters with a higher mass halo with $m_{2}>m_{1}$, we focus on minihalos that collapse at high redshifts from the white noise isocurvature density fluctuations, for which the scale radius goes as $r_{2}\propto m^{5/6}_{2}$. In this case, the tidal acceleration scales as $m_{2}\,r^{-3}_{2}\propto m_{2}^{-3/2}$, and the encounter time as $r_{2}$, so $\Delta v\propto m_{2}^{-2/3}$. Since even in the worst case $m_{2}\simeq m_{1}$ we have $\Delta v/\sigma_{1}\ll 1$, particle ejections by a single tidal shock are unlikely, and hence the quadratic scaling Eq. (17) is justified.

We now show that the cumulative heating from many encounters is still dominated by nearly equal encounters with $m_{2}\simeq m_{1}$. For the case $m_{2}<m_{1}$, dynamical heating per encounter scales has $m^{2}_{2}$, while the encounter rate will scale as $m_{2}^{-1}\,\mathrm{d}f_{\rm sub}(m_{2};\,M)/\mathrm{d}\ln m_{2}$. If the PS mass function Eq. (9) is applicable to subhalos, $\mathrm{d}f_{\rm sub}(m_{2};\,M)/\mathrm{d}\ln m_{2}\propto m_{2}^{1/2}$ for minihalos and the encounter rate goes as $m^{-1/2}_{2}$. In principle, subhalos of subhalos or multiple levels of subhalos also contribute to the disrupting tidal field, which are not included in Eq. (9). Without a detailed calculation, we argue that including those $\mathrm{d}f_{\rm sub}(m_{2};\,M)/\mathrm{d}\ln m_{2}$ should depend logarithmically on $m_{2}$ because minihalos are not easy to disrupt through hierarchical merging, and accordingly the encounter rate goes as $m^{-1}_{2}$. In either case, the increase in the encounter rate does not compensate for the decreased heating per encounter. For the case $m_{2}>m_{1}$, dynamical heating per encounter decreases as $m_{2}^{-4/3}$, but the encounter rate increases only as $r_{2}^{2}\,m^{-1}_{2}\,\mathrm{d}f_{\rm sub}/\mathrm{d}\ln m_{2}\propto m_{2}^{2/3}\,\mathrm{d}f_{\rm sub}/\mathrm{d}\ln m_{2}$, again insufficient to compensate for the decreased heating per encounter. The above analysis should be valid unless the more massive minihalos have somehow increased their density by dissipative processes above the value reached through virialization at the collapse epoch.

Encounters between halos of comparable masses with impact parameters comparable to their common scale radii have a rate inside the cluster halo,

[TABLE]

where $R$ is the scale radius of the cluster halo. Given the fractional heating per encounter Eq. (17), the timescale of dynamic disruption is

[TABLE]

where we have assumed virial equilibrium for both the minihalo and the host halo $R/r_{1}=(M/m_{1})\,(\sigma_{1}/\sigma)^{2}$. This result can also be expressed in terms of $t_{\rm df}$ of Eq. (16),

[TABLE]

Since the velocity dispersion for minihalos is much smaller than that of the host halo (see Eq. (D)), we have $t_{\rm dispt}<t_{\rm df}$. The dominant process for minihalo destruction is therefore dynamical heating by high-speed encounters among comparable minihalos, where they repetitively go through each other before they are disrupted by gradual heating and expansion. Numerically, the disruption timescale is (where we have inserted the relation $G\,M/R\sim\sigma^{2}$ for the host halo by the virial theorem)

[TABLE]

Even for the largest possible $\mathrm{d}f_{\rm sub}(m_{1};M)/\mathrm{d}\ln m_{1}$, this is much longer than the age of the Universe within cluster-size or galaxy-size hosts, and even true for minihalos orbiting small CDM halos that form from the adiabatic density fluctuations, say one with $M=10^{6}\,M_{\odot}/h$ and $\sigma=5\,{\rm km/s}$.

6.3 High-speed encounters with stars

After the gas in halos condenses into stars and other compact objects, the potential for tidal disruption of minihalos is greatly enhanced. The reason for perturber minihalos of similar masses ($m_{2}\simeq m_{1}$) dominating the destruction rate by high-speed encounters rather than perturber minihalos of larger masses has to do with the fact that the more massive minihalos are also less dense, with $\rho_{s}\propto m_{1}^{-3/2}$ in the regime that minihalos collapse and assemble from the white-noise initial isocurvature density fluctuations. By contrast, for a compact perturber $m_{2}$, the minimum impact parameter is fixed by the size of the minihalo $m_{1}$, then the amount of dynamical heating per encounter increases with the mass of the perturbing object as $m_{2}^{2}$. If the number density of perturbers only decreases approximately as $m_{2}^{-1}$, dynamical heating should be dominated by the most massive perturbers.

Tidal interactions between generic small halos and passing stars have been extensively studied (Goerdt et al., 2007; Green & Goodwin, 2007; Angus & Zhao, 2007; Zhao et al., 2007). Here we present a simple order-of-magnitude estimate applied to axion minihalos. A minihalo of mass $m_{1}$ and radius $r_{1}$ encountering a star with mass $m_{\star}$, at an impact parameter $b$ and a relative velocity on the order of the host halo’s internal velocity dispersion $\sigma$ is subject to a tidal acceleration $\sim G\,m_{\star}r_{1}/b^{3}$ during the encounter time $\sim b/\sigma$, and suffers an impulsive velocity perturbation,

[TABLE]

Unlike mutual encounters between minihalos, the velocity perturbation induced by a stellar perturber can be large enough that minihalo particles are unbound after a single encounter. This happens for $\Delta v/\sigma_{1}\simeq 1$, which sets a minimum impact parameter $b_{\rm min}$ for the minihalo to survive a complete disruption (Goerdt et al., 2007; Schneider et al., 2010; Tinyakov et al., 2016),

[TABLE]

For a range of minihalo parameters, $b_{\rm min}$ is comparable or larger than the minihalo scale radius $r_{1}$. Provided that the overall tidal disruption rate is dominated by encounters with small impact parameters, we therefore conclude that a minihalo is disrupted whenever it passes within $b_{\rm min}$ from a perturber star.

Let $\bar{n}_{\star}$ be the mean stellar number density within the host halo (including the stellar populations in all member galaxies). The stellar disruption timescale for minihalos orbiting within the host halo can be estimated as

[TABLE]

Here the fiducial value $\bar{n}_{\star}=10^{5}\,{\rm kpc}^{-3}$ roughly corresponds to enclosing $\sim 10^{13}\,M_{\odot}$ of stars within a sphere of $\sim 300\,{\rm kpc}$. Note that the internal velocity dispersion of the host halo $\sigma$ cancels out in this calculation. Eq. (24) suggests that the most massive minihalos are the most vulnerable because they are the puffiest. For minihalos in the most interesting mass range concerning this work $m_{1}<10^{-6}\,h^{-1}\,M_{\odot}$, and in axion cosmology scenarios with $M_{0}>10^{-10}\,M_{\odot}/h$, stellar disruption is not significant throughout the current age of the Universe. More massive minihalos however are expected to have been destroyed by passing stars in an intracluster environment.

By chance minihalos can pass the central regions of larger halos at high speed on radial orbits. Apart from stars, supermassive black holes residing at the center of member galaxies are probably the most disruptive perturbers, for which the previous analysis remains applicable. Eq. (24) implies that the disruption timescale is inversely proportional to the average mass density of the perturbers within the cluster halo. Despite being much more destructive individually, supermassive black holes still make a smaller contribution to the total mass budget than the stars do and hence are not the major sources of disruption.

Finally, we note that even if all minihalos orbiting within the intracluster space are dynamically destroyed, we still expect a remaining contribution from an intergalactic population of minihalos. Those minihalos are not subject to dynamic disruption from encounters with stars. Some of those are isolated, while others are bound to larger intergalactic halos. The latter are susceptible to additional dynamic disruption from mutual encounters inside the parent halo.

7 Impact on microlensing light curves of highly magnified stars

Despite several theoretical uncertainties, our analyses in the previous sections lead us to the following physical picture if a Peccei-Quinn phase transition in the early Universe induces isocurvature density fluctuations on small scales:

• •

Minihalos below the mass scale of planets form by gravitational instability starting at the epoch of radiation-matter equality, with characteristic densities much higher than in the standard CDM cosmology.

• •

Plenty of these minihalos should survive over the age of the Universe as subhalos orbiting inside the halos of galaxy clusters.

• •

The high area covering factor of minihalos implies that a line of sight traversing a cluster lens near a typical lensing critical curve (with an impact parameter of tens of kiloparsecs) probes nearly Gaussian random fluctuations in the lensing convergence $\kappa$ (and hence also in the shear). For the mass scale parameter $M_{0}\sim 10^{-13}$–$10^{-6}\,M_{\odot}/h$, which corresponds to the first gravitationally collapsed minihalos having massses $\sim 10^{-15}$–$10^{-8}\,M_{\odot}/h$, the level of fluctuations is $\Delta\kappa\gtrsim 10^{-4}$–$10^{-3}$, on scales $\sim 10$–$10^{4}\,{\rm AU}/h$.

In most applications of gravitational lensing, convergence fluctuations of $\sim 10^{-4}$–$10^{-3}$ on the lens plane are too small to lead to any interesting observational consequences. However, caustic transiting stars behind a galaxy cluster lens provide extreme situations in which these minuscule small-scale fluctuations in the convergence can leave observable imprints.

7.1 Detectable scales of surface density fluctuations under microlensing

We consider the effect of convergence and shear irregularities on the total flux of a point source. The flux magnification of one image equals the inverse determinant of the lensing Jacobian matrix at the image position, whose matrix elements are expressed in terms of the convergence and shear. When an image is highly magnified ($|\mu|\gg 1$), the corresponding Jacobian matrix is fine-tuned to a level $\sim 1/|\mu|$ to be nearly degenerate. Any irregularity in the convergence and/or shear at a level $\sim 1/|\mu|$ substantially perturbs the image flux.

The detected highly magnified stars behind MACS J1149 (Kelly et al., 2018) and behind MACS J0416 (Chen et al., 2019; Kaurov et al., 2019) are observed thanks to magnification factors of hundreds to thousands. If the projected mass distribution in the cluster lens were smooth, the magnification would simply vary as the inverse of the angular separation from the image to the critical curve, and could reach values above a million when a luminous star crosses the cluster caustic. However, intracluster stars introduce intermittent microlensing flux variations and micro-caustic crossings at which a pair of micro-images dominate the flux. This makes the highly magnified stars more easily identifiable from their variability, even though the maximum magnifications reached are reduced to $\sim 10^{4}$ (Venumadhav et al., 2017; Diego et al., 2018; Oguri et al., 2018). Flux variations of up to a factor of ten have been observed from microlensing events in the fields of MACS J1149 (Kelly et al., 2018) and of MACS J0416 (Chen et al., 2019). During micro-caustic crossings, with typical durations of days to weeks, the star flux becomes susceptible to even minuscule nonsmoothness of $\Delta\kappa\sim 10^{-4}$ in the lens surface density.

During a micro caustic transit event, the two dominant micro-images are highly elongated to a length $\approx 5000\,{\rm AU}\,(|\mu|/10^{4})\,(R_{\star}/100\,R_{\odot})$. Surface density irregularities are detectable down to these scales, which are interestingly comparable to the axion minihalo sizes we have discussed.

7.2 Micro-Fold Model

We now demonstrate how minuscule surface density fluctuations imprint irregularities in microlensing light curves. Appendix E shows that geometric optics is applicable to our problem.

For a concrete case study, we consider the highly magnified star LS1 detected in a lensed galaxy at $z_{S}=1.49$ behind the cluster lens MACS J1149 at $z_{L}=0.544$ (Kelly et al., 2018). This sets the angular diameter distances to the source $D_{S}=1.79\,$Gpc, to the lens $D_{L}=1.35\,$Gpc, and from the lens to the source $D_{LS}=0.95\,$Gpc. Two resolved macro-images of LS1 are observed at a separation of $0.13^{\prime\prime}$ from the cluster critical curve, symmetrically positioned around a point where the macro lens model GLAFIC (Kawamata et al., 2016) predicts a surface mass density corresponding to a convergence $\kappa_{0}=0.83$, and a magnification for each macro-image of $\sim 300$.

The presence of microlenses implies that each macro-image actually consists of several unresolved micro-images. During a transit of a fold caustic (the most common caustic type), a pair of micro-images usually become much brighter than the combined flux from all the rest (which is in any case nearly constant). This justifies considering only the two dominant micro-images; instead of simulating a network of microlensing caustics induced by many intracluster stars, constructing a simple fold model for the micro-caustic crossing suffices. The fold model is described by two parameters: the local convergence $\tilde{\kappa}_{0}$ and an eigenvalue gradient vector $\boldsymbol{\tilde{d}}$ at the point between the two images along the micro-critical curve. The total smooth convergence $\kappa_{0}$ is split into a smooth contribution $\tilde{\kappa}_{0}$, including the DM and any diffuse baryonic component, and the average contribution from stellar objects, which is $\kappa_{\star}\approx 0.005$ near the LS1 images in MACS J1149 (Oguri et al., 2018). We use $\tilde{\kappa}_{0}=\kappa_{0}-\kappa_{\star}\simeq 0.83$ as a good approximation.

The vector $\boldsymbol{\tilde{d}}$ is perpendicular to the micro-caustic and, in the fold approximation, is related to the peak magnification $\mu_{\rm pk}$ at the time when the source stellar disk grazes the micro-caustic on the source plane:

[TABLE]

where $\tilde{\alpha}$ is the angle between the micro-critical curve and the degenerate direction of the local Jacobian matrix, and the stellar radius is $R_{S}$. This fold approximation is accurate when the highly magnified images are much closer to the micro-critical curve than the typical angular scale of variation of the vector $\boldsymbol{\tilde{d}}$.

The presence of microlenses generates a region around the macro-critical curve of the cluster lens where the micro-critical curves interact and join together, forming a network of width $\sim\kappa_{\star}|\boldsymbol{\tilde{d}}|^{-1}$ (Venumadhav et al., 2017). Within this network, the scale of variation of $\boldsymbol{\tilde{d}}$ is about the mean separation between microlenses, $\sim\theta_{\star}/\kappa_{\star}^{1/2}$, where $\theta_{\star}$ is the Einstein angular radius of each individual microlens of mass $M_{\star}$, $\theta_{\star}=(4\,G\,M_{\star}/D_{\rm eff}\,c^{2})^{1/2}$, with $D_{\rm eff}=D_{L}\,D_{S}/D_{LS}$. At the same time, the typical magnitude of $\boldsymbol{\tilde{d}}$ in the corrugated network is

[TABLE]

For demonstrating typical cases, we shall set in the next subsection $|\boldsymbol{\tilde{d}}|=\kappa^{3/2}_{\star}/\theta_{\star}$, with fiducial values for $\kappa_{\star}$ and $M_{\star}$ as in the above equation. We shall focus on situations of high magnification, with $\mu_{\rm pk}\gtrsim 10^{3}$–$10^{4}$, which requires the lensed star to be close enough to a micro-caustic. In fact, the typical average magnification of an image at an angular separation $\theta_{\star}/\kappa^{1/2}_{\star}$ from a micro-critical curve is $\mu\sim|1-\tilde{\kappa}_{0}|^{-1}\,(|\boldsymbol{\tilde{d}}|\,|\sin\tilde{\alpha}|\,\theta_{\star}\,\kappa_{\star}^{-1/2})^{-1/2}\sim|1-\tilde{\kappa}_{0}|^{-1}\,\kappa^{-1/2}_{\star}\ll 10^{3}$–$10^{4}$, for $\kappa_{\star}\sim 10^{-2}$. For the high magnification of our interest, the two dominant micro images have a much smaller angular separation than $\theta_{\star}/\kappa^{1/2}_{\star}$, justifying our use of a fold model locally.

The timescale of caustic transients depends on an effective source velocity $v_{t}$. Roughly, this is the relative transverse velocity between the source star and the lens cluster (defined in Eq.(12) of Venumadhav et al., 2017), which is constant. Although intracluster stellar micro-lenses and DM substructure have random velocities of $\sim 1000\,{\rm km/s}$ within the host cluster, their impact is dwarfed by the even larger apparent velocity of the micro images on the image plane owing to the large magnification factor. Therefore, for the numerical simulations in the next subsection, the assumptions of a stationary micro-fold caustic model and time independent realizations of axion minihalos along the line of sight are justified.

7.3 Numerical Examples

The analytical estimates in Eq. (25) and Eq. (26) suggest that microlensing peak magnifications can reach $\sim 10^{3}$–$10^{4}$ for reasonable lens and source parameters. We employ a numerical simulation to verify that fractional fluctuations $\sim 10^{-4}$–$10^{-3}$ in the lens surface density on appropriate scales are sufficient to imprint observable effects in the total flux.

On top of the micro-fold model, we add a spatially varying convergence perturbation $\Delta\kappa$ due to minihalos. Without having to generate individual minihalos, we are justified to model $\Delta\kappa$ as a Gaussian random field as shown in Section 5. We generate random realizations of $\Delta\kappa$ according to the homogeneous and isotropic power spectrum in Eq. (14), converting the linear Fourier wave vector to the angular Fourier wave vector $\boldsymbol{\ell}=D_{L}\,\boldsymbol{q}_{\perp}$. In the case of MACS J1149, the critical surface density is computed to be $\Sigma_{\rm crit}=2.3\times 10^{9}\,M_{\odot}/{\rm kpc}^{2}$, and $\bar{\Sigma}_{\rm cl}=0.83\,\Sigma_{\rm crit}$. The perturbed deflection required for inverse ray-tracing can then be computed in the Fourier domain $\Delta\boldsymbol{\alpha}(\boldsymbol{\ell})=(-2\,i\,\boldsymbol{\ell}/\ell^{2})\,\Delta\kappa(\boldsymbol{\ell})$.

In Figure 6 we present a numerical example of how minihalos collectively induce irregularities in the light curve during a micro caustic transit event, assuming parameters appropriate for the case of MACS J1149 LS1. In the absence of minihalos, the total flux rises smoothly as $|t-t_{*}|^{-1/2}$ after the two highly elongated micro images become dominant, peaks when the finite source effect kicks in, and then plummets when the two micro images merge and disappear. Minihalos cause sizable brightening and fading ”bumps” in the light curve despite the minuscule surface density fluctuations they induce. As the microlensing peak is approached, these irregularities become more prominent. In this example, the irregularities occur on a timescale of several days for $v_{t}\sim 300\,{\rm km/s}$. Unlike the usual microlensing signature of compact lenses, these irregularities lack abrupt changes; due to caustic crossings, microlensing light curves usually show asymmetric peaks, which have a slowly varying wing followed by an abrupt cutoff, or the time reversed behavior.

In the perturbed light curve the time of micro image merger shifts relative to that in the unperturbed one. This is due to coherent deflection perturbation generated by long wavelength modes. Causing merely a uniform remapping of the image plane coordinates, this bears no observable significance.

As visualized in Figure 6, we observe that the perturbing effect is one dimensional in nature: only non-smoothness along the degenerate direction of the micro-fold matters as the two micro-images are highly elongated along this direction.

The perturbed light curve varies depending on a number of factors, as we show in Figure 7. Firstly, the length of the elongated micro-image scales with the source size, which smooths out irregularities below certain scales. For a smaller source star, therefore, irregularities are present on shorter timescales.

Secondly, larger convergence fluctuations will result in larger magnitude irregularities. This can be the case if more minihalos survive dynamic disruption and are bound to intermediate sized subhalos, since the amount of minihalos contributing to surface density clumpiness will then be more than what Eq. (14) predicts.

Moreover, susceptibility to lensing perturbations depends on the $\boldsymbol{\tilde{d}}$ vector of the micro-fold model. Although in Figure 6 we have set $|\boldsymbol{\tilde{d}}|=\kappa^{3/2}_{\star}/\theta_{\star}$, actually $|\boldsymbol{\tilde{d}}|$ varies among micro-caustic crossing events. At a fixed distance from the micro-caustic, a smaller $|\boldsymbol{\tilde{d}}|$ corresponds to a larger magnification of the two micro-images. This means that the flux is more prone to surface density perturbations and is subject to larger irregularities.

If the typical convergence fluctuation $\Delta\kappa$ greatly exceeds the inverse magnification factor, even a micro-image can be strongly perturbed and break into yet smaller images. For example, this is likely the case for the scenario with $M_{0}=10^{-6}\,M_{\odot}/h$ shown in the left panel of Figure 4. The underlying physical principle is analogous to the disruption of macro-images into micro-images due to a super-critical number density of intracluster microlenses (Venumadhav et al., 2017; Diego et al., 2018; Oguri et al., 2018), or due to a high abundance of CDM subhalos (Dai et al., 2018b). A detailed analysis of this regime of strong image disruption is beyond the scope of this work, but is certainly an important situation to consider. Excessively large surface density fluctuations on length scales being considered here should be constrained by microlensing observations because they would smooth out stellar microlensing peaks.

8 Discussion

Minihalos are predicted to form in many cases when the dark matter is composed of axion particles that can solve the strong CP problem of QCD. So far, no viable method for detecting these axion minihalos through their gravitational lensing effect has been proposed. A generic problem to detect them is their low surface density. This paper proposes that axion minihalos can be detected during extreme magnification events when luminous stars cross micro-caustics produced by intracluster stars near the macro-critical curves of massive galaxy clusters. Substantial deviations from the predicted lightcurve of a fold micro-caustic crossing over timescales from hours to days will the signature of large numbers of axion minihalos which should superpose along a line of sight to produce nearly Gaussian fluctuations in the surface density.

The predicted lightcurve will in general be affected by limb-darkening effects in the source star. These effects should nevertheless be accurately predictable once the spectral type of the source star is known. An irregular behavior of the flux might, however, arise for other physical reasons. We now discuss some of these other possible sources of deviations from a predicted lightcurve.

8.1 Microlensing by Planet-Sized Lenses

Planets or planetesimals inside the foreground galaxy cluster may imprint distortions onto stellar microlensing light curves. Inside the micro-critical curve network, and without sub-stellar lenses, the typical separation between a segment of micro-critical curve and nearby microlens stars is $\sim(\kappa_{\star}\,\Sigma_{\rm crit}/M_{\star})^{-1/2}\approx 0.1\,{\rm pc}\,(\kappa_{\star}/0.01)^{-1/2}\,(D_{\rm eff}/1\,{\rm Gpc})^{1/2}\,(M_{\star}/0.3\,M_{\odot})^{1/2}$. Bound planets orbiting their host stars much more closely than this do not induce independent critical curves.

If any rogue or wide-orbit planets lie far (in projection) from any stellar critical curves, the minimal surface density for them to form a finer network of critical curves is $\kappa_{\rm ffp}\gtrsim 1/\bar{\mu}$, where $\bar{\mu}\sim 100$–$1000$ is the magnification factor of the image being observed. This would require a comparable amount of mass in rogue planets and stars, which is ruled out by Milky Way microlensing surveys implying $\kappa_{\rm ffp}\lesssim 10^{-3}\,\kappa_{\star}$ (Sumi et al., 2011; Mróz et al., 2017).

Free-floating planets which happen to lie in the vicinity of the stellar micro-critical curve can have greatly enhanced microlensing effects. If these planets have mass $M_{\rm ffp}$, corresponding to an angular Einstein scale $\theta_{\rm ffp}$, one requires $\kappa_{\rm ffp}\gtrsim(\theta_{\rm ffp}\,\tilde{d})^{2/3}\sim\kappa_{\star}\,(M_{\rm ffp}/M_{\star})^{1/3}$ for breakup of stellar micro critical curves to be common (Venumadhav et al., 2017). This implies a much higher abundance of free-floating planets than main-sequence stars, which is unlikely.

If the stellar micro-critical curve is indeed substantially disrupted by sub-stellar lens objects (see, e.g., Dai & Guerras, 2018, for an analogous consideration of quasar microlensing), finer variability structures are expected in the light curve near a stellar microlensing peak. However, in this case, signatures of planets or planetesimals would likely be sharp minor peaks rather than smoothly-behaved irregularities.

8.2 Blended Source Stars

The observed flux of the magnified star within its point spread function is likely to be contaminated by blended fainter stars. Although lensing of these fainter stars (or intrinsic variability if they are Cepheids) could contaminate the lightcurve, this is unlikely to happen because the highly magnified stars that are observed are usually among the most luminous in a galaxy, and they are sufficiently rare to make the superposition of two detectable lensed stars highly improbable. Contaminating microlensing is expected to produce sharp minor peaks in the light curve rather than the more smoothly-behaved irregularities from axion minihalos, as can be seen in Figure 6 and Figure 7. If the contaminator is of a different stellar type, variability will be chromatic, which is not generally expected from lensing of a single star.

8.3 Inhomogeneous gas

The dominant baryonic component inside a cluster halo is the ionized gas, which is shock heated to $T\sim 10^{8}\,$K and has a number density $n_{e}\sim 10^{-3}\,{\rm cm}^{-3}$. What is of our concern is the surface density fluctuations in the intracluster gas along the line of sight. However, what matters to the inquiry of this paper is inhomogeneities on minuscule scales $\sim 10$–$10^{4}\,$AU, which is orders of magnitude smaller than even the mean free path $\lambda_{\rm mfp}\sim 1\,$pc of the dilute gas. Although X-ray surface brightness measurements have clearly shown sizable turbulent density fluctuations $\Delta\delta_{\rm gas,3d}\sim 10\,\%$ (fractional) on scales of $10$–$1000\,$kpc (Schuecker et al., 2004; Kawahara et al., 2008; Churazov et al., 2012), the characteristic amplitude of fluctuation is expected to decrease toward smaller scales following a Kolmogorov-like power-law $\Delta\delta_{\rm gas,3d}(k)\propto k^{\alpha}$ with $\alpha\simeq-1/3$ (Gaspari et al., 2014), and is likely to cut off well before scales of $1/k\sim 10$–$10^{4}\,$AU due to dissipation (Luan & Goldreich, 2014). For fixed volumetric density fluctuation, the surface density fluctuation is further suppressed by projection $\Delta\delta_{\rm gas,2d}(q_{\perp})\sim\Delta\delta_{\rm gas,3d}(q_{\perp})\,(q_{\perp}\,L)^{-1/2}$ (Appendix C), where $L\sim 1\,$Mpc is the path length traversing the gas halo. Given the huge hierarchy between our projected scale of interest $1/q_{\perp}$ and $L$, roughly $q_{\perp}\,L\sim 10^{7}$–$10^{10}$, we hardly expect any detectable contribution to the lens surface density fluctuations from the hot diffuse gas. Cold and compact ISM gas structures of tiny sizes, if existent, may possibly perturb microlensing light curves, although their abundance and physical properties in the intracluster environment or in the surroundings of the BCG are poorly known. The lensing phenomenon we have been considering will probe or constrain any such small-scale baryonic structures.

9 Conclusion

Recently discovered lensed individual stars behind cluster lenses have the largest magnification factors among all currently known gravitational lensing phenomena. Flux amplification by $\sim 10^{3}$–$10^{4}$ is expected when they transit microlensing caustics induced by intracluster stars. We have shown that, during these extreme magnification events, their lightcurves are susceptible to minuscule non-smoothness $\sim 10^{-3}$–$10^{-4}$ in the lens surface mass density, across very small projected scales $\sim 10$–$10^{4}\,{\rm AU}$, which are hardly accessible with other observational means.

Within the strongly motivated paradigm of axion DM particle, solar-system sized minihalos in the mass range $M\sim 10^{-10}$–$10^{-6}\,M_{\odot}$ are predicted to copiously orbit the DM halo of galaxy clusters. This prediction should be generic if the Peccei-Quinn symmetry breaking occurs after inflation. Our calculations have shown, for a promising range of axion mass, that a large number of minihalos along the line of sight are able to collectively induce detectable surface density fluctuations. This result applies to $M_{0}=10^{-13}$–$10^{-6}\,M_{\odot}/h$, translating into typical masses for the earliest gravitationally collapsed minihalos in the range $\sim 0.01\,M_{0}\approx 10^{-15}$–$10^{-8}\,M_{\odot}$, an interval which includes the predicted parameters for the QCD axion. To our knowledge, monitoring microlensing variability of highly magnified stars is the first practical lensing-based method to probe axion minihalos.

Primordial black holes formed in the early Universe are another DM candidate (Zel’dovich & Novikov, 1967; Hawking, 1971; Carr & Hawking, 1974; Carr et al., 2016; García-Bellido, 2017). The mass budget of stellar mass black holes (Bird et al., 2016; Clesse & García-Bellido, 2017; Sasaki et al., 2016) has been subject to stringent limits set by microlensing and dynamics. Recently, however, Montero-Camacho et al. (2019) showed that the majority of DM may be made of black holes in the mass range $10^{-17}\,M_{\odot}\lesssim M_{\rm PBH}\lesssim 10^{-12}\,M_{\odot}$ without contradicting known astrophysical constraints. Such primordial black holes of asteroid masses may cluster into minihalos in a hierarchical fashion (see e.g. Inman & Ali-Haïmoud, 2019), created from a white noise power spectrum in a similar way as axion minihalos, and they may produce similar lensing effects as the axion minihalos although with some detailed differences arising from dynamical relaxation and the discrete nature of black holes that may give rise to distinct observational signatures.

Accurate and high cadence light curves of highly magnified stars around the times of microlensing caustic transits have not been acquired so far. This will require dedicated monitoring using either space telescopes or large ground-based telescopes to reach the extremely faint fluxes at optical or near-infrared wavelengths, executed over days to weeks around the peak time with high cadence. In the case of stars lensed by MACS J1149 and MACS J0416, one or two microlensing peaks were mapped out only at crude levels. Since microlensing caustic transits occur randomly, an approximate time window bracketing the moment of flux culmination will need to be forecasted based on low cadence pre-monitoring that can detect a trend in which the star gradually brightens following a $f_{\nu}\propto|t-t_{0}|^{-1/2}$ law. Because any detected highly magnified star is generally expected to undergo intermittent micro-caustic transits for many years to follow, and the opportunity to discover DM minihalos or any other dark or baryonic small-scale lumpiness is unique and highly rewarding, there is a strong incentive to design and carry out these programs.

The authors would like to thank Matthew McQuinn, Scott Tremaine and Matias Zaldarriaga for insightful discussions, and Ken Van Tilburg for helping us clarify issues regarding the expected mass scale of QCD axion minihalos. We acknowledge the use of the publicly available Python package Colossus to carry out cosmological calculations (Diemer, 2018). LD acknowledges the support from the Raymond and Beverly Sackler Foundation Fund. JM has been partly supported by Spanish fellowship PRX18/00444, and by the Corning Glass Works Foundation Fellowship Fund.

Appendix A Mass scale of first gravitationally collapsed minihalos

We evaluate the mass parameter $M_{0}$ defined in Eq. (1) for the standard cosmological thermal history. In this Appendix, we adopt the natural units for which $c=\hbar=k_{B}=1$.

Let $T_{0}$ be the radiation temperature at time $t_{0}$ at the onset of axion field oscillation. The axion mean density at present is $\bar{\rho}_{a0}=(3\,H_{0}^{2}\,\Omega_{a})/(8\,\pi\,G)$. In Eq. (1), we write $k_{0}=a(t_{0})\,H(t_{0})$. The scale factor at $t_{0}$ is

[TABLE]

where the temperature of the cosmic microwave background (CMB) $T_{\gamma 0}=2.725\,$K, the temperature of the relic neutrinos $T_{\nu 0}=(4/11)^{1/3}\,T_{\gamma 0}$, and $g_{*s}\approx 78$ is the effective number of relativistic degrees of freedom for entropy when the temperature is $\sim 1\,$GeV (Husdal, 2016). We calculate the Hubble parameter $H(t_{0})$ using the Friedmann equation during the era of radiation domination $H(t_{0})=[(8\,\pi\,G/3)\,(\pi^{2}/30)\,g_{*\epsilon}\,T^{4}_{0}]^{1/2}$, where $g_{*\epsilon}\approx 78$ is the effective number of relativistic degrees of freedom for energy density at a temperature $\sim 1\,$GeV (Husdal, 2016). Putting all pieces together, we find

[TABLE]

We note that our definition of $M_{0}$ is adopted from Fairbairn et al. (2018), which is a factor $4\,\pi^{4}/3$ larger than what is used in several other references (e.g. Davidson & Schwetz, 2016; Hardy, 2017). Following the definition of the latter authors, the typical mass for the smallest minihalos that collapse the earliest from gravitational instability should be two orders of magnitude smaller than $M_{0}$. The exact number is not sharply defined because even these minihalos contain finer axion substructures as uncovered in numerical simulations.

Appendix B Linear growth of isocurvature modes

Under the assumption that the DM is made of axions, the isocurvature matter modes under consideration here are superhorizon modes at the onset of axion field oscillation, but enter the horizon during the era of radiaton domination. These perturbations do not grow during the radiation-dominated epoch. After the epoch of radiation-matter equality, they lock onto a growing mode that scales linearly with the scale factor.

To verify this picture, consider a toy universe composed of a radiation fluid and a matter fluid that does not interact with the radiation. We normalize the scale factor such that it is equal to unity $a_{\star}=1$ when the mean radiation density equals the mean matter density $\rho_{r}(a_{\star})=\rho_{m}(a_{\star})$. Define a synchronous-comoving coordinate system in which the peculiar velocity of the matter fluid is always zero. We can derive the following set of equations (Efstathiou & Bond, 1986; Ma & Bertschinger, 1995)

[TABLE]

Here $a_{\star}$ and $H_{\star}=H(a_{\star})$ are the scale factor and the Hubble parameter at the epoch of radiation-matter equality, respectively, $\delta_{m}$ and $\delta_{r}$ are matter and radiation overdensities, respectively, $h$ is the trace of the spatial metric perturbation in the synchronous gauge (to be distinguished from the dimensionless Hubble parameter $h$ elsewhere), and we have defined $\chi_{r}:=\theta_{r}/(a^{2}\,H)$ with $\theta_{r}$ being the velocity divergence of the radiation fluid. The notation $(\cdots)^{\prime}$ represents the derivative with respect to the scale factor $a$.

We consider modes that enter the horizon well before the epoch of radiation-matter equality $c\,k/(a_{\star}\,H_{\star})\gg 1$. The isocurvature initial condition is set at some initial time $a_{i}$ in the superhorizon regime $c\,k/(a_{i}\,H(a_{i}))\ll 1$. At the initial time, the isocurvature initial condition requires $\delta_{r}(a_{i})=-\delta_{m}(a_{i})\,a_{i}$, $\chi_{r}(a_{i})=h(a_{i})=h^{\prime}(a_{i})=0$.

In the limit of matter domination $a\rightarrow\infty$, the matter mode satisfies

[TABLE]

which has a growing solution $\delta_{m}\propto a$.

By numerically solving the above set of linear equations, we find that deep in the era of matter domination $a\gg 1$ the solution for the matter perturbation is roughly

[TABLE]

where $C$ is an order-unity coefficient with a logarithmic dependence on $c\,k/(a_{\star}\,H_{\star})$. The exact value of $C$ is cosmology dependent, but its value is not much different from unity.

We note that in principle the initial isocurvature overdensity can be much greater than unity. In this case, virialized structures can form well before the epoch of radiation-matter equality (Kolb & Tkachev, 1994; Enander et al., 2017). A distribution of large initial overdensity patches were suggested by Kolb & Tkachev (1996) according to numerical results. Microlensing constraints, e.g. suggested by Fairbairn et al. (2017) and Fairbairn et al. (2018), have relied on ruling out extremely compact mini structures that collapse from such patches. We caution that it is unclear whether an extrapolation to huge initial overdensity values is physically plausible in realistic models, or whether such patches are sufficiently common to be important. This study restricts to initial isocurvature fluctuations not much greater than unity.

Appendix C surface density power spectrum

In this Appendix, we consider a slab with homogeneous and isotropic clumpiness property. We derive the relation between the two-dimensional power spectrum of the surface density field and the three-dimensional power spectrum of the volumetric density field.

We focus in the regime where the slab length $L$ along the direction of projection (i.e. the line of sight) is much larger than the dimensions perpendicular to the direction of projection (i.e. parallel to the plane of the sky). The surface density field can be written as

[TABLE]

Here $r_{\parallel}$ and $\boldsymbol{r}_{\perp}$ are the coordinates, in proper units, along and perpendicular to the direction of projection, respectively. The three-dimensional wave vector $\boldsymbol{q}$ can be decomposed into a perpendicular component $\boldsymbol{q}_{\perp}$ and a parallel component. The Fourier transform of the surface density field is

[TABLE]

Using $\left\langle\tilde{\rho}(\boldsymbol{q})\,\tilde{\rho}^{*}(\boldsymbol{q}^{\prime})\right\rangle=(2\pi)^{3}\,\delta_{D}(\boldsymbol{q}-\boldsymbol{q}^{\prime})\,P_{\rho}(q)$, we can calculate

[TABLE]

Therefore, the projected power spectrum is

[TABLE]

The regime we are interested in is $q_{\perp}\,L\gg 1$.

We specialize to the case that $P_{\rho}(q)\propto q^{\gamma}$ with $\gamma<1$. In this case, the dominant contribution to the integral comes from modes with $q_{\parallel}\ll q_{\perp}$ but $q_{\parallel}\sim 1/L$. We therefore obtain $P_{\Sigma}(q_{\perp})\sim P_{\rho}(q_{\perp})\,L$. Carrying out the $q_{\parallel}$-integral explicitly, we can fix the normalization factor of order unity and obtain

[TABLE]

This relation can be verified for a halo model of density inhomogeneity. Under the assumption that halos are uniformly distributed in space, $P_{\rho}(q)$ and $P_{\Sigma}(q)$ can be explicitly computed in terms of the halo mass function and the halo density profiles, and Eq. (C5) indeed holds.

Appendix D Dynamic scales of minihalos

Minihalos that form from the white noise isocurvature density fluctuations have extremely small internal velocity dispersions,

[TABLE]

Generally, $\sigma_{v}(M)$ is much smaller than the internal velocities of large CDM halos of any galaxies or galaxy clusters, $\sim 1$–$1000\,{\rm km/s}$. The orbital timescale is

[TABLE]

These timescales are much shorter than the age of the present Universe but are usually not so short compared to galactic orbital timescales. Even for minihalos of $M=10^{-10}\,M_{\odot}/h$, the dynamic timescale $\sim 30\,$kyr is quite long.

Inside a galaxy cluster, minihalos are flying around at the huge internal velocities of clusters, $\sigma_{\rm cl}\sim 1000\,{\rm km/s}$, taking less than a few years to cross the size of even the largest axion minihalos under our consideration. Therefore, dynamic interactions among minihalos are always much faster than their internal orbital timescales, which is in the regime of high-speed encounters.

Appendix E Validity of geometric optics

We show here that in the context of lensing discussed in this paper, geometric optics is valid if observations are done at UV, optical, or near-infrared wavelengths. The angular size of a typical star producing the observed micro-caustic crossing events regulates the maximum magnification factor, up to a wavelength at which wave diffraction effects become important. Simple estimates are presented in this Appendix to justify that.

For a single lens plane, the chromatic amplification factor for the amplitude of an electromagnetic wave is given by the following diffraction integral (see, e.g., Dai et al., 2018a)

[TABLE]

where $\lambda$ is the wavelength and $\tau(\boldsymbol{x};\boldsymbol{y})$ is the ray travel time as a function of the source position $\boldsymbol{y}$ and the image position $\boldsymbol{x}$,

[TABLE]

Considering only the two dominant micro-images during a micro-caustic transit, the Shapiro term $\phi(\boldsymbol{x})$ can be written as the sum of the contribution from the micro-fold model and a stochastic contribution from surface density fluctuations, $\phi(\boldsymbol{x})=\phi_{\rm fold}(\boldsymbol{x})+\Delta\phi(\boldsymbol{x})$.

Geometric optics is valid if the diffraction integral in Eq. (E1) can be replaced by complex Gaussian integrals around extrema of $\tau(\boldsymbol{x};\boldsymbol{y})$. This requires that $\tau(\boldsymbol{x};\boldsymbol{y})$ is locally well approximated by a quadratic function of $\boldsymbol{x}$ in the lowest order Fresnel zones, which in the case of an elongated micro-image with magnification $\mu$, has angular dimensions

[TABLE]

where we denote $D_{0}=D_{L}\,D_{S}/D_{LS}$. The validity of the quadratic approximation for $\tau(\boldsymbol{x};\boldsymbol{y})$ around extrema is typically limited by the longer dimension $\Delta\,x_{1}$.

For the micro-fold model, the third derivative of $\tau(\boldsymbol{x};\boldsymbol{y})$ is typically comparable to $|\boldsymbol{\tilde{d}}|$. So the departure from the quadratic approximation for the phase of the integrand in Eq. (E1) is roughly

[TABLE]

This is negligibly small for reasonable parameters considered in this paper.

If included, small scale surface density fluctuations may be the dominant contribution to the third derivative of $\tau(\boldsymbol{x};\boldsymbol{y})$. The characteristic amplitude can be estimated as $\sim\Delta_{\kappa}/(r_{\perp}/D_{L})$, where the convergence fluctuation has a typical amplitude $\Delta_{\kappa}$ on a typical transverse proper length scale $r_{\perp}$. Similar to Eq. (E4), the correction to the phase of the diffraction integrand is

[TABLE]

For a range of parameters relevant to this work, this correction is also unimportant.