© The Authors 2026

1 Introduction

The central few hundred parsecs of the Milky Way (MW) contain a nuclear stellar disc (NSD): a dense, rapidly rotating stellar structure that dominates the gravitational potential at galactocentric radii 30 < R < 300 pc (Schultheis et al. 2025). The NSD has a radial scale length of R ≃ 90 pc, a vertical scale-height of H ≃ 30–45 pc, and a stellar mass of M_NSD = ${10.5}_{-1.0}^{+1.1}$ × 10⁸ M_⊙ (Launhardt et al. 2002; Nishiyama et al. 2013; Schödel et al. 2014; Gallego-Cano et al. 2020; Sormani et al. 2022).

Whether the MW’s NSD is purely axisymmetric remains an open question. Indeed, while a fair fraction of the current literature has treated the NSD as an axisymmetric structure (e.g. Nogueras-Lara 2022), some observational studies previously suggested the existence of a nuclear bar in the MW. Alard (2001) used H- and K-band star counts from 2MASS to derive the projected stellar surface density in the bulge region. After subtracting the main bar contribution, they found a residual nuclear component with higher star counts at negative longitudes than at positive longitudes, opposite to the orientation of the main bar, which they interpreted as a signature of a nuclear bar. Rodriguez-Fernandez & Combes (2008) derived the mass distribution in the Galactic Centre by fitting the 2MASS star count map to a model including a disc, bulge, and nuclear bar. Using this inferred gravitational potential, they performed hydrodynamical simulations to model the behaviour of the interstellar gas. The resulting gas flow patterns reproduced a characteristic parallelogram feature observed in the longitude-velocity diagram of molecular spectral lines, which they interpreted as evidence of a nuclear bar whose major axis is oriented at an angle of α ≈ 60–75° relative to the Sun-Galactic Centre line. Nishiyama et al. (2005) used red clump stars at a latitude of b = +1° to trace the main bar orientation as a function of longitude (−10.5° < ℓ < 10.5°). They reported a change in the inclination angle for the Galactic bar at fields with |ℓ| < 4° (corresponding to a nuclear bar radius R ~ 570 pc), which they interpreted as the effect of a distinct inner bar structure. Gonzalez et al. (2011) later confirmed this result using red clump stars from the VISTA Variables in the Via Lactea (VVV) survey to trace the orientation of the largescale bar at latitudes b = ±1°. While the above results seem compelling, Gerhard & Martinez-Valpuesta (2012) employed N-body simulations to show that these observed asymmetries and the change in the slope of the star count profiles can be explained without invoking a nuclear bar. Instead, they demonstrated that the observed features arise naturally from the intrinsic density distribution of a single, barred, boxy bulge that transitions from elongated to more axisymmetric within the inner few degrees. This transition causes an apparent flattening of the longitude profiles near |ℓ| ≈ 4°, similar to the one observed by Nishiyama et al. (2005) and Gonzalez et al. (2011), mimicking the expected signature of a nuclear bar.

A key challenge in determining whether a nuclear bar is present in the MW stems from our edge-on perspective of the Galactic Centre, which is heavily obscured by non-uniformly distributed dust (e.g. Bally et al. 1988; Gonzalez et al. 2012; Schultheis et al. 2014; Nogueras-Lara et al. 2021; Henshaw et al. 2023). Notably, most of the gas and dust (~3/4) within |ℓ| < 3.5° resides at positive longitudes (e.g. Bally et al. 1988; Henshaw et al. 2023), which complicates the use of star counts and photometric methods by producing an apparent deficit of stars at positive longitudes. This introduces significant biases in the results, challenging the reliability of photometric analyses.

An alternative approach to detecting a nuclear bar involves the use of stellar kinematics. Sormani et al. (2022) developed a dynamical equilibrium model fitted to 1D kinematic data from the KMOS NSD survey (Fritz et al. 2021). Although their results indicate that the model fits are consistent with an axisymmetric system, it remains uncertain whether the 1D velocity histograms employed in their analysis are adequate to uncover evidence of non-axisymmetry. Nonetheless, nuclear bars are relatively common features in MW–like galaxies, with roughly 50% of barred galaxies with stellar masses above log(M_⋆/M_⊙) > 10.5 being double-barred (Erwin 2024). Therefore, the possibility that our Galaxy hosts a nuclear bar warrants careful consideration.

Despite the current small sample sizes of existing 6D kinematic datasets, the anticipated availability of large and high-quality datasets from upcoming surveys such as the KMOS VVVX-GalCen¹ (Nogueras-Lara et al., in prep.), the Galactic Bulge Time Domain Survey with Roman (Terry et al. 2023), MOONS (Gonzalez et al. 2020), and JASMINE (Kawata et al. 2024) will enable the dissection of the NSD region in unprecedented detail, allowing for new insights into the formation and evolution of the Galactic Centre. The goal of this paper is to make testable kinematic predictions that can confirm or rule out the presence of a nuclear bar in the MW with future surveys.

This paper is organised as follows. In Sect. 2, we describe the models employed in this study. In Sects. 3 and 4, we present the two main kinematic diagnostics for detecting a nuclear bar:² the vertex deviation and the asymmetry observed in the μ_ℓ versus Galactic longitude, ℓ, distribution. In Sect. 5, we provide a discussion of the results and their implications. Lastly, Sect. 6 briefly outlines prospects for future work.

Fig. 1 Face-on (bottom) and side-on (top) stellar density distribution of the models. Left to right, the axisymmetric NSD, isolated strong NSB, isolated weak NSB, and large-scale bar. The large-scale bar model is oriented at −27° relative to the Solar-GC line (dotted line). The nuclear bar angle, α, is a free parameter. In all cases, the Solar position is (x, y, z) = (−8.2, 0, 0.012) kpc.

2 The models

We considered three models of nuclear structure:

• NSD: an axisymmetric nuclear stellar disc (Sect. 2.1).

• W-NSB: a weak nuclear bar (Sect. 2.2).

• S-NSB: a strong nuclear bar (Sect. 2.2).

In addition, we employed a large-scale MW model that includes the outer bar (Sect. 2.3). We constructed the full models by embedding the nuclear components at the centre of the largescale MW model, which includes an outer bar but no nuclear structure. Note that for the rest of the paper, when referring to the S-NSB and NSD models, these are always embedded in the large-scale bar. Otherwise, we refer to the ‘isolated’ NSD or NSB models. Fig. 1 shows edge-on (top) and face-on (bottom) views of the models. Fig. 2 shows the surface density in the plane of the sky of the large-scale MW model without a nuclear component (left), the surface density of the isolated NSD model (centre), and the ratio between the two (right). To create the density panels, we rescaled the large-scale bar model (see Sect. 2.3 for details) to achieve an NSD/bar surface density ratio consistent with Sormani et al. (2022). The combined models were then constructed by randomly sampling stars from each component (nuclear and large-scale bar), with the relative sampling probability at each position (ℓ, b) determined by the surface density ratio shown in the right panel of Fig. 2. This ensures the correct levels of contamination from large-scale bar stars. When sampling from the large-scale bar component, stars are drawn from the original, unscaled model to preserve the correct kinematic properties. The large-scale MW model is described in more detail in Sect. 2.3.

We note that vertical proper motions (μ_b) could in principle be used to reduce contamination from large-scale bar stars, as these exhibit systematically higher |μ_b| compared to the kinematically cooler NSD (e.g. Shahzamanian et al. 2022; Nogueras-Lara 2022). However, as we demonstrate in Appendix F, applying a conservative cut of |μ_b| > 5 mas yr⁻¹ reduces bar contamination by at most ~2%, which has a negligible impact on our kinematic diagnostics.

Lastly, although the large-scale MW model and the NSD/NSB models are each internally self-consistent, we do not evolve them together under their combined gravitational potential. Therefore, the combined models are not fully self-consistent. However, this has a negligible impact on the results given the very different spatial scales of the models.

Fig. 2 Surface density of the large-scale MW model with no nuclear component (left), the isolated NSD (middle), and the ratio between the two (right). The ratio illustrates how prominent the MW’s NSD would be if we were to see it from the Sun (see Fig. 1).

2.1 The axisymmetric NSD model

We adopt the axisymmetric, self-consistent dynamical model of the NSD from Sormani et al. (2022)³. This model was fitted to the line-of-sight and proper motion kinematic distributions obtained by combining the KMOS infrared (K-band) spectroscopic survey of Fritz et al. (2021) with a preliminary version of the VIRAC2 proper motion catalogue (Smith et al. 2025). The final fiducial NSD model has a total mass of M_NSD = 9.7 × 10⁸ M_⊙, an exponential radial scale length of R_disc = 74 pc and a vertical scale length of H_disc = 26 pc.

2.2 The strong and weak NSB models

To obtain isolated nuclear bar models, we evolve the isolated NSD model of Sormani et al. (2022) within the self-consistent MW model of Binney & Vasiliev (2023)⁴ using the N-body code GADGET-4 (Springel et al. 2021). We use softening lengths of 5,50, and 25 pc for the NSD, dark matter, and main stellar disc components, respectively. We find that the initially axisymmetric NSD is mildly unstable, and starts forming a nuclear bar after evolving for 0.6 Gyr. The left panel of Fig. 3, illustrates the evolution of the A₂ bar amplitude, computed via the usual m = 2 Fourier moment (e.g. Li et al. 2024), alongside the radius containing 98% of the NSD’s mass (right panel). The output cadence is sufficiently fine that the A₂ profile is well-sampled across the full evolution; a running average is applied in Fig. 3 to reduce shot noise between successive outputs. The bar amplitude peaks at t = 1.7 Gyr, after which it begins to weaken. We extract two snapshots from this simulation to construct our barred NSD models: the strong NSB (S-NSB) at t = 1.70 Gyr (black dashed line) when the bar is at peak strength, and the weak NSB (W-NSB) at t = 2.60 Gyr (blue dashed line) during the bar’s declining phase.

Fig. 3 Evolution of the NSD model. Left: bar strength, measured by the Fourier A2 parameter, as a function of time. The parameter A2 was calculated considering all stars within the radius Rmax that encloses 98% of the model’s mass. Right: evolution of Rmax with time. The initially stable NSD starts forming a bar at t ≈ 0.5 Gyr, which reaches maximum strength at t ≈ 1.7 Gyr. Black and blue dashed lines indicate the times at which the S-NSB and W-NSB models are extracted, respectively.

Fig. 4 Plots in galactocentric Cartesian coordinates of the isolated NSD (left), strong NSB (middle), and weak NSB (right) models oriented at −60° relative to the Sun-GC line. Overlaid are the velocity ellipses coloured by the vertex deviation, lv. We have placed an ellipse in each bin with N > 25 stars. The dashed lines reflect a 2° field of view centred on the Solar position at (x, y) = (−8.2, 0) kpc. The NSD shows radially oriented ellipses, while the S-NSB model ellipses show strong alignment with the bar’s major axis. The alignment is weaker for the W-NSB model.

2.3 The large-scale MW model

We employ a self-consistent N-body+Smooth Particle Hydrodynamic (SPH) model of a large-scale MW-like bar which corresponds to the model TG07v3 in Deg et al. (2025). Here, we provide only a brief overview.

The initial conditions consist of a bulgeless stellar disc with a Sérsic surface density profile (Sérsic index n = 1.75), a gas disc (7% of the stellar mass), and an NFW dark matter halo (Navarro et al. 1997) with mass M_h = 1.17 × 10¹² M_⊙ and scale radius R_h = 19 kpc. The stellar disc has a mass of M_d = 5.5 × 10¹⁰ M_⊙ and scale length R_d = 0.43 kpc (for the Sérsic profile), while the gas disc has an exponential scale length of R_g = 6.5 kpc. The model is evolved using CHANGA (Menon et al. 2015) for 10 Gyr using a base timestep of Δt = 0.5 Myr, which provides sufficient temporal resolution that our results are insensitive to the choice of output cadence. Star formation occurs when gas density exceeds 0.1 amu cm⁻³ and temperature falls below 15 000 K, with a star formation efficiency of 5% per dynamical time and Supernova feedback via the Keller et al. (2014) prescription.

The resulting model develops a bar that reaches a stable configuration by ~3.5 Gyr and remains approximately constant thereafter. At t = 3 Gyr, the bar has a semi-major radius of R_bar ~ 4.25 kpc and pattern speed of Ω ~ 42 km s⁻¹ kpc⁻¹, consistent with measurements of the MW (e.g. Wegg et al. 2015; Portail et al. 2017). The model also reproduces the characteristic bimodal distribution of red clump stars observed along various lines of sight, mimicking the observed morphology of the MW’s boxy/peanut bulge (Gonzalez et al. 2015). We note that this model does not include a metal-poor spheroidal stellar component in the bulge; we assess the potential impact of such a population in Appendix A and find it has negligible effect on our results.

To match the observed NSD/bar surface density ratio from Sormani et al. (2022), we rescale the model positions by α = 0.95 and masses by γ = 0.44. This rescaling is essential to prevent the large-scale bar from dominating the surface density in the Galactic centre region. To maintain consistency with the laws of motion, we also rescale time by β = 1.39 and velocities by δ = 0.68. The resulting scaled pattern speed is Ω_scaled ≈ 31 km s⁻¹ kpc⁻¹ and the bar semi-major axis is R_bar,scaled ≈ 4.03 kpc, both of which are consistent with observational constraints (e.g. Minchev et al. 2007; Portail et al. 2015; Pontzen et al. 2013; Sanders et al. 2019; Bovy et al. 2019; Chiba & Schönrich 2021; Clarke & Gerhard 2022).

3 Diagnostic 1: Vertex deviation

Given a set of stellar velocities, the velocity-dispersion tensor can be defined as $$\begin{array}{rl}{\sigma }_{ij}^2& =⟨(v_i-⟨v_i⟩)(v_j-⟨v_j⟩)⟩\\ & =⟨v_i{v}_j⟩-⟨v_i⟩⟨v_j⟩,\end{array}$$(1)

where i and j represent velocity components in the adopted orthogonal coordinate frame (Binney & Tremaine 2008). In the analysis in this section, we consider a 2D tensor where i and j can represent either the line-of-sight direction or the transverse direction parallel to the Galactic plane in the direction of increasing longitude ℓ. We denote the velocities in these two direction as v_los and v_ℓ respectively. The dispersion tensor is symmetric and can therefore always be diagonalised. The orthogonal eigenvectors (corresponding to eigenvalues σ_max and σ_min) define the principal axes of the velocity ellipse, with the major axis indicating the direction of maximum velocity dispersion (see Fig. B.1 and Appendix B).

The orientation of the velocity ellipse in the (v_los, v_ℓ) plane can be characterised by the vertex deviation, l_v, which is given by: $$l_{\mathrm{v}}=\frac{1}{2}\arctan \left(\frac{2{\sigma }_{12}^2}{{\sigma }_{11}^2-{\sigma }_{22}^2}\right),$$(2)

where 1 and 2 denote the v_los and v_ℓ directions, respectively (see Appendix B). By definition, l_v takes values −90° ≤ l_v ≤ 90°.

In an axisymmetric nuclear disc, since ${\sigma }_{12}^2\approx 0$, the velocity ellipsoid of all stars along a given line of sight is expected to be approximately aligned with one of the coordinate axes (i.e. l_v ≈ 0° or l_v ≈ ±90°) depending on the direction of highest velocity dispersion. For example, ellipses on opposite sides of the x = 0 line in the left panel Fig. 4 cancel out and their “average” orientation is aligned with the x and y axes. Note that strictly speaking this is valid only for an observer at infinite distance. For an observer at finite distance, the lines-of-sight that do not pass exactly through the Galactic centre (e.g. the dashed lines in Fig. 4) are not exactly parallel to the x axis. However, the nuclear disc is small enough that this effect is negligible.

In the presence of a bar instead, stars preferentially stream along the bar’s major axis, leading to a non-zero covariance ${\sigma }_{ij}^2$ and hence a non-vanishing l_v (see middle and right panels of Fig. 4). The vertex deviation has been successfully employed in both numerical (e.g. Debattista et al. 2019) and observational (e.g. Simion et al. 2021; San Martin Fernandez et al. 2025) studies of the MW bulge to detect and characterise the Galactic bar.

The aspect ratio of the velocity ellipse is characterised by the anisotropy parameter: $$\beta =1-\frac{{\sigma }_{min}^2}{{\sigma }_{max}^2},$$(3)

which takes values 0 ≤ β ≤ 1. For β = 0, the velocity distribution is isotropic (circular), and becomes increasingly elongated as β → 1. One important limitation of using the vertex deviation is that l_v becomes increasingly unstable as β → 0. In such cases, for a stellar population lacking a preferred direction along either coordinate axis (i.e. σ_min ≈ σ_max), the vertex deviation becomes poorly defined and lacks statistical significance when assessed via bootstrapping (Eq. (4)).

Fig. 5 Velocity ellipses for the models, with the major and minor axes shown in solid and dashed lines, respectively. We uniformly sampled N = 2400 stars from the region |ℓ| < 0.9° and |b| < 0.25°. In each panel we denote the vertex deviation (lv), the anisotropy (β), the bar angle α (except for the NSD case), and the velocity dispersions in the line of sight and longitude directions. The top-left panel shows the ellipse for the large-scale MW model with no nuclear component. As expected, the direction of highest dispersion is in the line of sight direction. The bottom left shows the ellipse for the axisymmetric NSD model. The other panels show the ellipses for the S-NSB model with various nuclear bar angles. In this case, the direction of highest dispersion depends on the nuclear bar angle.

3.1 Error handling

To quantify the uncertainties, we employ a bootstrap resampling method. Starting with a sample of N stars, we generate B resampled datasets, each maintaining the original sample size. For each re-sampling, we compute the vertex deviation $l_{\mathrm{v}}^{\ast }$. This process yields a distribution of $l_{\mathrm{v}}^{\ast }$ values, from which the uncertainty is derived as: $${\sigma }_{l_{\mathrm{v}}}=\sqrt{\frac{1}{B}\sum _{k=1}^B{(l_{\mathrm{v},k}^{\ast }-l_{\mathrm{v}})}^2},$$(4)

where l_v is the estimate from the original sample, and $l_{\mathrm{v},\mathrm{k}}^{\ast }$ are the bootstrap samples.

We also assess the impact of observational uncertainties in v_ℓ and v_los on the derived l_v measurements. Using typical KMOS uncertainties (~10 km s⁻¹ for v_los and ~23 km s⁻¹ for v_ℓ based on Fritz et al. 2021 and VIRAC2), we applied a Monte Carlo approach convolving model velocities with these errors. In agreement with San Martin Fernandez et al. (2025), we find that for a given sample size, N, these observational uncertainties are negligible (σ_ℓv < 5°); therefore, all uncertainties reported use bootstrapping only.

3.2 Vertex deviation in the models

We explore whether the vertex deviation can serve as a diagnostic for the presence of a nuclear bar. Specifically, we investigate: (1) whether a nuclear bar produces an l_v signature which is distinct from that of an NSD, (2) the number of stars required to reliably measure l_v, and (3) the impact of dust extinction on the measurement. We first consider results from the strong nuclear bar model (S-NSB). The impact of bar amplitude is assessed in Sect. 3.5.

To develop intuition for the spatial variation of l_v across the models, we present in Fig. 4 face-on density maps of each nuclear model, subdivided into a 10 × 10 spatial grid. At each grid cell, we overlay the corresponding velocity ellipse measured in (v_ℓ, v_los) space. The dashed lines indicate a 2° field of view centred on the Galactic Centre as seen from the Solar position at (x, y) = (−8.2, 0) kpc. In the isolated NSD case (left panel), the major axes of the velocity ellipses are oriented radially towards the Galactic Centre, reflecting the dominant contribution of line-of-sight velocity dispersion in this axisymmetric system. In contrast, for the strong NSB model (middle panel), the ellipses are strongly aligned with the major axis of the nuclear bar (α = −60°), tracing the streaming motions of stars on bar-supporting orbits. In the weak NSB case (right panel), the alignment of the ellipses along the bar major axis is substantially weaker, as expected from increased levels of symmetry.

In Fig. 5, we present the velocity ellipses derived from our models. The ellipses are constructed using N = 2400 stars sampled uniformly from the region |ℓ| < 0.9° and |b| < 0.25°. Within each panel, key parameters are indicated: the vertex deviation l_v, the anisotropy β, the bar angle α (where applicable), and the velocity dispersions along the line-of-sight (σ_los) and along the Galactic longitude (σ_ℓ). The top-left panel illustrates the velocity ellipse for the large-scale bar model without a nuclear component. In this configuration, the velocity dispersion is highest along the line-of-sight direction, consistent with observations of the MW of metal rich ([Fe/H] > 0.30) stars tracing the largescale bar (e.g. Babusiaux et al. 2010). The bottom-left panel shows the velocity ellipse for the axisymmetric NSD model, which shows velocities only weakly longitudinally biased (also reflected by the low anisotropy β = 0.15). As a result, the vertex deviation for the axisymmetric NSD model is l_v ≈ ±90°. The remaining panels display the velocity ellipses for the S-NSB model at various nuclear bar angles α. Compared to the NSD model, the anisotropy is high (β ≈ 0.50) across all angles. As expected, the value of the vertex deviation, l_v roughly matches the nuclear bar angle, α.

Fig. 6 presents 2D maps in (ℓ, b)-space with bins coloured by l_v. Overlaid are the fractional density contours from the right panel of Fig. 2 for reference. As a control, the top-left panel shows the large-scale bar model, which lacks a nuclear component. The inset shows the orientation of the large-scale bar (red line), which is −27° relative to the Sun-GC line. Consistent with expectations, l_v is negative across the region, reflecting the outer bar’s orientation (Zhao et al. 1994; Soto et al. 2007; Babusiaux et al. 2010; Debattista et al. 2019; Simion et al. 2021; San Martin Fernandez et al. 2025).

The top-right panel of Fig. 6 shows the axisymmetric NSD model (indicated by the blue circle in the inset). As shown in Fig. 5, the velocity ellipse of the axisymmetric NSD is characterised by low anisotropy with a slight bias towards higher dispersion in the Galactic longitude direction. From the view-point of the Sun, this results in l_v ≈ 90° in the region |ℓ| ≲ 1° and |b| ≲ 0.25°. The middle and bottom set of panels in Fig. 6 present the results for the S-NSB and W-NSB models, respectively, with the nuclear bar oriented at various angles α relative to the Sun-GC line. The insets show the NSB orientation (blue line), while the outer bar remains fixed at −27° in all cases. For both the S-NSB and W-NSB models, a nuclear bar at α = 60°, produces positive values of l_v. More generally, nuclear bar orientations in the range 0° ≤ α ≤ 90° yield positive l_v, whereas orientations in the range −90° ≤ α ≤ 0° produce negative l_v values. As we show later in Sect. 3.5, l_v varies smoothly with the NSB angle. In the limiting case of α = 90°, the resulting signal is difficult to distinguish from that produced by an axisymmetric NSD.

Fig. 7 shows the variation of l_v (top) and β (bottom) with Galactic latitude |b| for the NSD (black curve). The other curves show the variation for the S-NSB model with the nuclear bar at varying angles, α. Here, we only consider stars satisfying |ℓ| < 1°. As expected, l_v for the NSD model is l_v≈90° in the midplane, while this value changes for the S-NSB model depending on α. We find that in the midplane β ≈ 0.22 for the NSD model, while it takes on larger values (0.4 < β < 0.5) for the S NSB model. In all cases, as |b| increases, both l_v and β converge towards a single value. This is due to increased contamination from stars associated with the large-scale bar. At |b| ≈ 0.5°, the vertex deviation plateaus around l_v ≈ − 27°. This marks the transition point where l_v becomes fully dominated by the large-scale bar.

Fig. 6 Two-dimensional maps in (ℓ, b) space coloured by the vertex deviation, lv. Top row: maps for the large-scale Milky Way model without a nuclear component (left) and with an axisymmetric NSD model (right). The NSD is indicated by the blue circle in the inset. Middle panels: maps for the S-NSB model with the nuclear bar at various orientations, α. Insets show the geometric configuration of the outer bar (red line) and nuclear bar (blue line). White markers in the top-left panel indicate the fields used in the analysis of Sect. 3.3. Bottom panels: same as above but for the W-NSB model.

Fig. 7 Vertex deviation, (lv) and anisotropy, (β) as a function of Galactic latitude, |b|, for the NSD (black curves) and S-NSB model with varying values of the nuclear bar angle, α (coloured curves). As expected, the NSD has lv ≈90° in the midplane. In contrast, the midplane values for the S-NSB model are dependent on α. Across all model configurations, both lv and β converge to roughly similar values at |b| ≥ 0.5°. This marks the latitude at which contamination from large-scale bar stars dominate the density distribution.

3.3 The impact of sample size on vertex deviation

We now explore the dependence of the precision of l_v on the sample size. In Fig. 8, we show the distribution of $l_{\mathrm{v}}^{\ast }$ (defined in Eq. (4)) across four representative fields indicated by the white markers in the top-left S-NSB panel of Fig. 6. The sky Galactic coordinates are denoted in each panel. The left and right group of panels show the results obtained using the S-NSB and W-NSB models, respectively. In all cases, the models have a nuclear bar angle α = −60°; however, we have checked that the results are consistent across all α. We perform the bootstrap re-sampling method with B = 1000 iterations for each field (panel), considering sample sizes of 120 (dotted red), 500 (dashed green) and 1000 (solid blue) stars per field. For the S-NSB case, the average uncertainty (across all panels) in l_v is σ_lv = ±6.18° for a sample of 120 stars, which decreases to around σ_lv = ±2.97° for 500 stars. We find that for the case of a strong nuclear bar (but see Sect. 3.5 for the weak case), around 500 stars in each field may suffice for a reliable measurement of l_v. On the other hand only 120 stars per field may lead to spurious measurements. We also note that some variation in l_v is expected across the fields due to both the different levels of contamination from the large-scale bar, and due to the fields probing different kinematic regimes. Across the four fields in Fig. 8, l_v ranges from −66.4° to −62.8° (with σ_ℓv = ±1.13°) for N = 1000 (blue histogram).

Combining stars from multiple fields to increase statistical significance of l_v is a potential strategy for further reducing the overall uncertainty. To illustrate this, Fig. 9 shows the uncertainty in l_v as a function of sample size, N, for stars sampled from the fields in Fig. 8. The red and green curves show results for the S-NSB and W-NSB (discussed in Sect. 3.5) model, respectively. In both cases, the uncertainty decreases rapidly up to N ≈ 500, then more gradually thereafter. We extend the analysis to N = 3500, similar to the sample size in the current KMOS survey of Fritz et al. (2021).

3.4 The impact of observational biases on vertex deviation

To assess whether observational selection biases could systematically affect our kinematic measurements and lead to a false detection of a nuclear bar, we applied a selection function to the axisymmetric NSD model and re-measured the kinematic parameters. Our selection function incorporates three physical effects: distance-dependent detectability (favouring the near side of the NSD), midplane extinction and crowding (suppressing detection near b = 0°), and asymmetric dust extinction (suppressing detection at l > 0°). Full details of the selection function are provided in Appendix D. The probability distribution functions for each of these are show in the left column of Fig. 10.

The right column of Fig. 10 illustrates the impact of these selection effects on the stellar distributions. We compare the unbiased NSD model (black) to the ‘observed’ sample after applying the selection function (red) for distance (top panel), Galactic latitude (middle panel), and Galactic longitude (bottom panel). The selected sample shows the expected biases: stars are preferentially drawn from the near side (d < 8.2 kpc), higher latitudes (|b| > 0.1°), and negative longitudes where extinction is lower.

Fig. 11 shows the bootstrap distribution of the vertex deviation $l_v^{\ast }$ for the NSD model with (red) and without (blue) selection effects applied. Using a bootstrap approach, we resampled velocities for N = 12500 stars within the region (|ℓ| < 0.9° and |b| < 0.25°). For the NSD model without observational biases (blue histogram), we obtain a median value of l_v = −87.4 ± 5.2°. After applying the observational selection effects (red histogram), the vertex deviation measurement shifts to l_v = 82.2 ± 5.1°.

In general, we find that applying the selection function shifts the vertex deviation towards positive values. To understand this effect, consider once more the left panel of Fig. 4. The combined probability function preferentially samples stars which are on the near side of the NSD and towards positive Galactic longitudes. This translates to stars which are in the bottom left quadrant (x < 0, y < 0), where velocity ellipses with positive l_v dominate. Consequently, if extinction and observational biases are not carefully accounted for, this could lead to a spurious detection of a nuclear bar with its near side at negative Galactic longitudes.

We note that the distance bias component of our selection function (which preferentially samples foreground stars) has minimal impact on l_v when applied in isolation. This is because the distance bias affects stars on both sides of ℓ = 0° symmetrically, and their velocity ellipse orientations effectively cancel when combined (see left panel of Fig. 4). The dominant systematic effect arises from the asymmetric longitude-dependent extinction, which can be partially corrected observationally through extinction mapping.

3.5 A weak nuclear bar

We now investigate how nuclear bar amplitude influences l_v and β by repeating the vertex deviation analysis at t ≈ 2.60 Gyr, when the nuclear bar has declined to approximately 65% of its peak amplitude. Figure 12 shows l_v (left) and β (right) as a function of nuclear bar angle, α, for the S-NSB (red) and W-NSB (green) models. The primary bar angle is fixed at −27°. In each case, we sample 1800 stars within the region (|ℓ| < 0.9°, |b| < 0.25°), which we combine to compute the kinematic properties.

Bar amplitude has minimal effect on l_v (left panel), consistent with the findings of San Martin Fernandez et al. (2025), who showed that l_v is a relatively blunt probe of non-axisymmetry. Specifically, in the presence of a non-axisymmetric structure, l_v is determined by the nuclear bar orientation α, but it does not vary with the strength of the non-axisymmetry. In contrast, we find that β decreases significantly for the weaker bar, by approximately a factor of two, resulting in values comparable to those of the NSD model (β ≈ 0.15).

The more isotropic velocity distribution associated with the weaker bar (i.e. lower β) also leads to larger systematic uncertainties in l_v. In the right group of panels of Fig. 8, we repeat the error analysis of Sect. 3.3 for the weak bar case. We find greater variation in the bootstrapped l_v values, with uncertainties of 16.08°, 9.43°, and 6.52° for samples of 120, 500, and 1000 stars per field, respectively. Across the fields, the median l_v value ranges between −70.9° and −45.7°. This systematic increase in uncertainty is also evident in Fig. 9. The red curve shows the uncertainty in l_v as a function of sample size when combining stars across the central fields. Compared to the strong bar case, we find a systematic increase in the uncertainty by approximately 8° across all values of N_sample.

Fig. 8 Bootstrapped distributions of lv as a function of sample size. Using the S-NSB (left panels) and W-NSB (right panels) models with a nuclear bar angle α = −60°, we sample stars from fields shown with (ℓ, b) coordinates indicated in each panel. We also indicate the position of each field (white markers) in the top left panel for the S-NSB model in Fig. 6. For each field, we perform B = 1000 bootstrap iterations using sample sizes of 120 (dotted red), 500 (dashed green), and 1000 (solid blue) stars. For the S-NSB model, the mean uncertainty (using Eq. (4)) decreases from ±6.18° for 120 stars to ±2.97° for 500 stars, and ±2.01° for 1000 stars. We note that lv varies marginally across fields, ranging from −66.4° to −62.8° (with a standard deviation of σℓv ± 2.58°;). The uncertainty in lv for the W-NSB case is 16.08°, 9.43° and 6.52° for 120, 500 and 1000 stars per field respectively. Across the fields, we find that the median lv value ranges between −70.4° and −60.1°. This variation likely reflects both the intrinsic differences in stellar kinematics across fields and varying contamination from large-scale bar stars.

Fig. 9 Uncertainty in the vertex deviation σlv as a function of sample size, N. Stars are sampled from fields in Fig. 8. We show the S-NSB (red) and W-NSB (green) cases, both with α = −60°. The uncertainty decreases steeply up to N = 500, then more gradually thereafter. The W-NSB case yields systematically higher uncertainties (by ~8°) compared to S-NSB at all sample sizes due to its more isotropic velocity distribution (lower β).

Fig. 10 Left: observational selection functions applied to the models. (a) Distance-dependent detection probability favouring stars on the near side of the disc (d < 8.2 kpc). (b) Midplane extinction reducing completeness to 30% at b = 0°. (c) Asymmetric longitude extinction with steeper decline at positive longitudes. The total probability is the product Ptotal(d, ℓ, b) = P(d) × P(b) × P(ℓ). Right: impact of selection on the NSD model’s distance (top), latitude (middle), and longitude (bottom) distributions, showing the unbiased sample (black) versus the observed sample after applying selection effects (red).

Fig. 11 Bootstrapped distribution of $l_v^{\ast }$ for the axisymmetric NSD model without (blue) and with (red) the observational selection function applied. Each distribution results from 10 000 bootstrap iterations sampling N = 12500 stars from the region |ℓ| < 0.9° and |b| < 0.25°. The unbiased NSD yields lv = (−87.4 ± 5.2)°, consistent with an axisymmetric disc. However, applying realistic selection effects shifts the measurement to lv = (82.2 ± 5.1)° (towards positive values), demonstrating that observational biases could mimic a nuclear bar with its near side at negative longitudes if not properly corrected.

Fig. 12 Vertex deviation (left) and anisotropy (right) of 1800 stars sampled from the region |ℓ| < 0.9° and |b| < 0.25° as a function of the nuclear bar angle, α. Red curves show S-NSB results; green curves show W-NSB. While lv remains roughly identical between strong and weak nuclear bars, β decreases by approximately a factor of two for the weaker bar, approaching NSD-like values. This demonstrates that lv is a robust orientation diagnostic but a blunt probe of bar amplitude, while β may serve as a complementary bar amplitude indicator.

Fig. 13 Synthetic RGB image produced with PYNBODY (Pontzen et al. 2013) of the stellar disc of the N-body+SPH model. We show the model face-on with the main bar aligned roughly −27° to the Sun-GC line and with the Sun located at (x, y) = (−8.2, 0) kpc. The central region of this model hosts an orthogonal sub ~1 kpc NSD with its near side towards negative Galactic longitude. This configuration tests whether orthogonal NSDs can be distinguished from true nuclear bars using kinematic diagnostics.

3.6 An orthogonal nuclear disc

External galaxies frequently host elliptical NSDs (e.g. Bittner et al. 2020). In barred galaxies, NSDs are predominantly supported by the x₂ orbital family, with the NSD’s major axis oriented orthogonal to the primary bar (Cole et al. 2014). Since the Sormani et al. (2022) NSD model assumes axisymmetry, it cannot capture this orthogonal configuration. We therefore adopt a self-consistent N-body+SPH model (see Appendix C for details) to investigate the kinematic signature of an elliptical, orthogonally oriented NSD. We show a synthetic RGB rendering of the model in Fig. 13, which features a primary bar with radius R ≈ 3 kpc and a ~1 kpc-scale NSD. While not MW-scaled, the large NSD allows us assess how an orthogonal configuration impacts the l_v signal.

As before, the primary bar is oriented at −27° to the Sun–GC line, with its near side at positive longitudes. The orthogonal NSD has its near side towards negative longitudes, implying it should exhibit positive l_v if the large-scale bar shows negative l_v. We sample 2400 stars within the galactocentric radius R < 1 kpc which also satisfy (|ℓ| < 1°, |b| < 0.3°). The left panel of Fig. 14 shows the velocity ellipse for all stars satisfying the selection criterion. The resulting l_v and β values which are consistent with the axisymmetric NSD case (see Fig. 5). Since the NSD in this model formed late, we re-sample the model again (with N = 2400) only considering stars with ages τ_⋆ < 5 Gyr. This effectively reduces contamination from the large-scale bar. For this younger population, β increases by approximately a factor of ~3. We also observe that l_v shifts away from the coordinate axis, towards positive values.

Therefore, an orthogonally oriented NSD in the MW would produce positive l_v signal, mimicking the kinematic signature of a nuclear bar viewed nearly side-on with its semi-major axis at negative longitude (see Fig. 12). This degeneracy implies that weak positive l_v in the Galactic Centre is consistent with either scenario and does not exclude the presence of a mildly elliptical NSD.

Fig. 14 Velocity ellipses for the quasi-axisymmetric orthogonal NSD model, constructed from stars in the region |ℓ| < 0.9° and |b| < 0.25°. Left: all stars show lv ≈ −89.5° and low anisotropy (β = 0.22), appearing nearly axisymmetric due to contamination from older large-scale bar stars. Right: result when selecting younger stars (τ⋆ < 5 Gyr) that predominantly trace the late-forming NSD yields lv≈79.2° with higher anisotropy β = 0.62. This demonstrates that an orthogonal NSD can produce positive lv, similar to a nuclear bar with its near side at negative longitudes.

4 Diagnostic 2: Asymmetry in the (μ_ℓ, ℓ) plane

The (μ_ℓ, ℓ) plane provides a potential diagnostic for identifying an edge-on nuclear bar in the Galactic Centre. In an axisymmetric disc, the distribution of stars in the (μ_ℓ, ℓ) plane is symmetric about the ℓ = 0° line. Conversely, the presence of a nuclear bar can induce, depending on the orientation of the major axis, a measurable deviation from this symmetry.

To capture these asymmetries, we first bin stars by their Galactic longitude ℓ. Within each bin, a 1D histogram of μ_ℓ is constructed and smoothed using Gaussian kernel density estimation. We identify the peak of each histogram, denoted as μ_ℓ,max. Then we plot the relationship between μ_ℓ,max and ℓ.

Fig. 15 shows μ_ℓ,max versus ℓ for the S-NSB model (hollow red markers) with various nuclear bar angles, α. A smoothed fit to these points is shown by the red line. Blue markers and lines depict the corresponding distributions for the axisymmetric NSD case. The leftmost panels show the cases for α = 0° (top) and α = 90° (bottom), where the μ_ℓ,max − ℓ relationship is roughly consistent with that of the axisymmetric NSD. In contrast, for smaller α (second column), the distribution shows a ‘step’ feature in μ_ℓ,max which is not visible in the NSD case. As α increases (third/fourth columns) the prominence of the step features decreases quickly, eventually becoming flat.

The black open circles in Fig. 15 show the μ_l,max distribution for stars sampled from individual fields in the midplane in order to emulate observations. In each field, we sample N = 1500 stars from the S-NSB model, with uncertainties estimated via a bootstrapping approach. For each bin in ℓ, we performed resampling of the stars with replacement for 1000 iterations and the uncertainty determined via Eq. (4).

In cases where α is small (e.g. α ≲ 20°), the step feature in the distribution is observable. However, as α becomes larger, it is unlikely that this method will reveal the presence of a nuclear bar. In addition, we note that even with ~1500 sources per field the distribution can be very noisy. Nonetheless, the sample sizes available in surveys in the near future will allow this diagnostic to be employed.

In Fig. 16, we test the impact of extinction on the μ_1,max distribution. The blue markers/curve again show the distribution for the axisymmetric NSD. We then apply the selection function described in Sect. 3.4 (see also Fig. 10) to the NSD and re-compute the μ_ℓ,max − ℓ distribution (red markers/curve). The inset in Fig. 16 shows the density distribution of the NSD when the selection is applied. Extinction has a strong impact on the μ_l,max-ℓ distribution, creating a strong peak which is slightly shifted towards positive longitude. If not properly accounted for, extinction can easily produce a signal mimicking that of a nuclear bar; a major drawback of this diagnostic.

5 Discussion and conclusion

In this study, we have employed N-body simulations to systematically investigate kinematic diagnostics for detecting a nuclear bar in the Galactic centre. The primary objective was to develop kinematic diagnostics that are both theoretically robust and observationally practical given current and future data. The models were constructed to test the difference in kinematic signatures between an axisymmetric NSD and a nuclear bar. We also tested the role of extinction due to the presence of dust and gas in the Galactic centre using a toy model. We assessed three methods for detecting a nuclear bar: (1) Vertex deviation of the velocity ellipse (Sect. 3), (2) the asymmetry in the μ_ℓ − ℓ distribution (see Sect. 4), and (3) the correlation/anti-correlation between the Gauss-Hermite moment h₃ and the line-of-sight velocity, V_los (see Appendix E).

Our study indicates that the vertex deviation, l_v, provides the most robust method of detecting a nuclear bar in the Galactic centre (see Sect. 3.2). An NSD produces a velocity ellipse with its major axis aligned with the Galactic longitude velocity (i.e., l_v ≈ ± 90°), while a nuclear bar deviates from this, producing an l_v value which depends on its angle, α. A nuclear bar with its near side at negative longitudes produces positive l_v, which contrasts sharply with the background negative l_v of the largescale bar. Conversely, a nuclear bar with its near side towards positive longitudes produces a negative l_v. Importantly, contamination of roughly 20-25% from the main bar minimally affects l_v in the central region where the NSD density dominates, which spans roughly |ℓ| < 0.9° and |b| < 0.25°. However, in cases where α = ±90°, the l_v signals of an NSD and NSB are practically identical.

We also explored the strong possibility that the MW’s NSD is not fully axisymmetric, but elliptical, with its major axis orthogonal to the main bar, as is observed in external galaxies (e.g. Bittner et al. 2020). We find that in such cases, l_v is positive (see Sect. 3.6), and mimics the kinematic signature of a nuclear bar viewed nearly side-on with its semi-major axis at negative longitude. Therefore, should a positive l_v in the Galactic Centre be detected in future studies an NSD cannot be easily ruled out.

The uncertainties associated with l_v are a function of sample size and nuclear bar strength. In the weak bar case (W-NSB), stellar velocities are more isotropic (smaller β) (see Fig. 12), which systematically increases the uncertainty. Therefore, the sample size required for a statistically significant measurement of l_v changes with nuclear bar strength. In the strong NSB case (S-NSB), we find that roughly N = 500 is sufficient, whereas this increases to N = 1000 for a weaker NSB. In this context, one of the key advantages of using vertex deviation is the ability to combine stars from different observational fields, increasing the statistical significance of the l_v measurement. This suggests that the data available in the current KMOS survey of the NSD (Fritz et al. 2021) may allow for the detection of a nuclear bar via the vertex deviation diagnostic. This possibility will be addressed in a companion paper.

We find that l_v could be impacted by extinction; however, the effect is not strong. We applied a three-component probabilistic selection function that accounts for distance-dependent detectability, midplane extinction and crowding, and asymmetric dust extinction as a function of Galactic longitude (Appendix D). We find that this shifts l_v towards positive values, possibly mimicking an NSB with its near side towards negative Galactic longitudes.

Similar to San Martin Fernandez et al. (2025) for a largescale bar, we find that while l_v reveals non-axisymmetry, it should not be considered as a measure of the level of non-axisymmetry. Indeed, for a given nuclear bar angle, l_v will take on roughly the same value, regardless of bar strength (see Fig. 12). Instead, the anisotropy, β, of the velocity ellipsoid may be employed in order to characterise the bar amplitude, as it increases with bar amplitude.

The observed asymmetry in the μ_ℓ − ℓ plane provides another method of detecting a NSB. In Fig. 15 we show the distribution of μ_1,max versus ℓ for various orientations (panels) of the NSB (red markers/curves). The NSB will induce deviations from the otherwise flat distribution produced by an axisymmetric NSD (blue points). This detection method is most effective in cases where the NSB angle is small (−20° < α < 20°). This produces a ‘step’ feature in the inner |ℓ| < 1°. For NSB angles |α| > 30°, the asymmetry (and therefore the central slope) becomes weaker, and this method becomes unreliable. In such cases, using the vertex deviation is preferred. Two important caveats are that unlike the vertex deviation, this method is highly sensitive to the asymmetric extinction in the Galactic centre (see Fig. 16). Secondly, while the bootstrapping we applied reflects how μ_ℓ varies depending on sampling, individual sources in observational data will have their own μ_ℓ uncertainties, further diluting the signal.

We also tested the correlation or anti-correlation between the Gauss–Hermite moment h₃ and the line-of-sight velocity V_los. Although this method has proven effective for identifying NSDs in external galaxies (e.g, de Lorenzo-Cáceres et al. 2013), we find no observable distinction between the NSD and NSB in our models: both exhibit similar h₃–V_los anti-correlation for all NSB orientations (see Appendix E). Consequently, this diagnostic does not provide meaningful constraints on the presence of an NSB in the MW, and we therefore include the full analysis only in the Appendix for completeness.

Overall, our findings demonstrate that l_v is both robust against extinction and selection effects and sensitive across a wide range of nuclear-bar orientations. The μ_ℓ–ℓ asymmetry method plays a complementary but limited role, while Gauss–Hermite diagnostics do not provide discriminating power in this context.

Fig. 15 Distribution of μℓ,max versus ℓ distribution for the S-NSB model (red markers and smoothed fit) compared to the axisymmetric NSD (blue). Each panel shows a different nuclear bar angle, α. For small angles |α| ≲ 20°, the NSB produces a ‘step’ feature in the central region (|ℓ| < 1°) that clearly deviates from the flat NSD distribution. This asymmetry weakens substantially as |α| ≳ increases, and becomes indistinguishable from the NSD at α = 0° or ±90°, limiting this diagnostic’s applicability to a narrow range of nuclear bar orientations. Black open circles show measurements from individual midplane fields with N = 1500 stars each from the S-NSB model. Uncertainties are estimated via bootstrap resampling.

Fig. 16 Distribution of μℓ,max versus ℓ for the NSD model without selection effects (blue markers and curve), and with the observational selection function applied (red). The inset shows the resulting face-on density distribution after selection. The distribution is no longer constant, but produces a strong peak in the distribution which could be misinterpreted as evidence of a nuclear bar if observational biases are not carefully accounted for.

6 Future work

The number of stars required to reliably measure the vertex deviation is available in the KMOS survey (Fritz et al. 2021), which contains ≈3000 sources. This possibility shall be explored in a companion paper. Nevertheless, a detailed picture of the NSD region will likely require the larger sample sizes and precise measurements from future surveys such as the KMOS public survey, VVVX-GalCen (Nogueras-Lara et al. in prep). This will provide temperatures, velocities, spectral indices, and metallicities for ~40 000 sources in the NSD region. The MOONS REDdened Milky WAY (REDWAY) survey (Cirasuolo et al. 2020; Gonzalez et al. 2020), will also provide high resolution spectra for ~13 000 stars within the region |ℓ| < 1° and |b| < 0.5°, enabling measurements of radial velocities and metallicities. Similarly, the JASMINE Galactic Centre survey (GCS) (Kawata et al. 2024) will cover a region −0.7° < |ℓ| < 1.4° and |b| < 0.6°, and is expected to provide accurate positions and proper motion measurements for about 1.2 × 10⁵ stars. Lastly, the Roman Galactic Bulge Time Domain Survey (Terry et al. 2023) will provide high-precision photometric and astrometric measurements of ~3.3 million stars over a 0.28 deg² field centred on SgrA*, with proper motion precisions of 2.5–3.5 μas/yr. These datasets will facilitate the application of the diagnostic tools discussed in this paper, advancing our understanding of the Galactic Centre.

Acknowledgements

We thank the anonymous referee for providing valuable feedback which improved the manuscript. KF, MCS, ZF, XL, AV, and MD acknowledge financial support from the European Research Council under the ERC Starting Grant “GalFlow” (grant 101116226). MCS further acknowledges financial support from the Fondazione Cariplo under the grant ERC attrattività no. 2023-3014. JLS acknowledge support from the Royal Society (URF\R1\191555; URF\R\241030). The simulations in this paper were run both at the High Performance Computing Facility of the University of Lancashire and at the DiRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC High Performance Computing (HPC) Facility: www.dirac.ac.uk. This equipment was funded by the Department for Business, Innovation and Skills (BIS) National E-infrastructure capital grant ST/J005673/1, STFC capital grant ST/H008586/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure. FNL gratefully acknowledges financial support from grant PID2024-162148NA-I00, funded by MCIN/AEI/10.13039/501100011033 and the European Regional Development Fund (ERDF) “A way of making Europe”, from the Ramón y Cajal programme (RYC2023-044924-I) funded by MCIN/AEI/10.13039/501100011033 and FSE+, and from the Severo Ochoa grant CEX2021-001131-S, funded by MCIN/AEI/10.13039/501100011033.

References

Appendix A The metal-poor halo

Observations reveal that while the Milky Way’s boxy-peanut bulge is predominantly traced by metal-rich stars, the more spheroidal metal-poor ([Fe/H] < −0.5) population must include a contribution of ~ 10% from the halo (e.g. Ness et al. 2013; Rojas-Arriagada et al. 2020). This metal-poor halo is not included in the large-scale bar simulation employed in this work and could therefore represent an additional source of contamination in our analysis.

To assess the potential impact of this population, we construct a synthetic metal-poor component by sampling velocities from Gaussian distributions with dispersions σ_r = 99 km s⁻¹ and σ_ℓ = 138 km s⁻¹ (Table 1 in Babusiaux et al. (2010)), which characterises the kinematics of metal-poor stars in the bulge. When sampling stars from our models to compute the diagnostic metrics, we add a 10% contribution from this synthetic metal-poor population. We find that the resulting impact on our diagnostics is consistently within the uncertainties determined through bootstrapping (see Sect. 3.1). Given its negligible effect on our results, we exclude this component from our final modelling framework.

Appendix B Vertex deviation derivation

In 2 dimensions, the velocity dispersion tensor defined by Eq. (1) can be written as: $${\sigma }_{ij}=\left(\begin{array}{cc}{\sigma }_{11}^2& {\sigma }_{12}^2\\ {\sigma }_{21}^2& {\sigma }_{22}^2\end{array}\right),$$(B.1)

where 1 refers to the first coordinate (in our case, the line-of-sight direction v_los) and 2 to the second coordinate (in our case, the transverse direction v_ℓ), and the tensor is symmetric, ${\sigma }_{12}^2={\sigma }_{21}^2$. Because it is symmetric, the dispersion tensor is always diagonalisable. The eigenvalues are given by: $${\sigma }_{max/min}=\frac{1}{2}[{\sigma }_{11}^2+{\sigma }_{22}^2\pm D],$$(B.2)

where the + (−) sign corresponds to σ_max (σ_min), and we have defined $$D\equiv \sqrt{{({\sigma }_{11}^2-{\sigma }_{22}^2)}^2+4{\left({\sigma }_{12}^2\right)}^2}.$$(B.3)

The corresponding eigenvectors are: $$D\equiv \sqrt{{({\sigma }_{11}^2-{\sigma }_{22}^2)}^2+4{\left({\sigma }_{12}^2\right)}^2}.$$(B.4)

The vertex deviation, denoted as l_v and illustrated by the black arc in Fig. B.1, is defined as the angle that the eigenvector corresponding to the largest eigenvalue (V_max) makes with the v_los axis. Elementary trigonometry yields: $$\tan (2l_{\mathrm{v}})=\frac{2{\sigma }_{12}^2}{{\sigma }_{11}^2-{\sigma }_{22}^2},$$(B.5)

and therefore: $$\boxed{{\displaystyle l_{\mathrm{v}}=\frac{1}{2}\arctan \left(\frac{2{\sigma }_{12}^2}{{\sigma }_{11}^2-{\sigma }_{22}^2}\right)}}$$(B.6)

The vertex deviation l_v takes values within the range [−90°, 90°] and represents the angle that the semi-major axis of the velocity ellipse forms with the v_los axis as illustrated in Fig. B.1.

Fig. B.1 Schematic illustration of the velocity ellipse. The two eigenvectors (Eq. B.4) of the velocity dispersion tensor (Eq. B.1) are shown as green and red arrows. The vertex deviation (Eq. B.6) is the angle indicated by the black arc.

Note that sometimes a modified formulation is employed in the literature (e.g. Soto et al. 2007; Babusiaux et al. 2010; Debattista et al. 2019; Simion et al. 2021; San Martin Fernandez et al. 2025), which incorporates the absolute value of the denominator: $$l_{\mathrm{v}}=\frac{1}{2}\arctan \left(\frac{2{\sigma }_{12}^2}{|{\sigma }_{11}^2-{\sigma }_{22}^2|}\right).$$(B.7)

In this case, l_v takes values in the interval [−45°, 45°] and is the angle that the semi-major axis of the velocity ellipse makes with the nearest coordinate axis (either v_los or v_ℓ). We do not use this alternative definition in this paper.

Appendix C The SPH model

We employ an N-body+SPH model of an isolated disc galaxy to explore how an elliptical, orthogonal NSD impacts our results. This model is equivalent to model M3_nc_b in Fiteni et al. (2021), and has also been described extensively in Cole et al. (2014) and Debattista et al. (2017). Therefore, we give here only a brief overview of the model. Initially, the model comprises of a corona of hot gas embedded within a dark matter halo with a virial radius of r₂₀₀ ≃ 200 kpc and a virial mass of M₂₀₀ = 10¹² M_⊙. Both the gas and dark matter are represented by 5 × 10⁶ particles each, with softening parameters of ϵ = 50 pc and ϵ = 100 pc for the gas and dark matter respectively. All stars form from the gas, with no stars present at t = 0. The initial radial profile of the gas corona is identical to that of the dark matter halo, but contains only 10% of its mass. Gas particles are assigned a tangential velocity that results in a spin parameter $\lambda \equiv J|E{|}^{1/2}/(G{M}_{\text{vir\hspace{0.17em}}}^{5/2})=0.065$, where J and E denote the total angular momentum and the energy of the gas particles respectively, and G is the gravitational constant. The conditions for star-formation to occur require that the cold gas (T < 15 000 K) density exceeds 0.1 cm⁻³ and is part of a converging flow. Star particles forming from the cooling gas also have a softening parameter of ϵ = 50 pc. The star-formation efficiency is set to 5%, and we employ the feedback recipe described by Stinson et al. (2006). The supernova feedback model injects 10% of the 10⁵¹ erg per supernova back into the interstellar medium. Additionally, we implement turbulent diffusion for gas mixing as outlined by Shen et al. (2010). The model was evolved for 10 Gyr with GASOLINE (Wadsley et al. 2004, 2017), and develops a long-lived primary bar early in its evolution at roughly 3 Gyr. The orthogonal NSD structure forms later and is comprised mostly of young stars with ages τ_⋆ < 4 Gyr.

Appendix D Observational selection effects

To generate a realistic mock observation from the models, we applied a three-component probabilistic selection function that accounts for distance-dependent detectability, midplane extinction and crowding, and asymmetric dust extinction as a function of Galactic longitude. Each star in the model was assigned a detection probability based on these effects, and the final ‘observed’ sample was generated through Monte Carlo sampling. The three components of our selection function are illustrated in Fig. 10.

Due to extinction effects and crowding in the midplane, stars on the far side of the disc (d > 8.2 kpc) are less likely to be detected. We model this selection bias using a skewed normal distribution (Fig. 10a): $$P(d)=\frac{1}{N}\cdot {\varphi }_{\text{skew\hspace{0.17em}}}(d;\mu =8.2\mathrm{k}\mathrm{p}\mathrm{c},\sigma =1\mathrm{k}\mathrm{p}\mathrm{c},\alpha =-1),$$(D.1)

where ϕ_skew is the skewed normal probability density function with location parameter μ = 8.2 kpc, scale parameter σ = 1 kpc, and skewness parameter α = −1. The negative skewness produces an asymmetric distribution that favours stars on the near side of the NSD (d < 8.2 kpc), mimicking the observational reality that nearer stars are more readily detected than those on the far side of the Galactic Centre. The distribution is normalised such that max(p_dist) = 1.

Severe stellar crowding and extinction significantly reduce detection completeness near b = 0°. We also modelled the completeness with latitude (Fig. 10b): $$P(b)=1-0.7\exp (-\frac{|b|}{b_{\text{scale\hspace{0.17em}}}}),$$(D.2)

where b_scale = 0.1° is the characteristic latitude scale. This functional form imposes a maximum 70% incompleteness at the Galactic plane (b = 0°), with completeness recovering exponentially towards higher latitudes. As visible in Fig. 10b, the detection probability drops to only 30% at the midplane and recovers to nearly full completeness by |b| ~ 0.5°.

Observations of the inner Galaxy reveal a strong asymmetry in dust distribution, with approximately two-thirds of the obscuring material concentrated at positive Galactic longitudes. We model extinction using a hyperbolic tangent profile (Fig. 10c): $$P(\ell )=0.5+0.25[1-\tanh \left(\frac{\xi \cdot \ell }{{\ell }_{\text{scale\hspace{0.17em}}}}\right)],$$(D.3)

where the extinction strength coefficient is: $$\xi =\{\begin{array}{ll}3.0& \text{\hspace{0.17em}for\hspace{0.17em}}\ell >0\\ 1.0& \text{\hspace{0.17em}for\hspace{0.17em}}\ell \le 0\end{array}$$(D.4)

and ℓ_scale = 0.5 controls the transition sharpness. This formulation ensures that detection probability decreases monotonically as longitude increases. The total detection probability for each star is computed as the product of the three independent components shown in Fig. 10: $$P_{\text{total}}(d,\ell ,b)=P(d)\times P(b)\times P(\ell ).$$(D.5)

To achieve a target overall completeness C_target, we normalised the probability distribution: $$p_{\text{norm}}(d,\ell ,b)=\frac{p_{\text{total}}(d,\ell ,b)}{⟨p_{\text{total}}⟩}\times C_{\text{target}},$$(D.6)

where ⟨p_total⟩ is the mean probability across all model stars. The final observed sample is generated through Bernoulli trials: each star is retained if a uniform random variate u ~ 𝒰(0, 1) satisfies u < p_norm. For our fiducial model, we adopted C_target = 0.15%, yielding a sparse sample that reflects the extreme observational challenges of detecting individual stars in the crowded, heavily extincted Galactic Centre region. The actual achieved completeness may differ slightly from the target due to the stochastic sampling process and probability clipping.

Appendix E Gauss-Hermite moments

Line-of-sight velocity distributions have been widely used as kinematic diagnostics of galaxy features (e.g. Gerhard 1993; van der Marel & Franx 1993; Chung & Bureau 2004; Bureau & Athanassoula 2005; Debattista et al. 2005; Iannuzzi & Athanassoula 2015; Li et al. 2018; Zakharova et al. 2024; Fraser-McKelvie et al. 2025). Specifically, correlations and anti-correlations between h₃ (‘skewness’) and the mean line-of-sight velocity in edge-on discs provide insight into the underlying stellar orbital structure of different components. A positive h₃–V_los correlation typically indicates elongated orbits (indicative of a bar), whereas an anti-correlation suggests the presence of more circular orbits characteristic of a stellar disc (Bureau & Athanassoula 2005; Iannuzzi & Athanassoula 2015; Li et al. 2018; Zakharova et al. 2024). An anti-correlation is also associated with the presence of an NSD (Chung & Bureau 2004; Zakharova et al. 2024; Fraser-McKelvie et al. 2025).

In Fig. E.1, we plot 2D maps of the line-of-sight velocity (left), the h₃ Gauss-Hermite moment (middle), and the velocity dispersion (right). We show the distributions for the large-scale MW (top), NSD (middle), and S-NSB (α = 90°) models. In the S-NSB case, we tested the diagnostic for several values of α, and found no variation among the results. We only show the results for α = 90° as a representative case. As expected, we see positively correlated h₃ and V_los for the large-scale bar model (top row). Consistent with observations (e.g. Chung & Bureau 2004; Fraser-McKelvie et al. 2025), the presence of the NSD (middle row) produces an anti-correlation between the mean velocity and h₃. In the bottom row, we show V_los and h₃ for the S-NSB model. We again find that these parameters are anti-correlated, indistinguishable from the NSD model case. This is consistent with de Lorenzo-Cáceres et al. (2013), who carried out a detailed kinematical and stellar population analysis in a sample of four (SB0 – SBb) double-barred galaxies using integral-field spectroscopy. Their Figure 4 shows the face-on 2D maps of h₃, and reveals that the distribution is not dependent on the nuclear bar angle, and is thus most likely related to the NSD in which the nuclear bar is embedded. The velocity dispersion (right column in Fig. E.1) also shows no clear observable difference between the NSD and S-NSB cases, save for a greater central dispersion for the S-NSB. Therefore, this method is unlikely to reveal the presence of nuclear bar in the Galactic centre.

Fig. E.1 Left to right: Two-dimensional maps of the line-of-sight velocity, the h3 Gauss-Hermite moment, and the line-of-sight velocity dispersion. We show the distributions for the large-scale MW (top, note different scale), NSD (middle), S-NSB (α = 90°) models.

Appendix F Vertical motion kinematics

The left column of Fig. F.1 shows the μ_b distributions for stars within 7.8 < d [kpc] < 8.4 across various fields for the NSD, S-NSB, and large-scale bar models. The kinematically cooler NSD is characterised by low velocity dispersion, with typical values |μ_b| ≲ 5 mas yr⁻¹. The NSB shows moderately higher dispersion, reaching |μ_b| ≈ 8 mas yr⁻¹, while a subset of large-scale bar stars extends to significantly higher |μ_b| values. Applying a conservative cut of |μ_b| > 5 mas yr⁻¹ across all fields reduces contamination from large-scale bar interlopers by at most ~2%, having negligible impact on our kinematic diagnostics. Nonetheless, observers may wish to apply such cuts when constructing clean samples of NSD stars. The right column of Fig. F.1 shows the distance distribution of stars excluded by the |μ_b| ≲ 5 mas yr⁻¹ cut.

Fig. F.1 Left: Vertical proper motion, μb for several fields. We show the distribution for the isolated NSD (black), isolated S-NSB (blue) and large-scale MW (red) models for stars in the distance range 7.8 < d [kpc] < 8.4. Right: Distance distribution for stars excluded by the vertical cut |μb| < 5 mas/yr.

All Figures

	Fig. 1 Face-on (bottom) and side-on (top) stellar density distribution of the models. Left to right, the axisymmetric NSD, isolated strong NSB, isolated weak NSB, and large-scale bar. The large-scale bar model is oriented at −27° relative to the Solar-GC line (dotted line). The nuclear bar angle, α, is a free parameter. In all cases, the Solar position is (x, y, z) = (−8.2, 0, 0.012) kpc.
In the text

	Fig. 2 Surface density of the large-scale MW model with no nuclear component (left), the isolated NSD (middle), and the ratio between the two (right). The ratio illustrates how prominent the MW’s NSD would be if we were to see it from the Sun (see Fig. 1).
In the text

Fig. 3 Evolution of the NSD model. Left: bar strength, measured by the Fourier A2 parameter, as a function of time. The parameter A2 was calculated considering all stars within the radius Rmax that encloses 98% of the model’s mass. Right: evolution of Rmax with time. The initially stable NSD starts forming a bar at t ≈ 0.5 Gyr, which reaches maximum strength at t ≈ 1.7 Gyr. Black and blue dashed lines indicate the times at which the S-NSB and W-NSB models are extracted, respectively.

In the text

Fig. 4 Plots in galactocentric Cartesian coordinates of the isolated NSD (left), strong NSB (middle), and weak NSB (right) models oriented at −60° relative to the Sun-GC line. Overlaid are the velocity ellipses coloured by the vertex deviation, lv. We have placed an ellipse in each bin with N > 25 stars. The dashed lines reflect a 2° field of view centred on the Solar position at (x, y) = (−8.2, 0) kpc. The NSD shows radially oriented ellipses, while the S-NSB model ellipses show strong alignment with the bar’s major axis. The alignment is weaker for the W-NSB model.

In the text

Fig. 5 Velocity ellipses for the models, with the major and minor axes shown in solid and dashed lines, respectively. We uniformly sampled N = 2400 stars from the region |ℓ| < 0.9° and |b| < 0.25°. In each panel we denote the vertex deviation (lv), the anisotropy (β), the bar angle α (except for the NSD case), and the velocity dispersions in the line of sight and longitude directions. The top-left panel shows the ellipse for the large-scale MW model with no nuclear component. As expected, the direction of highest dispersion is in the line of sight direction. The bottom left shows the ellipse for the axisymmetric NSD model. The other panels show the ellipses for the S-NSB model with various nuclear bar angles. In this case, the direction of highest dispersion depends on the nuclear bar angle.

In the text

Fig. 6 Two-dimensional maps in (ℓ, b) space coloured by the vertex deviation, lv. Top row: maps for the large-scale Milky Way model without a nuclear component (left) and with an axisymmetric NSD model (right). The NSD is indicated by the blue circle in the inset. Middle panels: maps for the S-NSB model with the nuclear bar at various orientations, α. Insets show the geometric configuration of the outer bar (red line) and nuclear bar (blue line). White markers in the top-left panel indicate the fields used in the analysis of Sect. 3.3. Bottom panels: same as above but for the W-NSB model.

In the text

Fig. 7 Vertex deviation, (lv) and anisotropy, (β) as a function of Galactic latitude, |b|, for the NSD (black curves) and S-NSB model with varying values of the nuclear bar angle, α (coloured curves). As expected, the NSD has lv ≈90° in the midplane. In contrast, the midplane values for the S-NSB model are dependent on α. Across all model configurations, both lv and β converge to roughly similar values at |b| ≥ 0.5°. This marks the latitude at which contamination from large-scale bar stars dominate the density distribution.

In the text

Fig. 8 Bootstrapped distributions of lv as a function of sample size. Using the S-NSB (left panels) and W-NSB (right panels) models with a nuclear bar angle α = −60°, we sample stars from fields shown with (ℓ, b) coordinates indicated in each panel. We also indicate the position of each field (white markers) in the top left panel for the S-NSB model in Fig. 6. For each field, we perform B = 1000 bootstrap iterations using sample sizes of 120 (dotted red), 500 (dashed green), and 1000 (solid blue) stars. For the S-NSB model, the mean uncertainty (using Eq. (4)) decreases from ±6.18° for 120 stars to ±2.97° for 500 stars, and ±2.01° for 1000 stars. We note that lv varies marginally across fields, ranging from −66.4° to −62.8° (with a standard deviation of σℓv ± 2.58°;). The uncertainty in lv for the W-NSB case is 16.08°, 9.43° and 6.52° for 120, 500 and 1000 stars per field respectively. Across the fields, we find that the median lv value ranges between −70.4° and −60.1°. This variation likely reflects both the intrinsic differences in stellar kinematics across fields and varying contamination from large-scale bar stars.

In the text

Fig. 9 Uncertainty in the vertex deviation σlv as a function of sample size, N. Stars are sampled from fields in Fig. 8. We show the S-NSB (red) and W-NSB (green) cases, both with α = −60°. The uncertainty decreases steeply up to N = 500, then more gradually thereafter. The W-NSB case yields systematically higher uncertainties (by ~8°) compared to S-NSB at all sample sizes due to its more isotropic velocity distribution (lower β).

In the text

Fig. 10 Left: observational selection functions applied to the models. (a) Distance-dependent detection probability favouring stars on the near side of the disc (d < 8.2 kpc). (b) Midplane extinction reducing completeness to 30% at b = 0°. (c) Asymmetric longitude extinction with steeper decline at positive longitudes. The total probability is the product Ptotal(d, ℓ, b) = P(d) × P(b) × P(ℓ). Right: impact of selection on the NSD model’s distance (top), latitude (middle), and longitude (bottom) distributions, showing the unbiased sample (black) versus the observed sample after applying selection effects (red).

In the text

Fig. 11 Bootstrapped distribution of $l_v^{\ast }$ for the axisymmetric NSD model without (blue) and with (red) the observational selection function applied. Each distribution results from 10 000 bootstrap iterations sampling N = 12500 stars from the region |ℓ| < 0.9° and |b| < 0.25°. The unbiased NSD yields lv = (−87.4 ± 5.2)°, consistent with an axisymmetric disc. However, applying realistic selection effects shifts the measurement to lv = (82.2 ± 5.1)° (towards positive values), demonstrating that observational biases could mimic a nuclear bar with its near side at negative longitudes if not properly corrected.

In the text

Fig. 12 Vertex deviation (left) and anisotropy (right) of 1800 stars sampled from the region |ℓ| < 0.9° and |b| < 0.25° as a function of the nuclear bar angle, α. Red curves show S-NSB results; green curves show W-NSB. While lv remains roughly identical between strong and weak nuclear bars, β decreases by approximately a factor of two for the weaker bar, approaching NSD-like values. This demonstrates that lv is a robust orientation diagnostic but a blunt probe of bar amplitude, while β may serve as a complementary bar amplitude indicator.

In the text

Fig. 13 Synthetic RGB image produced with PYNBODY (Pontzen et al. 2013) of the stellar disc of the N-body+SPH model. We show the model face-on with the main bar aligned roughly −27° to the Sun-GC line and with the Sun located at (x, y) = (−8.2, 0) kpc. The central region of this model hosts an orthogonal sub ~1 kpc NSD with its near side towards negative Galactic longitude. This configuration tests whether orthogonal NSDs can be distinguished from true nuclear bars using kinematic diagnostics.

In the text

Fig. 14 Velocity ellipses for the quasi-axisymmetric orthogonal NSD model, constructed from stars in the region |ℓ| < 0.9° and |b| < 0.25°. Left: all stars show lv ≈ −89.5° and low anisotropy (β = 0.22), appearing nearly axisymmetric due to contamination from older large-scale bar stars. Right: result when selecting younger stars (τ⋆ < 5 Gyr) that predominantly trace the late-forming NSD yields lv≈79.2° with higher anisotropy β = 0.62. This demonstrates that an orthogonal NSD can produce positive lv, similar to a nuclear bar with its near side at negative longitudes.

In the text

Fig. 15 Distribution of μℓ,max versus ℓ distribution for the S-NSB model (red markers and smoothed fit) compared to the axisymmetric NSD (blue). Each panel shows a different nuclear bar angle, α. For small angles |α| ≲ 20°, the NSB produces a ‘step’ feature in the central region (|ℓ| < 1°) that clearly deviates from the flat NSD distribution. This asymmetry weakens substantially as |α| ≳ increases, and becomes indistinguishable from the NSD at α = 0° or ±90°, limiting this diagnostic’s applicability to a narrow range of nuclear bar orientations. Black open circles show measurements from individual midplane fields with N = 1500 stars each from the S-NSB model. Uncertainties are estimated via bootstrap resampling.

In the text

Fig. 16 Distribution of μℓ,max versus ℓ for the NSD model without selection effects (blue markers and curve), and with the observational selection function applied (red). The inset shows the resulting face-on density distribution after selection. The distribution is no longer constant, but produces a strong peak in the distribution which could be misinterpreted as evidence of a nuclear bar if observational biases are not carefully accounted for.

In the text

	Fig. B.1 Schematic illustration of the velocity ellipse. The two eigenvectors (Eq. B.4) of the velocity dispersion tensor (Eq. B.1) are shown as green and red arrows. The vertex deviation (Eq. B.6) is the angle indicated by the black arc.
In the text

	Fig. E.1 Left to right: Two-dimensional maps of the line-of-sight velocity, the h3 Gauss-Hermite moment, and the line-of-sight velocity dispersion. We show the distributions for the large-scale MW (top, note different scale), NSD (middle), S-NSB (α = 90°) models.
In the text

	Fig. F.1 Left: Vertical proper motion, μb for several fields. We show the distribution for the isolated NSD (black), isolated S-NSB (blue) and large-scale MW (red) models for stars in the distance range 7.8 < d [kpc] < 8.4. Right: Distance distribution for stars excluded by the vertical cut |μb| < 5 mas/yr.
In the text