General multipoles and their implications for dark matter inference

The flux ratios of strongly lensed quasars have previously been used to infer the properties of dark matter. In these analyses it is crucial to separate the effect of the main lensing galaxy and the low-mass dark matter halo population. In this work, we investigate flux-ratio perturbations resulting from general third-and fourth-order multipole perturbations to the main lensing galaxy’s mass profile. We simulate four lens systems, each with a different lensing configuration, without multipoles. The simulated flux ratios are perturbed by 10-40 per cent by a population of low-mass haloes consistent with CDM and, in one case, also a satellite galaxy. This level of perturbation is comparable to the magnitude of flux-ratio anomalies in real data that has been previously analyzed. We then attempt to fit the simulated systems using multipoles instead of low-mass haloes. We find that multipoles with amplitudes of 0.01 or less can produce flux-ratio perturbations in excess of 40 per cent. In all cases, third-or fourth-order multipoles can individually reduce the magnitude of, if not eliminate, flux-ratio anomalies. When both multipole orders are jointly included, all simulated flux ratios can be fit to within the observational uncertainty. Our results indicate that low-mass haloes and multipoles are highly degenerate when modelling quadruply-imaged quasars based just on image positions and flux ratios. In the presence of this degeneracy, flux-ratio anomalies in lensed quasars alone cannot be used to place strong constraints on the properties of dark matter without additional information that can inform our priors.


INTRODUCTION
The cold dark matter (CDM) paradigm, which posits that dark matter consists of non-relativistic, collisionless particles (e.g.Planck Collaboration et al. 2016), is successful at describing cosmic structure at scales larger than ∼ 1 Mpc (Springel et al. 2005;Planck Collaboration et al. 2020) and has been adopted as the standard in cosmology.CDM predicts a clumpy distribution of dark matter on sub-galactic scales and the existence of a large population of low-mass haloes (e.g., Vogelsberger et al. 2014;Schaye et al. 2015).On the other hand, warm dark matter models (WDM) predict a smaller amount of such objects with a less concentrated mass density profile (e.g.Bode et al. 2001;Viel et al. 2005;Lovell et al. 2014;Lovell 2020).The difference between CDM and still-viable WDM models is strongest at halo masses lower than 10 9  ⊙ (Hsueh et al. 2020), where most of these objects are predicted to be faint or even completely dark.Strong gravitational lensing allows us to detect them via their gravitational effect on the strongly lensed images.
In this paper, we focus on galaxy-scale strong lensing of unresolved sources, specifically quadruply-imaged quasars (quads).The image configuration is determined by the mass distribution of the lens and the position of the source.Low-mass haloes associated with the lens galaxy, called subhaloes, and haloes along the line of sight, called field haloes, can produce measurable changes in the relative fluxes of the lensed images due to the dependence of the image magnifications on the second derivative of the lensing potential.
★ E-mail: jascohen@ucdavis.eduThis method of investigating dark matter, known as flux-ratio analysis, cannot be used to precisely determine the masses and positions of individual haloes that cause the perturbations.Rather, the viability of a dark matter model is assessed based on the probability that its halo population could have produced the observed flux ratios, marginalised over the possible individual halo configurations.Dark matter models that produce larger numbers of low-mass haloes will lead to a higher incidence of lens systems that show so-called flux-ratio anomalies, in which a standard unperturbed smooth mass distribution cannot reproduce the observed flux ratios.In contrast, dark matter models that suppress the formation of low-mass haloes will lead to fewer and less significant flux-ratio anomalies in samples of lensed quasar systems.Single lens systems provide only a weak inference on dark matter models, since even in the presence of a large number of associated low-mass perturbers, these may be spatially distributed in a way such that no flux-ratio anomaly is produced.Thus, large samples of lens systems are needed.
Flux-ratio analysis was first proposed by Mao & Schneider (1998) and Metcalf & Madau (2001), and it was originally limited to subhaloes.Soon after, Dalal & Kochanek (2002) applied it to a sample of seven quads and reported results that were consistent with CDM simulations at 90 per cent confidence.
Follow-up studies argued for the inclusion of additional components that could also influence flux ratios, such as field haloes and stellar disks (Möller et al. 2003;Inoue & Takahashi 2012;Metcalf 2005;Despali et al. 2018).The contribution from field haloes is especially important given that they are often more numerous than substructures (Despali et al. 2018).Furthermore, field haloes should provide a cleaner test of dark matter models because their properties are not influenced by the tidal effects that can be so important for subhaloes.Baryonic components such as stellar disks have been discovered in real lens systems and, when included in the lens model, have successfully been able to reproduce the flux-ratio anomalies in those systems without having to resort to low-mass dark matter haloes (Hsueh et al. 2016(Hsueh et al. , 2017)).Similarly, more general explorations have shown that baryonic structures in lensing galaxies can mimic perturbations by low-mass haloes if not properly accounted for in the lens model (Gilman et al. 2017(Gilman et al. , 2018;;Hsueh et al. 2018).Sensitive high-resolution imaging can be used to estimate the contribution of baryonic structure in lensing galaxies, but other complexities may remain as confounding factors.
In this paper, we will investigate an important form of additional complexity for lens mass models, namely, the angular structure in the lensing galaxy, parameterised here as multipole perturbations.The most common mass profile used to model lens galaxies is the elliptical power law (EPL) with external shear (e.g., Tessore & Metcalf 2015).We will hereafter refer to this type of base model as the EPL model.The multipoles that we consider add Fourier-type perturbations to the angular part of the density profile, leaving the radial part unchanged.The use of multipoles is motivated by optical and infrared observations of elliptical galaxies, which show that the isophotes of many of them deviate from perfect ellipticity (e.g., Bender et al. 1988Bender et al. , 1989;;Cappellari 2016).These deviations can be modeled by simply-parameterized multipole components.While many treatments of elliptical galaxy isophotes focus only on fourthorder multipoles, Hao et al. (2006) present an extensive investigation of the surface brightness distribution in elliptical galaxies in which they fit both third-and fourth-order multipoles with a variety of orientation angles.
It is thus natural to consider multipole components in the mass distributions of galaxies as well.Demonstrating the impact of angular complexity in the lens galaxy on flux ratios, Evans & Witt (2003) and Congdon & Keeton (2005) showed that the joint inclusion of thirdand fourth-order multipoles with unrestricted orientation angles can reproduce many anomalies that had been observed at the time.We focus specifically on third-and fourth-order multipoles because they are expected to cause flux-ratio effects degenerate with those of perturbing haloes.Lower-order multipoles have effects that are analogous to changes in the macro-model parameters, and higher-order multipoles may introduce deviations from ellipticity that either produce greater than four images or are unphysical (Evans & Witt 2003;Congdon & Keeton 2005).Despite these findings, flux-ratio analyses since then have neglected to implement them completely.Hsueh et al. (2020) study some of the same lenses as Evans & Witt (2003) and Congdon & Keeton (2005), including angular structure in one lens in the form of an exponential disk, but they modelled all other lenses with only an EPL and shear.Recent flux-ratio analyses by Gilman et al. (2021Gilman et al. ( , 2022Gilman et al. ( , 2023) ) have included multipoles but in a specific and restrictive way.In those analyses, the orientation angle of the fourth order multipole is fixed to align with the EPL.Third order multipoles are not included.Gilman et al. (2024) inferred constraints on WDM from simulated lens systems including both third-and fourth-order multipoles, but the fourth order was fixed to align with the EPL.The effect of multipoles has also been considered in the context of the extended emission of lensed galaxies, where O' Riordan & Vegetti (2023) found that the inclusion of third-and fourth-order multipoles with unrestricted orientation angles could produce false substructure detections.
In this paper, we extend earlier work on lensed quasar flux-ratios (Evans & Witt 2003;Congdon & Keeton 2005) in several important ways.First, while those papers modeled real lenses with multipole components, our investigation uses simulated lenses so that we can directly compare the perturbative effects of low-mass haloes with those of multipoles.In addition, we consider the effects of thirdand fourth-order multipoles separately as well as jointly, and have more generality in our base models by allowing the power-law index to be different from the isothermal value.Our particular focus is an investigation of the potential for general third-and fourth-order multipoles to perturb the flux ratios of quadruply-imaged quasars in a way that is degenerate with perturbations from low-mass dark matter haloes.In Section 2, we describe our procedure for obtaining EPL base models for a sample of real lens systems.In Section 3, we detail the creation of four simulated lens systems from the combination of EPL base models and CDM low-mass halo populations plus, for one of the systems, a satellite galaxy.Section 4 describes how we model the simulated lenses in our sample using just an EPL model plus multipole components, and Section 5 presents the results.We discuss the implications of our results and future work in Section 6.

MODELLING REAL DATA: EPL𝛾
To quantify the effect of multipoles in realistic scenarios, we create simulated strong gravitational lens systems using image configurations from real lens systems taken from recent flux-ratio analysis studies (Hsueh et al. 2020;Gilman et al. 2020).To ensure applicability of our results across different image configurations, we select lens systems that fall into one of each of the general categories (see Figure 1 for visualizations): cross (WGD J0405-3308), fold (WFI 2026-4536) and cusp (B1422+231).We also select the cross configuration system PS J1606-2333, which has a luminous satellite associated with the main lensing galaxy, as the basis for a fourth simulated strong lens system.The satellite will allow us to investigate the degeneracy between multipoles and haloes beyond the low-mass range we otherwise consider (see Section 3.2).
For the mass model, we use an EPL.The corresponding dimensionless surface mass density (convergence) is given by where  is an ellitpical radius such that  2 = () 2 +  2 .The model parameters are the power-law slope, , axis ratio, , and scale length,  =   √ , where   is the Einstein radius.We define external shear with amplitude,  ext , and orientation angle,  ext .The satellite in J1606 is modelled as a singular isothermal sphere (SIS), an EPL profile with  = 1 and  = 1.Fitting only the observed image positions, we perform Markov chain Monte Carlo (MCMC) sampling to approximate posterior distributions of the mass model parameters for each lens system.The location of the luminous satellite in J1606 is fixed to the observed position, but we allow all other parameters, including the source position, to vary freely.For convenience, we will hereafter refer to the distributions generated in this step, including the system with the satellite, as EPL distributions.These sets of parameters describe the base models of our simulated data that will be perturbed either by low-mass haloes (EPL+CDM; Section 3) or by multipoles (EPL+MP; Section 4).
Table 1.Macro-model parameters used to create mock observations.The columns are   , the lens redshift,   , the source redshift,  E , the Einstein radius, dRA, the lens right ascension with respect to the observation center, dDec, the lens declination with respect to the observation center, , the power law slope,  , the ellipticity, , the position angle,  ext , the external shear strength,  ext , the shear angle, dRA src , the source right ascension with respect to the observation center, and dDec src , the source declination with respect to the observation center.

SIMULATED DATA: EPL𝛾+CDM
We create our four simulated lens systems by adding populations of CDM subhaloes and field haloes to the base models drawn from the EPL distributions described in Section 2. The results are mock quads with image positions that match those of a real lens system and flux ratios that are perturbed only by low-mass haloes.These lens systems do not contain any multipole components.We do not include both low-mass haloes and multipoles in any of our simulated data or models, as we intend to investigate whether the same observations can be produced with the presence of either individually.This represents a conservative approach to study their degeneracy.

Background source
Typically, flux-ratio investigations focus on emission from regions of the background objects that are large enough to avoid being affected by microlensing by stars in the primary lensing galaxy.These include mid-infrared emission from dust surrounding quasar accretion disks, which are typically smaller than 10 pc (Burtscher et al. 2013); emission from the narrow-line regions surrounding a quasar, which can extend up to 60 pc (Müller-Sánchez et al. 2011;Nierenberg et al. 2017); or radio emitting regions, for which individual observations give estimates of sizes smaller than 10 pc (Lee et al. 2017;Kim et al. 2022).Generally, as the size of the background source increases, it becomes less susceptible to flux perturbations from low-mass haloes (Dobler & Keeton 2006).In this paper, we want to quantify the degeneracy between low-mass haloes and multipoles in the scenario in which the effect of the former is maximal, hence the background sources in our mock observations and models are point-like.Their location in each realization is drawn from the MCMC chains associated with the modelling of the real data.In Section 6, we further investigate the effect of the source size and its implication for the degeneracy under study.

Low-mass halo population
To generate the CDM halo populations that we add to the EPL models, we largely follow the process described in Hsueh et al. (2020) with updated treatments of the mass-concentration relations for subhaloes and field haloes.All low-mass haloes are modelled as NFW profiles (Navarro et al. 1997;however, see Heinze et al. 2024).
For the field halo mass-concentration relation, we use that reported in Table 1 of Duffy et al. (2008) from N-body simulations.We use values derived using the virial radius definition of relaxed haloes between redshifts 0 and 2. Unlike Hsueh et al. (2020), we apply the associated scatter on the parameters.We follow the implementation of Despali et al. (2016), which is based on the approach introduced by Sheth & Tormen (1999), for the field halo mass function.We use their best-fitting parameters optimized over all considered redshifts and cosmologies.
We determine subhalo concentrations from a redshift-dependent mass-concentration relation extracted from the ShinUchuu N-body simulation (Ishiyama & Ando 2020;Moliné et al. 2023, also see O'Riordan et al. 2023 for more details).This relation is derived in terms of  max and  max , the radius of maximum tangential velocity and mass enclosed within it, as these more accurately describe the characteristics of haloes in simulations than the usual virial quantities.Our choice of mass-concentration relation results in more concentrated subhaloes than does the typical one from Duffy et al. (2008).We use a subhalo mass function that comes from fitting to the data in Lovell ( 2020) and has been reparameterised in terms of  max .After drawing the subhalo mass,  max , from the mass function, we draw the corresponding  max value from a log-normal distribution with mean and standard deviation Here,  =  +  and  =  + .Values for , , , , , , , and  are listed in Table 2, and  is the redshift.The NFW profile for a subhalo then has normalization where  s is the scale radius and  =  max / s = 2.16 (Bullock et al. 2001).We do not include tidal truncation or a dependence on the distance from the main lens centre of the mass-concentration relation because Despali et al. (2018) have shown the effects to be small compared to the scatter on the mass-concentration relation.
We generate populations consistent with predicted CDM subhalo and field halo mass functions down to a halo mass of 10 5  ⊙ , and we assume the total mass in substructure in the region of the lensed images to be ∼ 2 per cent of the total mass of the main lens in that region, which is roughly consistent with observational constraints (Dalal & Kochanek 2002;Hsueh et al. 2020;Gilman et al. 2020).This substructure fraction is higher than simulation predictions by Xu et al. (2015), but since we are testing the ability of multipoles to mimic low-mass haloes that strongly perturb the flux ratios, a bias towards models that have more perturbing haloes is a conservative choice.

Selecting realizations for simulated lens systems
Because we will proceed to stress test multipoles as they try to reproduce the flux ratios resulting from these EPL+CDM models, we generate 2000 realizations for each of the four main types of simulated lens in our sample (cross, fold, cusp, or satellite) and then select for each type the realization that produces the most extreme flux-ratio anomalies.These four EPL+CDM models should thus present the flux ratios that are most difficult for the EPL+MP models to reproduce.If models with multipoles but without low-mass haloes can fit perturbations produced by the most extreme halo populations, they should be able to do so in nearly all cases.We stress that though they are strong, the perturbations in our simulated lens systems are of comparable magnitude to flux-ratio anomalies in real observations (e.g.Nierenberg et al. 2020).Each simulated lens system contains total flux-ratio perturbations in excess of 10 per cent, and some images are perturbed beyond 20 per cent.Table 1 lists the macromodel parameters used in each of the four simulated systems, and Table 3 presents the image positions and flux ratios.We add to each flux in the models uncertainties that are based on the observations of the real lenses on which they were based.

MODELLING OF SIMULATED DATA: EPL𝛾+MP
In accordance with previous works examining the lensing effects of complex angular structure, we describe the convergence of multipoles in polar form Here,   and   are the standard multipole sine and cosine amplitudes, and  is the multipole order.For ease of interpretation, we also describe multipoles in terms of their overall amplitudes, and orientation angles, We restrict out focus in this paper to third-and fourth-order multipoles ( = 3, 4).To assess the possible degeneracy between the lensing effects of multipoles and haloes above our standard mass range, we do not include the SIS satellite galaxy in any of our models of that lens system.

RESULTS
In this section, we attempt to reproduce the simulated lens systems with multipoles in two distinct ways.First, we step through the parameter space of strengths and orientation angles for thirdand fourth-order multipoles separately to determine the range of flux ratio perturbations that each order, strength and angle can produce (Section 5.1).In this first part, the multipoles are not fit to the data.Rather, we simply show that configurations of multipoles exist that can produce flux ratio anomalies consistent with those produced by low-mass haloes.Second, we then use MCMC to fit the multipole and mass model parameters to the flux ratios of the simulated systems.(Section 5.2).

Independent investigation of third-and fourth-order multipoles
At this stage, we are not trying to fit to the flux ratios in our simulated sample, but rather to explore the dependence of the flux-ratio perturbations on the multipole amplitudes (  , where  is either 3 or 4) and position angles (  ).We do this by generating a grid of (  ,   ) pairs and, for each grid point, adding multipole components with these parameters to the 200 base EPL models.We do this exercise for  = 3 and  = 4 separately.While, judging from the isophotes, multipole amplitudes are not expected to be much larger than   ≈ 10 −2 in real galaxies (see Hao et al. 2006), we explore 11 multipole strengths ranging from 10 −3 −10 −1 with equal logarithmic spacing.For each multipole strength, we examine an evenly spaced set of 10 orientation angles encompassing the full range over which they are unique, i.e., from −60 to +60 degrees for the third-order multipoles and from −45 to +45 degrees for the fourth-order multipoles.To correct for any astrometric perturbations introduced by the addition of the multipole component, we optimize the macro-model parameters to fit to the simulated image positions after adding the multipole.
We show the results of this exploration in Figures 1 (third-order) and 2 (fourth-order).In each of the panels in the first two columns, the horizontal lines show the flux ratios of our simulated lenses while the points show the flux ratios produced by the EPL+MP models.In the left-hand columns we show how the perturbations change with multipole amplitude, showing only the points for values of   corresponding to the orientations of the highest-likelihood realizations.We calculate the likelihood of a realization from where x  and   are the image positions and flux ratios of a model realization (denoted  ) and simulated data (denoted  ).  and    are the simulated uncertainties.As expected, flux-ratio perturbations get larger with increasing multipole amplitude.However, the size of this effect is dependent on the configuration and the particular image in question.
The centre columns of Figures 1 and 2 show perturbations due to multipoles over the full range of orientation angles with   values fixed to those of the highest-likelihood realizations.All four lens systems show clear periodic behavior as the orientation angle changes.In all cases, some realizations bring the flux ratios closer to the simulated values than the macro-model alone.The simulated flux ratios in the fold and cusp lens systems can be reproduced within observational uncertainty by either third-or fourth-order multipoles, though the highest-likelihood realizations for the cusp system both have potentially unrealistic 1 amplitudes of   = 0.1.Though no thirdor fourth-order multipole perturbations can reproduce the simulated flux ratios in the cross and satellite systems, the magnitude of the discrepancy between model-predicted and simulated flux ratios can be significantly reduced by either order with reasonable amplitudes.
The flux-ratio perturbations induced by the SIS in conjunction with low-mass haloes in our simulated satellite system are comparable to those induced by low-mass haloes alone in other mocks.While satellite galaxies may be directly observable from their light, there may be a degeneracy between their inferred properties and multipole amplitudes.We leave the investigation of this potential degeneracy to a future work.

Joint investigation of third-and fourth-order multipoles
We now fit the mass model and multipole parameters simultaneously to both the flux ratios and image positions of our simulated lens systems.These data are the elements of a vector d.Similarly, the macro model parameters are the elements of a vector θ MM = {  , ...}, and the multipole amplitudes form a vector θ  = { 3 ,  3 ,  4 ,  4 }.From Bayes' theorem, the posterior distribution of the parameters given the data is The first term in the numerator Pr (d|θ  , θ MM ) depends only on 2 defined previously (see Equation 8).The prior probability Pr (θ  , θ MM ) = Pr (θ  ) Pr (θ MM ) is the probability of a given parameter value before the data is observed, based on other information.The normalisation of the posterior, or the evidence, Pr (d) can be ignored in this case as we only consider one model.The posterior we use in practice is then The calculation of this many-dimensional posterior is intractable, so we use MCMC to obtain samples of θ  and θ MM .The density of these samples represents the posterior probability distribution.To conduct the MCMC sampling, we use the emcee 2 ensemble sampler with 120 walkers and 10,000 burn-in steps that are discarded, followed by 20,000 recorded steps.We use broad uniform priors on the mass model parameters θ MM and use a normally distributed prior with mean zero and standard deviation of one per cent on the multipole amplitudes.
With both third-and fourth-order multipole parameters free, all flux ratios in the simulated lens systems can be fit within 1. Figure 3 displays flux ratio histograms from the joint MCMC sampling of both multipole orders and macro-model parameters (EPL+MP).For reference, we also include the flux ratio histograms that result from the MCMC sampling with only an EPL.These EPL-only models are unable to reproduce the simulated flux ratios within observational uncertainty in every case.Generally, the addition of third-and fourthorder multipoles to the model makes the model flux ratio distributions narrower and roughly centered around the simulated values.In other words, when multipoles with amplitudes that are consistent with the amplitudes of observed isophotes are included in the model, the simulated flux ratios do not appear anomalous despite the fact that they were perturbed by a population of CDM low-mass haloes, and, in one case, also a satellite galaxy.
Figures 4, 5, 6 and 7 show the resulting posterior distributions for the third-and fourth-order multipole parameters for all four lens systems.Non-zero third-and fourth-order multipole amplitudes are preferred in every case, and most orientation angles are constrained within 20 degrees.
Isodensity contours for the best-fitting models for each system and their underlying macro-models are shown in Figure 8.The shapes of the contours for the cross and fold systems generally become 'boxier'.However, offsets from perfect alignment between the fourthorder multipole and the EPL major axes along with the inclusion of third-order multipoles, makes this effect more irregular than pure 'diskiness' or 'boxiness'.Deviation between the perturbed and unperturbed contours for the satellite system occurs roughly in alignment with the location of the missing satellite galaxy.It is worth noting again that the simulated data were generated without any multipoles.

DISCUSSION AND CONCLUSIONS
We find that there is a substantial degeneracy between the perturbative lensing effects of low-mass dark matter haloes and multipoles when modelling quadruply-imaged quasars.Individually, third-and fourth-order multipoles with amplitudes around or below 0.01 can produce flux-ratio perturbations over 40 per cent.This is sufficient to reproduce perturbations induced not only by a population of lowmass CDM haloes, but also a satellite galaxy.Our results agree well with early investigations of this topic by Evans & Witt (2003) and The range of flux-ratios as a function of third-order multipole parameters  3 (left) and  3 (middle).Only flux ratios from models with position error less than 1 for each image are included.The horizontal lines and shaded regions represent the flux-ratios and uncertainties in the simulated data, is generated by an EPL+CDM model.The filled dots show the flux-ratios of each image for each multipole parameter value.Open circles represent the flux-ratios from the EPL+MP model that most closely matches the simulated flux ratios.The solid, dashed and dotted vertical black lines in the left column indicate the mean, 1 and 2 values, respectively, of the third-order multipole amplitude distributions from elliptical galaxy isophotes (Hao et al. 2006).The solid black line in the middle column indicates the orientation angle of the EPL major axis in each simulated lens system.The rightmost column shows the image configuration for each lens system, and the color of each image corresponds to data for its respective flux ratio in the other columns.The black point marks the image used in the denominator of the flux ratios.Congdon & Keeton (2005).From the analysis of real data, they find that all but the most extreme flux-ratio anomalies, such as those in the cusp systems B2045+265 and B1422+231, which differ from EPL predictions in excess of ∼50 per cent, can be fit using multipoles with amplitudes smaller than 0.01.

A more general treatment of multipoles
In our analysis of the ability of multipole components to reproduce flux ratio anomalies, we utilized two important general treatments of the multipoles: we included third-order multipoles in addition to the more standard fourth-order one, and we did not fix their position angle to the major axis of the underlying EPL model.
The freedom allowed to the position angle plays a key role in the degeneracy between multipoles and low-mass haloes.As we have shown, the orientations of multipoles that most closely reproduce the perturbed flux ratios in our simulated data do not always align with the EPL.Evans & Witt (2003), Kochanek & Dalal (2004) and Congdon & Keeton (2005) also found that a larger portion of observed fluxratio anomalies can be reduced or eliminated when multipole angles are allowed to vary freely as opposed to being fixed.
Recent flux-ratio analyses, where fourth-order multipoles have been included (e.g.Gilman et al. 2021Gilman et al. , 2022Gilman et al. , 2023Gilman et al. , 2024)), have fixed the orientation angles to align with the EPL major axis.However, neither the total mass density nor the light distribution of real galaxies seem to suggest this assumption to be generally valid.Using high-angular resolution observations of real gravitational lens systems with an extended source, Powell et al. (2022) and Stacey et al. (2024) have found that the multipoles can be misaligned with respect to the underlying EPL.Additionally, tilted dark matter haloes are commonly observed in numerical simulations (e.g Han et al. 2023).At the same time, the observed isophotes of elliptical galaxies only have order  = 4 perturbations aligned with the EPL for larger values of  4 , in excess of one per cent, i.e. for truly boxy isophotes (Hao et al. 2006, their fig. 4).The fourth-order multipole in the best-fitting model for our cusp system has an amplitude of 1.5 percent, and it does not align with the EPL.The best-fitting models for our other systems all have fourth-order multipoles with amplitudes ≲ 0.5 per cent.

Isophotes and iso-density contours
An important comparison to our work is the study of isophotes in elliptical galaxies, since multipole patterns in the former are often  Figure 3. Flux ratio histograms created using random samples of 10,000 models from MCMC chains for all four simulated lens systems: cross (a), fold (b), cusp (c) and satellite (d).For comparison, we show distributions resulting from using EPL models (blue) and EPL+MP models with third-and fourth-order multipoles included simultaneously (red).Vertical black lines show flux ratios from the simulated data, and the surrounding grey regions show the associated uncertainties.Corner plots from simultaneous MCMC sampling of EPL, shear and both multipole order parameters on the cross lens system.MCMC was conducted in the sine/cosine amplitude basis (left), and we also present the same data transformed into the overall amplitudes and angles (right) for ease of interpretation.
Only multipole parameters are displayed, but all EPL and shear parameters, along with the source position, were allowed to vary freely.The red-dashed contours and histograms show the Gaussian prior distributions used for the   and   parameters.used as motivation for including multipole components in the mass distribution of lens galaxies.
The amplitudes of the multipoles that we find to be highly degenerate with low-mass haloes are consistent with amplitudes of the isophotes observed in elliptical galaxies.Indeed, our best-fitting values fall within the central 1-2  of the distributions found by Hao et al. (2006), even though we use Gaussian priors that allow for more extreme values.Interestingly, despite isophote observations indicating that third-order multipoles may on average be weaker than fourth-order (Hao et al. 2006), our findings indicate that they can still have significant impacts on flux ratios.We repeated our analysis with much tighter priors on the multipole amplitudes, Gaussians with widths 0.35 per cent instead of 1 per cent, and found that the cross and fold simulated lens systems can still be fit with  2 < 2.  The best-fitting realizations for the cusp and satellite systems both had  2 < 4, much better than the best-fitting EPL models.Even with priors on the multipole amplitudes that are comparable to, and tighter in the case of  4 , those suggested from isophotes, multipoles can have a significant impact on the flux ratios.
Still, it is currently unclear how closely galaxy mass distributions resemble their observed isophotes.Recently, Stacey et al. (2024) analysed a sample of three strong gravitational lens systems and inferred multipole amplitudes in the total mass density distribution of up to ∼0.01 for both third and fourth orders.They found that in some cases the multipole amplitudes in the isophotes of the lens galaxies differ from those inferred in the total mass density distribution by 0.01 or more.Modelling very long baseline interferometric (VLBI) observations of a lensed radio jet, Powell et al. (2022) found the thirdorder multipole to have an amplitude that is roughly twice as large as that of the fourth order, in disagreement with expectations from ob-servations of isophotes in elliptical galaxies.Whether this is related to a genuine discrepancy between the light and mass distributions or to the fact that third-order multipoles in the mass density are more degenerate with low-mass haloes, at least for resolved lensing observations (O'Riordan & Vegetti 2023), remains to be determined.As already pointed out by Stacey et al. (2024), the environment within which the lens galaxies reside may also play a role in observed differences between the isophotes and the iso-density contours.In general, our lack of knowledge about how well the projected mass density profile at the Einstein radius is traced by light limits our capability to use isophote distributions as a strong prior for the properties of multipoles in the mass distribution of gravitational lens galaxies.

Flux ratio
Figure 9. Flux ratios as a function of source size for fixed EPL+CDM (circles) and EPL+M34 (crosses) models that predict the same image positions and point-source flux ratios.Source size denotes the standard deviation of a Gaussian brightness distribution.Note that the physical extents of the radio-emitting, mid-infrared-emitting and warm dust regions of quasars are typically expected to be ≲ 10pc.

Considerations on the source size
To understand the role played by the source on the relative effect of low-mass haloes and multipoles, we quantify the level of flux-ratio anomaly induced by the two contributions, i.e., low-mass haloes and multipoles, independently as a function of the source size (see Figure 9).We select two models, one perturbed by low-mass haloes and the other by multipoles, which result in the same image positions and flux ratios from a lensed point source.As we increase the source size from a point to 40 pc, we find that the flux ratios perturbed by low-mass haloes approach their unperturbed values (i.e. the values expected from an EPL-only mass distribution).Even when the source is only 10 pc in size, there is a ∼10 percent change of two of the flux ratios from their point-source values.On the other hand, the flux ratios perturbed by multipoles are hardly affected by increasing the source size up to 20 pc (i.e.0-3 per cent).While in both cases larger source sizes produce smaller flux ratio anomalies, the effect is much smaller for the multipoles than for low-mass haloes.As the multipoles are less sensitive to changes in the source size, the more extended the source, the lower the multipole amplitudes that are needed to reproduce the effect of a given halo population.We can conclude, therefore, that our assumption of a point source gives a conservative quantification of the degeneracy between low-mass haloes and multipoles.

Implication for dark matter inferences using flux ratios
Traditionally, flux-ratio constraints have been based on the concept that any anomaly that is seen with respect to the flux ratios predicted by the smooth EPL model (perhaps with the addition of a stellar disc; Hsueh et al. 2016Hsueh et al. , 2017Hsueh et al. , 2018)), is due to perturbations by clumpy dark matter.Thus, the underlying assumption is that the prevalence of observed flux-ratio anomalies can be used to infer the mass-function of low-mass dark matter haloes and, in turn, the properties of dark matter.If, instead, some unknown fraction of those anomalies are caused by multipole structure in the primary lensing galaxy, the constraining power of this method is diminished.
In this paper we have shown that, when presented with a simple set of observables, namely four image positions and four fluxes, it is generally possible to reproduce those observations with a model consisting of an EPL base model plus either (1) low-mass CDM haloes but no multipoles or (2) multipoles but no low-mass haloes.
The fluxes in real lens systems are likely to be affected by a combination of these two effects.Because we cannot identify from just the lensed image positions and fluxes whether observed anomalies are due to low-mass haloes or multipoles, any dark matter inference from real flux-ratio data will be affected by what weight is given to each of these two components, i.e., by the priors on the relative contributions of the perturbing haloes and multipole components.For example, neglecting the contribution of the multipole component or reducing its complexity, e.g., by not including third-order multipoles or by fixing the multipole position angle to match the major axis of the EPL, will lead to a bias in favour of colder dark matter models.Most extremely, an inference could be made under the assumption that all perturbations were due to low-mass haloes, as was done in Hsueh et al. (2020), in which case a CDM-like model, which predicts large numbers of low-mass haloes, is likely to be preferred.Similarly, an inference could be made assuming the maximum multipole contribution, in which case models that predict fewer low-mass haloes, such as WDM, are likely to be preferred.All but the most extreme flux-ratio anomalies seen in real lensed quasar systems are comparable to or smaller than the anomalies considered in this paper, indicating that our ability to fit the anomalies with multipole components is not due to selectively choosing easy systems to fit.Thus, one conclusion to be drawn from our results is that, in the absence of other data that can inform our priors, lensed quasar flux ratios alone are ineffective at placing constraints on dark matter models.
In a follow-up paper, we will quantify how the inclusion of multipole components in the lens mass distribution affect existing constraints on dark matter from the two recent investigations by Hsueh et al. (2020) and Gilman et al. (2020).Despite their different approach to multipoles, both works found nearly identical constraints on the WDM particle mass.However, neither work allowed for a general treatment of multipoles as presented in this paper.We expect that allowing more freedom in the angular structure of the lens galaxies will reduce the contribution from low-mass haloes to the observed flux ratios, and potentially weaken the constraints on dark matter.

Future prospects
As the prospects to derive meaningful constraints on the properties of dark matter with flux ratio only observations strongly depends on independent knowledge about the general mass density distribution of lens galaxies, it is fair to wonder where the required information could be obtained.
One option is deep observations that could allow one to identify and quantify the prevalence of multipoles in the light distribution of lens galaxies.However, it is currently unknown how well the total mass distribution follows that of the light.Though the best-fit multipole amplitudes that we have found in this paper are not in tension with those observed in galaxy isophotes by Hao et al. (2006), the current number of quadruply-imaged quasar observations is not sufficient to precisely quantify the distributions of inferred multipole parameters across all elliptical galaxies from lens modelling.
Another option is modelling high-resolution images of strongly lensed extended emission, which is expected to allow for the constraint of more complex macro model features including multipoles (e.g.Powell et al. 2022;Stacey et al. 2024).In principle, this information could be applied to the analysis of flux ratios when available (Gilman et al. 2024).However, high-resolution imaging of lensed extended emission is at the moment not available for the majority of the ∼ 15 lens systems in the current flux ratio sample.In addition, O'Riordan & Vegetti (2023) have shown that the degeneracy between multipoles and low-mass haloes persists in the modelling of lensed extended emission, leaving as reliably detectable only lowmass haloes which are located in projection close to the lensed images.Furthermore, Herle et al. (2023) have shown that automatic lens finding techniques may lead to the selection of gravitational lens galaxies that have different mass distributions when the sources are extended instead of unresolved.If confirmed, their result may place some doubts on the possibility of transferring knowledge from one type of lens system to another.Galan et al. (2022) have recently shown that it may be possible to distinguish between the lensing effects of multipoles and low-mass haloes using extended sources and a non-analytical description of the lensing potential, but this method has yet to be demonstrated on real data.If proven successful, in combination with the large number of gravitational lens systems expected from the James Webb Space Telescope, Euclid and the SKA (Nierenberg et al. 2023;Collett 2015;McKean et al. 2015), it may represent a path forward to obtain the necessary information on the mass density distribution of lens galaxies.Similarly, numerical simulations could provide useful insights on the mass and light distribution of (lens) galaxies (e.g.Gao & White 2006;Łokas 2022;Despali et al. 2022;Han et al. 2023).
To conclude, our results as well as previous similar works demonstrate that the analysis of flux-ratio anomalies, in its current form, cannot directly constrain the properties of the low-mass dark matter halo population.This follows from the main issue identified in this work, that the priors imposed on multipole parameters play a dominant role in the inference process.Though our focus is on CDM, we expect our results to be applicable to other dark matter models, but potentially with some differences.We expect, for example, the degeneracy here identified to be stronger for warm (WDM) and weaker for self-interacting (SIDM) dark matter.In the former case, the amplitude of the multipoles should not be significantly different than in CDM, but as low-mass haloes are less numerous and concentrated they result in weaker flux ratio anomalies.In SIDM models, galaxies are predicted to be rounder (e.g.Brinckmann et al. 2018;however, see Despali et al. 2022), while low-mass haloes, which experience core-collapse, have a stronger lensing effect (Gilman et al. 2023).We will quantify the degeneracy between multipoles and low-mass haloes in different dark-matter models in a follow-up publication.
Figure1.The range of flux-ratios as a function of third-order multipole parameters  3 (left) and  3 (middle).Only flux ratios from models with position error less than 1 for each image are included.The horizontal lines and shaded regions represent the flux-ratios and uncertainties in the simulated data, is generated by an EPL+CDM model.The filled dots show the flux-ratios of each image for each multipole parameter value.Open circles represent the flux-ratios from the EPL+MP model that most closely matches the simulated flux ratios.The solid, dashed and dotted vertical black lines in the left column indicate the mean, 1 and 2 values, respectively, of the third-order multipole amplitude distributions from elliptical galaxy isophotes(Hao et al. 2006).The solid black line in the middle column indicates the orientation angle of the EPL major axis in each simulated lens system.The rightmost column shows the image configuration for each lens system, and the color of each image corresponds to data for its respective flux ratio in the other columns.The black point marks the image used in the denominator of the flux ratios.

Figure 2 .
Figure 2. Same description as Figure 1, except with fourth-order multipoles.

Figure 4 .
Figure 4. Corner plots from simultaneous MCMC sampling of EPL, shear and both multipole order parameters on the cross lens system.MCMC was conducted in the sine/cosine amplitude basis (left), and we also present the same data transformed into the overall amplitudes and angles (right) for ease of interpretation.Only multipole parameters are displayed, but all EPL and shear parameters, along with the source position, were allowed to vary freely.The red-dashed contours and histograms show the Gaussian prior distributions used for the   and   parameters.

Figure 5 .Figure 6 .
Figure 5. Same description as Figure 4, except for the fold system.

Figure 7 .
Figure 7. Same description as Figure 4, except for the satellite system.

Table 2 .
Unitless constants used in equations 2 and 3.

Table 3 .
Image positions and flux ratios for each of the simulated lens systems along with the fractional degree of flux perturbation by low-mass haloes compared to the best-fitting EPL and shear-only fiducial model.All position units are in arcsec.These are generated using an EPL+shear macro-model in addition to a population of perturbing CDM subhaloes and field haloes.No multipoles are present in any of the simulated systems.Positions are given with respect to the observation centre, and all dRA and dDec uncertainties are 0.005.Flux ratio uncertainties are listed in parentheses next to each flux ratio.