A Comprehensive Kinematic Model of the Large Magellanic Cloud Disk from Star Clusters and Field Stars using Gaia DR3: Tracing the Disk Characteristics, Rotation, Bar, and Outliers

The internal kinematics of the Large Magellanic Cloud (LMC) disk have been modeled by several studies using different tracers with varying coverage, resulting in a range of parameters. Here, we model the LMC disk using 1705 star clusters and field stars, based on a robust Markov Chain Monte Carlo method, using Gaia DR3 data. The dependency of the model parameters on the age, coverage, and strength of the clusters are also presented. This is the first comprehensive 2D kinematic study using star clusters. Red clump (RC) stars and young main-sequence stars are also modeled for comparison. The clusters and field stars are found to have distinctly different kinematic centers, disk inclination, position angle of the line of nodes, and scale radius. We also note a significant radial variation of the disk parameters. Clusters and young stars are found to have a large residual proper motion and a relatively large velocity dispersion when compared to the RC field population, which could be due to perturbation from the bar and spiral arms. We trace the presence of the large residual proper motion and noncircular motion among clusters likely to be due to the bar and detect a decrease in the scale radius as a result of the possible evolution of the bar. The kinematically deviant clusters point to a spatiotemporal disturbance in the LMC disk, matching with the expected impact factor and time of the recent collision between the LMC and the Small Magellanic Cloud.


INTRODUCTION
The Magellanic Clouds (MCs) consist of two irregular dwarf galaxies observable from the Southern Hemisphere.The Large Magellanic Cloud (LMC), the larger of the pair, is situated at 49.59±0.09kpc away (Pietrzyński et al. 2019, hereafter, P20), while its smaller companion, the Small Magellanic Cloud (SMC), lies at a distance of 62.44±0.47kpc (Graczyk et al. 2020).The close proximity of the MCs to the Milky Way (MW) offers a unique opportunity to study them and gain valuable insights into the mechanisms that drive galaxy mergers and their consequences for galaxy evolution.
The LMC-SMC-MW is an interacting system (Putman et al. 1998;Diaz & Bekki 2012;Hammer et al. 2015).Traditionally, the kinematics of the LMC was thought to be modified by its interactions with the MW.In the last couple of decades, it was believed that the LMC was experiencing its first encounter with the MW, as indicated in the studies by Kallivayalil et al. (2006a,b), and Besla et al. (2007).However, the recent study by Vasiliev (2023) indicates the potential scenario where the MCs are in their second passage around the MW.Either way, the current morphology of the LMC is predominantly influenced by the interactions with the neighboring SMC rather than the MW (Besla et al. 2012, hereafter B12;Yozin & Bekki 2014;Choi et al. 2018a,b).
There have been significant efforts to understand the internal kinematics of the LMC.They were previously studied primarily using the H I gas and stars in the LMC (Kim et al. 1998;van der Marel et al. 2002;Olsen & Massey 2007; Table 1.The data sets used in the kinematic model of the LMC are summarized here.The length (L d ) and the number (N d ) of data sets are provided as well.Clusters and associated field regions are labeled based on the cluster richness in the 3 rd row, such as C10 representing clusters with 10 or more stars, and so on.A total of 48 data sets are used in our study to perform the kinematic modeling.2. We note that the average number of members in each cluster in our study is ∼ 40.To investigate potential biases in the data sets, we conducted modeling of groups based on the cluster members' strength across a range from 10 to 60 members, with intervals of 10.The field regions associated with these richness groups were also modeled to compare with the nature of clusters.We created 12 models for the cluster groups and field groups based on cluster richness.
The associated nearby field regions were also modeled with respect to each age group.It is to be noted that the age-wise grouping is only valid for the clusters, and the associated field regions are age-wise heterogeneous but analyzed for comparison.
4. We investigated the variations in kinematic parameters concerning the spatial extent of the LMC.Employing circular regions, we analyzed clusters and associated field regions, ranging from 2 • to 7 • in increments of 0.5 • from the kinematic center on the sky plane.We used the kinematic center that we calculated from the primary model encompassing clusters and field stars for the radially separated groups.
Table 1 summarizes the data sets that are modeled.The subsequent section covers the implementation of the kinematic model using the above data sets.

KINEMATIC MODEL OF THE LMC
We performed the kinematic model of the LMC corresponding to the observed median PM for the data sets referenced in Section 2. The following subsections provide a detailed explanation of the theoretical background we adopted in our modeling and the Bayesian methodology employed to estimate the optimal kinematic parameters for the data sets.

Analytical background of the modeling
The model PM for the candidates in the data sets is formulated based on the equations outlined by van der Marel et al. (2002, hereafter, V02).These equations define the directional vectors of local PM in the West and North (M W , M N ) in the sky plane using a series of 12 independent model parameters.The parameters selected for our model encompass the inclination of the LMC disc (i), the position angle of the line of nodes measured from West (θ), dynamic centers (α 0 , δ 0 ), the amplitude of the tangential velocity of the LMC's center of mass (v t ), the tangential angle made by v t (θ t ), scale radius (R 0 ), optimization factor (η), and the amplitude of the rotational velocity (v 0 ).Our modeling process is aimed at determining the most fitting values for these nine kinematic parameters.It's important to note that our modeling was confined to the LMC's sky plane due to the absence of line-of-sight velocity components in our data sets.
The velocity of sources at any given point (α, δ) within the LMC sky plane is primarily composed of three components.It consists of the velocity of the center of mass (COM) of the LMC, the velocity component associated with the precession and nutation of the galaxy, and the internal rotation component of the galaxy.We assumed that there is no precession and nutation of the LMC disc within the spatial coverage of the galaxy (less than 7 • ) considered in this study (V02, N22).We assumed the RA and DEC of the field regions associated with each cluster to have the same spatial coordinates as the cluster centers.Also, we assumed a fixed distance of 49.59 kpc (P20) to the LMC center, D 0 and a line-of-sight velocity of the COM, v sys , as 262.2 km s −1 (V02).The parametric form for the rotational velocity in the disc plane was adopted from similar kinematic models used in previous studies of the LMC (G21, N22).
The following subsection deals with the method and estimation of the best-fitting kinematic parameters for the data sets summarized in Table 1.

Modeling procedure
We found the best-fitting 48 models and the associated kinematic parameters for 48 data sets referenced in Table 1.We performed the fitting of the parameters using the MCMC serial stretch move sampling algorithm introduced by Goodman & Weare 2010.The code is developed and implemented in C language.A similar MCMC approach is used in the studies by W20, C22 and N22.The observed PM for the sources in data sets are along the RA and DEC directions, but the model formulation is for the PM in the West and North directions.For that, we took the negative of the source median PM in the RA direction.Now the model PM (M W,m , M N,m ) and the observed PM (M W,o , M N,o ) can be used to construct the log-likelihood function (ln L), which can be used with the associated West and North direction standard errors of observed data sets (σ W,o ,σ N,o ) to sample the best fitting parameter with MCMC.The equation for ln L is given by, ln L = −0.5 The priors for the model were uniformly chosen, except for the rotation velocity amplitude, v 0 , for which we used a Gaussian prior of 76 km s −1 (van der Marel & Kallivayalil 2013; G21; N22).When working with modeling datasets based on cluster age groups, we note that the spatial distribution of clusters in the sky plane is not homogeneous.This results in the non-convergence of η in its posterior sampling distributions while performing the MCMC with a uniform prior.To avoid this, we employed Gaussian weighting on η specifically for these age-based datasets, using the corresponding value of the parameter estimated from the primary model encompassing all cluster and field regions.
We executed 2000 steps of MCMC iteration involving 200 walkers evolving sequentially at each step.From the sampled posterior values for the nine parameters, we focused on the final 50%, calculating their median alongside the 16 th and 84 th percentile errors for estimation.In the following section, we present the comprehensive results obtained from the above modeling procedures.

RESULTS
This is the first 2D kinematic model of the LMC employing clusters and neighboring field regions with Gaia DR3 data.The data coverage of the LMC considered in this study is within ∼ 7 • from the LMC center.The primary modeling is performed with 1705 clusters and nearby field regions, estimating nine best-fitting model parameters defining the kinematic properties of the LMC.Figures 1 and 2 show the sampled posterior distribution of the parameters for clusters and nearby field regions, respectively.We also performed the parameter estimation for the rest of the data sets mentioned in Table 1.In the following subsections, we provide the estimated kinematic parameters obtained for all the modeled data sets.Additionally, we present the residual PM derived from comparing the model and observed data.

Cluster and Field kinematics
Table 2 provides the estimated kinematic parameters for the primary model involving clusters and field stars.The position angle of the line of nodes is conventionally measured from the North of the galaxy, whereas the modeling we performed involved measuring from the West.Therefore, we define the position angle of the line of nodes measured from North (PA) as Θ, which is (θ-90 • ) throughout the sections.The COM PM values in the West and the North directions (µ W,com and µ N,com ) can be obtained using the estimated tangential velocity vector estimated in our models (using v t and θ t and Equation 10 in V02).The orthographic projections of the sky coordinate and PM vectors are performed using equations 1, 2, and 3, as mentioned in G21. Figure 3 shows the bulk motion and rotation of the LMC obtained with clusters and field by subtracting the COM PM values estimated from the model.Below we compare the notable kinematic properties of clusters and field regions obtained from the modeling.
Similarly, for the field regions, the values are 1.851 ± 0.002 mas yr −1 and 0.344 ± 0.004 mas yr −1 .These values suggest no significant difference between the observed COM PM between cluster and field.• We note an offset of 28±8 arcmin between the dynamic centers (α 0 , δ 0 ) of the cluster and field.
• The modeled rotational parameters, R 0 and η 0 appear to be larger for field regions compared to clusters, while the v 0 remains almost similar.
Two major notable differences in the kinematic properties obtained are, (1) for the value of R 0 , where clusters have a relatively low value when compared to the field population, and (2) a significant offset in the dynamic centers of cluster and field population.In the next sub-section, we look into the control population we used in our model to understand the kinematic nature of the cluster and field population.

Comparison with control population
Table 3 provides the estimated kinematic parameters for the control population as mentioned in Section 2. Below we compare the cluster and field with the young MS stars and RC stars to understand their kinematic nature.
• µ W,com and µ N,com appear to be similar in all the data sets of the control population without any significant difference.
• The value of α 0 for the field and RC population are the same within errors, whereas the value of δ 0 shows a difference.The α 0 for clusters shows a significant offset in comparison to other populations.In general, the younger population tends to show a southern and eastward dynamic center.
• We note a significant shift in R 0 , |∆R 0 | ≈ 1.7 kpc between young MS and RC population.A corresponding change of η is observed as well.We also estimate a larger value for v 0 = 90.92+1.26 −1.18 km s −1 for the young MS population in comparison with the other data sets.As observed in the previous comparison (subsection 4.1), the value of Θ appears to be changing across the control population.Notably, the clusters and young MS population show a relatively small value in comparison to the field and RC population.Specifically, the R 0 for the cluster and young MS population falls within the range of 2 kpc, while the field and RC population have values exceeding 2 kpc.
The model parameters estimated for the cluster and the control samples are used to obtain the values of rotational velocity (V rot ) of the LMC. Figure 4 shows the spatial distribution of V rot among the clusters, field stars, Young stars, and RC population.Though the overall appearance of the plots is similar, we note that clusters (panel a) and the young MS stars (panel c) appear to have similar spatial distributions of V rot , whereas the field (panel b) and the RC stars (panel d) appear to have similar V rot distribution.We also note that the elongated bar feature is more prominent in the V rot maps of clusters and young MS stars.This suggests that the kinematic model of the LMC depends on the age of the population.In the subsequent subsection, we look into the kinematic model of the LMC with data sets of different ages, using the cluster age groups and nearby field regions as mentioned in Section 2.

Age dependent kinematics of the LMC
The age-dependent variation of the kinematic parameters of the LMC was estimated using the data sets mentioned in Section 2. Table 4 provides the estimated kinematic parameters for the cluster and nearby field data sets based on the cluster age groups.Below, we compare their kinematic properties.
• As noted in the previous subsections, the COM PM shows minimal variations across different age groups.Even slight variations in the parameter v t estimated as in Table 4 for different age groups do not result in substantial shifts in µ W,com and µ N,com .
• The clusters ranging in age from ∼ 100 Myr to 1.25 Gyr show a small offset in the (α 0 , δ 0 ) to the South-East with respect to the older clusters, as evident from C AG−2 and C AG−3 .Meanwhile, the nearby field population in these age groups shows variation only in the South.
• Θ attains its minimum value of 111 • .11+2.08 −2.04 for the C AG−3 , with field also showing a similar trend.Notably, the i (35 • .36+1.2  −1.31 ) for C AG−3 is the largest when compared to other age groups.In the case of the field regions, we do not detect any significant shift in i, but a smaller shift in Θ is noted.
• In the cluster data sets, the minimum value of R 0 is estimated to be 1.34 +0.17 −0.16 kpc for the C AG−3 , while its maximum value occurs at 3.14 +0.18  −0.22 kpc for the C AG−1 .Similarly, in the field data sets, there is a comparable pattern with the minimum and maximum values of 2.28 +0.06 −0.07 kpc and 2.84 +0.08 −0.08 kpc for F AG−3 and F AG−1 , respectively.η also shows a similar trend in the age-based data sets.
• The values of v 0 remain relatively consistent across all datasets, and any observed changes are not significant, as they are within the margin of errors.
As noted in the above comparison, there is a noticeable deviation of kinematic parameters for the cluster age groups, in contrast to the datasets of the field region associated with each age group.The comparison presented in this section has shown that the dependence of kinematic parameters on the age of the population is indeed present, and the variations are statistically significant.In the following subsections, we further investigate the dependence of estimated parameters on the cluster richness as mentioned in Section 2.

Influence of cluster richness on the kinematic model
Table 5 shows the estimated kinematic parameters for the data sets based on cluster richness, as outlined in Section 2. The cluster data sets and the nearby field regions are modeled based on cluster richness, which helped in checking for any kinematic changes in the model based on poor and rich clusters.
Figure 5 shows the variation of the kinematic parameters corresponding to different cluster and field groups based on cluster richness.We note an offset of dynamic centers from poor to rich clusters in the increasing RA directions for both clusters and field populations.Meanwhile, the shift of COM PM for poor to rich clusters is seen increasing in North and West directions.Similarly, estimates for Φ, Θ, R 0 , and v 0 show a decrease from C/F10 to C/F60 groups (as seen in Table 5).However, a notable deviation becomes apparent in the kinematic parameters when considering entire cluster groups compared to field groups.This once again indicates that the kinematic model for clusters and the field population indeed manifests distinct kinematic properties.We also note that the cluster richness groups do not have similar spatial coverage.Overall, we do detect variation in the COM PM, dynamic center, i, and R 0 , which may partially be due to the spatial coverage, and the variation in Θ is statistically insignificant.

Impact of the LMC's spatial coverage on kinematic properties
Table 6 shows the estimated kinematic parameters for the data sets based on the spatial coverage of the LMC.The circular regions of radii increasing from 2 • to 7 • with a step size of 0.5 • from the kinematic center (from the primary model, see Table 2) of the LMC were used to perform the modeling.
Figure 6 shows the radial variation of the kinematic parameters in the case of clusters and nearby field regions.Within the radii considered in this study, there is a notable offset in the amplitude of PM, of the order of ∼ 0.06 mas yr −1 in (µ W,com , µ N,com ) between the inner and outer cluster/field regions.Similarly, we note a positional offset of ∼ 0.4 • for clusters and ∼ 0.6 • for the field between the estimated values of α 0 and δ 0 .The value of i steeply decreases from the inner to outer radii, from ∼ 50 • to ∼ 32 • in the case of clusters, and ∼ 43 • to ∼ 31 • in the case of field population.Meanwhile, the value of Θ increases from inner to outer radii, from ∼ 108 • to ∼ 122 • in the case of clusters, and ∼ 110 • to ∼ 128 • in the case of field population.The rotational parameters, R 0 and v 0 after 5 • show convergence to the values estimated as in the primary model (see Table 2).Therefore, the estimated parameters show a significant radial dependence for both clusters and field populations.

Residual PM of the LMC: Cluster vs Field
The residual PM value for the clusters and nearby field populations is found by subtracting the net modeled PM from the net observed PM values.The residual PM vectors for clusters and field regions are plotted in the X-Y plane as shown in Figure 7, panels (a) and (b), respectively.Clusters show larger residual PM in the spatial plot when compared with the field population.The spatial residual plot for clusters shows relatively large residuals in the bar region and in the northern LMC, whereas such large residuals are not found in the corresponding plot for the field population.The residual PM amplitudes (|Residual PM|) of both the clusters and the nearby field are used to generate the probability distribution plot as in Figure 7, panel (c) depicting the distribution of their values.We obtained an RMS value of 0.146±0.002mas yr −1 for clusters and 0.069±0.001mas yr −1 for the nearby field population.The RMS distribution of clusters shows a broader profile, while the field shows a narrower one.

DISCUSSION
In this study, we performed the kinematic modeling of the LMC using 1705 star clusters and nearby field regions for the first time.We also created an extensive additional 46 models using two control samples (young MS stars and RC stars), different cluster age groups and nearby fields, cluster richness groups, and samples with varying spatial coverage of the LMC.The clusters were cleaned using field star contamination (D24) and found less than 2% of MW source contamination in each cluster while cross-matching with the MW/ LMC source catalog by J23.The median  Table 6.The kinematic best-fitting parameters obtained for the LMC based on the spatial coverage of the LMC with clusters and nearby field regions are tabulated below.

Data
value of the PM is used for the study for both clusters and fields and hence is unlikely to be affected by the foreground MW contamination.
In the following subsections, we discuss the important results we obtained from the modeling.We compare the estimated kinematic parameters with the previous studies, analyze the rotation curve of the LMC for the control samples, trace the spatial variation of residual PM across the LMC, analyze the kinematic outliers from the model, and comment on the kinematic variation in parameters based on cluster age.

Comparison of kinematic properties with previous studies
The estimated kinematic parameters of the LMC are compared with those from previous studies, mainly in which either the authors estimated the kinematic centers and COM PM values or adopted them.First we compare the estimations of (α 0 , δ 0 ), (µ W,com , µ N,com ), i, and Θ from studies by N22, C22, G21, W20, and van der Marel & Kallivayalil (2014, hereafter V14).The parameters estimated are tabulated in Table 7 along with estimations from this study for comparison.
Figure 8 shows the distribution of (α 0 , δ 0 ) and (µ W,com , µ N,com ) estimations from the literature studies as mentioned above, along with the estimations from this study.We can notice a significant variation in the values of (α 0 , δ 0 ) from various studies.Our estimations (shown in red) are located well within the range of other estimations, except the α 0 for the cluster population.Most clusters (∼ 75%) used for the modeling in this study are of ages younger than 1 Gyr (D24).Due to this, the estimated values of (α 0 , δ 0 ) for clusters are closer to the estimation by W20 using the young stellar population.Meanwhile, the estimation from the young MS population is much closer to the W20 estimations.The field population and RC population are closely in agreement with the recent studies by G21, V14, N22, and C22.The offset of (α 0 , δ 0 ) in the older and younger stellar population was previously noted by W20.
The µ W,com and µ N,com values for the clusters and field population also align with the recent studies by N22, C22, G21, and W20 in a similar manner.Notably, the clusters show a relatively larger PM in the North direction when compared to the field.However, the kinematic centers and COM PM are distinct for the clusters and field population.The values of i and Θ estimated in this study are compared with the previous studies (see Table 7) and are more or less in agreement.However, the kinematic parameters such as i and Θ are known to vary with respect to the coverage of the LMC (N22).The variations, as seen from the table, also contribute to the coverage.In this study, for the first time, we studied the variation of kinematic parameters as a function of the spatial coverage of the LMC using star clusters, as described in the subsection 4.5.The inner regions of the LMC show a larger inclination with respect to outer regions, which is already noted in the study by Saroon & Subramanian (2022) using RC population.
The following section discusses the rotational parameters (v 0 , R 0 ) and rotation curves of the LMC for various populations.

Rotation of the LMC
The rotation curves for the clusters, nearby fields, young MS stars, and RC stars from the control sample (as mentioned in subsection 4.2) in the LMC plane, after deprojection are shown in Figure 9 (a.Clusters, b.Nearby field, c. young MS stars, d.RC stars).The running average for the data points is made with a bin resolution of 0.25 kpc, and shown in red.The fitted model is shown with the blue curve.The rotational velocity amplitude (v 0 ) does not show a significant variation across the clusters and field population.The estimated v 0 for clusters, fields, and RC stars are consistent with the estimates by W20.The young MS stars in our study show a slightly larger v 0 (∼ 91 km s −1 ) when compared to the RC stars.However, it is consistent with the observed v 0 (∼ 90 km s −1 ) for the younger population studied by N22.
In the cluster age-dependent kinematic model (subsection 4.3), the variation of R 0 for cluster age groups is provided.The age group C AG−3 (log(age) = 8.0 -8.65) and young MS stars from the control population show reduced values of R 0 at 1.34 +0.17 −0.16 kpc and 1.29 +0.08 −0.08 kpc, respectively.This is suggestive of a steeper rise of the rotational velocity when compared to the older population.The decreased R 0 for these age groups may imply a redistribution of mass, leading to a higher mass density in the central regions.Also, in the youngest group C AG−4 (log(age) = 6.55 -8.0), the cluster density in the galaxy shifts towards the North-East regions (D24).In C AG−4 , there is a shift in the value of R 0 to a slightly higher value of 1.98 +0.17 −0.18 kpc.The field population near clusters is dominated by the older population, hence the R 0 value remains high.
Panel (a) of Figure 10 shows all the four model velocity profiles fitted to the control population provided in Figure 9.The v 0 for the control population is almost identical (∼ 81 to 84 kms −1 ), whereas the value for the young MS stars remains slightly higher (∼ 91 kms −1 ).It is clear from this figure that the clusters and young MS stars show a steeper rise of the V rot suggesting an increased mass density in the inner regions.A similar trend of small R 0 in the case of the younger population and larger R 0 in the case of the older population is noted in the study by W20.The variation in the value of R 0 with respect to population may be due to the effect of the evolution of the bar and its activity with time, resulting in the redistribution of mass in the inner regions.
Notably, the dispersion of the rotation velocity profile resides at ∼ 23 km s −1 and ∼ 11 km s −1 up to 6 kpc in the case of clusters and field population, respectively.In the case of the young MS population, there is a larger σ rot of ∼ 22 km s −1 up to 2 kpc, then a slight decline and further rises after 4 kpc.Meanwhile, the RC population has a σ rot of ∼ 20 km s −1 within 1 kpc and slowly declines and levels off at ∼ 12 km s −1 afterward.The larger dispersion values in the central 2 kpc might be due to the non-circular motions due to the bar.
Panel (b) of Figure 10 shows the radial variation of the V rot /σ rot ratio for the four control population.As all of the population have similar V rot (except slightly higher value for the young MS stars), the higher values observed for the field stars suggest a relatively low σ rot .The RC population has a similar profile, though with a slightly lower ratio (suggesting a relatively large σ rot ).The young MS stars have a more or less similar profile till ∼ 3 kpc, and we note a lower value of the ratio beyond this radius (suggesting a relatively large σ rot ).The cluster has a shallow profile and lowest value for the ratio, suggestive of the highest value of σ rot at all radii.Overall, the low value of σ rot the field population points to the minimal disk heating for this population.However, the ratios suggest that young clusters and young stars have relatively large velocity dispersion, which is suggestive of relatively recent perturbation(s) in the LMC disk.This may point to the heating of the gas in the disk (resulting in the formation of stars with similar kinematics), with minimal heating of the stellar disk.The heating may be due to internal perturbations, such as the bar and spiral arms, or external perturbation, which is the interaction with the SMC.
The rotational velocity maps shown in Figure 4 closely match with the velocity maps from G21.Additionally, the variation in slope between evolved and younger populations observed in Figure 10 is consistent with G21.However, there is a discrepancy in the rotation velocity dispersion profiles when compared to G21.

Spatial variation of residual PM: Clusters and Field
In order to trace regions with large variations in the residual PM, the mean |Residual PM| is traced for different radii (1 to 6 • with a bin size of 1 • using annular regions) as shown in Figure 11.We observe a similar trend in radial variation of mean |Residual PM| for both the clusters and field population, though the values are significantly different.For inner radii less than 2 • , the mean residual PM of the LMC is larger for both populations, likely to be due to the presence of the bar.Beyond 2 • , we note a decline in mean residual PM, and it increases after 4 • .This could be due to the presence of spiral arms in the galaxy, where dense star formation is noted.However, there is a clear signature of larger residual PM for clusters when compared to the field population.
The residual PM RMS profiles in Figure 7 can be compared with the study by C22 based on the numerical simulations by B12 with a population older than 1 Gyr. Figure 8 in C22 has provided the predicted logRMS |Residual PM| as a function of the impact parameter (P in kpc) and impact timing (T in Myr) of the recent LMC-SMC interaction.We adopted the value of T as 149 Myr, based on the significant CF peak in the SMC as observed in Figure 8b of D24 to estimate the likely value of P. We note the assumed value of T is also in good agreement with the estimate by Zivick et al. (2018), who suggests a recent direct collision between the LMC and the SMC at 147 Myr ago.The field logRMS |Residual PM| of −1.611 and −1.219 from the primary model and control sample support the Future 40 and Future 60 models provided in C22.We estimate that the impact parameter, P, must reside at less than 3 kpc, and most likely match with the impact parameter of 2 kpc as per B12 simulations.However, the RC stars in the control sample show a larger residual (0.08 mas yr −1 ) compared to the estimates of C22 (0.058 mas yr −1 ).We note that our model involves treating the COM PM and kinematic centers as free parameters, unlike C22.Hence a direct comparison of our results with C22 may be unrealistic.
The clusters and young stars show larger disc heating, with logRMS |Residual PM| of -0.853 and -0.957.This suggests that the young population is more perturbed, unlike the RC stars and other field stars.Also, the clusters, irrespective of age, show larger residual PM.The kinematic deviant clusters are analyzed in the next section.

Kinematically deviant clusters: Tracing the regions of larger residual PM
The residual PM distribution of clusters, as shown in Figure 7a suggests that some clusters have a large residual PM.We traced clusters with significantly large residuals (>3σ of |Residual PM|) and identified the region with the highest      We note several outlier clusters in the observed v rot profile (Figure 9) of the LMC.The clusters falling outside the 3σ margin from the running average of the v rot profile are treated as outliers.Their locations are traced, and their internal PM is plotted in Figure 12b.
The residual PM vectors are plotted in Figure 12c for clusters with the relative angular difference (θ r ) between the observed and modeled internal motion greater than 45 • and 90 • .The red vectors indicate clusters with θ r > 90 • , suggestive of counter-rotation.These are dominantly seen in the North-West of the galaxy.The deviant clusters also trace the non-circular motion of the bar.
The North-West region with the largest residual PM is located between a radial distance of 1-2 • .The same location also shows the presence of counter-rotating clusters.This region has PM deviation both in angle (panel c) and in value (panel a), pointing to the region experiencing an event of perturbation.As this region is located at an impact distance similar to that estimated by C22 and in this study, we speculate that this may be the area that experienced the largest impact due to the recent LMC-SMC collision.We also note that there could also be other reasons, such as perturbations that could occur at the ends of the bar, but a similar disturbance is not seen at the eastern end of the bar.Hence bar perturbation may not be the reason for this disturbance.

Kinematic variation of parameters in the younger clusters
The kinematic model of the LMC for different age groups is studied as in subsection 4.3.We note the kinematic parameters of the LMC, mainly the i, Θ, and R 0 , show variation in the younger age group C AG−3 (100-440 Myr).To further investigate the change in kinematic parameters, we split the C AG−3 into two groups (one from 100 to 200 Myr and the other from 200 to 440 Myr) and performed the kinematic modeling.The variation of i, Θ, and R 0 are shown in Figure 13.
The age group from 100 to 200 Myr shows the largest inclination (∼ 39 • ) compared to other age groups.Meanwhile, we note the smallest Θ value (∼ 105 • ) in this age range when compared to other age ranges.The cluster density shifted towards central regions of the LMC in the age group C AG−3 and towards North-East in C AG−4 (D24).We note this effect in the variation of R 0 to the inner radii of the galaxy, as shown in Figure 13c.The value of R 0 reaches 1.1 kpc for clusters of age between 100 to 200 Myr.A similar effect is seen in the nearby young population, where the value reaches 0.5 kpc.Meanwhile, the field population shows smaller i and lesser variation in R 0 compared to other age groups.The kinematic model of cluster groups based on richness also shows this trend of variation in R 0 .The younger and richer clusters with more than 60 members are located more toward the central regions of the galaxy as well.We note R 0 smaller than 1 kpc at this cluster group.In summary, we note that the i, Θ, and R 0 values for clusters and young MS stars in the age range 100-200 Myr are distinct from the other groups.A possible reason for this deviation is discussed in subsection 5.7.The spatial PM plots of clusters and field population, as shown in Figure 3, do not bring out the bar feature explicitly.The residual PM plots, as shown in Figure 7, on the other hand, show the bar feature, though more prominent in the case of clusters.The radial gradient of the residual PM (Figure 11, panel a) shows a large gradient in the bar region, pointing to residuals resulting from kinematics related to the bar.The residual PM plots shown in Figure 12 (panels a and c) show the possible presence of non-circular motions in the bar region.A significant reduction in the value of R 0 for clusters and young MS stars with respect to RC stars (see Figure 10a) is probably related to the mass distribution within the inner 2 kpc, where the bar is present.The fact that the younger population reaches the maximum value of V rot in a shorter radius when compared to the older population points to a redistribution of mass, most likely due to the evolution of the bar.

Kinematic signatures of the LMC bar
In summary, this study has traced three kinematic signatures of the LMC bar, relatively large residual PM in the inner 2 • radius for both clusters and field population, the presence of non-circular motion among star clusters, and a decrease of R 0 as a result of the possible evolution of the bar.Recently, J23 also found that the dynamics of the inner disk are dominated by the bar.

Effect of the recent LMC-SMC interaction on the LMC disk
The panels (a) and (c) of Figure 12, show that though there are kinematically deviant clusters spread through the LMC disk, there is a specific location where the deviation is maximum.Panel (c) indicates that the same location also shows the presence of cluster motion in all directions, suggesting a possible specific external disturbance.It is therefore important to identify the source of the disturbance identified in this study.
It is possible that this disturbance is caused by the LMC-SMC collision at ∼ 150 Myr ago as the radial distance of this region from the center also matches the value of the impact factor of the LMC-SMC collision as derived by C22.The age-dating of the disturbance, as shown in Figure 13 suggests a significant shift in inclination i, PA, and R 0 for clusters in the age range 100-200 Myr.The kinematic parameters of clusters younger than this period are more similar to the overall disk properties.This is an indication that the clusters formed during this period have significantly different kinematic parameters.We speculate that this is the first evidence of direct collision in the LMC disk and the spatial and temporal kinematic disturbance, matching with the predicted time and location of the LMC-SMC collision.We also speculate that the interaction mainly affected the gas, and the clusters and stars born from the disturbed gas bear the signature of the perturbation.

SUMMARY
We summarise the results and conclusions from this comprehensive study of the LMC disk kinematics using Gaia DR3 data.
1. We performed the kinematic model of the LMC that corresponds to the observed median PM using 1705 star clusters and field regions.This is a comprehensive study with a total of 48 data sets analyzed to create 48 models to study the dependence of the kinematic model on various parameters.
2. The model PM for the tracers in the data sets is formulated based on the equations outlined by V02.The model parameters estimated include inclination of the LMC disc (i), the position angle of the line of nodes (Θ), dynamic centers (α 0 , δ 0 ), the amplitude of the tangential velocity of the LMC's center of mass (v t ), the tangential angle made by v t (θ t ), scale radius (R 0 ), optimization factor (η), and the amplitude of the rotational velocity (v 0 ).We assumed a fixed distance to the LMC center and a line-of-sight velocity of the COM, v sys .The fitting of the parameters was performed by an MCMC technique, and model parameters were estimated for all datasets listed in Table 1.
3. This is the first 2D kinematic model of the LMC employing clusters and neighboring field regions with Gaia DR3 data.The data coverage of the LMC considered in this study is within ∼ 7 • from the LMC center.The parameters estimated in this study show good agreement with estimations in the literature when the comparison is made between similar populations.
4. There is no significant difference between the observed COM PM between cluster and field.We note an offset of 28±8 arcmin between the dynamic centers (α 0 , δ 0 ) of the cluster and field.Estimated Θ values of 122 for clusters and 128 • .02+0.79 −0.81 for field regions point to an offset of offset of 5 • .8±1.7 between them, while the inclination, i, remains almost similar.The modeled rotational parameters, R 0 and η 0 appear to be larger for field regions compared to clusters, while the v 0 remains almost similar.
5. We also estimated the kinematic model parameters for two control populations, RC stars and young MS stars.
We find that the younger population tends to show a southern and eastward dynamic center and a relatively smaller value of R 0 .We also note a varying value of Θ across the control population.
6.This study establishes that the kinematic model of the LMC disk varies with the age of the cluster population used for the estimation, in contrast to the surrounding field population.
7. The estimated parameters show a significant radial dependence for both cluster and field population.The value of i steeply decreases from the inner to outer radii, from ∼ 50 • to ∼ 32 • for clusters, and ∼ 43 • to ∼ 31 • for field population.The value of Θ increases from inner to outer radii, from ∼ 108 • to ∼ 122 • for clusters, and ∼ 110 • to ∼ 128 • for field population.The rotational parameters, R 0 and v 0 after 5 • show convergence to the values estimated from the full sample.
8. Clusters show larger PM residuals when compared with the field population.The RMS distribution of the residual PM for clusters shows a broader profile, while the corresponding distribution for the field shows a narrower one.9.The rotational velocity amplitude (v 0 ) for clusters, fields, and RC stars (∼ 81 to 84 km s −1 ) are consistent with the estimates by W20.The young MS stars in our study show a slightly larger v 0 (∼ 91 km s −1 ) when compared to the RC stars, as also noted by N22.
10.The modeled rotational parameters, R 0 and η 0 appear to be larger for field population when compared to clusters, while the v 0 remains almost similar.We note a significant shift in R 0 of ≈ 1.7 kpc between young MS and RC population.The variation in the value of R 0 with respect to population may be due to the redistribution of mass in the inner regions, and it may be due to the evolution of the bar and its activity over time.
11.The dispersion of the rotation velocity is found to be ∼ 23 km s −1 for clusters, and is likely to have contributions from the bar and spiral arms.The value of ∼ 11 km s −1 for field population is relatively low and suggests insignificant stellar disk heating.Young MS stars show a relatively large velocity dispersion, similar to clusters.
12. The residual PM for the cluster and the field decreases from the center up to 3 • , then increases beyond 4 • .The increased value in the inner 2 • region is likely to be the effect of the bar.We also detect evidence for non-circular motion in the bar region among the clusters.

Figure 1 .
Figure 1.The corner plot representing the sampled posterior distribution of nine kinematic parameters for the primary cluster data set is shown here.The vertical red lines represent the median values, and the black dashed lines represent the 16th and 84th percentiles.

Figure 2 .
Figure 2. The corner plot representing the sampled posterior distribution of nine kinematic parameters for the primary field data set is shown here.The vertical blue lines represent the median values, and the black dashed lines represent the 16th and 84th percentiles.

Figure 3 .
Figure 3.The observed PM plots for the clusters and field in the LMC sky plane are shown here.(a) Observed bulk motion of clusters; (b) Internal rotation PM of clusters; (c) Observed bulk motion of field; (b) Internal rotation PM of field.

Figure 4 .
Figure 4. Spatial distribution of the rotational velocity (Vrot) of the LMC for the clusters and the three control populations are depicted here.

Figure 5 .
Figure 5.The variation of kinematic parameters based on cluster richness is shown here.Estimated parameters for cluster groups (C10 to C60) and nearby field groups (F10 to F60) are shown in red and blue colors from Panel (a) to (f).(a) Variation in (α0, δ0); (b) Variation in i; (c) Variation in Θ; (d) Variation in (µW,com, µN,com); (e) Variation in R0; (f) Variation in v0.

Figure 6 .
Figure 6.The variation of kinematic parameters based on the spatial coverage of the LMC is shown here.Estimated parameters for different spatial coverage from 2 • to 7 • with a step size of 0.5 • from the center (based on the primary model) of the LMC are shown here.The corresponding cluster data sets (Cr=2 to Cr=7) and nearby field data sets (Fr=2 to Fr=7)are shown in red and blue colors from Panel (a) to (f).(a) Variation in (α0, δ0); (b) Variation in i; (c) Variation in Θ; (d) Variation in (µW,com, µN,com); (e) Variation in R0; (f) Variation in v0.

Figure 7 .
Figure 7. Residual PM of clusters and field population.(a) Spatial plot of residual PM vectors for the clusters; (b) Spatial distribution of residual PM vectors for the field population; (c) The Gaussian KDE showing the distribution of the |Residual PM|.

Figure 8 .
Figure8.The parameter space of the estimated (α0, δ0) and (µW,com, µN,com) are compared with the reference studies (Table7) are provided in panels (a) and (b), respectively.The corresponding estimated parameters for the clusters (CLS), field (FLS), young MS stars (YMS), and RC stars (RCS) are marked with red dots along with labels.The reference studies are marked with black dots along with labels.

Figure 9 .
Figure 9.The rotation curve of the LMC based on the parameter estimation for the control sample (Clusters, Nearby field, young MS stars, and RC stars) are plotted in panels (a) to (d).The magnitude of the rotational velocity in the LMC plane (Vrot) is shown with black dots in each panel.The red curve represents the running average over the observed Vrot with a bin size of 0.25 kpc, and the blue curve shows the best-fitting model depending on the rotational parameters (v0, R0, η) estimated for each data set (see Table3).
Figure 9.The rotation curve of the LMC based on the parameter estimation for the control sample (Clusters, Nearby field, young MS stars, and RC stars) are plotted in panels (a) to (d).The magnitude of the rotational velocity in the LMC plane (Vrot) is shown with black dots in each panel.The red curve represents the running average over the observed Vrot with a bin size of 0.25 kpc, and the blue curve shows the best-fitting model depending on the rotational parameters (v0, R0, η) estimated for each data set (see Table3).
density of such clusters, as illustrated in Figure12(a) with the 2D Gaussian KDE.We find a significant density of clusters in the North-West at ∼ 1-2 • with large residual PM.The bar region extending from North-West to South-East of the LMC shows relatively large residuals in PM.

Figure 10 .
Figure 10.Variataion of R0 and Vrot/σrot for the control population (Clusters, Nearby field, young MS stars, and RC stars).(a) Modeled velocity profiles showing the variation in R0 across the control sample; (b) Variation of Vrot/σrot across the control sample, estimated using a bin size of 0.4 kpc in radii.

Figure 11 .
Figure 11.The radial variation of mean |Residual PM| in the sky plane.Annular regions with a bin size of 1 • are used.Mean |Residual PM| of clusters and field population are marked with red and blue dots, respectively.

Figure 12 .
Figure 12.The clusters of larger residual PM are plotted here.(a) Spatial 2D Gaussian KDE of cluster locations with > 3σ of |Residual PM|; (b) Internal PM vectors of clusters > 3σ and < 3σ in Vrot profile of Figure 9(a); (c) Residual PM vectors of clusters with model and observed internal PM vectors differing by θr >45 • and θr >90 •.

Figure 13 .
Figure 13.The variation of i, Θ, and R0 with respect to the cluster age groups (subsection 5.5) are shown here.Red, Blue, and Green markers are used for clusters, fields, and young MS stars, respectively.(a) Variation of i; (b) Variation of Θ and (c) Variation of R0 are shown for 5 age groups.The age groups are valid only for clusters (Section 2).The field and the young MS stars are age-wise heterogeneous but located near the clusters within each age group.

Table 2 .
The kinematic best-fitting parameters obtained for the LMC from the primary model of clusters and nearby field regions are tabulated below.

Table 3 .
The kinematic best-fitting parameters obtained for the LMC based on the control sample (clusters, field population, young MS stars, RC stars) are tabulated below.

Table 4 .
The kinematic best-fitting parameters obtained for the LMC based on the cluster age groups and nearby field regions are tabulated below.As mentioned in Section 2, the associated field population for each cluster age group is age-wise heterogeneous.

Table 5 .
The kinematic best-fitting parameters obtained for the LMC based on the cluster richness and their nearby field regions are tabulated below.

Table 7 .
Comparison of the estimated kinematic parameters with previous studies.The first two rows list the estimations from this study.
b Estimated with fixed (α 0 , δ 0 ) based on G21, using VMC data c Estimated with fixed (α 0 , δ 0 ) and COM PM, using RC population with Gaia EDR3 d Estimated with Gaia EDR3 e Estimated with fixed (α 0 , δ 0 ), using Gaia EDR3 f Estimated using carbon stars, with SkyMapper DR1 g Estimated using RGB stars, with SkyMapper DR1 h Estimated using Young stars, with SkyMapper DR1 i Estimated with the third epoch of Hubble Space Telescope (HST) data, using PM.j Estimated with the third epoch of HST data, using PM + v LOS of old stars.k Estimated with the third epoch of HST data, using PM + v LOS of young stars.