© The Authors 2025

1. Introduction

Andromeda is an SA(s)b galaxy characterized by an unusual morphology (Arp 1964; Haas et al. 1998; Helfer et al. 2003). It exhibits weak spiral structures (e.g., Nieten et al. 2006; Richards et al. 2006), scarce star formation SFR ≃ 0.4 M_⊙ yr⁻¹ concentrated at 10 kiloparsecs (kpc) in the main disk (e.g., Barmby et al. 2006; Braun et al. 2009; Chemin et al. 2009; Tabatabaei & Berkhuijsen 2010; Azimlu et al. 2011), and two rings at 1 kpc and 10 kpc, as observed in the gas, dust, and star formation distributions (e.g., Barmby et al. 2006; Nieten et al. 2006; Tabatabaei & Berkhuijsen 2010; Azimlu et al. 2011). It also features a central gas hole, long suspected to be of dynamical origin (Josey & Arimoto 1992). In addition, its non-symmetry observed at all scales both in the gas and stellar components (e.g., Pellet 1976; Josey & Arimoto 1992) is usually understood as being due to an intense collisional past, as supported by the presence of numerous relics observed as spectacular giant stellar loops and tidal streams in the outskirts (e.g., Ibata et al. 2005; McConnachie et al. 2009, 2018).

Relying on the spectacular dust rings and spiral arms observed in the mid-infrared with the Infrared Array Camera (IRAC) on board the Spitzer Space Telescope (Barmby et al. 2006), Block et al. (2006) stressed the presence of an elongated off-center inner ring with projected diameters 1.5 kpc by 1 kpc. Given the fact that both rings at 1 kpc and 10 kpc are off-centered, Block et al. (2006) argued that the most likely scenario for the formation of these rings is a head-on collision of the Cartwheel type (Struck-Marcell & Higdon 1993; Horellou & Combes 2001). Unlike the Cartwheel, where the companion is about one-third of the mass of the target (major merger), in Andromeda, the collision can be called a minor merger, producing much less contrasted rings in the main disk. Block et al. (2006) proposed that the collision partner was a M-32 like galaxy, comprising about one-tenth of the mass (dark matter included) at the beginning. After stripping experienced in the collision, the M 32 mass is now 1/23 that of the main target, M 31. The M 32-like plunging head-on with an impact parameter of 4 kpc, at a relative velocity of 265 km s⁻¹ 210 Myr ago, has triggered the propagation of an annular wave, now identified with the 10 kpc ring; there is also a second wave propagating more slowly behind it (e.g., Appleton & Struck-Marcell 1996), which would correspond to the inner ring. In addition, the inner ring has formed and is propagating in a tilted and warped disk, which accounts for its almost face-on appearance, in contrast with the inclined main disk of M 31 (Jacoby et al. 1985). This scenario explains why the cold gas has been expelled from the central region and also why there are shocks and hot gas present in the inner kpc. We must note that in addition to the large inclination of the main disk (e.g., Ma 2001), there are various shells connected to the large scale structures superposed on this two-ring morphology (e.g., Escala et al. 2022), adding to the complexity.

A stellar bar, or more exactly a triaxial bulge, has long been detected in the center of M31 (e.g., Stark 1977; Athanassoula & Beaton 2006; Beaton et al. 2007), in particular, due to the isophote twist. According to Saglia et al. (2018), there is an axisymmetric classical bulge, with old age (> 10 Gyr) and high metallicity in the M31 center, showing a negative [Z/H] gradient. No signature of the bar or boxy bulge is seen in the age or [α/Fe] maps, which are approximately axisymmetric. However, in Fig. 14 of Saglia et al. (2018), we can see that the metallicity of the old stellar population is enhanced along the bar and boxy bulge, as seen in other barred galaxies. Blaña Díaz et al. (2017) showed that a triaxial bulge with zero pattern speed does not reproduce the measured stellar kinematics. A box-peanut component with a non-zero pattern speed must be the result of the buckling instability of a (previously thin) bar, although there is not presently any thin bar outside the boxy-bulge radius. The near-infrared isophotes outside this radius align onto the major axis of the galaxy. A more detailed discussion of a possible bar in the disk of M31 is reported in Sect. 6.

In addition to being a quiescent galaxy with little star formation, M31 has an ultra-weak nuclear activity (del Burgo et al. 2000). Melchior et al. (2000) first detected only a small amount of molecular gas (1.5 × 10⁴ M_⊙) within 1.3′ (305 pc in projection). Relying on NOEMA interferometric observations, Dassa-Terrier et al. (2019) detected small molecular gas clumps within 250 pc, corresponding to 8.4 ± 0.4 × 10⁴ M_⊙ (see also Melchior & Combes 2017). The molecular gas is very clumpy and no ongoing star formation has been detected in this central region (e.g., Olsen et al. 2006; Li et al. 2009; Azimlu et al. 2011). In Melchior & Combes (2011), we detected molecular gas with IRAM-30 m with very large line splittings up to 260 km s⁻¹ in the northwest (NW) side of the disk along the minor axis. Towards some lines of sight, there were even three distinct velocity components. We discussed the fact that these apparently counter-rotating components appeared to be compatible with the Block et al. (2006) scenario, implying an inner ring with a tilted inner disk superposed on the main disk. In this framework, the observed narrow component was interpreted as the ring, the second wider peak was associated with the inner disc, and the third weak mid-wide component at the systemic velocity was thought to be associated with the main disk seen in projection. In Melchior & Combes (2016), were presented new molecular observations of dense gas on both sides of the minor axis. This further supported the previous scenario, as the narrow and wide components were at opposite velocity on both sides of the major axis.

Some ionized gas has also been widely detected in the central field (e.g., Rubin et al. 1971; Ciardullo et al. 1988; Boulesteix et al. 1987; Bogdán & Gilfanov 2008; Liu et al. 2010). This ionized gas (approximately 1500 M_⊙ according to Jacoby et al. 1985) can be accounted for by mass lost from evolving stars. Recent works discussed a possible very weak star formation activity (Leahy et al. 2022) compatible with the clumpy molecular gas distribution.

In the present article, we further investigate this central kpc field, relying on different observations achieving some kinematic estimates; namely, optical ionized gas spectroscopy with the imaging Fourier transform spectrometer SITELLE at CFHT and molecular gas observations with the IRAM-30 m telescope. In a first step, we present a dynamical modeling based on tilted rings embedded in a global potential, mainly due to stars, as the gas mass is very weak in this region (∼1500 M_⊙). We relied on a kinematic map extracted from SITELLE observations with the SN3 filter targeting the Hα and [N II] lines (Martin et al. 2018). In a second paper, we will present the full exploitation of the three SITELLE data cubes with a proper stellar continuum subtraction with Nburst (Chilingarian et al. 2007) and excitation diagrams.

In the model, we follow the disk tearing methodology proposed by Raj et al. (2021) to account for misaligned disks spinning around galaxy nuclei. These authors showed that in the settling process to align angular momenta, thin and highly inclined precessing disks are more likely to yield large warp amplitudes, while low-viscosity disks are generally expected to become unstable and break in one or two rings. Although the origin of the torques are different, this scheme is quite adapted to the characteristics of the gas in the central region of the highly inclined disk of M31. In Sect. 2, we describe the data we obtained with SITELLE at CFHT, and we extract in Sect. 3 a perturbed kinematic map derived from the optical ionized gas ([N II], Hα, [S II]). In Sect. 4, we present the characteristics of the modeling that we developed to reproduces the velocity field. The main features can be well reproduced with a customized static modeling of a central tilted disk extending to a warp region connecting the main disk and an independent offset 1-kpc ring. In Sect. 5, we discuss the main results achieved, that is how the proposed modeling reproduce the data. In Sect. 6, we compare our results to previous works, in particular, focusing on the proposition that these features correspond to a bar. Our conclusions are presented in Sect. 7. Throughout this paper, we adopt a distance to M31 of 780 kpc, thus, 1 arcsec = 3.8 pc.

2. Data

Here, we present the observational data we used to constrain the dynamical modeling. In Sect. 2.1, we discuss the archival data from molecular gas observations at IRAM-30 m. In Sect. 2.2, we recall the characteristics of the SITELLE data cube used to extract the kinematics. In Sect. 2.3, we present the 2D maps from archival data, namely, the dust emission and upper limits on dust extinction.

2.1. Archival molecular gas data

To make a comparison with the molecular gas component, we again reduced the CO(2–1) data from IRAM-30 m, taken with the HEterodyne Receiver Array (HERA) focal plane array, between November 2011 and March 2012 (Melchior & Combes 2013). The spatial resolution is 11″ and spectral resolution 2.6 km s⁻¹, we refer to the above paper for more details. We have built a data cube, with pixel size of 5.33″, and 13 km/s. The field of view of the cube is 7.2 × 7.9′, corresponding to 1.63 × 1.79 kpc.

2.2. SITELLE SN3 data cube

We rely on a data cube obtained with SITELLE and installed at CFHT¹. The observations were performed on August 24, 2016, with the SN3 filter designed to detect the Hα-6563Å Balmer line and the [N II]-6548,6583 Å and [S II]-6716,6731 doublets. The 11′×11′ field of view has been centered on M31 optical nucleus (J2000: RA: 00h42m44.37s, Dec: 41° 16′08.34″; Crane et al. 1992). This field has been integrated for 4.1 hours with a pixel size of 0.321″. Such an interferometric technique enables us to reach a good spectral resolution (here R = 4800), but it suffers from the multiplexing disadvantage (Maillard et al. 2013): all the photons collected by each interferogram are redistributed over all the channels of each spectrum. This effect is important for the central region of M31, but has been limited with the use of the SN3 filter. The reduction of the data cube has been described in Martin et al. (2018). This region devoid of star-forming regions is dominated by diffuse ionised gas (DIG), as discussed in Sect. 3.

2.3. Archival dust extinction and emission maps

As displayed in Fig. 1, the central kpc field of view of Andromeda has been mapped in dust emission by Block et al. (2006). The map in the left corresponds to the 8 μm Spitzer map with the stellar continuum at 3.6 μm (Block et al. 2006) has been subtracted. It has been overlaid with contours corresponding to the Hα map (discussed in the next section). On the right panel, the A_B extinction map, computed assuming the dust is seen in front of the bulge, is presented, as discussed in Melchior et al. (2000). The apparent extinction is larger on the near part of the disk (NW), which explains why the Hα emission is weaker in the northwest region.

Fig. 1. Dust location in the central kpc field of view of Andromeda. Left panel: Dust emission map obtained from the 8 μm Spitzer map with a subtraction of the stellar continuum at 3.6 μm. (Block et al. 2006). The contours correspond to the ionized gas intensity displayed in the top left panel of Fig. 5. The dashed square corresponds to the inner ring region displayed in the right panel. Right panel: an upper limit estimate on the AB extinction (Melchior et al. 2000). The gray contours correspond to the dust emission at 8 μm displayed in the left panel. The black cross symbol represents the M31 center.

Both maps show the presence of an off-center inner ring at 1 kpc. These maps provide some constraints on the geometry of this ring detected in dust emission in the northwest part, while it is emitting in ionized gas more strongly in the southeast (SE). The overall dust and gas content in this region is weak, and we went on to assume that the gas is dominated by the stellar potential. In this paper, we further explore the modeling of this region with the support of the new SITELLE and new CO(2–1) velocity maps of gas observed along the minor axis, to test this scenario.

2.4. Off-centering of the multiple components

As discussed in Block et al. (2006), the inner and 10 kpc rings are off-centered. The galaxy itself is clearly perturbed and non symmetric. However, while the central velocity field looks a priori asymmetric (Fig. 5), the central circumnuclear disk is only slightly off-centered (Fig. 18). We assume that the disk and warp components are centered, but we do introduce extinction, depending on the depth of each component. As later discussed, this will introduce some asymmetry.

In parallel, the inner ring displayed in Fig. 1 is clearly off-centered, which makes it an improbable inner Lindblad resonant ring, in the case of a bar scenario. We rely on its apparent position in the Spitzer map (center, position angle, and axis-ratio) to constrain its position.

3. Extraction of the kinematics from the SITELLE SN3 data cube

While a complete analysis of the data cubes is in progress, we rely on a kinematic map extracted from the data for the purposes of this study, as follows. The idea is to optimize the extraction of the kinematic information and the small pixel size. Indeed, the brightness of the DIG is low when compared to the surface brightness of the old stellar population background in the center of M31. Even if its presence is detectable in some frames of the spectral cube, we must start by a careful removal of the other components of the background spectrum in order to reveal the full content of its emission. We can then detect its different components and estimate their flux as well as their velocity which can finally be used as an input for a fitting procedure. We describe these three steps with more details in the following.

In Sect. 3.1, we discuss how we modeled the stellar background. In Sect. 3.2, we discuss the automatic detection of the gas velocity relying on the five emission lines. In Sect. 3.3, we present the line fitting used to get the velocity dispersion map.

3.1. Background modeling: Sky lines, Galactic foreground, and stellar continuum

This background spectrum is indeed the superposition of the old stellar population continuum and the atmospheric OH lines, neither of which are trivial to estimate. As the bulge stellar emission is strong in the whole field of view, there is no region where the atmospheric emission can be directly estimated. We could nevertheless find two corners, northeast (NE) and southwest (SW), where the bulge emission (and especially its Hα absorption line) is small enough that its effect on the emission lines can be partially neglected. We combined these regions to obtain a spectrum which was fitted to compute a model of the foreground emission. This emission contains both the atmospheric sky lines and some Galactic emission in [N II] (6563, 6584), [S II] (6717, 6731) and Hα at a mean velocity of −38 ± 5 km/s. (see Fig. 2). The Galactic emission seems homogeneous both in terms of flux and velocity in the whole field of view so that once subtracted, it should have a negligible impact on our estimation of the emission of the DIG of Andromeda.

Fig. 2. Background subtraction. Top panel: Fit of the sky lines over a low brightness region of the field of view. The vertical gray lines indicate the position of the atmospheric sky lines. The five red vertical lines indicate the position of the foreground emission of the DIG of our own Galaxy in the [N II] (6563, 6584), [S II] (6717, 6731) and Hα lines at a velocity of −38 ± 5 km/s. Middle panel: Sky lines model and Galactic foreground model. During the fitting procedure the velocity of the sky lines and the velocity of the Galactic foreground were two independent parameters. Bottom panel: Analyzed spectrum and model of the background.

The old stellar population emission varies considerably in the field of view but at a sufficiently large scale that it can be approximated over large bins of data. We choose a binning surface of 200 × 200 pixels² because it is small enough to capture the variations of the old stellar population emission and large enough that the DIG emission varies sufficiently in velocity that a median spectrum taken over this region is not significantly contaminated by it (see Fig. 2). Another challenge comes from the fact that the emission level displays a strong pixel-to-pixel gradient, especially near the center of the galaxy. As the emission (and particularly the absorption lines) can be approximately considered to be proportional to the mean level of emission in the observed band, a naive mean spectrum of a large region would be dominated by the spectra of the brightest pixels. In order to get a more balanced estimation of the spectrum covering a binned region, the spectra were normalized by their mean emission level before being median-combined resulting in a normalized estimation of the background spectrum of the binned region. This last spectrum was then smoothed to remove the remaining emission-lines and obtain a normalized model of the galactic background in the 200 × 200 pixels² box (see Fig. 3). Using an old stellar population modeling instead of a simple smoothing is certainly a more precise approach, but the noise level of the individual spectra used for fitting (even binned 7 × 7) is high enough that we estimate the difference between both processes (smoothing vs. stellar population modeling) to be smaller than a few percent and thus negligible. The normalized model of the galactic background at a particular position is then linearly interpolated between the positions of the binned spectra and multiplied by the mean emission at the same position (see Fig. 2).

Fig. 3. Example of the old stellar population modeling in a 200 × 200 pixels2 box near the center of the galaxy. Top panel: Random spectra taken in the box to illustrate the strong brightness gradient. Bottom panel: Same spectra after a normalization by their mean flux. The resulting median spectrum of all the normalized spectra contained in the box is shown in blue, the smoothed version used as an estimation of the old stellar population is shown in red. The noise level of the individual 7 × 7 binned spectra used during the fitting procedure is shown as an orange envelope.

3.2. Automatic detection of the emission components of the DIG

The automatic detection of the different emission components of the DIG was achieved with the same algorithm as the one used in (Martin et al. 2021). This algorithm is based on the correlation of the spectra with a five emission lines comb ([N II] (6563, 6584), [S II] (6717, 6731), and Hα). We used it on 7 × 7 binned spectra and we obtained a map of the estimated velocity of the DIG at each binned position. In some cases, two velocity components are found along the same line of sight. As discussed in the next section, an example of such two components is displayed in the bottom panel of Fig. 4. The derived velocity map is displayed in the top right panel of Fig. 5, as well as the corresponding intensity (top left panel). Error bars on these velocities are estimated between 10 and 60 km s⁻¹, according to the flux signal-to-noise ratio. When two velocity components can be detected, we display their flux ratio, and velocity difference in the bottom panel of Fig. 5.

Fig. 4. Examples of emission lines fitting on background corrected spectra displaying only one (top panel) or two components (bottom panel). Top panel: Fit over a background subtracted spectrum revealing the pure emission of the DIG. Bottom panel: Spectrum displaying two resolved velocity components along the line of sight.

Fig. 5. Intensity and kinematic information extracted from the ionized gas with SITELLE SN3 data cube. Top panels: Ionized gas intensity (left) and the heliocentric velocity (right) corresponding to the main velocity component. The ionized gas intensity has a threshold of 6 × 10−18 erg cm−2 s−1 Å−1. The dashed squares correspond to the field of view displayed in the bottom panels. Bottom panels: Flux ratio ϕ2/ϕ1 (left) and velocity difference V2 − V1 (right) for the second kinematic component with respect to the first one. The contours correspond to the intensity levels (1.8 × 10−17, 3.6 × 10−17 erg cm−2 s−1 Å−1) of the flux map of the main component (top left). The black cross symbol representing the M31 optical center, 00h42m44.37s + 41d16m08.34 (Crane et al. 1992).

In parallel, the most remarkable features of the molecular gas observations are presented in Fig. 6. Two components are detected along the minor axis, where a single disk would be expected with a systemic velocity. In addition, one of the two has a large velocity dispersion, as expected for a large velocity gradient (typical of an inclined disk in a 12-arcsec beam). The second one is relatively weak and narrow. In Melchior & Combes (2011, 2013), we proposed that this narrow component could be accounted for by an offset ring as revealed by Block et al. (2006). On the northwest side, the broad (resp. narrow) component is blue-shifted (resp. red-shifted), while on the southeast side, it is the opposite.

Fig. 6. CO(2–1) spectra averaged over ∼60″ regions in diameter in the northwest (top) and the southeast (bottom) of the observed IRAM-30 m moment-zero map. The NW corresponds to (−160″,  120″) offsets, and the SE to (150,   − 110″) offsets relative to the M31 center RA = 00:42:44.35 Dec = 41:16:08.6. They correspond to the same regions as Fig. 7.

3.3. Line fit to extract the velocity dispersion

In addition to the positions where two ionized-gas components are clearly resolved along the line of sight and with comparable intensity (see bottom panel of Fig. 4), in most cases, the (possibly multiple) components could not be spectrally resolved, but the lines were broadened. We used a sincgauss model (Martin et al. 2016) to fit these lines and map their broadening (see Fig. 8). Note that this model works under the hypothesis that the observed emission-line is the convolution of a Gaussian (the real emission-line) and an instrumental line function (a sinc). In some cases, given the spectral resolution of our data, this hypothesis might be reasonable, but, when the real emission-line is the sum of two different velocity components, the measured broadening of the Gaussian model is a biased measure of the velocity separation of the emission-lines. Depending on the flux ratio of the two unresolved emission-lines, the measured broadening may be off by more than 30% from the real velocity separation (see Fig. 9). The measured broadening must thus be considered as an indication of the existence of two or more unresolved components along the line of sight giving only an approximate insight on their real velocity separation. It might however be useful when the S/R is low enough to prevent the modeling of two different components with unknown fluxes – which is the case in the present study.

Fig. 7. Hα spectra averaged over 12″ regions in the southeast (left) and the northwest (right) of the observed SITELLE map. The SE corresponds to (124″, −72″) offsets, and the NW to (−199, 72″) offsets relative to the M31 center RA = 00:42:44.35 Dec = 41:16:08.6.

Fig. 8. Velocity dispersion derived from the Gaussian fit on SITELLE sky-subtracted cube.

Fig. 9. Measured broadening when using a sincgauss model (broadened Gaussian emission-line) to fit two unresolved emission-lines. We have simulated the fit with a sincgauss model of two unresolved emission-lines with different flux ratios. We can clearly see that the measured broadening can be very different from the real velocity separation and depends strongly on the flux ratio of the unresolved lines. Note that, at the resolving power of our data (R = 4800), two components are considered resolved when separated by more than 63 km s−1. The dotted line is the bisectrix, meaning equality between the two quantities.

3.4. Velocity dispersion and double-peaked emission profiles

The previous kinematic map displayed in the right panel of Fig. 5 reveals a complex velocity field with discontinuities, which is capable of producing double-peaked features when smoothed in a larger beam. As displayed in Fig. 7, spectra integrated with a 12-arcsec beam on the SITELLE-SN3 data cube display double-peak features at positions along the minor axis and on the east side. Opitsch et al. (2018) also mapped two components in the velocity field. In CO observations along the minor axis (Melchior & Combes 2011, 2013), spectra exhibits also two components, that we interpreted as an inner ring (corresponding to the thin component) and a disk (corresponding to the thick component).

We also extracted the velocity dispersion from a single Sinc-Gaussian fit, as displayed in Fig. 8. For this figure, a threshold in flux has been required, with more than 3σ significance. The sensitivity is low in the very center due to the multiplex disadvantage. Interestingly, we can note that several regions exhibit a velocity dispersion larger than 100 km s⁻¹, in particular in the southwest region. However, taking into account the multiphase properties of this diffuse gas is beyond the goal of this paper.

We also observe, in the bottom panels of Fig. 5, that some positions present double peaks at the nominal resolution (1 arcsec). Interestingly, they lie at the edge of the densest regions, suggesting a warping, or a large inclination gradient along the line of sight. In the southern region, the velocity difference reaches 300 km s⁻¹ and seems typical of the southwest region where velocities range between −400 and −650 km s⁻¹, with patchy regions at different velocities.

3.5. Velocity map analysis

To extract geometrical constraints from the velocity map, we searched the maximal and minimal values of the velocity at different radii. We were then able to identify regions with possibly independent morphology, and to derive their respective position angle (PA =θ). In this Sect. 3.5.1, we present the methodology used to extract information from the map. In Sect. 3.5.2, we subsequently identify the main regions of interest (or main patterns of the velocity field).

3.5.1. Methodology of the map analysis

We first computed a set of N concentric circles of regularly spaced radii, expanding from the center of M31. The radius r_i (in kpc) of each circle is defined as $r_i=i\times \frac{R_{max}-R_{min}}{N}$. On each of these circles, we selected the minimal and maximal values V_min and V_max of the velocity field displayed in Fig. 5. Figure 10 illustrates this process for N = 9 circles. Each circle i with a radius r_i (in kpc) is separated into $N_i=\frac{r_i}{10\mathrm{kpc}}\times 360$ arcs. Each arc then has a constant length of π/180 kpc corresponding to 17 pc.

Fig. 10. Simplified illustration of the map processing, on N = 9 concentric circles, with radii regularly spaced between Rmin = 0 and Rmax = 0.8 kpc. Top panel: Velocity map overlaid by the concentric circles (in gray dotted lines). The red dots (resp. blue triangles) are the maximal (resp. minimal) values Vmax (resp. Vmin) for each circle, and have been computed with n = 7. Bottom panel: Previous extrema Vmax and Vmin of the velocity displayed as a function of the distance from the center (i.e., the radius of its corresponding circle). The negative distances corresponds to the blue shifted part of the map, where the velocities are negative. The gray lines represent the origin axes.

We then defined a running box of 2n × 2n pixels centered on each arc. For each box, we define an averaged velocity $V_{\mathrm{ave}}=\frac{1}{4n^2}{\mathrm{\Sigma }}_{i=0}^{i=2n}{\mathrm{\Sigma }}_{j=0}^{j=2n}{V}_{\mathit{ij}}$, with n = 7. In parallel, we eliminated pixels with a value is smaller than 10% of the maximum value in a similar box with n = 10. Hence, for each circle, we kept the following boxes:

• Pixels with a defined flux (i.e., above a threshold that depends on its surroundings), and boxes with less than 50% undefined-flux pixels;

• The standard deviation of the velocity inside the n = 7 box is smaller than 150 km s⁻¹;

• The sign of the velocity is the same as the sign of the median velocity over the neighboring points.

These criteria enabled us to remove most artifacts and to smooth sufficiently to get reliable results. The minimal and maximal velocities V_min and V_max are thus computed for each ring. Figure 11 displays the resulting map produced with N = 300 concentric circles, R_max = 1 kpc and R_min = 0.

Fig. 11. Minimal and maximal values, Vmin, $V_{min}^{2\mathrm{nd}}$, and Vmax, of the velocity on N = 300 concentric circles on the map, centered on the BH, with radii regularly spaced from Rmin = 0 to Rmax = 1 kpc. Left panel: Velocity map overlaid by the extreme points of the velocity. The red dots (resp. blue triangles) are the maximal (resp. minimal) values at each radii. As discussed in the text, the green points correspond to the $V_{min}^{2\mathrm{nd}}$ velocity found when considering the second component detected and presented in Fig. 5. Right panel: Minimal and maximal values of the velocity, as a function of the distance from the BH. The negative distances corresponds to the blue shifted part of the map.

The same analysis was conducted on the second velocity component, defined with a flux ratio above 30% (as displayed in the bottom panels of Fig. 5). Only the minimal points are computed for this second component, as the maxima are in patchy and noisy parts of the map. We superpose these secondary velocities $V_{min}^{2\mathrm{nd}}$ as green points in Fig. 11.

3.5.2. Separate regions of the velocity map

The extreme values of the velocity as a function of the galactocentric distance reveal well-defined independent structures in the velocity map (see Fig. 11). Firstly, this analysis stresses the very asymmetrical nature of the gas component of M31, detected both in the kinematics and in the structural characteristics. This asymmetry is mainly noticeable on both sides of the minor axis of the kinematics, as described in the following.

3.5.2.1. Northeast region.

This region hosts the maximal values of the velocity, which follow an almost-continuous line from the center to the edge of the map. The lines follow the winding visible at the center of the map. Around 0.4 kpc, the winding abruptly stops and the line roughly follows the axis of the position angle of Andromeda, around PA = 35° from the northern axis (cf. Hyper-Linked Extragalactic Databases and Archives (HyperLEDA)). Such a structure is compatible with the signature of a rotating gaseous disk, whose position angle varies with the galactocentric distance. This variation induces a warped gaseous disk, similar to that depicted in Raj et al. (2021). Moreover, the value of the velocity as a function of the galactocentric distance (right panel of Fig. 11) suggests a distorted rotation curve, hinting toward a variation of the inclination of the gaseous disk as well. We can already note at this stage that the few red points found at the eastern-most part of the map are not compatible with a warped gaseous disk, and may indicate non-circular motions.

3.5.2.2. Southwest region.

On the contrary, this region hosts several distinct structures rather than an almost continuous line. One of those structures, near the center, is the counterpart of the winding seen in the northeast region as shown in Fig. 12. As seen in Fig. 13, the position angle of the red and blue symbols in the central 350 pc radius are roughly symmetrical. Moreover, the values of the velocity for the first and second velocity component (blue and green triangles) are in the same order of magnitude as their northeast counterpart, with a probable asymmetrical warp of the gaseous disk in this region.

Fig. 12. Zoom-in on the winding in the very center of the map. The analysis is conducted on N = 100 radii from Rmin = 0 to Rmax = 0.35 kpc. Top panel: Vmin (blue), $V_{min}^{2\mathrm{nd}}$, and Vmax velocities superposed on the kinematic map, together with the maximal velocity (red bullets), minimal velocity (blue triangle) and minimal velocity of the second component (green triangle), for each radii. Bottom panel: Corresponding velocity values as a function of the galactocentric radius. The negative radii correspond to the distance of the points in the southwest region.

Fig. 13. Position angle (PA) of the main kinematic features as a function of the galactocentric center, computed for the maximal values of the velocity in the northeast region (red) and the southwest region (blue). We distinguish the winding below 300 pc. Beyond 300 pc, the position angle slowly tends toward approximately 40°.

Besides the winding pattern visible for galactocentric distances below 400 pc, the distribution of the minimal values on the map shows small segments, but no long and coherent structure emerges, neither in the primary nor in the secondary velocity component. This suggests disk tearing rather than disk warping, inducing discontinuity in the geometrical morphology of the gas in the southwest region. Beside the rift of the tearing, we expect double velocity components, as we would have the two gas structures that were teared apart really close to one another on the line of sight. This prediction is compatible with two bottom panels of Fig. 5, which shows a small crescent of high velocity difference (−300 km/s) where the rift would be. This could also account for the perturbed velocity dispersion pattern displayed in Fig. 8.

3.5.3. Constraints on the gas distribution

The distribution of the extreme values extracted from the kinematic maps at different galactocentric distances suggests:

1. A highly asymmetrical structure,

2. A warped disk in the northeast region, with continuously varying position angle (and inclination),

3. A teared disk in the southwest region, with discontinuously varying position angles (and inclination).

In addition, we define the position angle (PA) as the angle between the northern axis and the direction of the maximal point from the center. Figure 13 displays this PA value as a function of the galactocentric distance. While the redshifted points follow relatively regular pattern, the blueshifted points are irregular but for the inner 300 pc. The value of the inclination angle is discussed in the subsequent section, as it depends on the gravitational potential.

4. Dynamical modeling

The goal of the present work is to model the ionized gas in the inner kpc of M31 and to reproduce the peculiar features displayed in its velocity field in the right panel of Fig. 5. Several disk components, with different inclinations and position angles will be considered, to take into account the warping and misalignment of the gas in the central perturbed kpc zone. We identify four peculiar regions of interest:

1. The winding of the velocity field at the very center of the map;

2. The discontinuity of the velocity in the eastern part of the map, where we can see negative (blue) velocities superposed with positive (red) velocities (see zone between points A and B in Fig. 11);

3. A crescent-shaped region of high velocity (and high velocity gradient) in the southwest region;

4. The main disk of M31, visible through projection in the northeast part of the map.

As the gas content in the central kpc is weak in molecular (Melchior & Combes 2011, 2013, 2016, 2017; Dassa-Terrier et al. 2019), in atomic (Braun et al. 2009), and in ionized gas (Jacoby et al. 1985), we assumed the gravitational potential to be dominated by stars. We study the kinematics of the structure via static simulations of gas particles, located in a quasi-spherical potential. The positions of the particles are given by the geometry of the gas component, while their velocity vectors are derived from the total potential. We assume circular orbits corrected by an asymmetric drift and we choose to adapt our modeling inspired by Raj et al. (2021). Their work presented a simulation of misaligned disks, characterized with different position angles and inclinations, spinning around the galaxy nucleus, and their torque-driven differential precession triggers warps in the inner galaxy. We thus consider the central part of M31 to be a set of misaligned circles, each one with its own position angle θ and inclination angle i. The angles θ and i depend only on the distance r from the center of the galaxy.

In this section, we describe the principle and the main features of the modeling. In Sect. 4.1, we define the gravitational potential in the central region of M31. Then, in Sect. 4.2, we derive the density and velocity fields. Finally, in Sect. 4.3, we describe the geometry of the gas in these fields. We chose to use a symmetric modeling, and to adjust separately the NE and SW regions of the map.

4.1. Potentials

The gravitational potential results from the superposition of the potentials produced by each massive structure composing the galaxy. We represent the spherical components (bulge, dark matter) by Plummer potentials, and the disks by Miyamoto-Nagai potentials (Miyamoto & Nagai 1975). Contrary to previous models of the stellar component out to 6 kpc in radius (e.g., Blaña Díaz et al. 2017, 2018; Saglia et al. 2018), our goal is not to obtain a very precise model of the stellar bulge, but to model only the gaseous disk within 1 kpc, in terms of several components, warped and tilted. Our main point in the present paper is to focus on modeling the complex gas structure; for that we need an analytical model of the potential to derive easily frequencies and asymmetric drift, but not a perfect bulge model.

The potential Φ_Plummer of a Plummer sphere of scale R, at a distance r from its center, is given by:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{Plummer}}(r)=-\frac{\mathit{GM}}{\sqrt{r^2+R^2}},\end{array}$$(1)

where M is the total mass of the component, and G the gravitational constant. For a disk of scale length, R, thickness, h, and total mass, M, the value of the Miyamoto-Nagai potential Φ_MN, in cylindrical coordinates (r, z), is:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{MN}}(r,z)=-\frac{\mathit{GM}}{\sqrt{r^2+{(\sqrt{z^2+h^2}+R)}^2}}·\end{array}$$(2)

As a result, the spherical components have two parameters: their mass M and scale length R, whereas the disks have three: their mass, M, scale length, R, and thickness, h. We take into account five structures of different sizes, shapes, and natures:

1. The dark matter halo, characterized by M_DM,  R_DM;

2. The main stellar disk, characterized by M_SD,  R_SD, and h_SD;

3. The gaseous disk, characterized by M_GD,  R_GD, and h_GD;

4. A classical spherical galactic bulge, characterized by M_B,  R_B;

5. A nucleus, characterized by M_N,  R_N.

The total potential is the sum of these five potentials; hence, it is partly cylindrical and not completely spherical. However, the central kpc is dominated by the bulge, so the total potential can be considered almost spherical. This justifies our simplified assumption of circular orbits for the particles, depending only on the distance from the galactic center. Also, we did not take into account the inclination of the disk components when warped or tilted to compute their velocities, since the stellar disk brings a perturbation of less than 4% on the rotational velocity (within the central kpc, as can be seen in Fig. 19). The nucleus component accounts for the mass of the three stellar disks discussed in Bender et al. (2005) as well as the black hole.

4.2. Density and velocity fields

To generate particles to model the central region of M31, we define the density and velocity fields. We assume for the gas the density ρ(M) corresponding to a Miyamoto-Nagai disk, and we distribute 2.4 × 10⁶ gas particles inside a peculiar disk, whose thickness H depends on the radius r, following:

$$\begin{array}{c}\hfill H(r)=\alpha \times r+\beta ,\end{array}$$(3)

where α is fixed at 6%, and β is the thickness of the disk in the galactic center, fixed at 30 pc. The radius of the disk in which the particles are generated is 5 kpc, and its maximal thickness is 330 pc (see Fig. 14).

Fig. 14. Cross-section of the gaseous disk along a diameter. The crossing of the r and z axis represents the center of the galaxy.

Each of the gas particle is given a velocity vector, whose norm is derived from the total potential Φ as followed:

$$\begin{array}{c}\hfill ||\overrightarrow{v}||(r)=\sqrt{r\frac{d\mathrm{\Phi }}{\mathit{dr}}}·\end{array}$$(4)

We assumed circular orbits, with some epicyclic perturbations, accounting for the velocity dispersion, and therefore add an asymmetric drift. We adopt a Toomre parameter of Q = 1.3, as discussed in Melchior & Combes (2011), and compute the radial, and azimuthal dispersions according to the epicyclic approximation, and the z dispersion as a function of the plane thickness. To sum up, we have 2.4 million gas particles, each characterized by a coordinate vector (x, y, z, v_x, v_y, v_z) in phase space.

The result of this distribution is a flared disk, whose velocity field would be a classical spider diagram. In this work, we assume that a warping and tilting of the disk can result in the peculiar velocity field we have observed. The mass of the gas component is known to be weak, which allows us to consider its morphological transformation as negligible with respect to the total potential.

4.3. Geometry

We modeled the gas in the central part of the galaxy as a continuous set of concentric rings, whose radius range is from 0 to 5 kpc, to take into account the high inclination of the main disk. We generated 2.4 million particles, characterized by a set of six parameters (x, y, z, v_x, v_y, v_z) inferred from the Eqs (3) and (4). To model a warped disk, each concentric circle is tilted and offset, to reproduce the kinematic geometry, according to adjustable parameters. Hence, a circle is characterized by a set of five parameters: r its radius, θ its position angle, i its inclination angle, δN its offset in the northern direction, and δE its offset in the eastern direction. Each of the 2.4 million particles and their velocity vector are rotated, then translated, to model the misalignment and the offset of the circles. After these transformations, the new coordinates are (x′,y′,z′,v′_x, v′_y, v′_z). The absolute value of the velocity vector (corresponding to the circular velocity) is not modified by the transformation. The galactocentric distance is modified only if there is an offset. The characteristic of the transformation depends only on the initial galactocentric distance: $r=\sqrt{x^2+y^2+z^2}$. We do not take into account the coordinate shifts to modify the corresponding velocities. This simplifying approximation is justified by the complexity of optimizing for too many parameters.

We decompose the gaseous disk between three substructures, independent from one another and adjustable separately: (1) a main disk, belonging to the galactic plane, centered on the black hole; (2) a warped nuclear disk, tilted with respect to the galactic plane, centered on the black hole; and (3) an inner ring, tilted with respect to the galactic plane, with an offset center. The nuclear disk presents a warped zone and is teared apart from the main disk. The disk tearing phenomenon happens with a geometry very similar to that described in Raj et al. (2021), but for different physical reasons. The ring is also teared apart from the remainder of the gaseous disk. Apart from the nuclear disk, we consider the substructures to be stiff: inside the ring or inside the main disk, the values of the position angle, inclination angle, and offset are constants. A simplified schematics of the geometric modeling can be found in Figs. 15 and 16.

Fig. 15. Schematic representation of the geometry of the gas in the central kpc of M31. The graph represents a simplified cross section of the gas disk in the sky plane containing the central black hole, represented by the black cross. The blue part is the main disk, aligned with the galactic disk of M31, the green part represents the ring and the red part shows the nuclear disk. The green cross represents the center of the ring, offset with respect to the nucleus. The dashed part of the red line corresponds to the warped disk, whereas the full lines are stiff components. The dashed black line points to where the disk is supposedly teared. We also show a number of parameters of the modeling: rw is the distance between the central black hole and the end of the warped region; rmd indicates where the main disk begins; rr is the inner radius of the ring, as well as the end of the nuclear disk as they are supposed to be from a single structure that teared; wr is the width of the ring; and rid is the radius of the stiff disk inside the inner region. The substructures are not to scale for visibility reasons.

Fig. 16. Schematic representation of the nuclear warped disk (corresponding to the red substructure in Fig. 15). The black star is the center of M31. The x axis represents the line of sight, and the yOz plane the plane of sky. The filled pink color refers to the inner nuclear disk, and also to the end of the warped region. The dashed lines represent two perpendicular diameters of the substructure. One of the two coincides to the major axis of the projection. The plain lines show where the second diameter would be if the substructure was not inclined (i.e., fully in the plane of sky). The successive circles in the middle are a simple representation of the warped region, defined as a set of misaligned circles with continuously varying inclination and position angles. For visibility reasons, the different parts of the substructure are not to scale with respect to each other, and only a few of the misaligned circles have been drawn. Bright red line at the outer bound of the substructure shows where the warped disk supposedly ripped from the 1 kpc ring. The nuclear warped disk is located inside the 1 kpc radius of M31, superposed with the ∼1 kpc inner ring and the main disk.

4.3.1. Main disk

This disk has a hole in the middle, where we only see the inner and nuclear disk. The radius of the hole r_md is an adjustable parameter of the modeling. The other parameters’ values are fixed at the main orientation of the M31 disk: θ_md = 37°, i_md = 77°, and δE = δN = 0. The thickness of the main disk follows the Eq. (3). Even though the truncation radius we use (of 5 kpc) for the galactic disk is way beyond the field of view (approximately 1 kpc in radius), the ionized gas can still be seen in projection in the central area.

4.3.2. Ring

The inner ring is supposed to coincide with the dust ring seen in Spitzer maps and identified by Block et al. (2006). We note that a schematic representation of this inner ring is displayed in Fig. 17. It is tilted and offset with respect to the main disk, thus, it can be defined by six parameters:

Fig. 17. Schematic representation of the main disk and the inner ring (corresponding to the blue and green substructures in Fig. 15). The stars represent the center of the ring (in green), and the center of M31 (in black). In our modeling, we neglect the component of the offset parallel to the line of sight, and only consider the offset toward the east direction (here $-\overrightarrow{u_y}$) and north direction (here $\overrightarrow{u_z}$). The x axis represents the line of sight, and the yOz plane the plane of sky. The green dotted and dashed axis is parallel to the z axis. In each color, the dashed lines represent two perpendicular diameters of the substructures. One of the two coincides to the major axis of the projection. The plain lines show where the second diameter would be if the substructure was not inclined (i.e., fully in the plane of sky). The angular parameters of the modeling are (ir, θr) and (imd, θmd); the inclination and position angles of the ring and of the main disk.

1. r_r the inner radius (see Fig. 15);

2. w_r the width (see Fig. 15);

3. θ_r the position angle;

4. i_r the inclination;

5. δN_r the offset of the ring center in the northern direction;

6. δE_r the offset of the ring center in the eastern direction.

4.3.3. Warped nuclear disk

The goal of the nuclear disk is to reproduce the winding at the very center of the velocity map, as presented in Fig. 18, as well as the discontinuity in the velocity profile as seen in the northeast part of the map (see Fig. 5). The varied geometrical angles are defined in Fig. 16.

Fig. 18. Central 75 × 75 arcsec2 (285 × 285 pc) field of view of Andromeda. The white cross corresponds to the optical center. The intensity contours are the same as in the bottom panels of Fig. 5.

5. Fitting the kinematics

The parameters for the modeling are adjusted to fit the velocity map shown in the top right panel of Fig. 5. The modeling can be decomposed between the gravitational potential induced by the massive components of the galaxy, and the morphology of the gaseous disk immersed within this potential. Within 1 kpc scale, we expect the potential to be dominated by the stars, mainly apportioned in the bulge, the main stellar disk, and the nucleus. The parameters of the potentials are needed to determine a theoretical rotation curve for the galaxy (Sect. 5.1), which will then be used to infer the geometrical parameters of the gas (Sect. 5.2).

5.1. Theoretical rotation curve

The theoretical rotation curve is entirely determined by the morphology of the massive structures (as explained in Sect. 4.1). Several previous works give a good approximation of the majority of these parameters. The main stellar disk parameters M_SD, R_SD, and h_SD are set according to Dalcanton et al. (2023). The dark matter halo parameters M_DM and R_DM are set according to Melchior & Combes (2011), as well as the parameters for the galactic bulge M_B and R_B. The nucleus parameters M_N and R_N are set according to Bender et al. (2005) and Bacon et al. (2001). Finally, we assume the mass of the gaseous disk to be 10% of the stellar disk mass, and to extend to a 2 kpc-radius, with a maximal thickness of 200 pc.

In many of these works, the modeling chosen for the disks were exponential disks. Using a 2D curve fitting, we converted the scale and mass parameters to fit a Miyamoto-Nagai disk. A summary of the chosen parameters is displayed in Table 1.

The theoretical rotation curve is computed using Eq. (4). The contribution of each component is represented on Fig. 19. We computed the asymmetric drift, the difference between the rotational velocity V_rot and the circular velocity V_circ, according to the following formula, adapted to the mid-plane of an axisymmetric rotating disk of density ρ(R):

Fig. 19. Theoretical rotation curve derived from Eq. (4), using parameters from Table 1. The contributions are from the nucleus (green), the bulge (blue), the stellar disk (orange), the gaseous disk (red), and the dark matter halo (purple). The black line is the total rotation curve. The blue triangles correspond to the measured ionized gas velocities, averaged over the NE and SW side, and corrected from asymmetric drift.

$$\begin{array}{c}\hfill V_{\mathrm{circ}}-V_{\mathrm{rot}}=\frac{{\sigma }_r^2}{2V_{\mathrm{circ}}}[\frac{{\sigma }_{\theta }^2}{{\sigma }_r^2}-1-\frac{dln\rho {\sigma }_r^2}{dlnR}]\end{array}$$

where the radial and azimuthal velocity dispersions σ_r and σ_θ are computed from the epicyclic approximation, and the Toomre criterion discussed in Sect. 4.2.

5.2. Geometrical constraints

The parameters of the main disk are fixed. We let the position angle and inclination of the inner ring to vary, and also the geometrical parameters of the nuclear warped disk. Since the southwest region is affected by clear velocity perturbations, plausibly due to a large-scale shock, the geometrical parameters were essentially fitted on the northeast region. A summary of the geometrical parameters can be found in Table 2.

5.3. Fitting results

Four different scenarios have been explored for the complex geometry of the central kpc in M31. They are schematically described in Fig. 20. They differ mainly by the degree of decoupling of the three different dynamical structures considered in projection in the central region. The corresponding geometrical parameters are listed in Tables 2 and 3.

Fig. 20. Four scenarios proposed for the central region in M31, represented by their schematic transversal cut. Grey lines correspond to transition zones, blue lines to the main disk, red lines to the warped nuclear disk, and green lines to the inner ring. When two structures are decoupled, two distinct velocity components may be observed along the same line of sight, while when they are linked by a transition zone, there would be a single wide velocity component. Top-left: Scenario 1, where the nuclear and main disks are linked. The inner ring is decoupled, implying two velocity components along a line of sight. Top-right: Scenario 2, where the nuclear disk is linked to the inner ring. The main disk is decoupled, with possibly two velocity components on some line of sights. Bottom-left: Scenario 3, where the three structures are decoupled, implying distinct velocity components along a line of sight. Bottom-right: Scenario 4, which is the same as scenario 3, but there is now a link between the main disk and inner ring, limiting the occurrence of distinct velocity components.

According to these parameters, we predict the velocity field in the four scenarios, represented on Fig. 21 and Fig. 22, for the cube moments 1 and 2 respectively. The fiducial model (scenario 3) is represented at the bottom left. In scenario 1 (top left), the nuclear disk has an enhanced warp. The inclination of the ring is lower in scenario 2 (top right) and the radius of the inner warped disk smaller in scenario 4 (bottom right).

Fig. 21. Velocity field obtained in the four scenarios described in Fig. 20, in the same order. The maps have been masked with the velocity map of Fig. 5; namely, we do not plot any velocity when there is no signal in the data.

Fig. 22. Same as Fig. 21 for the second moment of the cube. The maps have been masked with the velocity dispersion map of Fig. 8; that is, we do not plot any velocity dispersion when there is no signal in the data.

The velocity maps of Fig. 21 show that the velocity twist in the very center is reproduced in the four scenarios, thanks to the warped nuclear disk. Its amplitude is better fitted in scenarios 3 and 4. The inner ring reproduces the blue-shifted region in the SE, and the best fit is scenario 3. The main disk is responsible for the red-shifted region in NE, and is optimally fit in scenarios 2 and 3. The high velocity dispersion in regions NW and SW is not perfectly reproduced in either scenario in Fig. 22, but scenario 3 is approaching the closest. Globally, scenario 3 corresponds to the best qualitative model.

6. Discussion

The recent kinematical data obtained with Hα and [N II] lines on SITELLE in the central kpc of M31 have allowed us to build a fully sampled velocity field of the ionized gas. Combined with the more patchy molecular gas velocity field, obtained with the CO lines at IRAM-30 m telescope, we identify three dynamical components, the main disk, a tilted ring, and a nuclear warped disk. The tilted ring is imposed by the discovery of two widely different velocity components (by ∼260 km/s) along the same line of sight, toward the center, around the minor axis (Melchior & Combes 2011, 2016). To account for the twisted velocity field around the very center (within 200 pc in radius; see Fig. 5), we show that a warped nuclear disk, and then an off-centered ring propagating through a warped disk reproduce most of the observations. The signature of an almost face-on nuclear disk surrounded by a warped area is a good candidate to explain the winding found around the black hole of Andromeda. The extent of the winding, and its characteristic shape can be explained by the tearing of a warped zone between the tilted ring and the nuclear disk. The analysis of the velocity dispersion map shows that the highest dispersions are due to a superposition of distinct sub-structures, rather than a broadening of the velocity. By superposing the geometrical structures on the dispersion map, we see that a superposition of substructures induces dispersions between 80 km s⁻¹ up to 120 km s⁻¹, a warped zone can induce a dispersion between 40 and 80 km s⁻¹ where the warped region is edge on. The average dispersion for a simple disk structure is around 40 km s⁻¹. We can interpret in that way the observed velocity dispersion map, to deduce where we might find several velocities along the line of sight. This interpretation is corroborated by the double sinc fitting of the cube.

We note that the warped structure of M31 has been reported by many authors in the past, namely in the H I gas (e.g., Chemin et al. 2009) and references therein. Two warps have been seen, one in the outer parts, at radii larger than 18 kpc, and one in the nuclear disk, inside 6 kpc. Chemin et al. (2009) found that the nuclear region shows a “warp in the warp” and the gas is almost face-on here, according to rotation velocities, inclinations, and position angles. Ciardullo et al. (1988) also noted that the ionized gas distribution is more circular in the nuclear region, implying a face-on disk. Over larger scales, an interpretation of the morphological and kinematical perturbations in the central 4 kpc of M31 has been proposed in terms of a bar, along with the accompanying non-circular motions and velocity jumps in shocks (Blaña Díaz et al. 2018; Opitsch et al. 2018; Feng et al. 2022, 2024).

Although the disk of M31 does not show evidence of a thin bar in the old stellar population, the triaxiality of the bulge has been recognized for a long time on red images (Lindblad 1956; Stark 1977; Athanassoula & Beaton 2006; Beaton et al. 2007). Stark (1977) showed that the triaxiality of the bulge is witnessed by the twist of isophotes, which are misaligned with the disk PA, by an angle of ∼10°. By modeling this morphology, Stark (1977) estimated that a star in the disk at 5 kpc from the center would feel a tangential force typically about 5% of the radial force. The triaxial bulge has axis ratios in the range 1:0.6:0.4, and extends to a radius of 2.6 kpc. The existence of an underlying bar in the disk is not obvious, due to the unfavorable inclination of the galaxy, but has been inferred from face-on models to be ∼4 kpc in radius (Athanassoula & Beaton 2006). However, the thin part of the bar, outside the boxy bulge, is not actually observed. The position angle of the near-infrared light outside the bulge has the same PA as the main disk (Beaton et al. 2007). The kinematics of the atomic gas shows some non-circular motions in the center, but nothing specific to the expected S-shape feature in the velocity field, due to a bar in the disk (Braun 1991; Chemin et al. 2009). There is no bar signature in the HI morphology, but H I is depleted within 2 kpc. Instead there is a highly contrasted ring structure at radius 10 kpc, which appears superposed to some pieces of spiral arms. In the PV-diagram, it is possible to distinguish two H I rings, at 2.5 and 5 kpc, which are not centered, but lopsided (Chemin et al. 2009).

The ionized gas is not completely depleted in the center, and although the morphology does not show a bar response, but more a ring with apparent shocks, the kinematics reveal non-circular motions (Opitsch et al. 2018). In the center, the velocity field of the [O III] line shows a kind of S-shape, mostly on the SE side. On the contrary, the stellar velocity field reveals regular rotation, showing no signature of bar and no S-shape feature.

Blaña Díaz et al. (2017) provided a density model of the M31 bulge, composed of one-third nitial classical bulge (ICB) and two-thirds box-peanut bulge (BPB). The BPB alone is not sufficiently concentrated to fit the data. The triaxial bulge rotates with a pattern speed of Ω_p ∼ 40 km/s/kpc, such that corotation is present at 5.8 kpc and OLR (outer Lindblad resonance) at 10.4 kpc. N-body simulations of only collisionless components (stars and dark matter) were performed and compared to the stellar data (morphology and velocities). The best fit is obtained from the morphology, position angle, boxiness, ellipticity, and asymmetry of the simulations within a radius if 4 kpc. Given their best fit model, and using the M2M (made-to-measure) technique to converge more accurately with the observations, Blaña Díaz et al. (2018) also fit the kinematical data obtained with the VIRUS-W IFU by Opitsch et al. (2018). This allowed them to find the best mass-to-light ratio to account for the rotation curve and fit the stellar velocity field together with the velocity dispersion. The final model is a refinement and a better fit of the Gauss-Hermite kinematic parameters h3 and h4. More information could be provided via the comparison with the gas kinematics.

This approach was undertaken by Feng et al. (2022), who analyzed the gas velocity field of the H I from Chemin et al. (2009) and the [O III] from Opitsch et al. (2018). They identified velocity jumps on a series of slits parallel to the minor axis. The width of the slits is 1.2′ = 274 pc, which is also the spacing between slit position. There are thus 17 slits from −9.6 to 9.6 arcmin from the minor axis. The position-velocity diagrams (PVD) are quite noisy in the [O III] line, and when jumps are seen, they do not involve a large flux. Those correspond to the S-shape already remarked by Opitsch et al. (2018) in their velocity field. In the H I gas, the PVD are rather regular and expected for a smoothly rotating gas, which is not surprising, since the H I morphology does not show any signature of a bar. Hydrodynamical simulations are then carried out to illustrate the comparison, but not meant to match all the details of shocks. The result of the gas simulation has a morphology quite classical for a barred spiral galaxy, but very far from the M1 observations: there is a high central gas concentration in a nuclear disk with almost circular orbits, slightly perpendicular to the bar (×2 orbits), and then very narrow nuclear spirals, developing in an empty region, until a radius of ∼3 kpc. This morphology is not seen either in H I nor in the ionized gas.

A more detailed analysis is provided by Feng et al. (2024), with more adapted hydrodynamical simulations, using the same best fit model for the potential, obtained by the above M2M method. To obtain a better fit with the dynamical model, and to reproduce the position and amplitude of the identified shocks, they need a twice lower pattern speed of the bar, this time Ω_p ∼ 20 km/s/kpc, such that corotation is now at 14 kpc. It is only if they allowed the nuclear disk to be tilted, with a more face-on orientation with respect to the large-scale disk, that the pattern speed could be increased to converge with the best fit of the stellar component. However, in this model, there is no reason for such a tilt. Again, the gas morphology in the model does not correspond to M31 observations, with a very high gas concentration in the center, with no lopsidedness. This disagreement of gas density is important, since the shocks are due in the model to gas with highly elliptical streamlines, spending a large time of their orbit in the center. However, this cannot be seen in the observations.

In summary, the bar model has serious problems when we are attempting to account for the gas morphology and kinematics in the central kpc of M31. Although there is a triaxial bulge, with significant pattern rotation, in addition to the classical bulge, which is able to reproduce the stellar kinematics, there might not be a thin bar present in the disk. One possibility is that M31 was a barred spiral galaxy before the encounter with an M32-like companion, as proposed by Block et al. (2006). This almost head-on collision would have destroyed the thin bar in the disk, while propagating a ring wave, explaining the 10 kpc ring feature in the interstellar medium. The dynamically hotter triaxial bulge would have remained unperturbed after the encounter.

7. Conclusion

We have built a fully sampled data cube of the ionized gas within a radius 1.2 kpc in the center of M31, through Hα and [N II] mapped with SITELLE, at CFHT. Combined with the more patchy molecular gas velocity field, previously obtained with the CO lines at IRAM-30 m telescope, and the dust photometry, we propose an interpretation of the complex central dynamics via a static modeling. We have introduced three dynamical components: the main disk with the inclination of 77° and PA = 37°, a tilted ring, required to explain spectra with double-velocity components, and a nuclear warped disk. The mass model of the central kpc is dominated mainly by the bulge and the nuclear disk, although we added the main disk and the dark matter halo to fit the rotation curve. The kinematics of the ionized and molecular gas is then computed in this potential and the velocity field confronted with the observations. The tilted ring component is allowed to be off-centered, with a spatial excursion reproducing the observations. The morphology and orientation angles of the tilted ring and warped nuclear disk were optimized to obtain the best qualitative compatibility with the observations. Several close configurations are possible and the complex dynamics of the central region is thought to be the result of a recent head-on collision with a M-32 like galaxy, as previously proposed by Block et al. (2006). The dynamical perturbations observed correspond to the nuclear disk tilted and warped in the collision which is now settling down back to equilibrium, keeping some lopsidedness, while the m = 1 waves are damping down. The central kpc of M31 is still embedded in a larger triaxial box-peanut-shaped bulge, as a remnant of the pre-collision bar. However, the gas disk within 1–2 kpc is too perturbed to show any bar signature.

Acknowledgments

We are most grateful to the anonymous referee for the very constructive comments that helped us to substantially improve the manuscript. This work is based on observations carried out under project numbers 067-11 and 221-11 with the IRAM-30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This paper is based on observations obtained with SITELLE, a joint project of Université Laval, ABB, Université de Montréal, and the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The authors wish to recognize and acknowledge the very significant cultural role that the summit of Mauna Kea has always had within the indigenous Hawaiian community.

References

All Tables

All Figures

Fig. 1. Dust location in the central kpc field of view of Andromeda. Left panel: Dust emission map obtained from the 8 μm Spitzer map with a subtraction of the stellar continuum at 3.6 μm. (Block et al. 2006). The contours correspond to the ionized gas intensity displayed in the top left panel of Fig. 5. The dashed square corresponds to the inner ring region displayed in the right panel. Right panel: an upper limit estimate on the AB extinction (Melchior et al. 2000). The gray contours correspond to the dust emission at 8 μm displayed in the left panel. The black cross symbol represents the M31 center.

In the text

Fig. 2. Background subtraction. Top panel: Fit of the sky lines over a low brightness region of the field of view. The vertical gray lines indicate the position of the atmospheric sky lines. The five red vertical lines indicate the position of the foreground emission of the DIG of our own Galaxy in the [N II] (6563, 6584), [S II] (6717, 6731) and Hα lines at a velocity of −38 ± 5 km/s. Middle panel: Sky lines model and Galactic foreground model. During the fitting procedure the velocity of the sky lines and the velocity of the Galactic foreground were two independent parameters. Bottom panel: Analyzed spectrum and model of the background.

In the text

Fig. 3. Example of the old stellar population modeling in a 200 × 200 pixels2 box near the center of the galaxy. Top panel: Random spectra taken in the box to illustrate the strong brightness gradient. Bottom panel: Same spectra after a normalization by their mean flux. The resulting median spectrum of all the normalized spectra contained in the box is shown in blue, the smoothed version used as an estimation of the old stellar population is shown in red. The noise level of the individual 7 × 7 binned spectra used during the fitting procedure is shown as an orange envelope.

In the text

	Fig. 4. Examples of emission lines fitting on background corrected spectra displaying only one (top panel) or two components (bottom panel). Top panel: Fit over a background subtracted spectrum revealing the pure emission of the DIG. Bottom panel: Spectrum displaying two resolved velocity components along the line of sight.
In the text

Fig. 5. Intensity and kinematic information extracted from the ionized gas with SITELLE SN3 data cube. Top panels: Ionized gas intensity (left) and the heliocentric velocity (right) corresponding to the main velocity component. The ionized gas intensity has a threshold of 6 × 10−18 erg cm−2 s−1 Å−1. The dashed squares correspond to the field of view displayed in the bottom panels. Bottom panels: Flux ratio ϕ2/ϕ1 (left) and velocity difference V2 − V1 (right) for the second kinematic component with respect to the first one. The contours correspond to the intensity levels (1.8 × 10−17, 3.6 × 10−17 erg cm−2 s−1 Å−1) of the flux map of the main component (top left). The black cross symbol representing the M31 optical center, 00h42m44.37s + 41d16m08.34 (Crane et al. 1992).

In the text

	Fig. 6. CO(2–1) spectra averaged over ∼60″ regions in diameter in the northwest (top) and the southeast (bottom) of the observed IRAM-30 m moment-zero map. The NW corresponds to (−160″,  120″) offsets, and the SE to (150,   − 110″) offsets relative to the M31 center RA = 00:42:44.35 Dec = 41:16:08.6. They correspond to the same regions as Fig. 7.
In the text

	Fig. 7. Hα spectra averaged over 12″ regions in the southeast (left) and the northwest (right) of the observed SITELLE map. The SE corresponds to (124″, −72″) offsets, and the NW to (−199, 72″) offsets relative to the M31 center RA = 00:42:44.35 Dec = 41:16:08.6.
In the text

	Fig. 8. Velocity dispersion derived from the Gaussian fit on SITELLE sky-subtracted cube.
In the text

Fig. 9. Measured broadening when using a sincgauss model (broadened Gaussian emission-line) to fit two unresolved emission-lines. We have simulated the fit with a sincgauss model of two unresolved emission-lines with different flux ratios. We can clearly see that the measured broadening can be very different from the real velocity separation and depends strongly on the flux ratio of the unresolved lines. Note that, at the resolving power of our data (R = 4800), two components are considered resolved when separated by more than 63 km s−1. The dotted line is the bisectrix, meaning equality between the two quantities.

In the text

Fig. 10. Simplified illustration of the map processing, on N = 9 concentric circles, with radii regularly spaced between Rmin = 0 and Rmax = 0.8 kpc. Top panel: Velocity map overlaid by the concentric circles (in gray dotted lines). The red dots (resp. blue triangles) are the maximal (resp. minimal) values Vmax (resp. Vmin) for each circle, and have been computed with n = 7. Bottom panel: Previous extrema Vmax and Vmin of the velocity displayed as a function of the distance from the center (i.e., the radius of its corresponding circle). The negative distances corresponds to the blue shifted part of the map, where the velocities are negative. The gray lines represent the origin axes.

In the text

Fig. 11. Minimal and maximal values, Vmin, $V_{min}^{2\mathrm{nd}}$, and Vmax, of the velocity on N = 300 concentric circles on the map, centered on the BH, with radii regularly spaced from Rmin = 0 to Rmax = 1 kpc. Left panel: Velocity map overlaid by the extreme points of the velocity. The red dots (resp. blue triangles) are the maximal (resp. minimal) values at each radii. As discussed in the text, the green points correspond to the $V_{min}^{2\mathrm{nd}}$ velocity found when considering the second component detected and presented in Fig. 5. Right panel: Minimal and maximal values of the velocity, as a function of the distance from the BH. The negative distances corresponds to the blue shifted part of the map.

In the text

Fig. 12. Zoom-in on the winding in the very center of the map. The analysis is conducted on N = 100 radii from Rmin = 0 to Rmax = 0.35 kpc. Top panel: Vmin (blue), $V_{min}^{2\mathrm{nd}}$, and Vmax velocities superposed on the kinematic map, together with the maximal velocity (red bullets), minimal velocity (blue triangle) and minimal velocity of the second component (green triangle), for each radii. Bottom panel: Corresponding velocity values as a function of the galactocentric radius. The negative radii correspond to the distance of the points in the southwest region.

In the text

	Fig. 13. Position angle (PA) of the main kinematic features as a function of the galactocentric center, computed for the maximal values of the velocity in the northeast region (red) and the southwest region (blue). We distinguish the winding below 300 pc. Beyond 300 pc, the position angle slowly tends toward approximately 40°.
In the text

	Fig. 14. Cross-section of the gaseous disk along a diameter. The crossing of the r and z axis represents the center of the galaxy.
In the text

Fig. 15. Schematic representation of the geometry of the gas in the central kpc of M31. The graph represents a simplified cross section of the gas disk in the sky plane containing the central black hole, represented by the black cross. The blue part is the main disk, aligned with the galactic disk of M31, the green part represents the ring and the red part shows the nuclear disk. The green cross represents the center of the ring, offset with respect to the nucleus. The dashed part of the red line corresponds to the warped disk, whereas the full lines are stiff components. The dashed black line points to where the disk is supposedly teared. We also show a number of parameters of the modeling: rw is the distance between the central black hole and the end of the warped region; rmd indicates where the main disk begins; rr is the inner radius of the ring, as well as the end of the nuclear disk as they are supposed to be from a single structure that teared; wr is the width of the ring; and rid is the radius of the stiff disk inside the inner region. The substructures are not to scale for visibility reasons.

In the text

Fig. 16. Schematic representation of the nuclear warped disk (corresponding to the red substructure in Fig. 15). The black star is the center of M31. The x axis represents the line of sight, and the yOz plane the plane of sky. The filled pink color refers to the inner nuclear disk, and also to the end of the warped region. The dashed lines represent two perpendicular diameters of the substructure. One of the two coincides to the major axis of the projection. The plain lines show where the second diameter would be if the substructure was not inclined (i.e., fully in the plane of sky). The successive circles in the middle are a simple representation of the warped region, defined as a set of misaligned circles with continuously varying inclination and position angles. For visibility reasons, the different parts of the substructure are not to scale with respect to each other, and only a few of the misaligned circles have been drawn. Bright red line at the outer bound of the substructure shows where the warped disk supposedly ripped from the 1 kpc ring. The nuclear warped disk is located inside the 1 kpc radius of M31, superposed with the ∼1 kpc inner ring and the main disk.

In the text

Fig. 17. Schematic representation of the main disk and the inner ring (corresponding to the blue and green substructures in Fig. 15). The stars represent the center of the ring (in green), and the center of M31 (in black). In our modeling, we neglect the component of the offset parallel to the line of sight, and only consider the offset toward the east direction (here $-\overrightarrow{u_y}$) and north direction (here $\overrightarrow{u_z}$). The x axis represents the line of sight, and the yOz plane the plane of sky. The green dotted and dashed axis is parallel to the z axis. In each color, the dashed lines represent two perpendicular diameters of the substructures. One of the two coincides to the major axis of the projection. The plain lines show where the second diameter would be if the substructure was not inclined (i.e., fully in the plane of sky). The angular parameters of the modeling are (ir, θr) and (imd, θmd); the inclination and position angles of the ring and of the main disk.

In the text

	Fig. 18. Central 75 × 75 arcsec2 (285 × 285 pc) field of view of Andromeda. The white cross corresponds to the optical center. The intensity contours are the same as in the bottom panels of Fig. 5.
In the text

Fig. 19. Theoretical rotation curve derived from Eq. (4), using parameters from Table 1. The contributions are from the nucleus (green), the bulge (blue), the stellar disk (orange), the gaseous disk (red), and the dark matter halo (purple). The black line is the total rotation curve. The blue triangles correspond to the measured ionized gas velocities, averaged over the NE and SW side, and corrected from asymmetric drift.

In the text

Fig. 20. Four scenarios proposed for the central region in M31, represented by their schematic transversal cut. Grey lines correspond to transition zones, blue lines to the main disk, red lines to the warped nuclear disk, and green lines to the inner ring. When two structures are decoupled, two distinct velocity components may be observed along the same line of sight, while when they are linked by a transition zone, there would be a single wide velocity component. Top-left: Scenario 1, where the nuclear and main disks are linked. The inner ring is decoupled, implying two velocity components along a line of sight. Top-right: Scenario 2, where the nuclear disk is linked to the inner ring. The main disk is decoupled, with possibly two velocity components on some line of sights. Bottom-left: Scenario 3, where the three structures are decoupled, implying distinct velocity components along a line of sight. Bottom-right: Scenario 4, which is the same as scenario 3, but there is now a link between the main disk and inner ring, limiting the occurrence of distinct velocity components.

In the text

	Fig. 21. Velocity field obtained in the four scenarios described in Fig. 20, in the same order. The maps have been masked with the velocity map of Fig. 5; namely, we do not plot any velocity when there is no signal in the data.
In the text

	Fig. 22. Same as Fig. 21 for the second moment of the cube. The maps have been masked with the velocity dispersion map of Fig. 8; that is, we do not plot any velocity dispersion when there is no signal in the data.
In the text