© The Authors 2024

1. Introduction

The Large Magellanic Cloud (LMC) is the most massive Milky Way satellite galaxy (McConnachie 2012). Situated at a distance of about 50 kpc, it lies well within the Milky Way halo alongside its lesser counterpart, the Small Magellanic Cloud (SMC; ∼62 kpc). Easily visible with the naked eye, the relative proximity of the LMC and SMC makes them prominent features in the southern sky, and they have been known to mankind since the dawn of the southern civilisations (Dennefeld 2020).

The past orbit of the LMC and the exact epoch of arrival of the clouds in such a proximity to the Milky Way is difficult to discern, but they are likely on their first approach (Besla et al. 2007), as evidenced by the LMC’s high tangential velocities approaching the Milky Way escape velocity (Kallivayalil et al. 2006). This statement is further supported by the fact that the LMC is still associated with a group of its own satellites that have not been tidally stripped (Jethwa et al. 2016; Sales et al. 2017; Kallivayalil et al. 2018). However, Vasiliev (2024) propose a model where the clouds are on their second approach, with the first being 5 − 10 Gyr ago and with a pericentric distance larger than 100 kpc, far enough to retain the current system of satellites.

The total mass of the LMC has also been a subject of debate for many years. The close interactions with the SMC and the Milky Way have led to significant perturbations in the kinematics in the LMC outer regions, where reliable tracers are also scarce, making direct mass measurements almost impossible. But a plethora of recent studies using different independent methods (perturbations of stellar streams, Erkal et al. 2019; Vasiliev et al. 2021; Shipp et al. 2021; Koposov et al. 2023; momentum balance of the Local Group, Peñarrubia et al. 2016; census of the LMC satellites, Erkal & Belokurov 2020) have firmly put it in the range of 1 − 2 × 10¹¹ M_⊙, which is only five to ten times lower than the mass of the Milky Way (see the review by Vasiliev 2023). Most recently, the virial mass of the LMC was estimated by Watkins et al. (2024) to be 1.8 × 10¹¹ M_⊙ using globular clusters as tracers going out to 13.2 kpc and assuming a Navarro et al. (1997, NFW) dark matter (DM) halo density.

The LMC is still star forming and thus classified as a dwarf irregular (dIrr). However, owing to its high mass and recent arrival in the vicinity of the Milky Way, it possesses surprisingly regular features. To a first approximation, the LMC can be described morphologically as a barred spiral disc galaxy showing a typical rotation pattern, as demonstrated by Luks & Rohlfs (1992), Kim et al. (1998), Olsen & Massey (2007) for the gaseous component and by van der Marel et al. (2002), Olsen et al. (2011), Gaia Collaboration (2021a) for the stellar disc. van der Marel & Kallivayalil (2014) modelled the LMC as a flat rotating disc in three dimensions, and Vasiliev (2018) used axisymmetric Jeans dynamical models to simultaneously fit its rotation and velocity dispersion using Gaia Data Release 2 (DR2; Gaia Collaboration 2018a) proper motions (PMs). Niederhofer et al. (2022) showed that the stars in the bar region follow elongated orbits.

A more detailed overview of the morphological and kinematic structure of the LMC, however, revealed the effects of the close interactions with the SMC and the Milky Way. The LMC-SMC interactions likely formed the Magellanic Stream ∼2 Gyr ago (Besla et al. 2012; Diaz & Bekki 2012). The Magellanic Stream is a large body of stripped ionised and neutral gas that follows the orbits of the clouds and spans more than 200° on the sky. In addition, the two galaxies are connected through the Magellanic Bridge, a younger gaseous and stellar feature that probably originated from a direct collision between the two galaxies a few hundred million years ago (Diaz & Bekki 2012; Besla et al. 2013; Wang et al. 2019; Zivick et al. 2019). Multiple morphological and kinematic substructures discovered in the periphery of the LMC attest to its perturbed nature (Mackey et al. 2016, 2018; Belokurov & Erkal 2019; El Youssoufi et al. 2021; Cullinane et al. 2022a,b).

The LMC disc appears elongated (van der Marel 2001) and also shows deviations from a simple planar structure, as it is warped and truncated in the west direction towards the SMC (Mackey et al. 2018; Choi et al. 2018). The bar appears off-centre (de Vaucouleurs & Freeman 1972) and might be a transient unvirialised structure formed through the interactions and close encounters between the clouds (Zhao & Evans 2000; Besla et al. 2012).

The dynamical centre of the LMC has been equally difficult to establish, with discrepant results from photometry (van der Marel 2001), stellar (van der Marel & Kallivayalil 2014; Wan et al. 2020; Gaia Collaboration 2021a; Niederhofer et al. 2022), and gas (Luks & Rohlfs 1992; Kim et al. 1998) kinematics. The determination of the dynamical centre from stellar kinematics also depends on the type of data, line-of-sight (LOS) velocities versus PMs, and the type of stellar populations used as tracers (e.g. young versus old stars). It is worth mentioning that the kinematic centre has been found very close to the centre of the bar, based on Gaia DR3 PMs, which is the most complete and precise stellar PM database in the LMC today (Gaia Collaboration 2021a; Jiménez-Arranz et al. 2023).

The literature is also not conclusive about the LMC figure rotation nor its bar pattern speed. Dottori et al. (1996) studied the spatial distributions of two young massive LMC clusters and concluded that their differences can be attributed to an off-centred bar induced perturbation that propagates with Ω = 13.7 ± 2 km s⁻¹ kpc⁻¹. Gardiner et al. (1998) used N-body simulations to study the effects from perturbations due to an off-centred bar to the global distribution of the gas and star formation activity in the LMC. They found that a spiral structure compatible with the LMC emerges for bar pattern speeds in the range of 40 − 50 km s⁻¹ kpc⁻¹ while noting that their findings are in tension with Dottori et al. (1996). Shimizu & Yoshii (2012) derived Ω = 21 ± 3 km s⁻¹ kpc⁻¹ under the assumption that the Shapley Constellation III star forming region (Shapley 1951) was found in the Lagrange point L₄ of a rotating non-axisymmetric bar potential. Most recently, Jiménez-Arranz et al. (2024a) used several independent methods to estimate the LMC bar pattern speed using Gaia DR3 PM and LOS velocity data. They conclude that the Tremaine & Weinberg (1984) method is not applicable to the LMC, owing to its strong dependency on the orientation of the galaxy frame and the viewing angle of the bar perturbation. Using the Dehnen et al. (2023) method, they found Ω = −1 ± 0.5 km s⁻¹ kpc⁻¹ (a bar that barely rotates or even shows a marginal counter rotation) but acknowledge that such a configuration is incompatible with the spiral arm structure of the stellar disc. Finally, Jiménez-Arranz et al. (2024a) find that fitting a bi-symmetric velocity model (Gaia Collaboration 2023a) to the tangential velocity field yields a co-rotation radius R_c = 4.20 ± 0.25 kpc and Ω = 18.5 ± 1 km s⁻¹ kpc⁻¹. With a co-rotation to bar radius ratio of R_c/R_b = 1.8 ± 0.1, they argue for a slow bar in the LMC. In a separate work, Jiménez-Arranz et al. (2024b) used the KRATOS suite of N-body simulations to analyse the bar pattern speed of the LMC, finding Ω = 10 − 20 km s⁻¹ kpc⁻¹. They concluded that recent close interactions between the LMC and SMC do not necessarily alter the frame rotation rate and bar size and that long bars are typically slow.

Tahmasebzadeh et al. (2022) introduced a new method to derive the bar pattern speed using Schwarzschild orbit super-position dynamical models (Schwarzschild 1979) to fit the stellar motions in a triaxial bar potential with figure rotation. Initially demonstrated on a simulated galaxy, Tahmasebzadeh et al. (2024) used the method to measure the bar pattern speed in NGC 4371. Motivated by the non-axisymmetric kinematics in the inner region of the LMC, we test the performance of the orbit super-position dynamical models in the inner region of the LMC utilising Gaia DR3 LOS velocities, thus deriving the LMC bar pattern speed independently.

In this study we present a series of equilibrium dynamical models of the LMC that take advantage of the most recent 3D kinematic data in the LMC from Gaia DR3 (Gaia Collaboration 2023b). Although the LMC is in a complex dynamical environment and likely significantly out of equilibrium, we concur that our approach allows us to learn about and highlight kinematic and morphological features inherently not accounted for in equilibrium models. The article is organised as follows: In Sect. 2 we describe our dataset. In Sect. 3, we present three-dimensional axisymmetric Jeans dynamical models of the LMC, whose results establish the basis for our Schwarzschild models. In Sect. 5, we show a bar-disc decomposition of the LMC and derive the galaxy’s viewing angles and de-projected morphology, assuming co-centric axisymmetric disc and triaxial bar. In Sect. 6, we present our Schwarzschild models and the results for the LMC bar pattern speed. Sect. 7 summarises our work and provides an outlook of our plans for future LMC dynamical studies.

2. Data

For this study we utilise Gaia Data Release 3 (DR3¹Gaia Collaboration 2023b). We select Gaia DR3 sources with full 6D phase-space information: sky position, parallax, PMs, and LOS velocities from the Gaia radial velocity spectrometer (RVS), crossmatched with the VISTA Magellanic Clouds (VMC²) survey DR6 (Cioni et al. 2011). All our sources have both Gaia and VMC entries. There are ∼72 000 sources in our selection, which include foreground and background contaminants. We use this complete catalogue without removing the contaminants for our Jeans modelling (see Sect. 3), where we evaluate the probability of each source of being a genuine LMC member. However, in this section we apply the LMC membership criteria from Gaia Collaboration (2021a), based on parallax and PMs to illustrate some properties of the selected Gaia data (∼28 000 LMC sources). Throughout this article we adopt the widely used PM convention μ_RA ≡ μ_RA cos δ.

In Fig. 1 we show the Gaia CMD of the LMC stars with 3D kinematics, compared to a larger sample of ∼760 000 sources with high precision PM measurements (μ_err < 0.1 mas yr⁻¹), but no LOS velocity information, that we analysed in Kacharov & Cioni (2023) applying 2D Jeans models. It is evident that the stars, which have Gaia DR3 LOS velocities include only the brightest populations – intermediate age asymptotic giant branch (AGB) stars and young red supergiants. There are very few main sequence stars. In the middle and right panels of Fig. 1 we show the PM vector diagram of the selected sources, colour coded by the LOS velocity measurements and the μ_RA μ_Dec correlation coefficient. The gradient of the LOS velocity field shows that the galaxy rotation is well resolved. On the other hand, we do not observe any noticeable gradient in the μ_RA μ_Dec correlation coefficient in the PM vector plot. It is derived from the full 5-dimensional astrometric solution for the source and defines the proper motion error ellipse. It depends on the Gaia scanning pattern. Highly correlated solutions could contribute additional systematic uncertainties in the analysis. The mean correlation in our dataset is −0.03 with a standard deviation of 0.11, which we deem reasonably small. We treat the μ_RA, μ_Dec, and V_LOS as independent in our 3D Jeans modelling.

Fig. 1. Data selection. Left panel: Gaia CMD of the selected LMC targets with 3D kinematics (blue). Middle panel: PM vector diagram of the LMC sources with 3D kinematics, colour-coded by the LOS velocity. Right panel: Same as the middle panel but colour-coded by the μRA μDec correlation coefficient. The grey points in all panels show additional LMC sources with good Gaia PM measurements (μerr < 0.1 mas yr−1) but no LOS velocity entries. These are the sources used in our 2D Jeans model.

In the left panel of Fig. 2 we show the density distribution of the selected LMC sources, where notable features of the galaxy morphology are recognisable like its dense stellar bar and the spiral structure of the disc. We also plot the spatial distribution of the μ_RA μ_Dec correlation in the middle panel that can be directly compared to the stellar density distribution. The PM correlation map can be directly traced to the Gaia scanning pattern (see Fig. 4 in Gaia Collaboration 2021a), but there is no dependence with the stellar crowdedness or any other LMC morphological features. Finally we show the renormalised unit weight error (RUWE) spatial distribution of the LMC sources in the right panel of Fig. 2. The RUWE value provides an indication of the quality of the astrometric solution, which affects both the positions and PMs of sources. Higher RUWE values suggest larger systematic errors in the astrometric measurements, potentially leading to less reliable PM estimates. The median RUWE index for the selected ∼28 000 LMC stars is 1.026, with 90% of the targets having a RUWE less than 1.18. However, we highlight the systematic increase of the RUWE index in the central most crowded bar region, which indicates possible elevated systematic uncertainties in the astrometric solution there. About 4% of the LMC member stars have RUWE > 1.4 and they are located predominantly in the bar.

Fig. 2. Spatial distribution and quality control. Left panel: Density distribution of the selected LMC targets with 3D kinematics (blue) with the larger sample of stars with good PM measurements plotted under with light grey dots. Middle panel: Similar to the left panel but the stars with 3D kinematics are colour-coded according to the μRA μDec correlation coefficient. Right panel: Stars with 3D kinematics colour-coded according to the RUWE index.

Throughout this work we assume a distance modulus to the LMC of 18.49 mag (de Grijs et al. 2014; Crandall & Ratra 2015), derived from a large compilation of different tracers, which corresponds to a distance of 49.9 kpc. This allowed us to convert the Gaia DR3 PMs and their uncertainties to linear velocities and make a fair comparison between the velocity errors across the three spatial dimensions. Figure 3 shows the individual velocity uncertainties along the LOS and in the plane of sky as a function of Gaia G magnitude. The median RV error in our LMC sample is 3.2 km s⁻¹ with 90% of the stars having RV errors less than 6.9 km s⁻¹. The PM uncertainties along the RA and Dec directions trace each other very closely and are noticeably larger than the RV uncertainties with a median μ_RA error 5.1 km s⁻¹ (0.021 mas yr⁻¹), μ_Dec error 5.3 km s⁻¹ (0.022 mas yr⁻¹) with 90-th percentiles, respectively 9.0 km s⁻¹ (0.038 mas yr⁻¹) and 9.6 km s⁻¹ (0.041 mas yr⁻¹).

Fig. 3. Velocity error distribution as a function of GaiaG magnitude along the LOS and in the plane of the sky. The PMs have been converted to linear velocities assuming an LMC distance of 49.9 kpc. The densest regions are represented only with contours due to the significant overlap between the points.

3. Axisymmetric Jeans Models

3.1. Model setup

In Kacharov & Cioni (2023) we presented a series of discrete axisymmetric dynamical models of the LMC with varying gravitational potentials, based on the Jeans (1915) hydrodynamical equations, using ∼1 × 10⁶ PM measurements in the same VMC field (the complete sample down to GaiaG = 18 mag shown in Fig. 1, but with contaminants included). We found that these models provide realistic estimates of the enclosed total mass within the extent of the data and describe well the rotation and velocity dispersion of the LMC disc. However, they do not capture well the kinematic properties of the central region of the galaxy, dominated by the LMC bar – an intrinsically non-axisymmetric feature. In addition, when using solely PM data, it is not feasible to fully constrain the orientation of the line of nodes and the pointing of the angular momentum vector, which leads to uncertainties in the overall shape and morphology of the galaxy.

In this article we present a similar setup of discrete Jeans axisymmetric models of the LMC, but using only stars, for which we have 3D Gaia kinematics. The inclusion of RV measurements allowed us to fully constrain the morphology and orientation of the LMC disc, purely from its kinematical properties, albeit still bound to the assumption for axisymmetry, which is intrinsic to the solution of the Jeans equations (Cappellari 2008).

The Jeans equations, derived from the collisionless Boltzmann equation, link the gravitational potential of a self-gravitating system to its stellar density distribution and stellar velocity moments. In order to solve them we need to make assumptions for the gravitational potential and the density distribution and parametrise them.

For the stellar density we assume that the galaxy surface brightness follows the projected stellar density and that it is exponential. We adopt the LMC exponential surface brightness profile from Gallart et al. (2004) with a scale radius r_e = 98.7 arcmin and scale it to a total luminosity L_LMC = 1.31 × 10⁹ L_⊙ from the NASA/IPAC Extragalactic Database (NED). While we do keep the above figures fixed in the dynamical model fit, we vary the projected flattening (q_proj), defined as the projected minor to major axis ratio, the orientation angle (ψ), which gives the direction of the LMC’s projected minor axis on the sky measured from north to east), and the inclination angle (θ), which is defined to be 0° if seen face-on and 90° if seen edge-on. They are free parameters in our setup. The intrinsic flattening can then be calculated as $q_{\mathrm{intr}}=\sqrt{q_{\mathrm{proj}}^2-cos{\theta }^2}/sin\theta$ in the axisymmetric case.

For the gravitational potential we have two cases. In one case, we assumed that the total mass of the LMC is contained within a spherical DM halo that follows a NFW mass density profile. This means that mass contributions from the stellar and gaseous component are also included in the DM halo and not treated separately. The NFW radial density profile (Eq. 1) is characterised by two parameters, a characteristic density (ρ_s) and a scale radius (r_s), that are free parameters in our Jeans dynamical model:

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_s}{\left(\frac{r}{r_s}\right){(1+\frac{r}{r_s})}^2}·\end{array}$$(1)

In the second case, we have a mixed potential that consists of an NFW DM halo and a visible mass component. The latter follows the exponential surface brightness distribution and an assumed stellar mass to light ratio (M/L_*). We explored models where M_*/L is a free parameter in the model or fixed to M_*/L = 1.5 M_⊙/L_⊙ for a total luminous mass of 2 × 10⁹ M_⊙. We represent both the mass density and surface brightness distributions as multi-Gaussian expansions (MGEs; Emsellem et al. 1994), which makes them convenient to project and de-project – operations necessary to solve the Jeans equations.

We used the PYTHON version of the Jeans Anisotropic Modelling (JAM) code by Cappellari (2008, 2012) and adopt a Bayesian framework with a discrete likelihood function, initially introduced by Watkins et al. (2013) (but also see Zhu et al. 2016; Kamann et al. 2020; Kacharov et al. 2022). In the discrete likelihood framework, the JAM code predicts the first ($\overline{v_z}$, $\overline{v_r}$, $\overline{v_{\varphi }}$) and second ($\overline{v_z^2}$, $\overline{v_r^2}$, $\overline{v_{\varphi }^2}$) velocity moments at the position of each star in the sample in the three cardinal directions – LOS (z), radial (r), and tangential (ϕ). We compute the probability of each star’s three observed velocity components to be drawn from a 3D Gaussian with the first moments as the mean values in the three cardinal directions of the model and corresponding variances ${\sigma }^2=\overline{v^2}-{\overline{v}}^2$, also in all three directions and call it the dynamical likelihood ($P_{\mathrm{dyn}}^{\mathrm{LMC}}$). We assume that there is no cross-term correlation between the two PMs and the RV. This setup requires two additional model parameters, which are set free in our fit – the velocity anisotropy (${\beta }_z=1-\frac{\overline{v_z^2}}{\overline{v_R^2}}$) and a systemic rotation parameter, defined as $\kappa =\frac{\overline{v_{\varphi }}}{{(\overline{v_{\varphi }^2}-\overline{v_R^2})}^{1/2}}$. When κ = ±1 and β_z = 0, the system reduces to an isotropic rotator and when κ = 0, there is no net angular momentum. Negative β_z values indicate tangential anisotropy, while positive β_z values indicate radial anisotropy (Cappellari 2008). Although both the anisotropy and the rotation parameter can vary with radius, here we assume that they are radially constant.

We also included in our Jeans models a combined foreground and background population component with a flat surface density designed to capture the foreground stellar and background galaxy contaminants in our Gaia DR3 data. This component is controlled by a single free parameter (ϵ). It is defined as a fraction of the central surface brightness of the LMC MGE, as in Watkins et al. (2013), Zhu et al. (2016), Kacharov et al. (2022). Hence, we can write the spatial probability distribution functions of the LMC and contaminant population components as

$$\begin{array}{c}\hfill P_{\mathrm{spa}}^{\mathrm{LMC}}=S/(S+ϵC)\end{array}$$(2)

and

$$\begin{array}{c}\hfill P_{\mathrm{spa}}^{\mathrm{cont}}=1-P_{\mathrm{spa}}^{\mathrm{LMC}},\end{array}$$(3)

where S is the MGE surface brightness profile at the position of each star and C is the central surface brightness of the LMC MGE. We also define $P_{\mathrm{dyn}}^{\mathrm{cont}}$ analogically to $P_{\mathrm{dyn}}^{\mathrm{LMC}}$, but with 0 mean motion and very large variance.

The posterior probability of a star (i) to belong to the LMC or to the foreground component of the model is then given by the joint probability of the above specified two probability distributions, and we maximised the log posterior function of the entire sample:

$$\begin{array}{c}\hfill lnP={\displaystyle \sum _i}ln(P_{\mathrm{spa},i}^{\mathrm{LMC}}{P}_{\mathrm{dyn},i}^{\mathrm{LMC}}+P_{\mathrm{spa},i}^{\mathrm{cont}}{P}_{\mathrm{dyn},i}^{\mathrm{cont}}).\end{array}$$(4)

Thus, when the best-fit model is obtained, we can easily compute the probability for each star to be a genuine LMC member or to belong to the contaminant population. Finally, we had five additional trivial model parameters that describe the kinematic centre of the LMC (α₀, δ₀) and its mean space motion (v_los⁰, μ_α⁰, μ_δ⁰), which brings the total number of free parameters in our Jeans model to 13.

We used the EMCEE affine invariant Markov chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010) in PYTHON (Foreman-Mackey et al. 2013). We ran our models on a CPU cluster engaging 96 cores and used 192 walkers and 3000 steps in the MCMC, which we confirmed to be enough for the fit to converge. We converted the stellar coordinates and PMs to an orthographic projection (α, δ, μ_α, μ_δ → X, Y, μ_X, μ_Y, see Eq. (1)–(3) in Gaia Collaboration 2021a); corrected the three velocity vectors for perspective effects, as in van de Ven et al. (2006, Eq. (6)); and rotated the coordinate system (X, Y, μ_X, μ_Y → X_p, Y_p, μ_Xp, μ_Yp) every time we called the likelihood function to account for the incremental random changes in the kinematic centre, mean space motion, and disc orientation during the fit.

3.2. Model results: LMC orientation

The best-fit parameters from our 3D Jeans model are summarised in Table 1. We also include there the results from our 2D Jeans model, based solely on Gaia PM data and described in Kacharov & Cioni (2023, this model includes the grey points in Fig. 1 as kinematic tracers). The way the 2D and 3D Jeans models are constructed and fit are essentially the same, which makes them directly comparable. The only difference (besides not including RV data in the 2D model) is that we kept the centre fixed to the best-fit coordinates from Gaia Collaboration (2021a) in the 2D model, while here we fit for it, using the added knowledge of the LOS velocities. The 2D model is also fit to significantly more and fainter stars with good PM data down to GaiaG 18 mag, compared to the 3D model (see Fig. 1). This may make the two models sensitive to different stellar populations in the LMC; the 2D model being more sensitive to the kinematics of old stars, represented by the numerous, but faint red giant branch (RGB) population.

The main difference between the 2D and 3D model results concerns the orientation of the LMC disc in space, its angular momentum vector, and the corresponding viewing angles from Earth. The difference in the angle ψ, which describes the orientation of the rotation axis in the plane of the sky is 20° (124° for the 3D model versus 104° for the 2D model). The difference in the best-fit inclination angle θ is 8.5° ±0.2°: θ = 34.0° ±0.1° for the 2D model, very close to the Gaia Collaboration (2021a) estimate of θ = 33.3°, versus θ = 25.5° ±0.1° for the 3D model. We note, however, that the Gaia Collaboration (2021a) study makes an assumption for a flat LMC disc, which makes their inclination angle an upper limit. They find a line of nodes orientation angle ψ ∼ 130°, closer to our 3D model result. The discrepancy between the 2D and 3D results also leads to a significant difference in the projected flattening of the galaxy – q_proj = 0.91 ± 0.01 for the 3D model versus q_proj = 0.84 ± 0.01 for the 2D model. Clearly, adding V_LOS data to the Jeans model fit of the LMC disc significantly changes how we perceive the galaxy from a purely kinematic perspective. These results can be compared to photometric constraints, which we do in Sect. 5. We note that the resulting intrinsic flattening of the LMC disc is essentially estimated to be the same according to the 2D and 3D fits q_intr ∼ 0.24. Previous literature studies put the orientation of the line of nodes at 120° −155° and the inclination angle in the range of 25° −40° (van der Marel & Kallivayalil 2014; Gaia Collaboration 2021a, and references therein). Our 3D kinematics analysis is consistent with these constraints, albeit on the lower end. On the other hand the 2D kinematic analysis yields a line of nodes orientation angle closer to the photometric major axis and an inclination angle close to the middle of the distribution, but clearly incompatible with the LOS velocity field. While the 3D velocities provide the most robust kinematic constraint of the galaxy orientation, PMs alone can also provide information about the line of nodes orientation and inclination assuming rotation of a thin axisymmetric disc. The formalism for this is given in Gaia Collaboration (2018b). In the case of the LMC, it appears that the line of nodes orientation constrained from PMs alone and the axis of maximum V_LOS gradient are offset by 20°. This discrepancy reflect the complex and perturbed nature of the LMC disc. Furthermore, the mismatch between the orientation of the LMC’s morphological major axis and the axis of maximum V_LOS gradient is a known issue, which prompted van der Marel & Cioni (2001), van der Marel & Kallivayalil (2014) to argue for an elongated non-axisymmetric disc in the LMC.

The left panel of Fig. 4 shows an infrared CMD from the VMC survey of all stars in our 3D velocity sample that have a high probability (P_mem > 0.8) of being genuine LMC member stars, according to our 3D axisymmetric Jeans dynamical model fit. Similarly, the right panel of Fig. 4 shows the infrared CMD of the high-probability contaminant sources (foreground Milky Way stars and background galaxies). They have a very flat spatial distribution, as seen from the inlet of the figure, showing that the discrete Jeans model is very successful at separating genuine LMC member stars from foreground/background contamination. We find 28 000 stars with LMC membership probability > 0.8. In comparison, applying the Gaia Collaboration (2021a) membership selection criteria, based on parallax and PM cuts, we are left with 27 653 LMC stars, which is an excellent agreement.

Fig. 4. LMC membership and population selection. Left panel: Infrared CMD from the VMC survey of the high-probability LMC member stars according to the 3D Jeans model fit. Right panel: Infrared CMD of the high-probability contaminant sources (foreground Milky Way stars and background galaxies) according to the 3D Jeans model fit. The inlet shows their spatial distribution in the VMC footprint. The dashed line is the same in both panels. It is used to separate the old and young LMC stellar populations.

In order to investigate whether the difference in the inferred orientation of the LMC, based on the 3D and 2D kinematics, could be alleviated by using a more homogeneous stellar population, we fitted the same Jeans model to 3D kinematic data from only the old stellar population and only the young stellar population, separately. We used the same photometric cuts as in El Youssoufi et al. (2019), Cioni et al. (2019) in the infrared VMC CMD to separate the stars into old and young populations (see Fig. 4). The young population includes main sequence stars, as well as blue, yellow, and red supergiants, while the old population consists of AGB and upper RGB stars. The results from these two additional models are also summarised in Table 1. They are essentially in line with the model that uses the 3D kinematic data of all stars.

Thus, the significant discrepancy of the LMC viewing angles inferred from the 3D and 2D Jeans axisymmetric models (the former are also in tension with photometric constraints) are not due to stellar population differences, but speak of a non-negligible deviation of the LMC disc from axisymmetry (see also van der Marel & Cioni 2001; van der Marel & Kallivayalil 2014).

Figure 5 shows the spatial distribution of the high-likelihood genuine LMC members with 3D kinematics in an orthogonal coordinates system, where the projected rotation axis is vertical and the line of nodes horizontal. We show the best-fit ellipse that describes the projection of the assumed axisymmetric LMC disc, together with its bulk PM direction according to the 3D Jeans model. We show the locations of the inferred kinematic centres from the three Jeans models (based on 3D kinematics of all stars, and the old, and young populations separately), as well as the 2D kinematic centre published in Gaia Collaboration (2021a), with respect to the photometric centre inferred from the Gaia DR3 density map of all sources (see Sect. 5). We find that the distance between the 3D kinematic centre inferred from all stars and the photometric centre is 16.6 arcmin. The 2D kinematic centre is closest to the photometric centre with an angular separation of 11.7 arcmin.

Fig. 5. Distribution of the high probability LMC member stars according to our 3D Jeans model (grey dots), rotated at angle ψ = 124°, so that the best-fit projected rotation axis is vertical and the major axis horizontal. The best-fit kinematic LMC centres are indicated with crosses, based on all stars with 3D velocities (blue cross), old stars (red cross), young stars (yellow cross), and the centre based on 2D kinematics from Gaia Collaboration (2021a, cyan cross). The large ellipse indicates the best-fit projected shape of the LMC disc, based on the axisymmetric Jeans model. The small ellipse indicates the shape of the bar, based on a photometric surface density fit (see Sect. 5). The dashed lines indicate the photometric centre, based on the surface density fit. The solid black line shows the bulk PM of the LMC.

3.3. Model results: LMC kinematics

While the Jeans model fits were performed on the discrete individual velocity measurements, we show the model performance using binned Voronoi maps of the Gaia DR3 kinematic data. We divided our data in Voronoi bins, each containing 150 stars with 3D kinematic information and converted the PMs to linear velocities using the same adopted distance to the LMC throughout this article. We also subtracted the systemic PM of the LMC, as fitted by our model. For each bin we inferred the mean velocity, velocity dispersion, and their uncertainties in the three cardinal directions using a maximum likelihood approach that takes into account the individual velocity errors.

Figure 6 shows the observed projected 3D rotation field of the LMC from the Gaia DR3 data, together with the results from our best-fit Jeans model, and the model residuals. We kept the velocity scale consistent in all maps for easier comparison of the results. There are several things worth noting from these kinematic maps. The Jeans model fits very well the projected rotation of the LMC disc in all three directions, as shown by the residual maps, which are very flat and do not show much systematic structure. In the central region, however, the data shows a noticeable twist in the rotation axis, aligned with the LMC bar, which is not caught by our axisymmetric model, by definition. Most of the rotation happens in the plane of the sky, as evident from the large amplitude difference between the LOS and the V_X, V_Y velocity fields, confirming the smaller inclination angle that we find in this analysis, compared to other studies, based only on PM data (e.g. Gaia Collaboration 2021a). We plot the projected rotation curve of the LMC in Fig. 7. The rotation amplitude in the plane of the sky reaches ∼75 km s⁻¹ at 6 kpc from the centre in the plane of the sky, while along the LOS direction it is only ∼35 km s⁻¹ at the same distance. We sample the rotation curve from the Jeans model posterior distribution (shown with red lines in Fig. 7) to find that the Gaia DR3 data and its uncertainties provide very tight constraints to the Jeans solution with almost no room for variation within the model limits.

Fig. 6. Voronoi maps showing the 3D Jeans model fit to the LMC rotation field from the Gaia DR3 kinematic data. The first column represents the data in the X, Y, and LOS directions, respectively. Each Voronoi bin is colour-coded according to the median velocity of all stars that fall into it. The second column represents the Jeans model results and the third column the model residuals. PMs are converted into linear velocities. The velocity scale is kept unchanged between all maps for consistency and easier comparison.

Fig. 7. Projected rotation curve of the LMC in the three cardinal directions. The blue points show the maximum likelihood mean velocities in each of the Voronoi bins defined in Fig. 6. The yellow points show the corresponding mean velocities of the best-fit 3D Jeans model in the same Voronoi bins. The read lines show 150 random draws from the model posterior distribution along the horizontal axis (vertical axis in the case of VX).

Analogically, Fig. 8 shows the Jeans model match to the projected velocity dispersion map of the LMC in the three cardinal directions and Fig. 9 shows the projected velocity dispersion profiles as a function of radius. Again we have kept the velocity dispersion range consistent in all panels of Fig. 8 for better clarity. The observational maps are quite anisotropic and show different velocity dispersion in the three directions with the X_p having the highest dispersion in the central region and the LOS direction – the lowest. The Jeans model tries to accommodate this with a relatively high anisotropy parameter ${\beta }_z\propto 1-\frac{{\sigma }_z^2}{{\sigma }_R^2}=0.7\pm 0.01$. The highest velocity dispersion in the plane of the sky is observed along the LMC bar, which is a non-axisymmetric feature and cannot be precisely reproduced by our Jeans model. The Gaia DR3 data shows that the velocity dispersion can reach above 50 km s⁻¹ in the plane of the sky, while only reaching ∼24 km s⁻¹ in the LOS direction. Our best-fit Jeans model predicts central velocity dispersions σ_X = 42 km s⁻¹, σ_Y = 40 km s⁻¹, and σ_LOS = 28 km s⁻¹.

Fig. 8. Voronoi maps showing the 3D Jeans model fit to the LMC velocity dispersion field from the Gaia DR3 kinematic data. The first column represents the data in the X, Y, and LOS directions, respectively. Each Voronoi bin is colour-coded according to the maximum likelihood velocity dispersion of all stars that fall into it, taking into account their individual velocity uncertainties. The second column represents the Jeans model results and the third column the model residuals. PMs are converted into linear velocities. The velocity scale is kept constant amongst all maps for consistency and easier comparison.

Fig. 9. Projected velocity dispersion profile of the LMC in the three cardinal directions. The blue points show the maximum likelihood velocity dispersion in each of the Voronoi bins defined in Fig. 8. The yellow points show the corresponding velocity dispersions of the best-fit 3D Jeans model in the same Voronoi bins. The read lines show 150 random draws from the model posterior distribution along the horizontal axis.

3.4. Model results: LMC potential

Finally, we discuss the gravitational potential of the LMC, as constrained by the best-fit NFW only gravitational potential, characterised by ${\rho }_s=1\pm 0.05\times {10}^{-2}{\mathrm{M}}_{\odot }$ pc⁻³ and r_s = 10.5 ± 0.4 kpc and the mixed potential (visible + DM mass), characterised by ${\rho }_s=3.4\pm 0.3\times {10}^{-3}{\mathrm{M}}_{\odot }$ pc⁻³, r_s = 20.0 ± 1.5 kpc, and a fixed M_*/L = 1.5 M_⊙/L_⊙. We note that if we let the mass-to-light ratio to be a free parameter in our MCMC fit, it tends to go to 0, so the overall best solution is the DM only model. Figure 10 shows the derived cumulative mass, mass-to-light (M/L) ratio, and the circular velocity curve according to our Jeans dynamical models with both types of gravitational potential. The most distant kinematic tracers in our data selection are at a projected distance of 6.2 kpc from the centre of the LMC and both potentials have virtually undistinguishable mass properties within this radial extent. We find an enclosed mass of $M_{\mathrm{encl}}=1.36\pm 0.01\times {10}^{10}{\mathrm{M}}_{\odot }$ for the DM only model versus $M_{\mathrm{encl}}=1.34\pm 0.01\times {10}^{10}{\mathrm{M}}_{\odot }$ for the model with a mixed potential (see Table 1).

Fig. 10. Cumulative mass (left panel), mass-to-light ratio (middle panel), and rotation curve (right panel) of the LMC according to best-fit NFW only (black) and mixed (magenta) gravitational potentials of the axisymmetric Jeans model. The vertical dashed lines denote the extent of the kinematic tracers and the viral radius (R200) and the horizontal lines denote the enclosed mass and mass-to-light ratio within these radii. In the rotation curve plot, we overlay the total mean velocity of each Voronoi bin defined in Fig. 6 as orange dots, as well as the total velocity of all LMC member stars as a transparent blue cloud.

However, the difference between both choices of gravitational potentials becomes noticeable if we extrapolate them beyond the extent of the kinematic tracers. The DM only model has a virial mass $M_{200}=1.81\pm 0.01\times {10}^{11}{\mathrm{M}}_{\odot }$ at a virial radius r₂₀₀ = 75 ± 2 kpc. This virial mass estimate is well in line with other independent LMC mass estimates from literature (see Vasiliev 2023), albeit on the higher end, and matches precisely the estimate by Watkins et al. (2024). The latter study and our models assume the same NFW shape of the gravitational potential, but use kinematic tracers at different radial scales. While we focus on the inner region stellar kinematics of the LMC disc and bar, Watkins et al. (2024) use halo GCs out to 13.2 kpc. The fact that both studies find the same estimate for the LMC virial mass is astounding and shows that the NFW DM density approximation is reasonable at these radial scales. But these total LMC mass estimates should not be taken at face value and must be treated with caution. The Jeans model with a mixed potential yields a significantly higher virial mass of $M_{200}=4.3\pm 0.5\times {10}^{11}{\mathrm{M}}_{\odot }$ at a virial radius r₂₀₀ = 156 ± 6 kpc, clearly showing the dangers of extrapolating these results beyond the extent of the kinematic tracers. To accommodate for the visible mass component that follows the adopted exponential distribution of the LMC surface brightness, the NFW halo in the mixed potential has a lower central density, but is significantly more extended than in the DM only case. Both virial radius estimates are significantly larger than the LMC-Milky Way distance. Our mass estimates are mostly driven by the strong assumption that the DM density of the LMC follows a strict NFW profile, which is most certainly not true that far out.

On the other hand, the virial mass derived from our DM only 2D Jeans model is about twice lower than from the DM only 3D Jeans model (see Table 1 and Kacharov & Cioni 2023), putting it on the lower end of LMC independent mass estimates from the literature. The enclosed mass within the maximum radial extent of our kinematic tracers is much more reliable and the numbers are in a reasonable agreement between all five Jeans models that are outlined in Table 1.

The unrealistically small statistical uncertainties are due to the good constraints that the kinematic tracers provide to the scale density and scale radius of the NFW profile, but they are also very model dependent and the systematic errors from the choice of the potential and surface brightness parametrisation will dominate. One could also see in Table 1, that changing the number and type of kinematic tracers will change the derived parameters more than their formal 1σ statistical uncertainties.

In the middle panel of Fig. 10 we plot the total M/L ratio as a function of radius for both types of gravitational potentials. The total M/L ratio is computed by dividing the cumulative total mass by the cumulative brightness. As mentioned in the previous section, we use an exponential surface brightness profile with a fixed total luminosity of 1.3 × 10⁹ L_⊙. We have a central M/L = 2.17 ± 0.03 M_⊙/L_⊙ for the DM only fit, which is typical for an old stellar population, but there is no distinction between stellar and DM contribution in this model. The mixed potential has a central M/L = 2.35 ± 0.03 M_⊙/L_⊙, which shows that the DM dominates the potential under the model assumption of a stellar M_*/L = 1.5 M_⊙/L_⊙. The M/L ratio gradually rises reaching ∼11 M_⊙/L_⊙ at 6.2 kpc – the maximum radial extent of our kinematic tracers. This radius contains 93% of the total LMC luminosity according to our exponential parametrisation. At the virial radius the M/L ratio is 140 ± 10 M_⊙/L_⊙ according to our fiducial DM only parametrisation, emphasising that the LMC resides in a very massive DM halo.

In the right panel of Fig. 10 we plot the circular velocity curve of the LMC, as computed from the best-fit NFW only profile, which reaches a maximum of 115 ± 2 km s⁻¹ and 90 km s⁻¹ at 5 kpc in excellent agreement with the result from Vasiliev (2018). The resulting rotation curve from the mixed potential follows very closely the fiducial DM only one within the extent of our kinematic tracers, but extrapolated outwards, it reaches a maximum of 126 ± 4 km s⁻¹ with the caveats discussed above. We compare the rotation curve to the total mean velocities of each of the Voronoi bins defined in Fig. 6. They all follow tightly the circular velocity curve and lie on, or under it if the Voronoi bin is far from the line of nodes (maximum velocity gradient). This assures us that the NFW parametrisation of the gravitational potential is adequate at least within the radial extent of the kinematic data.

Overall the 3D Jeans model fits provide a good description of LMC disc’s kinematics and the galaxy’s fundamental dynamical parameters, but to do better in the central region we need to work with triaxial models, which can describe the kinematics of the LMC bar.

4. Probing the axisymmetry of the LMC from 3D kinematics

We can probe the overall axisymmetry of the LMC by examining the correlation between the mean velocities in the LOS and Y directions. In an axisymmetric system they are related via the equality ⟨V_LOS⟩= − ⟨V_Y⟩ tanθ (van de Ven et al. 2006). We measure these quantities in Voronoi bins, similar to the ones shown in Figs. 6 and 7, which are for one particular orientation – the one derived from the 3D Jeans dynamical model. We estimate the residual variance between these two velocity components by fitting a linear function, assuming different line of node orientations (Fig. 11). The slope of the line corresponds to the inclination angle of the galaxy. We do this exercise for all LMC member stars, as well as for the old and young LMC populations separately, according to the age separation defined in Fig. 4. We find a best fit line of nodes orientation angle ψ = 122° with a corresponding inclination angle θ = 26°, when considering all LMC member stars. The old population yields ψ = 120° and θ = 25.5°, while for the young stars we find ψ = 126° and θ = 28°. These estimates are practically identical to the results obtained from the 3D Jeans model.

Fig. 11. LMC axisymmetry tests. Left panels: Residual variance of the mean VLOS versus the mean VY correlation for different line of nodes orientations, colour-coded by the inclination angle for all stars (top), old stars (middle), and young stars (bottom) rows, respectively. Middle panels: Line of nodes orientation versus inclination angle, colour-coded by the residual variance for the old stars and separated into old and young populations. The shaded regions in these panels correspond to the best fit values from the 3D Jeans model and their uncertainties. Right panels: Mean VLOS versus mean VYtan(θ) for the LMC orientation with minimum residual variance for each respective population. These plots are colour-coded by distance from the centre of each Voronoi bin, and the dashed lines indicate a perfect anti-correlation.

The resulting correlation between the mean V_LOS and the mean V_Y at the estimated line of nodes angle for the three cases that minimises the scatter is also shown in Fig. 11. We colour-coded the separate velocity Voronoi bins by their distance from the galaxy’s centre to check whether the inner bins, most influenced by the kinematic effects of the bar, are outliers in this plot. We do notice that the most central regions (dominated by the bar) are systematically offset from the line for all stars, while this effect is almost entirely mitigated, when considering only the old population. On the other hand, the young population shows the largest offset from the perfect correlation line, which indicates a more significant deviation from axisymmetry. This is not surprising, given the non-uniform and patchy spatial distribution of these young stars, which have likely not reached an equilibrium state.

In addition, the most distant Voronoi bins from the LMC centre in our footprint, located on the approaching side of the disc (and on the side of the SMC) also show an increased scatter in the old stellar population, likely an effect from the interaction between the three galaxies.

Overall, we did not see a distinct effect induced from the non-axisymmetry of the LMC bar in these plots. The bar is prominent in both the old and young stellar populations. The kinematic effect of the bar can mostly be seen when looking directly at the Voronoi maps in Fig. 6, as a twist in the velocity gradient direction in the bar region. This effect, however, appears to follow the same pattern in both the ⟨V_LOS⟩ and ⟨V_Y⟩ Voronoi maps, keeping the correlation between them.

5. De-projecting the LMC surface density distribution assuming a triaxial bar

The goal of this study is to at least partially abandon the common assumption for axisymmetry when it comes to the dynamical modelling of the LMC and include a triaxial bar, along with an axisymmetric disc. To achieve this, we need a new parametrisation of the LMC surface brightness distribution, which in the case of the Jeans models, was assumed to be purely exponential. Here we provide a new two-component photometric decomposition of the Gaia DR3 stellar density map of the LMC³ (Gaia Collaboration 2021b) using GALFIT⁴ (Peng et al. 2010), which we use to de-project the galaxy and calculate its viewing angles.

5.1. Photometric decomposition of the LMC

The Gaia DR3 density map of the LMC is shown in the left panels of Fig. 12. We have obtained it as a fits image with pixel resolution 1271 × 991 px and pixel scale 81.84 arcsec px⁻¹. The fits file contains a WCS coordinate solution, which allowed us to convert pixels to equatorial coordinates. We convolve this image with a Gaussian PSF kernel of 5 px and complement it with a corresponding Poissonian error image (each pixel has the square root of the value in the original image) for the GALFIT analysis.

Fig. 12. LMC disc – bar decomposition. From left to right: Gaia DR3 density map of the LMC, the best-fit GALFIT model, the GALFIT model residuals, and the Gaia DR3 LOS Voronoi map in the same orientation. The white dashed lines denote the best-fit major and minor photometric axes that cross at the best-fit photometric centre according to the GALFIT solution. The top row is the best-fit solution when all GALFIT parameters are left free; the bottom row shows the GALFIT solution when the orientation of the large photometric axis and projected flattening are set to the orientation of the line of nodes and projected flattening according to the 3D Jeans model and using the adopted Sérsic indices of the disc and bar nd = nb = 0.5.

Then we proceed with the GALFIT fit and morphological decomposition of the LMC. We look for a two component solution, employing two Sérsic profiles – one describing the LMC disc and the other one – the LMC bar. We also allowed for a background component with a spatial gradient. The main free parameters of the fit are summarised in Table 2. Each of the two Sérsic profiles is characterised by five parameters – integral magnitude (M), scale radius (R), Sérsic index (n), axis ratio (q), and position angle (ψ). In addition we fit for the centre of the LMC, requiring that the bar and the disc share the same centre. The latter requirement might be viewed as controversial, since many studies have reported an off-centred bar for the LMC (e.g. Besla et al. 2012), but it is necessary for running our triaxial dynamical models. Our study was especially motivated by the close coincidence of the PM kinematic centre of the LMC as inferred by Gaia Collaboration (2021a) with the centre of the bar.

GALFIT uses the minimum χ² approach to determine the best fit. We give the minimum reduced χ_R² = χ²/N_px value in the bottom of Table 2. Since the variance of the χ² distribution would be 2N_px, the reduced χ_R² has a standard deviation of $\sqrt{2}$. We did several fits with varying different additional constraints, discussed below. The results are shown in Fig. 12 with the best-fit centre indicated with white lines.

The result of the general photometric decomposition, where all Sérsic parameters for both components were free during the fit is shown in the top row of Fig. 12. The residual image shows that the disc and bar of the LMC are modelled reasonably well. The most noticeable components in the residuals are the LMC spiral arms and dust lanes, which we did not attempt to include in the fit. The on sky projected orientation and flattening of the LMC disc according to this photometric fit is close to what we obtain from the Jeans 2D kinematic fit with |Δψ_d| = 18.7° and |Δq_d| = 0.06. And similarly to the 2D Jeans kinematic fit, it does not agree with the direction of the LOS velocity gradient. The V_LOS rotation axis (as obtained from the 3D Jeans kinematic solution) appears to be inclined by 38.5° with respect to the photometric major axis. The LOS velocity map is also shown in the top row of Fig. 12 in the same orientation.

In the second GALFIT solution we decided to fix the orientation angle ψ_d = 124° to correspond to the orientation angle of the projected rotation axis, according to our best fit 3D Jeans dynamical model. All other GALFIT parameters were kept free as before. The projected flattening of the LMC disc that we find from this photometric decomposition is very close to what we find from the 3D kinematic model (|Δq_d| = 0.02). However, as one would expect, this photometric fit is significantly worse than the first one, having a minimum χ² 6.5 standard deviations higher than the original fit. This is a compromise that we have to make in order to find an axisymmetric disc orientation, compatible with the galaxy’s line of nodes orientation.

Furthermore, we notice that the best-fit Sérsic index values for both the disc and the bar components are equal or slightly below 0.5. A Sérsic index of 0.5 corresponds to a Gaussian distribution. Usually, galaxies have Sérsic indices between 1 and 4, with values 1 − 2 being typical for discs and bars, and n > 2 more typical for spheroidal components (Kormendy & Kennicutt 2004). But the fact that we find so low n values for the LMC may not be too surprising if the stellar component of the LMC is severely tidally limited by its proximity to the Milky Way. In addition, the surface density in the central region of the LMC remains approximately constant in an extended region, which is described better by a Gaussian distribution rather than an exponent. Hence, the third GALFIT solution that we present fixes n_d = n_b = 0.5, as well as q_d = 0.91, in addition to the fixed ψ_d angle (bottom row of Fig. 12). We found that adopting q_d = 0.91 from the dynamical solution of the 3D Jeans model, provides better constraints on the parameters of our Schwarzschild models, compared to using the photometric fit of q_d = 0.93. The GALFIT fits with free and fixed Sérsic indices are practically identical and within 1σ of the χ² distribution, but we need to impose this additional constraint, because we need to be able to represent the derived LMC surface brightness distribution with a MGE. In this case the MGE consists of two Gaussian components – one for the disc and one for the bar (Table 3). It is scaled to the same total LMC luminosity of 1.31 × 10⁹ L_⊙ used in our Jeans models (Table 1).

5.2. LMC de-projection and viewing angles

Once we defined the disc – bar morphological decomposition of the LMC on the plane of the sky – we could proceed with the 3D de-projection of the galaxy. We essentially determined the intrinsic shapes of the disc and bar as well as the viewing angles of the projected image.

The intrinsic shape of the LMC triaxial bar is determined by two axis ratios (q_b⁰ = c_b⁰/a_b⁰ and p_b⁰ = b_b⁰/a_b⁰), where a_b⁰, b_b⁰, c_b⁰ are the intrinsic sizes of the long, intermediate, and short bar axes, respectively. The intrinsic shape of the axisymmetric disc is determined by a single axis ratio (q_d⁰ = b_d⁰/a_d⁰), where a_d⁰, b_d⁰ are the intrinsic sizes of the long and short disc axes. There are also two unknown viewing angles – the inclination (θ) and the tip angle (ϕ). The latter is only relevant for triaxial ellipsoids. The third angle (ψ) can be fixed such that the projected short axis of the disc is vertical (see Fig. 13 for the definition of the viewing angles).

Fig. 13. Definition of the viewing angles of a triaxial ellipsoid. The LOS direction is along the Z axis.

In the general case, the de-projection of a triaxial ellipsoid is a degenerate problem, but we can break the degeneracy with the aid of the co-existing axisymmetric disc. The problem becomes solvable when we assume that the short axis of the bar is aligned with the short axis of the disc (Tahmasebzadeh et al. 2021).

We used an MCMC sampling approach to infer the intrinsic shapes of the LMC disc and bar from our morphological decomposition of the galaxy described in the previous section. This works as follows: We propose intrinsic shapes of the two morphological components and project them with a set of proposed viewing angles. The likelihood function maximises the similarity of the resulting projected morphology to the modelled image from GALFIT. We aim to match the projected flattening of both the disc and bar as well as the angle between the projected major axes of the two components.

The projection equation is from Binney (1985):

$$\begin{array}{c}\hfill \left(\begin{array}{c}{X}_p\\ Y_p\\ 0\end{array}\right)=\left(\begin{array}{ccc}-sin\varphi & cos\varphi & 0\\ -cos\theta cos\varphi & -cos\theta sin\varphi & sin\varphi \\ sin\theta cos\varphi & sin\theta sin\varphi & cos\theta \end{array}\right)·\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right).\end{array}$$(5)

The results of the de-projection routine are summarised in Fig. 14 for the two LMC decompositions shown in Fig. 12 in corresponding rows.

Fig. 14. LMC de-projection and viewing angles. Left panels: The corner plots show the MCMC distributions of the viable intrinsic axis ratios of the axisymmetric disc (qd0), the triaxial bar (qb0, pb0), and the viewing angles θ and ϕ in rad that result in the observed morphological decompositions of the LMC from the two cases shown in Fig. 12. Middle panels: Zoom-in of the (θ, ϕ) planes colour-coded by the disc intrinsic flattening. The dashed black lines show the median angles from the de-projection, while the red dashed line shows the best-fit inclination angle (θ) from the 2D Jeans dynamical model in the top panel and the 3D Jeans model in the bottom panel. Right panels: Embedded ellipsoids of the de-projected disc and bar with axis ratios as inferred from the de-projection of the LMC morphological decomposition from Fig. 12 and viewed from the three cardinal directions and from viewing angles θ and ϕ.

First we turn our attention to the top row of Fig. 14, corresponding to the general GALFIT solution with a mismatch between the projected long axis of the disc and the line of nodes orientation. The corner plot shows that the viewing angles are well constrained and we find an inclination angle θ = 39.0° ±2.3° and a tip angle ϕ = −17.8° ±0.6°. The inclination is slightly larger than what we inferred from the 2D Jeans dynamical model for the disc (θ = 34° ±0.1°), but in this de-projection routine we used a slightly lower projected axis ratio q_d = 0.78 from GALFIT, compared to q_d = 0.84, as inferred from the 2D Jeans model (see Tables 1 and 2). The lower angle inferred from the dynamical modelling of θ = 34° ±0.1° is still within the allowed range of angles for the de-projection (see the middle panel of the top row in Fig. 14). The intermediate-to-major axis ratio for the bar is also well constrained at p_b⁰ = 0.30 ± 0.02, while we are able to only put an upper limit to the minor-to-major axis ratio for the bar at q_b⁰ ≲ 0.3. The intrinsic flattening of the disc, however, is degenerate in the targeted projection, so we need to adopt the additional prior constraint that q_d⁰ < 0.3. The top right panel of Fig. 14 shows points sampled on the surfaces of two ellipsoids, representing the LMC disc and bar, using the median axis ratios, as inferred from our de-projection routine and how they project on the sky at the median viewing angles – a 2D geometry essentially the same as the one inferred from the GALFIT LMC decomposition shown in Fig. 12.

Now we look at the de-projection results shown in the bottom rows of Fig. 14, corresponding to the GALFIT solution, where the major axis of the disc and the line of nodes orientation are aligned. This is the de-projection that we adopt for our Schwarzschild models described in the next section. The viewing angles that we find for this targeted projection are θ = 25° ±3.0° and ϕ = −61.3° ±1.1°. The inclination angle is in line with the disc inclination that we infer from the 3D Jeans dynamical model (θ = 25.5° ±0.2°; Table 1) within the uncertainties. We note that the targeted projected flattening of the disc in this case is q_d = 0.91 – the same as inferred from the 3D Jeans model. We also manage to get a good constraint of the intrinsic flattening of the LMC disc (q_d⁰ = 0.23 ± 0.05) from this targeted projected geometry without the need for prior assumptions. The best-fit axis ratios for the bar are p_b⁰ = 0.38 ± 0.02 and again an upper limit for q_b⁰ ≲ 0.35. The bottom right panel of Fig. 14 illustrates the inferred 3D geometry of the LMC under our simplifying assumptions and the resulting projected morphology from the fitted median viewing angles, which is in excellent agreement with the GALFIT model shown in the bottom row of Fig. 12.

6. Schwarzschild models and the bar pattern speed

The orbit-superposition technique (Schwarzschild 1979) is a powerful tool to study galaxy dynamics (Rix et al. 1997; Gebhardt et al. 2003; Cappellari et al. 2007), not limited to the assumption of axisymmetry (van den Bosch et al. 2008; Neureiter et al. 2021). It relies on the idea of superposing a large number of individual stellar orbits to simulate the overall mass distribution and kinematics of galaxies. The method starts with an assumed gravitational potential. Then, a collection of orbits is numerically integrated in this potential. Each orbit spans a certain region of phase space and represents many stars. Individual orbits are weighted and superimposed to match the observed surface brightness and kinematic data, amongst other constraints. The method has been further developed to include the stellar age and metallicity, enabling chemo-dynamical decomposition of galaxy structures (Zhu et al. 2022; Ding et al. 2023; Jin et al. 2024).

Here we use a modified version of the Schwarzschild code originally developed by van den Bosch et al. (2008) and later published as DYNAMITE⁵ (Jethwa et al. 2020; Thater et al. 2022) to model the LMC as a two-component stellar system, consisting of an axisymmetric disc and a triaxial bar, following an approach outlined in Tahmasebzadeh et al. (2022). The most important modification is that the orbit integration is done in a rotating frame of reference, in which the barred galaxy potential is stationary. The frame rotation corresponds to the bar pattern speed (Ω) – a free parameter in our Schwarzschild models. As a consequence, the retrograde orbits cannot be derived by reversing the sign of the prograde orbit family and need to be integrated separately. In addition, the 8-fold orbit symmetry in non-rotating reference frames is reduced to 4-fold symmetrisation.

6.1. Data preparation

Currently DYNAMITE takes as input data Voronoi tessellation maps of the velocity field, velocity dispersion, and higher velocity moments with their corresponding uncertainty maps in the LOS direction. The kinematic data set is complemented with a MGE of the projected surface luminosity density distribution and viewing angles. We used the projected luminosity density distribution MGE derived with GALFIT from the Gaia DR3 density map of the LMC and shown in Table 3. The corresponding intrinsic shapes and viewing angles are as defined in Sect. 5.2 (bottom panel of Fig. 14). We allowed for 1% relative uncertainty in the projected and intrinsic luminosity distributions.

Here we focus on the preparation of the LOS kinematic data for the Schwarzschild modelling. The input Voronoi maps are based on the ones shown in the bottom left panels of Figs. 6 and 8 to illustrate the performance of the discrete 3D Jeans model in recovering the observed LOS velocity field and LOS velocity dispersion of the LMC, respectively. To construct the Voronoi maps, we consider only stars with high probability of being genuine LMC members according to our 3D Jeans model. To be on the conservative site, we impose a minimum LOS velocity error of 1 km s⁻¹ for individual stars (if the quoted Gaia DR3 LOS velocity uncertainty is less than 1 km s⁻¹, we replace it with 1 km s⁻¹) and remove stars with Gaia DR3 LOS velocity errors larger than 15 km s⁻¹. Some Gaia DR3 sources have unrealistically low LOS velocity uncertainty entries. The minimum LOS velocity error in our sample is only 0.11 km s⁻¹. In total ∼1900 sources (out of ∼28 000) have Gaia DR3 errors below 1 km s⁻¹. Imposing a minimum LOS velocity uncertainty allows for a bit more flexibility in our Schwarzschild models and prevents them from getting stuck in local χ² minima driven by spurious kinematic features. The resulting median LOS velocity error of the individual stars is 3 km s⁻¹. Each individual Voronoi cell contains about 100 stars with Gaia DR3 LOS velocity measurements and its mean velocity, velocity dispersion, and respective uncertainties are computed using a maximum likelihood approach, taking the individual errors of the stellar velocity measurements into account. We ignore higher velocity moments (h₃, h₄, …) because these limited statistical samples are insufficient to estimate them reliably.

The resulting LOS dispersion map is, however, quite noisy, and initial tests showed that this becomes an issue for the Schwarzschild models when looking for a global χ² minimum. We noticed that there is a significant difference in the velocity dispersion of the old and young LMC populations, with the young stars showing a considerably lower dispersion. We opted to use only the old stars of the LMC, as defined in Fig. 4, as kinematic tracers for both the LOS rotation pattern and velocity dispersion. We also applied a Gaussian smoothing filter with a kernel size of 25 arcmin to the velocity dispersion map of the old stellar population to further decrease the remaining irregularities. Such smoothing was not applied to the rotation field in order to not lose any subtle effects in the map caused by the bar kinematics. The median uncertainty of the Voronoi bins is 2 km s⁻¹ in the velocity field and 1.5 km s⁻¹ in the velocity dispersion field.

Since in these models the bar is centred by construction, we adopted the kinematic centre determined by Gaia Collaboration (2021a) in our input Voronoi maps, as it is closest to the photometric centre of the bar.

6.2. Gravitational potential

The gravitational potential consists of visible (mass-follows-light) and DM (a classical NFW halo) components, defined by three parameters – stellar mass-to-light ratio M_*/L, virial mass-to-stellar mass ratio (M₂₀₀/M_*), and NFW halo concentration (C = r₂₀₀/r_s), where r₂₀₀ is the virial radius and r_s – the NFW halo scale radius. We add a very light black hole in the centre with a mass of 10 M_⊙, which only purpose is to provide some additional stability to the orbit integration in the innermost region of the galaxy, but we expect that its overall effect on the dynamical properties of the modelled galaxy is negligible.

We examined a dense grid of possible potentials, varying log M₂₀₀/M_* and log C in step size of 0.05, and M_*/L in step size of 0.2 M_⊙/L_⊙. In addition the potential is stationary in a rotating reference frame defined by the bar pattern speed Ω. We tested a wide range of Ω ranging from −10 km s⁻¹ kpc⁻¹ (counter-rotating bar, motivated by the possible solutions found by Jiménez-Arranz et al. 2024a), up to 30 km s⁻¹ kpc⁻¹ with a step size of 1 km s⁻¹ kpc⁻¹.

6.3. Orbit library and model setup

Next we setup the orbit library used in our Schwarzschild dynamical models. We use a 30 × 15 × 10 three-dimensional main grid, defined by the three integrals of motion – orbit energy (E), I₂ (similar to the vertical angular momentum), and I₃. The orbits are sampled logarithmically in radius from r_min = 10^1.5 arcsec to r_max = 10^4.5 arcsec (7.6 pc to 7.6 kpc), covering the radial extent of our kinematic tracers. We also define a retrograde orbit grid, using the exact same settings and in addition, we use an orbit dithering parameter (d = 3), which defines a mini-grid around each of the main grid bins for a total of 2 × 30 × 15 × 10 = 9000 orbit bundles (each bundle consists of d³ = 27 orbits) to integrate. We impose a 10⁻⁶ relative energy conservation accuracy and integrate them in the selected gravitational potential for 200 periods. We then sample 50 000 points in the meridional plane for each orbit.

The Schwarzschild dynamical solution is a linear combination of all orbit bundles from which we can compute the resulting projected model surface brightness and model the LOS velocity distribution (LOSVD) in the same Voronoi tessellation as the input data. We find the orbit bundle weights, using a non-linear least squares methodology (NNLS, Lawson & Hanson 1974). The model LOSVD of each spatial bin is represented by a velocity histogram with 200 velocity bins 3 km s⁻¹ wide, that are matched to the input Gaussian LOSVDs (Fig. 15) and a χ² value is computed.

Fig. 15. Example LOSVD match between the input data (red) and the orbit superposition (black) for the 35 central Voronoi bins (out of 243). The input LOSVDs are represented by Gaussians, while the model LOSVDs are histograms with 200 velocity bins 3 km s−1 wide and constructed from all sampled points from the weighted orbit library that fall in the same projected Voronoi bins.

This procedure is repeated for different gravitational potentials and bar pattern speeds.

6.4. Model results

The upper left panel of Fig. 16 shows the full grid of Schwarzschild models that we computed varying the NFW parameters, the stellar M_*/L ratio, and the bar pattern speed. We find a global χ² minimum for Ω = 8 km s⁻¹ kpc⁻¹, log C = 0.9, log M₂₀₀/M_* = 4.0, and M_*/L = 0.3 M_⊙/L_⊙ (with an upper limit of 0.7 M_⊙/L_⊙ at the 3σ level). All model parameters appear to be well constrained by the explored χ² grid. We do expect a certain level of correlation between the mass parameters, given that our kinematic tracers are mostly concentrated in the inner region of the LMC and do not trace the full potential. There are multiple combinations of halo virial mass, concentration, and M/L_* that can result in the similar inner mass density distributions. While the gravitational potential parameters do have a physical meaning, we caution that here we treat them merely as pseudo-parameters that describe the potential shape within the extent of the kinematic tracers, as best as the chosen parametrisation allows. M_*/L = 0.3 M_⊙/L_⊙ is unrealistically low and the corresponding log M₂₀₀/M_* = 4.0 implies a virial mass M₂₀₀ ∼ 4 × 10¹² M_⊙, which is also non-physical. On the other hand, the measurement of the bar pattern speed is tightly constrained by the Schwarzschild model between 6 and 13 km s⁻¹ kpc⁻¹ at the 3σ level.

Fig. 16. Schwarzschild model results. Top-left panel: Full grid of Schwarzschild models (potential parameters and bar pattern speed) that we computed, colour-coded by χ2. The best-fit model is indicated with a cross, and the one-, two-, and three-σ contours are overplotted with red, blue, and cyan, respectively. Top-right panel: Performance of the best Schwarzschild model with Ω = 8 km s−1 kpc−1. Top-row: Model input. This includes surface brightness MGE, LOS velocity, and smoothed LOS velocity dispersion. Middle row: Model output. This contains the same quantities as computed from the weighted orbits superposition. Bottom row: Model residuals. Bottom-left panel: Enclosed cumulative mass profiles of the stellar (M*/L = 0.3 M⊙/L⊙), DM, and total mass distributions according to the best-fit Schwarzschild model. Bottom-middle panel: Circular velocity curve according to the same model together with the frame rotation with Ω = 8 km s−1 kpc−1. The shaded regions in both panels indicate the 1σ confidence intervals. Bottom-right panel: Orbit weights as a function of their circularity and radius for the best-fit model. The vertical line denotes the size of the bar, and the horizontal lines differentiate different types of orbits: λz ∼ 0 are radial orbits, λz ∼ 1 are highly circular orbits, and λz < 0 are retrograde.

If we restrict the explored Schwarzschild model grid to more realistic stellar M_*/L ratios ranging from 1 to 1.5 M_⊙/L_⊙, motivated from predictions of composite stellar population models for LMC type galaxies (van der Marel et al. 2002), the best fit DM halo parameters correspond to M₂₀₀ = 1 − 2 × 10¹² M_⊙, which is still an order of magnitude higher, than independent virial mass estimates for the LMC (1 − 2 × 10¹¹ M_⊙; Vasiliev 2023). The result is a manifestation of the same extrapolation issue that we discussed with our Jeans models. However, it is even more pronounced due to the steeper stellar mass distribution adopted here, which was motivated from the Gaia DR3 density image of the LMC (Fig. 12), versus the exponential profile adopted in the Jeans analysis.

The enclosed mass profile and the corresponding circular velocity curve of the LMC are shown in the bottom left and middle panels of Fig. 16. The circular velocity plots are also complemented with a curve showing the best fit rate of the bar frame rotation in each model grid. The shaded regions correspond to the 1σ confidence intervals.

In Fig. 17 we plot together the cumulative total mass profiles for select Schwarzschild models from our model grid and the best fit Jeans 3D models (DM only and with mixed gravitational potential), as well as their corresponding circular velocity curves. The selected Schwarzschild models include the one with the global χ² minimum (Fig. 16) and four additional ones with higher χ² values but more realistic LMC physical parameters. The first pair includes the best fit models at fixed stellar M_*/L = 1 and 1.5 M_⊙/L_⊙, corresponding to virial masses M₂₀₀ = 1 and 2 × 10¹² M_⊙, respectively. To select the second pair we additionally constrain the grid to only models with a fixed log(M₂₀₀/M_*) = 1.9, corresponding to a realistic LMC virial mass for the above stellar M_*/L ratios. The bar pattern speed In all five models is constrained between 8 and 15 km s⁻¹ kpc⁻¹. The estimated enclosed mass profiles within the radial extent of our kinematic tracers (6.2 kpc) agree within the uncertainties. The tight constraints on the mass profile that all models impose assures us that we have properly explored the entire parameter space and that there is not further room to expand the parameter grids for the Schwarzschild modelling, unless we change the potential parametrisation entirely. One could see that the over all best fit Schwarzschild model follows the 3D Jeans models most closely up to the radial extent of the kinematic tracers, but it does not have an extrapolative power at larger distances, as it becomes extremely steep. While all models agree well within a radial distance of ∼3 kpc, we observe a significant divergence going outwards in the models with a realistic virial mass (1 − 2 × 10¹¹ M_⊙), especially noticeable in the circular velocity plots, which can explain their worse χ² results.

Fig. 17. Dynamical models comparison. Left panel: Enclosed cumulative total mass profiles for select Schwarzschild dynamical models discussed in this work compared to the 3D Jeans models. Right panel: Corresponding circular velocity curves. We overplot the figure rotation rate for a bar pattern speed of 11 km s−1 kpc−1 (with a shaded region from 8 to 15 km s−1 kpc−1). The vertical lines correspond to the size of the LMC bar and the extent of our kinematic data. Everything beyond the second vertical line is an extrapolation of the model.

In the circular velocity plot we have also added the frame rotation rate for a bar pattern speed of ${11}_{-3}^{+4}$ km s⁻¹ kpc⁻¹, in line with the 1σ confidence intervals of the five Schwarzschild models combined. The point where the circular velocity curve and the frame rotation curve cross is the co-rotation radius, which according to the different dynamical models appears at roughly R_cor ∼ 10 kpc, depending on which rotation curve yields the most realistic results that far out in the LMC disc.

The resulting kinematic maps based on this best performing model are shown in the top right panel of Fig. 16. The best fit model captures well the surface brightness and morphology of both the LMC disc and bar, which we modelled with Gaussian surface brightness distributions (Table 3). The triaxial bar model is also successful in capturing the complexity of the LOS velocity field in the central region. There is a hint of a decoupled bar rotation pattern in the input data, that creates asymmetries and twists in the LOS velocity field in the central region, which are also present in the dynamical model with Ω ∼ 10 km s⁻¹ kpc⁻¹. All tested models, however, have trouble reproducing properly the input velocity dispersion in the inner region, predicting significantly higher values in the centre and a stronger gradient than the data. This could be due to the complexity of the stellar populations in the LMC (we already noted that the young stars have significantly lower velocity dispersion than the old stars, but the latter are also not a uniform population, spanning a wide age range), deviations of the real surface brightness distribution in the central region from our simple MGE input, or velocity bias in the LOS measurements of the Gaia DR3 catalogue due to crowding. Still qualitatively, we can find similarities between the input velocity dispersion field and the Schwarzschild model result. The model forecasts dips in dispersion at both extremities of the bar. These dips can also be discerned in our input velocity dispersion map, although they are somewhat less pronounced, likely due to the implemented smoothing filter. They also appear slightly displaced from their predicted locations according to the model, likely a consequence of the actual LMC bar being off-centred. Going outwards, the disc’s velocity dispersion aligns well with the predictions made by the Schwarzschild model.

Figure 18 shows how the projected LOS velocity and velocity dispersion fields change by varying Ω, when all other model parameters are kept fixed (M₂₀₀/M_* = 2.9, log C = 0.7, and M_*/L = 1.0 M_⊙/L_⊙). We note that the orbit weights are re-fitted for each case to the input data. By observing these maps, one can more easily understand the behaviour of the χ² distribution with respect to the bar pattern speed. Models with Ω = 0 to Ω ∼ 12 km s⁻¹ kpc⁻¹ are qualitatively most similar to the observed V_LOS field. However the telltale twist in the rotation axis that we observed in the data starts to become noticeable in the models only for Ω > 6 km s⁻¹ kpc⁻¹. Similarly, the model velocity dispersion maps are most similar to the input velocity dispersion for Ω in the range 9 − 15 km s⁻¹ kpc⁻¹. By increasing Ω above 15 km s⁻¹ kpc⁻¹, a strong kinematic effect from the bar develops, which shows a distinctive counter-rotating pattern, not visible in our input data.

Fig. 18. Projected LOS velocity and velocity dispersion maps for the best-fit Schwarzschild potential with log C = 0.7, log M200/M* = 2.9, and M*/L = 1.0 M⊙/L⊙ but varying Ω from 0 to 27 km s−1 kpc−1.

In the bottom right panel of Fig. 16 we show the orbit weights distribution for the best fit Schwarzschild model. We characterise each orbit by its circularity (λ_z ≡ J_z/J_max(E), where J_z is the angular momentum of the orbit along the z-axis and J_max(E) is the maximum angular momentum of a circular orbit with the same binding energy E) and radius (r) and bin them in a grid, where the weight of each bin is the sum of the weights of all orbits that fall into that bin. Not surprisingly, the disc dominated outer region of the LMC can be described with highly circular orbits (λ_z ∼ 1), while the bar motions require a wide range of different orbit families, including radial (λ_z ∼ 0) and retrograde orbits (λ_z < 0). One may notice that the modelled LOSVD in some central Voronoi bins (Fig. 15) is double peaked. This is owing to contribution from retrograde orbits. We did test model fits, where we gradually removed the most extreme retrograde orbits (λ_z < −0.75) from our library and confirmed that the double peaks disappear. However, we were not able to obtain a better χ² solution than using the full orbit library. By cutting retrograde orbits, even if the few affected bins are better fitted, the overall central velocity dispersion of the rest of the modelled bins increases, making them less compatible with the observational data.

Figure 19 shows the intrinsic axis ratios p₀, q₀, and a triaxiality parameter T = (1 − p₀²)/(1 − q₀²) as a function of radius as recovered from the weighted orbit superposition. We use the variations of these parameters to estimate the size of the LMC bar. As expected, the q₀ ratio is relatively flat around 0.25 in accordance with our input for the flattening of the LMC disc. The p₀ ratio is settled around 0.4 in the central region dominated by the bar and then rapidly increases to one when the bar relevance fades and the axisymmetric disc starts to dominate. Similarly the T parameter is around 0.9 (high triaxiality) in the bar dominated region and then falls to zero going outwards. We considered the crossing of these two curves at ∼2.2 kpc to indicate the intrinsic size of the LMC bar (R_bar). This size also corresponds to 2σ of the Gaussian surface brightness density profile that represents the bar in our MGE model for the LMC and is similar to what one would measure by eye from the Gaia DR3 density map of the LMC. While the surface brightness MGE and the Gaia DR3 density map are projected images of the LMC, the projected bar size is likely not too far from its intrinsic size due to the low inclination of the LMC that we adopted in our model.

Fig. 19. Triaxiality profile of the LMC as reconstructed from the orbit superposition. The vertical lines in all panels indicate the size of the LMC bar and the radial extent of our kinematic tracers.

Considering Ω ∼ 11 ± 4 km s⁻¹ kpc⁻¹ and a co-rotation radius R_cor ∼ 10 ± 4 kpc, we find a dimensionless bar rotation parameter R_cor/R_bar ≃ 4.5 ± 1.8. The main uncertainty here is the circular velocity of the LMC, since the co-rotation radius appears slightly outside of the extent of our kinematic tracers and the dynamical solution becomes entirely dependent on the parametrisation of the LMC gravitational potential (see Fig. 17). Galaxy bars can be classified as fast if 1 < R_cor/R_bar < 1.4 or slow if R_cor/R_bar > 1.4 (Debattista & Sellwood 2000). Hence, we argue that the LMC hosts a slow bar in agreement with the recent studies (Jiménez-Arranz et al. 2024a,b), based on independent methods. In addition, galaxy bars that are heavily influenced by tidal interactions tend to be slower than secular bars with a dimensionless bar rotation parameter in the two to three range (Gajda et al. 2017) and possibly larger closer to pericentric passages (Gajda et al. 2018), which fits well with our measurement for the LMC bar pattern speed. In such cases, the bar length and rotational velocity can change on a short timescale. Therefore, we can only report instantaneous (present-day) values that might be very different from genuine ones. The significant divergence in the LMC bar parameters existing in the literature (Jiménez-Arranz et al. 2024a) could be understood, as different techniques based on different tracers may capture different stages of the bar parameter fluctuations.

7. Summary and conclusions

In this work, we have presented a series of equilibrium dynamical models for the LMC utilising 3D kinematic data for ∼3 × 10⁴ genuine member stars from Gaia DR3. While the LMC is significantly out of equilibrium as a result of complex tidal interactions with the Milky Way and a recent collision with the SMC (e.g. Besla et al. 2012), which has undoubtedly disrupted the LMC’s kinematic and structural properties, equilibrium dynamical models still remain a valuable tool for analysing the galaxy’s kinematic properties. By applying these models, we could quantify the degree of the LMC’s deviation from equilibrium, offering crucial insights into its current dynamical state.

We utilised discrete 3D Jeans dynamical models to establish the basic dynamical properties of the LMC, its enclosed mass profile, and in particular its disc shape and orientation, under the assumption of axisymmetry. We also used these models to separate high probability LMC member stars from the Milky Way foreground contaminants and unresolved background galaxies based on their 3D kinematic properties and spatial distribution. We have shown that our Jeans model based Bayesian membership selection is as good as the selection method outlined in Gaia Collaboration (2021a), based on PM and parallax.

The fiducial Jeans model uses all member stars with 3D velocities from Gaia within the VMC footprint as kinematic tracers and assumes a spherical NFW gravitational potential. In addition, we computed models where only the kinematics of the old or young LMC stellar populations are used as tracers and found similar results, although we note that the young stellar population is significantly clumpier and has a noisier kinematic distribution than the old stars, which have smoother density and velocity profiles. We found an enclosed mass within 6.2 kpc of ∼1.4 × 10¹⁰ M_⊙ and a virial mass of M₂₀₀ ≃ 1.8 × 10¹¹ M_⊙ based on our fiducial model. The latter result is in excellent agreement with independent estimates of the total LMC mass (Vasiliev 2023), although we take it with caution, as we do not want to imply that a simple NFW potential is sufficient to describe the dynamical complexity of the LMC at all radii. We also tested a mixed gravitational potential that adds the scaled surface brightness profile by a fixed stellar M_*/L = 1.5 M_⊙/L_⊙ to the NFW halo. This model has the same performance as our fiducial model within the extent of the kinematic tracers, although the fitted NFW parameters define a more extended DM halo that extrapolates to an unrealistically high virial mass in order to accommodate the added visible mass. If we let the M_*/L be a free parameter in our model fit, it tends to go to zero, thus reducing the model to the fiducial one.

Our 3D Jeans models constrain the orientation and shape of the LMC disc based purely on its kinematic properties and without any prior knowledge of the surface brightness distribution, besides the assumption of a radially exponential decline. We found a projected flattening of 0.91 at an inclination angle of 25.5° and an orientation angle of 124°, in reasonable agreement with previous studies (e.g. van der Marel & Kallivayalil 2014; Gaia Collaboration 2021a). We also acknowledge the discrepancy between the best-fit kinematic axis and the photometric major axis, which we estimate at an orientation angle of ∼86° from the Gaia DR3 density map of the LMC. These discrepancies speak to the deviation of the LMC disc from axisymmetry, which has also been noted in previous studies (van der Marel & Cioni 2001; van der Marel & Kallivayalil 2014).

Despite these caveats, our axisymmetric Jeans models describe the LMC disc velocity field quite well, capturing the projected velocity gradients due to the disc rotation and the velocity dispersion profiles in the three cardinal directions. Additional tests exploring the velocity gradient correlations in the V_LOS and V_Y directions (van de Ven et al. 2006) confirmed the findings of the Jeans model and that, at least purely kinematically, the old LMC disc behaves as if it were axisymmetric. The young stars, on the other hand, show a significant discrepancy from this assumption.

We used the Gaia DR3 density image and the GALFIT software to decompose the LMC into the co-centric axisymmetric disc and triaxial bar components and to de-project them in order to obtain their intrinsic shapes, density distributions, and viewing angles. The de-projection of a triaxial ellipsoid is a degenerate problem; however, in this case, we can solve it if we assume that the short axis of the triaxial bar is aligned with the short axis of the axisymmetric disc (Tahmasebzadeh et al. 2021). We also opted to keep the disc shape and orientation as derived kinematically from the Jeans model. Although we fitted the Gaia density image of the LMC with a mixture of Sérsic density distribution to represent the two morphological components, we found Sérsic indices of n_d = 0.50 and n_b = 0.34 for the disc and bar, respectively. Since a Sérsic distribution with n = 0.5 reduces to the Gaussian distribution, we represented the surface brightness of the LMC as an MGE with two Gaussian components – one representing the disc and the other the bar, scaled to a total luminosity of 1.31 L_⊙ – as input for our Schwarzschild models together with the derived viewing angles. As kinematic input, we provided Voronoi binned maps of the LOS projected rotation and velocity dispersion fields of the LMC’s old stellar population from Gaia DR3.

We fit the Schwarzschild orbit superposition dynamical models to this dataset using a modified version of the DYNAMITE code (Jethwa et al. 2020; Thater et al. 2022), which allows for orbit integration in a rotating reference frame. The rotation rate that best fits the input data for a given gravitational potential corresponds to a stationary reference frame for the bar and thus gave us the bar pattern speed (Ω). Our gravitational potential consisted of the scaled de-projected MGE density distribution plus a spherical NFW DM halo. We explored a large grid of models, sampling the potential using three parameters – stellar mass to light ratio (M_*/L), dark mass to stellar mass ratio (log M₂₀₀/M_*), and DM concentration (log C) – and testing frame rotation rates from −10 to 30 km s⁻¹ kpc⁻¹. The Schwarzschild models recover a very similar mass distribution as our Jeans models and fit the LOSVD of the inner region of the LMC better, particularly its bar. We found a bar pattern speed Ω = 11 ± 4 km s⁻¹ kpc⁻¹, in line with other recent studies that use independent methodologies (Jiménez-Arranz et al. 2024a,b). With a bar size of ∼2.2 kpc and a co-rotation radius of ∼10 kpc, we argue for a slow bar in the LMC, in line with expectations for bars born from tidal interactions (Gajda et al. 2017, 2018).

This is the first study that attempts to include the LMC’s bar as a triaxial component in an equilibrium dynamical model based on the Schwarzschild orbit superposition method and fit directly to the observed data. Despite significant success in deriving the bar pattern speed of the LMC and describing reasonably well the LOSVD in the inner region of the LMC, there are important caveats to the current models. In the near future, we plan to expand our Schwarzschild models to work effectively with 3D velocity data. The incorporation of the full velocity field will advance our understanding of the dynamics of the LMC, introducing a new layer of depth and complexity to existing studies. Furthermore, our models are still hindered by the assumption of an axisymmetric disc, an assumption that may not hold in reality, as it is plausible that the disc of the LMC is not axisymmetric (van der Marel & Cioni 2001). This assumption therefore risks oversimplifying the actual configurations of the stellar system, and we plan to abandon it and present a fully triaxial Schwarzschild model of the LMC in the near future. Another factor to consider is that the system might be entirely out of equilibrium, and models assuming so may not provide a meaningful description of the LMC’s complex stellar kinematics. In light of this, it is worth noting that the LMC bar could be a transient feature born from a recent collision with the SMC (Besla et al. 2012). Furthermore, it is plausible that the bar is off-centred (de Vaucouleurs & Freeman 1972), thus introducing potential asymmetries in the stellar system’s morphology and dynamics that our models fail to account for.

There might also be problems in the Gaia data itself due to the high stellar density in the LMC’s innermost regions. Looking towards the future, upcoming projects such as the 4MOST 1001 Magellanic Fields survey (Cioni et al. 2019) and Gaia DR4 SIF data (Gaia Collaboration 2023c) promise to significantly improve the kinematic information in the centre of the galaxy. The 4MOST project, for instance, is set to deliver significantly more accurate LOS velocities in the crowded bar region of the LMC. Likewise, Gaia DR4 SIF data will feature considerably more accurate proper motions in the bar, thus providing a wealth of improved 2D kinematic information. Importantly, this influx of precise new data has the potential to greatly enhance the fidelity of our dynamical models and their predictions.

Acknowledgments

We thank Alice Zocchi, Sabine Thater, Edward Lilley, and Prashin Jethwa for insightful discussions and help with the Schwarzschild modelling technique using DYNAMITE. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of NASA’s Astrophysics Data System. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

References

All Tables

All Figures

Fig. 1. Data selection. Left panel: Gaia CMD of the selected LMC targets with 3D kinematics (blue). Middle panel: PM vector diagram of the LMC sources with 3D kinematics, colour-coded by the LOS velocity. Right panel: Same as the middle panel but colour-coded by the μRA μDec correlation coefficient. The grey points in all panels show additional LMC sources with good Gaia PM measurements (μerr < 0.1 mas yr−1) but no LOS velocity entries. These are the sources used in our 2D Jeans model.

In the text

Fig. 2. Spatial distribution and quality control. Left panel: Density distribution of the selected LMC targets with 3D kinematics (blue) with the larger sample of stars with good PM measurements plotted under with light grey dots. Middle panel: Similar to the left panel but the stars with 3D kinematics are colour-coded according to the μRA μDec correlation coefficient. Right panel: Stars with 3D kinematics colour-coded according to the RUWE index.

In the text

	Fig. 3. Velocity error distribution as a function of GaiaG magnitude along the LOS and in the plane of the sky. The PMs have been converted to linear velocities assuming an LMC distance of 49.9 kpc. The densest regions are represented only with contours due to the significant overlap between the points.
In the text

Fig. 4. LMC membership and population selection. Left panel: Infrared CMD from the VMC survey of the high-probability LMC member stars according to the 3D Jeans model fit. Right panel: Infrared CMD of the high-probability contaminant sources (foreground Milky Way stars and background galaxies) according to the 3D Jeans model fit. The inlet shows their spatial distribution in the VMC footprint. The dashed line is the same in both panels. It is used to separate the old and young LMC stellar populations.

In the text

Fig. 5. Distribution of the high probability LMC member stars according to our 3D Jeans model (grey dots), rotated at angle ψ = 124°, so that the best-fit projected rotation axis is vertical and the major axis horizontal. The best-fit kinematic LMC centres are indicated with crosses, based on all stars with 3D velocities (blue cross), old stars (red cross), young stars (yellow cross), and the centre based on 2D kinematics from Gaia Collaboration (2021a, cyan cross). The large ellipse indicates the best-fit projected shape of the LMC disc, based on the axisymmetric Jeans model. The small ellipse indicates the shape of the bar, based on a photometric surface density fit (see Sect. 5). The dashed lines indicate the photometric centre, based on the surface density fit. The solid black line shows the bulk PM of the LMC.

In the text

Fig. 6. Voronoi maps showing the 3D Jeans model fit to the LMC rotation field from the Gaia DR3 kinematic data. The first column represents the data in the X, Y, and LOS directions, respectively. Each Voronoi bin is colour-coded according to the median velocity of all stars that fall into it. The second column represents the Jeans model results and the third column the model residuals. PMs are converted into linear velocities. The velocity scale is kept unchanged between all maps for consistency and easier comparison.

In the text

Fig. 7. Projected rotation curve of the LMC in the three cardinal directions. The blue points show the maximum likelihood mean velocities in each of the Voronoi bins defined in Fig. 6. The yellow points show the corresponding mean velocities of the best-fit 3D Jeans model in the same Voronoi bins. The read lines show 150 random draws from the model posterior distribution along the horizontal axis (vertical axis in the case of VX).

In the text

Fig. 8. Voronoi maps showing the 3D Jeans model fit to the LMC velocity dispersion field from the Gaia DR3 kinematic data. The first column represents the data in the X, Y, and LOS directions, respectively. Each Voronoi bin is colour-coded according to the maximum likelihood velocity dispersion of all stars that fall into it, taking into account their individual velocity uncertainties. The second column represents the Jeans model results and the third column the model residuals. PMs are converted into linear velocities. The velocity scale is kept constant amongst all maps for consistency and easier comparison.

In the text

Fig. 9. Projected velocity dispersion profile of the LMC in the three cardinal directions. The blue points show the maximum likelihood velocity dispersion in each of the Voronoi bins defined in Fig. 8. The yellow points show the corresponding velocity dispersions of the best-fit 3D Jeans model in the same Voronoi bins. The read lines show 150 random draws from the model posterior distribution along the horizontal axis.

In the text

Fig. 10. Cumulative mass (left panel), mass-to-light ratio (middle panel), and rotation curve (right panel) of the LMC according to best-fit NFW only (black) and mixed (magenta) gravitational potentials of the axisymmetric Jeans model. The vertical dashed lines denote the extent of the kinematic tracers and the viral radius (R200) and the horizontal lines denote the enclosed mass and mass-to-light ratio within these radii. In the rotation curve plot, we overlay the total mean velocity of each Voronoi bin defined in Fig. 6 as orange dots, as well as the total velocity of all LMC member stars as a transparent blue cloud.

In the text

Fig. 11. LMC axisymmetry tests. Left panels: Residual variance of the mean VLOS versus the mean VY correlation for different line of nodes orientations, colour-coded by the inclination angle for all stars (top), old stars (middle), and young stars (bottom) rows, respectively. Middle panels: Line of nodes orientation versus inclination angle, colour-coded by the residual variance for the old stars and separated into old and young populations. The shaded regions in these panels correspond to the best fit values from the 3D Jeans model and their uncertainties. Right panels: Mean VLOS versus mean VYtan(θ) for the LMC orientation with minimum residual variance for each respective population. These plots are colour-coded by distance from the centre of each Voronoi bin, and the dashed lines indicate a perfect anti-correlation.

In the text

Fig. 12. LMC disc – bar decomposition. From left to right: Gaia DR3 density map of the LMC, the best-fit GALFIT model, the GALFIT model residuals, and the Gaia DR3 LOS Voronoi map in the same orientation. The white dashed lines denote the best-fit major and minor photometric axes that cross at the best-fit photometric centre according to the GALFIT solution. The top row is the best-fit solution when all GALFIT parameters are left free; the bottom row shows the GALFIT solution when the orientation of the large photometric axis and projected flattening are set to the orientation of the line of nodes and projected flattening according to the 3D Jeans model and using the adopted Sérsic indices of the disc and bar nd = nb = 0.5.

In the text

	Fig. 13. Definition of the viewing angles of a triaxial ellipsoid. The LOS direction is along the Z axis.
In the text

Fig. 14. LMC de-projection and viewing angles. Left panels: The corner plots show the MCMC distributions of the viable intrinsic axis ratios of the axisymmetric disc (qd0), the triaxial bar (qb0, pb0), and the viewing angles θ and ϕ in rad that result in the observed morphological decompositions of the LMC from the two cases shown in Fig. 12. Middle panels: Zoom-in of the (θ, ϕ) planes colour-coded by the disc intrinsic flattening. The dashed black lines show the median angles from the de-projection, while the red dashed line shows the best-fit inclination angle (θ) from the 2D Jeans dynamical model in the top panel and the 3D Jeans model in the bottom panel. Right panels: Embedded ellipsoids of the de-projected disc and bar with axis ratios as inferred from the de-projection of the LMC morphological decomposition from Fig. 12 and viewed from the three cardinal directions and from viewing angles θ and ϕ.

In the text

	Fig. 15. Example LOSVD match between the input data (red) and the orbit superposition (black) for the 35 central Voronoi bins (out of 243). The input LOSVDs are represented by Gaussians, while the model LOSVDs are histograms with 200 velocity bins 3 km s−1 wide and constructed from all sampled points from the weighted orbit library that fall in the same projected Voronoi bins.
In the text

Fig. 16. Schwarzschild model results. Top-left panel: Full grid of Schwarzschild models (potential parameters and bar pattern speed) that we computed, colour-coded by χ2. The best-fit model is indicated with a cross, and the one-, two-, and three-σ contours are overplotted with red, blue, and cyan, respectively. Top-right panel: Performance of the best Schwarzschild model with Ω = 8 km s−1 kpc−1. Top-row: Model input. This includes surface brightness MGE, LOS velocity, and smoothed LOS velocity dispersion. Middle row: Model output. This contains the same quantities as computed from the weighted orbits superposition. Bottom row: Model residuals. Bottom-left panel: Enclosed cumulative mass profiles of the stellar (M*/L = 0.3 M⊙/L⊙), DM, and total mass distributions according to the best-fit Schwarzschild model. Bottom-middle panel: Circular velocity curve according to the same model together with the frame rotation with Ω = 8 km s−1 kpc−1. The shaded regions in both panels indicate the 1σ confidence intervals. Bottom-right panel: Orbit weights as a function of their circularity and radius for the best-fit model. The vertical line denotes the size of the bar, and the horizontal lines differentiate different types of orbits: λz ∼ 0 are radial orbits, λz ∼ 1 are highly circular orbits, and λz < 0 are retrograde.

In the text

Fig. 17. Dynamical models comparison. Left panel: Enclosed cumulative total mass profiles for select Schwarzschild dynamical models discussed in this work compared to the 3D Jeans models. Right panel: Corresponding circular velocity curves. We overplot the figure rotation rate for a bar pattern speed of 11 km s−1 kpc−1 (with a shaded region from 8 to 15 km s−1 kpc−1). The vertical lines correspond to the size of the LMC bar and the extent of our kinematic data. Everything beyond the second vertical line is an extrapolation of the model.

In the text

	Fig. 18. Projected LOS velocity and velocity dispersion maps for the best-fit Schwarzschild potential with log C = 0.7, log M200/M* = 2.9, and M*/L = 1.0 M⊙/L⊙ but varying Ω from 0 to 27 km s−1 kpc−1.
In the text

	Fig. 19. Triaxiality profile of the LMC as reconstructed from the orbit superposition. The vertical lines in all panels indicate the size of the LMC bar and the radial extent of our kinematic tracers.
In the text