Testing Primordial Black Hole Dark Matter with ALMA Observations of the Gravitational Lens B1422+231

The nature of dark matter remains as a puzzle. It is plausible that a fraction of the dark matter in galaxies and their surrounding dark matter halos could come from primordial black holes (PBHs). The idea that a fraction of dark matter could be PBHs in the mass range $10-10^{3}M_{\odot}$ has become more popular (e.g. Bird et al., 2016; Carr et al., 2016) after the gravitational wave detections of merging binary black holes with masses above $10M_{\odot}$ (Abbott et al., 2016, 2019). Observational constraints on the fraction of dark matter in a wide range of compact object masses show that PBHs of mass $1-10^{3}M_{\odot}$ might be a candidate for dark matter (Carr et al., 2016). The mass and abundance of PBH dark matter are constrained by microlensing (Green, 2016; Mediavilla et al., 2017), star clusters (Brandt, 2016; Green, 2016), dwarf galaxies (Koushiappas & Loeb, 2017; Zhu et al., 2018), and wide binary stars (Quinn et al., 2009; Yoo et al., 2004). Recently, there has been strong observational evidence for a 150 $M_{\odot}$ intermediate-mass black hole (IMBH), the remnant of a binary black hole merger (Abbott et al., 2020). Such IMBHs could also be PBHs. At highest accuracy, constraints on the fraction of dark matter as PBHs have to be calculated based on an extended PBH mass function. PBHs with a unified multi-modal mass function that peaks at around $10^{-5}$, $2$, $30$ and $10^{6}M_{\odot}$ could provide the dark matter and satisfy a variety of observational constraints (Carr & Kühnel, 2020; Carr et al., 2021a). A recently proposed method calculates the allowed mass range of primordial black holes as dark matter based on published constraints for the black hole mass function (Carr et al., 2016). In particular, bounds from microlensing and accretion show $10-10^{3}M_{\odot}$ PBHs can constitute a fraction or all of dark matter (Carr & Kühnel, 2020; Carr et al., 2021b). Unfortunately, the mechanism for forming IMBHs, either from massive stars or primordial black holes, is not well understood (see Miller & Colbert, 2004, for review).

Gravitational lensing provides some of the most compelling and direct evidence for dark matter (e.g. Clowe et al., 2006). General relativity predicts that light rays are deflected by massive objects. The deflection of light due to a gravitational lens galaxy or galaxy cluster can distort the shapes of the source galaxy images in the case of weak gravitational lensing, or produce multiple distorted images of the source galaxy in the case of strong gravitational lensing (Schneider et al., 1992; Kochanek, 2006). Furthermore, different mass distribution in the gravitational lens produce different image positions, shapes and flux densities. This makes gravitational lensing one of the most powerful tools for probing the distribution of dark matter.

Numerical simulations have shown that in the cold dark matter paradigm hierarchical structure formation produces a large population of dark matter substructures (DMSs) within galaxy-sized and cluster-sized dark matter halos (Klypin et al., 1999; Moore et al., 1999). The DMSs should potentially host dwarf galaxies or satellite galaxies, but the observed abundance of dwarf galaxies does not agree with the predicted abundance of DMSs (Bullock, 2010; Bullock & Boylan-Kolchin, 2017). This ‘missing satellite problem’ has recently been resolved by taking observational selection function into account (Nadler et al., 2021). The distribution of dark matter in dark matter halos is clumpy because of the presence of DMSs. Their signatures on the lensed images make galaxy-galaxy strong gravitational lenses the best targets for detecting DMSs observationally. Several techniques have been developed for detecting individual DMS (Vegetti & Koopmans, 2009; Vegetti et al., 2012; Hezaveh et al., 2013a; MacLeod et al., 2013; Nierenberg et al., 2014; Hezaveh et al., 2016b), measuring the DMS power spectrum (Hezaveh et al., 2016a; Cyr-Racine et al., 2019), and extracting the halo mass function (Gilman et al., 2020; Ostdiek et al., 2022) with strong gravitational lensing. DMSs near a strongly lensed image, either inside the lens galaxy halo or along the line of sight, perturb the smooth gravitational potential of the lens halo and change the magnification and the flux densities of lensed images (Dalal & Kochanek, 2002; Chen et al., 2003; Metcalf, 2005; Wambsganss et al., 2005; Xu et al., 2012). This leads to anomalous flux density ratios that cannot be explained by a smooth lens. The presence of DMSs affects the flux density ratios observed in all wave bands.

DMSs may not be the only dark perturbers causing flux density ratio anomalies. A population of PBHs as a fraction of dark matter could be another non-smooth mass component of the lens galaxy. Compact perturbers with masses $10-10^{5}M_{\odot}$ are considered IMBHs, which can be of primordial origins. A stellar component on top of smoothly distributed dark matter in the lensing galaxy has been shown to enhance microlensing fluctuations and broadens the magnification probability distribution toward both magnification and demagnification (Schechter & Wambsganss, 2002; Schechter et al., 2004). The flux density ratio anomalies observed in optical and X-ray bands are usually caused by stellar microlensing. Similarly, microlensing by a population of IMBHs can produce additional magnification on top of a macro-lensed image, producing an uneven magnification distribution. The distortion pattern of an individual IMBH along the line of sight in the lensing galaxy has an angular scale of a few $\mu$as to mas. Even though individual microlens’ Einstein radius is much smaller than the source size, the net magnification due to microlensing by a population of stars converges to a constant low magnification for large source sizes, shown by Barvainis & Ivison (2002) for 30 lensed quasars at submillimeter bands. A population of IMBHs can in principle produce analogous effects and have a significant contribution to the mm-wave flux density ratio anomaly. IMBHs could also be the answer to the issue raised by Xu et al. (2015) that DMSs are unlikely to be the full reason for lensed flux density ratio anomalies in the radio band. Schechter et al. (2004) demonstrated using microlensing simulations that the magnification probability distribution depends on the mass spectrum of the compact objects. Assuming a fraction of dark matter is PBHs, a strongly lensed quasar could show flux density ratio anomalies caused by quasar microlensing where PBHs act as microlenses. The degree of flux density ratio anomaly observed in a multiply imaged, strong gravitational lens system could place new constraints on PBH dark matter.

In this paper, Section 2 describes the first flux density ratio measurement in the mm-wave band for the strong gravitational lens, B1422+231. Section 3 presents the lens and source models for the ALMA observations. In Section 4 and 5, we examine whether PBHs alone can produce the quadruply-lensed flux density ratio anomalies by performing PBH microlensing simulations.

B1422+231 is a quadruply-lensed radio-loud quasar that is known for its flux density ratio anomaly. It has previously been observed in radio, infrared, optical and X-ray wave bands. We observed B1422+231 for the first time in the mm-wave band with ALMA. Continuum observations in Band 6 at 233 GHz over a total bandwidth of 7.5 GHz were carried out in two execution blocks on 2019 June 19 and 27 with ALMA using a hybrid configuration of C43-9 and C43-10 (Project ID: 2018.1.00915.S). The ALMA array consisted of 43 long baseline antennas and baseline lengths ranged from 83.1 m to 16.2 km. The calibrators were J1256-0547, J1427+2348 and J1436+2321. The total on-source integration time was 69.02 minutes. The data collected in the first execution block passed quality-assurance stage 0 and 2. The data collected in the second execution block were classified “SemiPass” in quality-assurance stage 0 due to bad phase transfer from the phase calibrator to the target. The Common Astronomy Software Applications package (CASA; (McMullin et al., 2007)) and ALMA pipeline (Version 6.2.1.7; Pipeline Version 2021.2.0.128) were used for calibration and imaging for data from both execution blocks. The measurement set from each execution block was calibrated separately, first with the ALMA pipeline and then phase self-calibrated. The phase self-calibration used a 80s solution interval for the first execution block and a 70s solution interval for the second execution block. The image signal-to-noise (S/N) for the second execution block is 2.3 times lower than that of the first execution block, due to the large noise in phase. Amplitude and phase calibration were performed after combining the two phase self-calibrated measurement sets in order to align the flux scales measured on two observation dates, assuming the source was not variable. To improve sensitivity to the weakest demagnified lensed image, imaging was carried out using natural weighting of the visibilities. Figure 1 shows the final image of the self-calibrated combined measurement sets.

The total continuum flux density of B1422+231 at 233 GHz we measure is $6.305\pm 0.082$ mJy, summing the flux densities from all four lensed images. This is slightly lower than the measured flux density of $7.5\pm 0.8$ mJy at 238 GHz by the Submillimeter Array (Keating, 2018). The missing flux density is likely caused by the fact that our maximum recoverable scale (MRS) is smaller than the extended emission along the lensing arcs near the image triplet (see Figure 3). A MRS of $0.322^{\prime\prime}$ is derived for the hybrid ALMA configuration, based on the $5^{\text{th}}$ percentile baseline length 810.7 meters. The angular sizes of the arcs are expected to be comparable to the image separations ($\sim 0.50^{\prime\prime}$ between image A and B and $\sim 0.82^{\prime\prime}$ between image B and C). The delivered data were incomplete relative to our proposal request due to proposal scheduling priority. We requested hybrid observations including the more compact C43-6 configuration, but no data was taken with this configuration. Assuming the radio synchrotron spectral index $-1.35$ derived from measurements at 15.0 and 22.5 GHz (Tinti et al., 2005; Stacey et al., 2018), we expect a flux density of 6.13 mJy at 233 GHz and 5.96 mJy at 238 GHz. Our measured flux density is consistent with these values, suggesting the detected compact emission at in Band 6 is mostly dominated by the synchrotron emission from the quasar’s active galactic nucleus (AGN) with almost no thermal dust emission from the dusty torus expected to surround the accretion disk (Urry & Padovani, 1995).

The starting point of our mass model for B1422+231 is a singular isothermal ellipsoid (SIE) plus external shear $\gamma_{x}$ (SIE$+\gamma_{x}$) model. Optical and near-infrared observations have shown that B1422+231 belongs to a group consisting of five nearby galaxies that provide sufficient external shear to produce the observed image configuration (Kundic et al., 1997). The source galaxy is at redshift $z_{s}=3.62$ and the foreground lens galaxy is at redshift $z_{l}=0.34$ (Patnaik et al., 1992; Tonry, 1998). We modeled the ALMA observation of B1422+231 in three steps: 1) lens modeling with image positions and/or flux densities as constraints assuming a point source; 2) visibility-space parametric source modeling while holding the lens model fixed; and 3) visibility space joint parametric modeling of the lens and source. A Hubble constant $H_{0}=67.7$ km Mpc^-1 s^-1 and a matter density parameter $\Omega_{M}=0.307$ (Planck Collaboration et al., 2016) are assumed.

If a fraction of dark matter is in PBHs, microlensing by a population of PBHs creates a scatter in lensing magnifications of strongly lensed images. We perform PBH microlensing simulations to quantify the effects of PBH dark matter on image component flux density ratios. The brute force inverse ray shooting algorithm, GPU-D, is used to simulate microlensing at the three cusp image components of B1422+231 (Thompson et al., 2010; Bate et al., 2010; Vernardos & Fluke, 2014). Such ray tracing simulations have traditionally been used to produce magnification maps and light curves for a source lensed by a random field of stars. The total surface mass density $\kappa$ has a smooth component represented by $\kappa_{s}$ and a compact object componenet $\kappa_{c}$ consisting of point masses with $\kappa=\kappa_{s}+\kappa_{c}$. Following Vernardos & Fluke (2014), we introduce a smooth matter fraction

The three image components of B1422+231, namely, A, B and C, are simulated separately with the same smooth matter fraction $f_{s}$, but different surface density $\kappa$ and shear $\gamma$ determined by the macromodel for the lens galaxy. Within each simulated area that represents an image component, $\kappa$ and $\gamma$ are assumed to be constant. This assumption is valid because the three unresolved image components have small angular sizes. The variation of $\kappa$ and $\gamma$ within each image component is negligible compared to the uncertainties of the macromodel. For every image component, the differences between the $\kappa$ and $\gamma$ values predicted by different best-fit macromodels in the literature are much larger than their variations within the image component (e.g. Mediavilla et al., 2009; Schechter et al., 2014). The exact $\kappa$ and $\gamma$ values from the macromodel are not important, because here we are more interested in the scatter of magnification due to PBH microlensing rather than the exact magnification. To also account for magnification uncertainties from the macromodel, PBH microlensing simulations have to be performed for all the macromodels sampled during lens modeling, which is outside the scope of this study. We choose to adopt the ($\kappa$, $\gamma$) values at A (0.38, 0.473), B (0.492, 0.628), and C (0.365, 0.378) from Schechter et al. (2014) for our simulations, which have been used for stellar microlensing simulations. Our best-fit ‘SIE (${\bar{S}}_{\nu}$)’ model predicts similar ($\kappa$, $\gamma$) values. We assume that the fraction of matter in dark matter is equal to $\Omega_{DM}/\Omega_{m}=0.842$ (Planck Collaboration et al., 2016) and that smooth matter consists of non-PBH dark matter and baryonic matter. Our simulations explored the scenarios where the fraction of dark matter in PBHs, $f_{PBH}$, was $10\%$ and $50\%$. The PBH surface density $\kappa_{PBH}$ is given by

The number of microlenses (i.e. PBHs), $N_{PBH}$, is

where $A$ is the area where PBHs are randomly distributed and $\left<M\right>$ is the mean mass of the PBHs following a power law mass function given below with a mass function exponent $\gamma\in(-1,+1)$ (Carr et al., 2017),

The mass of a PBH, continuous random variable $M$, satisfying the above mass function can be mapped to a continuous uniform distribution $\mathcal{U}(0,1)$

The two extreme cases where the power law mass function exponent $\gamma=-1$ and $+1$ were explored in our simulations. We generated $N_{PBH}$ random PBH massess and uniformly distributed them in an circular region with area $A$ in the lens plane. This circular region was larger than the area in which light rays were shot. To minimize edge effects, a smaller area within the actual area where light rays were traced was used for estimating the magnification probability distribution. The chosen area must be at least a few times the size of the source in order for the magnification map produced to contain a large enough statistical sample of random source positions on the map. Figure 5 shows an example magnification map of image component B of B1422+231 after convolution with our best-fit source Sérsic profile in Table 3. Image components with larger surface density have more densely populated PBHs. The large-scale horizontal patterns reflect the shear direction. High resolution magnification maps showing individual caustics are useful for simulating microlensing light curves for a source size comparable to the Einstein radius of the microlenses. However, for the purpose of finding the magnification of a source much larger than the average PBH Einstein radius, the PBH caustics do not need be resolved. Our simulations produced low resolution magnification maps that may contain multiple PBHs per pixel. For each $f_{PBH}$ and $\gamma$ combination, 25 independent 2000$\times$2000 pixels magnification maps were produced. The width of each pixel is equivalent to 3.6 to 26.5 times of the Einstein radius of the average PBH mass, depending on $f_{PBH}$ and $\gamma$. The deflection due to each PBH was still solved exactly. Because it was not necessary to resolve the caustics, the number of light rays needed was drastically decreased, averaging 153.33 light rays per pixel. This significantly reduced the computational cost. To calculate the magnifications of A, B and C at 233 GHz, we convolved the magnification maps with the Sérsic profile in our best-fit source model in Table 3, where the source major axis $a_{s}=9.03$ mas. A smaller source size would allow larger fluctuations in the magnification probability distribution, hence a higher probability for flux density ratio anomaly. Fast Fourier transform implemented in Astropy (Astropy Collaboration et al., 2013, 2018) with periodic boundaries that wrap around the magnification map was used for convolution. Varying the source position angle has an negligible effect on the magnification probability distribution.

Figure 6 shows the probability distribution function of the flux density ratios and the magnification probability distributions of image component A, B and C as functions of the PBH dark matter fraction, $f_{PBH}$, and the power law index $\gamma$ of the PBH mass function. We find that for our assumed source major axis of 9.03 mas, the probability distribution function of the flux density ratio (A+C)/B and the magnification probability distributions of the three cusp image components are all well-described by Gaussian distributions. This is in contrast to the asymmetric multi-modal distribution functions of magnification for small source sizes (Schechter & Wambsganss, 2002; Schechter et al., 2004, 2014). Because the assumed source size at 233 GHz is much larger than the Einstein radius of an average PBH, the large fluctuations of magnification probability distribution as the source crosses the caustics are essentially removed when convoluted with the source. Nonetheless, the widths of the probability distribution functions shown in Figure 6 indicate that there is a still non-negligible scatter in magnifications and flux density ratios caused by the microlensing of PBH dark matter. Table 4 lists the P-values and confidence intervals of flux density ratio found by our simulations.

We find that a larger fraction of PBH dark matter, $f_{PBH}$, or a larger power index $\gamma$ for the PBH mass function both noticeably widens the probability distribution function of the flux density ratios and the magnification probability distributions. With a power law index $\gamma=-1$, the number of PBHs sharply declines as a function of PBH mass. With a power law index $\gamma=1$, the number of PBHs is approximately constant in each mass bin in the specified mass range $10-1000M_{\odot}$. The difference in the power law index $\gamma$ most directly affects the average PBH mass. The number density of PBHs is proportional to the fraction of PBH dark matter. The lower right panel of Figure 6 shows that for given $f_{PBH}$ and $\gamma$, the probability distributions of magnifications normalized by macromodel predicted magnification ($\mu/\mu_{th}$) are very close to Gaussian distribution and they does not depend on the surface mass density or shear. More clumpy dark matter, either due to a higher fraction of PBH dark matter or fewer but more massive PBHs, increases the probability of the flux density ratio anomaly. The scatters of flux density ratios and magnifications in Figure 6 are conservative, because stellar populations are considered part of the smooth matter and any scatter produced by stellar microlensing is ignored. Figure 6 shows that a change in the PBH mass function alters the probability of having a large flux density ratio anomaly more significantly than a change of PBH dark matter fraction from 10% to 50%. Without the inclusion of a massive DMS, PBH dark matter between $10-1000M_{\odot}$ with a flat mass function alone can produce the observed flux density ratio anomaly. This highlights the need to include the effects of PBH dark matter when considering the causes of flux density ratio anomalies of multiply imaged strong lenses.

Assuming PBH dark matter is the only cause for the flux density ratio anomaly of B1422+231 in our ALMA observation, our simulations show that the probability of observing a flux density ratio (A+C)/B$>$1.434 is 31.54% if 50% dark matter is $10-1000M_{\odot}$ PBHs with a flat mass function (power law index $\gamma=+1$), and 14.25% if 10% dark matter is $10-1000M_{\odot}$ PBHs with a flat mass function (see Table 4). The cases where the PBH mass function is highly skewed towards low masses (power law index $\gamma=-1$) have negligibly small probabilities, irregardless of the PBH dark matter fraction. Carr et al. (2017) claims that any power law index$\gamma$ outside $(-1,+1)$ is not permitted. The logarithmic PBH mass function proposed by Carr et al. (2017) also highly skews towards low mass PBHs, so the probabilities assuming a logarithmic PBH mass function should be similarly small compared to those from the power law when $\gamma=-1$. Hence, if PBH dark matter has a mass function highly skewed towards low masses, its contribution to the flux density ratio anomaly should be extremely small even if a significant fraction of dark matter is PBHs. Recent constraints on $f_{PBH}$ allow 100% dark matter being PBHs in the intermediate mass range $10M_{\odot}<M<10^{3}M_{\odot}$ (Carr & Kühnel, 2020). However, more updated constraints suggest that the fraction of PBH dark matter in this mass range is very small, while PBHs below $10^{-10}M_{\odot}$ may still make up the majority of dark matter (Carr et al., 2021b). Such low mass PBHs can be treated as smooth matter microlensing simulations. Given that the P-values are small for $f_{PBH}$, it is very unlikely that PBH dark matter is the only cause for our measured anomalous flux density ratio if $f_{PBH}<0.1$.

The P-values are expected to be higher if there is PBH dark matter outside the $10-1000M_{\odot}$ mass range. An additional fraction of PBH dark matter from other mass ranges will result in more microlensing events and widen the probability distributions of the flux density ratios. Treating stars as microlenses instead of a smooth mass component will have a similar effect. In our simulations, stars in the mass range of $10-1000M_{\odot}$ are part of the smooth baryonic matter surface density. Therefore, our statistics in Table 4 are conservative estimates of the effects of PBH dark matter on flux density ratios.

DMSs are likely a major contributor to the flux density ratio anomaly of B1422+231 according to previous studies (Chiba, 2002; Nierenberg et al., 2014; Xu et al., 2015). Their models add a DMS modeled by a singular isothermal sphere to the macrolens model near image A. It is possible that the primary cause for the flux density ratio anomaly in B1422+231 is a DMS and the secondary cause is PBH dark matter. To more conclusively disentangle of the effects of PBH dark matter from DMS, more comprehensive mass models for multi-frequency observations, including low mass perturbers in a wide mass range such as DMs, IMBHs and stars, need to be compared. One key difference between DMS and PBH in causing flux density ratio anomalies is that PBHs as a fraction of dark matter affect all the image components, whereas usually one DMS affects one image component at a time. In principle, there are many low-mass DMSs, but their number density is still much lower than that of PBHs. One can potentially distinguish DMS-caused flux density ratio anomaly from PBHs by checking if only one of the image components has a anomalous flux density.

The biggest uncertainty in constraining the fraction of dark matter in PBHs is the unknown PBH mass function. The uncertainties from the smooth mass model and any massive DMS are the main source of systematic errors of the flux density ratio probability distribution functions. Incorporating the results of an ensemble of microlensing simulations with varied PBH mass functions into a full Bayesian analysis of a lens model including a smooth halo component, DMSs, and PBHs has the potential to provide constraints on PBH dark matter fraction and mass function. Quadruply-lensed compact mm bright quasars with larger flux density ratio anomalies than B1422+231 will give tighter constraints on PBH dark matter.

The mass distribution in the foreground lens galaxy of multiply imaged strong gravitational lenses are often not well described by simple smooth lens models. Quadruple lenses such as B1422+231 show flux density ratio anomalies in multi-frequency observations. In radio wave bands, the flux density ratio anomaly is caused by DMSs. We argue that PBH dark matter can be the cause of flux density ratio anomalies in the mm-wave band. In addition to DMSs, PBH dark matter may also be a secondary contributor to the flux density ratio anomaly. In optical and X-ray bands, stellar microlensing is the primary cause for flux density ratio anomaly, but both DMSs and PBH dark matter should be taken into account for lens modeling.

We present the first continuum imaging of B1422+231 using ALMA at 233 GHz. Four unresolved point sources are detected with a total continuum flux density of $6.305\pm 0.082$ mJy. The detected compact emission is consistent with the spectral energy distribution of synchrotron emission from the quasar’s AGN with almost no thermal dust emission. The flux density ratio (A+B)/C is measured to be $1.434\pm 0.017$, consistent with the measurements in radio and mid-infrared frequencies. The smooth lens model is not well constrained by the visibility data directly. Our lens models from mm-wave image plane modeling are similar to those from narrow-line [O III] (Nierenberg et al., 2014). Assuming a Sérsic source profile, the source major axis is estimated to be $9.03^{+0.15}_{-0.12}$ mas, corresponding to 66.9 pc at the redshift of this source, similar to the intrinsic [O III] size (Nierenberg et al., 2014).

Our ray tracing microlensing simulations of the three cusp image components show that PBHs in the intermediate mass range $10M_{\odot}<M<10^{3}M_{\odot}$ as a fraction of dark matter can produce flux density ratio anomalies. The magnification probability distribution of image components are well described by Gaussian distributions. A larger fraction of dark matter in PBHs and a less negative power law index for the PBH mass function can both increase the probability of flux density ratio anomaly. Assuming PBHs are the only cause of flux density ratio anomaly for B1422+231, we determine an upper limit to the probability of (A+C)/B$>$1.434 to be 31.54% for up to 50% dark matter being PBHs, and 14.25% for up to 10% dark matter being PBHs. Our measured flux density ratio (A+B)/C=1.434 is within 0.4806$\sigma$ confidence interval of the prediction for 50% dark matter in PBHs with a flat mass function. PBHs with a highly skewed mass function towards low masses has very low probability of being the only cause for the observed flux density ratio anomaly. This is the first study that quantifies the effects of PBHs on the flux density ratio anomaly in mm-wave bands. Our analysis places new constraints on the fraction of PBH dark matter and the PBH mass function using a quadruple lens with flux density ratio anomaly and ray tracing simulations.