The Pandora project – I. The impact of radiation, magnetic fields, and cosmic rays on the baryonic and dark matter properties of dwarf galaxies Free

ABSTRACT

Enshrouded in several well-known controversies, dwarf galaxies have been extensively studied to learn about the underlying cosmology, notwithstanding that physical processes regulating their properties are poorly understood. To shed light on these processes, we introduce the Pandora suite of 17 high-resolution (3.5 parsec half-cell side) dwarf galaxy formation cosmological simulations. Commencing with magneto-thermo-turbulent star formation and mechanical supernova (SN) feedback, we gradually increase the complexity of physics incorporated, ultimately leading to our full-physics models combining magnetism, on-the-fly radiative transfer and the corresponding stellar photoheating, and SN-accelerated cosmic rays. We investigate multiple combinations of these processes, comparing them with observations to constrain what are the main mechanisms determining dwarf galaxy properties. We find hydrodynamical ‘SN feedback-only’ simulations struggle to produce realistic dwarf galaxies, leading either to overquenched or too centrally concentrated, dispersion-dominated systems when compared to observed field dwarfs. Accounting for radiation with cosmic rays results in extended and rotationally supported systems. Spatially ‘distributed’ feedback leads to realistic stellar and H i masses, galaxy sizes, and integrated kinematics. Furthermore, resolved kinematic maps of our full-physics models predict kinematically distinct clumps and kinematic misalignments of stars, H i, and H ii after star formation events. Episodic star formation combined with its associated feedback induces more core-like dark matter central profiles, which our ‘SN feedback-only’ models struggle to achieve. Our results demonstrate the complexity of physical processes required to capture realistic dwarf galaxy properties, making tangible predictions for integral field unit surveys, radio synchrotron emission, and for galaxy and multiphase interstellar medium properties that JWST will probe.

1 INTRODUCTION

Dwarf galaxies are intriguing dark-matter-dominated systems, subject to some of the most persistent controversies in the theory of galaxy formation. Most notable are the well-known missing satellites, cusp-core, and too-big-to-fail problems (Garrison-Kimmel et al. 2014; Oñorbe et al. 2015; Sawala et al. 2016; Bullock & Boylan-Kolchin 2017), which have prompted many studies challenging the standard cosmological Lambda cold dark matter model (Spergel & Steinhardt 2000; Libeskind et al. 2013). Some of these problems likely arise due to inaccurate treatment of the complex baryonic physics. ‘Missing’ satellites have been largely accounted for by considering a combination of their disruption due to photoionization (Efstathiou 1992; Benson et al. 2002; Bose, Deason & Frenk 2018), photoevaporation (Bullock, Kravtsov & Weinberg 2000) linked to a filtering mass (Gnedin 2000; Okamoto, Gao & Theuns 2008), photoheating starvation (Hoeft et al. 2006; Katz et al. 2020), and the impact of stellar feedback (Dekel & Woo 2003; Garrison-Kimmel et al. 2019). Furthermore, some problems may be alleviated by observational advances, which allow the detection of fainter previously ‘missing’ systems (e.g. Belokurov et al. 2008; Koposov et al. 2015) and by accounting for observational biases (Oman et al. 2016). While dwarf galaxies are the most abundant type of galaxy, their low masses and luminosities reduce their detectability and study to relatively modest redshifts (Tolstoy, Hill & Tosi 2009). Detailed observations are mostly limited to those found within the Local Group or not far beyond (e.g. Kirby et al. 2017; Read et al. 2017; Wheeler et al. 2017). However, excitingly well-resolved observations of dwarf galaxies at high redshift will very soon become available with JWST (Jeon & Bromm 2019).

Similarly, as for the missing satellites problem, baryonic physics may help to resolve the cusp-core controversy as well, with this line of research gaining traction over the years. The gravitational response of dwarf galaxy haloes to explosive supernova (SN) events is capable of producing cores in the central region of dwarf density profiles (Navarro, Eke & Frenk 1996). This process is particularly efficient when driven by resonantly cyclic SN bursts (Governato et al. 2012; Pontzen & Governato 2012), where systems featuring younger stellar populations could produce more cored profiles than those in which star formation ceased rapidly early on (Read, Walker & Steger 2019). In fact, because of their shallow potential wells and low masses, dwarf galaxies are particularly sensitive probes of baryonic physics (Bullock & Boylan-Kolchin 2017). Furthermore, due to their rapid quenching at high redshifts (Frebel, Simon & Kirby 2014; Rey et al. 2020; Chiti et al. 2021) and during the epoch of reionization (Barkana & Loeb 1999), dwarf galaxies are unique archaeological probes of physical processes shaping the formation of galaxies in the early Universe.

Various galaxy formation simulations generated with different codes and sub-grid prescriptions for star formation and stellar feedback have been able to reproduce the central scaling relations between the stellar component and the haloes of galaxies predicted by abundance matching (see e.g. Dubois et al. 2014; Vogelsberger et al. 2014; Crain et al. 2015; Oñorbe et al. 2015; Fattahi et al. 2016; Henden et al. 2018; Hopkins et al. 2018; Pillepich et al. 2018b). These scaling relations are rather uncertain in the dwarf regime, where very high-resolution simulations are needed to resolve the internal processes and the interstellar medium (ISM) of dwarfs (see e.g. Smith, Sijacki & Shen 2019; Wheeler et al. 2019; Agertz et al. 2020; Gutcke et al. 2021). While SN feedback has traditionally been invoked as the dominant physical process regulating the properties of dwarf galaxies (White & Rees 1978; Dekel & Silk 1986; White & Frenk 1991; Efstathiou 2000), it is now understood that other physical processes beyond SN feedback are also required to reproduce realistic dwarf properties (e.g. Rosdahl et al. 2018; Smith et al. 2019). Early stellar feedback in the form of e.g. winds from massive stars or stellar radiation provides a mechanism to rapidly suppress star formation, with the impact of photoheating shown particularly important within the mass regime of dwarf galaxies (Rosdahl et al. 2015; Emerick, Bryan & Mac Low 2018). Alternatively, warm or self-interacting dark matter have been explored as a means to alleviate the cusp-core problem (Villaescusa-Navarro & Dalal 2011; Lovell et al. 2014; Vogelsberger et al. 2016). Additionally, recent observational evidence supporting the presence of active galactic nuclei (AGNs) in dwarf galaxies has sparked interest in whether AGN feedback in these galaxies may have a significant effect on their evolution (e.g. Habouzit, Volonteri & Dubois 2017; Dashyan et al. 2018; Koudmani et al. 2019; Koudmani, Henden & Sijacki 2021; Koudmani, Sijacki & Smith 2022). However, well-known baryonic physics, such as magnetic fields, stellar radiation, and ∼GeV-energy cosmic rays, still need to be studied in detail in full cosmological simulations of dwarf galaxies.

Estimates of the magnetic energy in galaxies reveal the importance of magnetic fields for the ISM, with observations suggesting equipartition between the thermal, turbulent, and magnetic pressure components (Basu & Roy 2013; Beck 2015). Magnetic fields with strengths of |$\gtrsim \!\!\, \mu \mathrm{G}$| have been detected in dwarf galaxies (Chyzy et al. 2011). Such magnetic fields can affect the galaxies in multiple ways. At sub-ISM scales, magnetic fields are well known to be an important factor in regulating the dynamics of ISM turbulence (Padoan & Nordlund 2011). They affect thermal instabilities and thus regulate the different phases of the ISM (Iffrig & Hennebelle 2017; Körtgen et al. 2019) as well as influencing gas fragmentation (Inoue & Yoshida 2019). On galactic scales, magnetic fields have the potential to impact galactic outflows (Grønnow, Tepper-García & Bland-Hawthorn 2018; Steinwandel et al. 2019), halo gas mixing (Van De Voort et al. 2021), global galactic properties (Pillepich et al. 2018a; Martin-Alvarez et al. 2020), and possibly play a role during galaxy merging (Whittingham et al. 2021). Magnetic fields are none the less not expected to have a major effect on fundamental properties such as the final stellar mass of galaxies (Pakmor et al. 2017; Su et al. 2017; Martin-Alvarez et al. 2020).

Even though SN explosions may be the most obvious form of stellar feedback, stellar radiation is widely recognized as an important agent in galaxy evolution. Stellar radiation is particularly fundamental for low-mass, late-type dwarf galaxies (Rosdahl et al. 2015). Accounting for radiative feedback leads to a gentler self-regulation of dwarf galaxies and reduces the importance of explosive SN events (Agertz et al. 2020). Katz et al. (2020) showed how photoheating by stellar radiation can lead to high-redshift quenching of dwarf galaxies by evaporating the filaments responsible for their gas supply. Stellar radiation also augments the impact of individual SN events by pre-processing gas parcels where these explosions will take place (Geen et al. 2015). Through this effect, stellar radiation has been argued to support the driving of galactic outflows (Emerick et al. 2018). However, other studies have claimed opposite effects (Rosdahl et al. 2015; Smith et al. 2021), as photoionization feedback may disrupt local star-forming clouds much more rapidly than the first SN exploding, hence leading to a less bursty and less effective SN feedback.

Finally, cosmic rays have gained a lot of attention in recent years due to their ability to efficiently drive continuous and colder galactic outflows, especially when compared with winds driven solely by SNe (e.g. Booth et al. 2013; Salem & Bryan 2014; Girichidis et al. 2018; Samui, Subramanian & Srianand 2018; Buck et al. 2020; Dashyan & Dubois 2020; Farcy et al. 2022; Rodríguez Montero et al. 2023). Cosmic rays are able to reduce star formation rates in isolated (Pakmor et al. 2016; Dashyan & Dubois 2020) and cosmological (Jubelgas et al. 2008; Hopkins et al. 2020) galaxy formation simulations. Such reduction is found for most galaxy masses in isolated simulations. In a cosmological context, Hopkins et al. (2020) found this effect only for their larger galaxy masses, also varying based on the selected cosmic ray diffusion coefficient κ_∥. This disparity points towards some dependence on the employed physics of cosmic ray transport (diffusion, streaming, etc.), as well as on the κ_∥ coefficient, which remains poorly constrained. Depending on their implementation, cosmic rays may also alter the appearance of galaxies (Buck et al. 2020) and their ISM (Commerçon, Marcowith & Dubois 2019; Dashyan & Dubois 2020; Nuñez-Castiñeyra et al. 2022). When converting a small fraction of SN energy into cosmic rays (typically |$\sim \!\!10~{{\ \rm per\ cent}}$|⁠; e.g. Dashyan & Dubois 2020; Hopkins et al. 2020), this non-thermal energy has the potential to enhance the deposition of momentum by SNe (Diesing & Caprioli 2018; Rodríguez Montero et al. 2022). Most importantly, once they escape SN remnants, cosmic rays can establish a ‘smooth’ pressure gradient beyond galactic scales which is believed to accelerate gas at the edges of galactic discs (Hopkins et al. 2021).

In addition to their individual effects, due to the highly non-linear nature of galaxy formation physics, these different processes have a complex interplay when combined in the same simulation. For example, in the presence of magnetic fields the efficiency of star formation ramps up in molecular clouds (Federrath & Klessen 2012; Zamora-Avilés et al. 2018), but these clouds will then be rapidly dissipated by the radiation produced by the newly formed stars (Murray, Quataert & Thompson 2010). Likewise, radiation has the potential to puff up gas discs in galaxies (Roškar et al. 2014), whereas cosmic rays are capable of accelerating the diffuse gas located at high altitudes above discs (Breitschwerdt, McKenzie & Völk 1991; Girichidis et al. 2016). In combination, these effects may thus increase the amount of gas expelled by cosmic rays.

In this work, we aim to shed light on the role played by each of these additional physical processes, as well as on their interplay in the formation of dwarf galaxies. For this purpose, we perform a suite of cosmological zoom-in simulations of galaxy formation that, starting with our fiducial star formation and stellar feedback models, gradually includes additional and more sophisticated physics up to ‘full-physics’ simulations featuring magneto-thermo-turbulent (MTT) star formation and mechanical SN feedback with magnetic fields, on-the-fly radiative transfer, and cosmic rays. We investigate the evolution and resulting properties of the simulated dwarf galaxies and compare them with observations of similar galaxies in our local Universe to understand which physical processes are more likely to shape dwarf properties. To produce an agnostic and objective assessment of our different models, and to be able to better discriminate them against observations and empirical relations, we base the selection of all the model parameters in Pandora exclusively on physical considerations (see Section 2, with a summary of our simulations in Table 1), without any aim to match any specific observables. That is, we do not attempt to tune our model parameters and wherever a specific model fails to match expected quantities, we focus on investigating the physical mechanisms at play causing this result.

This work is organized as follows. The numerical methodology to generate and evolve our simulations is described in Section 2, with further descriptions of our magnetic fields (Section 2.2), radiative transfer (Section 2.3), and cosmic rays (Section 2.4) models. Section 2.5 describes our simulations suite and Section 2.6 discusses which observed dwarf galaxies provide the best comparison to Pandora. Section 3 reports our main results, commencing with the evolution of the stellar mass–halo mass relation (Section 3.2), followed by the study of stellar and gas morphology of our galaxy (Section 3.3), its resolved and integrated kinematics (Section 3.4), the colour–magnitude relation (Section 3.5), magnetic field and synchrotron synthetic observations (Section 3.6), and finally by analysing the impact of different physical processes on its dark matter distribution (Section 3.7). Finally, we conclude with a summary of our work and its main conclusions in Section 4.

2 METHODS

All galaxy formation simulations studied in this work are new cosmological zoom-in simulations generated with our own modified version of the publicly available ramses code (Teyssier 2002). ramses models the dark matter and stellar components as an ensemble of collisionless particles. These are coupled to each other and the baryonic gas through the gravity solver of the code. The evolution of the gaseous component is instead solved on an adaptive mesh refinement (AMR) grid that is recursively refined in regions of interest. We extend our version of ramses to simultaneously and self-consistently model radiative transfer (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), constrained transport (CT) magnetohydrodynamics (MHD; Fromang, Hennebelle & Teyssier 2006; Teyssier, Fromang & Dormy 2006), and cosmic rays (Dubois & Commerçon 2016; Dubois et al. 2019). Each of these physical components and the configurations we employ for them are described in Sections 2.2 (MHD), 2.3 (radiative transfer), and 2.4 (cosmic rays).

We re-simulate one of the dwarf galaxies studied by Smith et al. (2019), labelled as Dwarf1 in their work. Note that our simulations are generated with the ramses code instead of arepo, and that they feature different resolutions as well as physical models and implementations. An initial analysis of our set-up showed good agreement of our halo and galaxy masses with those reported by Smith et al. (2019). As we employ different sub-grid prescriptions for e.g. star formation and stellar feedback, we expect, however, multiple differences to emerge between our results and theirs. Our initial conditions are for a cubic box of ∼14.73 comoving Mpc (cMpc) per side, initialized at z = 127 and discretized with a uniform grid of 256³ cells. At the centre of our computational domain a dwarf galaxy forms in a halo of virial mass |$M_\text{vir} (z = 0) \sim 10^{10} \, \mathrm{M_\odot }$|⁠. We follow the formation of this dwarf galaxy using a convex hull zoom region. At z = 127, this region is approximately 2.5 cMpc across, and it is evolved using a passive refinement scalar advected with the fluid and the high-resolution dark matter particles. The dark matter and stellar particle mass resolutions in this region are |$m_\text{DM} = 1.5 \times10^{3} \, \mathrm{M_\odot }$| and |$m_\text{*} = 400 \, \mathrm{M_\odot }$|⁠, respectively. In the zoom region, we allow the AMR to refine the grid down to a resolution of Δx ∼ 7 physical pc (or 3.5 pc radius/half-cell size). Cells are marked for refinement when their total contained mass |$(\Omega _m/\Omega _b)\, m_\text{baryons} + m_\text{DM}$| surpasses 8 × m_DM or when their size is Δx > λ_J/4, with λ_J being the local Jeans length. All our simulations employ the Ade et al. (2016) cosmology. Furthermore, we impose an initial metallicity floor of |$10^{-4} \, \mathrm{Z_\odot }$|⁠, corresponding to the critical metallicity required for gas fragmentation to allow PopII stellar clusters to be formed (⁠|$\, \mathrm{Z_\odot }= 0.012$|⁠; Schneider et al. 2012). Due to their computational cost, we evolve the majority of our simulations only to z = 3.5, when the mass of the studied halo is |$M_\text{vir} (z = 3.5) \sim 5 \times10^{9} \, \mathrm{M_\odot }$|⁠. We only continue some of our hydrodynamical runs (HD, HD+Boost, NoFband NoFb + NoZ; listed in Table 1) down to z = 0.5.

Unless indicated in Section 2.5, all our simulations include metal cooling above and below a threshold of 10⁴ K interpolating cloudy cooling tables (Ferland et al. 1998) and following Rosen & Bregman (1995), respectively. We model the effects of ionizing radiation as an homogeneous ultraviolet (UV) background according to Haardt & Madau (1996), activated at z < 9. We always assume the baryonic gas to be monatomic and ideal (i.e. with a specific heat ratio γ = 5/3).

To determine the position of the main galaxy and its halo, we use halomaker (Tweed et al. 2009), employing a shrinking spheres algorithm (Power et al. 2003) to attain a higher centring accuracy. The application of the halo finder to the dark matter determines the location and properties of the studied dwarf galaxy halo. Whenever required, we obtain the centre of the galaxy and its angular momentum by applying halomaker to the baryons (gas and stars), with the centre of the galaxy being located inside the central region of the halo (r < 0.2 r_halo).

2.1 Star formation and stellar feedback

The process of star formation in our simulations is modelled by converting a fraction of the gas mass contained within a given cell into a new stellar particle. For this, we investigate two different prescriptions. The first is our fiducial prescription: an MTT star formation model, presented in Kimm et al. (2017), and Martin-Alvarez et al. (2020) in its MHD extension. In this model, we only allow cells at the highest level of refinement to form stars (Rasera & Teyssier 2006) as long as they fulfil the condition that the gravitational pull due to their gas content is higher than the support provided by the combined turbulent, thermal, and magnetic pressure. Star-forming cells convert gas into stars following a Schmidt law (Schmidt 1959),

with star formation efficiency, ϵ_ff, and free-fall time, t_ff. In this model, ϵ_ff is a local parameter determined by the MTT properties of the local gas cells. For its computation, we follow the multifree fall version of the Padoan & Nordlund (2011) model as described by Federrath & Klessen (2012). For further details, we refer to appendix B of Martin-Alvarez et al. (2020).

A small subset of our simulations assume an alternative star formation model, the more commonly used gas density threshold. This model assumes a fixed star formation efficiency of ϵ_ff = 0.015 and allows star formation to proceed according to equation (1) whenever a cell at the highest refinement level (Rasera & Teyssier 2006) has a gas density that exceeds |$\rho _g \gt \rho _\text{th} = 10\,\, \mathrm{m_H}\, \mathrm{cm}^{-3}$|⁠, with m_H the hydrogen atom mass. These values follow Smith et al. (2019), and are selected to produce a star formation model analogous to that by Krumholz & Tan (2006). These simulations using density threshold star formation are the most similar within this study to the simulations studied by Smith et al. (2019) in terms of configuration, but the caveats described above remain.

Whenever radiative transfer is included in our simulations, stellar particles emit radiation into specific energy bins. Radiative feedback from stars is described in more detail in Section 2.3. The majority of our simulations also feature SN feedback, employing the mechanical SN feedback prescription (Mech) by Kimm & Cen (2014) and Kimm et al. (2015). To determine when SN events occur, each stellar particle has its initial mass function (IMF) stochastically sampled during the first 50 Myr after its formation. Each stellar particle undergoing an SN event injects into its hosting and neighbouring cells mass, momentum, and energy, with a specific energy of ε_SN = E_SN/M_SN (except for the boosted feedback simulations – as indicated in Section 2.5), where |$E_\text{SN} = 10^{51} \, \mathrm{erg}$| and |$M_\text{SN} = 10\, \, \mathrm{M_\odot }$|⁠. We assume a Kroupa IMF (Kroupa 2001), returning a fraction η_SN = 0.213 of the total SN exploding mass to the ISM gas. A further fraction η_metals = 0.075 of this total mass corresponds to the gas returned as metal mass. In some of our simulations our SN feedback also injects magnetic and cosmic ray energy back to the ISM. A detailed description of these injections by the SN feedback are provided in Sections 2.2 and 2.4, respectively.

2.2 Magnetohydrodynamics

We model magnetic fields and their coupling to the gas fluid with the CT ideal MHD implementation of ramses by Fromang et al. (2006) and Teyssier et al. (2006). This implementation models the magnetic field |$\vec{\mathbf{ B}}$| as a face-centred quantity in each cell. The algorithm ensures that the divergence of the magnetic field fulfils the solenoidal constraint (⁠|$\vec{\nabla } \cdot \vec{\mathbf{ B}} = 0$|⁠) exactly to the numerical precision. This guarantees the absence of spurious MHD modifications and the preservation of conserved quantities, which is not ensured with other methods¹ (Tóth 2000). The time evolution of the magnetic field is computed solving the induction equation.

As magnetic diffusivity in most astrophysical environments is negligible, we set this quantity to zero in our simulations. This implies that all magnetic diffusive effects in our simulations will result from the numerical magnetic diffusivity emerging when the domain is discretized into a finite grid.

In the absence of battery terms in the ideal MHD induction equation (such as the implementation of a Biermann battery, e.g. Attia et al. 2021), a magnetic seed has to be introduced to obtain any |$\vec{B} \ne 0$| field. We investigate two different approaches: (a) an ab-initio magnetic field, and (b) an SN-injected magnetic field. In the first approach, we permeate the simulated domain with a uniform magnetic field along its z-axis with comoving strength B₀. This method is the most commonly employed in MHD simulations (e.g. Pakmor et al. 2016; Marinacci et al. 2018; Martin-Alvarez et al. 2018), and can be interpreted as a magnetic field of primordial origin coherent on large scales. Galaxy formation simulations seeded with sufficiently small B₀ retain negligible or no memory of the initial seed (Marinacci et al. 2015), with only the strongest primordial magnetic fields being able to affect the properties of galaxies (B₀ > 10⁻¹² G; Martin-Alvarez et al. 2020). The second method of seeding injects small-scale circular loops of magnetic field around SNe as the SN events take place. This guarantees |$\vec{\nabla } \cdot \vec{B} = 0$| when magnetizing the SN ejecta in our mechanical feedback. As the ejecta expand, they advect the injected magnetic field to larger scales in the ISM. Each SN explosion is assumed to inject |$E_\text{inj,mag} = 0.01 E_\text{SN} \sim 10^{49} \, \mathrm{erg}$|⁠. This corresponds to a magnetic field strength of ≳10⁻⁵ G when injected at scales of ∼10 pc, in reasonable agreement with the observed high magnetization of SN remnants (Parizot et al. 2006). This magnetic injection model is capable of reproducing the magnetic fields observed in galaxies (Martin-Alvarez et al. 2021). Additional details of the magnetic field SN injection implementation can be found in appendix A of Martin-Alvarez et al. (2021).

2.3 Radiative transfer

The implementation for ionizing emissivity in our simulations is similar to the one in the sphinx simulations (Rosdahl et al. 2018, 2022), which provide a good match to observational constraints on the reionization history of our Universe. We employ the ramses-rt implementation by Rosdahl et al. (2013) and Rosdahl & Teyssier (2015). Radiative transfer is particularly sensitive to the distribution of multiphase gas within the ISM. Kimm & Cen (2014) find that a resolution of ∼5 pc is required to have a converged escape of ionizing radiation into the ISM. Due to our similarly high spatial resolutions (Δx ∼ 7 pc) we expect our radiative transfer approach within the ISM as well as its escape from the galaxy to be reasonably well resolved.

We separate the radiation into three spectral bins: H i (13.6 eV ≤ ϵ_photon < 24.59 eV), He i (24.59 eV ≤ ϵ_photon < 54.42 eV), and He ii (54.42 eV ≤ ϵ_photon). We allow the radiation solver to subcycle over the hydrodynamical solver up to a maximum of 500 steps, and assume a reduced speed of light |$0.01\, c$|⁠. This is sufficient for modelling the propagation of ionization fronts through the ISM of galaxies. In our radiative transfer simulations, stellar particles are the only source of ionizing radiation. Each stellar particle radiates energy into its hosting cell with a spectral energy distribution given by the bpassv2.0 model (Eldridge, Izzard & Tout 2008; Stanway, Eldridge & Becker 2016) according to particle mass, metallicity, and age.

2.4 Cosmic rays

Our simulations including cosmic rays use the cosmic ray implementation in ramses by Dubois & Commerçon (2016) and Dubois et al. (2019). This implementation assumes cosmic rays to behave as a fluid characterized by an energy density which is evolved by an implicit solver. Unless explicitly indicated for a given model, our simulations with cosmic rays account both for anisotropic diffusion and streaming of cosmic rays. We assume the cosmic ray diffusion coefficient to be |$\kappa _\parallel = 3 \times 10^{28} \, \mathrm{cm}^2 \text{s}^{-1}$|⁠. This value has been found to be consistent with observations of γ-rays generated through cosmic ray hadronic losses (Ackermann et al. 2012; Salem, Bryan & Corlies 2016; Pfrommer et al. 2017b). This value is also in agreement with estimates for the isotropic coefficient in the Milky Way (Trotta et al. 2011; Cummings et al. 2016). We assume that the streaming velocity is equal to the Alfvén velocity. In addition to streaming and adiabatic energy losses, the cosmic ray implementation in ramses also accounts for hadronic and Coulumb cooling of cosmic rays (Guo & Oh 2008). Cosmic rays are assumed to be generated by SN explosions. In our simulations with cosmic rays, each SN event injects a fraction of its total energy as cosmic rays to its hosting cell, where this energy follows |$E_\text{CR} = f_\text{CR} E_\text{SN} = 10^{50} \, \mathrm{erg}$| and is extracted from the standard thermal injection. The selected fraction f_CR = 0.1 is in agreement with observational estimates for cosmic ray injection by SN remnants (Morlino & Caprioli 2012; Helder et al. 2013). Finally, we note that the selected values for κ_∥ and f_CR have been frequently studied in the literature (e.g. Wiener, Pfrommer & Peng Oh 2017; Pfrommer et al. 2017a; Butsky & Quinn 2018; Dashyan & Dubois 2020).

2.5 The simulation suite

Our suite of simulations builds up from a fiducial simulation employing hydrodynamics (HD) with star formation and stellar feedback physics up to a full-physics simulation with radiative transfer, cosmic rays, and magnetic fields. A compilation of all the simulations studied in this work is presented in Table 1. In this table, and throughout the rest of the manuscript, we label our simulations according to the physical processes included: magnetohydrodynamics (MHD), radiative transfer (RT), and cosmic rays (CR). According to the source of magnetic fields, we have simulations with a primordial magnetic field of a moderately high strength (B₀ = 3 × 10⁻¹³ G) simply labelled MHD; a strong primordial magnetic field (B₀ = 3 × 10⁻¹¹ G) labelled sMHD, and magnetized SN seeding (note that these simulations still include a negligible primordial field, B₀ = 3 × 10⁻²⁰ G) labelled iMHD. We do not need to discriminate different radiative transfer configurations, as all our simulations share the same radiative transfer implementation (Section 2.3). For cosmic rays, our simulations contain both cosmic ray streaming and diffusion (simply identified as CR) or exclusively diffusion with deactivated streaming (in the simulation labelled RTnsCRiMHD; identified by nsCR). Finally, the names of some of our simulations also contain information of their treatment for star formation and stellar feedback. Unless indicated otherwise, all our simulations include our fiducial MTT star formation implementation. The simulations using the standard density threshold star formation criterion are identified by adding thSf to their name. We also have some hydrodynamical simulations where SN feedback events do not return any energy to the ISM. These have HD replaced in their name by NoFb. Amongst these simulations, we distinguish those where SN feedback events do not take place at all (NoZ) from those where the SN feedback only returns mass (including metals) to the ISM. Finally, we include two simulations where we boost the specific energy of the mechanical feedback by a factor of 4 (Boost, with |$\varepsilon _\text{SN} = 2\, E_\text{SN} /\, 0.5\, M_\text{SN}$|⁠). This calibration is inspired by the sphinx simulations (Rosdahl et al. 2017), which found that such a boost factor for the energy injection by SN feedback² reproduces well the observations by Read et al. (2017) and the stellar mass to halo mass ratio predicted by abundance matching (Behroozi, Wechsler & Conroy 2013) in their low-mass galaxies at z = 6. With this calibration of the stellar mass–halo mass relation, sphinx obtains a remarkable match to observations of the reionization history.

In summary, our simulations can be separated into four groups:

• HD simulations, exploring different prescriptions for star formation and parameters for stellar feedback.

• MHD simulations, exploring different sources of magnetic fields.

• RT/RTMHD simulations, exploring the effects of stellar radiation and its combination with magnetic fields.

• CRMHD/RTCRMHD simulations, exploring the impact of cosmic rays and their combination with radiative transfer and magnetic fields.

Out of these, we focus most of our investigation on seven simulations: NoFb, HD, HD+Boost, RT, RTiMHD, RTnsCRiMHD, and RTCRiMHD. This subsample builds up from no SN feedback up to ‘full-physics’ simulations. We place special emphasis on three simulations that serve us to explore how standard modern simulations fare against the upcoming more complete models, and to study the importance of additional, well-known baryonic physics. These are stellar feedback HD, calibrated stellar feedback HD + Boost, and our ‘full-physics’ simulation RTCRiMHD. A detailed view of the full-physics dwarf galaxy at z = 3.5 is shown in Fig. 1. Gas density maps and synthetic optical views for other representative simulations are presented in Fig. 2.

2.6 Comparable dwarf galaxies in observations

In order to aid in reviewing our results and providing some context for galaxy properties and relations, we include observational data in various figures. The observational data included feature a broad population of low-mass and dwarf galaxies. While our simulated models may resemble the properties of many of these different systems across the multiple relations, we focus the comparison on those observed galaxies that have similar masses and are found in comparable environments. Our simulated dwarf galaxy is an isolated system that forms in a relatively quiet environment. Consequently, we favour comparison with the most isolated dwarf galaxies in the Local Group, with similar halo masses (⁠|$M_\text{vir} (z = 0) \sim 10^{10} \, \mathrm{M_\odot }$|⁠) and within a comparable stellar mass range (⁠|$M_*\sim 10^{6}\!-\!5 \times10^{7} \, \mathrm{M_\odot }$|⁠). The best analogues in the observational data are the dwarf galaxies WLM, LeoA, VV124, and SagDIG. On the other hand, dwarf galaxies that have lower halo and stellar masses at z = 0 are candidates for reionization quenching of star formation. Some of these systems, such as LeoT (⁠|$M_\text{halo} (z = 0) \sim 5 \times 10^8 \, \mathrm{M_\odot }$|⁠) have only recently re-ignited their star formation (Irwin et al. 2007; Rey et al. 2020). As these correspond to a different category of dwarf galaxies undergoing different evolutionary processes, we avoid comparing our models with them.

3 RESULTS

3.1 Dwarf galaxy morphology with different ISM physics

We evolve all our simulations from z = 127 to z = 3.5. In Fig. 2 we present face-on views of the dwarf galaxy at this lower redshift for our most representative simulations. The panels are projections of |$(5 \, \mathrm{kpc})^3$| cubic boxes centred on the galaxy with their line of sight oriented along the total baryonic angular momentum of the system. For each simulation, the two rows represent gas density (top) and a synthetic optical RGB image for the SDSS [u,g,r] filters in rest frame (bottom). We assume each stellar particle to be a single stellar population, with emission as predicted by Bruzual & Charlot (2003). All emission is integrated along the line of sight, accounting for dust obscuration by assuming a 3D absorption screen with an exponential extinction law of R_v ∼ 3.1. We model the dust content of each cell employing a simple metal-to-dust ratio of 0.4, found to be a reasonable approximation of a full post-processing radiative transfer treatment for this type of synthetic images by Kaviraj et al. (2016).

The physical processes included in the simulations and their different configurations significantly affect the properties and appearance of the galaxy. Focusing first on the simulations with no stellar feedback (top rows, three leftmost columns), they exhibit, unsurprisingly, extended gaseous discs and massive stellar components with an exceptionally pronounced stellar density peak (bulge) at their centre. Enabling the return of metals to the ISM has a substantial impact on the morphology of gaseous discs. They develop spiral structures and overdense ‘knots’, which correspond to the clearly visible stellar clusters in the optical emission, particularly prevalent in the simulation with the density threshold star formation model. By construction, the density threshold model allows for a greater fraction of cold, dense gas to become star-forming, whereas with the MTT model, gas must reach significantly higher densities before becoming star-forming. This leads to the density-threshold models converting a larger fraction of their gas into stars, consequently reaching higher stellar masses.

Including SN feedback (top rows, two rightmost columns and central rows, two leftmost columns) modifies the appearance of the dwarf galaxy considerably, becoming an irregular system with no clear disc. Furthermore, both the stellar and gas mass content is reduced, with the optical image showing a much fainter luminosity (note that we decreased the luminosity of the NoFb models in this figure by 1 dex to avoid image saturation). The galaxies produced with the density-threshold star formation model have notably more diffuse distribution of gas and stars, due to a more efficient impact from SN feedback. This is the result of gas being allowed to be converted into stars at earlier times in the density threshold models, as well as due to the formation of additional stars in diffuse regions. Hence, in the density threshold model, more SN events occur earlier than in the MTT model, and more events take place in regions of lower density, where SN feedback is more impactful. The simulations using the MTT model for star formation have a clumpier appearance, resembling local Universe dwarf irregulars as well as the more massive high-redshift analogues observed by GOODS and GEMS (Elmegreen et al. 2009). Finally, in addition to featuring even lower stellar masses and fainter luminosities, the boosted SN feedback simulations have depleted the system of the densest gas.

The remaining panels in Fig. 2 compare various of our simulations with additional physical processes.

• Magnetic fields: The simulation with a strong primordial magnetic field (central rows, central column, sMHD) shows a higher central concentration of the baryonic component (Martin-Alvarez et al. 2020) than the HD simulation. The presence of magnetic fields further produces a smoother gas structure of the ISM (Körtgen et al. 2019; Martin-Alvarez et al. 2020), also apparent in the iMHD simulation (central rows, fourth column). However, note that this latter simulation has had its star formation quenched for a considerable time, which is the aspect dominating its appearance.

• Radiative transfer: The addition of radiative transfer (central rows, last column) generally leads to a higher fraction of gas in the circumgalactic medium (CGM). We will show later how radiative transfer leads to a higher gas mass being retained in the galaxy, as the photoheating due to the stellar radiation prevents a fraction of gas to form stars and leads to gentler SN feedback. The stellar component is similar to that in the HD simulation. However, the spatial distribution of the stars is more extended and diffuse as well as bluer due to some recent star formation. The combination of radiative transfer and magnetic fields (bottom rows, first two columns) is particularly interesting. We will explore the interplay of these two physical processes in more detail next. For now, we note two important aspects in these simulations. First, the appearance of the gas and stars is somewhere in between the MHD-only and radiative transfer-only runs. Second, it is worth noting that the simulation with the strong primordial magnetic field (RTsMHD; bottom rows, first column) shows a dense concentration of stars, which we will later show corresponds to an overmassive and overdense stellar component.

• Cosmic rays: Finally, our cosmic ray simulations are shown in the last three columns of the bottom rows. The first of these three is the non-RT case, where we find an appearance similar to the iMHD simulation, but with a lower stellar mass and a more centrally concentrated gas-rich environment. We will later show that while both the iMHD and CRiMHD simulations follow similar star formation histories, the simulation with cosmic rays is nevertheless able to more efficiently eject gas and prevent star formation. Our most interesting simulations are the full-physics RTnsCRiMHD and RTCRiMHD simulations. These show a more significant gas reservoir as well as a denser CGM. The bottom rightmost panel in Fig. 1 illustrates how most of the gas in RTCRiMHD is outflowing.

3.2 The stellar mass to halo mass relation

One of the most studied relations in galaxy formation is the stellar mass M_* to halo mass M_halo relation. A physically insightful form of this relation employs the baryonic conversion ratio, M_*/(f_bM_halo), versus M_halo, where f_b = Ω_b/Ω_m is the baryon fraction, and Ω_b and Ω_m are the cosmological baryon and matter density parameters, respectively. This form assigns to a specific halo mass an efficiency for the conversion of its baryonic mass into stellar mass. In Fig. 3 we explore the evolution of the baryonic conversion efficiency for each simulation of our dwarf galaxy.

In order to contextualize the differences across simulations, we include shaded bands corresponding to various abundance matching relations. As many of these relations only reach halo mass limits as low as |$M_\text{halo} \sim 10^{10} \, \mathrm{M_\odot }$| we often resort to extrapolating them into the studied regime. Consequently, such extrapolations have to be examined with care. For example, the extrapolation of the Behroozi et al. (2013) and Read et al. (2017) relations predict considerably higher stellar masses than a variety of other models such as Garrison-Kimmel et al. (2014) or Moster, Naab & White (2018). We include models with stellar masses in between these two regimes, by showing the predictions by Moster et al. (2013) and Jethwa et al. (2018). Models predicting stellar masses below the presented range (e.g. Garrison-Kimmel et al. 2014; Moster et al. 2018) are not shown in Fig. 3 in order to facilitate visual examination of the differences across our models.

Furthermore, we compare our galaxy with the population of galaxies of similar halo masses (⁠|$M_\text{halo} (z = 6) \sim 1.5 \times10^{9} \, \mathrm{M_\odot }$|⁠) in the sphinx simulation (see cyan data point). sphinx is a high-resolution galaxy formation simulation including on-the-fly radiative transfer that reproduces well the reionization history of our Universe. This comparison at z = 6 may provide us with some insight regarding how our different models compare to those that reproduce observed reionization constraints. Finally, we include in the same figure observational data from Read et al. (2017) as black points for LeoT, Carina, and various isolated LG dwarf galaxies. We first discuss the general evolution of Pandora dwarfs in our no-feedback and SN-only simulations, and then move on to the discussion of the simulations with additional physics.

In the absence of stellar feedback, we find baryonic conversion ratios to be too large (>0.1). Out of the two star formation prescriptions studied, the MTT model yields somewhat lower albeit still unrealistically large values. The return of metals to the ISM by stars seems to only significantly increase the stellar mass of the MTT star formation prescription model. As expected, accounting for SN feedback dramatically reduces the stellar mass of our dwarf galaxy at all redshifts explored. While the effect is comparatively minor shortly after the onset of star formation (z ∼ 10), SN feedback leads to a large reduction in stellar mass at lower redshifts. In fact, at z ∼ 4, the HD simulation reaches a conversion ratio of |$\sim \!\!2~{{\ \rm per\ cent}}$|⁠, comparable with the extrapolation of the Behroozi et al. (2013) model to this halo mass. The simulations with standard SN feedback show no large differences between the density threshold and our MTT star formation implementations in terms of the baryonic conversion efficiency. With the classical density threshold star formation implementation, SN feedback has a slightly higher impact. We attribute this to SN events taking place at lower densities and thus providing a more efficient deposition of momentum (Iffrig & Hennebelle 2015). In practice, star formation in the MTT simulations typically occurs at densities |$\rho _\text{gas} \gtrsim 10^{-21} \, \, \mathrm{g}\, \, \mathrm{cm}^{-3}$|⁠, whereas in the density threshold case, we allow star formation only when |$\rho _\text{gas} \gt \rho _\text{th} \sim 10^{-23} \, \, \mathrm{g}\, \, \mathrm{cm}^{-3}$|⁠. Furthermore, in the density threshold simulations, the earlier onset of star formation suppresses the initial collapse of the gas. However, most gas is re-accreted later on and the simulation reaches baryonic conversion efficiency comparable to HD.

A common way of further reducing the stellar mass of simulated galaxies is to boost the strength of SN feedback (motivated by considering that other sources of stellar feedback such as binary stars, hypernovae, top heavy IMF, etc. have been neglected). This is what we do in our HD + Boost simulation. This simulation has, as expected, a lower baryonic conversion efficiency that continues to decrease until it reaches |$\sim 0.4~{{\ \rm per\ cent}}$| by z = 3.5. Here, the expulsion of gas is much more efficient and the galaxy maintains a significantly lower stellar mass. We further note that the simulation with boosted SN feedback and density threshold for star formation has the most extreme effect in reducing the stellar mass, due to the powerful SNe exploding in relatively low density environments, to a level that is likely unrealistic. Unfortunately, the high computational cost of the simulations including more complex physical processes impedes evolving all our simulations to z ∼ 0. However, we extend some of our hydrodynamical simulations to z = 0.5 to demonstrate that our dwarf galaxy and its halo do not undergo significant mass evolution after z ∼ 4. In particular, the halo mass evolves from |$M_\text{halo} (z = 3.5) \sim 4 \times 10^9 \, \mathrm{M_\odot }$| to |$M_\text{halo} (z = 0.5) \sim 8 \times 10^9 \, \mathrm{M_\odot }$|⁠, but most importantly, the stellar mass is barely changed, evolving from |$M_*(z = 3.5) \sim 1.0 \times 10^7 \, \mathrm{M_\odot }$| to |$M_*(z = 0.5) \sim 1.3 \times 10^7 \, \mathrm{M_\odot }$| for HD, and |$M_*(z = 3.5) \sim 2.2 \times 10^6 \, \mathrm{M_\odot }$| to |$M_*(z = 0.5) \sim 2.3 \times 10^6 \, \mathrm{M_\odot }$| for HD + Boost. Consequently, our dwarf galaxy is efficiently quenched by z ∼ 3.5. We therefore expect most of our results at z ∼ 3.5 to remain generally valid at later times, facilitating observational comparison with local dwarf galaxies.

• Magnetic fields: Although the simulations including magnetic fields look distinct in Fig. 2, by z = 4 their baryonic conversion efficiency is approximately equal to that of our’default’ SN feedback simulations, and only slightly lower for sMHD. The time evolution of the simulations with magnetic fields shows some minor differences, which are most notable at early times. In agreement with previous work, we find magnetic pressure delays the onset of star formation (occurring approximately at z ≳ 14; Martin-Alvarez et al. 2020; Koh, Abel & Jedamzik 2021), where the small-scale magnetic pressure support is also accounted for by our MTT star formation prescription (see appendix B in Martin-Alvarez et al. 2020). The stellar mass attained by z ∼ 10 depends thereby on our set-up for the magnetic field. For the intermediate and injected magnetic field simulations (iMHD and MHD), the additional magnetic support and star formation delay early on leads to a higher amount of gas accumulating in the galaxy. In the absence of any significant early stellar feedback due to, for example, stellar radiation or stellar winds, star formation proceeds unimpeded until a significant number of SN events take place. This leads to a higher stellar mass by z ∼ 10. Contrarily, for our simulation with a strong primordial magnetic field (sMHD) we find an initial reduction of the stellar mass. In this simulation, star formation is delayed even further, allowing the first SNe to commence regulating star formation. This sensitivity to the assumed initial magnetic field makes relic dwarf galaxies formed at extremely high redshifts a promising observational window to primordial magnetic fields.

• Radiative transfer: The inclusion of stellar radiation and its immediate feedback in the aftermath of star formation results in a small reduction of the stellar mass at z ∼ 10 for all our simulations with radiative transfer. In the no-magnetic field model (RT), this initial suppression does not significantly reduce the stellar mass at lower redshifts, eventually leading to a slightly higher stellar mass by z = 3.5. Such positive feedback due to stellar radiation has already been suggested previously by e.g. Smith et al. (2021). This is caused by a less efficient expulsion of gas from the dwarf galaxy when stellar radiation is included in our simulations.

  The combination of radiative transfer and magnetic fields is particularly relevant for star formation: magnetic fields affect the efficiency and spatial distribution of star formation, whereas radiation regulates its local suppression through the evaporation of molecular clouds. The sphinx-mhd simulation (Katz et al. 2021) pioneered the study of this combination of physical processes in high-resolution cosmological simulations, investigating its effects on reionization and as well as the on the properties of the simulated population of galaxies at the end of cosmic dawn (z ≳ 6). One of the results found by Katz et al. (2021) was that magnetic fields can affect the surface brightness of galaxies. Here, we delve deeper into exploring this combination of magnetic fields and radiation, although for a single dwarf galaxy.

  Our simulations with radiative transfer and magnetic fields show a stronger suppression of early star formation than the RT simulation. This is particularly notable for the RTsMHD simulation, which shows the delayed onset and suppression of star formation also observed in sMHD. This is expected, as the effects of radiation will not emerge until the formation of the first stars in our simulation. Interestingly, the suppression of the baryonic conversion efficiency in the RTiMHD and RTsMHD simulations is comparable to that of the HD + Boost simulation until z ∼ 7. However, they display a different time evolution at lower redshifts. This could be due to the more concentrated star formation in the presence of strong magnetic fields (Martin-Alvarez et al. 2020) leading to higher stellar radiation fluxes. However, this does not lead to an irreversible expulsion of gas from the system. Accretion around z ∼ 7–8 triggers an important peak of star formation at z ∼ 6. This burst of star formation drives the stellar masses of the RTiMHD and RTsMHD simulations to values comparable to the fiducial HD simulation. These simulations display an additional star formation burst at z ∼ 5. The positive feedback of the stellar radiation observed in the RT simulation disappears in the RTiMHD model, where |$\sim \!\!\, \mu \mathrm{G}$| magnetic fields are present in the galaxy from z ∼ 10. On the other hand, the dwarf galaxy in the RTsMHD simulation has an even higher final stellar mass, most likely due to the high magnetic pressure in the CGM further confining outflows and preventing outflowing gas from escaping the halo.

• Cosmic rays: In our simulation with cosmic rays and magnetic fields CRiMHD, cosmic rays only have a secondary effect at the onset of galaxy formation, behaving similarly to its non-cosmic ray counterpart (iMHD), thus reaching a higher initial stellar mass than HD due to the early effects of magnetic fields at z ≳ 10. The significant effects appear as cosmic ray energy density builds up and becomes comparable to the thermal energy density. As SN feedback events are the exclusive source of cosmic rays in our simulations, their impact is only found at later times. Indeed, from z ∼ 8 star formation is quenched in the CRiMHD simulation, and its final stellar mass ends up below the iMHD model and slightly below that of the HD simulation.

  Our most interesting models are the two ‘full-physics’ simulations, RTCRiMHD and RTnsCRiMHD. Both include magnetic fields, radiative transfer, and cosmic rays, with their only difference being whether they account for cosmic ray streaming. The addition of cosmic rays to radiation and magnetic fields further reduces the early gas content. The star formation peak observed in RTiMHD and RTsMHD occurs at z ∼ 6.3 in these simulations. RTCRiMHD and RTnsCRiMHD maintain instead a low baryonic conversion efficiency at z ≳ 5, (coincidentally) comparable to that of our HD+Boost simulation. Indeed, RTCRiMHD displays the lowest stellar mass across all our MTT-star formation simulations, even lower than HD+Boost. The HD + Boost simulation employs a comparable feedback boost to that used by the sphinx simulation. This may suggest that the inclusion of cosmic rays in this simulation may lead to the same halo mass–stellar mass relation and the same reionization history at z ≳ 5, in agreement with observations without requiring the calibration boost of SN feedback. On the other hand, Farcy et al. (2022) find the inclusion of cosmic rays to reduce the escape fraction of ionizing radiation in isolated galaxies. Furthermore, they also find cosmic rays generate a smaller stellar mass reduction than their calibrated SN feedback, with a similar boost to the one employed here. The main differences between our set-ups are the inclusion of radiative transfer in their boosted SN feedback model (as done in sphinx), and the fact that the Pandora simulations are evolved in a cosmological context, whereas the galaxies simulated by Farcy et al. (2022) are isolated. Consequently, it will be important to revisit the impact of cosmic rays on reionization using sphinx-like simulations.

  We note that RTnsCRiMHD and RTCRiMHD follow distinct evolutionary tracks in the time interval from z = 6 to z = 4. After z ∼ 5, our galaxy commences its final merger with a neighbouring dwarf galaxy. In the RTnsCRiMHD simulation, the dwarf galaxy star formation and associated SN feedback heats the neighbouring dwarf galaxy, preventing significant star formation inside it until both systems merge. On the contrary, SN feedback in RTCRiMHD occurs later, such that the star formation persists in the neighbouring dwarf galaxy for longer. This leads to the higher stellar mass in the RTCRiMHD simulation at around z ∼ 5. However, both models converge to a similar stellar mass at z ∼ 3.5.

  Finally, we note that in the presence of cosmic rays, the positive feedback effect from stellar radiation (Smith et al. 2021) in our RT simulation disappears, and in fact the stellar mass is reduced by ∼20  per cent compared to that in the HD simulation. We will later discuss the interplay of these two physical effects, as the efficiency of cosmic rays in ejecting gas is affected by radiative transfer.

The analysis of the stellar mass growth in all our simulations reveals that, due to their shallow potential, dwarf galaxies are extremely sensitive to the inclusion of different physical processes that affect their baryonic conversion efficiency. It is encouraging that a subset of ‘full-physics’ simulations converge to a similar stellar mass when the galaxy is ultimately quenched. However, we stress that all these models assume the same star formation and SN feedback model which may (in part) lead to the apparent final stellar mass ‘convergence’.

We conclude this section by comparing our simulations with the predictions by Behroozi et al. (2013), Moster et al. (2013), and Jethwa et al. (2018), as well as the observations by Read et al. (2017). First, we note that these relations (except Moster et al. 2013) correspond to z = 0 galaxies, whereas our simulations, due to their expensive computational cost, only reach z ∼ 3.5. Secondly, abundance matching relations are generally uncertain and subject to multiple caveats. In particular, the relations by Behroozi et al. (2013) and Read et al. (2017) yield higher stellar mass per halo mass. They have also been argued to be too high when considering reionization constraints (Graus et al. 2019) and observational incompleteness (Garrison-Kimmel et al. 2014). Finally, we also note that the relation obtained by Jethwa et al. (2018) is for satellite galaxies, whereas our simulation corresponds to an isolated system. With these considerations in mind, we find all our simulations with MTT star formation and standard SN feedback have |$M_*(z = 3.5) \sim 10^7 \, \mathrm{M_\odot }$|⁠, whereas those with boosted SN feedback have |$M_*(z = 3.5) \sim 2\times 10^6 \, \mathrm{M_\odot }$|⁠. These boosted SN feedback simulations have a stellar mass that may be too low when compared with both our extrapolation of the Behroozi et al. (2013) scaling relation and the estimates for LG dwarfs, but in better agreement by those predicted by Jethwa et al. (2018) from Milky Way satellites (note that our simulated galaxy is a field dwarf). None the less, we caution that there are still significant observational uncertainties and large data scatter in this low galaxy mass regime preventing us to draw any firm conclusions. Contrarily, our simulations with standard SN feedback only appear to provide a good match to both abundance matching and observations. Furthermore, due to its positive feedback, RT has a stellar mass that may be somewhat high for this halo mass range. RTsMHD appears to be in tension with both observations and predictions (Behroozi et al. 2013; Jethwa et al. 2018). We attribute this to the extremely strong primordial magnetic field used. While our results based on a single simulated galaxy should be taken with caution, we note that B₀ < 3 × 10⁻¹¹ G leads to more reasonable stellar masses. Finally, note that our RTCRiMHD and RTnsCRiMHD simulations provide the best match to comparable observed galaxies by Read et al. (2017) for all halo masses.

While reproducing the observational estimates of the M_*/(f_bM_halo) versus M_halo relation is an important test, this is not a sufficient condition to appropriately simulate dwarf galaxies as many other of their properties may be unrealistic. We proceed to review how these other properties are affected in the following sections.

3.3 Dynamical masses and sizes of stars and gas

In this section, we review how our different physical models affect the kinematics and morphology of our simulated galaxy and compare with the corresponding properties of observed galaxies. The panels of Fig. 4 show the dynamical mass (top row) and half-mass radius (bottom row) as a function of stellar (left) and gas mass (right) for a representative subset of our simulations that differ significantly from the HD simulation. We include evolutionary tracks as segmented lines, sampling the median values of each simulation during z = (7, 6, 5) ± 0.5 and z ∈ [3.4, 3.7], with symbols shown for the final redshift interval. Note that the dynamical mass is more accessible observationally than the ‘true’ halo mass. The dynamical mass is calculated separately for the stars and gas as follows:

where σ_{v, i} is the mass-weighted 3D velocity dispersion and R_{1/2, i} is the half-mass radius of either stars or gas (with i=‘*’ or ‘gas’, respectively). We limit the radial extent of all our mass measurements as well as the computation of σ_{v, i} to twice the half-mass radius of either stars or gas. The half-mass radius is in turn obtained for each component employing the mass within a 0.2 r_halo sphere centred on the galaxy.

The M_{dyn, *} versus M_* relation shows the expected correlation, with our simulations clustering in three main groups. These are in order of decreasing dynamical mass: no feedback simulations, standard SN feedback, and boosted SN feedback. Amongst the simulations without the boosted SN feedback, the inclusion of magnetic fields appears to have a relatively small effect, possibly driving a slight increase of M_{dyn, *} especially in the case with strong primordial magnetic fields. This is likely the consequence of a higher central density concentration, most pronounced for the sMHD simulation. Accounting for radiative stellar feedback generally increases the dynamical mass, and results in larger sizes and a larger baryonic mass. The most interesting effect is that of cosmic rays, which decrease the dynamical mass. Unsurprisingly, simulations with boosted SN feedback have a lower stellar mass as well as a lower dynamical mass. As HD and HD + Boost evolve to z = 0.5, they undergo a negligible increase of M_* and a small change in M_{dyn, *}, which appears to be set by the total mass of the galaxy (see Section 3.4.3).

We compare our simulations with the large sample of observations compiled by McConnachie (2012). Note that several of these observed systems correspond to satellites; they may have undergone significant dark matter-stripping while preserving their stellar mass (Penarrubia, Navarro & McConnachie 2008). On the other hand, our simulated dwarf galaxy is evolving in a ‘field’ environment without experiencing any substantial stripping. We also include data for Leo A, SagDIG, and Aquarius by Kirby et al. (2017), which may represent the closest present-day analogues of our simulated galaxy being also isolated systems. We use the mass model for WLM presented by Leung et al. (2021) to estimate a dynamical mass for that galaxy. An important caveat regarding our comparison with observations throughout this section is that our main simulated models only reach z ∼ 3.5, whereas the observations correspond to z ∼ 0. Therefore, to further reinforce our analysis, we include all simulation models that reach z = 0.5. These are labelled in Fig. 4 and the rest of the manuscript as encircled versions of the corresponding model, and show only small displacements of the dynamical masses. Overall, most of our simulations provide a good match to the observations within their considerable scatter. Bearing in mind that WLM and LeoA are the systems that best resemble the Pandora dwarf, the CRiMHD and boosted SN feedback simulations show the largest tension with these data.

The R_{1/2, *} versus M_* plot shows that the sizes of our simulated galaxy are in good agreement with the broad distribution of the observations (Wolf et al. 2010; McConnachie 2012; Kirby et al. 2013; Koposov et al. 2015)³ for most of our simulations with SN feedback. Despite a large stellar mass reduction, boosting SN feedback only has a minor impact on the stellar half-mass radius. Both for HD and HD + Boost we find negligible changes to their stellar half-mass radii when considering their evolution down to z ∼ 0.5. While inclusion of magnetic fields tends to slightly increase our R_{1/2, *} compared with their non-MHD counterparts, for stronger primordial fields than those explored here, Sanati et al. (in preparation) find a trend towards a mild size increase with increasing B₀, followed by a considerable shrinkage. However, when including radiation with our extreme primordial magnetic field (RTsMHD), we find the most considerable shrinkage. We attribute this to the combination of a lower efficiency of SN feedback due to stellar radiation (see also Rosdahl et al. 2015; Smith et al. 2021), and outflow confinement by the strong CGM magnetic pressure generated by the primordial fields. For the other models, the inclusion of radiation leads to higher R_{1/2, *} values, due to flattened and more extended stellar density radial profiles. We attribute these to a more extended distribution of the gas, as discussed later on in this section.

Furthermore, cosmic rays somewhat reduce R_{1/2, *}, particularly in the absence of radiation (CRiMHD). This is due to the non-thermal ISM support from cosmic rays (Dashyan & Dubois 2020), which inhibits the number and mass of star-forming clumps (Farcy et al. 2022). With less star-forming clumps where gravity dominates over the gas support at large distances from the galaxy centre, the galaxy features a lower R_{1/2, *} measurement. Such size reduction has also been reported for simulations of larger galaxy masses (Buck et al. 2020). The reduction of the stellar half-mass radius appears less significant for the simulations combining radiation and cosmic rays, with sizes comparable to that of RTiMHD. RTCRiMHD exhibits a decrease of R_{1/2, *} at z ∼ 3.5, which we attribute to its recent merger event⁴ (see Section 3.7).

The right column of Fig. 4 reviews the same relations for the gas component. We compare our dynamical masses with an observational estimate generated using the WLM mass profiles from Leung et al. (2021). Our simulations cluster in three groups. From higher to lower gas mass content these are no feedback, stellar radiation, and SN feedback runs. While, as expected, no feedback yields unrealistic M_gas and M_{dyn, gas} values for this halo, the inclusion of SN feedback drives our simulations to the lower values of the two quantities. In particular, the boosted feedback simulation has the lowest M_{dyn, gas} while maintaining a comparable gas content to the fiducial HD case. As we evolve HD and HD + Boost to z = 0.5, their dynamical masses measured through the gaseous component remain relatively unchanged, and their gas masses remain around |$M_\text{gas}\sim 2 \times 10^7 \, \mathrm{M_\odot }$|⁠. Simulations with magnetic fields do not significantly affect the gas mass content of our galaxies, but increase dynamical mass estimates based on the gas, consistent with higher concentrations (Martin-Alvarez et al. 2020). Including stellar radiation leads to a consistent increase of gas mass in the studied galaxy by ∼0.3–0.5 dex, in agreement with predictions from isolated galaxy simulations (Smith et al. 2021), but maintains the same approximate M_{dyn, gas}/M_gas ratio as the SN feedback simulations. Cosmic ray simulations have a mild reduction in gas mass, which we attribute to cosmic ray-driven galactic outflows. Additionally, the RTCRiMHD simulation shows an increased dynamical mass estimate from the gas dynamics at z = 3.5, likely due to a recent merger.

The gas mass–size relation shows that the boosted SN feedback produces the smallest sizes. This is likely due to the small amount of remaining gas, as the majority of gas is effectively removed from the system by the SN feedback already at high redshift. As the HD and HD + Boost simulations evolve to z ∼ 0.5, their gaseous extent decreases in size due to a reduction in the frequency of SN feedback. On the other hand, stellar radiation leads to a slightly more extended gas distribution. Due to their additional non-thermal support of the ISM, the inclusion of cosmic rays also increases R_{1/2, gas}. We include data by Valenzuela et al. (2007) and Leung et al. (2021) for comparison. All our simulations with stellar feedback have somewhat lower gas contents than observed galaxies with the same gas half-mass radius, favouring our models with radiative transfer. We note however that the stellar masses of some of the observed systems we compare to (Valenzuela et al. 2007) are higher than that of our simulated dwarf galaxy⁵, and that our system is not evolved to z = 0.

3.4 Integrated and resolved kinematics of stars and gas

3.4.1 Stellar kinematics

The study of the stellar kinematics serves as an important further diagnostic when comparing simulations with observations. Fig. 5 shows the integrated stellar kinematics of our simulated dwarf galaxy along with observational data for local dwarf galaxies as compiled by Wheeler et al. (2017).⁶ We include evolutionary tracks sampling median values during z = (7, 6, 5) ± 0.5 and z ∈ [3.4, 3.7] for the HD, RT, RTiMHD, RTsMHD, CRiMHD, RTnsCRiMHD, and RTCRiMHD simulations. Each of the kinematic quantities is computed as a mass-weighted average within 2R_{1/2, *}. From top to bottom, panels show the stellar rotational velocity (v_rot), stellar velocity dispersion (σ_{v, *}), and the ratio of rotational velocity to velocity dispersion (v_rot/σ_{v, *}). For all three quantities, our NoFb simulations are in tension with observations (not shown, due to their values being significantly above the ranges displayed for v_rot and σ_{v, *}). The inclusion of SN feedback provides better agreement. Rotational velocities with fiducial SN feedback are of the order of |$\sim \!\!1\!-\!5\, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠, and become somewhat lower when the strength of the SN feedback is boosted. During their time evolution to z = 0.5, HD and HD + Boost maintain low rotational velocities. Their evolution suggests that other models may maintain similar v_rot values or undergo a slight reduction. With additional physics included, rotational velocities of our simulated dwarf are in the |$\sim \!\!7\!-\!10\, \mathrm{km}\, \, \mathrm{s}^{-1}$| range, found in reasonable agreement with observations of isolated dwarfs (shown as cyan data points).

Most of our fiducial SN feedback simulations are located at the upper end or slightly above the σ_{v, *} observational data. The boosted SN feedback simulations, with their lower stellar masses, show a lower velocity dispersion that is in good agreement with observations. When evolved to z = 0.5, the stellar velocity dispersion of HD is reduced by |$\sim 5 \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠, which brings the model in closer agreement with observations. This suggests that simulations with additional physics may experience a similar reduction. HD + Boost maintains an approximately constant value for σ_{v, *} during its evolution. Simulations with stellar radiation have particularly large stellar velocity dispersion, whereas the models with cosmic rays (with or without radiation) have the lowest values of σ_{v, *} for a realistic M_*. RTnsCRiMHD in particular yields a reasonable match to observations of isolated dwarfs in the same stellar mass range. While RTCRiMHD is within the range of observations prior to its merger at z ∼ 5, this event is the cause of the higher final σ_{v, *} values than the RTnsCRiMHD case. The ratio of rotational velocity to velocity dispersion is within the scatter of observations for all our simulations. While some minor variation is observed when evolving HD and HD + Boostto z = 0.5, their qualitative properties remain unchanged. The only notable trend for this quantity is that the inclusion of stellar radiation leads to a more rotationally dominated support of the stellar component, which appears to best resemble the isolated dwarf galaxies from the Wheeler et al. (2017) sample.

3.4.2 H i properties

Including radiative transfer leads to a more self-consistent modelling of the H i content in our dwarf galaxy. This allows us to explore the subset of radiation-hydrodynamics simulations and compare their H i properties to observations. Fig. 6 revisits the kinematic properties of the studied galaxy, but now for the H i component, comparing with H i observations for similar mass galaxies. We compute each of these quantities adopting the same procedure as done for Figs 4 and 5, except for the maximal rotational velocity max(v_{rot, HI}), which corresponds to the maximum value of the rotational component of the velocity profiles within 2R_{1/2, HI}. We compare our models at z ∼ 3.5 with H i observational results for similar galaxies at z ∼ 0 by Chyzy et al. (2011), McConnachie (2012), Johnson et al. (2015), and Kirby et al. (2017).

In agreement with previous studies (e.g. Rey et al. 2022), we find a higher time-variability for measurements associated with the gaseous components, such as the H i values presented here. This variability is illustrated by the included error bars, where the median values remain relatively robust once the final merger event has taken place. We also find these shown quantities to remain relatively unchanged in the HD and HD + Boost simulations.⁷ Our measurements of the dynamical mass using H i (top panel) as a tracer are more robust than those using the total gas, with their scatter reduced by ∼0.5 dex and all our models matching well the comparable galaxies from Kirby et al. (2017). Our simulations spread around the expected values for max(v_{rot, HI}) (second panel), with the simulations combining radiative transfer and MHD falling slightly below the observations. For the H i velocity dispersion σ_{v, HI} (third panel), all our simulations have values of |$\sigma _{\text{v},\text{HI}} \sim 10 \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$| and the same approximate scatter, providing a good match to the distribution of observations. Finally, the ratio of max(v_{rot, HI})/σ_{v, HI} (fourth panel) is in best agreement with observations for the RT and RTnsCRiMHD simulations, dominated by the value measured for σ_{v, HI}. None the less, all our radiative transfer simulations predict reasonable H i content and properties of our simulated dwarf galaxy.

3.4.3 Comparing H i and stellar kinematics

With different components of the galaxy subject to different dynamics, we expect some variation in their kinematic properties. This is reflected in our overall results, where the stellar component displays larger velocity dispersion (v_rot/σ_{v, *} < 1), while the H i is dominated by rotation (max[v_{rot, HI}]/σ_{v, HI} > 1), in agreement with their collisional (H i) and collisionless (stellar) nature. Rotational velocities for the stars are |$\lesssim 10\, \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠, whereas the maximum H i rotational velocity tends to be within |$10\!-\!30\, \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠. In terms of velocity dispersion, H i displays |$\sigma _{\text{v},\text{HI}} \lesssim 10\, \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠, whereas σ_{v, *} tends to be within |$14\!-\!20 \, \, \, \mathrm{km}\, \, \mathrm{s}^{-1}$|⁠. Notably, employing the maximum velocity as done by the observations will bias the measurement towards higher rotational support estimates. It is also important to consider that the stellar velocity dispersion may be subject to some degree of dynamical heating due to particle mass resolution (Wheeler et al. 2019; Ludlow et al. 2023).

Another comparison worth exploring is that of the estimated dynamical masses from each component as the dynamical mass is expected to trace the total mass of the galaxy. We show in Fig. 7 the ratio of the stellar dynamical mass (M_{dyn, *}; standard symbols) and the gaseous dynamical mass (M_{dyn, HI}; lighter symbols) with respect to the total mass. Both components tend to provide relatively accurate estimates for the total mass of the galaxy. The H i measurement tends to yield ratios >1, whereas the stellar component tends towards somewhat lower values. However, combining the two estimates provides a good approximation to M_total, with the exception of the RTsMHD model with an extreme primordial magnetic field. The RTCRiMHD model shows the largest separation from the two components, likely as a consequence of its recent ‘wet’ merger.

3.4.4 Spatially resolved kinematics

Our review of various kinematic quantities revealed important differences between our models featuring different physical processes. These differences will be relevant for observations with access to spatial velocity information, employing integral field spectroscopy, such as CALIFA (García-Benito et al. 2015) or MaNGA (Penny et al. 2016). In Fig. 8 we show multiple kinematics projections of our studied dwarf system at z = 3.5, separated into its baryonic components. Each simulation is presented by a pair of columns. In the left column we show the line-of-sight stellar velocity dispersion (top row), the velocity along the line of sight for the stars (second row), the gas (third row), neutral hydrogen (fourth row), and ionized hydrogen (fifth row). These two latter components are obtained either using the ionization fractions computed self-consistently in simulations with radiative transfer, or using the gas thermodynamical properties for our simulations without. Stellar velocity projections are density-weighted, whereas gas velocity projections employ a ρ² weighting, to more closely mimic the emission-weighted maps. For the velocity projections, we display only gas cells with densities |$\rho \gt 0.1 \mathrm{m_H}\, \, \mathrm{cm}^{-3}$| to avoid a bias due to the motion of gas in the lowest density regions. We have verified that such a cut does not affect and of the main features of the projections. We smooth out gas velocity fields outside the central circle for clarity. In the right column, we show from top to bottom a synthetic SDSS mock image (following the same procedure employed in Fig. 2), stellar density, gas density, neutral hydrogen density, and ionized hydrogen density. We analyse the projections for the redshift interval z ∼ 5 to z ∼ 3.5 and describe below the main results that remain consistent through this period.

While significantly dispersion supported, as found from our integrated kinematics analysis, the HD simulation displays a level of coherent rotation in stars and gas, with the gas being more centrally concentrated and largely in the neutral state. Including radiative transfer in our simulations leads to a dwarf galaxy which is more gas-rich. This is exemplified by the RT simulation, which exhibits a more extended H i distribution and various clear dust bands. Our simulated dwarf for this particular simulation is embedded within an extended gaseous distribution. This model also exhibits a somewhat higher degree of coherence of rotation both in gas and stars, but it is worth noting that the line-of-sight velocity field and velocity dispersion are affected by a close-by gas-rich companion at this particular redshift. Gas within the dwarf contains a mixture of neutral and ionized hydrogen, and interestingly these two components have somewhat different line-of-sight velocity fields.

Turning now to our full physics run, RTCRiMHD, both the stellar and gaseous distribution is more centrally concentrated than in the RT model, with the dwarf being somewhat more gas-poor. In the synthetic SDSS image it is apparent that the morphology of the dwarf is more amorphous. Importantly, all line-of-sight velocity fields contain a number of kinematically distinct ‘clumps’, which are (partly) caused by significant cosmic ray-driven outflows, also causing kinematically misaligned and decoupled motions in stars, ionized hydrogen as well as neutral hydrogen. Performing the same analysis on the RTCRiMHD simulation at different times (not shown here) we further find that kinematic misalignment between the stellar component and ionized hydrogen are particularly prevalent during or shortly after significant star formation events that trigger cosmic ray-driven outflows.

3.5 Colour–magnitude relation

As part of our analysis, we review the optical emission properties of our simulated dwarf galaxy. Fig. 9 shows the g − i colour versus absolute magnitude G (top) and versus stellar mass within 2 R_{1/2, *} (bottom), compared with the observed properties of the ELVES dwarf galaxies (Carlsten et al. 2021, 2022). In order to measure the emitted magnitudes for each of the filters, we produce synthetic observations accounting for dust obscuration. For these, we assume each stellar particle to be a single stellar population with its emission following Bruzual & Charlot (2003) depending on its age and metallicity. We adopt a simple dust obscuration model through a 3D dust absorption screen (with metal-to-dust ratio 0.4; Kaviraj et al. 2016). Magnitudes are obtained by integrating the emission within the central 1.25 kpc. The measurement for our galaxies at z = 3.5 is shown with shaded symbols. Since for the majority of our simulations we only evolve the dwarf galaxies to z = 3.5, their stellar populations are still relatively young in comparison to the ELVES sample of local dwarfs. Hence, to provide a more adequate range for comparison, we also show the same measurement for our galaxies assuming their stellar populations have evolved passively down to z = 0.5 (i.e. we increase all particle ages by ∼6.8 Gyr) in their colour computation. The selected value of z = 0.5 redshift for this passive evolution is based on the availability of the HD and HD + Boost simulations to allow for direct comparison. Further evolution to z = 0 would increase the g − i colour by an additional |$(0.04 \pm 0.02)\, \text{mag}_\text{AB}$|⁠. Therefore, the resulting data points are obtained assuming that no additional star formation will occur during this period. Including these data points provides us with an upper limit for the colour attainable by the simulated galaxies. To further review these estimates, we also include the exact measurements at z = 0.5 for the simulations that reached that redshift (i.e. for the HD and HD + Boost models, at z = 0.5). These are shown as encircled data points. We compare this ‘manual’ ageing with the actual simulation data for the HD and HD + Boost runs which reached z = 0.5 to examine the accuracy of this approach. For HD, the approximation is remarkably accurate, with errors in the estimated colour and absolute magnitude of |$\sim 0.05\, \textrm {mag}_\text{AB}$| and |$0.05\, \textrm {mag}_\text{AB}$|⁠, respectively. These estimates experience larger deviations for the absolute magnitude in HD + Boost, with only |$\Delta (g-i) \sim 0.05\, \textrm {mag}_\text{AB}$| but |$\Delta G \sim 0.5\, \textrm {mag}_\text{AB}$|⁠.

We conclude that the ‘passive ageing’ procedure leads to reasonable estimates of the observed g − i colour, with a possible overestimation of the absolute magnitude G of up to |$\sim 1 \, \textrm {mag}_\text{AB}$|⁠, at least for those two models. If our other models followed reasonably similar evolutions, we would expect the majority of them to be relatively close to the estimated colour upper limit. None the less, we note that some degree of additional star formation may take place in Pandora models, and some could even undergo a mild re-ignition of star formation (Rey et al. 2020).

Unsurprisingly, at z ∼ 3.5, our galaxies are relatively blue compared with the overall observed ELVES population. By z ∼ 0.5 however, we expect most of our simulated dwarfs to be in good agreement with the ‘Early type’ observed sample, and possibly at the upper colour end of the ‘Late type’ if they underwent some star formation. By ‘manually’ ageing the simulated galaxies we have assumed that they remain completely quenched from z = 3.5 to z = 0.5, hence our g − i colours should be considered only as a conservative upper limit, and the stellar masses and absolute G magnitude as a conservative lower limit. If the trend displayed by the HD and HD + Boost models is reproduced by our other simulations, this could potentially imply that all of our feedback models would lead to a population of early-type dwarfs in the local Universe. Hence alternative feedback models, assembly histories, or star formation re-ignition are needed to produce the observed sample of ‘Late type’ dwarfs. To determine this, simulation to low redshifts and larger samples of simulated dwarfs are needed.

3.6 Constraints on magnetic field strength

Galaxy formation simulations simultaneously featuring radiative transfer and magnetic fields are still rare. The leftmost panels in Fig. 10 show how our simulations compare with observational relations between various H i properties and magnetic fields in local dwarf galaxies by Chyzy et al. (2011). Our average magnetic field at galactic scales is of the order of |$\sim \!2\!-\!5 \, \mu \mathrm{G}$|⁠, with slightly higher field strengths in the RTsMHD simulation due to a rather high value of its primordial magnetic seed. Encouragingly, for all of the H i quantities examined, our galaxies are either consistent with current observations or fall within the range where Chyzy et al. (2011) can only provide upper limits for the magnetic field. Simulated galaxies have magnetic fields which are also compatible with measurements at higher M_HI values, where these observations suggest |$B \sim 3 \!-\! 15 \, \mu \mathrm{G}$|⁠.

Detailed simulations of galaxies modelling the ionization state of hydrogen as well as MHD, as is the case for part of our Pandora suite studied here, will help our understanding of observations in the radio (Heald et al. 2021; Lopez-Rodriguez et al. 2022b) and far-infrared (Borlaff et al. 2021) wavelength ranges. These simulations will be crucial to aid both current and upcoming surveys such as the Survey of extragALactic magnetiSm with SOFIA (SALSA) using the Stratospheric Observatory for Infrared Astronomy (SOFIA; Lopez-Rodriguez et al. 2022a) or SKA. These surveys aim to probe magnetization down to dwarf galaxies. We showcase this capability for our simulations on the right-hand side of Fig. 10, where we present synthetic synchrotron observations for the simulations combining radiative transfer and magnetic fields.

To produce these synthetic observations we employ the publicly available code polaris⁸ (Reissl, Wolf & Brauer 2016; Reissl et al. 2019). We extract from the simulations an adaptive grid centred on the galaxy, resolved with double the local AMR resolution, interpolating all native ramses quantities required by polaris. Due to their limited spatial resolution, our simulations do not fully reproduce the entire magnetic energy inverse-cascade (Martin-Alvarez et al. 2018), thus failing to capture additional magnetic energy expected to reside below the local grid size (Schekochihin et al. 2002). As we are only interested in a qualitative comparison across our models, we account for the loss of additional magnetic energy by a simple boost of the magnetic field strength. To estimate the approximate amount of energy boost required, we expand a Kazantsev-like spectrum from our simulation resolution to the scales where turbulence should be converged. Martin-Alvarez et al. (2022) show the peak of the magnetic energy spectrum k_peak to scale with the finest simulation resolution Δx_min following |$k_\text{peak} \propto \Delta x_\text{min}^{-1}$|⁠. Körtgen, Federrath & Banerjee (2017) suggest resolutions |$\Delta x_\text{min} \lt 0.1 \, \mathrm{pc}$| are required to approach turbulence convergence, finding the average turbulence to stabilize at values of |$\Delta x_\text{min} \lesssim 0.03 \, \mathrm{pc}$|⁠. Extrapolating the Martin-Alvarez et al. (2022) scaling of k_peak from our suite resolution |$\Delta x_\text{min,sim} \sim 7 \, \mathrm{pc}$| to |$\Delta x_\text{min,exp} = 0.005 \, \mathrm{pc}$| yields a ratio of k_{peak, exp}/k_{peak, sim} ∼ 40.1. Expanding a Kazantsev-like spectrum (∝ k^3/2) by such ratio provides an increase of the average magnetic field strength of approximately 2 dex. As we are only interested in a qualitative exploration of these synthetic observations, we leave the review of alternative models to capture the sub-grid magnetic energy (e.g. Reissl et al. 2019) to future work.

The use of polaris requires additional quantities not modelled by our version of ramses. We follow a similar approach to that used by Reissl et al. (2019) to determine these quantities on the grid. We bound the Lorentz factor of electronic cosmic rays between a fixed γ_min = 4 minimum (Webber 1998) and a γ_max = 300 maximum (i.e. γ_max ≫ 1). We also follow Reissl et al. (2019) and fix the power-law index to p = 3 (Miville-Deschênes et al. 2008), related to the spectral index by α = (p − 1)/2. In order to obtain the distribution of thermal electrons, we compute their number density as a function of gas number density according to Pellegrini et al. (2020). Finally, we model the local energy density of electronic cosmic rays |$e_{e^{-} \text{CR}}$| following the CR2 model also presented by Reissl et al. (2019), which assumes energy equipartition between electronic cosmic rays and the local magnetic field (without any magnetic field increase). The produced synthetic radio observations show the resulting intensity map at λ = 6.2 cm. We overlay on these maps white dashes oriented perpendicularly (i.e. rotated 90 deg) to the local linear polarization (accounting for Faraday depolarization) in order to align with local magnetic field. These dashes are only shown in regions where the linearly polarized intensity is higher than 5  per cent of the average intensity of the entire map. Finally, we assume the observed galaxy to be at a distance of 2 Mpc, convolving our images with a 2D gaussian profile with a full width at half-maximum equal to the resolution of the corresponding telescope (12.6 arcsec for VLA and 0.1 arcsec for SKA).

In all maps, synchrotron emission traces regions of strong magnetic field. Due to its extreme primordial magnetic field, RTsMHD features the brightest emission, extending well beyond the galaxy. The observed polarization traces the large-scale structure of the magnetic field at the outskirts of the galaxy and into its halo. Contrarily, models with SN-injected magnetic fields have weaker and more turbulent magnetic fields outside the galaxy. In these simulations, the polarization vectors seemingly match only features at and below galactic scales. Our simulations with cosmic rays have a less extended emission, which is seen as a single gaussian in the VLA observations. While magnetic field lines are comparably turbulent in the three SN-injected scenarios, the polarization vectors appears more coherent in the models combining radiative transfer and cosmic rays, particularly at scales of the VLA resolution.

3.7 Halo density profiles

Dark matter density profiles in dwarf galaxies estimated from observations frequently suggest shallower central slopes than that of the theoretical expectation from dark matter assembly in a cold dark matter cosmology (Navarro, Frenk & White 1997). Such cores may be carved out of the high-density NFW cusps through the action of baryonic feedback sufficiently violent to induce rapid variations in the gravitational potential (see e.g. Navarro et al. 1996; Pontzen & Governato 2012). Cores may be particularly prevalent in dwarf galaxies due to their shallower gravitational potentials, which make their dark matter more susceptible to the SN-driven outflows (Read, Agertz & Collins 2016). Hence, in this section, we analyse the density distribution of our simulated dwarfs using their density profiles, and compare these with observations.

We centre our profiles in the galaxy rather than using the centre of the dark matter halo, motivated by our aim to resemble observational measurements. We review the effects of centring the measurements on the dark matter component in Appendix  B. We find an overall density increase of |$\sim 0.3 \cdot 10^8 \, \mathrm{M_\odot }$| for this choice, with the global trends for each set of models preserved. Density profiles may be at times sensitive to the selected galaxy centre, which is not always unambiguously defined. To obtain the galaxy centre for the radial profiles, we commence our calculation using the centre of the host dark matter halo, and select an initial sphere of the baryonic mass extending up to 0.2 r_halo. We recursively apply the shrinking spheres method (Power et al. 2003) to the baryonic mass inside this sphere until a centre is obtained.

Fig. 11 compares the average central dark matter density within the inner 150 pc of our dwarfs with data from Read et al. (2019). Following Read et al. (2016) and Read et al. (2019), we compute the bands of expected central densities for cusps and cores according to the cosmology of our simulations and employing the mass–concentration relation as a function of redshift provided by Macciò et al. (2007). These bands indicate their density expectation for haloes that have (blue bands) or have not (black bands) undergone a transition from cusps to cores, at redshifts z = 3 (shaded bands) and z = 0 (dashed lines). The vast majority of our simulations have central densities in agreement with cusp profiles (or even higher). At face value, the lack of cores could be attributed to the fact that our simulations have been evolved only to z ∼ 3.5, and thus cannot capture additional SN feedback cycles. While there is the possibility of additional star formation, based on the few simulations evolved to lower redshift, we do not expect significant evolution for the Pandora dwarf after z ∼ 4. When evolved to z = 0.5, both the HD and HD + Boost simulations remain within the z = 0 cuspy halo band (dashed black lines), which is approximately equal to the band at z = 0.5 (not shown for clarity). HD has no change in central density, whereas it is somewhat reduced in HD + Boost. Considering that Pandora does not undergo any additional mergers beyond z ∼ 5 (see Appendix  A), we expect all our other models with SN feedback to either maintain their central densities or further reduce them as they evolve to z = 0. Read et al. (2019) suggest that additional star formation at latter redshifts would lead to lower central densities. This is in qualitative agreement with our models displaying significant star formation after z ≲ 6. On the other hand, the models where most of the star formation takes place at earlier times (e.g. the iMHD, NoFb models) have considerably higher central densities. Although cusp destruction is hypothesized to be driven by SN feedback, we note that boosting the SN specific energy has only a moderate effect in our simulations by z = 3.5, even though this run is very efficient at quenching the dwarf and was driving a significant gas fraction out of the halo. The inclusion of radiative transfer leads to lower central dark matter densities. This effect is intensified when stellar radiation is combined with cosmic rays, as reflected by the tracks of RTCRiMHD and RTnsCRiMHD runs at z ∼ 5. While these lower densities persist in RTnsCRiMHD down to our final studied redshift, the central dark matter density in the RTCRiMHD simulation increases after z ∼ 4. We attribute this to a ‘wet’ merger, which erases the cored profile in the RTCRiMHD simulation and leads to a central density that resembles more the fiducial HD case. None the less, we note that the final central density measured for RTnsCRiMHD increases by a non-negligible amount when measured using the dark matter centre (reaching |$\sim 1.5 \times 10^8 \, \mathrm{M_\odot }$|⁠).

These trends are further quantified in Fig. 12, which shows density profiles for the gaseous, stellar, and dark matter components. The leftmost column shows some representative simulations (NoFb, HD, HD + Boostand RTiMHD) at z = 3.5. The vertical solid line indicates r = 0.15 kpc, used in Fig. 11 to calculate the average central dark matter density. The dotted lines spanning from |$r = 0.15 \, \mathrm{kpc}$| up to |$r = 2 \, \mathrm{kpc}$| show the NFW density scalings ρ_DM ∝ r⁻¹ (black) and ρ_DM ∝ r⁻³ (grey). The grey and blue bands within |$r = 0.15 \, \mathrm{kpc}$| denote the values of ‘cuspy’ and ‘cored’ dwarf profiles, respectively. These are computed as in Fig. 11, and are shown at z = 3 (shaded bands) and z = 0 (dashed horizontal lines). The y range shown corresponding to the halo mass that matches our simulated dwarfs at z = 3.5.

Compared to the no feedback simulations, the inclusion of SN feedback and radiative transfer (see bottom-left panel) leads to a significantly shallower central dark matter density profile that persists to up to a few 100 pc. However, this shallower dark matter density profile still falls within the ‘cusp’ category as defined by Read et al. (2019), which is indicated by the grey coloured bands. In terms of the baryonic contribution to the total density profile, it is interesting that due to SN feedback and local radiation feedback from young stars, the central stellar distribution is diffuse with gas dominating in the centre.

The central and right-hand columns of Fig. 12 show the redshift evolution of the density profiles from z = 5 to z = 3.5 for the ‘full-physics’ RTCRiMHD and RTnsCRiMHD simulations, respectively. Here, focusing first on the dark matter density profiles, the central profile is somewhat shallower than the ρ_DM ∝ r⁻¹ scaling and there is a range of redshifts where the central density distribution is more in line with a ‘core’. As discussed previously, due to the merging event taking place at z ∼ 4 both the baryonic and dark matter density profiles steepen in the central region for the RTCRiMHD simulation, while this is not the case for RTnsCRiMHD. This indicates that, while radiative transfer and cosmic rays contribute to the formation of dark matter cores, ‘wet’ galaxy mergers may transform these cores to cusps again. Furthermore, it is interesting to contrast the distribution of the baryonic component when the cosmic rays are included as well. While in the SN feedback only simulation as well as the RTiMHD model, gas is the dominant component in the centre of the galaxy, this is not the case in the simulations with cosmic rays. Now, at all redshifts explored, the stellar component dominates the total density profile within the central region with gas being significantly depleted (apart from the aftermath of the merger at z = 3.5 in the RTCRiMHD simulation). This is caused by efficient cosmic ray-driven outflows, pushing the gas out of our simulated dwarfs.

4 CONCLUSIONS

In this paper, we introduce the Pandora suite of high-resolution cosmological zoom-in simulations of a dwarf galaxy (with halo mass |$M_\text{vir} (z = 0) \approx 10^{10} \, \mathrm{M_\odot }$|⁠) combining magneto-hydrodynamics, radiative transfer, and cosmic rays. This set of simulations is generated using our own modified version of the ramses code (Teyssier 2002), employing multiple extensions by Fromang et al. (2006, MHD), Rosdahl & Teyssier (2015, RT), and Dubois et al. (2019, CRs). Our suite of simulations builds up from a fiducial model (HD) comprised of hydrodynamics, an MTT star formation model and a mechanical SN feedback scheme. For a subset of our simulations we also explore an alternative star formation model based on a gas density threshold and vary the SN feedback strength. We investigate multiple configurations for these models, gradually increasing the complexity of physical mechanisms, ultimately leading to two ‘full-physics’ simulations, RTnsCRiMHD and RTCRiMHD, which simultaneously account for magnetic fields, stellar radiation, and SN-generated cosmic rays.

The simulated dwarf galaxy is a gas-rich system at very high redshift, but does not evolve significantly after z ∼ 4, as it resides in a field environment. Consequently, we evolve all our simulations, summarized in Table 1, down to z ∼ 3.5. We also evolve a small subset of simulations all the way to z = 0.5 confirming the largely passive evolution since z ∼ 3.5, when the dwarf is quenched. We compare in detail the properties of our simulated dwarf with a wealth of observational data, with the aim to place constraints on the most likely physical mechanisms regulating the evolution of field dwarfs. Our main findings are as follows:

• Unsurprisingly, due to their shallow potential well, the properties of dwarfs are very sensitive to the physical processes included in the simulation. SN feedback remains one of the key processes regulating their final stellar masses. With our fiducial SN feedback, simulated dwarf galaxies are approximately located on the extrapolated stellar mass–halo mass relation of Behroozi et al. (2013), and agree with observations of isolated dwarf galaxies of similar halo masses (Read et al. 2017). While the combination of radiation, cosmic rays, and magnetism has only a moderate effect on the final stellar mass of our galaxy, both radiation and cosmic rays can significantly delay the growth of stellar mass in dwarf galaxies.

• In terms of the overall morphology and spatial gas distribution, models with strong SN feedback produce an overquenched, amorphous, and compact system. The inclusion of stellar radiation leads to a much more extended and gas-rich dwarf galaxy, with some central dust lanes observed in our mock optical images. Models that incorporate cosmic rays in addition to radiation lie somewhere in between of these two extremes, due to cosmic ray-driven outflows that deplete the gas reservoir in the ISM and can more effectively quench star formation.

• Consequently, strong SN feedback simulations are a poor match to observed dwarfs with a similar estimated halo mass (LeoA, WLM, SagDIG) in the dynamical mass–stellar mass relation, while our ‘full physics’ simulations provide much more reasonable predictions for the mass–size and dynamical mass–mass relations. All of this implies that ‘dialling’ up SN feedback strength in simulations to match (some) observational constraints cannot realistically account for the effects of radiation and cosmic rays.

• While a number of models with stellar radiation lead to integrated rotational velocities and velocity dispersion for stars and H i in good agreement with the kinematics of local isolated dwarf galaxies, spatially resolved kinematics reveals that cosmic ray-driven outflows can induce more realistic and diverse kinematics, with distinct ‘clumps’ and misaligned motions in stars, ionized, and neutral hydrogen as observed with IFU surveys in some dwarfs.

• Our fiducial SN models display cusp-like dark matter profiles, where increasing the SN feedback strength only leads to a minor reduction of profile cuspiness. This is the case even for the overquenched models, as the majority of the SN feedback takes place in a single burst. However, episodic removal of gas in our ‘full physics’ models leads to more core-like dark matter profiles. We note that this cusp-core transformation and its longevity is further compounded by the dwarf galaxy merger history, with a single ‘wet merger’ able to re-establish a cusp.

While our simulations explore additional physical processes frequently omitted in many numerical studies (due to their complexity and/or their computational cost), we note that our models are still far from complete. Some of the physical processes that we are missing and which may affect dwarf galaxy formation are stellar winds (Agertz et al. 2021), higher accuracy modelling of ISM turbulence either by alternative refinement strategies (Martin-Alvarez et al. 2022) or subgrid turbulence models (Semenov, Kravtsov & Gnedin 2016; Kretschmer & Teyssier 2020), and more realistic cooling prescriptions (Katz 2022). Furthermore, the possibility of AGN activity in dwarf galaxies has gained traction in recent years (e.g. Pardo et al. 2016), with numerical studies suggesting that AGN feedback has the potential to affect the formation and evolution of dwarf galaxies (Dashyan et al. 2018; Koudmani et al. 2019; Koudmani et al. 2022).

Our numerical simulations have unveiled the complexity of physical processes needed to more realistically model dwarf galaxies. ‘SN feedback-only’ models struggle to match realistic masses, sizes, and kinematics of observed dwarf galaxies, either leading to overquenched objects (for their halo mass) or too centrally concentrated baryons.

Inclusion of local stellar radiation sources and SN-driven cosmic rays leads to more extended, rotationally supported systems, where star formation and feedback is more spatially distributed and better able to regulate dwarf properties. Detailed resolved kinematics of dwarf galaxies from IFU surveys together with the upcoming JWST constraints on the multiphase nature of outflows in dwarf galaxies working in conjunction with detailed simulation models such as attempted here will help us to unravel the surprisingly complex nature of our Universe’s smallest building blocks.

ACKNOWLEDGEMENTS

The authors kindly thank the referee for their careful consideration of this manuscript, and their insightful comments and suggestions, which have highly improved the quality of this manuscript. The authors would like to thank M. Smith for sharing his initial conditions. This work was supported by the ERC Starting Grant 638707 ‘Black holes and their host galaxies: co-evolution across cosmic time’. SMA, DS, and MGH acknowledge support from the UKRI Science and Technology Facilities Council (grant number ST/N000927/1). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

DATA AVAILABILITY

The data employed in this manuscript are to be shared upon reasonable request contacting the corresponding author.

Footnotes

References

APPENDIX A: HALO MASS GROWTH AND THE RE-IGNITION OF STAR FORMATION

There is an integral connection between the evolution of galaxies and their hosting dark matter haloes. Exploring the evolution of the halo for the studied galaxy thus provides further insight on the fate of the galaxy. Events such as mergers, sustained gas accretion, or encounters with gaseous structures have the potential to re-ignite star formation (Wright et al. 2019; Rey et al. 2020; Gutcke et al. 2022). In order to further illustrate the lack of significant evolution for the Pandora galaxy after the simulated period (i.e. z = 3.5), we show in Fig. A1 the halo mass growth. This is shown for the HD simulation until z = 0.5, and down to z = 0 for a dark-matter-only simulation. We include mergers with systems down to 0.1 of the mass of the galaxy, as red (minor mergers) and blue (major mergers) circle markers, respectively. The studied halo and galaxy only undergo significant mergers prior to z = 4, and evolve secularly afterwards. We also verify the absence of any significant gas accretion or encounters with any gaseous structures in HD down to z = 0.5 (see Fig. A2). While there is some progressive halo growth during the z = 3.5 to z = 0 interval, the absence of any further galaxy mergers leads to no star formation re-ignition through a merger-driven channel. Combined with the lack of significant gas accretion, we predict a very low probability for star formation re-ignition in the presence of SN feedback. This is in agreement with the expected evolution of a very isolated dwarf galaxy.

APPENDIX B: CENTRAL DENSITIES OF THE GALAXY VERSUS THE HALO

In Section 3.7 we measured the central density of the Pandora galaxy to explore how different physics affect the formation of dark matter cores. In order to better reproduce observations, we performed this measurements centred on the galaxy. While galaxies such as Pandora are generally located at the centre of their hosting dark matter halo, we measure small separations of ∼20–30 pc, with somewhat larger deviations during disruptive events such as mergers. Studies that focus on dark matter core formation also opt to centre their measurements on the dark matter component (e.g. Orkney et al. 2021). To explore the impact of our centring choice, we show the central densities of Pandora in Fig. B1, now using the halo centre obtained with our shrinking spheres algorithm applied exclusively to the dark matter. This selection results in an approximate increase of the average central density by |$\sim 0.3 \times 10^8 \, \mathrm{M_\odot }\, \mathrm{kpc}^{-3}$|⁠, with a somewhat larger effect on RTnsCRiMHD and HD. Importantly, the described trends for the different physical models are preserved when employing this alternative dark-matter-only centring.