of phenomena from cosmic infrared background excess to black hole merger rates. See also Khlopov ( 2010 ) for an overview of cosmological implications of PBH formation mechanisms. PBH nondetection can itself constrain the scale of isocurvature perturbations in cold DM (CDM) in the early Univ erse (P assaglia & Sasaki 2022 ). A small DM fraction of high-mass PBH could seed supermassive black holes and galaxy formation (Carr & Silk 2018 ), and in turn the observed population of supermassive black holes can also be used to constrain the PBH mass function (Cai et al. 2023 ).

A powerful and direct way to probe the PBH contribution to DM is strong gravitational lensing (see e.g. Treu 2010 , and references therein). Intrinsically point masses, PBHs are particularly ef fecti ve deflectors. Their observational signature depends only on their mass, with the deflection angle in terms of impact parameter ξ modelled as α = 4 GMc -2 ξ -1 . The method of gravitational imaging (Koopmans 2005 ;Vegetti et al. 2010 ;He et al. 2022 ) could in principle be used to detect individual PBHs of masses greater than 10 3 M (Banik et al. 2019 ), and lensing constraints from compact radio sources can also be used to constrain high-mass PBHs (Zhou et al. 2022b ). Because constraints from lensing are completely independent of others that have been used to constrain PBHs in a similar mass range, such as dynamical constraints (Carr & Sakellariadou 1999 ;Quinn et al. 2009 ;Brandt 2016 ), X-ray background constraints on accretion rate (Inoue & Kusenko 2017 ), or Lyman-α forest enhancement constraints (Afshordi, McDonald & Spergel 2003 ;Mack, Ostriker & Ricotti 2007 ;Murgia et al. 2019 ), lensing provides a vital crosscheck.

Anomalies in the ratios of flux between images of the same lensed source can reveal substructure in the lensing mass distribution. This technique, suggested initially by Mao & Schneider ( 1998 ), can probe structure at lower mass scales than those accessible with gravitational imaging. Such flux ratio anomaly studies rely on observations of lensed sources that are large enough to a v oid being affected by stellar microlensing; see Dobler & Keeton ( 2006 ) for an analysis of the effect of source size on flux ratio analysis. Examples of such sources include radio emission in radio-loud quasars (Mao & Schneider 1998 ;Metcalf & Madau 2001 ;Dalal & Kochanek 2002 ;Hsueh et al. 2020 ), mid-infrared emission from the hot dust in active galactic nuclei (AGNs; Chiba et al. 2005 ), and the narrow-line region of AGNs (Moustakas & Metcalf 2003 ;Nierenberg et al. 2014Nierenberg et al. , 2017Nierenberg et al. , 2020 ) ). Gilman et al. ( 2020a , b ) presented an analysis framework that uses the flux ratios among images in quadruply imaged quasars (quads) to constrain the properties of DM structure in strong lens systems. These techniques can be adapted to constrain a variety of DM models, including CDM, warm DM (Gilman et al. 2020a ), self-interacting DM (Gilman, Zhong & Bovy 2023 ), and fuzzy DM (Laroche et al. 2022 ), given a prescription for the halo mass function and density profiles of haloes.

In this paper, we present new constraints on the PBH parameter space by analysing the flux ratio anomalies in a sample of strongly lensed quasars observed in the narrow-line regime (Nierenberg et al. 2014(Nierenberg et al. , 2017(Nierenberg et al. , 2020 ) ). In Section 2 , we explain our method of sampling the posterior distribution of our PBH parameters of interest using forward modelling. In Section 3 , we present the results of our modelling and comparison to real data, and in Section 4 we discuss further expansions on this study. When necessary, we use the cosmology parameters of Planck Collaboration VI ( 2020 ) throughout this analysis, although we stress that our results do not depend sensitively on this assumption.

In this section, we first describe the goal of this paper, to obtain a posterior distribution on the PBH parameters of interest, which we achieve using an approximate Bayesian computation (ABC) forward modelling method. We first model the lens substructure using the method developed in Gilman et al. ( 2019 ), and then we model the effect of a possible PBH population.

We are striving to measure the posterior probability of DM model parameters; here, our likelihood function L can be written as

where D i refers to the observed image positions and flux ratios for a certain lens, θ f , M represents our target model parameters, M PBH and f PBH , m r is a certain lens model realization, and θ r is the set of non-PBH model parameters that we marginalize o v er. We use the method described by Gilman et al. ( 2020a ) to sidestep e v aluating this integral directly, which would require a computationally intractable exploration of a vast parameter space. This sidestepping is accomplished by forward modelling data, generating flux ratios from many sets of model parameters, and then comparing the results to the observed data via a summary statistic; from this process, we can extract θ f , M that represents our posterior probability distribution. This is an ABC method (Rubin 1984 , see also, e.g.

Sisson, Fan & Beaumont 2018 ) of creating a large set of stochastically varying simulated data and accepting simulations close to the real data to sample a posterior. ABC has been used in astrophysical forward modelling problems where a direct calculation of the likelihood function is infeasible and data can be simulated (see e.g. Weyant, Schafer & Wood-Vasey 2013 ; Akeret et al. 2015 ; Birrer, Amara & Refregier 2017 ).

We use the sample of 11 quads selected for flux ratio analysis by Gilman et al. ( 2022 , section 2.2) because the size of the source, either observed as [O III ] emission from the narrow-line region or CO (10-11) radio emission, is larger than the scale that would be affected by microlensing or image arri v al time delay, and the main lensing galaxy does not require modelling for a known stellar disc component (Hsueh et al. 2016(Hsueh et al. , 2017 ; ;Gilman et al. 2017 ). Photometry data used for each lens are referenced in Table 1 .

We generate a lens model using LENSTRONOMYfoot_2 (Birrer & Amara 2018 ; Birrer et al. 2021 ). The lens model is optimized to match the observed image position, with the added astrometric uncertainty. Any draw of parameters that does not match the observed image positions would be rejected in the posterior, so we reduce computation time by requiring the lens model fit the positions.

We compute the magnification, and thus the flux, of each image in our lens system model realization, then obtain the three flux ratios r model between the four images. Only flux ratios are used because the intrinsic source brightness is not known. We compare the forwardmodelled flux ratios to the observed flux ratios r obs with the summary statistic

(2)

MNRAS 522, 5434-5441 ( 2023)

We generate O10 5 -10 6 lens model realizations sampling from our parameter space from which we choose the 1500 lowest summary statistics to represent a sample of the posterior distribution. We construct a continuous approximation of the likelihood function for each lens by applying a kernel density estimate to the accepted samples, and multiply the resulting likelihoods to obtain the final posterior.

The lens and halo substructure modelling process follows from Gilman et al. ( 2019Gilman et al. ( , 2020a , b ) , b ). The lensing galaxy, or main deflector, is modelled as a power-law ellipsoid with external shear. The properties of the main deflector that are optimized during initial lens model fitting are the Einstein radius, centroid, ellipticity, ellipticity angle, and shear angle. If the main deflector has any known satellite galaxies, they are included in the model as a singular isothermal sphere mass profile. The main deflector mass M host , log profile slope γ macro , and shear γ ext are sampled in the forward model. Subhaloes are rendered from 10 6 to 10 10 M , from the lowest mass we are sensitive to to the highest mass of halo we expect to be entirely DM. The projected mass density sub and power-law slope α parametrize the subhalo mass function (SHMF),

where F ( M halo , z) is a function to scale the number density of subhaloes with main lensing halo mass and redshift as described in Gilman et al. ( 2020a ). The pivot mass m 0 is set to 10 8 M (Fiacconi et al. 2016 ). For the line-of-sight haloes, we use the Sheth-Tormen halo mass function (Sheth, Mo & Tormen 2001 ) with two-halo term ξ 2halo as a scaling factor to account for correlated structure near the host halo (see Gilman et al. 2019 ) and δ los as an o v erall amplitude scaling factor:

Given the PBH mass is distributed along with CDM subhaloes, the Sheth-Tormen mass function should be broadly applicable, but there may be an enhancement of the power spectrum on small scales caused by isocurvature perturbations from PBH (Afshordi et al. 2003 ;Gong & Kitajima 2017 ). Our constraint is more conserv ati ve because we do not take this enhancement into account. The free model parameters and priors are as follows:

(i) M PBH (M ), the PBH monochromatic mass, with a prior of 10 4 to 10 6 M chosen to include PBHs that are large enough to affect the flux ratios given the background source size but not larger than the minimum rendered halo mass;

(ii) f PBH , the PBH mass fraction of total DM, with a prior of 0-50 per cent;

(iii) sub , the SHMF normalization, with a prior of 0-0.1 kpc -foot_4 . We allow broad uncertainty in the number of subhaloes to account for uncertainties associated with tidal stripping;

(iv) α, the log slope of SHMF, with a prior ranging from -1.85 to -1.95 as predicted by Lambda CDM N -body simulations (Springel et al. 2008, Fiacconi et al. 2016 );

(v) δ los , the line-of-sight halo mass function scaling factor, with a prior of 0.8-1.2 that accounts for differences between theoretical models of the halo mass function (e.g. Despali et al. 2016 ) and uncertainties in cosmological parameters;

(vi) γ macro , the log slope of main deflector mass profile, with a data-moti v ated prior of 1.9-2.2 (Auger et al. 2010 );

(vii) σ source , the background source size, differing depending on whether the source is observed in narrow-line (M üller-S ánchez et al. 2011 ) or other regions (Chiba et al. 2005 ;Stacey et al. 2020 ) surrounding the background quasar listed for each lens in Table 1 ; (viii) M host , the mass of the main lens host galaxy [see Table 1 and Gilman et al. ( 2020a ) for a discussion of these priors constructed from individual lens data];

(ix) γ ext , the external shear in the main lens plane [see Table 1 , with ranges based on the individual lens data determined in Gilman et al. ( 2022 )];

(x) δ xy (mas), the image position uncertainty;

(xi) δ f , the image flux uncertainty.

References to photometric measurement information are listed in Table 1 . Our target parameters are M PBH and f PBH , and we marginalize o v er the others when they are sampled together in the posterior. For each realization, the model parameters are drawn from a prior distribution and the halo placement is stochastic. Some lenses have photometrically estimated redshifts (Gilman et al. 2020a ), so we sample the redshift probability distribution function and marginalize o v er it for those lenses.

Lens RXJ1131 + 1231 was modelled with two Gaussian source components to match the data of Sugai et al. ( 2007 ). Lenses with an imaged satellite companion are modelled with the companion in the source plane as a single isothermal sphere with position uncertainty of 50 mas.

If some fraction f PBH of DM exists in the form of PBHs, the distribution of these objects should follow a population of haloes with Navarro-Frenk-White (NFW) profiles (Navarro, Frenk & White 1997 ). Thus, the first step in our analysis is to generate a population of NFW haloes and subhaloes throughout the lensing volume. We create a realization of DM haloes and subhaloes using the LENSTRONOMY affiliate package PYHALO 2 (Gilman et al. 2021 ). We can calculate the mass fraction of DM rendered in haloes f halo , which we use to determine the number of PBHs that are clustered with the halo mass, and a stochastic distribution of line-of-sight and lens-plane subhalo masses, which we use to determine the clustered PBH position.

To determine the spatial distribution of the black holes at each redshift plane along the line of sight, we first compute the projected mass in DM at the lens plane from the population of NFW haloes distributed throughout the lensing volume. If we were to render haloes down to the minimum halo mass in CDM, then all of the PBHs would track the density of DM in haloes. Ho we ver, as we only render a fraction of the total mass of DM in haloes, a number N clustered = ( f PBH )( f halo )( ρ DM ( z))( V / M PBH ) will track the DM density in haloes. We distribute this population of PBHs with a spatial probability density that varies in proportion with the project mass in DM at each lens plane, as illustrated by Fig. 1 . The mass added in PBH is remo v ed from the rendered particle-DM subhalo mass. We place the remaining N smooth point masses randomly across each redshift plane, tracking the smooth background distribution of DM that we do not place in haloes. For each image of the lens, we add PBH at discrete redshift steps along the line of sight within a circular aperture of 0.24 arcsec for 10 8 M , scaling with the root M PBH to a minimum of 0.15 arcsec. For images closer than 0.24 arcsec, we add half the distance between the two points to the aperture and centre it at the mid-point between the images, treating the rendering area for those images as one aperture as illustrated in Fig. 2 .

In Fig. 3 , we have plotted ef fecti ve multiplane convergence for a lens model with M PBH = 10 5.5 M and f PBH = 0.4 alongside that of a model with no PBH substructure. This ef fecti ve multplane convergence is defined as the multiplane convergence (half the divergence of the effective deflection through the lensing planes α eff ) from the lens model minus the macromodel convergence, thus κ eff ≡ 1 2 ∇ • α effκ macro .On a convergence map, which corresponds to surface density, the point masses produce markedly different lensing signatures to the less centrally concentrated NFW profiles.

After the PBHs have been distributed in the lens model, the deflection from the new point masses is accounted for by refitting the lens model to the observed image positions. We raytrace through this final lens model to get the simulated image flux ratios, and calculate our summary statistic for the realization.

We present our constraint on PBH DM from the posterior distribution of our target PBH parameters for 11 lenses (in Fig. 4 ), which were combined and marginalized o v er the main deflector and subhalo parameters described in Section 2.2 . We obtained a 95 per cent upper limit on f PBH of 0.17 across the probed mass range. We see that there is a tentative anticorrelation between M PBH and f PBH , as we would expect.

Our constraint is plotted along with others in the same mass range in Fig. 5 . The constraint is stronger than that placed by radial velocity measurements of three wide binary systems that could be disrupted by a PBH population (Quinn et al. 2009 ), but it is partially within the bounds of the other four constraints. Ho we ver, our method is totally independent of the other bounds, and thus provides an important cross-check of the assumptions of other methods and their potential systematic uncertainties.

The X-ray accretion background constraint depends sensitively on assumption about the physics of gas accretion on to PBH and the possible subsequent formation of an accretion disc, the density of the interstellar medium (ISM), and PBH motion through the ISM. The constraint shown from Brandt ( 2016 ), similar to that of Quinn et al. ( 2009 ), is from the survi v al of the Eridanus II star cluster that would be dynamically heated into dispersal by PBH DM. This assumes that the Eridanus II cluster formed in place. The dynamical constraint placed by Carr & Sakellariadou ( 1999 ) assumes that PBHs will drift to the centres of galaxies, but this has been argued to be a v oidable if PBHs are regularly dynamically ejected as well (Xu & Ostriker 1994 ). The LSS constraint in Fig. 5 is from the effect of PBHs on the matter power spectrum as probed by the Lyman-α forest, which in turn depends on assumptions and modelling of its thermodynamics (Viel et al. 2013 ;Villasenor et al. 2022 ).

Finally, as shown by Banik & Bovy ( 2021 ), N -body simulations with DM particle masses comparable to the PBH mass range we consider result in small-scale perturbations to stellar streams. Interpreting this result in the context of PBHs suggests that streams also can constrain the contribution of PBHs to the DM. . Red corresponds to a density higher than that of the mean DM density, while blue corresponds to an underdensity. Black circles are plotted at each of the four quad image positions, and the black curves are the critical curves, which follow the region of maximum image magnification. Small-scale features in the convergence map that appear to track towards the origin are associated with black holes rendered around the path followed by the lensed light rays. Deformation of the critical curve by the PBH population suggests they will strongly perturb image flux ratios.

We develop a new method for including PBH substructure in a lens model for flux ratio analysis, and present independent constraints on the fraction of DM that could be composed of relatively massive PBHs. We obtained a constraint on f PBH less than 0.17 for M PBH = 10 4 to 10 6 M (95 per cent CL). The mass distribution for the PBHs in this work is monochromatic as a conserv ati ve constraint, but the limit can be converted to an arbitrary extended mass distribution via the method presented in Carr et al. ( 2017 ). This constraint is totally independent of others in the same mass range. In the spirit of a first application of this method, we make several simplifying assumptions throughout this process. We do not account for the effect of PBH formation on the assembly history of subhaloes and how that could possibly affect the mass functions and density profiles that we are also assuming. We allow for a very general parametrized form of these functions and marginalize o v er the parameter space to reduce the rigidity of our models. As samples of quads impro v e and our method becomes more constraining, we will revisit the simplifying assumptions.

In the future, these constraints will be impro v ed by applying the method to larger samples of lenses that are currently being disco v ered (Schmidt et al. 2023 ) and will be disco v ered in wide field surv e ys such as the Vera C. Rubin, Euclid, and Roman Observatories (e.g. Oguri & Marshall 2010 ). Lens systems can also be followed up with adaptive-optics-assisted instruments from the ground (Wright et al. 2019 ;Wizinowich et al. 2022 ). Forthcoming data from the JWST (Nierenberg et al. 2021 ) will allow us to push to lower PBH mass scales because JWST will measure flux ratios in the mid-infrared. This emission comes from a more spatially compact ( ∼1-10 pc) region around the background source. The minimum deflection angle that impacts our data is determined by the size of the source, so the more compact source size will allow us to push to lower PBH mass scales than we can currently measure.

Using 50 000 simulated lens model realizations of B1422 + 231, we tested the performance of our method by applying it to simulated data. We chose a realization with a low target mass and mass fraction of PBH and used the simulated flux ratios as the 'true' flux ratios in the computation of the summary statistic. From this, we obtain the posterior distributions shown in Fig. A1 . We repeat the e x ercise using a high PBH mass and mass fraction, and show the resulting inference on the right of Fig. A1 . The other parameters described in Section 2 were fixed in the middle of their uniform prior ranges. This process was carried out similarly for the lenses PS J1606 -2333 and WGD J2038 -4008, and the marginalized joint posterior distribution of all three lenses is shown in Fig. A2 .