# Tidally Induced Morphology of M33 in Hydrodynamical Simulations of Its Recent Interaction with M31

Marcin Semczuk, Ewa Łokas, Jean-Baptiste Salomon, E. Athanassoula, Elena D'onghia

► To cite this version:

Marcin Semczuk, Ewa Łokas, Jean-Baptiste Salomon, E. Athanassoula, Elena D'onghia. Tidally Induced Morphology of M33 in Hydrodynamical Simulations of Its Recent Interaction with M31. The Astrophysical Journal, 2018, 864 (1), pp.34. 10.3847/1538-4357/aad4ae. hal-02021051

HAL Id: hal-02021051

https://hal.science/hal-02021051v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Tidally Induced Morphology of M33 in Hydrodynamical Simulations of Its Recent Interaction with M31

Marcin Semczuk<fnref n="1" />, Ewa L. Łokas<fnref n="1" />, Jean-Baptiste Salomon<fnref n="1" />,<fnref n="2" />, E. Athanassoula³, and Elena D'Onghia<fnref n="4" />,<fnref n="5" />,<fnref n="6" />

Received 2018 April 12; revised 2018 June 26; accepted 2018 July 17; published 2018 August 29

## Abstract

We present a hydrodynamical model of M33 and its recent interaction with M31. This scenario was previously proposed in the literature in order to explain the distorted gaseous and stellar disks of M33, as well as the increased star formation rate in both objects around 2 Gyr ago. We used an orbit integration scheme to find which estimate of the transverse velocity of M31 favors the interaction scenario more and then tried to reproduce it in our simulations. M33 was modeled as a stellar and gaseous disk embedded in a live dark matter halo, while M31 was approximated only with a live dark halo. In the simulations, the two galaxies passed each other with a pericenter distance of 37 kpc. Tides excited a two-armed spiral structure in the M33 disk, which is found to be the predominant spiral signal in the observed galaxy and has long been known as a feature easily induced by tidal interactions. We found that the gaseous warp produced by the interaction did not resemble enough the observed one, and we performed an additional simulation including the hot gas halo of M31 to show that this feature can be properly reproduced by tidal forces and ram pressure stripping acting simultaneously on the gaseous disk. In addition to the spiral arms, tidal forces produced a stellar stream similar to the one observed and triggered a star formation burst at radii similar to where it is observed.

**Key words:** galaxies: evolution – galaxies: individual (M33) – galaxies: interactions – galaxies: kinematics and dynamics – galaxies: structure – Local Group

## 1. Introduction

The Triangulum Galaxy (M33) is a late-type spiral and the third-largest galaxy in the Local Group (LG), behind the two more massive members of the LG, the Milky Way (MW) and the Andromeda Galaxy (M31). M31 and M33 form a pair that is rather separated from the MW, with their relative distance ~4 times smaller than their respective distances to the MW. Aside from their spatial proximity, which results in proximity on the sky, there are several other observational hints pointing toward the possibility that the galaxies are gravitationally bound and have interacted in the past.

### 1.1. Observational Evidence for the M33–M31 Interaction

The first morphological feature of M33 that might have been induced by the interaction is its gaseous warp. It was first found by Rogstad et al. (1976) and later, several works (Corbelli & Schneider 1997; Putman et al. 2009, hereafter P09; Corbelli et al. 2014, hereafter C14; Kam et al. 2017) confirmed that the H I disk of M33 extends farther than its stellar component and is strongly warped in the outer parts. This warp results in the continuous twist of the position angle (C14; Kam et al. 2017) and has a symmetric geometry, characteristic of tidally induced structures. However, the geometry of the warp of M33 is a little peculiar, because it is in the shape of the letter S, while the spiral arms, both in the stellar and the gaseous disk, have the reverse chirality, i.e., a Z-like shape. Initially, Rogstad et al. (1976) proposed that the origin of the warp could have been a primordial distortion of the disk, which has survived to

the present because the tilted rings would precess at different rates at different radii (Kahn & Woltjer 1959; Hunter & Toomre 1969). This scenario appeared more attractive than the interaction with M31, due to the angular separation between the galaxies, which was believed to be too large to produce sufficient tidal forces. However, a more recent study by P09 found that the gaseous features were very probably induced tidally by M31 in the past 1–3 Gyr. This finding was based on the orbit analysis that was constrained by the measurements of the radial velocities of M33 and M31 as well as the proper motion of M33 derived by Brunthaler et al. (2005).

While the vertical distortion of the gaseous disk of M33 has been known since the 1970s, the stellar disk was until very recently believed to be unperturbed in this dimension. McConnachie et al. (2009) reported finding a stellar structure extending northwest and south from the disk of M33. This feature stretches up to three times farther than the size of the disk and has an orientation similar to the H I warp. McConnachie et al. (2010) and Lewis et al. (2013) confirmed the alignment of the orientations; however, Lewis et al. (2013) also pointed out an offset between the two structures. The lack of a precise overlap of both features has been suggested to be due to the fact that the stellar component is only affected by tidal forces while the gas could also have experienced shocking and ram pressure stripping (RPS).

Both the gaseous warp and the stellar-stream-like distortion are features that strongly suggest some sort of past tidal interaction. Although they are not commonly found in spiral galaxies, M33 possesses another, more common morphological trait that may have been induced by tides, namely the dominant grand-design spiral structure. It is generally believed that M33

<footnote n="1">¹Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland</footnote>

<footnote n="2">²Institut UTINAM, CNRS UMR6213, Université Bourgogne Franche-Comté, OSU THETA, Observatoire de Besançon, BP 1615, F-25010 Besançon Cedex, France
³Aix-Marseille Université, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, F-13388 Marseille, France</footnote>

<footnote n="4">⁴Department of Astronomy, University of Wisconsin, 2535 Sterling Hall, 475 N. Charter Street, Madison, WI 53076, USA</footnote>

<footnote n="5">⁵Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="6">⁶Vilas Associate Professor.</footnote>
is a multi-armed spiral with a clear two-armed structure in the inner stellar disk (e.g., Humphreys & Sandage 1980; Newton 1980; Considere & Athanassoula 1988; Puerari 1993). It has been known for a long time that a grand-design spiral structure may be induced by tidal interactions. This view is supported by plenty of observational evidence on galaxies like M51, where one of its two spiral arms is pointing toward the companion flying by (Gunthardt et al. 2006). The scenario of triggering the two-armed spiral structure was also confirmed in several numerical simulations. It was shown that spiral arms in a disky galaxy may be induced either by a smaller companion passing by (e.g., Dobbs et al. 2010; Oh et al. 2015; D'Onghia et al. 2016; Pettitt et al. 2016) or during the orbital motion around a bigger perturber, even of the size of a galaxy cluster (Byrd & Valtonen 1990; Semczuk et al. 2017). In these simulations, the arms are triggered during the pericenter passages, and later, they dissolve and wind up over about 1–2 Gyr.

Besides two-fold spiral arms, the stellar disk of M33 is also known to host a small bar (Corbelli & Walterbros 2007; Hernández-López et al. 2009). While bars can be easily generated by global instabilities in cold disks (for a review on bar dynamics, see, e.g., Athanassoula 2013), it was shown that tidal interactions can also induce bar structures, e.g., in smaller companions perturbed by more massive hosts (Łokas et al. 2014, 2016; Gajda et al. 2017).

In addition to the morphological features present in both the gaseous and the stellar disk of M33, a gaseous bridge-like structure connecting M33 and M31 has been observed by several authors. First, Braun & Thilker (2004) reported finding an H I stream that seems to join M31 to M33. The presence of this structure was later confirmed by Lockman et al. (2012). Wolfe et al. (2013) found that about 50% of this structure is made of distinct clouds, while the rest is a diffuse component. While some authors, like Bekki (2008), argued that the H I bridge may have originated from an interaction between M33 and M31, Wolfe et al. (2013) objected, stating that the collapse timescale of these clouds (~400 Myr) is much shorter than the time since the hypothetical interaction (~1–3 Gyr). More recently, observations of Wolfe et al. (2016) have cast even more doubt on the presence of the bridge by confirming that the majority of H I material in that region takes the form of discrete clouds.

The last argument supporting the interaction scenario arises from the analysis of the star formation histories (SFHs) of both galaxies. Bernard et al. (2012) obtained the SFH of M31 in the outer disk and reprocessed the fields studied by Barker et al. (2011) in M33. They found that a rapid increase in star formation rates (SFRs) took place in both M33 and M31 about 2 Gyr ago. Bernard et al. (2012) concluded that these bursts could be triggered by a close passage of the galaxies and that a similar increase is consistent with state-of-the-art simulations of galaxy interactions and mergers.

## 1.2. Previous Works

Despite recent progress in high-resolution hydrodynamical and N-body simulations of galaxies, not many authors have attempted to reproduce particular observed interacting systems. Among the most popular objects are M51 (e.g., Salo & Laurikainen 2000; Theis & Spinneker 2003; Dobbs et al. 2010), the Antennae galaxies (e.g., Teyssier et al. 2010; Renaud et al. 2015), the Cartwheel galaxy (e.g., Horellou &

Combes 2001; Renaud et al. 2018), and the Magellanic Clouds (Besla et al. 2012; D'Onghia & Fox 2016; Pardy et al. 2018). Interestingly, M33 was not a popular target, in spite of (or perhaps due to) its proximity and the wealth of observational data.

Inspired by the discovery of Braun & Thilker (2004), Bekki (2008) carried out simple test-particle simulations to verify whether the discovered bridge-like structure could possibly have originated from the interaction between M33 and M31. The outcome of this numerical experiment was in favor of this scenario, and he proposed that the interaction might have happened 4–8 Gyr ago. He also suggested that the H I warp of M33 might be fossil evidence of such a past interaction; however, the resolution of his simulation was not sufficient to study the detailed structure of the M33 disk.

The second attempt at modeling the M33–M31 interaction was motivated by the discovery of the extended stellar stream by McConnachie et al. (2009). In the very same paper, McConnachie et al. (2009) presented results of high-resolution N-body simulations of M33 passing near M31. Both galaxies were modeled as exponential disks with bulges, embedded in dark matter halos. The relative orbit had a pericenter of 53 kpc and tidal forces excited a warp in the M33 disk that wound up and in projection closely resembled the observed distortion.

The simulations of Bekki (2008) and McConnachie et al. (2009) were only constrained by the 3D velocity vector of M33, since the radial and transverse velocities were known only for this galaxy. Brunthaler et al. (2005) obtained proper motions of M33 by water maser observations. For M31, van der Marel & Guhathakurta (2008) estimated the global transverse velocity by analyzing the line-of-sight (LOS) kinematics of its satellites. This estimate, however, was not used by McConnachie et al. (2009) to constrain the mutual orbit. Later, Sohn et al. (2012) used long-time *Hubble Space Telescope (HST)* observations of three fields in M31 to obtain its proper motions. Their measurements were corrected by van der Marel et al. (2012b) for the internal kinematics of M31. The obtained transverse velocity of M31 was found to be $17 \pm 17 \text{ km s}^{-1}$ and implied that M31 will merge with the MW in the future. This measurement also had implications for the possible past orbit between M31 and M33. Shaya & Tully (2013) used numerical action methods and integration backward in time to find the orbital history of the galaxies in the LG. They found that the results consistent with the measurements of the proper motions of M33 and M31 suggest that M33 is now at its closest approach to M31.

Recently, Patel et al. (2017b) used backward orbit integration and showed that assuming the transverse velocity of van der Marel et al. (2012b) yields orbits with an unlikely recent (~2 Gyr ago) and close (<100 kpc) pericenter passage. However, the robustness of the measurements by van der Marel et al. (2012b; and specifically the corrections for the internal motions of M31) was recently questioned by Salomon et al. (2016). Since the astrophysical implications of the measurements of van der Marel et al. (2012b) were enormous, Salomon et al. (2016) separately estimated the velocity vector of M31 by modeling the galaxy and its satellites as a system with cosmologically motivated velocity dispersion and density profiles. The resulting radial velocity was consistent with the observed one, while the tangential component was surprisingly high, ~149 km s⁻¹. This conflicting result created the problem of which transverse
velocity to adopt for M31 if one wanted to model the orbital history of galaxies in the LG.

### 1.3. This Study

In this paper, we aim to show that the observationally constrained structural parameters of M33 and M31 combined with a relative orbit that is consistent with measurements of the observed velocities (one or the other in the case of M31) can reproduce the following traits of the interaction found in M33: the gaseous warp, stellar stream, two-armed spiral structure, and an increase in SFH. In order to reach this goal, we used high-resolution N-body/hydrodynamical simulations. The model presented here does not attempt to reproduce exactly the history of the interaction between both galaxies; rather, our aim was to show that combining observables with numerical methods can result in structures similar to the observed ones, which in the past were often assigned to the interaction scenario.

This paper is organized as follows. In Section 2, we present the orbit integration method that helped us to answer which of the measurements of the transverse velocity of M31 favors more the interaction scenario. In Section 3, we give the details of the simulations (initial conditions, numerical methods, and the adopted orbit) that were carried out in order to reproduce the observed M33. Section 4 describes the properties of the simulated galaxy and its similarities to the observed M33. Section 5 provides the discussion of our results, and Section 6 summarizes them.

## 2. Orbit Integration

To construct a model of the interaction between M33 and M31, one requires mass models of both galaxies and their mutual orbit. The orbit is constrained by the final relative position and velocity of the galaxies that can be derived from the observed sky coordinates, distances, LOS velocities, and proper motions. The positions on the sky of both galaxies are known very precisely. Measurements of distances and LOS velocities have been carried out by different authors (e.g., van der Marel & Guhathakurta 2008; Gieren et al. 2013 and references therein), and while errors can be smaller or bigger, they converge to similar values. The proper motions of M33 were only measured once by Brunthaler et al. (2005). The proper motions of M31 were derived by several authors, and some of the results are in conflict, which creates problems for attempts to model the relative orbit of the M31–M33 system.

The two most recent and conflicting results for the proper motions of M31 were obtained by van der Marel et al. (2012b) and Salomon et al. (2016). In order to verify which of the measurements better favors the scenario with a recent pericenter passage, we performed a semi-analytic orbit integration similar to the one presented, e.g., in Patel et al. (2017b). The adopted values of the relative 3D positions and velocities between the two galaxies are described in Appendix A. Computational details of the orbit integration and parameters we used are briefly outlined in Appendix B.

Our approach to finding which one of the two transverse velocity estimates favors more the interaction scenario consisted of three steps. In the first step, we integrated the relative orbit of the M33–M31 system backward in time for 5 Gyr, starting from the central values of the two sets of phase-space coordinates given by Equations (2) and (3) (hereafter

called the vdM12 set) and then those given by Equations (2) and (4) (hereafter called the S16 set). The orbit for the vdM12 set had no pericenter in this time period, and the current position was the closest approach. The orbit for the S16 set had a pericenter passage about 2 Gyr ago, with a pericentric distance >100 kpc. Neither of these orbits was satisfying, i.e., neither had the pericenter <100 kpc to presumably reproduce the observed morphology of M33, so we performed additional orbit shooting.

In the second step, we integrated multiple orbits forward in time, starting from the modified values obtained in the first step for the vdM12 and S16 sets. The modification was made by varying the magnitude and the direction of the initial velocity vector. By changing these parameters (the magnification factor and the rotation angle), we obtained a grid of possible orbits. From the calculated sample of orbits, we selected those that had a pericenter passage closer than 100 kpc and checked how close they lie to the sets vdM12 and S16 after the pericenter passage. The proximity was defined in terms of the $\chi^2$ statistics, and lower values were found for the S16 set. We used the orbit found in this way in the following simulations; however, it turned out to be slightly divergent from the one obtained by the orbit integration. The differences appeared after including hydrodynamics and increasing the resolution of our simulations. The final orbit used in our fiducial model is then a result of many iterative corrections and is discussed in greater detail in Section 3.2.

One may argue that the method of selecting the orbit that favors the interaction scenario used in the second step is incomplete since there is no guarantee that the range of initial magnitudes and directions of the velocity vector can reproduce all possible values of the observed positions and velocities allowed by the observational errors. Because of that, in the third step, we tested our calculations with a simpler and more straightforward method and discussed its results instead of the statistics of $\chi^2$ obtained in the second step. The method is similar to the one presented in P09, and it consists of integrating orbits backward in time starting from values randomly selected from the range allowed by the vdM12 and S16 sets. We randomly selected 15,000 positions and velocity vectors for both vdM12 and S16 from the range enclosed by the $1\sigma$ error bars and integrated orbits backward in time for 5 Gyr.

Figure 1 shows the distributions of the pericenter distances and lookback times when the pericenters took place for both initial sets of values. What we call the pericentric distances here are in fact just distances at the closest approach. We did not integrate those orbits forward in time, or further backward, to really tell if it is a true pericenter or if the pericenter is about to take place (or took place more than 5 Gyr ago). The analysis of Figure 1 and especially of the horizontal histograms confirms our first findings that the S16 set favors the interaction scenario more strongly than vdM12. Around 20% of the orbits starting from vdM12 values are now at the closest approach, while the majority of the rest had their closest approach more than 2 Gyr ago with only one point having a pericenter recently (<2 Gyr ago) and closer than 100 kpc. Only around 0.5% of the orbits starting from S16 are now at the closest approach, while 82% have had a recent pericenter, i.e., less than 2 Gyr ago but not at the present time. Of all pericentric distances for S16, 65% are smaller than 100 kpc. Therefore, we conclude that the measurements of proper motions of M31 made by
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distributions of pericenter distances and lookback times for orbits obtained via semi-analytic orbit integration backward in time. The left panel shows the distribution for 15,000 orbits that started from randomly selected values from within the error bars of the relative position and velocity calculated by adopting the velocity of M31 derived by van der Marel et al. (2012b; Equations (1) and (2)). The right panel shows the same results but after adopting values for M31 derived by Salomon et al. (2016; Equations (1) and (3)). The pink square in both plots indicates the area of values ($r_{\text{peri}} < 100$ kpc and $t_{\text{peri}} < 2$ Gyr), suggested by Putman et al. (2009) and McConnachie et al. (2009), that would reproduce the observed morphological features of M33.</figcaption></figure>

Salomon et al. (2016) favor the interaction scenario more than do the results of van der Marel et al. (2012b). In the rest of this paper, we will aim to reproduce in simulations the S16 set of the relative position and velocity between M33 and M31.

This finding is not very surprising, since the error estimates derived by Salomon et al. (2016) are about twice as large as those of van der Marel et al. (2012b) and naturally allow for more orbital configurations. Our findings are also partially in agreement with the results of Patel et al. (2017a, 2017b), where one of the conclusions was that the interaction scenario is not very plausible once one adopts the orbital values obtained by van der Marel et al. (2012b).

# 3. The Simulations

## 3.1. Initial Conditions for the Two Individual Galaxies

Our aim when creating models of both galaxies is to reproduce their observationally derived rotation curves (C14 and Corbelli et al. 2010). This task is not straightforward since the interaction between the two galaxies changes their structural parameters. This paper is devoted to the investigation of the changes that M33 might have undergone due to such a hypothetical scenario, and because of that, we approximate M31 only as a Navarro–Frenk–White (NFW) dark matter halo with the virial mass $M_{\text{M31}} = 2 \times 10^{12} M_{\odot}$ and concentration $c_{\text{M31}} = 28$. These parameters were chosen to reproduce the rotation curve derived by Corbelli et al. (2010) only with the dark matter component (see the upper panel of Figure 2). The baryonic content of M31 contributes more to the rotation curve than the dark matter in the inner 10–15 kpc (see Figure 14 in Corbelli et al. 2010), and therefore our adopted parameters exceed the values of the parameters of the halo obtained by Corbelli et al. (2010), where the disk and the bulge were included in the modeling of the mass distribution. Our model of M31 consisted of $2 \times 10^5$ particles.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Upper panel: the rotation curve of the initial model of M31. The observed rotation curve (Corbelli et al. 2010) is represented by black dots. Lower panel: the rotation curve of the initial model of M33, after evolution in isolation for 10 Gyr without star formation, feedback, and cooling. The observed rotation curve (Corbelli et al. 2014) is represented by black dots.</figcaption></figure>
Including the disk of M31 in the simulation would be an interesting possibility given the plethora of substructure surrounding this galaxy (e.g., Lewis et al. 2013). However, it would also make it more difficult to match the inclinations of the two galaxies at the same time. The physical implication of not including the disk of M31 is the spherical symmetry of its potential. The axisymmetric disk component would slightly bend the relative orbit between M33 and M31 out of the initial plane, and this effect would have to be taken into account during iterative corrections of the orbit.

Our model of M33 is more detailed and consists of three components: a dark matter halo and stellar and gaseous disks. The initial dark matter halo had an NFW profile with $M_{M33} = 5.2 \times 10^{11} M_{\odot}$ and a concentration $c_{M33} = 11$. We increased the value of the halo mass with respect to what was estimated by C14 to take into account the mass loss that will occur during the tidal stripping in the simulations.

The best estimate of the stellar mass of M33 provided by C14 is $4.8 \times 10^9 M_{\odot}$. The combined gas mass of H, H$_2$, and helium is estimated to be $2.43 \times 10^9 M_{\odot}$ (C14). These values summed up give a total baryonic mass of $7.23 \times 10^9 M_{\odot}$, with the stars contributing 66% and the gas 34%. The surface density profiles of stars and baryons in total provided by C14 have a break at ~10 kpc. Because of this break, fitting one exponential stellar disk combined with the halo parameters and the gaseous disk will not reproduce the given rotation curve to a satisfactory degree. One must divide the stellar disk into two components with different slopes in order to have enough rotation in both the central and outer parts. In our simulations, however, we decided to start with the simplest model and keep the number of parameters of the problem as small as possible. In order to preserve this simplicity and still reproduce the rotation curve in the inner parts (and both inner and outer ones at the end of simulation, after tidal stripping of the halo), we adopt the initial baryonic disk mass to be $M_B = 8.49 \times 10^9 M_{\odot}$ with a scale length $R_B = 2.5$ kpc and scale height $z_B = R_B/5 = 0.5$ kpc.

The initial conditions for the described disk, which consisted of $10^6$ particles, were generated using the procedures provided by Widrow & Dubinski (2005) and Widrow et al. (2008). This method of generating initial conditions is suited for stellar disks but not gaseous ones. We therefore divided this initial baryonic disk into stellar (35%) and gaseous (65%) particles and evolved it in isolation so that the gaseous disk stabilized (similarly to Fouquet et al. 2017). The simulation in isolation was carried out with the GADGET-2 N-body/smoothed particle hydrodynamics (SPH) code (Springel et al. 2001; Springel 2005). The galaxy was evolved for 10 Gyr with the initial temperature set to 5000 K. The subgrid processes were turned off, because the purpose of this run was to avoid gas instabilities and produce initial conditions for the fiducial simulations, rather than to study the evolution of the galaxy in isolation. During these 10 Gyr, the gaseous disk thickened and stabilized. It became pressure-supported instead of being supported by velocity dispersion, as the stellar one is. The initial stellar and gas fractions (35% and 65%) that were kept constant in this run were adjusted to reach the values inferred from C14 (66% and 34%) after the evolution with star formation during the fiducial simulation. This tuning was done via a trial-and-error method by running simulations with subgrid physics and iteratively checking which initial values lead to the final values close to the desired ones.

<table><caption>**Table 1**
Parameters of Initial Conditions for Simulated Galaxies</caption><thead><tr><th>Components</th><th>Properties</th><th>Values</th></tr></thead><tbody><tr><td rowspan="3">M33 dark matter halo</td><td>Virial mass</td><td>5.2 × 10<sup>11</sup> M<sub>☉</sub></td></tr><tr><td>Concentration</td><td>11</td></tr><tr><td>Number of particles</td><td>10<sup>6</sup></td></tr><tr><td rowspan="3">M33 baryonic disk</td><td>Mass</td><td>8.49 × 10<sup>9</sup> M<sub>☉</sub></td></tr><tr><td>Scale length</td><td>2.5 kpc</td></tr><tr><td>Scale height</td><td>0.5 kpc</td></tr><tr><td>M33 gaseous disk</td><td>Number of particles</td><td>6.5 × 10<sup>5</sup></td></tr><tr><td>M33 stellar disk</td><td>Number of particles</td><td>3.5 × 10<sup>5</sup></td></tr><tr><td rowspan="3">M31 dark matter halo</td><td>Virial mass</td><td>2.0 × 10<sup>12</sup> M<sub>☉</sub></td></tr><tr><td>Concentration</td><td>28</td></tr><tr><td>Number of particles</td><td>2 × 10<sup>5</sup></td></tr></tbody></table>

The initial surface density profile of the gas follows an exponential distribution, with the same slope as the stellar disk. In general, the observed gas disk profiles are flatter in the center and do not follow an exponential law, which is also the case for M33 (C14). This discrepancy is an artifact of the simple method of generating the initial conditions that we used; nevertheless, it has little influence on the model, since as shown in Section 4.1 the gas profile in the center will be significantly lowered and will mimic the observed one well, due to the conversion of some of the gas into stars.

The initial rotation curve plotted in the lower panel of Figure 2 seems to match the data in the inner parts, while in the outer parts it overproduces the rotation. Note, however, that these values will be lowered as a result of tidal stripping of the dark matter halo. All of the parameters of the initial conditions for both galaxies are summarized in Table 1.

## 3.2. The Orbit

The orbits derived by the orbit integration scheme described in Section 2 happen to diverge from the orbits in the N-body simulations with the same initial parameters and mass profiles. This divergence is due to the fact that the orbit integration does not include tidal stripping and the resulting mass loss. In addition, the orbits found in N-body simulations diverge from the orbits in SPH simulations, and those differ among each other for runs with different particle numbers. All of these differences are of the order of a few tens of kiloparsecs (up to ~30 kpc), which makes it more difficult to accurately fit the adopted observed relative position and velocity (S16). It also changes the pericentric distance (by ~10 kpc), which can affect the tidally induced morphological features. Because of these numerical uncertainties, the orbit used in our fiducial simulation is not exactly the same as the one derived in Section 2, but instead it is a result of multiple iterative corrections that were made to best fit the morphological features (mostly the gaseous warp) and to end up as close as possible to the S16 set.

We applied a similar strategy to find the initial rotational angular momentum (spin) of our model of M33. The direction of this vector is the same as the vector normal to the disk's plane, which can be derived from the observed inclination and the position angle of M33. The first runs of our simulations were performed with the observationally derived values; we found, however, that this vector was changing during the orbital evolution. Therefore, we iteratively corrected it, so that at the present time it reproduces reasonably well the observed
<table><caption>**Table 2**
Initial Orbital Parameters and Inclination of M33</caption><thead><tr><th>Vector</th><th>X</th><th>Y</th><th>Z</th><th>Unit</th></tr></thead><tbody><tr><td>Position</td><td>-278.1</td><td>326.7</td><td>44.5</td><td>kpc</td></tr><tr><td>Velocity</td><td>56.5</td><td>-6.5</td><td>32.8</td><td>km s<sup>-1</sup></td></tr><tr><td>Normalized spin</td><td>0.048</td><td>0.998</td><td>-0.038</td><td>...</td></tr></tbody></table>

position angle and inclination as well as the observed morphology (we focused mostly on the gaseous warp).

The initial orientation of the disk of M33 with respect to its orbit around M31 can be parametrized by two angles: α, the angle between the spin of the disk and the orbital angular momentum, and β, the angle between the spin and the direction of the velocity of M33 on its orbit. We started with β derived from the normal vector from observations and varied it by ±45° to see whether both directions in the parameter space made the gaseous warp flatter. Hence, we assumed that for this parameter we are approximately in the best place. For α, we tried eight different values from the range between ~15° and ~80°. Values close to the prograde case, α = 0°, keep the tidal disturbance within the two dimensions of the plane of the disk, and the warp-like distortions are not present. Crossing α = 90° and approaching the retrograde case, α = 180°, decrease the effect of the tidal perturbation of the disk and make it hard to generate spiral arms (D'Onghia et al. 2009, 2010; Łokas et al. 2015). The final value of α turned out to lie only a few degrees away from what was derived from the observations, and we found that as long as it is not too close to 0°, it has very little influence on the 3D appearance of the galaxy. We also found that rather than changing the inclination, other parameters were much more influential in shaping the warp, namely the magnitude of the tidal perturbation (parametrized either by the pericenter distance or by the mass profile of M31) and the inclusion of the hot gas halo of M31, which is discussed later.

The values of the initial position, velocity, and spin that we used in our fiducial simulation are summarized in Table 2. M31 was placed at the center of the coordinate system, and its initial velocity was (0, 0, 0) km s⁻¹. The upper panel of Figure 3 shows how the relative distance and velocity between M33 and M31 was changing in time for the adopted orbit. The orbit has a pericenter of 37 kpc 2.7 Gyr after the beginning of the simulation. The pericenter distance is larger than the sum of the sizes of the disks of both galaxies, hence the approximation of the M31 potential as a pure dark matter halo does not imply that we ignore additional effects resulting from the crossing of the disks. We note that dynamical friction decreased the apocenter of the orbit by a factor of roughly 2.8, from 432 to 151 kpc after the pericenter passage.

The lower panel of Figure 3 shows the shape of the relative orbit in its plane. It is worth noting that the position of M31 was changing significantly due to the attraction from M33 (~100 kpc), which is not so surprising given that the adopted masses result in a mass ratio of M33 to M31 as large as ~26%. The orbit presented here was centered on M31 at each time step to make the image more clear.

### 3.3. The Code and the Fiducial Run

The fiducial simulation presented in this paper was carried out with a modified version of the GADGET-2 N-body/SPH code

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Upper panel: the time dependence of the relative distance *d* and velocity *v* between M33 and M31 for the adopted orbit. Lower panel: the shape of the relative orbit between M33 and M31 in the orbital plane. At each time step, it was centered on M31.</figcaption></figure>

(Springel et al. 2001; Springel 2005) that includes star formation and feedback and cooling processes added as described in Hammer et al. (2010) and Wang et al. (2012). The subgrid physics was implemented according to the recipes given by Cox et al. (2006). We chose to adopt the values of the parameters as advised by Cox et al. (2006), namely a star formation efficiency of 0.03, a feedback index of 2, a density threshold for star formation of 0.0171 $M_{\odot}$ pc⁻³, and a timescale of feedback thermalization of 8.3 Myr. Instead of tuning the subgrid parameters listed above, we were changing the initial gas fraction to match the observed one at the present time. If one intended to reproduce the simulations presented here with a code using a different subgrid model, this value would have to be tuned again to conform with the different subgrid model, as it is generally believed that subgrid schemes can change the outcome of simulations very significantly (e.g., Arthur et al. 2017).

The implementation of the subgrid physics using the prescriptions of Cox et al. (2006) includes cooling processes described by the procedures given by Katz et al. (1996). These procedures result in assigning every gas particle the neutral hydrogen mass fraction. We stress here that throughout this paper, we will use the name “gas” for gas particles in general while the term “neutral hydrogen” or $H I + H_2$ will be used for
masses of gas particles reduced by this fraction given by the cooling procedures in order to match better the observations that mostly probe the H I content of M33.

We followed the evolution of the system for 4.5 Gyr with outputs saved every 0.05 Gyr. We adopted the softening lengths of 0.09 kpc for the stellar and gaseous particles, 0.63 kpc for the dark matter of M33, and 10 kpc for the dark matter of M31.

### 3.4. The RPS Experiment

As described below, we managed to reproduce several features of M33 in the fiducial simulations. We found however that one of these features, namely the gaseous warp, cannot be reproduced to a satisfying degree of similarity (especially at larger radii) by tidal effects alone. In order to argue that the RPS originating from the hot gas halo of M31 would significantly improve the resemblance, we performed simulations with lower resolution that included the hot gas halo of M31. The existence of the circumgalactic medium around M31 was demonstrated by Lehner et al. (2015); however, there are no strong observational constraints on the mass of the hot gas. Other galaxies are believed to have hot gas halo masses roughly 100 times smaller than their dark matter masses (Miller & Bregman 2015). We started the modeling of the hot gas halo with the mass value close to this estimate; however, the final adopted mass was smaller than this, as a result of tuning the effect of RPS to obtain the most similar gaseous warp.

In the simulation including the RPS, M31 was modeled as a Hernquist (Hernquist 1990) dark matter and hot gas halo. Parameters of the dark matter component were fitted to reproduce the NFW profile discussed in Section 3.1 and took the values of $2 \times 10^{12} M_{\odot}$ for the mass and 41 kpc for the scale radius. The hot gas halo had the same scale radius and the mass of $3 \times 10^9 M_{\odot}$. The dark matter halo consisted of $2 \times 10^5$ particles with the softening length of 10 kpc, and the gas halo was made of 70,630 particles with the softening of 0.2 kpc. Initial conditions for M31 for this simulation were generated using the CLUSTEP code[^7] (Ruggiero & Lima Neto 2017). Initial conditions for M33 were generated in the same manner as described in Section 3.1 with very similar parameters, except for the particle numbers and softening lengths. Particle numbers were $2 \times 10^5$ for dark matter, 126,000 for the gas, and 74,000 for stars, while softenings were 1.4 kpc, 0.2 kpc, and 0.2 kpc, respectively.

## 4. Properties of the Simulated Galaxy

The initial model of the M33 galaxy was transformed during the simulated time by both the subgrid physics (star formation, cooling, and feedback) and tidal interactions with the M31-like dark matter halo. The most important morphological features have been induced by the tidal interaction with M31 during the pericenter passage. To verify this, we carried out simulations with the same initial conditions for M33 but in isolation for 4.5 Gyr. A short discussion of this case is included in the appendix.

### 4.1. General Properties

Star formation converts $2.61 \times 10^9 M_{\odot}$ of gas mass into stars (hereafter called young stars, i.e., those stars that were

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Surface density profiles for stars, neutral hydrogen, and baryons in total in our model of M33 at the time of the best match, compared with density profiles for stars, H I, and baryons from C14.</figcaption></figure>

born after the beginning of the simulation, unlike the old stars, which were present in the initial conditions), which resulted in a total stellar mass of $M_* = 5.59 \times 10^9 M_{\odot}$ (at the time of the best match; the selection of this time is explained in Section 4.2). This mass contributes 66% of the baryonic mass, which is the ratio we were aiming to obtain, consistent with the findings of C14.

The star formation created a central plateau in the surface density profile of the gas, which is also seen in observations (C14). It is clearly visible in Figure 4, where we compare the obtained surface density profiles for stars, hydrogen, and baryons in total with the observationally measured profiles given by C14. We find that the simulated profiles are in reasonably good agreement with the observed ones. We reproduced the previously mentioned gaseous “core” and also approximated the slopes of the stellar and total baryonic distributions.

Tidal stripping does not strongly affect the baryonic mass of M33; however, a significant fraction of the dark matter mass was stripped. We fitted the NFW profile to the dark matter density distribution at the time of the best match and obtained a mass of $M_{vir} = 4 \times 10^{11} M_{\odot}$ and concentration $c = 10.5$. These values are within $1\sigma$ from the best-fit model of C14, and the mass is ~77% of the initial halo mass.

The rotation curve of M33 at the time of the best match is presented in Figure 5. This curve was calculated directly from the sum of forces from the particles in simulations. We note that the initial curve presented in Figure 2 was calculated using analytical formulae and was plotted only to show what we aim to reproduce. The obtained final rotation curve exceeds the observed data by tens of km s⁻¹ in the inner ~3 kpc. This is due to the bar formation in this region, which affects the potential of the galaxy. In the outer parts, the simulated curve falls slightly under the observed one, which is a result of the tidal stripping of dark matter. However, despite these discrepancies on both ends, we find that in general the

<footnote>[^7]: https://github.com/ruggiero/clustep</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Rotation curve of our model of M33 at the time of the best match compared with the observed rotation curve provided by C14 (black dots).</figcaption></figure>

simulated curve is a reasonable approximation of M33, being flat at around 100 km s⁻¹ throughout the majority of radii. The resemblance of the curve to the observed one is better when the bar is excluded from the analysis and the disk component is calculated via the thin disk approximation from the fit to the surface density profile that excludes the steep inner part originating from the bar.

The tidal interaction during the pericenter passage induced spiral arms and a warp in both the gaseous and the stellar disk. These features survived for 0.75 Gyr after the pericenter and are present at the time of the best match (see Figure 6). The two-armed, grand-design spiral structure characteristic of tidal encounters is present in the gas, as well as the old and young stars. The off-planar tidal distortion (the warp) is mostly visible in the gaseous disk and the old stellar disk. This is due to the fact that the new stars are born in the inner region of the galaxy and the disk that they create is too small to exhibit a similar disturbance. The young stellar disk forms a small bar, which is also present in the observed M33 (e.g., Corbelli & Walterbros 2007; Hernández-López et al. 2009). However, as discussed in Section 4.5, the bar is not a tidal feature, and it also forms in isolation. Rather, the origin of the bar is related to the fact that the young stellar disk formed from the gas by star formation is not stable against bar formation. As shown in Athanassoula et al. (2013), this process can be weakened by increasing the gas fraction of the galaxy; however, in our case, the gas fraction is constrained by observations and subgrid physics.

## 4.2. The Warped Gaseous Disk

Figure 6 clearly demonstrates that both the gaseous and stellar disk have been warped by the tidal interaction. McConnachie et al. (2009) showed using N-body simulations that exciting such a distortion with similar parameters of M33 and M31 is possible. Their model, however, was not constrained by the proper motions of M31 and did not include gas physics. In this subsection, we briefly describe how we found the best orbital epoch and the viewing position that would result in the geometry of the warped disk similar to the observed one, and we compare our simulated images with those obtained by radio observations.

In order to find the projection and the orbital epoch that will result in a neutral hydrogen map resembling best the observed one, we visually inspected images of the rotated gas disk of

M33. The rotations were constrained by the demand for the relative 3D position and velocity of M31 with respect to M33 to be consistent within 2σ with the values of adopted Equations (2) and (4) (reversed, because Equations (2) and (4) are the relative coordinates of M33 with respect to M31 and here we look at the system from the perspective of M33). Using this procedure, we found the time of the best match to be 0.75 Gyr after the pericenter, which is 3.45 Gyr since the beginning of the simulation and just 0.05 Gyr before the apocenter. The relative position and velocity of M31 in the reference frame corresponding to the best-match projection at this time are $X_{\text{rel,M31}} = (56.1, 91.5, 106.2)$ kpc and $V_{\text{rel,M31}} = (7.1, -38.8, 41.1)$ km s⁻¹. One out of six of these coordinates is 1.75σ away from the adopted values. Another one lies within 1.25σ and the remaining four are within 1σ.

The neutral hydrogen density map for the best-matching time and projection from the fiducial simulation is presented in the middle panel of Figure 7 (coordinates (ξ, η) are distances in R.A. and decl. from the center of M33). The upper panel of the same figure shows the H I map published by P09. By comparing both images, we find that tidal interactions managed to only reproduce the inner warp of the gas disk, i.e., features marked 1 and 2 in red in both maps. At larger radii in the simulated image, strong spiral/tidal extensions are the most eye-catching feature, while in observations this is not the case. We identify several substructures in the image of P09 that slightly resemble the spiral features in the simulated image (marked 3–6); however, in observations they are merely composed of one or two differently shaped contours, while in the image from simulations they are the dominating signal at these radii. Another discrepancy is that in the image from the fiducial simulation, the gas structure is more or less symmetric, while in the P09 image there seems to be more gas northwest from the galaxy. Also, the northern part of the warp, marked 1, is longer than its southern counterpart, 2. This asymmetry in the observed gas is pointing toward M31 on the sky, and this coincidence motivated us to investigate whether RPS from the hot gas halo of M31 may improve the simulated image.

The best-matching map from the simulations with RPS (found in the same way as for the fiducial simulation) is presented in the lower panel of Figure 7. Strong spiral/tidal extensions present in the fiducial model (3 and 6) were adequately weakened by RPS, and the dominating signal is coming from the S-oriented warp. The resemblance to the observations is still far from perfect (i.e., the galaxy is too elongated in simulations in comparison with the more flattened distribution in observations), but qualitatively, the S-shaped structure is similar to the observed one, with the inner spiral arms having the reversed Z-chirality. In the RPS image, there is also too much gas at the lowest density, located southeast from the disk. We suspect that this material (and maybe other discrepancies) may arise due to the limitations of the SPH scheme used in the code to mimic the hydrodynamics of the gas. It is well known that SPH codes fail the blob test (Agertz et al. 2007; Hopkins 2015) and thus underestimate the RPS of the cold gas and its mixing with the hot component. Because of these numerical problems, we decided not to explore the RPS experiment in greater detail, since the physics of this process would not be properly modeled anyway, while the other M33 features can be reproduced simply with simulations of tidal evolution. The asymmetry of the neutral hydrogen in the RPS image at first sight seems to be reversed
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Surface density distributions seen from three different directions for the gas and —the old and young stellar particles for the simulated M33 at the time of the best match.</figcaption></figure>

with respect to the observed asymmetry that motivated us to perform that experiment; however, it is the case only for the lowest density material. If we focus on the higher density, e.g., the fourth or fifth contour, where the gas physics is better resolved, then the asymmetry is reproduced and the gas stretches up to 1° north, while in the south it only reaches ~0.6–0.7.

Corbelli & Schneider (1997), C14, and Kam et al. (2017) all reported that the gaseous warp quantitatively manifests itself in the continuous change of the position angle of the disk with increasing radius. The upper panel of Figure 8 presents the

radial dependence of the position angle obtained from images based on simulations and the comparison with the results from C14. The lines from the simulations were obtained by ellipse-fitting to column density distributions. The curves from both simulations exhibit a drop in values, of the order of 30°, with the one from the fiducial simulations being more important and resembling more the observed curve. The drop from the RPS case is not as big, about 20°, and corresponds well with the impression given by the map in Figure 7, i.e., that the overall shape of M33 is more elongated and the shift is milder than in the observations.
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Upper panel: the H I column density map of M33 published by P09 (© AAS. Reproduced with permission). Middle panel: best-matching neutral hydrogen column density map from our fiducial simulations. The isodensity contours were plotted at the same levels as in P09. The red numbers mark the tidally induced features that resemble the corresponding features found in the image of P09. Lower panel: best-matching neutral hydrogen column density map from simulations that included RPS from the hot gas halo of M31. The inclusion of RPS seems to be necessary to weaken the strong tidal arms (3–6) that are almost invisible in the observed galaxy.</figcaption></figure>

Supplementary quantitative information about the warp is given by the radial dependence of the inclination presented in the lower panel of Figure 8. The curves from the simulations show more structure than the ones obtained by C14. This is due to the fact that the extended gaseous features are less regular in

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Upper panel: the radial dependence of the position angle obtained for the best-matching neutral hydrogen density map from the fiducial simulation (red line) and the RPS experiment (blue line) in comparison with the measurements from C14 (black points with error bars). Lower panel: the radial dependence of the inclination for the same cases as in the upper panel.</figcaption></figure>

the simulations than in observations. Contrary to the visual impression, values from the simulations lie relatively close to the observed ones. This coincidence is probably just the result of a very simple method of deriving the inclination that we applied, and a bigger discrepancy would be seen if we have fitted a tilted ring model as in C14. Our models, however, are not very precise and do not require such sophisticated tools to analyze them; with the ellipse-fitting, we just wanted to show that both the inclination and the position angle have approximately similar values and radial dependence.

The warp of the M33 gaseous disk manifests itself also in the kinematics of the gas. Figure 9 presents the comparison between the kinematic map published by C14 and the one obtained from our fiducial simulation at the time and the projection of the best match. The simulated image reproduces reasonably well the observed one, and we find that the distortion of the disk is clearly visible in the twist of the zero-velocity line, as was also found by other authors (e.g., Kam et al. 2017). The image from fiducial simulations has a much thinner outline than the observed one, which is the result of the previously mentioned offset between the observed inclination and the one in the best-match projection. The best resemblance between the images is seen when one follows the lines that separate cyan from blue and red from yellow. These clearly show where the disk ends and tidal features start to be visible. For example, the tidal arm marked as 3 in Figure 7 is seen in the simulated kinematic map as a red arc southwest of the main
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Left: H I velocity map of M33 from C14 (reproduced with permission © ESO). Right: gas velocity map of M33 for the best-matching projection and time from our simulation. The simulated image has a similar resolution to the one from C14. While constructing our map, we took only those bins that contained more than five particles in order to cut off the extended tidal tails and better mimic the image of C14.</figcaption></figure>

disk. The corresponding arc is also present in the image from C14. The extended tidal arm is not clearly visible in the H I map; however, this red arc may be a hint that perhaps some sort of tidal arm is present there. We note that the image obtained from the RPS experiment possesses features similar to those of Figure 9.

### 4.3. The Stellar Stream

Unlike the gaseous warp, which has been known since the 1970s (Rogstad et al. 1976), its stellar counterpart was discovered only recently by McConnachie et al. (2009). Because this structure has been known for a much shorter time, it has been examined to a lesser extent than the warp. In Section 4.2, we selected the time and the best viewing position to match the warp. In this subsection, we briefly discuss the fact that for the same projection, the stellar stream is also present. In general, we could have repeated the whole procedure to find a stream that resembles the observed one more; however, its shape is not as well investigated as that of the warp, and it changes its orientation from C-like to S-like when stars with different metallicities are selected (McConnachie et al. 2010; Lewis et al. 2013).

Figure 10 presents the surface density distribution of all stars from our simulation. We decided to show the simplest possible image that can be obtained from the simulations and not translate it into surface brightness or apply complex processing to mimic images like Figure 13 of McConnachie et al. (2010), since performing such procedures on the simulation data would not be straightforward and require many assumptions and model-dependent methods. In this paper, we only want to quickly and qualitatively compare the resulting structure to observations.

We find that in general the simulated stream resembles the observed one. The two most evident features are in both cases the northwestern and southeastern extensions. The simulated streams span ~2°, which is similar to the observations. The

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Surface density distribution for all stars from the simulated M33 at the time of the best match and seen from the best viewing position.</figcaption></figure>

structure discussed in McConnachie et al. (2010) seems a bit more vertically elongated; however, this could also be accomplished in the simulations by choosing a different viewing position. The image from the simulation obviously has a better resolution, and some substructures that are visible here (for example, a small gap in the southeastern extension) would probably be smoothed out and not visible in observations. The surface density of the edge of the simulated disk is about 10 times bigger than the density of the stream. This agrees with the contours given by McConnachie et al. (2010), where the ratio of the surface brightness between the first and the last contour is 5.75. Apart from resembling the observed stellar structure, our simulated stream is also very similar to the simulated one presented by McConnachie et al. (2009), where the stream is made of wound up and projected tidal tails.
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Surface density maps of the neutral hydrogen (left), old stars (middle), and young stars (right) of the simulated M33 at the best-match time and projection.</figcaption></figure>

### 4.4. The Spiral Arms

As seen in Figure 6, both the gaseous and the stellar disks possess a grand-design spiral structure at the time of the best match. The two-armed spiral signal is characteristic of tidally induced spirals and is also present in the observed M33. Figure 11 shows how the spiral arms of our simulated galaxy appear in projection on the sky, seen from the best-match point of view. For the neutral hydrogen map, the similarity with the H I image of C14 is very good. The projection and the fragmented nature of the gaseous disk makes the arms look more flocculent, which is also the case for the C14 image where it is hard to identify the number of arms.

In the stellar images, the spiral arms are clearly more visible in young stars. This closely agrees with observations, since generally young stars are a better tracer of spiral structure. The spiral arms reach deeper into the young stellar disk than the old disk, a result of the fact that the young stellar disk was formed out of gas and is about twice as small as the old disk (see Figure 6). The spiral structure penetrates the young disk down to ~3–6 arcmin and is stopped there by the presence of a bar. The size of these spiral arms is similar to that seen in the B–V–I image of C14 made from the Massey et al. (2006) survey. The connection between the arms and the bar is a bit different in our simulated image compared to the one presented by C14. However, this results from the way the bar is formed in our simulations (which is not by tidal interaction; see Section 4.5), and it would be very difficult to match its position angle with the observed one.

Besides the qualitative inspection of the spiral structure, several quantitative properties may be measured for spiral arms, namely their number (i.e., multiplicity), strength, pitch angle, and pattern speed. First, we measured different Fourier modes $|A_m|$ (similarly to, e.g., Pettitt et al. 2016) for the three different components to see what their dominant spiral signal is. Their time evolution is presented in Figure 12. A grand-design ($m=2$) signal is clearly being tidally induced shortly after the pericenter passage and is stronger than other modes in all three components. Unfortunately, we do not reproduce higher modes, which are present in the observed M33 (and also in other tracers, e.g., in $H\alpha$; see Kam et al. 2015). Reproducing the multi-armed structure in M33, which is not a dominant one, would require better tuning of the initial model of M33, since the number of arms is mostly a function of a disk-to-halo mass ratio (Athanassoula et al. 1987; D’Onghia 2015). In our work, we focused mainly on tidally induced morphological features in

<figure><img src="image_18.png" /><figcaption>**Figure 12.** Time dependence of the $|A_m|$ Fourier coefficients measured in rings of 4 kpc $< R < 11$ kpc for the gas (upper panel) and old stars (middle panel) and 3 kpc $< R < 6$ kpc for young stars (lower panel). The first vertical dashed line indicates the time of the pericenter passage and the second one marks the time of the best match.</figcaption></figure>

M33, and we did not modify the model too much after we found the one that fits the rotation curve.

In order to measure the pitch angle and the strength of spiral arms, we expanded the surface distribution of gas and stellar
<figure><img src="image_19.png" /><figcaption>**Figure 13.** Upper panel: the time dependence of the pitch angle for the gas and old and young stars measured in the same rings as in Figure 12. Lower panel: the time dependence of the arm strength for the same components measured in the same rings. The first vertical dashed line indicates the time of the pericenter passage and the second one marks the time of the best match.</figcaption></figure>

particles in logarithmic spirals as discussed in, e.g., Sellwood & Athanassoula (1986). A more detailed description of the method we used can be found in Section 3.2 of Semczuk et al. (2017). The time evolution of these parameters after the pericenter passage is presented in Figure 13. The pitch angle starts from high values during the pericenter, ~40°–60° for all components. Later, it exponentially decreases, reaching 19°6 for the gas, 21°8 for old stars, and 17°9 for young stars, at the time of the best match. The value for the best tracer of the spiral structure, i.e., young stars, agrees very well with the pitch angle for the *m* = 2 structure, 16°5 obtained by Considere & Athanassoula (1988) and 17°1 found by Puerari (1993). The values for the other components also lie close to these measurements. This finding argues for the possibility that the grand-design component of the spiral structure of M33 was indeed induced tidally, and the pitch angle had enough time after the encounter to wind up to the observed value. We note that higher values for the pitch angle of the *m* > 2 structure may be found in the literature (e.g., Sandage & Humphreys (1980) and also measurements may reach ~40° for *m* = 2 (references in Considere & Athanassoula 1988). However, as described by Considere & Athanassoula (1988), these values originate from measurements in the outer, looser spiral structure of M33, and the best fit for the brightest inner two-armed structure is indeed ~16°. We conclude that the degree of winding up of the brightest component of the spiral structure in our simulations agrees very well with the brightest component in observations.

<figure><img src="image_20.png" /><figcaption>**Figure 14.** Radial dependence of the pattern speed measured between the time of the best match and 0.15 Gyr before this time.</figcaption></figure>

The lower panel of Figure 13 shows the time evolution of the arm strength. It peaks at approximately the same time as $|A_m=2|$ in Figure 12, and later it slowly decreases with time, which means that the spirals are dissolving. At the time of the best match, the highest value of the arm strength is in the young stellar component, which closely corresponds with the fact that the spiral structure is most pronounced in this component.

The last property that we measured for the spiral arms of the simulated M33 is the pattern speed. We used the method discussed, e.g., by Dobbs (2011), which relies on tracking the maximum of the surface density $\Sigma_{max}$ in polar coordinates at a given radius between two time epochs. The pattern speed can then be calculated as the difference between the angular positions of the maxima, divided by the time difference. The radial dependence of the mean pattern speed of the two-armed structure at the time of the best match is presented in Figure 14. We find that for both the gaseous and the old stellar arms, the pattern speed decreases radially and tightly follows the inner Lindblad resonance. This indicates that the spiral arms are kinematic density waves, as expected for tidally induced arms (Dobbs & Baba 2014). Estimates of the pattern speed of M33 found in the literature (e.g., 15 km s⁻¹ kpc⁻¹ by Courtes & Dubout-Crillon 1971 or 28 km s⁻¹ kpc⁻¹ by Puerari 1993) highly exceed the values obtained from our simulation. This discrepancy probably arises from the fact that those estimates were based on the quasistationary density wave theory, where the pattern speed of the arms is constant. This is not the case in our simulation, and it was only recently reported (Saha & Elmegreen 2016) that a constant pattern speed can, in fact, be obtained in an *N*-body case.

### 4.5. The Bar

The observed M33 is known to possess a small bar in the central parts of its disk (Elmegreen et al. 1992; Regan & Vogel 1994; Corbelli & Walterbros 2007; Hernández-López et al. 2009). As can be seen in Figures 6 and 11, our simulated M33 also forms a bar, mostly in the young stellar disk. In order to verify whether the origin of this bar is related to the tidal interaction or due to the secular instability of the disk, we measured the time evolution of the Fourier mode $|A_2| \equiv |A(2, p=0)|$ in the inner 3 kpc of the face-on surface distribution of all stars for the fiducial simulation and the simulation of M33 in isolation. The upper panel of Figure 15
<figure><img src="image_21.png" /><figcaption>**Figure 15.** Upper panel: the time evolution of the $|A_2|$ Fourier mode calculated within the radius of 3 kpc for our fiducial simulation and for the case in isolation. The first vertical dashed line indicates the time of the pericenter passage and the second one marks the time of the best match. Lower panel: the radial profiles of $|A_2|$ measured at the time of the best match for the fiducial simulation and at the corresponding time for the isolated case. Dashed lines indicate estimates of the bar length.</figcaption></figure>

shows that the tidal evolution has very little influence on the bar formation and even slightly suppresses it.

The fact that the bar in isolation grows bigger is also confirmed by the measurement of its length in both cases. We estimate the bar length by calculating the radial dependence of $|A_2|$ and finding where it drops below a fixed threshold (here we assume it to be 0.5 of the maximum value of the inner peak). This method yields a length of 1.56 kpc at the time of the best match for the simulation with M31 and 1.9 kpc for the isolated case at the corresponding time (lower panel of Figure 15). While this estimate (1.56 kpc) agrees quite well with one of the measurements in the literature (1.5 kpc by Elmegreen et al. 1992), it is at least two times larger than other measurements (e.g., 0.4 kpc by Regan & Vogel 1994 or 0.7 kpc by Corbelli & Walterbros 2007). The value obtained by Elmegreen et al. (1992) has been suspected of being biased by the influence of the very tight spiral arms surrounding the bar in that region (Hernández-López et al. 2009), and the agreement with it seems to be coincidental.

The bar in our simulation is formed in the young stellar disk, only because this disk is formed from the gas and is very unstable. If we started our simulation later (or earlier) the bar length would be different, since the bar would have had enough time to grow bigger (as in isolation) or would not have had enough time to grow to the present size. In this paper, we aim to investigate the very recent history of M33, although we cannot determine when exactly the processes responsible for bar formation started. Our results only confirm that the influence of tidal interactions on bars can be ambiguous, and

they do not always accelerate bar formation (Pettitt & Wadsley 2018).

## 4.6. The SFH

Bernard et al. (2012) reprocessed the data from Barker et al. (2011) for two fields in M33 at radii 9.1 and 11.6 kpc and found a peak in the SFH that occurred ~2 Gyr ago. The density of the SFR had a value of $0.6 \times 10^{-9} M_{\odot} \text{yr}^{-1} \text{pc}^{-2}$ at that time, which is approximately three times larger than the average values found for earlier times. Bernard et al. (2012) also found similar activity in the SFH of M31 and suggested that the past mutual interaction of both galaxies might have caused those synchronized peaks.

We checked whether M33 in our simulation also exhibits a rapid increase in SFR at similar radii. Figure 16 presents the time evolution of the densities of SFR measured in rings (and one circle) of sizes of 1 kpc. At the beginning of the simulation, the values of the SFR were obviously high in the central regions, where most of the star formation happened. As the star formation slowly quenches, these values decrease and stabilize. At radii 6 to 12 kpc, the SFR rapidly increases again after the pericenter passage (marked as the first vertical dashed line in Figure 16). The values are from 2 to 10 times greater than the average ones before. The values closest to those found by Bernard et al. (2012) occur at distances from 6 to 9 kpc. However, it is difficult to compare the absolute values between simulations and observations, because the SFR was measured in a different way and in different kind of fields (here in rings, in observations in two square fields). The peak of the SFR in the simulation is separated by around 0.75 Gyr from the time of the best match, a period ~2.6 times smaller than the one found by Bernard et al. (2012). When comparing these times of the peak, it is worth keeping in mind that deriving the SFH is strongly model dependent. For example, for the same fields, Barker et al. (2011) obtained the increased SFR at ~3 and ~6 Gyr ago. We conclude that the pericenter passage of ~37 kpc can induce an increase in the SFR of the same relative order and at similar radii as observed. We also found that introducing the hot gas halo of M31 does not change these results regarding SFHs at different radii.

## 4.7. The Extended Gaseous Structure

More than a decade ago, Braun & Thilker (2004) found a faint H I stream that seems to connect M31 and M33. Using test-particle simulations, Bekki (2008) showed that such a bridge-like structure may have been tidally induced by the interaction between the two galaxies. To see whether such a feature would be present in the model of interaction presented here, we used the simulations that included the RPS, since as we showed in Section 4.2, this additional physical process is crucial in reproducing better the outer regions of the cold gas of M33. We found no cold gas component stretching in a bridge-like form between M33 and M31. This finding partially agrees with the results of Wolfe et al. (2016), who showed that the majority of H I between the two galaxies is contained in small discrete clouds rather than in a bridge. In our simulations, we do not have enough resolution to model such objects.

In our fiducial simulation, we found cold gas stretching from M33 in the direction of M31 and also on the other side of M33 (similarly to the stars as seen in Figure 10). This structure formed as a result of tidal stripping; however, as discussed in
<figure><img src="image_22.png" /><figcaption>**Figure 16.** Time evolution of the SFR normalized to the area of the ring (or circle) in which it was measured. The first dashed vertical line indicates the time of the pericenter passage; the second marks the time of the best match.</figcaption></figure>

Section 4.2, we had to include RPS to weaken these features, since their signal was too strong inside 1° from the disk of M33, which is not seen in observations.

Despite not finding an extended structure of neutral hydrogen northwest from M33, in the RPS experiment, we

find low-density gas material southeast from M33 (partially visible in the lower panel of Figure 7). The position of this material corresponds well with the H I clouds recently found by Keenan et al. (2016). The LOS velocities of the material from the simulation and the clouds of Keenan et al. (2016) are consistently oriented toward the observer; however, the observed material is approaching faster than M33 itself, while in our simulations this velocity is lower than the one of M33. We did not investigate the properties of this material in greater detail, since one of the conclusions of Keenan et al. (2016) was that at the moment we cannot rule out the possibility that this gas is in fact some contamination from the Magellanic stream.

# 5. Discussion

## 5.1. Accuracy and Uniqueness of the Model

As stated in the introduction, the model presented here is not aimed to precisely reproduce the Triangulum Galaxy and its present relative position and velocity with respect to Andromeda. The main goal of this study was to show that once we adopt the values of the structural parameters that are similar to the ones derived from observations and use them in numerical simulations, the observed morphology of M33, the burst in its SFH, and the relative position and velocity with respect to M31 may be reproduced with reasonable accuracy. The purpose of this paper is to show that the observed parameters allow for the interaction with M31 to sufficiently disturb the disk of M33, and that such scenario utilizes an orbit that is consistent with measurements of distances and velocities of both galaxies.

We noted several times that discrepancies between our model and the observations still persist. The biggest discrepancy seems to be in the relative position between M31 and M33 at the best-match time and projection, specifically in the X-direction, where it lies only within $1.75\sigma$ of the adopted error bar. This corresponds to a difference of 41.1 kpc. The errors we adopted here as $1\sigma$ originated from the assumption that measurements have Gaussian error distributions, and we can easily transform them into different reference frames. These errors are not the maximal possible errors. The maximal possible errors of the relative positions are the same in all directions and equal to the sum of errors of both distances (assuming that the sky coordinates have exact values), i.e., 63 kpc. This value is greater than the discrepancy of our results. However, the adopted errors that we used seem slightly underestimated when we compare them with a rather large range of observationally derived distances to both galaxies, especially M33. For example, in the NASA/IPAC Extragalactic Database (NED), the measurements of the distance to M33 range from ~620 kpc to ~960 kpc (with the mean value of 849 kpc and the standard deviation of 212 kpc). This sample may be biased by very old measurements; nevertheless, Gieren et al. (2013) considered a sample of more recent attempts and the range of distances was still quite extended, from ~720 kpc to ~970 kpc. Comparing this range with our adopted error of 23 kpc suggests that taking such a small error may be recognized as a little too optimistic. Reviewing the results presented in the literature and their discrepancies shows that the astrophysics of the M33–M31 system is a branch of science with a distance accuracy of ~50–100 kpc, rather than ~25 kpc, as it is for the projected position errors.

This huge range of possible values that are only weakly observationally constrained gives modelers a variety of
parameters to explore. The distance to M33 is not the only parameter with such a big allowed range. For example, the mass of M31 may vary from ~7 × 10¹¹ M☉ (Evans et al. 2000) up to ~2.5 × 10¹² M☉ (according to the introduction in Patel et al. 2017a). In general, one could run simulations with different halo models of M33 and M31, gaseous and stellar disks of M33 and M31, the bulges of both galaxies included or not, the gaseous halo of M31, as well as different relative orbits and inclinations of both galaxies, particle resolutions, hydrodynamical solvers, and subgrid models and their parameters. The simplest strategy would be to perform simulations probing a multidimensional grid of parameters; however, such a task would be computationally expensive and the analysis of the simulations would be very difficult, given the abundance of observational properties to reproduce.

While constructing the model described in this paper, we applied a different strategy, given the limited resources that we had. In the first step, we found the initial orbit using an orbit integration scheme. This orbit was later tested and corrected in N-body runs. These runs were used as a basis for the following SPH simulations, which were also iteratively corrected after comparing the obtained properties of the galaxy with observations. The model presented here is a result of about 80 such iterations. Some of these intermediate steps reproduced certain observables better than the fiducial model (e.g., the peak in the SFH occurred a longer time ago); however, this was obtained at a cost of less similarity to other observables. Reproducing a particular galaxy is therefore a minimalization issue, and the fiducial model we present here seemed to be a compromise in terms of the similarities with the most important characteristics of M33.

It is certainly possible to obtain a better model. The simulation including the hot gas halo of M31 looks very promising, and exploring this scenario further with a code that resolves the RPS better (Agertz et al. 2007; Hopkins 2015) would be very interesting. Probing better the parameter space of orbits and inclinations in this setup perhaps could also give a better result. Given the huge uncertainty of the proper motions of M31, we do not exclude the possibility that the true velocity vector lies somewhere between the sets of vdM12 and S16.

Most features of M33 that were induced due to the interaction with M31 in our simulation are of tidal origin and therefore depend strongly on the orbit, inclinations, and mass models. The subgrid physics seems not to be very important, and our results are rather robust against different models of star formation, etc. The peak in SFH after the pericenter passage might be the only observable that could have been changed by subgrid physics. We tested our model with a code with different subgrid prescriptions, and a similar effect of the increase in SFR during the pericenter passage was also found. Some extreme model with a very low density threshold and very high feedback could probably alter this result and prevent a strong increase in star-forming activity.

## 5.2. Comparison with Previous Works

The first model of the interaction between M33 and M31 was presented by Bekki (2008). The test-particle simulations presented there were tailored to reproduce the H I bridge found earlier by Braun & Thilker (2004). In our model, we do not find such a bridge-like structure, as was the case in Bekki (2008). In the fiducial simulations without RPS, we obtained instead tidal arms that were extended up to ~4°–5° in both directions

(similarly to the stellar tails in Figure 10). These features, however, were dissolved in the run with RPS that reproduced better the outer gaseous structures of M33. This discrepancy with the results of Bekki (2008) comes from many differences between our works; however, the most important is probably the inclusion of hydrodynamics and hot gas halo in our case, since the run without it produced similar features. Our finding of no bridge-like structure is in agreement with recent results (e.g., Wolfe et al. 2013) that question the possibility that H I between M31 and M33 is due to the past interaction between the two galaxies. Apart from the bridge, Bekki (2008) does not discuss other features associated with the interaction scenario, hence it is difficult to further compare with his work.

The second model of the interaction was presented by McConnellie et al. (2009). They performed high-resolution N-body simulations in order to reproduce the stellar distortion of the disk, the discovery of which was reported in the very same paper. Similarly to Bekki (2008), they adopted a small (8 × 10¹⁰ M☉) halo mass of M33 and a relatively large halo mass of M31 (2.47 × 10¹² M☉). This assumption of a small mass ratio of only ~3% between M33 and M31 makes the modeling of the orbit much easier, since such a system can be approximated as a small satellite orbiting a massive spiral. We assumed such masses of both galaxies with the aim to reproduce their rotation curves (C14; Corbelli et al. 2010), and this resulted in the mass ratio of ~26%, for which the system has to be treated more like a two-body interaction of similar-sized objects.

Another difference between our work and McConnellie et al. (2009) is that their orbit, just like Bekki’s, was not constrained by the measurement of the M31 transverse velocity, which was unknown at that time. Despite these different assumptions, we find that our model agrees with the one presented by McConnellie et al. (2009). The pericenter distance they used was similar (~50 kpc), and the warp in their disk looks very similar to ours, especially for the stellar component (Figure 10). The stream-like structure, just as in our case, is made of tidally distorted material that wound up and, in a specific projection, has a different orientation than the disk. After the pericenter passage, their M33 ends up more distant from M31 than ours (the apocenter is ~260 kpc in their model and ~151 kpc in ours), but this is a result of the smaller mass ratio discussed earlier. Unfortunately, any quantitative results regarding, e.g., the spiral arms or SFH (due to the lack of hydrodynamics in their simulations) were not discussed by McConnellie et al. (2009) so we cannot compare our model with theirs in more detail.

Besides the work in the literature that aimed to reproduce the hypothetical past interaction between M33 and M31, several authors have argued that such a scenario was unlikely using different orbital studies. First, Shaya & Tully (2013) used the numerical action method and backward orbit integration to conclude that M33 and M31 are currently at the closest approach. According to the information in the appendix of their paper, their scheme did not include the effect of dynamical friction. However, as shown in Appendix A of van der Marel et al. (2012a), orbital studies with and without implementation of dynamical friction may vary a lot (i.e., the apocenter may change from 800 kpc to 100 kpc), with the former agreeing much better with the simulations. After finding that M33 is now at its closest approach to M31, Shaya & Tully (2013) suggested that Andromeda XXII (And XXII), the only candidate for a
satellite of M33 (Chapman et al. 2013), could have disturbed the disk of M33. We will discuss this alternative scenario more in Section 5.3.

More recently, Patel et al. (2017b, 2017a) explored the past orbital history of two systems, the MW–Large Magellanic Cloud and M31–M33. In the first paper (Patel et al. 2017b), they applied backward orbit integration to find possible orbits of both systems and compared them to orbits from dark matter only cosmological simulations. In Section 7.1 of their paper, they state that from backward orbit integration of 10,000 initial velocities drawn from $4\sigma$ proper-motion error space, they found that less than 1% of orbits had a pericenter smaller than 100 kpc in the last 3 Gyr. As can be seen in Figure 13 of their paper, these 10,000 velocities cover well the area surrounding the velocity that results from adopting the van der Marel et al. (2012b) proper motions, while the area near the results of Salomon et al. (2016) is not probed.

In Section 2 of our paper, we extend these results of Patel et al. (2017b) to also probe the region surrounding estimates of the velocity resulting from adopting the values of Salomon et al. (2016). Our findings agree with the conclusion of Patel et al. (2017b) that once one adopts the transverse velocity of M31 from van der Marel et al. (2012b), the recent and close pericenter passage is not very plausible. This agreement also holds if we take into account the fact that we used different mass models of the galaxies (in particular, M33 is almost twice as massive in our work). The conclusions about the recent pericenter change once the estimates of van der Marel et al. (2012b) are replaced by those of Salomon et al. (2016). Since the publication of the most recent estimates of the M31 transverse velocity, no other measurements were made to support one or the other of the conflicting estimates. Furthermore, cosmological simulations (Carlesi et al. 2016) do not strongly exclude any of these estimates. Moreover, such simulations do not exclude a close pericenter passage, as discussed by Patel et al. (2017b), where the analysis of orbits from the Illustris-Dark (Vogelsberger et al. 2014) simulation showed that orbits with pericenters even lower than 55 kpc are very common (~30%) for similar pairs of halos. The probability of orbits with pericenters smaller than 50 kpc was found to be even higher in a recent study of orbits of similar pairs in cosmological simulations that included baryonic physics (Shao et al. 2018).

In the second paper, Patel et al. (2017a) used a Bayesian inference scheme and the Illustris-Dark simulations to estimate the masses of MW and M31. One of the conclusions of this paper was that after adopting the criteria for the recent interaction scenario, the mass of M31 would be ~$10^{12} M_{\odot}$. This value is not far from the estimates of the M31 mass found in the literature; for example, it lies within the error bars of the result of Watkins et al. (2010), $1.4 \pm 0.4 \times 10^{12} M_{\odot}$. We argue that the precision of the estimates of the mass of M31 is not high enough to rule out ~$10^{12} M_{\odot}$ and abandon the interaction scenario.

## 5.3. Alternative Scenarios

Both Shaya & Tully (2013) and Patel et al. (2017b) proposed a passage of And XXII as a new alternative scenario that could disturb the gaseous and stellar disk of M33. Here we try to roughly estimate the impact of the tidal force that such a passage would have had and compare it with the impact in the scenario presented in this paper. Elmegreen et al. (1991)

defined the parameter $S$ that quantifies the strength of the tidal interaction. Their definition makes use of the masses of the perturber and the perturbed galaxy, the pericentric distance, the size of the perturbed galaxy as well as the timescales of the interaction and of the motions of the stars in the perturbed disk. Since constraining the orbit of And XXII around M33 would give too many free parameters, we decided to apply instead a parameter that does not use the timescales and was discussed, e.g., in Oh et al. (2015):

$$P = \left( \frac{M_{\text{ptb}}}{M_{\text{g}}} \right) \left( \frac{R_{\text{g}}}{d} \right)^{3} \qquad (1)$$

Here, $M_{\text{ptb}}$ and $M_{\text{g}}$ denote the masses of the perturber and the perturbed galaxy, $R_{\text{g}}$ is the size of the perturbed galaxy, which we take to be five times the disk scale length, and $d$ is the pericentric distance. For the case of the simulation described in this paper, the tidal parameter takes the value of $P = 0.03$, where we assumed that the mass of the perturber is the mass of M31 enclosed within the pericentric distance. We know from our simulation that the tidal interaction characterized by this number can disturb the disk of M33 and the induced morphology mimics the observed one to a certain degree.

We now want to find out what size of pericenter And XXII would need to have in order to provide the same tidal impact. In order to estimate this, we need the mass of And XXII. Adopting the value of Shaya & Tully (2013), $1.3 \times 10^7 M_{\odot}$, we obtain the pericenter of 1.3 kpc. Tollerud et al. (2012) estimated the mass in And XXII enclosed within the half-light radius to be $log(M_{1/2}/M_{\odot}) = 6.67 \pm 1.08$. If we take the upper limit of this measurement and assume that the half-light radius encompasses only 10% of the total mass, the estimated pericenter will still be only 4.5 kpc. Such small pericenter values suggest that the low-mass And XXII would have to pass extremely close to M33 in order to have the tidal impact similar to the one induced by M31 in our simulation. We expect that such a close passage would perturb more the inner parts of the disk of M33, which seems to be undisturbed. The warp and the stellar stream occur at much greater radii. We find that the scenario in which And XXII induced the warp in M33 is less likely, due to the small mass of this satellite.

Another possibility is that the more massive satellites of M33 disturbed its outer disk, but they are yet to be discovered. This scenario could perhaps explain the warped gaseous feature, since gas is less massive and hence less bound and easier to distort. A different possible scenario for the gas disk distortion could also be an asymmetric gas accretion. However, these two explanations would not be very convincing in the case of the stellar component. Studying MW analogs, Gómez et al. (2017) found that a satellite should have a mass of at least 1% of the host to induce some sort of vertical distortion in the disk. It would be very surprising if a satellite of M33 of a corresponding brightness existed and has not yet been discovered. Perhaps the distorted stellar and gaseous disks are remnants of a major merger that M33 underwent in the past, but this scenario would raise another question of why the inner disk remained unaffected by such an event. Hopefully, future observations and simulations will shed more light on the peculiar history of M33.

Very recently, van der Marel et al. (2018) utilized Gaia DR2 results to estimate the proper motions of both M31 and M33. Based on the agreement of the new result for M31 with their
previous *HST* estimate and on the orbit integration sample from Patel et al. (2017b), they argue that M33 may be on its first infall into M31. We find this interpretation unlikely for the following three reasons. First, the error bars of the van der Marel et al. (2018) estimates are very large, and the resulting 3D relative velocity between M31 and M33, based solely on *Gaia* DR2 measurements, is in $1.75\sigma$ agreement with the relative velocity for the best match of our model (Section 4.2), in which the two galaxies recently interacted. Second, the orbit integration sample from Patel et al. (2017b) was selected from the $4\sigma$ space around the van der Marel et al. (2012b) measurements and naturally will fall far from the new measurements (in Figure 4 of van der Marel et al. 2018), as that region was not fully probed and therefore the distance between the new results and points with a pericenter might just be an artifact of the probing of the phase space. Finally and most importantly, we do not consider a semi-analytic orbit integration method to be a robust computational tool that can rule out with any confidence one scenario or another for the past history of M33 and M31. Orbits obtained by this method do not reproduce the simulations accurately, since it does not include mass loss originating from tidal stripping and parameters of the dynamical friction have to be fitted to a particular setup.

Moreover, in the close proximity to M31, there are some massive satellite galaxies that would most likely have some influence on the M33–M31 orbit. The mass of M32 was estimated to be $8 \times 10^{10} M_{\odot}$ (Dierickx et al. 2014), which is around 20% of the mass of M33, and M110 has a similar magnitude to M32 (McConnachie 2012). Including these factors in the orbit integration would create a lot of degeneracy, since their proper motions are unknown and finding the true orbit between M33 and M31 would be even more complicated. The simple scenario of the interaction that we discussed here favored more the Salomon et al. (2016) estimates, but perhaps the influence of M32 and M110 would change this conclusion, and other estimates of the transverse velocity of M31 would be in agreement with the M33–M31 interaction scenario.

## 6. Summary

In this work, we revisited the scenario discussed in P09 and McConnachie et al. (2009), which was proposed to explain the stellar and gaseous disks of M33 which were disturbed by a recent passage close to M31. We used the orbit integration method to verify which of the measurements of the transverse velocity of M31 favors this scenario more. We found that while the estimates of van der Marel et al. (2012b) do not support this common history of the two galaxies, the estimates of Salomon et al. (2016) allow for it. We performed N-body/SPH simulations, aiming to reproduce the observed disturbed morphology of M33 and at the same time to be consistent with the 3D relative position and velocity of the galaxies resulting from the estimates of Salomon et al. (2016). The fiducial setup presented here fulfills the orbital conditions and had the pericenter at the distance of 37 kpc, which was close enough to tidally disturb the disk of M33. Mass models of both galaxies in our simulations were constructed to roughly reproduce the observed rotation curves.

We found that the tidal impulse originating from such an interaction is sufficient to excite a two-armed spiral structure similar to the one found to be the dominant spiral component in

the observed M33. Tides also induced distortion in the stellar and gaseous disks at larger radii, with the former having the shape and the extent similar to the observed one. The disturbance of the gaseous disk, however, was found to be similar in the inner parts while in the outer parts, it was dominated by strong tidal/spiral features that are not present in the observed data. We showed by performing an additional run including the hot gas halo of M31 that the RPS is a crucial component in modeling the gaseous warp of M33 in greater detail.

Finally, we also found that the tidal forces in our simulations were sufficient to compress the gas in M33 during the pericenter passage and trigger a burst of star formation at similar radii to those found by Bernard et al. (2012) and hypothesized to be due to the passage near M31, since similar activity was found in its SFH at approximately the same time. The model presented here did not aim to reproduce the observed M33 in great detail, and this was not achieved; rather, its aim is to demonstrate that observationally constrained structural and orbital parameters of the system allow for the interaction to trigger in M33 features similar to the observed ones.

This work was supported in part by the Polish National Science Centre under grant 2013/10/A/ST9/00023. We are grateful to an anonymous referee for useful comments that helped to improve the paper. We thank L. Widrow for providing procedures to generate N-body models of galaxies for initial conditions. We are grateful to S. Fouquet for sharing the modified version of the GADGET-2 code and to his collaborators F. Hammer and J. Wang for incorporating the modifications. We thank S. Z. Kam and L. Chemin for sharing the H I data of M33 as well as E. Corbelli, M. Putman, and the editorial board of A&A for granting us the permission to reprint certain figures. Acknowledgments are also due to J.-C. Lambert for his tutorial about the tool for the visualization of the simulations gnemo2 and for making it public. We appreciate insightful discussions with G. Gajda, K. Kowalczyk, I. Ebrova, N. Peschken, and J. S. Gallagher III that contributed to this paper. M.S. is grateful for the hospitality of Laboratoire d'Astrophysique de Marseille and the University of Wisconsin at the time of his visits. E.A. acknowledges financial assistance from the CNES and access to the HPC resources of CINES under the allocation 2017-[A0040407665] made by GENCI. E. D. acknowledges the support and the hospitality of the Center for Computational Astrophysics (CCA) at the Flatiron Institute during the preparation of this work and the Vilas Associate Professor Fellowship at the University of Wisconsin, Madison. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

## Appendix A
### Relative Positions and Velocities

Following van der Marel et al. (2012b), we adopted the distance to M33 and its error as $D_{M33} = 794 \pm 23$ kpc (McConnachie et al. 2004) and for M31, $D_{M31} = 770 \pm 40$ kpc (van der Marel & Guhathakurta 2008 and references therein). We combined these values with the distance of the Sun from the Galactic center $R_0 = 8.29 \pm 0.16$ kpc (McMillan 2011) to obtain the relative position vector between M33 and
M31 in the so-called Galactocentric rest frame. This Cartesian reference frame (X, Y, Z) is defined as follows: the center lies on the Galactic center, the X-axis points in the direction from the Sun to the Galactic center, the Y-axis points in the direction of the Sun's rotation in the Galaxy, and the Z-axis points toward the Galactic north pole. The resulting relative position vector between M33 and M31 with error bars is

$$X_{\text{rel}} = (-97.2 \pm 23.5, -121.6 \pm 34.8, \\ -129.8 \pm 19.0) \text{ kpc.} \qquad (2)$$

In order to calculate the relative velocity between the two galaxies in the same reference frame, first we used the values of the proper motions of M33 from Brunthaler et al. (2005), the LOS velocity as used in van der Marel et al. (2012b), $v_{\text{LOS,M33}} = -180 \pm 1 \text{ km s}^{-1}$ (van der Marel & Guhathakurta 2008), and applied the correction for the motion of the Sun ($U_{\odot}, V_{\odot} + V_{\text{LSR}}, W_{\odot}$) = (11.1 $\pm$ 0.7, 255.2 $\pm$ 5.1, 7.25 $\pm$ 0.36) km s$^{-1}$ (Schönrich et al. 2010; Reid et al. 2014). Then we considered the proper motions and LOS velocity of M31 as given by van der Marel et al. (2012b), corrected for the same solar values. This yields the first option for the relative velocity,

$$V_{\text{rel,vdM12}} = (-23.2 \pm 34.3, 177.4 \pm 29.2, \\ 93.7 \pm 38.5) \text{ km s}^{-1}. \qquad (3)$$

The second possible value of the relative velocity was obtained by replacing the values for M31 with those given by Salomon et al. (2016) and correcting for the same solar values. This yields

$$V_{\text{rel,S16}} = (-72.0 \pm 64.4, 86.4 \pm 48.0, 10.6 \pm 62.1) \text{ km s}^{-1}. \qquad (4)$$

For all calculations done here, we assumed that error bars follow a Gaussian distribution. We have adopted this approximation in order to be able to easily transform error bars from one reference frame to another.

## Appendix B: Orbit Integration Scheme

In order to quickly find possible orbits of the two galaxies without running computationally costly simulations, we performed semi-analytic orbit integrations, similar to those described in Patel et al. (2017b) and Dierickx & Loeb (2017a, 2017b). We treated both galaxies as two interacting NFW halos (Navarro et al. 1997), and we included the dynamical friction from M31 acting on M33. We thus solved the following equations of motion:

$$\begin{aligned} \ddot{x}_{\text{M33}} &= -\nabla\psi_{\text{M31}} + f_{\text{DF}}, \\ \ddot{x}_{\text{M31}} &= -\nabla\psi_{\text{M33}}. \end{aligned} \qquad (5)$$

Here, $x_{\text{M33}}$ and $x_{\text{M31}}$ denote the positions of the M33 and M31 galaxies. The functions $\psi_{\text{M33}}$ and $\psi_{\text{M31}}$ are the NFW potentials of M33 and M31, respectively, given by

$$\psi_i = -\frac{GM_i}{r[\ln(1+c_i) - c_i/(1+c_i)]} \ln\left(1 + \frac{r}{r_{s,i}}\right), \qquad (6)$$

where $M_i$ are the virial masses of the galaxies, $c_i$ are their concentration parameters, and $r_{s,i} = r_{v,i}/c_i$ are the scale radii, i.e., the virial radii $r_{v,i}$ divided by concentrations. For both

galaxies, we adopted halo parameters that were estimated by modeling the observed rotation curves. For M33, we assumed $M_{\text{M33}} = 4.38 \times 10^{11} M_{\odot}$ and $c_{\text{M33}} = 11$ (C14). The halo mass includes the estimate of the mass of baryons in M33, discussed in greater detail in Section 3.1. For M31, we took $M_{\text{M31}} = 2 \times 10^{12} M_{\odot}$ and $c_{\text{M31}} = 28$ to fit the rotation curve given by Corbelli et al. (2010) only with the potential of the halo (see Section 3.1). For the supplementary parameters necessary to compute the NFW profile, we adopted the virial overdensity to be $\Delta_c = 102$ and the critical density $\rho_c^0 = 136 M_{\odot} \text{ kpc}^{-3}$.

Finally, $f_{\text{DF}}$ in Equation (5) denotes the acceleration caused by the dynamical friction, approximated in our calculations according to the Chandrasekhar formula (Chandrasekhar 1943):

$$f_{\text{DF}} = -\frac{4\pi G^2 M_{\text{M33}} \ln \Lambda \rho(r)}{\nu^2} \left[ \text{erf}(X) - \frac{2X}{\sqrt{\pi}} \exp(-X^2) \right] \frac{\nu}{\nu}, \qquad (7)$$

where $\rho(r)$ is the density of the M31 halo at a given radius $r$, $\nu$ is the relative velocity vector with its magnitude $|\nu| = \nu$, and $X = \nu/\sqrt{2\sigma}$, where $\sigma$ is the 1D velocity dispersion, which we approximated according to the formula derived by Zentner & Bullock (2003). For the Coulomb logarithm $\ln \Lambda$, we adopted the formula presented in Hashimoto et al. (2003) and used by Dierickx & Loeb (2017a, 2017b):

$$\ln \Lambda = \ln \left( \frac{r}{1.4\epsilon} \right), \qquad (8)$$

where $\epsilon$ is the softening length of the galaxy subject to dynamical friction. We fit this parameter for M33 by comparing the integrated orbits with orbits in preliminary, collisionless, low-resolution simulations. We estimated the value $\epsilon = 28.5 \text{ kpc}$ to be the best match. To integrate Equations (5), we used a symplectic leapfrog integration method as described in Springel et al. (2001).

## Appendix C: Simulations in Isolation

In order to confirm that the morphological features like the gaseous warp, the stellar stream, and the grand-design spiral arms originate from the tidal interaction with M31 and are not induced by secular processes (e.g., feedback or star formation) in the disk of M33, we evolved this galaxy in isolation for the same time as in the case with the perturber. Figure 17 presents surface density maps for all three components of M33 evolved in isolation at 3.45 Gyr, the time of the best match for the fiducial model. A quick look at these images and the comparison with Figure 6 is enough to confirm that the aforementioned features are indeed tidally induced, since none of the components of the isolated galaxy shows signs of similar distortions. Particles present above and below the gaseous disk in Figure 17 that are distributed in an hourglass-like shape have been ejected there by feedback and contribute a very small fraction of the total gaseous content. They are absent in the images of the fiducial model, most likely because of tidal stripping. The young stellar disk also reveals a small bar which is not a tidal feature as discussed in Section 4.5.
<figure><img src="image_23.png" /><figcaption>**Figure 17.** Surface density distributions seen from three different directions for the gas and the old and young stellar particles for the M33 galaxy simulated in isolation at the time corresponding to the time of the best match for the fiducial model.</figcaption></figure>

## ORCID iDs

Marcin Semczuk https://orcid.org/0000-0002-8191-8918

Ewa L. Łokas https://orcid.org/0000-0001-7138-8899

E. Athanassoula https://orcid.org/0000-0001-6079-1332

## References

* Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963

* Arthur, J., Pearce, F. R., Gray, M. E., et al. 2017, MNRAS, 464, 2027

* Athanassoula, E. 2013, in Bars and Secular Evolution in Disk Galaxies: Theoretical Input, ed. J. Falcón-Barroso & J. H. Knapen (Cambridge: Cambridge Univ. Press), 305

* Athanassoula, E., Bosma, A., & Papaioannou, S. 1987, A&A, 179, 23

* Athanassoula, E., Machado, R. E. G., & Rodionov, S. A. 2013, MNRAS, 429, 1949

* Barker, M. K., Ferguson, A. M. N., Cole, A. A., et al. 2011, MNRAS, 410, 504

* Bekki, K. 2008, MNRAS, 390, L24

* Bernard, E. J., Ferguson, A. M. N., Barker, M. K., et al. 2012, MNRAS, 420, 2625

* Besla, G., Kallivayalil, N., Hernquist, L., et al. 2012, MNRAS, 421, 2109

* Braun, R., & Thilker, D. A. 2004, A&A, 417, 421

* Brunthaler, A., Reid, M. J., Falcke, H., et al. 2005, Sci, 307, 1440

* Byrd, G., & Valtonen, M. 1990, ApJ, 350, 89

* Carlesi, E., Hoffman, Y., Sorce, J. G., et al. 2016, MNRAS, 460, L5

* Chandrasekhar, S. 1943, ApJ, 97, 255

* Chapman, S. C., Widrow, L., Collins, M. L. M., et al. 2013, MNRAS, 430, 37
Considere, S., & Athanassoula, E. 1988, A&AS, 76, 365

Corbelli, E., Lorenzoni, S., Walterbros, R., et al. 2010, A&A, 511, 89

Corbelli, E., & Schneider, S. E. 1997, ApJ, 479, 244

Corbelli, E., Thilker, D., Zibetti, S., et al. 2014, A&A, 572, 23

Corbelli, E., & Walterbros, R. A. M. 2007, ApJ, 669, 315

Courtes, G., & Dubout-Crillon, R. 1971, A&A, 11, 468

Cox, T. J., Jonsson, P., Primack, J. R., & Somerville, R. S. 2006, MNRAS, 373, 1013

D'Onghia, E. 2015, ApJL, 808, L8

D'Onghia, E., Besla, G., Cox, T. J., & Hernquist, L. 2009, Natur, 460, 605

D'Onghia, E., & Fox, A. J. 2016, ARA&A, 54, 363

D'Onghia, E., Madau, P., Vera-Ciro, C., Quillen, A., & Hernquist, L. 2016, ApJ, 823, 4

D'Onghia, E., Vogelsberger, M., Faucher-Giguere, C.-A., & Hernquist, L. 2010, ApJ, 725, 353

Dierickx, M., Blecha, L., & Loeb, A. 2014, ApJL, 788, L38

Dierickx, M., & Loeb, A. 2017a, ApJ, 836, 92

Dierickx, M., & Loeb, A. 2017b, ApJ, 847, 42

Dobbs, C., & Baba, J. 2014, PASA, 31, e035

Dobbs, C. L. 2011, MSAIS, 18, 109

Dobbs, C. L., Theis, C., Pringle, J. E., & Bate, M. R. 2010, MNRAS, 403, 625

Elmegreen, B. G., Elmegreen, D. M., & Montenegro, L. 1992, ApJS, 79, 37

Elmegreen, D. M., Sundin, M., Sundelius, B., & Elmegreen, B. 1991, A&A, 244, 52

Evans, N. W., Wilkinson, M. I., Guhathakurta, P., Grebel, E. K., & Vogt, S. S. 2000, ApJL, 540, L9

Fouquet, S., Lokas, E. L., del Pino, A., & Ebrová, I. 2017, MNRAS, 464, 2717

Gajda, G., Łokas, E. L., & Athanassoula, E. 2017, ApJ, 842, 56

Gieren, W., Górski, M., Pietrzyński, G., et al. 2013, ApJ, 773, 69

Gómez, F. A., White, S. D. M., Grand, R. J. J., et al. 2017, MNRAS, 465, 3446

Gunthardt, G., Diaz, R. J., & Aguero, M. P. 2006, AJ, 152, 150

Hammer, F., Yang, Y. B., Wang, J. L., et al. 2010, ApJ, 725, 542

Hashimoto, Y., Funato, Y., & Makino, J. 2003, ApJ, 582, 196

Hernández-López, I., Athanassoula, E., Mújica, R., & Bosma, A. 2009, RMxAC, 37, 160

Hernquist, L. 1990, ApJ, 356, 359

Hopkins, P. F. 2015, MNRAS, 450, 53

Horellou, C., & Combes, F. 2001, Ap&SS, 276, 1141

Humphreys, R. A., & Sandage, A. 1980, ApJS, 44, 319

Hunter, C., & Toomre, A. 1969, ApJ, 155, 747

Kahn, F. D., & Woltjer, L. 1959, ApJ, 130, 705

Kam, Z. S., Carignan, C., Chemin, L., et al. 2017, AJ, in press (arXiv:1706.04248)

Kam, Z. S., Carignan, C., Chemin, L., Amram, P., & Epinat, B. 2015, MNRAS, 449, 4048

Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 105, 19

Keenan, O. C., Davies, J. I., Taylor, R., & Minchin, R. F. 2016, MNRAS, 456, 951

Lehner, N., Howk, J. C., & Wakker, B. P. 2015, ApJ, 804, 79

Lewis, G. F., Braun, R., McConnachie, A. W., et al. 2013, ApJ, 763, 4

Lockman, F. J., Free, N. L., & Shields, J. C. 2012, AJ, 144, 52

Łokas, E. L., Athanassoula, E., Debattista, V. P., et al. 2014, MNRAS, 445, 1339

Łokas, E. L., Ebrová, I., del Pino, A., et al. 2016, ApJ, 826, 227

Łokas, E. L., Semczuk, M., Gajda, G., & D'Onghia, E. 2015, ApJ, 810, 100

Massey, P., Olsen, K. A. G., Hodge, P. W., et al. 2006, AJ, 131, 2478

McConnachie, A. W. 2012, AJ, 144, 4

McConnachie, A. W., Ferguson, A. M. N., Irwin, M. J., et al. 2010, ApJ, 723, 1038

McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., et al. 2004, MNRAS, 350, 243

McConnachie, A. W., Irwin, M. J., Ibata, R. A., et al. 2009, Natur, 461, 66

McMillan, P. J. 2011, MNRAS, 414, 2446

Miller, M. J., & Bregman, J. N. 2015, ApJ, 800, 14

Navarro, J. F., Frenk, C. S., & White, S. D. 1997, ApJ, 490, 493

Newton, K. 1980, MNRAS, 190, 689

Oh, S. H., Kim, W.-T., & Mok, H. 2015, ApJ, 807, 73

Pardy, S. A., D'Onghia, E., & Fox, A. J. 2018, arXiv:1802.01600

Patel, E., Besla, G., & Mandel, K. 2017a, MNRAS, 468, 3428

Patel, E., Besla, G., & Sohn, S. T. 2017b, MNRAS, 464, 3825

Pettitt, A. R., Tasker, E. J., & Wadsley, J. W. 2016, MNRAS, 458, 3990

Pettitt, A. R., & Wadsley, J. W. 2018, MNRAS, 474, 5645

Puerari, I. 1993, PASP, 105, 1290

Putman, M. E., Peek, J. E. G., Muratov, A., et al. 2009, ApJ, 703, 1486

Regan, M. W., & Vogel, S. N. 1994, ApJ, 434, 536

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Renaud, F., Athanassoula, E., Amram, P., et al. 2018, MNRAS, 473, 585

Renaud, F., Bournaud, F., & Duc, P.-A. 2015, MNRAS, 446, 2038

Rogstad, D. H., Wright, M. C. H., & Lockhart, I. A. 1976, ApJ, 204, 703

Ruggiero, R., & Lima Neto, G. B. 2017, MNRAS, 468, 4107

Saha, K., & Elmegreen, B. 2016, ApJL, 826, L21

Salo, H., & Laurikainen, E. 2000, MNRAS, 319, 393

Salomon, J.-B., Ibata, R. A., Famaey, B., et al. 2016, MNRAS, 456, 4432

Sandage, A., & Humphreys, R. M. 1980, ApJL, 236, L1

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Sellwood, J. A., & Athanassoula, E. 1986, MNRAS, 221, 195

Semczuk, M., Łokas, E. L., & del Pino, A. 2017, ApJ, 834, 7

Shao, S., Cautun, M., Deason, A. J., & Frenk, C. S. 2018, MNRAS, 479, 284

Shaya, E. J., & Tully, R. B. 2013, MNRAS, 436, 2096

Sohn, S. T., Anderson, J., & van der Marel, R. P. 2012, ApJ, 753, 7

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Yoshida, N., & White, S. D. M. 2001, NewA, 6, 79

Teyssier, R., Chapon, D., & Bournaud, F. 2010, ApJL, 720, L149

Theis, C., & Spinneker, C. 2003, Ap&SS, 284, 495

Tollerud, E. J., Beaton, R. L., Geha, M. C., et al. 2012, ApJ, 752, 45

van der Marel, R. P., Besla, G., Cox, T. J., Sohn, S. T., & Anderson, J. 2012a, ApJ, 753, 9

van der Marel, R. P., Fardal, M., Besla, G., et al. 2012b, ApJ, 753, 8

van der Marel, R. P., Fardal, M. A., Sohn, S. T., et al. 2018, arXiv:1805.04079

van der Marel, R. P., & Guhathakurta, P. 2008, ApJ, 678, 187

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Wang, J. L., Hammer, F., Athanassoula, E., et al. 2012, A&A, 538, A121

Watkins, L. L., Evans, N. W., & An, J. H. 2010, MNRAS, 406, 264

Widrow, L. M., & Dubinski, J. 2005, ApJ, 631, 838

Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239

Wolfe, S. A., Lockman, F. J., & Pisano, D. J. 2016, ApJ, 816, 81

Wolfe, S. A., Pisano, D. J., Lockman, F. J., McGaugh, S. S., & Shaya, E. J. 2013, Natur, 497, 224

Zentner, A. R., & Bullock, J. S. 2003, ApJ, 598, 49