Sergio Martin-Alvarez

Debora Sijacki

Martin G Haehnelt

Marion Farcy

Yohan Dubois

Vasily Belokurov

Joakim Rosdahl

Enrique Lopez-Rodriguez

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

Keywords: magnetic fields -radiative transfer -cosmic rays -galaxies, dwarf -galaxies, formation -methods, numerical

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 chal-★ E-mail: martin-alvarez@stanford.edu (SMA)

lenging the standard cosmological ΛCDM 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 et al. 2018), photoevaporation (Bullock et al. 2000) linked to a filtering mass (Gnedin 2000;Okamoto et al. 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 et al. 2009).

Detailed observations are mostly limited to those found within the Local Group or not far beyond (e.g. Read et al. 2017;Wheeler et al. 2017;Kirby 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 et al. 1996). This process is particularly efficient when driven by resonantly cyclic SN bursts (Pontzen & Governato 2012;Governato et al. 2012), where systems featuring younger stellar populations could produce more cored profiles than those in which star formation ceased rapidly early on (Read et al. 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 et al. 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;Pillepich et al. 2018b;Hopkins et al. 2018). These scaling relations are rather uncertain in the dwarf regime, where very high resolution simulations are needed to resolve the internal processes and the ISM of dwarfs (see e.g. Wheeler et al. 2019;Smith 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 et al. 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 (AGN) in dwarf galaxies has sparked interest in whether AGN feedback in these galaxies may have a significant effect on their evolution (e.g. Habouzit et al. 2017;Dashyan et al. 2018;Koudmani et al. 2019Koudmani et al. , 2021Koudmani et al. , 2022)). However, wellknown 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 impor-tance 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 ≳ 𝜇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 et al. 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 nonetheless not expected to have a major effect on fundamental properties such as the final stellar mass of galaxies (Su et al. 2017;Pakmor 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 recognised 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 selfregulation 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 et al. 2018;Dashyan & Dubois 2020;Buck et al. 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 et al. 2019;Dashyan & Dubois 2020;Nuñez-Castiñeyra et al. 2022). When converting a small fraction of SN energy into cosmic rays (typically ∼10%; e.g. Hopkins et al. 2020;Dashyan & Dubois 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 nonlinear 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 et al. 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 [START_REF] Breitschwerdt | Astronomy and astrophysics[END_REF]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 star formation and mechanical SN feedback with magnetic fields, onthe-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.

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 evolu-tion 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; Teyssier et al. 2006;Fromang et al. 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 setup 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, initialised at 𝑧 = 127 and discretised with a uniform grid of 256 3 cells. At the centre of our computational domain a dwarf galaxy forms in a halo of virial mass 𝑀 vir (𝑧 = 0) ∼ 10 10 M ⊙ . We follow the formation of this dwarf galaxy using a convex hull zoom region. At 𝑧 = 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 𝑚 DM = 1.5 • 10 3 M ⊙ and 𝑚 * = 400 M ⊙ , respectively. In the zoom region, we allow the AMR to refine the grid down to a resolution of Δ𝑥 ∼ 7 physical pc (or 3.5 pc radius/half-cell size). Cells are marked for refinement when their total contained mass (Ω 𝑚 /Ω 𝑏 ) 𝑚 baryons + 𝑚 DM surpasses 8 × 𝑚 DM or when their size is Δ𝑥 > 𝜆 𝐽 /4, with 𝜆 𝐽 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 Z ⊙ , corresponding to the critical metallicity required for gas fragmentation to allow PopII stellar clusters to be formed ( Z ⊙ = 0.012; Schneider et al. 2012). Due to their computational cost, we evolve the majority of our simulations only to 𝑧 = 3.5, when the mass of the studied halo is 𝑀 vir (𝑧 = 3.5) ∼ 5 • 10 9 M ⊙ . We only continue some of our hydrodynamical runs (HD, HD+Boost, NoFb and NoFb+NoZ; listed in Table 1) down to 𝑧 = 0.5. Unless indicated in Section 2.5, all our simulations include metal cooling above and below a threshold of 10 4 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 𝑧 < 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 (𝑟 < 0.2 𝑟 halo ). (Right-hand panels) Synthetic observations of the galaxy in a (2.5 kpc) 3 box. The top panel shows radio synchrotron (𝜆 = 6.2 cm, dashes show 90º-rotated polarisation to align with the magnetic field) as would be observed by VLA in its D configuration in the top right corner and for SKA-like resolutions at the bottom left. The middle and bottom panels, show respectively the optical emission as would be observed by SDSS, and the near infrared as would be observed by JWST. For these synthetic observations we artificially position the galaxy at a distance of 2 Mpc from Earth in order to emulate Local Group dwarf galaxies analogue distances. Each image is convolved with the PSF of the corresponding telescope. For the optical and near-IR panels, we provide the apparent magnitude within the circular aperture displayed by the white circle wedges. (Bottom row) From left to right, these four panels are projections of a (7.5 kpc) 3 box showing magnetic field, hydrogen ionization fraction, cosmic ray energy density, and gas density (gray) separated into inflowing (blue gas; 𝑣 gas,radial < -10 km s -1 ) and outflowing (red gas; 𝑣 gas,radial > 10 km s -1 ) gas. The galaxy in our full-physics model appears particularly extended (compared to the fiducial hydrodynamical case, discussed in Fig. 2), with an envelope of outflowing gas that correlates spatially with the high energy density region of cosmic rays and strong magnetic fields, extending to approximately 3 kpc.

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: a magneto-thermoturbulent (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),

𝜌 star = 𝜖 ff 𝜌 𝑔 𝑡 ff , (1)
with star formation efficiency, 𝜖 ff , and free-fall time, 𝑡 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 multi-free 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 𝜌 𝑔 > 𝜌 th = 10 m H / 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 a SN event injects into its hosting and neighbouring cells mass, momentum and energy, with a specific energy of 𝜀 SN = 𝐸 SN /𝑀 SN (except for the boosted feedback simulations -as indicated in Section 2.5), where 𝐸 SN = 10 51 erg and 𝑀 SN = 10 M ⊙ . We assume a Kroupa initial mass function (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.

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 ì 𝐵 as a face-centred quantity in each cell. The algorithm ensures that the divergence of the magnetic field fulfils the solenoidal constraint ( ì ∇ • ì 𝐵 = 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 methods1 (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 discretised 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 ì 𝐵 ≠ 0 field. We investigate two different approaches: a) an ab-initio magnetic field, and b) a SN-injected magnetic field. In the first approach, we permeate the simulated domain with a uniform magnetic field along its 𝑧 axis with comoving strength 𝐵 0 . This method is the most commonly employed in MHD simulations (e.g. Pakmor et al. 2016;Martin-Alvarez et al. 2018;Marinacci 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 𝐵 0 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 (𝐵 0 > 10 -12 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 ì ∇ • ì 𝐵 = 0 when magnetising 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 𝐸 inj,mag = 0.01𝐸 SN ∼ 10 49 erg. This corresponds to a magnetic field strength of ≳ 10 -5 G when injected at scales of ∼ 10 pc, in reasonable agreement with the observed high magnetisation of supernova 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).

Radiative Transfer

The implementation for ionizing emissivity in our simulations is similar to the one in the sphinx simulations (Rosdahl et al. 2018(Rosdahl et al. , 2022)), which provide a good match to observational constraints on the reionization history of our Universe. We employ the ramsesrt implementation by Rosdahl et al. (2013); 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 (Δ𝑥 ∼ 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: HI (13.6 eV ≤ 𝜖 photon < 24.59 eV), HeI (24.59 eV ≤ 𝜖 photon < 54.42 eV), and HeII (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 𝑐. 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 et al. 2008;Stanway et al. 2016) according to particle mass, metallicity and age.

Cosmic rays

Our simulations including cosmic rays use the cosmic ray implementation in ramses by Dubois & Commerçon (2016); 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 𝜅 ∥ = 3 • 10 28 cm2 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 et al. 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 𝐸 CR = 𝑓 CR 𝐸 SN = 10 50 erg and is extracted from the standard thermal injection. The selected fraction 𝑓 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 𝑓 CR have been frequently studied in the literature (e.g. Pfrommer et al. 2017a;Wiener et al. 2017;Butsky & Quinn 2018;Dashyan & Dubois 2020).

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 (𝐵 0 = 3•10 -13 G) simply labelled MHD; a strong primordial magnetic field (𝐵 0 = 3•10 -11 G) labelled sMHD, and magnetised SN seeding (note that these simulations still include a negligible primordial field, 𝐵 0 = 3 • 10 -20 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 𝜀 SN = 2 𝐸 SN / 0.5 𝑀 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 2 reproduces well the observations by Read et al. (2017) and the stellar mass to halo mass ratio predicted by abundance matching (Behroozi et al. 2013) in their low-mass galaxies at 𝑧 = 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:

(i) HD simulations, exploring different prescriptions for star formation and parameters for stellar feedback.

(ii) MHD simulations, exploring different sources of magnetic fields.

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

(iv) 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 7 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 fullphysics dwarf galaxy at 𝑧 = 3.5 is shown in Fig. 1. Gas density maps and synthetic optical views for other representative simulations are presented in Fig. 2.

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 features 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 (𝑀 vir (𝑧 = 0) ∼ 10 10 M ⊙ ) and within a comparable stellar mass range (𝑀 * ∼ 10 6 -5•10 7 M ⊙ ). 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 𝑧 = 0 are candidates for reionization quenching of star formation. Some of these Table 1. Suite of simulations studied in this work. Columns indicate from left to right the unique symbol and ID label of the simulation, the solver employed, the initial (uniform) magnetic field seed 𝐵 0 , whether the simulation accounts for radiative transfer (Section 2.3) and cosmic rays (Section 2.4), the star formation prescription, the stellar feedback configuration, and further details regarding the configuration of the simulation. From top to bottom, simulations increase in complexity, accounting for additional baryonic physics. systems, such as LeoT (𝑀 halo (𝑧 = 0) ∼ 5 • 10 8 M ⊙ ) 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.

Simulation

RESULTS

Dwarf galaxy morphology with different ISM physics

We evolve all our simulations from 𝑧 = 127 to 𝑧 = 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 kpc) 3 cubic boxes centred on the galaxy with their line-ofsight 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 [START_REF] Bruzual | [END_REF]. 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 𝑅 𝑣 ∼ 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, 3 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 over-dense '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, 2 rightmost columns and central rows, 2 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 below. 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 over-massive and over-dense 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.

The stellar mass to halo mass relation

One of the most studied relations in galaxy formation is the stellar mass 𝑀 * to halo mass 𝑀 halo relation. A physically insightful form of this relation employs the baryonic conversion ratio, 𝑀 * /( 𝑓 𝑏 𝑀 halo ), versus 𝑀 halo , where 𝑓 𝑏 = Ω b /Ω m is the baryon fraction, and Ω b , Ω 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 contextualise 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 𝑀 halo ∼ 10 10 M ⊙ 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 et al. (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 (𝑀 halo (𝑧 = 6) ∼ 1.5•10 9 M ⊙ ) 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 𝑧 = 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 onto 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 (𝑧 ∼ 10), SN feedback leads to a large reduction in stellar mass at lower redshifts. In fact, at 𝑧 ∼ 4, the HD simulation reaches a conversion ratio of ∼2%, 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 𝜌 gas ≳ 10 -21 g cm -3 , whereas in the density threshold case, we allow star formation only when 𝜌 gas > 𝜌 th ∼ 10 -23 g 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 ∼ 0.4% by 𝑧 = 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 𝑧 ∼ 0. However, we extend some of our hydrodynamical simulations to 𝑧 = 0.5 to demonstrate that our dwarf galaxy and its halo do not undergo significant mass evolution after 𝑧 ∼ 4. In particular, the halo mass evolves from 𝑀 halo (𝑧 = 3.5) ∼ 4 • 10 9 M ⊙ to 𝑀 halo (𝑧 = 0.5) ∼ 8 • 10 9 M ⊙ , but most importantly, the stellar mass is barely changed, evolving from

�� -� �� -� �� -� �� � �� � � * /(� � � ���� ) � ���� ( � ☉ ) ��� � ������ ��� ��� ������� ��� ���� ��� ��� ��� J e
�� � �� -� �� -� �� -� �� � � ���� ( � ☉ ) � * /(� � � ���� ) ��� � ������ ��� ��� ������� ��� ���� ��� ��� ��� ��� J e th w a + 2 0 1 8 M os te r+ 20 13 (z = 3. 5) Be hro ozi +20 13 (z = 0)
Re ad +2 01 7

z = 10 z = 8 z = 7 z = 6 z = 5 z = 3.5 M ha
�� � �� -� �� -� �� -� �� � �� � ���� ( � ☉ ) � * /(� � � ���� ) ��� � ������ ��� ��� ������� ��� ���� ��� ��� ��� J e th w a + 2 0 1 8 M os te r+ 20 13 (z = 3. 5) Be hro ozi +20 13 (z = 0)
Re ad +2 01 7 For each simulation, we show the conversion efficiency at several redshifts. For the NoFb, HD, and HD+Boost simulations, we include four additional evolution markers corresponding to redshifts 𝑧 = 0.5, 1, 2, 3 moving from rightmost towards the left of the panel. These illustrate how simulations including SN feedback have negligible star formation between 𝑧 = 3.5 and 𝑧 = 0.5. The blue band shows an extrapolation of the baryonic conversion efficiency predicted by Behroozi et al. (2013) for 𝑧 = 0, the yellow band indicates the 1𝜎 prediction by Moster et al. (2013) at 𝑧 = 3.5 and the green band corresponds to the abundance matching prediction from Read et al. (2017). We also include the stellar mass-halo mass relation predicted by Jethwa et al. (2018) for the Milky Way satellites population as the red band in the bottom right corner. We provide additional details regarding their comparison and about additional relations at the end of Section 3.2. The dotted black line indicates an evolutionary track assuming a constant 𝑀 * . We include as the cyan data point the result from the sphinx simulation (Rosdahl et al. 2017) for 𝑀 halo ∼ 1.5 • 10 9 M ⊙ at 𝑧 = 6. Black data points show estimates for LeoT, Carina and isolated dwarf galaxies in the LG (Read et al. 2017). Baryonic conversion efficiency changes considerably depending on the adopted physical models. The 'full-physics' simulations at 𝑧 = 3.5 provide the best match to LG dwarf galaxies observed at 𝑧 = 0, the sphinx population and evolve along the scalings of Behroozi et al. (2013) and Read et al. (2017) (both for 𝑧 = 0) after 𝑧 ∼ 6.

z = 10 z = 8 z = 7 z = 6 z = 5 z = 3.5 M ha
�� � �� �� -� �� -� �� -� �� � � ���� ( � ☉ ) � * /(� � � ���� ) ���� ����+��� ����+���� ����+������� �� ��+����� ��+���� ��+��������� ���� ��� ���� �� ������ ������ ������ ���������� ��������
𝑀 * (𝑧 = 3.5) ∼ 1.0 • 10 7 M ⊙ to 𝑀 * (𝑧 = 0.5) ∼ 1.3 • 10 7 M ⊙ for HD, and 𝑀 * (𝑧 = 3.5) ∼ 2.2 • 10 6 M ⊙ to 𝑀 * (𝑧 = 0.5) ∼ 2.3 • 10 6 M ⊙ for
HD+Boost. Consequently, our dwarf galaxy is efficiently quenched by 𝑧 ∼ 3.5. We therefore expect most of our results at 𝑧 ∼ 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 𝑧 = 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 𝑧 ≳ 14; Martin-Alvarez et al. 2020;Koh et al. 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 𝑧 ∼ 10 depends thereby on our setup 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 𝑧 ∼ 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 𝑧 ∼ 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 𝑧 = 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 (𝑧 ≳ 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 𝑧 ∼ 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 𝑧 ∼ 7 -8 triggers an important peak of star formation at 𝑧 ∼ 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 𝑧 ∼ 5. The positive feedback of the stellar radiation observed in the RT simulation disappears in the RTiMHD model, where ∼ 𝜇G magnetic fields are present in the galaxy from 𝑧 ∼ 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 𝑧 ≳ 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 𝑧 ∼ 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 𝑧 ∼ 6.3 in these simulations. RTCRiMHD and RTnsCRiMHD maintain instead a low baryonic conversion efficiency at 𝑧 ≳ 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 𝑧 ≳ 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 setups 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 𝑧 = 6 to 𝑧 = 4. After 𝑧 ∼ 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 𝑧 ∼ 5. However, both models converge to a similar stellar mass at 𝑧 ∼ 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% 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). Firstly, we note that these relations (except Moster et al. 2013) correspond to 𝑧 = 0 galaxies, whereas our simulations, due to their expensive computational cost, only reach 𝑧 ∼ 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 𝑀 * (𝑧 = 3.5) ∼ 10 7 M ⊙ , whereas those with boosted SN feedback have 𝑀 * (𝑧 = 3.5) ∼ 2 • 10 6 M ⊙ . 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). Nonetheless, 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 𝐵 0 < 3 • 10 -11 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 𝑀 * /( 𝑓 𝑏 𝑀 halo ) vs 𝑀 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.

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 𝑧 = (7, 6, 5) ± 0.5 and 𝑧 ∈ [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

𝑀 dyn,i = 𝜎 2 v,i (2𝑅 1/2,i ) 𝐺 , (2)
where 𝜎 v,i is the mass-weighted 3D velocity dispersion and 𝑅 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 𝑟 halo sphere centred on the galaxy. The 𝑀 dyn, * vs 𝑀 * 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 𝑀 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 𝑧 = 0.5, they undergo a negligible increase of 𝑀 * and a small change in 𝑀 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 et al. 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 𝑧 ∼ 3.5, whereas the observations correspond to 𝑧 ∼ 0. Therefore to further reinforce our analysis, we include all simulation models that reach 𝑧 = 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 this data.

The 𝑅 1/2,* versus 𝑀 * 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) 3 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 𝑧 ∼ 0.5. While inclusion of magnetic fields tends to slightly increase our 𝑅 1/2,* compared with their non-MHD counterparts, for stronger primordial fields than those explored here, Sanati et al. (in prep) find a trend towards a mild size increase with increasing 𝐵 0 , 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 We only show a subset of our simulations: the representative models that differ significantly from the HD run, with segmented lines representing their evolution from redshift 𝑧 = 7 down to 𝑧 = 3.5 (see text for more details). Additional physical processes affect non-linearly the dynamical mass and galaxy size, with most notable 'outliers' being runs with strong primordial field (both with and without radiation) and boosted SN simulations. In general, simulations with radiation lead to more gas rich and extended galaxies. Overall, a number of simulations provide a good match to recent observations of isolated dwarfs, which are reasonable analogues to our simulated system. 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 𝑅 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.

� �� � �� � �� � �� � �� � � ����* (� ☉ ) � * (� ☉ ) ��� �����+���� �� � �� � �� � �� � �� � � ��� (� ☉ ) � ������� (� ☉ ) Behroozi+2013 Read+2017 �+������+���� ����������� ���� �������+���� �����+���� �� � �� � �� � �� � �� -� �� � � * (� ☉ ) � �/��* (���) ��� ������� ������� ����������+���� �����+���� �� � �� � �� � �� -� �� � � ��� (� ☉ ) � �/����� (���) ���� �� ��+����� ���� ���� �� ������ ������ ������ ���������� ��������
Furthermore, cosmic rays somewhat reduce 𝑅 1/2,* , particularly in the absence of radiation (CRiMHD). This is due to the nonthermal 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 𝑅 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 𝑅 1/2,* at 𝑧 ∼ 3.5, which we attribute to its recent merger event4 (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 𝑀 gas and 𝑀 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 𝑀 dyn,gas while maintaining a comparable gas content to the fiducial HD case. As we evolve HD and HD+Boost to 𝑧 = 0.5, their dynamical masses measured through the gaseous component remain relatively unchanged, and their gas masses remain around 𝑀 gas ∼ 2 • 10 7 M ⊙ . Simulations with magnetic fields do not sig-nificantly 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 𝑀 dyn,gas /𝑀 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 𝑧 = 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 𝑧 ∼ 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 𝑅 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 5 , and that our system is not evolved to 𝑧 = 0.

Integrated and resolved kinematics of stars and gas

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) 6 . We include evolutionary tracks sampling median values during 𝑧 = (7, 6, 5) ± 0.5 and 𝑧 ∈ [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 2𝑅 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 order ∼ 1 -5 km s -1 , and become somewhat lower when the strength of the SN feedback is boosted. During their time evolution to 𝑧 = 0.5, HD and

� �� � �� � � ��� (�� � -� ) � * (� ☉ ) ��� �� ��� ��� � �������+���� (��/���/��) �������+���� (���) �������+���� (��) � �� �� �� �� σ � (�� � -� ) ��� �� ��� ��� � �� � �� � �� � �� -� �� � � * (� ☉ ) � ��� / σ � �� ��+����� ��+���� ��+��������� ���� ��� ���� �� ������ ������ ������ ���������� �������� Figure 5
. Integrated stellar kinematic properties versus 𝑀 * for our lated galaxy at 𝑧 ∼ 3.5, compared against available observations of local dwarf galaxies (Wheeler et al. 2017) at 𝑧 = 0. From top to bottom, panels show mass-weighted rotational velocity 𝑣 rot , velocity dispersion 𝜎 v,* and 𝑣 rot /𝜎 v,* , all measured within 2𝑅 1/2,* . Vertical shaded bands show, as in Fig. 4, the expected stellar mass range for our simulated dwarf galaxy based on abundance matching. We only show evolution tracks for some representative models, displaying the evolution from 𝑧 = 7 until 𝑧 = 3.5. Most of our simulations attain 𝜎 v,* values generally located at the upper end of the observations distribution or above it. Amongst them, our models with cosmic rays have the lowest values of 𝜎 v,* at a given 𝑀 * .

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 ∼7 -10 km 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 𝑧 = 0.5, the stellar velocity dispersion of HD is reduced by ∼ 5 km 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 𝑀 * . 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 𝑧 ∼ 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+Boost to 𝑧 = 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.

HI properties

Including radiative transfer leads to a more self-consistent modelling of the HI content in our dwarf galaxy. This allows us to explore the subset of radiation-hydrodynamics simulations and compare their HI properties to observations. Fig. 6 revisits the kinematic properties of the studied galaxy, but now for the HI component, comparing with HI observations for similar mass galaxies. We compute each of these quantities adopting the same procedure as done for Figs. 4 and5, except for the maximal rotational velocity max(𝑣 rot,HI ), which corresponds to the maximum value of the rotational component of the velocity profiles within 2𝑅 1/2,HI . We compare our models at 𝑧 ∼ 3.5 with HI observational results for similar galaxies at 𝑧 ∼ 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 HI 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 7 . Our measurements of the dynamical mass using HI (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(𝑣 rot,HI ) (second panel), with the simulations combining radiative transfer and MHD falling slightly below the observations. For the HI velocity dispersion 𝜎 v,HI (third panel), all our simulations have values of 𝜎 v,HI ∼ 10 km s -1 and the same approximate scatter, providing a good match to the 7 We have quantified the changes in the kinematic quantities of the HD and HD+Boost models between 𝑧 ∼ 3.5 and 𝑧 = 0.5. Note however that this measurement is done for the gas component: with no radiative transfer, these simulations do not provide a self-consistent modelling of hydrogen ionization. Keeping this caveat in mind, we find that their maximal rotational velocity remains within ∼ 20 -30 km s -1 , the 𝜎 v,gas within ∼ 5 -10 km s -1 , and their ratio remains at ∼ 3. The evolution of 𝑀 dyn,gas is reviewed in Fig. 4. distribution of observations. Finally, the ratio of max(𝑣 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 . Nonetheless, all our radiative transfer simulations predict reasonable HI content and properties of our simulated dwarf galaxy.

������ �������� ��� � ��� � ��� ��� � ��� � ��� ��� ��� ��� ������ �������� ����������� �����+���� �������+���� �� � �� � �� � �� � �� � �� � �� � �� � � ������ (� ☉ ) � �� (� ☉ ) �����+���� �������+���� �� � �� � ���[� ������ (�� � -� )] ����������� ���� �����+���� �����+���� �������+���� �� � σ ���� (�� � -� ) ��� ��� ��� � ������ ��� �������� �����+���� �������+���� �� � �� � �� � �� � �� � �� � �� � � �� (� ☉ ) ���[� ������ ]/σ ����

Comparing HI 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 HI is dominated by rotation (𝑚𝑎𝑥 [v rot,HI ]/𝜎 v,HI > 1), in agreement with their col-

�� -� �� � ��� ��� ��� ��� ��� ��� � ����� / σ � � ��� / � ����� Figure 7.
Ratio of the dynamical mass with respect to the total mass, inferred through the stellar (standard symbols) and HI (lighter symbols) components.

The two estimates of the dynamical mass encompass well the total mass, except for the model with an extreme primordial magnetic field (RTsMHD).

lisional (HI) and collisionless (stellar) nature. Rotational velocities for the stars are ≲ 10 km s -1 whereas the maximum HI rotational velocity tends to be within 10 -30 km s -1 . In terms of velocity dispersion, HI displays 𝜎 v,HI ≲ 10 km s -1 , whereas 𝜎 v,* tends to be within 14-20 km 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 (𝑀 dyn, * ; standard symbols) and the gaseous dynamical mass (𝑀 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 HI 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 𝑀 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.

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 𝑧 = 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 𝜌 2 weighting, to more closely mimic the emission-weighted maps. For the velocity projections, we display only gas cells with densities 𝜌 > 0.1m H 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 𝑧 ∼ 5 to 𝑧 ∼ 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 HI 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.

Colour-magnitude relation

As part of our analysis, we review the optical emission properties of our simulated dwarf galaxy. Fig. 9 shows the 𝑔 -𝑖 colour versus absolute magnitude 𝐺 (top) and versus stellar mass within 2 𝑅 1/2,* (bottom), compared with the observed properties of the ELVES dwarf galaxies (Carlsten et al. 2021(Carlsten et al. , 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 [START_REF] Bruzual | [END_REF] 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 𝑧 = 3.5 is shown with shaded symbols. Since for the majority of our simulations we only evolve the dwarf galaxies to 𝑧 = 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 𝑧 = 0.5 (i.e. we increase all particle ages by ∼ 6.8 Gyr) in their colour computation. The selected value of 𝑧 = 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 𝑧 = 0 would increase the 𝑔 -𝑖 colour by an additional (0.04 ± 0.02) mag 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 𝑧 = 0.5 for the simulations that reached that redshift (i.e. for the HD and HD+Boost models, at 𝑧 = 0.5). These are shown as encircled datapoints. We compare this 'manual' ageing with the actual simulation data for the HD and HD+Boost runs which reached 𝑧 = 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 ∼ 0.05 mag AB and 0.05 mag AB , respectively. These estimates experience larger deviations for the ab-solute magnitude in HD+Boost, with only Δ(𝑔 -𝑖) ∼ 0.05 mag AB but Δ𝐺 ∼ 0.5 mag AB .

We conclude that the 'passive ageing' procedure leads to reasonable estimates of the observed 𝑔 -𝑖 colour, with a possible overestimation of the absolute magnitude 𝐺 of up to ∼ 1 mag 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. Nonetheless, 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 𝑧 ∼ 3.5, our galaxies are relatively blue compared with the overall observed ELVES population. By 𝑧 ∼ 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 Observational data:

Carlsten+2021 Carlsten+2022 z=0.5 z=0.5 �� � �� � ��� ��� ��� ��� ��� � * (� ☉ ) � -� (��� �� ) �� ��+����� ��+���� ��+��������� ���� ��� ���� �� ������ ������ ������ ���������� �������� Figure 9
. 𝑔 -𝑖 colour vs 𝐺 absolute magnitude (top) and stellar mass (bottom) for the simulated dwarf galaxies accounting for dust obscuration compared with ELVES dwarfs at 𝑧 ∼ 0 (Carlsten et al. 2021(Carlsten et al. , 2022)). These observations are divided into 'Early type' (red), 'Late type' (blue) and those with an ambiguous classification (green). Points show the galaxy at 𝑧 = 3.5 (shaded symbols) and the colour upper limit estimate for 𝑧 ∼ 0.5 (standard symbols; see text for details), respectively. We also include as encircled symbols the exact measurements at 𝑧 ∼ 0.5 for the HD and HD+Boost simulations, which we evolved down to 𝑧 = 0.5. The coloured bands in the bottom panel indicate the range of stellar masses predicted from the abundance matching analysis as labelled.

have assumed that they remain completely quenched from 𝑧 = 3.5 to 𝑧 = 0.5, hence our 𝑔 -𝑖 colours should be considered only as a conservative upper limit, and the stellar masses and absolute 𝐺 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 reignition 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.

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 re-lations between various HI properties and magnetic fields in local dwarf galaxies by Chyzy et al. (2011). Our average magnetic field at galactic scales is on the order of ∼ 2 -5 𝜇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 HI 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 𝑀 HI values, where these observations suggest 𝐵 ∼ 3 -15 𝜇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 farinfrared (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 magnetisation 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 polaris8 (Reissl et al. 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. [START_REF] Martin-Alvarez | Monthly Notices of the Royal Astronomical Society Martin N[END_REF] show the peak of the magnetic energy spectrum 𝑘 peak to scale with the finest simulation resolution Δ𝑥 min following 𝑘 peak ∝ Δ𝑥 -1 min . Körtgen et al. (2017) suggest resolutions Δ𝑥 min < 0.1 pc are required to approach turbulence convergence, finding the average turbulence to stabilize at values of Δ𝑥 min ≲ 0.03 pc. Extrapolating the [START_REF] Martin-Alvarez | Monthly Notices of the Royal Astronomical Society Martin N[END_REF] scaling of 𝑘 peak from our suite resolution Δ𝑥 min,sim ∼ 7 pc to Δ𝑥 min,exp = 0.005 pc yields a ratio of 𝑘 peak,exp /𝑘 peak,sim ∼ 40.1. Expanding a Kazantsev-like spectrum (∝ 𝑘 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) spectral index by 𝛼 = ( 𝑝 -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 𝑒 𝑒 -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 degrees) to the local linear polarisation (accounting for Faraday depolarisation) in order to align with local magnetic field. These dashes are only shown in regions where the linearly polarised intensity is higher than 5% 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 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 polarisation 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 polarisation 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 polarisation vectors appears more coherent in the models combining radiative transfer and cosmic rays, particularly at scales of the VLA resolution.

Halo density profiles

Dark matter density profiles in dwarf galaxies estimated from observations frequently suggest shallower central slopes than that of We follow Read et al. (2019) and include bands for the expected density at 𝑧 = 3 (shaded bands) and 𝑧 = 0 (dashed bands). The latter provide a good approximation to the bands expected at 𝑧 = 0.5, to be compared with the encircled symbols. We only show evolutionary tracks for a number of representative models, with coloured lines displaying the evolution from 𝑧 = 7 until 𝑧 = 3.5, and symbols showing the simulated data at 𝑧 = 3.5.

Simulations with very high central densities are indicated with arrowheads at the top of the graph, separated by artificial displacements along 𝑥 axis for visibility. Following Read et al. (2019), we colour the observational points corresponding to galaxies that stopped forming stars more than 6 Gyr ago as black, between 6 and 3 Gyr as purple, and less than 3 Gyr ago as blue. Colour bands show the prediction by Read et al. (2019) for the central density in haloes with cusps (black) or cores (blue). While the majority of our runs have cuspy central dark matter distributions, some of our simulations with radiation and cosmic rays fall within the 'core' band.

the theoretical expectation from dark matter assembly in a cold dark matter cosmology (Navarro et al. 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 susceptible to the SN-driven outflows (Read et al. 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 ∼ 0.3 • 10 8 M ⊙ 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 𝑟 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 𝑧 = 3 (shaded bands) and 𝑧 = 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 𝑧 ∼ 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 𝑧 ∼ 4. When evolved to 𝑧 = 0.5, both the HD and HD+Boost simulations remain within the 𝑧 = 0 cuspy halo band (dashed black lines), which is approximately equal to the band at 𝑧 = 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 𝑧 ∼ 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 𝑧 = 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 𝑧 ≲ 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 𝑧 = 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 𝑧 ∼ 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 𝑧 ∼ 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. Nonetheless, we note that the final central density measured for RTnsCRiMHD increases by a non-negligible amount when measured using the dark matter centre (reaching ∼ 1.5 • 10 8 M ⊙ ).

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+Boost and RTiMHD) at 𝑧 = 3.5. The vertical solid line indicates 𝑟 = 0.15 kpc, used in Fig. 11 to calculate the average central dark matter density. Dotted lines spanning from 𝑟 = 0.15 kpc up to 𝑟 = 2 kpc show the NFW density scalings 𝜌 DM ∝ 𝑟 -1 (black) and 𝜌 DM ∝ 𝑟 -3 (gray). Gray and blue bands within 𝑟 = 0.15 kpc denote the values of 'cuspy' and 'cored' dwarf profiles, respectively. These are computed as in Fig. 11, and are shown at 𝑧 = 3 (shaded bands) and 𝑧 = 0 (dashed horizontal lines). The 𝑦 range shown corresponding to the halo mass that matches our simulated dwarfs at 𝑧 = 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 𝑧 = 5 to 𝑧 = 3.5 for the 'fullphysics' RTCRiMHD and RTnsCRiMHD simulations, respectively. Here, focusing first on the dark matter density profiles, the central profile is somewhat shallower than the 𝜌 DM ∝ 𝑟 -1 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 𝑧 ∼ 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 𝑧 = 3.5 in the RTCRiMHD simulation). This is caused by efficient cosmic ray-driven outflows, pushing the gas out of our simulated dwarfs.

CONCLUSIONS

In this paper, we introduce the Pandora suite of high resolution cosmological zoom-in simulations of a dwarf galaxy (with halo mass 𝑀 vir (𝑧 = 0) ≈ 10 10 M ⊙ ) 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, a magneto-thermo-turbulent 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 𝑧 ∼ 4, as it resides in a field environment. Consequently, we evolve all our simulations, summarised in Table 1, down to 𝑧 ∼ 3.5. We also evolve a small subset of simulations all the way to 𝑧 = 0.5 confirming the largely passive evolution since 𝑧 ∼ 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:

(i) 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.

(ii) In terms of the overall morphology and spatial gas distribution, models with strong SN feedback produce an over-quenched, 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 these two extremes, due to cosmic ray-driven outflows that deplete the gas reservoir in the ISM and can more effectively quench star formation.

(iii) 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.

(iv) While a number of models with stellar radiation lead to integrated rotational velocities and velocity dispersion for stars and HI 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.

(v) 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 [START_REF] Martin-Alvarez | Monthly Notices of the Royal Astronomical Society Martin N[END_REF] or subgrid turbulence models (Semenov et al. 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. 2019Koudmani 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 over-quenched 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 multi-phase 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.

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. 𝑧 = 3.5), we show in Fig. A1 the halo mass growth. This is shown for the HD simulation until 𝑧 = 0.5, and down to 𝑧 = 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 𝑧 = 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 𝑧 = 0.5 (see Fig. A2). While there is some progressive halo growth during the 𝑧 = 3.5 to 𝑧 = 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 VS 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 ∼ 0.3 • 10 8 M ⊙ 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 centering.

The

  Pandora project. I: the impact of radiation and cosmic rays on baryonic and dark matter properties of dwarf galaxies Sergio Martin-Alvarez, Debora Sijacki, Martin G. Haehnelt, Marion Farcy, Yohan Dubois, Vasily Belokurov, Joakim Rosdahl, Enrique Lopez-Rodriguez To cite this version: Sergio Martin-Alvarez, Debora Sijacki, Martin G. Haehnelt, Marion Farcy, Yohan Dubois, et al.. The Pandora project. I: the impact of radiation and cosmic rays on baryonic and dark matter properties of dwarf galaxies. Monthly Notices of the Royal Astronomical Society, 2023, 525 (3), pp.3806-3830. 10.1093/mnras/stad2559. hal-03892311 1 INTRODUCTION

Figure 1 .

 1 Figure 1. 'Full-physics' simulation (RTCRiMHD) projections centred on the dwarf galaxy at 𝑧 = 3.5. (Inset box) Mass projection of the entire simulated box, with the zoom-in circle encompassing approximately twice the virial radius of the galaxy studied. (Large panel) RGB projections of a (20 kpc) 3 box showing gas density (blue), gas temperature (orange) and stellar density (yellow). Square labels indicate the projected regions shown in the bottom row (7.5 kpc) and the right-hand column (2.5 kpc). The circular panel in the central zoom region shows a gas density close-up view of the inner 500 pc of the dwarf galaxy.(Right-hand panels) Synthetic observations of the galaxy in a (2.5 kpc) 3 box. The top panel shows radio synchrotron (𝜆 = 6.2 cm, dashes show 90º-rotated polarisation to align with the magnetic field) as would be observed by VLA in its D configuration in the top right corner and for SKA-like resolutions at the bottom left. The middle and bottom panels, show respectively the optical emission as would be observed by SDSS, and the near infrared as would be observed by JWST. For these synthetic observations we artificially position the galaxy at a distance of 2 Mpc from Earth in order to emulate Local Group dwarf galaxies analogue distances. Each image is convolved with the PSF of the corresponding telescope. For the optical and near-IR panels, we provide the apparent magnitude within the circular aperture displayed by the white circle wedges. (Bottom row) From left to right, these four panels are projections of a (7.5 kpc) 3 box showing magnetic field, hydrogen ionization fraction, cosmic ray energy density, and gas density (gray) separated into inflowing (blue gas; 𝑣 gas,radial < -10 km s -1 ) and outflowing (red gas; 𝑣 gas,radial > 10 km s -1 ) gas. The galaxy in our full-physics model appears particularly extended (compared to the fiducial hydrodynamical case, discussed in Fig.2), with an envelope of outflowing gas that correlates spatially with the high energy density region of cosmic rays and strong magnetic fields, extending to approximately 3 kpc.

Figure 2 .

 2 Figure2. Face-on projection centred on the studied dwarf galaxy at 𝑧 = 3.5 for different simulation models as indicated on the panels. Panels show in each column gas density (top) and a synthetic optical observation (bottom). The synthetic observations employ the SDSS [u,g,r] filters in the galaxy rest-frame, accounting for dust extinction but not convolved with any telescope PSF. We reduce the luminosity of the NoFb images by 1 dex to avoid saturated images.

  Hydro Hydro Hydro Hydro Hydro Hydro Hydro Hydro

  th w a + 2 0 1 8 M os te r+ 20 13 (z = 3. 5) Be hro ozi +20 13 (z = 0) Re ad +2 01 7 z = 10 z = 8 z = 7 z = 6 z = 5 z = 3.5 M ha lo g ro w th a t c o n s ta n t M * ������ (� = �� �������+��) ����+���� MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD MHD

Figure 3 .

 3 Figure3. Baryonic conversion efficiency (𝑀 * /( 𝑓 𝑏 𝑀 halo )) versus halo mass 𝑀 halo for all our simulations of the dwarf galaxy studied here. Simulations are separated into hydrodynamical simulations comparing star formation and SN feedback models (Hydro, top left), MHD (MHD, top right), radiative transfer with and without MHD (RT, bottom left) and CRMHD with and without radiative transfer (CRs, bottom right), with the HD model included in all panels as reference. For each simulation, we show the conversion efficiency at several redshifts. For the NoFb, HD, and HD+Boost simulations, we include four additional evolution markers corresponding to redshifts 𝑧 = 0.5, 1, 2, 3 moving from rightmost towards the left of the panel. These illustrate how simulations including SN feedback have negligible star formation between 𝑧 = 3.5 and 𝑧 = 0.5. The blue band shows an extrapolation of the baryonic conversion efficiency predicted byBehroozi et al. (2013) for 𝑧 = 0, the yellow band indicates the 1𝜎 prediction byMoster et al. (2013) at 𝑧 = 3.5 and the green band corresponds to the abundance matching prediction fromRead et al. (2017). We also include the stellar mass-halo mass relation predicted byJethwa et al. (2018) for the Milky Way satellites population as the red band in the bottom right corner. We provide additional details regarding their comparison and about additional relations at the end of Section 3.2. The dotted black line indicates an evolutionary track assuming a constant 𝑀 * . We include as the cyan data point the result from the sphinx simulation(Rosdahl et al. 2017) for 𝑀 halo ∼ 1.5 • 10 9 M ⊙ at 𝑧 = 6. Black data points show estimates for LeoT, Carina and isolated dwarf galaxies in the LG(Read et al. 2017). Baryonic conversion efficiency changes considerably depending on the adopted physical models. The 'full-physics' simulations at 𝑧 = 3.5 provide the best match to LG dwarf galaxies observed at 𝑧 = 0, the sphinx population and evolve along the scalings ofBehroozi et al. (2013) andRead et al. (2017) (both for 𝑧 = 0) after 𝑧 ∼ 6.

Figure 4 .

 4 Figure 4. Dynamical mass (top row) and half-mass radius (bottom row) as a function of stellar mass (left column) and gas mass (right column). Dynamical masses are estimated from stars (left) and gas (right), and computed according to Equation (2). Data from different observations at 𝑧 = 0 is included for comparison (see legends). In order to aid the comparison of our 𝑧 ∼ 3.5 results with observations, we include 𝑧 = 0.5 values where available as encircled symbols of the corresponding models. Dashed gray lines (top row) correspond to the one-to-one relation between dynamical and component mass. On the left column, vertical shaded bands indicate the stellar mass range predicted for the simulated galaxy halo mass (𝑀 halo (𝑧 = 3.5) ∼ 4.1 • 10 9 M ⊙ ) according toBehroozi et al. (2013, shown in blue) andRead et al. (2017, shown in green). We only show a subset of our simulations: the representative models that differ significantly from the HD run, with segmented lines representing their evolution from redshift 𝑧 = 7 down to 𝑧 = 3.5 (see text for more details). Additional physical processes affect non-linearly the dynamical mass and galaxy size, with most notable 'outliers' being runs with strong primordial field (both with and without radiation) and boosted SN simulations. In general, simulations with radiation lead to more gas rich and extended galaxies. Overall, a number of simulations provide a good match to recent observations of isolated dwarfs, which are reasonable analogues to our simulated system.

5

  Masses are 𝑀 * ∼ 7.6 • 10 7 M ⊙ for NGC6822 (Read et al. 2017)) and 𝑀 * ∼ 7.6 • 10 7 M ⊙ for NGC3109 (McConnachie 2012), while our dwarf galaxy has 𝑀 * ∼ 10 7 M ⊙ . 6 Observational source manuscripts are: Mateo (1998); Koch et al. (2007); Simon & Geha (2007); Walker et al. (2009); Fraternali et al. (2009); Leaman et al. (2009); Geha et al. (2010); Koposov et al. (2011); Leaman et al. (2012); Tollerud et al. (2012); McConnachie (2012); Frinchaboy et al. (2012); Ho et al. (2012); Collins et al. (2013); Martin et al. (2013); Kirby et al. (2014); Martin et al. (2014); Salomon et al. (2015); Walker et al. (2015).

Figure 6 .

 6 Figure 6. HI dynamical mass and kinematics vs 𝑀 HI for our simulations featuring radiative transfer (i.e. self-consistently modelling HI). From top to bottom, panels display the dynamical mass estimate using HI (𝑀 dyn,HI ; Equation (2)), maximum HI rotational velocity (max(𝑣 rot,HI )), HI velocity dispersion (𝜎 v,HI ) and the max(𝑣 rot,HI )/𝜎 v,HI ratio. Error bars show median quartiles during the redshift interval 𝑧 ∈ [3.4, 3.7]. Overall, all our radiative transfer simulations models at 𝑧 = 3.5 are in good agreement with HI observations at 𝑧 ∼ 0, with the RT, RTnsCRiMHD and RTCRiMHD simulations providing the best match to the local isolated dwarfs reported by Kirby et al. (2017).

Figure 8 .

 8 Figure 8. For the subset of simulations shown (as labelled), the left column displays line-of-sight stellar velocity dispersion (top row) and line-of-sight velocities for stars (second row), gas (third row), neutral hydrogen (fourth row) and ionized hydrogen (bottom row). Panels in right column show a synthetic SDSS mock image (top row), and the corresponding density maps next to each line-of-sight velocity field. All projections are computed along the x-axis of the box. Only gas cells with total densities higher than 0.1 m H cm -3 are employed in the line-of-sight computation.

  Figure10. (Leftmost panels) From top to bottom, the 𝑦-axis of various panels displays maximum rotational velocity of HI, max(𝑣 rot,HI ), total HI mass, 𝑀 HI , HI velocity dispersion, 𝜎 v,HI , and the star formation rate averaged over the last 100 Myr, SFR 100 Myr versus the average magnetic field strength 𝐵. Our simulated dwarf galaxies have magnetic fields at 𝑧 ∼ 3.5 that are in agreement with observations byChyzy et al. (2011) of local dwarf galaxies at 𝑧 ∼ 0. (Right panels) Comparison of the synchrotron synthetic views (for 𝜆 = 6.2 cm) for our simulated dwarf galaxy at 𝑧 = 3.5 generated using polaris (see text for details). We assume the galaxy to be observed at a distance of 2 Mpc, convolving the images with the corresponding telescope resolution. White dashes display the local linear polarisation orientation rotated 90 degrees to align with the local magnetic field. Each row displays from left to right: VLA-like synthetic map, SKA-like synthetic map and magnetic field strength (streamlines represent magnetic field lines). From top to bottom, rows show the RTsMHD, RTiMHD, RTnsCRiMHD and RTCRiMHD simulations. Emission traces regions of stronger magnetic field. Due to the turbulent nature of magnetic field lines, polarisation only captures large-scale magnetic field orientation features. Our simulations with cosmic rays and radiative transfer have a more concentrated synchrotron emission.

Figure 11 .

 11 Figure11. Average dark matter density within the inner 150 pc of the simulated dwarf galaxy versus 𝑀 halo , where the models HD and HD+Boost are also shown at their later redshift 𝑧 = 0.5 with encircled symbols. Our results are compared with Local Group dwarf galaxies data byRead et al. (2019). We followRead et al. (2019) and include bands for the expected density at 𝑧 = 3 (shaded bands) and 𝑧 = 0 (dashed bands). The latter provide a good approximation to the bands expected at 𝑧 = 0.5, to be compared with the encircled symbols. We only show evolutionary tracks for a number of representative models, with coloured lines displaying the evolution from 𝑧 = 7 until 𝑧 = 3.5, and symbols showing the simulated data at 𝑧 = 3.5. Simulations with very high central densities are indicated with arrowheads at the top of the graph, separated by artificial displacements along 𝑥 axis for visibility. FollowingRead et al. (2019), we colour the observational points corresponding to galaxies that stopped forming stars more than 6 Gyr ago as black, between 6 and 3 Gyr as purple, and less than 3 Gyr ago as blue. Colour bands show the prediction byRead et al. (2019) for the central density in haloes with cusps (black) or cores (blue). While the majority of our runs have cuspy central dark matter distributions, some of our simulations with radiation and cosmic rays fall within the 'core' band.

Figure 12 .

 12 Figure12. Radial density profiles of the gaseous (dashed), stellar (dotted) and dark matter (solid) components for various representative models. The vertical solid line indicates 𝑟 = 0.15 kpc, used in Fig.11andRead et al. (2019). Colour bands show the prediction for the central density in haloes with cusps (gray) or cores (blue) for 𝑀 halo = 4.5 • 10 9 M ⊙ at 𝑧 = 3 (shaded bands) and 𝑧 = 0 (dashed lines) computed in an analogous manner to that of(Read et al. 2019). Finally, we include two additional lines spanning from 𝑟 = 0.15 kpc up to 𝑟 = 2 kpc showing the NFW density scalings, 𝜌 DM ∝ 𝑟 -1 (gray) and 𝜌 DM ∝ 𝑟 -3 (black). (Left column) From top to bottom, the panels show NoFb, HD, HD+Boost and RTiMHDat 𝑧 = 3.5, illustrating how SN feedback and radiation decrease the central dark matter density and flatten the radial slope of the dark matter profile. (Central column) Radial density profile evolution for the full-physics RTCRiMHD simulation from 𝑧 = 5 (top) to 𝑧 = 3.5 (bottom). (Right column) Same as the central column but for the RTnsCRiMHD simulation. These two columns illustrate that cosmic ray feedback together with radiation lead to a central dark matter profile more in line with 'cores', which is transformed back to a cusp at 𝑧 ∼ 3.5 in RTCRiMHD, due to a 'wet merger'.

Figure A1 .

 A1 Figure A1. Pandora halo mass growth down to 𝑧 = 0 for the HD simulation parsed with the halo mass growth from a dark matter-only simulation beyond 𝑧 = 0.5. We overplot all major (blue circles) and minor (red circle) mergers experienced by the halo. The studied halo forms rapidly at high redfshift, with no mergers beyond 𝑧 ∼ 4. After that time, the halo only grows in a gradual and secular manner.

Figure A2 .Figure B1 .

 A2B1 Figure A2. Gas density evolution of the studied halo in the HD simulation. No significant gas accretion nor encounters with gaseous structures occur between 𝑧 = 3.5 and 𝑧 = 0.5 due to the considerable spatial isolation of the simulated galaxy.

MNRAS 000,1-26 (2023)

See e.g.Hopkins (2016) for a performance comparison of various MHD solver methods in simple and complex astrophysical problems.MNRAS 000,1-26 (2023)

Note that our boosted SN feedback models do not include radiative transfer, whereas the sphinx simulations do.MNRAS 000,1-26 (2023)

Note that here observational data shows the half-light radius rather than the half-mass radius we measure in our simulations, and that observational measurements are projected.MNRAS 000,1-26 (2023)

Note that this merger event takes place at slightly different times in each of our models, featuring also somewhat different mass ratios and impact parameters.MNRAS 000,1-26 (2023)

https://portia.astrophysik.uni-kiel.de/polaris/ MNRAS 000,1-26 (2023)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-26 (2023)

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

DATA AVAILABILITY

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

DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).