2.1 Circular Model

The simplest model included in XookSuut is the circular rotation model, which is the most frequently adopted for describing the rotation of galaxies. It assumes no other movements than pure circular motions in the plane of the disk and describes the rotation curve of disk galaxies.

Assuming that particles follow circular orbits on the disk, the circular model is given by the projection of the velocity vector V→ along the line-of-sight direction:

V _t is the circular rotation or azimuthal velocity and is a function of the galactocentric distance; Vsys is the systemic velocity and is assumed constant for all points in the galaxy. In this equation and in the following, r is the radius of a circle in the disk plane, which projects into an ellipse in the sky plane. The angle θ is the azimuthal angle relative to the disk major axis, and i is the disk inclination angle.

2.2 Radial Model

When radial motions are not negligible, the disk circular velocity is described by two components of the velocity vector: the tangential velocity (Vt) and the radial one (Vr). In this way, the model including radial velocities is described by the following expression:

Comparing with equation 1, the only difference is the addition of the Vrsinθ term. This term accounts for axisymmetric radial flows (inflow or outflow) on the disk plane.

2.3 Bisymmetric Model

The bisymmetric model describes an oval distortion on the velocity field, such as that produced by stellar bars (e.g., ^{Spekkens & Sellwood 2007}; ^{Sellwood & Spekkens 2015}), or by a triaxial halo potential. In the presence of an oval distortion particles follow elliptical orbits elongated towards an angle that in general differs from that of the disk position angle (e.g., ^{Spekkens & Sellwood 2007}). This kinematic distortion shows a characteristic “S” shape in the projected velocity field that makes the minor and major axes not orthogonal (e.g., ^{Kormendy 1983}). Given that this pattern has been mostly observed in the velocity field of barred galaxies, we will refer to the origin of the oval distortion to stellar bars, although it is not necessarily the case, as mentioned before. The model that intends to describe this pattern is called bisymmetric model (e.g., ^{Spekkens & Sellwood 2007}) since most of the perturbation is kept in the second order of an harmonic decomposition on the disk plane. The bisymmetric model is described by following the expression:

V2,t and V2,r are the nonaxisymmetric velocities induced by the oval distortion and represent, respectively, the tangential and radial deviations from V _t , where the latter describes the disk circular rotation. The angular variable θbar is the location relative to the position angle of the bar ( ϕbar ), in this way³:

Note that in this expression both angles are measured on the disk plane. If ϕbar represents the major (minor)-axis position angle of a bar, then both V2,t(r) and V2,r(r) have positive (negative) values. Unlike the disk position angle ϕdisk' , ϕbar is not a variable than can be easily recognized from the velocity field of barred galaxies; however, its projection on the sky plane is related to ϕdisk' and the disk inclination angle, as follows:

where ϕbar' is the position angle of the bar on the sky plane. Computationally, it is more practical to estimate ϕbar instead of ϕbar'. When the oval distortion is produced by a stellar bar, ϕbar' is expected to be aligned with the photometric position angle of the bar, while the radial profile of V2r(r) and V2t(r) should extend to the length of the bar.

2.4 Harmonic Decomposition

Similar to the photometric decomposition of galaxy images into light profiles via Fourier expansions, the line of sight (LoS) velocity field of a galaxy can be expressed as a sum of harmonic terms as follows:

where c _m and s _m are the harmonic velocities, m is the harmonic number, and θ and r have the same meaning as before. For convenience we have taken the inclination angle out of the Fourier expansion; also note that the 0^th order of the expansion c0(r) is assumed a constant equal to the systemic velocity. However, in addition to the expansion up to M = 1, where we recover the radial model, note that c1≈Vt and s1≈Vr; the expansion to higher orders does not offer a direct interpretation of c _m and s _m since these terms represent a mere decomposition of the LoS velocities, although it is possible to assign to these velocities a physical meaning. The harmonic number is closely related to perturbations of the gravitational potential; under the epicycle theory, such perturbations will induce the appearance of harmonic sectors in the LoS velocities in such a way that if the gravitational potential contains a perturbation of order m, the LoS velocities contain the m + 1 and m - 1 harmonic terms of the Fourier expansion (see ^{Schoenmakers et al. 1997}, for a detailed description). For instance, a bar-potential can be described by a 2^nd order perturbation, which means that the LoS velocity field will contain the 1^st and 3^rd harmonic terms of equation 6 (e.g., ^{Wong et al. 2004}; ^{Fathi et al. 2005}). Similarly, this analysis can be extended for the case of spiral arms (e.g., ^{van de Ven & Fathi 2010}).

In XooSuut the harmonic model, (equation 6), can be expanded to any harmonic order, although, most of the non-circular motions induced by spiral arms or bars are captured by a third order expansion (e.g., ^{Trachternach et al. 2008}).

The harmonic model was first included in the GIPSY task RESWRI I (e.g., ^{Schoenmakers 1999}) RESWRI under the assumption of a thin disk. Afterwards the harmonic decomposition was generalized in KINEMETRY (e.g., ^{Krajnovi´c et al. 2006}) Kinemetry including not only disks but also triaxial structures. The major difference between RESWRI and XooSuut is the assumption of a flat disk. While RESWRI and KINEMETRY allow varying the disk position angle and inclination during the fitting analysis, XookSuut keeps these angles fixed to allow the residual velocities of a circular model to be adjusted with non-circular motions and not absorbed by the variations of these angles, as explained before. However, when large variations of ϕdisk' or i are present throughout the disk, XookSuut will fail in the interpretation of the harmonic velocities, even when the fit is successful.

3.1 χ² Minimization Technique

As we will see in further sections, Bayesian methods like MCMC require to start sampling around the maximum a posteriori, or maximum likelihood, to generate new samples also known as chains; this requires necessarily to find those parameters that minimize the residuals from a given kinematic model and the data. Therefore, in the following we describe the method to solve for each of the different kinematic components of the models, and the constant parameters. The first part corresponds to the analysis adopted in DiskFit (e.g., ^{Spekkens & Sellwood 2007}; ^{Sellwood & Spekkens 2015}) with minimum changes.

A given set of initial conditions for the geometric parameters defines the projected disk with an elliptical shape on the sky plane. Ideally, the initial conditions for the gaseous disk geometry should be close to that of the stellar disk. This geometry will be the starting configuration for the minimization analysis; then, the field of view is divided into K concentric rings of fixed width that follow the same orientation as before. The maximum length of the ellipse semi-major axis can be easily set-up as described in Appendix A; this will create a 2D mask and only those pixels inside this maximum ellipse will be considered for the analysis. The geometry of this mask will be adapted in subsequent iterations until reaching the orientation that better describes the observed velocity field.

The algorithm will solve for each ring a set of velocities that will depend on the kinematic model considered, namely equations 1-3 or equation 6 . Thus, the number of different velocity components to derive will be K velocities in the circular model (Vt,K); 2K in the radial model (Vt,K, Vr,K); 3K in the bisymmetric model (Vt,K, V2r,K, V2t,K) and 2MK in the harmonic model (c1,K,...,cM,K,s1,K,...,sM,K).

The velocity map consists of a two-dimensional image of size nx×ny, with N observed data points Dn with individual errors σn. Let V← be the set of velocities that describe the corresponding kinematic model (namely, V←=〈V←t, V←2,t,V←2,r〉 for the bisymmetric model and similarly for other models). Frequentist methods adopt the chi-square χ2, to derive from a model the set of parameters that describes the data. In this case, the reduced χr2 for the different kinematic models is given by:

Here v is the total number of degrees of freedom (i.e., ν=N-Nvarys, and Nvarys is the number of parameters to estimate from the model ); Wk,n are a set of weights that depend on the pixel position, and will serve to define an interpolated model; and V←k is the set of velocities in the k-th ring that describes the considered kinematic model.

Each kinematic component from V←k would require different weights. For instance, for the circular model the weights adopt the following expression:

where the super-index t makes reference to the circular rotation component. The radial model requires two different weights for the different kinematic components:

Similarly, the bisymmetric model would require three different weights:

Finally, the harmonic decomposition model will have 2M weights given by:

Note that for M = 1 it reduces to the radial model.

The wk,n terms define the interpolation method to be performed between the K-rings. As in DiskFit, these weights adopt the form of a simple linear interpolation given by the usual expression:

where r _k and r _k+1 are the position of the k ^th and (k+1)th rings respectively, and δrk=rk+1-rk is the spacing between rings. As the first ring (k = 1) cannot be placed at the kinematic centre, XookSuut implements different strategies for assigning velocities to pixels down the first ring. Depending on the spatial resolution of the data or the signal-to-noise ratio (S/N), one may opt for one of the following extrapolation options.

The first method is to assume that velocities grow linearly from zero to the velocities derived in the first ring (V←1). This implies that V←0=0 at r = 0; therefore, the kinematic center does not rotate. In the second approach, the set of velocities and positions (V←1,r1) and (V←2,r2) are used to extrapolate velocities to pixels down r ₁; in this way V←0≠0 at r = 0. The third option allows the user to fix the velocity at the origin to some value. Then V←0 and V←1 are linearly interpolated for sampling pixels down r ₁.

As V←k is linear in equations 1-3 and equation 6, we can set the derivative with respect to V←j in equation 12, giving as result:

Rearranging this expression we obtain:

The minimization technique from equation 24 was first introduced by ^{Barnes & Sellwood (2003)}, and subsequently incorporated into DiskFit (^{Sellwood & Spekkens 2015}). Here, equation 24 is generalized for the harmonic decomposition model. The latter expression is a system of linear equations for the V←k unknowns; thus V←k values are solved arithmetically. As mentioned before, the number of velocity components V←k depends on the number of rings and the adopted kinematic model; therefore the dimensions of the matrix to solve will increase as more rings are included in the analysis, and as the kinematic model becomes more complex.

Given a set values for the constant parameters and K rings positioned at r _k on the disk plane, we can solve for V←k in equation 24 by assigning uniform weighting factors (wk,n=1). This way we are obtaining a row-stacked velocities (r _k vs. V _k ). In essence the row-stacked velocities represent the average velocity of each ring; then, we use these velocities as initial conditions to perform an iteratively least squares analysis (LS) through equation 12, but now with the proper weighting factors (namely equations 13-19). Note that if the initial geometric parameters passed to the algorithm are close the true ones, then the arithmetic solutions to V←k should be close the true velocities. This can speedup the MCMC sampling as we will see in further sections.

The minimization procedure in equation 12 is performed by constructing a 2D model from the interpolation weights of equation 21. In each χ2 iteration a new set of velocities V←k' are obtained, together with a new set of constant parameters. The latter will define a new geometry for the mask, and new row-stacked velocities will be obtained for the next iteration. Multiple rounds of χ2 minimization will be performed up to some maximum iteration defined by the user, or until the difference in χ2 evaluations varies less that 10%. Commonly after three iterations the disk geometry becomes stable and simultaneously V←K. Figure 1 summarizes all the fitting procedure in a flowchart.

Figure 1: Flowchart of the fitting procedure to derive the best kinematic model. 

Rings not proper sampled with data may give absurd values of V←k when solving equation 24. To avoid this problem, we define a covering factor to guarantee a minimum number of data per ring. If the covering factor is 1, it means that rings must be 100% occupied by data to estimate V←k, as described in Appendix A. In addition, isolated pixels in the image due to low S/N may not be desired during the analysis XookSuut allows to remove these pixels by excluding those with velocity errors greater than certain threshold defined by the user.

To perform the LS analysis in equation 12, XookSuut adopts the Levenberg-Marquardt (LM) algorithm included in the lmfit package (e.g. ^{Newville et al. 2014}). This algorithm has the advantage that it is fast, although it is widely known to be susceptible to getting trapped at a local minimum. Note that DiskFit adopts the Powell method since this method only performs evaluation of functions with no derivatives performed.

So far the algorithm adopts an LS method for deriving the best parameters defined by the kinematic models. In the following we use sampling methods to infer the posterior distributions of the parameters.

3.2 Bayesian Analysis

The novelty of XookSuut resides in the estimation of the posterior distribution of the non-circular motions and the model parameters. Given the high-dimension of the models, it is desired to perform a thorough analysis of the prior space to obtain the most likely solutions to the problem for each kinematic model regardless of its complexity. For this purpose we adopt Bayesian inference methods for sampling their posterior distributions. Other packages like 2DBAT and KinMSpy (i.e., ^{Oh et al. 2018}; ^{Davis et al. 2013}) also use Bayesian approaches for extracting rotational curves of galaxies. The difference is that KinMSpy is able to fit non-circular motions (radial and bisymmetric).

According to Bayes’ theorem, given a set of data D described by a model function M with parameters α→=α1, α2, .., αn, the posterior distribution of α→ given D follows the expression:

where p(α→|D,M) is the joint posterior distribution of the whole set of parameters; p(D|α→,M) is the probability density of the data given the parameters and the assumed model; p(α→) is the prior probability distribution of the parameters and p(D,M) is a normalization constant also know as marginal evidence or evidence. It is common to find equation 25 expressed in terms of the likelihood function L, as follows:

with the evidence defined as:

where the integral is computed over all the parameter space defined by the priors, Ωα. The evidence can be interpreted as the likelihood of the observed data under the model assumptions; in other words, it is the average of the likelihood over the priors.

The final goal of Bayesian inference is to obtain the posterior distribution p(α→|D,M) of all parameters α→ describing the model function M. Multiple methods have been developed for this purpose. For instance, Markov-Chain Monte Carlo (MCMC) methods evaluate the unnormalized posterior distribution (i.e., p(α→|D,M)∝L(α→)p(α→)), by generating samples or chains from the likelihood function. One of the main characteristics of Markov chains is that the position of a point in the chain depends only on the position of the previous step. Different algorithms with automating chain tunning have been developed to efficiently sample the posterior distribution. Among the most popular MCMC samplers are those who implement affine-invariant ensemble sampling and ensamble slice sampling (e.g., ^{Foreman-Mackey et al. 2013}; ^{Karamanis et al. 2021}).

Other methods, such as nested sampling (NS, ^{Skilling 2006}), are designed to compute the evidence by numerical integration of equation 27, which often makes them computationally more expensive than MCMC methods. This integral is performed from the priors space (or prior volume), and unlike MCMC, does not require an initialization point. Nevertheless, computing the evidence is crucial for model comparison as it represents the degree to which the data are in agreement with the model. Although the main goal of NS is to compute the evidence, the posterior distribution is obtained as a by-product; because of that, NS methods are becoming popular for the inference of parameters in astronomy (see ^{Ashton et al. 2022}, for a thorough description of the method).

One of the advantages of nested sampling with respect to MCMC methods is regarding the convergence criteria. There is no defined convergence criteria among MCMC algorithms, although some of them are based on the number of independent samples in the chains, the so-called integrated autocorrelation time (IAT); however this is often evaluated a posteriori. If the whole chain contains between 10-50 times the IAT, then it is a good indicator that chains are converging (^{Foreman-Mackey et al. 2013}; ^{Karamanis et al. 2021}). In contrast, in nested sampling the stopping evaluation criterion is well defined, since sampling stops after the whole prior space has been integrated.

A detailed discussion of these two sampling methods is, however, beyond the scope of this paper. Following we show the implementation of MCMC and nested sampling methods for the parameter extraction of the kinematic models presented in Sec. 2.

4.1 Toy Model Example

As an example of its use, we run XookSuut on a simulated velocity map. This is the velocity field of a galaxy at 31.4 Mpc with a 32 optical radius. We model a velocity field with an oval distortion described by equation 3. For the rotation curve we adopt the parameterization from ^{Bertola et al. (1991}).. The non-circular motions were modeled using the Gamma probability density function; Gamma(2,3.5) for describing the V2,t component and Gamma(2,3) for V2,r. The constant parameters were set to ϕdisk'=77∘, ϕbar=35∘, i=55∘, x0=76.5, y0=75.2 and Vsys=2142 km s^-1 . The field of view (FoV) is defined as 64×74 and the pixel scale was set to 0.5. Finally we convolved the image for decreasing its spatial resolution. We simulate a circular PSF with a 2D Gaussian function with a 1 full width at half maximum (FWHM). We perturb the velocity profiles by adding Gaussian noise centered at zero and with a standard deviation of 5 km s^-1 .

We started XookSuut by assigning random values to each of the constant parameters, except for ϕbar which is initialized around its maximum (i.e., ϕbar=45∘). We set the initial and last ring exploration at 2.5 and 40 respectively; we also estimate the radial velocities each 2.5. A LS analysis was performed with 3 round iterations before starting the MCMC run. For comparison with Bayesian methods we adopt 1000 bootstraps during the LS analysis. For MCMC sampling, we run a total of 4000 iterations with 60 different chains, which represent twice the number of free variables for this case. We discarded 50% of the joint chain, for a total of 120k posterior samples on each parameter. The IAT for this run resulted in 127 which is superior to 50.

On the other hand, nested sampling only requires the prior information, for which we adopted the uniform priors from Table 1. No initial LS was performed for this case. We stopped the sampling procedure only when the remaining evidence to be integrated was ≤0.1.

We test all the different kinematic models, i.e., circular, radial, bisymmetric models and we arbitrarily expand the harmonic series up to M = 3. MCMC and NS derive the posterior distribution for each variable of the kinematic model; thus, we can take advantage of corner plots to represent their marginalized distributions and explore possible correlations between parameters. The median values and 2σ errors for each parameter are shown in Table 2, while in Figure 2 we show the marginalized posteriors; for simplicity we only show the constant parameters for the bisymmetric model, although note this should be a 30×30 dimensions plot.

Fig. 2 Marginalized distributions of the parameters describing the bisymmetric model for our toy model described in § 4.1. This corner plot shows only the parameters describing the disk geometry. Contours enclose 68% and 95% of the data. Histograms of individual distributions are shown on top, together with the median values and 2σ credible intervals for each parameter. The orange straight lines represent the true values. As observed, all parameters but y0, are recovered within the ±1σ region. The color figure can be viewed online.  

MCMC and nested sampling methods converge to the same solutions found by the LS method; this represents a great success for Bayesian methods to derive kinematic parameters from an input velocity map given the large dimension of the likelihood function. We note that the input parameters are recovered in the bisymmetric model, which is not the case for the circular, radial and harmonic, as expected; this is better appreciated in the corner plot from Figure 2. MCMC and nested sampling recover the input parameters within the 1σ credible interval, except for y0 which lies within 2σ; among the constant parameters, the position angle of the oval distortion shows the largest uncertainty. From Table 2, we notice that the uncertainties derived with Bayesian methods and bootstraps are of the same order.

The different velocities, Vt(r), V2,r(r) and V2,t(r), are also recovered within the 2σ errors as observed from the rightmost panel of Figure 3. For consistence, in Appendix C we include the results using DiskFit; we notice that the XookSuut results are in total agreement with those obtained with DiskFit.

Fig. 3 XookSuut results for the toy model example, for the bisymmetric model case. Figures from left to right: simulated velocity field; best two dimensional interpolated model from MCMC; residual map (input minus output). Overlaid on these maps are iso-velocity contours starting at 0 km s^-1 with steps of ±50 km s^-1 . The fourth column shows the radial variation of the different kinematic components included in the model. Here the input velocities are shown with discontinuous lines while continuous lines represent the velocities derived by XookSuut using MCMC (green), nested sampling (blue) and LM + bootstrap (red). The shadowed regions represent 2σ errors. The color figure can be viewed online. 

Table 2 also shows the root mean square (RMS) for each kinematic model; models including non-circular rotation show a RMS value around 8.5 km s^-1 . This leads to the question of which model is the preferred one for describing a particular velocity field. In a statistical sense, when comparing different models one should choose the one with fewer parameters, since more variables in a model often reduce further the RMS, which does not necessarily provide the best physical interpretation of the data. Statistical tests such as the Bayesian information criterion (BIC) penalizes over the variables from the model; BIC is defined in terms of the likelihood (or the chi-square) as, BIC=-2lnLα^+Nvarysln(N), where α^ represents the parameters that maximize the likelihood function and N is the number of data; thus, the model with the lowest BIC should be preferred. Additionally, the evidence Z computed from NS, is a measure of the agreement of the data with the priors; in this way large (small) Z values are more (less) compatible with the priors.

However, when there is little information about the data, or only the data itself, it is difficult to select a model description of the data based on any information other than statistical tests. Regardless of the statistical method adopted, it is important to have a physical motivation for accepting or rejecting a model; otherwise, erroneous interpretations of the velocity field could arise.

For the toy model example, Table 2 shows that non-circular models have similar BIC values. Even when the harmonic decomposition model seems to perform a good fitting based on the residuals, the physical interpretation of the m = 2 components are meaningless for this example. Thus in a real scenario the radial and bisymmetric model should be compared. The Bayesian evidence for the radial and bisymmetric models results in lnZ=-36136 and - 36135, respectively. In a statistically sense both solutions are equally probable; in other words, the data are insufficient for making an informed judgment. This is not surprising given the simplicity of our velocity field model.

4.2 Simulations

We carried out a set of 1000 simulations with different inclination angles ranging from 30∘<i<70∘, disk position angle 0∘<ϕdisk'<360∘, and to avoid degeneracy, the bar position angle varies from 5∘<ϕbar<85∘; velocity profiles and kinematic center have the same values as in the toy example. We also adopt the same sampling configurations as before. We notice that sometimes XookSuut detects the minor-axis bar position angle instead of the major one; in such cases, ϕbar is found 90° away from the minor axis. This result is also a totally acceptable model since the difference resides only in the sign of the bisymmetric components, V2,r and V2,t, which for these cases both have negative values. XookSuut computes the projected major axis position angle of the bar via equation 5, while the projected minor axis position angle is computed by shifting ϕbar by 90°.

In Figure 4 we show the results of the analysis in corner plots; the derived values are shown relative to the true ones, namely Δα=αrecovered-αtrue; for the radial dependent velocities, we subtract the velocity profile of each component from the derived velocities. These results show that the median values of Δα lie around zero for all parameters describing the bisymmetric model; furthermore, the scatter of the differences is contained within the average value of the 2σ credible interval for each parameter. Results from this analysis demonstrate that MCMC and NS methods are able to recover the true parameters of our simulated velocity maps, even when each model is described by +30 free variables.

Fig. 4 XookSuut results for the bisymmetric model on 1000 synthetic velocity maps with oval distortions. Values are reported with respect their true value, namely Δα=αrecovered-αtrue, where α is any of the parameters considered in this example. Straight orange lines represent the true values (i.e., Δα=0). Blue and lime colors show the results for NS and MCMC methods, respectively. The inner and outer density contours contain 68% and 95% of the data, respectively. The upper histograms represent the 1D distributions of the parameters on the x-axis, while the size of the error bars represents the average value of the 2σ credible interval for each parameter. Note that all parameters are recovered within the reported error bars. Values on-top the histograms represent the 50% percentile of Δα, together with the ±2σ dispersion. The color figure can be viewed online. 

The final accuracy of the recovered parameters would depend on the details of data themselves (resolution, S/N , spatial coverage etc.). Therefore, ad-hoc simulations are encouraged.

4.3 NGC 7321

We proceed to evaluate XookSuut over the velocity field of a galaxy hosting a stellar bar. For this purpose we adopt the galaxy NGC 7321 observed as part of the CALIFA galaxy survey (e.g., ^{Sánchez et al. 2012}). This object has been previously analyzed by ^{Holmes et al. (2015)} using DiskFit. The Hα velocity field of this object represents a good example of a galaxy with a strong kinematic distortion, most probably caused by the stellar bar. ^{Holmes et al. (2015)} found the bisymmetric model as the best model for reproducing the inner distortion observed in this object; they found best fit values and 1σ errors for the constant parameters given by ϕdisk' = 12±1∘, i=46±2∘, Vsys = 7123±3 km s^-1 and a bar position angle oriented at ϕbar =47±6∘. We implemented XookSuut on the Hα velocity map taken from the CALIFA data products (e.g., ^{Sánchez et al. 2016}). We adopted the same ring configurations as before, excluding pixels from the error map with values larger than 25 km s^-1 ; we proceed to explore the non-circular motions up to r = 18, and set the maximum radius for the circular velocities up to 40; this leads to a total of 36 free variables that will be estimated with Bayesian inference. We adopt 3 rounds of iterations for the LS method, and also compute the errors on the parameters with 1000 bootstraps. For MCMC, we adopt 5000 steps and drop half of the total samples to let the joint chain stabilize. For NS we stop the sampling when the remaining evidence to be integrated is 0.1.

Figure 5 shows the marginalized distributions of the constant parameters. From the 1D histograms we obtain ϕdisk' = 11±0∘, i=46±1∘, Vsys = 7123±1 km s^-1 and ϕbar =46±6∘; additionally we compute the intrinsic scatter of the data as ≈16 km s -1 . Table 3 shows a summary of these results. As can be read from this table, the constant parameters derived by XookSuut are in concordance with those previously reported by ^{Holmeset al. (2015)}, although our uncertainties are smaller when comparing the errors at 2σ, probably due to differences in methods. The bottom figure shows the best model and residual map obtained from nested sampling. The kinematic distortion observed in the central region is well reproduced with the bisymmetric model. The bisymmetric motions, i.e., the bar-like flows, are oriented at 46±8∘ on the sky plane. The rightmost panel shows the radial profile of the different velocity components derived from NS, MCMC and LS+bootstrap methods. Again, the uncertainties reported from NS and MCMC are of similar magnitude, and these are larger than those obtained with bootstrap methods.

Fig. 5 XookSuut results for NGC 7321 for the bisymmetric model. Top figure shows the marginalized distribution for the constant parameters with MCMC methods in green colors and NS in blue. As in Figure 2, the quoted values on top the histograms represent the median and 2σ errors obtained from the posterior distributions. Bottom figure shows from left to right the Hα velocity map; the two-dimensional model from NS; the residual map; and the radial profile of the different velocities. Shadow regions represent the 2σ credible intervals obtained from each method (bootstrap in red colors). The color figure can be viewed online. 

In Appendix D we show the implementation of XookSuut on other data with different instrumental configurations.