3D-Barolo: 3D Galaxy Kinematics Algorithm

• 3D-Barolo is a computational algorithm that fits 3D models to galaxy spectroscopic data-cubes, accurately deriving rotation curves and intrinsic velocity dispersions.
• It employs tilted-ring modeling and explicit convolution with instrumental responses to mitigate beam-smearing and low signal-to-noise effects.
• The algorithm automates source detection and parameter initialization, enabling robust analysis of large survey datasets from instruments like ALMA and SKA pathfinders.

3D-Barolo is a computational algorithm designed for extracting the rotation curves and intrinsic velocity dispersions of galaxies from emission-line spectroscopic data-cubes. Unlike classical 2D techniques relying on velocity fields derived from first moment or Gaussian fits, 3D-Barolo models the observational data in its native three-dimensional format (two spatial, one spectral dimension) using tilted-ring modeling directly in 3D space. This approach allows accurate recovery of kinematic parameters in environments affected by beam-smearing and in galaxies with minimal spatial resolution or low signal-to-noise ratios. The algorithm incorporates explicit convolution with instrumental spatial and spectral response and employs automated modules for source detection and initial parameter estimation, facilitating the analysis of data from instruments including ALMA, SKA pathfinders, and new-generation integral field units (Teodoro et al., 2015).

1. Motivations and Theoretical Foundations

Traditional 2D fitting methods extract velocity fields by collapsing data-cubes to intensity-weighted or best-fit velocity maps, then apply a tilted-ring analysis. While computationally efficient (fast analytic solutions), these approaches are limited by beam-smearing, which artificially flattens velocity gradients and induces spurious velocity dispersions. Moreover, ambiguities arise in velocity assignment for cases with multi-component emission or thick, asymmetric discs and high-inclination systems.

3D-Barolo’s methodology circumvents these limitations by fitting synthetic 3D models to the full data-cube, reproducing instrumental effects via convolution. The key advantage is direct modeling of beam-smearing and the disentanglement of rotation kinematics from intrinsic velocity dispersion. This enables recovery of accurate rotation curves $V_{\rm rot}(R)$ and intrinsic dispersions $\sigma_{\text{int}}(R)$ in scenarios where traditional methods fail—specifically, for galaxies barely resolved across two beams and for data with global S/N as low as 2–3 (Teodoro et al., 2015).

2. Algorithmic Architecture and Workflow

3D-Barolo operates through a multi-stage pipeline, starting with input data and proceeding through parameter estimation, modeling, convolution, fitting, and regularization:

1. Preprocessing and Mask Generation: Data-cubes in FITS format are ingested, and a mask is constructed via a Duchamp-derived source finder, typically by thresholding (>3σ post-smoothing) to isolate emission regions.
2. Parameter Initialization: The algorithm determines initial guesses for key parameters by analyzing the masked data:
   • Center: flux-weighted centroid.
   • Systemic velocity: midpoint between the 20% peaks of the spectral line profile.
   • Position angle and inclination: extracted from velocity gradients and axisymmetry fits.
   • Initial $V_{\rm rot}$ and $\sigma_v$: using analytic relations and/or defaults.
3. Tilted-Ring Model Construction: The disc is divided into rings; for each, Monte Carlo sampling generates clouds with coordinates $(R_c, \theta_c, z_c)$ according to chosen vertical profiles (Gaussian, $\text{sech}^2$, exponential). Sky-plane projection uses inclination $i$ and PA. Line-of-sight velocities are assigned as

$$v_{\rm LOS} = V_{\rm sys} + V_{\rm rot}(R) \sin i \cos \phi + v_{\rm random}$$

with $$v_{\rm random} \sim \mathcal{N}(0, \sigma_{\text{total}}^2)$$ and $$\sigma_{\text{total}}^2 = \sigma_v(R)^2 + \sigma_{\text{instr}}^2$$.

1. Instrumental Convolution: Model cubes are convolved with the instrument’s PSF (2D Gaussian, FFT-based) and spectral response (incorporated via $\sigma_{\text{instr}}$ for the velocity channels).
2. Residual Minimization: A set of objective functions (χ²-like, absolute, relative) measures voxel-wise discrepancies between model and data, weighted by azimuth to emphasize the major axis. Nelder–Mead simplex optimization updates ring parameters so as to minimize the global residual.
3. Ring-by-Ring and Multi-Step Strategy: Rings are fitted sequentially (outwards from $R_{\min}$), freezing inner rings’ parameters after convergence. A second pass regularizes geometric parameters ($i(R)$, PA(R), $V_{\rm sys}$, center) via polynomial or Bézier smoothing, refitting $V_{\rm rot}$ and $\sigma_v$ only.
4. Error Estimation: Monte Carlo routines explore the parameter space around the best fit, quantifying uncertainties via the criterion that the residual increases by a fixed percentage (default 5%, approximating $1\sigma$).

3. Mathematical Formalism

3D-Barolo’s modeling is governed by the following formal equations:

• Observed Intensity:

$$I_{\rm obs}(x, y, v) = \iiint I_{\rm em}(R, \phi, v') \cdot \text{PSF}(x-x', y-y') \cdot \text{LSF}(v-v') \, d^2x' \, dv'$$

• Line-of-Sight Velocity:

$$v_{\rm LOS}(R, \phi) = V_{\rm sys} + V_{\rm rot}(R) \sin i \cos \phi$$

• Velocity Dispersion:

$$\sigma_{\text{total}}^2 = \sigma_{\text{int}}^2 + \sigma_{\text{instr}}^2$$

• Residual Functions:

$$\Delta r = \begin{cases} \frac{(M - D)^2}{\sqrt{D}} & \text{(χ²-like)} \ M - D & \text{(absolute)} \ \frac{M - D}{M + D} & \text{(relative)} \end{cases}$$

Here, $M$ and $D$ denote model and data voxel intensities, respectively.

4. Application Scope and Performance Validation

Empirical assessments demonstrate that 3D-Barolo robustly recovers rotation curves and velocity dispersions across a range of resolutions:

• High-Resolution Cubes (e.g., THINGS): Output agrees with 2D fits; runs complete in approximately one day for cubes of size $1024 \times 1024 \times 87$ on dual-core laptops.
• Mid/Low Resolution (WHISP dwarfs, 30″ beam, 15″ ring width): Success rate up to ∼94%, typical runtime under one minute per galaxy; results are consistent with classical tilted-ring analyses.
• Very Low Resolution (single-dish, ∼2–3 beams across): 2D methods yield inaccurate velocity fields and inflated dispersions; in contrast, 3D-Barolo maintains physical accuracy for $V_{\rm rot}(R)$ and $\sigma_{\text{int}}(R)$ when S/N exceeds 2–3.
• Mock Testing: Reliable kinematic extraction down to two resolution elements; inclination accurate for $45^\circ \lesssim i \lesssim 75^\circ$; robust for spectral channel counts of $\gtrsim 8$ per line profile.

This suggests that, even under severe observational constraints, the three-dimensional approach mitigates biases inherent to velocity field extraction (Teodoro et al., 2015).

5. Data Requirements, Automation, and Practical Considerations

3D-Barolo’s applicability presupposes certain minimum data characteristics:

• Spatial: $\geq 2$ beams across the galaxy (≥2 resolution elements per spatial side).
• Spectral: $\geq 8$–12 channels sampled across the emission profile.
• Signal-to-Noise: S/N$\geq$2–3 per ring post-masking.

Ring width $\Delta R \gtrsim 0.5$ FWHM$_{\rm beam}$ is advised to balance independence against noise fitting. Automation features allow batch processing for large survey datasets (ASKAP, MeerKAT, APERTIF, ALMA, VLT-IFUs). The Duchamp-derived finder, mask smoothing, and thresholding facilitate robust source segmentation. Outputs such as position–velocity diagrams support post-fit diagnostics; manual fixes are occasionally needed for uncertain inclination or center estimation, especially when fits diverge ring-to-ring.

A plausible implication is that these design features position 3D-Barolo as a scalable tool for heterogeneous survey data, with broad utility across local and high-redshift galaxy kinematics (Teodoro et al., 2015).

6. Limitations and Regularization Strategies

Parameter estimation in 3D-Barolo may be affected by ring-to-ring jitter due to noise and fit degeneracies. Regularization is implemented as a two-step smoothing procedure:

1. Initial Fit: All parameters are optimized independently for each ring.
2. Smoothed Geometry: Low-order polynomials or Bézier curves are fitted to geometric parameters (inclination, position angle, systemic velocity, center), then $V_{\rm rot}(R)$ and $\sigma_v(R)$ are refitted with geometry terms held fixed.

Uncertainty quantification is achieved by randomly perturbing parameters from their best-fit values, yielding error bars where the residual function increases by the adopted threshold.

This structured regularization is crucial for large datasets and for cases with uncertain geometric parameters due to low resolution or marginal S/N, ensuring physical consistency in kinematic profiles extracted by the algorithm (Teodoro et al., 2015).