Kinematic Decomposition of the H i Gaseous Component in the Large Magellanic Cloud

We perform a profile analysis of the combined H i data cube of the Large Magellanic Cloud (LMC) from observations with the Australia Telescope Compact Array and the Parkes radio telescope. For the profile analysis, we use a newly developed algorithm that decomposes individual line profiles into an optimal number of Gaussian components based on a Bayesian nested sampling. The decomposed Gaussian components are then classified into kinematically cold, warm, and hot gas components based on their velocity dispersion. The estimated masses of the kinematically cold, warm, and hot gas components are ∼12.2%, ∼58.3%, and ∼29.5% of the total H i mass of the LMC, respectively. Our analysis reveals the highly complex H i structure and kinematics of the LMC that are seen in previous studies but in a more quantitative manner. We also extract the undisturbed H i gas bulk motions and derive new H i gas bulk rotation curves of the LMC by applying a 2D tilted-ring analysis. In contrast to previously derived H i rotation curves, the newly derived bulk rotation curves are much more consistent with the carbon star kinematics, with rotation velocity linearly increasing in the inner part and reaching a maximum of ∼60 km s−1 at the outermost measured radius. By comparing the lower bulk rotation curves with previous studies, we conclude that there is a lower dynamical contribution of dark matter in the central part of the LMC.


Introduction
The Large Magellanic Cloud (LMC), one of the closest galaxies to the Milky Way (∼50 kpc; Pietrzyński et al. 2013), is an ideal laboratory for allowing us to examine the state of its interstellar medium (ISM) and the combined effects of energy deposition in great detail. The enormous gaseous features are evidence of tidal and dynamical interactions with the Small Magellanic Cloud (SMC) and the Milky Way  and references therein). Due to its close proximity, we can identify gas clouds that are moving at anomalous velocities in the galaxy as part of the gaseous feature of the Magellanic System at high angular resolutions. Detailed study of the interplay between star formation and the ISM in the LMC has benefited greatly from its proximity to the Milky Way. There have been extensive multiwavelength observational studies of the LMC and theoretical modeling to investigate the detailed properties of the galaxy's ISM and its link to star formation processes. These investigations include, but are not limited to, the following: (1) star formation processes and their impact on the structure and physical properties of the ambient ISM (Meixner et al. 2006(Meixner et al. , 2013Book et al. 2009;Olsen et al. 2011;Nidever et al. 2017), (2) the search for invisible massive compact halo objects (MACHOs) through monitoring of gravitational microlensing events (Alcock et al. 2000), (3) the study of stellar populations and star formation history (Harris & Zaritsky 2001;Cioni et al. 2014), and (4) the origin and structure formation of the Magellanic System (Putman et al. 2003;Mastropietro et al. 2004;Diaz & Bekki 2012). We refer to van der Marel (2006) and D' Onghia & Fox (2016) for a review of the structure and kinematics of the LMC and Magellanic Stream, respectively.
Observational studies of the LMC at high angular and spectral resolutions show highly complex structure and ISM kinematics. For example: shell-like structures and holes are seen to be associated with the star-forming regions (Kim et al. 1999); H I clouds/clumps extend much further than the H I scale height that could be formed during thermally unstable cooling ; there is evidence of tidally induced H I flows near the 30 Doradus region (Fukui et al. 2017); possible stripped gas forms the Leading Arm IV ; and H I arms extend toward the SMC that form part of the Magellanic Bridge . To study the dynamics of the LMC and its interaction with the Milky Way, studies of proper motion have also proved to be valuable. They show that the Gaia DR2 measurements of proper motion for the LMC are more affected by the assumed dynamical centers of stars (e.g., the center of the H I gas disk or photometric centers) than by systematic uncertainties (Gaia Collaboration et al. 2018). The dynamical center of the galaxy derived using the H I kinematics shows an offset with respect to those of the proper motions from Hubble Space Telescope (HST) and Gaia observations (Kallivayalil et al. 2013;van der Marel & Sahlmann 2016; but see van der Marel & Kallivayalil 2014).
In conjunction with the H I kinematics and results from proper motion (Piatek et al. 2008;Kallivayalil et al. 2013), the study of Indu & Subramaniam (2015) examines the highvelocity gas in the western disk of the LMC and southwest and southern parts of the Magellanic Bridge. Another approach includes the use of kinematics of carbon stars, from which general expressions for the LMC velocity field are derived (van der Marel et al. 2002;van der Marel & Kallivayalil 2014). This method mitigates the large uncertainty of transverse velocity when using the LMC line-of-sight kinematics and models of the Magellanic Stream.
In a qualitative sense, previous studies have shown that modeling of three-dimensional spectral line data cubes is not straightforward when classical moment analysis is applied to velocity profiles that are asymmetric, non-Gaussian, and/or have multiple components. Parameters such as amplitude (A), central velocity (V 0 ), and velocity dispersion (σ) become sensitive to the estimation methodology.  performs a 2D tiltedring analysis of an H I MOMENT1 velocity field of the LMC using intensity-weighted mean (IWM) velocities of individual line profiles and compiling them into a 2D map. As noted, such analysis can yield biased measurements for highly asymmetric non-Gaussian profiles, as is the case for the LMC (Gentile et al. 2004;Oh et al. 2011). This is also true for other model-based parametric profile analysis, for which a single Gaussian or Hermite H 3 polynomial fit is used to determine central velocities of line profiles.
In addition to the methodological limitation of profile analysis, the H I gas in the LMC has been found to significantly deviate from the stellar kinematics in several studies due to the tidal interaction with the SMC and the Milky Way (Weinberg 2000;van der Marel & Cioni 2001;Bekki & Chiba 2007;Mastropietro 2010;Besla et al. 2012). One of the unambiguous deviations is the kinematic center derived from the carbon star kinematics, which is offset by more than ∼1°from the kinematic center of the H I kinematics. Otherwise, the kinematic center of the carbon stars is consistent with the centers of both the stellar bar and the outer isophotes of the LMC (van der Marel et al. 2002). As such, H I data alone would be insufficient to investigate the global ISM structure and kinematics of the galaxy owing to its highly disturbed kinematic nature. Furthermore, this can cause large uncertainties in deriving the underlying global kinematics, and subsequently in mass modeling (e.g., disk-halo decomposition, dark matter distribution) of the LMC.
To circumvent the limited capability of previous profile analysis on modeling the highly complex ISM of the LMC and to improve the analysis of the H I kinematics, we perform a multiple Gaussian decomposition on the combined Australia Telescope Compact Array (ATCA) and single-dish Parkes telescope H I data cube. A newly developed algorithm, BAYGAUD , is used. This algorithm allows us to decompose individual line profiles into an optimal number of Gaussian components based on a Bayesian nested sampling algorithm, practically using the MULTINEST (Feroz & Hobson 2008;Feroz et al. 2009bFeroz et al. , 2019 library, which implements a nested sampling. In this paper, we classify the decomposed Gaussian components of the H I spectral line data cube into kinematically cold, warm, and hot gas components according to their velocity dispersion, and derive their physical properties. Of the decomposed Gaussian components, we also extract H I gas bulk motions of the LMC, which are consistent with a reference model velocity field constructed using the carbon star kinematics. Lastly, we derive an H I gas bulk rotation curve of the LMC, which minimizes the effects of noncircular motions in the gas disk. In the course of deriving the rotation curves, following the general expressions for the LMC line-of-sight velocity distribution described in van der Marel et al. (2002), we correct for the transverse, nutation, and precession motions of the galaxy caused by its proximity and the gravitational effect of the Milky Way.
The structure of this paper is as follows: Section 2 discusses the complex ISM structure and kinematics of the LMC from our Bayesian profile decomposition applied to the combined H I data cube of the galaxy. Section 3 discusses a model reference velocity field constructed using the carbon star kinematics, which represents the global kinematics of the LMC. Sections 4 and 5 discuss the rotation curve analysis of the extracted H I gas bulk motions and the shape of bulk rotation curves compared with the previous ones derived using the H I intensity-weighted velocity field and carbon star kinematics. Section 6 summarizes the main results of this paper.

The Combined ATCA and Parkes H I Data Cube
We use the combined H I data cube from ATCA and Parkes single-dish observations presented in Kim et al. (2003) to perform the kinematic analysis of H I gaseous components in the LMC. The data cube is sensitive to scales from parsecs to kiloparsecs. The data cube has an angular resolution of 60″ and a spectral resolution of ∼1.6 km s −1 . It has a dimension of 1998 × 2230 × 120 (20″ per pixel) and an rms flux of ∼15 mJy beam −1 , which corresponds to ∼2.4 K. We refer the reader to Staveley-Smith et al. (2003) and Kim et al. (2003) for full descriptions of the calibration and imaging of the data cube as well as the method used for combining the interferometric and single-dish observations.
To make all the profiles of the cube independently while doing profile decomposition, we use a Python task of the scipy. ndimage package to regrid the cube to a pixel size of one beam resolution (60″). The pixel size of 60″ corresponds to a physical size of ∼14.6 pc assuming the LMC distance to be 50.1 kpc. This is still smaller than the typical size of giant molecular clouds in the Milky Way (∼50 pc; Blitz et al. 2007). The resulting cube has a dimension of 666 × 743 × 120.
In the left panel of Figure 1, we show a multiwavelength three-color composite image of the stellar and gaseous components of the LMC. It is constructed using the ATCA +Parkes H I integrated intensity map (H I gas, blue; Kim et al. 2003), Spitzer 3.6 μm (old stars, red; Meixner et al. 2006), and Hα (star formation over a period of tens of millions of years, yellow; Gaustad et al. 2001). The composite image reveals the localized gas disturbance presumably influenced by hydrodynamical processes in the LMC. Such disturbed gas kinematics in the regions are represented by highly asymmetric non-Gaussian H I gas profiles as shown in the right panel of Figure 1.
To quantify the kinematic properties of the H I gas component of the LMC, we extract the moment maps from the H I data cube of the galaxy using standard moment analysis (see Figure 2). When deriving the MOMENT2, we use two different clipping options where the MOMENT2 values are computed using the channels whose flux is higher than 1σ or 3σ where σ is the mean standard deviation of the line profiles. This is for checking the effect of the clipping option on the extracted MOMENT2 maps given the noise characteristics of the cube. For the MOMENT0 and MOMENT1 maps, we choose the 1σ and 3σ clipping options, respectively, after testing with several clipping options.
As discussed earlier, conventional moment analysis is limited in modeling non-Gaussian line profiles. The estimated intensity, central velocity, and velocity dispersion of a line profile are often biased by the data noise and complex velocity structure. This latter is imprinted in the distorted velocity pattern of the MOMENT1 map seen in Figure 2. In the following sections, we apply a new profile decomposition technique to the H I data cube of the LMC, in order to improve the kinematic analysis of the galaxy's complex ISM.

The Optimal Decomposition of H I Gas Profiles
As seen in the right panel of Figure 1, it is common to observe multiple velocity components along different lines of sight and asymmetric non-Gaussian H I velocity profiles (see also Kim et al. 2003;Staveley-Smith et al. 2003). Analysis of these complex profiles in an appropriate and quantitative manner is an essential prerequisite for addressing the structure and kinematics of the ISM and the link to other hydrodynamical processes in the galaxy.  presented a new profile decomposition method, BAYGAUD, 5 based on a Bayesian nested sampling statistics to model a line-of-sight velocity profile into an optimal number of Gaussian components. The Bayesian approach circumvents the common issue that is found in the χ 2 minimization technique, in which the algorithm gets trapped in local minima because of poorly defined initial estimates. It also allows us to classify the decomposed Gaussian components into two categories: (1) kinematically cold, warm, and hot components against their velocity dispersions, and (2) bulk and noncircular random motions with regard to a user-defined underlying global rotation of the disk of a galaxy.
We run BAYGAUD on the combined H I data cube and obtain a set of best-fitting parameters of individual line profiles. When carrying out the analysis, we let BAYGAUD fit a given velocity profile with up to five Gaussian components simultaneously. From visual inspection of line profiles of the data cube, this is likely to be good enough to model the non-Gaussianity in the profiles. As described in , BAYGAUD makes five Bayesian fits to each line profile with m = 1 to 5 using the following equation: where G(x) represents the line-of-sight velocity profile, m is the number of Gaussian components. a i , σ i , and μ i are the Gaussian parameters for area (integrated intensity), standard deviation (velocity dispersion), and central velocity, respectively. b j are coefficients of the nth-order polynomial for the baseline fit. Among the five different fits, the most appropriate model is chosen based on the Bayesian evidence. The efficiency of calculating the Bayesian evidence is one of the computational challenges in Bayesian inference because it requires a large amount of computational time. An efficient method to calculate the Bayesian evidence at high accuracy through an efficient sampling of the parameter space has been developed. It is known as the MULTINEST algorithm, which makes use of importance nested sampling (INS, Liu 2008;Feroz et al. 2019). The INS algorithm and its performance test for highly multimodal problems are fully described in Feroz & Hobson (2008; see also Feroz et al. 2009bFeroz et al. , 2011Feroz et al. , 2019. The algorithm has been successfully implemented in solving several astrophysical problems with multimodal and/or degenerate distributions (e.g., Shaw et al. 2007;Feroz et al. 2009a).
For comparing the Bayesian evidences of two models selected out of the five model fits, we use the ratio (R) of posterior probabilities, where Pr(M 0 ) and Pr(M 1 ) are the prior probabilities of the models, D is the data, Pr(D|M 0 ) and Pr(D|M 1 ) are their Bayesian evidences ( 0  , 1  ) or marginal likelihoods. This ratio can be further reduced to the so-called Bayes factor if the prior probability ratio of the two models, Pr(M 1 )/Pr(M 0 ), is set to unity. We use uniform priors of the individual model parameters, which are set to be the same for the five model fits in the BAYGAUD analysis.
To select the best model, we adopt Bayes factor >10, the socalled "strong" model selection criterion as suggested by Jeffreys (1961) and Raftery (1995). This implies that one model has at least 10 times more evidence than the other in explaining the observed data. From this, among the five different models fitted to each line-of-sight H I velocity profile of the LMC, we find the most appropriate model that has a Bayesian evidence more than 10 times that of the second best model.

Quantification of the Non-Gaussianity of H I Gas Profiles
In Figure 3, we present all the optimally decomposed Gaussian components of each velocity profile of the LMC H I data cube. The integrated signal-to-noise ratio (S/N) of all the decomposed Gaussian components is larger than 4. We derive the integrated flux (F int ) in units of K km s −1 of a decomposed Gaussian component from its Gaussian fit as follows: is a Gaussian function with gas velocity v, central velocity v 0 , and velocity dispersion σ. The velocity range over which the line profile is integrated is quantified by n in units of σ. We set n = 3 (i.e., v 0 − 3σ v v 0 + 3σ), over which ∼99.7% of the total flux of the profile is covered. The integrated noise (σ int ) in units of K km s −1 is estimated as where σ chan is the rms noise of each channel, ΔV is the channel width of the cube, and N is the number of channels corresponding to 6σ velocity width (i.e., 6σ/ΔV ).
In addition, a map of the optimal number of Gaussian components extracted is shown in the N_Gauss panel of Figure 4. Non-blank pixel values of the map range from 1 to 5. A significant number of H I velocity profiles are modeled as a set of multiple Gaussian components.
This demonstrates the complex structure and kinematics of the H I component in the LMC in a quantitative manner. The kinematic disturbances could be induced by tidal interactions with the SMC and the Milky Way as well as star formation activity in the LMC as discussed in previous studies (see For et al. 2014 and references therein). The Hα image in panel (d1) of Figure 4 shows several distinct star-forming regions, as Figure 2. H I moment maps extracted from the combined H I data cube of the LMC described in Section 2.1: (a) MOMENT0 (integrated intensity in units of K km s −1 , clipping the channels whose fluxes are below 1σ), (b) MOMENT1 (centroid velocity in units of km s −1 , clipping the channels whose fluxes are below 3σ; velocity contours are spaced by 20 km s −1 ), (c1) MOMENT2 (velocity dispersion in units of km s −1 , clipping the channels whose fluxes are below 1σ), and (c2) MOMENT2 (same as in (c1) but 3σ clipped). See Section 2.1 for details. indicated by circles, where star clusters are located (Carlson et al. 2012;Anderson et al. 2014). This star-forming activity interacts with the ISM in a way that clearly disturbs the ambient ISM, resulting in such non-Gaussian H I gas profiles. Together, tidally disturbed gas motions that are associated with the SMC and the Milky Way could contribute to part of the gas disturbance. In Figure 4, active star-forming regions including the Tarantula nebula (30 Doradus), identified with Hα emission, and H I gas streaming motions in the outer region of the galaxy associated with the Magellanic Stream and the Magellanic Bridge are identified by the large number (e.g., >3) of optimal Gaussian components, which even reaches 5.

Kinematic Classification of H I Gas Profiles
In the left panel of Figure 5, we show the distribution of all the decomposed Gaussian components from the H I data cube with respect to their velocity dispersions. The profiles with a peak intensity S/N less than 2 are omitted. We fit the skewed histogram of the H I velocity dispersions with multiple Gaussians and find that a fit using five Gaussian components appears to be plausible in explaining the distribution. A fit with a smaller number of Gaussian components (2, 3, or 4) could not explain the distribution well, particularly failing to fit the tail region with σ > 40 km s −1 .
We classify the decomposed Gaussian components of the H I gas in the LMC into three main kinematic groups with respect to their velocity dispersion, i.e., "kinematically" (1) cold (ΔV < σ < 4.8 km s −1 ), (2) warm (4.8 km s −1 < σ < 18.1 km s −1 ), and (3) hot (σ > 18.1 km s −1 ), where ΔV is the channel resolution of the cube (1.6 km s −1 ) as indicated by the vertical dotted lines in the left panel of Figure 5. For comparison, previous studies show that the mean velocity dispersion values of narrow and broad H I gas components of galaxies in the local universe are 5 and 8 km s −1 , respectively (Young & Lo 1996, 1997de Blok & Walter 2006;Ianjamasimanana et al. 2015Ianjamasimanana et al. , 2017Mogotsi et al. 2016).
In the right panel of Figure 5, we show the H I mass histogram of all the decomposed Gaussian components of the cube with respect to their velocity dispersions. We also superimpose the ones derived from single-component Gaussian fits and the MOMENT2 velocity dispersions. As described in Section 2.1, the MOMENT2 is derived applying two different clipping options, 1σ and 3σ clippings, to the cube, while the MOMENT0 that is used for measuring the corresponding H I mass of the velocity dispersion histogram in Figure 5 is derived using the 1σ clipping option.
In Figure 5, the histograms denoted as "Moment (1σ)" and "Moment (3σ)" represent the results derived using the MOMENT2 maps extracted applying the 1σ and 3σ clipping options, respectively. The MOMENT2 velocity dispersion values measured using the 1σ clipping option are predominantly higher than the others with a mode value of ∼21 km s −1 in the mass-velocity dispersion histogram. This is mainly attributed to the non-Gaussian or multicomponent profiles as discussed in the previous sections (see Figures 1 and 3). For the same reason, compared to the single Gaussian fits, the optimal Gaussian fit was able to retrieve more H I gas with narrow velocity dispersions less than ∼5 km s −1 , which corresponds to ∼1 × 10 7 M e .
The intensities, velocity fields, and velocity dispersions of the kinematically cold, warm, and hot H I components with respect to their velocity dispersions are mapped in Figures 4(a) -(c). By comparing these velocity fields (panels (a2), (b2), and (c2)), they do not appear to be much different except for the outer and some localized regions such as the 30 Doradus region. The localized regions can be largely associated with star-forming regions as shown in the Hα image in Figure 4(d1). The star formation being traced by Hα emission represents star formation over the last few tens of millions of years (Harris & Zaritsky 2009).
As shown in panels (b1) and (c1) of Figure 4, high-intensity regions of the kinematically warm and hot components tend to be located around the active star-forming regions, particularly  (Gaustad et al. 2001). N_Gauss: the optimal number of Gaussian components derived for each spaxel whose integrated S/N is larger than 4. See Section 2.4 for more details.
the 30 Doradus region. In addition, a spiral-like structure in the kinematically warm and hot components is found in the outer region, while the cold component is mainly distributed in the disk. The spiral-like streaming motions of the gas can be associated with the Magellanic Stream and the Magellanic Bridge. The radiative and mechanical stellar feedback could impact on the structure and kinematics of the ambient interstellar medium, giving rise to turbulent gas motions with larger velocity dispersions as found in the kinematically warm and hot gas components (e.g., Lee et al. 2019).
On the other hand, the kinematically cold component appears to have a more regularly rotating velocity field and is distributed more uniformly in intensity than the warm and hot components. Some of the kinematically cold gas component can be either leftover ISM from the recent star formation or a cold gas reservoir for future star formation.
From the Bayesian profile decomposition, we have been able to better deblend the multiple gas components with an optimal number of Gaussian components, deriving an H I mass of 4.73 × 10 8 M e , which is comparable to that measured from the 1σ clipped MOMENT0 map. The fractions of the kinematically cold, warm, and hot components are 12.2%, 58.3%, and 29.5% of the total H I gas in the galaxy, corresponding to 5.77 × 10 7 M e , 2.76 × 10 8 M e , and 1.40 × 10 8 M e , respectively. Notably, compared to both the MOMENT0 and single-component Gaussian fits with integrated S/N larger than 4, we extract a higher mass of ∼1 × 10 7 M e in the kinematically cold components with σ < 5 km s −1 . Some of the kinematically cold (σ < 5 km s −1 ) component will be directly associated with star formation because H I gas should pass through this phase before turning into H 2 (Ianjamasimanana et al. 2012). The fractions of the decomposed Gaussian components and their estimated H I masses are summarized in Table 1.

The Global Kinematics of the LMC
Gas components in the LMC are found to be significantly disturbed in both morphology and kinematics (van der Marel et al. 2002). This is mainly caused by two effects: (1) the tidal and hydrodynamical interactions with the SMC and the Milky Way, and (2) the internal stellar feedback in the LMC. To first order, they result in a clear offset (1°) for the kinematic center between the H I and stellar kinematics. The kinematic center of the LMC derived from the carbon star kinematics is consistent with the morphological center of the stellar bar, which is determined from the outer isophote fits to the stellar component. On the other hand, the H I kinematics may be more affected by recent violent events (van der Marel et al. 2002).
It has also been found that the dynamical centers for young stellar populations are significantly offset from those of older populations (Wan et al. 2020). H I gas and young stars are likely to be more affected by hydrodynamic and gravitational forces from recent baryonic processes in the LMC and interactions with the Milky Way and the SMC than old stars, which formed before such events. Examining the kinematically decomposed gas components as shown in Figure 3, we find that part of the gas bulk motions remains alongside the disturbed  gas motions. We could potentially extract this component more accurately by selecting those components moving at close to the velocities expected from the global kinematics to represent the global kinematics of the LMC. To this end, we need a reference velocity field.

The Model Line-of-sight Velocity Distribution: Correcting for Transverse, Precession, and Nutation Motions
For the reference global kinematics of the LMC, we construct a model using the carbon star kinematics presented in van der Marel et al. (2002). The stellar velocity field is not directly subject to modification by stellar feedback and hydrodynamic processes, and may be better suited for tracing the undisturbed global kinematics of the LMC. However, owing to the large angular extent (>10°) of the LMC on the sky, the "flat sky" assumption is not valid. A correction for the transverse velocity is therefore required. In addition, the nutation and precession of the LMC that are induced by the tidal torques from the Milky Way affect its line-of-sight velocities.
According to van der Marel et al. (2002), the model for lineof-sight velocity is given as  where W ts and V tc are given as follows: We refer the readers to van der Marel et al. (2002) for a complete description of the model for line-of-sight velocity and its application to the kinematic analysis of carbon stars in the LMC (see also Di Teodoro et al. 2019 for the kinematic analysis of the SMC). Using Equation (9), we can correct the observed line-of-sight velocity field, v LOS , of the LMC disk for the precession and nutation as well as the transverse motions, which yields the internal kinematics of the LMC disk.

A Global Model Velocity Field Based on the Kinematics of Carbon Stars
In Figure 6, we show the rotation velocity derived from the carbon star kinematics that is corrected for the nutation, precession, and transverse motions of the LMC disk, and the associated kinematic parameters including the time derivative of the inclination angle estimated in van der Marel et al. (2002). This represents the global kinematics of the LMC because the carbon star kinematics is expected to be dynamically less affected by any hydrodynamical effect of the LMC than the gaseous component. To suppress any unphysical small-scale fluctuations in the rotation velocity of the carbon stars, we smooth the curve using a model circular velocity caused by a pseudo-isothermal halo (e.g., Begeman et al. 1991) given by where ρ 0 and R C are the core density and core radius of a halo, and G is the gravitational constant. The fit result is shown by the dashed line in the left panel of Figure 6. We note that the reference velocity field being used for profile decomposition does not need to be an accurate prediction of the global kinematics. It is sufficient if it represents the velocities of the H I components left behind close to the extracted global kinematics.
Using the smoothed rotation curve and the kinematic parameters of carbon stars in van der Marel et al. (2002), we construct a velocity field model after taking the effect of nutation, precession, and transverse motions of the LMC into account as given in Equation (9) for the model velocity. The resulting reference model velocity field that is used for classifying the decomposed Gaussian components of the LMC H I data cube is presented in the right panel of Figure 6.

Extraction of H I Bulk Gas Motions
Given the reference velocity field constructed in Section 3.1, we extract corotating H I bulk gas components whose kinematics is closest to the global kinematics predicted from the LMC carbon star kinematics within a velocity limit. The velocity limit is chosen by considering the several effects, including the gravitational potential of the LMC, stellar feedback, bar potential, geometry, and tidal interactions with the SMC and the Milky Way etc. In this work, we use the H I velocity dispersion map derived from the single Gaussian fit (SGFIT) as a velocity limit to the LMC H I data cube.
To check the adopted sensitivity of the velocity limit when extracting the H I gas bulk motions, we obtain a set of velocity limit maps by multiplying the SGFIT velocity dispersion map by three factors (1.0, 2.0, and 3.0). While the velocity limit with 3.0 times the SGFIT velocity dispersion map appears to be wide enough to cover the velocities of any leftover corotating H I gas disk associated with the global kinematics of the LMC, we restrict the velocity limit to 1.0 times the SGFIT velocity dispersion map. This is a conservative choice but provides sufficient filling factor for the resulting velocity field to perform the 2D tilted-ring analysis in Section 4.2.
Of the decomposed Gaussian components in Figure 3, we select the ones whose central velocities are the closest to the reference velocity field in the right panel of Figure 6 within the velocity limit defined using the SGFIT velocity dispersion map as H I gas bulk motions. If no significant (i.e., integrated S/ N < 5) H I gas is found at the predicted global velocities within the velocity limit, we put a blank value into the corresponding pixel of the H I gas bulk velocity field. We then make 2D maps of integrated intensity, central velocity, and dispersion of the resulting bulk motions. The derived H I mass of the bulk motion components is 1.21 × 10 8 M e , which is about 26% of the total H I mass of 4.73 × 10 8 M e as given in Table. 1. A selection bias introduced by the reference velocity field can arise in the line profiles that include multiple kinematic components within the range of velocity limit. In this case, the reference velocity field used may affect the extracted bulk H I gas by forcing it be linked to the other kinematic components.
To check this possibility, we derive the fraction of pixels in the extracted bulk H I velocity field that have more than one kinematic component within the velocity range out of all the pixels available. It is ∼27%, which is relatively small. The effect of such a bias on the bulk H I velocity field is likely to be insignificant.
For illustration purposes, we only show the maps of the bulk motion extracted using the velocity limits of the SGFIT velocity dispersion map in Figure 7. A distinct difference is shown in the centroid velocities of the bulk motions compared to the ones determined from both the moment analysis, MOMENT1 in Figure 2, and the SGFIT result in Figure 7.
Primarily, the LOS velocity difference is most significant near the northeastern region of the stellar bar and the outer regions of the galaxy, corresponding to the regions where signs of significant gas disturbance are found, as indicated by the large velocity dispersion of ∼20 km s −1 in the SGFIT velocity dispersion map.
The velocity difference is mainly attributed to the multiple kinematic components projected along the LOS. In practice, this biases the measurement of the centroid velocity of a line profile estimated in an intensity-weighted manner. In contrast, the centroid velocity map of the bulk motions shows better consistency with the one expected from the carbon star kinematics. The gas bulk motions less affected by hydrodynamical processes in the galaxy can be extracted from the optimally decomposed Gaussian components based on the reference global kinematics from the kinematics of carbon stars. The smaller filling factor of the bulk motion map could be due to the relatively narrow velocity limit defined by the SGFIT velocity dispersion map when extracting the bulk motions. As mentioned above, we note that all the extracted bulk velocity fields have integrated S/N values larger than 5.

Rotation Curve Analysis
We derive bulk gas rotation curves of the LMC using the bulk velocity field extracted in Section 4.1. To this end, we fit the model LOS velocity in Equation (9) to the bulk velocity fields. In addition, we need to correct for the transverse, nutation, and precession motions. For this, we use the 2D Bayesian automated tilted-ring fitter (2DBAT), which fits 2D tilted-ring models to a velocity field using a Bayesian nested sampling tool (Oh et al. 2018). For the full description of 2DBAT including its performance test and practical application to a sample of galaxies from H I interferometric observations, we refer to Oh et al. (2018).
In this paper, we modify the likelihood function based on a Student-t distribution in 2DBAT to account for the additional parameters of the model LOS velocity in Equation (9), which are for the transverse, precession, and nutation motions as well where v LOS (ρ, Φ) is the observed line-of-sight velocity at (ρ, Φ), ν is the normality parameter set to 3 in this analysis, and Γ is the gamma function given as The distribution is scaled by a free parameter, ò in Equation (13), which is largely proportional to the standard deviation of the residual map between the observed and model velocity fields. We regularize the radial variations of rotation velocity v ROT (R), position angle f(R), and inclination i(R) with B-spline functions to suppress any sudden change or unphysical fluctuations when deriving their most representative global shapes. From fitting with different sets of uniform priors for the kinematic center (α 0 , δ 0 ), we find that it could not be reliably determined from the H I 2D tilted-ring analysis alone as there is a large difference among them. In all cases, the best-fitting parameters of the kinematic center converge at the edge of the given priors, which deviate significantly from the one derived from the carbon star kinematics. This is mainly due to the shallow velocity gradient in the bulk velocity field along the kinematic major axis of the galaxy, as indicated by almost parallel isovelocity contours (see Figure 7). This results in a degeneracy between the solution for the kinematic center and the other tilted-ring parameters. Together, the fitting can also be affected by the distorted velocity field in the outer part, which is presumably caused by the tidal interaction with the SMC and the Milky Way. Having more outer pixels would contribute significantly to the likelihood estimation in Equation (13). We therefore fit the parameters of the likelihood function in Equation (13) after fixing the kinematic center (α 0 , δ 0 ) derived from the carbon star kinematics (van der Marel et al. 2002).
We fit the 2D tilted-ring models that incorporate the kinematic effects of the nutation, precession, and transverse motions of the LMC to the bulk velocity field extracted using the velocity limit map of 1.0 times the SGFIT velocity dispersion map. In the course of the nonparametric 2D tiltedring analysis for the bulk velocity field, we set various orders of B-spline curves for v ROT (R), f(R), and i(R) in Equation (9) as summarized in Table 2. This is for checking the sensitivity of the B-spline orders adopted to the regularization of the parameters. Given the regularization setup of f, i, and v ROT in Table 2, we derive eight rotation curves in total. All the fits were made without any failure. We find that there is a clear radial variation of the position angle f of the galaxy, which is well described by the cubic spline function. No significant radial variation is found for the inclination i in the cubicsplined fit. For the rest of this paper, we discuss the fitting

Asymmetric Drift Correction
We correct for the dynamical effect of pressure support caused by random gas motions in the LMC on its rotation curves. This so-called asymmetric drift correction is usually required for dwarf galaxies whose maximum rotation velocity is comparable to gas velocity dispersions in the disk. As shown in Figure 8, the H I gas velocity dispersion in the outer region of the LMC is ∼10 km s −1 , which is much lower than the rotation velocity of ∼40 km s −1 . Following the method described in Bureau & Carignan (2002) (see also Oh et al. 2011), we correct for the asymmetric drift as follows: where v ROT is the rotation velocity derived from the 2D tiltedring analysis in Section 4.2 and v COR is the asymmetric driftcorrected velocity. The asymmetric drift correction σ D is given as where σ and ρ are the velocity dispersion and volume density of H I, and R is the radius of a galaxy. In particular, ρ can be converted to the H I surface density Σ by assuming an exponential distribution in the vertical direction and a constant scale height (for a first approximation). For the surface density Σ and velocity dispersion σ, we use the integrated H I (MOMENT0) and velocity dispersion (MOMENT2) maps, respectively.
Using the geometrical parameters of the LMC disk (center, position angle, and inclination) derived in Section 4.2, we derive the corrected radial profiles of Σ and σ, and the profile of Σσ 2 is given in Figure 8(c) as a black line. To avoid large fluctuations in the derivative in Equation (16), we fit Σσ 2 with an analytical function, where I 0 and R 0 are the fitted values in units of M e pc −2 km 2 s −2 and arcseconds, respectively. α is given in units of arcsec −1 . The analytic fit Σσ 2 (R) and the resulting asymmetric drift correction σ D are shown as the square symbols in Figures 8(c) and (d).
Lastly, the asymmetric drift-corrected bulk rotation curve is presented as the connected pentagon symbols in Figure 8(e). The asymmetric correction increases the rotation curve in the outer parts by ∼10 km s −1 .

Bulk Rotation Curves of the LMC
We present the H I bulk rotation curve of the LMC that has been corrected for the asymmetric drift as well as the transverse, nutation, and precession motions of the LMC in Figure 9. We also overplot the rotation curves derived from the H I IWM velocity field (MOMENT1) ) and from the carbon star kinematics (van der Marel et al. 2002). Compared to the former, the bulk rotation curve is more solidbody-like in the inner part, lower in amplitude by up to ∼30 km s −1 , and extends further in radius to ∼4.2 kpc. This significant velocity difference between the two curves, which is larger than the mean H I velocity dispersion value (3σ clipped MOMENT2) of ∼8.7 km s −1 , is mainly caused by the larger velocity gradient in the H I IWM velocity field than in the bulk velocity field. Additionally, the lack of correction for the nutation and precession of the H I IWM rotation velocities could contribute to the difference.
As discussed in van der Marel et al. (2002), the large velocity difference between the H I MOMENT1 rotation curve and the new bulk one could also be attributed to the large offset of ∼1°.2 in the kinematic center between the carbon stars and H I gas component. As discussed earlier, the kinematic center of the carbon stars is more consistent with the morphological center of the bar (van der Marel et al. 2002) than that of the H I IWM velocity field adopted in . This is likely due to the recent dynamical history of the LMC, the different ages for the two tracers, and their differing energy and angular momentum dissipation rates. In this respect, the global kinematic parameters of the LMC could not be reliably derived from the H I kinematic analysis alone. This is more serious for the H I velocity fields extracted in a conventional manner like the single Gaussian fit or IWM velocity field, which is highly affected by any noncircular disturbed gas motions in the galaxy.
On the other hand, as shown in Figure 9, the rotation curve derived from the carbon star kinematics is largely in fair agreement with the bulk rotation curve. In detail, the resulting bulk rotation velocities are well consistent with the carbon star rotation velocity out to ∼3.7 kpc from the kinematic center with discrepancies less than ∼10 km s −1 , which is comparable to the mean H I MOMENT2 (3σ clipped) value of ∼8.7 km s −1 . This indicates that the new H I bulk velocity field extracted using the profile decomposition method is able to minimize the effect of random noncircular gas motions and better trace the undisturbed global gas motions of the LMC disk that are expected from the carbon star kinematics. Unlike the sparsely distributed carbon stars in the disk of the LMC, the H I gas is more uniformly distributed in the disk and thus provides a much higher filling factor in the extracted velocity field. This allows us to derive more reliable rotation velocities in the inner part of the galaxy, reducing uncertainties. Instead, carbon stars are better at tracing the galaxy kinematics in the outer region beyond ∼5 kpc, where H I intensities become weak. As shown in Figure 9, the hybrid rotation curve of the LMC consisting of the H I bulk and carbon star rotation velocities is solid-body-like, linearly rising in the inner part of the galaxy. This is much lower than the previous H I tilted-ring analysis down to ∼30 km s −1 , indicating an even lower central gravitational potential for the dark matter halo component. Except for the last measured rotation velocity at ∼9 kpc, whose uncertainty is large, the hybrid rotation curve is likely to continue to rise at the radius of ∼7 kpc, reaching ∼60 km s −1 .

Summary
In this paper, we make a practical application of a newly developed profile decomposition technique-BAYGAUD, described by -to an H I data cube of the LMC from combined ATCA and Parkes observations in order to quantify the complexity of its ISM structure in terms of morphology and kinematics. For each line profile, using the new profile decomposition method based on Bayesian nested sampling, we fit a series of models constructed with different numbers of Gaussian components starting from 1 to a userspecified number. The model is then selected from among the competing models with respect to their Bayesian evidence, which is estimated from an efficient sampling of the parameter space using MULTINEST (Feroz & Hobson 2008), a Bayesian inference library.
We fit multiple Gaussian components (up to five) to individual H I line profiles of the LMC, and separate the bulk motion of the H I gas, which is corotating with the stellar disk as traced by carbon stars in van der Marel et al. (2002), from noncircular random or streaming motions. These disturbed gas motions are stirred up by star formation activity in the galaxy or interactions with the SMC and/or the Milky Way.
Together, the optimally decomposed H I gas components are classified into kinematically cold (1.6 < σ < 4.8 km s −1 ), warm (4.8 < σ < 18.1 km s −1 ), and hot (σ > 18.1 km s −1 ) with respect to their velocity dispersions. The corresponding H I gas masses of the kinematically decomposed components are about 12.2%, 58.3%, and 29.5% of the total H I mass, respectively. We find that they do not much differ in their velocity fields except for some localized and outer regions, which are presumably associated with active star formation (e.g., 30 Doradus) and interactions with the SMC and the Milky Way; the kinematically cold component appears to be distributed in the inner disk of the LMC, displaying a more regular rotation pattern. These interactions enhance the kinematically warm and hot H I gas components in the LMC in a way that results in larger velocity dispersions.
We derive the bulk rotation curves of the corotating H I gas disk associated with the global kinematics of the stellar disk traced by carbon stars presented in van der Marel et al. (2002). For this, we perform a 2D tilted-ring analysis using the H I bulk velocity field using a modified version of 2DBAT (Oh et al. 2018) to account for the effect of the LMC's transverse, nutation, and precession motions on the observed velocity field. These nonparametric fits of the 2D tilted-ring models are regularized using B-spline functions such as constant and cubic splines for radially varying ring parameters like position angle (f), inclination angle (i), and rotation velocity (v ROT ), except for the kinematic center (α 0 , δ 0 ) and systemic velocity (v SYS ), which are constrained with a constant single value over the galaxy radii.
The bulk rotation curves with different regularizations are in general consistent with each other except for the one with the cubic-splined position angle. The radial variation of position angle, which is well described by a cubic spline, is most likely to take into account the systematic streaming gas motions that are affected by the oval potential of the stellar bar of the galaxy. The change in position angle is pronounced within the stellar bar but the resulting rotation velocity is in good agreement with other bulk rotation curves derived with different regularizations within the estimated errors of rotation velocities.
We find that the velocity gradient of the extracted bulk velocity field is consistent with that of carbon star kinematics along the kinematic major axis but shallower than that computed from the previous classical intensity-weighted mean H I analysis. Localized deviations in the H I MOMENT1 velocity field are largely removed in the bulk velocity field, which shows more ordered motions. The H I gas bulk rotation curves are well consistent with those of carbon stars, which are sparsely distributed in the LMC disk but are not directly affected by hydrodynamical effects in the galaxy because of their relatively larger gravitational potential than the gaseous component. The carbon stars are expected to better trace the global disk kinematics of the galaxy than gaseous components, which tend to be vulnerable to such hydrodynamical effects in the galaxy.
The new H I bulk rotation curves are much shallower down to ∼30 km s −1 than the previous one derived from an H I MOMENT1 analysis ), which is not corrected for the random gas motions and the effect of precession and nutation of the LMC. The lower bulk rotation curves than the previous one indicate a lower dynamical contribution of dark matter in the inner part of the galaxy. Being in fair agreement with the carbon star kinematics of the disk, the newly derived bulk rotation curves are rather close to those of a solid body Figure 9. Upper panel: rotation curves of the LMC. Gray filled circles: the best fitting result of the 2D tilted-ring analysis derived using the bulk velocity field, corrected for the transverse, nutation, and precession motions as described in Section 4.2. Upward and downward triangles: the fit results derived using the receding or approaching side of the bulk velocity field, respectively. Filled pentagons: asymmetric drift-corrected H I bulk rotation curve derived in Section 4.2. Open circles connected by the dashed line: the rotation curve derived from the carbon star kinematics in van der Marel et al. (2002). Tripods: the rotation curve derived using the H I MOMENT1 map in . Filled circles: the rotation curve derived from the H I moment analysis in Indu & Subramaniam (2015). Lower panels: filled circles indicate the position angle (PA), inclination (INCL), and systemic velocity (VSYS) derived from the 2D tilted-ring analysis of the H I bulk velocity field. The dashed lines and tripods show the ones derived from the carbon star kinematics and the H I MOMENT1 analysis in .