Cosmic clocks: a tight radius–velocity relationship for H i-selected galaxies

Abstract

H i-selected galaxies obey a linear relationship between their maximum detected radius R_max and rotational velocity. This result covers measurements in the optical, ultraviolet, and H i emission in galaxies spanning a factor of 30 in size and velocity, from small dwarf irregulars to the largest spirals. Hence, galaxies behave as clocks, rotating once a Gyr at the very outskirts of their discs. Observations of a large optically selected sample are consistent, implying this relationship is generic to disc galaxies in the low redshift Universe. A linear radius–velocity relationship is expected from simple models of galaxy formation and evolution. The total mass within R_max has collapsed by a factor of 37 compared to the present mean density of the Universe. Adopting standard assumptions, we find a mean halo spin parameter λ in the range 0.020–0.035. The dispersion in λ, 0.16 dex, is smaller than expected from simulations. This may be due to the biases in our selection of disc galaxies rather than all haloes. The estimated mass densities of stars and atomic gas at R_max are similar (∼0.5 M_⊙ pc⁻²), indicating outer discs are highly evolved. The gas consumption and stellar population build time-scales are hundreds of Gyr, hence star formation is not driving the current evolution of outer discs. The estimated ratio between R_max and disc scalelength is consistent with long-standing predictions from monolithic collapse models. Hence, it remains unclear whether disc extent results from continual accretion, a rapid initial collapse, secular evolution, or a combination thereof.

galaxies: dwarf, galaxies: fundamental parameters, galaxies: kinematics and dynamics, galaxies: spiral, galaxies: structure

1 INTRODUCTION

Based on the cold dark matter (CDM) scenario for galaxy evolution, the main structural and dynamical properties of galaxies’ haloes and discs are expected to obey simple virial scaling relations (Fall & Efstathiou 1980; Mo, Mao & White 1998; Dutton et al. 2007). These properties are typically specified as a radius R, rotation velocity amplitude V, and a mass M, or alternatively luminosity L as a proxy for mass. For haloes in virial equilibrium, we expect V ∝ R ∝ M^1/3 (Mo et al. 1998, hereafter MMW98). Although the dark matter (DM) is not directly observable, scaling relations are observed in the properties of the baryons, although the slopes (power-law exponents) of the relations are not exactly as predicted for the haloes (e.g. Courteau et al. 2007).

The most used scaling relation is the velocity–luminosity relation, better known as the Tully–Fisher relation (hereafter TFR; Tully & Fisher 1977), and similarly the Baryonic TFR (McGaugh et al. 2000), which is a velocity–mass relationship. Baryonic physics is messy. The scaling between luminosity and baryonic mass depends on the star formation history, which varies between galaxies (Grebel 1997; Tolstoy, Hill & Tosi 2009; Weisz et al. 2011; Williams et al. 2011), the initial mass function (IMF), which also apparently varies between galaxies whether they are dominated by young stellar populations (Hoversten & Glazebrook 2008; Lee et al. 2009; Meurer et al. 2009; Gunawardhana et al. 2011) or old ones (Treu et al. 2010; Cappellari et al. 2012; Conroy & van Dokkum 2012; Dutton, Mendel & Simard 2012; Smith, Lucey & Carter 2012; van Dokkum & Conroy 2012), and the dust content and distribution (Calzetti, Kinney & Storchi-Bergmann 1994; Gordon et al. 2001; Tuffs et al. 2004). Theory and observations indicate that feedback from star formation (Governato et al. 2010; Oh et al. 2011) or active galactic nuclei (e.g. Bonoli et al. 2016) can rearrange the distribution of baryons, and in the process drag along the DM into an altered distribution, affecting all scaling relations.

The radius–velocity (RV) relationship has received somewhat less attention. Courteau et al. (2007) and Dutton et al. (2007) fit scaling relations to R, V, and L in a sample of luminous spiral galaxies having optical spectroscopic observations. They found that the scatter in the RV relationship was the highest compared to the LV and RL relationships. Some of the scatter in the RV relation is due to the uncertainties and ambiguities of measuring R. This includes lack of a uniform definition of radial scalelength (cf. Pohlen & Trujillo 2006), contamination by the bulge component, selection effects (especially with surface brightness), and errors due to dust. However, if instead of using a scalelength to characterize R we consider an outer radius, then some of these concerns (e.g. bulge and dust) are minimized and a tighter relationship can be found. This will allow a better measurement of the intrinsic scatter in the RV relationship, which is very sensitive to the spin of the haloes in which galaxies lie (e.g. Mo et al. 1998; Courteau et al. 2007; Dutton et al. 2007; Obreschkow & Glazebrook 2014).

Here, we demonstrate a nearly linear RV relationship in various measurements of H i-selected galaxy samples. In Section 2, we present our primary samples and detail the measurements we use. Section 3 shows the observed correlations and quantifies the slopes and scatters; we also test the results on a large comparison sample selected and measured in the optical giving consistent results. In Section 4, we show what a linear RV relation means in the context of CDM-dominated galaxy evolution models. Our results are discussed further in Section 5, where we estimate the spin parameter of galaxies, the properties of discs near their outer extents, and then discuss how these results relate to ideas on what limits the extent of galaxy discs. Our conclusions are presented in Section 6.

2 SAMPLES AND MEASUREMENTS

We measure the RV relationship in three primary samples. The first uses optical data from the Survey of Ionization in Neutral Gas Galaxies (SINGG; Meurer et al. 2006), which is an H α and R-band follow-up survey to the H i Parkes All Sky Survey (HiPASS; Koribalski et al. 2004; Meyer et al. 2004; Zwaan et al. 2004) combined with single dish H i data from HiPASS. The second uses data from the Survey of Ultraviolet emission in Neutral Gas Galaxies (SUNGG; Wong 2007), which observed HiPASS-selected galaxies in the ultraviolet (UV) with GALEX, for a sample largely overlapping with SINGG. Here, we use sub-samples of SINGG and SUNGG designed to ensure that reasonable rotation amplitudes can be derived from the HiPASS H i data. Specifically, both samples are selected to have major to minor axial ratios a/b ≥ 2 and to be the only apparent star-forming galaxy in the system. The a/b cut guarantees a minimum inclination of about 60°, thus limiting projection errors in calculating orbital velocities. For those galaxies observed by SINGG, the isolation criterion was determined using the H α images, which are roughly the same size as the HiPASS beam. For those SUNGG galaxies not observed by SINGG, isolation was determined morphologically; systems with companions of similar angular size, obvious signs of interaction, or noted as interacting with another galaxy in the literature were excluded. These selection criteria result in 71 and 87 galaxies from SINGG and SUNGG, respectively, with an overlap of 47 galaxies in common. The third sample uses H i imaging data of the 20 galaxies studied by Meurer, Zheng & de Blok (2013, hereafter MZD13).

In all three samples, H i data are used to infer the maximum rotation amplitude. The implicit assumption is that the H i in these galaxies is dominated by a rotating disc. It is important to bear in mind that the selection of the samples requires detectable amounts of H i, and thus is biased against gas-poor disc galaxies (e.g. S0 galaxies and ellipticals). Note that the a/b cut applied to the SINGG and SUNGG samples also is likely to remove early type and S0 galaxies from our samples. As pointed out by Meurer et al. (2006), very few such galaxies are found in the SINGG sample. The selection against early type galaxies may have implications on the types of haloes they are associated with, as discussed in Section 5.1. As we show below, the implied rotational amplitudes range from ∼10 km s⁻¹ (dwarf galaxies) to ∼300 km s⁻¹ (the largest spirals).

The radii used for the SINGG and SUNGG samples depend on the maximum extent of the galaxies observed in the optical and UV, respectively. Both surveys are designed to measure the total light of extended nearby galaxies using a series of concentric elliptical annuli. For SINGG, the apertures are set in a manner slightly modified from that given in Meurer et al. (2006). As noted there, the aperture shape (a/b and position angle) and centre are set by eye to include all the apparent optical emission. In most cases, this shape matches well the apparent shape of the galaxy in the R band, i.e. a tilted disc. We then grow the apertures to an arbitrarily large size, and determine, by eye, where the raw (before sky subtraction) radial surface brightness profile levels off. The surface brightness of the galaxy at that radius is on the order of 1 per cent of the sky brightness.

The radius where the raw surface brightness profiles flatten is called the maximum radius R_max. Since the R-band light almost always can be traced further than H α, R_max typically measures the maximum detectable extent in the optical continuum. Most exceptions are dwarf galaxies with strong minor-axis outflows. Optical sizes were estimated in this manner by two of us. First by DH and then by GRM, who ‘tweaked’ the size estimates in about half of the SINGG sample. Typically those that were adjusted were made larger because the raw profiles indicated that a small amount of additional flux could be gained doing so. Here, we use the tweaked aperture radii. Compared to using the initial estimates, the use of the tweaked apertures increases R_max by 0.06 dex on average and also reduces the scatter in the residuals of the fits described below by 0.06 dex (when taken in quadrature). The SUNGG maximum radius is set in a similar manner; it is determined separately in near-ultraviolet (NUV) and far-ultraviolet (FUV), and the maximum of the two is taken as R_max.

For both the SINGG and SUNGG samples, we interpolate enclosed flux versus aperture semi-major axis profiles to determine the radii containing 50 per cent (R₅₀) and 90 per cent (R₉₀) of the flux in the R band and UV, respectively.

For the MZD13 sample, we use three radii: R_max is the maximum extent of the H i radial profiles as given in the original studies used by MZD13, whereas R₁ and R₂ represent the extent of the region of the H i surface mass density profile |$\Sigma _{\rm H\, \small {I}}$| fitted with a power law by MZD13. These radii are set by eye to mark kinks in the H i radial profiles, indicating changes of slope in |$\log (\Sigma _{\rm H\, \small {I}}(R))$|⁠. On average, they are close to the radii that contain 25 per cent and 75 per cent of the H i flux, respectively (MZD13). Unfortunately, neither MZD13 nor the studies they employed calculated R₅₀ and R₉₀ for the H i data.

The shape of the rotation curve (hereafter RC) V(R) of galaxies varies systematically with mass, or peak rotational velocity, from nearly solid body (linearly rising), for the lowest mass galaxies, to RCs that are flat at nearly all radii, or even slightly declining at large R for the most massive galaxies (Persic & Salucci 1991; Persic, Salucci & Stel 1996; Catinella, Giovanelli & Haynes 2006). Unless stated otherwise, we take V to be the maximum rotational amplitude. For most cases, this will be the amplitude at the flat part of the RC. In the majority of other cases, it will be the farthest measured point of the RC. We take these definitions to be synonymous. For the SINGG and SUNGG samples, we derive V from the full width at half-maximum of the H i spectrum from HiPASS, assuming a flat RC over all relevant radii. We follow the method of Meyer et al. (2008) and correct the line widths for inclination, and broadening resulting from turbulence, relativity, instrumental effects, and data smoothing. As with Meyer et al. (2008) the inclinations are derived from a/b. For the H i sample we interpolate the RCs, from the various original studies used by MZD13 to arrive at rotation amplitudes at R₁, R₂, and R_max separately (i.e. V(R₁), V(R₂), and V(R_max)).

3 RESULTS

3.1 Observed correlations

We show the observed correlations separately for each data set in three figures. Fig. 1 shows the RV relationship for the SINGG optical data. In the left-hand panel, the y-axis gives the radius as R₅₀(R), i.e. the radius containing 50 per cent of the R-band light; similarly the middle panel shows R₉₀(R) as the radius; whereas the right-hand panel uses R_max as defined from the SINGG optical data. The velocity in all panels is the circular velocity V defined from the HiPASS H i line widths (Section 2). Similarly, Fig. 2 shows R₅₀(NUV), R₉₀(NUV), and R_max from the SUNGG UV data in the left-hand, middle, and right-hand panels, respectively, against V derived from HiPASS. Fig. 3 shows the H i radii R₁, R₂, and R_max plotted against the circular velocities interpolated at those radii V(R₁), V(R₂), and V(R_max) in the left-hand, middle, and right-hand panels, respectively.

Figure 1.

Radius R plotted against circular velocity V on a logarithmic scale as derived from SINGG optical data. The radii plotted in the left-hand, middle, and right-hand panels are R₅₀(R), R₉₀(R), and R_max(optical), respectively, whereas all three panels plot the same HiPASS H i based V. The solid line shows the iteratively clipped ordinary least squares bisector fit to the plotted quantities, whereas the dotted line shows the fit offset by ±3σ_log(R), where σ_log(R) is the dispersion of the R residuals. Each panel is annotated with the mean log orbital time, 〈log (t_orb)〉, and its dispersion (after clipping), the fitted slope β and its error, and the rms of the residuals in the ordinate. The parallel grey dashed lines from the bottom to top show where t_orb = 10⁸,  10^8.5,  10^9.0,  10^9.5, and 10¹⁰ years, respectively.

Figure 2.

Radius R plotted against circular velocity V on a logarithmic scale as derived from SUNGG UV data. The radii used here are R₅₀(NUV), R₉₀(NUV), and R_max(UV) in the left-hand, middle and right-hand panels, respectively. The meanings of the various lines and annotations are the same as in Fig. 1.

Figure 3.

Radius R plotted against circular velocity V on a logarithmic scale for the H i sample of MZD13. The radii R₁ (left-hand panel) and R₂ (centre panel) delimit the region where the H i surface brightness profile is a power law with index γ ≈ −1 (see MZD13 for details), whereas R_max (right-hand panel) is the maximum detected extent of H i. The meanings of the various lines and annotations are the same as in Fig. 1.

We fit the RV relations in log–log space as

\begin{eqnarray} \log (R) = \alpha + \beta \log (V) \end{eqnarray}

(1)

using an ordinary linear least squares bisector algorithm (Isobe et al. 1990) weighting each point equally, and iteratively clipping points that deviate from the fit by more than three times the dispersion in R. Table 1 reports the results of the fits, giving the coefficients α, β, the dispersion of the residuals σ_log(R), σ_log(V), and Pearson's correlation coefficient r_xy. Some of these quantities are also listed in Figs 1–3.

Table 1.

Fit parameters.

Sample	Radius	Nused	Nrej	α	β	log (torb)	σlog(R)	σlog(V)	rxy
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
SINGG	R50(R)	71	0	1.42 ± 0.12	1.07 ± 0.06	8.34 ± 0.15	0.15	0.14	0.899
SINGG	R90(R)	71	0	1.71 ± 0.12	1.09 ± 0.06	8.68 ± 0.14	0.15	0.14	0.908
SINGG	Rmax(opt)	71	0	1.95 ± 0.12	1.13 ± 0.06	9.00 ± 0.14	0.15	0.13	0.914
SUNGG	R50(NUV)	88	0	1.11 ± 0.11	1.33 ± 0.08	8.54 ± 0.22	0.18	0.24	0.797
SUNGG	R90(NUV)	88	0	1.56 ± 0.15	1.24 ± 0.08	8.81 ± 0.21	0.18	0.22	0.805
SUNGG	Rmax(UV)	87	1	2.16 ± 0.12	1.04 ± 0.06	9.02 ± 0.16	0.16	0.16	0.830
H i	|$R_{\rm 1}({\rm H\, \small {I}})$|	20	0	1.49 ± 0.15	1.09 ± 0.08	8.45 ± 0.23	0.24	0.22	0.859
H i	|$R_{\rm 2}({\rm H\, \small {I}})$|	19	1	1.64 ± 0.19	1.18 ± 0.10	8.81 ± 0.15	0.16	0.14	0.876
H i	|$R_{\rm max}({\rm H\, \small {I}})$|	20	0	1.99 ± 0.24	1.10 ± 0.11	8.98 ± 0.17	0.18	0.16	0.825
PS1	R50(r)	692	6	1.32 ± 0.07	1.07 ± 0.03	8.26 ± 0.13	0.14	0.13	0.770
PS1	Rb(r)	694	4	1.41 ± 0.07	1.14 ± 0.03	8.52 ± 0.14	0.15	0.13	0.761
PS1	R90(r)	689	9	1.68 ± 0.06	1.07 ± 0.03	8.62 ± 0.12	0.12	0.11	0.818

Column (1): the galaxy sample fitted. Column (2): the radius measured. Column (3): the number of data points used in the fit. Column (4): the number of data points rejected from the fit. Column (5): the zeropoint of the fit. Column (6): the slope of the fit. Column (7): average log of the orbital time of the fitted data points. Column (8): dispersion in the log of the residuals in radius R of the fitted points. Column (9): dispersion in the log of the residuals in orbital velocity V of the fitted points (or implied orbital velocity V^΄ for the PS1 sample). Column (10): Pearson's correlation coefficient using all data points.

For the optical and UV samples, the fits are the ‘best’ at R_max, where best is defined as having the highest r_xy and lowest σ_log(R) and σ_log(V). For the H i sample, the fit at R₁, marking where the H i profiles flatten, is much worse than the other two fits. The flattening is likely to be due to the increasing dominance of molecular gas at small radii (Bigiel et al. 2008; Leroy et al. 2008). The fits at R₂ and R_max have similar scatters, indistinguishable statistically. In summary, the fits are their best, or close to it, at R_max, where β is close to but slightly greater than unity, that is, a linear relationship.

A linear RV implies that the orbital time

\begin{eqnarray} t_{\rm orb}= \frac{2\pi R}{V} \end{eqnarray}

(2)

is constant (assuming the orbit shape is well approximated by a circle). We list the mean log (t_orb) in Table 1 and in the panels of Figs 1–3. The RV relation at R_max is nearly identical in the three figures even though R_max is defined at very different wavelengths, which are sensitive to different physical processes. Fig. 4 a overplots the three samples at R_max, showing the excellent correspondence in the RV relationships. They all imply that t_orb ≈ 1 Gyr, with a scatter of 0.14–0.18 dex (38 per cent to 51 per cent). Thus, H i-selected disc galaxies behave like clocks and rotate once in a Gyr at their outermost detected radii. This holds for galaxies which range in radius from R_max ∼ 1.5 kpc, having V ∼ 10 km s^{− 1} to those with R_max ∼ 50 kpc and V ∼ 300 km s^{− 1}. The RV relationship for t_orb = 1 Gyr is shown with the dashed line in Fig. 4a.

Figure 4.

Maximum radius R_max plotted against circular velocity V on a logarithmic scale for our combined primary samples is shown in the left-hand panel. Here, orange circles, purple diamonds, and green triangles show the SINGG, SUNGG, and H i samples displayed in Figs 1c, 2c, and 3c, respectively. In the right-hand panel, we show other estimates of the maximum disc size compared to rotational amplitude. In the right-hand panel, the blue four pointed stars show the optically determined R_max of edge-on galaxies for the sample of Kregel, van der Kruit & de Grijs (2002) against their V. The red pentagons show the H ii region truncation radius R_HII plotted against V for galaxies in common between the samples of Martin & Kennicutt (2001) and MZD13. The grey dashed line in both panels shows the relationship expected for t_orb = 1 Gyr.

The SINGG–RV relation is equally well defined at R₅₀, R₉₀, and R_max. However, the meaning is less clear when using R₅₀ and R₉₀. The velocity used, V, is determined from the line width of integrated H i velocity profiles of galaxies that are spatially unresolved. The H i in galaxies typically is weighted to larger radius than the easily observed optical emission (e.g. Leroy et al. 2008), hence the derived V is also applicable to large radii. The rotation amplitude at R₉₀ and R₅₀ will be systematically overestimated using V as one goes to lower rotation amplitudes and shorter radii (i.e. the effect will be stronger for R₅₀ than R₉₀). Hence, if we used the true V values at R₉₀ and R₅₀ then we should see shallower β values than shown in Fig. 1.

The RV relations in the UV also are defined using H i velocity profiles. Here, we see significantly larger σ_log (R) residuals when using R₅₀ and R₉₀ as well as steeper β values compared to the RV relation at R_max. We posit that the worse fits are due to whether or not galaxies have a central starburst, and the degree to which they are affected by dust. These will have more of an impact on the distribution of the UV luminosity at small radii than in the determination of R_max.

3.2 Results for a comparison sample

Our primary results are for three samples having H i-based V measurements, one has H i-based R measurements, and two have overlapping selections based on H i properties. In order to address whether our results may be a byproduct of working with H i data, we now consider a sample that is selected and measured in the optical; the sample of 698 disc galaxies of Zheng et al. (2015). The sample was selected from the Pan-STARRS1 (PS1) Medium Deep Survey (Chambers et al. 2016) fields, and measured from the stacked survey images. The galaxies were selected to have images in all PS1 bands (g, r, i, z, y), spectroscopic redshifts from SDSS-III, ¹ to be fairly face-on (a/b < 2), and to be well resolved with a Petrosian (1976) radius²R_P > 5 arcsec. Galaxy profiles are then measured to 2R_P. This algorithm recovers ≥ 94 per cent of the total light for galaxies having a Sérsic (1963) index n ≤ 2 typical of disc galaxies (Graham et al. 2005). Zheng et al. (2015) found that radial surface brightness profiles typically show a ‘break’, or change in slope, in the bluer bands with the break less apparent towards longer wavelengths. They fitted stellar population models to annular photometry in the five bands to derive stellar mass density profiles and integrated to yield the total stellar mass. They recorded R₅₀, R₉₀, and the break radius R_b all measured in the r band. Hence, as with the other samples, we have three fiducial radii to work with. Instead of using a measured rotational velocity, we use |$V^{\prime }_{\rm c}$|⁠, the circular velocity estimated from the stellar mass based TFR of Reyes et al. (2011). This fit to the TFR has been shown to well represent the kinematics of an SDSS based sample (Simons et al. 2015) that has redshifts similar to this PS1 sample.

The resulting RV^΄ relationships are shown in Fig. 5. We fit this sample in the same manner as done for the other samples (Section 3.1). The fit parameters are tabulated in Table 1. In all three cases, the fits to the PS1 data are nearly linear (β ≈ 1), with the fit at R₉₀ being closest to linear and having the smallest scatter |$\sigma _{\log (R_{\rm 90})} = 0.12$| dex of any of our RV fits.

Figure 5.

Radius R plotted against estimated circular velocity V^΄ on a logarithmic scale for the Pan-STARRS1 sample of Zheng et al. (2015). Here, V^΄ is not a direct measurement of circular velocity, but estimated from the stellar mass based TFR of Reyes et al. (2011). The radii used here are the radius containing 50 per cent of the PS1 r-band light R₅₀(r)), the break radius R_b(PS1), and the radius containing 90 per cent of the r-band light R₉₀(r) in the left-hand, middle and right-hand panels, respectively. The meanings of the various lines and annotations are the same as in Fig. 1.

Amongst the PS1 sample fits, the one at R_b has the largest scatter, 0.15 dex in |$\sigma _{{\rm log}(R_{\rm b})}$|⁠, and deviates the farthest from linear (β = 1.12). Nevertheless, the scatter about the mean orbital time of 0.13 dex compares well with the other RV fits. The larger scatter compared to the fit at R₉₀ may arise because the strength of the radial profiles breaks is highly variable with some galaxies ‘breaking down’ (more typical), others ‘breaking up’, and some showing no discernible breaks (Freeman 1970; Pohlen & Trujillo 2006; Zheng et al. 2015).

The (logarithmic) mean t_orb at R₅₀ and R₉₀ for the PS1 sample 0.18, 0.42 Gyr, respectively, is close to that for the SINGG sample 0.22, 0.48 Gyr. Meanwhile, at R_b, the mean t_orb = 0.36 Gyr, intermediate between those at R₅₀ and R₉₀. Hence, R_b ≈ 1.8R₅₀ ≈ 0.8R₉₀ for the PS1 sample. Although Zheng et al. (2015) do not measure R_max, they note that typically R₉₀ = R_P and that most of the light is recovered at 2R_P, hence we expect R_max ≈ 2R₉₀ for the PS1 sample. For flat RCs, then we infer that t_orb(R_max) = 2t_orb(R₉₀) = 0.83 Gyr for the PS1 sample, within 0.08 dex of the SINGG sample. The t_orb estimates for the PS1 sample at R₅₀, R₉₀ and that implied at R_max are all lower than those for the SINGG sample, suggesting a more general optical selection of galaxies may result in smaller galaxies than an H i selection. However, the differences are all close to or about equal to the 0.06 dex systematic error noted in Section 2. Hence, to that level of accuracy, the same RV relationship for H i selected galaxies applies to all disc galaxies at low redshifts.

The somewhat tighter fits to the PS1 sample does not necessarily mean that the intrinsic scatter in the RV relations is lower than for the SINGG sample. This is because an inferred rather than measured velocity is used. Since the V^΄ in Fig. 5 is derived from luminosities, these are essentially RL or RM_⋆ correlations we are showing. Saintonge & Spekkens (2011) find a very tight RL relation (having a scatter of 0.05 dex in R) using their SFI++ sample of spiral galaxies, and employing isophotal radii and I-band luminosities. Similarly, both Courteau et al. (2007) and Hall et al. (2012) find smaller scatters in their RL relations than their VR relations. In part, this is because errors in L are effectively reduced by a factor of 3–4 due to the TFR scaling, making them competitive or better than velocity errors (Saintonge & Spekkens 2011). Velocities are also more prone to errors in inclination, position angle, and non-circular motions. Working with stellar mass (fitted to photometry), as done with our PS1 sample, also improves the accuracy by effectively spreading any error over five bands and weighting the results to the reddest bands. But improved accuracy may not be the only cause for the tight fits in Fig. 5. Saintonge & Spekkens (2011) performed a careful error analysis of the scatter in their scaling relations and estimate the intrinsic scatter in their RL relations (∼0.034 dex in R) is less than half of that in their RV relations (∼0.084 dex in R).

4 EXPECTATIONS FROM SIMPLE GALAXY EVOLUTION MODELS

A constant t_orb at R_max implies a constant spherically averaged mean mass (baryons and DM) density ρ interior to R_max since

\begin{eqnarray} \rho = \frac{3\pi }{Gt_{\rm orb}^2}, \end{eqnarray}

(3)

where G is the gravitational constant. Our adopted t_orb(R_max) = 1 Gyr yields the mean mass density interior to R_max:

\begin{eqnarray} \langle \rho (R_{\rm max})\rangle = 2.1\times 10^{-3}\, {\cal M}_{\odot }\, {\rm pc^{-3}}. \end{eqnarray}

(4)

The closure density of the Universe, ρ_c, is given by

\begin{eqnarray} \rho _c = \frac{3H_z^2}{8\pi G}, \end{eqnarray}

(5)

where H_z is the redshift (z) dependent Hubble constant. This can be used to estimate the collapse factor of matter within R_max. Adopting H₀ = 70 km s⁻¹ Mpc⁻¹ and the results of the Planck Collaboration (Planck Collaboration XVI 2014) that the ratio of the cosmic matter density to the closure density Ω_M = 0.315, we have

\begin{eqnarray} \frac{\langle \rho (R_{\rm max})\rangle }{\Omega _{\rm M}\rho _c} = 49000. \end{eqnarray}

(6)

The third root of this is the average collapse factor f_c(R_max) of the matter within R_max compared to the present day matter density of the Universe:

\begin{eqnarray} \langle f_{\rm c}(R_{\rm max})\rangle = 36.6. \end{eqnarray}

(7)

The ‘virial radius’ R₂₀₀ is usually defined as the radius where the mean density of the enclosed mass is 200 times larger than ρ_c. From equation (5), it is apparent that ρ_c depends only on redshift, and thus, by definition, a linear RV relationship is expected at R₂₀₀ at any given epoch given by equation (2) of MMW98:

\begin{eqnarray} R_{200} = \frac{V}{10 H_z}. \end{eqnarray}

(8)

From equation (2), the orbital time at the virial radius at the present epoch is t_orb = 8.8 Gyr.

The RC interior to R₂₀₀ depends on the distribution of DM and baryons. MMW98 used an analytical approach to examine the expected structure of disc galaxies within DM haloes under a variety of plausible assumptions about cosmogeny and distribution of the baryons and DM. They adopt a simple isothermal sphere to parametrize DM haloes, and show that this framework is convenient for understanding how the galaxy scaling relations are influenced by the properties of their haloes. Adopting this approach, MMW98 derived equation (8). An isothermal sphere has a flat (constant) RC and a density profile

\begin{eqnarray} \rho (R) = \frac{V^2}{4\pi G R^2}. \end{eqnarray}

(9)

A flat RC is well supported observationally in most disc galaxies, especially at large radii (e.g. Rubin, Thonnard & Ford 1978; Bosma 1981; Mathewson, Ford & Buchhorn 1992; de Blok et al. 2008; Epinat et al. 2008), whereas the shape of the inner part of the RC varies systematically with mass, or V (e.g. Persic & Salucci 1991; Catinella et al. 2006). Since our results are the most consistent at R_max, where RCs are typically flat, the isothermal approximation suffices for our purposes. For a pure exponential disc galaxy in a dominant isothermal halo, MMW98 derive the disc scalelength R_d relative to R₂₀₀ in their equation (12). Combining that with equation (8) yields

\begin{eqnarray} R_{\rm d} = \frac{V}{\sqrt{200}H_z}\left(\frac{j_{\rm d}}{m_{\rm d}}\right)\lambda , \end{eqnarray}

(10)

where j_d is the fraction of the total angular momentum in the disc, m_d is the fraction of the total mass in the disc, and λ is the spin parameter. For systems that are not purely exponential discs in an isothermal halo, the scalefactor (⁠|$1/[\sqrt{200}H_z]$|⁠) will vary depending on the detailed distributions of DM and baryons (MMW98). Thus, a linear relationship between R_d and V should exist if (j_d/m_d)λ is constant.

If the DM and baryons are well mixed when galaxies collapse, one would naively expect j_d = m_d and thus a constant j_d/m_d (MMW98). This is also the working assumption of Fall & Efstathiou (1980), whose simple models were consistent with the observations of the time. Although it is impossible to observationally confirm this expectation because of the invisible nature of DM, simulations allow it to be tested. Posti et al. (2018), using a similar approach to ours, and matching of galaxy properties to halo properties from a variety of recent cosmological N-body simulations find that |$j_{\rm d}/m_{\rm d}^{2/3}$| is approximately constant, close to what all we require. We note that a constant j_d/m_d is difficult to reproduce in more detailed numerical simulations (Governato et al. 2010). The dependence of λ on mass and other parameters for dark haloes, as measured in numerous simulations, is also weak (e.g. MMW98; Cole & Lacey 1996; Bett et al. 2007; Macciò et al. 2007). Thus, naive considerations tell us that we expect a linear RV relation when a disc scalelength, or anything proportional to it, is used to measure size. This would be the case for the isophotal sizes of pure exponential discs that have constant central surface brightness as originally proposed by Freeman (1970).

Disc galaxies, however, are not that simple. They typically contain a bulge increasingly apparent with morphological type (e.g. Hubble 1926). Since Freeman's landmark work, it has become apparent that discs obey a surface brightness–luminosity relationship (e.g. Kauffmann et al. 2003), and that the radial profiles frequently show breaks from being pure exponential (Freeman 1970; Pohlen & Trujillo 2006; Zheng et al. 2015). However, allowing for these complications may not necessarily cause major changes to the RV relation. MMW98 derive the behaviour of an exponential disc in a halo having the typical profile found in CDM-only simulations (Navarro, Frenk & White 1997, hereafter NFW), which is allowed to respond to the disc's mass. They find relationships for R_d and the maximum rotational velocity that differ from the isothermal case by form functions that depend on j_d, m_d, λ, and halo concentration c. Of these, c is the parameter that is expected to have the largest impact (Dutton et al. 2007). For example, MMW98 considered the case of a bulge plus disc embedded in an NFW halo, and found that disc size depends on assumptions about angular momentum transfer between the bulge, disc, and halo. They found that disc sizes can vary by a factor of about 2, whereas maximum velocities only vary by ≲ 20 per cent.

5 DISCUSSION

The formalism presented in Section 4 allows us to place our results in a cosmological context. We continue with this approach in Section 5.1 by examining the constraints on the spin parameter λ and its dispersion implied by our results. Section 5.2 discusses what our results imply for the properties at the disc outskirts. Section 5.3 argues that our results are best explained by a true physical truncation of discs. Although the formalism presented thus far implies continual accretion limits the extent of discs, Section 5.4 considers other scenarios for limiting the extent of discs. Finally, we present some ancillary implications of our results in Section 5.5.

5.1 Spin parameter

Equation (10) is readily re-arranged to be

\begin{eqnarray} \lambda = \frac{\sqrt{50}}{\pi }\frac{t_{\rm orb}(R)}{t_{\rm H}}\frac{R_{\rm d}}{R}, \end{eqnarray}

(11)

where t_orb(R) is the orbital time at radius R, and |$t_{\rm H} = H_z^{-1}$| is the Hubble time (13.96 Gyr for our H₀), and assuming j_d/m_d = 1. Thus, spin parameter can be estimated from the orbital time at a given radius and the scaling of that radius with disc scalelength. Since R_d was not measured in our samples, indirect estimates of this scaling must be made. We do this using two approaches.

First, if all baryons are in an un-truncated exponential disc, we can use the SINGG results shown in Fig. 1 and Table 1 to estimate λ. Noting that the radius containing 50 per cent and 90 per cent of the light of such a disc corresponds to 1.68 and 3.89 times R_d, and converting the mean orbital times in the log from Table 1, then we have λ = 0.021 and 0.020 estimated from t_orb(R₅₀) and t_orb(R₉₀), respectively. Being virtually identical, we take λ = 0.020 to be the average spin under the pure exponential disc assumption.

Secondly, we estimate R_max/R_d, and thus λ by scaling from the sample of Kregel et al. (2002) shown in Fig. 4b. They fit models including both an exponential disc and bulge to the light distribution of edge-on galaxies. Their disc model is truncated, yielding a maximum radius R_max, which they find to be on average a factor 3.6 ± 0.6 times larger than R_d. Their sample yields a significantly shorter average t_orb = 0.76 Gyr than what we find, probably due to systematic differences in how R_max is determined. If so, then we scale their results to estimate

\begin{eqnarray} R_{\rm max}\sim (3.6\pm 0.6) \frac{1}{0.76}R_{\rm d} = (4.7\pm 0.8)R_{\rm d}. \end{eqnarray}

(12)

Following equation (11), we have λ = 0.034 for t_orb = 1 Gyr. Since this scales from an estimate that avoids bulges, it produces a longer R_d scalelength and hence higher λ value than assuming all the light comes from an exponential disc.

In comparison, measurements of typical average spin parameters of haloes created in cosmological simulations range from λ = 0.03 to 0.055, with each simulation set producing a broad distribution that is close to log normal and consistently found to have a width in the range σ_log (λ) ≈ 0.21–0.23 dex (Cole & Lacey 1996; Bullock et al. 2001; Bett et al. 2007; Macciò et al. 2007). Our first estimate, λ = 0.020 (assuming pure exponential discs), is somewhat below the expectations of cosmological simulations, whereas the second estimate λ = 0.034 (from scaling the results of Kregel et al. 2002) is at the low end of the expectations from simulations. The 0.23-dex difference in these estimates is indicative of the systematic uncertainty involved. In addition, neglection of the gaseous disc, or equivalently assuming the same distribution for it as the stars, will underestimate the angular momentum of the baryons, and thus λ. Improved estimates of λ can be made with better modelling of the baryonic mass and angular momentum distribution in galaxies (e.g. Obreschkow & Glazebrook 2014; Butler, Obreschkow & Oh 2017), but would still require assumptions about the coupling with the unseen DM halo. Our approach using equation (11) assumes a singular isothermal sphere and j_d/m_d = 1, both of which might introduce additional systematic bias in our estimate of λ.

Despite the likely systematic offset between our observational estimate of λ (via equation 11) and its true value, we can none the less discuss the relation between the relative scatter in t_orb and λ, or, equivalently, the absolute scatter in logt_orb and log λ. The observed scatter in t_orb(R_max) of ∼0.16 dex has several sources: (1) the physical dispersion in λ, (2) measurement uncertainties in R_max and V_max, (3) variations in the ratio R_max/R_d, (4) variations in j_d/m_d, and (5) deviations from the assumed iso-thermal density profile. ³ We assume that (1) is the dominant source, but expect that the other sources make a non-negligible addition to the scatter of t_orb. In Section 2, we (crudely) estimated (2) the scatter in t_orb due to measurement errors as 0.06 dex. Removing this (in quadrature), but neglecting (3), (4), and (5), the scatter in λ is σ_log λ ∼ 0.15 dex. This is somewhat lower than predicted by CDM models (0.23 dex; Macciò, Dutton & van den Bosch 2008). This is remarkable, given that we have not even accounted for the scatter of t_orb(R_max) coming from the sources (3), (4), and (5).

The explanation for the relatively low scatter in t_orb(R_max) is likely two-fold. First, to the extent that the disc surface mass density at R_max is similar in all galaxies, ⁴ high spin systems are likely to have their disc truncated at smaller radii relative to R_d than low-spin systems. Therefore, the scatter in R_max/R_d (source 3) is negatively correlated to that of λ (source 1), hence reducing the scatter in t_orb(R_max), relative to the scatter in t_orb(R_d), which would be similar to the scatter in λ. R₅₀ is the closest proxy we have to R_d and we do indeed find that the scatter in t_orb(R_max) is less than that of t_orb(R₅₀) in both the optical and UV samples (Table 1). The effect is more prominent in the UV sample. Secondly, our sample is likely biased towards a more narrow range of spin parameters than present in a volume-complete sample of all DM haloes. The lowest λ haloes have little angular momentum and are more likely to be bulge dominated (i.e. S0 and elliptical galaxies); hence, they will have little H i, and not make it in to our samples. A rationale for our H i-selection producing a bias against high λ systems is less obvious. Effectively, all H i-selected galaxies are detected in the optical and UV (Meurer et al. 2006; Wong 2007); the detection limits are not biasing the samples. More speculatively, there may be a bias against high λ systems if the baryons they contain have not been able to cool enough for H i or stars to form.

5.2 Properties at disc galaxy outskirts

In Section 2, we defined R_max as the radius of readily detectable emission. It is largely determined by the amplitude of large-scale ‘sky variations’ in the R band (optical) and NUV (ultraviolet). These variations represent how well we can flat-field our data. The surface brightness of these variations provides a crude estimate of the limiting surface brightness at, or just interior to, R_max. The situation is slightly different for R_max(H i) – the limiting surface brightness is the measured |$\Sigma _{\rm H\, \small {I}}$| at the last measured point in published H i profiles. Of course, a galaxy may extend beyond R_max at fainter levels than the limiting surface brightness. Histograms of limiting surface brightness are shown in Fig. 6.

Figure 6.

Histograms of limiting surface brightness in the R band and NUV for the SINGG and SUNGG samples are shown in the top and middle panels, respectively. The bottom panel shows the Σ_gas at the farthest point in the radial profiles of the galaxies in the H i sample, hence they correspond to the faintest H i recorded for each galaxy. The dotted vertical line in each panel shows the median of the distribution. The top axis on the top and middle panels gives the limiting surface brightness in observed units. The bottom axes are calibrated in physically meaningful quantities: stellar mass density (Σ_⋆), star formation intensity (Σ_SFR), and H i surface mass density (Σ_gas) for the top, middle, and bottom panels, respectively.

The bottom axes of Fig. 6 show the limiting surface brightnesses converted to physically meaningful quantities. The R-band surface brightness μ_R is converted to the stellar mass density Σ_⋆, assuming a mass-to-light ratio M/L_R = 2 M_⊙/L_{R, ⊙}. For standard IMF assumptions, the adopted M/L_R is reasonable for a star-forming population, but probably will result in an underestimate for stellar mass densities if the relevant stellar population is old (Bell et al. 2003). To convert μ_NUV to star formation intensity, we adopt the FUV conversion factor of Leroy et al. (2008) and assume an intrinsic colour (FUV − NUV)₀ = 0 ABmag, which is reasonable for the outer discs of galaxies (Gil de Paz et al. 2007; Thilker et al. 2007; Zaritsky & Christlein 2007; Boissier et al. 2008; Goddard, Kennicutt & Ryan-Weber 2010; Hunter, Elmegreen & Ludka 2010; Werk et al. 2010a; Lee et al. 2011). The ISM density Σ_g assumes that the ISM is dominated by H i and is corrected by a factor of 1.3 to account for heavier elements.

The medians of the distributions are marked in Fig. 6 and correspond to Σ_⋆ = 0.42 M_⊙ pc^{− 2}, Σ_SFR = 2.8 × 10^{− 12}M_⊙ yr^{− 1} pc^{− 2}, and Σ_g = 0.58 M_⊙ pc^{− 2}. These may be considered typical conditions at or near R_max. The star formation intensity at R_max is weak compared to the stellar and gaseous contents. The time needed to form the observed stellar populations at the current star formation intensity is t_build = Σ_⋆/Σ_SFR = 150 Gyr, whereas the time required to process the gas through star formation is t_g = Σ_g/Σ_SFR ∼ 200 Gyr. Equivalently, both the specific star formation rate (⁠|$t_{\rm build}^{-1}$|⁠) and the star formation efficiency (⁠|$t_{\rm g}^{-1}$|⁠) are low in outer discs. Thus, at R_max, the current in situ star formation is too feeble to either create the stellar populations or transform the accumulated ISM into stars in a reasonable amount of time.

For a galaxy to have the same R_max in the R band, the NUV implies that the NUV-R colour at R_max is similar to the -‘colour-’ of the large-scale sky fluctuations, i.e. NUV-R ≈ μ_sky(NUV) − μ_sky(R) ≈ 3.8 ABmag. This is a ‘green-valley’ colour, i.e. intermediate between the blue and red sequences (Schiminovich et al. 2007), validating the long t_build we derive above.

The slope of the RV relationship in the R band (β_optical = 1.13 ± 0.06) is slightly steeper than in the NUV (β_UV = 1.04 ± 0.06). A comparison between Figs 1 and 2 shows that the values of R_max in the optical and UV are nearly equal to that at the high end, where V ≳ 200 km s⁻¹, whereas for V ≲ 50km s⁻¹, we find that, on average, galaxies have R_max(R) ≲ 0.7R_max(NUV). Hence, at R_max(NUV), galaxies are redder for large spirals than dwarfs. This may be due to the relative importance of an old component in the disc or halo for spirals compared to dwarfs. It may also be a sign of ‘downsizing’ that lower mass galaxies are less evolved into stars than high-mass galaxies.

5.3 The edge of the disc

Our results imply a distinct physical edge in the light distribution corresponding to R_max. The first line of evidence for this is the result that the three tracers give nearly identical estimates of t_orb; for a given V_max they yield the same radius in the distribution of stars, star formation, and atomic hydrogen. If discs were purely exponential, then the equality in t_orb would be remarkably coincidental, since the different measurements of R_max are set by independent observational limits for each tracer. If the observational limits were consistent within each band, then one could argue that R_max is effectively an isophotal radius. Previous studies (Saintonge et al. 2008; Saintonge & Spekkens 2011; Hall et al. 2012) have shown that isophotal radii produce tighter RV relationships than those using an exponential scalelength, perhaps because of the difficulty in consistently measuring R_d in the face of contamination from the bulge, breaks in radial profiles, and biases in setting the range of radii to fit with an exponential (Freeman 1970; Pohlen & Trujillo 2006; Hall et al. 2012; Zheng et al. 2015). However, as shown by Fig. 6, the limiting surface brightness is not consistent between galaxies, hence R_max is not an isophotal radius.

Indeed, the observed scatter in the RV relationship provides a second line of evidence that we are dealing with a truncation in the disc. The dispersion in limiting surface brightnesses shown in Fig. 6 is 0.40, 1.05, and 0.44 dex in the R band, NUV, and H i, respectively. Assuming a pure exponential disc and adopting equation (12) for the mean scaling between the disc scalelength and R_max, these dispersions should contribute 0.08, 0.27, and 0.09 dex to the respective scatter in the RV relationships, whereas the corresponding observed scatters are 0.15, 0.16, and 0.18 dex. Thus, the expected induced scatter, in this scenario, is larger than the observed scatter in the NUV, whereas it would make a considerable fraction ( ∼ 25 per cent in quadrature) of the observed scatter in the R band and H i.

That we are seeing a real edge to the disc is most apparent in the H i sample. Using the data from MZD13, we find an average power-law index γ = −4.6 ± 0.5 for |$\Sigma _{\rm H\, \small {I}}(R)$| profiles between R₂ and (H i) (the uncertainty is the standard error on the mean). If this slope is maintained towards larger R, the total H i content is well constrained. This is unlike the region R₁ to R₂, where the H i traces DM well but γ ≈ −1, which cannot be maintained indefinitely. Modern H i observations are sufficiently deep that large improvements in sensitivity of observations do not result in large changes to the H i content. For example, Gentile et al. (2013) present HALOGAS survey data of NGC 3198 with an H i surface brightness sensitivity 10 times fainter than the THINGS observations used by MZD13. Those improved observations result in an increase of 6 per cent in the H i flux, and 21 per cent (0.08 dex) in maximum radius compared to the THINGS data.

We conclude that discs are not purely exponential all the way to R_max, but must have a steep fall-off in surface brightness near R_max. An edge, or steep fall-off, in surface brightness has been noted in the optical by van der Kruit and collaborators (van der Kruit 2007; Kregel et al. 2002) and in H i by van Gorkom (1993). Our results are similar to theirs (Fig. 4), indicating that disc has nearly identical truncations at R_max in stars, star formation, and atomic hydrogen, and that it is this physical disc truncation that we are observing.

Baryons clearly exist beyond R_max in galaxies. For example, at the rotational amplitude of our Galaxy V = 220 km s⁻¹, then t_orb = 1 Gyr corresponds to R_max = 33 kpc. The RC of the Milky Way disc can be traced out to R ≈ 20 kpc (Sofue, Honma & Omodaka 2009; Burch & Cowsik 2013; Bhattacharjee, Chaudhury & Kundu 2014), whereas halo blue horizontal branch stars can be detected out to R ≈ 60 kpc (Xue et al. 2008) and globular clusters out to R ≈ 100 kpc (e.g. Pal 3; Koch, Côté & McWilliam 2009). In M31, the stellar disc can be traced to at least R = 40 kpc as shown by Ibata et al. (2005). Using their adopted RC (Klypin, Zhao & Somerville 2002), t_orb = 1.04 Gyr at this radius, nicely consistent with our average t_orb at R_max. Ibata et al. (2005) point out that additional fainter disc material may be detected out to 70 kpc, whereas Ibata et al. (2014) show faint but prominent features at larger radii relate to the halo, which extends to at least 150 kpc, about half the virial radius of R_vir ≈ 290 kpc (Klypin et al. 2002). Clearly there are stars well beyond where t_orb is 1 Gyr in both the Milky Way and M31. But they are primarily located in their host's halo, rather than disc.

5.4 Alternative mechanisms to truncate discs

The cosmological approach we adopted in Section 5.1 implies that discs grow with cosmic time (the |$H_z^{-1}$| dependence) due to accretion. Disc growth is also predicted in simple semi-analytic model extensions to cosmological N-body simulations, albeit with weaker growth (Dutton et al. 2011). However, other mechanisms may also be at play in setting the extent of galactic discs. These include the limitations in the angular momentum in an initial protogalactic collapse (van der Kruit 1987), truncation in star formation due to disc stabilization (Kennicutt 1989; Martin & Kennicutt 2001), ionization by the ultraviolet background (UVB; van Gorkom 1993), and spreading of the disc due to internal angular momentum transfer (Roškar et al. 2008a,b).

The fact that we see the linear RV expected for the cosmological accretion scenario is a strong argument in its favour. Likewise, simple semi-analytic models of galaxy evolution that incorporate accretion can account for the redshift evolution of the RV relationship and other virial scaling relations (Dutton et al. 2011). However, ‘smoking-gun’ observations of intergalactic gas being ‘caught in the act’ of accreting on to galaxies have been elusive. In a naive interpretation of the accretion scenario, one would expect outer discs to be largely gaseous. Instead, the equality of t_orb in the R band and H i combined with the near equality of |$\Sigma _{\rm H\, \small {I}}$| and Σ_* in the outer discs implied by Fig. 6 suggests that they are well evolved in to stars (albeit typically less so than inner disc). This conclusion should be considered tentative since our methods for estimating the various R_max values as well as |$\Sigma _{\rm H\, \small {I}}$| and Σ_* are crude.

An older scenario for producing a truncated but evolved outer disc is the concept of a rapid initial collapse of galaxies including their discs (Eggen, Lynden-Bell & Sandage 1962; Freeman 1970). van der Kruit (1987) shows that an initial uniformly rotating spherical gas cloud in a potential with a flat RC that collapses while conserving angular momentum will produce an exponential disc that truncates at 4.5 times the disc scalelength. In practise, Kregel et al. (2002) found that stellar discs truncate at R_max ∼ 3.6R_d, i.e. somewhat smaller. However, as argued in Section 5.2, they are likely measuring shorter R_max values than we do. Indeed, van der Kruit (2007) notes that the truncations examined by Kregel et al. (2002) correspond to μ_V ∼ 26.5–27.5 mag arcsec⁻², brighter than our estimates of the surface brightness at R_max (Fig. 6). The fact that van der Kruit (2007) often find H i beyond their optical truncation radii is consistent with them underestimating R_max compared to us, since our H i and optical R_max values are consistent. When we scale their results to our t_orb(equation 12), we find that the ratio between R_d and R_max is a factor of 4.7 ± 0.8, consistent with what is expected from a monolithic early collapse.

Kennicutt (1989) note that star formation, as traced by H ii regions in spiral galaxies, typically cuts off at a radius R_HII, beyond which few bright H ii regions are detected. Martin & Kennicutt (2001) confirmed this result with improved observations of more galaxies. Fig. 4 b plots R_HII for five galaxies from Martin & Kennicutt (2001), which are also in the sample of MZD13. The R_HII values (which correspond to t_orb = 48 to 390 Myr) are considerably smaller than the R_max values in our primary samples, but similar to the break radius R_b of the PS1 sample. Christlein, Zaritsky & Bland-Hawthorn (2010) find that the H α distribution of edge-on spirals typically has a downward break at 0.7R₉₀(R), which, they note, may correspond to the R_HII break. This scaling is very close to the R_b ≈ 0.8R₉₀ scaling we find for the PS1 sample, strengthening the notion that R_b and R_HII are related. The fact that Christlein et al. (2010) find H α emission beyond their break radius and the SUNGG UV measurements continue out to ∼2.2R₉₀(R) demonstrates that the limits of galaxies traced by prominent H ii regions does not measure the full extent of star formation in galactic discs.

Instead, UV emission is a better tracer of star formation in outer discs. The existence of extended UV (XUV) discs (Gil de Paz et al. 2005; Thilker et al. 2005, 2007) demonstrates that star formation can extend beyond the portion of the disc readily observed in the optical. These outer discs can also be probed using resolved stellar populations from the ground (Cuillandre et al. 2001; Ibata et al. 2005) or space (e.g. Bruzzese et al. 2015). The close match in the RV relationships at R_max shown in Fig. 4 implies that star formation extends to the limits of the H i disc.

One mechanism that has been promoted for limiting the extent of galaxy discs is ionization by the UVB posited by van Gorkom (1993) to explain the steep decline in |$\Sigma _{\rm H\, \small {I}}$| profiles at large R in NGC 3198 and other galaxies. The scenario was consistent with modelling of the time (Maloney 1993). The column densities he considered are similar to or somewhat smaller than the typical |$\Sigma _{\rm H\, \small {I}} \sim 0.1\, {\rm {\cal M}_{\odot }\, pc^{-2}}$| we find at R_max(H i). If ionized by the background is setting R_max(H i) then one should be able to detect emission from the ionized disc beyond R_max. Bland-Hawthorn, Freeman & Quinn (1997) present evidence for finding this emission in the outer disc of NGC 253. However, other searches for ionized disc gas beyond the H i edges of galaxies have not been successful (e.g. Madsen et al. 2001; Dicaire et al. 2008; Adams et al. 2011; Hlavacek-Larrondo et al. 2011). Recent very deep integral field spectroscopy of the outermost disc of UGC 7321 finds very low surface brightness H α, consistent with ionization by the UVB, but this emission does not extend beyond the |$\Sigma _{\rm H\, \small {I}} \sim 0.1\, {\rm {\cal M}_{\odot }\, pc^{-2}}$| contour (Fumagalli et al. 2017). Although UVB may ionize the ‘skin’ of H i discs, ionization gas does not extend much beyond the observable H i disc, which marks the true maximum extent of the cool ISM disc.

Roškar et al. (2008a,b) model the interplay between star formation and disc dynamics in isolated spiral galaxies. Their simulations produce star formation edges like that seen by Kennicutt (1989) and Martin & Kennicutt (2001), beyond which the gaseous part of the disc has a high Toomre (1964) disc stability parameter Q and thus produces little in situ star formation. Instead, most of the old stars at large radii formed at smaller radii and ‘migrated’ outwards due to resonances with transient spiral features. Such a process can explain the downward breaking surface brightness profiles, ‘U’ shaped age, and colour profiles commonly seen in spiral galaxies (e.g. Pohlen & Trujillo 2006; Bakos & Trujillo 2013; Zheng et al. 2015). The material in discs between R_HII and R_max may then be a combination of weak XUV disc star formation in the Q stable portion of the disc combined with outwardly migrating older stars. Although this scenario is appealing, it is not obvious how it would result in a linear RV relation largely consistent at different wavelengths down to the dwarf galaxy regime. Low-mass galaxies are also a concern because they do not have spiral density waves that are likely to drive radial migrations.

5.5 Other implications

There is a strong relationship between the H i radius and H i mass in galaxies of the form

\begin{eqnarray} M_{\rm H\, \small {I}} \propto R_{\rm H\, \small {I}}^2. \end{eqnarray}

(13)

This was emphasized recently by Wang et al. (2016), who note that it has been found for samples selected in a wide variety of ways (Broeils & Rhee 1997; Verheijen & Sancisi 2001; Swaters et al. 2002; Noordermeer et al. 2005; Wang et al. 2013). The correlation implies that the average H i surface brightness within |$R_{\rm H\, \small {I}}$| is constant. The scatter in this relationship is ∼0.06 dex, tighter than our RV relationship. The RV relationship at R_max is peripherally related to this result. It has long been known that a maximum |$\Sigma _{\rm H\, \small {I}} \sim 10\, {\rm {\cal M}_{\odot }\, pc^{-2}}$| is set by the conversion of the interstellar medium into a molecular form (e.g. Bigiel et al. 2008), and many galaxies reach this saturation in their central regions. The outer radius adopted by Wang et al. (2016) is where |$\Sigma _{\rm H\, \small {I}} = 1\, {\rm {\cal M}_{\odot }\, pc^{-2}}$| brighter than adopted for our H i sample (MZD13). This effectively limits the range of allowed average surface brightness. Within galaxies, H i has a predictable distribution giving a power-law fall-off in |$\Sigma _{\rm H\, \small {I}}$|⁠, which is apparently set by the disc maintaining a constant stability parameter (Meurer et al. 2013; Wong et al. 2016). The limited dynamic range of |$\Sigma _{\rm H\, \small {I}}$|⁠, combined with the shallow power-law radial profile, results in the narrow range of average |$\Sigma _{\rm H\, \small {I}}$|⁠.

Since radial density profiles typically decrease monotonically with R, the central density should not be less than the average density at R_max. This corresponds to a constraint on the slope of the inner RC of galaxies – the gradient should not be less than that implied by the RV relation at R_max, hence the orbital time should be less than or equal to ∼1 Gyr in the central parts of galaxies. Fig. 7 tests this assertion by plotting the histogram of central velocity gradients, d_RV(0), for 57 galaxies comprising the final sample of Lelli et al. (2013) with valid measurements. The dashed line shows the gradient d_RV(0) = 6.1 km s⁻¹ kpc⁻¹ corresponding to t_orb = 1 Gyr. There are no galaxies with a shallower gradient. The shallowest d_RV(0) = 9.0 km s⁻¹ kpc⁻¹ in their sample corresponds to the irregular galaxy IC 2574 (de Blok et al. 2008). Following the discussion in Sections 4 and 5, a galaxy with a central density less than 〈ρ(R_max)〉 would have had to have collapsed less than our samples, and that would mean they either have a higher λ or j_d/m_d > 1 (i.e. they have a larger fraction of the spin in the disc than the fraction of mass in the disc) or some combination of the two. Apparently, such galaxies have not (yet) formed.

Figure 7.

Histogram of central velocities gradients of a wide range of disc galaxies from the sample of Lelli, Fraternali & Verheijen (2013). The solid vertical line marks the velocity gradient equivalent to the mean orbital time t_orb = 1 Gyr at R_max of our samples, whereas the dashed lines indicate the dispersion about this mean time of 0.16 dex. Galaxies do not have central densities less than |$\rho = 2.1 \times 10^{-3}\, {\cal M}_{\odot }\, {\rm pc^{-3}}$|⁠.

6 CONCLUSIONS

We have shown that disc galaxies display a nearly linear radius–velocity (RV) relationship at the outermost radius R_max observed in the optical, UV, and H i emission at 21cm. The RV relationship is consistent between data sets and implies a constant orbital time of ∼1 Gyr at this radius. A comparison of our H i-selected and optically measured SINGG sample with the much larger optically selected and measured PS1 sample of Zheng et al. (2015) shows nearly identical RV relations at two fiducial radii. This suggests that our results are robust against the vagaries of sample selection and may be generic to disc galaxies in the low-redshift Universe. Within R_max, matter has collapsed by a factor of 37 to |$\rho _M = 2.1\times 10^{-3}\, {\cal M}_{\odot }\, {\rm pc^{-2}}$|⁠, a factor of 4.9 × 10⁴ times higher than the present day average matter density in the Universe.

We argue that R_max in our data sets corresponds to the edge of the disc. Recent studies indicate that R_max(H i) is limited by the available ISM in the disc rather than external ionization by the UVB. The star formation intensity at R_max is an order of magnitude too weak to build up the existing stellar populations or consume the available gas within a Hubble time. Hence, star formation at its current rate is not solely responsible for setting this radius. Although R_max appears to mark a sharp truncation in the disc of galaxies, it does not enclose all baryons. Stars in the halo are distributed to much larger radii, and their kinematics indicate the DM also extends further, likely to the virial radius.

Instead, R_max must be set by other processes such as accretion (e.g. Sancisi et al. 2008; Brook et al. 2012; Mollá et al. 2016). Continuous cosmic accretion provides a natural explanation for the RV relation. In that scenario, the RV relation gives a constraint on the average spin parameter, which we estimate to be in the range λ = 0.02–0.035. This estimate is likely to be biased due to the crudeness of our estimates and the requirement for H i in our samples, which will bias them against the typically gas poor and low spin elliptical and S0 galaxies. The scatter in the orbital times provides a constraint on the dispersion of spin parameters σ_log (λ) ≈ 0.16 dex, somewhat smaller than expected by theory (∼0.22 dex), probably also due in part to the previously mentioned biases. The scatter in orbital times may also underestimate that in λ if R_max corresponds to a consistent disc surface brightness or mass density.

An older theory, consistent with our results, is that R_max is set by a rapid early collapse. Unfortunately, this scenario does not make a prediction on the RV relationship. However, a crude estimate of the scaling of R_max with disc scalelength R_d (following the results of Kregel et al. 2002) is consistent with long-standing theoretical predictions for this scenario (R_max/R_d ∼ 4.5 van der Kruit 1987). Our estimates of the gas and stellar surface densities near R_max are very similar, indicating a high degree of evolution of outer discs. The relatively flat metallicity gradients in the outskirts of galaxies also indicates a high degree of chemical evolution in disc outskirts (Werk et al. 2010b, 2011). Hence, an early rapid collapse model is nominally consistent with our results. However, our estimates of R_max/R_d and the surface densities at R_max are crude. Better estimates are needed to test this interpretation.

The RV relationship has some practical implications. Since the conversion of angular to physical radius is distance dependent, whereas the conversion of velocities is not (to first order), then one could use our RV relationship to estimate distances. However, since the observed relationship is linear it is not as powerful as the TFR where luminosity goes as orbital velocity to a power of three to four (e.g. Meyer et al. 2008). Furthermore, due to the likely evolution in this relationship (Dutton et al. 2011), one must take care to limit its use to the local Universe.

A simple RV scaling relation provides a convenient tool to estimate the extent of galaxy discs. We used the R_max found here in the model we developed to explain the nearly constant ratio of star formation rate (as traced in the UV) to the H i mass (Wong et al. 2016). Further development of this model would be useful for determining a wide range of properties along the star-forming main-sequence of galaxies. Of particular relevance would be using such an approach, combined with observed column density distributions within galaxies, to model the likely cross-section of H i absorbers (e.g. Rao & Briggs 1993; Ryan-Weber, Webster & Staveley-Smith 2003, 2005; Zwaan et al. 2005; Braun 2012). Similarly a realistic disc truncation radius could be usefully employed in setting the initial conditions for detailed dynamical simulations of local galaxies, or for modelling the inclusion of baryons in semi-analytic models.

Acknowledgements

GRM acknowledges useful discussions with Tim Heckman, Brent Tully, Ken Freeman, Joss Bland-Hawthorn, Emma Ryan-Weber, and Martin Zwaan. GRM also thanks the Center for Astrophysical Sciences of the Johns Hopkins University, and the National Astronomical Observatories of the Chinese Academy of Sciences for their hospitality during visits while this paper was being developed. ZZ is supported by the National Natural Science Foundation of China Grant No. 11703036. Partial funding for the SINGG and SUNGG surveys came from National Aeronautics and Space Administration (NASA) grants NAG5-13083 (Long-Term Space Astrophysics program), GALEX GI04-0105-0009 (NASA GALEX Guest Investigator grant), and NNX09AF85G (GALEX archival grant) to GRM. The SINGG observations were made possible by a generous allocation of time from the Survey Program of the National Optical Astronomy Observatory (NOAO), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under a cooperative agreement with the National Science Foundation. GALEX is a NASA Small Explorer, launched in 2003 April. We gratefully acknowledge NASA's support for construction, operation, and science analysis for the GALEX mission, developed in cooperation with the Centre National d'Etudes Spatiales of France and the Korean Ministry of Science and Technology. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Footnotes

1

http://skyserver.sdss.org/

2

where the local surface brightness is 20 per cent of the average interior surface brightness.

3

This includes variations in the ratio of circular velocity measured over the disc to that at the virial radius, V/V₂₀₀.

4

In Section 5.2, we show that the surface brightness limits at R_max varies greatly, but this does not preclude the corresponding mass densities just interior to where this limit is found to be similar.