Cosmicflows-4: The Baryonic Tully-Fisher Relation Providing $\sim 10,000$ Distances

The BTFR invokes a correlation between the combined contributions of the stars and gas in a galaxy with rotation rate, a measure of the total mass of the system. The observational extension on the TFR is the incorporation of H i fluxes. The sources of this information are discussed in §2.1. Issues regarding the inference of stellar mass from observed luminosities are discussed in §2.2. Luminosities are acquired at both optical and infrared bands and consistency is required in mass inferences across bands. The formulations derived in this section are tentative. Refined relations are given by the Bayesian analysis discussed in §2.6 and §2.7. The dominant source of bias in the BTFR procedure arises from a faintness detection limit on H i fluxes. The nature of this bias is discussed in §2.3. The impact of this bias on the summation of gas and stellar contributions to total baryonic mass are discussed in §2.4. Our treatment of the bias is validated through tests with a mock catalog in §2.5. With the general properties of the BTFR now established, we carry out a Bayesian analysis to refine parameters, first at optical bands in §2.6, and then the coupling to the infrared in §2.7.

The core product of §2 is the formulation of the BTFR across three optical bands and one infrared band with an arbitrary absolute scale set by assuming Hubble expansion with $H_{0}=75$ km s^-1 Mpc^-1. In §3, the scale is shifted (slightly) to match constraints imposed by BTFR values for samples of galaxies with Cepheid and Tip of the Red Giant Branch distance measurements. Our method for determining distances to individual galaxies is discussed in §4, with uncertainties the topic of §4.1. The message in this latter section is that the sum of error estimates on observed and derivative parameters does not fully explain the BTFR scatter; ie, an unexplained or intrinsic component contributing to scatter is significant.

Given distance estimates for the $\sim 10,000$ galaxies in our sample, in §5 we derive a best value for $H_{0}$ consistent with our data. Various tests are carried out in §6 that the sample must reasonably pass if there is to be confidence in the relative (not absolute) distances that are derived. A tabulation of our results is presented in §7. The BTFR distances derived here are a component of the ensemble of distances that will be assembled in an upcoming Cosmicflows-4 publication. The absolute scale of the BTFR distances and associated determination of $H_{0}$ given here should be considered preliminary.

The mass in neutral Hydrogen gas is (Haynes & Giovanelli, 1984)

where $F_{21}$ is the integrated flux in the 21cm H i signal in Jy ${\rm km~{}s^{-1}}$ and $d$ is the distance of the target in Mpc. Taking the factor $K_{g}=1.33$ in Equation 2 accounts for the contribution from interstellar Helium that must accompany Hydrogen. In the literature (Bradford et al., 2016) the value $K_{g}=1.4$ has been used to account for molecular, ionized, and metal enriched constituents of interstellar gas. We have evaluated results based on these alternative gas multipliers and not found a significant effect. We assume $K_{g}=1.33$.

Values for the integrated H i flux $F_{21}$ are drawn from two sources. One is the All Digital H i file posted in the Extragalactic Distance Database¹¹1https://edd.ifa.hawaii.edu (Courtois et al., 2009; Dupuy et al., 2021). This compilation for more than 20,000 galaxies draws on observations from the Green Bank 100-meter, 300-foot, and 140-foot telescopes, the Arecibo Telescope, the Parkes Telescope, the Nancay Telescope, and the Effelsberg Telescope. All observations summarized in this file, whether originating from our own observations or from others, have been analyzed in a coherent fashion by drawing raw material from the various archives. The other source is the ALFALFA 100% compilation (Haynes et al., 2018), available from the Arecibo Legacy Fast ALFA Survey website²²2http://egg.astro.cornell.edu or at the Extragalactic Distance Database.³³3EDD fully includes material from the ALFALFA 40% compilation (Haynes et al., 2011) analyzed in a consistent way with data from other sources. An unweighted average is taken if fluxes are available from both sources. In all cases, the observations are made with single beam radio telescopes. In the vast majority of cases, the beams are substantially larger than the optical images of the targeted galaxies. Flux could be lost with particularly extended gas-rich galaxies or with very nearby galaxies of large angular extent.

The relationship between stellar mass and luminosity has been extensively studied at infrared bands where flux is dominated by evolved stars on the red giant branch and optimal correlations are anticipated between light and mass. Particular attention has been given to the L-band where observations have been carried out with the space telescopes Spitzer ([3.6] filter) and WISE (adjacent W1 filter).

Stellar population synthesis or star formation history models can inform the value and scatter of the parameter $\Upsilon^{*}=M^{*}/L$ at the closely similar [3.6] and W1 bands (Eskew et al., 2012; McGaugh & Schombert, 2014; Cluver et al., 2014; Meidt et al., 2014; Querejeta et al., 2015; Hall et al., 2018; Schombert et al., 2019; Lelli et al., 2019). There is general agreement from models suggesting $\Upsilon^{*}_{W1}\simeq 0.5~{}M_{\odot}/L_{\odot}$. The choice of stellar initial mass function is a dominant source of systematic uncertainty. As McGaugh & Schombert (2015) have pointed out, constraints on $\Upsilon^{*}_{W1}$ based on the minimization of scatter of the BTFR provides constraints on systematics in the stellar mass distribution models.

With an initial assumption $\Upsilon^{*}_{W1}=0.5~{}M_{\odot}/L_{\odot}$, with observed H i and WISE $W1$ fluxes, and taking distances⁴⁴4Hubble flow distances are derived from velocities in the Local Sheet rest frame (Tully et al., 2008), assuming group mean velocities for galaxies identified to be in a group. The Local Sheet frame is preferred over the cosmic microwave background frame nearby because of a coherence in galaxy flows that extends over at least 40 Mpc. Velocities in alternate reference frames converge at large distances. from the Hubble law assuming $H_{0}=75$ km s^-1 Mpc^-1, we generate Figure 2. The value of $\Upsilon^{*}_{W1}$ is the single free parameter in this plot. Scatter arises from variance in $\Upsilon^{*}_{W1}$ and a multitude of other sources (uncertainties in H i fluxes, line widths, inclinations, photometry, and obscuration, distance deviations from the Hubble law, and intrinsic scatter).

At optical bands, where Population I components compete with Population II, $\Upsilon^{*}_{\lambda}$ have color terms (Taylor et al., 2011). The dependence with SDSS $g-i$ is shown in Figure 3 for the cases of the $r$, $i$ and $z$ bands. On the ordinate, we plot stellar mass to light ratios (assuming Hubble flow distances based on $H_{0}=75$ km s^-1 Mpc^-1). Here, where we require both optical and infrared information to be available, the stellar mass is based on photometry at the WISE $W1$ band, as discussed in the previous paragraph, while the absolute luminosity, $L_{\lambda}$, is derived from SDSS $\lambda=r$, $i$ or $z$ photometry.⁵⁵5All $SDSS$ optical and $WISE$ infrared photometry used in this study has been carried out from raw archival material and inclinations and obscuration issues have been given detailed attention as discussed exhaustively by Kourkchi et al. (2019); Kourkchi et al. (2020a, b). This mixed photometry provides a linkage between the infrared and optical, assuring that $<\Upsilon^{*}_{\lambda}L_{\lambda}>\simeq<\Upsilon^{*}_{W1}L_{W1}>$. The linkage is accomplished through the relations demonstrated by the blue lines in Figure 3

where $M_{*}$ is the stellar mass, and $\lambda$ is any of $r$, $i$ and $z$ bands. In this study, for the absolute magnitude of the Sun we adopt the values $4.64$, $4.53$, and $4.50$ mag in the AB system at SDSS $r$, $i$ and $z$ bands, respectively (Willmer, 2018). Our derived color dependencies are in a good agreement with the predictions of stellar population models (Bell et al., 2003; McGaugh & Schombert, 2014).

A tentative distance to a galaxy is given by its deviation from the fit to the data in Figure 2. Scatter about the fitted line arises for many reasons beyond unknowns that could be called intrinsic. Horizontal scatter arises from uncertainties in line width measurements and inclination adjustments. Vertical scatter arises from photometry errors (small), obscuration (small in the infrared), variations in $\Upsilon^{*}_{\lambda}$, and uncertainties in H i flux. Vertical scatter also arises from departures in distance from the expectations of uniform cosmic expansion. These distance components have a unique signature. Galaxies at distance $d$ that are greater than given by the Hubble law are intrinsically brighter than assumed, whence $M^{d}_{b}>M^{75}_{b}$ where these galaxies are given representation by $M^{75}_{b}$ in Figure 2. Such galaxies lie lower in this figure than they would if baryonic mass on the ordinate was based on their correct distance. The inverse is the case for galaxies nearer than distances given by the Hubble law; the baryonic masses of such galaxies will lie high with respect to what their placement would be if they were given correct distances.

If the only reason for scatter in the BTFR were variations in distance at a given systemic velocity then the distance for a galaxy, $j$, with respect to the arbitrary $H_{0}=75$ scale would be

where $d_{75}$ is the distance given by the Hubble law, and $M^{fid}_{b,j}$ is the baryon mass on the fiducial relation at the observed linewidth of the galaxy. Of course, only a small part of the scatter in the BTFR is attributable to distance deviations; indeed, is a component that diminishes roughly linearly with redshift. Assuming the scatter that is not associated with distance variations is randomly distributed about the fiducial relation, then measures of distances are embedded in deviations from the fiducial relation as systematics but buried in the substantial random uncertainties of non-distance causes. A more sophisticated recovery of distance estimates than offered by the simple application of Eq. 5 will be described in §4.

Although the BTFR seen in Figure 2 looks promising, it contains a bias that is revealed by considering the run of the Hubble parameter

given systemic velocity, $V_{ls,j}$, for galaxy $j$ with distance $d_{j}$. The function $f_{j}$ accounts for cosmic curvature⁶⁶6$f_{j}=1+1/2(1-q_{0})z_{j}-1/6(1-q_{0}-3q_{0}^{2}+j_{0})z_{j}^{2}$ where $z_{j}$ is the redshift of the galaxy, $j_{o}\simeq 1$, and $q_{0}=0.5(\Omega_{m}-2\Omega_{\Lambda})=-0.595$ assuming $\Omega_{m}=0.27$ and $\Omega_{\Lambda}=0.73$. (Visser, 2004) which is slight within the range of the current sample. The bias is in the sense that $H_{j}$ values drift lower (distances tending to overestimation) with increasing redshift. The bias arises because the primary selection of targets comes from H i surveys that are flux limited, particularly the dominant ALFALFA survey (Haynes et al., 2018). See Figure 4 and Dupuy et al. (2021) Figure 3. Photometry is comfortably available for any galaxy that meets the H i detection criterion. Nearby, H i is easily detected in any galaxy relevant to our interests. However, at increasingly great redshifts galaxies with lesser H i flux drop out of the sample. The consequence is manifested in the top panel of Figure 5 which plots the differential stellar and gas components vs. linewidth overlain with binned differences in increasing redshift intervals. It is clear that the measured H i component of mass becomes increasingly important at larger redshifts. The observed ratio $M_{\star}/M_{HI}$ decreases with increasing systemic velocity.

The bias against candidates with lesser H i flux is given quantitative expression in the middle panel of Figure 5. The difference, $\Delta$, gives the departure from the dashed straight line at a given ${\rm log}W^{i}_{mx}$ in the top panel. Letting $\Delta\equiv{\rm log}\widetilde{M}_{HI}-{\rm log}{M}_{HI}$, the linear relation in the middle panel is expressed as

What we will call the pseudo-gas mass parameter, $\widetilde{M}_{HI}$, is reduced from observed $M_{HI}$ values such that $<M_{\star}/\widetilde{M}_{HI}>\simeq constant$ with redshift. With this new parameter we can formulate a pseudo-baryonic mass, $\widetilde{M}_{b}$, that obeys a pseudo-BTFR

In §2.5, it will be demonstrated with mock data that unbiased distance estimates can be recovered based on Equation 8.

The separate mass contributions from stars and gas are shown as a function of H i profile linewidths in the top two panels of Figure 6. The bottom panel shows the ratio of the two for selections of the sample. The stellar mass has by far the tighter correlation with rotation. Gas masses peak abruptly at $M_{g}\simeq 4\times 10^{11}~{}M_{\odot}$ while stellar masses reach $M_{*}\simeq 2\times 10^{12}~{}M_{\odot}$. Gas masses are relatively unimportant for the largest galaxies.

In the bottom panel, a variant of the top panel of Figure 5, stellar to total gas mass ratios are shown for galaxies at the velocity extremes of our sample: those with velocities less than 4,000 ${\rm km~{}s^{-1}}$ are identified in blue while those with velocities greater than 10,000 ${\rm km~{}s^{-1}}$ are identified in red. The displacement toward lower $M_{*}/M_{g}$ values at higher velocities, the consequence of the H i flux bias, is evident. However another point can be made. The horizontal dotted line separates between an upper zone of stellar mass dominance and a lower zone of gas mass dominance. Our sample includes all manor of galaxies within the low velocity (distance) cutoff of 4,000 ${\rm km~{}s^{-1}}$. Galaxies with large systemic velocities (distances) such as those represented in red in the figure are massive galaxies with the stellar component dominant. This situation somewhat mitigates the H i flux bias problem. The H i detection bias is aggravated at large distances but most sample targets at large distances have only minor gas contributions. Nearby, the gas contributions can be substantial but well formulated because the bias is small.

The redshift dependency of the stellar to H i gas ratios is also evident in the top row panels of Figure 7. Several trends are evident. For one, there is an upward shift in stellar to H i mass with increasing linewidths. Second, the preponderance of samples shifts to higher velocities with increasing linewidth intervals. But third and a direct consequence of the H i flux bias, there are downward trends in stellar to H i mass ratios with systemic velocity across all the upper panels. The bottom row panels of Figure 7 are generated in a similar fashion using the H i pseudo-mass, $\widetilde{M}_{HI}$. The consequence of the adjustments is a flattening of the slope of the stellar to H i gas mass ratios with redshift. Utilizing $\widetilde{M}_{HI}$ minimizes the redshift dependent scatter about the pseudo-BTFR and results in distances free of the H i flux bias.

The adjustments made to compensate for the bias caused by the H i flux cutoff can be evaluated through the study of a mock catalog of objects with known distances and otherwise similar observed properties as our real sample. We construct a mock catalog with the following recipe. (1) Create a sample that populates the fiducial BTFR assuming distances based on $H_{0}=75$ km s^-1 Mpc^-1 and assuming a Schechter function that is appropriate for H i selection (Kourkchi et al., 2020b). The construction assumes the gas to stellar mass fraction given by Hall et al. (2012) with the scatter prescribed in that paper. (2) Create Gaussian scatter in both log linewidth and log baryon mass. (3) Translate this sample to a multitude of distance shells so now the operational parameters are apparent magnitudes and observed H i fluxes. Populate the shells such that after imposing an H i flux limit there is rough agreement with the histogram with redshift of the observed sample. (4) Only retain candidates that exceed an H i flux limit. (5) Determine the distances that would be measure for the candidates that are retained and compare those distances with the known input distances.

The analysis of the mock material follows what was done with the real data. The top panel of Figure 8 shows the stellar to gas mass fractions as a function of linewidths with superimposed symbols giving binned averages within systemic velocity intervals. In the bottom panel of this figure, individual points record the vertical position of a candidate with respect to the dashed line in the top panel as a function of systemic velocity. Mean values in velocity intervals are represented by open symbols and a linear fit is given to the distribution. The slope of this line of $-5.79\pm 0.05\times 10^{-5}$ is very close to that of Eq. 7 describing the real data.

The fit in the bottom panel of Figure 8 provides the information needed for adjustments to mock distance moduli. Figure 9 illustrates the results with respect to input values scattered around $H_{0}=75$ km s^-1 Mpc^-1 for the raw observed moduli in red and the adjusted moduli in green. The raw observed moduli drift to higher distances at larger velocities as expected from the bias. The adjustment returns distances that are close to the mean values that would be recovered if there were no H i flux cutoff. These same results are seen in Figure 10 where the ordinate is now the Hubble parameter $V_{j}/d_{j}$, with means in velocity bins and horizontal fits represented by the dashed lines for values at $V>4000$ ${\rm km~{}s^{-1}}$. In the lower panel adjusted case, the input value of $H_{0}$ is recovered and beyond 4,000 ${\rm km~{}s^{-1}}$ the averaged values are stable about constant $H_{0}$.

We carried out SDSS photometry at $u,g,r,i,z$ bands (Kourkchi et al., 2020a, b). The $u$ material will not be considered further because measurement uncertainties are larger in that band. The TFR at $g$ has increased scatter because of stochastic contributions from young populations. The utility of $g$ photometry is its service in identifying such contributions, as seen in Figure 3. On the other hand, the TFR at $r$, $i$ and $z$ bands have comparable scatter, although color terms are slightly different.

Our aim in this section is to explore the space of the model free parameters to find the best set of parameters that minimize the scatter about the linear relation between baryonic mass and H i linewidth in Equation 8 where $\widetilde{M}_{b}$ is the pseudo-baryonic mass and is derived using the adjusted H i mass, ${\rm log}\widetilde{M}_{HI}={\rm log}{M}_{HI}-C_{adj}V_{ls}$. Here, $C_{adj}$ is the adjustment factor and one of the model parameters. In this section, we relax the assumption we made in Section 2 about the non-variability of $\Upsilon^{*}_{W1}$ across all of our sample. We attempt to initially determine the functional form of the $\Upsilon^{*}_{\lambda}$ at optical wavebands and then find an equivalent version at infrared bands. Accordingly, we define $\Upsilon^{*}_{\lambda}$ as a linear function of the optical colors $g-i$ given by Equation 4.

In our optimization process, the objective is to construct the posterior distribution of the model parameters, $\mathcal{P}(\Theta|\mathcal{D})$, where $\Theta$ is the vector that holds all model parameters (i.e. $Slope$ and $ZP$ of the BTFR, $\alpha_{\lambda}$ and $\beta_{\lambda}$, and the bias factor $C_{adj}$. The observed data, $\mathcal{D}$, consist of linewidths and galaxy luminosities at different wavebands.

The posterior distribution can be determined following the conditional probability law, $\mathcal{P}(\Theta|\mathcal{D})\propto\mathcal{P}(\mathcal{D}|\Theta)\mathcal{P}(\Theta)$. We set $\mathcal{P}(\Theta)=1$ given a lack of any prior knowledge on the distribution of the model parameters. This assumption reduces our problem to finding a set of optimum parameters that maximize the likelihood function.

Having adopted the likelihood function we take advantage of the Markov Chain Monte Carlo (MCMC) method to explore the parameter space and sample the posterior distribution. This process allows us to incorporate all uncertainties of the observables upon which our model is constructed. Data points are mutually independent since they represent different galaxies. Moreover, it is a reasonable assumption that all uncertainties on the measured observables are almost Gaussian implying the following likelihood function

where $n$ is the galaxy index and $N$ is the total number of galaxies in our analysis. To minimize the effect of the Malmquist bias, the data$-$model residuals ($r_{n}=\mathcal{D}_{n}-\mathcal{M}_{n}(\Theta)$) are calculated along the linewidth parameter (Strauss & Willick, 1995; Tully & Pierce, 2000). Here, $\sigma_{n}$ is the uncertainty of the residual parameter, $r_{n}$, which penalizes a model based on the uncertainty of the real measurements and predictions, and is calculated using $\sigma_{n}^{2}=\sigma_{D_{n}}^{2}+\sigma_{\mathcal{M}_{n}}^{2}$. For each galaxy $\mathcal{D}_{n}={\rm log}W^{i}_{mx,n}$ and $\sigma_{D_{n}}$ is the uncertainty on the measured linewidths. For any model with given parameters, $\Theta$, the model uncertainty $\sigma_{\mathcal{M}}$ is calculated through a propagation of uncertainties of the model observables, i.e. H i flux, magnitudes, and color indices.

In this study, we sample the posterior distribution using Python package emcee (Foreman-Mackey et al., 2013), emcee leverages fast linear algebra routines to accelerate the process and takes the logarithm of the posterior likelihood

where ${\bf r}$ is the $N\times 1$ data$-$model residual as define earlier, and $\Sigma$ is the $N\times N$ diagonal covariance matrix where $\Sigma_{n,n}=\sigma^{2}_{n}$.

We adopt this likelihood for each of $r$, $i$ and $z$ bands and generate 64 chains each with the length of 10,000 samples. Each chain initially starts from a random location in the parameter space. This allows to expand the exploration area and obtain more robust results. Mostly after the first 500 burn in steps, each chain converges and follows the Markov statistics. However, to be more conservative, we remove the first 1,000 steps of each MCMC sample and combine all 64 walkers to generate the posterior distribution of model parameters.

The corner plot for the posterior distribution of model parameters based on the $i$ band photometry is illustrated in Figure 11. In this diagram, the top-most panel of each column displays the one dimensional distribution of the corresponding sampled parameter, with the overplotted red solid line representing the median values, and the lower/upper bounds corresponding to 16/84 percentiles (black dashed lines). Horizontal and vertical red lines are drawn at the location of the median values that we adopt as the optimum values of our model parameters. The shape of the posterior distribution looks almost the same if either $r$ or $z$ band photometry is used to build the BTFR.

Table 1 lists the optimized parameters of our model with the use of different optical bands and color combinations to calculate the stellar to mass ratios (Equation 4). The most probable $Slope$ and $ZP$ that we find through the MCMC simulations are in good agreement with each other at different wavebands. We tested that adding an extra linear dependency on the $r-z$ color term does not significantly alter the BTFR parameters.

The top row panels of Figure 12 illustrate the BTFR after adjustments for the H i flux selection bias, for $r$, $i$ and $z$ bands, where for each waveband the model parameters are given in the corresponding row for the band in Table 1.

Our tools for measuring galaxy distances can be slightly sharpened by using a combination of $r,i,z$ material. This unification can be accomplished in two ways. One is to extract the stellar mass of galaxies from the best fits of $\Upsilon^{*}_{\lambda}$ at different wavebands (Table 1). The other way is define an average linear relation

where $L_{riz}$ is the average of the absolute luminosities of a galaxy at $r$, $i$ and $z$ bands. We fit the parameters of this relation together with the BTFR coefficients through the same MCMC procedure presented in Section 2.7. Table 1 lists the output of our analysis and the bottom left panel of Figure 12 illustrates the fitted straight line of the corresponding BTFR model.

The same methodology can be used to find model parameters at the infrared WISE $W1$ band. However, optical color indices are not available for galaxies in our sample that only have infrared photometry coverage. Therefore, we need a formulation that only uses infrared information to establish empirical functions involving $\Upsilon^{*}_{W1}$. At the same time, though, we need to maintain a coherence between the optical and infrared measurements. Consequently, in deriving the functional form of the infrared parameter $\Upsilon^{*}_{W1}$, we restrict to galaxies with both optical and infrared data coverage. We derive their stellar masses, $M^{*}_{riz}$, following Equation 4 and the optimized parameters listed in Table 1 using combined $<riz>$ optical measurements. The optical$-$infrared coupling is achieved by requiring

which results in $M^{*}_{riz}/L_{W1}\simeq\Upsilon^{*}_{W1}$.

Figure 13 plots the optical stellar mass to $W1$-band luminosity versus the $W1-W2$ color. Compared to $g-i$ color with a dynamical range of $\sim 1.5$ mag, $W1-W2$ of most galaxies falls in a very short range $\sim 0.2$ mag. Still, a reasonable correlation pattern is identified when we bin data points. $M^{*}_{riz}/L_{W1}$ decreases toward redder $W1-W2$ colors at intermediate $W1-W$ values. Statistics at blue and red extremes are poor where dependencies are taken to be approximately constant. The purple dashed line in Figure 13 shows the best fit given by the maximum likelihood method. $M_{*}/L_{W1}$ takes the constant values of $0.66\pm 0.09$ on the blue side of $W1-W2=-0.63\pm 0.01$ mag. As the $W1-W2$ color gets redder, $M_{*}/L_{W1}$ linearly decreases with color until reaching the value of $0.51\pm 0.08$ at $W1-W2=-0.46\pm 0.04$ mag. Redward, the value is taken to be constant at $0.52\pm 0.08$.

The scatter about this formulation is shown in the lower panel of Figure 13. The rms scatter in $M^{*}/L_{W1}$ of $0.11$ about a median value of $0.59$ corresponds to $19\%$ uncertainty on the stellar mass to light ratios.

Figure 14 shows our model of $M_{*}/L_{W1}$ together with three theoretical models that are constructed by adopting various initial mass functions and stellar evolution scenarios (Hall et al., 2018). These models describe different galactic environments and are only valid for restricted ranges of $W1-W2$. As seen, our fitted function for the $M_{*}/L_{W1}$ color dependency (purple line) is in general consistent with those advocated by Cluver et al. (2014), and Eskew et al. (2012). Most of our data points are located on the right side of the applicability zone of the Meidt et al. (2014) curve. Our constraints on $M_{*}/L_{W1}$ at both $W1-W2$ color extremes are weak.

Adopting this multi-linear model for $\Upsilon^{*}_{W1}$, we follow the same MCMC procedure to map the posterior distributions of the BTFR parameters using $W1$ band photometry. Results are summarized in Figure 15. Here, only Slope and ZP are adopted as free parameters, because the parameters describing $\Upsilon^{*}_{W1}$ have been already found by the analysis of the galaxies that are in common between our optical and infrared sub-samples. Fixing the $\Upsilon^{*}_{W1}$ parameters guarantees that the estimated stellar mass to light ratios are in agreement with those estimated using the optical data and improves the consistency between the BTFR parameters we find in different scenarios.

The parameters of the fitted BTFR with the $W1$ band photometry are included in Table 1. The infrared correlation is illustrated in the bottom right panel of Figure 12. The WISE infrared and SDSS optical versions are consistent with each other thanks to the coupling imposed through Equation 12.

In the previous section, we described how a galaxy distance is calculated for a given set of observables and model parameters. Here, prior to calculating the uncertainty of the measured distances, we study the contribution of each error term in the scatter of data points about the BTFR and subsequent distance measurement uncertainties. The 1$\sigma$ scatter of our sample galaxies is $0.18$ dex about the direction of baryonic mass in the BTFR diagram. This scatter is roughly equivalent to an uncertainty of $0.45$ mag in the distance modulus measurements. We want to evaluate the contribution of different sources of uncertainty to the ensemble value.

In our study, uncertainties have mainly four sources: (1) uncertainties associated to the model parameters that describe correlations between different quantities, i.e. the errors on the fitted coefficients, (2) errors that are related to the scatter of the model residuals that basically monitors how well the constructed model performs, (3) uncertainties of the observables that are used in the model, and (4) the unknown uncertainties that are related to the underlying physical processes that govern the dynamic and composition of our sample galaxies.

To estimate the effect of each error term, first we generate an ensemble of 10,000 galaxies by randomly drawing from the original sample. In this process, a galaxy may appear multiple times in the ensemble, however its parameters would be dispersed differently. For each galaxy in the ensemble, we adopt the corresponding baryonic mass and distance from the original sample. Assuming that the simulated galaxies completely obey the BTFR gives us a way to compute their linewidths. By introducing uncertainties of different natures to the parameters of the ensemble galaxies, we measure how galaxies are dispersed along the vertical axis of BTFR which is then translated to the error in the distance modulus by $d\/DM=2.5\/d{\rm log}M_{b}$.

Below, we calculate the contribution of various uncertainties in the total scatter of the BTFR. For the sake of discussion, the BTFR and other relevant parameters are taken from our $riz$ model. The error budget is roughly the same if we consider other BTFRs at different wavebands. In this analysis the values and uncertainties of the H i flux ($F_{21}$) and linewidth (${\rm log}W^{i}_{mx}$) are taken from Table 1 of Kourkchi et al. (2020b) and given again in §7 Table 3. The uncertainties on the apparent magnitudes at SDSS optical ($r$, $i$, and $z$) and WISE infrared ($W_{1}$ and $W_{2}$) bands are less than $0.05$ mag, the value we conservatively adopt in this analysis.

• •

  Uncertainties in the calculated stellar mass to light ratios that have three sources. One originates from the errors in the coefficients $\alpha_{\lambda}$ and $\beta_{\lambda}$ in Equation 4. The other is related to the uncertainty of the $g-i$ color, which is $\sim 0.07$ mag considering that the uncertainties on the apparent magnitudes at SDSS optical ($r$, $i$, and $z$) and WISE infrared ($W_{1}$ and $W_{2}$) bands are less than $0.05$ mag. Given the rms scatter of the residuals of the $\Upsilon^{*}_{\lambda}$ models, we associate an uncertainty of $0.1$ dex on the logarithm of the calculated stellar mass to light ratios, $\rm{log}\Upsilon^{*}_{\lambda}$ $[M_{\odot}/L_{\odot}]$ (see Figures 3 and 13). This level of uncertainty is equivalent to an average of $\sim 20\%$ error on the derived $\Upsilon^{*}_{\lambda}$. Dispersing the stellar mass to light ratios of the ensemble test galaxies and recalculating their total baryonic mass, results in an scatter of $0.08$ dex on the ${\rm log}\widetilde{M}_{b}$ which implies an uncertainty of $0.2$ mag on the distance modulus.

• •

  Uncertainty of the gas mass. Assuming a Gaussian error of $0.1$ on $K_{g}$ and in calculation of the baryonic mass from Equation 2 generates a scatter of 0.01 dex on ${\rm log}\widetilde{M}_{b}$ which is equivalent to the moduli uncertainty of $0.02$ mag. Another error factor in the calculation of the gas mass is the error of the measured H i flux. Dispersing the data points based on the error of the H i flux generates a dispersion of less than $0.01$ dex along ${\rm log}\widetilde{M}_{b}$, equivalent to an error of $0.02$ mag in distance modulus. The combined error of 0.03 is minor.

• •

  The error of the velocity multiplier $C_{adj}$ in Equation 7 and the scatter of residuals about the linear correlation shown in the middle panel of Figure 5. This scatter has been modeled by the dashed line in the bottom panel of Figure 5. Considering these two types of errors introduces an scatter of $0.08$ dex along ${\rm log}\widetilde{M}_{b}$, an error of $0.2$ mag on distance modulus.

• •

  Uncertainty of the H i line width, ${\rm log}W^{i}_{mx}$. We randomly disperse the linewidth of the simulated galaxies according to the error on the ${\rm log}W^{i}_{mx}$ parameter which already includes the uncertainty of the linewidths and inclinations measurements. A dispersion of $0.1$ dex on ${\rm log}\widetilde{M}_{b}$ is implied which is translated to $\sim 0.25$ mag error on the distance modulus.

• •

  Adding all mentioned known uncertainties we end up with a scatter of $0.15$ dex in ${\rm log}\widetilde{M}_{b}$ which is smaller than the 0.18 dex that we observe overall. Following the quadratic rule of addition of Gaussian errors, this unexplained scatter is equivalent to $0.1$ dex in ${\rm log}\widetilde{M}_{b}$ that translates to $23\%$ uncertainty on the calculated baryonic mass from the BTFR correlation. We attribute this residual to the intrinsic uncertainty of the BTFR correlation which could originate from unknown physical processes such as variations between the distributions of baryonic and dark matter. We include this unknown uncertainty in our analysis to get more reasonable errors on the measured distances.

Now that we have estimated the sources and ensemble contributions of different uncertainties in the final measurements, we can focus on each galaxy individually. In practice, we measure the distance of each galaxy using an ensemble of 10,000 sets of parameters that are generated by randomly dispersing the input observables based on their reported Gaussian errors. In addition to accounting for the uncertainties of the galaxy observables, each measurement adheres to the BTFR and $\Upsilon^{*}_{\lambda}$ models whose parameters are drawn from either their posterior distributions or the covariance matrices that are generated in the maximum likelihood fitting processes. The scatter about the fitted correlations is considered in the same process as explained above. In the end, we add an additional uncertainty of $0.1$ dex on ${\rm log}\widetilde{M}_{b}$ when measuring distances of the ensemble galaxies. For each galaxy, we perform 10,000 measurements using the simulated parameters and we report the median of all deduced distances as the final galaxy distance and adopt the $1\sigma$ scatter as the statistical error.