Dwarf galaxies have historically posed challenges to the cold dark matter CDM model and while many of the socalled dwarf galaxy problems have been mitigated by incorporating baryonic processes the observed diversity of dwarf galaxy rotation curves remains a contentious topic. Meanwhile the growing observational samples of active galactic nuclei AGN in dwarf galaxies have prompted a paradigm shift in our understanding of dwarf galaxy evolution traditionally thought to be regulated by stellar feedback. In this study we explore the potential role of AGN feedback in shaping dark matter distributions and increasing the diversity of dwarf galaxy rotation curves using a new suite of cosmological zoomin simulations of dwarf galaxies with the FIRE3 model. Our findings indicate that the presence of active black holes BHs in dwarf galaxies can lead to diverse outcomes ranging from cuspier to more corelike profiles. This variability arises from the dual role of BHs in providing additional feedback and regulating the extent of stellar feedback. Consistent with previous research we find that AGN feedback is most impactful when cosmic ray CR modelling is included with CRs from any source significantly influencing dark matter profiles. Overall our results highlight that the interplay between stellar feedback BHs and CRs produces a broad spectrum of dark matter density profiles which align with observed correlations between rotation curve shapes and baryonic dominance. This underscores the importance of including the full range of baryonic processes in dwarf galaxy simulations to address the persistent smallscale challenges to the CDM paradigm.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - High Energy Astrophysical Phenomena']
Diverse dark matter profiles in FIRE dwarfs black holes cosmic rays and the cuspcore enigma
{'No Header': ', 1-22 (2024)', 'Diverse dark matter profiles in /f.pc/i.pc/r.pc/e.pc dwarfs: black holes, cosmic rays and the cusp-core enigma': "Sophie Koudmani, 1 , 2 , 3 ★ Douglas Rennehan, 3 Rachel S. Somerville, 3 Christopher C. Hayward, 3 Daniel Anglés-Alcázar, 4 , 3 Matthew E. Orr, 3 Isabel S. Sands 5 and Sarah Wellons 6  \n- 1 St Catharine's College, University of Cambridge, Trumpington Street, Cambridge CB2 1RL, UK\n- 2 Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 3 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA\n- 4 Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT, 06269, USA\n- 5 TAPIR, California Institute of Technology, Mailcode 350-17, Pasadena, CA 91125, USA\n- 6 Department of Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA  \nSubmitted to MNRAS", 'ABSTRACT': "Dwarf galaxies have historically posed challenges to the cold dark matter (CDM) model and, while many of the so-called 'dwarf galaxy problems' have been mitigated by incorporating baryonic processes, the observed diversity of dwarf galaxy rotation curves remains a contentious topic. Meanwhile, the growing observational samples of active galactic nuclei (AGN) in dwarf galaxies have prompted a paradigm shift in our understanding of dwarf galaxy evolution, traditionally thought to be regulated by stellar feedback. In this study, we explore the potential role of AGN feedback in shaping dark matter distributions and increasing the diversity of dwarf galaxy rotation curves, using a new suite of cosmological zoom-in simulations of dwarf galaxies with the /f.pc/i.pc/r.pc/e.pc-3 model. Our findings indicate that the presence of active black holes (BHs) in dwarf galaxies can lead to diverse outcomes, ranging from cuspier to more core-like profiles. This variability arises from the dual role of BHs in providing additional feedback and regulating the extent of stellar feedback. Consistent with previous research, we find that AGN feedback is most impactful when cosmic ray (CR) modelling is included, with CRs from any source significantly influencing dark matter profiles. Overall, our results highlight that the interplay between stellar feedback, BHs, and CRs produces a broad spectrum of dark matter density profiles, which align with observed correlations between rotation curve shapes and baryonic dominance. This underscores the importance of including the full range of baryonic processes in dwarf galaxy simulations to address the persistent 'small-scale challenges' to the CDM paradigm.  \nKey words: methods: numerical - galaxies: active - galaxies: dwarf - galaxies: evolution - galaxies: formation - dark matter", '1 INTRODUCTION': "Dwarf galaxies are pivotal in the Lambda - Cold Dark Matter ( Λ CDM) model, serving as the building blocks of hierarchical structure formation. Defined as galaxies with stellar masses below 3 × 10 9 M ⊙ (similar to the mass of the Large Magellanic Cloud), dwarf galaxies are crucial astrophysical probes in near-field cosmology and galaxy formation studies.  \nDue to their shallow potential wells, dwarf galaxies are ideal testbeds for studying galactic feedback processes ('baryonic physics'). They are also highly dark-matter dominated, making them valuable for testing dark matter models. However, their sensitivity to both dark matter and baryonic processes has led to several controversies, collectively known as the 'dwarf galaxy problems', which highlight discrepancies between observations and Λ CDM predictions from N-body simulations (see Sales et al. 2022, for a recent review). Key issues include the missing satellites problem, which  \nnotes a deficit in the observed number of Milky Way satellites compared to theoretical predictions (e.g. Kauffmann et al. 1993; Klypin et al. 1999; Moore et al. 1999), and the too-big-to-fail problem, where the central masses of the largest observed satellites do not match those of the most massive simulated subhaloes (e.g. BoylanKolchin et al. 2011). Additionally, the cusp-core problem refers to the observation that dwarf galaxy dark matter halo profiles are less centrally peaked or 'cuspy' than predicted by the CDM model with some dwarf galaxies instead exhibiting core-like profiles (e.g. Flores &Primack 1994; Moore 1994).  \nThese discrepancies have led to explorations of alternative dark matter models, such as warm dark matter (WDM) (Blumenthal et al. 1982). WDMmodels,whichinclude a free-streaming cut-off at dwarf galaxy scales, can address the missing satellite problem and reduce central halo densities, thus partially resolving the too-big-to-fail problem(Lovell et al. 2012). However, WDM struggles with the cuspcore problem, as the resulting cores are not large enough (e.g. Shao et al. 2013). Another approach involves self-interacting dark matter (SIDM), where particles scatter elastically, creating cored density  \nprofiles and reducing halo ellipticity (see Tulin & Yu 2018, for a review). Constraints from strong lensing measurements of galaxy clusters suggest a velocity-dependent cross section that decreases from dwarf to cluster scales (e.g. Meneghetti et al. 2001; Randall et al. 2008; Vogelsberger et al. 2012). Additionally, fuzzy dark matter (FDM), which proposes extremely low-mass particles like ultralight axions, produces cored central profiles and suppresses small-scale structures in the dwarf regime due to their long de Broglie wavelengths (see Niemeyer 2020, for a review).  \nHowever, there has also been a large body of theoretical work focusing on improving the baryonic physics modelling which may largely alleviate the 'dwarf galaxy problems' observed in pure CDM simulations. In particular, star formation suppression from reionization could significantly decrease the number of bright dwarf satellites (e.g. Efstathiou 1992; Okamoto et al. 2008; Fitts et al. 2017; Katz et al. 2020) thereby resolving the missing satellite problem. Furthermore, stellar feedback (e.g. Navarro et al. 1996; Governato et al. 2010; Parry et al. 2012; Brooks & Zolotov 2014; Hopkins et al. 2014, 2018; Chan et al. 2015; Kimm et al. 2015; Emerick et al. 2018; Smith et al. 2018, 2019; Gutcke et al. 2021; Jackson et al. 2024) may significantly decrease the central dark matter densities of dwarf galaxies. In particular, it has been demonstrated that cyclic supernova (SN) bursts play an important role in creating cores in dwarf galaxies (e.g. Pontzen & Governato 2012) and that late-time star formation is crucial for maintaining these cores (e.g. Read et al. 2019; Muni et al. 2024).  \nNotwithstanding this progress in alternative dark matter models and improved baryonic physics, some prominent discrepancies between simulations and observations still appear difficult to resolve without fine-tuning the models. In particular, theoretical models significantly struggle with explaining the observed diversity of dwarf rotation curves (see e.g. Santos-Santos et al. 2020, for a recent compilation). Due to the self-similar nature of structure formation in the Λ CDMuniverse, the full mass profile of a halo can be characterised by a single parameter such as the virial mass (Navarro et al. 1997). Hence dwarf galaxy rotation curve shapes are predicted to be almost identical for similar maximum circular velocities - and yet a great diversity is observed ranging from extremely cuspy profiles to large cores for the same inferred dark matter halo mass (see Oman et al. 2015). Galaxy formation simulations based on CDM have traditionally struggled to reproduce this observed diversity of dwarf rotation curves (Dutton et al. 2016; Sawala et al. 2016; Garrison-Kimmel et al. 2019; Applebaum et al. 2021). Some groups find that SIDM may be able to improve these discrepancies (e.g. Ren et al. 2019; Kaplinghat et al. 2020), whilst others find no preference between SIDM and CDM simulations with stellar feedback (e.g. Santos-Santos et al. 2020; Zentner et al. 2022).  \nThe most fundamental issue for CDM simulations lies in obtaining a balance in the burstiness of star formation histories, with some simulation set-ups predominantly producing smooth (or quenched) star formation histories and cusps (e.g. Vogelsberger et al. 2014; Sawala et al. 2016; Gutcke et al. 2022) whilst others predominantly produce bursty star formation histories and cores (e.g. Navarro et al. 1996; Read & Gilmore 2005; Governato et al. 2010; Pontzen & Governato 2012; Chan et al. 2015; Garrison-Kimmel et al. 2019).  \nOne potential solution could be additional feedback processes in the dwarf galaxy regime interacting with SNe and thereby modulating the star formation histories, leading to a natural diversity in rotation curves (see discussion in Garrison-Kimmel et al. 2019). Recently, the role of active galactic nuclei (AGN) in resolving the remaining dwarf galaxy problems has garnered significant interest. Traditionally, AGN feedback was only considered in massive galax-  \nies, with stellar feedback thought to dominate in low-mass galaxies. However, these models have been called into question by the growing observational samples of AGN in dwarfs, with detections spanning the whole electromagnetic spectrum from X-ray (e.g. Schramm et al. 2013; Baldassare et al. 2015, 2017; Lemons et al. 2015; Miller et al. 2015; Mezcua et al. 2016, 2018a; Pardo et al. 2016; Aird et al. 2018; Birchall et al. 2020, 2022; Latimer et al. 2021) to optical (e.g. Greene &Ho2004,2007;Reineset al. 2013; Chilingarian et al. 2018; Molina et al. 2021; Polimera et al. 2022) to IR (e.g. Satyapal et al. 2014; Sartori et al. 2015; Marleau et al. 2017; Kaviraj et al. 2019) and to radio observations (e.g. Greene et al. 2006; Wrobel & Ho 2006; Wrobel et al. 2008; Nyland et al. 2012, 2017; Reines & Deller 2012; Reines et al. 2014; Mezcua et al. 2018b, 2019; Reines et al. 2020; Davis et al. 2022). Intriguingly, recent targeted surveys indicate high AGN occupation fractions, on the order of 10 per cent, in dwarf galaxies, suggesting there may not be a drastic decline of AGN activity in the dwarf regime (e.g. Baldassare et al. 2015; Birchall et al. 2022; Bichang'a et al. 2024; Mezcua & Domínguez Sánchez 2024).  \nThese tantalising observations have motivated theorists to investigate this largely unexplored regime. Analytical models point towards favourable AGN energetics in the dwarf regime (compared to stellar feedback, see e.g. Dashyan et al. 2018) and suggest that AGN activity may be able resolve all of the remaining 'dwarf galaxy problems' (Silk 2017). However, hydrodynamical simulations have painted a more complex picture. AGN activity in dwarf galaxies may be significantly suppressed by stellar feedback evacuating gas from the central region (e.g. Dubois et al. 2015; Anglés-Alcázar et al. 2017b; Habouzit et al. 2017; Trebitsch et al. 2018; Hopkins et al. 2022a; Byrne et al. 2023a) and black holes (BHs) may be wandering in dwarf galaxies due to their shallow potential wells (e.g. Bellovary et al. 2021; Ma et al. 2021; Sharma et al. 2022; Beckmann et al. 2023) further decreasing their accretion efficiency. The modelling of BH growth is another crucial aspect as the fiducial Bondi model suppresses the growth of low-mass BHs compared to gas-supplylimited or torque-driven schemes due to its strong dependence on BH mass (e.g. Anglés-Alcázar et al. 2013, 2015; Koudmani et al. 2022; Wellons et al. 2023; Gordon et al. 2024). The cosmic environment also has a crucial influence on AGN activity in the dwarf regime with minor mergers triggering AGN episodes whilst longterm residence in dense environments is detrimental to BH growth (e.g. Kristensen et al. 2021).  \nDespite these caveats, various groups have identified physical regimes where AGN in dwarfs can accrete efficiently in accordance with observed samples. These dwarf AGN significantly influence their host galaxies by driving powerful outflows (e.g. Koudmani et al. 2019, 2022) and suppressing star formation (e.g. Barai & de Gouveia Dal Pino 2019; Sharma et al. 2020, 2022, 2023; Koudmani et al. 2021, 2022; Wellons et al. 2023; Arjona-Galvez et al. 2024). Given this potentially significant impact by AGN on the baryon cycle in dwarf galaxies, could AGN feedback also impact the dark matter distributions in dwarfs as predicted by analytical models?  \nThe interaction between AGN feedback, star formation, and dark matter in dwarf galaxies has not been thoroughly explored, though sometrends have been noted, such as AGN mildly suppressing central dark matter densities (see Koudmani et al. 2022; Arjona-Galvez et al. 2024), hinting at a possible role of AGN in driving cusp-to-core transformations in dwarfs. However, the central density suppressions by the AGN in these simulations were mostly minor and did not change the qualitative nature of the central profiles. This is likely due to the effective equation of state (Springel & Hernquist 2003) employed for modelling the interstellar medium (ISM) in these simulations, which leads to smoother star formation histories alongside less bursty AGN  \nfeedback at low redshifts, often operating as 'maintenance-mode' feedback (also see discussion in Koudmani et al. 2022). Recently, Arora et al. (2024) examined the role of BH feedback in core formation within the NIHAO simulations and found negligible effects in low-mass systems, largely due to inefficient BH growth in this mass regime, constrained by the Bondi accretion scheme. In contrast, they found significant impacts of BH feedback on core formation in more massive galaxies, aligning with earlier studies that show efficient AGN activity can create cores in these galaxies (e.g. Martizzi et al. 2013; Peirani et al. 2017; Macciò et al. 2020).  \nIt is therefore timely to investigate whether efficient AGN activity can induce cusp-to-core transformations in dwarf galaxies and potentially account for the observed diversity in dwarf rotation curves. Clearly, within this context it is crucial to model the multi-phase nature of the ISM but also to consider the role of non-thermal components such as magnetic fields or cosmic rays (CRs) (e.g. Uhlig et al. 2012; Booth et al. 2013; Hanasz et al. 2013; Salem & Bryan 2014; Pakmor et al. 2016; Farber et al. 2018; Holguin et al. 2019; Dashyan & Dubois 2020; Martin-Alvarez et al. 2020; Steinwandel et al. 2020). In particular, it has been shown that CRs may drive core formation (e.g. Hopkins et al. 2020b; Martin-Alvarez et al. 2023) though their impact in the dwarf regime remains controversial as gas-rich mergers may re-establish cores in these systems (see discussion in Martin-Alvarez et al. 2023). What is more, it has been found that CRs may significantly boost the efficiency of AGN feedback (e.g. Su et al. 2021, 2024; Byrne et al. 2023b; Wellons et al. 2023) by suppressing cold gas accretion from the CGM due to enhanced CR pressure support. Therefore, it is important to also assess the role of CRs when investigating cusp-core transformations with AGN.  \nThe aim of this paper is to systematically investigate the impact of AGNfeedback and CRs on dark matter profiles in dwarf galaxies and to assess whether incorporating these additional baryonic processes may be the key to resolving the diversity of rotation curve problem. To this end, we perform simulations with and without CRs and/or AGN of three different dwarf haloes with explicit multi-phase ISM physics within the /f.pc/i.pc/r.pc/e.pc-3 galaxy formation framework and contrast the resulting dark matter profiles, their cosmic evolution and their connection with different feedback processes.  \nThe remainder of this paper is structured as follows. In Section 2, we introduce our three dwarf simulation set-ups and describe the different variations of the /f.pc/i.pc/r.pc/e.pc-3 galaxy formation model explored in this work. The results of our analysis are presented in Section 3. Firstly, we focus on the stellar and BH properties of our simulated dwarfs, including visualisations of the stellar distributions in Section 3.1, comparisons with observed nearby low-mass galaxies (see Section 3.2) and observed BH masses in dwarfs (see Section 3.3). We then analyse how our different feedback set-ups influence the dark matter profiles at 𝑧 = 0 in Section 3.4 and examine the driving forces behind the diversity of dark matter profiles in 3.5. We investigate the cosmic evolution of the central dark matter densities in Section 3.6 and assess whether this could explain the observed diversity of rotation curves in Section 3.7. We discuss our results and caveats to our modelling in Section 4 and conclude in Section 5.", '2.1 General set-up': "Our simulation suite is based on a set of 3 initial conditions (ICs) of dwarf haloes selected from the /f.pc/i.pc/r.pc/e.pc ICs set, spanning a range of halo masses in the 'classical dwarf' regime. We select two low-mass  \ndwarfs: m10q and m10y , with 𝑧 = 0 halo masses of 𝑀 vir = 8 × 10 9 M ⊙ and 𝑀 vir = 1 . 4 × 10 10 M ⊙ , respectively. Our third halo, m10z , is an intermediate-mass dwarf with 𝑧 = 0 halo mass 𝑀 vir = 3 . 5 × 10 10 M ⊙ . All of these ICs were generated using the /m.pc/u.pc/s.pc/i.pc/c.pc IC generator (Hahn & Abel 2011) and have been thoroughly explored with the /f.pc/i.pc/r.pc/e.pc-2 model (Hopkins et al. 2018). In the standard /f.pc/i.pc/r.pc/e.pc-2 runs, m10q has a quiescent growth history and is relatively isolated, m10y forms early and has a large core, and m10z is an ultra-diffuse galaxy.  \nTo simulate these systems, we use the /g.pc/i.pc/z.pc/m.pc/o.pc hydrodynamics+gravity code, employing the mesh-free finite mass method (Lanson & Vila 2008a,b; Gaburov & Nitadori 2011; Hopkins 2015, 2017) that tracks constant-mass gas resolution elements in a Lagrangian fashion. Our simulations include both baryons and dark matter, with mass resolutions of 𝑀 gas = 262 M ⊙ and 𝑀 DM = 1303 M ⊙ , respectively. The gravitational softening for the dark matter, stars and BH particles is set to 𝜖 soft = 0 . 0625 ckpc ℎ -1 , whereas the softening for the gas resolution elements is set adaptively. This allows us to easily resolve core formation as demonstrated in previous work (e.g. Chan et al. 2015; Hopkins et al. 2018).", '2.2 Galaxy formation model': 'For the galaxy formation physics, we use the newly updated /f.pc/i.pc/r.pc/e.pc-3 model(Hopkins et al. 2023a) which tracks the state of the multi-phase ISM and multiple forms of stellar feedback, including feedback from SNee 1 of Type Ia and II, stellar winds from OB and AGB stars as well as multi-wavelength photo-heating and radiation pressure. Star formation occurs in gas that is molecular, dense, self-gravitating, self-shielding, and Jeans-unstable. All simulations include magnetic fields, using the magneto-hydrodynamics methods from Hopkins & Raives (2016) and Hopkins (2016). The /f.pc/i.pc/r.pc/e.pc-3 model also has new optional physics modules compared to the previous /f.pc/i.pc/r.pc/e.pc (Hopkins et al. 2014) and /f.pc/i.pc/r.pc/e.pc-2 (Hopkins et al. 2018) models. Two of the most important optional features are a novel sub-grid BH model, including seeding, accretion and feedback, as well as the inclusion of CR physics. To separate the impacts of CRs and BHs, we perform four simulations for each set of ICs, including runs without CRs or BHs, with both CRs and BHs and with only CRs or BHs, respectively. For the runs without BHs but with CRs, no BHs are seeded and the CRs only stem from SNe. For the runs without CRs and with BHs, we incorporate all AGN feedback channels except for CR injection and SNe do not inject CRs either. See Table 1 for an overview of the simulation suite.  \nThe modelling of BHs and CRs remains highly uncertain. Several studies have investigated the impact of different modelling assumptions for CRs (e.g. Hopkins et al. 2020b, 2023b; Butsky et al. 2023) and BHs (Anglés-Alcázar et al. 2017b, 2021; Su et al. 2021, 2024; Byrne et al. 2023a,b; Wellons et al. 2023) within the /f.pc/i.pc/r.pc/e.pc model. For our study, we base our BH model on the extensive study by Wellons et al. (2023) who identified a modelling space that allows for efficient BH feedback whilst still matching observational constraints, in particular with regards to star formation distributions. We summarise these BH modelling choices in Section 2.2.1 and point the interested reader to Wellons et al. (2023) for more details and an in-depth parameter study of alternatives to our choices. For the CR modelling, we employ the numerically-efficient subgrid CR model introduced  \n1 We couple the SN momentum following Hopkins et al. (2023a) and do not yet include the updates to the momentum coupling from Hopkins (2024). Though we note that these changes mainly affect high-mass galaxies at low resolution, whilst the impact on high-resolution dwarfs is very limited.  \nTable 1. Overview of dwarf zoom-in simulation suite. We list the simulation names (first column), ICs (second column), and whether BHs (third column) and CRs (fourth column) are included. We also provide the central dark matter densities at 𝑟 = 150 pc (fifth column), the ratios between the central dark matter density for the hydro runs and dark matter only simulations (sixth column), and the shape of the central dark matter density profile (seventh column) as shown in Figure 8. All density estimators are rounded to two significant figures.  \nin Hopkins et al. (2023b) allowing us to explore several different IC set-ups to 𝑧 = 0 whilst minimising the computational overhead costs that would arise from explicit CR transport. We summarise the key characteristics of the CR subgrid model in Section 2.2.2 and refer the interest reader to Hopkins et al. (2023b) for a more detailed overview of this approach.', '2.2.1 BH modelling': "For our BH modelling, we largely adopt the default /f.pc/i.pc/r.pc/e.pc-3 BH implementation (see Hopkins et al. 2023a), except for the AGN feedback injection where we use the continuous wind mode (rather than particle spawning) and for the CR injection where we employ the subgrid model from Hopkins et al. (2023b). Below we summarise the salient details as well as our modifications to the /f.pc/i.pc/r.pc/e.pc-3 BH model and refer the interested reader to Hopkins et al. (2023a) and Wellons et al. (2023) for a more in-depth description and parameter studies.  \nBHs may be seeded into the simulation from any star-forming gas particles that satisfy physically-motivated properties for the formation of low-mass seeds. In particular, the gas particles must have extremely high densities ( Σ gas ≳ 5000 M ⊙ pc 2 ) and low metallicities ( 𝑍 gas ≲ 0 . 001 Z ⊙ ) in order to spawn BH particles of mass 𝑀 BH = 100 M ⊙ . With a DM particle mass of 𝑀 DM ∼ 10 3 M ⊙ , we therefore do not fully resolve the dynamical friction forces that would drive the orbits of our seed BHs towards the centre of the galaxy. Hence we follow Wellons et al. (2023) and employ a subgrid prescription where the BHs are 'drifted' towards the local binding energy extremum of the stars, dark matter and BH particles so that they are not artificially ejected (also see Hopkins et al. 2023a). Similarly, since we cannot resolve the dynamics of binary BHs, two BHs are merged if their interaction kernels and force softenings overlap and they are gravitationally bound to one another.  \nOnce a BH particle is present, it grows via accretion by draining gas from the nearest 𝑁 ∼ 256 gas cell neighbours. Wellons et al. (2023) parameterized the accretion rate as / 𝑀 acc ≡ 𝜂 acc 𝑀 gas Ω , where 𝜂 acc is the accretion efficiency, 𝑀 gas is the gas mass within the BH kernel, and Ω = √︁ 𝐺𝑀 tot / 𝑅 is the orbital frequency determined at the kernel size 𝑅 , based on the total mass enclosed with 𝑅 , 𝑀 tot . /f.pc/i.pc/r.pc/e.pc-3 modulates the BH accretion rate via a gas reservoir that represents an accretion disc from which the BH particle grows at a mass growth rate / 𝑀 BH = 𝑀 𝛼 disc / 𝑡 dep , where 𝑀 𝛼 disc is the mass of the accretion disc and 𝑡 dep = 42 Myr ( 1 + 𝑀 BH / 𝑀 𝛼 disc ) 0 . 4 is the depletion time.  \nWhile Wellons et al. (2023) tested various distinct accretion models (i.e. variations in 𝜂 acc ), their best-fitting models favoured accretion that was powered by gravitational torques. Physically, these torques result from asymmetries in a galaxy's gravitational potential, or from interactions between (a) the dark matter and gas, (b) the stellar component and gas, or (c) the self-interaction of the gas with itself via shocks and dissipation (Hopkins & Quataert 2010, 2011). Theaccretion efficiency in the case of gravitational torque-dominated accretion comes from the work in Hopkins & Quataert (2011),  \n𝜂 acc = 𝐶 ([ 𝑀 BH + 𝑀 d ]/ 𝑀 d ) 1 / 6 1 + 3 𝑀 1 / 3 d , 9 ( 𝑀 gas / 𝑀 d ) , (1)  \nwhere 𝐶 is a resolution-dependent constant, 𝑀 BH is the BH mass, 𝑀 d is the mass of angular-momentum support material, and 𝑀 d , 9 = 𝑀 d / 10 9 M ⊙ . The model has been studied extensively in the context of cosmological zoom-in simulations in Anglés-Alcázar et al. (2013, 2017b,a, 2021) and Hopkins et al. (2016). It is important to note that we choose a value of 𝐶 = 0 . 1 in order to have reasonable BH masses within our three simulated dwarf galaxies. However, the parameter 𝐶 cannot take into account all sub-grid processes that may occur below the resolution scale.  \nOne important piece of (partially) unresolved physics is the impact of stellar feedback on the accretion flow in the vicinity of the BH. Wellons et al. (2023) introduce a scaling factor 𝑓 acc (such that 𝜂 acc → 𝑓 acc 𝜂 acc ) that depends on the local gravitational acceleration, 𝑎 g , eff = 𝐺𝑀 ( < 𝑅 )/ 𝑅 , as  \n𝑓 acc = 𝑎 g , eff 𝑎 g , eff + 𝑎 g , crit , (2)  \nwhere 𝑎 g , crit ≈ 10 -7 cms -2 . Below 𝑎 g , crit (an effective surface density of Σ tot ∼ 3000 M ⊙ pc 2 ), a significant fraction of gas is expelled due to stellar feedback (Grudić et al. 2019; Hopkins et al. 2022a). However, we emphasise that in addition to this stellar feedback effect on unresolved scales, our simulations explicitly capture the impact of stellar feedback on larger resolved scales (comparable to the size of the BH kernel).  \nThe accretion process itself drives powerful outflows, jets and radiation that influence the properties of the galaxies hosting the BHs, collectively coined BH feedback. The /f.pc/i.pc/r.pc/e.pc-3 BH feedback model includes three channels: radiative feedback, mechanical feedback, and CR feedback. Each mode injects energy at a rate proportional to the BH accretion rate / 𝑀 BH .  \nAsgasfalls into a BH, it loses gravitational energy that is converted into radiation, with an efficiency that is dependent on the location of the innermost stable circular orbit of the BH. While the location of that orbit depends on the spin of the BH, the typical efficiency used in the literature is 𝜖 r = 0 . 1, implying that 10% of the rest massenergy accretion rate is converted into radiation. Given that we do not track the BH spin evolution in our simulations, we also adopt this fiducial radiative efficiency value in our work. In /f.pc/i.pc/r.pc/e.pc-3, that radiation is treated in the same manner as the stellar feedback with multiband radiation transport and metallicity-dependent opacities using the LEBRON method (see Hopkins et al. 2020a), except using a quasar template spectrum (Shen et al. 2020). The radiation momentum flux depends on the luminosity absorbed ( 𝐿 abs ) in the gas as / 𝑝 = 𝐿 abs / 𝑐 .  \nThere are other outflows generated from the accretion process besides radiation-driven winds. For example, there may be jets or hydromagnetic winds from the accretion disc itself. The physics that drives these processes occurs on scales that are much below the resolution of our simulations and, therefore, need to be parameterized into a sub-grid model. The /f.pc/i.pc/r.pc/e.pc-3 BH model assumes that these winds have mass outflow rates proportional to the accretion rate. We assume a mass loading of one, such that / 𝑀 out = / 𝑀 BH . The kinetic energy rate in the wind itself depends on the wind velocity at launch, 𝑣 wind , such that / 𝐸 mech = / 𝑀 BH 𝑣 2 wind / 2. Wellons et al. (2023) recast the mechanical luminosity in terms of the new mechanical efficiency, 𝜂 mech ≡ 𝑣 2 wind /( 2 𝑐 2 ) , such that / 𝐸 mech = 𝜂 mech / 𝑀 BH 𝑐 2 . In this work, we choose 𝑣 wind = 10 , 000 km s -1 (see e.g. Fiore et al. 2017, compiled in fig. 3 of Rennehan et al. 2024) and, hence, have 𝜂 mech = 5 . 6 × 10 -4 .  \nThe /f.pc/i.pc/r.pc/e.pc-3 model also includes a CR model and, therefore, it also includes a BH feedback channel that accounts for relativistic massive particles accelerated by the AGN jets. We describe the CR propagation model in Section 2.2.2, but note here that the CRs are injected at a rate / 𝐸 CR = 𝜂 CR / 𝑀 BH 𝑐 2 , where 𝜂 CR is the efficiency of CRfeedback. In this work, we use an efficiency of 𝜂 CR = 0 . 01 which was identified as yielding realistic galaxy properties in the Wellons et al. (2023) parameter study.  \nFor the mechanical and radiative feedback, we use the solid angle weighting to inject the AGN feedback (akin to SN feedback injection, see Hopkins et al. 2018, for details). Following Su et al. (2021), we inject the majority of the AGN energy within a relatively narrow opening angle by weighting the AGN wind energy and radiation pressure received by each gas cell as  \n𝑤 ( 𝜃 ) = 𝜖 jet ( 𝜖 jet + cos 2 𝜃 ) ( 𝜖 jet + 1 )( 𝜖 jet + ( 1 -cos 2 𝜃 )) . (3)  \nWe set 𝜖 jet = 0 . 35 which leads to an effective base opening angle of ∼ 7 ° . The angle weighting is applied with respect to the total gas angular momentum accreted by the BH, 𝐽 gas , which is measured in the centre-of-mass frame of the BH - accretion disc sink particle. This allows for AGN energy injection without significantly disrupting the small-scale gas supply in the central region which favours repeated bursts and/or continuous AGN activity allowing us to probe the 'physically interesting' regime for AGN influencing cusp to core transformations. As discussed in Su et al. (2021), the large-scale collimation properties of the jet are not hugely sensitive to the small-scale opening angle since the feedback energy chooses the path of least resistance effectively recollimating even for large opening angles, also see Koudmani et al. (2019). Long-range radiation transport and CRs are treated via the LEBRON approximation and injected isotropically (see Hopkins et al. 2020a, 2023b).  \n2.2.2 CR modelling  \nThe /f.pc/i.pc/r.pc/e.pc-3 model includes a comprehensive CR model that explicitly follows the multi-species/multi-spectral CR dynamics (Hopkins et al. 2021a,b,c, 2023a) with the entire distribution function (Hopkins et al. 2022b) integrated with the magnetohydrodynamics (MHD) solver. However, there are also optional, simplified models that make a series of approximations valid in our regime of interest - dwarf galaxies. In particular, we use the sub-grid CR model described in (Hopkins et al. 2023b) and briefly describe our motivation for this model below.  \nTheprimary assumption is that the CR energy spectrum is spatially constant, and that CRs above ≳ 1 GeV dominate the CR pressure (e.g. Chanet al. 2019). This 'single-bin' model only evolves the integrated CRenergy density, drastically reducing the computational time while providing a reasonable approximation of the CR pressure, which is the dominant variable that impacts galaxy evolution. The study in Chan et al. (2019) showed that higher values of the isotropic diffusion coefficient 𝜅 are necessary to explain 𝛾 -ray luminosities in dwarf and 𝐿 ∗ -galaxies. Importantly, the values of the diffusion coefficient that best match the observations leads to most CRs escaping from galaxies with low gas density (i.e. dwarfs), while in starbursting galaxies with dense gas, the CR energy is lost almost entirely to collisions. Once the CRs escape the dense gas, the adiabatic gain/loss terms from the 𝑃 d 𝑉 work on the low density gas also becomes negligible.  \nGiven that several of the CR physical processes are seemingly subdominant in the dwarf regime, Hopkins et al. (2023b) made a series of further simplifying assumptions that are applicable in a 'streaming+diffusion' limit. In particular, the reduced sub-grid model assumes that the CR energy equation is in steady state ( 𝜕 𝑡 𝑒 cr → 0) and that the magnetic fields below the resolution scale are isotropically 'tangled'. The latter assumption allows the anisotropic diffusion tensor 𝜿 | | to be replaced with an isotropic diffusion coefficient 𝜅 , averaged over the resolution scale (and spatio-temporally constant). The dominant loss term for CR energy is assumed to be a combination of hadronic/pionic, Coulomb, and ionization loses - therefore, neglecting losses from diffusive reacceleration (Hopkins et al. 2022b). Streaming contributions from advective and Alfvén velocities are also ignored, leading to a two parameter model depending on the isotropic diffusion rate 𝜅 and a streaming velocity 𝑣 𝜅 . An approximate solution for the spatial distribution of CR energy density for such a model is presented in Hopkins et al. (2023b), and depends on an integral of an exponentially-declining function of the loss function.  \nThe resultant CR model is solved using a LEBRON-type approximation (see Hopkins et al. 2018, for details). The 'optical depth' of the CRs is equivalent the loss function we describe above, and the CR energy at a given distance from a source (i.e. SNe or AGN) is solved using the gravity tree in GIZMO . In /f.pc/i.pc/r.pc/e.pc-3, 10 per cent of the SN energy is converted to CRs and, as we mention above, 10 per cent of the BH luminosity is in the form of CRs. In this work, we use the values from Hopkins et al. (2023b); 𝜅 = 5 × 10 28 cm 2 s -1 and 𝑣 𝜅 = 20 km s -1 .", '3.1 Overview of the simulation suite': "We begin our analysis of the simulation suite by inspecting visualisations of our simulated dwarfs' stellar distributions. Figure 1 shows mock u/g/r Hubble images of the stars in each dwarf set-up, made by ray tracing through the gas: for each pixel, we sum up the luminosities  \nFigure 1. Mock u/g/r Hubble images showing the stellar distributions of our simulated dwarf galaxies. The side length of each projection is 10 kpc. The rows depict our three different IC set-ups, whilst the columns represent the different galaxy formation physics variations explored. The stellar sizes and colours are clearly impacted by the addition of CRs and/or AGN feedback hinting at a potential influence by these processes on dark matter distributions.  \n<!-- image -->  \nof the star particles in each band along the line of sight. Dust associated with the gas elements attenuates the light, assuming a Milky Way-like attenuation curve in each band and a fixed dust-to-metals ratio. The side length of each projection is 10 kpc.  \nThe three rows correspond to our three different ICs whilst the different columns represent our four physics variations based on the /f.pc/i.pc/r.pc/e.pc-3 physics modules: the fiducial /f.pc/i.pc/r.pc/e.pc-3 set-up (without CRs or AGN), /f.pc/i.pc/r.pc/e.pc-3 with AGN, /f.pc/i.pc/r.pc/e.pc-3 with CRs and /f.pc/i.pc/r.pc/e.pc-3 with both CRs and AGN. We note that there are several uncertainties associated with modelling AGNandCRsingalaxyformationsimulationsandherewe only explore a small subset of this parameter space that was identified by Wellons et al. (2023) as producing realistic galaxy properties. Furthermore, our simulations likely represent an upper limit for the impact of CRs since we are using an effective subgrid model for the CR population (Hopkins et al. 2023b) which only accounts for local losses and assumes a constant diffusion and effective streaming speed. See Section 2.2.2 and Hopkins et al. (2023b) for more details.  \nKeeping these caveats in mind, we note that both CRs and AGN have a significant impact on the stellar distributions. Generally, the addition of CRs leads to more compact stellar distributions, whereas the simulations without CRs also have an extended diffuse component. This is partly due to CRs suppressing star formation as the additional pressure support prevents gas accretion from the CGM onto the galaxy (also see Hopkins et al. 2020b, 2023b; Farcy et al. 2022; DeFelippis et al. 2024; Thomas et al. 2024). For the runs with AGN feedback, the picture becomes even more complex. The runs without CRs but with AGN appear qualitatively similar to their non-  \nAGN counterparts if not slightly more diffuse. The runs with CR, however, display systematic differences with added AGN feedback. This is consistent with the analysis from Wellons et al. (2023) where they find that CR-driven feedback is among the most efficient AGN feedback channels. In all cases, the runs with AGN are more compact and redder, pointing towards older stellar populations. Overall, these visualisations demonstrates that both AGN and CRs have a significant non-linear impact on dwarf galaxy evolution.", '3.2 Stellar assembly': 'Before assessing the impact of our baryonic modelling choices on dark matter distributions, it is vital to confirm that our feedback implementations produce realistic dwarf galaxy properties. To this end, we assess the stellar properties of our simulated dwarfs and compare these with observational constraints in Fig. 2, which shows the stellar mass - halo mass (SMHM) relation, galaxy colours and stellar half mass radii. Note that in this and later plots the symbol style indicates the model variations explored in our simulation suite, as indicated by the legend.', '3.2.1 SMHM relation': "Firstly, we examine if the simulated dwarf galaxies align with observational expectations for the SMHM relation. The left panel of Fig. 2 shows our 12 dwarf zoom-in simulations in 𝑀 halo -𝑀 stellar space at redshift 𝑧 = 0.  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Integrated stellar properties of our simulations in comparison with observational constraints. The symbol style in this and later plots indicates the model variations. Left panel: Stellar mass - halo mass (SMHM) relation. SMHM relations based on empirical models as well as individual observed galaxies are plotted for comparison as indicated by the legend. Middle panel: (Intrinsic) simulated galaxy colours against stellar mass. The colours from the SDSS galaxies are indicated as grey contours and the integrated photometry of local dwarfs is plotted as grey squares (Mateo 1998) and grey diamonds (Carlsten et al. 2022), for comparison. Right panel: Stellar half-mass radii against stellar mass with the data from local dwarfs (Carlsten et al. 2022) shown as grey diamonds for comparison. Overall , though the addition of AGN and/or CRs introduces significant scatter in the simulated distributions, our simulated stellar properties are in good agreement with the observations for all simulation set-ups explored.  \n<!-- image -->  \nFor comparison, we plot several SMHM relations from the literature, including Behroozi et al. (2013, 2019) and Moster et al. (2018). The low-mass end of the SMHM relation remains largely unconstrained by observations, evident from the significant differences between these models. The Behroozi et al. (2013, 2019) models extend down to log ( 𝑀 halo / M ⊙ ) = 10 . 0, and the Moster et al. (2018) relation was shown to also be valid down to the dwarf regime, at least to log ( 𝑀 halo / M ⊙ ) = 10 . 0, by O'Leary et al. (2023). We plot the extrapolated relations by Behroozi et al. (2013, 2019) and Moster et al. (2018) as dashed and dotted lines, respectively. Extrapolating SMHM relations below the minimum mass considered has pitfalls, especially regarding physical processes like reionization suppression, which are effective at low masses but not included in most of the models. The scatter is especially difficult to model in this regime since more bursty star formation histories at lower stellar masses are expected to increase the scatter whilst high-redshift quenching due to reionization suppression may act to substantially reduce the scatter in the dwarf regime (O'Leary et al. 2023).  \nOur simulations with the m10q ICs, the least massive halo, are firmly in the extrapolated regime, whilst m10y and m10z are at the lower end of validity for the Behroozi et al. (2013, 2019) and Moster et al. (2018) relations. Hence for our simulations comparing to SMHM relations constructed for low-mass galaxy samples, rather than the fiducial SMHM relations for the general galaxy population, is more instructive. We include SMHM relations derived for dwarf galaxies by Read et al. (2017), Jethwa et al. (2018), and Nadler et al. (2020) as shaded regions. The differences in these relations arise due to different assumptions for surface-brightness incompleteness correction, see discussion in Behroozi et al. (2019). We also show various observed data points from the SPARC data base (Li et al. 2020) to indicate the scatter in the observations.  \nWe note that the empirical relations as well as the SPARC data base use different definitions for the halo mass, either employing the virial mass as defined by the mass enclosed in a sphere whose mean density is 200 times the critical density of the Universe at 𝑧 = 0 or the halo mass based on the collapse of a spherical top-hat perturbation following Bryan & Norman (1998). We have checked for our  \nsimulated dwarfs that the difference between these two definitions is at most 0.1 dex, and plot the simulated halo masses as the average value between these two definitions. The simulated stellar mass is calculated within twice the stellar half mass radius. Given the significant uncertainty and scatter in SMHM relations, we do not adjust the stellar masses for different radial cut-offs or initial mass functions.  \nFrom this comparison with the empirical relations and SPARC observations, none of our models can be ruled out based on their integrated stellar masses. All our simulations fall within the range of the empirical SMHM relations and are well within the scatter of the observed data. The simulation set-ups with CRs (indicated by squares and diamonds) are generally less efficient at forming stars than the simulations without CRs (indicated by stars and crosses). Most simulations fall towards the upper end of the predicted SMHM relations, with the m10zCR+AGN simulation the only set-up resulting in an 'undermassive' galaxy. Though we note that even this set-up is still within the scatter of the most-up-to-date dwarf galaxy SMHM relation by Nadler et al. (2020).", '3.2.2 Galaxy colours': 'Next, we compare the galaxy colours of our simulated dwarfs with observational data. The middle panel of Fig. 2 presents the 𝑔 -𝑟 colour versus stellar mass for our simulated dwarfs. For comparison, we display the colours of SDSS galaxies from the NASA Sloan Atlas (NSA) as a grey contours, the colours of Milky-Way dwarf satellites from Mateo (1998) as grey squares and the ELVES dwarf satellites, a nearly volume-limited sample of Milky Way-like hosts in the Local Volume, from Carlsten et al. (2022) as grey diamonds.  \nWe use the Bruzual & Charlot (2003) stellar population synthesis model to calculate galaxy colours as a function of stellar age and metallicity, assuming a Kroupa (2001) IMF, including all star particles within twice the stellar half-mass radius. For the MW dwarf satellites, we convert the 𝐵 -𝑉 colours using the transformation equations from Jester et al. (2005), noting that while these equations are for stars, they should approximate galaxies well unless strong emission lines are present. One important caveat is that our simu-  \ngalaxy colours do not account for dust attenuation. However, we checked that, as expected for dwarf galaxies, our simulations have relatively low metallicities and therefore dust attenuation is expected to be minimal.  \nAll simulated dwarfs are overall in good agreement with both the local dwarf galaxy constraints and the colour distribution of SDSS galaxies. Although we note that SDSS is very incomplete in the mass regime we consider here and the ELVES dwarfs were selected to be satellites while our dwarfs are field objects. The SDSS galaxies cluster around 𝑔 -𝑟 ∼ 0 . 35, whilst most of our dwarfs are offset towards somewhat redder colours, clustering around 𝑔 -𝑟 ∼ 0 . 45. Accounting for dust attenuation would likely make this offset more severe. However, this discrepancy is also inherently linked with completeness bounds in SDSS in this mass range, with blue galaxies likely overrepresented.  \nInterestingly, the m10qCR set-up has a similar stellar mass to its AGNcounterpart m10qCR+AGN , but is shifted towards bluer colours due to a late-time burst in star formation which is also visible in the stellar projections in Figure 1. Conversely, the m10q set-up has a lower stellar mass and is redder than its AGN counterpart, indicating a complex interplay between CRs and AGN activity. For m10y , neither CRs nor the AGN have a strong impact on the integrated galaxy colours, though there are small, yet systematic, offsets with the CR runs being redder than the no CR runs, and the AGN runs being redder than their counterparts without AGN. Finally, for m10z , the set-up with neither CRs nor AGN yields the bluest galaxy colour, m10z , whilst the three other set-ups all have relatively similar colours, despite spanning an order of magnitude in stellar mass, indicating that all three of these formed the majority of their stars early in cosmic history.', '3.2.3 Stellar sizes': 'We now turn towards examining the stellar sizes which provide crucial clues to both the star formation and dynamical histories of our dwarfs. We plot the stellar half mass radii of our simulated dwarfs at 𝑧 = 0 as colour-coded symbols. For comparison, we also plot the observed effective radii of local dwarf galaxies as grey diamonds based on the ELVES catalogue from Carlsten et al. (2022), which also includes the Milky Way dwarf measurements from McConnachie (2012) as a subsample. To make this comparison more consistent, we calculate the simulated stellar half mass radii based on the averages of 2D projections along random sightlines. Again there is appreciable scatter in the observations and our simulated dwarfs fall well within this scatter.  \nAs we found from our visual inspection of the stellar distributions in Fig. 1, the stellar distributions are generally more compact with CRs. The impact of AGN feedback on the stellar sizes is more complex, leading to both more extended distributions (e.g., compare m10yCR with m10yCR+AGN ) or more compact distributions (e.g., compare m10qCR with m10qCR+AGN ). The differences between the AGN and no-AGN set-ups are more pronounced if CRs are present, in agreement with the findings from Wellons et al. (2023) that CRs enhance the efficacy of the AGN feedback in /f.pc/i.pc/r.pc/e.pc.  \nOverall our stellar sizes are in good agreement with the observations, though the m10qCR+AGN simulation is slightly more compact than local dwarfs at similar stellar masses. However, we note that the observed radii are based on stellar light whilst we are calculating the simulated radii based on the stellar mass distribution. Observational effects can add significant scatter leading to both more compact half-light radii (e.g. Klein et al. 2024) or more extended observed distributions (e.g. Parsotan et al. 2021), with the stellar sizes also  \nFigure 3. BHmass - stellar mass scaling relations compared to observational constraints. We show the redshift evolution of the dwarf zoom-in simulations from 𝑧 = 4 to 𝑧 = 0 and plot data from observed dwarfs and globular clusters as indicated by the figure legends. We also show the observational BH mass - stellar scaling relations from the literature by Reines & Volonteri (2015) and Greene et al. (2020), greyed out in the extrapolated regime, as well as the high-redshift relation from Pacucci et al. (2023) and NSC mass - stellar mass relation from Neumayer et al. (2020), which can be taken as an upper limit. Our simulated BHs are generally overmassive compared to the (heavily) extrapolated local relations though they are well within the limits set by the NSC relation and their offsets are much smaller than those uncovered by JWST at high redshifts.  \n<!-- image -->  \nbeing significantly dependent on the wavelength (e.g. Cochrane et al. 2023).', '3.3 BH assembly': "Having analysed the stellar properties, we turn to analyse the evolution of the BHs in our simulations in the context of observational constraints.  \nFig. 3 shows our simulated dwarfs in stellar mass - BH mass space. The evolutionary tracks of the dwarf simulations from 𝑧 = 4 to 𝑧 = 0 are plotted as dotted lines, colour-coded by simulation setup. Note that the BHs are generally seeded around 𝑧 ∼ 8-10 with the exact seeding time depending on the metallicity and density conditions. However, we here focus on later redshifts as by 𝑧 = 4, the most massive progenitor dwarf hosts a BH for all model variations so that we can track these alongside the stellar mass growth of the main progenitor. Integer redshifts are marked with the marker symbols for the respective simulation set-ups as indicated in the legend.  \nWe note that there is substantial uncertainty at the low-mass end of the BH - galaxy scaling relations due to the difficulty in reliably measuring the masses of intermediate-mass BHs (IMBHs) which have a much smaller gravitational sphere of influence and generally lower luminosities. Indeed there are no BH mass measurements for dwarf galaxies in the stellar mass range covered by our simulations, 10 6 M ⊙ ≲ 𝑀 stellar ≲ 10 8 M ⊙ . For more massive dwarfs there are only very limited measurements and we show these observed BH masses in massive dwarfs for reference, including virial BH mass  \nestimates for active dwarf galaxies from Baldassare et al. (2020) and Wasleske & Baldassare (2024) as dark-grey squares and stars, respectively, and measurements from Greene et al. (2020) (based on dynamics, reverberation mapping or scaling from the stellar velocity dispersion in the specific case of TDE events that seem to have reliable determinations of the host galaxies) as dark-grey hexagons and triangles for mean values and upper limits, respectively. Error bars are omitted for clarity. We note that there are no dynamical nuclear IMBH mass measurements below 𝑀 stellar ∼ 10 9 M ⊙ and, excluding upper limits, the lowest BH mass estimate based on TDEs stems from a host galaxy of stellar mass 𝑀 stellar = 2 . 5 × 10 8 M ⊙ . Belowthis mass, we are firmly in the extrapolated regime and indicate this by greying out the local BH - galaxy scaling relations below this mass.  \nTo extend our comparison to lower-mass systems, we include IMBH candidates from globular clusters around the Milky Way and M31, compiled by Greene et al. (2020), plotted in light-grey. Some globular clusters, like 𝜔 Cen, are hypothesised to be remnants of dwarf galaxies disrupted by the Milky Way's tidal field (e.g. Freeman 1993; Bekki & Freeman 2003; Meza et al. 2005). The presence of IMBHs in globular clusters is controversial, with significant uncertainty in BH mass values. Thus, we show a range for each cluster, indicating the lowest and highest BH mass values reported in the literature, with crosses for mean values and triangles for lower/upper limits. In addition, we plot an updated firm lower limit for the BH mass in 𝜔 Cen which was recently determined by Häberle et al. (2024) and a new IMBH detection in M31's most massive globular cluster which may also originate from a stripped dwarf galaxy (Pechetti et al. 2022), both indicated in olive-green.  \nWe also show several observed BH mass - stellar mass scaling relations from the literature to give an indication of the scatter in the extrapolated relations. The Greene et al. (2020) relations are based on dynamical measurements by Kormendy & Ho (2013) as well as more recent dynamical measurements and upper limits in the lowmass regime. We plot the relation for the total galaxy sample as a solid black line, and the early-type and late-type galaxy relations as loosely and densely dashed black lines, respectively. The latetype relation has a significantly lower normalization, likely due to a weaker correlation between disc mass and BH mass (e.g. Kormendy & Richstone 1995). We also show the Reines & Volonteri (2015) scaling relation, based on local AGN, including dwarfs, as a dasheddotted black line. This relation is significantly shallower than that of Greene et al. (2020), leading to larger BH masses in the low-mass dwarf regime.  \nFor the classical dwarf mass range covered by our simulations, 𝑀 stellar ≲ 10 8 M ⊙ , we would expect a flattening of the scaling relations, with the transition mass dependent on the seeding mechanism (see e.g. Greene et al. 2020), which sets a lower limit for BH masses. The extrapolated BH - galaxy scaling relations may hence be considered lower limits for BH masses in the dwarf regime. We also plot the high-redshift scaling relation derived by Pacucci et al. (2023) from overmassive BHs uncovered by JWST to indicate the high BH mass to stellar mass ratios that can be reached with optimal growth conditions in the early Universe (modulo observational uncertainties). To estimate upper limits, we show the nuclear star cluster (NSC) - galaxy scaling relations. NSCs and massive BHs often coexist in galaxies, with NSCs being the dominant central mass component in dwarfs 2 . Neumayer et al. (2020) compile NSC mass measurements  \nextending to low-mass dwarfs with stellar masses ≲ 10 6 M ⊙ , making the relation valid for the entire 𝑧 = 0 galaxy stellar mass range considered here.  \nKeeping these observational caveats in mind, we find that all of our dwarf AGN end up being overmassive compared to the extrapolated local BH - galaxy scaling relations. Most runs are offset by about an order of magnitude with the exception of m10zAGN which closely follows the early-type relation from Greene et al. (2020) and is only minimally offset from the Reines & Volonteri (2015) relation. Again, we emphasise that these offsets are not necessarily unexpected given that we are in the heavily extrapolated regime for these relations. Indeed the offsets are not as severe as the intrinsic high-redshift offsets inferred from recent JWST observations by Pacucci et al. (2023) and well within the bounds provided by the nuclear star cluster relation.  \nFor the redshift evolution, we plot the BH mass of the most massive progenitor for redshifts 𝑧 = 1 , 2 , 3 , 4. We note that in all cases the set-ups without CRs have more significant stellar and BH mass growth due to the weaker AGN feedback. The difference is especially significant at low redshifts. At late times, between 𝑧 = 1 and 𝑧 = 0, all simulations with CRs only have minimal stellar mass growth. For m10yCR+AGN and m10zCR+AGN , there is still significant BH growth during this time indicating that whilst the CR pressure has shut down star formation, the BH is still able to accrete. For m10qCR+AGN BH growth is already shut down from 𝑧 = 3, however, between 𝑧 = 4 and 𝑧 = 3, this run similarly experiences a rapid BH growth phase whilst star formation is already quenched.  \nWe also investigated the role of BH mergers versus accretion in growing the BHs in our dwarf simulations and find that in all cases BH growth is predominantly driven by gas accretion. In particular, we verified that merging with other BH seeds represents a negligible growth channel with a maximum of two, five and nine seed BHs forming in the m10q , m10y and m10z simulations, respectively. With a seed mass of 100 M ⊙ , this represents only a very small fraction of the final BH mass in all cases.", '3.4 Dark matter distributions': "Fig. 4 shows the spherically averaged dark matter density profiles as a function of radius. The three panels display these profiles for our three sets of ICs, respectively. For these profiles, we centre the dark matter distributions on the minimum potential of the main halo 3 . In each case, we plot the mean 𝑧 = 0 profiles, averaged over Δ 𝑇 = 500 Myr, for our four different galaxy formation model variations as colourcoded solid lines as indicated by the legend. We also performed additional dark-matter only (DMO) runs and we plot these as solid black lines. In each cases the temporal dispersion of these profiles is indicated by shaded regions which show the 1 𝜎 scatter in the densities over Δ 𝑇 = 500 Myr. We note that for the majority of the simulation runs, the central density is not very sensitive to the temporal bin width as the evolution of the central density is relatively steady over the last Gyr. The m10zCR simulation, however, experiences significant fluctuations in the central density from ∼ 12-13 Gyr due to mergers with substructures, so for our 𝑧 = 0 profiles we choose the bin width  \nmass, and BHs only start to dominate at stellar masses above a few times 10 10 M ⊙ .  \nFigure 4. Spherically averaged dark matter density profiles as a function of radius. The panels display simulated profiles based on our three sets of ICs, with colour-coding indicating the galaxy formation model variations as listed in the legend. The runs without CRs are shown as dotted lines for clarity. Density profiles from dark-matter-only (DMO) simulations are shown as solid black lines. The simulated profiles are averaged over Δ 𝑇 = 500 Myr and the 1 𝜎 temporal scatter is shown by the colour-coded shaded regions. The expected profiles for an NFW halo are plotted as grey-shaded regions indicating the 2 𝜎 scatter due to varying halo concentrations. The dark matter softening length is marked as a grey dashed line, and the Power radius (enclosing the central 2000 dark matter particles, see Power et al. 2003) is indicated by colour-coded dashed lines for each of the simulation set-ups. Flattening of the profiles within these radii may be partly due to numerical artefacts. CRs and/or AGN feedback cause significant diversity in dark matter profiles, ranging from cuspy to various degrees of cored profiles. Central densities at 𝑟 = 150 pc, a characteristic radius for distinguishing cusps from cores in observations (Read et al. 2019), are highlighted by marker symbols corresponding to the different galaxy formation physics configurations (with the circle symbol representing the DMO simulations).  \n<!-- image -->  \nto be small enough as to not be affected by these transient features. For all of the simulated profiles, we highlight the central dark matter density at 𝑟 = 150 pc, a characteristic radius for distinguishing cusps from cores in observations (Read et al. 2019), with colour-coded symbols.  \nFurthermore, we indicate the gravitational softening length at 𝑟 = 62 . 5 pc and the Power radius 4 (enclosing the nearest 2000 neighbours, see Power et al. 2003) as grey and colour-coded dashed vertical lines, respectively. Within these radii, the dark matter profiles may be impacted by numerical heating effects. In addition, we show the expected NFW profile based on the virial mass as grey shaded regions indicating the 2 𝜎 scatter from the halo mass - concentration relation (Dutton & Macciò 2014) to indicate the typical variations in profiles due to varying halo concentrations. The DMO runs are in good agreement with the NFW expectations, with m10y having slightly higher and m10z having slightly lower concentration than the mean relation.  \nFirst of all, we find that the introduction of AGN and CRs introduces a significant amount of scatter. The simulations with the default /f.pc/i.pc/r.pc/e.pc-3 model, without CRs or AGN, all lead to distinct cores at 𝑟 = 150 pc. The addition of CRs and/or AGN, however, leads to a variety of outcomes from stronger cores to weakened cores to cusps at the same inner radius 5 .  \nFor simulations without CRs, the addition of an AGN enhances the core size and further suppresses the central density for all three cases. Though it should be noted that the magnitude of the additional density suppression varies and is significantly weaker for m10y at 𝑟 = 150 pc.  \nFor simulations with CRs but without AGN, the profiles become cuspier than their no-CRs counterpart for m10q and m10y , whilst for m10z there is small suppression of the central densities with CRs. The addition of an AGN to these CR-based set-ups leads to a further increase of the central densities for m10q , whilst the central densities for m10y are minimally decreased. For m10z there is a much stronger impact, with the central density notable enhanced, yet still in the core regime.  \nOverall, there is a variety of outcomes when introducing CRs and/or AGN to dwarf simulations which mostly do not follow clear trends due to the dual role of these processes as an additional source of feedback that may both suppress star formation (and therefore suppress core formation) or enhance core formation via the additional energy input. Hence great care has to be taken to untangle the contribution from these additional baryonic processes and their interplay with star formation. A further consideration is that small perturbations in hydrodynamical simulations (e.g. floating-point round-off or random number generators) can lead to chaotic-like behaviour with the small perturbations growing over time and manifesting as macroscopic differences in galaxy properties (e.g. Genel et al. 2019; Keller et al. 2019; Borrow et al. 2023). However, the variability we observe here is significantly larger than the variability expected from this 'numerical butterfly effect' which typically leads to central density fluctuations of at most 10 to 15 per cent (see Keller et al. 2019) whereas we observe variations ranging from 70 to 90 per cent for our different sets of ICs and galaxy formation physics configurations (see Table 1). Furthermore, as we show in the next Section, these central density trends are relatively steady as a function of cosmic  \nFigure 5. Central dark matter densities as a function of cosmic time. The dark matter densities are binned over Δ 𝑇 = 500 Myr and the temporal 1 𝜎 scatter in each bin is shown by the shaded regions. The colour-coding indicates the respective simulation set-ups, as listed in the legend. For clarity, we plot the simulations without CRs as dotted lines. We measure the central densities at the inner radius 𝑟 = 500 pc to ensure that the central densities are well-resolved for all set-ups considered. For the m10q and m10y runs, the central densities are relatively steady at low redshifts, whilst the m10z runs still shows significant fluctuations due to a higher number of mergers and strong feedback events at late times.  \n<!-- image -->  \ntime which would not be expected if the variability resulted from the amplification of small numerical perturbations. Overall, this then points to a physical origin of the diversity of dark matter profiles that we observe in our simulations.", '3.5 Driving forces of core formation': 'Our analysis from Section 3.4 indicates that our different feedback physics configurations significantly affect the central dark matter distributions of dwarf galaxies. In this Section, we explore the driving forces behind these trends by contrasting the cosmic evolution of the central dark matter densities with SN and AGN energy injection histories.  \nFig. 5 shows the central dark matter densities at 𝑟 = 500 pc as a function of cosmic time. We choose this radius as it is small enough to showcase the differences between the models explored whilst being large enough to ensure we are not being dominated by numerical heating effects, i.e. we choose a radius that is at least the size of the Power radius for all set-ups explored. Again we average the dark matter densities over time bins of Δ 𝑇 = 500 Myr and indicate the mean densities as solid lines and the standard deviation in each bin by the shaded regions.  \nThe m10q set-up is shown in the first panel. Reflecting its early formation history, the dark matter density trends we found at 𝑧 = 0 in Fig. 4 have persisted for the past 6 Gyr. The central densities are relatively steady with only m10qAGN showing significant evolution after redshift 𝑧 = 1 with the central density of this run showing a slow decline.  \nSimilarly for m10y (middle panel), the 𝑧 = 0 trends have persisted for the past ∼ 5 Gyr. Both of the CR-based runs, m10yCR and m10yCR+AGN , do not show any significant evolution at low redshifts, whilst the runs without CRs, m10y and m10yAGN show a steady decline from 𝑧 ∼ 1 onwards. Note that the density fluctuations for these runs are also much stronger, reflecting the potential perturbations characteristic of core formation (e.g. Pontzen & Gov-  \n2; Martizzi et al. 2013). Focusing on early times between redshifts 𝑧 = 4 and 𝑧 = 1, we find that most of the runs are quite similar, except for m10yCR which exhibits a strong core during that period. We also note that at 𝑡 ∼ 5 Gyr, the m10y system experiences a major merger. For the DMO run this significantly raises the central density, and the strong early core for the m10yCR set-up is erased by this merger.  \nFor the m10z system (right panel), there is significantly more latetime evolution as well as stronger fluctuations as reflected in its more irregular morphology, see Fig. 1. Indeed, the general central density trends we observe at 𝑧 = 0 have only been in place for the past ∼ 2 Gyr and the m10zCR run experiences strong fluctuations for most of this period. We have verified that these potential fluctuations are driven by late-time mergers with substructures that pass through the centre and are independent of the slice width chosen or the method used to determine the halo centre when calculating the central density. Between redshifts 𝑧 = 2 and 𝑧 = 0 . 5, the set-ups m10zAGN , m10zCR , and m10zCR+AGN , all show a relatively similar redshift evolution of their central densities. The m10z set-up, however, closely follows the DMO run and even surpasses it from 𝑡 = 7-9 Gyr, indicating strong baryonic cooling at the centre. At 𝑧 ∼ 0 . 5 a merger induces density fluctuations for all of our m10z set-ups, but most prominently for the CR-based runs m10zCR and m10zCR+AGN . After this merger, the two CR-based set-ups have relatively steady central densities. On the other hand, the runs without CRs, m10z and m10zAGN experience significant density suppressions at late times.  \nIn the following two subsections, we turn to investigate the origin of these cosmic evolution patterns in the dark matter central densities, examining SN and AGN energy injection as a function of time.', '3.5.1 SN feedback': 'Fig. 6 shows the injected SN energy as a function of cosmic time with the instantaneous energy shown as dashed lines and the cumulative energy shown as solid lines, both binned over 250 Myr for clarity.  \nFigure 6. Instantaneous (dashed lines) and cumulative (solid lines) SN energy injected as function of cosmic time. The three columns display simulations based on our three sets of ICs, with colour-coding indicating the galaxy formation physics configurations as listed in the legends. We also indicate the 𝑧 = 0 binding energy of the halo as dotted lines. The cumulative SN energy injected compared to the binding energy is a strong predictor of core formation as well as the central dark matter density more generally. Some deviations from this pattern are introduced due to the important role of late-time star formation bursts in preserving cores.  \n<!-- image -->  \nThe three columns correspond to our three sets of ICs. The no-AGN and AGN runs are shown in the top and bottom row, respectively.  \nWe focus on the SN energy injected in the central region of the halo within a fixed 1 kpc aperture (in physical coordinates). This choice is motivated by the fiducial inner radius considered by observers for determining the shapes of dwarf galaxy rotation curves (see Section 3.7) and it also corresponds to the scale of the half-mass radius for most of our systems (see Fig. 2).  \nThe SN energy injection rates are recalculated from the stellar ages and metallicities based on S/t.pc/a.pc/r.pc/b.pc/u.pc/r.pc/s.pc/t.pc99 tables for a Kroupa (2001) IMF. We focus here on the SN II rates, which constitute the majority of SN events (see Hopkins et al. 2023a), for simplicity. For reference, we also indicate the 𝑧 = 0 binding energy of the respective haloes as dotted lines, which we estimate as 𝑓 bar 𝑀 vir 𝑉 2 vir , where 𝑓 bar is the universal baryon fraction and 𝑉 vir is the virial velocity of the halo (also see Wellons et al. 2023).  \nFirstly, we note that both CRs and AGN have a significant impact on star formation, and therefore the SN energy histories. As discussed in previous works (Sparre et al. 2017; Su et al. 2018; Genel et al. 2019; Keller et al. 2019; Hopkins et al. 2023c), due to their shallow  \npotential wells, dwarf systems are very sensitive to baryonic feedback processes so that even small changes can significantly affect star formation histories.  \nThis highly variable behaviour is quite clearly illustrated by the m10q set-up in the first column. With CRs, the system experiences a higher initial burst which is powerful enough to quench the dwarf system from 𝑧 ∼ 6. Hence both m10qCR and m10qCR+AGN end up with cuspy profiles, even though the cumulative injected SN energy for m10qCR+AGN slighly exceeds the 𝑧 = 0 binding energy, as the absence of late-time star formation means that a cusp can be reestablished quite easily after the initial burst. Note that m10qCR has a slighly reduced central density at 𝑧 = 0 compared to m10qCR+AGN due to low-level late-time star formation activity; however, these bursts are not powerful enough to induce a core. Without CRs, the initial star formation burst is less powerful, allowing for continued star formation activity. For the m10q set-up, there are a few more powerful bursts at intermediate redshifts, which then quench the systems at much later times from 𝑧 ∼ 1 . 5. For the m10qAGN setup, the AGN regulation means that the intermediate redshift star formation bursts are less intense so that star formation can continue  \nuntil 𝑧 = 0. Both m10q and m10qAGN end up with cumulative SN energy almost an order of magnitude above the 𝑧 = 0 binding energy leading to cored profiles; however, m10qAGN has a more extreme core due to the continued star formation activity (also see Muni et al. 2024).  \nFor m10y , there is a similar picture albeit not as extreme as for the m10q halo. With CRs, high-redshift star formation is somewhat burstier, and the more powerful bursts then also mean that most of the SN activity is restricted to 𝑧 ≳ 1. This is still sufficient to (cumulatively) exceed the 𝑧 = 0 binding energy and induce a core for both of the CR-based runs; however, the cores are significantly weaker and less extended than for the runs without CRs. Both m10y and m10yAGN have continued star formation activity until 𝑧 = 0, higher cumulative injected SN energies, and therefore much stronger cores. Adding in AGN only has a small impact for both m10yAGN and m10yCR+AGN compared to their no-AGN counterparts. Again the bursts at high redshift are somewhat weaker, allowing for more SNactivity at lower redshifts and therefore somewhat stronger cores. However, as discussed in Sections 3.4 and 3.6, the differences induced by AGN feedback are only minimal for the m10y profiles at 𝑧 = 0.  \nFor m10z , we obtain the most complex interplay between star formation, AGN activity and CRs. Here, the set-ups without CRs have more powerful initial star formation bursts (underlining that the initial intensity of the burst is not necessarily directly dictated by the galaxy formation physics configuration). Overall, m10z , m10zCR , and m10zAGN all experience episodes of star formation for most of cosmic time, although for m10z , those are shifted to lower redshifts, whilst for the m10zCR simulation SN activity is more concentrated at higher redshifts. m10zAGN experiences bursty star formation throughout cosmic time, leading to the highest cumulative injected SN energy and the most pronounced core. Likewise, the cumulative SN energy in the m10z and m10zCR set-ups also exceeds the 𝑧 = 0 binding energy and both of these result in cored profiles. Interestingly, however, star formation in the m10zCR+AGN set-up is significantly suppressed, and powerful bursts are mostly restricted to high redshifts ( 𝑧 > 2), so that the total injected SN energy remains well below the 𝑧 = 0 binding energy by a factor of three. This indicates that for the m10zCR+AGN set-up, the core formation may be driven by the AGN.  \nThis analysis underlines that CRs and AGN primarily affect cuspto-core transformation, and the diversity of dark matter profiles, by regulating star formation histories in our simulations. The only setup where the AGN feedback is actively driving core formation is m10zCR+AGN , where star formation is too heavily suppressed to drive core formation. In the cases of m10q and m10y , the AGN is not efficient at suppressing star formation, leading to similar or higher levels of SN activity compared to the equivalent no-AGN setups. However, in these scenarios, the AGN may still enhance core formation; we examine this possibility in the next section.', '3.5.2 AGN feedback': "In the /f.pc/i.pc/r.pc/e.pc-3 model, AGN feedback is injected in several channels including mechanical feedback, radiative feedback and CRs (see Section 2.2.1 for details). Comparing this multi-channel AGN energy to the SN energy that is predominantly injected mechanically is not straightforward, especially since we are most concerned with feedback-induced gravitational potential perturbations. Therefore, we need to assess how the AGN energy will couple dynamically. Fully tracking and analysing coupling efficiencies of different AGN feedback modes is beyond the scope of this paper, so we make some simplifying assumptions based on previous literature.  \nFor the AGN, only a very small fraction of the injected energy is directly coupled in the form of mechanical winds, with only ∼ 0 . 6 per cent of the AGN luminosity injected as wind energy (see Section 2.2.1).  \nFor the radiative feedback, on the other hand, the full luminosity is injected into the simulation domain using the LEBRON method (see Hopkins et al. 2020a) for radiation transport. This radiation transport then predicts the radiation pressure forces, i.e. momentum flux, from the absorbed photon luminosity. Overall only a small fraction of the luminosity will therefore be energetically coupled to the gas as a radiatively driven wind. The LEBRON method has been most thoroughly tested for stellar feedback. Whilst AGN radiation has a different spectrum, we can still draw some tentative conclusions from these studies for our estimates. For the stellar feedback, only 0 . 4 𝐿 / c of the photon momentum couples to the gas in dwarfs, whilst the remainder escapes the galaxy, and the vast majority of this photon coupling occurs in the single-scattering regime (see analysis in Hopkins et al. 2020a). We here assume that a similar fraction of the luminosity would couple to the gas for radiative AGN feedback. Following (King & Pounds 2015), a fraction of approximately 5 per cent of that availabile AGN luminosity is then coupled as mechanical energy in the form of a radiatively driven wind in the single scattering regime, yielding an overall mechanical energy efficiency of 2 per cent compared to the AGN luminosity for the radiative feedback channel.  \nFor the AGN-driven CR feedback, 1 per cent of the accreted restmass energy, i.e. 10 per cent of the luminosity with the assumed radiative efficiency of 10 per cent, is injected as CRs. For CR transport, we employ the subgrid model from Hopkins et al. (2023b), which interpolates between the limit in which CRs escape the galaxies with negligible losses and that in which CRs lose most of their energy catastrophically before escaping, using a formalism akin to the LEBRON method. For dwarf galaxies, we are generally far from the proton calorimetric limit where most energy is lost before CRs escape the dense gas, and we assume that (as with the radiation) 40 per cent of the CRs couple before escaping - which likely provides an upper limit. This then yields an overall mechanical energy efficiency of 4 per cent compared to the AGN luminosity for the CR feedback channel.  \nKeeping in mind that these assumed efficiencies may be overestimating the actual coupled AGN energy, we show the 'effective' injected AGN energy based on our assumed efficiencies for the different channels as a function of cosmic time in Fig. 7. Again the cumulative energy is shown as the solid lines and the instantaneous energy is shown as dashed lines binned over 250 Myr for clarity. The AGN energy is based on the accretion rate of the most massive BH in the main halo at 𝑧 = 0 and its most massive progenitor halo at higher redshifts. The three different panels represent our three different ICs and the colour coding of the lines indicates the different physics configurations as listed in the legends.  \nOverall, we note that AGN energy injection is more steady than SN energy injection, with the AGN in most set-ups being active continuously, except for m10qCR+AGN where the AGN activity, just like the star formation, is quenched at early times and is then mostly quiescent at late times apart from a relatively weak late-time burst. In all cases, the injected AGN energy matches or exceeds the binding energy though we caution that due to the extremely uncertain mechanical coupling efficiencies this has to be taken with a grain of salt and the AGN energy that is actually coupled mechanically to the gas may be lower. Furthermore, we note that whilst the SN feedback is distributed throughout the galaxy the AGN energy injection is only focussed at the centre, with radiation and mechanical winds injected  \nFigure 7. Instantaneous (dashed lines) and cumulative (solid lines) effective AGN energy injected as a function of cosmic time. The three panels display simulations based on our three sets of ICs, with colour-coding indicating the galaxy formation physics configurations as listed in the legends. The 𝑧 = 0 binding energy of the halos is indicated by colour-coded dotted lines. AGN feedback is injected through three channels: radiation, mechanical winds, and, for some, CRs. The overall 'effective' energy is calculated based on reported couplings from the literature. Runs with CRs show more efficient AGN feedback regulation, with m10zCR+AGN experiencing efficient bursts that induce a core despite low star formation activity. Without CRs, BHs grow more efficiently, leading to higher energy injection rates and increased core sizes for all set-ups.  \n<!-- image -->  \nin a collimated fashion, and therefore AGN feedback is likely less efficient than SN feedback in affecting the ISM at a fixed energy injection rate. Indeed, past simulations have demonstrated that a significant fraction of the AGN wind energy escapes the ISM once the winds have opened up a central cavity (e.g. Bourne et al. 2015; Koudmaniet al. 2019; Torrey et al. 2020; Mercedes-Feliz et al. 2023).  \nFor m10q , in the CR-based run m10qCR+AGN , the AGN does not have a significant impact on the central densities since the late-time activity, which is crucial for core formation, is mostly suppressed. Without CRs, for the m10qAGN set-up, there is a sharp AGN burst at 𝑡 ∼ 8 Gyr which is also associated with pronounced central dark matter density fluctuations and followed by a steady decline in the central dark matter density suggesting that the AGN may be contributing to the formation of the core with this set-up.  \nFor m10y , the m10yAGN set-up has a strong AGN burst at 𝑡 ∼ 5 . 5 Gyr, which is again associated with strong central density fluctuations indicating gravitational perturbations as well as a drop in the star formation rate. With the CRs, the AGN energy injection is again significantly suppressed (though not completely shut down like in the m10qCR+AGN set-up) and the relatively low levels of activity do not have a significant impact on core formation.  \nFor m10z , contrary to the other two sets of ICs, the energy injected from the CR-based run, m10zCR+AGN , exceeds the AGN energy injected with the equivalent no-CR set-up m10zAGN . This is mostly driven by a significant burst in the m10zCR+AGN set-up at 𝑡 ∼ 9 Gyr, which takes the overall injected AGN energy budget above the binding energy of the halo and is followed by strong density fluctuations. Most notably, even though central star formation is completely suppressed from 𝑡 ∼ 10 Gyr in the m10zCR+AGN simulation, the AGN activity persists at relatively high levels, maintaining the core in the absence of SN-driven winds. This is driven by two factors in our AGN modelling. Firstly, at late times the ratio between BH mass and subgrid accretion disc mass in the m10zAGN run is quite large ( 𝑀 BH / 𝑀 d ∼ 10 4 ), which leads to depletion timescales of the order of ∼ 1 Gyr (also see Section 2.2.1), allowing for persistent AGN activity even after the inflows onto the BH - accretion disc particle have sub-  \nlst the gas density is significantly decreased in the central region (and therefore the gas is not star-forming), the BH is still able to accrete at low rates from this supply since there is no explicit density criterion for the BH accretion. This likely represents an optimistic estimate of the BH accretion rate and we discuss the implications and caveats of these modelling choices in more detail in Section 4.3.2. For the m10zAGN set-up the AGN has a significant burst at 𝑡 ∼ 12 Gyr which again is associated with strong central density fluctuations suggesting that here the AGN may also enhance the formation of the core.  \nOverall, we note for most of our AGN set-ups, including m10qAGN , m10yAGN , m10yCR+AGN , m10zAGN , and m10zCR+AGN , the equivalent no-AGN set-up already has a SNinduced core. For all of these cases, apart from m10zCR+AGN , the addition of an AGN further decreases the central densities enhancing the existing core. However, it is difficult to disentangle whether this additional decrease is driven directly by AGN energy injection or indirectly by suppressing high-redshift star formation and 'delaying' SNactivity to later times favouring core formation. Indeed, it is likely that both of these factors contribute to the enhancement of cores.  \nThe m10zCR+AGN set-up is our only simulation where the AGN is efficient in globally suppressing star formation compared to its no-AGN counterpart (also see Fig. 2). Here, there is no SN activity for the last ∼ 3 . 5 Gyr from 𝑧 ∼ 0 . 5 and only weak SN activity between 𝑧 = 2 and 𝑧 = 0 . 5, yet the core is maintained. Indeed, the central density significantly decreases from 𝑧 = 2 to 𝑧 = 1 and is then maintained at low levels despite the late-time merger at 𝑧 ∼ 0 . 5, providing strong evidence for AGN-driven core formation.  \nOverall, we conclude that AGN feedback can induce significant scatter in dark matter density profiles, predominantly indirectly by regulating star formation histories but also in some cases directly by driving powerful potential fluctuations.  \nFigure 8. Central dark matter densities at 𝑟 = 150 pc as a function of virial mass ( left panel ), stellar mass to virial mass ratio ( middle panel ), and BH mass to virial mass ratio ( right panel ). Simulated 𝑧 = 0 data points are highlighted by marker symbols corresponding to the respective galaxy formation physics configurations as indicated by the legend. The binned redshift evolution from 𝑧 = 3 to 𝑧 = 0 is represented by dotted lines. The left and middle panels include various observed central dark matter densities from Read et al. (2019) for comparison. In the middle panel, a dashed grey line marks the stellar mass to virial mass ratio where previous studies suggest stellar-feedback-driven core formation may become inefficient. Our simulations broadly agree with this theory, except for m10zCR+AGN , which forms a core despite a decreasing stellar mass to virial mass ratio, though has an increasing BH mass to virial mass ratio, indicating AGN-driven core formation.  \n<!-- image -->", '3.6 Cosmic evolution of central dark matter densities in observational context': 'We further analyse the origin of the diversity of our simulated dark matter profiles by examining the central densities in the context of observations. Fig. 8 shows the central dark matter densities at 𝑟 = 150pc as a function of virial mass (left panel), stellar mass to virial mass ratio (middle panel) and BH mass to virial mass ratio (right panel). The dotted lines indicate the (binned) redshift evolution of our simulated dwarfs from 𝑧 = 3 to 𝑧 = 0, again binned over 500 Myr. For clarity we only highlight the 𝑧 = 0 data point with the marker symbol corresponding to our different physics configurations, as indicated by the legend. In the left panel we also show the expected central densities for an NFW profile and a cored profile following the /c.pc/o.pc/r.pc/e.pcNFW profile from Read et al. (2016) as grey lines with the 2 𝜎 concentration scatter indicated by dark-grey and light-grey shaded regions, respectively. Furthermore, we show the observed central densities of 16 observed nearby dwarf galaxies, including classical dwarf spheroidals and gas-rich dwarf irregulars from Read et al. (2019).  \nThe first panel demonstrates that our simulations are generally in good agreement with the observed central densities of nearby dwarfs, keeping in mind that the sample sizes are limited in both cases (16 observed dwarfs versus 12 simulated dwarfs). Two of the m10z simulations, m10zCR and m10zAGN , lead to very strong cores, with the central densities being more severely suppressed than in the observations. Our simulations do not reproduce the observed strong cusps with central densities above the NFW mean relation. However, we note that the vast majority of our dwarfs have late star formation histories, with only the CR set-ups for m10q and m10y resulting in an early truncation of star formation, as defined by Read et al. (2019), with no activity for the past 6 Gyr. These dwarfs with early truncation are associated with cuspy profiles hinting that we would require more early-forming set-ups to reproduce these strong cores.  \nThe middle panel shows the simulated and observed central densities as a function of stellar mass to virial mass ratio 𝑀 stellar / 𝑀 vir . As discussed in Read et al. (2019), the observed dwarfs show a clear trend with higher 𝑀 stellar / 𝑀 vir being associated with lower central dark matter densities. Moreover, in the observations, as well as in previous theoretical studies, it is found that no stellar-feedback-driven  \ncores form below a critical ratio of ( 𝑀 stellar / 𝑀 vir ) crit = 5 × 10 -4 since the integrated SN energy is insufficient to drive the required potential fluctuations (see discussion in Di Cintio et al. 2014). We also find a strong correlation between 𝑀 stellar / 𝑀 vir and central density suppression for our simulated dwarfs. Indeed, the correlation is significantly tighter than for the observed dwarfs likely due to increased scatter induced by measurement uncertainties in the observations. Looking at the redshift evolution tracks of our simulated dwarfs, we find that most systems also follow these trends in a temporal sense with two notable exceptions. The CR-based runs with the m10q ICs show a slight decrease as the 𝑀 stellar / 𝑀 vir decreases. However, this is mostly a by-product of these runs being quenched at very high redshifts so that the central dark matter density evolution is dominated by the general assembly history of the halo with the virial mass increasing and the halo concentration decreasing (whilst the stellar mass remains constant) at late times. The second exception to the 𝑀 stellar / 𝑀 vir trend, however, is more notable: for the m10zCR+AGN set-up the central densities as a function of 𝑀 stellar / 𝑀 vir are steadily decreasing so that a core forms despite the 𝑧 = 0 𝑀 stellar / 𝑀 vir ratio being below the critical ratio identified by Di Cintio et al. (2014) and Read et al. (2019). This provides further evidence that this run could be experiencing AGN-driven core formation.  \nWe further examine this possibility in the third panel where we plot the simulated central densities (as well as their binned redshift evolution from 𝑧 = 3) for the AGN-based simulation set-ups as a function of BH mass to virial mass ratio 𝑀 BH / 𝑀 vir . For the runs without CRs, m10qAGN , m10yAGN , and m10zAGN , the central densities again steadily decrease as a function of 𝑀 BH / 𝑀 vir . However, we caution that this does not necessarily indicate that the BHs are contributing to core formation due to the strong correlation between 𝑀 BH and 𝑀 stellar (see Fig. 3 for scaling relations of our simulated dwarfs). The m10qCR+AGN set-up is the only simulation where the central density ratio decreases as 𝑀 BH / 𝑀 vir decreases, as with the stellar mass to halo mass ratio, and this mainly reflects that the central density is unaffected by baryonic processes at late times. The m10yCR+AGN set-up has a very shallow gradient again reflective of the relatively quiescent late-time evolution. Finally, the m10zCR+AGN set-up exhibits a significant decrease in central densities as the 𝑀 BH / 𝑀 vir increases (opposite trend from 𝑀 stellar / 𝑀 vir ) providing further support for the BH-driven core formation scenario.  \nFigure 9. Circular velocity profiles at 𝑧 = 0 as a function of radius including all matter (solid lines) and baryons only (dashed-dotted lines). The profiles are averaged over Δ 𝑇 = 500 Myr with the shaded regions indicating the 1 𝜎 scatter. The first three columns display simulations based on our three sets of ICs with the no-AGN runs (top row) and AGN runs (bottom row) plotted separately for clarity. The markers highlight the circular velocity at the fiducial radius which is used as a reference point in observational studies of rotation curves. The fourth column translates these velocity curves to the rotation curve shape - baryon dominance parameter space, also showing observed dwarf galaxies in a similar mass range from the SPARC sample as compiled in Santos-Santos et al. (2020). Cosmic environment, CRs and BH feedback can all influence the rotation curve shape parameter, with our simulations broadly reproducing the observed trends with baryonic dominance.  \n<!-- image -->', '3.7 The diversity of rotation curves': "As discussed in Section 1, one of the most prominent remaining dwarf galaxy 'problems' pertains to the observed diversity of dwarf galaxy rotation curves. From our analysis in the previous sections, we have found that AGN and CRs may significantly increase the scatter in dark matter density profiles and hence may contribute to resolving this on-going controversy. Ultimately, to make statistical predictions for dwarf galaxy rotation curve shapes, we would need a much larger simulation sample spanning a larger range of environments and halo masses. Nevertheless, we can use our simulation suite to examine the broad trends for rotation curve shapes and assess the potential for more sophisticated galaxy formation models including the impact of AGN and CRs at the low-mass end to resolve the long-standing diversity of rotation curves problem.  \nIn Fig. 9, we present circular velocity curves for all our dwarf simulations. As for the dark matter profiles, we average the velocity curves over Δ 𝑇 = 500Myrandindicatethe1 𝜎 scatter with the shaded regions. Following the methodology presented in Santos-Santos et al. (2020), we plot circular velocity profiles for the total matter content of our dwarfs (dark matter, gas, stars, and BHs) as solid lines and for baryons (gas, stars, and BHs) as dotted-dashed lines.  \nThe total circular velocity curves for each set of ICs converge between 1 to 10 kpc, with the extent of the matter deficits in the central regions reflecting the core sizes from Fig. 4. For the baryonic circular velocity curves, however, there are significant differences. We also examine the gas and stellar distributions (not shown here) and find that both of these are significantly modified for galaxies with pronounced cores, which tend to have more extended stellar distributions (also see Figure 1). For the cored dwarfs this then leads to slowly rising rotation curves for both dark matter and stars. Dwarfs that are quiescent at 𝑧 = 0, and generally have cuspier profiles, are more gas poor (especially at the centre), leading to overall lower baryonic circular velocities in the central region.  \nFollowing Santos-Santos et al. (2020), we quantify the cusp versus core nature of the profiles by calculating characteristic velocity ratios, comparing the inner 'fiducial' rotation velocity 𝑉 fid with the asymptotic flat rotation velocity 𝑉 max . The fiducial rotation velocity  \nis measured at the fiducial inner radius defined as:  \n𝑟 fid = GLYPH<18> 𝑉 max 35 km s -1 GLYPH<19> kpc . (4)  \nThe rotation curve shape parameter 𝜂 rot can then be defined as  \n𝜂 rot = 𝑉 fid 𝑉 max = 𝑉 ( 𝑟 fid ) 𝑉 max . (5)  \nRapidly rising rotation curves have 𝜂 rot ≲ 1, with the NFW profile yielding a rotation curve shape parameter of 𝜂 rot ∼ 0 . 65. Cored profiles have 𝜂 rot ≪ 1; see Santos-Santos et al. (2020) for a detailed discussion of this parameter.  \nAsin Santos-Santos et al. (2020), we also inspect a second velocity ratio - the baryonic importance parameter, which is defined as  \n𝜂 bar = GLYPH<18> 𝑉 b , fid 𝑉 fid GLYPH<19> 2 = GLYPH<18> 𝑉 b ( 𝑟 fid ) 𝑉 ( 𝑟 fid ) GLYPH<19> 2 , (6)  \nwhere 𝑉 b , fid is the contribution from baryons (gas, stars, BHs) to the circular velocity curve. Under the assumption of spherical symmetry, 𝜂 bar represents the mass fraction of baryons within the fiducial radius.  \nThe two velocity ratios 𝜂 rot and 𝜂 bar can be used to plot the 'cuspiness' of the rotation curves against the contribution from baryons in the inner region of the galaxy. We plot this distribution in the rightmost panel in Fig. 9. For comparison, we also show observed data as collated by Santos-Santos et al. (2020). We restrict the observed galaxy sample to the classical dwarf regime with 𝑉 max ≤ 60 km s -1 to compare this more easily and directly with our simulated dwarfs. The observational data stems mainly from the SPARC sample and, in all cases, the rotation curves were inferred from high-resolution HI and/or H 𝛼 velocity fields. Given that all of the observational rotation curves are inferred from the baryonic contributions, the dwarfs with extremely low baryonic contributions are missing from the observed samples. Similarly to Santos-Santos et al. (2020) we therefore show our low𝜂 bar simulated dwarfs at 𝜂 bar = 0 . 05 for clarity.  \nFocussing first on the observations, we find that, as discussed in Santos-Santos et al. (2020), the baryonic importance parameter appears to be a poor predictor of rotation curve shapes with a significant amount of scatter in the data. Specifically, this plot highlights  \nthe appearance of profiles ranging from cusps to weak cores at low baryonic dominance and cusps to strong cores at high baryonic dominance in the observations. This large amount of scatter may naively seem at odds with the theory that baryons are driving cusp to core transformations and has been extremely difficult to reproduce with simulations. Despite the high scatter, there are nevertheless some interesting observational trends where strong cores in dwarfs tend to be associated with higher baryonic importance and strong cusps tend to be associated with lower baryonic importance. However, both of these aspects seem counter-intuitive if cores are associated with strong outflows and cusps with baryons cooling at the centre.  \nOur simulations cover almost all observed scenarios: cusps and weak cores at low baryonic dominance (e.g., m10qCR and m10zCR+AGN ) and cusps and strong cores at high baryonic dominance (e.g., m10y and m10zAGN ). The only missing scenario are the 'ultra cuspy' dwarfs at both low and high baryonic dominance, which might be influenced by environmental factors, suggesting our sample might be too small to cover the required range of assembly histories. To reproduce these extremely cuspy dwarfs with CDM-based simulations, we would need strong cosmic inflows so that feedback becomes inefficient leading to overcooling (see e.g. Smith et al. 2019). Another possible pathway for inferring ultra-cuspy dwarfs from gas rotation curves constitutes strong, transient perturbation to the gas, usually resulting in an offset between the gas, stars and dark matter. However, in the fiducial /f.pc/i.pc/r.pc/e.pc-3 model (as well as in the SPARC sample), these transient ultra-cuspy dwarfs are generally associated with more massive systems ( 𝑀 halo ≳ 10 11 M ⊙ and 𝑉 max ≳ 50 kms -1 ), see Sands et al. (2024) for details.  \nNevertheless, even with our small sample, we obtain a range of profiles at fixed baryonic importance ranging from cusps to cores, a feat that has been historically difficult to reproduce in simulations (e.g. Oman et al. 2015; Garrison-Kimmel et al. 2019; Santos-Santos et al. 2020). We are able to reproduce these trends due to the interplay of the different feedback processes resulting in varying overall feedback efficiencies for the same baryon content, either enhancing or reducing the overall feedback impact compared to simulations that do not include CRs or BHs. Cusps are generally associated with highly concentrated stellar distributions in gas-poor dwarfs, whilst cores are created by galactic fountain activity in gas-rich dwarfs that produces a core through cyclic feedback events but does not expel the baryons. In this manner, we are able to match both the large observed scatter and the weak observed trends with baryonic dominance.", '4.1 Comparison with previous theoretical work': "We have investigated the impact of BH feedback and CRs on the cusp versus core problem in dwarf galaxies with a special focus on whether including additional baryonic feedback channels may help explain the observed diversity of rotation curves in dwarf galaxies.  \nSeveral groups have investigated the role of stellar feedback in driving cusp to core transformations in galaxy formation simulations with mixed results (e.g. Governato et al. 2010; Parry et al. 2012; Pontzen & Governato 2012; Hopkins et al. 2014, 2018; Kimm et al. 2015; Emerick et al. 2018; Smith et al. 2019; Gutcke et al. 2021). Most simulations tend to either predominantly produce cusps or predominantly produce cores which appears to be closely related to the star formation threshold employed and, more crucially, whether the cold, dense ISM phase is resolved in the simulations (see discussion in Jahn et al. 2023). Nevertheless, for a given ISM model it has proven  \nvery challenging for simulations to reproduce the observed diversity of dwarf rotation curves (e.g. Oman et al. 2015; Garrison-Kimmel et al. 2019; Santos-Santos et al. 2020). The past decade has unveiled growing observational samples of dwarf galaxies with AGN. These observations raise the question of whether active BHs may increase the diversity in star formation histories as well as directly impact central dark matter distributions and thereby increase the diversity of dwarf rotation curves in simulations.  \nPrevious studies found systematic yet small-scale suppressions of central dark matter densities by AGN feedback in dwarfs (Koudmani et al. 2022; Arjona-Galvez et al. 2024), yet in both cases, the simulations did not explicitly resolve the multiphase ISM. Cusp-to-core transformations with AGN feedback were also explored by Arora et al. (2024). However, they employed the fiducial Bondi accretion model for BH growth, which suppresses the growth of low-mass BHs such that AGN activity in dwarfs is insignificant in these simulations and does not affect dark matter distributions.  \nIn our study, we investigate, for the first time, the impact of AGN feedback on dark matter profiles in dwarfs with a resolved multiphase ISM model and a BH growth scheme that does not suppress AGNactivity in dwarf galaxies. Several works exploring BH physics within the /f.pc/i.pc/r.pc/e.pc model have found that CRs appear to be a key ingredient for effective AGN feedback (e.g. Su et al. 2021; Wellons et al. 2023; Byrne et al. 2023b), so we also investigate the impact of CR feedback on the interplay between baryons and dark matter in dwarfs. We find that efficient AGN feedback can actively drive core formation and, together with CRs, significantly enhances the diversity of dwarf rotation curves by leading to more varied star formation histories. Indeed, in many cases our runs including CRs (with or without AGN) lead to cuspier profiles due to suppressed star formation activity. This is in contrast with recent findings from Martin-Alvarez et al. (2023), who find that CRs are prone to enhance core formation. We note that for one of our set-ups (compare m10z to m10zCR ), we also observe this behaviour, and indeed it is plausible that CRs could act to either suppress or enhance cores depending on their interaction with the ISM and the assembly history of the given dwarf galaxy. Our most important finding is that BHs and CRs as additional baryonic processes lead to varying overall feedback efficiencies that can lead to both cuspier or more cored profiles for a given assembly history.  \nWe also demonstrate the promising potential for these additional baryonic processes to resolve the diversity of dwarf galaxy rotation curves problem: our simulations exhibit a wide variety of profiles, including cuspy profiles at low baryonic dominance and strong cores at high baryonic dominance, thereby reproducing the puzzling relationship between rotation curve shapes and the gravitational importance of baryons from the observations. The former is associated with highly concentrated stellar distributions in gas-poor dwarfs, whilst the latter is a signature of galactic fountain activity in gas-rich dwarfs that produces a core through cyclic feedback events but does not expel the baryons.  \nOne remaining discrepancy with the observations lies in the 'extremely' cuspy dwarfs which are even denser than NFW expectations. Within our relatively modest sample size of simulations we are unable to reproduce these. Apart from the obvious need to explore more halo masses and environments, as these objects may be linked to strong cosmic inflows which lead to overcooling (see e.g. Smith et al. 2019), in future work, it would also be important to explore observational uncertainties in rotation curve measurements as well as theoretical uncertainties in the modelling of BHs and CRs, and we discuss both of these aspects in the following sections.", '4.2 Observational uncertainties': "In observations, the dark matter distribution of dwarf galaxies may be inferred from the rotation curves of gas or stars, yet it has been shown that in practice these rotation curves may be highly unreliable tracers of the 'true' underlying circular velocity distribution. In particular, various dynamical perturbations may lead to large discrepancies of 50 per cent or more (see Pineda et al. 2017; Roper et al. 2023; Downing & Oman 2023; Sands et al. 2024, for some recent studies). Processes that are associated with inducing AGN activity, including mergers and strong cosmic gas inflows, lead to inaccurate estimates of the matter distributions in dwarfs. Interestingly, all of these studies find that observational uncertainties are more likely to lead to an underestimation of the true circular velocity, i.e. overestimation of the occurrence of cored profiles in observations, whilst the inverse error is less common though may still occur in the central few kiloparsecs due to dynamical phenomena, especially for more massive dwarfs (see discussion in Sands et al. 2024).  \nApart from the observational uncertainties in inferring dark matter mass distributions, it also remains extremely challenging to determine the drivers behind core formation from observations. In particular, there are no observational constraints on the potential role of AGNincusp-to-core transformations. This is mainly due to the number of observed AGN in nearby dwarf galaxies being very limited. To further complicate matters, AGN are associated with disturbed gas rotation curves - see e.g. Manzano-King & Canalizo (2020), who find a strong association between AGN activity and disturbed gas kinematics in observed dwarf rotation curves, making it even more difficult to establish a link between AGN activity and central dark matter densities. JWST may partly alleviate these issues by disentangling rotational and outflow components with high-resolution NIRSpec-IFU observations (e.g. Bohn et al. 2023). From our simulations, we predict that if cored dwarf galaxies were found below the critical mass ratio of GLYPH<16> 𝑀 stellar 𝑀 vir GLYPH<17> crit = 5 × 10 4 this would be a strong indicator of AGN-driven core formation in dwarfs.", '4.3 Theoretical uncertainties': "The diversity of models presented in this work (with or without CRs and with or without AGN) are meant to represent different outcomes of 'the same underlying physics' under different conditions. In particular, the impact of AGN feedback is expected (and observed) to be highly variable in the dwarf regime. The AGN activity levels in dwarfs are highly dependent on the seeding model assumed (with some dwarfs possibly not hosting any BHs at all) and whether the BH is off-centre for the majority of the dwarf's history, i.e. the no-AGN model here may represent the same underlying physics as the AGN model but in a situation where no BH seeded or only an extremely weakly accreting off-centre BH is present. Hence just seeding and BH dynamics could drive substantial variability in addition to other effects such as varying feedback efficiencies depending on BH spin evolution. Similarly, the CR subgrid model we employ here represents a case where losses within the galaxy are minimal and hence the impact of CRs in our dwarfs likely represents an upper limit. With full CR transport the outcome will likely lie somewhere inbetween the no-CR and CR simulations presented in this work. It would require significantly higher resolution simulations (and correspondingly more detailed BH modelling) and explicit CR transport to depict the full variability in dark matter distributions introduced by baryonic physics within one galaxy formation model. Our simulations provide strong motivation for exploring such detailed models  \nin future work. Below we discuss the caveats of our CR and BH modelling in more detail.", '4.3.1 CR modelling': 'Whilst the focus of our paper lies on the interplay between AGN feedback and dark matter distributions in dwarfs, we also assess the impact of CRs on the cusp versus core problem since CR injection has been identified as an important AGN feedback channel. For this, we take advantage of the CR subgrid model by Hopkins et al. (2023b), which allows for the inclusion of the impact of CRs on galaxy formation without the significant computational overhead associated with full CR transport.  \nAs discussed in previous sections, for dwarf galaxies, this subgrid model presents a good approximation since, in this mass regime, we are far away from the calorimetric limit for protons. Consequently, we can safely assume that losses within the galaxy only play a secondary role. Nevertheless, with this subgrid model, we may still be somewhat overestimating the effects of CR feedback at late cosmic times, as we are not fully capturing the inhomogeneity of CRs in the CGM, which may modify thermal instabilities (Butsky et al. 2020). In particular, we find that this implementation of CR feedback suppresses star formation so strongly that the addition of CRs generally leads to cuspier profiles. Indeed Hopkins et al. (2023b) also find in their validation simulations with a 10 11 M ⊙ halo that the subgrid CR model leads to suppressed star formation and a somewhat weaker core compared to an equivalent simulation with full CR transport, though these differences are small ( ∼ 1 per cent) and both CR models lead to a weaker core compared to their no-CR set-up. Our results should therefore be regarded as an upper limit for the impact of CR feedback and, whilst exploring alternative CR implementations is beyond the scope of this paper, our work provides strong motivation for investigating cusp-to-core transformations with AGN and full CR transport in future work.', '4.3.2 BH modelling': "The modelling of AGN feedback in cosmological simulations remains subject to significant theoretical uncertainties due to the vast dynamic range of relevant scales spanning (at least) 14 orders of magnitude from the event horizon ( ∼ 10 -6 pc for Sgr A*) to the cosmic web ( ∼ 10 8 pc). Hence BH growth, feedback, and dynamics cannot be modelled ab-initio and need to be included as 'subgrid' processes in cosmological simulations.  \nHere, we model the BH growth based on the torque-based accretion model (Hopkins & Quataert 2011), which allows for efficient BH growth in the dwarf regime because it does not suffer from the strong BH mass dependence inherent to the widely used Bondi model (also see Anglés-Alcázar et al. 2013, 2015; Koudmani et al. 2022; Wellons et al. 2023). This leads to relatively high BH masses compared to the extrapolated scaling relations - although we emphasise again that this comparison should be interpreted with caution because we would expect a flattening of the scaling relations in this mass regime (Greene et al. 2020). Nevertheless, from an accretion perspective, our simulations may be exploring an upper limit on the impact of BHs in the dwarf regime.  \nThe inflow rates provided by the torque-based model are then coupled to a gas reservoir representing the accretion disc, with the depletion time of this reservoir set by the thin 𝛼 -disc model. However, there are a few caveats to this approach. Firstly, the 𝛼 -disc is generally only applicable in the radiatively efficient accretion regime.  \nSecondly, this approach only tracks the mass flow rates through the disc and does not follow angular momentum flow rates. This in turn means that the BH spin evolution cannot be tracked self-consistently; therefore we have to assume constant efficiencies for the disc and jet luminosities, which otherwise could be directly inferred from the disc and BH properties (see e.g. Beckmann et al. 2019; Sala et al. 2021; Talbot et al. 2021, 2022, 2024; Huško et al. 2022, 2024; Koudmani et al. 2024).  \nFor the AGN feedback, we only explore one set of parameters for the three AGN channels (mechanical winds, radiation and CRs) that allows us to match the main observed characteristics of nearby dwarf galaxies; see Fig. 2. However, as discussed in Wellons et al. (2023), there are several other 'plausible' AGN feedback set-ups within the /f.pc/i.pc/r.pc/e.pc model, leaving a large parameter space to be explored. In particular, we assume collimated injection for winds and radiation; however, the cusp versus core problem may be sensitive to alternative injection geometries, as discussed in Zhang et al. (2024) for the case of stellar feedback. This may allow us to also have a direct impact on dwarf core formation at lower halo masses, as in the m10zCR+AGN set-up, without compromising in terms of producing realistic stellar properties.  \nThe final caveat in our model is that we do not have sufficient resolution to self-consistently model BH dynamical friction. Hence we employ a drift force to model the unresolved dynamical friction, which generally ensures that the BHs remain close to the centre of the halo. However, in dwarf galaxies, off-centre BHs may actually be physical (see e.g. Pfister et al. 2019; Bellovary et al. 2021; Ma et al. 2021) and in particular for massive BH binaries, the BHs orbiting around the centre could be an important contribution to core formation via core scouring. This has mainly been investigated in the context of massive elliptical galaxies (see e.g. Rantala et al. 2018, 2019), however, these mechanisms may also translate to the dwarf regime and, in turn, the central dark matter densities of dwarfs are expected to affect binary shrinking times (see Tamfal et al. 2018). To also account for the dynamical effects of massive BHs on core formation in future work, it will be important to include the unresolved effects of dynamical friction as accurately as possible following dynamical friction estimators based on high-resolution N-body simulations (e.g. Ma et al. 2023; Genina et al. 2024; Partmann et al. 2024).", '5 CONCLUSIONS': "In this work, we have investigated the impact of feedback from massive BHs on cusp-to-core transformations and the associated 'diversity of rotation curves problem' in dwarf galaxies with a new suite of high-resolution cosmological zoom-in simulations based on the /f.pc/i.pc/r.pc/e.pc-3 galaxy formation model. Our suite includes three different dwarf haloes spanning a range of masses and environments within the classical dwarf regime of 8 × 10 9 M ⊙ < 𝑀 halo < 4 × 10 10 M ⊙ . For each of these haloes, we investigate simulations with and without AGN as well as with and without CRs yielding four physics variations for each of the three haloes. We note that these different set-ups may represent different outcomes of the same underlying galaxy formation physics under different conditions, with the efficiency of both AGN and CRs expected to be highly variable in dwarf galaxies. Our main conclusions are as follows:  \n- (i) AGNmaydrive core formation directly as an additional source of feedback dynamically heating the central dark matter distribution; see the m10zCR+AGN set-up.  \n- (ii) AGNmayalsoenhance core formation indirectly by suppressing high-redshift star formation and shifting SN activity to later times, which favours core formation at 𝑧 = 0; see e.g. the m10qAGN and m10zAGN set-ups.\n- (iii) CRs may also indirectly affect cusp-to-core transformations by regulating star formation histories - for most of our simulations, this leads to suppressed star formation and therefore suppressed core formation; see e.g. the m10qCR and m10yCR set-ups.\n- (iv) Our simulation suite is in good agreement with observed dark matter central density distributions, broadly following the trends from the nearby dwarfs in Read et al. (2019). In particular, cores are generally more pronounced for higher stellar-to-halo mass ratios, 𝑀 stellar / 𝑀 vir . The only exception to this trend comes from the AGN-driven core in the m10zCR+AGN set-up, which has a strongly suppressed 𝑀 stellar / 𝑀 vir ratio and a cored profile. This parameter space is a key target for observational searches for cores induced by AGN feedback.\n- (v) BH feedback and CRs create a variety of cusps and cores in circular velocity profiles, with correlations between rotation curve shapes and baryonic influence that align with observations, potentially helping to resolve the diversity of rotation curves problem in dwarf galaxies.  \nOverall, our findings suggest that AGN in dwarf galaxies can significantly impact central dark matter densities. BHs influence dwarf galaxy rotation curves through two key mechanisms: directly driving core formation via AGN feedback and indirectly regulating core formation by altering star formation histories and SN activity. CRs further contribute to this regulation by enhancing AGN feedback and suppressing star formation, thereby influencing core formation. The combined effects of AGN and CRs lead to a spectrum of central dark matter densities, reflecting the varying levels of AGN activity, star formation, CR influence, and their non-linear interaction. Due to the bursty star formation and accretion histories in dwarfs, AGN feedback is likely to occur in a highly stochastic manner; this may explain the wide range of dark matter central densities inferred from dwarf observations at fixed halo mass. Linking surveys of dwarf rotation curves with searches for AGN signatures will be crucial for future observations to elucidate the impact of BHs on dark matter profiles and to determine whether the diversity of rotation curve problem may be explained by baryonic processes.", 'ACKNOWLEDGEMENTS': "The authors are grateful for helpful discussions with Vasily Belokurov, Martin Bourne, Jenny Greene, Vid Irsic and Sergio MartinAlvarez. The authors would also like to thank Isabel Santos-Santos for providing the observational rotation curve data. The simulations presented in this work were run on the Flatiron Institute's research computing facilities (the Iron compute cluster), supported by the Simons Foundation. SK has been supported by a Flatiron Research Fellowship and a Junior Research Fellowship from St Catharine's College, Cambridge. The Flatiron Institute is supported by the Simons Foundation. DAA acknowledges support by NSF grant AST-2108944, NASA grant ATP23-0156, STScI grants JWSTGO-01712.009-A and JWST-AR-04357.001-A, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023028 by the Research Corporation for Science Advancement. Support for ISS was provided by NSF Collaborative Research Grant 2108318. SW received support from the NASA RIA grant 80NSSC24K0838.", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.'}