4.1. Disc and Halo Models

In this work, we use the method described by ^{Sofue (2012)}, where the total rotation velocity at radius 𝑅 is defined as:

where 𝑉_𝑡 is the total rotation velocity and 𝑉_𝑏,𝑉_𝑑,𝑉_ℎ,𝑉_𝐻𝐼 are the rotation velocity components for the bulge, disc, DM halo and HI disc, respectively. For the present study, since our sample is composed of late spiral and irregular galaxies, we will consider the velocity contribution of a bulge component as negligible.

The stellar disc velocity component, 𝑉_𝑑, was determined with the analytical expression obtained from the Poisson equation for a thin disc following an exponential surface density profile and a constant M/Lalong the radius (^{Freeman 1970}; ^{Binney & Tremaine 1987}).

where 𝑀 _𝑑 is the total stellar disc mass, 𝑎 _𝑑 the exponential scale radius and 𝐷 (𝑋) is defined as:

being 𝐼 _𝑖 y 𝐾 _𝑖 the modified Bessel functions and X =Rαd . The two free parameters are the total stellar disc mass, 𝑀 _𝑑 , and the exponential scale radius, 𝑎 _𝑑 . In the present investigation, the DM halo was modelled assuming three different density profiles: Navarro-Frenk-White (hereafter, NFW) (^{Navarro et al. 1997}), pseudo-isothermal (ISO) (^{Kent 1986}), and Burkert-type (BURK) (^{Burkert 1995}). They are described as:

where 𝜌 ₀ is the core density and ℎ is the scale length of the DM halo. Considering spherical symmetry, the total halo mass is obtained by integrating along the radius, giving the following expressions for the NFW, ISO and BURK profiles, respectively:

The rotation velocity component for the dark matter halo at radius 𝑅 is:

DM density profiles are often expressed in terms of the rotation velocity at the virial radius, 𝑉 ₂₀₀, and the concentration parameter, 𝑐, defined as:

where 𝑅 ₂₀₀ is the radius of a spherical halo with a mass equal to 𝑀 ₂₀₀ and the virial mass 𝑀 ₂₀₀ is defined as the total mass contained within a sphere whose mean density equals 200 times the critical density of the Universe (𝜌 _𝑐 ) at 𝑧 = 0

We have considered a Hubble constant value of 𝐻 ₀ = 73 𝑘𝑚 𝑠 ⁻¹ 𝑀𝑝𝑐 ⁻¹ (^{Riess et al. 2022}) throughout this investigation. This value is the same as the one assumed in the published fits used as comparison (^{Li et al. 2020}).

4.2. HI Disc Velocity Component

Since all the elements in our sample are classified as late-type galaxies (mostly Sm and Irr), it should be necessary to take into account the dynamical effects of a separate HI disc component because such galaxies tend to have an important fraction of their total baryonic mass in the form of HI (^{Roberts 1962}; ^{Huchtmeier & Richter 1989}; ^{Korsagaetal.2019}). In this work, the HI disc velocity component is obtained from the SPARC database HI mass model for each galaxy. This component tis calculated with the ^{Casertano (1983)} formula. In this investigation, in order to isolate the contributions of the stellar disc and DM halo, the HI disc velocity component for each point in the RC is subtracted from the observed rotation velocity using the following expression:

where 𝑉 _𝑜𝑏𝑠 is the observed rotation velocity , 𝑉 _𝐻𝐼 the HI disc velocity component and 𝑉 _{𝑐𝑜𝑟𝑟} the corrected velocity. This latter component is then inserted in the fitting algorithm and analysed.

4.3. Photometry and (𝑀/𝐿) _★ Ratio

In all the problems that involve the fitting of observational data using a certain theoretical model it is important to select a suitable set of constraints to ensure that our algorithm does not give unrealistic results. One important constraint in RC analysis is the surface brightness photometry (^{Katz et al. 2017}; ^{Korsaga et al. 2019}; ^{Li et al. 2020}). The surface brightness profile of the galaxy provides information about how the stars are distributed throughout the disc. In recent years, many studies have shown that near-infrared (NIR) surface photometry can be used as an excellent stellar disc mass tracer, since these bands are less susceptible to systematic uncertainties caused by dust emission that afflict the optical bands (^{Kennicutt et al. 2003}; ^{Leroy et al. 2008}; ^{Meidt et al. 2014}). Another reason is that most of the stellar disc total mass is composed of low-mass stars, whose main emission is at NIR bands (^{Kroupa & Jerabkova 2021}). This is seen as a low dispersion on the (𝑀/𝐿)-color correlation at NIR bands compared to the optical bands (^{Bell & de Jong 2001}). In this work, we used the exponential disc scale length (𝑎 _𝑑 ) and the 3.6 𝜇m total luminosity (LT3.6μm), extracted from the SPARC database (Table 1). It is important to note that, although the thin disk model (^{Freeman 1970}; ^{Binney & Tremaine 1987}) represents a simplification with respect to models used in other studies, it has the advantage of only depending on two photometric parameters (𝐿 _𝑇 , 𝑎 _𝑑 ) to estimate the stellar disk contribution, instead of requiring a complete luminosity profile obtained with high spatial resolution observations (^{Casertano 1983}). In the case of low-mass galaxies, the exponential profile is usually a good approximation for a large interval of radii (^{Lelli et al. 2016}b). However, it is common to find small deviations on the disc caused by secondary structures, such as bars or spiral arms.

Another important constraint in our model is the (𝑀/𝐿) _★ ratio for the stellar disc. As mentioned in the Introduction, the uncertainty in this parameter is mainly responsible for the disc-halo degeneracy. A way to counter this degeneracy is by considering the maximal disc hypothesis, where the contribution of the disc mass to the RC is considered to be as large as the RC itself allows it. Qualitatively, it is defined as a disc whose contribution to the total potential exceeds the contribution of the dark matter halo in the inner regions. However, in recent years evidence was found suggesting that the application of the maximum disc hypothesis usually resulted in (𝑀/𝐿) _★ values that were too high compared to those predicted by the stellar population synthesis (SPS) models (^{McGaugh & Schombert 2014}; ^{Schombert & McGaugh 2014}). This effect is more noticeable in late type galaxies (Sc, Sd, Irr), where the contribution of a bulge is small or nonexistent. Thus, when a plausible range for (𝑀/𝐿) _★ is considered, the resulting fits show DM halos whose contributions to the total rotation velocity are dominant at all radii. The discs that show this property are classified as submaximal (^{Sackett 1997}).

For this first implementation and testing of our procedure, we decided to fix the (𝑀/𝐿) _★ ratio instead of considering it as a free parameter. This was done in order to avoid adding another free parameter to the algorithm, hampering the interpretation of our results. We fixed this ratio at (𝑀/𝐿) _★ = 0.5, constant along the disc, which many studies have converged to take as the fiducial valueat 3.6 𝜇m (^{McGaugh & Schombert 2014}; ^{Schombert & McGaugh 2014};^{Lellietal.2016a}), consistent with the results of stellar population synthesis (SPS) models (^{Bruzual & Charlot 2003}; ^{Meidt et al. 2014}).

4.4. AGA Selection Criterion, the χR2 Function

As we briefly mentioned in the Introduction, in the present work we use the χR2 function as the selection criterion for AGA, which is defined as:

where 𝑁 is the number of RC data points, 𝑃 the number of free parameters, Viobs the observed rotation velocity, Vit the calculated total rotation velocity and 𝜎 _𝑖 the observational uncertainties. After establishing the photometry constraints, the remaining free parameters are the ones corresponding to the DM halo (𝜌 ₀ and ℎ), thus 𝑃 = 2.

4.5. Uncertainty Estimation

In order to estimate the uncertainty range of our results, we use the methodology described in the AGA main paper (^{Cantó et al. 2009}), based in the generation of synthetic data sets. The process is quite simple: for each data point of the measured curve, a random data point is generated inside the observational velocity uncertainty interval 2𝜎. The procedure is repeated until a set of synthetic RCs is obtained. These synthetic curves serve as statistical equivalents on the uncertainty estimation. Each one is treated and analysed in the same way as the measured RCs. The formula used for the generation of synthetic RCs is:

where 𝑉 ^′ (𝑅) is the synthetic rotation velocity value, 𝑉 _𝑜𝑏𝑠 (𝑅) the measured rotation velocity, 𝜎 the observational error and 𝜉 is a random variable in the interval [0, 1]. In the present work, the uncertainties were estimated using sets with an average of 6 synthetic RCs per galaxy. This number was chosen after some tests designed to find the best compromise between time and accuracy. The uncertainties are defined as the standard deviation of the parameter distribution for the synthetic data set.

4.6. Python Subroutine

Adapting AGA to our procedure required to implement some modifications in the code. As mentioned in section 4, the stellar disc contribution is calculated from equation (2), which contains the first and second order modified Bessel functions. The language in which AGA is compiled (FORTRAN90) does not support these functions implicitly, so it was necessary to add them as an external module. These functions were adapted from the FORTRAN special function package, SPECFUN (^{Cody 1993}). Individual tests were performed for each function to ensure they worked correctly.

The fitting procedure was compiled as a Python subroutine, whose main purpose was to provide AGA with all the information needed from a set of machine-readable tables containing the SPARC data for each galaxy in the sample. These data consisted of the measured RC, the disc exponential scale length (𝑎 _𝑑 ), the galaxy total luminosity at 3.6 𝜇m (L⨀3.6μm) and the HI velocity component. Once a fit was obtained for the specified RC, the subroutine then estimated the corresponding 𝑐, 𝑉 ₂₀₀, 𝑀 ₂₀₀ and 𝑅 ₂₀₀ values by using the numerical solver included in the SYMPY package (^{Meurer et al. 2017}). In order to guarantee a good coverage of the parameter space, we established a set of loose constraints for our fitting variables. In the case of the DM profile core density, 𝜌 ₀, we set a range between 0 ≤ 𝜌 ₀ ≤ 1 𝑀 _⊙ pc⁻³, while the scale parameter ℎ could vary between 0 ≤ ℎ ≤ 100 kpc. The typical values of 𝜌 ₀ and ℎ tend to be of the order of 𝜌 ₀ ≈ 10⁻² 𝑀 _⊙ pc⁻³ and ℎ ≈ 10³ − 10⁴ pc (^{Sofue 2012}, ²⁰¹³).

Each AGA run consisted of 150 generations, containing 1000 individual points in the parameter space (200 parent + 800 offspring), with a box size factor of 0.5 between generations. The box size factor reduces the parameter-space volume from which new points can be generated around each individual between one generation and the next. Every AGA run starts with a box size spanning the full specified parameter range. The result is then stored and used as a starting point for the next run. This procedure may seem redundant at first, but it is important to note that for each run, the sizes of the sampling boxes are reset to their initial values, so each run starts searching for the solution using boxes of the same size as the ones used in the first run, but centered on improved initial values (^{Canto et al. 2009}). This enhances the accuracy of´ the algorithm and reduces the probability of converging on a local minimum. Finally, after 6 consecutive runs, the results are stored as the best-fit parameters.

5.1. Comparison with Previous Results

Figures 1 and 2 compare our obtained best-fit parameter values against their corresponding values found in the literature (^{Li et al. 2020}) assuming each DM density profile for the parameters 𝑐 and𝑉 ₂₀₀, respectively. A qualitative analysis shows that most galaxies in these diagrams lie near the red 𝑦 = 𝑥 line, within uncertainties, implying a great degree of consistency between studies. However, some deviations were found (e.g. Figure 2, top), which will be addressed below.

Fig. 1 Comparison diagrams of the best-fit 𝑐 values obtained with our procedure against their corresponding published fits in the SPARC database (^{Li et al. 2020}). The 𝑦 = 𝑥 line is indicated in red. The corresponding residual graphs, described in terms of the deviation parameter 𝛿 _𝑐 (Equation 17) are located below each comparison diagram. On all DM profiles, most galaxies in our sample (≳ 90%) lie within 𝛿 _𝑐 ≤ 1. The color figure can be viewed online. 

Fig. 2 Comparison diagrams of the best-fit𝑉200 values obtained with our procedure against their corresponding published fits in the SPARC database (^{Li et al. 2020}). The 𝑦 = 𝑥 line is indicated in red. The corresponding residual graphs, described in terms of the deviation parameter 𝛿 _𝑉200 (equation 18) are located below each comparison diagram. On all DM profiles, most galaxies in our sample (≳ 90%) lie within 𝛿 _𝑉200 ≤ 1. The color figure can be viewed online. 

In order to verify in a quantitative manner if our results are consistent with those found in the literature, we compared our obtained best-fit𝑐 and𝑉 ₂₀₀ samples against their corresponding published values through a set of paired 𝑡tests. According to these tests, for a sample size of 100 galaxies with a 95% confidence level, the critical 𝑡 value, 𝑇 (from which we can reject the same-mean hypothesis) is equal to 1.98. If the absolute values of 𝑡 _𝑐 or 𝑡 _𝑉200 are higher than𝑇, we can consider that the difference between the means of our compared samples is statistically significant. The 𝑡-test results, along with their corresponding Pearson correlation values for the parameters 𝑐 (𝑟 _𝑐 ) and 𝑉 ₂₀₀ (𝑟 _𝑉200 ) are condensed in Table 2.

The 𝑡-test analysis shows that our fits pass the equal mean hypothesis in 5 out of the 6 possible tests (2 fit parameters, 3 DM profiles). The only 𝑡-test that did not pass corresponds to the comparison between our obtained vs the published V200NFW samples. This is seen in the V200NFW comparison diagram (Figure 2, top) as a divergence from the 𝑦 = 𝑥 line at V200NFW ≳ 200 km s⁻¹ . It is interesting to note that most of the galaxies that deviate from the y = x line show ℎ values that lie near the upper limit imposed for our procedure (ℎ ≈ 100 kpc), indicating that the global minimum probably did not fall inside our specified ℎ range. In order to explore this, several runs were performed assuming different upper range values for ℎ. We found that, although the χR2 global minimum was reached, the resulting parameter values would not be physical (i.e. 𝑉 ₂₀₀ > 500 km s−1) (^{Li et al. 2020}). This behaviour was found exclusively when a NFW profile was considered.

Along with the 𝑡-tests, to quantify the deviation between the obtained and the published parameter values for each individual galaxy, we define a deviation parameter (𝛿 _𝑐 and 𝛿 _𝑉200 for 𝑐 and 𝑉 ₂₀₀, respectively) by using the following expressions:

where 𝑐 ^𝑏𝑖𝑏 and V200bib are the published 𝑐 and 𝑉 ₂₀₀ values from SPARC, respectively. The number of galaxies in our sample with 𝛿 _𝑐 and 𝛿 _𝑉200 below different thresholds is described in Table 3.

Table 3 shows that in the case of the ISO and BURK profiles the fraction of galaxies in our sample with |𝛿 _𝑐 | ^, |𝛿 _𝑉200 | ≤ 1 is 95% and 94%, respectively, while the NFW profile only reached 85%. These results indicate that the core-type profiles, such as ISO and BURK, are closer to the respective published values than those obtained assuming a NFW profile.

Regarding the Pearson correlation (𝑟 _𝑐 ,𝑟 _𝑉200 ) between datasets (Table 2), high values (𝑟 _𝑐 ,𝑟 _𝑉200 ≈ 0.8 − 0.94) are observed in 4 out of the 6 comparison diagrams (Figures 1 and 2), implying a good compliance of our results with the literature. In contrast, two of these cases show low correlation values (𝑟 _𝑐 ,𝑟 _𝑉200 ≈ 0.3), specifically in the 𝑐 _𝑁𝐹𝑊 comparison diagram (Figure 1, top) and the V200BURK diagram (Figure 2, bottom). These low values are consequence of a handful of outlying galaxies with exceptionally high 𝛿 _𝑐 or 𝛿 _𝑉200 : one in the 𝑐 _𝑁𝐹𝑊 diagram (F571-8) and three in the V200BURK diagram (NGC0247, UGC11557 and D631-7). The main cause for the observed inconsistencies in these galaxies may be attributed to multiple factors, e.g. a divergence from a pure exponential profile in the stellar disk produced by secondary features (i.e. a bar or bulge), or artifacts produced by the use of different numerical solver routines between studies. Without taking into account the outlying galaxies, the correlation values between the samples are 𝑟 _𝑐 = 0.93 and 𝑟 _𝑉200 = 0.92 for the 𝑐 _𝑁𝐹𝑊 and V200BURK comparison diagrams, respectively. It is important to emphasize that these galaxies represent a small fraction of the total sample. Most of the galaxies analysed with our procedure (≳ 85%), show consistency with the literature. The influence of other factors, such as the assumed stellar disc (M/L)* value and the consideration of a separate HI disc component in the resulting mass profile fits, will be explored in a future work.

5.2. Comparison Between Our Results

A comparison of our best-fit parameter samples assuming different DM models serves as an internal consistency check for our procedure. Figures 3 and 4 show all the comparison diagrams (ISO vs NFW, BURK vs NFW and BURK vs ISO) between our best-fit parameter samples 𝑐 and 𝑉 ₂₀₀, respectively. A qualitative analysis shows that all comparison diagrams display a positive trend, with varying degrees of dispersion. In the case of the 𝑐 parameter (Figure 3), we found high correlation values (𝑟 ≈ 0.9) in all comparisons, displaying a nearly linear trend. For the case of 𝑉 ₂₀₀, the diagrams show lower correlation values (0.43 ≤ 𝑟 ≤ 0.71) and in general, a greater dispersion. This could imply a difference in the sensibility between fitting parameters. From this comparison, we conclude that in general our results show consistency between different models.

Fig. 3 Comparison diagrams between our obtained best-fit 𝑐 samples assuming different DM profiles. 

Fig. 4 Comparison diagrams between our obtained best-fit𝑉200 samples assuming different DM profiles.