Evidence for non-axisymmetry in M31 from wide-field kinematics of stars and gas

As the nearest large spiral galaxy, M31 provides a unique opportunity to learn about the structure and evolutionary history of this galaxy type in great detail. Among the many observing programs aimed at M31 are microlensing studies, which require good three-dimensional models of the stellar mass distribution. Possible non-axisymmetric structures like a bar need to be taken into account. Due to M31's high inclination, the bar is difficult to detect in photometry alone. Therefore, detailed kinematic measurements are needed to constrain the possible existence and position of a bar in M31. We obtained $\approx$ 220 separate fields with the optical IFU spectrograph VIRUS-W, covering the whole bulge region of M31 and parts of the disk. We derive stellar line-of-sight velocity distributions from the stellar absorption lines, as well as velocity distributions and line fluxes of the emission lines H$\beta$, [OIII] and [NI]. Our data supersede any previous study in terms of spacial coverage and spectral resolution. We find several features that are indicative of a bar in the kinematics of the stars, we see intermediate plateaus in the velocity and the velocity dispersion, and correlation between the higher moment $h3$ and the velocity. The gas kinematics is highly irregular, but is consistent with non-triaxial streaming motions caused by a bar. The morphology of the gas shows a spiral pattern, with seemingly lower inclination than the stellar disk. We also look at the ionization mechanisms of the gas, which happens mostly through shocks and not through starbursts.


Introduction
Because of its proximity, M 31 is the best case after the Milky Way to learn about the detailed evolutionary history of a large spiral galaxy. Therefore, M 31 has been studied by several major surveys in recent years, combining large-scale photometry with pointed spectroscopic observations. Two of these large programs are the Spectroscopic and Photometric Landscape of Andromeda's stellar halo survey (SPLASH; Guhathakurta et al. 2005Guhathakurta et al. , 2006Gilbert et al. 2006) and the Pan-Andromeda archaeological survey (PAndAS; McConnachie et al. 2009). These surveys studied the stellar halo of M 31 in great detail, to measure its global properties (Gilbert et al. 2012(Gilbert et al. , 2014 and look at structures within the halo like the Giant Southern Stream (Gilbert et al. 2009). This provided information about the formation history of M 31. One of the results of the PAndAS survey is that white dwarfs galaxies around M 31 all lie in a thin plane (Ibata et al. 2013), which poses problems for the current understanding of galaxy formation. The Panchromatic Hubble Andromeda Treasury (PHAT) program (Dalcanton et al. 2012) This paper includes data taken at The McDonald Observatory of The University of Texas at Austin. This research was supported by the DFG cluster of excellence "Origin and Structure of the Universe".
Full  are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/611/A38 looked at the disk of M 31 and obtained photometry for 117 million individual stars (Williams et al. 2014). The Herschel Exploitation of Local Galaxy Andromeda (HELGA; Fritz et al. 2012) observed M 31 in the far infrared and sub-millimeter wavelengths, measuring the distributions of dust ) and molecular clouds (Kirk et al. 2015). Microlensing events and variable stars were monitored by the Pan-STARRS 1 survey of Andromeda (PAndromeda; Lee et al. 2012Lee et al. , 2014a. All these surveys help to understand the structure, the accretion history and the history of star formation within M 31. M 31 is classified as an unbarred spiral SA(s)b by de Vaucouleurs et al. (1991), however, there has been evidence from photometry and kinematics that it is barred, see Sect. 1. In this paper we present new spectroscopic observations of M 31. The paper focusses on the description of the data and the hints for a bar that can be seen there in. A discussion of other possible explanations, like a superposition of disks and rings, goes beyond the scope of this paper and will be presented in future papers based on the dataset described here.
The current Lambda cold dark matter (ΛCDM) paradigm (Planck Collaboration XIII 2016) has been very successful in explaining both the large scale structure of the universe and the observed properties of galaxies. One component of this model is the cold dark matter, which resides in halos around galaxies. The nature of this dark matter is not yet understood. The currently preferred candidates for dark matter are non-baryonic elementary particles, so-called weakly interacting massive particles (WIMPs), see for example Bertone (2010) and references therein. Baryonic candidates for dark matter are mostly ruled out, because the fraction of baryonic matter in the universe is only 15% of the total matter (Planck Collaboration XVI 2014). This baryonic dark matter fraction could be composed of large astrophysical objects, like brown dwarfs, Jupiter-sized planets or black holes, collectively known as massive astrophysical compact halo objects (MACHOs; Griest 1991). In order to place firm constraints on the existence of these massive objects, microlensing studies have been conducted. In such a microlensing event, a MACHO passes between a bright object and the earth, the light from the source objects is deflected by the gravity of the MA-CHO, which leads to a perceived increase in brightness of the source object. The observation of the frequency of such events allows to determine the number density of MACHOs (Paczynski 1986). The MACHO survey (Alcock et al. 1993) finds that the contribution of MACHOs to the total halo mass of the Milky Way is 20% (Alcock et al. 2000), with an average MACHO mass of ≈0.4 M . The Expérience pour la Recherche d'Objets Sombres (EROS) and EROS-2 projects (Aubourg et al. 1993;Afonso et al. 2003) measure a significantly lower fraction for the same masses, at less than 8% (Tisserand et al. 2007). The optical gravitational lensing experiment (OGLE; Udalski et al. 1992) finds a fraction that is comparable to the one measured by the EROS-survey (Wyrzykowski et al. 2011), less than 7% for MACHO masses lower than 1 M . In addition to these results on the galactic halo, the galactic bulge microlensing surveys within OGLE-III (Wyrzykowski et al. 2015), microlensing observations in astrophysics (MOA) II (Sumi et al. 2013) and EROS-2 (Hamadache et al. 2006) have become major tools for understanding the structure of the Milky Way, see for example Wegg et al. (2016) for a new analysis.
Several microlensing surveys have also been focused toward M 31, like the PAndromeda survey (Lee et al. 2012) and the Wendelstein Calar Alto pixellensing project (WeCAPP; Lee et al. 2015). A total of 56 events have been detected in M 31 (Lee et al. 2015). These events do not have to be caused by MA-CHOs, they can also happen due to lensing by other stars in M 31, so-called self-lensing (Riffeser et al. 2006). To constrain the number of events that are caused by self-lensing, a proper understanding of the three-dimensional distribution of the stars is needed. Models with different mass distributions of the galaxy result in widely varying predictions of the event rate of selflensing and lensing through MACHOs. It is therefore important to model the stellar mass distribution in the galaxy as accurately as possible. Based on the data presented in this paper and numerical models from Blaña Díaz et al. (2017, hereafter B17) and Blaña et al. (in prep.), new predictions for microlensing events will be presented in a future paper (Riffeser, in prep.).
A large fraction of disk galaxies in the local universe is barred, ranging from about 50% in the optical (Barazza et al. 2008) to about 60% to 70% in the infrared (Eskridge et al. 2000;Menéndez-Delmestre et al. 2007). It is now thought that global instabilities in the disk lead to the quick formation of bars. In this process, the m = 2 mode grows strongly by swing-amplification and forms a long-lasting bar non-linearly (Sellwood & Wilkinson 1993;Sellwood 2013). Over time, the inner part of the bar goes through a buckling phase, which is a short but violent vertical instability not long after bar formation (Combes & Sanders 1981;Combes et al. 1990; Raha et al. 1991;Merritt & Sellwood 1994;Athanassoula & Misiriotis 2002;Debattista et al. 2005). The instability bends out of the plane of the disk, then settles back to the plane, redistributing energy to smaller spatial scales and to higher stellar velocity dispersion, thereby thickening the bar (Raha et al. 1991). The buckled part of the bar becomes a threedimensional so-called boxy/peanut shaped (B/P) bulge, the part that has not buckled is referred to as the thin or flat bar. While this buckling phase is frequently seen in simulations, it has only recently been detected in observations by Erwin & Debattista (2016) for two local spiral galaxies.
While the Milky Way was originally thought of as unbarred, it is now widely accepted that it contains a bar. For a review, see Gerhard & Wegg (2015), and references therein. Recently, signs for a bar have also been detected in the innermost parts of the third large spiral galaxy in the Local Group, M 33 (Hernández-López et al. 2009).
The question if M 31 also contains a bar or not has not been settled, for example Widrow et al. (2003) refer to the bulge as "barlike", while Kormendy et al. (2010) classify M 31 as containing a classical bulge. If present, a bar is not easily detected in images of M 31 because of its high inclination of 77 • (Walterbos & Kennicutt 1987). This is too high to see a bar directly in the image, but too low to recognize its shape above and below the stellar disk, as is possible in an edge-on view (Athanassoula & Beaton 2006). Nevertheless, the boxy appearance of the isophotes in near-infrared images already hint at the existence of a bar (Beaton et al. 2007). However, boxy isophotes do not need to be caused by bars, there are numerous examples of early-type galaxies that are boxy without having a bar, see for example Kormendy et al. (2009).
It is possible to detect a bar in a galaxy with an inclination similar to the one of M 31, Kuzio de Naray et al. (2009) investigate the galaxy NGC 2683 with a similar inclination to M 31, by studying ionized gas velocities and the overall morphology. By combining the photometric and kinematic data, they find evidence for the presence of a bar in NGC 2683 and constrain its orientation and strength.
According to Stark & Binney (1994), there are three arguments for a triaxial structure in M 31: 1. There is a twist in the inner isophotes in the bulge with respect to the outer disk, first seen by Lindblad (1956). He was subsequently the first one to claim that M 31 has a bar. These twists cannot be reproduced by a rotationally symmetric distribution of stars (Stark 1977). 2. The velocities of the H i gas are not symmetric about the minor axis (Rubin & Ford 1971). 3. The ionized gas has the appearance of a spiral pattern, which is rounder than the appearance of the disk, as seen by Jacoby et al. (1985), Boulesteix et al. (1987) and Ciardullo et al. (1988). Stark (1977) shows that the features measured by Lindblad (1956) can be explained by a family of triaxial bulge models. Stark & Binney (1994) narrow these models down by simulating the velocities of the gas in this potential. Berman (2001) and Berman & Loinard (2002) simulate the gas velocities in the triaxial bulge potential that was derived using the method of Stark (1977) and they are in agreement with the non-circular gas velocities in the inner disk. Their model has a fast pattern speed of 53.7 km s −1 kpc −1 . Therefore, it would be more fitting to call their triaxial bulge a bar. In our understanding, a triaxial bulge would be a non-rotating structure. Most observed bars have after the buckling phase a threedimensional inner part of the bar, the B/P bulge, and a flat outer part (see e.g., Athanassoula 2005or Martinez-Valpuesta et al. 2006. In this nomenclature, which we adopt throughout this paper, the "triaxial bulge" of the models by Berman (2001) and Berman & Loinard (2002) is a B/P bulge. According to Gordon et al. (2006), the model by Berman (2001) explains the morphology of dust in M 31, with spiral arms emerging from the bar. However, the fact that the two prominent dust rings do not share the same center, which also do not coincide with the optical center of M 31, lead Block et al. (2006) to propose a different scenario, where these rings are not created by a bar, but instead are shock waves due to the collision of the small companion galaxy M 32 with M 31. Athanassoula & Beaton (2006) test four different bar models and qualitatively compare the velocities to H i kinematics from Rubin & Ford (1970), Brinks & Shane (1984) and Brinks & Burton (1984), and the overall morphology to observations in the near infrared by Beaton et al. (2007). They find that in order to explain the boxy appearance of the isophotes in Beaton et al. (2007), a classical bulge needs to be present. The triaxial bulge seen by Lindblad (1956), Stark (1977) and Stark & Binney (1994), corresponds to the B/P bulge from Athanassoula & Beaton (2006). The fact that the boxy isophotes in Beaton et al. (2007) do not coincide with the disk argues for a misalignment of the bar and disk.
While the arguments for a bar in Athanassoula & Beaton (2006) are mostly qualitative, a more quantitative result is obtained by B17, who test 84 different models and compare them to 3.6 µm infrared photometry from Barmby et al. (2006), H i kinematics from Chemin et al. (2009) andCorbelli et al. (2010), as well as stellar kinematics from Saglia et al. (2010, hereafter S10) and data from the work presented in this paper. Again, they rule out solutions which do not contain a classical bulge component embedded within the B/P bulge, finding for the mass of the classical bulge M class = 1.0−1.4 × 10 10 M and for the half mass radius r class = 0.5−1.1 kpc. In the preferred model in B17, the mass of the classical bulge component is M class,best = 1.1×10 10 M and the half-mass radius is r class,best = 0.53 kpc, for the B/P component, the parameters are M B/P,best = 2.2 × 10 10 M and r B/P,best = 1.3 kpc. The projected position angle of the bar is PA bar = 55.7 • , which is 17.7 • more than the disk position angle of PA = 38 • (de Vaucouleurs 1958). The bar in the model by B17 entered the buckling phase about 2 Gyr after it formed, the buckling phase ended about 1 Gyr later. This happened at least 0.8 Gyr ago, but it could have happened earlier, because once the buckling phase is over, the galaxy evolves very slowly with time, so it can't be said when these events happened exactly. The intrinsic length of the bar is 1000 , which projected onto the sky with M 31's orientation and inclination becomes 600 .
According to the models by Athanassoula & Beaton (2006) and B17, there are at least the following separate structural stellar components in M 31, from the innermost to the outermost: 1. A classical spherical bulge in the center, 2. a B/P bulge, which is the inner thicker part of the bar, 3. a thin bar, this is the outer part of the bar, 4. a disk and 5. a halo. This is the first paper in a series of papers about our observations of M 31, this one covers the description of the data, the kinematics, and the gas fluxes, in an accompanying paper we will present results on the stellar populations. The data will also be made available at the CDS. This paper is structured as follows: In Sect. 2, our observations of M 31 are described, before we present the methods used to fit the kinematics in Sect. 3. Section 4 then presents the results for the stellar and gas kinematics, as well as the gas morphology.
In Sect. 5, we search for arguments for the bar in the data, before summarizing our findings in Sect. 6.

The IFU spectrograph VIRUS-W
The observations were carried out with the integral-field unit (IFU) spectrograph VIRUS-W (Fabricius et al. 2012a) at the 2.7 m telescope at the McDonald observatory. The IFU consists of 267 fibers which are arranged in a rectangular hexagonal dense-pack scheme with a filling factor of 1/3. The field-of-view of the instrument is 105 × 55 at the 2.7 m telescope, with the long edge of the fiberhead aligned along the east-west axis. Each fiber covers a circle with diameter 3.2 on sky. The actual spectrograph has two different resolution modes, each realized with a volume phase holographic (VPH) grating. We used the high-resolution mode, where the grating has a line frequency of 3300 lines per millimeter and a resolution of R ≈ 9000, which corresponds to an instrumental dispersion of σ inst = 15 km s −1 . For a VPH a change of the grating angle results in an effective change of the blaze function, such that the throughput for a specific wavelength range is optimized. For our observations, we adjusted the grating angle to 353 • which puts the highest throughput between 4900 Å and 5100 Å and presented a good compromise between covering the Hβ absorption feature at 4861 Å and the Mgb region around 5180 Å. The complete wavelength range is 4802 Å to 5470 Å. In this range, we see the emission lines Hβ at 4861 Å, the doublet of forbidden lines of doubly ionized oxygen [O iii]λλ4959, 5007 and the doublet of the forbidden nitrogen lines [N i]λλ5198, 5200.  (Kormendy & Bender 1999). Therefore, we call all pointings in this area "bulge pointings" and the ones in the stripes "disk pointings". Along the major axis, we reached approximately one disk scalelength of r d = 24 or 5.3 kpc (Courteau et al. 2011). We did not dither our observations, because we were primarily interested in covering a large area of M 31 efficiently. We observed each galaxy pointing with an exposure time of ten minutes, except for the central pointing, which was only observed for five minutes, because it covers the bright nucleus of M 31, where sufficient signal-to-noise values were already reached with this shorter exposure time. Before and after each galaxy pointing, we nodded the telescope away from the galaxy to a sky position, which was exposed for five minutes. The seeing varied between 1.3 and 3.0 during the observations. At the beginning and at the end of each observation night, we took bias, flat and arc images for calibration.

Data reduction
The data reduction follows the standard procedure for VIRUS-W as described in Fabricius et al. (2014). We used the fitstools package (Gössl & Riffeser 2002) and the Cure pipeline developed for HETDEX (Hill et al. 2004). First, we created master biases, flats and arcs by taking the mean of the individual images for each morning and evening with appropriate clipping of spurious events. The master bias frames were then subtracted from all other frames. Cure traces the fiber positions on the master flat frames and then extracts the positions of the spectral line peaks along these traces in the master arc frames. To model the distortion and the spectral dispersion, a two-dimensional seventh degree Chebyshev polynomial is used. The resulting model transforms between pixel positions on the detector and fiberwavelength pairs and vice-versa. Foreground stars, which appear as brightly illuminated fibers, were removed using a κ-σclipping method.
We used 27 lines for the wavelength calibration. They are listed in Table B.1. After the tracing of the fiber positions and the calibration of the wavelengths, we extracted spectra from the science frames by moving along the trace center positions and averaging the values in a 7 pixel wide aperture. The extraction was performed in ln(λ)-space, the step width corresponds to 10 km s −1 . The flatfield frames were extracted in the same way as the data frames. We corrected the fiber to fiber throughput variation and the vignetting by dividing each observation by the corresponding flatfield observation. The flatfield data are very stable over the course of one observing night, with the deviations being on the order of 0.03%. However, the resulting flat-fielded spectra still exhibit the rather strong variation of sensitivity as function of wavelength that is due to the strongly peaked diffraction 1 The image has been taken from http://archive.stsci.edu/ cgi-bin/dss_form. The digitized sky surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. c 1995 by the Association of Universities for Research in Astronomy, Inc. efficiency of the VPH grating. This would have complicated the later throughput calibration. Therefore, in the next step we divided the spectra at each wavelength by the mean value of the flat field spectrum at that wavelength, where the mean is taken across all fibers in the flatfield observation.
To remove the contribution to the signal by the sky background, we pointed the telescope to a sky position before and after each galaxy pointing. The sky frames were fiber-extracted and flat-fielded in the same way as the galaxy observations. They are quite homogeneous after being flat-fielded, with a root-meansquare of the order of 1%. We averaged the two sky pointings observed just before and after each galaxy pointing into one. The difference between the two sky pointings is very small, about 3%. We then subtracted this flat-fielded sky image from the flat-fielded galaxy image.
Because the different observing runs took place in different months of the year, we also had to correct during the extraction for the relative motion of the earth around the sun. We used the web-tool by Edward Murphy 2 based on an algorithm described in Meeks (1976) to calculate the relative velocity of the Earth toward M 31 at the time of the observation. For each observing run, we used the value for the mean date of the run. The correction for the run in October 2012 is c Oct12 = −3.6 km s −1 , for the run in February 2013 it is c Feb13 = 19.1 km s −1 and for the run in August 2013 it is c Aug13 = −28.1 km s −1 . For the run in October 2011, the correction is c Oct11 = 0 km s −1 .
The sky position for every fiber in each pointing of M 31 was determined following Adams et al. (2011). The accuracy of this method is estimated by Adams et al. (2011) to be 0.2 , much less than the fiber diameter of VIRUS-W of 3.2 . The coordinates were converted to distance in arcseconds relative to RA = 10:41:07.04 and Dec = 41:16:09.41 (J2000). This is the coordinate of the pointing covering the center of M 31 in our observations.
In the central region, the signal-to-noise ratio (S/N) of each fiber spectrum lies well above 30. This remains true out to a radius of about 140 along the major axis and 100 along the minor axis. At larger radii m fiber spectra were then spatially binned to reach a minimum S/N of 30 using the Voronoi-binning method by Cappellari & Copin (2003), resulting in a total of 7563 binned spectra.

Measuring the stellar kinematics with pPXF
The kinematics were measured using the penalized PiXel Fitting (pPXF) routine by Cappellari & Emsellem (2004) and Gas AND Absorption Line Fitting (GANDALF) code by Sarzi et al. (2006). GANDALF uses pPXF as its first step. pPXF measures the stellar kinematics by broadening a weighted sum of template star spectra with a trial line-of-sight velocity distribution (LOSVD) and subsequently changing the parameters of the LOSVD until the residuals between the measured and the model spectrum are minimized. We used spectra from 41 kinematic standard stars obtained with VIRUS-W. They are listed in Table B.2.
In pPXF, the LOSVD is expanded as a Gauss-Hermite series following van der Marel & Franx (1993) and Gerhard (1993): H m are the Hermite polynomials and hm the Gauss-Hermite coefficients, the sum is broken off after M entries. We only looked at the Gauss-Hermite moments h3 and h4.

Fit of the emission lines with GANDALF
The kinematics of the ionized emission lines were fitted with GANDALF (Sarzi et al. 2006). It fits the kinematics of the [O iii]λ5007 line and ties the kinematics of the other emission lines to that line. In Fig. 2, a fit with GANDALF to a spectrum is plotted. Close examination of the spectra showed that in many regions multiple gas components at different radial velocities exist. In order to properly treat these multiple peaks, we added another [O iii] component to be fitted with GANDALF, as well as a second Hβ component and a second [N i] component. The initial guesses for the gas velocities have to be slightly different for the two components, otherwise GANDALF does not fit separate components. An example of a fit with two lines is shown in Fig. 3. Here, the two lines have almost the same amplitude and are clearly separated.
This was not always seen so clearly, there were also cases where one line was stronger than the other or where the two lines were almost blended together or where there was only one line clearly visible with a skewed line shape, that could, however, be described by the combination of two lines. In Fig. 4, we show where we found one line and two lines. We saw only one line in about 3000 bins, while we saw two lines in about 3500.   Kormendy & Bender (1999) with an ellipse fit to a K-band image by S10.
In order to reliably fit the different components, we had to feed GANDALF information whether only one component was present or if there were really two. In order to get these initial guesses, we applied the following method: First, we cross-correlated a model spectrum only consisting of the [O iii]λ5007 line with each spectrum. The program fitted the resulting cross-correlation function with a set of gaussians. These gaussians all had the same dispersion of σ = 20 km s −1 and their peak velocities were 40 km s −1 apart. The program changed the amplitudes of the individual Gaussians to get the best approximation of the input cross-correlated spectrum. We told the program to only pick the Gaussian with the largest amplitude to have an estimate for the one-component fit and the one with the largest amplitude plus the one with the second largest amplitude for the two component estimate. If the amplitude of the second component was less than 0.25% of the amplitude of the first one, we decided to take the initial guess with only one component. We also discarded unrealistically high velocity dispersions of σ gas > 100 km s −1 . We used the central velocities of the Gaussians as the initial guesses for GANDALF, letting it fit one line for the cases where we had found only one line and letting it fit two lines where we had found two lines. After a first iteration, we checked all fits manually, updated the initial guesses for the spectra where the fit failed and let GANDALF fit a second iteration. This second iteration resulted in 85% of the spectra being fitted correctly, the rest was left out of the analysis.
In order to check if the spectra could be contaminated by a contribution of Planetary Nebulae (PNe), we looked at the catalog of PNe positions from Merrett et al. (2006). If the position of a PNe is closer than 1.6 (the radius of a fiber), we took this bin as affected. Overall, 166 bins were affected, which is 2.5%. All of these spectra showed [O iii] lines, about two thirds of them showed two lines, while one third showed only one line. The spectra did not look systematically different from the unaffected spectra, there was for example no exceptionally bright line in any of them. We therefore concluded that contamination by PNe is negligible.

Flux calibration
In each observation night, we observed photometric standard stars, they are listed in Table B.3. These spectra were compared to photometrically calibrated spectra from the literature (Oke 1990;Le Borgne et al. 2003) to get the throughput for the particular observation night as a function of wavelength. The throughput curves were calculated with a program from Müller (2014), the error for the throughput curve is about 4%.
The atmospheric extinction was also corrected for each individual observation. Differences in observing conditions between the individual pointings were taken into account by comparing the integrated flux f tot in one spectrum to the flux f V of a photometric model image. This image was constructed by combining the bulge-disk decomposition by Kormendy & Bender (1999) with an ellipse fit to a K-band image performed by S10. The model image is shown in Fig. 5.
GANDALF calculates 1-σ errors in the routine, which are consistent with what we obtained from Monte-Carlo simulations of fitting a noiseless spectrum with random added noise. These 1-σ errors are the ones tabulated in Table B.5.

Stellar kinematics
In this section, the stellar kinematics measured with GANDALF is presented. The data will be made available at the CDS, an example is given in Tables B.4 and B.5. The former table lists the individual on-sky fiber positions and the corresponding bin numbers that the fiber signals were assigned to. The latter table gives the actual stellar and gas kinematic data and the gas flux values with respect to the on-sky aperture size of the VIRUS-W fibers.
We compare our data with the measurements by S10 in the in Appendix A. In general, both datasets agree very well.
The stellar velocity map is plotted in Fig. 6, the heliocentric velocity of M 31 (v sys = −300 km s −1 ; de Vaucouleurs et al. 1991) has been subtracted.
A schematic view of M 31 based on Henderson (1979) is plotted in Fig. 7, with the naming of the receding and approaching sides of M 31 taken from the stellar velocity map in Fig. 6. Overall, the stellar velocity field is regular, rotation is clearly visible. The velocities increase strongest along the disk major axis, A38, page 6 of 22 The magenta contours are the surface brightness of the IRAC 3.6 µm image by Barmby et al. (2006). The line in the colorbar is the median errorbar of the velocities. 78 is plotted, this corresponds to an inclination angle of i = 77 • , the disk major axis has position angle PA = 38 • . The half of the galaxy to the north-west of the major axis is the near side (Henderson 1979). with the highest velocity in the bulge region being v bulge,max = 136 ± 4 km s −1 in the outermost bulge pointing in the receding side and the lowest being v bulge,min = −157 ± 4 km s −1 on the opposite side. The "bulge region" is the region of M 31 where the bulge-to-total ratio of M 31 is higher than 0.5, taken from a decomposition by Kormendy & Bender (1999). The absolute maximum of the velocities is reached with v max = 208 ± 3 km s −1 in the outermost disk major axis disk pointing in the receding side, and the lowest value of v min = −193 ± 2 km s −1 already reached in the disk pointing at a radius of 930 on the approaching side. On that side, the velocity roughly remains at that value to the outermost radii. The median velocity error is dv = 3.8 km s −1 . When plotting the velocity map with contours, a twist for the line of zero velocity becomes apparent which not aligned with the photometric minor axis (see Fig. 8).
The stellar velocity dispersion of M 31 is plotted in Fig. 9. The velocity dispersion drops faster along the disk minor axis than the major axis. The red regions are two lobes with high values of velocity dispersion of over 165 km s −1 , they are separated by a line of lower velocity dispersion along the major axis of the bar. Around this we see a region with slightly lower values (≈160 km s −1 ) that has the shape of a small yellow parallelogram. Further out, we see a light blue region, which has even lower values of σ (≈130 km s −1 ), which is elongated along the disk major axis. The maximum value is σ max = 188 ± 5 km s −1 , located at a distance of 46 from the center, the minimum value is σ min = 55 ± 4 km s −1 , in the outermost northern disk pointing at x = −100 and y = 1100 . The mean velocity dispersion for the whole dataset is σ mean = 116 ± 4 km s −1 , in the central 20 it is σ central = 160 ± 5 km s −1 , considering only the bulge region it is σ mean,bulge = 137±4 km s −1 , in the disk it is considerably lower with σ mean,disk = 103 ± 4 km s −1 . The disk velocity dispersion is still quite high, which is in agreement with what Fabricius et al. (2012b) find for similar disk galaxies.
In Fig. 10, we plot the Gauss-Hermite moment h3. In the very central 100 , h3 has high positive and negative values, while further out, the absolute values become lower. Even further out in the bulge, the absolute values become higher again. In the disk along the major axis, the signs of h3 flip and even higher absolute values are reached. Corresponding regions of high and low values are relatively symmetric with respect to the center.
In the central 50 , h3 has values of about ±0.05. Further out, the values drop to about ±0.02, before rising at the edges of the bulge region to ±0.08. In the bulge, the mean absolute value is |h3 bulge | = 0.03 ± 0.02. In the disk, the values are higher, the mean absolute disk value is |h3 disk | = 0.04 ± 0.03. Along the disk major axis, the values are highest, here, the mean absolute value is |h3 major | = 0.07 ± 0.02. The maximum and minimum values of the whole map are also found along the disk major axis, being  h3 max = 0.14 ± 0.02 and h3 min = −0.16 ± 0.02. Over most of the bulge region, we see a correlation between h3 and v, this will be discussed in Sect. 5.
In Fig. 11, the Gauss-Hermite moment h4 is presented. This map shows that along the disk major axis, the absolute values of h4 are very low. At the ends of the bulge region, the values of h4 become negative. Perpendicular to the disk major axis, the values of h4 become positive. The mean value over the whole dataset is h4 = 0.02 ± 0.02. In a 240 wide and 350 long strip along the disk major axis, the mean value is lower, at h4 outer bulge = 0.002. Along the disk minor axis, the values of h4 are generally higher, in regions that have a distance to the disk major axis of more than 120 , the mean value is h4 minor = 0.03. In the disk pointings along the disk major axis, h4 is higher in the northern half of the galaxy, having a mean value of h4 major,north = 0.04, while in the southern half, it is h4 major,south = −0.01.

Gas kinematics
In this section, the kinematics of the ionized gas is presented. These data will be made available at the CDS, an example is given in Table B.5. The motions of the gas are more complicated than the ones of the stars from Sect. 4. In the center two lobe-like regions with high values of velocity dispersion of over 165 km s −1 can be seen. They are separated by a narrow ridge of lower velocity dispersion that is roughly aligned with the major axis. At larger radii one sees a region with slightly lower values (≈160 km s −1 ) that has the shape of a parallelogram. Further out towards the edge of the bulge region, the velocity dispersion takes values of about 130 km s −1 . This region forms an elongated rim around the bulge region and is aligned with the disk major axis. More complicated gas kinematics are often observed in disk galaxies (e.g., Falcón-Barroso et al. 2006;Ganda et al. 2006), because contrary to the dissipationless stars, the gas can also lose energy through radiation. The dense gas traced by ground-state CO and H i transitions is most likely to have settled onto closed orbits via hydrodynamic interactions. Associated with this dense gas are regions of ionized gas (Stark & Binney 1994), which is then seen in the optical emission lines. As mentioned in Sect. 2, we observed in about half of the investigated binned spectra two gas components separated in velocity. We tried several ways to sort the two components, either by flux or proximity in position space. In the end, we settled on sorting by velocity, because it produced the smoothest maps. The velocity map of the first component, which is the faster one of the two, is plotted in Fig. 12, the one of the second component in Fig. 13. The first gas component has a median absolute value of |v [OIII],1 | = 162 ± 5 km s −1 . The maximum value is v [OIII],1,max = 294.7 ± 4.5 km s −1 at the coordinates (−220 , 281 ), at 360 from the center along the disk major axis in the receding side. The minimum value is v [OIII],1,min = −340 ± 3.0 km s −1 at the coordinates (25 , −100 ), which is at about 100 south of the center in the approaching side. The minimum is closer to the center than the maximum on the other side, leading to an asymmetric appearance with respect to the center. Apart from that, there is a "spiral-like signature" in the innermost 100 × 100 and a large "S-shape" between the approaching and receding gas velocities. These structures are marked in the magnified inset of Fig. 12. For the second component, the median absolute value is significantly lower than for the first component, with |v [OIII],2 | = 73.2 ± 5.5 km s −1 . The maximum value is v [OIII],2,max = 183.8 ± 5.2 km s −1 at the coordinates (−190 , 280 ), which is in the same region as the maximum for the first component. The minimum value is v [OIII],2,min = −240.2 ± 3.8 km s −1 , at the coordinates (65 , −170 ), at 182 from the center along the disk major axis in the receding side. This again corresponds to a region of low velocities in the first component. The overall shape of the velocity field for the second component is similar to the first, also with an S-shape along the line of zero velocity. This structure has a slightly different shape than the one for the first component, an arm of approaching velocities is extending further into the region of receding velocities. Additionally, on the western side of the bulge, at about 500 along the disk major axis from the center, there is an arm of velocities of about zero. This arm seems to be unconnected to the kinematics of the rest of this component. These structures are marked in the magnified inset of Fig. 13.
The velocity dispersion of the first component is plotted in Fig. 14. The general appearance is noisy, overall, the mean value is σ [OIII],1 = 35 ± 5 km s −1 . The velocity dispersion is low near the northern end of the bulge along the major axis, where the velocities are high (σ ≈ 50−60 km s −1 ).
In the velocity dispersion values of the second component, a ring-like structure of high values is visible, which increases the mean value of the velocity dispersion to σ [OIII],2 = 37 ± 6 km s −1 . The values in the ring are σ [OIII],2 = 57 ± 10 km s −1 . This ring of high velocity dispersions corresponds to a similar ring in Fig. 4, which is where two peaks are measured. In this ring, a second peak with a lower amplitude than the first peak is present. However, this second peak is broader, therefore the velocity dispersion is larger.

Gas fluxes
The fluxes of Hβ, [O iii]λ5007 and [N i]λ5198 are plotted in Figs. 16 to 24. The values will be made available at the CDS, Fig. 13. Velocity map of the more slowly rotating second gas component. The lines are analogous to Fig. 9. The inset shows a zoom into the inner regions, with the S-shape and the low-velocity arm highlighted, which are mentioned in the text. There are regions where the first component has higher flux than the second one, as well as regions where the opposite is true. When averaging the fluxes along ellipses, however, it becomes apparent that the first component is stronger than the second one. This is true for all different lines. In the Hβ and [O iii] maps, a spiral structure is visible over most of the bulge region, with an incomplete elliptical ring inside, which is oriented roughly along the minor axis, with a "spoke" along the ring's short axis. The ring and spoke are highlighted in the zoomed-in inset of Fig. 18. This ring is at a smaller radius than the one mentioned above in the velocity dispersion of the second gas component. In the southwest, one of the arms of the spiral pattern is prominent in the Hβ and the [O iii]. The fluxes of Hβ are lower than the     , no clear pattern can be seen there apart from the fact that it becomes brighter in the center. The overall filamentary appearance of the gas morphology could be either due to heating by shocks or supernovae of type Ia (Jacoby et al. 1985).
We fitted ellipses to the maps of the [O iii] gas and compared them to the stellar surface brightness. The comparison is shown in Fig. 25. From about 30 to about 200 , the decline of the [O iii] is similar to the one of the stars, while further out, it drops significantly compared to the light of the stars. This corresponds to the edge of the spiral structure visible in the maps. The second [O iii] component has lower flux values, its profile is more irregular, but it also declines fast from 200 outwards.   different strengths and viewing angles. Stark & Binney (1994), Athanassoula & Beaton (2006) and B17 all find that the bar in M 31 is neither viewed end-on nor side-on, but instead at an intermediate angle. The special signatures of a bar are plateaus in the velocity and the velocity dispersion, and minima in the higher moments h3 and h4. These are theoretically predicted by Bureau & Athanassoula (2005) and measured on several barred galaxies by Chung & Bureau (2004). We adopted the model of B17 as the one against which we compared our results. In their model, the bar has a position angle of PA bar = 55.7 • , extending out to 600 in projection. However, in the direction of the bar, our coverage only extends to 500 , we are therefore missing the crucial regions of the end of the bar. Pointings covering these regions are being observed and the results will be presented in a future paper. The co-rotation radius of B17 and the outer Lindblad resonance are also outside our observed regions.

Bar signatures in the stellar kinematics
In Fig. 26, cuts through the kinematic fields along the disk major axis are plotted. Left corresponds to the north-eastern re-  ceding side, right to the south-western approaching side. The velocity rises to 70 km s −1 at 35 , before remaining relatively constant until 350 on the left and 460 on the right, before reaching 160 km s −1 at −750 on the left and at 660 on the right. An axisymmetric stellar disk rises rapidly, but smoothly and remains flat at large radii. Bureau & Athanassoula (2005) find that a bar seen end-on produces a "double-hump" profile, where the velocity rises rapidly to a local maximum, then drops slightly to create a local minimum and then rises again slowly to the flat section at the end. When the bar is not seen end-on, but at an intermediate angle, this behavior is weakened and the local maximum and minimum disappear and form a single plateau at moderate radii, which is what we see in the profiles of M 31. This is expected from the models by B17, where the bar is seen at an intermediate angle. The rising part of the velocity curve comes from the orbits of the stars that are parallel to the bar, while the slower growth afterwards, or the local minimum in the case of an end-on bar, is caused by an inner ring structure caused by the bar (Bureau & Athanassoula 2005).  The velocity dispersion profile in Fig. 26 shows that σ has two off-centered maxima, with a slight drop of about 8 km s −1 in between. This drop was first observed by Kormendy (1988), who sees it as evidence for a central disk and subsequently the supermassive black hole. Bureau & Athanassoula (2005) find that such minima can also be caused by the bar, because the orbits parallel to the bar become more circular and thus lower the dispersion, however, the minimum in models by Bureau & Athanassoula (2005) is usually much stronger than what we measure here. One should keep in mind that such minima are not uniquely related to bars. Comerón et al. (2008) suggest that cold gas is accreted to the center and a subsequent starburst creates then stars with low velocity dispersion. They claim that the gas is most probably driven into the center by spiral arms, but also, albeit less probable, bars. An inner disk can also lead to this minimum, such a structure has been found in the very center of M 31 (Tremaine 1995;Peiris & Tremaine 2003), however, with a scalelength of only about 1 , this structure is not resolved by our observations. The drop could also be caused by the classical bulge (B17). σ drops to 140 km s −1 at 400 , before staying roughly constant out to 600 and then reaching 80 km s −1 at 950 . These plateaus, which are not completely flat, but only less inclined than the rest of the profile, are also frequently seen in barred galaxies (Bureau & Athanassoula 2005).
The h3 profile in Fig. 26 shows that the slope of h3 changes sign several times. The maxima and minima in h3 correspond to the points where the slope of the velocity profile changes, this is similar to the behavior of the simulations by Bureau & Athanassoula (2005) and has also been observed in barred galaxies (Fisher 1997;Chung & Bureau 2004;Falcón-Barroso et al. 2006;Ganda et al. 2006).
In Fig. 27, we show the correlation between h3 and v. In the outer regions, h3 is anti-correlated to v, which is the expected behavior for disky structures in a galaxy (Binney & Tremaine 1987;Bender et al. 1994;Fisher 1997;Binney & Tremaine 2008). In the central 30 , h3 is again anti-correlated to the velocity. This could mean that a disky structure is also present at the center, potentially explaining the slight drop in velocity dispersion. The radial extent of this central anticorrelation corresponds roughly to the rapidly rising part of the rotation curve, a behavior that is also seen in other disk galaxies (Chung & Bureau 2004) and interpreted as a decoupled inner disk. In between the two anti-correlated regions, h3 becomes correlated with the velocity, which is taken by Bureau & Athanassoula (2005) to be a sign for a bar. The correlated region is oriented more along the bar major axis than the disk major axis. It is more prominent in the south-west than in the north, because the northern edge of the bulge is more affected A38, page 12 of 22 Fig. 27. Correlation between the stellar velocity and h3. Plotted is sign(h3 · v), with the disk and bar major axis from Fig. 9. The color blue represents where the two quantities are anti-correlated, while in orange areas, they are correlated. by dust (Draine et al. 2014). The correlation means that there are more stars rotating faster than the circular velocity in projection, which may be a consequence of elongated motions. However, the correlation does not necessarily have to be caused by a bar, it can also be caused by the superposition of an axisymmetric bulge component embedded in a rotating disk, depending on the bulge-to-disk ratio. If the bulge is brighter than the disk, the observed mean velocity is going to be biased towards the bulge stars, with the disk creating a tail of high-velocity material (van der Marel & Franx 1993;Bureau & Athanassoula 2005).
The h4 profile in Fig. 26 is relatively constant, with the exception of a minimum at 670 , where h4 drops from 0 to −0.07. This minimum corresponds to the radius where h3 reaches its maximum. On the opposite side, the drop in h4 at −750 is not as pronounced. Outside of 1000 , the values of h4 stay roughly constant at larger values of ≈0.01 for positive radii and ≈0.07 for negative radii. B/P bulges often show dips in the very center in h4 (Debattista et al. 2005;Méndez-Abreu et al. 2008), however, this only applies to low inclinations of i < 30 • . It is therefore not surprising that we do not see a central drop in h4 in M 31.

Bar signatures in the gas kinematics and morphology
The simulations of Athanassoula & Beaton (2006) and B17 do not take into account the gas content in M 31. However, the gas content in M 31 is estimated to be only ≈ 7% of the stellar mass (Corbelli et al. 2010), so we can safely assume that the gas will follow the potential set by the stellar and dark matter components. A detailed discussion of the gas behavior in this potential is provided in Blaña et al. (in prep.). As mentioned above, the velocity maps for both gas components display an S-shape in the line of zero velocity. This S-shape is stronger than the twist in the stellar velocity field. Such S-shaped twists in the gas velocity are characteristic of velocity fields with oval distortions or bars (Bosma 1981). Many barred galaxies show them, like NGC 1068 (Emsellem et al. 2006), NGC 1300 (Peterson & Huntley 1980) Fig. 28. The orange and red points between 0 < r < 500 , 0 km s −1 < v < 100 km −1 and 100 < d < 200 correspond to the arm of zero velocity visible in Fig. 13. To the left of the center, there is an almost flat band of negative velocities, this is the zone of approaching velocities on the eastern side of the bulge in Fig. 13.
The region of very low velocities at the edge of the bulge on the near side of the galaxy could be produced by a large scale warp in the gas. Such a warp can project small velocities from further out into the line-of-sight (Melchior & Combes 2011).
In Figs. 28 and 29, position-velocity diagrams are plotted for the gas. We show the gas kinematics in this way to compare with similar diagrams from measurements and simulations of neutral gas, for example Chemin et al. (2009) and Athanassoula & Beaton (2006). Figure 28 shows the first component and Fig. 29 the second one. In this way, differences between the two components become immediately apparent. Both components occupy different regions in the PV-diagrams. The two components are real as they have different velocities. We also compared these PV-diagrams to the ones of the planetary nebulae from Merrett et al. (2006). As expected, neither gas Fig. 30. Cut through the velocity maps along the disk major axis within an aperture of 40 . The red line is a cut through the stellar velocity (Fig. 6), the blue one is a cut through the first velocity component (Fig. 12), the green one the second velocity component (Fig. 13). The black dashed line is the velocity for the triaxial bar model by Berman (2001). component corresponds to the PNe, since the PNe follow the kinematics of stars.
In the PV-diagram of the first gas component, there is a prominent steep band of velocities in the center. On the receding side, the velocities stay more or less constant once they reach the plateau. The second component also shows a similar steep band of velocities as the first component, but it is less pronounced and much wider.
The position-velocity diagram of the first gas component agrees very well with the overall shape of the model of Berman (2001), who simulates a triaxial bulge in M 31: the velocity reaches a peak, then stays constant and rises to a second peak. The comparison plot of his data to cuts through our maps is shown in Fig. 30. While generally similar, our position-velocity diagrams appear slightly more asymmetric: for the approaching side of M 31, the peak is closer to the center than for the other side.
The overall shape of the position-velocity diagram is very similar to what is seen in simulations ) and observations of barred galaxies Merrifield & Kuijken 1999), however, it could also be compatible with the superposition of other structures, like rings and disks, but this discussion goes beyond the scope of this paper.
The pattern that is visible in the flux maps in [O iii] and Hβ is similar to the one seen by Jacoby et al. (1985) in a Hα + [N ii] filter and in [O iii], by Boulesteix et al. (1987) in [N ii] and by Ciardullo et al. (1988) in Hα + [N i]. As seen in Sect. 4.3, the Hβ and [O iii] maps have the overall appearance of a spiral structure, with an incomplete elliptical ring inside, which is oriented roughly along the minor axis, with a spoke along the ring's short axis. The 100 long spoke along the disk major axis could be an inner disk that is projected to the observed elongated structure at M 31's inclination (Jacoby et al. 1985). The inner spiral pattern seems to be tipped to a lower inclination with respect to the outer part, which according to Jacoby et al. (1985) can be caused by a non-axisymmetry, like a bar, and cannot be explained by axisymmetric features alone.
While the ring structure lies in the region where B17 find two inner Lindblad resonances, it is hard to deproject the structures we see in our images. A detailed analysis of the structures is postponed to a later paper, where gas will be taken into account in the dynamical model.
The fact that the incomplete ring is aligned along the minor axis, almost perpendicular to the position of the assumed bar, leads Block et al. (2006) to conclude that this ring is not caused by a bar, but by the collision of M 31 with M 32, which also caused a split of the so-called 10 kpc ring, which is a structure appearing further out in the gas. While the collision model is in better agreement with the ring structure, especially with the fact that the rings are off-centered, it has low flux inside the ring, whereas in our flux maps, as well as the ones by Jacoby et al. (1985), Boulesteix et al. (1987) and Ciardullo et al. (1988), there is flux present there. We therefore think that the bar is a more likely explanation for the rings. However, as mentioned above, to fully reproduce the visible gas rings, we need to include gas into our bar models, which will be done in the future.

Bar or triaxial bulge?
The first models that tried to explain the triaxiality in M 31 (Stark 1977;Stark & Binney 1994;Berman 2001;Berman & Loinard 2002) called their structures a "triaxial bulge" instead of a "bar". These two concepts are fundamentally different: a bar is a rotating structure that results from an instability in the disk, while a triaxial bulge is a consequence of an anisotropy in velocity dispersion, generated, for example, by a merging processes. However, as already stated in Sect. 1, the models in the aforementioned papers actually describe rapidly rotating systems, which in our understanding are bars.
After buckling, most observed bars consist of a threedimensional inner part of the bar, the B/P bulge, and a flat outer part (see e.g., Athanassoula 2005or Martinez-Valpuesta et al. 2006. In this nomenclature, the triaxial bulge of the models by Berman (2001) and Berman & Loinard (2002) is a B/P bulge. Nevertheless, it is still interesting to see if our findings are also consistent with a completely non-rotating structure: Firstly, we attempt to quantify the level of cylindrical rotation in M 31, as cylindrical rotation is generally attributed to the presence of a bar ). Based on Saha & Gerhard (2013), Molaeinezhad et al. (2016, hereafter M16) define the parameter m cyl = m avg + 1 where m avg is a measure of the change in rotational velocity in the direction perpendicular to the major axis. m avg averages over the slopes of the decrease in rotational velocity within the bulge region at various distances to the minor axis. M16 first test the discriminatory power of this quantity on data from the Giraffe Inner Bulge Survey (Zoccali et al. 2014) of the Milky Way and on a barred N-body model of the Milky Way by Martinez-Valpuesta & Gerhard (2011), and then apply it to 12 further disk galaxies with inclinations ranging from 70 • to 90 • . They find that the B/P bulges in their sample as well as the Milky Way bulge show levels of cylindrical rotation of m cyl 0.7.
As M16 we derived rotation curves parallel and at varying distances to the major axis of M 31 first and then turned these into one-sided rotation curves parallel to the minor axis at varying distances to the minor axis. We performed this analysis on the south-eastern half of M 31 as it is this half where the disk is located on the far side, such that the bulge kinematics are not affected by dust screening of the disk. For the extent of the bulge along the major axis and along the minor axis we adopted values of 460 and 270 respectively from the bulge to disk decomposition in Kormendy & Bender (1999) and used the photometric disk position angle of 38 • (de Vaucouleurs 1958). The resulting  Molaeinezhad et al. (2016). Top: stellar line-of-sight velocity parallel to the major axis of galaxies at different heights from the disk plane. Bottom: line-of-sight velocity gradient parallel to the minor axis at different distances from the minor axis. Profiles are color-coded according to the distance from the minor axis. The solid black profile represents the average value m avg . The cylindrical rotation is then m cyl = m avg + 1 as a quantity to express the level of cylindrical rotation. plot is shown in Fig. 31. The parameter m cyl for the south-eastern side is 0.72. Changing the exact values of the position angle and also varying the radii for the bulge in major axis direction and in minor axis direction, and also the area of exclusion around the major axis (we use 50 ) had a minor effect on m cyl . Tests showed that for reasonable parameter ranges the value never becomes smaller than 0.68. Hence the level of cylindrical rotationas quantified by this method -seems to be consistent to what is reported in M16 for their B/P bulges.
Secondly, as described in Athanassoula & Beaton (2006) and B17, the spurs in the photometry outside the bulge argue for a planar bar, the "thin" or "long bar" component described in both papers. Although both papers give different lengths for the thinbar components, the fact still remains that such a component is required to explain the photometric appearance of the galaxy.
Thirdly, the outer 10 kpc ring of M 31 is associated with the outer Lindblad resonance of the bar by B17. This structure is only possible in a rotating barred potential, as there are no characteristic resonance radii in non-rotating triaxial potentials. In a system with a triaxial bulge, this ring would be an ad-hoc feature which is difficult to explain.
Finally, the correlation between v and h3 (see Fig. 27) is a clear signature for a bar and not a triaxial bulge.

Explanations for multiple components
Several scenarios have been suggested by various authors for the origin of the multiple gas components. When gas is subject to a non-axisymmetric potential like a bar, it leads to the formation of gas streams and rings, see for example Kim et al. (2012) for numeric models. The line-of-sight passes through several of these streams, therefore leading to multiple peaks in the gas lines. A detailed comparison will only be possible once gas is included into the dynamical model, which will be done in a future paper.
There is also an alternative explanation for such streams and rings in M 31, which is that they were created by the collision of M 31 with its small companion galaxy M 32 (Block et al. 2006). A galaxy with similar inclination as M 31, which has also a B/P bulge, is NGC 2683, which has been observed by Kuzio de Naray et al. (2009). Investigating the Hα velocity field, they see S-shaped twists and argue for the presence of a bar at a position angle of 5 • higher than the disk position angle. Differences in their PV-diagram and the ones for M 31 lead Melchior & Combes (2011) to conclude that the interpretation of Athanassoula & Beaton (2006) in terms of a bar in M 31 is wrong and that the ring structures are only due to the collision of M 31 with M 32 suggested by Block et al. (2006). However, the bar in NGC 2683 is only at an angle of 5 • to the disk major axis, so it is seen quite side-on, whereas the bar in M 31 is at an angle of 17 • to the disk major axis and it is seen more end-on, which can explain the differences according to Kuzio de Naray et al. (2009). Melchior & Combes (2011) also measure molecular emission lines in CO with two peaks, in Melchior & Combes (2016) sometimes even three. While the molecular gas is not expected to follow the ionized gas, it is still interesting that they also observe a line splitting. Melchior & Combes (2011) try to explain the double CO components with a tilted ring model coming from the collision model by Block et al. (2006). A small disk with inclination of 43 • -more face-on than the large-scale disk of M 31is surrounded by a ring. The velocity profiles extracted from the simulated velocity fields by Melchior & Combes (2011) consist of a very broad component, which is blueshifted and a narrow one, which is redshifted, the blueshifted part comes from the inner disk and the redshifted part from the ring. Our [O iii] measurements taken at the same positions look quite different, both components roughly have the same width, so the scenario of Melchior & Combes (2011) does not predict the actual shape of our measured [O iii] spectrum.
Another scenario for multiple components is that large-scale warps in the outer galaxy disk project lower velocities from the outer regions into the inner disk. In the neutral H i gas, Chemin et al. (2009) find up to five different gas components in M 31. The velocity maps for all these components look very similar and have basically the same regular kinematic pattern. A component with a steep slope in velocity in the center is the main H i component, while components with flatter slopes are due to slower gas in the outer H i disk that is projected into the central areas due to warps in the H i disk. However Athanassoula & Beaton (2006) conclude that the different slopes can be explained by the streaming gas motions caused by the bar. Comparing the position-velocity diagrams of Chemin et al. (2009) with our own, we found that their main H i component agrees with the first gas component that we observed while the The starburst grids are from Dopita et al. (2000), and assume a gas density of N e = 350 cm −3 and use a spectral energy distribution obtained from models from Starburst99 (Leitherer et al. 1999) for an instantaneous star formation episode. The grids assume a range of metallicity Z for both stars and gas in the starburst, shown by the solid lines for Z = 0.2, 0.4, 1.0, 2.0 Z , and different values for the ionization parameter q, shown by the dashed lines for q = 0.5, 1, 2, 4, 8, 15, 30×10 7 cm s −1 . The shock grids, without precursor H ii region, are from Allen et al. (2008) and assume an electron density of N e = 1 cm −3 and a solar value for the gas metallicity. The solid lines show models with increasing shock velocity v = 150, 200, 300, 500, 750, 1000 km s −1 , and the dashed lines models with magnetic parameter b = 0.5, 1.0, 2.0, 4.0. The grids were obtained using the program ITERA (Groves & Allen 2010). We assumed solar metallicity for the shocks and twice solar metallicity for the dusty active galactic nucleus (AGN) grids, because we measured these values in the bulge in our stellar population measurements, which will be presented in Saglia et al. (2018). zero velocity spiral arm in the second component agrees with warped components in H i. However, it should be kept in mind that the spatial resolution in Chemin et al. (2009) in the center is lower than our own.

Ionization mechanisms of the gas
In order to investigate which mechanisms are responsible for ionizing the gas, diagnostic diagrams are used, which compare the ratios of line fluxes of different emission lines. The most widely utilized of these diagrams compares [O iii]λ5007/Hβ to [O i]λ6300/Hα (Veilleux & Osterbrock 1987). Since Hα or [O i]λ6300 lie outside our observed wavelength range, we could not use this standard diagnostic diagram. Sarzi et al. (2010) devises alternative diagnostic diagrams for the SAURON spectrograph (Bacon et al. 2001 Fig. 32, the one for the second component in Fig. 33. There is not much difference between the two diagrams, the data points occupy similar regions in both. The ionization does not happen via starbursts, instead most of the points lie in the region where the ionization is due to shocks in the gas. That shocks are the primary source for ionization is supported by the fact that the emission seems to be strengthened along the lines of the velocity discontinuities. These shocks could be triggered by the bar. About half the points, however, also lie in between the Shock region and the AGN region, which is occupied by Seyfert galaxies in Sarzi et al. (2010), but since the AGN in M 31 is very weak (del Burgo et al. 2000), excitation by the AGN is unlikely. These datapoints lying between the regions can be partially explained by the fact that the errors of the datapoints are large.  Sarzi et al. (2006). Contamination by planetary nebulae is also not an issue, we checked the 166 bins that are affected by PNe. They lie in the upper parts of the shock grid and between the shock and the AGN grid, there is no systematic offset of the PNe positions with respect to the rest of our spectra.

Main results
For several decades, the question has been posed if M 31 is a barred galaxy or not (Lindblad 1956;Stark 1977; A38, page 16 of 22 Stark & Binney 1994;Athanassoula & Beaton 2006;B17). The high inclination of M 31 made it difficult for a long time to see the bar directly in optical photometry. Twists in the inner isophotes were the first hints for its existence (Lindblad 1956;Stark 1977), while more recently, the boxy appearance of the isophotes in infrared images was used to constrain its orientation and length (Beaton et al. 2007;B17). However, a kinematic confirmation of the bar has been missing so far.
In this paper, we have presented results from observations of M 31 with the optical integral field unit spectrograph VIRUS-W. We derived the kinematic properties of the stars and two different ionized gas components and measure the gas fluxes. We found hints in the kinematics that M 31 contains a bar. Our main results are: 1. Cuts through the stellar velocity field reveal a plateau at moderate radii, which is a signature that can be created by a bar that is seen at an intermediate angle (Bureau & Athanassoula 2005). The line of zero velocity does not coincide with the photometric minor axis and shows a twist toward the edge of the bulge. 2. The stellar velocity dispersion field shows a drop in the center, with two maxima aligned along the minor axis. While this central drop can be caused by inner disky structures or dust in the center (Falcón-  or the classical bulge (B17), it can also be due to the bar, since the stellar orbits that make up the bar become more circular in the center (Bureau & Athanassoula 2005). There are also two inclined plateaus at intermediate radii in the σ profile, a behavior that is also often seen in barred potentials (Bureau & Athanassoula 2005). 3. The Gauss-Hermite moment h3 is anti-correlated with the velocity v in the disk regions and the very center, which is the expected behavior for a disk component (Bender et al. 1994).
In the majority of the bulge region, h3 is correlated with the velocity v, which can be achieved by elongated motions which occur along the bar direction (Bureau & Athanassoula 2005). 4. The gas kinematics, measured using the [O iii]λ5007 line, is more complicated than the stellar kinematics. Many spectra exhibit two separate peaks, resulting in two kinematically distinct components. One of the two components has faster velocities than the other one. The fast component has disklike rotation with a very steep gradient in the center. The line of zero velocity is S-shaped, again pointing to a bar. The overall shape of the rotation curve of the first component is qualitatively in agreement with gas in a barred dynamical model of M 31 by Berman (2001). 5. The position-velocity diagrams of the gas components look similar to what is expected from simulations (Athanassoula & Beaton 2006), with arms of steep and flatter slopes. The overall shape of the position-velocity diagram agrees with observations in H i by Chemin et al. (2009). They measure up to 5 different components, the most luminous of which is the main H i disk, which coincides with the steep slope in the position-velocity diagram. The other components belong to the branches with lower slope in the position-velocity diagrams. According to Chemin et al. (2009), they are low velocities from the outer regions of the H i disk, which have been projected to the center due to warps. However, they can also agree with the lower velocity branches that exist in a barred potential (Athanassoula & Beaton 2006). 6. When looking at the morphology of the gas, we see a spiral pattern, similar to what is described by Jacoby et al. (1985) in Hα + [N ii] and [O iii], Boulesteix et al. (1987) in [N ii] and Ciardullo et al. (1988) in Hα + [N ii]. This spiral pattern has a lower inclination than the disk, which means it could have been tilted by a non-axisymmetry, like a bar (Jacoby et al. 1985). 7. Block et al. (2006) find that the morphology of M 31 in the far infrared is not caused by a bar, but is instead the result of a density wave caused by the collision of M 32 with M 31. Comparing their gas morphology with our measured one, the ring somewhat corresponds to what we see, but there is no high emission inside the ring, which we observe in [O iii]. Building on this model, Melchior & Combes (2011) propose a scenario of a tilted ring in the center over another rotating disk, to explain line splittings they measure in CO observations. Their model predicts that the component in the ring has a very narrow line, while the disk leads to a very broad line, which is not what we see, our two components in that region have the same width. We therefore think that this model is inconsistent with our observations. 8. From diagnostic diagrams for the ionization of the gas, we find that the gas is mostly ionized by shocks and not by starbursts, which is consistent with the low star formation rate of M 31 (Davidge et al. 2012). This is also in agreement by the filamentary appearance of the gas morphology, which could be either due to shocks or supernovae of type Ia (Jacoby et al. 1985).
To summarize, there are several hints in the kinematics of the stars and the ionized gas and its morphology for a bar in M 31.

Future work
From the analysis of the absorption line strengths, we measured the age, metallicity and α/Fe-overabundance of the stellar populations. These results will be presented in an accompanying paper (Saglia et al. 2018). These measurements, combined with the kinematics presented in this paper, are used as a basis for a madeto-measure dynamical model (Blaña et al., in prep.). This model will allow us to truly test the bar hypothesis and to check if our observations are caused by the bar, especially once gas is taken into account. They will also provide new estimates on the number of microlensing events caused by MACHOs or self-lensing, which will be used in future microlensing surveys toward M 31. We also plan to observe further pointings along the bar major axis, to calculate the pattern speed of the bar.

Appendix A: Comparison with data from S10
We compared the stellar and gas kinematics measured with VIRUS-W with the ones from S10. S10 observed M 31 with a longslit spectrograph along 6 directions, the slit positions are plotted in Fig. A.1. We generated cuts through our own stellar velocity maps along the slit directions of S10 and compared them with their data in Fig. A.2. Our velocities agree within 5% with the ones by S10, except for the minor axis on the near side (positive radii in Fig. A.2), where the deviation is about 30%. The comparison plot for the velocity dispersion is shown in Fig. A.3. Overall, our data agree well with the measurements of S10, but we do not reproduce the central spike in σ that is caused by the supermassive black hole, because we lack resolution in the very central regions. The overall agreement is within 4%. The comparison for h3 can be seen in Fig. A.4. The standard deviation of the difference between our values and the ones by S10 is 0.020-0.027, which agrees with the root mean square value of the error of the VIRUS-W values. In Fig. A.5, the corresponding cuts for h4 are shown. Again, we find good agreement, the standard deviation of the difference between our values and the ones by S10 is 0.02, which agrees with the root mean square of the error of the VIRUS-W values.
We also compared our gas kinematics to the data from S10. The comparison plot is shown in Fig. A.6. Generally, the first component agrees with the values from S10. For large radii, the velocity from S10 often lies between both velocity components. This is expected, the spectral resolution in S10 (≈69 km s −1 ) does not resolve the two gas components, so consequently they observe a single broader line at intermediate velocity. Similarly, the single component gas velocity dispersions of S10 are affected by the relatively lower spectral resolution and tend to be systematically larger than our values. Still, overall our data agree very well with S10.  Comparison of our velocities with the ones from S10. Black are the velocities measured by S10, red are cuts through our velocity maps along the same directions.   Comparison of our gas velocities with the ones from S10. Black are the velocities measured by S10, red are cuts through the first gas component from Fig. 12 and blue the cuts through the second component from Fig. 13. A38, page 20 of 22 Notes. The information about the stars was taken from the ELODIE (Prugniel et al. 2007) and LICK (Worthey et al. 1994) catalogs. The coordinates came either from the ELODIE catalog or van Leeuwen (2007), using the SIMBAD interface (Wenger et al. 2000).

Notes.
The full table is available at the CDS. Flux values are given with respect to the aperture of a single VIRUS-W fiber. The projected fiber diameter on sky is 3.2 . Hence to convert flux to per arcsec 2 quantities, the values have to be divided by the corresponding on-sky fiber aperture area of 8.04 arcsec 2 . Where multiple fibers are binned together, average flux values are reported such that the same normalisation holds true. To obtain the integrated flux for a specific bin, the values in this table must be multiplied by the number of respective fibers that are listed for that bin in Table B.4.