Exploring stellar and ionized gas non--circular motions in barred galaxies with MUSE

We present MUSE integral field stellar and ionized velocity maps for a sample of 14 barred galaxies. Most of these objects exhibit"S"-shape iso-velocities in the bar region indicative of the presence of streaming motions in the velocity fields. % By applying circular rotation models we observe that bars leave symmetric structures in the residual maps of the stellar velocity. %which demonstrates the capabilities of the MUSE instrument for detecting kinematic bar signatures. % We built non-circular rotation models using the \xs~tool to characterize the observed velocity fields. In particular we adopt bisymmetric models and a harmonic decomposition for a bar potential for describing the non-axisymmetric velocities. We find that both models reproduce the observed kinematic features. % The position angle of the oval distortion estimated from the bisymmetric model correlates with the photometric bar position angle $(\rho_{pearson} = 0.95)$, which suggest that non-circular velocities are caused by the bar. However because of the low amplitudes of the $s_3$ harmonic we can not rule out radial flows as possible source. % Because of the weak detection of \ha~in our objects we are not able to compare gas to stellar non-circular motions in our sample, although we show that when galaxies are gas rich the oval distortion is also observed but with larger amplitudes. % Finally, we do not find evidence that the amplitude of the non-circular motions is dependent on the bar size, stellar mass or the global SFR.


INTRODUCTION
Stellar bars are among the multiple non-axisymmetric components observed in galaxies (de Vaucouleurs et al. 1991), and it is estimated that nearly 1 3 of disk galaxies exhibit a stellar bar (e.g., Sellwood & Wilkinson 1993), with the number increasing when observing in infrared bands (e.g., Knapen et al. 2000). Commonly, bars are thought to play multiple roles in the evolution of galaxies including radial migration of stars (e.g., Minchev & Famaey 2010), triggering the star-formation (SF) and nuclear activity (e.g., Combes 2001) and for shaping the metallicity gradients in galaxies (e.g., Sánchez-Blázquez et al. 2011), among others. All these processes are intrinsically related with the dynamics of the bar. Therefore, understand their kinematic properties is crucial to unreveal their real influence in galaxy evolution. While bars are relatively easy to recognize in continuum images (Wozniak & Pierce 1991), their kinematic counterpart is not evident until detailed examination of their velocity fields. Even so, non-axisymmetric motions induced by bars have been identified in a wide variety of objects. These bar-driven motions when projected into the line of sight manifest in the form of oval distortions which have been observed in the velocity field of molecular (e.g., Weliachew et al. 1988), neutral (e.g., Bosma et al. 1977;Peterson et al. 1978) and ionized gas (e.g., Fathi et al. 2005;Holmes et al. 2015), as well as in stellar velocity maps (e.g., Kormendy 1983;. Such distortion recognized for producing a Sshaped pattern in the velocity field is less evident when the bar axis is closely aligned parallel or perpendicular to the line of nodes; these viewing angles disfavor the detection of bar stream motions (e.g., van Albada & Roberts 1981;Pence & Blackman 1984;Athanassoula & Misirio-tis 2002), making difficult to identify barlike-flows from the line-of-sight velocities. Hence it is common that barflow studies are biased toward galaxies with bars lying at intermediate orientations from the disk major/minor axis (e.g., Pence & Blackman 1984).
With the advent of recent observational techniques such as the integral field spectroscopy (IFS) we are able to spatially resolve their ionized and stellar kinematic properties (e.g., Fathi et al. 2005; Barrera-Ballesteros et al. 2014;Holmes et al. 2015;Fraser-McKelvie et al. 2020;Gadotti et al. 2020). In the optical, their kinematic counterpart have been identified mostly using Hα. This, however, tends to bias the studies towards gas-rich systems. In addition that Hα traces in most cases the location of young stars. Conversely, stellar bars are dominated by old stellar populations (e.g., Sánchez-Blázquez et al. 2014), whose ionization is mostly dominated by old stars (e.g., Stasińska et al. 2008;Gomes et al. 2016;Lacerda et al. 2018;Sánchez 2020).
On the other hand, studies of bars using the stellar velocity maps in most IFS-galaxy surveys are limited by inherent spatial resolution effects, precluding the identification of bar-driven flows (e.g., Barrera-Ballesteros et al. 2014).
In the present work, we address the kinematic study of bars by taking full advantage of integral field spectroscopic data from the multi-unit spectroscopic explorer instrument to detect bar-like flows on the stellar and Hα velocity maps of 14 galaxies. We built kinematic models to describe the bar flows and try to relate the amplitude of the bar non-circular motions with some global properties of galaxies.
The paper is structured as follows. In section 2 we describe the data; in section 3 the analysis of the velocity maps and the kinematic models adopted; in section 4 we describe the results and the conclusions are presented in section 5. Throughout this paper we adopted a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7.

DATA
In this study we use the dataproducts from the AMUSING++ galaxy compilation (Galbany et al. 2016;López-Cobá et al. 2020). AMUSING++ is a compilation of ∼ 600 nearby galaxies observed with the multi-unit spectroscopic explorer (MUSE) instrument (e.g., Bacon et al. 2010). With a field-of-view (FoV) of 1 × 1 , delivering ∼90 K spectra per data-cube, MUSE is the most advanced integral field spectrocopy (IFS) instrument in the optical range covering from 4800Å to 9300Å. It combines both a moderate spectral resolution (R ∼ 3000), and a high spatial resolution limited by the atmospheric seeing.
The AMUSING++ cubes were analyzed using the pipe3d package (e.g., Sánchez et al. 2016). This tool has been extensively used on a number of IFS galaxy surveys such as CALIFA (e.g., Sánchez et al. 2012), MaNGA (e.g., Bundy et al. 2015) and SAMI (e.g., Allen et al. 2015). pipe3d models the stellar population using a combination of simple stellar populations (SSP), and derives the properties of the emission lines once subtracted this model to the original spectra.
On the other hand, the recovery of the stellar population properties, such as the stellar kinematics, is performed over a tesselation map that follows the galaxy light distribution, typically around the V-band. In general, the size of these tesselations or voxels depend on the target signal-to-noise (S/N) of the continuum, here chosen to reach S/N≥ 30. For the particular case of the MUSE data, these correspond to ∼ 1 for regions with high S/N, and ∼ 2 for regions with low S/N, typically for the outskirt of the disks. We emphasize that this is one of the major advantage of the current data compared to other IFS-surveys with coarser spatial resolutions. Then the spectra within each voxel are coadded and over this spectrum is performed a simple stellar population analysis (SSP) to recover the velocity, age of the SSPs, metallicity, among other properties. This procedure is performed over each voxel of the tesselation map.
In addition to the above analysis, the ionized gas properties such as the flux, velocity and velocity dispersion, are recovered spaxel-by-spaxel by implementing a moment analysis on a set of ∼ 50 emission-lines. We refer the reader to the AMUSING++ presentation paper (e.g., López-Cobá et al. 2020) for more details of this procedure. The final outputs of the pipe3d analysis are two dimensional maps of the main properties of ionized gas and the underlying stellar populations.
For the current analysis, we selected a sample of barred galaxies from the AMUSING++ dataset. Galaxies were required to exhibit clear stellar bars in the MUSE gri images. Barred galaxies in interaction were excluded since their kinematics is dominated by the interaction process and not by internal ones; in addition, we exclude those objects where the higher S/N voxels have sizes larger than 3 times the seeing of the observing night (based on the ESO DIMM archive). This condition excludes galaxies observed with non-optimal atmospheric conditions. Finally the apparent bar length must be resolved in the stellar velocity map and it must fit in the MUSE FoV. These criteria lead to a sample of 14 galaxies. Since AMUSING++ is a compilation of Yellow contours represent the isophotes. Bottom panel shows the radial variation of the disk position angle and ellipticity of the isophotes. The vertical red line shows the approximate length and position angle of the bar, estimated in r bar ∼ 9.5± 1.5 and φ = 114 ± 3 • respectively. The red shadow region represents the bar length standard deviation. The green line ontop of the continuum image represents these values. The galaxy is oriented North-East with North pointing up and East to the left. objects, our sample of bars is not an statistical representation of the population of barred galaxies in the Local Universe. Therefore, the results of this analysis are not statistical significant.

Photometric properties of bars
Stellar bars are dominated in general by old stellar populations (e.g., Sánchez-Blázquez et al. 2014;Vera et al. 2016;Fraser-McKelvie et al. 2019;Neumann et al. 2020); therefore, it is common to study these structures in the reddest bands of the optical spectrum or in infrared bands (e.g., Díaz-García et al. 2016), where their morphological properties are enhanced. Among those properties is the bar-length, albeit there is no unambiguous method for estimating it (see Athanassoula & Misiriotis 2002, for a thorough description of methods). Parametric methods often decompose the galaxy light distribution into axisymmetric and non-axisymmetric components to recover the disk, bulge and bar properties (e.g., Laurikainen et al. 2018;Méndez-Abreu et al. 2019). An easier way consists in analyzing the object light distribution by tracing isophotes. Abrupt changes in the disk position angle (φ disk ) and the axial ratio ( ) have been commonly used to estimate the bar length, as well as its orientation in the sky (e.g., Kormendy 1983;Pence & Blackman 1984;Wozniak & Pierce 1991;Pérez et al. 2009). Then, the true length of the bar along the major axis can be estimated by pure trigonometric relations as follows: where r bar is the apparent length of the bar in the sky, φ = φ disk − φ bar,phot , with φ bar,phot representing the bar position angle, and i is the disk inclination. Instead of adopting narrow band images from the MUSE cubes, we use r−, i− or z− band images from DES (e.g., Abbott et al. 2018) or Pan-STARRS (e.g., Chambers et al. 2016) whenever available. These images have a larger FoV and they are deeper than our MUSE data. For extracting the isophotes we use the photoutils Python package (Bradley et al. 2016), which relies on the ellipse fitting analysis introduced by Jedrzejewski (1987). After background subtraction, we adopt the last faintest isophote to derive the galaxy position angle φ disk and disk inclination. Then we trace their radial variations, and measure the bar position angle and bar length based on the three consecutive isophotes that maximize the difference ∆P = P i+1 − P i , where P is φ disk or the ellipticity ( ). Although this analysis is just a proxy for the bar size and its orientation, it offers a first order estimation of these properties. As a sanity check, we performed a visual comparison between our estimations of φ bar,phot and r bar with their apparent projections in the continuum images. In all cases a mutual agreement is reached. Figure 1 shows the implementation of this procedure for one object in our sample, ESO476-16. The estimated photometric bar position angle for this object is φ bar,phot ∼ 114 • , which is consistent with the apparent orientation of the bar in the sky. We finally applied this analysis to all our galaxies in our sample.

Kinematic signatures of bars
A common and noticeable effect of bars on the velocity field is the change in orientation of the dynamical major and minor axes (e.g., Barrera-Ballesteros et al. 2014), a signature of bar stream flows. As mentioned before, oval distortions on the velocity field of barred galaxies are well known from the gas velocity fields (CO, H I, Hα). These kind of perturbations induce large deviation of the circular motions which can be observed in residual maps of simple axisymmetric models of circular rotation (e.g., . Thus, we adopt kinematic models to identify non-circular motions most likely induced by bars.

Kinematic models
Multiple algorithms have been developed for describing the velocity field of galaxies. The vast majority of them are based on the tilted ring model (e.g., Begeman 1989), in which the velocity field is divided in concentric rings each one rotating with a different velocity around a fixed kinematic centre. In this work we use XookSuut (e.g., Lopez-Coba et al. 2021), which relies on the DiskFit (e.g., Spekkens & Sellwood 2007) and RESWRI (e.g., Schoenmakers et al. 1997) algorithms. XookSuut performs non-parametric circular and noncircular rotation models for describing the line-of-sight velocity (LoS). It adopts Markov Chain Monte Carlo (MCMC) methods for sampling the posterior distribution of the different kinematic components and creates an interpolated model based on the highest probability states.
XookSuut adopts the following log likelihood with flat priors on the parameters: where α represents all the parameters describing the considered model; V obs is the observed velocity map; σ is the error map of the measured velocities; W k,n is a series of weights computed at k independent rings and will serve for creating a two-dimensional interpolated model V model ; N is the total number of pixels included in the model. All XookSuut models adopt the thin disk approximation; i.e., it assumes the disk is intrinsically flat, with a constant disk ellipticity, constant position angle and fixed kinematic center. These conditions make it possible to represent kinematic models on twodimensional maps (hereafter 2D model).
The circular model is the simplest model included in XookSuut, which is described by the following equation: In this expression V t is the circular rotation, which is function of the radius. V sys is the systemic velocity assumed constant for all points in the galaxy. The azimuthal angle θ is measured on the galaxy plane and is related to the sky coordinates through the inclination angle (i), kinematic position angle (φ disk ) measured from North to East from the galaxy receding side, and the kinematic centre (x c , y c ).
In order to characterize the non-circular motions in our sample, we adopt the following non-circular models included in XookSuut. The first one is the bisymmetric model proposed by Spekkens & Sellwood (2007), and is suitable for describing non-circular motions induced by a bisymmetric distortion to an axisymmetric potential; for instance that caused by a stellar bar. This model fits elliptical stream lines around a fixed angle to the LoS velocities; the equation describing this model is given by: In this expression the V 2,r and V 2,t terms represent the bisymmetric deviations (radial and tangential respectively) from the circular velocity represented by the V t term. The phase φ bar , represents the position angle of the oval distortion with respect to the azimuthal angle. Its sky projection is given by the following expression (e.g., Bettoni & Galletta 1997;Spekkens & Sellwood 2007) : where φ disk is the kinematic position angle major axis of the receding side, and i is the galaxy inclination. The description of the non-circular motions in our second model is based on the epicycle theory, where the radial and tangential components of the non-circular motions are taken into account by inducing small perturbations to circular orbits. As shown by Schoenmakers et al. (1997), the LoS velocity V los , can be expressed as a Fourier series as follows: where c 0 is the 0 th harmonic term, assumed to be constant here (∼ V sys ). The harmonic coefficients c m and s m describe different velocity components as a function of the galactocentric radius, and θ is the azimuthal angle measured in the galaxy plane. The black rectangle on top of (iii) shows the region with the largest residuals. Iso-velocity contours spaced each 50 km s −1 are superimposed on each map. Bottom panels have the same meaning but this time for the Hα velocity map. Note the zero velocity line shows an "S"-shape distortion. Schoenmakers et al. (1997) showed that a perturbation to the potential of order m affects only the m = m ± 1 harmonics in the Fourier expansion (Eq. 6). For the case of a bar m = 2 and the harmonics describing the LoS velocity adopt the following expression (e.g., Franx et al. 1994): V los = c 0 + sin i c 1 cos θ + c 3 cos 3θ + s 1 sin θ + s 3 sin 3θ (7) Thus, this expression represents the expected behavior of the LoS velocity for an elongated potential. Note that a bisymmetric perturbation also includes the m = 3 and m = 1 harmonics (e.g., Spekkens & Sellwood 2007;Oman et al. 2019); however, as we will see in further sections the amplitude of the harmonic coefficients allow us to assess the source of the non-circular motions.
We use XookSuut to built circular and non-circular flow models of the stellar and Hα velocity maps of our sample of galaxies.
Unlike RESWRI (Schoenmakers et al. 1997), XookSuut does not allow the kinematic center and projection angles to change across radii. In fact, this behavior is not desired in our analysis since we expect that non-circular motions arise due to the bar potential, rather than being induced by the presence of a twisted disk.
XookSuut requires guess values for the disk position angle, inclination and kinematic cen-ter. For the first two we use the values estimated from the isophotal analysis described in Sec. 3.1, while the kinematic center was estimated by eye in each velocity map. Since interpolated models are created over rings on the disk plane, we choose the rings to be spaced each 2 . This value is grater than the spatial resolution of our objects (FWHM DIMM ∼ 1.5 on average). Finally, we use the corresponding error maps for discarding spaxels with low S/N; that is, we remove spaxels with errors larger than 10 km s −1 .
In the following we use the same galaxy ESO476-16 as a showcase to go through the different kinematic models. Figure 2 shows the circular rotation model on the stellar and Hα velocity maps. This figure highlights the differences in spatial resolution in both maps. Given that the recovery of the gas kinematics does not involve the binning procedure described in Sec. 2, the Hα velocity map shows a better spatial resolution than the stellar velocity.
The stellar velocity reveals a central distortion induced most probably by the presence of the bar, while the iso-velocity contours appear twisted near the minor axis. The stellar circular model shows a typical rotating disk with orthogonal major-minor axis. The superimposed iso-velocities show that the outermost regions in the galaxy are compatible with a pure rotating disk. The residual map in both cases (stars and Hα) show several important features in the inner parts of the disk. A close-up around the minor axis shows symmetric structures with large residuals of the order ±40 km s −1 . On the other side, the Hα velocity field shows a central ring with no data because of the low S/N of Hα in this region. Despite of that, XookSuut is able to predict values in these regions by linear interpolation. We notice that the previous structure with high residuals is also observed in the gas velocity field, albeit is affected by the low S/N data. This is more evident in the close-up to this region in the third panel. Such symmetric structures with blueshifted and redshifted components have been observed in residual maps of barred galaxies with similar amplitudes (e.g., Fathi et al. 2005;Castillo-Morales et al. 2007). Thus, we may conclude from figure 2 that a non-axisymmetric component is present in the residual velocities from both the ionized gas and stars.
Whether the velocity field of a barred galaxy is better described by a bisymmetric flow or by a pure axisymmetric radial flow is in general not straightforward to determine (e.g., Wong et al. 2004). In this work we will not make any a priori assumption on the kinematic model. Instead, we first adopt a bisymmetric model constraining the radial extension of the non-circular motions up to the the deprojected bar size plus an additional value ranging from 1 -2.5 , which takes into account the spatial resolution of the data. Additionally, we decompose the LoS velocity with a Fourier analysis including only the m = 1 and m = 3 harmonics (Eq. 7), since the ratio of the s 3 and s 1 coefficients can be used to distinguish between bar-like flows and radial flows. (e.g., Wong et al. 2004;Fathi et al. 2005;Elson et al. 2011).
The inclusion of complex models such as the bisymmetric model adds extra variables in the fitting procedure that may reduce the residuals, but at the expense of overfiting the data. In order to assess the bisymmetric model we compute the Bayesian information criterion (BIC, Kass & Raftery 1995), defined as BIC = N ln(χ 2 /N ) + ln(N )N params , where χ 2 is the quadratic sum of the residuals and N params the number of parameters to estimate from the model . Unlike χ 2 , BIC is more sensitive to the number of parameters to estimate from the model, favoring less complex ones.
We use the outputs of the circular rotation model (namely x c , y c , i, φ disk , V sys ), as input for the harmonic and bisymmetric models. We note that there is no significant differences in the results if we fix these parameters or we let them vary during fitting. In any case, we let XookSuut to derive the best set of parameters. Results adopting the bisymmetric model on the stellar velocity map of the showcase object is shown in Figure 3. This figure is a corner plot for the parameters describing the disk geometry and the bar orientation. The bisymmetric model finds an oval distortion oriented at φ bar = 43 • , with its corresponding sky projection (Eq. 5) oriented at φ bar,kin = 289 ± 1 • . The position angle of the oval distortion is in agreement with our previous estimation for the photometric bar position angle (114 The 2D representation of the bisymmetric and harmonic models together with the radial profile of the different velocity components are shown in Figure 4. In the top panels, the bisymmetric model reproduce successfully the "S"shape distortion observed in the inner regions, furthermore the residual map no longer exhibits the symmetrical patterns observed in Figure 2. On the other hand, the non-circular velocities, V 2,r and V 2,t , show maximum amplitudes of the same of order as the residuals in the circular rotation model, i.e., ∼ 40 km s −1 , with smooth profiles as those observed in gas velocity fields (e.g., Sellwood & Sánchez 2010;Holmes et al. 2015).
The bottom panels of figure 4 show the results of the harmonic decomposition. Again the 2D model seems to be a good representation of the LoS velocities, and no residual structures are observed in the central regions. The similar behavior of the c 1 harmonic with V t in the bisymmetric model, suggesting that the disk circular rotation is indistinguishable between both models and that the differences must reside only in the amplitudes of the non-circular components. The inclusion of the m = 3 terms reproduce much better the features observed in the velocity field, however the residual map is, at first sight, indistinguishable from the bisymmetric model indicating that the harmonic decomposition produces results as good as the bisymmetric one. Nevertheless, the interpretation of the harmonic coefficients is not straightforward. Based on the epicyclic theory, Franx et al. (1994) and Wong et al. (2004) computed the expected behavior of the m = 1 and m = 3 harmonics for an elliptical potential (i.e., an m = 2 perturbation). They found that for an elliptical potential, the slope between the s 3 and s 1 coefficients (hereafter ds 3 /ds 1 ), is found to be negative, i.e. ds 3 /ds 1 < 0, meanwhile for an axisymmetric radial flow |ds 3 /ds 1 | < 0.1, which means that s 3 s 1 . The bottom rightmost panel of figure 4 shows the s 1 /c 1 vs. s 3 /c 1 ratio together with the best fit line to the points. For the considered galaxy, the slope of the non-circular velocities is found in −0.11; there- fore, the s 1 and s 3 harmonics are compatible with being produced by a bar potential.

RESULTS
So far we have described our methodology for a single object. In the following we describe our results when applied to our sample of galaxies. Prior the kinematic modeling we removed spaxels with errors larger than 10 km s −1 to exclude low signal-to-noise data that could affect the final models. For all objects we first create circular models with initial values taken from the isophotal analysis. The 2D models for each individual object are shown in Appendix A. Then we use the best fit values of the geometric parameters (namely φ , i, x c , y c and also V sys ) as inputs for the bisymmetric and harmonic decomposition models.
It is expected that bar-like flows are not present across the entire disk, but only in the region influenced by the bar as observed from the residuals of Figure 3. Therefore, we fit non-circular motions up to the size of the photometric bar, plus a constant value that goes from 1 − 2.5 to compensate for the binning on the stellar velocity maps. The two-dimensional non-circular models for the stellar and gas velocity maps are shown in Figs. 10 and 11 from Appendix B. Table 1 shows XookSuut results for the constant parameters namely φ disk , i and V sys . The phometric esti-mations of the bar position angle and true bar length (i.e., Eq. 1) are shown in the third and fourth columns, respectively. In general, we do not find large differences in the constant parameters, when adopting a bisymmetric or an harmonic decomposition model. This means that the flat disk assumption is adequate for our objects, besides that both fits are numerically stable despite the underlying assumptions. Exception to this might be IC 0004 where a constant ellipticity by XookSuut does not seem to reproduce the kinematic orientation of the disk. Table 1 also shows the BIC value of the bisymmetric model, but normalized to the BIC circular model. Thus, models with BIC > 1 would not favor the bisymmetric model. We note that only in one object (PGC 055442) BIC definitely favors the circular model over the bisymmetric one. Figure 9 shows the 2D circular models of the stellar and Hα velocity maps, as well as the residual maps of the models. We notice that in the vast majority of cases Hα is not detected in the inner disks; however, a careful analysis of the residual velocities in both stellar and Hα, large scale residuals are observed around the bar region in most cases, for instance in NGC 692. Furthermore when Hα is detected, symmetric residual structures with opposite velocities are observed in the inner regions; for instance IC 2160, IC 0004, PGC 055442 and ESO 18-18. Other objects although show high residual towards the center, they appear to be affected by other sources of non-circular motions, such as spiral arms as in NGC 289 and NGC 3464.

Circular rotation models
On the other hand, the residuals of the stellar velocity maps show many examples of galaxies with "S"-shaped iso-velocities together with symmetric residual structures around the bar. Such residuals are not compatible with an axisymmetric rotating disk and are most probably induced by the presence of the stellar bar. The fact we are clearly detecting oval distortions in most of our stellar velocity maps is likely a consequence of the high spatial resolution of our data, since similar IFS studies with coarser resolution do not observe such structures (e.g., Barrera-Ballesteros et al. 2014). Figure 10 shows the bisymmetric and harmonic decomposition models for the stellar velocity maps. Bisymmetric models tend to reproduce much better the observed velocity map than circular rotation models do. For most objects, V 2,r and V 2,t show smooth behaviors, with velocities decaying at large radius, a possible signal of the kinematic ending of the bar.

Stellar non-circular motions
NGC 289 is the only object where the bisymmetric model, for both gas and stars, does not seem to be a reasonable physical model given the erratic behavior of the circular and non-circular velocities. As observed in figure 10 and table 1, the bar in NGC 289 is closely aligned with the major axis, being separated by less than 10 • ; and it is well known that bar stream lines in the velocity field become less evident when the bar lies close the galaxy major-axis (e.g., van Albada & Roberts 1981). Moreover, when φ bar = 90 • or 0 • , V 2r and V 2t from Eq. 4 become degenerate with the disk circular rotation, making it impossible to separate the non-circular velocities from the circular rotation. As a result, V t , V 2r and V 2t can lead to unphysical absurd values as is probably the case in NGC 289; for instance, note that for r < 40 the three velocities shows the same uncommon radial profiles. On the other hand the harmonic model does not depend on the bar phase; thus, the disk circular rotation related to c 1 seems to be offer a better representation of the rotation curve in this case.
When comparing the BIC ratio from the bisymmetric to circular models we note that, in the case of the stellar maps, these values are mostly close to 1, with the aforementioned exception. We expect BIC ratios lower than 1 would favour the bisymmetric model. While this is observed in our sample, in some cases the differences with the circular model are only marginal. It should be noted that because of the binning procedure during the SSP analysis, spaxels are correlated in the stellar maps, affecting the estimation of the BIC value. Therefore, we favour the bisymmetric models in those cases where we detect kinematic oval distortions, regardless of the BIC ratio.
The bar position angle is the main parameter in the bisymmetric model since it defines the axis where the non-circular motions are expected. To visualize better this parameter, Figure 5 shows the photometric bar position angle vs. the sky projected kinematic bar position angle for the objects in the sample. It can be seen that there is a good agreement between both angles; in deed, the Pearson correlation coefficient for this relation is extremely high, being ρ pearson = 0.95.
On the other side, the harmonic decomposition models show good representations of the LoS velocities. Table 1 shows that at least in one velocity map (stellar or ionized), the slopes in the harmonic coefficients ds 3 /ds 1 are negative, which is in favor of the hypothesis that the non-circular motions are induced by the bar. However there are objects where |ds 3 /ds 1 | < 0.1, which means that radial flows could be contributing to the observed non-circular motions. Figure 6. Comparison between ionized gas and stellar bisymmetric models for IC 2160. Black, blue and red curves represent the Vt, V2,r and V2,t components for the Hα (continuous line) and stellar (dashed line) velocity maps.

Ionized gas non-circular motions
As bar-like flows are expected along the bar, the lack of gas in the inner regions affects directly the estimation of φ bar , and hence the amplitude of the non-circular motions. Because of the weak emission of Hα in the inner disks in our objects, the bar region provides little or null information of the behavior of the non-circular motions. Only in 6 objects is observed plenty of gas in the inner disk to be able to perform reliable non-circular flow models. These objects are shown in Fig. 11. However, the bars in NGC 289 and NGC 3464 are closely aligned parallel to the disk position angle, preventing the bar streams from being separated from the disk circular rotation. Even so, the large residuals observed in fig. 9 indicate the presence of strong non-circular motions, al-though a large fraction of them could be attributed to the prominent spiral arms in NGC 289 (e.g., Pence & Blackman 1984).
We notice however that oval distortions appear in the Hα velocity maps when galaxies are gas rich; and such distortions are pronounced when the bar is elongated at an intermediate position angle from φ disk . A remarkable example is observed in IC 2160 where the bar is oriented ∼ 40 • away the kinematic position angle. As consequence, a strong twist in the kinematic major and minor axes is observed.

DISCUSSION
5.1. Bar-like or radial flows? Wong et al. (2004) categorized different source of noncircular motions based on the slope of the s 1 and s 3 coefficients. For radial flow dominated systems, they found that |ds 3 /ds 1 | < 0.1. Table 1 show objects that fall in this category; however, for describing barred galaxies the bisymmetric model is usually preferred over radial flows. In this model the perturbation is driven along a fixed axis (φ bar ) and the non-circular velocities correspond to those of elliptical orbits. In the epicyclic theory, radial stellar motions are not expected to contribute significantly to the non-circular motions (e.g., Sellwood & Binney 2002). In addition, the large speeds on the observed non-circular velocities could not persist for a very long time without rearranging the mass distribution of the disk (e.g., Wong et al. 2004;Spekkens & Sellwood 2007); therefore, invoking pure axisymmetric radial flows (inflows/outflows) brings considerable consequences in the stability of stellar disks.
In addition, numerical simulations also show that inflows in bars are expected to be < 5 km s −1 (e.g., Athanassoula 1992), but not of tens of kilometers per second as observed in most objects. Recent studies however, observe radial motions in spirals of the order of 10-30 km s −1 (e.g., Di Teodoro & Peek 2021). Such large speeds of gas inflow/outflow would invoke necessary efficient mechanisms of gas depletion through star formation processes or galactic scale winds to remove gas out of the galaxy. The former scenario, although it has been widely observed (e.g., López-Cobá et al. 2019, 2020, the link of axisymmetric radial flows with out of plane outflows is a subject that has not yet been addressed. For all the above reasons pure axisymmetric radial flows are infrequently considered for describing the non-circular flows in barred galaxies.

Deficit of ionized gas in bars
The absence of ionized gas in the central regions observed in halve of our sample, suggests it could be related with the presence of the bar itself. In fact, in most barred galaxies there is observed a lack of gas in the inner disk (e.g., Erroz-Ferrer et al. 2015;Fraser-McKelvie et al. 2020). However there is no clear explanations in the literature for why bars tend to not exhibit ionized gas. Two scenarios are related with the absence of molecular gas. If there is absence of molecular gas along the bars, no new stars are formed. Closely related might be the presence of, indeed, radial flows that could have transported cold gas towards the center. In such case one should expect a higher star formation rate (SFR) towards the center, which is not yet observed (e.g., Erroz-

Ferrer et al. 2015)
. It is therefore important to trace the cold gas abundance in bars to help to distinguish between both scenarios.

Connection between non-circular motions and galaxy properties
As noted in Figure 10 and 11, the amplitudes of the non-circular velocities V 2,r and V 2,t vary widely from galaxy to galaxy. Thus, we ask whether the amplitude of the non-circular motions depend on some galaxy properties or they are local processes. For instance recent studies show a two-slope relation between the stellar mass and the bar length Erwin (2019). Since bars are composed largely of old stellar populations, the bar must contribute to a large fraction of the total stellar mass, as observed in some barred galaxies (e.g, Sánchez-Blázquez et al. 2011). Thus the relation between r bar and the host galaxy stellar mass, should be expected as well with the bar-mass. Figure 7 shows the deprojected stellar bar length (Eq. 1) versus the stellar mass obtained from the SSP analysis. We notice that our objects fall around the steeper slope of this relation where galaxies with larger bars are located. Indeed most of our objects have deprojected bar sizes larger than 2.5 kpc.
We investigate if the stellar mass, size of the bar and the star formation rate (SFR) are coupled with the amplitude of the non-circular motions. Instead of adopting the average value of the circular rotation residuals, we characterize the non-circular motions with the amplitude of the bisymmetric components V 2,r and V 2,t . As these terms are function of radius we define their am-   Table 2. Pearson and Spearem correlation coefficients between non-circular motions and galaxy global properties. The first column from top to bottom correspond to the fraction of non-circular over the local circular rotation, the maximum amplitude of the non-circular motions and the average velocity of non-circular motions. These parameters are compared with the galaxy stellar mass, the size of the bar and the SFR integrated across the MUSE FoV. For each set of parameters we have computed the Pearson correlation coefficient and the Spearman rank correlation with p-values shown in parenthesis.
plitud, A bis , as the quadratic sum of both components: A bis (r) = V 2 2,t (r) + V 2 2,r (r) and the fraction of non-circular over circular motions as: f nc (r) = A bis (r)/V t (r) for r ≤ r bar (9) where V t is the tangential velocity component of the bisymmetric model, which gives a better description of the circular rotation than if we consider the pure circular model. Note again that these two parameters depend on the galactocentric distance. Thus in the second expression we are comparing at each distance the non-circular motions with the local circular velocity. Figure 8 shows A bis and f nc computed for the stellar velocity maps against the stellar mass and the bar lengths of our objects. In this figure we also include the integrated SFR derived with Hα, with the dust attenuation correction using the Cardelli et al. (1989) extinction law and case B of recombination (e.g., Osterbrock 1989), and use the distances reported in the AMUSING++ paper (e.g., López-Cobá et al. 2020). The first thing to notice is the strength of non-circular motions ranges between 10-50% of the local circular rotation on average, which is similar to what other studies have found (e.g., Bettoni & Galletta 1997).
In figure 8 we tried to find possible correlations between kinematic properties and global properties, while Table 2 shows the Pearson and Spearman coefficients for those relations. We find however only weak correlations between these parameters. Even more, the large p-values show that such correlations are not statistically significant. These results are in line with recent studies (e.g., Erroz-Ferrer et al. 2015), where no obvious relation is observed with the residuals of circular rotation and morphological properties of bars. However, we are aware of the low statistic provided by this sample.

CONCLUSIONS
In this study we analyzed the incidence of non-circular flows in the stellar and Hα velocity field of a sample of non-interacting barred galaxies observed with the MUSE spectrograph. The exquisite resolution of the data allowed us to detect oval distortions in the stellar velocity maps most probably associated to the presence of the bar. Such perturbations are recognized in the residual maps of pure circular rotation models as symmetric structures with blueshifted and red-shifted velocities, clearly evidencing the presence of a non-axisymmetric potential affecting the expected disk rotation. Despite not observing ionized gas in the inner disk in most of the objects, the bar leave imprints in the circular residual velocities in regions close to the bar. When ionized gas is observed along the bar, the oval distortion is also revealed in the Hα velocity field. This evidence that both stars and gas are affected by the bar potential.
We characterize the kinematics of these galaxies with models that include non-circular motions, in particular the bisymmetric model and a harmonic decomposition for a m = 2 potential perturbation. We use the slope of the harmonic coefficients, the ds 3 /ds 1 ratio, to determine the source of the observed non-circular motions.
Based on this parameter we find that the harmonic decomposition on the LoS velocity is compatible in most cases with being produced by an elongated potential. However, some radial flows could be present in some objects since |ds 3 /ds 1 | < 0.1 is observed in 4/14 cases. In these cases, we attribute to the dominance of other sources of non-circular motions such as spiral arms.
We find that the bisymmetric model produces successful fittings in the velocity maps within the bar region. We also find that the position angle of the oval distortion correlates with the photometric position angle of the bar. This supports the scenario of a stellar bar potential inducing deviations of the circular rotation in the observed velocity fields. We find that the lack of ionized gas along the bars restricts the robustness of a modeling of the non-circular motions, in which case the stellar velocity could be considered for studying bar-like flows when no gas is detected. However we notice that not only the stellar circular rotation is affected by the asymmetric drift, but also the non-circular motions since gas velocities appear rotate faster.
We find the average amplitude of the non-circular motions in our sample is ∼ 30 km s −1 for stars and ionized gas, while the strength of the non-circular motions reaches values of up to 50% of the local circular velocity, although this fraction varies among galaxies.
When trying to relate the non-circular motions with galaxy properties we do not find any clear correlation with the stellar mass, the bar size or the global SFR. These results point that bar-flows are rather local process; however, we stress that larger statistical samples of barred galaxies with high spatial resolution are required to reveal possible correlations between kinematic properties of bars and global properties of galaxies.
We would like to thank the anonymous referee for their comments and suggestions that contribute to improve the quality of this paper.
L Following are shown the maps of the circular and non-circular kinematic models for the sample of galaxies shown in Table 1. In each figure, only one object is shown as an example, the remaining figures are available in the online journal. Figure 9 shows the circular rotation models derived by XookSuut for the ionized-gas and stellar velocity maps of each object. The three figures on the left (right) correspond to the observed velocity map, best kinematic model and the residual map for the ionized (stellar) velocity maps.  Figure 10 shows the results from the isophotal analysis together with the kinematic modeling of XookSuut adopting bisymmetric and harmonic decomposition models. In this set of figures, only one object is shown as an example, the remaining figures are available in the online journal. This figure shows results only for the stellar velocity maps.