A little FABLE: exploring AGN feedback in dwarf galaxies with cosmological simulations

Contrary to the standard lore, there is mounting observational evidence that feedback from active galactic nuclei (AGN) may also play a role at the low-mass end of the galaxy population. We investigate this using the cosmological simulation suite FABLE, with a particular focus on the dwarf regime ($M_\mathrm{stellar}<10^{9.5} \ \mathrm{M_{\odot}}$). We find that overmassive black holes (BHs), with respect to the mean scaling relations with their host galaxies, drive hotter and faster outflows and lead to significantly reduced gas mass fractions. They are also more likely to display a kinematically misaligned ionized gas component in our mock MaNGA velocity maps, although we caution that cosmic inflows and mergers contribute to misalignments as well. While in the local Universe the majority of AGN in dwarfs are much dimmer than the stellar component, for $z \geq 2$ there is a significant population that outshines their hosts. These high-redshift overmassive BHs contribute to the quenching of dwarfs, whereas at late cosmic times supernova (SN) feedback is more efficient. While our results are overall in good agreement with X-ray observations of AGN in dwarfs, the lack of high-luminosity X-ray AGN in FABLE at low redshifts highlights an interesting possibility that SN feedback could be too strong in FABLE's dwarfs, curtailing AGN growth and feedback. We predict that future observations may uncover many more AGN in dwarfs with lower luminosities and at higher redshifts.


INTRODUCTION
In the ΛCDM model of structure formation, primordial density fluctuations grow and collapse into dark matter (DM) haloes. The gas collapses within these haloes, cools and condenses to form stars. The dark and baryonic components then build up together by accreting and merging with other structures. However, when comparing the DM halo mass function to the galaxy stellar mass function (GSMF), both the low-mass end and the high-mass end of the GSMF are found to be significantly suppressed.
Feedback processes are usually invoked to explain this mismatch and to regulate the overall efficiency of star formation. In the common theoretical model, the low-mass end of the GSMF is suppressed by reionization (e.g. Efstathiou 1992;Okamoto et al. 2008;Fitts et al. 2017) and SN feedback (e.g. Larson 1974;Dekel & Silk 1986;Mori et al. 2002), whilst AGN feedback regulates the high-mass end (e.g. Binney & Tabor 1995;Di Matteo et al. 2005;Bower et al. 2006;Croton 2006;Sijacki et al. 2007;Cattaneo et al. 2009).
At the high-mass end of the galaxy population, simulations have shown that AGN feedback can bring stellar properties into E-mail: skoudmani@ast.cam.ac.uk agreement with observations (Puchwein & Springel 2012;Dubois et al. 2014;Vogelsberger et al. 2014a;Khandai et al. 2015;Schaye et al. 2015;Weinberger et al. 2017;Henden et al. 2018). However, whether reionization and SN feedback are sufficient to regulate star formation at the low-mass end is still unclear. Whilst some find that these two processes can suppress star formation in low-mass haloes (e.g. Governato et al. 2007;Sawala et al. 2016), other groups find that additional physical processes (e.g. photoionization, radiation pressure, or stellar winds) are required for SNe to couple to the interstellar medium (ISM) effectively and drive outflows (see e.g. Hopkins et al. 2014;Kimm et al. 2015;Emerick et al. 2018;Smith et al. 2019). Recently, it has been speculated that AGN may also contribute to regulating star formation at the low-mass end (see e.g. Silk 2017). This idea has been motivated by several observational studies which have identified a population of dwarf galaxies hosting AGN.
Some of the first systematic searches for AGN in dwarf galaxies were carried out by searching SDSS dwarf galaxies for optical AGN signatures (see e.g. Greene & Ho 2004Reines et al. 2013;Chilingarian et al. 2018). These studies use the BPT (Baldwin et al. 1981) diagram to confirm the AGN nature of these sources and broad optical emission lines to estimate the BH mass, where available, measuring BH masses down to M BH ∼ 3 × 10 4 M .
Note that Cann et al. (2019) find that standard optical spectral classification schemes used to identify higher mass BHs do not apply when the BH mass falls below M BH ∼ 10 4 M , so BPT diagnostics may become less useful for exploring AGN activity in the lowmass regime. Moreover, low-mass AGN can easily be outshone by star formation in the optical band. Therefore only bright AGN are recovered resulting in low occupation fractions, e.g. Reines et al. (2013) obtain an AGN occupation fraction of ∼ 0.5 per cent for their optically-selected sample of AGN in dwarf galaxies.
These two issues underline the importance of multi-wavelength studies. For example, X-ray observations have the advantage that they are less biased towards type 1 AGN (especially at higher redshifts), and that they suffer less from contamination from star formation. Follow-up X-ray observations can then be used to confirm optically-selected AGN (Baldassare et al. , 2017Chilingarian et al. 2018). Furthermore, wide-field X-ray surveys have identified large samples of AGN in dwarf galaxies via their X-ray emission (e.g. Schramm & Silverman 2013;Lemons et al. 2015;Miller et al. 2015;Mezcua et al. 2016Mezcua et al. , 2018Pardo et al. 2016;Aird et al. 2018;Birchall et al. 2020). These X-ray surveys also include AGN in dwarf galaxies that would have been missed by the BPT diagnostics due to the optical AGN contribution being less than that from galaxy emission (see e.g. Birchall et al. 2020) resulting in higher AGN occupation fractions than for the optically-selected samples. For example, Pardo et al. (2016) find an overall AGN occupation fraction of 0.6−3 per cent, whilst the occupation for high X-ray luminosities (log(L X [erg s −1 ]) ≥ 41.5) is estimated to be ∼ 0.2 − 0.4 per cent (Aird et al. 2018;Mezcua et al. 2018;Birchall et al. 2020), more akin to the optically-selected samples which are likely dominated by bright sources.
Furthermore, infrared (IR) surveys have been able to identify 'optically-hidden' AGN in dwarf galaxies (e.g. Satyapal et al. 2014;Sartori et al. 2015;Marleau et al. 2017;Kaviraj et al. 2019), resulting in significantly higher AGN fractions than from optical samples, e.g. Kaviraj et al. (2019) find AGN occupation fractions for dwarf galaxies in the range ∼ 10 − 30 per cent. Though see Lupi et al. (2020) who caution that resolution effects and source overlapping may contaminate IR-selected AGN samples in the low-mass regime.
Finally, there are several future facilities that are expected to find large samples of BHs in dwarf galaxies. Long-term optical variability searches are a promising tool for finding AGN missed by other selection techniques, especially in the era of LSST (e.g. Baldassare et al. 2020). The gravitational wave observatory LISA will be aimed at detecting gravitational waves from the coalescence of BHs with masses M BH ∼ 10 4 −10 7 M and is therefore optimally placed to study the BH population in the low-mass regime. Recent simulation studies have highlighted the multimessenger signatures of BHs in dwarfs (Bellovary et al. 2019;Volonteri et al. 2020). JWST will also provide important insights into AGN in dwarf galaxies and alleviate the above-mentioned resolution issues of current IR studies (e.g. Cann et al. 2018). The combination of JWST with future X-ray missions (Athena, Lynx) will be a powerful tool for mapping the AGN population in dwarf galaxies at high redshifts (e.g. Civano et al. 2019;Haiman et al. 2019). The ngVLA will map the weakly-accreting AGN in dwarf galaxies and provide insights into the jet-driven feedback mechanisms in the low-mass regime as well as BH seeding mechanisms (e.g. Nyland & Alatalo 2018;Plotkin & Reines 2018).
Beyond constraining the AGN demographics at the low-mass end, recent observations have started studying whether AGN in dwarf galaxies could have a significant impact on their hosts -with a particular focus on quenching and outflows. Penny et al. (2018) study quiescent dwarf galaxies in the MaNGA survey and find that six out of 69 quiescent MaNGA dwarfs have optical AGN signatures. Interestingly, five out of six of these dwarf have an ionized gas component that is kinematically offset from the stellar component which could be a sign of AGN-driven outflows. Manzano-King et al. (2019) find direct evidence of spatially extended high-velocity ionized gas outflows in nine dwarf galaxies hosting AGN, for six of these the emission-line ratios of the outflow itself are consistent with AGN ionization. Dickey et al. (2019) study quiescent dwarfs in isolated environments and find that 16 out of 20 quiescent isolated dwarfs have AGN-like line ratios.
These recent observations of AGN activity in dwarf galaxies have motivated theorists to investigate this largely unexplored regime -with mixed results. On the one hand, analytical models of AGN feedback in dwarf galaxies look promising, suggesting that AGN could provide an alternative and more successful source of negative feedback than SNe (Silk 2017;Dashyan et al. 2018). On the other hand, in numerical simulations, strong SN feedback can hinder the growth of BHs in low-mass galaxies rendering AGN feedback ineffective (see e.g. Dubois et al. 2015;Anglés-Alcázar et al. 2017b;Habouzit et al. 2017;Trebitsch et al. 2018). However, if BHs in simulations are able to grow efficiently, they can have a significant impact on their hosts (Barai & De Gouveia Dal Pino 2019;Sharma et al. 2020). We have systematically tested the effect of AGN activity in a series of isolated simulations of dwarf galaxies (Koudmani et al. 2019). We found that the outflows reach significantly higher velocities and temperatures with added AGN feedback. However, AGN activity did not have a significant effect on the global star formation rate (SFR) in our isolated set-up.
The considerations above suggest an interesting possibility that AGN feedback may play an indirect role in regulating star formation at the low-mass end. For example, the AGN-boosted outflows might be able to hinder cosmic gas inflows. Cosmological simulations are needed to investigate these possibilities.
In this paper, we use the (Feedback Acting on Baryons in Large-scale Environments) simulation suite (Henden et al. 2018(Henden et al. , 2019a to explore the role of AGN feedback in the low-mass regime within a realistic cosmological environment. Keeping in mind that the parameters have been tuned so that SN feedback is strong to regulate the low-mass end of the GSMF, we investigate whether low-mass galaxies with efficient BH growth have significantly different physical properties than the ones without AGN activity -with a particular focus on outflow properties and SFRs. Though we note that any effect we find from AGN feedback in the low-mass regime is likely to be a lower limit. Furthermore, we aim to compare the galaxies to observations of dwarf galaxies with AGN and to make predictions for future observations. The outline of this paper is as follows. We describe the characteristics of the simulation suite and the sub-grid models in Sections 2.1 and 2.2, respectively. The methods for identifying haloes and galaxies in and the definition of the low-mass galaxy sample are given in Section 2.3 and we describe how we calculate the outflow properties of these galaxies in Section 2.4. Our methods for creating mock MaNGA observations and X-ray luminosities are outlined in Sections 2.5 and 2.6. In Section 3, we present our results on the AGN population (Section 3.1), the scal-ing relations (Section 3.2), the properties of low-mass galaxies with overmassive BHs (Section 3.3), mock MaNGA observations (Section 3.4), and mock X-ray observations (Section 3.5). We discuss our results in Section 4 and summarise our main findings in Section 5.

Basic simulation properties
We use the cosmological simulation suite for our analysis of AGN feedback in low-mass galaxies. The simulation suite is extensively described in Henden et al. (2018), and only a brief summary is given below.
The simulations were carried out with the massively parallel code (Springel 2010), where fluid dynamics is discretized on a moving mesh using a quasi-Lagrangian finite volume technique. The unstructured mesh is based on the Voronoi tessellation of a set of discrete points that cover the whole computational domain. These mesh-generating points are allowed to move with the local flow velocity, with minor corrections to avoid excessive distortion of the gas cells. Gravitational interactions are modelled using the TreePM approach, with stars and DM modelled as collisionless particles.
The suite consists of a large-volume cosmological box simulation as well as a series of zoom-in simulations of individual groups and clusters. As we are interested in low-mass galaxies in representative regions of the Universe, we focus only on the results from the cosmological box for this study.
The comoving 40 h −1 Mpc box was evolved using initial conditions for a uniformly-sampled cosmological volume based on a Planck cosmology (see Planck Collaboration XIII 2016). The volume contains 512 3 DM particles with masses m DM = 3.4 × 10 7 h −1 M and initially 512 3 gas resolution elements with target mass m gas = 6.4 × 10 6 h −1 M . The gravitational softening length is set to 2.393 h −1 kpc in physical coordinates below z = 5 and fixed in comoving coordinates at higher redshift.
The sub-grid models are largely based on the Illustris galaxy formation model (Vogelsberger et al. 2013(Vogelsberger et al. , 2014aGenel et al. 2014;Torrey et al. 2014;Sijacki et al. 2015). Whilst the models for star formation (Springel & Hernquist 2003), radiative cooling (Katz et al. 1996;Wiersma et al. 2009a), and chemical enrichment (Wiersma et al. 2009b) are unchanged from Illustris, the stellar feedback (Vogelsberger et al. 2013) and AGN feedback Springel et al. 2005;Sijacki et al. 2007Sijacki et al. , 2015 models have been updated (see Section 2.2). Both the stellar and the AGN feedback models were calibrated to reproduce observations of the local GSMF and the gas mass fractions of observed massive haloes.

Stellar feedback
Stellar feedback is included as a sub-grid model for galactic winds and outflows (see Vogelsberger et al. 2013, for details). In the Illustris simulation, these winds are purely kinetic. In , this is modified so that one third of the wind energy is imparted as thermal energy, which slows down the cooling of the gas thereby increasing the overall effectiveness of the stellar feedback (see Marinacci et al. 2014;Henden et al. 2018).

BH seeding and growth
BHs are modelled as collisionless sink particles. BH seeds of mass M BH,seed = 10 5 h −1 M are placed into every DM halo above 5 × 10 10 h −1 M . Subsequently, these seeds can grow via mergers with other BHs and via gas accretion (see Vogelsberger et al. 2013, for details). In brief, BH particles accrete at the Eddington-limited Bondi-Hoyle-Lyttleton-like rate boosted by a factor of α = 100 (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944;Springel et al. 2005). The radiative efficiency is set to be a constant r = 0.1. A fraction of (1 − r ) of the accreted mass is added to the BH mass, whilst the rest is available as feedback energy.

AGN feedback
In analogy to Illustris, the AGN feedback operates in either one of two states depending on the Eddington fraction f Edd : the high accretion rate quasar mode Springel et al. 2005) or the low accretion rate radio mode (Sijacki et al. 2007). The BH particles switch to the quasar mode once they exceed a critical Eddington fraction f Edd,QM . In the model, this threshold is set to f Edd,QM = 0.01. This AGN feedback model has been modified in by introducing a duty cycle to the quasar mode. The feedback energy is accumulated over δt = 25 Myr, before being released in a single event, to reduce overcooling (see Booth & Schaye 2009;Henden et al. 2018).
The remainder of the feedback energy that is not coupled to the surrounding gas by the quasar or radio mode goes into radiative electromagnetic feedback, which is approximated as an additional radiation field around the BH superposed with the redshift-dependent ultraviolet background (Vogelsberger et al. 2013).

Galaxy identification
Haloes and subhaloes are identified via friends-of-friends (FoF) and algorithms (Davis et al. 1985;Springel et al. 2001;Dolag et al. 2009). For the FoF search, we choose a linking length of 0.2 multiplied by the mean inter-particle separation. We then use to identify gravitationally self-bound subhaloes within each FoF halo. The 'central subhalo' is the subhalo at the minimum gravitational potential of the FoF halo. All other subhaloes in the FoF halo are categorised as satellite subhaloes.
We consider each subhalo with at least 100 star particles to be a resolved galaxy, yielding a minimum stellar mass of approximately 10 9 M . For this investigation, we only consider central subhaloes, as our focus is on intrinsic feedback processes rather than environmental effects.
At each redshift, our low-mass galaxy sample consists of all central subhaloes that host a BH and have a stellar mass in the range 9.0 ≤ log(M stellar [M ]) < 10.5 at the given redshift. Stellar mass here is defined as the mass of all star particles within twice the stellar half mass radius, r eff . The minimum mass is set by the above resolution requirements. In addition, we restrict our sample to galaxies with masses below the knee of the GSMF (M * = 10 10.5 M ), where AGN feedback had previously been assumed to be ineffective.
To investigate the impact of AGN feedback as a function of stellar mass, we divide up the low-mass galaxy sample into three mass bins of equal widths in log space: dwarf galaxies (9.0 ≤ log(M stellar [M ]) < 9.5), massive dwarf galaxies (9.5 ≤ log(M stellar [M ]) < 10.0), and M * galaxies (10.0 ≤ log(M stellar [M ]) < 10.5).

Outflow properties
We would like to assess whether AGN feedback can have an impact on the overall outflow properties of low-mass galaxies. To this end, we extract the outflow properties from the simulation data for the whole gas (i.e. not separating out the different gas phases). Though note that does not include cooling below 10 4 K, so these outflow properties mainly correspond to the warm and to the hot gas phase.
We obtain the outflow properties for each galaxy as follows: Firstly, we centre the gas coordinates on the potential minimum of the galaxy. Then we subtract the peculiar velocity (here defined as the mass-weighted mean velocity of all gas cells within twice the virial radius 1 , 2R vir ) from the gas velocities. We then calculate the outflow properties at two different target radii r t = 0.2R vir , R vir within shells of widths ∆r = 10, 20 kpc, respectively 2 . We select all gas cells within the r t ± ∆r 2 shell and calculate the radial velocities v r for all of these cells. We then mark gas cells with v r < 0 as inflowing and gas cells with v r > 0 as outflowing. Note that we only calculate the outflow properties for galaxies that have at least ten gas cells (at least five inflowing and at least five outflowing) within the shell considered.
We then use the inflowing and outflowing gas cells to calculate the mass flow rates through the r t ± ∆r 2 shell as follows: Here m cell is the gas cell mass, M out is the outflow rate, M in is the inflow rate, and M tot is the total mass flow rate.
We also calculate the outflow velocities and temperatures. We define the outflow velocity, v out , as the mass-weighted mean velocity of all gas cells within the shell centred at r t with v r > 0. Similarly, we define the outflow temperature, T out , as the mass-weighted mean temperature of all cells within the shell centred at r t with v r > 0.
Using the above quantities, we then also calculate the momentum outflow rate, P out , and the kinetic energy outflow rate, E kin,out , as: Finally, we also calculate the mass loading factor, β, as where SFR(r ≤ r t ) is the summed up SFR of all gas cells within r t .
1 Throughout this paper, the virial quantities of a system are defined by a mean density of 200× the critical density of the Universe. 2 Note that in both cases we tested various shell widths and found that the results are well converged and do not depend on the exact choice of ∆r, as long as the width is large enough to have a sufficient number of gas cells within the shell.

Mock MaNGA maps
The MaNGA survey found kinematically misaligned gas in quenched dwarf galaxies with AGN signatures hinting at AGNdriven outflows (see Penny et al. 2018). In Koudmani et al. (2019), we created mock MaNGA maps of isolated simulations of dwarf galaxies to explore these kinematic offsets. Here we adapt this procedure to cosmological simulations, as described below.

Sample selection
To accurately determine the kinematic position angle (PA), the galaxies need to have clear rotational features. To automatically select galaxies that are rotationally supported, we only keep galaxies with V max /σ > 1.1, where V max is the maximum value of the respective subhalo's spherically averaged rotation curve and σ is the 3D velocity dispersion of the same subhalo. Visual inspection of our sample galaxies confirms that this is an efficient criterion for ensuring sufficient rotational support.
In addition, we require that there should be at least 50 star particles and 50 ionized gas cells within the 1.5r eff 2D aperture (corresponding to the spatial coverage of the primary MaNGA sample). Here r eff is the stellar half-mass radius. This ensures that the kinematics are well resolved by the simulation.

Line-of-sight velocity maps
We match the resolution of our mock MaNGA maps to the pixel size of the MaNGA final reduced data cubes (0.5 arcsec). We analyse three different redshift (z = 0, 0.2, 0.4). Note that for the local z = 0 sample, we calculate the spatial resolution assuming that these galaxies are observed at the mean redshift of the primary MaNGA sample,z = 0.03. We do not consider higher redshifts than z = 0.4, as the resolution of the MaNGA survey would prohibit us from confidently identifying kinematic PAs.
We create separate line-of-sight (l.o.s) velocity maps for the star particles and the ionized gas cells. All star particles or ionized gas cells that belong to the galaxy's FoF group are included along the projection direction. We keep all gas cells with neutral hydrogen abundance, fraction of the hydrogen cell mass (or density) in neutral hydrogen, n H,neutral < 10 −2 to isolate the ionized gas 3 .
For each subhalo, we centre the velocity map on the minimum gravitational potential and subtract the systemic velocity, here taken as the mass-weighted mean velocity within 1.5r eff .
To separate out the effect of the inclination angle on kinematic misalignments, we rotate all galaxies into the same orientation for our analysis. To this end, we calculate the stellar total angular momentum L * within 1.5r eff and rotate each galaxy (both stars and gas) so that L * is aligned with the z-axis. We then project the stellar/ gas kinematics onto the y-z plane so that all galaxies are viewed edgeon. As an optional additional step, the disc can then be inclined by rotating about the y-axis to test the effect of the inclination angle.
We then project the ionized gas/ stellar velocity distribution onto a grid corresponding to the spatial extent and resolution of the MaNGA data. We define the stellar smoothing length as the radius enclosing the 64 nearest neighbours and the gas smoothing length as 2.5 times the cell radius. We set a maximum smoothing length as 50 times the pixel size, though in practice that limit is rarely reached.
We weight the velocity contributions to each pixel by the cubic spline smoothing kernel. In addition, we impose a luminosity-like weighting: the stars by their mass and the gas cells by the square of their density. Note that when creating these smoothed maps we consider all resolution elements within a 2.5r eff aperture as some might contribute to pixels within the 1.5r eff aperture even though their centres lie outside that aperture.
Finally, we convolve the resulting velocity map with a Gaussian filter of FWHM = 2.5 arcsec to model the effect of the MaNGA PSF.

Kinematic position angles
The kinematic PA is related to the projected angular momentum vector and traces the position of the maximum change in velocity on the map (see e.g. Krajnović et al. 2006). Following Penny et al. (2018), we determine the kinematic PAs of both ionized gas and stars using the fit_kinematic_pa routine (see Krajnović et al. 2006, for a detailed description). Note that this routine only returns PAs in the range 0°≤ PA < 180°, which does not take the sense of the rotation into account. We identify the redshifted part of the velocity map to convert the PAs to the 0°≤ PA < 360°range.
The kinematic offset between gas and stars is then defined as ∆PA = |PA stars − PA gas |. Gas and stars are considered to be in dynamical equilibrium if ∆PA = 0°or 180°. As in the observational study by Penny et al. (2018), we then define the gas as kinematically misaligned if 30°≤ ∆PA < 150°.

Mock X-ray luminosities
X-ray studies by Chandra and XMM-Newton have resulted in large samples of AGN candidates in low-mass galaxies. As X-ray surveys push to higher redshifts, it has become possible to start studying the redshift evolution of AGN activity in these galaxies. To compare the AGN population in to these observational findings, we calculate the X-ray properties of the galaxies. We define the bolometric BH luminosity as L BH = r M BH c 2 , where r is the radiative efficiency, M BH is the BH accretion rate and c is the speed of light (see Section 2.2.3). To obtain the BH luminosities in the soft (0.5 -2 keV) and hard (2 -10 keV) X-ray band, we use the bolometric corrections from Shen et al. (2020). To assess the detectability of these low-luminosity AGN, we also estimate the X-ray luminosities of X-ray binaries (XRBs) and the hot gas in the galaxy.
To obtain an estimate for the contribution from XRBs, we use the redshift-dependent relation from Lehmer et al. (2016) which relates the total X-ray luminosity L X,XRB to the SFR and stellar mass (M stellar ): The parameter values in the soft X-ray band are log(α 0 ) = 29.04, log(β 0 ) = 39.38, γ = 3.78, and δ = 0.99. In the hard X-ray band, the parameter values are given by log(α 0 ) = 29.37, log(β 0 ) = 39.28, γ = 2.03, and δ = 1.31. To estimate the luminosities, we evaluate both the SFR and M stellar within 2r eff .
We calculate the hot gas contribution using the relationship between the X-ray luminosity of the diffuse ISM and the SFR from Mineo et al. (2012): Again, we calculate the SFR within 2r eff here. Note that the above relation only applies to soft X-ray luminosities. To obtain the contribution in the hard X-ray band, we assume a photon index Γ = 3, which is a good representation of a thermal model with temperature 0.7 -1 keV (see Mezcua et al. 2018).

The AGN population
We start our analysis by inspecting the large-scale distribution of AGN in low-mass galaxies in our simulation. Figure 1 shows projections of the total gas surface density for the wholesimulation box at z = 0 (left panel) and z = 1 (right panel). The markers show the locations of active low-mass galaxies (with Eddington fractions f Edd > 10 −3 ). The galaxies' stellar masses are indicated by the marker style with the sample split into dwarfs with 9.0 ≤ log(M stellar [M ]) < 9.5, massive dwarfs with 9.5 ≤ log(M stellar [M ]) < 10.0, and M * galaxies with 10.0 ≤ log(M stellar [M ]) < 10.5. Focusing on the spatial distribution, there are low-mass galaxies with AGN in all types of cosmic environments: voids, filaments, and knots. There are significantly more active low-mass galaxies at redshift z = 1 than at z = 0.
The marker colour shows the galactic outflow rates at the virial radius, R vir for the total gas (see Section 2.4). Across all mass bins (and at both redshifts), we find active galaxies with large-scale outflows at the virial radius ( M out 10 M yr −1 ). There are higher outflow rates and velocities at redshift z = 1 compared to z = 0. Unsurprisingly, more massive galaxies tend to be associated with stronger outflows, given the larger SFRs and BH masses in those systems. Though the more massive galaxies also have larger escape velocities than the low-mass objects.
For each of the three mass bins at both redshifts, we select a representative example galaxy with a high outflow rate and show 2R vir × 2R vir × 2R vir slices of the outflows in the zoomed-in insets.
The inset colour-coding shows projections of the radial velocity of the total gas with the streamlines overlaid on the projection. Note that the galaxies shown in these projections have been rotated so that they are viewed edge-on. The more massive galaxies have collimated outflows, whilst the outflows for the example dwarf galaxies are more isotropic. For all of the example galaxies shown here, the outflows are able to propagate to the virial radius.
The outflow velocities of the total gas at R vir , as defined in Section 2.4, range from ∼ 50 to 150 km s −1 . In most cases the outflows decelerate between 0.2R vir and R vir . To estimate the escape velocity for these haloes we use the virial velocity, V vir , for simplicity, which ranges from ∼ 80 to 210 km s −1 . For most of these example galaxies, the mean outflow velocity at R vir does not exceed the virial velocity, indicating galactic fountain nature of the outflows. The exception to this is the dwarf example galaxy at z = 0. For this galaxy, the outflows accelerate from 42 km s −1 at 0.2R vir to 94 km s −1 at R vir , exceeding V vir = 81 km s −1 . Note that this galaxy has an overmassive BH relative to its stellar mass with M BH = 4 × 10 5 M and M stellar = 10 9 M (see scaling relations in Figure 3) and is in a crowded region with strong cosmic inflows.
Even though the outflow velocities of the total gas do not exceed V vir in the other cases shown, the fast phase of the outflow is able to reach velocities higher than V vir in all cases. The maximum ). The marker colour shows the mass outflow rate of the total gas at the virial radius, R vir . For each of the three mass bins at both redshifts, we select an example galaxy with a high outflow rate and show the projected outflow kinematics of the total gas in the zoomed-in insets (projection dimensions 2R vir × 2R vir × 2R vir ). The inset colour-coding indicates the radial velocity and the velocity streamlines are shown in black. Note that these insets have been rotated so that the galaxies are viewed edge-on.
velocities of the outflowing gas range from ∼ 150 to 500 km s −1 . We also checked the temperature distributions of these outflows and find that the outflows at R vir are multiphase with temperatures ranging from 8 × 10 3 K to 3 × 10 6 K. The fast phase of the outflow corresponds to the hot phase, with temperatures between ∼ 10 5 K and 10 6 K.
Having established that a significant number of low-mass galaxies with AGN in have appreciable outflows, we move on to investigate the radiative properties of this population. Figure 2 shows the distribution of the low-mass AGN population in L BH − f Edd space, where L BH is the bolometric BH luminosity. We show this distribution at a range of redshifts between z = 0 and z = 4, focusing on the high Eddington fraction end of this distribution ( f Edd ≥ 10 −4 ). The threshold Eddington fraction for transition to the quasar mode ( f Edd,QM = 10 −2 ) is shown as a blue dashed line. Furthermore, the fraction of BHs in the quasar mode (out of all BHs hosted by low-mass galaxies) is given in the upper left corner for each redshift.
The colour-coding in the upper panel indicates the average stellar mass in each pixel. Dwarf galaxies (with log(M stellar [M ]) < 9.5) dominate the galaxy sample at high redshifts and the high-z distribution is quite narrow. Moving to lower redshifts, the distribution broadens as we get a wider range of stellar masses.
There is a clear stellar mass gradient where at a fixed BH luminosity, BHs with higher Eddington fractions tend to be hosted by less massive galaxies. There are very few highly-accreting BHs at low redshifts (z < 1) rendering these AGN difficult to detect. Though note that may underproduce high-luminosity sources at low redshifts (see Section 3.5).
At z = 0, only one per cent of BHs in low-mass galaxies are in the quasar mode. This increases to 10 per cent at z = 1, and goes up to 51 per cent at z = 4. Future facilities (such as Athena, Lynx, JWST or the Roman Space Telescope) will be able to probe this early phase of significant activity, as predicted by our simulations. Note that we have compared the evolution of the number density of AGN for different luminosity bins in with Hopkins et al. (2007). We found that this is in good agreement and reproduces the cosmic downsizing effect, so the overall evolution of AGN luminosities in appears to be realistic.
At low redshifts, observational searches according to our model are likely only uncovering the tip of the iceberg. This issue is explored in more detail in the lower panel where the colour coding shows the average ratio between the (bolometric) BH luminosity L BH , and the central (bolometric) stellar luminosity L stellar . Here we take the central stellar luminosity to be the luminosity generated by all of the star particles within the stellar half mass radius. We estimate the luminosity of each star particle using the up-to-date Bruzual & Charlot (2003) stellar population synthesis models for a Chabrier (2003) IMF, and calculate the bolometric luminosity as a function of stellar age and metallicity.
The majority of AGN in low-mass galaxies are outshone by the stars in their host galaxies -even for highly accreting BHs. Therefore resolved studies as well as multi-wavelength surveys, such as probing the IR and the X-ray regime, are crucial tools for mapping this elusive population (see Section 1). There are very few bright AGN in low-mass galaxies at low redshifts, rendering these AGN hard to detect. The high-redshift population is much more promising with ∼ 25 − 50 per cent in the quasar mode. The bottom panel illustrates how in the majority of cases, the stars in the galaxy outshine the AGN, demonstrating the importance of multi-wavelength studies. Again this problem is alleviated for high Eddington fraction accretors at higher redshifts.

Scaling relations
We move on to investigate the relationship between the BHs and the stellar and gaseous components of their host galaxies. Figure 3 shows the scaling relations of BH mass, M BH , against stellar mass, M stellar , (left panels) and BH mass against gas mass, M gas , (right panels) at redshifts z = 0 and z = 1.
Both M stellar and M gas are calculated within twice the stellar half mass radius, so these are scaling relations for a proxy for the total galaxy mass, rather than comparing just to the bulge component. We use the total galaxy mass since bulge decompositions are rather difficult for observations at the low-mass end (e.g. MacArthur et al. 2003;Reines & Volonteri 2015;Greene et al. 2019) due to the high resolution requirements (although see Schutte et al. 2019).
To investigate whether there is a link between the scatter in the scaling relations and feedback activity in , we define 2D bins (bin dimensions: 0.1 dex × 0.2 dex) with at least ten objects and colour-code them by their mean outflow rate at R vir (using the outflow rate definition from Section 2.4 for the total gas). For more sparsely-populated bins, we plot the individual objects colour-coded by their respective outflow rates. This allows us to show the mean outflow properties of the low-mass galaxy population, whilst also indicating the behaviour of outliers.
The binned mean scaling relations (bin width 0.1 dex, at least ten objects per bin) for the galaxies are shown by the solid black lines with the filled black circles indicating the bin midpoint values. For comparison, we plot some of the observed M BH −M stellar relations. Here, we mainly focus on relations derived for late-type galaxies (marked by 'LT'). These are most applicable to our lowmass galaxy sample as late-type (star-forming) galaxies dominate the local GSMF at the low-mass end (e.g. Muzzin et al. 2013;Vogelsberger et al. 2014a). This distinction is important as scaling relations obtained from late-type galaxies have a lower normaliza-tion than scaling relations derived from early-type galaxies (Greene et al. 2010Reines & Volonteri 2015;Läsker et al. 2016;Davis et al. 2018).
The relation shown from Reines & Volonteri (2015) is derived from a sample of local AGN (a significant number of which are spirals or discs), mostly containing broad-line AGN but also a subsample of dwarfs and reverberation-mapped AGN. This relation also has a significantly lower normalization than for early-type galaxies (see Reines & Volonteri 2015, for a detailed discussion). Moreover, we include the scaling relation derived by Davis et al. (2018) from a sample of 40 spiral galaxies with directly measured BH masses 4 . Finally, we show the relations obtained by Greene et al. (2019). These are based on the dynamical sample from Kormendy & Ho (2013) as well as additional dynamical measurements published since then Saglia et al. 2016;Krajnović et al. 2018;Thater et al. 2019), particularly at low masses (Den Brok et al. 2015;Nguyen et al. 2018Nguyen et al. , 2019, as well as upper limits (Barth et al. 2009;Neumayer & Walcher 2012;De Lorenzi et al. 2013). In addition to the relation derived from the late-type Greene et al. (2019) sample, we also show the relation derived from the whole Greene et al. (2019) sample, including both early-type and late-type galaxies. This relation has a greater than 0.7 dex higher normalization, demonstrating the significant effect of including early-type galaxies in the analysis.
It should be noted that there is a scarcity of dynamical BH mass measurements below M stellar = 3 × 10 10 M , so care should be taken when comparing any of the above observed relations with simulated low-mass galaxies. Broad-line AGN can be used as an , though the BHs are undermassive compared to most of the observed scaling relations. We also show the binned distribution of the whole low-mass galaxy sample with colour-coded histograms. 2D bins with at least ten objects are colour-coded by the mean outflow rate of the total gas at R vir , and otherwise we plot the individual objects colour-coded by their respective outflow rates to indicate the behaviour of outliers. Higher outflow rates are obtained for BHs with larger positive offsets from the M BH − M gas relation. At z = 1, the M BH /M stellar and M BH /M gas ratios increase and higher outflow rates occur. alternative means to establish a scaling relation for observed lowmass galaxies, however this comes with its own difficulties, as the virial factor f vir still carries high uncertainties.
An additional complication is introduced by the different stellar mass-to-light ratios used in observational studies (see discussion in Reines & Volonteri 2015;Davis et al. 2018Davis et al. , 2019. Note that this has no effect on the slope of the scaling relation, just on the normalization. Fully adjusting for this effect would be beyond the scope of this paper, however, we reduce the stellar masses from Greene et al. (2019), which were obtained from the Bell et al. (2003) fitting functions, by -0.093 dex (Gallazzi et al. 2008;Zibetti et al. 2009), so that they are appropriate for a Chabrier IMF, as adopted in . For the Reines & Volonteri (2015) and Davis et al. (2018) relations, we make no adjustments, as these stellar masses were calculated assuming a Chabrier IMF and we do not attempt to harmonise the stellar mass measurements beyond a common IMF.
Keeping the above caveats in mind, we note that the local BH population in is undermassive compared to most of the observed scaling relations, with the exception of the Davis et al. (2018) relation. This relation has a significantly lower normalization than the other observed relations and brackets the lower end of the distribution. The closest agreement is reached with the Reines & Volonteri (2015) relation, though the mean relation lies consistently below this observed relation with discrepancies from 0.1 − 0.4 dex. The late-type relation from Greene et al. (2019) brackets the upper end of the distribution, whilst the overall relation from Greene et al. (2019), including both early-type and late-type galaxies, lies significantly above the mean relation derived from and has hardly any overlap with the simulation data. We also caution that at the low-mass end of the BH mass distribution, where BH masses have not grown considerably from the seed mass of M BH,seed = 10 5 h −1 M , alternative seeding models could affect the simulated scaling relations.
Despite the uncertainties regarding the normalization of the observed relations, Figure 3 suggests that may underproduce the masses of AGN in low-mass galaxies at z = 0 and that any effect we find from AGN on low-mass galaxies in might constitute a lower limit of the possible effect of AGN feedback in that mass range. This is also supported by the absence of high-luminosity X-ray dwarfs from at low redshifts (see Section 3.5). As we move to z = 1, we find a significantly higher M BH /M stellar ratio, suggesting that the BHs grew faster than their host galaxies were assembled in the simulation, similarly to what has been found in Illustris . Whether the observed scaling relations evolve at higher redshifts is still controversial. While some studies find that the M BH /M stellar , or equivalently the M BH /L stellar , ratio increases towards higher redshifts (e.g. Kormendy & Ho 2013;Bongiorno et al. 2014;Park et al. 2015;Yang et al. 2018;Ding et al. 2020), others find that there is no evolution up to z ∼ 1 − 2.5 (e.g. Schramm & Silverman 2013;Sun et al. 2015;Suh et al. 2020).
However note that, as with the local galaxies, these results are mainly just extrapolated from the high-mass end of the galaxy population. High-redshift studies of BHs in low-mass galaxies are only just at the beginning (see Sections 1 and 3.5). The next generation X-ray missions (Athena, Lynx) and deep spectroscopic surveys with the Roman Space Telescope and JWST will make it possible to obtain significant samples of AGN in low-mass galaxies at intermediate redshifts (see Greene et al. 2019). This will allow us to test our prediction that the M BH /M stellar ratio should increase with redshift.
Focusing now on the relationship between M BH and M gas , we note that the gas mass scaling relations have a similar slope to the stellar mass scaling relations. Furthermore, the M BH /M gas ratio likewise increases from z = 0 to z = 1. However, the M BH /M gas ratio has a significantly higher amount of scatter than the quite tightly correlated M BH − M stellar relation. Interestingly, BHs with large positive offsets from the mean M BH − M gas relation tend to have higher mass outflow rates, indicating that these overmassive AGN are able to drive more powerful outflows. For the stellar mass relation, the trend is weaker as the outflow activity also increases for higher stellar masses due to increased stellar feedback.
The increased outflow rates for overmassive BHs at a fixed gas mass provide an interesting hint that AGN can enhance galactic outflows once enough gas mass is available. The relationship between overmassive BHs and host galaxy properties is investigated more closely in the next section. Figure 4 shows various outflow quantities for the total gas at R vir (mass outflow rate M out , total mass flow rate M tot , outflow velocity . Outflow characteristics of the total gas and galaxy properties against total gas mass at z = 0 (upper panel) and z = 1 (lower panel) for the low-mass galaxy sample. In each case the binned mean relation is shown as a solid black curve with the filled black circles indicating the bin midpoints. The colour coding of the distribution indicates the offset from the M BH − M gas scaling relation from Figure 3, with blue for undermassive and red for overmassive BHs. 2D bins with at least ten objects are colour-coded according to the mean BH mass offset, and otherwise we plot the individual objects colour-coded by their respective BH mass offsets. Galaxies with overmassive BHs have increased outflow and inflow rates, leading to a bimodal distribution for the total mass flow rate. The outflows are faster and hotter in the overmassive regime. Overmassive BHs are also associated with reduced gas fractions and increased BH luminosities. Furthermore, the mass loading factor β is increased, whilst the sSFRs are suppressed.

Overmassive black holes and galaxy properties
v out , outflow temperature T out , mass loading factor β) as well as integrated galaxy properties (gas mass to stellar mass ratio f gas , BH luminosity L BH , specific SFR (sSFR)) plotted against M gas , both at z = 0 and at z = 1. The outflow properties at R vir are calculated within a spherical shell of width ∆r = 20 kpc. For reference, we also checked the outflow properties at 0.2R vir , where we used a shell of width ∆r = 10 kpc. See Section 2.4 for more details on the calculation of the outflow properties for the total gas.
We include the whole low-mass galaxy sample and colour-code the distribution by the offset from the mean M BH − M gas relation from Figure 3. Overmassive BHs, with BH masses above the mean BH mass <M BH > for a given gas mass are shown in red, whilst undermassive BHs are shown in blue. As in Figure 3, we show binned mean quantities where there are at least ten objects per bin, and otherwise the individual objects are plotted. In each case, we also show the binned mean relation, where there are at least ten objects per 0.15 dex gas mass bin, as a solid black line, with the bin midpoint values shown as filled black circles.
We plot the variation of outflow and galaxy properties with M BH at fixed gas mass as the gas supply is a crucial quantity for the effectiveness of the AGN. On the one hand, the AGN needs a sufficient amount of gas to be able to accrete efficiently. On the other hand, AGN feedback can only drive large-scale outflows if there is a sufficient amount of gas to couple to in the host galaxy.
To have a fair comparison, it is therefore important to keep M gas constant. Figure 4 shows that there are strong correlations between the offset from the M BH − M gas relation and all of the quantities considered at fixed gas mass. However, note that we would expect these correlations to be amplified by the tight relation between M BH and M stellar (see Figure 3), i.e. overmassive BHs at fixed gas mass will also tend to be associated with higher stellar masses. To ensure that the trends from Figure 4 are not just driven by increased stellar feedback activity, we show the same plot for overmassive BHs at In the lower panel, the colour coding indicates the offset from the mean M BH /M gas ratio of the respective stellar mass bin, with brown for above-average and green for below-average M BH /M gas ratios. 2D bins with at least ten objects are colour-coded according to the bins' mean offsets, and otherwise we plot the individual objects. At high redshifts, overmassive BHs are correlated with suppressed sSFRs across the whole stellar mass range (contrary to the low-redshift case, see Figure A1). The correlation is even stronger when the M BH /M gas ratio is considered, indicating that at high-redshift, the BHs are able to drive the gas out of dwarf galaxies and suppress star formation.
fixed stellar mass (see Appendix A). Broadly, the outflow trends are recovered at fixed stellar mass, though the separation of the sample into undermassive and overmassive BHs is less clear-cut. This is partly due to the small scatter in M BH at fixed M stellar . In addition, as discussed above, the AGN's ability to influence its host galaxy is tightly correlated with gas availability. Therefore not keeping the gas mass fixed further blurs the correlation, as there will be a wide range of gas masses for a given stellar mass. We will point out differences between the M gas -dependent and M stellar -dependent plots as we discuss the different quantities. The first column of Figure 4 shows the mass outflow rates as well as the total mass flow rates. At both redshifts, overmassive BHs are associated with increased outflow as well as increased inflow rates. For the total mass flow rate this leads to a bimodal distribution with both extreme inflow and outflow rates correlated with overmassive BHs. This bimodality is to be expected as on the one hand large gas inflow rates will allow BHs to grow more quickly and become overmassive, and on the other hand these overmassive BHs are then able to drive more powerful outflows than their undermassive counterparts. Outflow rates are higher at z = 1 than at z = 0, and in both cases increase towards higher gas masses. At z = 0, overmassive BHs reach outflow rates up to ∼ 50 M yr −1 , whilst at z = 1 these have outflow rates of up to ∼ 160 M yr −1 . Note that we obtain the same trends at fixed stellar mass, however the overmassive BH population is not as clearly separated and the variations in BH masses at fixed M stellar are smaller.
The second column shows the outflow velocity and temperature of the total gas at R vir . Outflows are on average faster and reach higher temperatures for overmassive BHs. Furthermore, the outflow velocities and temperatures are on average higher at z = 1 than at z = 0. This is particularly noteworthy considering that, at fixed halo mass, the average gas density within R vir will be higher at higher redshifts.
We note that when we inspect the same outflow properties at 0.2R vir , whilst the trends for the mass flow rates remain the same, the trends for v out and T out with BH mass offset are much more distinct. This suggests that the AGN mainly accelerates and heats the gas at a local level, while for the large-scale outflows the difference is somewhat less significant. This is especially true for the distribution at fixed M stellar where there is not a strong link between outflow temperatures and BH mass offsets at R vir , however, there is a clear correlation at 0.2R vir .
At z = 0, the mean outflow temperature stays relatively constant with increasing gas mass at ∼ 2 × 10 5 K. At z = 1, the outflow temperature increases with gas mass from ∼ 2 × 10 5 to ∼ 7 × 10 5 K. The outliers at both the low-temperature and the high-temperature end are dominated by undermassive BHs, suggesting that the link between BH mass and outflow temperature at R vir is slightly weaker. We found that the high-temperature outliers at 0.2R vir , on the other hand, are overmassive giving further support to the hypothesis that the AGN heating mainly affects the small-scale outflows.
For the z = 0 galaxies, the mean outflow velocity ranges from ∼ 10 − 40 km s −1 and for the z = 1 galaxies from ∼ 20 − 70 km s −1 . Galaxies with overmassive BHs at both redshifts reach typical outflow velocities of ∼ 100 km s −1 , with the outliers reaching up to ∼ 200 km s −1 .
Using the virial velocity V vir as a proxy for the escape velocity (as in Section 3.1), we find that at z = 0, the mean escape velocity increases with gas mass from ∼ 80 − 140 km s −1 , and at z = 1 from ∼ 100−220 km s −1 . Therefore, at both redshifts, the fastest outflows should be able to escape the halo, whilst for the systems that lie on the mean relation, the outflow velocities are not high enough to escape the halo. Though note that, as discussed in Section 3.1, the outflows are multiphase and that even if the outflow velocity for the total gas, as considered here, is too slow to escape the halo, the hot and fast outflow phase may still be able to escape.
At 0.2R vir , the outflow velocities are significantly higher with the mean outflow velocity ranging from ∼ 20 − 90 km s −1 at z = 0 and from ∼ 70 − 140 km s −1 at z = 1, implying that the outflows are decelerating. At both redshifts, typical outflow velocities for overmassive BHs are ∼ 200 − 300 km s −1 . The properties of AGN-driven outflows in the low-mass regime are still relatively unconstrained by observations due to the high resolution requirements. Manzano-King et al. (2019) presented the first direct detection and measurements of AGN-driven outflows in dwarf galaxies. They find outflow velocities from 375 − 1090 kms −1 for galaxies with masses ∼ 4 × 10 8 − 9 × 10 9 M . However, they measure these outflows in a relatively small region (1.5 − 3.0 kpc from the centre), whilst our outflow measurements at 0.2R vir correspond to length scales of ∼ 16 − 44 kpc in the dwarf regime. As we find that the outflows are decelerating, velocities in the above range could be plausible for the total gas in the central region. Though we cannot test this due to the resolution of , as the gravitational softening length is ∼ 3.5 kpc so the galaxies' central regions are not resolved in our simulation.
Next, we investigate the relation between overmassive BHs and gas fractions as well BH luminosities in the third column of Figure 4. f gas is reduced for large offsets from the mean BH mass, with overmassive BHs mostly found in galaxies with gas fractions f gas ≤ 1. This is in agreement with the observed suppressed gas fractions of low-mass galaxies with AGN activity (Bradford et al. 2018;Ellison et al. 2019). Furthermore, the overmassive BHs are  Figure 3, with blue for undermassive and red for overmassive BHs. Where there are fewer than ten objects per 2D bin, the individual objects are plotted. In cgs units, both momentum and energy outflow rates are clearly enhanced for overmassive BHs. However normalizing by the BH luminosity washes out this correlation at z = 1 due to strong coupling between BH masses and luminosities (see Figure 4). also correlated with high BH luminosities suggesting that these overmassive BHs are still actively accreting.
Finally, we investigate the properties related to star formation. The last column shows the mass loading factor at the virial radius β(R vir ) = M out (R vir )/SFR(≤ R vir ) and the sSFR = SFR/M stellar (the latter is evaluated within twice the stellar half mass radius). Note that at z = 0, 0.5 per cent of objects have a sSFR below 10 −11 (for 0.3 per cent of objects the sSFR is zero). These objects are also extremely gas-depleted and are not plotted here. Overmassive BHs are associated with higher mass loading factors and suppressed sSFRs for the whole gas mass range. For overmassive BHs at z = 0, we obtain sSFRs in the range ∼ 10 −11 − 2 × 10 10 yr −1 . Whilst at z = 1, where the mean sSFR is significantly higher, the majority of overmassive BHs have sSFRs in the range 10 −10 − 6 × 10 −10 yr −1 .
However, there are clear differences when we plot these two (star-formation related) quantities at fixed stellar mass instead (see Figure A1). For both β and the sSFR there is a turnover in the massive dwarf regime at ∼ 6 × 10 9 M at z = 0 and z = 1. Below this mass, the trends from Figure 4 are either reversed or washed out. This suggests that star formation in dwarfs is not affected by AGN feedback at low redshifts (z < 2).
For z ≥ 2, however, as the AGN reach higher luminosities (see Figure 2), there is a clear correlation between overmassive BHs and suppressed sSFRs across the whole stellar mass range. This is illustrated in Figure 5, where we plot the sSFR against stellar mass for z = 2, 3 and 4. We show the distribution of the sSFRs colourcoded by both the offset from the mean M BH − M stellar relation (upper panel) as well as by the offset from the mean M BH /M gas ratio for the respective stellar mass bin (lower panel). At high redshifts, overmassive BHs at fixed M stellar are associated with suppressed sSFRs down to the dwarf regime -contrary to the low-redshift case. The correlation is even stronger when considering the offset from the mean M BH /M gas ratio. This implies that at z 2, the AGN feedback can drive the gas out of dwarf galaxies and help regulate star formation, whilst at low redshifts, stellar feedback is the dominant quenching process. However, we note that based on the undermassive BHs at z = 0 compared to the observed scaling relations (see Figure 3) and the lack of high-luminosity X-ray in (see Figure 10), this quenching transition at z ∼ 2 could also be delayed to lower redshifts with more efficient BH growth at late times.
In observations, low-mass galaxies (with a typical stellar mass of M stellar ∼ 9 × 10 9 M ) with overmassive BHs (with respect to the M BH − σ relation) have been found to also have reduced star formation -although this correlation is less significant than for massive galaxies (Martín- Navarro & Mezcua 2018). This trend is also found in Illustris, across the whole galaxy mass range, whilst for IllustrisTNG the correlation is much weaker (Li et al. 2020). Furthermore, Sharma et al. (2020) find that isolated dwarf galaxies with overmassive BHs in R 25 experience star formation suppression starting at around z = 2.
Nevertheless, the enhancement of outflows by AGN is observed throughout cosmic time. We investigate the outflow properties more closely in Figure 6 where we plot the momentum and kinetic energy outflow rates for the total gas, as defined in Section 2.4, against gas mass. In the top row, we show these rates in units of the BH luminosity and in the bottom row in cgs units. We use the same colour-coding as in Figure 4 with red indicating overmassive BHs and blue indicating undermassive BHs with respect to the mean M BH − M gas relation.
Overmassive BHs have both higher momentum and higher energy outflow rates. However, once we normalise these rates by the BH luminosity the trends become weaker, as overmassive BHs also have higher BH luminosities (see Figure 4). This is particularly true for z = 1 where the BHs reach significantly higher luminosities and the BH offsets are correlated with the luminosities across a wider gas mass range.
We also investigated the momentum and energy outflow rates plotted against stellar mass (see Figure A2 in Appendix A). As with the previous plot, we recover similar trends at fixed stellar mass but the trends are weaker due to the variation in gas supply as well as the small scatter in BH masses at fixed M stellar . When we normalise the trends at fixed stellar mass by the BH luminosity, the correlation is reversed. This is because the correlation with BH luminosity at fixed stellar mass is much stronger (see Figure A1) and therefore able to invert the outflow trends.

Kinematic properties: mock MaNGA observations
In Section 3.3, we demonstrated that overmassive BHs drive more powerful outflows reaching higher velocities, temperatures, and mass loading factors. Next, we assess the potential impact of overmassive BHs on observations of low-mass galaxies. Penny et al. (2018) find that dwarf galaxies (which they define as galaxies with M stellar 5 × 10 9 M ) with AGN signatures are more likely to have an ionized gas component with large kinematic offsets with respect to the stellar component (difference between global kinematic PAs is 30°≤ ∆PA < 150°). They interpret this as indirect evidence for AGN-driven outflows which may be responsible for the kinematically misaligned ionized gas.
Here we test this hypothesis by producing mock MaNGA l.o.s. velocity maps for both the stars and the ionized gas following the procedure described in Section 2.5. We create these maps at three different redshifts (z = 0.0, 0.2, 0.4) with a pixel size of 0.5 arcsec (size of the MaNGA square spaxels) and we convolve the maps with a Gaussian filter of 2.5 arcsec, matching the MaNGA PSF. Further, we only include stars and gas within a 2D aperture with radius 1.5r eff , which corresponds to the spatial coverage of the primary MaNGA sample 5 . To check for the effect of orientation, we create the l.o.s. velocity projections for each galaxy at two fixed orientations: edge-on (inclination angle θ = 0°) and inclined by θ = 45°. We rotate the galaxies automatically by aligning L stellar with the vertical axis. Note that we only include galaxies that are rotationally supported. See Section 2.5 for more details on the mock l.o.s. velocity maps.
To assess the impact of the AGN-boosted outflows in galaxies with overmassive BHs, we select the galaxies below the 20 th percentile (P 20 ) and above the 80 th percentile (P 80 ) of the M BH /M gas distribution for each of the three stellar mass bins (dwarfs with 9.0 ≤ log(M stellar [M ]) < 9.5, massive dwarfs with 9.5 ≤ log(M stellar [M ]) < 10.0, and M * galaxies with 10.0 ≤ log(M stellar [M ]) < 10.5) at each redshift. Note that we only consider galaxies that are resolved by at least 50 star particles and 50 ionized gas cells within the 1.5r eff MaNGA aperture (see Section 2.5 for details). We then create mock MaNGA maps for these galaxies and determine the kinematic offsets between the ionized gas and the stars to assess whether heightened AGN activity produces a significantly higher percentage of kinematically misaligned galaxies. Figure 7 shows example mock MaNGA maps for galaxies from each of the three mass stellar bins at z = 0.0, 0.2, 0.4, demonstrating some of the qualitative features of interest. Note that all of the projections shown here are at an inclination of θ = 0°( edge-on). All of the example galaxies have overmassive BHs and would be classified as kinematically misaligned. In addition to the mock MaNGA maps, we also show the l.o.s. velocity maps of the whole gas at the original resolution of the simulation (projection dimensions: 0.2R vir × 0.2R vir × 2R vir ). The 1.5r eff MaNGA aperture is indicated by a dashed black circle. Furthermore, we show the gas density contours in light grey and the velocity streamlines in dark grey to highlight the inflowing and outflowing streams around the gas disc.
The kinematic PA is measured from the vertical axis and traces the position of the maximum change in velocity. We determine 5 Note that Henden et al. (2019a) find that galaxies with M stellar < 10 11 M have larger half-mass radii r eff than observed ones by roughly a factor of two. This means that we can better resolve galaxies with the MaNGA PSF than may be the case for the observed galaxies. the kinematic PA from our mock MaNGA maps for both ionized gas and stars using the fit_kinematic_pa routine (see Krajnović et al. 2006). The stellar l.o.s. velocity maps are dominated by the rotational motion of the disc. Since we have aligned L stellar with the vertical axis, we would therefore expect the kinematic PA of the stellar component to be PA ∼ 90°.
For the example galaxies shown here, this is recovered within ±22°. For the ionized gas component, however, we obtain large kinematic offsets, with respect to the stellar component, from ∆PA = 34°to ∆PA = 80°. This means that the method for determining kinematic PAs is sufficiently robust and the large kinematic offsets between stars and ionized gas are real features.
The gas kinematics at the original resolution of the simulation (third column of each panel) reveal complex inflow-outflow structures. Fast outflows can cause kinematic misalignment, see e.g. the example M * galaxy at z = 0.0 (Figure 7, left panel, third row). However, we identified a number of different physical processes other than outflows that can also result in kinematic misalignment. These include galactic fountains, as for the example dwarf galaxy at z = 0.2 (see Figure 7, middle panel, first row). Furthermore inflowing gas from mergers can also offset the rotational motion in projection, see e.g. the gas inflow from the lower left corner for the example M * galaxy at z = 0.4 (Figure 7, right panel, third row). Similarly, cosmic gas inflows can also cause misalignment, as for the example massive dwarf galaxy at z = 0.0 (Figure 7, left panel, second row). Finally, the entire gas disc can be misaligned so that misalignment stems from the rotational motion rather than gas inflows or outflows, as in the example dwarf galaxy at z = 0.0 (Figure 7, left panel, first row). This shows that kinematically misaligned ionized gas is not necessarily due to feedback activity alone, but misalignment still correlates with enhanced feedback (since inflows and mergers can also trigger feedback events).
To assess the relationship between overmassive BHs and kinematic misalignment more quantitatively, we systematically compare kinematic offsets for the undermassive (M BH /M gas below P 20 ) and the overmassive sample (M BH /M gas above P 80 ) in Figure 8. We plot the distributions of ∆PA for the edge-on, θ = 0°case (solid lines) and the θ = 45°case (dashed lines), with the shaded grey area indicating the misaligned regime. From left to right, the panels show the three different redshifts and from top to bottom, the panels show the three different stellar mass bins (as indicated by the plot titles). The plot titles also give the size of each sample for θ = 0°as well as θ = 45°, respectively.
We note that the massive dwarfs and M * galaxies sample sizes increase with decreasing redshift, whilst the size of the dwarf galaxies sample decreases. This is due to the above mentioned resolution criterion (at least 50 star particles and 50 ionized gas cells within 1.5r eff ) which forces us to discard an increasing number of dwarf galaxies as the amount of ionized gas decreases with redshift. Comparing the sample sizes for the two inclination angles, we can see that the θ = 45°sample size is greater than or equal to the θ = 0°s ample size in all cases. Again this is caused by the above mentioned resolution issues. When the galaxy is inclined towards the observer, we are likely to have a higher amount of gas in the l.o.s. due to the presence of outflows, so a higher number of galaxies will fulfill the resolution criterion at higher inclination.
Focusing first on the θ = 0°projections, we can see that in all cases, the galaxies above P 80 are either more likely (or for the M * galaxies at z = 0.4 equally likely, though note the small sample size) to be categorised as kinematically misaligned. Furthermore, these galaxies are also more likely to have extreme misalignments (60°≤ ∆PA < 120°). ) < 10.5, third row). All of the example galaxies shown here have a M BH /M gas ratio above the 80th percentile for their respective stellar mass bin and have been rotated to be viewed as edge-on. In some cases the kinematic misalignment is caused by fast outflows, however other factors, such as inflows or mergers, also play a role.
Moving on to the θ = 45°projections, we first note that in all cases these are more likely to be misaligned (or for the P 20 M * galaxies at z = 0.4 equally likely) compared to their θ = 0°c ounterparts, as it is easier for inflows or outflows to offset the projected gas kinematics when the galaxy is inclined towards the observer.
Comparing the P 20 sample to the P 80 sample at θ = 45°yields a more complex picture. In the majority of cases (six out of nine), the overmassive BHs are still more likely to be misaligned. However, in the remaining three cases, the P 20 galaxies are more or equally likely to be misaligned.
We interpret these distributions by noting that overmassive BHs are associated with both extreme mass inflow and outflow rates (see Figure 4) as well as higher outflow velocities. The higher outflow velocities make it easier to offset the rotational motion in projection leading to significant kinematic offsets. As we incline the galaxy disc towards the observer, the l.o.s. component of the rotational velocity decreases so outflows moving at slower velocities can also offset the rotational motion in projection, and the distinction between undermassive and overmassive BHs becomes less clear. We tested this hypothesis by inspecting the distribution of kinematic offsets split by v out instead of M BH /M gas and found that we obtain a similar distribution, demonstrating that outflow velocity is an important factor for shaping the ∆PA distribution -though outflows are not the only process generating kinematically misaligned gas (see Figure 7). Furthermore, observed differences between the distributions of ∆PA for AGN and non-AGN galaxies are also influenced by the inclination angle. We estimated inclination angles of the MaNGA galaxies from Penny et al. (2018) using the b/a axial ratios at r eff from the NASA-Sloan Atlas catalogue. We found that the median inclination angle (measured with respect the edge-on configuration) for the dwarfs without AGN signatures is θ ∼ 46°, whilst the median inclination angle for the dwarfs with AGN signatures is θ ∼ 55°, so this might slightly increase the difference in kinematic offsets between non-AGN and AGN dwarfs (though note that this difference in inclination angles is much smaller than the difference between the two angles we consider here). Also note that Penny et al. (2018) only considered quiescent dwarf galaxies. We do not impose a star formation based cut on our sample as this would leave us with too few objects, so this is a potential disparity as kinematic misalignment due to stellar feedback is likely more prominent in our sample.
The observed higher incidence of kinematically misaligned ionized gas for dwarf galaxies with AGN signatures could be caused by AGN-boosted outflow velocities (see also Duckworth et al. 2020), but also by other physical processes like cosmic inflows or mergers. Moreover, spectroscopic AGN signatures are easier to identify for galaxies which are inclined towards the observer, further increasing the likelihood of kinematic misalignment. To sum up, we find that kinematically misaligned gas is correlated with overmassive BHs but not always caused by AGN feedback itself.

Radiative properties: mock X-ray luminosities
As discussed in Section 1, there have been a number of systematic X-ray searches for AGN in dwarf galaxies. Mezcua et al. (2018)  Within each mass bin, we select the galaxies that fall below the 20 th percentile or above the 80 th percentile of the M BH /M gas ratio distribution. We show the PDFs of ∆PA for two fixed inclination angles: θ = 0°(edge-on) and θ = 45°. Galaxies with overmassive BHs are more likely to be offset across mass bins and redshifts. This difference is even more pronounced for θ = 0°, as rotational features are stronger and therefore only fast inflows or outflows can offset the rotational motion in projection.
(z 0.25) using the common footprint of the MPA-JHU catalogue and XMM DR7. Here we compare the X-ray properties of the dwarf galaxies to the these observations, focusing on the redshift evolution and the detectability of these BHs in dwarf galaxies using X-ray surveys.
Current X-ray studies searching for AGN in dwarf galaxies face two main challenges. Since BHs in dwarf galaxies have relatively low masses (M BH ∼ 10 5 M ), their luminosities are also correspondingly lower. This makes it difficult to identify AGN in dwarf galaxies, unless they have a high Eddington fraction. For example, Mezcua et al. (2018) estimate that 95 per cent of their AGN candidates have near to super Eddington accretion rates ( f Edd > 10 −2 ). Contrast this with Section 3.1, where we found that only one to ten per cent of low-mass galaxies have BHs with f Edd > 10 −2 between z = 0 and z = 1, though due to the efficient SN feedback at the low-mass end in , this might be a lower limit. Birchall et al. (2020), on the other hand, have a higher proportion of low-luminosity AGN in their sample, with a typical Eddington fraction of f Edd ∼ 10 −3 , as they focus on local galaxies. Furthermore, their technique for separating AGN from other Xray sources might allow them to identify more AGN sources (see below).
The second main issue for X-ray surveys is contamination from XRBs as well as hot gas emission. Due to the lower AGN luminosities, there is an overlap between the luminosity distributions of the contaminants and the AGN. In observational studies, it is therefore necessary to estimate the contributions from XRBs and hot gas to check that the integrated luminosity is significantly higher than the expected contribution from non-AGN sources. Birchall et al. (2020) use the relation from Lehmer et al. (2016) to estimate the XRB contribution based on M stellar and the SFR. For the hot gas contribution, they use the relation from Mineo et al. (2012) which gives the X-ray luminosity of the gas based on the SFR. They require that the total X-ray luminosity has to be three times higher than the estimated luminosity of the non-AGN sources (i.e. the putative AGN luminosity, L AGN , has to be twice as high as the luminosity of the contaminants, L XRB+Gas ). This leaves them with 61 AGN candidates out of the 86 X-ray active dwarf galaxies in their original sample. Note that in their sample, galaxies with significantly higher X-ray luminosities than expected from XRBs also exceed the expected contribution from hot gas, so the latter is less of an issue when identifying AGN.
The AGN candidates from Mezcua et al. (2018) are drawn from a parent sample of ∼ 2300 X-ray-selected type-2 AGN out . Distribution of the X-ray luminosities of BHs, XRBs, and hot gas across different redshifts. The top row shows the X-ray luminosity distributions for all dwarf galaxies (9.0 ≤ log(M stellar [M ]) < 9.5). The second row only includes the dwarfs that fulfil the AGN selection criterion from Birchall et al. (2020): L AGN > 2 × L XRB+Gas . We also examine the X-ray luminosity distributions using an even stricter criterion (L AGN > 10 × L XRB+Gas ) in the third row. In both cases, we give the fraction of dwarf galaxies with BHs that would get categorised as AGN using these criteria in the upper left hand corner. With the Birchall et al. (2020) criterion, a high fraction of BHs in dwarf galaxies are categorised as AGN, demonstrating that contamination from XRBs and hot gas should likely not pose a significant issue, even at the low-mass end. For comparison, we also show the X-ray distributions of the observed dwarf galaxies with AGN from Mezcua et al. (2018) and Birchall et al. (2020) in the bottom row. The local detections are in broad agreement with the distributions, though the observed AGN are shifted towards higher luminosities. From the luminosity distributions, we would expect many more AGN in dwarf galaxies to be detected by future X-ray surveys with higher sensitivities. to z ∼ 3 (see Suh et al. 2017). These source were classified by their spectroscopic type when available, or the photometric type. Mezcua et al. (2018) also check the contamination from XRBs using the Lehmer et al. (2016) relation. They find that the observed X-ray luminosities are more than ∼ 6σ larger than expected from XRBs, and using the Mineo et al. (2012) relation, they find that the X-ray luminosities of all their sources are more than ∼ 34σ above the luminosity expected from hot gas. So there is an extremely low chance of contamination -which is to be expected given the high Eddington fractions of the AGN candidates in this sample.
To study the redshift evolution of AGN and contaminant luminosities in , we plot the PDFs of the AGN, XRB, and hot gas Aird+18 Mezcua+18 Birchall+20 Figure 10. AGN fraction in dwarf galaxies (9.0 ≤ log(M stellar [M ]) < 9.5) as a function of redshift for two X-ray luminosity bins: 3.7 × 10 41 erg s −1 ≤ L X < 2.4 × 10 42 erg s −1 (upper panel) and L X ≥ 2.4 × 10 42 erg s −1 (lower panel). The AGN fractions for the fiducial run are shown as dark orange and dark blue hexagons, respectively. Furthermore, we plot the AGN fractions of three of the calibration runs, which tested different variants of the AGN feedback model, as grey hexagons. Observational data are also shown for comparison, as indicated by the legend. At low redshifts there are no objects in either of the X-ray luminosity bins, suggesting that the feedback model might suppress the occurrence of these highlyaccreting objects in the local Universe. At intermediate redshifts, is in agreement with observed X-ray AGN fractions. For z > 1, the simulation predicts that the AGN fraction should rapidly increase in both of the luminosity bins. X-ray luminosities in Figure 9. We use the bolometric corrections from Shen et al. (2020) to estimate the X-ray AGN luminosities in the 0.5 − 10 keV band, matching the band used by Mezcua et al. (2018). In analogy to the observational studies, we estimate the XRB and hot gas contributions from the relations by Lehmer et al. (2016) and Mineo et al. (2012), respectively. See Section 2.6 for more details on the calculation of the different X-ray luminosities.
The top row of Figure 9 shows the luminosity distributions of all dwarfs. There is a significant overlap between the luminosity distributions of AGN and XRBs. For the local galaxies, the contamination from XRBs is less significant as the sSFR is significantly lower (see Figure 4). Note that for all redshifts, the hot gas and the AGN contribution are well separated, with the hot gas distribution shifted towards significantly lower luminosities than the XRB contribution. In agreement with the observational studies, we therefore conclude that the hot gas contribution should not pose a significant issue for the identification of AGN in dwarf galaxies.
In the second row, we show the distribution of the X-ray luminosities for all dwarf galaxies that fulfil the Birchall et al. (2020) criterion (L AGN > 2 × L XRB+Gas ). In the top left corner, we give the fraction of dwarf galaxies with BHs, f detected , that fulfil this criterion. As expected from the full distribution, the highest fraction of AGN in dwarf galaxies are recovered for z = 0 with f detected = 0.43. This then decreases with redshift to f detected = 0.23 at z = 0.8 and then increases again, as the mean AGN luminosities move towards higher values with higher redshift, to f detected = 0.29 at z = 1.6. Note that since Birchall et al. (2020) identified 61 out of 86 X-ray active dwarf galaxies as AGN, they obtained an even higher AGN fraction using this criterion. This difference is likely linked to the harder X-ray band used in their study (2 − 12 keV) as well as completeness limits. Note that their original sample is shifted towards higher X-ray luminosities compared to and that it only has a few dwarf galaxies with X-ray luminosities below 10 39 erg s −1 .
In the third row, we show the distribution of the galaxies with an even stricter criterion of L AGN > 10 × L XRB+Gas . Whilst this efficiently separates the distributions of the AGN and the contaminants, it only classifies between three and seven per cent of sources as AGN.
In the bottom row, we plot the observed distributions from Mezcua et al. (2018) and Birchall et al. (2020). We show the Mezcua et al. (2018) galaxies in five different redshift bins, approximately corresponding to the redshifts considered. The last bin also includes a few higher redshift objects, as there is only a small number of AGN candidates at z ≥ 1.4. The Birchall et al. (2020) galaxies are plotted in the z < 0.2 panel, as all of their sources are at low redshifts, again note that the X-ray luminosities from Birchall et al. (2020) are based on a harder X-ray band (2 − 12 keV), so these distributions are not directly comparable.
We only include observed dwarf galaxies with M stellar ≥ 10 9 M in this plot to match the dwarf sample. The Mezcua et al. (2018) stellar masses were calculated from galaxy model templates based on the Chabrier IMF so we do not adjust these stellar masses. The stellar masses from Birchall et al. (2020) are derived using the Kroupa IMF, so we subtract 0.05 dex to adjust these to a Chabrier IMF.
The Mezcua et al. (2018) AGN luminosities are clearly separated from the contaminant distributions and more akin to the distributions using the stricter separation criterion of L AGN > 10 × L XRB+Gas . Though note that towards higher redshifts (z 1.0) the observed AGN luminosities have very high values (L AGN ∼ 10 43 erg s −1 ), demonstrating that the observations mainly pick up a few very bright objects. Interestingly, we have no objects with L AGN 10 43 erg s −1 at any of the redshifts considered here, indicating again that the observed dwarfs may be accreting more efficiently than in our sub-grid model. From our simulated X-ray luminosity distributions, we would expect many more candidates to be identified by future X-ray surveys with lower flux limits.
The local Birchall et al. (2020) AGN distribution is less clearly separated from the contaminants than the Mezcua et al. (2018) dwarfs. Though the AGN luminosity distribution is still shifted towards higher luminosities than the equivalent distribution. We also note that for both observational studies, there are no AGN sources with luminosities below 10 39 erg s −1 . This again suggests that AGN accreting at low Eddington fractions are missed by these surveys.
We therefore conclude that sensitivity limits will be the main issue for further X-ray searches for AGN in dwarf galaxies. Whilst emission from XRBs and hot gas is less of an issue and a high fraction of BHs in dwarf galaxies should be identified despite contamination from these sources (between ∼ 20 and ∼ 40 per cent).
Though note that the relations used here to estimate the XRB and hot gas contribution are derived from massive galaxies. As dwarf galaxies tend to have sub-solar metallicities they might have an enhanced high-mass XRB population (see discussion in e.g. Mezcua et al. 2018;Birchall et al. 2020). Lehmer et al. (2019) observe this enhancement for four dwarf galaxies, however, more data is needed to reliably estimate the XRB population in low-mass galaxies.
As X-ray searches are relatively unaffected by contamination from other sources, they can be used to establish an AGN fraction in different X-ray luminosity bins above the completeness limit of the respective survey. Next we compare these observed X-ray AGN fractions with AGN fractions from . In Figure 10, we show the AGN fraction against redshift in two different X-ray luminosity bins: 3.7 × 10 41 erg s −1 ≤ L X < 2.4 × 10 42 erg s −1 (upper panel) and L X ≥ 2.4 × 10 42 erg s −1 (lower panel).
For both of these luminosity bins, we plot the observed AGN fractions obtained by Aird et al. (2018), Mezcua et al. (2018), and Birchall et al. (2020). Though note again that the Birchall et al. (2020) luminosities were measured in a slightly different X-ray band (2 -12 keV band) whilst the other two studies use the 0.5 -10 keV band. As in the previous section, the luminosities are based on the 0.5 -10 keV band.
In addition, we show the AGN fraction obtained by Pardo et al. (2016) in the upper panel. Though this AGN fraction should be taken as an upper limit due to the relatively small sample, with over half of the detected AGN having X-ray luminosities below 10 41 erg s −1 . We also show the AGN fraction from Reines et al. (2013). However, this fraction should also be taken as an upper limit as it was derived from optically-selected dwarf galaxies. We plot this fraction as a green down-pointing triangle for comparison in both panels.
We show the simulation data from the fiducial run for the two luminosity bins as dark orange and dark blue hexagons, respectively. Here the AGN fraction is defined as the number of dwarfs with an overall X-ray luminosity L X = L AGN + L XRB+Gas in the respective luminosity bin normalised by the total number of dwarfs at the given redshift. We include the luminosity of the contaminants for consistency with the observations, though for these high-luminosity sources this makes little difference as the AGN is by far the dominant contribution (see Figure 9). Note that, as in the previous sections, we only consider central dwarf galaxies here so we might slightly overestimate the AGN fraction.
To get a sense of the uncertainty in the simulated AGN fractions, we also plot the AGN fractions obtained by three of the calibrations runs, where multiple variations of the AGN feedback parameters were tested, including variants without a duty cycle and weaker radio feedback (see Appendix A of Henden et al. 2018, for details). The AGN fractions from the runs with the alternative AGN feedback models are shown as grey hexagons.
For the 3.7 × 10 41 erg s −1 ≤ L X < 2.4 × 10 42 erg s −1 luminosity bin, there are no dwarfs at z = 0.0 for any of the AGN feedback variations. At z = 0.2, all of the three variant AGN feedback models produce AGN in this bin, with AGN fractions from ∼ 0.1 to 0.4 per cent. Although note that in absolute numbers, this corresponds to one to three AGN in this luminosity bin, so these AGN fractions will be affected by small-number statistics. The fiducial run only has AGN in this luminosity bin from z = 0.4 on-wards. We shade the plot area corresponding to z ≤ 0.1 grey to indicate that we have no objects in this redshift range. The absence of these objects could be due the outputs not being frequent enough 6 as accretion rates are highly variable and AGN might only spend a short time in the high accretion state. Furthermore, it could be that we do not find any bright objects as the sub-grid models prevent the AGN from reaching these high luminosities.
For the redshifts where there are both simulated and observed data points, 0.2 ≤ z ≤ 1, the AGN fractions from are in broad agreement with the observed fractions. Both the observed and the simulated data points are consistent with no significant evolution of the AGN fraction until z ∼ 1. From z = 1 on-wards, where we have no observational constraints, the simulations predict a significant rise in the AGN fraction from ∼ 1.5 per cent at z = 1 to ∼ 43 per cent at z = 4.
Next, we focus on the highest luminosity bin L X ≥ 2.4 × 10 42 erg s −1 in the lower panel. Here there are no objects below z = 0.4. At z = 0.4, only one of the calibration runs has any AGN in this luminosity bins, with an AGN fraction of ∼ 0.1 per cent, whilst the fiducial run only has AGN in this bin from z = 0.6. For 0.4 ≤ z ≤ 1.0 where we have both observational and simulated data points, these are in good agreement and both consistent with no significant evolution in the AGN fraction until z ∼ 1. Similar to the other luminosity bin, however, we see a significant increase in bright AGN at high redshifts, increasing from ∼ 0.2 per cent at z = 1 to ∼ 15 per cent at z = 4. We caution, however, that these AGN fractions are sensitive to the seeding model and the assumption that every DM halo above 5 × 10 10 h −1 M should host a BH (see Section 2.2.2 for details on the seeding model). The prediction that the high-luminosity AGN fraction in dwarf galaxies will significantly increase with redshift can be tested with upcoming X-ray missions. We estimate the X-ray luminosity limits across the redshift range using the flux limits for Athena 7 , and Lynx 8 . For a given redshift, we calculate the luminosity limit as L X = 4πd 2 L f X , where d L is the luminosity distance and f X is the flux limit. Both of the above instruments are more sensitive in the soft band (0.5 -2 keV) than in the hard band (2 -10 keV). Therefore, we use the soft and hard band flux limits to obtain lower and upper limits for detectable X-ray luminosities, respectively.
We find that Athena should be able to measure the AGN fraction for L X ≥ 2.4 × 10 42 erg s −1 up to z ∼ 1.5 − 3.2. The AGN fraction for 3.7 × 10 41 erg s −1 ≤ L X < 2.4 × 10 42 erg s −1 should be measurable up to z ∼ 0.7 − 1.5. Therefore, Athena should be able to detect the upturn in the AGN fraction at z 1.0 for the high-luminosity AGN, and it might also provide further constraints for the other luminosity bin.
Lynx, however, will be able to constrain the AGN fraction 6 The simulation outputs are spaced in redshift by ∆z = 0.2 for z ≤ 3 and by ∆z = 0.5 for 3 < z ≤ 4. 7 https://www.cosmos.esa.int/documents/400752/507693/ Athena_SciRd_iss1v5.pdf 8 https://wwwastro.msfc.nasa.gov/lynx/docs/ LynxInterimReport.pdf across the whole redshift range (0 ≤ z ≤ 4) for both of the luminosity bins considered here. At z = 4, we find that the soft X-ray luminosity limit is L X,soft ∼ 1.3 × 10 40 erg s −1 and the hard X-ray luminosity limit is L X,hard ∼ 1.6×10 40 erg s −1 , so Lynx will be able to observe the evolution of the AGN fraction even for low-luminosity AGN.

DISCUSSION
Dwarf galaxies are an important testbed of galaxy formation and cosmology. Notwithstanding massive theoretical effort, it remains challenging to match the observed abundances and structural properties of dwarfs, requiring either to revisit our physical models or to even think outside of the standard ΛCDM framework. For example, the majority of large-scale galaxy formation simulations employ strong SN feedback to reproduce the low-mass end of the observed GSMF (see e.g. Bower et al. 2012;Dubois et al. 2014;Vogelsberger et al. 2014a;Schaye et al. 2015;Pillepich et al. 2018). However, this approach is called into question by high-resolution (zoom-in) galaxy formation simulations (as well as small scale simulations of the ISM), which find it difficult to regulate star formation with SN feedback alone (e.g. Hopkins et al. 2014;Kim & Ostriker 2015;Kimm et al. 2015;Hu et al. 2016Hu et al. , 2017Marinacci et al. 2019;Smith et al. 2019) so that other sources of star formation regulation need to be invoked such as stellar winds, photo-heating, radiation pressure effects, cosmic rays and (MHD) turbulence.
Interestingly, when some of these additional physical processes are implemented alongside SN feedback, stellar feedback still struggles to regulate star formation in dwarf simulations (e.g. Kimm et al. 2018). To obtain efficient feedback in a cosmological setting, stellar feedback then needs to be enhanced, e.g. by boosting SN rates  or via high star formation efficiencies (Hopkins et al. 2018).
While strong stellar feedback seems required to match the observed dwarfs, it is likely to stunt AGN growth and feedback in these low-mass systems if a considerable gas reservoir is removed from the innermost regions (see e.g. Anglés-Alcázar et al. 2017b;Habouzit et al. 2017;Trebitsch et al. 2018;Habouzit et al. 2020). Moreover, it is worth noting that most galaxy formation simulations use the Bondi-Hoyle-Lyttleton rate to model BH accretion, which suppresses BH growth close to the seed mass due to its quadratic dependency on M BH . In cosmological simulations, seed masses of ∼ 10 5 − 10 6 M are typically employed, so that the slow growth phase then by construction of these models corresponds to the dwarf regime.
The suppression of BH growth in dwarf galaxies (both by strong SN feedback and the Bondi-like accretion prescription) has meant that so far AGN feedback in dwarfs has largely been neglected by simulators. Recent observations of AGN in dwarf galaxies question these assumptions, as contrary to the common theoretical models, there is a population of dwarf galaxies where AGN can accrete efficiently. AGN feedback in dwarf galaxies is an exciting possibility as it is under-explored theoretically, but very promising energetically (see Dashyan et al. 2018).
What is more, several BH seeding mechanisms should naturally predict BHs in dwarfs, for example remnants of Population III stars (e.g. Madau & Rees 2001;Heger et al. 2003;Volonteri et al. 2003;Heger & Woosley 2010;Whalen & Fryer 2012;Karlsson et al. 2013) Katz et al. 2015). Recent work has started exploring alternative seeding models in cosmological simulations, growing BHs from smaller seeds (see e.g. Habouzit et al. 2017;DeGraf & Sijacki 2020). However, observations cannot yet clearly distinguish between different seeding models (e.g. see discussion in Volonteri et al. 2008;Greene et al. 2019) and the true occupation fraction of BHs in dwarf galaxies is still unknown (see Section 1).
To investigate AGN in dwarf galaxies, in addition to improved seeding models, it will be necessary to explore accretion mechanisms beyond the widely adopted (yet very simplistic) Bondi model, with several models considered, such as a supply-limited accretion scheme (Beckmann et al. 2018(Beckmann et al. , 2019, torque-driven BH growth (Anglés-Alcázar et al. 2015, 2017a or accretion disc models (e.g. Power et al. 2011;Fiacconi et al. 2018;Bustamante & Springel 2019).
The modelling of the different AGN feedback channels is equally important. Here, with , we only consider simple, isotropic feedback, however observations find evidence for bipolar outflows from AGN (e.g. Rupke & Veilleux 2011;Maiolino et al. 2012). Note that AGN feedback in is modelled as mechanical feedback and a radiation field around the BH (see Section 2.2.3), but other processes that are not included here could also play an important role, such as outflows driven by radiation pressure (e.g. Bieri et al. 2017;Costa et al. 2018a,b) or jet-driven outflows. With a more realistic AGN feedback model, underestimating the BH luminosity would have an impact beyond just underestimating the mechanical AGN feedback and photoionization, since increasing the BH luminosity would also e.g. increase radiation pressure.
Even though under-predicts AGN activity in dwarf galaxies, as evidenced by the comparison with observed scaling relations and high-luminosity X-ray AGN fractions, there is still an effect by AGN feedback on dwarf galaxy properties in . This suggests that in reality the impact of AGN feedback on dwarfs could be even more significant (with the caveat that whilst the -BH growth is inefficient, the seeding is very efficient, see e.g. DeGraf & Sijacki 2020). Taking these results at face value, it is then inevitable to conclude that the effects of AGN feedback on the evolutionary history of dwarfs need to be explored, which may lead to re-evaluation of stellar feedback models as well.
To aid the development of theoretical models, it will be crucial to assess how many dwarfs harbour central BHs and how these BHs grow through cosmic time. Future observational facilities like Lynx or LISA will provide important observational constraints, as these instruments will probe the evolutionary history mapping both BH seeding and BH growth.
AGN feedback in dwarf galaxies might also contribute to the resolution of an on-going debate about the so-called 'cusp vs. core' problem. This controversy stems from the observed rotation curves of some dwarf galaxies, which seem to suggest cored DM density profiles rather than the cuspy profiles predicted by ΛCDM (Flores & Primack 1994;Moore 1994;de Blok et al. 2001;Gentile et al. 2004;Walker & Pẽarrubia 2011). Currently there is no consensus about the origin of the diversity of inferred DM density profiles. Some suggest that the inferred cored profiles could be due to uncertainties in circular velocity measurements (e.g. Marasco et al. 2018;Oman et al. 2019). Others have questioned the validity of ΛCDM and come up with alternative DM models -such as warm DM (e.g. Lovell et al. 2012) or self-interacting DM (e.g. Vogelsberger et al. 2014b). Yet others have invoked baryonic processes transforming cusps to cores via strong feedback (see e.g. Navarro et al. 1996;Governato et al. 2010;Di Cintio et al. 2013;Oñorbe et al. 2015;Fitts et al. 2017).
In , we find no significant differences in DM profiles for dwarf galaxies with overmassive BHs (not shown here), despite the clear impact on the baryonic component. This indicates that the AGN feedback is not strong enough to influence the DM component. Furthermore, this is likely also compounded with the resolution of the simulations. However, given that likely underestimates the luminosities of AGN in dwarf galaxies, AGN feedback may still play a role in transforming cusps into cores in dwarfs.

CONCLUSIONS
Recent observations have uncovered a population of dwarf galaxies hosting AGN, providing tantalizing hints that AGN feedback could also play a role in low-mass galaxies (see Section 1). We used the cosmological simulation suite to investigate the impact of AGN feedback on (central) low-mass galaxies (10 9.0 M ≤ M stellar < 10 10.5 M ), with a particular focus on dwarf galaxies (10 9.0 M ≤ M stellar < 10 9.5 M ). We examined the distribution of bolometric BH luminosities (Section 3.1) as well as the scaling relations (Section 3.2). Furthermore, we studied the correlations between overmassive BHs and host galaxy properties, in particular outflows, in Section 3.3. We also compared the low-mass galaxies to observational studies by creating mock MaNGA l.o.s. velocity maps (Section 3.4). Moreover, we estimated the X-ray luminosities for AGN, XRBs and hot gas to assess the detectability of AGN in dwarf galaxies with wide-field X-ray surveys (Section 3.5). Our most important findings are the following: (i) The majority of AGN in low-mass galaxies, and in particular in dwarf galaxies, are outshone by the stellar component. This renders these AGN difficult to detect. However, the high-redshift regime is more promising as more low-mass galaxies reach high Eddington fractions and luminosities.
(ii) The M BH −M stellar scaling relation for low-mass galaxies is broadly in agreement with observed late-type scaling relations, though note that the relation is undermassive compared to most of the observed relations. This indicates that SN feedback may be too strong in the low-mass systems stunting BH growth.
(iii) Low-mass galaxies with large positive offsets from the M BH − M gas or the M BH − M stellar relation have increased mass outflow and inflow rates. These overmassive BHs cause higher outflow velocities and hotter (albeit multiphase) outflows. While warm outflows typically generate galactic fountains, the hot outflow component is able to escape the host haloes, leading to a reduced gas reservoir. Future MUSE observations should be able to test our predicted outflow properties.
(iv) We find that quenching of dwarfs proceeds via two channels. At lower redshifts (z 2) where AGN accretion rates are low, SN feedback mainly regulates star formation in dwarf galaxies, whilst above z ∼ 2, AGN feedback is strong enough to suppress star formation in dwarfs as well. We note that these conclusions are sensitive to the strong SN model adopted, with AGN quenching possibly more important at lower redshifts, too.
(v) Galaxies with overmassive BHs are more likely to have an ionized gas component that is kinematically misaligned from the stellar component in mock MaNGA l.o.s. velocity maps in agreement with observations. While fast AGN-boosted outflows are partly responsible for this misalignment, we caution that other factors such as cosmic inflows or mergers are also important.
(vi) X-ray surveys are a promising tool for identifying AGN in dwarf galaxies. Even at the low-mass end, the AGN, XRB and hot gas X-ray luminosity distributions are sufficiently separated to identify a high fraction of AGN in dwarfs. Sensitivity remains the main issue, with the observed X-ray distributions shifted towards higher luminosities than . We predict that future X-ray surveys should uncover many more dwarf galaxies with AGN with lower Eddington fractions and at higher redshifts.
(vii) By comparing the occupation fraction of luminous AGN in dwarf galaxies in to observations, we find that luminous AGN are missing from at low redshifts, again indicating SN feedback-starved BH growth. Good agreement is, however, reached for intermediate redshifts. For high redshifts (z > 1), where there are currently no observational constraints, we predict that the fraction of luminous AGN in dwarf galaxies should rise rapidly.
Like the majority of galaxy formation simulations, has been designed such that SN feedback regulates the low-mass end of the galaxy mass function. Notwithstanding this, we find that AGN feedback in has a clear impact on the outflow properties of dwarf galaxies across cosmic time and contributes to quenching at high redshifts. Ongoing and upcoming state-of-the-art observations with e.g. JWST, MUSE, Athena, Lynx, ngVLA, and LISA will be able to probe the elusive AGN population in dwarfs to much lower luminosities and at higher redshifts. It is hence of paramount importance to develop next generation galaxy formation models with more realistic seeding models and more efficient AGN accretion in low-mass objects to explore this currently unknown territory, with the ultimate goal of elucidating the role of AGN in dwarf galaxies.

APPENDIX A: OVERMASSIVE BLACK HOLES AT FIXED STELLAR MASS
In Section 3.3, we illustrated how the offset from the M BH − M gas relation is correlated with the properties of the host galaxy, in particular with the outflow properties. We found that low-mass galaxies with overmassive BHs are associated with more powerful outflows, lower sSFRs and gas fractions, as well as higher BH luminosities.
We focused on the effect of varying BH mass at fixed gas mass, as gas abundance is the key driver for BH activity. On the one hand gas provides the fuel for the BH to grow through accretion, and on the other hand a certain amount of gas is needed to be able to drive outflows.
However note that overmassive BHs at fixed gas mass tend to also have larger stellar masses than their undermassive counterparts due to the tight relation between M BH and M stellar (see Figure 3). To ensure that the increased outflow activity is not just driven by stellar feedback, we present plots analogous to the ones in Section 3.3, but at fixed stellar mass. Figure A1, in analogy to Figure 4, shows outflow quantities at R vir (mass outflow rate M out , total mass flow rate M tot , outflow velocity v out , outflow temperature T out , mass loading factor β) as well as integrated galaxy properties (gas mass to stellar mass ratio f gas , BH luminosity L BH , sSFR) at z = 0 (top panel) and at z = 1 (bottom panel) plotted against stellar mass. . Outflow characteristics and galaxy properties against total stellar mass at z = 0 (upper panel) and z = 1 (lower panel) for the low-mass galaxy sample. In each case the binned mean relation is shown as a solid black curve with the filled black circles indicating the bin midpoints. The colour coding of the distribution indicates the offset from the M BH − M stellar scaling relation from Figure 3, with blue for undermassive and red for overmassive BHs. 2D bins with at least ten objects are colour-coded according to the mean BH mass offset, and otherwise we plot the individual objects colour-coded by their respective BH mass offsets. Galaxies with overmassive BHs have increased outflow and inflow rates, leading to a bimodal distribution for the total mass flow rate. The outflows are faster and hotter in the overmassive regime, similar to the fixed gas mass case. Overmassive BHs at fixed M stellar are also associated with reduced gas fractions and increased BH luminosities. The trends for the star-formation related properties (mass loading factor and sSFR), however are either washed out or inverted below M stellar 6 × 10 9 M compared to the fixed gas mass plots. See Figure 4 for corresponding plots at fixed gas mass.
We plot the distribution of the whole low-mass galaxy population. The colour-coding indicates the offset from the M BH − M stellar relation from Figure 3, with red for overmassive and blue for undermassive BHs. Where there are at least ten objects per 2D bin we show the mean BH mass offset of the respective bin, and otherwise we plot the individual objects. Note that the scatter in BH masses at fixed stellar mass is much smaller than at fixed gas mass (see Figure  3). We also plot the mean relation as a solid black curve, with the bin midpoints marked by filled black circles.
At fixed stellar mass, we recover the same outflow trends as at fixed gas mass, albeit with a weaker correlation due to the smaller scatter and differing gas reservoirs at fixed stellar mass. We find that overmassive BHs are correlated with both increased mass outflow rates and mass inflow rates, leading to a bimodal distribution for the total mass outflow rate (see Figure A1, first column).
Moving on to the second column, we find that overmassive BHs at fixed stellar mass tend to be associated with faster and hotter outflows, but the correlation is weaker than at fixed gas mass. We also inspected the outflow properties at smaller scales and found that there is a much stronger correlation between overmassive BHs and increased outflow velocities and temperatures at 0.2R vir . Similarly, we found that at fixed gas mass, the strength of the correlation between overmassive BHs and v out and T out is stronger at 0.2R vir . This suggests that the hot phase driven by the AGN is more dominant at smaller radii, indicating that some of the hot gas falls back as a galactic fountain.
The third column of Figure A1 shows that there is a clear correlation between overmassive BHs and suppressed gas fractions, in agreement with the results at fixed gas mass. Furthermore, overmassive BHs are generally associated with higher BH luminosities, though this trend is less significant for dwarf galaxies at z = 0.
The last column shows galaxy properties related to star formation: the mass loading factor β and the sSFR. Here we find the most significant difference between the trends at fixed gas mass and fixed . Momentum and kinetic energy outflow rates against total stellar mass at z = 0 (left panel) and z = 1 (right panel) for the low-mass galaxy sample. The top row is normalised by the BH luminosity whilst the bottom row shows the outflow quantities in cgs units. The binned mean relations are shown as solid black curves, with the bin midpoints indicated by the filled black circles. Colour coding indicates the offset from the M BH − M stellar scaling relation from Figure 3, with blue for undermassive and red for overmassive BHs. Where there are fewer than ten objects per 2D bin, the individual objects are plotted. In cgs units, both momentum and energy outflow rates are clearly enhanced for overmassive BHs. However normalizing by the BH luminosity washes out this correlation at both redshifts due to strong coupling between BH masses and luminosities (see Figure A1). See Figure 6 for corresponding plots at fixed gas mass. stellar mass. Whilst at fixed gas mass, we found that overmassive BHs are associated with increased β and suppressed sSFRs for the whole galaxy gas mass range, the trends are less clear at fixed stellar mass. Below M stellar 6×10 9 M , the correlation is either washed out or inverted. This suggests that for dwarf galaxies in there is no strong connection between AGN activity and quenching at low redshifts. Though note that the sub-grid model might underestimate AGN luminosities in this mass range (see Figure 10 in Section 3.5).
Finally, we also consider the momentum and kinetic energy outflow rates at fixed stellar mass in Figure A2 (see Figure 6 for the equivalent figure at fixed gas mass). We show these outflow rates at z = 0 (left panel) and z = 1 (right panel). The top row shows the momentum and energy outflow rates in units of the BH luminosity, whilst the bottom row shows the same quantities in cgs units. Again, the colour-coding indicates the offset from the M BH -M stellar relation. We show the mean bin values where there are at least ten objects per 2D bin, and otherwise plot the individual objects. The mean outflow rates are shown as solid black curves, with the bin midpoints marked by filled black circles.
Focusing first on the bottom row, we find that galaxies with overmassive BHs have increased momentum and energy outflow rates, at both redshifts. The correlation is less strong than at fixed gas mass -which we attribute again to the range of gas reservoirs at fixed M stellar and the small scatter in the M BH -M stellar relation.
However, when we normalise the momentum and energy outflow rates by the BH luminosity, the relationship is inverted (top row). This is because overmassive BHs are correlated with higher BH luminosities (see Figure A1), and this correlation is significantly stronger than the one for the outflow rates.
Overall, we find that the trends at fixed gas mass from Section 3.3 are still recovered at fixed stellar mass. An exception to this are the mass loading factor and the sSFR. At low redshifts, we only see a clear trend for the upper end of the massive dwarf regime (M stellar 6 × 10 9 M ) and the M * galaxies. This difference is most likely due to the strong SN feedback in which clears the gas out of dwarf galaxies, preventing the AGN from accreting efficiently. However, at higher redshifts (z ≥ 2), where accretion rates are higher and a significant fraction of BHs is in the quasar mode (see Figure 2), there is a clear correlation between overmassive BHs and suppressed sSFRs across the whole stellar mass range (see Figure 5). This indicates that if the AGN are able to accrete efficiently, they are able to regulate star formation in dwarf galaxies. This paper has been typeset from a T E X/L A T E X file prepared by the author.