SEAGLE--II: Constraints on feedback models in galaxy formation from massive early type strong lens galaxies

TL;DR

This study compares various galaxy formation models with observed strong lens galaxies, revealing that models with weaker feedback mechanisms best match the observed mass density slopes.

Contribution

It provides constraints on feedback models in galaxy formation by analyzing their impact on strong lens galaxy properties in simulations versus observations.

Findings

Models with weak feedback match observed mass density slopes.

Strong feedback models fail to reproduce the isothermal slopes.

Weak feedback models align better with SLACS and BELLS data.

Abstract

We use nine different galaxy formation scenarios in ten cosmological simulation boxes from the EAGLE suite of {\Lambda}CDM hydrodynamical simulations to assess the impact of feedback mechanisms in galaxy formation and compare these to observed strong gravitational lenses. To compare observations with simulations, we create strong lenses with $M_\star$ > $10^{11}$ $M_\odot$ with the appropriate resolution and noise level, and model them with an elliptical power-law mass model to constrain their total mass density slope. We also obtain the mass-size relation of the simulated lens-galaxy sample. We find significant variation in the total mass density slope at the Einstein radius and in the projected stellar mass-size relation, mainly due to different implementations of stellar and AGN feedback. We find that for lens selected galaxies, models with either too weak or too strong stellar and/or AGN feedback fail to explain the distribution of observed mass-density slopes, with the counter-intuitive trend that increasing the feedback steepens the mass density slope around the Einstein radius ($\approx$ 3-10 kpc). Models in which stellar feedback becomes inefficient at high gas densities, or weaker AGN feedback with a higher duty cycle, produce strong lenses with total mass density slopes close to isothermal (i.e. -d log({\rho})/d log(r) $\approx$ 2.0) and slope distributions statistically agreeing with observed strong lens galaxies in SLACS and BELLS. Agreement is only slightly worse with the more heterogeneous SL2S lens galaxy sample. Observations of strong-lens selected galaxies thus appear to favor models with relatively weak feedback in massive galaxies.

Figures

Tables

Table 1. Table 1: Main sub-grid parameters of the EAGLE simulations used in this work. Columns (left to right) are the name of the simulation, L 𝐿 L : the comoving side length of the volume, N 𝑁 N : the number of particles for individual type i.e. gas and DM, γ eos subscript 𝛾 eos \gamma_{\rm eos} : the power-law slope of the polytropic equation of state, n H ⋆ superscript subscript 𝑛 H ⋆ n_{\rm H}^{\star} : the star formation density threshold, f th subscript 𝑓 th f_{\rm th} : the star formation feedback efficiency, f th , max subscript 𝑓 th max f_{\rm th,max} : the asymptotic maximum and f th , min subscript 𝑓 th min f_{\rm th,min} : minimum values of f th subscript 𝑓 th f_{\rm th} , n H , 0 subscript 𝑛 H 0 n_{\rm H,0} : the density-term denominator for the Reference model , n n subscript 𝑛 n n_{\rm n} : the Reference model density-term exponent (from equation 4 ), C visc subscript 𝐶 visc C_{\rm visc} : the sub-grid accretion disc viscosity parameter (from equation 7 in Crain et al. 2015 ), and Δ ​ T AGN Δ subscript 𝑇 AGN \Delta T_{\rm AGN} : the temperature increment of stochastic AGN heating. The calibrated models reproduce the GSMF at z = 0.1 𝑧 0.1 z=0.1 . The reference variation models adopt a single-parameter variation of the Reference simulation (varied parameters are highlighted in bold). Except for FB σ 𝜎 \sigma (which uses the parameter n T subscript 𝑛 T n_{\rm T} ), all other models have n Z = 2 / ln ⁡ 10 subscript 𝑛 Z 2 10 n_{\rm Z}=2/\ln{10} with the same numerical value (see equation 2 ). For FBconst, this parameter is not applicable. This Table is partially reproduced from Crain et al. ( 2015 ) .

Identifier	Side length	$N$	${\gamma }_{\mathrm{eos}}$	$n_H^{\star }$	$f_{\mathrm{th}}$-scaling	$f_{\mathrm{th},\max }$	$f_{\mathrm{th},\min }$	$n_{H,0}$	$n_n$	$C_{\mathrm{visc}}/2\pi$	$\Delta T_{\mathrm{AGN}}$
	$L$ [cMpc]			[${\mathrm{cm}}^{-3}$]				[${\mathrm{cm}}^{-3}$]			${\log }_{10}$ [K]
Calibrated models
FBconst	50	752	$4/3$	Eq. 1	$-$	$1.0$	$1.0$	$-$	$-$	${10}^3$	$8.5$
FB$\sigma$	50	752	$4/3$	Eq. 1	${\sigma }_{\mathrm{DM}}^2$	$3.0$	$0.3$	$-$	$-$	${10}^2$	$8.5$
FBZ	50	752	$4/3$	Eq. 1	$Z$	$3.0$	$0.3$	$-$	$-$	${10}^2$	$8.5$
Ref (FBZ$\rho$)	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
Ref-100 (FBZ$\rho$)	100	1504	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
Reference-variations
ViscLo	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟐}}$	$8.5$
ViscHi	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$8.5$
AGNdT8	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.0$
AGNdT9	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$9.0$
NOAGN	50	752	$4/3$	Eq. 1	$Z,\rho$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$-$

Table 2. Table 2: Summary of the simulation settings and output products.

Galaxy Selection
Observable	Value	Name		Comments
$M_{\star }$	$\ge 1.76\times {10}^{10}{M}_{\odot }$	Stellar mass threshold		Taken from Auger et al. (2010a)
$\sigma$	> 120 km/sec	Stellar velocity dispersion		Kept lower than SLACS
$R_{50}$	> 1 kpc	Projected half-mass radius
Lens Candidates
	$M_{\star }$ threshold	$M_{\star }$ threshold
	for follow-up work	for this work
	- - - - - - - - - - - -	- - - - - - - - - - -	After 3
Simulation	$\ge 1.76\times {10}^{10}{M}_{\odot }$	>${10}^{11}{M}_{\odot }$	projections	Comments
Reference-100cMpc	-	67	201	100 cMpc box.
Reference-50 (FBZ$\rho$)	252	25	75	50 cMpcbox
FBconst	279	22	66	′′
FB$\sigma$	259	22	66	′′
FBZ	312	19	57	′′
ViscLo	289	29	87	′′
ViscHi	188	14	42	′′
AGNdT8	276	27	81	′′
AGNdT9	194	8	24	′′
NOAGN	312	37	111	′′
Object-properties	Value	Type		Comments
Orientation	3	x, y, z		Projected surface density maps
Redshift	$z_l=0.271$	-		Consistent with SLACS’
				mean lens-redshift of 0.3
Source Properties
Parameters	Value	Unit		Comments
Source Type	Sérsic	-		Consistent with SLACS lenses
				(Newton et al. 2011)
Brightness	23	apparent mag.		′′
Size ($R_{\mathrm{eff}}$)	0.2	arcsec		′′
Axis ratio ($q_s$)	0.6	-		′′
Sérsic Index	1	-		′′
Redshift	$z_s$=0.6	-		′′
Position	Random	Within caustics		Producing rings and arcs lens systems
				consistent with SLACS
Instrumental Settings
Parameters	Type	Value	Comments
PSF	Gaussian	FWHM=0.1 arcsec	-
Noise	HST ACS-F814W	2400 sec	-
Image Properties
Map used	Properties	Value
Surface density	(a) Size	512$\times$512 pixels
	(b) Units	pkpc
$\kappa$, Inv. mag. map and Lens	(a) Size	161$\times$161 pixels
	(b) Units	degrees

Table 3. Table 3: The median values of mass density slopes, t 𝑡 t , of the simulated lenses in different galaxy formation models subjected to weighting scheme with α 𝛼 \alpha =0.5 and 1.5 and their respective fractional change. Table 5 have the value for α 𝛼 \alpha =1.0.

Simulation	$\alpha$=0.5	$\alpha$=1.5	|$\Delta t$|/$t$
Ref-50	2.16	2.20	0.02
FBconst	1.98	2.08	0.05
FB$\sigma$	1.68	1.75	0.04
FBZ	1.61	1.81	0.12
ViscLo	1.68	1.88	0.12
ViscHi	2.10	2.18	0.04
AGNdT8	1.96	2.08	0.06
AGNdT9	2.11	2.17	0.03
NOAGN	1.61	1.79	0.11

Table 4. Table 4: The mean values of effective radius, R eff subscript 𝑅 eff R_{\rm eff} , of the lensing galaxies in different galaxy formation models and their respective mean Einstein radius, R Ein subscript 𝑅 Ein R_{\rm Ein} . The ratio R Ein / R eff subscript 𝑅 Ein subscript 𝑅 eff R_{\rm Ein}/R_{\rm eff} gives a good estimate of the type of strong lenses simulated from EAGLE and observations.

Simulation	<$\log \left(R_{\mathrm{eff}}\right)$>	rms	<$\log \left(R_{\mathrm{Ein}}\right)$>	rms	$R_{\mathrm{Ein}}/R_{\mathrm{eff}}$	rms
Ref-50	0.91	0.21	0.65	0.34	0.71	0.23
FBconst	0.84	0.26	0.68	0.35	0.83	0.22
FB$\sigma$	0.82	0.23	0.77	0.36	0.94	0.20
FBZ	0.72	0.28	0.81	0.33	1.13	0.21
ViscLo	0.83	0.20	0.77	0.30	0.93	0.25
ViscHi	1.08	0.13	0.52	0.27	0.46	0.27
AGNdT8	0.84	0.19	0.64	0.28	0.76	0.22
AGNdT9	1.13	0.16	0.71	0.43	0.63	0.21
NOAGN	0.56	0.23	0.75	0.35	1.33	0.24
\hdashlineSLACS	0.86	0.51	0.59	0.11	0.69	0.18
SL2S	0.83	0.49	0.95	0.60	1.52	0.70
BELLS	1.03	0.76	1.05	0.62	1.03	0.88

Table 5. Table 5: Mean, standard deviation and median values of mass density slopes inferred from lens modelling, t 𝑡 t , of the simulated lenses in different galaxy formation models. The KS test results for the mass density slopes (1D) and mass-size relation (2D) compared to SLACS, BELLS and SL2S, are also listed. The p-values that exceed 0.05, and hence indicate an acceptable agreement between the simulations and observations, are shown in bold.

$\log M_{\star }/M_{\odot }=11.0-12.0$
	Mass density slope ($t$)			Mass density slope KS test						Mass-size KS test
	. . . . . . . . . . . . . . . . .			. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .						. . . . . . . . . . . . . . . . .
Simulation	Mean	Std.	Median	SLACS		SL2S		BELLS		SLACS
				D-value	p-value	D-value	p-value	D-value	p-value	D-value	p-value
Ref-100	2.09	0.26	2.24	0.26	0.53e-2	0.43	0.46e-3	0.42	0.17e-2	0.44	0.57e-2
Ref-50	2.19	0.25	2.20	0.35	0.15e-5	0.51	0.27e-5	0.48	0.59e-5	0.41	0.29
FBconst	2.00	0.22	2.06	0.15	0.39	0.36	0.005	0.17	0.63	0.47	0.15
FB$\sigma$	1.62	0.22	1.60	0.76	1.25e-26	0.77	4.44e-13	0.99	2.52e-19	0.48	0.11
FBZ	1.60	0.21	1.65	0.82	5.08e-27	0.84	2.23e-14	0.63	1.24e-7	0.53	0.02
ViscLo	1.64	0.25	1.61	0.68	1.2e-22	0.65	0.9e-10	0.46	0.001	0.52	0.002
ViscHi	2.09	0.23	2.24	0.17	0.09	0.22	0.15	0.21	0.26	0.77	1.95e-7
AGNdT8	1.95	0.22	2.00	0.38	0.12	0.36	0.003	0.21	0.26	0.44	0.24
AGNdT9	2.18	0.24	2.25	0.23	0.01	0.24	0.10	0.22	0.23	0.82	1.17e-5
NOAGN	1.67	0.20	1.47	0.78	5.06e-20	0.78	1.38e-11	0.51	0.11e-3	0.58	5.12e-6

Table 6. Table 6: The priors used in the modeling with an EPL plus shear mass model, using LENSED .

Parameter	Prior type⋆⋆	Prior range				Description
		$\mu$	$\sigma$	min	max
$x_L$	norm	80.0	5.0	-	-	Lens position: x coordinate
$y_L$	norm	80.0	5.0	-	-	Lens position: y coordinate
$r_L$	unif	-	-	5.0	70.0	Einstein radius in pixel units
$t_L$	norm	1.1	0.1	-	-	Surface mass density slope
$q_L$	unif	-	-	0.2	0.99	Lens axis ratio
${\varphi }_L$	unif	-	-	0.0	180.0	Lens position angle in degrees, wrapped around
$\gamma _1^{}{}_L{}^{}$	norm	0.0	0.01	-	-	Shear vector
$\gamma _2^{}{}_L{}^{}$	norm	0.0	0.01	-	-	Shear vector
$x_S$	norm	80.0	30.0	-	-	Source position: x coordinate
$y_S$	norm	80.0	30.0	-	-	Source position: y coordinate
$r_S$	unif	-	-	0.1	10.0	Source size in pixel units
$ma{g}_S$	unif	-	-	-5.0	0.0	Source magnitude, adjusted with the background magnitude#
$n_S$	norm	1.0	0.1	-	-	Sérsic index
$q_S$	norm	0.5	0.1	-	-	Source axis ratio
${\varphi }_S$	unif	-	-	0.0	180.0	Source position angle in degrees, wrapped around
$\star$ All values are in pixels except $q$, $\gamma$, $t_L$, $ma{g}_S$, $n_S$, and $\varphi$. $\star \star$ norm = Gaussian (with mean $\mu$ and standard dev. $\sigma$), unif = Uniform
$\#$ Source’s real magnitude = Background magnitude - $ma{g}_s$, where background magnitude is flux due to background in mag/arcsec2

Table 7. Table 7: Comparison of the modeled density slopes and key parameters using different source structures for a typical EAGLE lens. The source redshift is 0.6 and lens redshift is 0.271. The remaining settings are the same as mentioned in the modeling section of the paper.

Input Sérsic source parameters
$R_{\mathrm{Sersic}}$=0.1 arcsec, $n_{\mathrm{Sersic}}$= 1 and 3
Modelled output
Parameters	SIE	EPL
$R_{\mathrm{Ein}}\left(\mathrm{arcsec}\right)$	2.62	2.62
$\mathrm{Density}\mathrm{Slope}\left(t\right)$	$\equiv$2.00	2.049
${\mathrm{Ellipticity}}_{\mathrm{Lens}}$	0.158	0.158
${\mathrm{PA}}_{\mathrm{Lens}}$	78.17	77.98
$\mathrm{Source}\left(n_{\mathrm{Sersic}}\right)$	1.637	1.617
$\mathrm{Source}{R}_{\mathrm{Sersic}}\left(\mathrm{arcsec}\right)$	0.097	0.118
${\chi }^2$	0.961	0.959
Input Sérsic source parameters
$R_{\mathrm{Sersic}}$=0.1 arcsec, $n_{\mathrm{Sersic}}$= 3
Modelled output
Parameters	SIE	EPL
$R_{\mathrm{Ein}}\left(\mathrm{arcsec}\right)$	2.62	2.63
$\mathrm{Density}\mathrm{Slope}\left(t\right)$	$\equiv$2.00	2.045
${\mathrm{Ellipticity}}_{\mathrm{Lens}}$	0.157	0.158
${\mathrm{PA}}_{\mathrm{Lens}}$	78.04	78.63
$\mathrm{Source}\left(n_{\mathrm{Sersic}}\right)$	2.754	2.771
$\mathrm{Source}{R}_{\mathrm{Sersic}}\left(\mathrm{arcsec}\right)$	0.107	0.095
${\chi }^2$	0.925	0.916
Input Sérsic source parameters
$R_{\mathrm{Sersic}}$=0.1 arcsec, $n_{\mathrm{Sersic}}$= 1 and 3, random noise
Modelled output
Parameters	SIE	EPL
$R_{\mathrm{Ein}}\left(\mathrm{arcsec}\right)$	2.61	2.61
$\mathrm{Density}\mathrm{Slope}\left(t\right)$	$\equiv$2.00	1.988
${\mathrm{Ellipticity}}_{\mathrm{Lens}}$	0.158	0.156
${\mathrm{PA}}_{\mathrm{Lens}}$	78.95	78.63
$\mathrm{Source}\left(n_{\mathrm{Sersic}}\right)$	1.672	1.634
$\mathrm{Source}{R}_{\mathrm{Sersic}}\left(\mathrm{arcsec}\right)$	0.105	0.102
${\chi }^2$	0.922	0.914

Table 8. Table 8: Comparison of the modeled density slopes from different model variations of EAGLE than Reference with different source sizes from that used in the main analysis.

Source sizes†	FBconst	FBZ	FB$\sigma$	AGNdT8	AGNdT9	ViscLo	ViscHi	NOAGN
- - - - - - - -	- - - - - -	- - - -	- - - -	- - - - - -	- - - - - - -	- - - - -	- - - - -	- - - - -
(arcsec)	Mass density slope ($t$)
0.05	2.00	1.66	1.62	1.99	2.21	1.62	2.18	1.49
0.1	2.01	1.66	1.62	1.97	2.19	1.60	2.16	1.55
0.3	2.01	1.65	1.60	1.98	2.19	1.61	2.16	1.57
0.4	2.01	1.65	1.63	1.98	2.18	1.61	2.16	1.60
$†$ All other source parameters have been kept unchanged.

Equations

$$n_{\mathrm{H}}^{\star }(Z)=\min \left[0.1\mathrm{c}{\mathrm{m}}^3{\left(\frac{\mathrm{Z}}{0.002}\right)}^{-0.64},10\mathrm{c}{\mathrm{m}}^3\right],$$

$$f_{\mathrm{th}}=f_{\mathrm{th},\min }+\frac{f_{\mathrm{th},\max }-f_{\mathrm{th},\min }}{1+{\left(\frac{T_{\mathrm{DM}}}{10^5\mathrm{K}}\right)}^{n_{\mathrm{T}}}}.$$

$$f_{\mathrm{th}}=f_{\mathrm{th},\min }+\frac{f_{\mathrm{th},\max }-f_{\mathrm{th},\min }}{1+{\left(\frac{Z}{0.1Z_{\odot }}\right)}^{n_Z}},$$

$$f_{\mathrm{th}}=f_{\mathrm{th},\min }+\frac{f_{\mathrm{th},\max }-f_{\mathrm{th},\min }}{1+{\left(\frac{Z}{0.1Z_{\odot }}\right)}^{n_Z}{\left(\frac{n_{\mathrm{H},\mathrm{birth}}}{n_{\mathrm{H},0}}\right)}^{-n_n}},$$

$$\kappa (R)=\frac{b}{2R},$$

$$\kappa (R)=\frac{(2-t_L)}{2}{\left(\frac{b}{R}\right)}^{t_L},$$

$$\text{centering}W(M_{\star })\equiv {\left(\frac{M_{\star }}{⟨M_{\star }⟩}\right)}^{\alpha },\text{@add@centering}$$

Full text

††pubyear: 2018††pagerange: SEAGLE–II: Constraints on feedback models in galaxy formation from massive early-type strong lens galaxies–D

1 Introduction

Large-scale numerical simulations have established the Cold Dark Matter (CDM) paradigm as a viable framework for galaxy formation (e.g. Davis et al. 1985; Frenk et al. 1988). The CDM model predicts that galaxies form in dark matter halos having a Navarro-Frenk-White (NFW) density profile (Navarro et al. 1996, 1997) and predict the abundance and distribution of substructures within these halos (e.g. Gao et al. 2004; Springel 2010). The physics of galaxy formation, however, complicates the description of the matter distribution on small (several kpc) scales. Moreover, the central regions of CDM halos can also be strongly modified by baryonic matter and their associated physical processes. Baryons settle into the centers of density concentrations due to dissipation, thereby modifying the inner DM slopes (e.g. Duffy et al. 2010; Sonnenfeld et al. 2012; Grillo 2012; Remus et al. 2013; Cappellari et al. 2013; Tortora et al. 2009, 2010; Tortora et al. 2014a; Pontzen & Governato 2014). Because a complete analytic theory of baryonic physics is lacking, hydrodynamic simulations that include many physical processes have emerged as the dominant tool to study the complex non-linear interactions taking place during galaxy formation (e.g. Schaye et al. 2010; Vogelsberger et al. 2014; Schaye et al. 2015; Dubois et al. 2016; Hopkins et al. 2016). State-of-the art hydrodynamical simulations with improved stellar and AGN feedback, for example, can reproduce the cosmic star formation history of the Universe and the galaxy stellar mass function. Hydrodynamic simulations are currently working only above certain mass and spatial resolutions, however, and physical processes on smaller scales are implemented via analytic prescriptions known as ‘sub-grid physics’. The impact of varying sub-grid physics prescriptions on large representative populations of stellar systems was first systematically explored in the ‘OverWhelmingly Large Simulations’ project (OWLS; Schaye et al. 2010), a suite of over fifty large cosmological hydrodynamical simulations with varying sub-grid physics. Calibration of sub-grid prescriptions to reproduce a limited number of observables has been explored extensively (Vogelsberger et al. 2014; Schaye et al. 2015; Crain et al. 2015; McCarthy et al. 2017), showing that their exact parameterizations are very important.

Strong gravitational lensing is one of the most robust and powerful techniques to measure the total mass and its distribution in galaxies on kpc scales (Kochanek 1991; Koopmans et al. 2006), allowing their inner structure and evolution over cosmic time to be studied in detail (Treu et al. 2006; Treu et al. 2009; Koopmans et al. 2006; Koopmans et al. 2009; Dutton & Treu 2014), independently of the nature of the matter or its dynamical state. In particular, the mass density profile of massive lensing galaxies at $z>0.1$ can trace their formation and evolution mechanisms (e.g. Barnabè et al. 2009, 2011). The last two decades have seen major progress in observational studies of strong lensing thanks to surveys such as the Lenses Structure and Dynamics survey (LSD; Treu & Koopmans 2004), the Sloan Lens ACS Survey (SLACS; Bolton et al. 2006; Koopmans et al. 2006; Bolton et al. 2008a, c; Koopmans et al. 2009; Auger et al. 2010a; Auger et al. 2010b; Shu et al. 2015, 2017), the Strong Lensing Legacy Survey (SL2S; Cabanac et al. 2007; Ruff et al. 2011; Gavazzi et al. 2012; Sonnenfeld et al. 2013a; Sonnenfeld et al. 2013b; Sonnenfeld et al. 2015) and the BOSS Emission-Line Lens Survey (BELLS; Brownstein et al. 2012). Future surveys such as the Euclid (Laureijs et al. 2011) and the Large Synoptic Survey Telescope (LSST; Ivezić et al. 2008), as well as the ongoing Kilo Degree Survey (KiDS; de Jong et al. 2015) and the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005), are expected to increase the number of known strong lenses by several orders of magnitude (Petrillo et al. 2017; Metcalf et al. 2018; Treu et al. 2018) and revolutionize strong lensing studies.

The paper is structured as follows. In Section 2, we summarize the EAGLE galaxy formation simulations and the relevant codes that are used in this paper. Section 3 describes the simulation and analysis pipeline. The mass models used are described in Section 4. We give a brief description of the strong lensing observations in Section 5. In Section 6, we compare mock lens samples with observations, in terms of their mass-size relations and the total matter density slopes. The implications of our results are discussed and summarized in Section 7. Throughout the paper, we use EAGLE simulations that assume a Chabrier stellar Initial Mass Function (IMF, Chabrier 2003) and compare these to observables derived under the same IMF assumption. The values of the cosmological parameters are $\mathrm{\Omega_{\Lambda}}$ = 0.693, $\mathrm{\Omega_{b}}$ = 0.0482519, $\mathrm{\Omega_{m}}$ = 0.307, ${h=H_{0}/(100\;{\rm km\;s^{-1}\;Mpc^{-1}})}$ = 0.6777 and ${\sigma_{8}}$ = 0.8288. These are taken from the Planck satellite data release (Planck Collaboration et al. 2014).

2 EAGLE Simulations

Although there have been simulation studies of strong lensing focusing on the mass-size relations, the total density slope and other observables (e.g. Remus et al. 2017; Peirani et al. 2017; Xu et al. 2017), the impacts of varying sub-grid physics (in particular baryonic feedback) on lensing statistics, their mass density slopes and stellar masses and sizes have not been studied comprehensively yet (Peirani et al. 2018). Duffy et al. (2010) analyzed the impact of baryon physics on dark matter structure but only had low-resolution models at low redshift.

Mukherjee et al. (2018) (hereafter M18), introduced the SEAGLE pipeline to systematically study galaxy formation via simulated strong lenses. SEAGLE aims to investigate and possibly disentangle galaxy formation and evolution mechanisms by comparing strong lens early-type galaxies (ETGs) from hydrodynamic simulations with those observed, analyzing them in a similar manner (although this is not always exactly possible).

As in M18, we make use of the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations (Schaye et al. 2015; Crain et al. 2015; McAlpine et al. 2016) – a suite of state-of-the-art hydrodynamical simulations – to create, model and analyze simulated strong lens-galaxies and compare them with observations. Throughout this study, we use ten selected galaxy formation scenarios (i.e. having different sub-grid physics prescriptions; Schaye et al. 2015; Crain et al. 2015), the GLAMER ray-tracing package (Metcalf & Petkova 2014; Petkova et al. 2014), and the LENSED lens-modeling code (Tessore et al. 2016). We preselect potential strong lenses based on their stellar masses and create projected mass maps for three different orientations. We calculate the half-mass radius from the simulated mass maps. We create mock lenses by ray tracing through the mass maps, placing an analytic Sersic (1968) source, at a higher redshift, having observationally motivated parameters. We ignore line-of-sight effects, which for massive ETGs is expected to be a good approximation (see e.g., Koopmans et al., 2006). We use a single-orbit HST-ACS F814W noise level and PSF to mimic strong lenses found in SLACS and BELLS observations (Auger et al., 2010a; Bolton et al., 2012a).

Throughout this work, we also discuss possible observational systematics (e.g. differences in model-fitting methodologies, differences in filters/bands of the observational surveys, possible lens selection biases, etc.) as well as resolution effects in the simulations, that might affect their comparison. The main aim of this study, however, is to illustrate the effects of the sub-grid physics parametrization adopted by the EAGLE models, and the strong sensitivity of a number of strong lens observables e.g. total mass density slope, mass-size relation, and Einstein radius to the variation of the key sub-grid physics. In future work, we will analyze other properties such as the dark matter fractions and the stellar Initial Mass Function (IMF). Although we assume a Chabrier IMF in this work, the impact of assuming a different IMF (e.g. stellar mass and feedback) is partially removed during the process of calibration (see Section 2.3). The impact of a changing IMF should therefore be very carefully examined and will be done in a future publication for the Reference model for which these models are available (see e.g. Barber et al., 2018). A full analysis is currently not possible for the other models and well beyond the scope of this work, where we focus on the impact of galaxy-formation models.

In this section we describe the EAGLE simulations used in this study. In Section 2.1, we broadly describe the types of model-variations that have been chosen and in Section 2.2, we describe the simulation setup and the sub-grid physics recipes that are used in those model variations. Section 2.3 describes the calibrated simulations and reference models variation are summarized in Section 2.4. The details presented here are kept concise, yet informative, to make this paper self-contained.

2.1 EAGLE model variations

The simulations explored in this paper are taken from Crain et al. (2015) plus the 100cMpc-Reference run from Schaye et al. (2015). Crain et al. (2015) divided the simulations into two categories. The first comprises four simulations calibrated to yield the $z=0.1$ galaxy stellar mass function (GSMF) and central black hole (BH) masses as a function of galaxy stellar mass. The second category comprises simulations that each vary a single sub-grid physics parameter with respect to the Reference model but without considering whether they match the GSMF (i.e. they are not calibrated). In the calibrated simulations, the models differ in terms of their adopted efficiency of feedback associated with star formation, and how this efficiency depends upon the local environment. In the Reference variation simulations, the sensitivity of the resulting galaxies to these variations are assessed. We note that similar variations have previously been done in the OWLS project (Schaye et al., 2010). The general conclusion from previous work has been that the properties of simulated galaxies are most sensitive to the efficiency of baryonic feedback (see e.g., Schaye et al., 2010; Scannapieco et al., 2012; Haas et al., 2013a, b; Vogelsberger et al., 2013). This has motivated us to largely focus in this study on the effect of baryonic feedback on lensing observables, in particular on the total mass density profile in the inner regions of massive ETGs ($\sim$5 kpc), which was not considered during the calibration process and thus is a more reliable tracer of various formation processes.

2.2 Subgrid physics

Any simulation has a certain resolution limit below which the physical processes cannot be simulated via the dynamics of the particles. Similarly, the physical processes on scales smaller than the resolution of the EAGLE simulations are incorporated via analytic prescriptions. In EAGLE eleven chemical elements have been considered in the simulations. The calculations of radiative cooling and heating rates using the CLOUDY (version 07.02) code of Ferland et al. (1998), account for variations in metallicity and for variations in the relative abundances of individual elements. The cooling rates are specified as a function of density, temperature and redshift. While implementing the cooling in EAGLE simulations, it is assumed that the optically thin gas is in a state of ionization equilibrium and is exposed to the CMB and an instantaneous, spatially uniform, temporally-evolving (Haardt & Madau, 2001) UV/X-ray background (Wiersma et al. 2009a). Stochastic star formation, as formulated by Schaye & Dalla Vecchia (2008), has been implemented, but with the metallicity-dependent density threshold of Schaye (2004). A density threshold for star formation, $n_{\rm H}^{\star}$, was imposed because star formation occurs only in cold ($T\ll 10^{4}\rm K$), dense gas. Because the transition from a warm, neutral phase to a cold, molecular phase only occurs at lower densities and pressures in more metal-rich (and hence dust-rich) gas, the metallicity-dependent star formation threshold put forward by Schaye 2004 (see his equations 19 and 24) was adopted:

[TABLE]

where Z is the gas metallicity. Every star particle constitutes a stellar population with a fixed Chabrier (2003) IMF. The mass-to-light (M/L) ratio includes all the stellar remnants. The stellar evolution and mass loss implemented in EAGLE, is based on the prescription proposed in Wiersma et al. (2009b). The simulations adopt the stochastic thermal stellar feedback scheme of Dalla Vecchia & Schaye (2012), in which the temperature increment, $\bigtriangleup T_{\rm SF}$, of heated resolution elements is specified. The fraction of the supernova energy budget that is available for feedback determines the probability that a resolution element neighboring a young star particle is heated. This fraction is referred to as $f_{\rm th}$ (Dalla Vecchia & Schaye, 2012). According to the convention, $f_{\rm th}=1$ equates to $1.736\times 10^{49}\rm erg\ M^{-1}_{\odot}$, being the level of injected energy per stellar mass formed. Lastly, AGN feedback has been implemented via a single mode, where energy is injected thermally and stochastically, analogous to energy feedback from star formation.

2.3 Calibrated simulations

In EAGLE model variations, the efficiency of the stellar feedback and the BH accretion were calibrated to broadly match the observed local ($z\approx$ 0) GSMF, subject to the constraint that galaxy sizes must be in agreement with observations. We explain why calibration was needed and then we briefly describe the calibrated simulations of Crain et al. (2015), that are also used in this paper. Table 1 provides a concise overview of all the important parameters and a brief description of the four calibrated EAGLE simulations, adapted from the above-mentioned work.

2.3.1 The necessity of calibration

The choice of sub-grid routines and the adjustment of their parameters can result in substantial alterations of the simulation outcomes. Schaye et al. (2015) argued that the appropriate methodology for cosmological simulations is to calibrate the parameters of the uncertain sub-grid routines for feedback with a small number of key observations, in order that simulations reproduce those representative observables, and then compare properties (between simulations and observations) whose quantities that are not considered during the calibration. The total mass density slope, examined in this paper, is one of those which was not used in calibration. The results thus obtained can reasonably be considered being a consequence of the implemented astrophysics. On the other hand, the impact of changing the IMF (e.g. Barber et al., 2018) is partly calibrated out, and will be more carefully considered in a separate paper for the Reference model.

2.3.2 A Constant Feedback (FBconst)

This is the simplest feedback model where, independently from the local conditions, a fixed amount of energy per unit stellar mass is injected into the ISM. This fixed energy corresponds to the total energy discharged by type-II SNe ($f_{\rm th}=1$). The thermal stellar feedback prescription employed in EAGLE becomes inefficient at high gas densities due to resolution effects (Dalla Vecchia & Schaye 2012). Model Reference (see Section 2.3.5) compensates for this known artifact by injecting more energy at higher gas density. Because this is not done in FBconst, the stellar feedback will be less effective in high-mass galaxies (where the gas tends to have higher densities) (Crain et al. 2015).

2.3.3 Velocity dispersion dependent feedback (FB$\sigma$)

This model prescribes stellar feedback based on the local conditions, inferred from neighboring DM particles. The efficiency, $f_{\rm th}$, is calibrated as a function of the square of the 3-dimensional velocity dispersion of the DM particles within a stellar particle’s smoothing kernel at the time of its birth ($\sigma_{\rm DM}^{2}$).

The prescription of $f_{\rm th}$ in its functional form, is a logistic (sigmoid) given by,

[TABLE]

$T_{\rm DM}$ is the temparature of the characteristic virial scale of environment of the star particle. The parameter $n_{\rm T}>0$ controls how rapidly $f_{\rm th}$ transitions as the dark matter ‘temperature’ scale deviates from $10^{5}\;\rm{K}$.

2.3.4 Metallicity dependent feedback (FBZ)

This model makes the radiative losses, $f_{\rm th}$, a function of the metallicity of the ISM. Energy dissipation associated with star formation feedback are likely to be more significant when the metallicity is sufficient for cooling from metal lines to dominate over the cooling contribution from H and He. The transition of outflowing gas in the simulations is expected to occur at $Z\sim 0.1Z_{\odot}$ for a temperature range $10^{5}\,K<T<10^{7}\,K$ (Wiersma et al., 2009a). This phenomenon can be numerically depicted by equation 2, but only after replacing ($T_{\rm DM}$,$n_{\rm T}$,$10^{5}\;\rm K$) with ($Z$,$n_{\rm Z}$,$0.1Z_{\odot}$) to obtain,

[TABLE]

where $Z_{\odot}=0.0127$ is the solar metallicity and $n_{\rm Z}=n_{\rm T}=2/\ln{10}$.

2.3.5 Reference (FB$Z\rho$)

The feedback associated with FB$\sigma$ and FBZ becomes numerically inefficient in the centers of high-mass galaxies because a significant fraction of the star particles form at densities more than the resolution-dependent critical density ($n_{\rm{H},t_{\rm c}}$) above which radiation loss of the feedback energy is quick (Dalla Vecchia & Schaye 2012). These spurious energy losses can be partly compensated when a density dependence is introduced in the expression for $f_{\rm th}$:

[TABLE]

where $n_{\rm H,birth}$ is the gas particle’s density at the time when it gets converted into a star particle. Hence, at a fixed metallicity $f_{\rm th}$ increases with density. Such a density dependence may have a physical basis, because the star formation law and hence the feedback energy injection rate per unit volume, has a supra-linear dependence on surface density, which may result in smaller radiative losses at higher densities. In this work we use both the 50 and 100 cMpc boxes of the Reference model. The 100 cMpc box has a much larger number of massive galaxies for comparison to strong lens observations, whereas we use the Reference-50 boxes to compare with other model variations.

2.4 Variations of the reference model

Schaye et al. (2015) demonstrated that it is possible to calibrate the Reference model to reproduce the Galaxy Stellar Mass Function (GSMF) and the observed sizes (in different bands) of galaxies at $z$ = 0.1. However, a systematic study of the model’s key sub-grid parameters and sensitivity of this model to the variations of sub-grid parameters are critical. In order to quantify these effects, Crain et al. (2015) conducted a series of simulations (listed in the lower section of Table 1) for which the value of a single parameter was varied from that adopted in the Reference model. Here, we briefly summarize the five Reference model variations that are used in this work. There are five more Reference-model variations available, but those have a smaller box size (25 cMpc) that provide insufficient numbers of high-mass galaxies for comparisons to observed strong lens galaxies.

2.4.1 Viscosity variations (ViscLo and ViscHi)

The viscosity parameter $C_{\rm visc}$ governs two important physical processes: (a) the angular momentum scale at which gas accretion onto BHs switches from the relatively inefficient viscosity-limited regime to the Bondi-limited regime, and (b) the rate (only during the viscosity-limited regime) at which gas transits through the accretion disc (Rosas-Guevara et al., 2015). It is important to note that in both cases (viscosity-limited and the Bondi-limited regime) are subjected to the Eddington limit. A lower value of the viscosity parameter $C_{\rm visc}$, corresponding to a higher sub-grid viscosity. When the sub-grid viscosity is high, an earlier onset of the dominance of AGN feedback is triggered at a larger energy injection rate during the viscosity-limited regime. The viscosity parameter could thus affect the efficiency of galaxy formation and the scale of the halo mass at the peak of the stellar fraction. Lower (higher) values for the viscosity increase (decrease) both of them. However, we note that Bower et al. (2017) showed that the transition from slow to fast black-hole growth, which leads to the quenching of star formation, occurs when the halo is sufficiently massive to make stellar feedback inefficient and depends only very weakly on $C_{\rm visc}$.

2.4.2 Temperature variations in AGN heating (AGNdT8 and AGNdT9)

Schaye et al. (2015) have examined the role of the AGN heating temperature in EAGLE by adopting $\Delta T_{\rm AGN}=10^{8.5}{\rm K\;and\;}10^{9}{\rm K}$. They demonstrated that a higher heating temperature produces less frequent but more energetic AGN feedback episodes. They concluded it is necessary to reproduce the gas fractions and X-ray luminosities of galaxy groups. Le Brun et al. (2014) also concluded that a higher heating temperature yields more efficient AGN feedback. We analyze two Reference-model variation simulations with $\Delta T_{\rm AGN}=10^{8}\rm K$ (AGNdT8) and $\Delta T_{\rm AGN}=10^{9}\rm K$ (AGNdT9), besides the Reference model itself which adopted $\Delta T_{\rm AGN}=10^{8.5}\rm K$. In massive galaxies, the heating events (less frequent but more energetic) are more effective at regulating star formation due to a higher heating temperature. AGNdT8 (AGNdT9) model has higher (lower) peak star fraction compared to the Reference model. The reduced efficiency of AGN feedback, when a lower heating temperature is adopted, leads to the formation of more compact galaxies, because gas can more easily accrete onto the centers of galaxies and form stars.

2.4.3 No AGN feedback (NOAGN)

The final model that we consider has no AGN feedback and is the most extreme EAGLE model variation for massive galaxies. It appears unrealistic because the lack of AGN feedback is expected to dramatically increase the baryon concentration in the inner regions of galaxies, producing overly massive and concentrated galaxies. The reason that this variation is included, is to clearly demonstrate the effect of the absence of AGN activity. All other parameters are kept the same as in the Reference run.

3 Creating Mock Lens Data

Here, we explain the SEAGLE (Simulating EAGLE LEnses) pipeline in more detail. We briefly summarize the selection criteria of the (lens) galaxies, the extraction of the galaxies from the simulations, the impact of projection on the lens galaxy convergence map (Section 3.1), ray-tracing with GLAMER to create mock lensed images (Section 3.2), and finally the automatic process to create masks around the lensed images used in the lens modeling (Section 3.3). The flow diagram shown in Figure 1 of M18 describes the SEAGLE pipeline and the resulting data products. The reader is referred to M18 for more details on the pipeline.

3.1 Galaxy selection and post-processing

The initial down-selection of (lens) galaxies is based on the lens redshift ($z_{\rm l}$) and stellar mass ($\rm M_{\star}$) range from SLACS. Auger et al. (2010a) find a broad lens redshift range of $0.075<z_{\rm l}<0.513$ and a lower limit on the total stellar mass of $\rm M_{\star}\geq 1.76\times 10^{10}\;M_{\odot}$. The luminosities and effective radii of SLACS lens galaxies are based on a de Vaucouleurs profile fit to the galaxy brightness distribution as observed with Hubble Space Telescope (HST). We choose their I-band filter value, assuming it is closest to the bulk of the stellar mass. These are turned into stellar masses assuming either a Chabrier or Salpeter stellar IMF (Salpeter, 1955). We use the former in this paper to remain consistent with EAGLE. We also use a lower limit on both the line-of-sight stellar velocity dispersion ($\sigma>120{\rm\leavevmode\nobreak\ km\,s}^{-1}$) and the half stellar mass radius ($\mathrm{R_{50}}>1{\rm\leavevmode\nobreak\ kpc}$) from the EAGLE snapshot catalogs to avoid blatant outliers e.g., due to mergers. Table 2 summarizes these initial selection criteria.

We select all sub-haloes that match these selection criteria and extract all their particles from the snapshot. We do this for a single redshift roughly in the middle of the SLACS redshift range, i.e. $z_{\rm l}=0.271$. We reiterate, as in M18, that the lens redshift is fixed at $z$=0.271 for all mock lenses, despite having a range of observed lens redshifts. This redshift is intermediate between that of SLACS at somewhat lower redshifts, and SL2S plus BELLS at somewhat higher redshifts. Choosing simulation boxes at different redshifts for all lenses, to account for the minor effect of evolution, is computationally not feasible. We expect the effect of evolution to be small around this redshift (Furlong et al., 2015, 2017) and to be smaller than the observed scatter in the inferred quantities for all galaxies. Although this neglects the effect of evolution in the simulated sample, this redshift is roughly in the middle of the bulk of the redshifts of the combined set of SLACS, BELLS and SL2S lenses. For more details on the galaxy extraction we refer to Section 3.2 of M18. We finally rotate the particle position vectors in several directions around the center of the lens galaxy. In the current paper, each galaxy is projected along the three simulation box axes. The particles using the same SPH kernel as used in the simulation are exported into projected surface density maps (for more specifications see Trayford et al. 2017). For each galaxy, we separately calculate the surface density maps for the individual particle types (DM, stars and gas), as well as their total surface density map. Stellar remnants are included in the star particles.

3.2 Creating mock lens systems

The surface density maps are created in units of solar masses per pixel on a square-pixel grid of 512 $\times$ 512 (Table 2). They form the input to the ray-tracing lensing code GLAMER (Metcalf & Petkova 2014; Petkova et al. 2014). The size in proper kilo persec (pkpc) (100 pkpc) and pixel scale ($\approx$0.2 pkpc) of the grid ensure that the surface density map and the corresponding convergence map are well resolved in the inner regions of the galaxy (see Tagore et al. 2018), down to the simulation softening length, and are consistent with the SLACS pixel scale of 0.05 arcsec at $z$=0.271, corresponding to $\approx$ 0.2 pkpc.

We choose a source redshift for GLAMER to convert these mass maps into convergence maps, by dividing the surface density maps by the critical surface density which is set by the lens and source redshifts (Meylan et al., 2006). We choose a fixed redshift of $z_{\rm s}=0.6$, typical for SLACS lenses. Similar to the lens redshift, we choose a fixed source redshift to reduce computational overhead, although this restriction can be let go in the future. The dependence of the Einstein radius on source redshift is weak, however, increasing by $<20\%$ from $z_{s}=0.6$ to 1.0. Since all quantities in this work are determined inside fractions of the effective radius, the impact of the choice of the source redshift is very small. To describe the source, we use an elliptical Sérsic brightness profile with an index $n=1$, apparent magnitude = 23 in the HST-ACS F814W filter (AB system), an effective radius of 0.2 arcsec, a position angle $\phi_{s}=0\leavevmode\nobreak\ \deg$, and a constant axis ratio $q_{s}$=0.6. We set the parameters as such to keep close resemblance to sources found in SLACS (see Figure 4 in Newton et al. 2011). As shown in M18, the choice of the source size has negligible influence on the quantities of interest in this analysis. Furthermore, in Section 4 of Tessore et al. (2016), it is shown that there is only a negligible impact on the recovered parameters when using a realistic source as opposed to using a pixellated or parametric source model. They also show that LENSED recovers the source parameters well for both an exact model (i.e. the truth is part of the model family) and an inexact model. Thus our constant-size analytic source model is expected to have a negligible impact on our conclusions related to the mass density slopes as is further motivated in Appendix B.

For each convergence map, the critical curves and caustics are calculated, using GLAMER. We then randomly put the Sérsic source inside the diamond caustics of the lens to create multiple lensed images. This helps to maximize the number of arc and ring-like systems in the simulations (this roughly mimics the large magnification bias in the observations). The pixel scale of the grid – representing the lensed images – is set to 0.05 arcsec with the PSF and noise corresponding to an HST-ACS F814W exposure of typically 2400 s. The final resulting images have sizes of 161 $\times$ 161 pixels with 8.0 arcsec side length. The images are exported in standard fits-file format. All parameter values are again listed in Table 2 and motivated mostly by the range of values inferred from SLACS lenses (e.g. Koopmans et al., 2006; Newton et al., 2011; Bandara et al., 2013).

3.3 Mask creation

To mask large areas of noisy pixels in the image and include only regions around the lensed images in the lens modeling (see Figure 4 in M18), we automatically create a mask for each lens system. The noisy lensed images are convolved with a Gaussian having FWHM of 0.25 arcsec to decrease the noise by about a factor of 5 and obtain a slightly larger footprint of the lensed images. A surface brightness threshold is set at typically 2.5–5 times below the original noise level. This threshold defines the edge of the mask, faithfully traces the lensed images below the noise, and sufficiently extends outside the lensed images to include some noise-dominated pixels in the original image (see e.g. middle panel of Figure 4 in M18). The central 7$\times$7 pixels of images (i.e. centered on where the lens galaxy is) are also masked, similar to what is done for real lenses. This removes any artificially bright central images that are purely the result of a too low central convergence due to the small, but still finite, size of the SPH kernel. Whereas in real lenses the central surface density in general is extremely high (i.e. leading to large gradients in the potential), thereby de-magnifying the central lensed image, in the mock lenses it leads to a too bright central image. To avoid a bias in the lens model, we mask this central region. This artificial core has however little impact on the outer images near the Einstein radius. The resulting mask is used in all subsequent modeling and only image data inside the mask are used for the lens modeling.

4 Modeling of the Lens Sample

In this section we describe the selection of the final mock lens sample (Section 4.1), and the subsequent gravitational lens modeling and convergence-map fitting, i.e. the modeling of the surface mass density as directly obtained from the simulations (Section 4.2).

4.1 The lens sample

Implementing an automated recipe for the lens modeling of galaxies with stellar masses $\rm M_{\star}<10^{11}\;M_{\odot}$ has proven difficult due to the finite resolution effect of the particles during projection causing an artificial ‘core’ in the inner density profile, which in turn creates prominent but artificial images in the central regions of the lenses during ray tracing. These artificial images are not observed in real lens systems and are particularly pronounced in lower-mass galaxies that are more affected by the finite resolution of the simulations. As in M18, we therefore restrict ourselves to galaxies with total stellar masses $\rm M_{\star}>10^{11}\;M_{\odot}$. These galaxies are far less affected by any resolution effects and still significantly overlap with the massive lensing galaxies of SLACS and SL2S. Moreover, the disc to total ratio (D/T) distributions also matches well between SLACS and EAGLE (Reference-100) and thus we should statistically select comparable ETG candidates. Of these massive galaxies, about 80% are central galaxies (the most massive subhalo of a given halo) and about 20% are satellites (subhalos other than the main subhalo) in the 100 cMpc box. For the 50 cMpc boxes they are mostly (>90%) central galaxies. Table 2 summarizes the selection of this restricted and more massive sub-sample, used for all comparisons with observed lenses in this work. Table 3 lists the total mass density slope (t) values and the effect of the selection bias that we introduce via a mass-weighting scheme. Table 4 lists the average Einstein radius and several other relevant quantities of the simulated strong lenses from the different model variations of EAGLE.

4.2 Gravitational lens modeling

Having created the mock lens systems, we model each lens system with the open source, publicly available lens modeling code: LENSED (Tessore et al. 2016; Bellagamba et al. 2017). We use either an Elliptic Power Law (EPL; Tessore & Metcalf 2015) or a Singular Isothermal Ellipsoid (SIE; Kormann et al. 1994) mass model, including external shear. We use the corresponding mask, noise level and PSF for each system. A total of 14 or 15 parameters are sampled using a Nested Sampling MCMC method for the SIE or EPL models, respectively. The EPL mass model (which includes the SIE) has been utilised in several previous studies and has proven to describe very well the underlying mass model of strong gravitational lenses in various observational studies (Treu & Koopmans 2004; Koopmans et al. 2006; Koopmans et al. 2009; Barnabè et al. 2009; Newton et al. 2011; Barnabè et al. 2011). When modeling with an SIE plus external shear, we use the prior settings tabulated in Table 3 of M18. The SIE model’s (dimensionless) surface mass density can be numerically stated as:

[TABLE]

where $b$ equates to the measured radius of the Einstein ring (formally only for $q$=1) and $R$ is the elliptical radius defined by: $R=\sqrt{qx^{2}+y^{2}/q},$ where $q$ is the axis ratio (minor over major axis) and $x$, $y$ are Cartesian coordinates on the image plane. The lens is allowed to vary in position angle and mass centroid as well. We perform the lens modelling on the lenses with an EPL mass model. From Tessore & Metcalf (2015) we write the convergence as

[TABLE]

where 0 < $t_{L}$ < 2 is the power-law surface mass density slope. The other parameters (e.g. ellipticity, position angle etc) are kept same as for the SIE model. EPL can emerge from an oblate 3D density distribution, with $\rho(r)\propto r^{-t}$, where $t=t_{L}+1$. Both models also include external shear parameters. Statistically we aim to compare the SLACS, BELLS and SL2S lenses with those from the simulated lenses via the ensemble of density slopes obtained from the EPL implemented lens-modelling technique.

However, many of the SLACS density slopes were obtained from a joint lensing and dynamics analysis, rather than only from lensing (Koopmans et al. 2009; Barnabè et al. 2009; Auger et al. 2010b; Barnabè et al. 2011). We assume here that there is no significant bias between the lensing and lensing plus dynamics analyses (Tortora et al. 2014a; Xu et al. 2017). A direct comparison of the model parameters with the convergence map fitting can be performed with the same model, which we do not discuss further in this work but was extensively studied in M18. As in the creation of the mock lenses, we use a Sérsic profile to model the source. Even though some of the SLACS, BELLS and SL2S sources show irregular morphologies, our main objective is to calculate the global properties of the galaxies acting as lenses, and the exact choice of the source model does not bias the lens parameters for different (and inexact) source models (see section 4.4 of Tessore et al. 2016). We also compare the recovered source size between SIE and EPL and found negligible difference that does not bias our results (see Appendix B). In Figure 11 we also demonstrate that there is no such correlation between source-size density slope which might bias our analyses. Additional tests were carried out in M18, where we found no change in the distribution or the value of the model parameters when changing the source model parameters between lens systems (see Appendix A of M18). The priors used in the lens and source modeling are listed in Appendix A (see also Table 3 of M18). The priors were chosen such that the convergence of lens modeling parameters occurs faster in the Nested Sampling optimization and leads to minimal biases. We note that the priors are generally much wider than the inferred errors, hence they mostly guide the convergence rather than impact the parameter errors. The overall modelled parameters give considerably good fit to the lens and optimised residuals (for details see SEAGLE-I).

5 Observations

Here we summarize the strong lensing observational surveys that we use to compare with our results. In Section 5.1 we briefly describe the observations. Section 5.2 describes the weighting scheme that is used to compare simulated lens ensemble properties with observation. We note that in our comparison between simulated and observed lenses, we show all of the SLACS lens galaxies for visual purposes, but only quantitatively compare these galaxies with simulated galaxies for the restricted range $\rm M_{\star}>10^{11}\;M_{\odot}$.

5.1 SLACS, SL2S and BELLS

In the SLACS survey, Bolton et al. (2006) selected potential lens candidates spectroscopically from SDSS. Since then SLACS has successfully identified more than a hundred confirmed strong lens systems, with HST follow-up. The SLACS galaxies are massive ETGs, specifically Luminous Red Galaxies (LRGs) with star-forming background sources emitting strong emission lines. The advantage of the SLACS survey is that for all lenses spectroscopic lens and source galaxy redshifts are available. The mean Einstein radius of SLACS lenses is 1.2 arcsec (Koopmans et al., 2006; Auger et al., 2010a) with sources having a typical size of about 0.2 arcsec (Koopmans et al., 2006; Newton et al., 2011) and typically being at $z_{\rm s}\approx 0.6$. Although it is the largest complete strong lens sample, SLACS has a relatively limited lens redshift range with the bulk of the lenses in the range of $z_{\rm l}\approx 0.1-0.3$.

The SL2S survey was initiated to increase the number of known lenses by a different methodology than SLACS. In SL2S, Cabanac et al. (2007) performed a dedicated search in the CFHTLS to find strong gravitational lenses. They focused on mostly galaxy-scale and group-scale lenses. SL2S aimed at providing a larger sample of strong lenses at higher redshift. RingFinder (Gavazzi et al. 2014), an automated software was used in SL2S to find lenses by searching 170 square degrees of the sky. RingFinder performed a search for blue arcs that are elongated tangentially and ring-like structures around red galaxies to select lens candidates. The most promising systems were followed up with HST and spectroscopy (Gavazzi et al., 2012). Even though SL2S lenses combined with SLACS provided evidence of structural evolution (Ruff et al. 2011), the SL2S sample is limited by a lack of source-galaxy redshifts for a considerable number of systems.

In BELLS, Brownstein et al. (2012) utilized the same spectroscopic methodology implemented in SLACS, to select the strong lenses, but they used Baryon Oscillation Spectroscopic Survey (BOSS; Eisenstein et al. 2011) spectra. BELLS discovered a sample of strong galaxy-galaxy lenses, at somewhat higher redshifts, that is of comparable size and homogeneity to that of SLACS at lower redshift. BELLS is also comparable in stellar mass to the SLACS lens galaxies. Both the BELLS and SLACS samples are complete in their spectroscopic lens and source galaxy redshifts. The lens redshifts of the three lens samples are within a similar range of 0.1–0.65, but the source redshifts cover a wide range from 0.3 to 3.5. Bolton et al. (2012b) have reported evidence for mild evolution in the mass density slope between BELLS and SLACS. We ignore this in the sample of mock lenses and compare observations with simulations only at $z=0.271$, in between the two samples, as discussed earlier.

5.2 Lens selection bias

Differences in lens-galaxy selection and follow-up can lead to differences in the population of lenses in the SLACS, BELLS and SL2S samples. For example, due to the relatively small fiber opening used in SDSS spectroscopic observations ($1.5^{\prime\prime}$ radius), the SLACS spectroscopic survey typically limits the search to lenses with an equivalent or smaller Einstein radius (although larger lenses could be found if one of the lensed images is inside the fiber and bright enough), and finite source effects play a role as well. SL2S on the other hand can select lenses directly from images and over a larger Einstein radius range, i.e. mass scale, typically yielding Einstein radii greater than $1^{\prime\prime}$, because they are less well resolved in ground-based images. These selection effects are hard to quantify though (see e.g., Dobler et al., 2008, for SLACS).

Observational selection biases often hinder a proper comparison between simulations and lens surveys, strong lensing being no exception. In this work, we assume that lens selection biases are dominated by the lens-galaxy mass and correlate sub-dominantly with the lens and source redshifts, and with the lens-galaxy mass density profile and ellipticity. This is a reasonable assumption if the lens mass models are close to isothermal (i.e., the caustics are shape invariant as a function of redshift and only scale in cross-section) and the source size is small compared to the Einstein radius (Dobler et al., 2008). Massive ellipticals also do not vary strongly in their ellipticity. The observed lens sample properties are then mainly affected by the lensing cross-section (Marshall et al. 2007), which is mass dependent, and by the magnification bias, which can be different between surveys. A precise analysis is difficult to implement and beyond the scope of this paper. We therefore take an empirical approach and only correct for the lens cross-section and we assume that the magnification bias does not correlate with galaxy mass111This holds exactly for SIE models if the source is a point source and the galaxy mass model (i.e. ellipticity for the SIE) does not vary with galaxy mass.. The square of the Einstein radius varies proportionally with the cross-section of lensing for the EPL model for a fixed ellipticity (generally also close to the SIE model). Assuming the Faber-Jackson relation (Faber & Jackson 1976) and a constant mass-to-light ratio, the Einstein radius is again proportional to the stellar mass of the respective galaxy. Hence we arrive at a direct observable (i.e. the stellar mass) in both the simulations and observations.

Motivated by the above arguments, we propose the following weighting scheme per lens:

[TABLE]

with $\langle M_{\star}\rangle$ being the average stellar mass of the galaxies in the sample and $\alpha$ being the exponent of the weight function. We re-weight each simulated strong lens (which we assume to be volume limited) when comparing distributions (i.e. histograms) of the mass-density slopes between observed lenses from SLACS, BELLS, SL2S and simulated lenses. Hence a weight $W_{i}$ for simulated lens $i$ implies it counts as $W_{i}$ galaxies (note that the weights are non-integers). Most of the lenses are massive systems, and in general drawn from the exponential tail of the mass function. Hence re-weighting should have a limited impact on the massive end of the distribution functions, but it does strongly affect the low-mass end. We test values of $\alpha=0.5,1.0$, and 1.5 to show that the weighting scheme does not affect the conclusions and are only to mimic the observation selection bias of the lenses depending on their stellar mass. Other options for re-weighting the lens galaxies, to account for their lens cross-section, are using either their Einstein radii or their stellar velocity dispersions, which we have not done in this work and leave for future improvements in the analysis when we study the redshift evolution of these lenses.

6 Results

In this section we compare the simulated EAGLE lenses with those from SLACS, BELLS and SL2S, in terms of their surface mass density profiles. Even though SL2S and BELLS lenses are typically at somewhat higher redshifts, we compare the simulated lenses at $z_{\rm{l}}$=0.271 assuming limited ETG evolution within the redshift range 0<$z$<1, as discussed earlier. This assumption is reasonable as it was pointed out by both Sonnenfeld et al. (2013b) and Koopmans et al. (2006), that the total mass density slopes (which are close to isothermal) do not strongly evolve with time in observed ETG lenses (although see Bolton et al., 2012a). We compare the mass-size relation in Section 6.1, the total density slopes in Section 6.2, and the Einstein radius in Section 6.3. We compare our results with OWLS simulation in Section 6.4. Table 2 summarizes the number of galaxies, lenses and projected mass maps. Tables 3, 4 and 5 give the effect of magnification bias (mimicked by a weighting scheme) on the total mass density slope ($t$) values, the average Einstein radii, the ratio of Einstein radius over effective radius and several other relevant quantities of the simulated strong lenses from different model variations of EAGLE.

6.1 The mass-size relation

Observationally the stellar mass (or luminosity to be precise) of an ETG correlates with its size (e.g. Baldry et al. 2012). Similarly, in our simulations the stellar masses of galaxies correlate with their sizes (Furlong et al. 2017). To assess whether a similar relation holds for the mock lenses at $z_{\rm l}=0.271$, we define the effective radius ($R_{\rm eff}$) as the stellar projected half-mass radius in the simulations, hence assuming a constant mass-to-light ratio. As demonstrated by Remus et al. (2013, 2017), this might lead to a slight overestimation of the actual size of the galaxy compared to observations (e.g. in the case of SLACS the effective radius is derived from a de Vaucouleur fit to the galaxy brightness distribution), but we ignore this minor (<0.05 dex) effect rather than fit a profile to the projected stellar mass for all simulated galaxies. We assume a constant Chabrier IMF for both the observed and simulated galaxies.

Figure 2 shows the mass-size relations for the nine selected EAGLE model variations, overlaid on SLACS. We find that the Reference model (REF) which was calibrated on the GSMF and galaxy sizes, yields somewhat larger effective radii compared to similarly massive SLACS galaxies. On the other hand, the models FBconst, FB$\sigma$ and FBZ, which (except for FBconst) were calibrated on the GSMF but not on galaxy sizes, all have similar effective radii as SLACS, except for two outliers around the lowest stellar mass end and above the relation that have unusually large effective radii222We note that each mock lens is shown three times (once for each principle-axis orientation), as discussed earlier, and hence the number of independent outliers is very small.. Due to the relatively low efficiency of stellar feedback in the FBconst, FB$\sigma$, FBZ models and the absence of AGN feedback in the NOAGN model, stars tend to form somewhat closer to the center of the galaxy (see Crain et al. 2015). The NOAGN model, however, leads to much more compact galaxies, with some systems even straddling the resolution limit of the simulations. The galaxies from the AGN model variations (AGNdT8 and AGNdT9) both have larger effective radii than the NOAGN model. When $\Delta T=10^{8}\rm K$ (AGNdT8) about half of the galaxies are more compact in size and in good agreement with SLACS, whereas for $\Delta T=10^{9}\rm K$ (AGNdT9) hardly any galaxy matches the observations. The higher temperature in the AGNdT9 model leads to more effective AGN feedback, keeping gas away from the center and increasing the size of the galaxy. For comparison, the Reference model assumes $\Delta T=10^{8.5}\rm K$, explaining its position halfway between AGNdT8 and AGNdT9 in mass-size relation (Figure 2). A low black hole accretion disc viscosity (ViscLo), i.e a high viscosity parameter ($C_{\rm visc}$), delays the onset of AGN feedback, allowing gas to settle closer to the galaxy center before star formation. The ViscHi model has the opposite effect, increasing the size of the galaxy.

Overall, we conclude that simulated galaxies from EAGLE better match the mass-size relation of SLACS lens galaxies when there is moderately low AGN activity or stellar feedback driving the galaxy formation, with only a mild impact from variations in the type of stellar feedback model. This trend is consistent with previous studies (Remus et al. 2017; Figure 1 in Peirani et al. 2018). Finally, we find that changes in the viscosity have a stronger impact by indirectly affecting AGN feedback.

6.1.1 Comparison with earlier EAGLE results

We now compare the inferred mass-size relation with the results by Schaye et al. (2015), Crain et al. (2015) and Furlong et al. (2017). This comparison is necessary to assess any selection bias within the simulations. If we are selecting an ETG population that is significantly different than the total galaxy sample analyzed in other aforementioned EAGLE works, this might invoke a bias in our lensing ETG sample and their properties. Moreover, our calculations are performed on mass maps and not directly on the cataloged particles. Schaye et al. (2015) compared the mass-size relation of the Reference model by Shen et al. (2003) and Baldry et al. (2012), and found excellent agreement. Similarly, Crain et al. (2015) compared the mass-size relation from the calibrated models (Figure 3 therein) and found $\approx$0.2-0.3 dex difference from the Ref-50 model at the higher-mass end. This result is consistent ($\approx$0.2 dex difference) with our findings in Figure 2 for our strong lensing sample. Figure 11 and 12 from Crain et al. (2015) (3rd panel from right) show a comparison of mass-size relation of Ref-50 model variations, from which it is concluded that AGNdT9 and ViscHi models yield larger galaxy sizes compared to the AGNdT8 and ViscLo models, respectively, consistent with our findings. In Figure 3, we compare the mass size relation of the Ref-100 cMpc model obtained in our analysis with that by Furlong et al. (2017). We find excellent agreement, within 0.1 dex. We also compare with Shen et al. (2003) and Baldry et al. (2012) and found that our results are consistent with them. The mass-size trends in this paper are thus consistent with the findings of other EAGLE studies showing no bias due to our selection or methodology in calculating the sizes. As for EAGLE, part of the difference lies in the fact that the Ref-50 simulations provide larger sizes than the Ref-100 simulation at $M_{\star}>10^{11}\rm M_{\odot}$, due to small number statistics for Ref-50 (see also Crain et al. 2015). However, some systematic differences are still present with strong-lens galaxies tending to be more compact than non-lensing galaxies. SLACS galaxies, therefore, appear about 0.2 dex smaller in size than non-lensing galaxies of similar mass (see right panel of Figure 4). In paper III of SEAGLE series we will explore the systematics and compare with non-lensing ETGs from SPIDER survey (La Barbera et al., 2010; Tortora et al., 2014a), which we will show, have sizes that agree much better with EAGLE, and we point out the methodological differences (e.g., measurements with different bands of observations, different fitting algorithm, etc.) that could potentially bias the analysis.

6.2 The total mass density slope

Keeping the mass-size results discussed in the the previous section in mind, in this section we assess whether the same galaxy formation models that (visually) reproduce the mass-size relation of SLACS lens galaxies also reproduce their mass density slopes, which is not an observable used in the calibration of the EAGLE simulations. We use the EPL surface mass density profile to model the simulated strong lenses with LENSED, closely mimicking real lens observations (see Section 4 for details). This allows for a more unbiased and systematic comparison with strong lens observations.

6.2.1 Calibrated simulations

As a first check, we confirm that the lens galaxies from the Reference-100 cMpc model show a similar distribution of density slopes as presented in M18 where the smaller 50 cMpc box was used. The latter has a much smaller number of massive galaxies. We confirm that EAGLE galaxies from the Reference model tend to have slightly steeper density slopes than SLACS, BELLS and SL2S (see left panel in Figure 4 and also Figure 12 in M18). However, the ratio of $R_{\rm Ein}$/$R_{\rm eff}$ can play a crucial role in this respect because the lens modeling is mainly dependent on information obtained near the Einstein radius. Since the total mass density can be sensitive to the radial scale at which it is measured (Xu et al. 2017), we will explore this aspect in Section 6.3.

In Figure 5 the density slopes for all EAGLE model variations are shown for the smaller 50 cMpc boxes. The FBconst model appears to yield galaxies most similar to SLACS with the total mass density profile being very close to isothermal. This can be attributed to its less efficient stellar feedback, which yields a mass profile, different than the Reference model. The FBZ and FB$\sigma$ models have more dark matter in the center of the galaxy compared to the FBconst and Reference models, leading to a shallower total density slope in their central regions. Hence, whereas the FBZ and FB$\sigma$ models visually reproduce the mass-size relation of SLACS rather well, they fail to reproduce their mass density slopes. We find the rather counterintuitive trend that when feedback efficiency increases from the FBZ, FB$\sigma$, FBconst to Reference models, the average total mass density slope steepens. We will see that variations in AGN feedback show the same trend and we will discuss the cause in the next section. In Section 6.3 we will also study the correlation of the ratio of $R_{\rm Ein}$/$R_{\rm eff}$ with the total mass density slope for different model variations.

6.2.2 Reference-model variations

There is a clear dependence of the total mass density slope on AGN feedback. As the stochastic temperature increment in AGN models increases from $\Delta T=10^{8}$K (AGNdT8) to $\Delta T=10^{8.5}$K (Reference) and $\Delta T=10^{9}$K (AGNdT9) the total density slope steepens. Generally, we would expect the opposite, since stronger AGN activity (i.e. temperature increments) should move or keep gas particles away from the galaxy center, preventing star formation. As mentioned in Le Brun et al. (2014), more energetic heating events associated with a higher heating temperature, even-though less frequent, are more effective at regulating star formation in massive galaxies. Crain et al. (2015) also pointed out that the peak galaxy formation efficiency decreases with increasing AGN temperature. The reduced efficiency of AGN feedback thus, counter-intuitively, manifests itself in a steeper total mass density slope. A similar trend is found when the viscosity parameter is increased, which impacts AGN feedback at fixed mass as discussed earlier. This trend is consistent with previous simulation studies (e.g. Remus et al., 2017; Xu et al., 2017). In short, the AGNdT8 model with its weaker AGN feedback (compared to the Reference model) produces lensing galaxies that are closer to isothermal and in better agreement with the results from SLACS, BELLS and SL2S lens galaxies. Table 5 summarizes the mean, median and standard deviation of the density slopes for all EAGLE model variations used in this work. The evolutionary trends will be studied in details in a forthcoming paper.

6.2.3 Correlations of slope and galaxy stellar mass

We correlate the total mass density slope and the stellar mass of the three prominent simulation models compared in our analysis, namely, Ref-50, FBconst and AGNdT8. Figure 6 shows the distribution of the density slopes calculated from lens modeling from both simulations and SLACS (Koopmans et al., 2009). We find at most a very mild trend in the total mass density slope with the stellar mass, consistent with strong lensing observations of SLACS (Koopmans et al., 2009). More massive galaxies tend to have a slightly lower total density slope than less massive galaxies in all three model variations (see also Tortora et al. 2014a, where this trend, with shallower (isothermal) profile at high mass and steeper profiles at lower masses are found). However, the intrinsic scatter in the distribution in each of the model variations is too large to draw any significant conclusion, especially since the high-mass end of the distribution contains very few galaxies in the simulations. This very mild trend is also consistent with theoretical work by Remus et al. (2017) and Xu et al. (2017).

6.2.4 Dependence on weighting scheme

We test different values of the $\alpha$ parameter in our weighting scheme to demonstrate the robustness of our results against the selection effects in the observations. In Figure 5, we show the variations in the median total mass density slope for three different values of $\alpha$=0.5,1.0, and 1.5. Although, the median density slope is sensitive to the weighting scheme, relative changes are well within the spread calculated for each of the model variations. This result implies that our conclusions do not strongly depend on the observational selection bias. We note that we do not separately compensate for the magnification bias, as a function of galaxy mass, but assume this effect is folded into the weighting scheme. The at-most mild trend of the density slope with galaxy mass, however, suggest that any re-weighting based on galaxy mass will make little difference in the conclusions. Tables 3 and 5 list the median values of the total mass density slope for different values of $\alpha$ parameter, and their relative change compared to the nominal model with $\alpha=1$. We note that we have not considered the errors on the measured slope in Figure 5. The errors on the measured slopes will slightly broaden the distributions. However, the rms error on the slopes is typically well below 0.2 (see Auger et al. 2010b), i.e. inside our chosen bin-size, and considerably smaller than the spread in the distribution. In addition, the slope measurements from the simulations have a similar spread, mimicking partly this broadening effect, thus reducing its impact. The changes in galaxy-formation processes is by far the most prominent source differences in the distributions.

6.3 Einstein radius comparison and correlation with the total mass density slope

The Einstein radius ($R_{\rm Ein}$) is a fundamental observable in strong gravitational lensing. However, to compare between strong lenses simulated from different model variations of EAGLE having a range of effective radii and strong lensing surveys having different observing strategies (e.g. SLACS, SL2S and BELLS), we need to compare the ratio of $R_{\rm Ein}$/$R_{\rm eff}$ (see Li et al. 2018). For SLACS, the values of $R_{\rm Ein}$/$R_{\rm eff}$ ratios populate $\approx$ 0.7 with very little scatter due to the small redshift-range for both the source and the lens (Koopmans et al. 2006; Koopmans et al. 2009). Whereas SL2S yields larger values of $R_{\rm Ein}$ with similar sized lensing ETGs as SLACS, due to the large spread in redshift-range of the lensing galaxies ($z_{\rm l}=0.2-0.8$) and the lensed sources ($z_{\rm s}=1-3.5$) (Sonnenfeld et al., 2013a). In BELLS, the $R_{\rm Ein}$/$R_{\rm eff}$ values mainly range from 0.5 to 1.5 with a sharp drop below 0.5, primarily due to a wide range of the source redshift from $z_{\rm s}$=0.8 to 3.5 with mean lens redshift of $z_{\rm l}$=0.52 (Li et al., 2018). We find that our best models, FBconst and AGNdT8, are closest in their $R_{\rm Ein}$/$R_{\rm eff}$ to the mean value of SLACS. Table 4 gives a complete overview of the mean of $R_{\rm eff}$, $R_{\rm Ein}$, the ratio $R_{\rm Ein}$/$R_{\rm eff}$ and their respective rms values for different model variations of EAGLE and observations (e.g. SLACS, SL2S and BELLS).

Figure 7 shows the correlation between the average total mass-density slope ($t$) and $R_{\rm Ein}$/$R_{\rm eff}$ ratios from different model variations of EAGLE. We find that as the feedback becomes stronger, the effective radius increases (consistent with Sales et al. 2010). This in-turn decreases the ratio $R_{\rm Ein}$/$R_{\rm eff}$ and steepens the total density slope since $t$ is calculated at the $R_{\rm Ein}$. The larger sizes of Einstein radius for strong lenses in SL2S are primarily due to the difference in observing strategy from SLACS. SLACS and BELLS selected lens candidates from spectroscopic signatures coming from two objects at different redshifts on the same line of sight in the SDSS spectra. The relatively small fiber used in SDSS spectroscopic observations, 1.5${}^{{}^{\prime\prime}}$ for SLACS and 1${}^{{}^{\prime\prime}}$ for BELLS in radius, they tend to select strong lenses with small Einstein radii. SL2S finds considerably more strong lenses with Einstein radii above 1${}^{{}^{\prime\prime}}$, since they can be more clearly resolved in ground-based images. For similar comparison of $R_{\rm Ein}$/$R_{\rm eff}$ in SLACS, BELLS and SL2S, readers are refereed to Figure 1 in Sonnenfeld et al. (2013a) and Li et al. (2018).

6.3.1 The assumption of a power-law density profile

Koopmans et al. (2006) tested the assumption of the shape of the density-profile itself, i.e. the power-law model. If the density profiles of lens galaxies are different from a power-law, but have the same shape for each galaxy (scaled to a common scale), or, if they are different from a power-law and different between lens galaxies, the power-law assumption might give biased results. In either case, it is expected that the inferred (average) logarithmic total mass density slope inside $R_{\rm Ein}$ will change with the ratio ($R_{\rm Ein}$/$R_{\rm eff}$) for a particular model variation. In the case of the total mass density slope is a broken power-law with a change in slope inside $R_{\rm Ein}$, one expects $t$ to change depending on where the change in slope occurs with respect to $R_{\rm eff}$. Thus one is expected to find some “average” slope weighed by the luminosity and kinematic profile, varying as function of ($R_{\rm Ein}$/$R_{\rm eff}$). This is due to the dependence of $R_{\rm Ein}$ mostly on the relative distances of the lens and the source and is not a physical scale of the lens galaxy itself. Koopmans et al. (2006) found no evidence of any clear systematic correlation between $t$ and $R_{\rm Ein}$/$R_{\rm eff}$ ratio (see Figure 5 therein). Figure 8 shows the trend in the total mass density slope and the ratio $R_{\rm Ein}$/$R_{\rm eff}$ for individual lenses. We also find no evidence of any correlation between $t$ and $R_{\rm Ein}$/$R_{\rm eff}$ ratio for both FBconst and AGNdT8 models, thus showing that our results are not biased by the power-law assumption. The small deviations of $t$ from 2.0 further support this. We conclude that our assumption of a single power-law shape for the total density profile is valid and reliable, consistent with the finding of Koopmans et al. (2006).

6.4 Comparison with OWLS simulations

In a previous study using five model variations from OWLS (Schaye et al. 2010) and also the DM-only simulation, Duffy et al. (2010) probed the mass density slope at $z$=2 and compared the results with SLACS lenses (Figure 3 therein). They found that implementation of AGN feedback, or extremely efficient feedback from massive stars, is necessary to match the observed stellar-mass fractions in groups and clusters. However, that made the inner density profiles shallower than isothermal. They concluded that a weak or no feedback produces galaxies with an isothermal profile. This is consistent with the results in this work, where we also found that weaker feedback leads to better agreement of the total mass density slope with SLACS, BELLS and SL2S observations. However, they also conclude that other observables, such as the stellar fractions, rule out those weak feedback models (e.g. see Crain et al. 2015). One way to explain this conundrum is that all the models miss something critical, which may well be the case. Another explanation could be that the strong lenses are a biased sample of the total ETG population in a volume limited sample. Previously, Sales et al. (2010) explored different feedback models in OWLS (Schaye et al. 2010) and found large variations in the abundance and structural properties of bright galaxies at $z$ = 2. They showed that models with inefficient or no feedback lead to the formation of overly massive and compact galaxies with a large fraction (upwards of 50 percent) of all available baryons (gas, stars, and stellar remnants) being retained in each halo. Increasing the efficiency of stellar or AGN feedback reduces the baryonic mass fraction fraction and increases the size of the simulated galaxies. This trend is also consistent with our findings.

The conclusion in Duffy et al. (2010) that NOAGN feedback produces an isothermal profile is in contradiction with our analysis. One reason could be that our analysis is carried out at a redshift of $z=0.271$, however, closer in redshift to where these lens galaxies are observed and is consistent with the results of several other simulation studies (Xu et al. 2017; Remus et al. 2017). Analysis of Duffy et al. (2010) is done at a significantly higher redshift of $z$=2. In the next section, we will discuss the possible reasons for these differences in light of potential systematics.

6.5 Potential systematics

There could be several effects that play a role in the comparison between observations and simulations. We describe three of these below.

Evolution of the density profile:

The inclusion of baryons results in differences in the total density profiles that depend on the efficiency of the radiative cooling and feedback. As pointed out in Remus et al. (2017) and Xu et al. (2017), there could be a significant steepening of the total mass density slope in the simulations at lower redshifts which might affect the density-slope analysis. Even though Koopmans et al. (2006) have shown that there is no strong evidence for evolution in the total mass density slope in SLACS with redshift, this only holds for the redshift range of $0.1\lesssim z\lesssim 0.3$ where the bulk SLACS lenses are found. Evolution might exist as we move to higher redshifts (Bolton et al., 2012a; Sonnenfeld et al., 2013b). This potentially could explain the differences between this work ($z$=0.271) and the analysis in Duffy et al. (2010) which was carried out at a higher redshift ($z$=2). Moreover, the galaxies analyzed in Duffy et al. (2010) are less massive than those used in our analysis, mostly due to the significant difference in the redshifts of both the analysis. Also, for a random lens system, we measured the density profile with the lensing galaxy at three lens redshifts of $z_{l}$=0.101, $z_{l}$=0.271 and z=0.474, with the source redshift remaining at $z_{s}$=0.6. We found the difference in the slope parameter to be 0.02 and 0.03 respectively i.e. much below the rms error. So we assume the effects are currently not significant in our case. A similar result is also reported recently by Wang et al. (2019) where they find the density slope to remain nearly invariant after $z$=1 with a mild increase towards $z=0$. However, in our case, a full-scale redshift evolution study is beyond the scope of this work.

Simulation resolution bias:

Duffy et al. (2010) found that the resolution of the simulations can strongly affect the region where the mass density slope is measured. Their density slope measurement, however, was typically done around an Einstein radius of $\sim 3$ kpc, only just above the resolution limit in the highest-resolution OWLS run at $z=2$. Similarly, Schaller et al. (2015) showed that below a radius of roughly $\sim$ 2-3 kpc, the matter density slope due to the resolution, is increasingly less reliable. This is not directly due to the softening length, but rather due to the radius enclosing a certain number of particles needed for the circular velocity to converge to within $\sim$10% (i.e.  the convergence radius) and the enclosed mass to within $\sim$20%. At radii smaller than the convergence radius, the mass profile becomes increasingly less reliable and typically displays a too shallow density profiles. The impact of baryons, especially a large number of stars dominating the potential in these regions, also becomes more uncertain. In our work, however, we analyze galaxies at much a lower redshift and at a much higher resolution, similar to Xu et al. (2017) and Remus et al. (2017) (i.e. to Illustris and Magneticum, respectively). In these lower-redshift and higher-resolution simulations, massive galaxies have a larger Einstein radius, in the range of 3-10 kpc, well above the resolution limit and also above the convergence radius in the simulations. We therefore expect these effects to play a minor role in the current EAGLE simulations around the Einstein radius of massive early types galaxies with $M_{\star}>10^{11}$ M⊙.

Intrinsic degeneracies:

There exists an intrinsic degeneracy between the source size and the density slope. This degeneracy is intrinsic to lensing and independent of the lens modelling technique adopted (see Appendix C & D). The errors on the slope have a dependence on the size of the source for that particular lens-system. However, we have simulated the lenses with exactly the same source parameters (see Table 2) for uniformity. Also, in SEAGLE-I, we have tested for a small sub-sample of simulated lenses with a range of source sizes and found no evidence that a varying source size has an influence on the main modelled quantities. We do note that the error on the inferred source depends on the noise and angular resolution of the data that the referee is suggesting. In a recent published study by Van de Vyvere et al. (2020), it is shown that size of the mass maps could also potentially affect the image reconstruction. They also tested for several realistic slope types and finally showed that density slope and the Einstein ring measurements remained completely unaffected regardless of the change in source orientation, lensing configuration and size (see Conclusions therein). Thus when a statistical inference is drawn from the analyses presented here, it remains robust when comparing these key quantities to SLACS/BELLS/SL2S. Moreover, comparison with the fitting results coming directly from the simulations demonstrates that the recovered values for slopes via lensing are trustful to draw conclusions of their feedback mechanism. Even though, an argument maybe made that SLACS lenses have a range of source sizes, Newton et al. (2011) showed along with comparison from the GEMS and HUDF samples that the values chosen in this analysis are realistic enough to simulate the lens sample regardless of that it is obtained from lenses or non-lenses. Nevertheless, as mentioned in Appendix B that the mean source sizes obtained are 0.218, 0.217 and 0.213 for Ref, FBconst and AGNdT8 repectively which corresponds to 9% , 8.5% and 6.5% mean errors respectively. Asserting that the assumption of an SIE model for the lens galaxy mass profile is the most significant source of systematic uncertainty, Marshall et al. (2007) estimated the systematic errors in the source size to be 12%. The mean errors in this analysis for key model variations are well within the expected limit for these surveys. Thus in either way the results presented in this work is robust in the scope of the parameter range that we chose to simulate the lenses.

Observational biases:

Dobler et al. (2008) found that the most significant instrumental selection effect is the finite size of the spectroscopic fiber, which selects against large separation lenses and results in a non-monotonic dependence of the rogue line probability (defined as the probability that a given luminous red galaxy (LRG) has a rogue [O II] line in its spectrum) on velocity dispersion. The situation is further complicated by the effects of atmospheric seeing, which can add flux from images outside or remove flux from images inside the fiber. Dobler et al. (2008) also reported that the lensing probability has a fairly weak dependence on the size of the source (see also the appendix of M18). Hence, whereas it is clear that lens galaxies are mass-selected and biased to higher-mass galaxies, some of the most massive lenses might not have been found in SLACS due to the above-mentioned effects. These massive systems are already rare to begin with and their absence would not bias the bulk of the lens population which peaks around $M_{\star}=10^{11.35}\rm M_{\odot}$ (Auger et al. 2010b). As was shown by Bolton et al. (2008b), SLACS lens galaxies also appear in all observational aspects to be similar to their LRG parent population, suggesting that they are not a biased LRG sub-sample. Also, BELLS is very similar to SLACS in the type of lens galaxy, given the more heterogeneous nature of the lenses and their environments in the SL2S survey (which were morphologically and not spectroscopically selected), a lesser agreement with SL2S is maybe not entirely unexpected. Nevertheless, previous observational (e.g. Auger et al. 2010b; Sonnenfeld et al. 2013b; Li et al. 2018) and simulation analyses (Xu et al. 2016; Remus et al. 2017) these surveys have been compared among each other with the assumption that different observational selections does not hinder a fair comparison.

Moreover, as pointed out in Tortora et al. (2014b) (Table 1), strong lensing galaxies tend to be more compact than non-lensing galaxies (e.g. SPIDER sample). However, SPIDER uses K-band data and Sérsic fitted values of $R_{\rm eff}$, while SLACS uses V-band and de Vaucouleurs fit. This can give different results. But Auger et al. (2010b) showed that using different fitting profiles gives negligible difference in $R_{\rm eff}$ values. Even though this is consistent with the argument that strong lensing prefers weaker feedback which in turn forms galaxies with relatively smaller sizes at fixed stellar mass compared with more efficient feedback models, it might bias correlations between galaxy properties. It could be that LRGs are a biased sub-sample of galaxies with respect to volume limited samples. We will explore this trend of galaxy sizes in light of dark matter fraction and explore possible systematics that might be causing the differences in a forthcoming work.

6.6 Kolmogorov-Smirnov statistics

Even though we find qualitatively and visually quite similar distributions between some of the model variations (i.e., FBconst, AGNdT8) and observations, we have not quantified this (dis)agreement. The Kolmogorov–Smirnov (KS) test (Kolmogorov, 1933) is a nonparametric test of the equality of continuous, one-dimensional probability distributions that can be used to compare a sample with a reference probability distribution, or be used to compare two samples. The KS statistics (D-value) quantifies the maximum probability difference between the cumulative probability distribution functions of two samples. A KS test also yields a p-value, being the probability that two distributions are in fact drawn from the same underlying distribution and are dissimilar at the current level (D) or larger, by random chance. In this work, we use the standard 1D KS test to compare the mass density slopes and we use the 2D KS test of Peacock (1983) to compare the mass-size relations. Table 5 summarizes the KS D- and p-values by comparing the results from the EAGLE model variations with those of SLACS, BELLS and SL2S, respectively.

We indeed find that the FBconst, AGNdT8 and ViscHi models which visually appeared most consistent with the observations, also have consistently high p-values (we assume a lower limit of acceptance of $p>0.05$). When we combine our analysis with the p-values from the 2D KS test for the mass-size relation, we find that only the FBconst and AGNdT8 model variations remain viable. The Reference model, even though displaying similarity to observations of the mass-size relation from SLACS, performs poorly in the mass density slope KS test. In addition, we can clearly rule out the NOAGN, ViscLo, FBZ and FB$\sigma$ model variations based on their failure to reproduce the observed strong lens distributions in slope, mass and size. This confirms our earlier visual inspection.

7 Summary and Conclusions

In this paper, we have systematically explored the impact of different galaxy formation processes used in the EAGLE hydrodynamical simulations – in particular stellar and AGN feedback – on strong lens observables in massive early-type galaxies with $\rm M_{\star}>10^{11}M_{\odot}$. Simulations of various mock-lens ensembles with the SEAGLE pipeline (M18) allow us to quantify in particular the (dis)agreement between the total mass density slopes around the Einstein radius and the stellar mass-size relation between these mock lens ensembles and observations from the SLACS, BELLS and SL2S lens surveys. We compared these observables with the outcome of a range of EAGLE model variations, varying stellar & AGN feedback and black hole accretion disc viscosity parameters (Schaye et al. 2015; Crain et al. 2015).

We select potential strong lenses based on the stellar mass ($\rm M_{\star}>10^{11}\rm M_{\odot}$) at a redshift of $z_{\rm l}=0.271$ and create projected mass maps for three different orientations. We create mock lenses by ray tracing through the mass maps, placing an analytic Sersic (1968) source with observationally motivated parameters at a higher redshift ($z_{\rm s}=0.6$). We add realistic HST noise and PSF to mimic strong lenses found in observations. We calculate the projected half-mass radius for each individual mass map. We also model these lenses with an elliptic power-law model (EPL) and obtain their mass density slopes around their respective Einstein radii. Their strikingly similar visual appearance (see Figure 1) and similar stellar mass function to SLACS, SL2S and BELLS, motivates us to compare these observed lens samples to the simulated lens systems. This allows us to compare our findings with observations and draw the following main conclusions:

(1) The stellar mass-size relation and total mass density slope of strong lens galaxies from SLACS, BELLS and SL2S agree best with EAGLE galaxy formation models that have weak or mild AGN activity or in which stellar feedback becomes inefficient at high gas densities (FBconst). In particular, the AGN model with a moderate temperature increment during active periods, $\Delta T=10^{8}\rm K$ (AGNdT8), shows excellent agreement with the observations. Models with no or high-temperature increments agree considerably less well in statistical KS tests. Similarly, the stellar-feedback model with a constant supernova energy injection per unit stellar mass into the surrounding medium (i.e. FBconst) also shows excellent agreement with the observations. Our finding that more efficient feedback yields larger galaxy sizes for a fixed galaxy mass is consistent with previous work by Sales et al. (2010), based on OWLS (Schaye et al., 2010).

(2) Models in which the energy injection per unit stellar mass formed depends either on metalicity or local environment perform less well. Models with a high viscosity also reproduce the total mass density slopes of observed lens galaxies, but perform poorly in reproducing the mass-size relation. The EAGLE Reference model (the benchmark model) also does not perform well, most likely due to a too efficient AGN feedback model. We note that agreement with SL2S is in general worse for all models, which we expect is due to its more heterogeneous selection (as opposed to SLACS and BELLS, they were not selected to be lenses).

(3) Quantitatively, we find that if the simulated lensed images are modeled using an elliptical power law (EPL) profile plus external shear, then the median total mass density slope of galaxies from the AGNdT8 and FBconst models, which have the highest $p$-values in the KS tests, are $t$=2.01 and $t$=2.07, respectively, in good agreement with the observations of SLACS, SL2S and BELLS. Galaxies in the EAGLE Reference model, however, tend to have a steeper median total mass density slope ($t$=2.24) than observed lens galaxies (i.e. $t$ =2.08 for SLACS, $t$=2.11 for BELLS and $t$=2.18 for SL2S). This trend in mass density slope agrees well with the results from other independent analyses (e.g. Remus et al., 2017; Peirani et al., 2018).

(4) We also assess whether in the best model variations that emerged in our analysis (FBconst and AGNdT8) and the benchmark model (Reference), $t$ correlates with stellar mass and found only a mild trend of slopes being shallower than isothermal at higher stellar mass. This is consistent with observations (Auger et al. 2010b; Tortora et al. 2014a) and simulations (Remus et al. 2017; Xu et al. 2017). However, we find no evidence of correlation at any significant level between $R_{\rm Ein}$/$R_{\rm eff}$ ratios and $t$. This is consistent with Koopmans et al. (2006), Auger et al. (2009), Koopmans et al. (2009) and Treu et al. (2009). Thus any selection bias based on mass should therefore not affect the conclusions.

(5) We also find that the mean $R_{\rm Ein}$/$R_{\rm eff}$ ratios in Reference, FBconst and AGNdT8 models are the closest to SLACS. We see a trend in the total mass density slope and $R_{\rm Ein}$/$R_{\rm eff}$ ratio where increasing the feedback efficiency, increases the $R_{\rm eff}$ thereby decreasing the value of $R_{\rm Ein}$/$R_{\rm eff}$ and steepening the total mass density slope ($t$) as in the lens modeling $t$ is calculated around $R_{\rm Ein}$.

Overall we conclude that more efficient feedback in massive galaxies yields steeper total mass density slopes at a radius of $\approx$ 3-10 kpc and that strong lens galaxies appear more consistent with galaxy formation models with somewhat more limited or weaker stellar and/or AGN feedback. Our findings are consistent with the work by Remus et al. (2017) and Peirani et al. (2018) using different simulations. Remus et al. (2017) used the Magneticum Pathfinder (Hirschmann et al., 2014) and two samples, taken from zoom-in re-simulations of Oser simulations (Oser et al., 2010) and Wind simulations (Hirschmann et al., 2013) differing in their baryonic feedback processes. Whereas Peirani et al. (2018) used two varying AGN feedback models of HORIZON-AGN simulations (Peirani et al., 2017).

Duffy et al. (2010) who looked at inner density slopes in the OWLS models found a similar trend, that a weaker feedback is preferred by strong lensing. However, NOAGN feedback does not produce an isothermal profile in our analysis and disagrees with Duffy et al. (2010). These differences may be due to the fact that their mass density slope was obtained at a much higher redshift ($z$=2) and for lower-mass galaxies. Also, they did not create simulated lenses and model them with an EPL model, as done in this work, which might lead to some additional biases. We note that LRGs could have other observational selection biases and might not represent volume limited samples. Our conclusions are not biased by this trend as the evolution of $R_{\rm eff}$ is considerably small (Furlong et al. 2017) in EAGLE.

Our results prefer galaxy-formation models that have been ruled out in Crain et al. (2015) after comparison with non-lensing observations. Furlong et al. (2017) found that the Reference model agrees well with the observed mass-size relation when compared to non-lensing galaxies. This finding is also seen in Duffy et al. (2010), who found that weak feedback is required to match the lensing observations (consistent with our work) but also pointed out that other observables, such as the stellar fractions, rule out those weak feedback models. These seemingly opposing conclusions could be due to either differences in the precise methodologies adopted in the strong-lensing (Duffy et al. 2010, this work) and their non-lensing studies (Crain et al. 2015; Furlong et al. 2017), or additional observational selection biases in the galaxy samples, or even from missing crucial physics. This also might indicate that LRGs that acts as lensing galaxy, might have different formation history than the rest. A complete redshift evolution study of the total mass density slope will be addressed in a forthcoming work.

In this work, we have demonstrated that observables of strong lens galaxies, in particular their total mass density profiles in the inner 3-10 kpc radial range, are very sensitive to variations in the feedback in galaxy formation models. Although strong lensing analysis could have systematical difference from non-lensing analysis in the methods of the modeling. We stress again that SLACS lens galaxies are not different from the parent population of non-lens galaxies from which they were drawn (Treu et al. 2006; Bolton et al. 2008b). In paper III of the SEAGLE series we will explore the systematic errors and compare simulated lenses to non-lensing ETGs from SPIDER survey (La Barbera et al., 2010; Tortora et al., 2012; Tortora et al., 2014a) and, we will show that mass-size relation of EAGLE matches very well with it. Whereas in this paper we have concentrated on the mass-size and mass density slopes, in forthcoming papers we will investigate the inner mass regions in more detail, focusing in particular on the effects of the dark matter distribution and the stellar IMF.

Data Availability

EAGLE data for all the snapshots are publicly available at http://icc.dur.ac.uk/Eagle/database.php. All the simulated lenses and other relevant quantities are available on request to SM/LVEK.

Acknowledgements

We thank the anonymous referee for her/his constructive suggestions to bring the paper to its present form. SM, LVEK, CT and GV are supported through an NWO-VICI grant (project number 639.043.308). SM also acknowledges the funding from COSMICLENS: ERC-2017-ADG, Grant agreement ID: 787886. SM thanks SURFSARA network in Amsterdam. SM thanks Nicolas Tessore and Lyne Van De Vyvere for all the help and support in modeling the lenses. MS is supported by VENI grant 639.041.749. CT also acknowledges funding from the INAF PRIN-SKA 2017 program 1.05.01.88.04. JS is supported by VICI grant 639.043.409. FB thanks the support from the grants ASI n.I/023/12/0 “Attività relative alla fase B2/C per la missione Euclid” and PRIN MIUR 2015 “Cosmology and Fundamental Physics: Illuminating the Dark Universe with Euclid”. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

Appendix A Prior used for density slope

Here in Table 6, we give the prior values used for modelling the simulated strong lenses with LENSED. We use a combination of uniform and Gaussian priors. In Mukherjee et al. (2018), we have explained for the motivation of priors used and also demonstrated the tests that we performed with different prior combinations for the EPL model.

Appendix B Source-size related tests

Here we present some results to demonstrate that the recovered source-sizes do not bias our conclusions. We compare the source-sizes between SIE and EPL and assess the source-size versus slope correlation in those models. These tests are in addition to those carried out in Mukherjee et al. (2018). Readers can consult the Appendix in the latter paper.

In Figure 9, we present the histograms of source-size comparison between SIE and EPL for Reference-100 simulation. We show that the recovered source-sizes agree with the input ones within the error limits for both the models. Also the SIE and EPL modelling provide consistent results. The difference between source-sizes from these two different models is on average 0.008 arcsec, i.e. only 0.4% of the source size. In Figure 10 We also compared the source-sizes between Reference, FBconst and AGNdT8 yielding a mean value of 0.218, 0.217 and 0.213 arcsec, respectively. Thus there is an overall perfect agreement. In Figure 11, furthermore, we compare the source-size of SLACS and the EAGLE Reference-100 model against the density slope. We fitted a liner function to both SLACS and EPL model. We have used the same range of values for SLACS that the EPL models are covering in their sample space. We find that the EPL and SLACS slope and source-size values correlate well with each other. In Figure 12, we also show the relative difference in source-size and EPL density slope from the SIE models. We find a mild anti-correlation having a slope of -0.09 with a Spearman rank of -0.24. No obvious bias is found in our analyses between the EPL and SIE model values and hence we believe the conclusions to be robust.

Finally, even if there were a small bias, such biases would occurs in real lenses as well (see Newton et al. 2011), and hence would broaden both the observed and simulated slope distributions and not impact the inference on the formation scenarios.

Appendix C Comparison with direct fitting

Previously we performed this test between density slopes inferred via convergence fitting, $t_{NM}$ and LENSED, $t_{LENSED}$ in SEAGLE-I and reported (Figure. 8 therein) that there could be a difference of 10% in Einstein radius (see also Küng et al. 2015) and demonstrated that we find a mean ratio of 0.91 for $t_{NM}$/$t_{LENSED}$, with a standard deviation of 0.17 (Figure 9 therein). Eventhough, the lens modeling fits the density profile (more precisely that of the potential) near the lensed images, whereas the direct fit is mostly fitting the higher density regions inside the mask, we do not find any biased results from these two different methods. In Figure 13 we have shown the mean density slope comparison between Reference, FBconst and AGNdT8 models.

Appendix D In-exact modeling: Source parameter variation

We have carried out several tests with a representative combination of source structures. Table 7 summarizes the results. We use three types of complex sources: (a) two Sérsic profiles with Sérsic indices 1 and 3, (b) one Sérsic profile with index n=3, and (c) two Sérsic profiles with indices 1 and 3 plus random noise/perturbation. We fit the source with one Sersic profile. We find that the density slope, Einstein radius and ellipticity (key parameters in this work) are obtained fairly consistently. Thus we believe the conclusions in this work to be robust against modest differences between the source model and the true source structure. We do acknowledge there could be slight differences on a case by case basis, but in a statistical sense the inferences from this analysis are not affected by the source structure, and our modelling behaves similarly between real and simulated lenses. We also tested if the impact of source size have any selective effects on model variations. Table 8 summarizes the results. We chose a random lensing galaxy from each model variation in the stellar mass range $M_{*}\sim 10^{11.0-11.1}M_{\odot}$. We use different source sizes other than the one used in this analysis. However, we find consistent result for mass density slope with variations that are appreciably lower than the average r.m.s. error that ranges from 0.1–0.2.