The Global Stability of M33 in MOND

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 ${{\rm{sech}} }^{2}$ vertical profile with characteristic thickness h_g given by

$$\begin{eqnarray}&&{{c}_{s}}^{2}=\pi {{Gh}}_{g}{{\rm{\Sigma }}}_{g}=\displaystyle \frac{{g}_{N,z}{h}_{g}}{2},\end{eqnarray} \tag{ 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 ${{\rm{sech}} }^{2}$ vertical profile remains a good description if we modify Equation (8) to

$$\begin{eqnarray}&&{{c}_{s}}^{2}=\displaystyle \frac{{g}_{N,z}{h}_{g}}{2}\nu \left(\sqrt{{g}_{N,r}^{2}+{g}_{N,z}^{2}}\right).\end{eqnarray} \tag{ 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 ${\widetilde{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

$$\begin{eqnarray}&&{{c}_{s}}^{2}\equiv \displaystyle \frac{{\widetilde{g}}_{N,z}{h}_{g}}{2}\nu \left(\sqrt{{{g}_{N,r}}^{2}+{{\widetilde{g}}_{N,z}}^{2}}\right).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{\widetilde{g}}_{N,z}=2\pi G\left({{\rm{\Sigma }}}_{g}+{{\rm{\Sigma }}}_{* }\right),\quad {h}_{g}={h}_{* }.\end{eqnarray} \tag{ 11 }$$

In general, h_g  ≠ h_*. We make the assumption that each component still follows a ${{\rm{sech}} }^{2}$ vertical profile. We then find what fraction of the stellar disk is enclosed at $\left|z\right|\lt {h}_{g}$. For the gas disk, this fraction is $\tanh \left(1\right)$. These fractions are unequal if h_g  ≠ h_*. In this case, we assume that the appropriate generalization of Equation (11) is

$$\begin{eqnarray}&&{\widetilde{g}}_{N,z}=2\pi G\left({{\rm{\Sigma }}}_{g}+{{\rm{\Sigma }}}_{* }\displaystyle \frac{\tanh \left({h}_{g}/{h}_{* }\right)}{\tanh \left(1\right)}\right).\end{eqnarray} \tag{ 12 }$$

This approach can easily be generalized to stellar disks with multiple components that have different scale heights and/or vertical profiles different from ${{\rm{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 center of the galaxy. Following our previous approach, the important quantity is the resulting vertical Newtonian gravity at a height h_g . This is approximately given by ${g}_{N,r}\left({h}_{g}/r\right)$ 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 \left(r/{h}_{* }\right)$, leading to our final approximation that

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{g}}_{N,z} & = & 2\pi G\left({{\rm{\Sigma }}}_{g}+{{\rm{\Sigma }}}_{* }\displaystyle \frac{\tanh \left({h}_{g}/{h}_{* }\right)}{\tanh \left(1\right)}\right)\\ & & +{g}_{N,r}\displaystyle \frac{{h}_{g}}{\sqrt{{r}^{2}+{{h}_{g}}^{2}}}\tanh \left(\displaystyle \frac{r}{{h}_{* }}\right).\end{array}\end{eqnarray} \tag{ 13 }$$

We can now solve Equation (10) by varying h_g using a Newton–Raphson algorithm, each time recalculating ${\widetilde{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}\left(r\right)$ shown in Figure 3, which we incorporate into our hydro simulations as described next.

Zoom In Zoom Out Reset image size

We set up the gas density and velocity field in por using our estimated gas disk thickness profile ${h}_{g}\left(r\right)$ 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 ${{\rm{sech}} }^{2}$ vertical profile with characteristic height ${h}_{g}\left(r\right)$. As a result, the gas density is

$$\begin{eqnarray}&&{\rho }_{g}\left(r\right)=\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{2\,{h}_{g}}{{\rm{sech}} }^{2}\left(\displaystyle \frac{z}{{h}_{g}}\right).\end{eqnarray} \tag{ 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}\equiv \sqrt{-{{rg}}_{r}}$ due to pressure gradients. The equation of radial hydrostatic equilibrium is

$$\begin{eqnarray}&&{{v}_{c}}^{2}={{v}_{a}}^{2}-{{c}_{s}}^{2}\left(\displaystyle \frac{\partial \mathrm{ln}{\rho }_{g}}{\partial \mathrm{ln}r}\right).\end{eqnarray} \tag{ 15 }$$

Because h_g depends on r, it is nontrivial to ensure even radial hydrostatic equilibrium of the gas. We approximately achieve this in the disk midplane by realizing that ρ_g  ∝ Σ_g /h_g (Equation (14)), implying that

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}{\rho }_{g}}{\partial \mathrm{ln}r}=\displaystyle \frac{\partial \mathrm{ln}{{\rm{\Sigma }}}_{g}}{\partial \mathrm{ln}r}-\displaystyle \frac{\partial \mathrm{ln}{h}_{g}}{\partial \mathrm{ln}r},\quad z=0.\end{eqnarray} \tag{ 16 }$$

The radial gradient of h_g is found by centered differencing. We avoid numerical difficulties at the origin by exploiting the fact that

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}{h}_{g}}{\partial \mathrm{ln}r}=\displaystyle \frac{r}{{h}_{g}}\displaystyle \frac{\partial {h}_{g}}{\partial r}.\end{eqnarray} \tag{ 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.

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 toward 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).

Zoom In Zoom Out Reset image size

The difference in behavior due to the gas temperature can be understood by relating the above values to the 1D velocity dispersion of the gas. This is ${\sigma }_{g}=\sqrt{{kT}/\left(\mu {m}_{p}\right)}=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 kiloparsecs 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 $\left(\sin i,0,\cos i\right)$. 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 center of M33 coincident with that of the central pixel.

Since each pixel contains a finite number of particles, there is a nonzero correlation between their mean LOS velocity ${\overline{v}}_{\mathrm{LOS}}$ and that of each particle. The real M33 has many more particles, so this is a numerical artifact that should be corrected. The standard approach is to enhance the measured variance by a factor of $N/\left(N-1\right)$ 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

$$\begin{eqnarray}&&{\sigma }_{\mathrm{LOS}}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{m}_{i}{\displaystyle \sum }_{i=1}^{N}{m}_{i}{v}_{{}_{\mathrm{LOS},i}}^{2}-{\left({\displaystyle \sum }_{i=1}^{N}{m}_{i}{v}_{\mathrm{LOS},i}\right)}^{2}}{{\left({\displaystyle \sum }_{i=1}^{N}{m}_{i}\right)}^{2}-{\displaystyle \sum }_{i=1}^{N}{{m}_{i}}^{2}}}.\end{eqnarray} \tag{ 22 }$$

Once we have obtained ${\overline{v}}_{\mathrm{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 center 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).

Zoom In Zoom Out Reset image size

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 noncircular 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. (2018a) 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 midplane. 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 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 T2_star (Section 2.4).

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

$$\begin{eqnarray}\begin{array}{rcl}{A}_{0} & \equiv & \displaystyle \frac{1}{2}\displaystyle \sum _{i}{m}_{i}\,\left(i\,{\rm{l}}{\rm{a}}{\rm{b}}{\rm{e}}{\rm{l}}{\rm{s}}\,{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}{\rm{e}}{\rm{r}}{\rm{e}}{\rm{n}}{\rm{t}}\,{\rm{p}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{i}}{\rm{c}}{\rm{l}}{\rm{e}}{\rm{s}}\right),\\ {A}_{m} & \equiv & \sqrt{{\left({\displaystyle \sum }_{i}{m}_{i}\cos m{\phi }_{i}\right)}^{2}+{\left({\displaystyle \sum }_{i}{m}_{i}\sin m{\phi }_{i}\right)}^{2}},\,m\geqslant 1.\end{array}\end{eqnarray} \tag{ 23 }$$

This quantifies the presence of non-axisymmetry, i.e., dependence on the azimuthal angle $\phi \equiv {\sin }^{-1}\left(y/r\right)$. We then find the strength of the mth Fourier mode using

$$\begin{eqnarray}&&{\textstyle \text{Strength of Fourier mode}}\,m\equiv \displaystyle \frac{{A}_{m}}{{A}_{0}}.\end{eqnarray} \tag{ 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).

Zoom In Zoom Out Reset image size

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 behavior 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₂/A₀ ≈ 0.2, similar to that observed in M33. There are extended periods with a very weak bar (A₂/A₀ ≲ 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₂/A₀ ≈ 0.3, see Figure 10). Including gas at 100 kK leads to similar overall behavior to the stellar-only model, but with the final A₂/A₀ ≈ 0.2 (not shown). In both cases, the A₂/A₀ 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 semiminor axis comparable to the 3 kpc aperture used to calculate A₂/A₀ (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.

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

$$\begin{eqnarray}&&\tan \left(m{\theta }_{m}\right)\equiv \displaystyle \frac{{\sum }_{i}{m}_{i}\sin \left(m{\phi }_{i}\right)}{{\sum }_{i}{m}_{i}\cos \left(m{\phi }_{i}\right)}.\end{eqnarray} \tag{ 25 }$$

We then consider how θ₂ 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 scale lengths of the stellar and gas disks (Table 1). Since our calculation of Ω_p 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.

Zoom In Zoom Out Reset image size

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₂₅, with their table 1 indicating ${r}_{25}=30^{\prime}$ 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\pm 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⁵ 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\,\pm 2$ km s⁻¹ kpc⁻¹ given in Section 4.3 of Corbelli & Walterbos (2007). Those authors assumed that 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. (2018b) 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 estimate of pattern speed, 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.

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.

Zoom In Zoom Out Reset image size

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, 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 two-armed spiral in the final snapshot of the gas. In a more advanced simulation, this may well cause a prominent 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 analog to the EFE, which violates the strong equivalence principle.

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 Newtonian-equivalent 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)

$$\begin{eqnarray}&&{g}_{N,\mathrm{ext}}=\displaystyle \frac{{{g}_{\mathrm{ext}}}^{2}}{{g}_{\mathrm{ext}}+{a}_{0}}=4.6\times {10}^{-3}\,{a}_{0}.\end{eqnarray} \tag{ 27 }$$

To retain our previous notation that g _N refers to the Newtonian gravity of M33 alone, we generalize Equation (18) to

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{g}}={\rm{\nabla }}\cdot \left[{{\boldsymbol{g}}}_{N,\mathrm{tot}}\nu \left({{\boldsymbol{g}}}_{N,\mathrm{tot}}\right)\right],\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{N,\mathrm{tot}}\equiv {{\boldsymbol{g}}}_{N}+{{\boldsymbol{g}}}_{N,\mathrm{ext}}.\end{eqnarray} \tag{ 29 }$$

In a more detailed model, g _N and g _N,ext would be directed toward 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 behavior 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 dominated by the external field, 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

$$\begin{eqnarray}&&{\rm{\Phi }}=-\displaystyle \frac{{GM}{\nu }_{\mathrm{ext}}}{R}\left(1+\displaystyle \frac{{K}_{0}}{2}{\sin }^{2}\theta \right),\end{eqnarray} \tag{ 30 }$$

where θ is the angle between g _N,ext and the position R , while K₀ was defined earlier (Equation (7)). This solution was confirmed numerically by direct integration of Equation (18) (Section 2.2 of Banik & Zhao 2018a). The external field dominates for distances R ≫ R_EFE, where the transition radius in the DML is

$$\begin{eqnarray}&&{R}_{\mathrm{EFE}}=\displaystyle \frac{\sqrt{{{GMa}}_{0}}}{{g}_{\mathrm{ext}}}.\end{eqnarray} \tag{ 31 }$$

Since the mass of M33 is M = 6.5 × 10⁹ M_⊙ (Table 1), the external field from M31 dominates beyond R_EFE = 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_EFE 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.

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 dominated by the external field. 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¹² 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⁻³ a₀. The results of the two 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⁻⁴ a₀. This is probably due to the larger box size making Equation (30) more accurate on the boundary. ¹⁴

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

Zoom In Zoom Out Reset image size

Since the EFE is not very significant at r ≲ 2 kpc (Figure 1), changes in the RC arise mainly from differences 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 that 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 toward lower r. This is not much dependent on the choice of direction for the EFE, which we vary in Section 4.4.

Zoom In Zoom Out Reset image size

Our simulation with the "fiducial" g _ext (Equation (26)) still displays a small RC discrepancy at r ≲ 4 kpc. In addition to changes in the distance of M33 at the level of a few per cent, 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 gigayears 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 recognizable simply by comparing its approaching and receding halves, though doing so for M33 does reveal differences of this order (Figure 10 of Kam et al. 2015). Other issues include a slight r-dependent inclination, i.e., a warp. Koch et al. (2018a) obtained i = 5508 ± 156 assuming fixed i for the whole disk (see their Table 2). If i = 50° at the center (as in the earlier work of Corbelli & Salucci 2007), the simulated RC would appear to exceed the observed RC by 7% even if the simulation matches the actual RC. 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 the inclination could be 56° at r = 5 kpc, even though it is only 50° at the center. 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 subdominant 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).

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

Zoom In Zoom Out Reset image size

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, ${\sigma }_{{\rm{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⁻¹.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 nonzero component toward + x (Equation (26)). The potential of a point mass is symmetric with respect to g _ext when g _ext dominates (Equation (30)), but there is an asymmetry between ± g _ext if the 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, 2017). At large distances, g _ext becomes dominant and this asymmetry vanishes. The asymmetry is discussed further in Section 4.4.

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.

Zoom In Zoom Out Reset image size

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 noninteracting 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)).

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{v}_{t}\equiv \displaystyle \frac{{{xv}}_{y}-{{yv}}_{x}}{r}.\end{eqnarray} \tag{ 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 $\overline{r}$ and tangential velocity ${\overline{v}}_{t}$. We interpolate this to obtain the expected ${\overline{v}}_{t}$ at any r ≲ 25 kpc. We then subtract ${\overline{v}}_{t}$ to obtain the excess tangential velocity:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{t}\,\equiv \,{v}_{t}-{\overline{v}}_{t}.\end{eqnarray} \tag{ 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, noncircular motions should mostly be restricted to the central 3 kpc, which is approximately the case. The noncircular 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.

Zoom In Zoom Out Reset image size

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 $\left|z\right|\lt 20\,\mathrm{kpc}$. The direction of this vector is the disk spin axis, which we denote $\widehat{{\boldsymbol{h}}}$. The time evolution of $\widehat{{\boldsymbol{h}}}$ is shown in Figure 21, which reveals a very slow precession by just over 1° Gyr⁻¹.

Zoom In Zoom Out Reset image size

Since the external field picks out only one preferred direction (Equation (26)), we expect $\widehat{{\boldsymbol{h}}}\cdot {{\boldsymbol{g}}}_{\mathrm{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 $\widehat{{\boldsymbol{h}}}$. Thus, we can consider $\widehat{{\boldsymbol{h}}}$ as precessing around g _ext. Since g _ext is within the x z plane (Equation (26)), this implies that $\widehat{{\boldsymbol{h}}}$ mostly precesses within the y z plane. In our simulation, the precession is toward  − y . Our isolated model of M33 exhibits only a very small amount of precession (<001), 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, the potential is deepest along the g _ext axis. As a result, stars at x > 0 would feel an extra force toward  + z , while stars at x < 0 would feel a force toward  − 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 x z 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 $\widehat{{\boldsymbol{h}}}$ to precess toward  − 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 (30) is 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 nonzero only if we consider the EFE. Their result is based on direct numerical integration of Equation (28) for 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

$$\begin{eqnarray}&&\dot{{\rm{\Omega }}}=\displaystyle \frac{{f}_{\mathrm{geom}}{{GMa}}_{0}r\sin \theta \cos \theta }{2{v}_{c}\,{{r}_{\mathrm{ext}}}^{3}},\end{eqnarray} \tag{ 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 $\left(x=r,y=z=0\right)$ is not representative over all ϕ. Note that only precession toward  − y is considered, because other components cancel out by symmetry. Equation (34) suggests a precession rate of 166 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 11 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 $\widehat{{\boldsymbol{h}}}$ before doing further analyses. Fortunately, the effect is rather small because $\cos 10^\circ =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). This is because the mean height of the disk is $\overline{z}\propto \sin \phi$, 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 $\overline{z}$ is stationary with respect to ϕ, which occurs when ϕ = ±π/2 and $\left|z\right|$ 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.

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₀, 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 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 $\overline{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 coordinate system with $\widehat{{\boldsymbol{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 $\widehat{{\boldsymbol{h}}}$ and g _ext, we expect the x z 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.

Zoom In Zoom Out Reset image size

The mild warping evident in Figure 22 is not sufficient to explain the  ≈5° warp inferred by Corbelli & Salucci (2007)—a change in $\overline{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 from M31. Tidal effects might also have been important—tides naturally pull down on one side of M33 while pulling up on the other, 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. 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). The planes of their satellite galaxies (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).

Our simulation with the EFE uses a lower resolution than 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 high-resolution hydro simulation with the EFE for 6.1 Gyr due to significant numerical drift of M33, which is much 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, the two 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).

Zoom In Zoom Out Reset image size

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 affects 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).

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), our work 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 toward 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 toward 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 satellite of very low mass, 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 edge-on 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 $\overline{z}$ maps in Figures 22 and 26, 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. (2018a) and used in the analysis of Chae et al. (2020). While the analysis of Desmond et al. (2018a) 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 toward 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 bowl-shaped 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.

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₀, 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 $\overline{z}\left(r\right)$ 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 $\overline{z}\left(r\right)$ 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 high-resolution photographs, with a 0.1 kpc warp in a galaxy 50 Mpc away requiring an angular resolution better than 0.4''.

Zoom In Zoom Out Reset image size

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 that those authors and Desmond et al. (2018c) 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 that 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 MOND-predicted warp should be nearly axisymmetric (Figures 22 and 26), but a random perturbation such as 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 because the warp material would still orbit the galactic center.