Radial acceleration relation of galaxies with joint kinematic and weak-lensing data

We combine kinematic and gravitational lensing data to construct the Radial Acceleration Relation (RAR) of galaxies over a large dynamic range. We improve on previous weak-lensing studies in two ways. First, we compute stellar masses using the same stellar population model as for the kinematic data. Second, we introduce a new method for converting excess surface density profiles to radial accelerations. This method is based on a new deprojection formula which is exact, computationally efficient, and gives smaller systematic uncertainties than previous methods. We find that the RAR inferred from weak-lensing data smoothly continues that inferred from kinematic data by about 2.5 dex in acceleration. Contrary to previous studies, we find that early- and late-type galaxies lie on the same joint RAR when a sufficiently strict isolation criterion is adopted and their stellar and gas masses are estimated consistently with the kinematic RAR.


Introduction
A well-known property of spiral galaxies is that their rotation curves are approximately flat at large radii [2,3].This flatness is, however, far from the only regularity contained in the dynamical properties of spiral galaxies.Indeed, a lot of additional regularities are captured by scaling relations like the Baryonic Tully-Fisher Relation [BTFR, [4][5][6][7][8], the Central Density Relation [CDR, 9,10], and the Radial Acceleration Relation [RAR, [11][12][13][14][15][16].The BTFR links the baryonic mass M b (stars and gas) of a galaxy to its flat rotation speed 1 .The CDR links the baryonic surface density at the galaxy center (R → 0) with the dynamical surface density inferred from the inner steepness of the rotation curve.The RAR links the observed centripetal acceleration at each radius (g obs = V 2 c /R) with that predicted from the observed distribution of stars and gas assuming Newtonian gravity (g bar = |∂Φ bar /∂R|).In particular, one finds g obs ≈ g bar at large accelerations, g bar ≫ a 0 , and g obs ≈ √ a 0 g bar at small accelerations, g bar ≪ a 0 .The acceleration scale a 0 occurs in all three scaling relations (BTFR, CDR, RAR).It has a universal value of about 10 −10 m/s 2 despite playing a different physical role in each relation.Scaling relations like the RAR are commonly measured using radio interferometry of the 21 cm spin-flip transition of atomic hydrogen.Such observations probe over many tens of kpc [21], sometimes reaching 100 kpc [22,23], without revealing any credible deviation from a universal RAR.Another probe is offered by weak gravitational lensing [e.g., [24][25][26].Recently, it was shown in Ref. [1] that lensing observations probe to at least 300 kpc and perhaps to Mpc scales.
Here, we rederive the RAR from the KiDS DR4 weak-lensing data [27][28][29] that was also used in Ref. [1], improving key aspects of the analysis.In particular, we estimate stellar masses using the same stellar population synthesis model [30] that was used for previous determinations of the RAR [15,21].We further introduce a new method for converting weak-lensing data to radial accelerations.This method is based on a new, exact deprojection formula that relates excess surface densities and radial accelerations.
We first introduce our new deprojection formula and the corresponding method for obtaining radial accelerations from weak-lensing data in Sec. 2. In Sec. 3 and Sec. 4, we then describe, respectively, the data set and baryonic mass estimates we use.We investigate the resulting radial acceleration in Sec. 5.After a brief discussion in Sec.6 we conclude in Sec. 7.

From excess surface density to acceleration
Weak-lensing observations measure the distortion of source images due to massive objects, called lenses, along the line-of-sight.Here, we consider galaxy-galaxy lensing where both the sources and lenses are galaxies.The goal is to measure the radial accelerations g obs around the lens galaxies.The main ingredients for this approach are the ellipticities of the source galaxies.These contain contributions from both intrinsic ellipticities and distortions induced by gravitational lensing.To isolate the effect of gravitational lensing, one averages over many source-lens pairs, after which the intrinsic ellipticities average out.In this way, one obtains an averaged tangential shear profile γ t (R) as a function of the projected distance from the lens R, which is related to the so-called excess-surface density (ESD) ∆Σ(R) [31,32].From such ESD profiles, one can then obtain averaged gravitational accelerations.
In this section, we introduce a new method for obtaining the radial acceleration g obs of an individual galaxy given its ESD profile ∆Σ.Following Ref. [1], we assume General Relativity (GR) and the usual thin-lens approximation [33], allowing for both baryonic and dark matter.Nevertheless, our results have a straightforward interpretation also in some modified gravity models; namely in models where lensing works as in GR, just with the lenses having a specific total mass M (r) (see below).Examples are discussed in Refs.[25,34,35].
The ESD profile ∆Σ(R) of an individual galaxy is defined in terms of its surface density Σ(R), Galactic scaling relations like the RAR are based on accelerations.Thus, to test these relations using weak-lensing data, we need to convert ESD profiles to accelerations.Here, we assume spherical symmetry which is a reasonable approximation at the large radii involved (up to a few Mpc, see below).The relevant acceleration is then the radial acceleration g obs (r) = GM (r)/r2 , where M (r) is the cumulative 3D mass calculated from the total density profile ρ(r), including both baryons and dark matter, and r is the spherical radius.In Ref. [1], two methods for converting ESD profiles to accelerations were proposed.The first is called the "SIS" method which assumes that the total density profile ρ of the galaxy is a singular isothermal sphere (SIS).With this assumption, one obtains the simple relation (2. 2) The second method is called the "PPL" method and is more elaborate.One fits a piecewise power-law density profile ρ to the observed ESD profile and, from the density profile obtained in this way, one finds the acceleration by integration.The advantage of the PPL method is that it does not assume a specific density profile.It does, however, come with a much higher computational cost and is more complicated to implement.
Here, we introduce a new method for converting ESD profiles to accelerations.We use a new, exact deprojection formula which relates the ESD profile ∆Σ and the radial acceleration g obs .We make no assumptions about the density profile ρ except that it is spherically symmetric and asymptotically falls off faster than 1/r with the spherical radius 2 r.Despite this generality, the method is simple to implement and computationally efficient.The basic relation is derived in Appendix A and reads It is instructive to compare this new relation to the "SIS" relation Eq. (2.2).We see that the SIS method obtains the radial acceleration g obs at a radius R from the ESD profile at the same radius R. Similarly, Eq. (2.3) places most weight on radii near R, but also averages over radii larger than R.
In particular, the integral in Eq. (2.3) involves ∆Σ(R) at arbitrarily large radii R.This is entirely expected: The ESD ∆Σ is defined in terms of the surface density Σ which is an integral over the density ρ that extends to infinity.In contrast, g obs (r) does not depend on ρ at arbitrarily large radii.Thus, getting g obs from ∆Σ(R) requires information about ∆Σ(R) out to infinity.
In practice, observations of the ESD profile ∆Σ extend only to a certain maximum radius R max .So we need to make a choice about what to do at radii larger than R max .Here, we assume that, at radii larger than R max , the ESD profile ∆Σ behaves as a singular isothermal sphere, i.e. ∆Σ ∝ 1/R, Other choices are possible.But without appealing to a specific theoretical model of how the density ρ should typically look like in our universe, we have no good reason to choose one over the other.Thus, there is a systematic uncertainty corresponding to how ∆Σ behaves3 at radii larger than R max .We quantify it by calculating g obs using two opposite, extreme assumptions and then taking the difference between these two versions of g obs .Here, we consider the two assumptions that ∆Σ drops to zero at the last data point and that ∆Σ keeps the value it has at R max up to infinity, That is, written symbolically, we adopt the following systematic error due to this uncertainty at large radii, In practice, this systematic error is often negligible.It only becomes important close to R max .This is for two reasons.First, the integrand ∆Σ in the integral Eq. ( 2.3) drops off with radius.Second, irrespective of the integrand, the integral Eq. ( 2.3) gives more weight to radii close to R and only relatively little weight to arbitrarily large radii.
Another systematic uncertainty arises because, in practice, we know the ESD profile ∆Σ(R) only in discrete radial bins (Appendix B.2), while the integral in the deprojection formula Eq. (2.3) is continuous.So we have to interpolate between the discrete radial bins.This also applies to the integral Eq. (2.8) for the statistical errors we introduce below.Here, for simplicity, we linearly interpolate between the discrete radial bins.To estimate the systematic uncertainty associated with this interpolation, we calculate g obs with both linear and quadratic interpolation and take the difference as an additional systematic error.Symbolically, Below, we will see that these systematic uncertainties are typically both small compared to the statistical uncertainties (except close to the last data point at R max ).This should be contrasted with the systematic uncertainties of the SIS and PPL methods.For them, a systematic uncertainty of 0.05 dex on g obs was estimated in Ref. [1]. 4 This is comparable to the statistical uncertainties found in Ref. [1].
For the statistical errors, we assume that the measurements of ∆Σ at different radii are independent.We further discuss this assumption in Appendix C.1.The statistical errors on g obs are then given by (see Eq. (2.3)), where σ ∆Σ denotes the statistical error on ∆Σ.This formula requires the statistical uncertainties of the ESD profile ∆Σ at radii beyond the last data point.However, these uncertainties are systematic ones that we treat separately.Thus, beyond the last radial bin R max , we adopt More specifically, we set σ ∆Σ to zero beyond the last bin edge, not beyond the last bin center.Otherwise, the θ integral would give exactly zero already at the last bin center.Between the last bin center and the last bin edge we interpolate linearly.So far we have considered ESD profiles and radial accelerations as a function of radius R. Below, we will mainly be interested in the RAR that requires g obs as a function of the Newtonian baryonic acceleration g bar ≡ GM b (R)/R 2 .For an individual galaxy, obtaining g obs (g bar ) from g obs (R) is straightforward, at least as long there is a one-to-one mapping between radii R and baryonic accelerations g bar .Indeed, instead of g obs (R) one can simply use g obs (R(g bar )) where the function R(g bar ) maps values of g bar to values of R for that galaxy.
Below, we will not be interested in the lensing signal of an individual galaxy but in the stacked signal of a large number of galaxies.In Appendix B, we explain how to adapt the method discussed above to this practically relevant case of stacking a large number of galaxies.We also explain in Appendix B how we obtain the required ESD profiles from observational data.
There are some subtleties when adapting our method to work with stacked data.Basically, there is a choice between first deprojecting individual galaxies using Eq.(2.3) and then stacking (which is slightly more susceptible to observational systematics) versus first stacking individual galaxies and then deprojecting (in which case the deprojection is not exact).In practice, both choices give almost identical results (see Appendix C.2) despite their different systematics.Thus, we trust our results despite these subtleties.The results shown below are obtained using the first method, i.e. we first deproject and then stack.We discuss this in detail in Appendix B.
We note that the deprojection formula Eq. ( 2.3) can also be used to infer cumulative 3D mass profiles M (< R) from weak-lensing observations.Indeed, in spherical symmetry we have M (< R) = (R 2 /G) g obs (R).Thus, Eq. ( 2.3) may also be useful when one is not directly interested in radial accelerations.

Data
We adopt H 0 = 73 km s −1 Mpc −1 for consistency with the RAR derived from kinematic measurements [15] and the Hubble constant measured with the BTFR [8].Below we sometimes use the notation h 70 , defined as H 0 /(70 km s −1 Mpc −1 ).Our choice corresponds to h 70 = 73/70.We adjust our masses and radii accordingly.
We use mostly the same data as Ref. [1].Specifically, we use source galaxies from the KiDS-1000 SOM-gold catalog [27,28,38,39] and lens galaxies from the KiDS-bright sample [29].We use the best-fit photometric redshifts z B,s for the sources s, and the redshifts z ANN,l from the machine-learning method "ANNz2" [40] for the lenses l.We include only lenses with 'masked = 0' and with redshift 0.1 < z ANN,l < 0.5.
In addition, we restrict the lens galaxies to be isolated using the criterion introduced by Ref. [1].That is, for each lens, we enforce a lower bound R isol on the 3D distance to the closest neighboring galaxy with at least 10% of its stellar mass.We use R isol = 4 Mpc/h 70 unless stated otherwise.This is larger than the value R isol = 3 Mpc/h 70 used by Ref. [1].We explore the effects of making this criterion stricter or laxer in Sec.5.2 below.We use the stellar masses of Ref. [1] for the isolation criterion.We do this for simplicity and to have a direct comparison to Ref. [1].Adopting our own mass estimates from Sec. 4 instead would not affect our conclusions, as we explain there.
Note that the isolation criterion is of great practical importance.Non-isolated galaxies pick up signals from surrounding structure which violates the assumption of spherical symmetry that we make; see also the discussion of the so-called two-halo term in Ref. [1].Even if we ignore this assumption of spherical symmetry, our analysis still requires isolated galaxies because we are not interested in the structure surrounding lens galaxies but in the intrinsic properties of the lens galaxies themselves.
Since the isolation criterion operates on 3D distances it depends on the redshifts of the KiDS-bright lenses.There are no spectroscopic redshifts available for this sample, so we use the photometry-based redshifts z ANN,l which, as all photometric redshifts, have non-negligible uncertainties.We discuss this more in Sec.6 below.An alternative would be to base the isolation criterion on projected distances R isol,proj instead of 3D distances R isol .This would give a stricter criterion which does not rely on photometric redshifts.We nevertheless prefer the isolation criterion based on 3D distances for two reasons.First, this allows for a direct comparison to Ref. [1].Second, the sample size becomes very small even for moderate values of R isol,proj .For example, for R isol,proj = 1 Mpc/h 70 , there are only 196 lenses left.Following Ref. [1], we impose an upper limit on the stellar masses of the lens galaxies because massive galaxies have significantly more satellites [see also Sect.2.2.3 of Ref. 26] and because our adopted light to stellar-mass conversions [30] (see Sec. 4) may require modest revision for the most massive ETGs [41] (we discuss this more in Sec.5.3).We adopt log 10 M * /M ⊙ < 11.1.For this, we use our stellar mass estimates from Sec. 4. This leaves 106,843 lenses in our sample which is less than the 259,383 used by Ref. [1].The biggest difference comes from the stricter isolation criterion we apply.Since our stellar masses are, on average, larger than those of Ref. [1], see Sec. 4, we use a slightly larger cut-off for log 10 M * /M ⊙ , namely of 11.1 instead of 11.0.
For the sources, we restrict the SOM-gold catalog to 'SG_FLAG = 1', 'SG2DPHOT = 0', and 'CLASS_STAR < 0.5' in order to remove stars from the sample [1,27].We also apply the quality cut 'IMAFLAGS_ISO = 0' and the recommended mask 'MASK & 28668 = 0' [1,27].In addition, when estimating the tangential shear around a lens, we follow Ref. [1] and consider only sources s that are not too close to a given lens l in redshift to avoid including sources in front of the lens l.Specifically, we consider only sources with z B,s > z ANN,l + 0.2.
For the critical surface density (see Appendix B.2 for where this is needed), we use the method of Refs.[1,42] to take into account uncertainties in the redshifts.That is, we adopt Here, D(z l ) and D(z s ) are the angular diameter distances to the lens and source, respectively, and D(z l , z s ) is the angular diameter distance between the lens and the source.Following Ref. [1], we calculate these distances assuming a flat FLRW cosmology with Ω m = 0.2793.Further, p z ANN,l (z) is a normal distribution centered on the lens redshift z ANN,l with standard deviation σ = 0.02 • (1 + z ANN,l ).The function n z l ,z B,s (z) is determined as follows.For a given value of the integration variable z l and a given best-fit photometric source redshift z B,s , we first find out in which of the five tomographic bins from Ref. [39] the z B,s value belongs and get the corresponding redshift distribution function from Ref. [39].Then we normalize this distribution to unity in the interval [z l , ∞).The stellar mass histogram assuming our stellar masses for the whole sample (gray), for early type galaxies (ETGs in red), for late type galaxies (LTGs in blue), and for the overlap of ETGs and LTGs (green).

Masses and mass distributions
In order to properly compare kinematic and lensing data, it is important that the stellar mass informing g bar is consistent for both samples.We do not observe stellar mass directly, but rather the apparent magnitude of each galaxy.This is converted to a luminosity for the assumed distance scale, and then to a stellar mass using a mass-to-light ratio indicated by a stellar population model informed by the observed colors.Stellar population models have converged so that most give similar results, but agreement is never perfect and systematic differences still occur at the factor of two level [43].
To insure a uniform stellar mass scale, we have reanalyzed the KiDS data used here with the same stellar population synthesis model [30] used for the kinematic data [15].Comparing our results to those of Ref. [1], we find good agreement overall, albeit with a slight morphology dependence.Following Ref. [1], we define early and late type galaxies (ETGs and LTGs) based on the color 5 split u − r ≷ 2.5.Indeed, the agreement is uniform enough that the bulk of the sample can be described as a simple multiplicative scaling factor, which we parameterize by defining the ratio Q such that 6 M * = QM KiDS * .The agreement in the average mass-to-light ratio is excellent for LTGs: Q LTG = 1.0.For ETGs, we find a modest offset: Q ETG = 1.4.This goes in the sense that we assess ETGs to be slightly more massive than do Ref.[1].We apply these average correction factors to each individual galaxy depending on which type they are.We show the resulting stellar masses in Fig. 1.
An offset like the one we find for ETGs is to be expected when comparing different stellar population synthesis models.Indeed, the exact same type-dependent factor is found in Ref. [1] when comparing stellar masses for ETGs in their KiDS and GAMA7 samples.The building block of any composite stellar population model is a single stellar population (SSP) of a given age: we use the SSPs from Ref. [45], while the KiDS team used the older SSPs from Ref. [46].In addition, the composite stellar population models we employ [30] assume physically-motivated differences in the star-formation and chemical-enrichment histories of early and late types: ETGs form most of their stars early and quickly reach high metallicities, whereas LTGs have a more continuous star formation and a more gradual chemical enrichment [see also 41,47].Unfortunately, Ref. [1] does not provide details about the star-formation and chemical-enrichment histories of their stellar population models, so we cannot discern what is the dominant cause of disparity in the stellar masses of ETGs.Two possible culprits are the treatment of AGB stars, for which we use an empirical calibration [30], and metallicity, for which we take care to use a distribution of stellar metallicities that is consistent with chemical evolution models rather than the usual approximation that all stars have the same metallicity.The metallicity plays an important role in the color-magnitude relation of ETGs [48] that may induce larger mass-to-light ratios in brighter ETGs [41].Though the difference is modest (∼40%), it is important to whether or not ETGs follow the same RAR as LTGs.
We follow Ref. [1] in using scaling relations to account for the baryonic mass in gas.Massive ETGs, such as those considered here (Fig. 1), are known to have X-ray-emitting coronae of hot ionized gas [e.g., 49].Thus, we add hot gas to the total baryonic mass according to the scaling relation [50], Star-forming LTGs have a non-negligible interstellar medium of atomic and molecular gas, so we add cold gas according to the scaling relation where The first term in Eq. (4.2) represents atomic gas according to the scaling relation from Ref. [21] with M * = 0.5L [3.6] .The second term takes into account molecular gas.Eq. ( 4.3) accounts for the variation of the hydrogen fraction X as metallicity varies with stellar mass [51].This accounting of the cold gas is practically indistinguishable from the relation adopted by Ref. [1] over the relevant mass range.Thus, our baryonic mass estimate for LTGs is close to that of Ref. [1].For ETGs, we use the hot gas mass estimate Eq. (4.1), while in Ref. [1] the same cold gas mass estimate as for LTGs is used.At large stellar masses, our hot gas mass estimate is larger than the cold gas mass estimate of Ref. [1].Thus, since ETGs tend to have large stellar masses (Fig. 1), our baryonic mass estimate for most ETGs is larger than that of Ref. [1].On average, our ETG baryonic masses are larger by a factor of about 1.7, which includes the effects of our different stellar M * /L, our different gas mass estimate, and our different choice of H 0 .Our results below mainly depend on the total baryonic mass, not on how this mass is distributed between gas and stars.Thus, we would expect to find similar results with Q ETG = 1.7 and the gas mass estimate of Ref. [1].
For simplicity, we treat both the stellar mass and the gas mass as point masses.This assumption is sensible because the bulk of the observed baryonic mass is typically contained within tens of kpc, which is comparable to the smallest radial bins we consider while the largest exceed 3 Mpc.Atomic gas disks in LTGs and hot gas halos in ETGs are sometimes observed to extend out to ∼ 100 kpc, but both components typically give a minor gravitational contribution out to these radii.Hot gas around galaxies may be more extended than that, forming the so-called circumgalactic medium, but its amount and distribution are poorly constrained, so they will be neglected in our work.In Appendix C.6, we explicitly show that modeling the hot gas from Eq. (4.1) as an extended distribution instead of a point mass does not significantly change our results.
As discussed in Sec. 3, we reuse the stellar masses of Ref. [1] for the isolation criterion defined there.We expect that using our own stellar masses instead would not significantly affect our results, because the isolation criterion considers only ratios of stellar masses which are mostly unchanged by the simple rescaling we adopt here.The only reason our rescaled stellar masses would make a difference at all is that we rescale masses of ETGs and LTGs differently.Thus, when an LTG lens has an ETG neighbor or vice versa, the isolation criterion would change when adopting our stellar masses.In practice, however, this does not change any of our conclusions for the following reason.Our stellar mass estimates agree with those from Ref. [1] for LTGs and are a bit higher than those from Ref. [1] for ETGs.Thus, for ETGs we are actually adopting a stricter isolation criterion by sticking with the stellar masses from Ref. [1].For LTGs, our isolation criterion would be stricter if we used our own stellar masses instead.But in Sec.5.2, we will see that our LTG sample is sufficiently isolated anyway, so this is not a problem.We also note that the precise quantitative meaning of R isol is not important to any of our results.It is only important that by adjusting R isol we can make the isolation criterion more or less strict.

The radial acceleration relation
Fig. 2 and Table 1 show the weak-lensing RAR obtained using our new method for converting excess surface densities to radial accelerations from Sec. 2 and our stellar and gas mass estimates from Sec. 4. The weak-lensing RAR smoothly continues the RAR obtained from kinematics [15].We also show the fit function used in Ref. [15], extrapolated to the smaller accelerations probed by weak lensing, adopting a 0 = 1.24 • 10 −10 m/s 2 [52].This extrapolated fit function matches the weak lensing data well, except at the last few data points at g bar ≲ 10 −14 m/s 2 where systematic uncertainties are large.We further discuss these last few data points in Sec. 6.In Appendix D.1, we compare the RAR obtained using our new method to that obtained using the SIS method used in Ref. [1].
The RAR from weak lensing data (yellow diamonds), derived using our new deprojection formula Eq. ( 2.3) and assuming our stellar and gas mass estimates described in Sec. 4. The binned kinematic RAR from Ref. [15] is shown as gray circles.The error bars indicate the statistical uncertainty on the mean in each bin (not the galaxy-by-galaxy variation).The colored band indicates the systematic uncertainties of the lensing result from interpolation between the discrete bins and from extrapolation beyond the last data point, see Sec. 2. The gray band indicates an additional systematic uncertainty of about 0.2 dex in stellar mass [1].For simplicity, we translate this 0.2 dex uncertainty into a ∼ 0.1 dex uncertainty on g obs using the fact that g obs scales as √ g bar .The dashed gray line shows the fitting function Eq. (5.1) from Ref. [15].The shaded region at g bar < 10 −13 m/s 2 indicates where the isolation criterion may be less reliable according to the estimate from Ref. [1].Our results from Sec. 5.2 suggest that LTGs may be sufficiently isolated down to g bar ≈ 10 −14 m/s 2 .We shade this region where LTGs may still be reliable in a lighter color.
The reason for the large systematic errors below g bar ≈ 10 −14 m/s 2 is the need to extrapolate beyond the last ESD data point, as explained in Sec. 2. Thus, in principle, we could likely get formally smaller systematic errors at g bar ≲ 10 −14 m/s 2 by simply extending our analysis to values of g bar below our current limit of g bar = 10 −15 m/s 2 .However, in practice, this is not very useful for two reasons.First, observational systematics become more than just a small correction beyond a few Mpc, see Appendix C.2.If we extended our analysis to, say, g bar = 10 −16 m/s 2 , a typical lens galaxy with M b = 5 • 10 10 M ⊙ would be probed up to R ∼ 8 Mpc.As described in Appendix C.1, we deal with such systematics by subtraction.That procedure may, however, be suspect when these systematics are no longer just a small correction.Second, the isolation criterion becomes less reliable at large radii.Indeed, we have verified that the last few data points shown in Fig. 2 are sensitive to the precise choice of R isol (see also Sec. 5.2).At even smaller values of g bar , such effects would become even more important.Thus, the systematic uncertainties shown in Fig. 2 do indicate where systematic uncertainties become important, at least qualitatively.
Table 1.The weak-lensing radial acceleration relation from Fig. 2. The kinematic data shown in Fig. 2 is available from Ref. [15], see also http://astroweb.cwru.edu/SPARC/.The accelerations g bar and g obs are measured in m/s 2 .Uncertainties on g obs are converted to log 10 space using linear error propagation.The systematic uncertainties listed here come from interpolating between discrete radial bins and from extrapolating beyond the last radial bin (see Sec. 2, colored band in Fig. 2).For brevity, the additional, fixed 0.1 dex systematic uncertainties on g obs from uncertainties in the stellar masses (shown as a gray band in Fig. 2) are not listed here.The shaded region in Fig. 2 indicates where the isolation criterion may be less reliable.This region corresponds to g bar < 10 −13 m/s 2 where the isolation criterion from Ref. [1] may no longer be reliable according to their estimate.We conservatively adapt this estimate despite using a stricter isolation criterion, namely R isol = 4 Mpc/h 70 instead of R isol = 3 Mpc/h 70 .Our results from Sec. 5.2 indicate that LTGs are sufficiently isolated at even smaller g bar , down to g bar ≈ 10 −14 m/s 2 .It's only the ETGs that we are concerned about already at g bar ≈ 10 −13 m/s 2 .This is why we shade the region above 10 −14 m/s 2 , where LTGs may still be reliable, in a lighter color.
Our method for converting ESD profiles to radial accelerations is based on the deprojection formula Eq. (2.3) which involves arbitrarily large radii.Thus, if the isolation criterion fails at large radii (small g bar ), this may in principle affect even the radial accelerations we infer at small radii (large g bar ).However, as we demonstrate in Appendix C.4, this is not a problem in practice.We find almost identical results if we artificially cut off the ESD profiles that enter the integral formula Eq. ( 2.3) at g bar = 10 −14 m/s 2 .

Similar RAR for early-and late-type galaxies
In Ref. [1], it is found that the RAR of ETGs deviates from that of LTGs.In particular, according to Ref. [1], LTGs follow the fit function Eq. (5.1) of Ref. [15] even at extremely small accelerations g bar , while ETGs deviate from it.Here, we argue that -with our stricter isolation criterion and our mass estimates from Sec. 4 -weak-lensing data does not indicate The acceleration implied by weak lensing for ETGs (red triangles) and LTGs (blue diamonds), relative to that of the fit function Eq. (5.1) from Ref. [15], here denoted by g RAR .The kinematic data from Ref. [15] is shown as gray circles.Error bars and bands are as in Fig. 2, except we do not show the stellar mass systematic uncertainty for clarity.In contrast to Ref. [1], we find good agreement between ETGs and LTGs down to about g bar = 10 −14 m/s 2 .The reason is that we use a stricter isolation criterion, R isol = 4 Mpc/h 70 , and larger baryonic masses for ETGs (see Sec. 4).
a difference in the RAR of ETGs and LTGs. 8This fits with kinematic data which indicates that ETGs and LTGs follow the same RAR at accelerations g bar above 10 −12 m/s 2 [15,16].
Specifically, we find no significant difference between ETGs and LTGs down to g bar ≈ 10 −14 m/s 2 , see Fig. 3.As we argue below, the difference to Ref. [1] is in the baryonic masses of ETGs (Sec.4) and the isolation criterion.Following Ref. [1], we here restrict the ETG and LTG subsamples to have the same stellar mass distribution (see the "overlap" shown in Fig. 1) by randomly removing an appropriate number of galaxies at each stellar mass value. 9n Appendix D.2 we show that we can reproduce the result of Ref. [1] if we adopt their baryonic masses and isolation criterion, i.e. if we adopt R isol = 3 Mpc/h 70 , Q ETG = 1.0, and their gas mass estimate instead of our preferred values R isol = 4 Mpc/h 70 , Q ETG = 1.4,and our gas mass estimate (see Sec. 4).Thus, the difference between our result and that of Ref. [1] is to be found in these differences in the masses and isolation criterion and not in the difference in the method for converting ESD profiles to radial accelerations (see Sec. 2).
We further show in Appendix D.2 that using our baryonic masses while keeping the isolation criterion of Ref. [1], R isol = 3 Mpc/h 70 , helps reduce the difference between ETGs and LTGs, but does not eliminate it.Thus, both the isolation criterion and the baryonic masses of the lenses are important.In the rest of this section, we will look into this more quantitatively.We will show that ETGs and LTGs indeed agree reasonably well when using our default choices for the isolation criterion and the baryonic masses of the lenses.In addition, we will argue that LTGs are sufficiently isolated down to about g bar ≈ 10 −14 m/s 2 , i.e. further down than the validity limit of the isolation criterion from Ref. [1] according to their estimate, namely g bar ≈ 10 −13 m/s 2 .The same does not hold for ETGs, which are much more sensitive to the precise choice of isolation criterion.
To understand the effect of the isolation criterion, we try various values of R isol , see Fig. 4. Following Ref. [1] we quantify the difference between the lensing-inferred radial accelerations and those implied by the fit function Eq. (5.1) by first calculating a χ 2 value and then converting this into a number of σ.For simplicity, we treat different radii as independent, i.e. we leave out the small off-diagonal elements of the covariance matrix [1].Concretely, where i runs over the g bar bins, g obs,lensing is the lensing-inferred radial acceleration, and σ 2 g obs statistical is the statistical uncertainty of g obs,lensing .Appendix D.2 suggests that the difference between ETGs and LTGs is larger at relatively small accelerations and smaller at relatively large accelerations.To keep track of this, we separately consider large and small accelerations, small g bar : 10 −14 m/s 2 − 10 −13 m/s 2 , large g bar : 10 −13 m/s 2 − 10 −11 m/s 2 . (5.3) We do not consider data below g bar = 10 −14 m/s 2 because there the systematic uncertainties become important, while χ 2 takes into account only the statistical uncertainties.We can now discuss the effect of R isol shown in Fig. 4. We first note that LTGs are barely affected by changes in R isol .This holds for both small and large accelerations g bar .LTGs follow the fit function Eq. (5.1) quite well down to g bar ≈ 10 −14 m/s 2 .That there is no dependence on R isol indicates that the LTG sample is sufficiently isolated down to g bar ≈ 10 −14 m/s 2 , as we will discuss in more detail below.At even smaller g bar -in the very tail of our lensing RAR -systematic uncertainties become important and even the LTGs are sensitive to R isol , so one should be cautious about this tail.
The situation is different for ETGs.With R isol = 3 Mpc/h 70 , they deviate significantly from the fit function Eq. (5.1), especially at relatively small accelerations g bar , but this deviation becomes smaller as we increase R isol .That the trend with R isol is stronger for small accelerations fits with the fact that the isolation criterion is expected to be less reliable at large radii [1].In any case, this result shows that for ETGs one must be careful with the isolation criterion.A too small R isol may result in an artificially increased signal, likely due to the two-halo term [1].This is why our default choice is R isol = 4 Mpc/h 70 , i.e. a stricter isolation criterion compared to Ref. [1].
That ETGs are more sensitive to the isolation criterion than LTGs makes sense from the perspective of the morphology-density relation, i.e. considering that ETGs are more clustered than LTGs [53].That is, if the isolation criterion fails, we would expect it to fail earlier for ETGs.This is consistent with what we see here.
Compared to Ref. [1], we here consider a stricter isolation criterion, i.e. we make R isol larger than 3 Mpc/h 70 .As a result, the lens sample becomes smaller and the statistical uncertainties increase.Thus, one may worry that the difference between LTGs and ETGs becomes statistically insignificant simply because of the larger uncertainties.In Appendix C.5, we demonstrate that this is not the case; we find the same qualitative behavior if we artificially keep the error bars fixed at the values they have for R isol = 3 Mpc/h 70 .
Fig. 4 also shows that ETGs do not flatten-out as a function of R isol in the same way LTGs do.This is important because, if our results depend on R isol , this means we are not (just) measuring intrinsic properties of the lens galaxies but (also) something about their environment.In contrast, when the trend flattens out -as for the LTGs -we are plausibly measuring an intrinsic property of the galaxies. 10Thus, ETGs -in contrast to LTGs -may not be sufficiently isolated down to g bar ≈ 10 −14 m/s 2 .ETGs look more reliable down to 10 −13 m/s 2 , which matches where the isolation criterion from Ref. [1] is reliable according to their estimate.
It would be interesting to see if the trend of ETGs flattens out at R isol larger than what we show in Fig. 4, i.e. at R isol > 4 Mpc/h 70 .We have verified that -within the statistical uncertainties -it does.However, beyond R isol = 4 Mpc/h 70 , statistical fluctuations due to the small sample size take over if we artificially keep the error bars fixed at the values they have for R isol = 3 Mpc/h 70 .So at R isol > 4 Mpc/h 70 , we cannot be sure if the trend flattens out just due to the larger error bars (i.e. the argument from Appendix C.5 no longer works).
We note that, at the strictest isolation criterion shown in Fig. 4, R isol = 4 Mpc/h 70 , ETGs match the fit function Eq. (5.1) almost perfectly down to g bar ≈ 10 −14 m/s 2 .However, since we are not seeing the trend flatten out, this almost perfect match may be a statistical fluctuation.This leaves open the possibility that other M * /L values for ETGs, i.e. other values of Q ETG , may provide a better fit.We have tried various values and find that Q ETG = 1.8 works well. 11However, because we cannot be sure that ETGs are sufficiently isolated, these numbers should not be overinterpreted.
To sum up, our LTG sample seems to be sufficiently isolated to derive a robust lensing RAR down to g bar ≈ 10 −14 m/s 2 .Indeed, both ETGs and LTGs follow the fitting function Eq. (5.1) quite well down to g bar ≈ 10 −14 m/s 2 when we impose a sufficiently strict interpolation criterion and use our mass estimate from Sec. 4.However, the results for ETGs may not be reliable all the way down to g bar ≈ 10 −14 m/s 2 because these depend quite a bit on the details of the isolation criterion that one adopts.In any case, we find no significant difference between the RAR for LTGs and ETGs.5.1) down to g bar ≈ 10 −14 m/s 2 .This is especially true for larger baryonic masses.For smaller masses, there are larger deviations.A part of this is likely due to statistical fluctuations.Indeed, the sample is smaller for smaller masses, as indicated by the larger error bars in Fig. 5.Still, there may be a small systematic shift to larger accelerations g obs .

Mass bins
If there is such a systematic shift, a possible explanation is a failure of the isolation criterion: In contrast to larger galaxies, smaller galaxies can be affected even by a relatively small neighbor.This is exacerbated by the m = 20 magnitude limit of the KiDS bright sample [1,29] which means that many of the neighbors that affect small lenses are not detected.
The effect of the magnitude limit was discussed in Appendix A of Ref. [1] who estimated that the effect on their results is small.However, only the effect on the whole lens sample was considered in Ref. [1], not the effect on a small-mass subsample.Indeed, small-mass galaxies are only a small part of the whole lens sample so, overall, they have only a small effect.If one 10 Strictly speaking, one may find a flattening-out at 0 σ even with a non-isolated sample simply because the error bars increase with R isol .However, Appendix C.5 shows that we find the same qualitative trends with R isol even when we artificially keep the size of the error bars fixed (until fluctuations take over at large R isol , as we discuss below).
11 Alternatively, a similar result can be obtained by keeping the M * /L fixed and adjusting the hot gas mass estimate.This is because, as already mentioned in Sec. 4, our results are mainly sensitive to the total baryonic mass.

Figure 5.
The RAR implied by weak lensing for four baryonic mass bins with bin edges log 10 M b /M ⊙ = [9.0,10.5, 10.8, 11.1, 11.5] assuming our mass estimates from Sec. 4. The error bars and bands are as in Fig. 2. All mass bins generally agree with the RAR fitting function Eq. (5.1) from Ref. [15].In the smallest mass bin, the radial accelerations g obs tend to fall above this fitting function's prediction more than for the other mass bins.This may be because the isolation criterion is less reliable for small masses.The last data point in the highest-mass bin is not shown because the inferred g obs is negative there.splits the sample by mass, however, the small-mass bins may be affected by this systematic effect.Thus, we advise some caution when considering small-mass subsamples. 12e have explored lifting our upper cutoff in stellar mass and considered a fifth mass bin including lenses up to log 10 M b /M ⊙ = 11.9.We find that the overall shape of the lensing RAR remains the same but with a slight offset of g obs towards larger values.One possible reason is that the isolation criterion fails for such high-mass galaxies because they have more satellites.Indeed, this is the reason the stellar mass cut was originally introduced in Ref. [1].However, the offset we find is basically constant across all g bar values which is in conflict with the expectation that the isolation criterion should break down gradually towards larger radii.In contrast, a roughly constant offset is expected if we underestimate the M * /L for high-mass galaxies.The high-mass end of our sample is dominated by ETGs and the offset we are seeing in their RAR may be an indication of a mass-dependent M * /L, which is a possibility we did not consider when calculating our Q ETG .A value of Q ETG = 1.8 roughly compensates the offset in g obs , implying that our baseline stellar masses for the most massive ETGs might need to be increased by a factor 1.8/1.4= 1.29.The required change seems plausible in terms of stellar population models [41], especially given the possible need for a mass-dependent IMF [54][55][56], a topic beyond the scope of this paper.

Discussion
In Sec.5.1 we saw that the RAR shows signs of a downturn in the last few data points at g bar ≲ 10 −14 m/s 2 .Due to the systematic uncertainties discussed there, one should be careful not to read too much into these last few data points.Still, if the downturn is real, it would be an interesting subject of further study.
For example, in the context of particle dark matter, a downturn of the RAR is expected due to the finite size of the dark matter halo.In fact, such a downturn should happen already well before g bar = 10 −14 m/s 2 [e.g.57,58].Thus, the extended part of the weak-lensing RAR before the potential downturn is likely to provide strong constraints on particle dark matter models.Indeed, the density of an NFW tail falls off as 1/r 3 while, before the potential downturn, the weak-lensing RAR requires a 1/r 2 profile.
In the context of modified gravity, a downturn at large radii may be expected for different reasons, related to the lenses no longer being isolated.Most relevant here are models that reduce to Modified Newtonian Dynamics [MOND, [59][60][61] in the non-relativistic limit.Many of these have the attractive property that they naturally explain why the weak-lensing RAR follows the MOND-inspired fit function Eq. (5.1) so well [25,62].Importantly, many of these models also feature a so-called External Field Effect (EFE) which leads to a downturn at large radii where galaxies can no longer be treated as isolated [62,63]. 13This may explain the possible downturn of the lensing RAR at g bar ≲ 10 −14 m/s 2 .
Unfortunately, quantitatively testing the EFE is complicated.One reason is that the EFE works differently in different models of MOND.Another is that one would actually expect two competing effects in such models.One is the EFE which one expects to lead to a downturn.The other is a version of the so-called two-halo term which tends to lead to an upturn instead.Thus, one would have to quantitatively work out the net effect of these two competing effects.
A final complication is relevant even beyond testing modified gravity models.Namely, the isolation criterion we apply to the lens galaxies is ultimately based on photometric redshifts.These redshifts are estimated to have an uncertainty of σ z = 0.02 (1 + z) [1].Since the relevant redshifts z are in the range 0.1 to 0.5, these uncertainties translate into a distance uncertainty much larger than the value R isol = 4 Mpc/h 70 we use for our isolation criterion.This complicates any quantitative analysis of the EFE and other effects related to non-isolation, and is a reminder of the importance of accurate spectroscopic redshifts.
More generally, this means our lens sample is unlikely to be as isolated as the formal definition of R isol suggests (no other galaxies with at least 10% of a lenses' mass within R isol ).We nevertheless trust our results for the following reasons.First, the isolation criterion from Ref. [1] was validated against simulations and found to be trustworthy down to about g bar = 10 −13 m/s 2 .Since we use a stricter isolation criterion, our results should also be trustworthy down to at least g bar = 10 −13 m/s 2 .Second, we have seen in Sec.5.2 that, down to g bar ≈ 10 −14 m/s 2 , the RAR for LTGs is almost independent of the value of R isol .This indicates that LTGs are sufficiently isolated down to g bar ≈ 10 −14 m/s 2 .ETGs, in contrast, are quite sensitive to the precise value of R isol .This makes sense considering that ETGs are more clustered [53], but means we cannot trust our results for ETGs quite as far down in g bar .In particular, as discussed above, we trust our results for ETGs down to g bar ≈ 10 −13 m/s 2 and those for LTGs down to g bar ≈ 10 −14 m/s 2 .
At even smaller accelerations, g bar ≲ 10 −14 m/s 2 , both ETGs and LTGs are unlikely to be sufficiently isolated.This fits well with what we suggested above in the context of modified gravity models; namely that the possible downturn of the lensing RAR at g bar ≲ 10 −14 m/s 2 may be due to the isolation criterion breaking down, and not due to an intrinsic property of isolated galaxies.Indeed, there is a limit to how isolated galaxies can be due to the large-scale structure of the universe [50].Investigating this in more detail is left for future work, as is the analysis of larger spectroscopic surveys as data improve.

Conclusion
We have combined weak-lensing and kinematic data to construct the RAR over a large dynamic range in acceleration.We have estimated the lens galaxies' stellar and gas masses consistently with previous kinematic determinations of the RAR.We have further employed a new deprojection formula that converts excess surface densities to radial accelerations for spherically symmetric lenses.
We find that the RAR inferred from weak lensing smoothly continues the RAR from kinematic data by about 2.5 dex, implying a dark matter density profile that scales as 1/r 2 out to large radii.In contrast to previous studies we find no significant difference between the RAR for ETGs and LTGs.This is partly due to our somewhat larger baryonic masses for ETGs and partly due to the fact that we impose a stricter isolation criterion.
At the last few data points -where systematic uncertainties are important -we find hints of a downturn of the RAR.We speculate that this may be related to a failure of the isolation criterion rather than an intrinsic property of the lenses.We will further investigate this in future work.This leaves the second integral.Plugging in our result Eq. (A.4) gives where we introduced the short-hand notation By plugging the expressions Eq. (A.8) and Eq.(A.9) for I 1 and I 2 into Eq.(A.5) we obtain an expression for g obs in terms of ∆Σ.It remains to simplify this expression.
Our first simplification step is to reduce the double integral I 22 to a single integral.To facilitate this, we change the order of integration, where Θ denotes the Heaviside step function.The lower integration boundary of the r integral and the Θ function ensure that we always have r < R < R ′ .This implies we can set the lower integration boundary of the R ′ integral to r, (A.12) The x integral can be done by Mathematica [65], All integrals can now be combined into a single integral with integrand proportional to ∆Σ.We have, with the auxiliary function h(R) given by Some of these terms cancel and we end up with Our final simplification step is to substitute r R ≡ sin θ , (A.17) with θ in the interval [0, π/2] which finally gives B RAR from stacking

B.1 Stacking a large number of galaxies
In Sec. 2, we have considered individual galaxies.However, the weak-lensing signal for any individual lens galaxy is small.Thus, one usually considers the stacked signal from a large sample of lens galaxies.
The stacked ESD profile used by Ref. [1] can be written in the form where the sum runs over the lens galaxies l, the w l (R) are unnormalized weights, N (R) = l w l (R) is a normalization factor, and ∆Σ l (R) is an estimate for the ESD profile of the lens l at the projected radius R. We will explain how w l (R) and ∆Σ l (R) are related to actual weak-lensing data below in Appendix B.2.
But first we discuss how to go from ESD profiles to accelerations with such stacked data.For this, we would like to use the exact deprojection formula Eq. (2.3).The question is how to go from an individual lens to a large sample of stacked lenses, i.e. how to define a reasonable stacked radial acceleration g stacked obs , Note that the stacked ESD profile ∆Σ stacked is, at each radius R, a weighted average of the ESD profiles of the individual lenses.So it has a straightforward physical interpretation as an average ESD profile.We would like to have a similarly straightforward interpretation for the stacked acceleration g stacked obs .
Unfortunately, the simplest idea one might have for a stacked radial acceleration does not in general have such a straightforward interpretation.Namely, one might want to simply apply Eq. ( 2.3) to the stacked ESD profile ∆Σ stacked (i.e.first stack, then deproject), But this is not a weighted average of the radial accelerations of the individual galaxies.The reason is that the weights w l (R) in the stacked ESD profile depend on R. Indeed, if we write out this definition in full, we have In contrast, any weighted average of the radial accelerations of the individual lenses can be written in the form, for some weights wl (R) with normalization factor N = l wl (R).Eq. (B.4) is of this form only when the combination N −1 (R) w l (R) does not depend on R, or if the ∆Σ l satisfy special properties.But this is not generally the case. 15 In practice, for the KiDS data we use (see Appendix B.2), the weights w l (R) of the lenses l all have roughly the same scaling with R. Specifically, they all roughly scale as the projected area covered by the radial bin R. Thus, the normalized weights w l (R)/ l ′ w l ′ (R) are roughly independent of R (but they do depend on l) so that Eq. (B.3) works quite well.Indeed, in practice, Eq. (B.3) gives results that are very close to what we get using the more general method we now propose, see Fig. 8 in Appendix C.2.
Here, we want to use a stacked radial acceleration that is an average of the radial accelerations of the individual lenses.That is, we want our definition of g stacked obs to be of the form Eq. (B.5), i.e. g stacked obs (R) = g averaged obs (R).This can equivalently be written as where g obs,l (R) denotes the radial acceleration of each individual lens l.The question then becomes how to choose the weights wl (R).Our choice is 15 Another disadvantage of Eq. (B.3) is that it cannot easily be generalized to stacking in g bar space instead of position space.For the special case of baryonic point particles, this is simple enough; Eq. (2.3) becomes g obs (g bar ) = 4G π/2 0 dθ∆Σ(g bar sin 2 θ).For an individual galaxy, this can be generalized to other baryonic mass distributions as well.However, in general, the integral then depends on properties other than g bar of the galaxy.This is a problem when applied to a stacked ESD profile ∆Σ stacked (g bar ) because then one no longer has access to properties of the individual galaxies other than g bar .In the following, we mostly consider baryonic point masses so this is only a secondary concern for us.We consider non-point masses only in Appendix C.6.acceleration g bar .If we approximate the baryonic mass as a point mass, this is where M b,l is the baryonic mass of the lens l.Other mass distributions require a different form of the function R l (g bar ), see for example Appendix C.6.For the stacked radial acceleration, we similarly have,

B.2 Stacking from observational data
We now explain how we obtain stacked ESD profiles and stacked radial accelerations from actual weak-lensing observations.Essentially, we define the weights w l (R) and the ESD estimates ∆Σ l (R) introduced in the previous section in terms of observational data.We start with the stacked ESD profile written in the form given by Ref. [1], Here, the sum over the source galaxies s runs over all sources within the radial bin R of the lens l.We indicate this in our notation by "D ls = R" next to the summation sign.The W ls are weights for each lens-source pair given by where w s estimates the precision of the ellipticity measurement of the source s and Σ crit,ls is the critical surface density that we define below.The tangential ellipticity ϵ t,ls of the source s is an estimate of the tangential shear γ t at a projected radius R from the lens l based on the second brightness moments of the source [33].It is given by where ϵ 1,s and ϵ 2,s are the ellipticity components of the source s with respect to an equatorial coordinate system and ϕ ls is the angle between the x-axis and the lens-source separation vector [33,42].The tangential shear is related to the ESD profile of the lens [31,32].In particular, the quantity Σ crit,ls ϵ t,ls is an estimate for the ESD of the lens l at radius R based on the tangential ellipticity ϵ t,ls of the source s, By taking a weighted average over the sources s, we get an estimate for the ESD of the lens l in a radial bin R that includes all relevant sources, If we now define we can write the original formula Eq. (B.12) for the stacked ESD in the following form, with the normalization factor N (R) = l w l (R).This is the form we already used in the previous section in Eq. (B.1).Thus, Eq. (B.16) and Eq.(B.17) are the definitions of ∆Σ l and w l in terms of observational quantities.
It remains to give the definition of the critical surface density Σ crit,ls .If one knows the exact redshifts of both the source, z s , and the lens, z l , one has In practice, however, we use a more complicated definition of Σ crit,ls that takes observational uncertainties into account.We discuss this in Sec.3; see Eq. (3.1).
For the continuous integrals in the deprojection formula Eq. (2.3) (see also Eq. (B.5)), we need to know ∆Σ l (R) for all values of R, not just for discrete radial bins, which is what we have discussed so far.For simplicity, we linearly interpolate ∆Σ l between these discrete bins (see Appendix C.1 for how we estimate the uncertainty associated with this interpolation).When there are no sources in some radial bin, we do not have an estimate of ∆Σ l there.In such cases, we linearly interpolate between the bins that do have sources.

C.1 Statistical and systematic uncertainties
Here, we discuss the statistical and systematic uncertainties of our stacked radial accelerations g obs and stacked ESD profiles ∆Σ.In the following, we show explicit expressions only for quantities stacked in position space.We use the same procedure when stacking in g bar space.We first repeat the formula Eq. (B.1) for our stacked ESD profiles ∆Σ(R) to make it easier for the reader to follow the rest of this section.We have where N (R) = l w l (R) normalizes the weights w l (R).The weights w l (R) and the ESD profiles ∆Σ l (R) of the lens l are defined in terms of observational data in Eq. (B.16) and Eq.(B.17).We similarly repeat the formula Eq. (B.6) for our stacked radial accelerations g obs (R).We have where N (R) = l wl (R) normalizes the weights wl (R) that are defined as the inverse of the squared statistical uncertainty of g obs,l , i.e. wl (R) = 1/σ 2 g obs ,l (R).We will define this statistical uncertainty in Eq. (C.7) below.The radial acceleration g obs,l of the lens l is given by We now consider the statistical uncertainties of our stacked ESD profiles and stacked radial accelerations.We start with the ESD profile estimates ∆Σ l for individual lenses l from Eq. (B.16).For these, we use the analytical error estimate used by Refs.[1,66], but adjust it for the fact that we are dealing with only a single lens for now.We find (see Appendix B.2 for the definition of W ls and Σ crit,ls ), where we take the ellipticity dispersion σ ϵ,s from Table I in Ref. [28] which lists values for five tomographic redshift bins.We choose the one that corresponds to the redshift of the source s.This expression follows directly from Eq. (B.16) and is simpler than the one derived in Ref. [66], even when considering only the diagonal elements of their covariance matrix.This is because each source occurs at most once in Eq. (B.16) since we consider only a single lens for now.
To obtain statistical uncertainties for the stacked ESD profiles and the stacked radial accelerations, we linearly propagate these σ ∆Σ l uncertainties, assuming that the ∆Σ l (R) are independent for different lenses l and different radii R. In particular, this assumes that only very few sources contribute to multiple lenses.In contrast, Refs.[1,66] take such cases into account.However, since our analysis -like that of Ref. [1] -relies on lenses being isolated, such cases should not be important.Indeed, we have verified that our simplified procedure reproduces the statistical uncertainties of Ref. [1] almost perfectly.
Concretely, we adopt the following statistical uncertainties for our stacked ESD profiles from Eq. (B.1), And similarly for our stacked radial accelerations from Eq. (B.6), The θ integral involves σ ∆Σ l at arbitrarily large radii and, as discussed in Sec. 2, we set σ ∆Σ l to zero beyond the last ∆Σ l data point.
In practice, we know ∆Σ l (R) only in a finite number of discrete radial bins, see Appendix B.2.As discussed in Sec. 2, this means there are systematic uncertainties related to both extrapolating beyond the last bin and interpolating between the discrete bins.For our Figure 6.The tangential and cross components of the ESD profile (left) and the radial acceleration (right) for about 45 million random coordinates for R < 3 Mpc (see Fig. 7 for larger radii), assuming h 70 = 1.These should ideally be close to zero.As discussed in the text, the statistical uncertainties shown here are only a rough estimate and should not be taken too seriously.That the tangential radial acceleration starts to systematically deviate from zero already at relatively small radii is likely due to the effect discussed in Appendix B.1, namely that one cannot downweigh bad data as optimally as for the ESD profile when using our preferred method (i.e.deproject first, then stack).This is supported by the fact that using Eq. ( 2.3) to convert the tangential ESD profile from the left panel to tangential radial accelerations (i.e.stack first, then deproject) does not reproduce the systematic trend in the right panel.We nevertheless trust our results because, after subtraction, the radial accelerations calculated using our "deproject first, then stack" method agree very well with those obtained using the "stack first, then deproject" method which has very different systematics (see Fig. 8).Thus, we assume that ∆Σ l,× is zero beyond the last data point rather than continuing with an SIS profile.Second, the cross profiles are obtained from the cross components ϵ ×,ls of the ellipticities instead of the tangential components ϵ t,ls .Using the notation of Appendix B.2, we have ϵ ×,ls = sin(2ϕ ls )ϵ 1,s − cos(2ϕ ls )ϵ 2,s . (C.11) Both tangential and cross profiles should ideally be close to zero for random coordinates.From Fig. 6 and Fig. 7, we see that there are systematic non-zero residuals at large radii.Such systematic effects at large radii are not uncommon, see for example Appendix A of Ref. [42].This effect is strongest for the tangential radial acceleration, for which the systematic Figure 8.The weak-lensing RAR obtained using our method of choice as described in Appendix B (i.e.deproject first, then stack, yellow diamonds) and using Eq.(B.3) (i.e.stack first, then deproject, gray diamonds).Since Eq. (B.3) uses the stacked ESD profile as its input, it is not affected by the systematic trend for the tangential radial acceleration of random coordinates in Fig. 6, right, but comes with its own systematic uncertainties instead (see Appendix B.1).That these two procedures with different systematics agree so well indicates that our results can be trusted.Error bars and bands are as in Fig. 2, except we do not show the stellar mass systematic uncertainty for clarity.
deviation from zero starts at smaller radii than for the other profiles shown.As discussed in Appendix B.1, the radial acceleration is particularly prone to such effects because one cannot downweigh unreliable data as well as for the ESD profiles.We suggest that this effect explains why the tangential radial acceleration shows a systematic deviation from zero already at relatively small radii.
This interpretation is supported by the following observation.Our radial accelerations are based on the integral formula Eq. (2.3).We have verified that applying this formula to the stacked ESD profile obtained from random coordinates (Fig. 6, left) does not reproduce the systematic trend we see in the tangential radial acceleration (Fig. 6, right).Thus, neither the data nor the integral formula by themselves are responsible for this systematic trend.But it may well be that the restrictions on downweighing bad data discussed in Appendix B.1 are.
We trust our results despite this systematic trend because -after subtracting the profile obtained from the random coordinates -the g obs obtained using our method of choice (i.e.deproject first, then stack) agrees very well with that obtained using Eq.(B.3) 17 (i.e.stack first, then deproject), see Fig. 8. Importantly, as explained in Appendix B, the "stack first, then deproject" method from Eq. (B.3) is not prone to the systematic effect we consider here because it uses the stacked and subtracted ESD profile as input which does not show such a systematic trend at relatively small radii, see Fig. 6.Fig. 8 shows some difference at the very last data point.But this is where systematic uncertainties are anyway large, so in practice the difference in this last data point is not important.
The statistical uncertainties shown in Fig. 6 and Fig. 7 are calculated as described in Appendix C.1.In particular, the calculation assumes that sources do not contribute to multiple lenses.That is a good approximation for the isolated lenses we consider in the main text.But it is not justified for the large sample of random coordinates we consider here.Thus, the numerical values of the statistical uncertainties shown in Fig. 6 and Fig. 7 give only a rough indication of the order of magnitude of the uncertainties.Note that these uncertainties do not enter any of our results in the main text.They are used only to guide the eye of the reader in Fig. 6 and Fig. 7.

C.3 Cross component from actual lenses
As discussed in Appendix C.2, both the cross and tangential components of the ESD profiles and radial accelerations should be zero for random coordinates.For actual lenses, the tangential component carries the lensing signal and should not be zero, while the cross component should still be zero.This can be used to validate the lensing data and the method used to analyse this data.In the following, we subtract the cross components obtained from random coordinates (see Appendix C.2) from the cross components of the stacked ESD profiles and stacked radial accelerations, just as we do for the tangential components (see Appendix C.1).
We first consider the cross-components of the stacked ESD profiles.Since the version of these cross ESD profiles from Ref. [1] is publicly available, we can directly compare our results to theirs.For this, we adopt the stellar and gas masses, stellar mass bins, choice of h 70 , and cuts on the KiDS bright sample from Ref. [1].We show the result in Fig. 9.We see that our results match those of Ref. [1] very closely.There are some differences, but these are much smaller than the error bars and may be due to any number of minor numerical differences.We also see that these cross ESD profiles are consistent with zero, as they should be.
In Fig. 10 we show the cross components of the stacked radial accelerations obtained using our method for converting ESD profiles to accelerations described in Sec. 2 (i.e.deproject first, then stack).In contrast to Fig. 9, we here use our own masses, mass bins, choice of h 70 , and cuts on the KiDS bright sample (see Sec. 4, Sec. 3 and Sec.5.3).For comparison, we also show the cross components of the radial accelerations obtained using Eq.(B.3) (i.e.stack first, then deproject).We see that both methods agree well with each other (see also Fig. 8) and produce a cross component consistent with zero.
There is a slight tendency for the cross components in Fig. 10 to fall below zero more often than above zero.However, this may well just be random fluctuations.Indeed, the cross components do generally stay within the error bars.In addition, we have verified that if we derive the cross-component of random coordinates as in Appendix C.2 but with a much smaller number of random coordinates matching our sample of actual lenses, it is not hard to find similar behavior, even after subtracting the cross components inferred from the (much larger) full random sample from Appendix C.2.
Figure 9.The cross components of the stacked ESD profiles for the four stellar mass bins defined by Ref. [1] as derived in this work (yellow) and in Ref. [1] (red).To allow a direct comparison to Ref. [1], for this plot we adopt their h 70 = 1, their stellar and gas masses, and their cuts on the KiDS bright sample, namely R isol = 3 Mpc/h 70 and log 10 M * /M ⊙ < 11.

C.4 Influence of larger radii on smaller radii
As discussed in Sec.5.1, the deprojection formula Eq. (2.3) involves data at arbitrarily large radii.Thus, one may worry that, if the isolation criterion fails at large radii, radial accelerations at all radii will be affected, not just those at large radii.Indeed, our results from Sec. 5.2 suggest that our isolation criterion is probably not reliable below g bar ≈ 10 −14 m/s 2 .Fig. 11 shows that this is not a problem in practice.We find very similar results when we artificially cut off the ESD profiles ∆Σ l of the individual lenses at g bar = 10 −14 m/s 2 instead of using data down to g bar = 10 −15 m/s 2 .We have verified that the same is true if we artificially cut off the data at g bar = 10 −13 m/s 2 .Mathematically, this is because most of the integration volume of the integral in Eq. (2.3) is close to R, i.e. large radii are downweighted relative to radii close to R.

C.5 R isol dependence with artificially fixed error bars
In Sec.5.2, we saw that ETGs are quite sensitive to the details of the isolation criterion as quantified by R isol , while, at least down to g bar ≈ 10 −14 m/s 2 , LTGs are almost independent of R isol .In particular, we considered R isol up to 4 Mpc/h 70 , which gives a stricter isolation criterion than that used in Ref. [1], namely R isol = 3 Mpc/h 70 .For large values of R isol , we found no significant difference between ETGs and LTGs.
Figure 10.The cross components of the stacked radial accelerations for the four baryonic mass bins defined in Sec.5.3 derived using our method described in Sec. 2 (i.e.deproject first, then stack, yellow).For comparison, we also show the cross components of the radial accelerations derived from the stacked ESD profile and Eq.(B.3) (i.e.stack first, then deproject, blue, see also Fig. 8).
A stricter isolation criterion implies a smaller lens sample and thus larger statistical uncertainties.Thus, in principle, the reason we find no difference between ETGs and LTGs may simply be that it becomes statistically insignificant because of the larger uncertainties.To counter this, Fig. 12 shows what happens when we artificially keep the error bars fixed at the level they are at for R isol = 3 Mpc/h 70 .
More specifically, we have verified that the statistical uncertainties scale as 1/ √ N to an excellent approximation.Here, N is the number of lenses in the sample.Thus, to keep the error bars at the level they are at for R isol = 3 Mpc/h 70 , we simply adjust the statistical uncertainties in the following way before calculating χ 2 , where N is the number of lenses that satisfy the isolation criterion with the value of R isol that is currently under consideration and N 3.0 is the number of lenses that satisfy the isolation criterion with R isol = 3 Mpc/h 70 .Comparing Fig. 4 and Fig. 12 shows that our results remain unchanged even with these artificially small statistical uncertainties.

C.6 Extended hot gas distribution
As discussed in Sec. 4, we assume that ETGs are surrounded by hot gas and we model this hot gas as a point mass.In reality, however, this hot gas may be quite extended.To illustrate Figure 11.The weak-lensing RAR derived using the deprojection formula Eq. (2.3) when using the full ESD profiles of the individual lenses (as in Fig. 2, yellow diamonds) and when artificially cutting off the data at g bar = 10 −14 m/s 2 (gray diamonds).We find very similar results in both cases.
Thus, if the isolation criterion fails below g bar = 10 −14 m/s 2 , this has only very little effect on g obs at larger values of g bar .Error bars and bands are as in Fig. 2, except we do not show the stellar mass systematic uncertainty for clarity.
the effect of this, here we follow Ref. [1] and model the hot gas around ETGs as having an SIS profile cut off at 100 kpc.
To obtain radial accelerations stacked in g bar space using this extended mass profile, we must modify the functions R l (g bar ) that map between radii R l and baryonic Newtonian accelerations g bar for each lens l, see Eq. (B.9) and Eq.(B.10).In particular, the point mass relation Eq. remains valid for LTGs, but for ETGs we now have where M * ,l is the stellar mass of the lens l, M gas,l is the hot gas mass of the lens l, and R c = 100 kpc is where we cut off the SIS profile.Fig. 13 shows that modeling the hot gas of ETGs as an SIS cut off at 100 kpc instead of a point mass has a small effect at relatively large g bar and leaves the results unchanged at small g bar .These results are due the SIS profile being less concentrated towards the center than a point particle.There would be a larger effect at large g bar for gas profiles that are  4, but with artificially rescaled error bars.The idea is that, at larger R isol , the lens sample is smaller so the error bars are larger.Thus, the number of sigmas at larger R isol may be small simply because of the larger error bars.To counter this, we artificially rescale the error bars to the level they were at for R isol = 3.0 Mpc/h 70 , see Eq. (C.12).Even with these artificially small error bars, the qualitative trend with R isol remains the same as in Fig. 4.
Figure 13.Same as Fig. 3, but we additionally show the result for ETGs with their surrounding hot gas modeled as an SIS cut off at 100 kpc instead of a point mass (gray color).This makes a small difference at relatively large g bar and does not change the result at small g bar .For clarity, we do not show lensing data for LTGs and kinematic data.  2 but with g obs derived using the SIS approximation proposed in Ref. [1] instead of our method based on the exact deprojection formula Eq. (2.3).The SIS method gives a similar result as our method, but the resulting RAR is less smooth and generally has a larger systematic uncertainty of 0.05 dex from converting ESD profiles to radial accelerations (see Sec. 2).even less concentrated towards the center.For profiles that are more concentrated towards the center, such as an NFW profile, we would expect a smaller effect.At sufficiently small g bar , where all the gas mass is enclosed so that g bar can be approximated to fall off as 1/r 2 , the effect will always be negligible, irrespective of the gas profile at larger g bar .

D.1 SIS method
Here, we compare the weak-lensing RAR derived using the SIS method proposed by Ref. [1] to that derived using our method based on the exact deprojection formula Eq. (2.3) (see Sec. 2).Fig. 14 shows the RAR derived using our lens sample but with the SIS method (i.e. using Eq.(2.2)).Comparing to Fig. 2 we see that both methods produce similar results, but the RAR produced using our new method is smoother and generally comes with a smaller systematic uncertainty.As discussed in Sec. 2, the systematic uncertainty using our new method becomes significant only close to the last data point.
Our lens sample satisfies a stricter isolation criterion than that of Ref. [1], namely R isol = 4 Mpc/h 70 rather than R isol = 3 Mpc/h 70 .This means a cleaner lens sample, but also a smaller sample and thus larger statistical fluctuations.Due to the integral in Eq. (2.3), our improved method still produces quite smooth results even with such reduced statistics.3, but not using our preferred baryonic masses and isolation criterion.Instead, we here adopt R isol = 3 Mpc/h 70 , Q ETG = 1.0, and the gas mass estimate of Ref. [1] in order to show that we can reproduce the finding of Ref. [1].

D.2 ETGs vs LTGs
In Fig. 15, we demonstrate that we can reproduce the finding of Ref. [1] that the RAR of ETGs and LTGs differs significantly when we adopt the same baryonic masses and isolation criterion.Specifically, for a reasonably direct comparison with Ref. [1], Fig. 15 uses R isol = 3 Mpc/h 70 , Q ETG = 1.0, and the gas mass estimate from Ref. [1] rather than our default choices R isol = 4 Mpc/h 70 , Q ETG = 1.4,and our gas mass estimate from Sec. 4.
In Sec.5.2 we discuss in detail how this result changes when adopting our baryonic masses and our isolation criterion.In Fig. 16, we show that adopting our baryonic masses (Q ETG = 1.4 and the gas mass estimates Eq. (4.2) and Eq.(4.1)) while keeping the isolation criterion of Ref. [1] (R isol = 3 Mpc/h 70 ) reduces the difference between ETGs and LTGs but does not eliminate it.

Figure 1 .
Figure1.The stellar mass histogram assuming our stellar masses for the whole sample (gray), for early type galaxies (ETGs in red), for late type galaxies (LTGs in blue), and for the overlap of ETGs and LTGs (green).

Figure 4 .
Figure 4. Top: The difference between the radial accelerations inferred from weak lensing and the RAR fitting function Eq. (5.1), measured in sigmas, as a function of how isolated the lenses are, quantified by R isol .We separately show the result for ETGs (red) and LTGs (blue) as well as for small (triangles with dashed lines) and large accelerations (diamonds with solid lines), see Eq.(5.3).LTGs are mostly unaffected by making the isolation criterion stricter.In contrast, ETGs do depend on R isol , but, with increasing R isol , tend towards what the fitting function Eq. (5.1) predicts.Middle and bottom: The actual accelerations behind these sigma numbers for R isol = 3 Mpc/h 70 and R isol = 4 Mpc/h 70 .

Fig. 5
Fig.5shows the RAR separately for four baryonic mass bins, with log 10 M b /M ⊙ bin edges[9.0,10.5, 10.8, 11.1, 11.5].All mass bins generally show similar behavior, following the fit function Eq. (5.1) down to g bar ≈ 10 −14 m/s 2 .This is especially true for larger baryonic masses.For smaller masses, there are larger deviations.A part of this is likely due to statistical fluctuations.Indeed, the sample is smaller for smaller masses, as indicated by the larger error bars in Fig.5.Still, there may be a small systematic shift to larger accelerations g obs .If there is such a systematic shift, a possible explanation is a failure of the isolation criterion: In contrast to larger galaxies, smaller galaxies can be affected even by a relatively small neighbor.This is exacerbated by the m = 20 magnitude limit of the KiDS bright sample[1,29] which means that many of the neighbors that affect small lenses are not detected.The effect of the magnitude limit was discussed in Appendix A of Ref.[1] who estimated that the effect on their results is small.However, only the effect on the whole lens sample was considered in Ref.[1], not the effect on a small-mass subsample.Indeed, small-mass galaxies are only a small part of the whole lens sample so, overall, they have only a small effect.If one

RM 2 c
l (g bar )| ETG,SIS =    G(M * ,l +M gas,l ) g bar , for g bar < a c,l , gas,l , for g bar ≥ a c,l .(C.13) where a g,l ≡ GM gas,l R , a c,l ≡ G(M gas,l + M * ,l )

Figure 12 .
Figure 12.Same as the top panel of Fig.4, but with artificially rescaled error bars.The idea is that, at larger R isol , the lens sample is smaller so the error bars are larger.Thus, the number of sigmas at larger R isol may be small simply because of the larger error bars.To counter this, we artificially rescale the error bars to the level they were at for R isol = 3.0 Mpc/h 70 , see Eq. (C.12).Even with these artificially small error bars, the qualitative trend with R isol remains the same as in Fig.4.

Figure 14 .
Figure 14. as Fig.2but with g obs derived using the SIS approximation proposed in Ref.[1] instead of our method based on the exact deprojection formula Eq. (2.3).The SIS method gives a similar result as our method, but the resulting RAR is less smooth and generally has a larger systematic uncertainty of 0.05 dex from converting ESD profiles to radial accelerations (see Sec. 2).

Figure 15 .
Figure 15.Same as Fig.3, but not using our preferred baryonic masses and isolation criterion.Instead, we here adopt R isol = 3 Mpc/h 70 , Q ETG = 1.0, and the gas mass estimate of Ref.[1] in order to show that we can reproduce the finding of Ref.[1].

Figure 16 .
Figure 16.As Fig.15but using our baryonic mass estimates from Sec. 4 (i.e. using Q ETG = 1.4 and the gas mass estimates Eq. (4.2) and Eq.(4.1)).This brings the RAR of ETGs and LTGs closer together, but does not eliminate the difference.
log 10 g bar log 10 g obs σ statistical log 10 g obs