Indranil Banik

Ingo Thies

Benoit Famaey

Graeme Candlish

Pavel Kroupa

Typeset using L A T E X twocolumn style in AASTeX63

Keywords:

The dynamical stability of disk galaxies is sensitive to whether their anomalous rotation curves are caused by dark matter halos or Milgromian Dynamics (MOND). We investigate this by setting up a MOND model of M33. We first simulate it in isolation for 6 Gyr, starting from an initial good match to the rotation curve (RC). Too large a bar and bulge form when the gas is too hot, but this is avoided by reducing the gas temperature. A strong bar still forms in 1 Gyr, but rapidly weakens and becomes consistent with the observed weak bar. Previous work showed this to be challenging in Newtonian models with a live dark matter halo, which developed strong bars. The bar pattern speed implies a realistic corotation radius of 3 kpc. However, the RC still rises too steeply, and the central line of sight velocity dispersion (LOSVD) is too high. We then add a constant external acceleration field of 8.4 × 10 -12 m/s 2 at 30 • to the disk as a first order estimate for the gravity exerted by M31. This suppresses buildup of material at the centre, causing the RC to rise more slowly and reducing the central LOSVD. Overall, this simulation bears good resemblance to several global properties of M33, and highlights the importance of including even a weak external field on the stability and evolution of disk galaxies. Further simulations with a time-varying external field, modeling the full orbit of M33, will be needed to confirm its resemblance to observations.

INTRODUCTION

The ubiquity of thin disk galaxies belies their central role in a major astronomical mystery -why do their rotation curves (RCs) go flat in their outskirts instead of following the expected Keplerian decline beyond the extent of their luminous matter (e.g. Babcock 1939;Rubin & Ford 1970;Rogstad & Shostak 1972;Roberts & Whitehurst 1975;Bosma 1981)? These acceleration discrepancies are historically related to the issue of disk galaxy stability. Hohl (1971) showed that a selfgravitating Newtonian disk is very unstable. Given the age of galaxy disks like that of our own Milky Way (Knox et al. 1999), it is clear that something fundamental is missing from these simulations.

ibanik@astro.uni-bonn.de * Humboldt fellow Ostriker & Peebles (1973) suggested the addition of a dominant dark matter (DM) halo around the disk, which can stabilize it. Such DM halos around galaxies are now considered an essential component of the Lambda-Cold Dark Matter (ΛCDM) cosmological paradigm (Efstathiou et al. 1990;Ostriker & Steinhardt 1995). In this model, disks are not self-gravitating, but are mostly held together by a DM halo. The extra gravity from the halo elevates the disk RC, and can make it asymptotically flat in the outskirts. If the total gravity g consists of a halo contribution g h and a disk contribution g d , then the fractional change in g from some disk surface density perturbation is reduced by a factor of (g d /g) due to the halo (Banik et al. 2018a). In this way, the halo mass interior to the radius of the disk can make the disk more stable and elevate its RC.

1.1. MOND arXiv:2011.12293v2 [astro-ph.GA] 1 Dec 2020

However, several problems remain to this day in matching detailed properties of RCs (e.g. Oman et al. 2015;Desmond 2017a,b) and the long-term stability of some disk galaxies (e.g. Sellwood et al. 2019, hereafter S19). Even if they were solved, agreement between theory and observations does not prove a theory correct -there may be alternative explanations for the same data. In particular, the RC anomalies could be a sign that Newtonian dynamics breaks down on galaxy scales. It would after all not be too surprising if a theory developed exclusively using Solar System data cannot be reliably extrapolated by many orders of magnitude in distance and acceleration. The best developed proposal for how a breakdown might occur is Milgromian Dynamics (MOND, Milgrom 1983;Bekenstein & Milgrom 1984;Milgrom 2014). In MOND, the gravitational field strength g at distance R from an isolated point mass M transitions from the Newtonian GM/R 2 law at short range to:

g = GM a 0 R , R GM a 0 . (1)
MOND introduces a 0 as a fundamental acceleration scale of nature below which the deviation from Newtonian dynamics becomes significant. When g a 0 , a system is said to be in the deep-MOND limit (DML) in which the point mass gravity declines only inversely with R, making the dynamics scale invariant (Milgrom 2009). Empirically, a 0 ≈ 1.2×10 -10 m/s 2 to match galaxy RCs (Begeman et al. 1991;Gentile et al. 2011). With this value of a 0 , MOND continues to this day to fit galaxy RCs very well using only their directly observed baryonic matter (e.g. Kroupa et al. 2018a;Li et al. 2018;Sanders 2019). The properties of polar ring galaxies (Lüghausen et al. 2013) and shell galaxies (Bílek et al. 2015) are also consistent with MOND. It has been suggested that MOND is a manifestation of the quantum vacuum, and thus holds important clues on how to unify quantum mechanics with gravity (Milgrom 1999;Verlinde 2016;Smolin 2017).

Equation 1 predicts that galaxy RCs should become asymptotically flat at a plateau with level v f , which depends only on the baryonic mass M of the galaxy according to:

v f = 4 GM a 0 .
(2)

The observational counterpart to this prediction is commonly known as the baryonic Tully-Fisher relation (BTFR, McGaugh et al. 2000), a generalization of the original luminous Tully-Fisher relation (Tully & Fisher 1977). Although this played an important role in the development of MOND, it was certainly not clear in the 1980s that the BTFR would continue to remain very tight once the gas mass was included in low mass galaxies, and more generally once galaxies with rather different properties were observed (McGaugh 2012;Lelli et al. 2016b). MOND can also be applied to systems which do not have spherical symmetry. This requires the use of a generalized Poisson equation, which can be derived from a Lagrangian (Bekenstein & Milgrom 1984;Milgrom 2010). We discuss this further in Section 2.4 since we will need this for our numerical simulations. Unlike the standard Poisson equation, the MOND version is nonlinear, as is readily apparent from Equation 2. The fact that MOND gravity from different sources cannot be superposed creates an external field effect (EFE, Milgrom 1986). As a result, the gravity from surrounding structures induces a weakening of the galaxy's gravitational field, but without affecting its inertial mass. This can cause the RC to decline at large distance (Haghi et al. 2016;Hees et al. 2016). In principle, it is possible to estimate this effect based on the environment of the galaxy in question (Wu & Kroupa 2015). Since large isolated galaxies are typically used for RC analyses, deviations from the BTFR are generally expected to be small. Nonetheless, Chae et al. (2020) recently reported a highly significant detection of the predicted correlation between external gravitational environment and RC deviations from the isolated MOND prediction.

The EFE is more important for dwarfs near a massive host. It is particularly important to explain the diverse velocity dispersions of galaxies such as Dragonfly 2 (DF2, Famaey et al. 2018;Kroupa et al. 2018b), DF4 (Haghi et al. 2019b), or DF44 (Bílek et al. 2019;Haghi et al. 2019a). The predicted dispersions for DF2 and DF4 are significantly affected by the EFE. However, the more isolated DF44 does indeed have a higher internal velocity dispersion despite a similar baryonic content. Some evidence for the EFE was also found in the dwarf satellite galaxies of M31 (McGaugh & Milgrom 2013) and in the Milky Way satellite Crater 2 (McGaugh 2016;Caldwell et al. 2017). Without the EFE, it is not possible to escape from a mass as its potential is logarithmically divergent (Equation 1). Thus, the fact that MOND with the EFE can accommodate the Galactic escape velocity curve also argues in favour of the EFE (Famaey et al. 2007;Banik & Zhao 2018b).

Recently, the second data release of the Gaia mission (Gaia Collaboration 2018) has been used to convincingly demonstrate that MOND cannot function without some form of EFE. The evidence comes from the velocity distribution of wide binary stars in the Solar neighbourhood -at separations 5 kAU, MOND should affect their orbits (Hernandez et al. 2012). The predicted boost to their orbital velocity is rather large without the EFE, but a more modest 20% if the Galactic EFE is included (Banik & Zhao 2018a). Observations of wide binaries completely rule out MOND without the EFE (Pittordis & Sutherland 2019). The more subtle effect with the EFE is still allowed observationally, and should be testable with future data releases and/or further analysis. The main source of contamination is probably hierarchical triples, which should be accounted for statistically (Belokurov et al. 2020;Clarke 2020).

In addition to these small scale tests, MOND might also be strongly constrained by the success of ΛCDM in fitting large scale observations like the cosmic microwave background (CMB, Planck Collaboration VI 2020). However, the excellent fit they obtain is based on an expansion rate history that requires a local Hubble constant below the observed value (e.g. Riess 2020). This could be due to our position within a large local supervoid underdense by ≈ 30% out to a radius of ≈ 300 Mpc (Keenan et al. 2013). Though observed in multiple surveys, such a large and deep underdensity is incompatible with ΛCDM at 6.04σ (Haslbauer, Banik, & Kroupa 2020). Those authors showed that such a void could arise naturally in a MOND cosmology supplemented by light sterile neutrinos, as proposed by Angus (2009). This model would have a nearly standard expansion history, primordial light element abundances, and CMB anisotropies. However, it would produce more structure than ΛCDM, in line with earlier analytic expectations (Sanders 1998). This allows for the existence of large 'Hubble bubbles' with enhanced apparent expansion rate, sufficient to explain several key local Universe observables at only 2.53σ tension using a background expansion history fixed to the baseline assumption of Planck Collaboration VI (2020) and without placing us at a special location in the void. MOND with sterile neutrinos therefore appears to account for astronomical observations ranging from the kpc scales of galaxies all the way to the Gpc scale of the local supervoid, without causing any obvious problems in the early Universe. The large scale implications of MOND should therefore be considered in future work, especially in light of upcoming surveys and ongoing tensions explaining current observations within ΛCDM.

Disk galaxy stability

In this contribution, we deal with numerical simulations of disk galaxies, focusing in particular on the Local Group galaxy M33. In principle, this regime allows us to avoid the use of a relativistic MOND theory, whose weak-field limit should give a classical modification of gravity. * Previous analytic and numerical studies indicate that MOND can stabilize a self-gravitating disk without DM (Milgrom 1989;Brada & Milgrom 1999). Essentially, this is because the deep-MOND force ∝ √ M rather than M (Equation 1), limiting the extra gravity created by an overdensity. This effect saturates in the deep-MOND limit (DML), so MOND can provide only a limited amount of additional stability even for a galaxy * A relativistic MOND theory was recently developed in which gravitational waves travel at the speed of light (Skordis & Z lośnik 2019).

with arbitrarily low surface density. This might explain why LSBs often have spiral features (McGaugh et al. 1995), usually taken as a sign that disk self-gravity is important (Lin & Shu 1964;Das et al. 2020).

The global stability of disks thus offers a promising way to break the CDM-modified gravity degeneracy in galaxies. M33 is one of the nearest galaxies with a rather low surface brightness. S19 recently conducted isolated N -body and hydrodynamical (hydro) simulations of M33 in the CDM context. As with many other galaxies, its observed RC required them to use a cored DM profile. They found it rather difficult to explain the weak observed bar together with the two-armed spiral (Corbelli & Walterbos 2007) in the presence of a live CDM halo. Weak bars are also observed in many other galaxies (Cuomo et al. 2019). They are problematic in the CDM picture because a bar in the baryons creates a response in the live DM halo, which would cause a resonant effect (Athanassoula 2002). This was demonstrated explicitly by S19, who compared their live halo simulations to a 'frozen halo' simulation in which the DM particles provide a fixed gravitational potential. † S19 showed that M33 would be more stable in a frozen halo, leading to a weak bar similar to that in the observed M33. Such a model is of course unphysical in the ΛCDM context, but useful to better understand the dynamics.

Interestingly, a frozen halo is somewhat similar to gravity theories like MOND where galaxies lack hypothetical DM. This is because features in the disk very rapidly become irrelevant to the potential as one moves away from its plane. Consequently, perturbations in the disk have very little effect on the Milgromian 'phantom DM' (PDM) halo as its density depends only on the local potential. ‡ Indeed, this is why a bare disk can be stable in MOND. The lack of an actual DM halo automatically removes the possibility of bar-halo angular momentum exchange. However, disk self-gravity is enhanced in MOND compared to the case of a frozen DM halo. The presence of an EFE would slightly decrease this enhancement, but would also decrease the spheroidal potential support from the PDM. In general, it would also break the axisymmetry and up-down symmetry of the problem. Clearly, numerical simulations are needed to properly understand how disk stability is influenced by a fundamental modification to the gravity law.

Several MOND N -body simulations have considered isolated disks, starting with the seminal work of Brada & Milgrom (1999) which used a custom potential solver but restricted particle motion to a plane. Tiret & Combes (2007) conducted a series of 3D simulations, finding rapid growth of the bar strength due to the enhanced role of disk self-gravity. However, the bar then weakened with time as the vertical velocity dispersion increased (see their section 5.1). The basic picture remained the same when those authors included gas as sticky particles in their simulations, though including dissipation in this way weakened the bar somewhat (Tiret & Combes 2008a). It was also shown that it is difficult to form a significant bulge through secular processes in clumpy disks (Combes 2014), potentially explaining the abundance of disks with a low bulge fraction (Kormendy et al. 2010). MOND simulations of disk galaxies have also been used to understand how they interact with each other (Tiret & Combes 2008b;Renaud et al. 2016;Bílek et al. 2018), and with the gravitational field of a galaxy cluster (Candlish et al. 2018). These works leave no doubt that MOND can stabilize a bare self-gravitating disk provided its surface density is low enough that its dynamics are significantly non-Newtonian, i.e. the central surface density

Σ 0 Σ M ≡ a 0 2πG , (3)
corresponding to a vertical Newtonian gravity a 0 . Moreover, the formation of disk galaxies from collapsing gas clouds was recently investigated for the first time in the MOND context (Wittenburg et al. 2020). This naturally led to the formation of exponential disks after 10 Gyr of evolution. Somewhat too compact stellar bulges were also formed, but this problem could perhaps be avoided with a more realistic formation scenario including gas accretion. In particular, it has previously been shown that gas-rich clumpy galaxies in the early universe do not secularly form bulges from the coalescence of clumps in the MOND context (Combes 2014). Gas accretion would be naturally included in cosmological MOND simulations, which are currently under way (N. Wittenburg et al., in preparation).

In this contribution, we follow up the work of S19 in conventional gravity by conducting MOND hydro simulations of a galaxy with initial parameters chosen to match the photometry and RC of M33, with and without the EFE from M31. These global properties can be matched in MOND (Sanders 1996) despite claims by Corbelli & Salucci (2007), in particular in the presence of a small gradient in the stellar mass-to-light ratio (see figure 22 of Famaey & McGaugh 2012). To investigate the dynamical stability of our models, we simulate M33 for at least 6 Gyr using the publicly available code Phantom of RAMSES (por, Lüghausen et al. 2015). por adapts the potential solver of the ramses code (Teyssier 2002) to solve the rather computer-friendly quasilinear formulation of MOND (QUMOND, Milgrom 2010). We analyze our M33 simulations similarly to S19.

In the following, we describe the initial conditions and setup of our isolated simulations (Section 2). We then present our results and analyses in Section 3, where we also compare with observations. We extend our model by including the EFE from M31 in an approximate manner (Section 4). Our conclusions are presented in Section 5. Videos of our isolated simulations with frames every 10 Myr are publicly available, along with the algorithm we use to extract data on gas cells. *

INITIAL CONDITIONS & SIMULATION SETUP

We first conduct two isolated hydro simulations of M33 in MOND, varying the initial gas temperature T to assess the impact of how we model the gas component (Section 2.3). We also run a stellar-only comparison simulation with the same total surface density profile. Our initial setup has the M33 disk spin vector pointing along the positive z-axis of our Cartesian (x, y, z) coordinate system. We frequently make use of cylindrical polar co-ordinates (r, z), where r ≡ x 2 + y 2 .

In MOND, encounters between M33 and other galaxies are less likely to end in a damaging merger due to the absence of dynamical friction from DM halos. A merger is possible for an almost head-on encounter, but geometrically a flyby is more likely. In this case, the other galaxy should still be visible. The nearest large galaxy to M33 is M31, but combining the orbit modeling of Patel et al. (2017) with the latest proper motion measurements of both galaxies indicates that a recent close interaction is very unlikely (figure 4 of van der Marel et al. 2019). There may have been an interaction ≈ 6.5 Gyr ago (Tepper-García et al. 2020), but M33 should still be stable in its current configuration as its dynamical time is much shorter (Section 2.4). We also ignore the EFE at first, but we consider its possible influence in subsequent simulations (Section 4).

Our simulations use QUMOND (Milgrom 2010), which is based on a Lagrangian that obeys the usual symmetries and thus conservation laws. Its numerical implementation requires only the standard linear Poisson solver (Equation 18). The non-linearity of MOND is handled using a non-linear algebraic step, but crucially there is no need for a non-linear grid relaxation stage, a key part of the older aquadratic Lagrangian formulation (AQUAL, Bekenstein & Milgrom 1984). Both AQUAL and QUMOND can be implemented using a publicly available numerical N -body solver called raymond (Candlish et al. 2015), which has clarified that the two formulations give rather similar results (Candlish 2016). Not only do the two approaches agree exactly in spherical symmetry, the differences would be very small even in quite non-spherical situations like a point mass in a dominant external field, as can be shown analytically Sanders (1996) which was based on nearinfrared photometry obtained by Regan & Vogel (1994). The stellar disk scale length agrees with the determination of Kam et al. (2015). Part of this disk is modeled as gas, giving a total gas fraction of 0.3. This leads to the gas having a component distributed similarly to the stars, which we expect because of stellar evolution. The gas is assumed to be isothermal at temperature T = 10 5 K (100 kK) or 25 kK, though T includes the turbulent motions. Radially, each component is modeled as an exponential disk with sech 2 vertical density profile. The aspect ratio is the ratio of their characteristic scales, e.g. the stellar disk density ∝ sech 2 (z/0.178 kpc).

In our hydro simulations, we then adjust the gas disk scale height at each radius using a Newton-Raphson algorithm in order to best achieve equilibrium (Section 2.3). Our stellaronly model includes the gas disk as collisionless particles, so all our simulations have the same total surface density profile.

(section 7.2 of Banik & Zhao 2018a). Therefore, our results can probably be generalized to AQUAL as well.

M33 mass distribution

We treat M33 as a combination of two radially exponential disks with different scale lengths. No bulge is present initially, though one could form if the initial setup is unstable. The more extended disk is purely gas while the less extended disk is mostly stars, with a small proportion of gas to account for gas lost from stars as they evolve. Both disks have a sech 2 vertical profile. Table 1 summarizes the parameters of our M33 model, whose initial surface density profile is shown in Figure 1. All our simulations start with the same total surface density profile and RC, with our hydro models all having the same initial mixture of stars and gas.

Initializing a Milgromian disk

We set up a MOND disk by adapting the Newtonian code Disk Initial Conditions Environment (dice, Perret et al. 2014) -our modified version is publicly available. * dice offers the advantage that the Jeans equations are not solved using the potential, which is difficult to define for an isolated system in MOND. Rather, dice uses only the Newtonian gravity g N , which it calculates using the principle of superposition accelerated by a fast Fourier transform. To get the true gravity g, we use the alge- * https://bitbucket.org/SrikanthTN/bonnpor/src/master/ Separate folders are used for the stellar-only and hydro versions. 3). M33 has a lower surface density almost everywhere, so MOND should greatly impact its dynamical properties. The lower dotdashed magenta line shows the vertical component of the Newtonian-equivalent external field from M31 if we put this at 30 • to the disk plane (Equation 26). Notice that the EFE could have a percent-level effect on the gravitational field in the central regions of M33 (Section 4).

braic MOND approximation:

g ≈ νg N , where (4)

ν = 1 2 + 1 4 + a 0 g N . ( 5
)
We employ the notation v ≡ |v| for any vector v. The 'simple' form of the ν function is used to interpolate between the Newtonian and MOND regimes (Famaey & Binney 2005) since it provides a good fit to a variety of data on galactic and extragalactic dynamics (Gentile et al. 2011;Banik & Zhao 2018a). It is rather similar to the function used by McGaugh et al. (2016) to fit the Spitzer Photometry and Accurate Rotation Curve dataset (SPARC, Lelli et al. 2016a). In the QUMOND approach, ν depends only on g N and is thus readily computable once standard techniques are used to obtain g N .

The algebraic MOND approximation is exactly correct in spherical symmetry and works rather well in axisymmetric problems (Angus et al. 2012;Jones-Smith et al. 2018). It is expected to work particularly well just outside the disk (Banik et al. 2018a). However, it becomes inaccurate within the disk due to the steep vertical gradient in ν caused by that in g N,z . Naively applying Equation 5 would imply a rapid change in g r with z, something that is physically unrealistic as it would cause ∇ × g = 0, allowing energy to be gained by a particle going around a closed loop. * To avoid this, we fix the value of g N,z entering the calculation of ν (Equation 5) to 2 tanh (2) πGΣ if |z| is small enough that the fraction of the local column density Σ (r) at even smaller |z| falls below tanh (2). This is based on the assumption that g N,z = 2πGΣ at the disk 'surface', which is valid for a thin disk. Once ν is calculated in this revised way, we set g = νg N .

To adjust the gravity law, we also need to change how the Toomre Q parameter is calculated -this is used to ensure the velocity dispersion is consistent with local stability (Toomre 1964). Classically, the Toomre stability condition (his equation 65) is that:

Q ≡ σ r Ω r 3.36 GΣ ≥ 1 , (6)
where σ r is the radial velocity dispersion, Ω r is the radial epicyclic frequency, and Σ is the disk surface density.

The QUMOND generalization of the Toomre condition was derived analytically in Banik et al. (2018a). Briefly, the effective value of the gravitational constant G should be enhanced according to:

G → Gν 1 + K 0 2 , K 0 ≡ ∂ ln ν ∂ ln g N . (7)
Their derivation requires g N to include a vertical component of 2πGΣ along with whatever radial component exists at the desired r. Since the Toomre condition is only concerned with local stability, we use the Q lim option in dice to require that Q ≥ 1.25 everywhere, thus leaving a modest safety margin. Results are similar if a slightly lower floor is imposed on Q (Section 4.3).

Figure 2 shows the M33 RC resulting from our modified dice algorithm. We also show the observed RC and the MOND fit based on a much more detailed surface density profile than we assumed (Figure 1). As expected, the MOND fit works rather well despite lacking the flexibility afforded by a DM halo. Since M33 is in the DML almost everywhere, its RC rises along with the enclosed mass (Equation 1) and lacks a region with an approximately Keplerian decline. Due to the disk geometry, our calculated RC nonetheless rises slightly above the asymptotic value of 101 km/s (Equation 2) before gently declining toward it. Despite our rather simplified mass model for M33, its observed RC is rather similar to that of our dice template.

Including gas

We model the gas component of M33 as initially isothermal at T = 10 5 K (≡ 100 kK) or 25 kK. We will see that the cooler model is much more realistic, so we use it for further analyses and simulations. Since we are only interested in the large scale behaviour of M33, we suppress star formation and metallicity-dependent cooling. We do allow gas in our simulations to cool, but only down to its initial temperature. Allowing further cooling would require us to allow star formation, which is beyond the scope of this project. Our rather high T partly compensates for the fact that feedback is not included in our simulations. The use of a constant T simplifies our initial setup because it is already in thermal equilibrium. Even so, ensuring dynamical equilibrium is non-trivial and we only approximately achieve this, as discussed next.

Thickness profile

The inclusion of gas is complicated because there are also stars. We fix the thickness of the stellar component and iteratively adjust that of the gas, following a similar approach to equation 11 of Corbelli (2003).

We start with the Newtonian result of Spitzer (1942) that an isothermal gas slab with sound speed c s has a sech 2 vertical profile with characteristic thickness h g given by:

c s 2 = πGh g Σ g = g N,z h g 2 , ( 8
)
where the gas surface density Σ g gives rise to a vertical Newtonian gravity of g N,z = 2πGΣ g at the disk 'surface', i.e. as z → ∞. We generalize this to MOND by supposing that a sech 2 vertical profile remains a good description if we modify Equation 8 to:

c s 2 = g N,z h g 2 ν g N,r 2 + g N,z 2 . ( 9
)
The issue now arises of what value to use for g N,z in the presence of a stellar disk with surface density Σ * and scale height h * . Since this will be found in a different way to evaluating g N,z at a particular position, we use g N,z to denote the appropriate value of g N,z in Equation 9. Thus, the generalization of Equation 8 to a stellar + gas disk is:

c s 2 ≡ g N,z h g 2 ν g N,r 2 + g N,z 2 . ( 10
)
If h g = h * and the stellar disk is also isothermal, we can treat Equation 8 as applying to the combined stellar + gas disk since there is no practical distinction between the individual components. Thus, it is clear that

g N,z = 2πG (Σ g + Σ * ) , h g = h * . (11)
In general, h g = h * . We make the assumption that each component still follows a sech 2 vertical profile. We then find what fraction of the stellar disk is enclosed at |z| < h g . For the gas disk, this fraction is tanh (1). These fractions are unequal if h g = h * . In this case, we assume that the appropriate generalization of Equation 11 is:

g N,z = 2πG Σ g + Σ * tanh (h g /h * ) tanh (1) . ( 12
)
This approach can easily be generalized to stellar disks with multiple components that have different scale heights and/or vertical profiles different from sech 2 . At large radii, the surface density of an exponential disk becomes very small. Since the gas disk is isothermal, this weakening of the vertical gravity causes the disk to flare (as observed for the Milky Way, Sánchez-Salcedo et al. 2008). As a result, it is no longer reasonable to assume that the vertical gravity comes mostly from disk material at similar r. Instead, it mostly comes from material at much smaller r. For an observer at large r, this material can be approximated as a point mass at the centre of the galaxy. Following our previous approach, the important quantity is the resulting vertical Newtonian gravity at a height of h g . This is approximately given by g N,r (h g /r) for r h g . To prevent our correction term diverging at small r, we put in a softening length of h g . We are then left with the vertical gravity receiving an extra contribution of g N,r at small r, contrary to the idea of applying a 'geometric correction' that only becomes significant many disk scale lengths out. We fix this with an extra factor of tanh (r/h * ), . Initial sech 2 scale height of our simulated M33 gas disk as a function of r, with the adopted temperature shown in the legend. This thickness is much larger than assumed when calculating the RC in dice (Section 1). This weakens the radial gravity at intermediate radii, implying that our hydro simulations are not set up exactly in radial equilibrium (see text).

leading to our final approximation that:

g N,z = 2πG Σ g + Σ * tanh (h g /h * ) tanh (1) + g N,r h g r 2 + h g 2 tanh r h * . ( 13
)
We can now solve Equation 10 by varying h g using a Newton-Raphson algorithm, each time recalculating g N,z as above. Once we have obtained h g at some r, we use it as an initial guess for h g at the next r in our grid, minimizing the number of steps required to reach convergence. In this way, we build up the gas thickness profiles h g (r) shown in Figure 3, which we incorporate into our hydro simulations as described next.

Density and velocity field

We set up the gas density and velocity field in por using our estimated gas disk thickness profile h g (r) from our modified dice algorithm (Section 2.3.1). Practically, this involves adjustments to por, which we now describe. In principle, these adjustments are unrelated to the gravity law -this only enters the algorithm by way of the RC and thickness profile, which are provided externally by dice. Other codes could in principle be used to prepare these input files.

At each r, we assume the gas follows a sech 2 vertical profile with characteristic height h g (r). As a result, the gas density is:

ρ g (r) = Σ g 2 h g sech 2 z h g . ( 14
)
We next assign the gas velocity field. Consider a gas parcel with density ρ and isothermal sound speed c s rotating around a galaxy at cylindrical radius r with speed v a , which in general can differ from the circular velocity v c ≡ √ -rg r due to pressure gradients. The equation of radial hydrostatic equilibrium is:

v c 2 = v a 2 -c s 2 ∂ ln ρ g ∂ ln r . ( 15
)
Because h g depends on r, it is non-trivial to ensure even radial hydrostatic equilibrium of the gas. We approximately achieve this in the disk mid-plane by realizing that ρ g ∝ Σ g /h g (Equation 14), implying that:

∂ ln ρ g ∂ ln r = ∂ ln Σ g ∂ ln r - ∂ ln h g ∂ ln r , z = 0 . ( 16
)
The radial gradient of h g is found by centred differencing. We avoid numerical difficulties at the origin by exploiting the fact that:

∂ ln h g ∂ ln r = r h g ∂h g ∂r .
(17)

The gas velocity field is set up by applying the above pressure correction to the RC calculated by dice based on a constant gas disk thickness of 0.4 kpc (Table 1). Since the actual thickness is generally larger (Figure 3), our algorithm is slightly out of equilibrium even in the radial direction. However, we expect the disk to settle down within a few dynamical times.

Simulation setup

We simulate M33 using por (Lüghausen et al. 2015) in non-cosmological particle-in-cell mode. por solves the field equation of QUMOND, which involves obtaining the Newtonian gravity g N and then solving:

∇ • g = ∇ • (νg N ) . (18)
The generalization to non-zero EFE is given in Equation 28, with the simple form of the interpolating ν function (Equation 5) used in all cases. Our adopted cubic box size has sides of 512 kpc, much larger than M33. This is necessary because por assumes that the gravitational potential Φ satisfies the boundary condition appropriate for a point mass M , namely:

Φ = GM a 0 ln R , ( 19
)
where R is the position relative to the centre of mass in any unit of length. The simulations are advanced for at least 6 Gyr, which represents ≈ 50 revolutions for a typical extent of 2 kpc (Table 1) and rotation velocity of 100 km/s (Equation 2). This should provide ample time for our M33 models to settle into dynamical equilibrium, though results for the first ≈ 1 Gyr might depend on details of our initialization procedure, especially for our hydro simulations (Section 2.3).

To provide adequate spatial resolution, we refine the simulation volume into levelmin = 7 up to levelmax = 13 refinement levels, i.e. there are at least 2 7 pixels across each of the three spatial dimensions. ramses uses adaptive mesh refinement (AMR) to improve the resolution in regions with a high density. We allow up to 13-7 = 6 levels of refinement, giving a maximum spatial resolution of 512/2 13 kpc = 62.5 pc, sufficient to resolve each disk scale length into many pixels (Table 1). Our chosen refinement conditions are mass sph = 10 3 M and m ref ine = 20. Refined regions operate at a reduced timestep as defined by the runtime parameter nsubcycle, which we set to (1, 1, 2, 2). We also set the Poisson parameter epsilon = 10 -4 , defining the convergence condition for the iterative Poisson solver. Note that a standard Poisson solver is sufficient for QUMOND simulations as the non-linearity of MOND is handled in an algebraic step (Equation 18). Further details of ramses can be found in Teyssier (2002), which also gives the default values for parameters that we do not adjust.

Our simulations use 10 6 particles for M33, which we divide amongst its inner and outer exponential disks in proportion to their mass. Since all particles have the same mass within each component, this ensures that masses are also equal between the components. In our stellar-only simulation, we do not implement the procedures described in Section 2.3 since these relate to the gas component. In this case only, we advance our por simulation in N -body mode (the runtime parameter hydro is set to false).

To include gas with a gas fraction of 0.3, we remove all particles in the more extended component and apply a uniform fractional reduction to the masses of particles in the less extended component. We then put this removed mass back in as gas using a modified condinit routine, which sets the properties of each gas cell using the procedures described in Section 2.3. The RC it uses is obtained from dice, which is slightly different to that generated by the mass distribution. However, the difference should be small because it arises only from the initial gas disk being thicker than assumed for the dice RC calculation.

As we are only interested in the large scale behaviour of M33, we suppress star formation and metallicitydependent cooling. We allow the gas to cool, but impose a temperature floor of T 2 star = T 2 ISM , where T 2 ISM is the initial temperature (100 kK or 25 kK). Feedback processes are not included in our simulations, differences between which are summarized in Table 2.

RESULTS OF ISOLATED SIMULATIONS

We extract por data on the simulated particles into plain text format using an algorithm that we make pub- Table 2. Summary of the simulations presented in this contribution, and the section in which each is mostly discussed.

If included, the EFE has a strength of 0.07 a 0 . In our stellaronly simulation, the gas is treated as stellar particles, thus preserving the total surface density profile in all models. We also ran an isolated simulation at 10 kK, but abandoned it as it became very unstable (not shown).

licly available. * For extracting the gas distribution, we use a modified version of the rdramses algorithm † to create a list of positions, velocities, and masses for all gas cells. This lets us analyze particles and gas in the same way, with gas cells treated as particles at the cell centres. We then subtract the combined centre of mass position and velocity in all analyses, correcting for a very small numerical drift over the course of our simulations. This drift is more significant once we include the EFE (Section 4). Apart from images of our simulations, all our quantitative analyses focus only on stars and gas at |z| < 20 kpc to ensure robustness against material ejected to large distance.

Rotation curve

The evolution of our M33 model can be seen in its RC, which we calculate based on the radial gravity g r acting on particles within 1 kpc of the disk plane. The circular velocity is calculated as:

v c = √ -rg r . ( 20
)
The simulated RCs of our isolated hydro models are shown in Figure 4. The mass distribution becomes more centrally concentrated, boosting the RC at low r. The * https://github.com/GFThomas/MOND/tree/master/extract por † http://www.astro.lu.se/∼florent/rdramses.php Note that we obtain reliable results only if rdramses is operated on just one core, though the simulation it analyses need not be.

Our modified version is available at https://seafile.unistra.fr/d/843b0b8ba5a648c2bd05/ decrease at intermediate r is caused by outward spreading of the disk, especially of its gas component. These changes mostly occur within the first Gyr, after which the RC changes very little.

For both gas temperatures, the final simulated RC rises much more steeply than observed in the central regions of M33 (r 4 kpc). This is not much dependent on which observations are used -figure 5 of Koch et al. (2018b) shows that their study gives similar results to the previous studies of Corbelli et al. (2014) and Kam et al. (2017), with the RC similar also to the more recent study of Utomo et al. (2019). This RC mismatch is undoubtedly a serious shortcoming of our isolated models, so we investigate it in more detail to try and understand the reason behind the discrepancy. We will see that it is related to the radial redistribution of material being too efficient, a process which is disrupted by including the EFE such that better agreement is obtained in this case (Section 4).

Cylindrically projected view and the bulge

To better understand the behaviour of our simulated M33, we show it in a cylindrical rz projection where the intensity of each pixel is:

I (r, z) = 1 r 2π 0 ρ (r, φ, z) dφ . (21)
The density ρ actually consists of many discrete particles and gas cells. We find the contribution to each pixel by summing up the total mass in a finite range of r and z regardless of the cylindrical polar angle φ. This is much better than a traditional edge-on (e.g. xz)

projection because flaring at large r can make it very difficult to see the crucially important central regions of the disk. While this may indeed occur in an external galaxy viewed close to edge-on, we can get a better view of our simulations.

Figure 5 shows the stellar particles in our isolated hydro simulations using this cylindrical projection. The initially very thin stellar disk remains rather thin for ≈ 1 Gyr before starting to thicken. By the end of our 100 kK simulation (6 Gyr), though the disk remains rather thin further out, an X-shaped bulge is apparent in the central 2 kpc (top panel). This is related to the rather thick gas disk (Figure 3) being unable to stabilize the much thinner stellar disk. We obtain similar results in our stellar-only simulation, where the gas component is treated as collisionless particles.

The situation is different at a lower temperature of 25 kK, which is closer to the 12 kK adopted in section 3.2 of S19. Our 25 kK model retains a thin stellar disk for at least 7 Gyr (bottom panel of Figure 5), demonstrating the importance of dissipation in the gas component despite its sub-dominant contribution to the central surface density (Figure 1).

Although Corbelli & Walterbos (2007) could obtain acceptable RC fits with a very small bulge having a scale length of only 0.15 kpc, their section 4.1 indicates that they "could not exclude a larger bulge using dynamical arguments". This is also apparent from our results -the RC appears similar in both our 25 kK and 100 kK models (Figure 4), even though only the 100 kK model has a significant bulge (Figure 5). The possibility of a central bulge can also be addressed using photometry. The surface brightness of M33 does in fact show a rise toward its central regions, more so than an inward extrapolation of an exponential law fitted to larger radii (Bothun 1992). This suggests the presence of a 2 kpc bulge (Regan & Vogel 1994). However, other studies found that M33 has only a very small bulge (Kormendy & McClure 1993;Gebhardt et al. 2001). More recent Spitzer photometry at 3.6 µm suggests a bulge scale length of only 0.4 kpc and a bulge fraction of just 4% (section 5.1.2 of Kam et al. 2015). Thus, M33 does not have the sort of bulge evident in Figure 5, indicating that our 100 kK model is unable to correctly reproduce this important aspect of the observations. However, qualitative agreement is gained using a temperature of 25 kK.

The importance of hydro effects in the gas is in line with some previous CDM simulations (Shlosman & Noguchi 1993;Berentzen et al. 1998). Including hydro effects apparently did not much influence the simulations conducted by S19, at least not to the point of obtaining a good match to observations. Nonetheless, the gas stabilizes their ΛCDM simulations somewhat -their section 3.2 indicates that after 2 Gyr, the bar strength A 2 /A 0 (defined in Section 3.3.1) is ≈ 0.2 with hydro treatment of the gas, but their nominal collisionless model reaches A 2 /A 0 ≈ 0.3 by this time (see the 'full mass' curve in their figure 5). Although dissipation in the gas was insufficient on its own to stabilize the M33 disk and prevent it forming a strong bar, hydro effects could well be sufficient when combined with a different gravity law and the lack of angular momentum exchange with a live DM halo. A reduction in A 2 /A 0 upon hydro treatment of the gas is also evident in the MOND model shown in figure 4 of Tiret & Combes (2008a), especially in their Sc galaxy. Thus, it is not surprising that the gas should help to stabilize the disk to some extent, especially in a model where it constitutes a larger fraction of the mass because no CDM component is assumed.

Velocity dispersion

In both models, there is at least a gradual increase in thickness. This is mirrored in the mass-weighted vertical velocity dispersion σ z in the central region of M33 (Figure 6). Since this region is dominated by stars (Figure 1), σ z is not much affected by whether we consider the gas. In our purely stellar simulation, σ z almost doubles (from 20 km/s to 40 km/s) over a 6 Gyr period, though it is rising much more slowly towards the end. The same trend is apparent if gas is included at 100 kK, though the final σ z is reduced slightly to ≈ 35 km/s. A more significant change arises in our model at 25 kK -the stellar σ z now stabilizes at ≈ 30 km/s after only ≈ 2 Gyr, suggesting the model is now stable. This is also approximately when the stellar disk stops thickening (Figure 5). A mild amount of disk heating at this early stage in the simulation could well be due to difficulties setting up the gas and stellar disks in equilibrium with each other (Section 2.3).

The difference in behaviour due to the gas temperature can be understood by relating the above values to the 1D velocity dispersion of the gas. This is σ g = kT / (µm p ) = 21.7 km/s for T = 100 kK, where k is the Boltzmann constant, m p is the mass of a hydrogen atom, and the mean molecular weight for neutral gas with a 25% primordial helium mass fraction is µ = 7/4. We see that σ g is similar to the initial stellar σ z , so the gas is unable to help stabilize the stellar component. However, reducing T to 25 kK reduces σ g to 10.9 km/s, which is much less than the stellar σ z . Therefore, gas in our cooler model can in principle help to stabilize the stellar disk.

To facilitate a more direct comparison with observations, we find the line of sight (LOS) velocity dispersion σ LOS of the stars. We assume that M33 is inclined by i = 50 • with respect to the plane of our sky (figure 1 of Corbelli & Salucci 2007). Although mild warping is evident, this inclination is particularly accurate for the central few kpc of M33. It is also consistent with the i = 52 • ± 2 • obtained in section 4.2.1 of Kam et al. (2015). We project the position and velocity of each stellar particle onto our LOS, which we take to be along the direction (sin i, 0, cos i). This lets us construct an 'image' of M33 based on the sky-projected position of each particle. We divide this image into squares of side 1 kpc, with the centre of M33 coincident with that of the central pixel.

Since each pixel contains a finite number of particles, there is a non-zero correlation between their mean LOS velocity v LOS and that of each particle. The real M33 has many more particles, so this is a numerical artefact which should be corrected. The standard approach is to enhance the measured variance by a factor of N/ (N -1) for a pixel with N particles. For particles with arbitrary masses m i and LOS velocity v LOS,i , the appropriate generalization of this result is:

σ LOS = N i=1 m i N i=1 m i v 2 LOS,i - N i=1 m i v LOS,i 2 N i=1 m i 2 - N i=1 m i 2 . ( 22
)
Once we have obtained v LOS and σ LOS , we perform iterative 5σ outlier rejection with stringent convergence conditions, one of which is that the number of 'accepted' particles must be the same as on the previous iteration and must exceed 150. We do not report σ LOS for pixels with fewer accepted particles as these results are more affected by Poisson noise. We found that the σ-clipping gives very similar results for any threshold number of standard deviations ≥ 2.5, including in the case of not applying outlier rejection at all. We nonetheless apply a conservative 5σ outlier rejection to ensure the robustness of our procedure. We also check in all cases that the v LOS values in the central pixel closely follow a Gaussian of width σ LOS . An example is given in Section 4.2.2.

Figure 7 shows our simulated stellar σ LOS map for our isolated 25 kK model. We include magenta contours bracketing the observational range of (28 -35) km/s, which does not consider the very central regions of M33 (section 3.2 of Corbelli & Walterbos 2007). Assuming only a very small region near the centre is excluded observationally, we expect to get a central σ LOS approximately within this range. However, the central pixel in Figure 7 has a much higher σ LOS = 57 km/s. Thus, the simulated stellar σ LOS is higher than observed. The situation is worse for the 100 kK model (not shown), where σ LOS = 62 km/s. This is related to the simulation developing a significant central bulge (Figure 5) and higher σ z (Figure 6).

Since the stellar and gas RCs trace the same potential, any mismatch between them is partly caused by asymmetric drift, which in turn depends on the size and shape of the stellar velocity dispersion tensor. This technique could yield constraints on σ z , as attempted by Corbelli & Walterbos (2007) in their section 4.1. They did not obtain conclusive results because they found that part of the mismatch must be due to non-circular motions caused by the bar, which we discuss next. Their only conclusion was that σ z 20 km/s, consistent with our simulations (Figure 6).

The stars in our 25 kK simulation form a rather thin disk. The gas has a thicker distribution (Appendix A). This is mainly due to the difficulties we encountered setting up the gas component, which forced us to make it rather hot (Section 2.3). Nonetheless, Koch et al. (2018b) mention in their section 5.1 that there is evidence for a rotationally lagging gas component, which could well be a thick gas disk. Directly detecting this would be difficult due to the inclination of M33, which makes it hard to know how far any detected gas lies from the disk mid-plane. A better comparison with observations would require a simulation at even lower T . In this respect, we attempted to reduce T further, but found that the disk becomes very unstable at 10 kK. This might be addressed by raising levelmax, but ultimately one would need to include the associated star formation and stellar feedback. A higher initial gas fraction might then be required to match the observed value for M33, depending on the resulting star formation rate.

It is thus unclear what effect an even thinner gas disk would have. On the one hand, it would increase the surface density enclosed within the stellar disk, making it less stable. On the other hand, since reducing T from 100 kK to 25 kK helps to stabilize our model, further reductions might reduce σ z further and make the disk even thinner. Any developing instabilities in the stellar disk could be more easily absorbed by the gas disk if it is allowed to cool further, which numerically would involve reducing T 2 star (Section 2.4).

Face-on view and the M33 bar

Figure 8 shows face-on (xy) views of the stellar particles in our 25 kK hydro simulation every 500 Myr. A weak two-armed spiral is apparent after ≈ 3 Gyr, similar to the actual M33 (figure 11 of Corbelli & Walterbos 2007). A fairly strong bar develops quite rapidly, but then weakens as the simulation progresses. This is also evident when viewing the gas face-on (Figure 9).

Results for the 100 kK simulation are shown in Appendix B. There is a very strong bar even after 6 Gyr, to the extent that the whole galaxy is essentially one giant bar. An accurate treatment of the gas component is therefore important to get a simulated M33 resembling the observed one in a MOND context. In what follows, we focus mostly on our 25 kK models.

Bar strength

To quantify the strength of our simulated bar, we follow the approach of S19 and perform a Fourier decomposition of the mass distribution within cylindrical radii of r = (0.5 -3) kpc. We begin by calculating:

A 0 ≡ 1 2 i m i (i labels different particles) , (23)
A m ≡ i m i cos mφ i 2 + i m i sin mφ i 2 , m ≥ 1 .
This quantifies the presence of non-axisymmetry, i.e. dependence on the azimuthal angle φ ≡ sin -1 (y/r). We then find the strength of the m th Fourier mode using:

Strength of Fourier mode m ≡ A m A 0 . ( 24
)
Of particular importance is the m = 2 Fourier mode, whose strength we show for the stellar distribution (Figure 10). We use a magenta line to overplot the nominal model of S19, which they labeled 'full mass' in their figure 5. This simulation contains a live DM halo, but only ran for 2.7 Gyr compared to our 7.4 Gyr. In a ΛCDM context, it is rather unlikely that the bar would get weaker if evolved for longer (e.g. Tiret & Combes 2007). Their figure 7 indicates that a temporary dip in bar strength is possible after a few dynamical times, but the bar then rapidly recovers its strength. A temporary dip is indeed evident in the S19 results reproduced here, but the bar then strengthens after ≈ 1 Gyr, suggesting that there is no scope for subsequent weakening. This led those authors to conclude that their model cannot match the observed rather weak bar of M33, whose strength has been estimated at 0.2 (section 4.3 of Corbelli & Walterbos 2007).

Our simulation initially forms a much stronger bar than the simulations of S19. However, it reaches peak strength after 1 Gyr and starts getting weaker. This is similar to the behaviour evident in figure 4 of Tiret & Combes (2008a), who ran their simulations for (7 -8) Gyr. By the end of our 25 kK simulation, A 2 /A 0 ≈ 0.2, similar to that observed in M33. There are extended periods with a very weak bar (A 2 /A 0 0.1). While this is mitigated to some extent by the generally more important m = 1 mode (Section 3.3.3), it is clear from the face-on view that there are indeed extended periods with very small departures from axisymmetry (Figure 8).

The bar is stronger in our purely stellar simulation (A 2 /A 0 ≈ 0.3, see Figure 10). Including gas at 100 kK leads to similar overall behaviour to the stellar-only model, but with the final A 2 /A 0 ≈ 0.2 (not shown). In both cases, the A 2 /A 0 ratio in our analysis does not indicate a genuinely weak bar comparable to observations. Rather, the face-on view shows that the entire disk is essentially one giant bar, with semi-minor axis comparable to the 3 kpc aperture used to calculate A 2 /A 0 (Appendix B). The very long bar it reveals implies a large corotation radius, since bars should be shorter than corotation (Contopoulos 1980). We discuss this next by finding the bar pattern speed.

Bar pattern speed

Fourier decomposing the stellar distribution allows us to extract the bar pattern speed Ω p . For this purpose, we first define an angle θ m for each Fourier mode, where ). In all hydro simulations, results are very similar when using the total mass distribution instead (not shown).

tan (mθ m ) ≡ i m i sin (mφ i ) i m i cos (mφ i ) . ( 25
We then consider how θ 2 changes with time, adding multiples of π using the method of exact fractions in order to minimize the change between successive snapshots. The bar rotates in a prograde sense at the speed shown in Figure 11. Our hydro models at 25 kK yield Ω p ≈ 30 km/s/kpc. Since galactic bars are typically fast and end close to their corotation radius (e.g. Guo et al. 2019), we use Figure 4 to show a line at this gradient. It intersects the RC at r ≈ 3 kpc regardless of whether we use the simulated or observed RC, suggesting a bar intermediate in length between the stellar and gas disk scale lengths (Table 1). Since our Ω p calculation is based on material at r = (0.5 -3) kpc, the estimated corotation radius of 3 kpc for our 25 kK models should be rather reliable.

This estimate can be compared with the observational result given in section 3.2.7 of Elmegreen et al. (1992) that the corotation radius is 0.4 r 25 , with their table 1 indicating r 25 = 30 for NGC 598 ≡ M33. For our adopted distance of 840 kpc (Kam et al. 2015, and references therein), this corresponds to an observed corotation radius of 2.9 kpc. This is based on the length of one star formation ridge, but the presence of other ridges that extend out to larger r means the result of Elmegreen et al. (1992) is not very secure. Section 5.2 of Corbelli et al. (2019) gives a corotation radius of (4.7 ± 0.3) kpc for the spiral pattern in M33. Since bars generally have a larger Ω p , it could well be that their estimated corotation radius exceeds that of Elmegreen et al. (1992) because the latter result corresponds to the bar. Although a comparison with observations would be quite valuable, corotation radii are rather difficult to measure for a bar as weak as that in M33.

Turning to our less realistic isolated models with only stars or with T = 10 5 K, both yield Ω p ≈ 10 km/s/kpc (Figure 11). A line at this gradient would intersect the flat portion of the RC at r ≈ 10 kpc, suggesting a very long bar (Figure 4). This is actually quite consistent with Appendix B, where we see that the whole galaxy is essentially one giant bar when viewed face-on. However, this is very different to the observed M33 (e.g. figure 11 of Corbelli & Walterbos 2007).

Despite these serious problems, the simulated Ω p in these models is quite consistent with the observational estimate of (10 ± 2) km/s/kpc given in section 4.3 of Corbelli & Walterbos (2007). Those authors assumed the bar ends at its own corotation radius and has a length of (0.6 -1.5) kpc. However, this agreement is not a success of the models -the RC of Koch et al. (2018a) reaches an amplitude of 72 km/s at 1.5 kpc, so a 1.5 kpc long bar ending at corotation has Ω p = v c /r = 48 km/s/kpc. Under this assumption, a shorter bar would have an even higher Ω p (Figure 4). Clearly, the Ω p estimate of Corbelli & Walterbos (2007) is seriously discrepant with subsequent analyses. Its main problem is that it seems to equate Ω p with the local slope of the RC instead of the angular frequency implied by the RC. In other words, in their pattern speed estimate, Corbelli & Walterbos (2007) probably assumed that a bar ending at corotation has Ω p = dv c /dr, instead of the correct Ω p = v c /r. Thus, the agreement between their estimated Ω p and that of our stellar-only or 100 kK models appears to be coincidental.

Other harmonics

We can use Equation 23 to find the strengths of modes with any m ≥ 1. Following Tiret & Combes (2007), we repeat our mode strength calculations for m = 1, 2, 3, 4, and 8. Our results are shown in Figure 12 based on the distribution of stellar particles with r = (0.5 -3) kpc in our 25 kK isolated model. Results are similar if the total mass distribution is used instead (not shown). It is clear that the m = 1 mode is dominant except for brief periods when it is overtaken by the m = 2 mode.

It is difficult to know whether the strong m = 1 mode is a problem for our model, since Elmegreen et al. (1992) stated in their section 4 that background intensity gradients make it difficult to reliably obtain the m = 1 mode strength observationally, causing them to not plot this at all. Some asymmetry between different sides of M33 is required to explain the ≈ 10 km/s differences in the RC inferred from its approaching and receding halves (Kam et al. 2015(Kam et al. , 2017)). In any case, the dominant modes of non-axisymmetry in the stellar distribution are clearly m ≤ 2, with very little power in higher harmonics. Moreover, Figure 9 shows a very symmetric twoarmed spiral in the final snapshot of the gas. In a more advanced simulation, this may well cause a prominent 23) for stellar particles at r = (0.5 -3) kpc in our isolated 25 kK hydro model. Some modes can have strength > 1 if multiple modes are present and their phases (Equation 25) are fortuitously aligned. Results remain rather similar when including the gas (not shown).

spiral arm traced by newly formed stars, perhaps resembling the observed two-armed spiral structure of M33 (e.g. figure 11 of Corbelli & Walterbos 2007). Therefore, our isolated 25 kK model seems to resemble the observed M33 in many respects. It is also possible to avoid a dominant m = 1 mode in one of our models with the EFE (Section 4.4).

The number of spiral arms is related to the wavelength most unstable to amplification by disk self-gravity. In general, we expect this to be shorter in Newtonian models with an appropriate DM halo, leading to a larger expected number of arms (section 3.3 of Banik et al. 2018a). While this is clearly a positive point for our isolated MOND simulation of M33 at 25 kK, the presence of M31 and its possible EFE on M33 motivate us to consider relaxing the assumption of isolation in the next section. This is especially relevant in light of the theoretical requirement for the EFE in MOND (Milgrom 1986), the clear disagreement between data on local wide binaries and MOND expectations if the Galactic EFE were artificially removed (Pittordis & Sutherland 2019), and the recent detection of the EFE through accurate RC measurements of galaxies in different environments (Chae et al. 2020). There is no classical analogue to the EFE, which violates the strong equivalence principle. 2 /d = 0.07 a 0 , which lies well within the DML. This is much less than the central surface gravity of M33, which slightly exceeds a 0 (Section 2.1). Therefore, the a priori expectation would be that M31 does not have a big effect on M33 via the EFE. This is nevertheless worth investigating for the above reasons.

To estimate the possible impact of the EFE on M33, we conduct another hydro simulation with the same disk template (and thus parameters) as our isolated model, but this time including a constant external field of 0.07 a 0 directed at 30 • to the disk plane:

g ext 0.07 a 0 =    cos 30 • 0 sin 30 •    . ( 26
)
This can be thought of as placing M31 at (173.2, 0, 100) kpc relative to M33. We use 30 • as a representative angle -this is the median inclination of a randomly chosen vector to a fixed plane. Other g ext orientations are considered in Section 4.4. Note however that the magnitude and direction of g ext would have been different at earlier times due to orbital motion of M33 around M31. It would be interesting to consider time variation of g ext by directly including M31 in the simulation volume as a collection of static particles, similar to the technique used in Thomas et al. (2017) to minimize the computational cost. However, this goes beyond the scope of the present work, which intends to check if the EFE can influence disk stability. Given the relatively slow orbit of M33 around M31 compared to the rotation period of M33 (Section 2.4), time variation of g ext should have only a small effect on our results. This is discussed further in Section 4.5.1.

To simplify our analysis, we neglect the effect of tides from M31 even though we consider its EFE. This is because the tidal stress on M33 is of order g ext (r d /d), where r d ≈ 2 kpc is the M33 disk scale length (Table 1), and d ≈ 200 kpc is the distance to M31 (Patel et al. 2017). The tidal stress from M31 is thus two orders of magnitude less significant than the EFE. Since even the EFE is rather sub-dominant in the central regions of M33 (Figure 1), neglecting the tidal stress from M31 should be a very good approximation at the present epoch.

Numerical implementation

Including the EFE in our por simulations requires adjusting the PDM calculation and the boundary condition for the potential, which we now describe. The simulation setup is otherwise the same as before (Section 2), apart from some changes to the resolution settings (Section 4.1.2). We also use the same disk template as the RC should not be affected very much in the region of interest (see below). Given our results in Section 3, we use T = 25 kK for all our models with the EFE.

Changes to the gravity solver

The gravity solver is adjusted to impose a constant external gravitational field on the simulated domain. Thus, g ext adds an additional source of gravity to Equation 18, altering ν. In the QUMOND approach, the precise way in which this occurs depends on the Newtonianequivalent external field g N,ext , which is what the external field would have been in Newtonian gravity. Since M31 can be thought of as a distant point mass, we assume g N,ext is parallel to g ext , with their magnitudes related by Equation 5:

g N,ext = g ext 2 g ext + a 0 = 4.6 × 10 -3 a 0 . (27)
To retain our previous notation that g N refers to the Newtonian gravity of M33 alone, we generalize Equation 18to:

∇ • g = ∇ • g N,tot ν g N,tot , (28) g N,tot ≡ g N + g N,ext . (29)
In a more detailed model, g N and g N,ext would be directed towards M31, making their magnitude and direction time-dependent. The EFE generally has the effect of reducing ν, making the system more Newtonian. However, partial cancellation between internal and external fields is also possible, which would raise ν for a weak EFE.

In addition to changing ν inside the simulation box, the EFE also changes the boundary condition -Equation 19 is only valid for a point mass if the EFE can be neglected. Φ has a rather complicated behaviour if internal and external fields are comparable (e.g. Thomas et al. 2018;Banik & Kroupa 2019). However, it is possible to obtain Φ analytically if the boundary is external field dominated, i.e. if g N,ext

g N (section 3 of Banik & Zhao 2018c). Treating the matter distribution as a point mass M , they obtained the solution:

Φ = - GM ν ext R 1 + K 0 2 sin 2 θ , ( 30
)
where θ is the angle between g N,ext and the position R, while K 0 was defined earlier (Equation 7). This solution was confirmed numerically by direct integration of Equation 18 (section 2.2 of Banik & Zhao 2018b). The external field dominates for distances R R EF E , where the transition radius in the DML is:

R EF E = GM a 0 g ext . ( 31
)
Since the M33 mass M = 6.5×10 9 M (Table 1), the external field from M31 dominates beyond R EF E = 41 kpc, much larger than M33 (Figure 1). We therefore increased the box size to 1024 kpc, ensuring a minimum distance of 512 kpc R EF E from M33 to the edge. At this distance, it would be quite accurate to treat M33 as a point mass. Therefore, we replace Equation 19 with Equation 30 for our simulations with the EFE.

Changes to resolution settings

Our preceding discussion shows that including a weak external field is quite computationally expensive because the box size must be large enough for the boundary to be external field dominated. Due to the high computational cost and memory requirement of our simulation, we reduce the maximum refinement level to levelmax = 12. The most resolved regions thus have a resolution of 1024/2 12 kpc = 250 pc, which is still much smaller than the scale length of the M33 disk (Table 1). We maintain the previous setting that the minimum refinement level is levelmin = 7, so even the most poorly resolved regions have a cell size of 8 kpc.

To check how these changes affect our results, we then run a higher resolution simulation for 6.1 Gyr with the EFE (Section 4.2.6). The simulation was not advanced further due to significant numerical drift of M33 in roughly the direction opposite to g ext , with the barycentre position changing approximately quadratically over time at an implied acceleration of 3.4 × 10 -3 a 0 . The results of both simulations are sufficiently similar that it is clear the lower resolution simulation has reached numerical convergence. Our higher resolution simulation with the EFE has exactly the same resolution settings and box size as our isolated 25 kK simulation and also uses the same disk template, so any differences can be attributed solely to the EFE. However, to discuss the evolution of M33 over a longer timespan and to minimize the role of edge effects, we focus mainly on our lower resolution EFE simulation, which we advance for 9.9 Gyr. The numerical drift of M33 is much smaller in this case, with an implied acceleration of just 3.9 × 10 -4 a 0 . This is probably due to the larger box size making Equation 30 more accurate on the boundary. *

Results

Rotation curve and surface density

One of the main problems with our isolated hydro simulations is that the RC rises too steeply at r 4 kpc (Figure 4). This problem is largely resolved once we include the EFE (Figure 13). In this model, the inner RC is still evolving somewhat by the end of the simulation, which further reduces the already small difference with the observed RC. Clearly, even a sub-dominant EFE can have an important effect on the stability and secular evolution of galactic disks in MOND. We discuss this further in Section 4.5.

Since the EFE is not very significant at r 2 kpc (Figure 1), changes in the RC arise mainly from dif- * The numerical issue of barycentre drift can be explored with fewer computations using the ring library procedure described in Banik & Zhao (2018b) -a point mass embedded in a uniform gext experiences a non-zero acceleration despite the use of a spherical domain. This can be reduced to a very small level by improving the resolution and enlarging the simulated region. ferences in the matter distribution. We therefore investigate how the surface density profile of the stellar component differs between our hydro models with and without the EFE (Figure 14). In the isolated models, the central Σ * is lower when T is reduced from 100 kK to 25 kK. This is expected given the cooler model prevents the formation of a substantial central bulge (Figure 5). However, the difference is not enough to reconcile the RC with observations (Figure 4). Including the EFE further suppresses the radial inflow of stars, causing Σ * to rise less steeply towards lower r. This is not much dependent on the choice of direction for the EFE, which we vary in Section 4.4.

Our simulation with the 'fiducial' g ext (Equation 26) still displays a small RC discrepancy at r 4 kpc. In addition to percent-level changes in the M33 distance, this could in principle be addressed by changing the initial configuration. Our adopted disk scale lengths (Table 1) are based on M33 today, but its surface density profile was probably different several Gyr ago. To get a shallower RC, the disk should have been more extended. For the same disk mass, this would reduce Σ and thus the disk self-gravity, putting it deeper into the MOND regime and making the disk more stable (Section 1.2). A time-varying external field also needs to be considered in further works.

Moreover, the RC itself could have small observational issues in the central region. The empirically determined RC of a galaxy is mainly sensitive to the radial velocities along its major axis. Even with rather weak spiral arms, it is quite possible for the inferred RC to differ from the true one at the 10 km/s level (figures 3 and 4 of Martinez-Medina et al. 2020). If the galaxy has a 180 • rotational symmetry, this would not be recognisable simply by comparing its approaching and receding halves, though doing so for M33 does reveal differences of this The initial distribution (Figure 1) is shown with a featureless pink line.

order (figure 10 of Kam et al. 2015). Other issues include a slight r-dependent inclination, i.e. a warp. This would still leave a mild disagreement at r ≈ 5 kpc, which could be due to a warp. The tilted ring fit shown in figure 2 of Corbelli & Salucci (2007) shows that at r = 5 kpc, the inclination could be 56 • even though it is only 50 • at the centre. We postulate that a warp might partially mask a rapid flattening of the RC at r ≈ 3 kpc. At larger radii, figure 16 of Kam et al. (2017) shows reasonably good agreement with our simulated RC, which is essentially flat at v f = 101 km/s (Equation 2). This prediction also applies when including the EFE, since g ext is sub-dominant at distances 40 kpc (Equation 31). Thus, the RC at (10 -25) kpc is not such a strong test of our model -most of the mass lies closer in. As a result, our isolated models also yield a similar RC in this distance range despite behaving rather differently at low r (Figure 4). In a different galaxy, the EFE could be more significant in the observationally accessible part of the RC (Section 1.1, see also Chae et al. 2020).

Bulge and velocity dispersion

A bulge is readily apparent after just 3 Gyr in our isolated 100 kK simulation, but the formation of a bulge is suppressed at 25 kK (Figure 5). As might be expected from Figure 14, a bulge also does not arise in our models with the EFE. This is evident in Figure 15, which shows the stellar component in a cylindrical rz projection (Equation 21).

The stellar σ z is similar between our 25 kK models with and without the EFE (Figure 6), but the slightly lower σ z when including the EFE suggests that M33 is more stable in this case. Although σ z is not directly observable, σ LOS is. Observationally, σ LOS = (28 -35) km/s (section 3.2 of Corbelli & Walterbos 2007). In our isolated hydro simulation, the central pixel has σ LOS = 57 km/s, well above this range (Figure 7). When we include the EFE, this decreases to 48 km/s (Figure 16). While still rather high, the agreement is much better. The exact choice of snapshot and pixel size would also have some effect, with the observing direction at fixed i contributing an uncertainty of ≈ 1 km/s.

In addition to σ LOS , observations can also tell us the distribution of v LOS , which we show in Figure 17 for the central pixel. A Gaussian of width σ LOS provides a very good description of the results. This also demonstrates that a mass-weighted variance calculated using standard techniques is sufficient for the central pixel.

The discrepancy in σ LOS could again probably be addressed by a small adjustment of the initial conditions or a time-varying EFE, though this still needs to be demonstrated. In addition, including stellar feedback processes would help to limit the central accumulation of material, which is probably also responsible for the mild tension with the RC (Figure 13). As in the RC case, the comparison with observations could be problematic. One important reason is that our calculated σ LOS is a mass-weighted velocity dispersion (Equation 22), whereas the observational estimate is based on the broadening of spectral lines. Since stars do not have an exactly linear mass-luminosity relation, this can create a slight mismatch between mass-and luminosity-weighted velocity dispersions (Aniyan et al. 2016). This issue was also discussed in section 4.3 of S19. * Since the luminosity of a star generally rises faster than its mass, the bias is likely to be in the right direction, especially in a galaxy that is still forming stars (Verley et al. 2009). This is different to our simulations where star formation is disabled, so they are not directly comparable to observations even if these give a true mass-weighted σ LOS for the stars.

Another interesting feature of our σ LOS map is the slight left-right asymmetry in our model with the EFE (Figure 16), something not evident in our analogous isolated model (Figure 7). This is because the underlying gravitational physics is symmetric with respect to ±x without any EFE, but this symmetry is broken by an external field with non-zero component towards +x (Equation 26). The potential of a point mass is symmetric with respect to g ext when g ext dominates (Equation 22). To ensure the robustness of our results, we use 5σ outlier rejection, but this has almost no effect on σ LOS .

external field is weak. This effect can induce asymmetry in tidal streams of satellite galaxies (Thomas et al. 2018), and in galaxies orbiting within a galaxy cluster (Wu et al. 2010(Wu et al. , 2017)). At large distances, g ext becomes dominant and this asymmetry vanishes. The asymmetry is discussed further in Section 4.4.

Face-on view and the bar

The face-on appearance of M33 remains rather circular once the EFE is included (Figure 18). A weak bar is evident, whose length is related to Ω p . This is not much affected by the EFE (Figure 11). Regardless of whether we use the simulated or observed RC, Ω p = 30 km/s/kpc corresponds to a corotation radius of ≈ 3 kpc (Figure 13). As argued in Section 3.3.2, this is quite reasonable.

In both our isolated and non-isolated models, the bar strength is similar to the observed value of 0.2 by the end of the simulation. The main difference is that the typical strength is larger when the EFE is included (Figure 19). This is more in line with the observed population of spiral galaxies, which often have strong bars (Laurikainen & Salo 2002;Garcia-Gómez et al. 2017). However, a meaningful comparison would need to restrict the observational sample to non-interacting galaxies with a similar surface density to M33, since the stability properties of Milgromian disks depend on their central surface density relative to a particular threshold (Equation 3).

Non-circular motion

Our isolated model of M33 appears mildly non-circular when viewed face-on, especially if we consider just the stars (Figure 8). Non-circular motions are necessary to sustain a non-circular shape. We quantify this using the in-plane tangential velocity

v t ≡ xv y -yv x r . ( 32
)
Most of this is just circular motion, so we first determine the mean v t for particles in different radial bins. This yields a list of mass-weighted mean radii r and tangential velocity v t . We interpolate this to obtain the expected v t at any r 25 kpc. We then subtract v t to obtain the excess tangential velocity:

∆v t ≡ v t -v t . ( 33
)
We find ∆v t for different parts of M33 using the same binning procedure as described in Section 3.2.1. To better illustrate the dynamics of M33, we now use a face-on view. The resulting map of ∆v t is shown in Figure 20. We previously estimated that the corotation radius is 3 kpc in our 25 kK models regardless of the EFE (Section 4.2.3). Thus, non-circular motions should mostly be restricted to the central 3 kpc, which is approximately the case. The non-circular motions as quantified by ∆v t are much smaller when the EFE is included (bottom panel).

The EFE is therefore able to suppress the bar instability somewhat. Another difference is that the isolated model is close to symmetric with respect to φ → φ + π, indicating a dominant m = 2 mode. The EFE breaks the symmetry between ±x but preserves that between ±y, at least for r 6 kpc. This is no doubt a consequence of the particular direction we adopted for g ext (Equation 26). Indeed, the pattern of ∆v t appears to rotate continuously with time in the isolated model, but undergoes rather little variation once the EFE is included.

Disk precession and warping

The EFE breaks the axisymmetry in our isolated models of M33. One consequence is that the whole disk precesses by a small amount over the course of the simulation. We quantify this based on the total angular momentum of all material at |z| < 20 kpc. The direction of this vector is the disk spin axis, which we denote h. The time evolution of h is shown in Figure 21, which reveals a very slow precession by just over 1 • /Gyr.

Since the external field picks out only one preferred direction (Equation 26), we expect h•g ext to remain constant. Indeed, the angle between these vectors (which is initially 60 • ) changes by < 1 • over 10 Gyr, much less than the change in h. Thus, we can consider h as precessing around g ext . Since g ext is within the xz plane (Equation 26), this implies h mostly precesses within the yz plane. In our simulation, the precession is towards -y. Our isolated model of M33 exhibits only a very small amount of precession (< 0.01 • ), which we attribute to numerical noise. Thus, the more significant precession evident in Figure 21 is caused by the EFE, and remains if the simulation is run at higher resolution (Section 4.2.6).

We can understand this by considering Equation 30 and approximating the disk as a point mass. Since K 0 < 0, the potential is deepest along the g ext axis. As a result, stars at x > 0 would feel an extra force towards +z, while stars at x < 0 would feel a force towards -z. Stars at ±y in general also feel an extra force out of the disk plane (Section 4.4), but this would be in the same direction due to the EFE preserving reflection symmetry with respect to the xz plane containing g ext . Thus, the net torque would come only from the asymmetry between ±x (there is none between ±y). Our preceding argument suggests that the torque would cause h to precess towards -y, which is indeed what happens in our simulation. This is a manifestation of the 'external field torque', which arises because the force between two masses is in general misaligned with their separation once the EFE is considered (e.g. Equation 30is not spherically symmetric).

We can obtain a rough estimate of the precession rate using appendix A of Banik et al. (2018b), which provides a general equation for the tangential force on a test particle due to a point mass -this is non-zero only if we consider the EFE. Their result is based on direct numerical integration of Equation 28for a point mass in an external field, with all accelerations in the DML. *

The results are reliably known only for the case of a point mass and a test particle, but this will be sufficient for our purposes. In this case, the precession rate for a ring of material at some radius r is approximately:

Ω = f geom GM a 0 r sin θ cos θ 2v c r ext 3 , ( 34
)
where θ is the minimum angle between g ext and a vector within the ring plane, while f geom = 1/2 is a geometric factor that accounts for azimuthal averaging, i.e. the fact that the torque on a particle at (x = r, y = z = 0) is not representative over all φ. Note that only precession towards -y is considered, because other components cancel out by symmetry. Equation 34 suggests a precession rate of 1.66 • /Gyr for a ring at r = 4 kpc with v c = 100 km/s. This is a rather rough estimate because we have neglected disk self-gravity and the choice of r is arbitrary (though it should be similar to the disk scale length). Nonetheless, the precession rate of 1.1 • /Gyr in our por simulation is not too far off, and is in the predicted direction (Figure 21). The precession of the M33 disk would have a small impact on our previous results, since ideally we should rotate M33 into a reference frame aligned with h before doing further analyses. Fortunately, the effect is rather small because cos 10 • = 0.985, so even after 10 Gyr, derived quantities like the RC would be only slightly affected by the precession. An interesting situation where disk precession does matter is the bifurcation apparent at large r in our cylindrical rz projection (Figure 15). * Banik & Kroupa (2019) provides a visualization of the point mass gravitational field when gext = 1.78 a 0 , the value appropriate to the Solar neighbourhood of the Milky Way (McMillan 2017). This is because the mean height of the disk is z ∝ sin φ, where φ is the azimuthal angle measured from the line of nodes with respect to the initial disk plane (in this case, the x-axis). The intensity of the rz projection is highest when z is stationary with respect to φ, which occurs when φ = ±π/2 and |z| is maximal. As a result, sharp ridges appear in the rz projection, with both ridges having the same angle to the r-axis (Figure 15). However, as the main purpose of this work is to consider the stability of M33's central regions, this does not affect our results very much. shows the total change in h from its initial orientation (along +z), while the dotted green line shows how much of this is a nutation rather than precession around gext. The behaviour of h is very similar in our higher resolution simulation with the EFE (Section 4.2.6), so we show only the total precession angle in this case (solid blue). The precession is almost entirely towards -y in both cases (solid black and dot-dashed magenta curves coincide). The direction and approximate magnitude of the precession are explained in the text using analytic arguments.

Since different parts of the disk are subject to different amounts of internal gravity, we expect that not all parts precess at exactly the same rate, which could cause the disk to warp. This possibility was explored previously by Brada & Milgrom (2000) using integration of test particle orbits in a potential generated by an exponential disk. * Their work neglected self-gravity of the disk, so was unable to explore how its stability might be influenced by the EFE. Since then, no other works have used self-consistent simulations to explore how a MONDian thin disk galaxy would be influenced by a weak external field. The effect on elliptical galaxies was explored in Wu et al. (2017), while Candlish et al. (2018) investigated the impact of using g ext ≈ a 0 , thereby addressing how disk galaxies might evolve in a cluster environment. Such a strong EFE is seriously damaging to the disk, as might be expected -MOND effects are greatly suppressed in these circumstances, essentially reducing the problem to a purely Newtonian disk, which is very unstable (Hohl 1971). However, our results suggest that a weak external field can actually make the disk more * They put gext at 45 • to the disk rather than our adopted 30 • . 26).

stable, reducing the concentration of material at low r (Figure 14).

To explore how the EFE might achieve this, we use the same binning and outlier rejection procedure described in Section 3.2.1 to show z, the mass-weighted mean z of the stellar particles in each pixel. The effects of disk precession are first removed by an appropriate rotation to align the co-ordinate system with h, allowing us to see what M33 would look like face-on if we had reliable depth perception. Its disk is indeed mildly warped (Figure 22). Since the only two vectors in the problem are h and g ext , we expect the xz plane to be a plane of symmetry (Equation 26). This is why the results are nearly symmetric with respect to ±y. However, they are not symmetric with respect to ±z because this symmetry is broken by the z-component of g ext . The resulting asymmetry in the potential is discussed further in Section 4.4.

The mild warping evident in Figure 22 is not sufficient to explain the ≈ 5 • warp inferred by Corbelli & Salucci (2007) -a change in z of 0.1 kpc between r = 0 and 3 kpc implies a slope of only 2 • . The more significant warp identified by Corbelli & Salucci (2007) could be due to recent gas accretion or because of a stronger historical EFE on M33 arising from a smaller separation with M31. Tidal effects might also have been important -tides naturally pull down on one side of M33 while pulling the other side up, which is difficult to do with a constant EFE. However, if tides were able to significantly affect the M33 disk at r = 3 kpc, then they would have much more significant effects further out, contradicting the regular appearance of the M33 disk. 26. The latter has the same resolution settings and box size as our isolated simulations, but was not advanced further due to significant numerical drift of M33. The slightly smaller disk precession in our higher resolution simulation (Figure 21) causes the disk to appear thinner at large r.

Thus, gas accretion or a stronger EFE in the past appear to be more promising explanations for its warping. * This is related to the M31-M33 orbit, which could be constrained with better proper motion data -especially for M31. In addition, it would also be important to consider the Milky Way, which in a MOND context must have previously interacted with M31 (Zhao et al. 2013). Their satellite galaxy planes (Ibata et al. 2013;Pawlowski & Kroupa 2013, 2020;Sohn et al. 2020) might have condensed out of tidal debris expelled during this interaction (Banik et al. 2018b;Bílek et al. 2018).

Numerical convergence

Our simulation with the EFE uses a lower resolution compared to our isolated simulations, which may affect our results (Section 4.1.2). To check if they are numerically converged, we run a higher resolution simulation with a box size of only 512 kpc and allow up to levelmax = 13 levels of refinement. The box size and resolution settings are thus the same as in our isolated simulations (Section 2.4). We only advance our highresolution hydro simulation with the EFE for 6.1 Gyr due to significant numerical drift of M33, which is much * It is of course not possible to disentangle tides from the EFE when following the orbit of a single satellite, but this can be done using numerical experiments, or through comparison with other systems.

less pronounced in our lower resolution model with a larger box size (Section 4.1.2).

We checked that the increased resolution has only a small effect on the RC at 6.1 Gyr. At r 3 kpc, both models including the EFE behave rather similarly, and have a lower v c than our isolated model. Another important consideration is whether a central bulge develops. We test this by using Figure 23 to show rz projections of our M33 simulations with the EFE. The results look rather similar, with no sign of a central bulge in either case. This is also true in our isolated model at the same T = 25 kK (Figure 5).

We therefore conclude that the lack of development of a central bulge in our hydro models at 25 kK is a numerically converged result. If anything, the disk appears slightly thinner in our higher resolution model with the EFE (bottom panel of Figure 23). Thus, neither the resolution nor the EFE affect our ability to stabilize a thin bulge-free Milgromian disk initialized to the M33 surface density profile. We therefore use the lower resolution settings when varying the direction of g ext (Section 4.4).

Reducing the Toomre parameter

Our simulations use a disk template where the Toomre parameter Q (generalized to QUMOND in Banik et al. 2018a) has a lower limit of Q lim = 1.25 (Section 2.2). Using Q lim > 1 ensures some margin of safety, but the exact choice is arbitrary. To explore the impact of Q lim, we run two more versions of our hydro simulation with the EFE at 30 • to the disk (Equation 26). In these simulations, Q lim is reduced to 1.1 and then 1.

The effect of varying Q lim is apparent in Figure 24, which shows σ z for stars with r = (0.5 -3) kpc. Reducing Q lim from 1.25 to 1.1 reduces σ z somewhat, as might be expected. In the model with Q lim = 1, σ z is never constant for an extended period, though it rises more slowly when t ≈ (2 -4) Gyr. The model appears to be marginally stable, and starts rapidly heating when t ≈ 5 Gyr. This suggests that our analytic Toomre criterion (Equation 7) is actually a rather good estimate of what initial conditions would be stable.

These results on σ z are mirrored in σ LOS , which for the central kpc square is 48 km/s for our fiducial model with Q lim = 1.25 (Figure 16). This drops to 47 km/s when Q lim = 1.1, but rises to 53 km/s when Q lim is reduced to 1. The m = 2 mode strength also drops somewhat when Q lim is reduced from 1.25 to 1.1, but then becomes higher if Q lim is reduced further to 1. Therefore, changes to Q lim do not have a significant impact on our results for a reasonable choice of Q lim, i.e. for a value slightly above 1 to ensure local stability.

Different external field orientations

To gain a deeper insight into how the EFE affects our M33 model, we run two more simulations where we put g ext in two special orientations. We use the same simulation duration, numerical settings, and external field 26) in which the floor on the Toomre Q parameter is varied to Q lim = 1 (solid red), Q lim = 1.1 (dot-dashed blue), or has the fiducial value of Q lim = 1.25 (black) used in all other simulations. The simulation at Q lim = 1 appears to become unstable after ≈ 5 Gyr, while the other cases remain stable. Results are very similar if the full mass distribution is used instead of just the stars (not shown).

strength of 0.07 a 0 , but set it to point either within the disk plane (along +x) or along the disk spin axis (along +z). We call these the disk-aligned and axially aligned models, respectively. The barycentre has a numerical acceleration of 2.3 × 10 -3 a 0 in the axially aligned case and 3.9 × 10 -4 a 0 in the disk-aligned case, with the drift being very nearly in the opposite direction to g ext . The disk-aligned model has a numerical drift similar in magnitude to the model with intermediate g ext discussed in the preceding section, but the numerical drift is greater for the axially aligned case. The overall appearance of the disk and its RC are quite similar in all our models with the EFE. They all have a stellar σ z that is ≈ (3 -5) km/s lower than in the isolated case, signifying somewhat more stability. As a result, the central σ LOS differs little between our simulations with the EFE -the axially aligned case yields 49 km/s while the disk-aligned case gives 47 km/s, bracketing the 48 km/s result for the intermediate orientation.

Some differences between these simulations are revealed upon Fourier analysis. The relatively strong m = 1 mode evident in Figure 12 also arises in our models with the EFE -but not for the axially aligned case. Figure 25 shows the corresponding results for this model. Depending on how well A 1 /A 0 is constrained observationally, this could mean the axially aligned g ext case is more realistic. Due to the different symmetry properties of our models with differently aligned g ext , the behaviour of the warp differs significantly, as discussed next. Warping of the disk is negligible ( 0.03 kpc) in the disk-aligned case, which is expected as the problem is symmetric with respect to ±z. Figure 26 shows the disk warp for the axially aligned case. Since the problem is axisymmetric, it is not surprising that the warp also retains approximate axisymmetry. Relative to the centre, the outer parts of the disk become warped in the direction opposite to g ext . The potential is symmetric with respect to ±g ext in both the analytically tractable isolated and g ext -dominated cases (Equations 19 and 30,respectively). This is not true for the intermediate case, as pointed out in earlier works (e.g. Thomas et al. 2018).

To gain insight into the potential, we make use of the DML numerical force library described toward the end of section 2 in Banik & Zhao (2018b). This provides calculations of the gravitational field from a point mass in the DML, with the EFE rigorously accounted for by direct numerical integration of Equation 28. The force evaluations relevant to our discussion here are for points at θ = π/2, since all points in the xy plane are at right angles to g ext as perceived from the centre of M33. Their numerical calculations show that at these points, the force is mostly toward the centre, but there is a small additional component along -g ext . The strength of this symmetry-breaking component is shown in Figure 27 . The gravity in the direction parallel to gext experienced by a test particle due to a point mass located in the orthogonal direction. This is a measure of asymmetry in the gravitational field between ±gext, since results would be zero in the symmetric case -which arises in the EFE-dominated limit (Equation 30). The results shown in solid red are from numerical force calculations in the DML of a point mass embedded in a uniform external field (section 2 of Banik & Zhao 2018b). The units are such that G = M = a 0 = gext = 1, so R EF E = 1 (Equation 31). Results shown here can be scaled to other gext if the problem remains in the DML (Milgrom 2009). The dot-dashed blue line shows the semi-analytic fit given in appendix A of Banik et al. (2018b), with the -sign changed to + in their equation A2 to fix an error in the original. The gravity in the radial direction is not shown here, but is the dominant component at all radii.

The DML force library in Banik & Zhao (2018b) was later fit using the fitting function given in appendix A of Banik et al. (2018b), bearing in mind the gravitational field at all locations. We noticed that the -sign in their equation A2 must be replaced by + to get the correct results. With this correction and restricting to the case of a test particle at θ = π/2 relative to a point mass M embedded in an external field along +z, we expect a vertical force of:

g z g ext = - 2 r 2 5 (1 + r 2 ) (1 + r 3 ) , ( 35
)
where r ≡ r/R EF E is the radius scaled to R EF E (Equation 31), the only physical scale in this DML problem.

Using M = 6.5 × 10 9 M (Table 1) and g ext = 0.07 a 0 as before, we get that g z = -0.0017 a 0 at r = 10 kpc.

We can estimate the height of the induced warp by approximating that the EFE-induced vertical gravity from M33 must be balanced by a geometric term similar to that in Section 2.3.1. Treating the disk as having an extent 10 kpc, the gravity it exerts directly towards itself is very nearly v f 2 /r, with v f given by Equation 2. We expect this to be approximately valid if g ext is sub-dominant to the disk gravity ( r 1), which is true for r R EF E = 41 kpc (Equation 31). With this approximation, the vertical component of the disk gravity near the z = 0 plane is -v f 2 /r (z/r) even without the EFE. We assume that in equilibrium, this approximately balances the additional g z induced by the EFE (Equation 35). As a result, the equilibrium height of the warp can be estimated as:

z eq ≈ r 2 g z v f 2 < 0 . ( 36
)
For the case of M33, we get that z eq = -0.061 kpc. This is fairly close to the actual warp of ≈ -0.15 kpc evident in Figure 26, comparing z at the disk centre with its value at r = 10 kpc. Our semi-analytic approach does not capture the full complexity of an extended disk. If we focus on a point P within the disk mid-plane (not necessarily z = 0) at r = 10 kpc, regions close to P do not exert a net gravitational force along z. The regions at lower r do, but all the material is not concentrated at r = 0. Moreover, even the inner regions of the disk become warped to some extent, making the z co-ordinate of P closer to the barycentric z of the whole system. This reduces the geometric (z/r) factor assumed in Equation 36, so the outer regions of the disk must warp more to compensate. As a result, we expect the true warping of the disk to exceed our semi-analytic estimate, which is indeed the case. Nonetheless, it is possible to determine the approximate extent of disk warping without N -body or hydro simulations.

Although the EFE is a natural part of MOND (Section 1.1) and follows inevitably from its governing equations (Milgrom 1986), its impact on our simulations is rather surprising given that it is sub-dominant at r 41 kpc. In QUMOND, any possible effect on the dynamics must come through the value of ν, which depends solely on the Newtonian gravity (Equation 28). Since g N ∝ ∼ 1/r 2 over the range r = (4 -40) kpc, we see that g N,ext should be sub-dominant by a factor of 100 at 4 kpc. In particular, a point mass of M = 6.5 × 10 9 M creates a Newtonian gravity of g N = 0.472 a 0 = 103 g N,ext at this distance. A similar level of isolation is evident if we consider the z-component of the internal and external Newtonian gravity for our intermediately aligned g ext : at 4 kpc, the former dominates by a factor of 101 if we assume g N,z = 2πGΣ (Figure 1). Thus, the EFE could in principle have a significant impact on the secular evolution of M33 after 100 dynamical times, i.e. after about 100 (r d /v f ) ≈ 2 Gyr, well within the duration of our simulations.

The EFE alters the resonant structure of the disk by breaking axisymmetry. This likely inhibits the usual process where a whole ring of material moves inward and ends up feeling even more radial gravity, with build-up of random motions eventually halting the inward migration. Our simulations are initially stable according to the QUMOND Toomre condition (Section 2.2), but radial migration always occurs to a limited extent. In addition to breaking the axisymmetry of the disk, the EFE generally also breaks its up-down symmetry and causes the disk to warp (Figure 22). This makes it even more difficult to form resonances as different rings of material are in different planes. Restoring axisymmetry or up-down symmetry with a special alignment of g ext has little effect on our results (Section 4.4), suggesting that the more efficient radial migration in our isolated simulation is related to it having both these symmetries. The EFE necessarily breaks at least one of these symmetries. Therefore, the EFE has very important implications for the development of instabilities, especially the bar instability (Figure 20). This in turn affects the efficient redistribution of disk material, which is apparent in Figure 14.

While our model may provide a viable scenario for the weak bar and bulge in M33, bars are often much stronger (Menéndez-Delmestre et al. 2007). Using isophotal ellipticity as a proxy for bar strength, they identified that only ≈ 40% of galaxies have a very weak bar that might be undetected in their analysis. Thus, it is necessary to get a mixture of both strong and weak bars. Our results suggest that this could partly arise from differences in the EFE.

Some evidence for the EFE has already been found in that galaxies with declining RCs in their outskirts generally have an identifiable perturber (Haghi et al. 2016;Chae et al. 2020). Our results show that even a weak perturber would also affect the overall level of dynamical stability, an issue which has not previously been discussed. Our work is the first to conduct a selfconsistent hydro simulation of how a weak external field affects the secular evolution of a Milgromian disk. The implications for the bar strength are somewhat ambiguous and degenerate with other properties like the gas temperature (Section 3). A much more powerful test could come from the warp induced by the EFE, as discussed next.

The external field warp

Our results on disk warping in Section 4.4 highlight that when the internal and external gravitational fields are comparable, the gravity due to a mass distribution becomes asymmetric with respect to ±g ext . Though this was known in prior work and is evident from the point mass potential (Figure 27), ours is the first to demonstrate the impact on a self-consistent numerical simulation of a disk galaxy. The asymmetry also has implications for tidal streams, as discussed further in Thomas et al. (2018) using the example of Palomar 5. Consider first the leading arm, which consists of material lost from the satellite in the direction towards the host galaxy. Suppose that a particle in the leading arm is almost directly ahead of the satellite in its orbit, so that as perceived from the satellite, the particle and the host are at right angles. As shown in Figure 27, the asymmetric potential causes the gravity on the particle to receive a slight excess contribution in the direction opposite to g ext , which in this case is provided by the host galaxy. Thus, the particle feels an extra force in the radially outward direction as viewed from the host. Due to angular momentum conservation, the angular velocity of the particle then decreases, causing it to move back towards the satellite. This has the effect of shortening the leading arm. The above logic remains exactly the same for the trailing arm, except for the last step -since the material is already trailing, a further decrease in its angular velocity relative to the host increases the length of the trailing arm. This is critical to the results shown in figure 5 of Thomas et al. (2018), and the putative comparison they draw with observations. Note however that the asymmetry relies on having comparable gravity from the satellite and the host. For a very low mass satellite, the host gravity would completely dominate, causing the satellite to experience a dominant external field. This would not only weaken the self-gravity of the satellite, but would also make it much more symmetric with respect to ±g ext . Therefore, tidal tails are not guaranteed to be asymmetric in MOND. Semi-analytic arguments can be used to narrow down the most promising targets to search for such asymmetry.

Tidal tails would be difficult to observe in this level of detail beyond the Local Group, limiting the statistics. Our work suggests another promising line of investigation -external disk galaxies viewed nearly edgeon should curve away from g ext . The recent results of Chae et al. (2020) already indicate that disk galaxies are sensitive to the magnitude of g ext , so a natural extension would be to test whether they are also sensitive to its direction. Comparing the z maps in Figures 22 and26, we see that even when g ext is significantly misaligned with the disk spin axis, the warp is still nearly axisymmetric and mostly in the opposite hemisphere to g ext , with approximately the same strength. * For a much more significant misalignment, g ext would by definition lie very close to the disk plane, making the situation similar to our disk-aligned EFE model. As expected on symmetry grounds, there is no warp in this case, so we do not expect a warp to develop in all disks experiencing a non-negligible EFE.

In cases with a suitable g ext , the main expected signal would be for the disk to warp away from the nearest major galaxy. This would be quite unusual conventionally, since any linear gravity theory is not compatible with a star orbiting around a mass located outside the orbital plane of the star. Even if a galactic disk were initialized to have a bowl-shaped warp, the usual expectation would be for the warp to disappear on a dynamical timescale due to the imbalance of forces along the spin axis. Any observation of such a feature would therefore be quite remarkable, especially if the direction follows expectations based on the large scale gravitational field. This was mapped by Desmond et al. (2018c) and used in the analysis of Chae et al. (2020). While the Desmond et al. (2018c) analysis used conventional methods to estimate g ext on individual galaxies, this was still sufficient for Chae et al. (2020) because an important part of the analysis was comparing deviations from the radial acceleration relation between galaxies experiencing a weak or strong EFE. The latter generally involved a relatively massive nearby galaxy, in which case g ext would be closely aligned with the direction towards it. Moreover, the critical quantity in QUMOND is the Newtonian external field (Equation 28). It should therefore be quite feasible to estimate the direction of g ext on SPARC galaxies by adapting existing work. In principle, this could lead to a detection of the EFE through three main channels:

1. The existence of a nearly axisymmetric bowlshaped warp in edge-on disk galaxies.

2. A correlation between the direction of curvature and that of the external field.

3. The expected magnitude of the curvature could be estimated based on g ext and parameters of the disk, allowing a comparison with observations. * The only difference is a slight central depression for the misaligned case, which could be incorporated into a more sophisticated analysis.

To address the last issue, we use Figure 28 to show the warp profile as a function of r, using annular bins aligned with the disk spin axis. The results could be scaled to different galaxies based on Equations 35 and 36, exploiting the scale invariance of dynamics in the DML (which should be valid as typically g ext 0.1 a 0 , Chae et al. 2020). Since nearby galaxies move over long periods, an important issue is the time required for the warp to develop as a result of the EFE. We address this by showing the warp profile z (r) at different times (Figure 28). From this, it is clear that the warp is already very apparent after just 1 Gyr, before then settling into what is presumably the equilibrium configuration after ≈ 2.5 Gyr. At r = 10 kpc, the orbital period is 2πr/v f = 0.61 Gyr, so the warp is apparent after only a few revolutions. This short timescale means it should be sufficient to predict z (r) using the present g ext . Indeed, MOND simulations of DF2 showed that memory effects due to even a quite strongly time-varying g ext play a rather small role in its internal velocity dispersion, which can be estimated quite accurately by considering it to be in equilibrium with the present g ext (figure 5 of Haghi et al. 2019b). Therefore, orbital motion of galaxies and the resulting changes in g ext would not much affect the above-mentioned disk warp effect. Another very useful feature is that it does not require very precise kinematics of the disk -a rough estimate of the RC amplitude would be helpful, but the test is based mainly on the shape of the galaxy rather than subtle features in its RC. The main observational input would thus be highresolution photographs, with a 0.1 kpc warp in a galaxy 50 Mpc away requiring an angular resolution better than 0.4 .

Searches for warps of this sort have previously been conducted to test fifth-force theories in which galaxies are almost Newtonian and reside in CDM halos (Ferreira 2019). These models generically violate the weak equivalence principle, creating an offset between the baryonic disk and the CDM halo (Desmond et al. 2018b). This would cause the disk to warp, an effect which those authors and Desmond et al. (2018a) claimed to detect. Recent results suggest that the signal was driven by only a small number of galaxies, so an appropriate choice of quality cuts would eliminate them -and therewith the signal (Desmond & Ferreira 2020).

It would be very interesting to revisit these analyses in a MOND context, since the environmental dependence would be rather different. In addition, the same analysis should be applied to mock galaxy images from a hydro ΛCDM cosmological simulation (e.g. Illustris, Pillepich et al. 2018). This would clarify whether ΛCDM contains processes which can create a warp correlated with galaxy and environmental properties in a manner analogous to MOND. In both cases, it would be necessary to apply careful quality control measures, e.g. excluding interacting galaxies significantly affected by tidal forces. One important outcome of our simulations is that the The 1D projection used here loses little information because the warp is nearly axisymmetric (Figure 26). Notice that the warp changes very little after 2.5 Gyr, with the main change before this being the development of a 'knee' at r ≈ (4 -6) kpc. There is no warp initially, so results at this time give an idea of the numerical noise (thin dotted red line).

MOND-predicted warp should be nearly axisymmetric (Figures 22 and 26), but a random perturbation like a satellite would generally affect one side much more than the other. Therefore, systematic errors could be reduced by focusing only on galaxies whose image is nearly symmetric about the sky-projected minor axis. Conventional warps are not symmetric in this sense as the warp material would still orbit the galactic centre.

5. CONCLUSIONS In MOND, disk galaxies are self-gravitating and lack DM halos, leading to different stability properties compared to ΛCDM (Milgrom 1989;Banik et al. 2018a). Our main objective was to see if these differences could help alleviate reported difficulties reproducing the observed dynamical properties of M33 with a live DM halo, where a strong bar tends to form in disagreement with observations (S19). Those authors could stabilize M33 to have a realistic morphology over at least 3 Gyr if the Toomre parameter Q = 2 (Equation 6), which would make it difficult to explain why M33 is currently forming stars (Verley et al. 2009).

To investigate the global stability of M33 in MOND, we set up N -body and hydro simulations of M33 in MOND using a modified version of dice (Perret et al. 2014), which we make publicly available (Section 2.2). We ran them with and without the estimated EFE from M31 using the publicly available por code (Lüghausen et al. 2015). por implements the rather computerfriendly QUMOND formulation of MOND (Milgrom 2010).

We ran isolated stellar-only and hydro simulations for just over 6 Gyr using the same total surface density profile initially. The stellar-only and T = 100 kK simulations become substantially non-circular even in the outer regions, indicating that the whole galaxy is essentially a very long bar (Appendix B) with pattern speed Ω p ≈ 10 km/s/kpc, implying a corotation radius of 10 kpc. These models also develop a substantial central bulge, unlike the observed M33 (Figure 5). This is related to σ LOS being 62 km/s at the centre of M33, much higher than observed.

Most of these problems can be resolved by using a lower gas temperature of 25 kK, which is more in line with the 12 kK adopted in section 3.2 of S19. In this cooler isolated model, no bulge forms. The stellar disk appears fairly circular in the outer regions (Figure 8), with a clear bisymmetric spiral evident in the gas when viewed face-on (Figure 9). A strong bar forms in the first Gyr, but it then loses strength for a variety of reasons, including buckling, radial gas inflow, and chaos (see also Tiret & Combes 2008a). As a result, the observed weak bar of M33 is recovered by the end of this simulation. Our Fourier analysis indicates that the dominant mode of non-axisymmetry is m = 1, with a significant m = 2 component but very little strength in higher harmonics (Figure 12). The A 2 /A 0 ratio for material with r = (0.5 -3) kpc is sometimes similar to the observed 0.2 (section 4.3 of Corbelli & Walterbos 2007), though it is generally smaller. Since the outskirts appear nearly circular, the bar is genuinely weak. Its pattern speed of 30 km/s/kpc implies a corotation radius of 3 kpc, so the bar is fairly short. Despite these promising results, there is still significant radial redistribution of material, causing the inner RC of M33 to rise more steeply than observed (Figure 4). The central σ LOS of 57 km/s is also above the observational range of (28 -35) km/s (section 3.2 of Corbelli & Walterbos 2007). We conclude that allowing a more dissipative gas component greatly helps to stabilize the thin observed disk of M33 against the formation of a very strong bulge and bar, but further improvements to the model are still necessary.

In MOND, the EFE from M31 can have a percentlevel effect on the dynamics of M33. We included the EFE with g ext = 0.07 a 0 at 30 • to the M33 disk. We demonstrated for the first time that this affects its secular evolution after a few Gyr, with the inward radial migration of material significantly suppressed (Figure 14). This makes the RC much closer to observations (Figure 13), with the minor ( 10 km/s) differences possibly being due to a warp that is not included when deriving the observed RC. The central σ LOS is also much lower at 48 km/s by the end of our simulation (9.9 Gyr), with the LOS velocity distribution closely following a Gaussian of this width (Figure 17). This is much more consistent with the observational range. The non-circular motions are also much smaller (Figure 20).

To test the numerical convergence of these results, we ran a higher resolution simulation for 6.1 Gyr. No central bulge formed (Figure 23), consistent with our lower resolution model and our isolated 25 kK simulation at this time (Figure 5). However, the isolated 100 kK model develops a substantial bulge. We therefore argue that the observed configuration of M33 is inherently stable in MOND if a realistic gas temperature is adopted, though the EFE helps to improve the agreement with some observables. Despite the low observed gas fraction near the centre, its dissipative nature is seemingly important to our main result that no central bulge forms, in agreement with observations (e.g. Kam et al. 2015). Since an important reason for the success of our model is a reduction in T from 100 kK to 25 kK, a colder gas component may help to further stabilize the stellar disk. However, our results from a very cold 10 kK model indicate that the disk becomes unstable once T is reduced this much. It is unclear whether stellar feedback would change this picture, though the ejection of material should make it even more difficult to form a substantial bulge. A reduced central surface density would cause the RC to rise more gradually in the central few kpc, which can also be achieved by starting with a larger scale length than observed. This would make the disk more stable by pushing it deeper into the MOND regime (Equation 3).

In general, the EFE breaks both the axisymmetry and up-down symmetry of the underlying gravitational physics. To isolate the effect of each in turn, we ran simulations where g ext is aligned with the disk spin axis or put in the disk plane (Section 4.4). The suppression of radial migration is similar for all three considered external field orientations, so we conclude that this arises from breaking the combination of axisymmetry and updown symmetry present only in the isolated case. In models with the EFE where one of these symmetries is preserved, the restoration of some symmetry prevents the disk as a whole from precessing, which it does with intermediately aligned g ext (Figure 21). The disk also becomes warped if the up-down symmetry is broken. We relate the warp to previous DML calculations of the gravity from a point mass in a weak external field, focusing on the extent to which results are asymmetric with respect to ±g ext (Figure 27). Possible observational signatures of this asymmetry are discussed in Section 4.5, focusing in particular on the induced warp in disk galaxies viewed close to edge-on (Section 4.5.1). The warp develops in only a few dynamical times (Figure 28) in the sense that the outskirts curve towards -g ext or away from the mass sourcing the EFE, potentially offering a clear and highly specific observational signature of the EFE. This would extend the previous result of Chae et al. (2020) that its magnitude has a detectable impact on galaxy RCs, violating the strong equivalence principle because the reported effect is not due to tides in conventional gravity (see their section 4).

MOND has enjoyed a great deal of success predicting the RCs of spiral galaxies, including M33 (Sanders 1996;Famaey & McGaugh 2012;Kam et al. 2017). It can also form such galaxies out of a collapsing gas cloud (Wittenburg et al. 2020), which may well fit into a broader cosmological context that better explains the large scale structure and expansion rate of the local Universe without violating early Universe constraints (Haslbauer et al. 2020). Using a numerical implementation of MOND, we were able to match the leading-order non-axisymmetric features of M33 reasonably well once its gas is assigned a not too high temperature. Our model produces no significant central bulge, but the simulated RC still rises more steeply than observed due to efficient inward radial migration. This process is hampered when we include the MOND-predicted weak EFE representative of the estimated gravity exerted by M31. This yields a much better representation of the observed RC. The central velocity dispersion σ LOS is however still too high, but this could be due to possible systematics in the measurements and their comparison to our simulations (Section 4.2.2). Further improvements to our models should be possible by e.g. including stellar feedback to reduce the central surface density, or forcing this to a lower initial value by starting with a more extended mass distribution.

Our work highlights for the first time the role of a weak external field on the stability and evolution of disk galaxies in MOND. Further simulations with a time-varying external field, modeling the full orbit of M33, will be needed to confirm, and perhaps refine, its resemblance to observations.

Figure 1 .

 1 Figure1. The total surface density profile of our M33 model (black) and its component parts -stars (red) and gas (blue). The horizontal magenta solid line near the top represents the MOND critical surface density (Equation3). M33 has a lower surface density almost everywhere, so MOND should greatly impact its dynamical properties. The lower dotdashed magenta line shows the vertical component of the Newtonian-equivalent external field from M31 if we put this at 30 • to the disk plane (Equation26). Notice that the EFE could have a percent-level effect on the gravitational field in the central regions of M33 (Section 4).

  Figure3. Initial sech 2 scale height of our simulated M33 gas disk as a function of r, with the adopted temperature shown in the legend. This thickness is much larger than assumed when calculating the RC in dice (Section 1). This weakens the radial gravity at intermediate radii, implying that our hydro simulations are not set up exactly in radial equilibrium (see text).

Figure 4 .

 4 Figure4. The simulated RC of M33 in our isolated hydro models based on the radial gravity felt by particles within 1 kpc of its disk plane (Equation20). The top (bottom) panel shows results for T = 100 kK (25 kK). Different curves correspond to different times, which we indicate in the legend. The data points are fromKoch et al. (2018a) based on an M33 distance of 840 kpc(Kam et al. 2015, and references therein). The grey line represents an angular frequency of 10 km/s/kpc (30 km/s/kpc) for our 100 kK (25 kK) model, which is roughly the bar pattern speed in each case (Section 3.3.2). Although the 100 kK and 25 kK RCs are quite similar, the disks have very different morphologies (Section 3.2).

Figure 5 .

 5 Figure 5. The stellar particles in our isolated hydro simulation of M33 with T = 100 kK (top) and 25 kK (bottom), shown using a cylindrical rz projection (Equation 21). Values are shown in units of 1000 M /pc 2 . The time elapsed since the start of the simulation is shown in the top left corner of each frame. Analogous results for the gas are shown in Appendix A. Notice the significant central bulge in the hotter model, and its absence in the cooler model.

Figure 6 .

 6 Figure 6. Time evolution of the simulated vertical velocity dispersion σz of material with r = (0.5 -3) kpc. 'Stellaronly' refers to our purely stellar N -body simulation. Other results are for hydro simulations at 25 kK, except for the pink dotted line (100 kK). The dotted black curve shows results when including the EFE (Section 4), but other models are isolated. The 1D velocity dispersion of the gas is σg ≈ 22 km/s at T = 100 kK but only σg ≈ 11 km/s at T = 25 kK, possibly explaining the difference in stability properties of the stellar disk (see text).

Figure 7 .

 7 Figure 7. The line of sight velocity dispersion of stars in different parts of the sky-projected image of M33, shown for the final snapshot at 7.4 Gyr for our isolated 25 kK simulation. Thin ellipses show the projected appearance of circles in the disk plane with radii of 2, 4, and 6 kpc, approximately integer multiples of the stellar disk scale length (Table 1). The magenta dotted (dot-dashed) contours show σ LOS = 35 (28) km/s, bracketing the observational range (section 3.2 of Corbelli & Walterbos 2007). The central σ LOS is 57 km/s, well above this range. Similar results are obtained if viewing M33 from within the yz plane, with a central σ LOS differing by < 0.1 km/s (not shown). Values are ≈ 5 km/s higher for the final snapshot of our isolated 100 kK simulation (not shown). In both cases, the v LOS values in the central pixel closely follow a Gaussian with the calculated σ LOS , as shown explicitly for a better fitting model in Section 4.2.2.

Figure 8 .

 8 Figure 8. Face-on view of the stellar mass distribution in our isolated hydro model of M33 at 25 kK, with the surface density Σ shown in 1000 M /pc 2 . Results for the analogous model at 100 kK are shown in Appendix B. Despite the significant difference in morphology, the RCs are nearly identical (Figure 4).

Figure 9 .Figure 10 .

 910 Figure 9. Similar to Figure 8, but for the gas. Notice the bisymmetric spiral towards the end. Results for the 100 kK model are shown in Appendix B.

Figure 11 .

 11 Figure11. Bar pattern speed for stars with r = (0.5 -3) kpc in our isolated stellar-only and 25 kK hydro simulations (solid red and dotted blue curves, respectively). Our 100 kK simulation (dotted pink) behaves similarly to the stellar-only case. The black curve shows the pattern speed of the stars in our simulation with the EFE (Section 4). In all hydro simulations, results are very similar when using the total mass distribution instead (not shown).

Figure 12 .

 12 Figure 12. Strengths of different Fourier modes (Equation23) for stellar particles at r = (0.5 -3) kpc in our isolated 25 kK hydro model. Some modes can have strength > 1 if multiple modes are present and their phases (Equation25) are fortuitously aligned. Results remain rather similar when including the gas (not shown).

4.

  INCLUDING THE EFE FROM M31 In MOND, M31 can affect M33 through the EFE (Section 1.1). Since the M31 RC has v f ≈ 225 km/s (Carignan et al. 2006) and its distance to M33 is d ≈ 200 kpc (table 1 of Patel et al. 2017), the external field on M33 is mainly the M31 gravity g ext = v f

Figure 13 .

 13 Figure 13. Similar to Figure 4, but now showing the RC of M33 in our model with the intermediately aligned EFE. The final RC (dotted black) is now quite close to the observations.

Figure 14 .

 14 Figure 14. The stellar surface density distribution in our isolated models at 100 kK (red dotted line) and 25 kK (blue dotted line). Results with the EFE included according to Equation 26 are shown in solid black, with other black lines indicating special external field orientations (Section 4.4).The initial distribution (Figure1) is shown with a featureless pink line.

Figure 15 .Figure 16 .Figure 17 .

 151617 Figure15. Similar stellar rz view to Figure5, but for our hydro simulation with the EFE. Notice the absence of a central bulge even after 10 Gyr, similarly to our analogous isolated model (bottom panel of Figure5). The bifurcation at large r is caused by precession of the whole disk (Section 4.2.5).

Figure 18 .Figure 19 .

 1819 Figure18. Similar stellar xy view to Figure8, but for our hydro model with the EFE. In both cases, the outer regions look nearly circular after 10 Gyr, indicating a genuinely weak bar and rather little non-circular motion (Section 4.2.4).

Figure 20 .

 20 Figure20. Face-on view of M33 showing the mean tangential velocity within its disk plane less the average for particles at similar r (Equation33) at the end of our 25 kK hydro simulation in isolation (top) and with the EFE set by Equation 26 (bottom). The curves show circles of radius 2, 4, and 6 kpc. Our binning and outlier rejection procedures are described in Section 3.2.1. The disks have been rotated into the plane defined by their overall angular momentum, which is important for our model with the EFE (Section 4.2.5). Both panels use the same colour scale, highlighting the much smaller non-circular motions in our model with the EFE. If we use snapshots at different times (not shown), the pattern in the isolated model appears to rotate, while results the EFE change very little and do not circulate.

Figure 21 .

 21 Figure 21. Time evolution of the disk spin axis h in our model with the fiducial EFE (changes are negligible in other models, not shown). The solid red and dot-dashed magenta curves show the individual components of h. The solid black lineshows the total change in h from its initial orientation (along +z), while the dotted green line shows how much of this is a nutation rather than precession around gext. The behaviour of h is very similar in our higher resolution simulation with the EFE (Section 4.2.6), so we show only the total precession angle in this case (solid blue). The precession is almost entirely towards -y in both cases (solid black and dot-dashed magenta curves coincide). The direction and approximate magnitude of the precession are explained in the text using analytic arguments.

Figure 22 .

 22 Figure22. The mean height z of stars in the M33 disk after 9.9 Gyr, found after applying a rotation as shown in Figure21. The circles have radii of 2, 4, and 6 kpc. The external field is in the initial xz plane (Equation26).

Figure 23 .

 23 Figure 23. Similar to Figure 8, but now showing the stellar rz projection at 6.1 Gyr for our lower resolution (top) and higher resolution (bottom) simulation with the EFE from Equation26. The latter has the same resolution settings and box size as our isolated simulations, but was not advanced further due to significant numerical drift of M33. The slightly smaller disk precession in our higher resolution simulation (Figure21) causes the disk to appear thinner at large r.

Figure 24 .

 24 Figure24. Similar to Figure6, but now showing σz for the stars at r = (0.5 -3) kpc in our hydro simulation with the EFE at 30 • to the disk (Equation26) in which the floor on the Toomre Q parameter is varied to Q lim = 1 (solid red), Q lim = 1.1 (dot-dashed blue), or has the fiducial value of Q lim = 1.25 (black) used in all other simulations. The simulation at Q lim = 1 appears to become unstable after ≈ 5 Gyr, while the other cases remain stable. Results are very similar if the full mass distribution is used instead of just the stars (not shown).

Figure 25 .

 25 Figure25. Similar to Figure12, but now showing Fourier analysis of the stars at the end of our hydro simulation with axially aligned gext. This is the only model to avoid a dominant m = 1 mode at all times. Its strength in the gas is also rather weak (not shown).

Figure 26 .

 26 Figure 26. Similar to Figure 22, but now showing z for our simulation with axially aligned gext ∝ +z. Notice that the warping is now axisymmetric, and occurs in the direction opposite to gext. Results remain fairly similar despite the orientation of gext differing by 60 • .

Figure 28 .

 28 Figure 28. The mean height z of different annuli in the M33 disk in an axially aligned external field. Different colours and linestyles show different times, as indicated in the legend.The 1D projection used here loses little information because the warp is nearly axisymmetric (Figure26). Notice that the warp changes very little after 2.5 Gyr, with the main change before this being the development of a 'knee' at r ≈ (4 -6) kpc. There is no warp initially, so results at this time give an idea of the numerical noise (thin dotted red line).

Figure 29 .

 29 Figure 29. Similar to Figure5, but now showing the gas in our isolated simulation at 100 kK (top) and 25 kK (bottom). The thinner gas disk in the cooler model is expected from the initial thickness profile (Figure3), though some thickening occurs in both models.

Figure 30 .Figure 31 .

 3031 Figure 30. Similar to Figure 8, but now showing the particles in our isolated simulation at 100 kK. Notice the significantly more non-circular shape of M33 in its outskirts, signifying a very strong and extended bar quite inconsistent with observations.

Table 1 .

 1 Parameters of our M33 model, adopted from the MOND RC fit in

	Component	Mass	Scale length Aspect ratio
	Stellar disk 4.7 × 10 9 M	1.78 kpc	0.1
	Gas disk	1.8 × 10 9 M	4 kpc	0.1

* This is not true of QUMOND if applied rigorously since it is derivable from a Lagrangian, therewith obeying the usual conservation laws(Milgrom 2010).

* They concluded that the uncertainties were probably not enough to accommodate σ LOS = 43 km/s.

* The fact that the tangential component of the gravity is antialigned with gext is also evident in figure2ofBanik & Kroupa (2019), though their results are for a much stronger gext = 1.78 a 0 as appropriate for the Solar neighbourhood.

ACKNOWLEDGEMENTS

IB is supported by an Alexander von Humboldt Foundation postdoctoral research fellowship. BF and RI acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 834148). GC acknowledges support from FONDECYT Regular No. 1181708. IB thanks Stacy McGaugh for providing the M33 baryonic mass model and RC measurements. He also thanks Jerry Sellwood for providing bar strength time series in several Newtonian simulations of M33 embedded in a DM halo. The authors are very grateful to the referee for helping to significantly improve this manuscript and the simulations it presents.

APPENDIX

A. EFFECT OF TEMPERATURE ON THE GAS DISK Figure 29 shows the gas in our isolated hydro simulations at different times using the cylindrical rz projection (Section 3.2). In the outskirts, the gas expands outward to a small extent in both models. Meanwhile, the density in the central pixel rises slightly, especially in the hotter (100 kK) model. Notice that the appearance changes little after the first Gyr, which is already enough to discern a substantial difference between the models. The retention of a thin gas disk in the cooler model is important to the global stability of M33, as discussed in the main article.

B. FACE-ON VIEWS OF THE 100 KK ISOLATED SIMULATION Figure 30 shows the stars in our isolated 100 kK simulation as viewed face-on at different times. As discussed in Section 3, the outer parts appear significantly non-circular, to the extent that the whole galaxy is essentially one giant bar. This is consistent with the corotation radius of 10 kpc implied by the bar pattern speed (Figure 11).

The gas in this simulation is shown in Figure 31 using the same face-on view. The high temperature causes a distinct lack of small-scale features, making the appearance quite different to our cooler isolated model at 25 kK (Figure 9).