Diverse dark matter profiles in FIRE dwarfs: black holes, cosmic rays and the cusp-core enigma
Sophie Koudmani Douglas Rennehan Rachel S. Somerville Christopher C. Hayward Daniel Anglés-Alcázar Matthew E. Orr Isabel S. Sands Sarah Wellons


firstpage–lastpage
An efficient observational strategy for the detection of the Oort cloud
    
Eran O. Ofek1,⋆,
Sarah A. Spitzer1,
Guy Nir2

    Received ; accepted 
=======================================================================




§ 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 fire-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.





methods: numerical – galaxies: active – galaxies: dwarf – galaxies: evolution – galaxies: formation – dark matter








§ 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 (similar to the mass of the Large Magellanic Cloud), dwarf galaxies are crucial astrophysical probes in near-field cosmology and galaxy formation studies.

Due 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 <cit.>. Key issues include the missing satellites problem, which notes a deficit in the observed number of Milky Way satellites compared to theoretical predictions <cit.>, 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 <cit.>. 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 <cit.>.

These discrepancies have led to explorations of alternative dark matter models, such as warm dark matter (WDM) <cit.>. WDM models, which include 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 <cit.>. However, WDM struggles with the cusp-core problem, as the resulting cores are not large enough <cit.>. Another approach involves self-interacting dark matter (SIDM), where particles scatter elastically, creating cored density profiles and reducing halo ellipticity <cit.>. Constraints from strong lensing measurements of galaxy clusters suggest a velocity-dependent cross section that decreases from dwarf to cluster scales <cit.>. 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 <cit.>.

However, 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 <cit.> thereby resolving the missing satellite problem. Furthermore, stellar feedback <cit.> 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 <cit.> and that late-time star formation is crucial for maintaining these cores <cit.>. 

Notwithstanding 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 <cit.>. Due to the self-similar nature of structure formation in the ΛCDM universe, the full mass profile of a halo can be characterised by a single parameter such as the virial mass <cit.>. 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 <cit.>. Galaxy formation simulations based on CDM have traditionally struggled to reproduce this observed diversity of dwarf rotation curves <cit.>. Some groups find that SIDM may be able to improve these discrepancies <cit.>, whilst others find no preference between SIDM and CDM simulations with stellar feedback <cit.>.

The 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 <cit.> whilst others predominantly produce bursty star formation histories and cores <cit.>. 

One 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 <cit.>. 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 galaxies, 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 <cit.> to optical <cit.> to IR <cit.> and to radio observations <cit.>. 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 <cit.>.

These tantalising observations have motivated theorists to investigate this largely unexplored regime. Analytical models point towards favourable AGN energetics in the dwarf regime <cit.> and suggest that AGN activity may be able resolve all of the remaining `dwarf galaxy problems' <cit.>. 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 <cit.> and black holes (BHs) may be wandering in dwarf galaxies due to their shallow potential wells <cit.> 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-supply-limited or torque-driven schemes due to its strong dependence on BH mass <cit.>. The cosmic environment also has a crucial influence on AGN activity in the dwarf regime with minor mergers triggering AGN episodes whilst long-term residence in dense environments is detrimental to BH growth <cit.>.

Despite 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 <cit.> and suppressing star formation <cit.>. 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? 

The interaction between AGN feedback, star formation, and dark matter in dwarf galaxies has not been thoroughly explored, though some trends have been noted, such as AGN mildly suppressing central dark matter densities <cit.>, 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 <cit.> employed for modelling the interstellar medium (ISM) in these simulations, which leads to smoother star formation histories alongside less bursty AGN feedback at low redshifts, often operating as `maintenance-mode' feedback <cit.>. Recently, <cit.> 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 <cit.>.

It 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) <cit.>. In particular, it has been shown that CRs may drive core formation <cit.> though their impact in the dwarf regime remains controversial as gas-rich mergers may re-establish cores in these systems <cit.>. What is more, it has been found that CRs may significantly boost the efficiency of AGN feedback <cit.> 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.

The aim of this paper is to systematically investigate the impact of AGN feedback 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 fire-3 galaxy formation framework and contrast the resulting dark matter profiles, their cosmic evolution and their connection with different feedback processes.

The remainder of this paper is structured as follows. In Section <ref>, we introduce our three dwarf simulation set-ups and describe the different variations of the fire-3 galaxy formation model explored in this work. The results of our analysis are presented in Section <ref>. Firstly, we focus on the stellar and BH properties of our simulated dwarfs, including visualisations of the stellar distributions in Section <ref>, comparisons with observed nearby low-mass galaxies (see Section <ref>) and observed BH masses in dwarfs (see Section <ref>). We then analyse how our different feedback set-ups influence the dark matter profiles at z=0 in Section <ref> and examine the driving forces behind the diversity of dark matter profiles in <ref>. We investigate the cosmic evolution of the central dark matter densities in Section <ref> and assess whether this could explain the observed diversity of rotation curves in Section <ref>. We discuss our results and caveats to our modelling in Section <ref> and conclude in Section <ref>.



§ METHODOLOGY
 



 §.§ General set-up


Our simulation suite is based on a set of 3 initial conditions (ICs) of dwarf haloes selected from the fire ICs set, spanning a range of halo masses in the `classical dwarf' regime.  We select two low-mass dwarfs: m10q and m10y, with z=0 halo masses of =8×10^9  and =1.4×10^10 , respectively.  Our third halo, m10z, is an intermediate-mass dwarf with z=0 halo mass = 3.5×10^10 . All of these ICs were generated using the music IC generator <cit.> and have been thoroughly explored with the fire-2 model <cit.>. In the standard fire-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. 

To simulate these systems, we use the gizmo hydrodynamics+gravity code, employing the mesh-free finite mass method <cit.> that tracks constant-mass gas resolution elements in a Lagrangian fashion.  Our simulations include both baryons and dark matter, with mass resolutions of = 262  and = 1303 , respectively. The gravitational softening for the dark matter, stars and BH particles is set to ϵ_soft = 0.0625  ckpc  h^-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 <cit.>.



 §.§ Galaxy formation model


For the galaxy formation physics, we use the newly updated fire-3 model <cit.> which tracks the state of the multi-phase ISM and multiple forms of stellar feedback, including feedback from SNee[We couple the SN momentum following <cit.> and do not yet include the updates to the momentum coupling from <cit.>. Though we note that these changes mainly affect high-mass galaxies at low resolution, whilst the impact on high-resolution dwarfs is very limited.] 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 <cit.> and <cit.>. The fire-3 model also has new optional physics modules compared to the previous fire <cit.> and fire-2 <cit.> 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 <ref> for an overview of the simulation suite.




The modelling of BHs and CRs remains highly uncertain. Several studies have investigated the impact of different modelling assumptions for CRs <cit.> and BHs <cit.> within the fire model. For our study, we base our BH model on the extensive study by <cit.> 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 <ref> and point the interested reader to <cit.> 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 in <cit.> allowing us to explore several different IC set-ups to z=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 <ref> and refer the interest reader to <cit.> for a more detailed overview of this approach.



  §.§.§ BH modelling
 
For our BH modelling, we largely adopt the default fire-3 BH implementation <cit.>, 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 <cit.>. Below we summarise the salient details as well as our modifications to the fire-3 BH model and refer the interested reader to <cit.>  and <cit.> for a more in-depth description and parameter studies.

BHs 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  pc^2) and low metallicities (Z_gas≲0.001) in order to spawn BH particles of mass =100. With a DM particle mass of M_DM∼ 10^3, 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 <cit.> 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 <cit.>. 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.

Once a BH particle is present, it grows via accretion by draining gas from the nearest N∼256 gas cell neighbours.  <cit.> parameterized the accretion rate as ≡η_accΩ, where η_acc is the accretion efficiency,  is the gas mass within the BH kernel, and Ω = √(GM_tot/R) is the orbital frequency determined at the kernel size R, based on the total mass enclosed with R, M_tot. fire-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 = M_α disc / t_dep, where M_α disc is the mass of the accretion disc and t_dep = 42 Myr (1 +  / M_α disc)^0.4 is the depletion time.

While <cit.> 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 <cit.>. The accretion efficiency in the case of gravitational torque-dominated accretion comes from the work in <cit.>,


    η_acc = C ([ + M_d] / M_d)^1/6/1 + 3M_d,9^1/3(/M_d),


where C is a resolution-dependent constant,  is the BH mass, M_d is the mass of angular-momentum support material, and M_d,9 = M_d / 10^9.  The model has been studied extensively in the context of cosmological zoom-in simulations in <cit.> and <cit.>.  It is important to note that we choose a value of C = 0.1 in order to have reasonable BH masses within our three simulated dwarf galaxies.  However, the parameter C cannot take into account all sub-grid processes that may occur below the resolution scale.

One important piece of (partially) unresolved physics is the impact of stellar feedback on the accretion flow in the vicinity of the BH.  <cit.> introduce a scaling factor f_acc (such that η_acc→ f_accη_acc) that depends on the local gravitational acceleration, a_g,eff = GM(<R)/R, as 


    f_acc = a_g,eff/a_g,eff + a_g,crit,


where a_g,crit≈ 10^-7 cm s^-2.  Below a_g,crit (an effective surface density of Σ_tot∼ 3000  pc^2), a significant fraction of gas is expelled due to stellar feedback <cit.>. 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).

The accretion process itself drives powerful outflows, jets and radiation that influence the properties of the galaxies hosting the BHs, collectively coined BH feedback. The fire-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 .

As gas falls 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 mass-energy 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 fire-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 <cit.>, except using a quasar template spectrum <cit.>.  The radiation momentum flux depends on the luminosity absorbed (L_abs) in the gas as ṗ = L_abs / c.  

There 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 fire-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 =. The kinetic energy rate in the wind itself depends on the wind velocity at launch, , such that Ė_mech = ^2 / 2.  <cit.> recast the mechanical luminosity in terms of the new mechanical efficiency, η_mech≡^2 / (2c^2), such that Ė_mech = η_mech c^2.  In this work, we choose = 10,000 (see e.g. , compiled in fig. 3 of ) and, hence, have η_mech = 5.6×10^-4. 

The fire-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 <ref>, but note here that the CRs are injected at a rate Ė_CR = η_CR c^2, where η_CR is the efficiency of CR feedback. In this work, we use an efficiency of η_CR = 0.01 which was identified as yielding realistic galaxy properties in the <cit.> parameter study.

For the mechanical and radiative feedback, we use the solid angle weighting to inject the AGN feedback <cit.>. Following <cit.>, 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


    w(θ) = ϵ_jet (ϵ_jet + cos^2θ)/(ϵ_jet+1)(ϵ_jet+(1-cos^2θ)).


We 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, J_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 <cit.>, 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 <cit.>. Long-range radiation transport and CRs are treated via the LEBRON approximation and injected isotropically <cit.>.





  §.§.§ CR modelling
 

The fire-3 model includes a comprehensive CR model that explicitly follows the multi-species/multi-spectral CR dynamics <cit.> with the entire distribution function <cit.> 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 <cit.> and briefly describe our motivation for this model below. 

The primary assumption is that the CR energy spectrum is spatially constant, and that CRs above ≳ 1 GeV dominate the CR pressure <cit.>.  This “single-bin” model only evolves the integrated CR energy 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 <cit.> showed that higher values of the isotropic diffusion coefficient κ are necessary to explain γ-ray luminosities in dwarf and L^*-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 P dV work on the low density gas also becomes negligible.

Given that several of the CR physical processes are seemingly subdominant in the dwarf regime, <cit.> 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 (∂_t e_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 <cit.>.  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 v_κ.  An approximate solution for the spatial distribution of CR energy density for such a model is presented in <cit.>, and depends on an integral of an exponentially-declining function of the loss function.

The resultant CR model is solved using a LEBRON-type approximation <cit.>. 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 fire-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 <cit.>; κ = 5×10^28 cm^2 s^-1 and v_κ = 20 km s^-1.
 


§ RESULTS
 


 §.§ Overview of the simulation suite
 

We begin our analysis of the simulation suite by inspecting visualisations of our simulated dwarfs' stellar distributions. Figure <ref> 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 of 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. 

The three rows correspond to our three different ICs whilst the different columns represent our four physics variations based on the fire-3 physics modules: the fiducial fire-3 set-up (without CRs or AGN), fire-3 with AGN, fire-3 with CRs and fire-3 with both CRs and AGN. We note that there are several uncertainties associated with modelling AGN and CRs in galaxy formation simulations and here we only explore a small subset of this parameter space that was identified by <cit.> 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 <cit.> which only accounts for local losses and assumes a constant diffusion and effective streaming speed. See Section <ref> and <cit.> for more details.

Keeping 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 <cit.>. 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-AGN 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 <cit.> 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. 



 §.§ 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. <ref>, 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.



  §.§.§ SMHM relation


Firstly, we examine if the simulated dwarf galaxies align with observational expectations for the SMHM relation. The left panel of Fig. <ref> shows our 12 dwarf zoom-in simulations in M_halo – M_stellar space at redshift z=0.

For comparison, we plot several SMHM relations from the literature, including <cit.> and <cit.>. The low-mass end of the SMHM relation remains largely unconstrained by observations, evident from the significant differences between these models. The <cit.> models extend down to log(M_halo / ) = 10.0, and the <cit.> relation was shown to also be valid down to the dwarf regime, at least to log(M_halo / ) = 10.0, by <cit.>. We plot the extrapolated relations by <cit.> and <cit.> 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 <cit.>.

Our 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 <cit.> and <cit.> 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 <cit.>, <cit.>, and <cit.> as shaded regions. The differences in these relations arise due to different assumptions for surface-brightness incompleteness correction, see discussion in <cit.>. We also show various observed data points from the SPARC data base <cit.> to indicate the scatter in the observations.

We 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 z=0 or the halo mass based on the collapse of a spherical top-hat perturbation following <cit.>. We have checked for our simulated 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.

From 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  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 <cit.>.



  §.§.§ Galaxy colours


Next, we compare the galaxy colours of our simulated dwarfs with observational data. The middle panel of Fig. <ref> presents the g - r 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 <cit.> as grey squares and the ELVES dwarf satellites, a nearly volume-limited sample of Milky Way-like hosts in the Local Volume, from <cit.> as grey diamonds.

We use the <cit.> stellar population synthesis model to calculate galaxy colours as a function of stellar age and metallicity, assuming a <cit.> IMF, including all star particles within twice the stellar half-mass radius. For the MW dwarf satellites, we convert the B - V colours using the transformation equations from <cit.>, 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 simulated galaxy 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.

All 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 g-r∼0.35, whilst most of our dwarfs are offset towards somewhat redder colours, clustering around g-r∼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. 

Interestingly, the  set-up has a similar stellar mass to its AGN counterpart , 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 <ref>. Conversely, the  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, , 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. 



  §.§.§ 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 z=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 <cit.>, which also includes the Milky Way dwarf measurements from <cit.> 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. 

As we found from our visual inspection of the stellar distributions in Fig. <ref>, 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  with ) or more compact distributions (e.g., compare  with ). The differences between the AGN and no-AGN set-ups are more pronounced if CRs are present, in agreement with the findings from <cit.> that CRs enhance the efficacy of the AGN feedback in fire.

Overall our stellar sizes are in good agreement with the observations, though the  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 <cit.> or more extended observed distributions <cit.>, with the stellar sizes also being significantly dependent on the wavelength <cit.>.



 §.§ 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.



Fig. <ref> shows our simulated dwarfs in stellar mass – BH mass space. The evolutionary tracks of the dwarf simulations from z=4 to z=0 are plotted as dotted lines, colour-coded by simulation setup. Note that the BHs are generally seeded around z ∼ 8–10 with the exact seeding time depending on the metallicity and density conditions. However, we here focus on later redshifts as by z=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.

We 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. 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 estimates for active dwarf galaxies from <cit.> and <cit.> as dark-grey squares and stars, respectively, and measurements from <cit.> (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 M_stellar∼ 10^9 and, excluding upper limits, the lowest BH mass estimate based on TDEs stems from a host galaxy of stellar mass M_stellar = 2.5 × 10^8. Below this mass, we are firmly in the extrapolated regime and indicate this by greying out the local BH – galaxy scaling relations below this mass.

To extend our comparison to lower-mass systems, we include IMBH candidates from globular clusters around the Milky Way and M31, compiled by <cit.>, 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 <cit.>. 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 <cit.> and a new IMBH detection in M31's most massive globular cluster which may also originate from a stripped dwarf galaxy <cit.>, both indicated in olive-green.

We 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 <cit.> relations are based on dynamical measurements by <cit.> as well as more recent dynamical measurements and upper limits in the low-mass 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 late-type relation has a significantly lower normalization, likely due to a weaker correlation between disc mass and BH mass <cit.>. We also show the <cit.> scaling relation, based on local AGN, including dwarfs, as a dashed-dotted black line. This relation is significantly shallower than that of <cit.>, leading to larger BH masses in the low-mass dwarf regime. 

For the classical dwarf mass range covered by our simulations, M_stellar≲ 10^8, we would expect a flattening of the scaling relations, with the transition mass dependent on the seeding mechanism <cit.>, 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 <cit.> 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[The scatter for the M_BH/M_NSC ratio at a given galaxy mass spans over three orders of magnitude, showing a trend of increasing ratio with galaxy mass, and BHs only start to dominate at stellar masses above a few times 10^10.]. <cit.> compile NSC mass measurements extending to low-mass dwarfs with stellar masses ≲ 10^6, making the relation valid for the entire z=0 galaxy stellar mass range considered here.



Keeping 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  which closely follows the early-type relation from <cit.> and is only minimally offset from the <cit.> 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 <cit.> and well within the bounds provided by the nuclear star cluster relation. 

For the redshift evolution, we plot the BH mass of the most massive progenitor for redshifts z=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 z=1 and z=0, all simulations with CRs only have minimal stellar mass growth.  For  and , 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  BH growth is already shut down from z=3, however, between z=4 and z=3, this run similarly experiences a rapid BH growth phase whilst star formation is already quenched.

We 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, this represents only a very small fraction of the final BH mass in all cases.



 §.§ Dark matter distributions
 

Fig. <ref> 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[Note that we also investigated alternative centering methods, including the shrinking sphere method from <cit.>, however, we find that in practice the differences are minimal as the common halo centre definitions are in good agreement for our dwarfs.]. In each case, we plot the mean z=0 profiles, averaged over Δ T=500 Myr, for our four different galaxy formation model variations as colour-coded 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 Δ T=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  simulation, however, experiences significant fluctuations in the central density from ∼ 12–13 Gyr due to mergers with substructures, so for our z=0 profiles we choose the bin width to 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 r=150 pc, a characteristic radius for distinguishing cusps from cores in observations <cit.>, with colour-coded symbols. 

Furthermore, we indicate the gravitational softening length at r=62.5 pc and the Power radius[We note that the required particle number for resolving the central region depends on the radius and density contrast in question but generally ranges between N=1000–3000; here we follow <cit.> and set N=2000.] <cit.> 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 <cit.> 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.
 
First of all, we find that the introduction of AGN and CRs introduces a significant amount of scatter. The simulations with the default fire-3 model, without CRs or AGN, all lead to distinct cores at r=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[We note that the central densities at r=150 pc are within the Power radius in most cases, however, for all cases, the cusp versus core categorisations based on the 2σ scatter around NFW expectations are unchanged between r=150 pc and the Power radius, demonstrating that the qualitative results are robust.].

For 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 r=150 pc.

For 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. 



Overall, 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 <cit.>. 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 <cit.> whereas we observe variations ranging from 70 to 90 per cent for our different sets of ICs and galaxy formation physics configurations (see Table <ref>). Furthermore, as we show in the next Section, these central density trends are relatively steady as a function of cosmic time 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.



 §.§ Driving forces of core formation
 

Our analysis from Section <ref> 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. 

Fig. <ref> shows the central dark matter densities at r=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 Δ T=500 Myr and indicate the mean densities as solid lines and the standard deviation in each bin by the shaded regions. 

The m10q set-up is shown in the first panel. Reflecting its early formation history, the dark matter density trends we found at z=0 in Fig. <ref> have persisted for the past 6 Gyr. The central densities are relatively steady with only  showing significant evolution after redshift z=1 with the central density of this run showing a slow decline.

Similarly for m10y (middle panel), the z=0 trends have persisted for the past ∼5 Gyr. Both of the CR-based runs,  and , do not show any significant evolution at low redshifts, whilst the runs without CRs,  and  show a steady decline from z∼1 onwards. Note that the density fluctuations for these runs are also much stronger, reflecting the potential perturbations characteristic of core formation <cit.>. Focusing on early times between redshifts z=4 and z=1, we find that most of the runs are quite similar, except for  which exhibits a strong core during that period. We also note that at t ∼ 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  set-up is erased by this merger.

For the m10z system (right panel), there is significantly more late-time evolution as well as stronger fluctuations as reflected in its more irregular morphology, see Fig. <ref>. Indeed, the general central density trends we observe at z=0 have only been in place for the past ∼2 Gyr and the  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 z=2 and z=0.5, the set-ups , , and , all show a relatively similar redshift evolution of their central densities. The  set-up, however, closely follows the DMO run and even surpasses it from t=7–9 Gyr, indicating strong baryonic cooling at the centre. At z∼ 0.5 a merger induces density fluctuations for all of our m10z set-ups, but most prominently for the CR-based runs  and . After this merger, the two CR-based set-ups have relatively steady central densities. On the other hand, the runs without CRs,  and  experience significant density suppressions at late times.



In 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.



  §.§.§ SN feedback
 

Fig. <ref> 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. The three columns correspond to our three sets of ICs. The no-AGN and AGN runs are shown in the top and bottom row, respectively. 

We 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 <ref>) and it also corresponds to the scale of the half-mass radius for most of our systems (see Fig. <ref>).

The SN energy injection rates are recalculated from the stellar ages and metallicities based on Starburst99 tables for a <cit.> IMF. We focus here on the SN II rates, which constitute the majority of SN events <cit.>, for simplicity. For reference, we also indicate the z=0 binding energy of the respective haloes as dotted lines, which we estimate as f_bar  M_vir  V_vir^2, where f_bar is the universal baryon fraction and V_vir is the virial velocity of the halo <cit.>.

Firstly, 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 <cit.>, due to their shallow potential wells, dwarf systems are very sensitive to baryonic feedback processes so that even small changes can significantly affect star formation histories.

This 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 z ∼ 6. Hence both  and  end up with cuspy profiles, even though the cumulative injected SN energy for  slighly exceeds the z=0 binding energy, as the absence of late-time star formation means that a cusp can be re-established quite easily after the initial burst. Note that  has a slighly reduced central density at z=0 compared to  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  set-up, there are a few more powerful bursts at intermediate redshifts, which then quench the systems at much later times from z∼ 1.5. For the  set-up, the AGN regulation means that the intermediate redshift star formation bursts are less intense so that star formation can continue until z=0. Both  and  end up with cumulative SN energy almost an order of magnitude above the z=0 binding energy leading to cored profiles; however,  has a more extreme core due to the continued star formation activity <cit.>.

For 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 z ≳ 1. This is still sufficient to (cumulatively) exceed the z=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  and  have continued star formation activity until z=0, higher cumulative injected SN energies, and therefore much stronger cores. Adding in AGN only has a small impact for both  and  compared to their no-AGN counterparts. Again the bursts at high redshift are somewhat weaker, allowing for more SN activity at lower redshifts and therefore somewhat stronger cores. However, as discussed in Sections <ref> and <ref>, the differences induced by AGN feedback are only minimal for the m10y profiles at z=0.

For 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, , , and  all experience episodes of star formation for most of cosmic time, although for , those are shifted to lower redshifts, whilst for the  simulation SN activity is more concentrated at higher redshifts.  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  and  set-ups also exceeds the z=0 binding energy and both of these result in cored profiles. Interestingly, however, star formation in the  set-up is significantly suppressed, and powerful bursts are mostly restricted to high redshifts (z>2), so that the total injected SN energy remains well below the z=0 binding energy by a factor of three. This indicates that for the   set-up, the core formation may be driven by the AGN.

This analysis underlines that CRs and AGN primarily affect cusp-to-core transformation, and the diversity of dark matter profiles, by regulating star formation histories in our simulations. The only set-up where the AGN feedback is actively driving core formation is , 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 set-ups. However, in these scenarios, the AGN may still enhance core formation; we examine this possibility in the next section.



  §.§.§ AGN feedback
 



In the fire-3 model, AGN feedback is injected in several channels including mechanical feedback, radiative feedback and CRs (see Section <ref> 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.

For 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 <ref>). 

For the radiative feedback, on the other hand, the full luminosity is injected into the simulation domain using the LEBRON method <cit.> 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 L/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 <cit.>. We here assume that a similar fraction of the luminosity would couple to the gas for radiative AGN feedback. Following <cit.>, 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.

For the AGN-driven CR feedback, 1 per cent of the accreted rest-mass 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 <cit.>, 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.

Keeping 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. <ref>. 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 z=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.

Overall, 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  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 in 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 <cit.>.

For m10q, in the CR-based run , 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  set-up, there is a sharp AGN burst at t ∼ 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.

For m10y, the  set-up has a strong AGN burst at t∼ 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  set-up) and the relatively low levels of activity do not have a significant impact on core formation.

For m10z, contrary to the other two sets of ICs, the energy injected from the CR-based run, , exceeds the AGN energy injected with the equivalent no-CR set-up . This is mostly driven by a significant burst in the  set-up at t∼ 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 t∼ 10 Gyr in the  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  run is quite large (M_BH/M_d∼ 10^4), which leads to depletion timescales of the order of ∼ 1 Gyr (also see Section <ref>), allowing for persistent AGN activity even after the inflows onto the BH – accretion disc particle have subsided. Secondly, whilst 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 <ref>. For the  set-up the AGN has a significant burst at t∼12 Gyr which again is associated with strong central density fluctuations suggesting that here the AGN may also enhance the formation of the core.



Overall, we note for most of our AGN set-ups, including , , , , and , the equivalent no-AGN set-up already has a SN-induced core. For all of these cases, apart from , 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' SN activity to later times favouring core formation. Indeed, it is likely that both of these factors contribute to the enhancement of cores.

The  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. <ref>). Here, there is no SN activity for the last ∼3.5 Gyr from z∼ 0.5 and only weak SN activity between z=2 and z=0.5, yet the core is maintained. Indeed, the central density significantly decreases from z=2 to z=1 and is then maintained at low levels despite the late-time merger at z ∼ 0.5, providing strong evidence for AGN-driven core formation.

Overall, 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.



 §.§ 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. <ref> shows the central dark matter densities at r=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). The dotted lines indicate the (binned) redshift evolution of our simulated dwarfs from z=3 to z=0, again binned over 500 Myr. For clarity we only highlight the z=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 coreNFW profile from <cit.> 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 <cit.>.

The 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,  and , 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 <cit.>, 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.

The middle panel shows the simulated and observed central densities as a function of stellar mass to virial mass ratio M_stellar/M_vir. As discussed in <cit.>, the observed dwarfs show a clear trend with higher M_stellar/M_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 cores form below a critical ratio of ( M_stellar/M_vir)_crit= 5 × 10^-4 since the integrated SN energy is insufficient to drive the required potential fluctuations <cit.>. We also find a strong correlation between M_stellar/M_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 M_stellar/M_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 M_stellar/M_vir trend, however, is more notable: for the  set-up the central densities as a function of M_stellar/M_vir are steadily decreasing so that a core forms despite the z=0 M_stellar/M_vir ratio being below the critical ratio identified by <cit.> and <cit.>. This provides further evidence that this run could be experiencing AGN-driven core formation.

We further examine this possibility in the third panel where we plot the simulated central densities (as well as their binned redshift evolution from z=3) for the AGN-based simulation set-ups as a function of BH mass to virial mass ratio M_BH/M_vir. For the runs without CRs, , , and , the central densities again steadily decrease as a function of M_BH/M_vir. However, we caution that this does not necessarily indicate that the BHs are contributing to core formation due to the strong correlation between M_BH and M_stellar (see Fig. <ref> for scaling relations of our simulated dwarfs). The  set-up is the only simulation where the central density ratio decreases as M_BH/M_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  set-up has a very shallow gradient again reflective of the relatively quiescent late-time evolution. Finally, the  set-up exhibits a significant decrease in central densities as the M_BH/M_vir increases (opposite trend from M_stellar/M_vir) providing further support for the BH-driven core formation scenario.






 §.§ The diversity of rotation curves
 

As discussed in Section <ref>, 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.

In Fig. <ref>, we present circular velocity curves for all our dwarf simulations. As for the dark matter profiles, we average the velocity curves over Δ T=500 Myr and indicate the 1σ scatter with the shaded regions. Following the methodology presented in <cit.>, 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.

The 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. <ref>. 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 <ref>). For the cored dwarfs this then leads to slowly rising rotation curves for both dark matter and stars. Dwarfs that are quiescent at z=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.

Following <cit.>, we quantify the cusp versus core nature of the profiles by calculating characteristic velocity ratios, comparing the inner `fiducial' rotation velocity V_fid with the asymptotic flat rotation velocity V_max. The fiducial rotation velocity is measured at the fiducial inner radius defined as:


    r_fid = (V_max/35  km  s^-1)  kpc.


The rotation curve shape parameter η_rot can then be defined as 


    η_rot = V_fid/V_max = V(r_fid)/V_max.


Rapidly 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 <cit.> for a detailed discussion of this parameter.

As in <cit.>, we also inspect a second velocity ratio – the baryonic importance parameter, which is defined as


    η_bar = ( V_b,fid/V_fid)^2 = ( V_b(r_fid)/V (r_fid))^2,


where V_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.

The 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. <ref>. For comparison, we also show observed data as collated by <cit.>. We restrict the observed galaxy sample to the classical dwarf regime with V_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 <cit.> we therefore show our low-η_bar simulated dwarfs at η_bar = 0.05 for clarity. 

Focussing first on the observations, we find that, as discussed in <cit.>, 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 the 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. 

Our simulations cover almost all observed scenarios: cusps and weak cores at low baryonic dominance (e.g.,  and ) and cusps and strong cores at high baryonic dominance (e.g.,  and ). 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 <cit.>. 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 fire-3 model (as well as in the SPARC sample), these transient ultra-cuspy dwarfs are generally associated with more massive systems (M_halo≳ 10^11 and V_max≳ 50  km   s^-1), see <cit.> for details.

Nevertheless, 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 <cit.>. 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.



§ DISCUSSION
 



 §.§ 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.

Several groups have investigated the role of stellar feedback in driving cusp to core transformations in galaxy formation simulations with mixed results <cit.>. 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 <cit.>. Nevertheless, for a given ISM model it has proven very challenging for simulations to reproduce the observed diversity of dwarf rotation curves <cit.>. 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.

Previous studies found systematic yet small-scale suppressions of central dark matter densities by AGN feedback in dwarfs <cit.>, yet in both cases, the simulations did not explicitly resolve the multiphase ISM. Cusp-to-core transformations with AGN feedback were also explored by <cit.>. 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.

In our study, we investigate, for the first time, the impact of AGN feedback on dark matter profiles in dwarfs with a resolved multi-phase ISM model and a BH growth scheme that does not suppress AGN activity in dwarf galaxies. Several works exploring BH physics within the fire model have found that CRs appear to be a key ingredient for effective AGN feedback <cit.>, 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 <cit.>, who find that CRs are prone to enhance core formation. We note that for one of our set-ups (compare  to ), 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. 

We 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.

One 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 <cit.>, 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.



 §.§ 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 <cit.>. 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 <cit.>.

Apart 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 AGN in cusp-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. <cit.>, 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 <cit.>. From our simulations, we predict that if cored dwarf galaxies were found below the critical mass ratio of (M_stellar/M_vir)_crit = 5 × 10^4 this would be a strong indicator of AGN-driven core formation in dwarfs. 



 §.§ 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 in-between 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 in future work. Below we discuss the caveats of our CR and BH modelling in more detail.



  §.§.§ 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 <cit.>, which allows for the inclusion of the impact of CRs on galaxy formation without the significant computational overhead associated with full CR transport.

As 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 <cit.>. 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 <cit.> also find in their validation simulations with a 10^11 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.



  §.§.§ 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. 

Here, we model the BH growth based on the torque-based accretion model <cit.>, 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 <cit.>. 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 <cit.>. Nevertheless, from an accretion perspective, our simulations may be exploring an upper limit on the impact of BHs in the dwarf regime.

The 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. Secondly, 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 <cit.>.

For 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. <ref>. However, as discussed in <cit.>, there are several other `plausible' AGN feedback set-ups within the fire 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 <cit.> 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  set-up, without compromising in terms of producing realistic stellar properties.

The 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 <cit.> 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 <cit.>, 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 <cit.>. 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 <cit.>.



§ 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 fire-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. 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:

    
  * AGN may drive core formation directly as an additional source of feedback dynamically heating the central dark matter distribution; see the  set-up.
    
  * AGN may also enhance core formation indirectly by suppressing high-redshift star formation and shifting SN activity to later times, which favours core formation at z=0; see e.g. the  and  set-ups.
    
  * 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  and  set-ups.
    
  * Our simulation suite is in good agreement with observed dark matter central density distributions, broadly following the trends from the nearby dwarfs in <cit.>. In particular, cores are generally more pronounced for higher stellar-to-halo mass ratios, M_stellar / M_vir. The only exception to this trend comes from the AGN-driven core in the  set-up, which has a strongly suppressed M_stellar / M_vir ratio and a cored profile. This parameter space is a key target for observational searches for cores induced by AGN feedback.
    
  * 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.


Overall, 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 Martin-Alvarez. 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 JWST-GO-01712.009-A and JWST-AR-04357.001-A, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 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.




mnras
 
























