Large-scale Hydrodynamical Shocks as the Smoking Gun Evidence for a Bar in M31

The formation and evolutionary history of M31 are closely related to its dynamical structures, which remain unclear due to its high inclination. Gas kinematics could provide crucial evidence for the existence of a rotating bar in M31. Using the position-velocity diagram of [OIII] and HI, we are able to identify clear sharp velocity jump (shock) features with a typical amplitude over 100 km/s in the central region of M31 (4.6 kpc X 2.3 kpc, or 20 arcmin X 10 arcmin). We also simulate gas morphology and kinematics in barred M31 potentials and find that the bar-induced shocks can produce velocity jumps similar to those in [OIII]. The identified shock features in both [OIII] and HI are broadly consistent, and they are found mainly on the leading sides of the bar/bulge, following a hallmark pattern expected from the bar-driven gas inflow. Shock features on the far side of the disk are clearer than those on the near side, possibly due to limited data coverage on the near side, as well as obscuration by the warped gas and dust layers. Further hydrodynamical simulations with more sophisticated physics are desired to fully understand the observed gas features and to better constrain the parameters of the bar in M31.


INTRODUCTION
Although M31 is the nearest large spiral galaxy to the Milky Way at a distance of 785 ± 25 kpc (McConnachie et al. 2005), its exact location on the Hubble tuning fork diagram is still unclear, which is key to understand its formation and evolution. The debates have a quite long history that probably started from Lindblad (1956) who first claimed that M31 is a barred galaxy based on the twist of central isophotes. The isophotal twist cannot be reproduced by an axisymmetric stellar distribution (Stark 1977). Later studies argued that unbarred galaxies can also have twisted inner isophotes as they can originate from triaxial bulges, not necessarily bars (Stark 1977;Zaritsky & Lo 1986;Bertola et al. 1988;Gerhard et al. 1989;Méndez-Abreu et al. 2010;Costantin et al. 2018), which left the true morphology of M31 a puzzle. In addition, the strong enhancement of star formation between 2 − 4 Gyr ago (Williams et al. 2015) could be cal bulge and the box/peanut bulge (BPB) in M31 contribute about one-third and two-thirds of the total stellar mass in the central part, respectively. Blaña Díaz et al. (2018) further constructed made-to-measure (m2m) models for M31, using constraints from 3.6 µm photometry (Barmby et al. 2006) and stellar kinematics . The models evolved from N−body buckled bar models. The BPB results from a buckled bar in their best model, with a half-length of ∼ 4 kpc and a pattern speed of 40 ± 5 km s −1 kpc −1 . They also checked a model with a very low pattern speed and found it does not match the data as nicely as the more rapidly rotating models. Moreover, Saglia et al. (2018) found that the stellar metallicity is enhanced along the proposed bar structure. Based on the m2m models in Blaña Díaz et al. (2018), Gajda et al. (2021) constrained the 3D distribution of metallicity and α-enrichment using observed [Z/H] and [α/Fe] maps , and an X-shaped metallicity distribution is found in the bulge region. These results imply that M31 is a barred galaxy.
On the other hand, evidence from gas kinematics in M31 also hints for a bar. Large non-circular motions of gas have been identified in H I (Brinks & Burton 1984;Chemin et al. 2009) and CO (Loinard et al. 1995(Loinard et al. , 1999Nieten et al. 2006). In addition, [OIII] emission from ionized gas shows a twisted zero-velocity curve in the central bulge region  , and its morphology suggests a tilted spiral pattern with a lower inclination angle compared to the stellar disk ). All of these are typical features seen in barred galaxies (Jacoby et al. 1985;Emsellem et al. 2006;Kuzio de Naray et al. 2009;Fathi et al. 2005), although there are alternative explanations (e.g. a head-on collision between M31 and M32 proposed by Block et al. 2006, non-circular motions of gas caused by a triaxial bulge). Nevertheless, the overall shape of position-velocity diagrams (PVDs) of [O III] in Opitsch et al. (2018) is similar to what has been observed in barred galaxies Merrifield & Kuijken 1999).
One of the most characteristic gaseous features in typical barred galaxies is a pair of dust-lanes on the leading side of the bar (Athanassoula 1992). Dust-lanes are generally associated with shocks, which result in sharp velocity jumps in gas kinematical maps, and thus can be identified by integralfield unit (IFU) data (e.g. Opitsch et al. 2018).
The goal of the current paper is to search and identify such shock (velocity jump) features from the observations. It should be noted that the shocks we are searching for are "large-scale" shocks (e.g. of the type proposed by Roberts 1969;Roberts et al. 1979) which are generally extended over a few kiloparsecs, rather than those "small-scale" shocks due to local turbulence of interstellar medium or supernova feedbacks. If the positions and properties of velocity jumps are similar to those expected from bar-driven shocks, this would provide independent, strong evidence for the existence of a bar in M31.
The paper is organized as follows: §2 describes the data of [ with the bar-driven shocks in hydrodynamical simulations. §6 mainly discusses the potential physical reasons for the asymmetry in the shock features between the far side and near side of M31. We briefly summarize our findings in §7.

OBSERVATIONAL DATA
The main results of the current paper are summarized in Fig. 1, which shows the positions of the possible bar-driven shocks superposed on the optical image of M31. The distance to M31 is adopted to be 785 ± 25 kpc (McConnachie et al. 2005), so 1 corresponds to 228 pc.
We use data from Opitsch et al. (2018) and Chemin et al. (2009) to study the gas kinematics in M31. Opitsch et al. (2018)  , respectively. Multiple gas components are expected when the line of sight passes through different gas streams, which could be caused by a bar (Kim et al. 2012) or a collision between M31 and its satellite galaxy M32 (Block et al. 2006). We use the main (higher-velocity) component of [O III] 1 to present the PVDs. We discuss the possible origin of the two components in §6.5. Similar to [O III], H I observations by Chemin et al. (2009) show multiple velocity components in their spectra. The H I disk is extended and starts to warp beyond ∼ 20 kpc (Newton & Emerson 1977;Henderson 1979). The H I component with a lower velocity is likely from the warped H I layer in the outer disk (Brinks & Burton 1984;Brinks & Shane 1984). Chemin et al. (2009) shows that gas from the warped region dominates the H I emission in the central 20 , producing shallow linear structures on their PVDs. The warp may obscure the inner H I disk with similar velocities on PVDs. To find the main H I component that best represents the disk rotation, Chemin et al. (2009) selected the component with the largest velocity relative to the galactic systemic velocity while re- The background image is taken from https://apod.nasa.gov/apod/ap220119.html. The image is credited to Subaru (NAOJ), Hubble (NASA/ESA), Mayall (NSF) and processed by R. Gendler and R. Croman. X-axis and Y-axis are taken along the major and minor axes of the stellar disk of M31, respectively. Solid, open, and dashed markers represent the Class I, Class II, and Class III shock features, respectively (see Table.

CONSTRUCTION OF PSEUDO-SLITS AND SHOCK IDENTIFICATION
Shocks on the leading side of the bar are commonly found in simulations (Roberts et al. 1979;Athanassoula 1992) and observations (Sandage 1961;Sandage & Bedke 1994). The gas velocity component perpendicular to shocks has an abrupt change, producing sharp velocity jump features on the PVDs. In Fig. 2, we plot the schematic view of the possible shocks under the viewing angles of M31. The position angle of the projected bar major axis (on the sky frame) from Blaña Díaz et al. (2018) is 55.7 • , deviating from PA disk = 38 • by 17.7 • . If M31 is a barred galaxy as suggested by Blaña Díaz et al. (2018), we expect that shocks appear on the leading side of the bar and extend roughly to the bar ends. Thus, we position pseudo-slits perpendicular to the disk major axis so that the slits cut through the shocks nearly perpendicularly. The standard width of the slit is chosen to be 1.2 which is about twice larger than the spatial resolution of H I (43.75 ). The slits cover the inner 20 where shocks are expected.
The large-scale bar-driven shocks generally produce sharp velocity jumps on PVDs, for which we try to search in this work. For each PVD, the goal is to identify the positions and amplitudes of shock features. Canny (1986) detected edges in signals by convolving data with a derivative of Gaussian, and it is now commonly used in 2D image edge detection.
We modify their algorithm to detect step-function shaped and δ−function shaped velocity jumps. The steps and instructions of our procedures are described as follows.
Requirement of sharpness. The amplitude of velocity jumps ∆V as well as the spatial extent ∆Y determines the sharpness of a shock feature. Simulations by  and Kim et al. (2012) have shown that typical bar-driven shocks have ∆V > 100 − 150 km s −1 within ∼ 100 − 200 pc (after projection 200 × cos 77 • = 45 pc). We expect the spatial extents of observed velocity jumps to be wider than those in simulations for several reasons: 1. Sharpest velocity jumps are expected when the slits cut shocks perpendicularly, which may not be the case shown in Fig. 2. 2. Resolution of the observation is usually lower than that of simulations. 3. Dust extinction and local turbulence could blur the velocity field and produce less clear shock features. Considering these effects, we aim to find shock features showing ∆V > 100 km s −1 , ∆Y < 0.6 on PVDs, which corresponds to a velocity gradient over 100 km s −1 /0.6 = 167 km s −1 arcmin −1 . We use a larger ∆Y of 2 to identify H I shock features, which is around three times the spatial resolution (43.75 ) of H I. Compared to those of [O III], the shock features of H I have a smaller velocity gradient due to lower spatial resolution.
Smoothing. We use boxcar smoothing to reduce noises on PVDs. Each observed data point is replaced by the median * of its adjacent 15 points for PVDs of [O III]. Our tests show that the number of points of 15 is good enough for show- * We also tested that using median gives sharper shock features than using mean, and the results of the flux density-weighted mean do not differ much from the non-weighted ones.  Locating a jump feature. Given a curve showing several velocity jumps, we would like to detect jump features automatically. Our procedure first creates a window of 0.6 (2 for HI) wide and use it to cover a selected part of the curve. Then it convolves the data within the window using a derivative of Gaussian to highlight the position of the velocity jump, sim-ilar to the edge-detection algorithm in Canny (1986). The derivative of the Gaussian operator has the form: Here s = 3 is a scaling factor and σ = 4 represents the dispersion of the Gaussian. Vertical lines in Figs. 3 and 4 indicate the positions of selected shock features using this method. Tests of the edge-detection algorithm in identifying step-function shaped and δ-function shaped jump features are given in Appendix A. Identification. We calculate the velocity difference within 0.6 (2 for H I) for each point on the boxcar smoothed curve. Since we want to find sharp velocity jumps, we focus on regions showing velocity differences over 100 km s −1 . Our algorithm targets the region with the largest velocity difference, and returns one shock position using the method above. Then it moves to regions with smaller velocity differences, which departs from the identified shock features by at least 1.5 . We repeat this process several times until all of the velocity jumps over 100 km s −1 are found. The criteria to identify shock features is summarized in Table 1. Classification of shock features. We classify the gas velocity jump features into three classes according to the likeliness that a bar-driven shock is present based on ∆V .
− Class I: if the feature shows a velocity jump over 170 km s −1 , it is classified as a gas feature "very likely" being a shock.
− Class II: if the feature shows a velocity jump between 125 km s −1 and 170 km s −1 , it is classified as a gas feature "likely" being a shock.
− Class III: if the feature shows a velocity jump between 100 km s −1 and 125 km s −1 , it is classified as a feature "possibly" being a shock.
The ∆V values in Table 1 are empirically chosen. Previous observations and simulations of barred galaxies give a rough reference value of shock velocity jump ∆V 100 km s −1 (e.g. , which we choose to be the threshold ∆V of shock features. We have tested that minor variations of the ∆V ranges do not change the main result. We expect that positions of Class I and Class II shock features in [O III] 1 are similar to those of H I, except for central regions due to the lack of H I.

Shock features on PVDs in [OIII]
Figs. 3 and 4 present the PVDs of [O III] 1 and the boxcar smoothed result (black curve) on the receding (X < 0) and approaching (X > 0) side of M31. Negative and positive Y represent the far and near side of M31. Overall, shock features are clearer on the receding side (i.e. Fig. 3), showing Class I shock features (thick red lines) in most panels. Class I shock features appear in S0 at Y ≈ −2.7 and −0.8 , together with a Class III shock feature (red dashed line) at Y ≈ 2.9 (near side). In Fig. 3, the clearest shock features are shown in slits (S-1, S-2, S-3, S-4). Positions of Class I shock features shift downwards from Y ≈ −3.1 to −3.6 as we go from S-1 to S-6. Further out, ∆V of Class I shock features decreases and the shock features turn into Class II (thin red line) in S-7. On the approaching side, Class I and Class II shock features show up in most panels, mainly between Y ≈ 3.2 and 4.6 . Using the method described in §3, we extract the clearest shock features of [O III] 1 on the receding side, which is shown in Fig. 5. Color indicates the X positions of the slits. Each curve corresponds to one panel in Fig. 3. We do not show the shock features along S-6 and S-7 slits in Fig. 5 for they are too close to the boundary of data coverage. The shock features are found between Y ≈ −3.1 and −3.5 , and shift upwards by ≈ 150 km s −1 as the slit moves from S-1 to S-5.

Potential shock features in the HI data
We also try to cross-validate the shock features in [ Fig. 6. Positions of Class II shock features are found between Y = −2.7 and Y = −3.9 on the far side of M31. Further out, ∆V of Class II shock features decreases, and the shock features turn into Class III (blue dashed lines) in S-9, S-10. On the approaching side, a Class I shock feature (thick blue line) and several weak shock features show up in panels (S+3, S+5, S+8) of Fig. 7, but in other regions the gas features with large velocities are quite clumpy and do not show clear shock features. Further out, Class II and III shock features appear near the disk major axis ranging from S+9 to S+13.
We make a cross-comparison of the PVDs of H I and [O III] 1 . The identified shocks distribution of these two tracers are quite similar, hinting for a common origin that is probably due to large-scale bar dynamics. PVDs of S-3, S-4, S-5, S-6 in Fig. 3 Fig. 9 we show the same comparison but on a smaller spatial scale. Overall shock features of H I coincide with [O III] 1 , especially at X = −1.2 and X = −3.6 . The H I features are quite clumpy, therefore determining the exact shock positions of H I is not easy. It is possible that the shocks in H I are not shifted from those in [O III] 1 , but obscured or hidden by some missing clumps instead. For example, the shock profile appears to be incomplete at X = −2.4 . If there were a clump near Y ∼ −3.5 with a velocity around −100 km s −1 , the shock feature would have been clearer to see.  Shock features are, however, absent near X ∼ 5 . The observation of [O III] 1 mainly covers the bulge region with |Y | < 5 , but it covers only |Y | < 4 at X = 5 . Class I and Class II shock features in [O III] 1 are found at larger Y distances on the near side than the far side by ∼ 0.5 . It is possible that shock features do exist, but extend beyond the data coverage near X ∼ 5 .  We also make simple isothermal 2D gas models in a constrained M31 potential to compare with the observed shock features. The simulations here are for illustrative purposes only, but not meant to match all the details of shocks.
We solve Euler equations in the initial frame using the gridbased MHD code Athena++ (Stone et al. 2020). We adopt a uniform Cartesian grid with a resolution of 4096 × 4096 covering the simulation domain with length L = 15 kpc in each direction. The setting corresponds to a grid spacing of ∆x = ∆y = 7.3 pc.
For simplicity, we use the isothermal equation of state P = Σ c 2 s , here P, Σ denote the gas pressure and gas surface density, respectively. The effective sound speed c s describes the turbulent properties of the gas. Kim et al. (2012) systematically explores the effects of c s on gas substructures in barred potentials. They found that models with larger c s are We set the initial surface density distribution of the gas disk to an exponential profile Σ gas (R) = Σ 0 exp(−R/R d ), here Σ 0 = 76 M pc −2 and R d = 6.0 kpc, which results in a total mass of 1.2×10 10 M . The initial azimuthal rotation velocity of the gas disk is set to balance the azimuthally averaged gravitational force. We input linear growth of the nonaxisymmetric force and gradually ramp up the barred potential in one bar rotation period. We accomplish this by increasing the fraction of bar potential from 0 to 1.0 and decreasing the fraction of axisymmetrized bar potential from 1.0 to 0 linearly with time in T grow = 2π/Ω b , similar to previous studies (Athanassoula 1992;Kim et al. 2012, e.g.).

Gravitational potential
The galactic potential is based on the best fit m2m N-body model constructed by Blaña Díaz et al. (2018). The m2m  in the dynamical model consists of a classical bulge component with mass of 1.18 +0.06 −0.07 × 10 10 M and a BPB component with mass of 1.91 ± 0.06 × 10 10 M . The bar in their model has a pattern speed of Ω b = 40 ± 5 km s −1 kpc −1 and a length of semi-major axis of 4 kpc. The major axis of the bar is oriented at a position angle of 54.7 • (in the galactic plane) with respect to the line of nodes. The dark matter halo in their model follows an Einasto profile and the mass of dark matter within the bulge region (R < 3.2 kpc) is 1.2 +0.2 −0.4 × 10 10 M . The authors obtained similar mass values using models with NFW dark matter profiles. The massto-light ratio of stellar component in 3.6 µm is found to be Υ 3.6 = 0.72 ± 0.02 M L −1 . We use their best fit models JR804 and KR241 as the basis of our galactic potential. The former includes an Einasto dark matter halo and the latter includes an NFW dark matter halo. We add a Plummer sphere in the center to represent the supermassive black hole with a mass of M BH = 2 × 10 8 M (Bender et al. 2005): Here a = 10 pc.

Simulated bar-driven shocks
For simplicity, we only compare simulated bar-driven shocks with the clearest shock features of [O III] 1 on the far side (slits S-1, S-2, S-3, S-4 in Fig. 3). We discuss the possible reasons for weaker shocks on the near side in §6.1. We start with models using a bar pattern speed around Ω b = 40 km s −1 kpc −1 given in Blaña Díaz et al. (2018). After projection with an inclination of 77 • , the shocks are found to be too close to the bar major axis and cannot extend as far as shock positions in [O III] 1 . Other parameters that could be important are the uncertain shape and strength of the thin bar. Hammer et al. (2018) suggested that a single merger event ∼ 2 Gyr ago could explain the recent active star formation (Williams et al. 2015). More recent work studied the age-velocity dispersion relation in M31's stellar disk (Bhattacharya et al. 2019), which suggests a merger with a mass ratio of around 1:5 3-4 Gyr ago. The bar could have been weakened by the merger (Ghosh et al. 2021) and left its shape and strength uncertain. One way to produce shocks with a larger distance from the bar major axis is to decrease the bar pattern speed (Li et al. 2015). We tested models with lower bar pattern speeds and their shocks have a better match with [O III] 1 data. When the bar rotates with a lower pattern speed around Ω b = 20 km s −1 kpc −1 , the nuclear ring size becomes very large and the effects of thermal pressure need to be considered to reduce the size of the nuclear ring. It requires a high c s of around 30 km s −1 to make a reasonably sized nuclear ring. Models with low Ω b and high c s produce less curved shocks, larger velocity jumps, and a nuclear ring with a reasonable size. Using an Einasto or NFW dark matter halo does not affect much the substructures. We present the shock pattern in two models: (1) Model 1 with a bar pattern speed Ω b = 20 km s −1 kpc −1 and sound speed c s = 30 km s −1 based on JR804 potential; (2) Model 2 with a bar pattern speed Ω b = 33 km s −1 kpc −1 and sound speed of c s = 15 km s −1 based on KR241 potential. The bar rotation period 2π/Ω b for Model 1 and Model 2 is 307 Myr and 186 Myr, respectively. The ratio of co-rotation radius to semi-major axis of the bar is R ≡ R CR /a bar ∼ 3.2 and R ∼ 1.8 for Model 1 and Model 2, respectively. Although simulations have reached quasi-steady after two bar rotation periods, small transient changes still appear on PVDs. We choose snapshots such that the shock features of the models are most similar to those of [O III] 1 .
The x-axis of the simulation grid is along the major axis of the disk. The major axis of the bar in Model 1 is positioned at an angle of 54.7 • to the x-axis in the face-on case. The upper panel of Fig. 11 illustrates the gas surface density of Model 1 at T = 749 Myr projected with an inclination of 77 • . Pink circles and brown triangles represent the shock positions of [O III] and H I, respectively. Large, medium-sized, and small markers represent the Class I, Class II, and Class III shock features. We present the face-on view of Model 1 in the lower panel of Fig. 11. In Fig. 12 Fig. 13 that a non-rotating bar (i.e. similar to a triaxial bulge) cannot reproduce the shock features. The black dashed lines in Fig. 13 represent the PVDs of a model that has the exact same setups with Model 1 but with zero bar pattern speed. It is clear that no sharp velocity jumps can be found in this non-rotating bar model.
Adjusting the inclination and the bar angle of models can fine-tune the match with the observed shock positions of [O III] 1 . As the inclination decreases, the angle of the bar (after projection) in respect of the line of nodes becomes larger, resulting in shock positions further away from the disk major axis. In the process of decreasing inclination, a smaller bar angle is needed to keep the shape of shocks.
We run Model 2 to test the effects of inclination and bar angle. Model 2 produces shocks much closer to the bar major axis than Model 1, which shows a different shock pattern from that in [O III] 1 . However, with a smaller inclination of 67 • and a bar angle of 50 • , shock positions in Model 2 can still be similar to those in [O III] 1 , especially on the far side. In the upper panel of Fig. 14 we present the gas surface density of Model 2 at T = 799 Myr projected with an inclination of 67 • . The lower panel of Fig. 14 shows the face-on view of Model 2. Fig. 15 illustrates the PVDs of Model 2 on the receding side and its comparison with [O III] 1 and H I. The average shock velocity jump ∆V is smaller than [O III] 1 by ≈ 30%. These results show that there may be a degeneracy between the inclination and the bar pattern speed. We also present the PVDs in a non-rotating Model 2 that has zero bar pattern speed in Fig. 16. There are no shocks in this model   also tested different pattern speeds within the range of 25 -50 km s −1 kpc −1 and sound speeds within the range of 1 -50 km s −1 . Shock positions move closer to the x-axis with a larger bar pattern speed and/or a larger sound speed. The line-of-sight velocity jump at shocks becomes larger as bar pattern speed decreases and sound speed increases. A detailed comparison including more bar parameters will be presented in a follow-up paper. Berman (2001) and Berman & Loinard (2002) performed hydrodynamical simulations in a M31 potential and found a much higher bar pattern speed of Ω b = 53.7 km s −1 by fitting the line-of-sight velocity of CO (Loinard et al. 1995(Loinard et al. , 1999 along the disk major axis. They adopted a simple analytical bar model, and the spatial resolution of their simulations We find that the shock features are generally weaker on the near side than on the far side. In Fig. 4 it is clear that shock features are absent near Y ∼ 5 in slit S+4. In Fig. 17 we overlay the positions of the [O III] 1 shocks on the observed 2D velocity map. Class I shock features (large crosses) are found mainly at −7.2 < X < 0 on the far side (bottom left). On the near side (top right), there are mostly Class II shock features (medium-sized crosses) in regions with small velocity gradients. Overall shock features are closer to the boundary of data coverage on the near side than the far side. This implies the shock features of [O III] 1 defined in the current study might be incomplete due to the limited coverage of observation.
Clear asymmetry in the position of shocks and gas kinematics has been observed in the barred galaxy NGC 1365 (Zánmar Sánchez et al. 2008). They suggest that the asymmetry in positions of the dust lane is likely caused by a minor merger event. They also provided an alternative explanation that the ram pressure of a gas stream has moved the shock to offset its original position.
Another possible scenario for the asymmetry in the shock feature is that the warp in the outer H I disk has a stronger extinction effect towards the near side. The H I disk is extended and has a large-scale warp in the outer region (Newton & Emerson 1977;Henderson 1979;Brinks & Burton 1984;Chemin et al. 2009;Corbelli et al. 2010). Fig. 18 presents a schematic diagram of the gas warp structure. The gas warp appears as a foreground on the near side and becomes a background on the far side. The asymmetric extinction by the gas warp could result in shock features different between the near side and the far side. Furthermore, abundant dust in the H I warp may also help blur the [O III] 1 shock features on the near side. The distribution of dust is found to follow that of H I in the outer disk, and the strong reddening effect suggests a large amount of dust (Cuillandre et al. 2001;Bernard et al. 2012). More recent detection of dust by Ruoyi & Haibo (2020) found that dust in the disk has an exponential distribution and extends over 2.5 times its optical radius (around 54 kpc).

Effects of changing slit orientations on shock features
Maximum velocity jumps may be expected when slits are positioned perpendicular to the shock fronts (Athanassoula 1992). We checked PVDs of [O III] 1 and H I in slits at different orientations and compared them to Fig. 2 to see if the velocity jump features can be shown even more clearly. Our findings are roughly consistent with the theoretical expectation. When the slits are positioned parallel to the disk major axis, the shock features still appear, but with a smoother profile. The shock features become sharper as we reposition slits to be nearly perpendicular to the shock fronts, with an increase of ∆V around 30 km s −1 compared to the more parallel slits. In appendix B we tested repositioning the slit to make it more perpendicular to the shock fronts and found that the shock properties are similar to §4.1.

Comparison of shock features with dust and CO morphology
In Fig. 19 we present the map of our identified shock features in [O III] 1 (red circles) and H I (blue triangles) on a background of dust (color-coded with surface density) and CO (white contours) morphology. We use data of dust from Draine et al. (2014) and CO from Nieten et al. (2006) Such coincidence also appears in the CO distribution, which generally matches the dust morphology. The shock features of H I, extending to the bar ends in Blaña Díaz et al. (2018), connect to a possible spiral structure. For the near side, shock features are weaker as we discussed in §6.1. There are not enough Class II shock features to show the relation between shock positions and dust morphology. Another reason might be that the amount of H I, CO, and dust in the bulge is scarce on the near side. Near center there is a lack of H I shock features and CO emission. This could be due to the low gas density in the central region (Brinks & Shane 1984;Li et al. 2009;Dong et al. 2016;Li et al. 2019Li et al. , 2020.

Comparison of shock features in other nearby barred galaxies
Other barred galaxies also show qualitatively similar barinduced shock features as in M31. Mundell & Shone (1999) detected velocity jumps of H I with ∆V ≈ 130 km s −1 on the leading side of bar of NGC 4151. H I in NGC 4151 is not interfered by the gas warp, showing a clearer view of shock features. The Hα velocity field of the barred galaxy NGC 4123 illustrates similar shock features with ∆V > 100 km s −1 (Weiner et al. 2001b). Velocity jumps in NGC 4151 and NGC 4123 have a smaller amplitude than our Class I shock features, probably due to the lower inclinations and lower masses of galaxies. A large velocity gradient near dust lanes has been observed in other barred galaxies as well, e.g. NGC 1530 (Reynaud & Downes 1998;Zurita et al. 2004), NGC 7479 (Laine et al. 1999), NGC 5448 (Fathi et al. 2005), NGC 1365 (Zánmar Sánchez et al. 2008). The recent highresolution PHANGS-ALMA survey presents 2D gas kinematics of CO for nearby spiral galaxies (Leroy et al. 2021). The velocity gradient near dust lanes can even be visually identified in many barred galaxies from their sample, e.g. NGC 2903, NGC 3627, NGC 4536 and NGC 4945. Apart from bar-driven shocks, spiral arms in non-barred galaxies could also produce velocity jumps, but with a much smaller amplitude (usually around 40 km s −1 in simulations Roberts 1969;Pettitt et al. 2020 Opitsch et al. (2018). The gas flow forms streams and ring structures in barred potentials (Kim et al. 2012). Multiple gas components could be found when the line of sight passes through several of such gas substructures. A collision between M31 and its satellite galaxy M32 (Block et al. 2006) could also produce gas streams and rings, leading to multiple observed gas components along the line of sight. It is also possible that part of [ (2018). The spatial location of the shocks and the amplitude of shock velocity jumps are qualitatively consistent with our preliminary simulations of bar-induced gas inflow in M31. This result provides independent, strong evidence that M31 hosts a large bar. A detailed comparison with more hydrodynamical simulations will be presented in a follow-up study to provide a better understanding of the gas features in the center of M31, and hopefully determine better the main bar parameters of M31.
Software: Athena++ (Stone et al. 2020), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007 Our procedure mainly follows the edge-detection algorithm in Canny (1986), modified to detect step-function shaped and δ-function shaped velocity jumps on PVDs. The main idea is to convolve data with the derivative of Gaussian operator Here s = 3 is a scaling factor and σ = 4 represents the dispersion of Gaussian. The maximum (or minimum) of the convolution result of the data and the derivative of Gaussian operator highlights the position of the velocity jump features. In Figs. 20 and 21 we test the derivative of Gaussian operator with different σ to identify jump features in cases of four different amplitude of noises. f represents the ratio of the Gaussian noise amplitude to the amplitude of velocity jumps. Here we test the cases of f = 0.5, 1, 1.5 and 2. Dashed lines indicate the identified position of jump features. In general operators with a small σ tend to highlight more local fluctuations of the sample (more peaks on the convolution result). For stepfunction shaped jump features, the edge-detection algorithm works well for all σ values. A small shift of the identified position is found for larger σ at small noises. However, for noisier data it seems larger σ have better performance. For the δ−function shaped jump features a smaller σ is preferred for it gives more precise positions of jump features. As the noise amplitude increases, the peaks caused by noises become larger for σ = 2 and one could mistake them as the jump positions. For a balance between the stability and accuracy of edge-detection, we adopt σ = 4 in this work.