{'fraction_non_alphanumeric': 0.03021971091996065, 'fraction_numerical': 0.030774663875085134, 'mean_word_length': 4.608941709111488, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Testing Primordial Black Hole Dark Matter with ALMA Observations of the Gravitational Lens B1422+231

We examine the ﬂux density ratio anomaly in the quadruply-imaged strong gravitational lens, B1422+231, and consider the contribution of 10 − 10 3 M (cid:12) primordial black holes (PBHs) as a potential dark matter constituent. We describe the ﬁrst ﬂux density ratio measurement of B1422+231 in the millimeter-wave band using the Atacama Large Millimeter Array (ALMA). This ﬁlls an im-portant multi-wavelength gap in our knowledge of this key lensed system. The ﬂux density of the quasar at 233 GHz is dominated by synchrotron emission and the source size is estimated to be 66.9 pc. The observed ﬂux density ratios at 233 GHz are similar to those measured in radio, mid-infrared and optical bands, which cannot be explained by a simple smooth mass model of the lens galaxy. We examine the probability of the ﬂux density ratio anomaly arising from PBH microlensing using ray tracing simulations. The simulations consider the cases where 10% and 50% of dark matter are 10 − 10 3 M (cid:12) PBHs with a power law mass function. Our analysis shows that the anomalous ﬂux density ratio for B1422+231 can be explained by a lens model with a signiﬁcant fraction of dark matter being PBHs. This study demonstrates the potential for new constraints on PBH dark matter using ALMA observations of multiply imaged strong gravitational lenses.


INTRODUCTION
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 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 (Abbott et al. 2016(Abbott et al. , 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 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 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 multimodal mass function that peaks at around 10 −5 , 2, 30 and 10 6 M 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 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 multi-ple 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 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 mi-crolensing. 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 µ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.

ALMA OBSERVATIONS OF B1422+231
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 ± 0.082 mJy, summing the flux densities from all four lensed images. This is slightly lower than the measured flux density of 7.5 ± 0.8 mJy 10 10 10 12 10 14 10 16 10 18 Obs.-frame frequency (Hz)   Lawrence et al. 1992;Impey et al. 1996;Chiba et al. 2005;Sluse et al. 2012;Pooley et al. 2012). Our ALMA Band 6 observation labeled with an arrow shows a flux density ratio of (A+C)/B = 1.434±0.017, consistent with measurements in radio and mid-infrared frequencies. The dotted and dashed horizontal lines show the predictions of the best-fit smooth lens model by Xu et al. (2015), (A+C)/B=1.247 and by Schechter et al. (2014), (A+C)/B=1.372. Data points with large deviation from the horizontal lines exhibit large flux density ratio anomalies.
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 is derived for the hybrid ALMA configuration, based on the 5 th percentile baseline length 810.7 meters. The angular sizes of the arcs are expected to be comparable to the image separations (∼ 0.50 between image A and B and ∼ 0.82 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).

Flux Density Ratio of Image Components
Four unresolved point sources are detected and labeled as image A, B, C, and D (see Figure 1). Image A, B, and C correspond to the brightest three image components in a cusp configuration of the quadruply-imaged B1422+231. The foreground lens galaxy is not detected, which is expected at this frequency. Image A, B, and C are individually fitted with two-dimensional elliptical Gaussian image components parametrized by peak intensity, peak pixel coordinates, major and minor axes, and position angle using CASA. The best-fit image properties are presented in Table 1. The image positions in our observation are consistent with those from radio, infrared and optical observations Lawrence et al. 1992;Impey et al. 1996;Nierenberg et al. 2014). The apparent alignment of position angles among A, B and C in Figure 1 only reflects the position angle of the beam, which has higher angular resolution along the tangential direction of the triplet's arc than the radial direction. The uncertainties of the fitted position angles for the three image components after devolution from the beam (in the last column of Table 1) are likely underestimated. The error estimates of the integrated flux to image components are calculated during the Gaussian fitting implemented in CASA. The error estimates to the fitting parameters are based on the formulation by Condon (1997). The image component positions from Gaussian fitting before self-calibration are the same as those in Table 1. Image component D has the largest position uncertainties from Gaussian fitting due to its low S/N, 1.3 mas in right ascension and 2.6 mas in declination. Figure 2 shows the flux density ratio of B1422+231 measured across the electromagnetic spectrum. We find that the flux density ratio (A+C)/B for B1422+231 at 233 GHz to be 1.434 ± 0.017, consistent with the values from radio, mid-infrared and most optical observations. The frequency-independent flux density ratios in radio, mm and mid-infrared observations suggest that the emissions at these frequencies likely originate close to the accre-tion disk of the AGN. Most flux density ratios across all frequencies in Figure 2 deviate significantly from the prediction from the best smooth lens model by Xu et al. (2015), which indicates that the mass distribution in the foreground lens galaxy is not smooth. The presence of DMSs is likely the primary contributor to the observed flux density ratio anomalies at low frequencies. The large scatter and deviation from the smooth model prediction in optical and X-ray observations suggests that stellar microlensing likely causes additional flux density ratio anomaly contributions.

LENS AND SOURCE MODEL OF B1422+231
The starting point of our mass model for B1422+231 is a singular isothermal ellipsoid (SIE) plus external shear γ x (SIE+γ 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 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 Ω M = 0.307 (Planck Collaboration et al. 2016) are assumed.

Lens Model
We first used only the image positions in Table 1 as constraints and modeled the lens system with the soft- Table 2. Best-fit lens model parameters. The first two lines are fitted using glafic. Note that ∆xs and ∆x l have opposite signs compared to ∆α cos(δ) in Table 1. The position angle θe is defined in degrees east of north in glafic, where +x direction is west and +y direction is north.  Table 3. Best-fit Sérsic source model parameters using visilens. The source position angle φs is defined in degrees counterclockwise from the lens major axis (see Figure 3), where +x direction is west and +y direction is north. ware package glafic (version 2; Oguri 2010), which minimizes χ 2 in the source plane. The SIE+γ x model includes seven parameters (columns 4-10 in Table 2): relative positions of the dark matter halo center to the phase center in RA and DEC directions ∆x l and ∆y l , Einstein radius θ Ein , ellipticity of the lens halo e, position angle of the lens dark matter halo θ e , external shear γ x and position angle of the external shear θ x . The source model is a point source with coordinates ∆x s and ∆x y relative to the phase center. We show the best-fit lens model with both the image positions and flux densities as constraints. However, because of the known flux density ratio anomaly of this system, the best-fit lens model using only image position as constraints is more reliable than that with flux densities as constraints. Table 2 shows that our model with only position constraints is similar to the published model obtained for [O III] images using only image positions as constraints (Nierenberg et al. 2014).

Source Model
Secondly, we reconstructed the source by holding the 'SIE (S ν )' lens model fixed and fitting a Sérsic profile using the software package visilens (Hezaveh et al. 2013b;Spilker et al. 2016), a Python package for model-ing interferometric observations of strong gravitational   lensing systems. The source emission was modeled with a Sérsic profile (Sérsic 1963) with seven parameters (columns 2-8 in Table 3): relative position of the source center to the lens center in the RA and DEC directions ∆X s and ∆Y s , total integrated flux density of the source F s , source major axis a s , Sérsic profile index n s , source minor to major axis ratio b s /a s and source position angle φ s . The software visilens was modified to make use of the importance nested sampling algorithm Multinest (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019 and its Python interface PyMultiNest (Buchner et al. 2014) for improved fitting efficiency and more precise estimation of the parameter uncertainties for the best-fit model in a multi-model parameter space. The calculation of the log-likelihood function in visilens remains unchanged. A constant efficiency mode of MultiNest with a sam-pling efficiency of 0.1 and 1000 live points was used. To reduce the number of visibilities to model, the continuum visibility data were averaged over 128 frequency channels in each spectral window spanning 1.875 GHz and 6.048 seconds in time. This is a prudent choice because the response R a to a point source resulting from this averaging is close to 1 (Equation (6.80) in Thompson et al. 2017). In comparison, the source model assuming the 'SIE (S ν )' lens model has higher log-evidence log Z than that assuming the 'SIE ([O III])' lens model in the last column of Table 3. The model image from our best lens and source model 'SIE (S ν )' is shown in Figure 3.1. All four image positions can be reproduced. Our model predicts an intrinsic source flux density of 0.2963 µJy and total magnification of 24.22, corresponding a total flux density of 7.176 mJy after lensing. The model-predicted extended arcs near the triplet images were not observed, due to the low sensitivity of the hybrid array to the angular scales of these arcs discussed in Section 2, resulting in a lower total observed flux density of 6.305 mJy. The model predicted Sérsic source major axis of 9.03 mas corresponds to 66.9 pc at the redshift of this source. This physical size is similar to the narrow-line [O III] size of 60 pc (Nierenberg et al. 2014). The fitted Sérsic index n s much smaller than 0.5 describes a sharp cutoff to the source surface profile. This is likely because of the lack of sensitivity to extended emission and the source not resembling the Sérsic profile.

Joint Modeling
Thirdly, we attempted modeling the visibility data directly by jointly fitting fourteen model parameters for the lens, external shear, and a Sérsic source using both the original version of visilens and our modified visilens with MultiNest sampler. The best-fit model image produced by the original visilens could not reproduce the relative image positions. The fitted image positions were more widely separated than those in the observation. The poor fits were caused by the low S/N of image D in the visibility data, resulting in large uncertainties in lens mass, lens position and source position. Large uncertainties for the source model parameters are also expected, given the under-constrained lens model and unresolved source. Jointly fitting the lens, external shear, and the source using the modified visilens with MultiNest sampler did not converge to produce a bestfit model. Therefore, we conclude that the joint lens and source model cannot be well constrained by the current data.
The failure of a smooth mass model is not surprising for this source. B1422+231 has been chosen for this study because of its known flux density ratio anomaly. Modeling visibilities directly makes use of both the position and flux density information of the image components. For a system with a flux density ratio anomaly, the best-fit smooth lens model constrained by the image component flux density information is bound to be inaccurate by design. Therefore, lens modeling for systems with flux density ratio anomalies are usually only based on the positions of image components. Previous work has shown that the best-fit SIE+γ x models with and without B1422+231 image flux densities are noticeably different (Nierenberg et al. 2014). Previous SIE+γ x models for B1422+231 observed in other frequencies can mostly reproduce the image component positions, but they also have difficulties matching the observed component flux densities (Kormann et al. 1994;Bradač et al. 2002;Chiba 2002;Nierenberg et al. 2014;Schechter et al. 2014). Their best-fit SIE+γ

PRIMORDIAL BLACK HOLE MICROLENSING SIMULATION
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 κ has a smooth component represented by κ s and a compact object componenet κ c consisting of point masses with κ = κ s + κ c . Following Vernardos & Fluke (2014), we introduce a smooth matter fraction (1) 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 κ and shear γ determined by the macromodel for the lens galaxy. Within each simulated area that represents an image component, κ and γ are assumed to be constant. This assumption is valid because the three unresolved image components have small angular sizes. The variation of κ and γ within each image component is negligible compared to the uncertainties of the macromodel. For every image component, the differences between the κ and γ 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 κ and γ 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 (κ, γ) 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 (S ν )' model predicts similar (κ, γ) values. We assume that the fraction of matter in dark matter is equal to Ω DM /Ω 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 P BH , was 10% and 50%. The PBH surface density κ P BH is given by The number of microlenses (i.e. PBHs), N P BH , is where A is the area where PBHs are randomly distributed and M is the mean mass of the PBHs following a power law mass function given below with a mass function exponent γ ∈ (−1, +1) (Carr et al. 2017 Table 3. The shape of the source Sérsic profile is illustrated in the bottom left corner. continuous uniform distribution U(0, 1) The two extreme cases where the power law mass function exponent γ = −1 and +1 were explored in our simulations. We generated N P BH 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 P BH and γ combination, 25 independent 2000×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 P BH and γ. 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. 2013Collaboration et al. , 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 P BH , and the power law index γ 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 welldescribed 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. 2004Schechter et al. , 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.

SIMULATION RESULTS AND DISCUSSION
We find that a larger fraction of PBH dark matter, f P BH , or a larger power index γ 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 γ = −1, the number of PBHs sharply declines as a func-tion of PBH mass. With a power law index γ = 1, the number of PBHs is approximately constant in each mass bin in the specified mass range 10−1000M . The difference in the power law index γ 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 P BH and γ, the probability distributions of magnifications normalized by macromodel predicted magnification (µ/µ 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 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 PBHs with a flat mass function (power law index γ = +1), and 14.25% if 10% dark matter is 10 − 1000M 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 γ = −1) have negligibly small probabilities, irregardless of the PBH dark matter fraction. Carr et al. (2017) claims that any power law indexγ 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 γ = −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 P BH allow 100% dark matter being PBHs in the intermediate mass range . The grey color bands show our measured flux density ratios within uncertainties from Table 1. The lower right panel shows the probability distribution functions of image components' magnification when fP BH = 0.5 and γ = 1, normalized by their respective theoretical magnification predicted by the macromodel. A Gaussian distribution is plotted for comparison. All plotted distributions are histograms with the same bin width 0.001. Table 4. P-value and Confidence Interval (CI) of the probability distribution of flux density ratios plotted in Figure 6. The P-values represent the cumulative probabilities of flux density ratio deviating away from the peak of the distribution more than the values measured by ALMA in Table 1. The confidence interval zσ is the proportion of flux density ratios that lie within the smallest symmetric interval (−zσ, +zσ) that includes the value measured by ALMA in Table 1, assuming a Gaussian distribution with standard deviation σ. more updated constraints suggest that the fraction of PBH dark matter in this mass range is very small, while PBHs below 10 −10 M 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 P BH , it is very unlikely that PBH dark matter is the only cause for our measured anomalous flux density ratio if f P BH < 0.1. The P-values are expected to be higher if there is PBH dark matter outside the 10 − 1000M 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 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.

Flux Density Ratio
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 mat-ter 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.

CONCLUSIONS
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 ± 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 ± 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 < M < 10 3 M 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σ 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.