A massive stellar bulge in a regularly rotating galaxy 1.2 billion years after the Big Bang
Federico Lelli, Enrico M. Di Teodoro, Filippo Fraternali, Allison W. S. Man, Zhi-Yu Zhang, Carlos De Breuck, Timothy A. Davis, Roberto Maiolino
AA massive stellar bulge in a regularly rotating galaxy1.2 billion years after the Big Bang
Federico Lelli, , , ∗ , Enrico M. Di Teodoro , Filippo Fraternali , Allison W.S. Man ,Zhi-Yu Zhang , Carlos De Breuck , Timothy A. Davis , Roberto Maiolino , School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, UK Arcetri Astrophysical Observatory, INAF, Florence 50125, Italy; ∗ e-mail: [email protected] Department of Physics & Astronomy, Johns Hopkins University, Baltimore MD 21218, USA Kapteyn Astronomical Institute, University of Groningen, Groningen 9700 AV, The Netherlands Dunlap Institute for Astronomy & Astrophysics, University of Toronto, Toronto ON M5S 3H4, Canada School of Astronomy and Space Science, Nanjing University, Nanjing 210023, P.R. China European Southern Observatory, Germany Headquarters, Garching bei M¨unchen 85748, Germany Kavli Institute for Cosmology, University of Cambridge, Cambridge CB3 0HA, UK Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK1 a r X i v : . [ a s t r o - ph . GA ] F e b osmological models predict that galaxies forming in the early Universe expe-rience a chaotic phase of gas accretion and star formation, followed by gas ejec-tion due to feedback processes. Galaxy bulges may assemble later via mergersor internal evolution. Here we present submillimeter observations (with spa-tial resolution of 700 parsecs) of ALESS 073.1, a starburst galaxy at z (cid:39) whenthe Universe was 1.2 billion years old. This galaxy’s cold gas forms a regularlyrotating disk with negligible noncircular motions. The galaxy rotation curverequires the presence of a central bulge in addition to a star-forming disk.We conclude that massive bulges and regularly rotating disks can form morerapidly in the early Universe than predicted by models of galaxy formation. In the standard Λ cold dark matter ( Λ CDM) cosmological model, galaxies form inside darkmatter (DM) halos when primordial gas cools, collapses, and begins star formation ( ). At earlytimes, when the Universe was only 10 % of its current age (redshift z (cid:39) ), Λ CDM modelspredict that galaxies undergo a turbulent phase of gas accretion from the cosmic web, leadingto violent bursts of star formation and black hole growth ( ), rapidly followed by massive gasoutflows due to feedback from supernovae and active galactic nuclei (AGN). Strong stellar andAGN feedback are needed in Λ CDM models to reproduce the observed number of galaxiesper unit volume per unit mass (the galaxy stellar mass function) and to quench star formationin the most massive DM halos, reproducing the observed population of quiescent early-typegalaxies at z = 0 ( ). In the hierarchical Λ CDM scenario, galaxy bulges (the spheroidal stellarcomponents in the central parts of massive galaxies) are expected to form either after the mergerof two galaxies of similar mass, or via secular dynamical processes within the stellar disk ( ).Observation of galaxies at high redshifts enables us to test these scenarios by searching for thesephysical processes in action.The 158- µ m emission line of singly ionized carbon, [C II ], is a major coolant of the interstel-2ar medium of galaxies and traces a combination of atomic and molecular hydrogen phases. The[C II ] line has been observed up to z (cid:39) . ( ) and several rotating disks have been identifiedat z > ( ). Most existing observations only marginally resolve the [C II ] distribution andkinematics, and thus cannot distinguish between rotating disks, galaxy mergers, or gas inflowsand outflows. Limited spatial resolution has prevented the measurement of well-sampled rota-tion curves (with more than six independent elements), which can be used to measure differentmass components by fitting dynamical models.We used the Atacama Large Millimeter/submillimeter Array (ALMA) to observe a star-forming galaxy at z = 4 . ( ) at a spatial resolution of ∼ (cid:48)(cid:48) , corresponding to 0.7 kpcin a Λ CDM cosmology ( ). The target, ALESS 073.1 (LESS J033229.3-275619), is a strongsub-millimeter source ( ). Its large far-infrared luminosity suggests an ongoing starburst withan estimated star-formation rate of ∼ (cid:12) ) per year (
10, 11 ). ALESS 073.1contains a dust-obscured AGN ( ). These extreme properties indicate that the galaxy maydrive a massive gas outflow ( ). Previous [C II ] observations with ALMA indicated a rotatingdisk ( ), but the spatial resolution (0.5 (cid:48)(cid:48) corresponding to ∼ µ m, which traces dust heated by thestar-formation activity. Figure 1B shows the integrated [C II ] intensity map, which traces thecold gas distribution. The [C II ] emission is about twice as extended as the dust emission andshows an asymmetry to the northeast. Similar lopsided gas disks are common in the nearbyUniverse, forming nearly half of the local population of atomic gas disks ( ).The [C II ] velocity field (moment-one map; Fig. 1C) shows a regularly rotating disk. Ro-tating disks and galaxy mergers can be difficult to distinguish when observed at low spatialresolutions ( ). This is not the case for ALESS 073.1 because the area with detectable [C II ]emission is covered by ∼
45 resolution elements. The kinematic major axis is perpendicular to3igure 1:
ALMA observations of ALESS 073.1. ( A ) Continuum emission at 160 µ m (rest-frame) tracing dust heated by young stars, ( B ) [C II ] intensity map tracing cold gas, and ( C )[C II ] velocity field showing a rotating disk. North is up and east is left. The kinematic center,located at a Right Ascension (R.A.) of 03 h m s and Declination (Dec.) of − ◦ (cid:48) (cid:48)(cid:48) , is represented by a white star. The beam size is plotted as the grey ellipse in the bottom-left corner. The physical scale is indicated by the scale bar in the bottom-right corner. In (A),iso-emission contours range from . to 1 mJy beam − (where 1 mJy = 10 − W m − Hz − )in steps of 0.11 mJy beam − . In (B), iso-emission contours range from 0.35 to 6 mJy beam − km s − in steps of 0.7 mJy beam − km s − . In (A) and (B), the lowest iso-emission contourcorresponds to a signal-to-noise ratio of ∼
3. In (C), iso-velocity contours range from − to +120 km s − in steps of 30 km s − , the bold contour indicates the systemic velocity (set tozero), and the grey line shows the kinematic major axis.4he kinematic minor axis, implying that noncircular motions are negligible in the inner partsof the gas disk ( ). The gas kinematics are less regular in the northeastern extension: Thismay be due to recently accreted gas, a warped outer disk, or streaming motions along a spiralarm ( ). The limited signal-to-noise ratio in the outer regions does not allow us to distinguishbetween these possibilities. Overall, the kinematic regularity of the [C II ] disk is unexpected fora starburst AGN-host galaxy at z (cid:39) .We model the [C II ] kinematics using the software B AROLO ( ), which fits the observa-tional data with a rotating disk model ( ). This determines the [C II ] surface brightness profile,the rotation curve, and the intrinsic velocity dispersion profile. The best-fitting rotation curve(Fig. 2A) rises steeply in the central parts and declines across the disk. Similar rotation curvesare observed in bulge-dominated disk galaxies in the nearby Universe ( ), so there mightbe a bulge component in ALESS 073.1. The normalization of the rotation curve depends on thedisk inclination, which we treat as a free parameter below.We use the derived rotation curve to constrain the mass distribution within the galaxy. Ourfiducial mass model has four components: cold gas disk, stellar disk, stellar bulge, and DMhalo. We also investigate mass models with three components, which we find either cannot fitthe observations or are less plausible in the Λ CDM paradigm ( ). We assume that the [C II ]density profile traces the distribution of the cold gas disk, whereas the dust density profiletraces the distribution of the stellar disk (Fig. 2B), because dust absorbs ultraviolet radiationfrom young stars and re-emits it at far-infrared wavelengths. The two additional components(bulge and DM halo) are modeled with analytic functions ( ). The bulge component implicitlyincludes the contribution of the central super-massive black hole, which contributes < % ofthe central mass in high- z galaxies ( ).For the best-fitting inclination of 22 ◦ , the [C II ] disk rotates at velocities ( V rot ) of 300 to400 km s − (Fig. 2A). Similar rotational speeds are observed in the most massive early-type5igure 2: Mass model of ALESS 073.1. ( A ): The observed rotation curve (black circles) adopt-ing the best-fitting inclination of 22 ◦ . The model rotation curve (black solid line) is the total ofcontributions from the stellar bulge (red dot-dashed line), stellar disk (blue dashed line), cold gasdisk (green dotted line), and DM halo (magenta dash-dotted line). ( B ): Surface density profilesof the bulge (red squares), stellar disk (blue stars), and cold gas disk (green circles) adoptingthe best-fitting masses ( ). The bulge profile assumes a de Vaucouleurs’ distribution ( ). Thestellar and gas disk profiles are obtained from azimuthal averages of the dust continuum and[C II ] intensity maps, respectively. The uncertainties on these profiles are dominated by the sys-tematic uncertainty on the normalizing masses ( ). ( C ): [C II ] velocity dispersion profile. ( D ):Disk rotational support versus radius. In (C) and (D), the dashed line and grey band show themedian value and median absolute deviation, respectively. In (A), (C), and (D), the error barscorrespond to ± σ uncertainties ( ). 6alaxies at z = 0 that host molecular gas disks ( ), so we conclude that they are the likelydescendants of galaxies similar to ALESS 073.1. The intrinsic velocity dispersion ( σ V ) reaches50 to 60 km s − at small radii and decreases to ∼
10 km s − in the outer parts (Fig. 2C). Theratio of rotation velocity to velocity dispersion (Fig. 2D) is ∼
10 within 2.5 kpc, similar to thatof cold gas disks at z = 0 ( ). The [C II ] disk of ALESS 073.1 is supported by rotation, notturbulence.The bulge component of the model is necessary to reproduce the observed high rotationspeeds at small radii. We obtain a bulge-to-total mass ratio M bul /M baryon = 0 . , where M baryon is given by the sum of all baryonic (non-DM) components within 3.5 kpc. This ratio is nearlyindependent of the disk inclination because it is largely driven by the shapes of the observedrotation curve and of the baryonic gravitational contributions, not by their absolute normaliza-tions. Baryons dominate the gravitational potential within 3.5 kpc, similar to massive galaxiesat z = 0 ( ) and z = 1 to 3 (
23, 24 ). In Λ CDM cosmology, DM halos are expected to have lowconcentrations at high z , whereas baryons can efficiently cool and collapse to the bottom of thepotential well, reaching high concentrations. Thus, DM halos become gravitationally dominantat large radii.Gas disks in the early Universe are expected to be more turbulent than their local analogues:the gas velocity dispersion is thought to increase with redshift, whereas the degree of rotationalsupport ( V rot /σ V ) decreases (
22, 25 ). The enhanced turbulence may be driven by gravitationalinstabilities from cold gas flows, and/or by stellar and AGN feedback ( ). At z ≥ , however,the limited spatial resolution can bias the measurements of V rot and σ V (beam smearing effects),leading to systematic overestimates of σ V and underestimates of V rot /σ V ( ). The gas velocitydispersion of ALESS 073.1 is consistent with the extrapolations from data at z < . (Fig. 3A)but lies at the low σ V boundary of the prediction of a semi-empirical model based on Toomre’sdisk instability parameter ( ), despite the galaxy hosts a starburst and a dust-obscured AGN.7igure 3: Turbulence and rotational support of galaxies as a function of cosmic time. (A)
Gas velocity dispersion and (B) V rot /σ V versus redshift. Grey bands show a semi-empiricalmodel based on Toomre’s disk instability parameter ( ). In panel A, the black dot-dashedand the blue dashed lines show the velocity dispersion evolution inferred from warm ionizedgas data (black circles; averages from galaxy samples) and cold neutral gas data (blue squares;averages for z = 0 and individual galaxies for z > . ), respectively ( ). In panel B, blackdiamonds and red bars show, respectively, the median and 90 % distribution of ionized-gas sur-veys ( ). The red star represents ALESS 073.1. Error bars correspond to ± σ uncertainties.8he rotational support of the gas disk of ALESS 073.1 is much higher than the value predictedby the same model (Fig. 3B), which has been used to describe kinematic measurements at z < . ( ). The V rot /σ V -ratio is similar to that of gas disks in the local Universe (Fig.3B). This suggests that both starburst and AGN feedback have only a gentle effect on the coldinterstellar medium of ALESS 073.1 because the [C II ] disk is regular and unperturbed, contraryto galaxy formation models with violent feedback.It remains unknown whether bulges exist in the early Universe ( z > ) given the lack ofhigh-resolution infrared images probing the stellar mass distribution. In ALESS 073.1 we usethe gas kinematics to infer the presence of a central massive component that is not traced by thedust. This is likely a stellar bulge hosting a super-massive black hole. Bulges form either fromthe secular evolution of stellar disks or after major galaxy mergers ( ). Both mechanisms couldbe at play in ALESS 073.1, but they must act on short timescales because the Universe was only1.2 billion years old at z (cid:39) . . Any major merger must have happened at even earlier cosmicepochs to provide time for the gas kinematics to relax to regular rotation. The orbital time atthe last measured point of the rotation curve is about × years, and about five revolutionsmay be required to relax the disk, so any major merger must have happened at z > . . ALMAobservations of a quasar-host galaxy at z (cid:39) . with similar spatial resolution show that the[C II ] kinematics can be heavily disturbed at this younger epoch ( ∼ ). Bulges may also develop from starformation inside AGN-driven gas outflows (
29, 30 ), which could occur on short timescales.Although we observe only a single object, we conclude that the Universe produced regularlyrotating galaxy disks with prominent bulges at < % of its current age. This implies that theformation of massive galaxies and their central bulges must be a fast and efficient process.9 cknowledgments We thank L. Legrand and H. ¨Ubler for providing tabular data from their publications. F.L.thanks P. Li for technical support with the Markov-Chain-Monte-Carlo fitting code.
Funding:
F.L. was supported by the ESO fellowship during the initial stages of this project.E.M.D.T. was supported by the U.S. National Science Foundation under grant 1616177. F.F.acknowledges support from the Friedrich Wilhelm Bessel Research Award Programme of theAlexander von Humboldt Foundation. A.W.S.M. is supported by a Dunlap Fellowship at theDunlap Institute for Astronomy & Astrophysics, funded through an endowment established bythe David Dunlap family and the University of Toronto. R.M. acknowledges ERC AdvancedGrant 695671 “QUENCH” and support from the Science and Technology Facilities Council(STFC).
Authors contributions:
F.L.: conceptualization, data curation, formal analysis, investigation,methodology, project administration, software, validation, visualization, writing - original draft.E.M.D.T.: conceptualization, formal analysis, methodology, software, visualization, writing -review & editing. F.F.: conceptualization, methodology, software, writing - review & editing.A.W.S.M.: formal analysis, writing - review & editing. Z.-Y.Z.: formal analysis, writing -review & editing. CdB: data curation, writing - review & editing. T.A.D.: resources, writing -review & editing. R.M.: resources, writing - review & editing.
Competing interests:
The authors have no competing interests.
Data and materials availability:
This paper makes use of the following ALMA data:ADS/JAO.ALMA II ] data cube (Data S1), the[C II ] moment maps (Data S2 and S3), the model outputs from B AROLO (Data S4-S6), themodel outputs from G
ALFIT (Data S7), and the best-fitting mass model (Data S8) are availableas supplementary data files.
Materials and Methods
Observations and data reduction.
ALESS 073.1 was observed on 18 October 2018 withALMA Band-7 and 47 working antennas, giving minimum and maximum baselines of 15and 2517 m, respectively. The weather conditions were good with precipitable water vaporbetween 0.35 and 0.40 mm. The flux, bandpass, and pointing calibrator was ALMA J0522-3627 (B 0521-365), while the phase calibrator was ALMA J0348-2749 (D 0346-279). A high-resolution spectral window with a bandwidth of 1.875 GHz was covered with 480 channels andcentered at 329.638 GHz to target the redshifted [C II ] emission line at a rest-frame frequencyof 1900.539 GHz. Three low-resolution spectral windows with bandwidths of 2 GHz were cov-ered with 128 channels each and centered at 331.201, 317.280, and 319.156 GHz to target thedust-emitting continuum. Two execution blocks of about 45 minutes each gave a total on-sourceintegration time of 1.5 hours. Two additional execution blocks of 45 minutes were flagged as“semi-pass” by the ALMA quality-assurance team. We attempted to calibrate these additionaldata but concluded that they were unavoidably damaged due to time-dependent phase jumps.We also analysed archival data of ALESS 073.1 from 2012 ( ) and 2015 ( ). The 2012 datawere combined with the 2018 data, increasing the overall signal-to-noise ratio with their short11aselines between 18 and 402 m. We used the 2015 data for the continuum imaging only, be-cause the [C II ] emission was filtered out by the very long baselines (up to 16.2 km) of thoseobservations ( ).The data reduction was performed using the Common Astronomy Software Applications(C ASA ) package ( ). The Fourier-plane data were flagged and calibrated using the standardC ASA pipeline. Data from different ALMA cycles were taken with different spectral setups anddefinition of visibility weights. We first split the science target from the Fourier-plane data andused the C
ASA task oldstatwt on the line-free channels to set consistent visibility weightsamong the different datasets. Next, we regridded the line-data to a common bandwidth of about1.6 GHz using 102 channels with a spectral width of 15.626 MHz, corresponding to a velocitywidth of ∼
14 km s − .The combined dataset was imaged using a H¨ogbom deconvolver with Briggs’ robust pa-rameter of 0.5 ( ). A continuum image was obtained combining all four spectral windows,excluding channels with line emission, reaching a sensitivity of 18 µ Jy beam − (where 1 µ Jy =10 − W m − Hz − ). The synthesized beam of the continuum map is . (cid:48)(cid:48) × . (cid:48)(cid:48) with aposition angle (PA) of − . ◦ . A dirty [C II ] cube was obtained without any continuum subtrac-tion and cleaning: these reduction steps were performed in the image plane using the GroningenImaging Processing System (G IPSY ) package ( ). The synthesized beam of the [C II ] cube is . (cid:48)(cid:48) × . (cid:48)(cid:48) with a PA of − . . The root-mean square (rms) noise of the dirty [C II ] cubedepends on frequency because a different number of visibilities from different ALMA cyclesare combined in each rebinned channel. This is taken into account during the cleaning processwith G IPSY .The continuum was subtracted from the [C II ] cube in the image plane, after fitting a first-order polynomial to the line-free channels. The cleaning was performed within a Boolean maskthat follows the kinematic structure of the [C II ] emission. The Boolean mask was constructed12y smoothing the continuum-subtracted [C II ] cube to a resolution of 0.4 (cid:48)(cid:48) and clipping at 3 σ s ,where σ s is the rms noise of the smoothed cube. The cleaning was done down to 3 σ d , where σ d is the rms noise of the continuum-subtracted dirty cube. The cleaned cube was spectrallysmoothed using the Hanning scheme, giving a velocity resolution of ∼
28 km s − (twice thechannel width). This increased the signal-to-noise ratio, giving a mean rms noise of 0.18 mJybeam − per channel. Gas distribution and kinematics.
We assume a Λ CDM cosmology adopting the Planck 2018cosmological parameters ( ): matter density Ω m = 0 . , cosmological constant density Ω Λ =0 . , and Hubble constant H = 67 . km s − Mpc − . In this cosmology, 0.1 (cid:48)(cid:48) corresponds to657.7 pc at the redshift of ALESS 073.1 ( z = 4 . ). Coordinates of Right Ascension (R.A.)and Declination (Dec.) are given using the International Celestial Reference System (ICRS).The [C II ] cube was analysed using the B AROLO package ( ). Moment maps were con-structed using the “smooth & search” option, setting the threshVelocity parameter to 1. Apseudo-3 σ surface brightness contour was estimated following standard procedures ( ). The[C II ] intensity map (moment zero) and velocity field (moment one) are shown in Figure S1(A and D). The [C II ] line is broadened by unresolved differences in rotation within the res-olution element (the beam-smearing effect). B AROLO directly models the 3D cube takingbeam-smearing effects into account, so does not rely on 2D moment maps. B AROLO divides the galaxy into a series of rings with each ring described by nine param-eters: coordinates of the kinematic center ( x , y ), position angle (PA), inclination ( i ), verticalthickness ( z ), systemic velocity ( V sys ), rotation velocity ( V rot ), radial velocity ( V rad ), and ve-locity dispersion ( σ V ). We use nine rings with a width of 0.06 (cid:48)(cid:48) , corresponding to half of theresolution element (cid:113) θ x × θ y , where θ x and θ y are the major and minor axes of the synthesizedbeam. For simplicity, we adopt a fully axisymmetric disk model, so the surface density of13igure S1: Moment Maps: integrated [C II ] intensity ( I [CII] ) map from the observed ( A ), model( B ), and residual data cube ( C ); [C II ] line-of-sight velocity ( V l . o . s . ) map from the observed ( D ),model ( E ), and residual data cube ( F ). North is up and east is left. The kinematic centre (R.A. =03 h m s ; Dec. = − ◦ (cid:48) (cid:48)(cid:48) ) is shown with a white star. The beam size is plottedas the grey ellipse in the bottom-left corner. In panels A and B, iso-emission contours rangefrom 0.35 mJy km s − ( ∼ σ ) to 6 mJy km s − in steps of 0.7 mJy km s − . In panel C, residualcontours are at − σ (blue) and + σ (red); non-axisymmetric, arm-like features are detected. Inpanels D and E, iso-velocity contours range from − to +120 km s − in steps of 30 km s − ;the bold contour highlights the systemic velocity (set to zero). In panel F, the residual contoursare at − and − km s − (blue) as well as +30 and +60 km s − (red). Non-circular motionsabove the velocity resolution ( ∼
30 km s − ) are detected along an arm-like feature towards thenorthwest (see also Fig. S4). 14ach ring is directly computed from the observed [C II ] intensity map using azimuthal averages. B AROLO can also produce disk models with non-axisymmetric gas distributions by renor-malizing the flux density of each spatial pixel to the observed intensity map. This techniquewas introduced to account for rotating disks with massive gas clumps and/or feedback-drivencavities (
37, 38 ) but it appears unnecessary for ALESS 073.1. The non-axisymmetric approachwould be nearly insensitive to the inclination angle because the observed intensity map (mo-ment zero) is always reproduced by construction. We therefore stick to the axisymmetric modelin which the gas surface density depends on radius ( R ) only.We fix the vertical thickness of the gas disk to 300 pc; the value of z has little effect onour results because the disk is nearly face-on. The radial velocity is negligible because thekinematic major axis of the disk is perpendicular to the kinematic minor axis (Fig. S1D), so wefix it to zero. If we leave V rad as a free parameter, we find values consistent with zero within theuncertainties, so we conclude that large-scale radial motions in the disk must be smaller than ∼
30 km s − (approximately the velocity resolution).A first run of B AROLO is used to estimate ( x , y ), PA, and V sys , considering the arithmeticaverage of the best-fitting values across the rings. When all parameters are left free, B AROLO occasionally converges to local rather than global minima giving unphysical solutions. Weverified that the best-fitting values of ( x , y ), V sys , and PA are plausible by visually inspecting[C II ] channel maps, position-velocity (PV) diagrams, and moment maps. The PA is then fixedto 40 ◦ because the velocity field does not show any strong indication for a warped disk within R (cid:39) . (cid:48)(cid:48) . At larger radii, the asymmetric [C II ] extension to the north may potentially bedue to a warp: this would only affect the outermost 2-3 points of the rotation curve, so wekeep the PA constant at all radii rather than modeling an uncertain, asymmetric warp. Thisanomalous kinematic component corresponds to non-circular motions of about 30-60 km s − along an extended arm-like feature (Figs. S1 and S4). Our conclusion of a central massive15ulge is driven by the rotation velocities at small radii, where there is no evidence for a warpeddisk, nor for non-circular motions.A second run of B AROLO is used to estimate an average inclination of ◦ ± ◦ ; theuncertainty is taken to be the standard deviation across the rings. This value of i is confirmedby comparing the model [C II ] intensity map with the observed one, but its formal uncertaintyshould be taken with caution because the [C II ] disk is close to face-on. The inclination angle hasa large effect on the normalization of the best-fitting rotation curve, so we build mass modelstreating i as a free parameter in a Bayesian context, scaling the rotation velocities by sin( i ) and using the value of i = 25 ◦ ± ◦ as a Gaussian prior. A previous kinematic analysis ofALESS 073.1 using the ALMA 2012 data (with a lower spatial resolution of about 0.5 (cid:48)(cid:48) ) gavedifferent inclination values depending on the fitting technique and modeling assumptions ( ): ◦ ± ◦ fitting 2D maps with an exponential disk model, ◦ ± ◦ fitting 2D maps with anarctangent model, and ◦ ± ◦ fitting the 3D cube with the K INMS software ( ). Our 2018observations resolve the [C II ] distribution and kinematics, so we can directly infer the diskinclination without a-priori assumptions on the gas density profile and rotation curve shape.The resulting inclination of ◦ ± ◦ is fully consistent with the previous 3D determination, butis discordant at the 2.5 σ level with the 2D determinations. When analysing poorly resolved data,it is preferable to fit the 3D cube over the 2D moment maps ( ). For example, the asymmetric[C II ] extension to the north of ALESS 073.1 may mimic a more inclined disk in moment mapsat lower spatial resolutions of ∼ (cid:48)(cid:48) (as in the ALMA 2012 data), leading the previous 2Dfitting to infer higher inclinations.The final run of B AROLO uses only V rot and σ V as free parameters. For each parameter, B AROLO provides asymmetric uncertainties ( δ + , δ − ) corresponding to a variation of 5 % ofthe residuals from the global minimum. We consider these errors as 2 σ deviations from thebest-fitting value and compute symmetric 1 σ uncertainties as √ δ + δ − / , which are used in our16igure S2: Kinematic Modeling: position-velocity diagrams along the disk major axis fromthe observed ( A ), model ( B ), and residual ( C ) data cube; position-velocity diagrams along thedisk minor axis from the observed ( D ), model ( E ), and residual data cube ( F ). Horizontal andvertical lines correspond to the galaxy systemic velocity and kinematic center, respectively.Blue contours are at 2 σ , 4 σ , and 8 σ ; grey contours are at -2 σ and -4 σ . In panels A and B, reddots with ± σ error bars show the projected rotation curve V rot sin( i ) .17ayesian mass models (see Eq. 1). Panels B and E of Fig. S1 show moment-zero and moment-one maps from our best-fitting 3D model, while panels C and F show the residual maps. Anaxisymmetric, rotating disk model provides a good match to the data. The residual intensitymap shows significant ( > σ ) and spatially resolved ( > beam) non-axisymmetric structureswith an arm-like morphology. The non-axisymmetric structures to the south are also seen in thedust distribution, as we discuss below (see Fig. S4). The residual velocity field shows that non-circular motions are generally negligible (below our resolution limit of ∼
30 km s − ), except inthe anomalous extension to the northwest. This may be due to recently accreted gas, streamingmotions along spiral arms, or to a warped gas disk. The low signal-to-noise ratio in the outerregions prevents us from determining which.To assess the reliability of our kinematic model, Figure S2 shows PV diagrams along themajor and minor axis of the disk. PV diagrams are a direct representation of the data becausethey do not require masking or clipping, unlike moment maps. The PV diagrams from the modelcube reproduce the observations along both kinematic axes. The majority of the residuals arewithin ± σ , implying that any non-circular motion must be smaller than the velocity resolution( < km s − ), or localized to spatial scales smaller than the ALMA beam ( < <
40 M (cid:12) pc − ). Inspection of the contours on the residual PVdiagrams reveal three small features at − σ along the major axis and one small feature at − σ along the minor axis. The size of these features is similar to our resolution element (both inspace and velocity), thus they are consistent with being noise peaks.Figure S3 shows the [C II ] channel maps from the observed and model cubes, confirming thatthe overall kinematics is reproduced by a simple rotating disk model. The arm-like feature tothe northwest is again an exception, visible at velocities between +103 and − km s − . Similardeviations due to non-circular motions are common in galaxies at z = 0 (
13, 15, 37, 38, 40 ), sodo no represent some peculiar feature of high- z galaxies.18igure S3: Channel Maps: observed [C II ] emission (colorscale and blue contours) overlaidwith the best-fitting disk model (red contours). Coloured contours are at 2 σ , 4 σ , 8 σ . Greycontours are at -4 σ and -2 σ . The beam size is plotted as the grey ellipse in the bottom-leftcorner. Each panel shows data at a different line-of-sight velocity, given in the top-left corner(with respect to the galaxy systemic velocity). The white star indicates the galaxy center.19 ust distribution. The dust distribution (Fig. S4A) appears elongated along a PA of about110 ◦ which is different from the kinematic PA of about 40 ◦ . The asymmetric shape of theALMA beam, however, may play a role because the dust emission is more compact thanthe [C II ] emission, so it is less well resolved. Inspection of the continuum map shows non-axisymmetric structures to the south, distorting the iso-emission contours. We use the G ALFIT package ( ) to fit 2D analytic functions to the dust distribution taking the point-spread function(PSF) into account. Because the ALMA continuum map was cleaned down to the noise level,we consider a Gaussian PSF with a full-width half-maximum equal to the synthesized beam( . (cid:48)(cid:48) × . (cid:48)(cid:48) with PA of − . ◦ ).If we run G ALFIT leaving all parameters free, we do not find satisfactory results due toparameter degeneracies. Our goal is to check whether the dust orientation is consistent withthe [C II ] orientation. Thus, we fix several G ALFIT parameters, adopting an exponential disk(S´ersic index n = 1 ) with the same center and inclination as the [C II ] disk (projected axis ratioof . ). The free parameters are the intrinsic PA, the effective radius, and the integrated flux.Because the dust distribution has substructures, we use an iterative approach masking residualpixels above 5 σ until the fitting results converge ( ). G ALFIT returns
PA = 41 ◦ ± ◦ at the firstiteration and PA = 36 ◦ ± ◦ at the fifth, final iteration (Fig. S4B). These values are consistentwith the kinematic PA. However, if we allow the axis ratio and/or the S´ersic index to vary, awide range of PA values are found, indicating that the intrinsic geometry cannot be constrainedwithout any prior information. We conclude that the orientation of the dust disk is consistentwith that of the [C II ] disk.The residual map (Fig. S4C) includes a prominent arm-like feature to the southwest whichmay connect to a broader feature to the southeast. Similar structures are also visible in theresidual [C II ] intensity map (Fig. S4D). It is unlikely that these features are due to a disruptingmerging companion or to a gas accretion event because the [C II ] kinematics appears largely un-20igure S4: Dust Distribution: observed continuum intensity ( I cont ) map at rest-frame 160 µ m( A ), model dust map ( B ), residual dust map ( C ), and residual dust map overimposed on the [C II ]residual intensity map ( D ). North is up and east is left. The kinematic centre is shown with awhite star. The beam size is plotted as the grey ellipse in the bottom-left corner. In panels A andB, iso-emission contours range from 0.055 ( ∼ σ ) to 1 mJy beam − in steps of 0.11 mJ beam − .In panel C, residual contours are at − σ (blue), + σ and + σ (red); an arm-like feature isvisible. In panel D, the same negative and positive contours are plotted with dashed and solidblack lines, respectively; the residual [C II ] contours are at − σ (blue) and + σ (red); the dustand gaseous arm-like features to the southwest are approximately at the same location.21erturbed at these locations. We interpret these features as a spiral arm embedded in a rotatingdisk. Similar structures have been observed in high-resolution dust maps of other submillime-ter galaxies at high redshifts ( ). Single-arm spiral structures can be excited by weak tidalperturbations due to, e.g., a nearby low-mass companion ( ) or a past galaxy encounter ( ).The [C II ] distribution displays a longer arm-like feature to the northwest (Fig. S4D), where thedust emission is weak or absent and the gas kinematics are less regular (Fig. S1). This maypotentially represent a second arm. Mass Models.
We build a mass model with four components: a cold gas disk, a stellar disk,a stellar bulge, and a DM halo. For each mass component, we compute the expected circularvelocity of a test particle using the
Rotmod task in G
IPSY . We also present alternative modelswith fewer mass components, which have the advantage of having fewer free parameters butare less plausible within the standard galaxy-formation scenario (Figure S6). The fiducial massmodel has six free parameters: disk inclination ( i ), gas mass ( M gas ), stellar disk mass ( M disk ),bulge mass ( M bul ), halo mass ( M ), and halo concentration ( C ).The gravitational contribution of the cold gas is computed by solving Poisson’s equation incylindrical symmetry and considering a disk of finite thickness ( ). We use the azimuthally-averaged radial density profile from the observed [C II ] intensity map (Fig. 2B). This is moreaccurate than assuming simplified functional forms like exponential radial density profiles ( ).For the vertical density distribution, we assume an exponential law with a fixed scale heightof 300 pc ( ). The vertical geometry has a minor effect on the resulting circular velocity inthe disk mid-plane: the difference between a zero-thickness disk and the equivalent sphericaldistribution is at most 20 % ( ). For numerical convenience, the gravitational contribution isnormalized to an arbitrary mass of 10 M (cid:12) , which is then scaled up or down during the fittingprocess using a dimensionless parameter Υ gas .22he gravitational contribution of the stellar disk is computed in the same way as for the gasdisk. ALESS 073.1 has been detected at multiple wavelengths but the stellar component hasnot been spatially resolved. We assume that the young stars are distributed in a similar wayas the dust because the ultraviolet emission from young stars is absorbed by dust grains andre-emitted in the far-infrared. Dust grains at high redshifts, when the Universe was ∼ (cid:12) after 10-100 Myr of stellar evolution and by core-collapse supernovae with lifetimes of about10 Myr ( ). The dynamical timescales ( R/V rot ) in ALESS 073.1 are ∼
50 Myr, so it is likelythat dust and stars did not have enough time to alter their respective radial distributions. Theazimuthally-averaged radial density profile from the ALMA continuum map can be describedby an exponential model with a scale length of 700 pc (Fig. 2B). This scale length is similar tocompact stellar disks in high-surface-brightness galaxies in the nearby Universe ( ), althoughit is likely that the stellar disk of ALESS 073.1 will grow in size from z (cid:39) to z (cid:39) followingthe inside-out galaxy formation scenario ( ). We use the dust radial density profile to solvePoisson’s equation in cylindrical symmetry, adopting an exponential vertical distribution withscale-height of 300 pc. The gravitational contribution is again normalized to a mass of 10 M (cid:12) and rescaled using a dimensionless parameter Υ disk .The gravitational contribution of the stellar bulge is computed assuming spherical symmetry( ) and a De Vaucouleurs’ projected density profile ( ). To constrain the bulge effectiveradius ( R e ), we analyze an image of ALESS 073.1 from the Hubble Space Telescope (HST)in the F160W filter ( ). We use G ALFIT ( ) and construct a PSF using seven nearby stars.The galaxy is unresolved, but we obtain an upper limit on R e of about 300 pc (less than oneHST pixel). The HST image may be affected by the AGN contribution because it measures therest-frame near-ultraviolet emission, so the value of R e is uncertain. Different values of R e donot change our general result that there is a central bulge component, but slightly modify the23est-fitting stellar masses. For example, a smaller value of R e = 200 pc would decrease thebulge mass by about 20 % , while increasing the disk mass by a similar amount. In submillimetergalaxies at z (cid:39) − . , the stellar component is sometimes observed to be more extended thanthe dust component (
53, 54 ). If we assume that the stellar disk has a purely exponential radialprofile with a scale-length larger than that of the dust ( > pc), the stellar-disk contributionto the rotation curve would peak at larger radii, leaving a more massive bulge component inthe central parts. Our fiducial mass model (assuming that the stellar disk is traced by the dustand the stellar bulge follows a De Vaucouleurs’ profile) is conservative regarding the presenceand mass of a central bulge component. Similarly to the gas and stellar disks, we introduce adimensionless scaling parameter Υ bul normalized to a mass of 10 M (cid:12) .The DM gravitational contribution is computed assuming a spherical halo with a Navarro-Frenk-White (NFW) density profile ( ). The halo mass M is measured at R , the radiuswithin which the mass volume density equals 200 times the critical density of the Universe. Forconvenience, we perform our model fitting by converting the halo mass to the circular speedat R , given by V = [10 G H ( z ) M ] / where G is Newton’s gravitational constant and H ( z ) is the Hubble parameter. In our adopted cosmology, H ( z ) is about 7.8 H at z = 4 . .The halo concentration is defined as C = R /R h , where R h is the characteristic scale-length of the NFW profile. MCMC fitting.
The best-fitting parameters are determined using a Markov-Chain-Monte-Carlo (MCMC) method in a Bayesian context. We define the likelihood L = exp( − . χ ) with χ = N (cid:88) i [ V obs ( R i ) − V mod ( R i )] δ V obs ( R i ) , (S1)where V obs is the observed rotation curve at radius R i , δ V obs is the associated error, and V mod isthe model rotation curve obtained by summing in quadrature the velocity contributions of each24ass component, namely: V = Υ gas V + Υ disk V + Υ bul V + V ( C , V ) . (S2)We adopt a Gaussian prior on i with a central value of 25 ◦ and a standard deviation of 3 ◦ ,taken from the modeling with B AROLO . The disk inclination changes the observed rotationvelocities as V obs sin( i ) . For the other parameters, we adopt Gaussian or log-normal priors thatare motivated by previous observations of ALESS 073.1 and theoretical expectations, as wedetail below. Gaussian and log-normal priors in Bayesian statistics ensure posterior propriety( ), but in our case they are also necessary because uninformative, flat priors would lead toparameter degeneracies. Our priors are summarized in Table S1. We impose boundaries tokeep the parameters in a physical range: inclinations between 0 and 90 degrees; gas, disk andbulge masses between 10 and 10 M (cid:12) ; halo velocity between 10 and 1000 km s − ; and haloconcentration between 0.01 and 100.Table S1: MCMC priors.Quantity Prior Type Center Standard Deviation i ( ◦ ) Gaussian 25.00 3.00 log(Υ gas ) Lognormal 0.32 0.50 log(Υ disk ) Lognormal 0.70 0.50 log(Υ bul ) Lognormal 0.70 0.50 log( V ) Lognormal Eq. S3 0.45 C Gaussian 3.40 0.85The cold gas mass ( M gas ) of ALESS 073.1 is constrained by previous observations. Thegalaxy has been detected in the CO( J = 2 → ) line giving a molecular gas mass of (1 . ± . × M (cid:12) ( ). This value assumes thermalized gas with L CO( J =2 → = L CO( J =1 → andthe CO-to-H conversion factor α CO = 0 . M (cid:12) (K km s − pc ) − , appropriate for starburstgalaxies ( ). The mass of atomic gas (hydrogen and helium) can be inferred from the CO( J =2 → ) and [C II ] fluxes using models of photo-dissociation regions (PDR), giving a value of25 . ± . × M (cid:12) (
6, 59 ). Thus, ALESS 073.1 has an estimated cold gas mass (atomic plusmolecular) of . × M (cid:12) . Considering the large uncertainties in both the α CO factor andPDR modeling, we center the prior at log(Υ gas ) = 0 . with a standard deviation of 0.5 dex.The total stellar mass ( M (cid:63) ) of ALESS 073.1 has been estimated by six different teams byfitting the spectral energy distribution (SED) with stellar population synthesis models (
10, 11,60–63 ). ALESS 073.1 has a well-sampled SED from the far-ultraviolet to the far-infrared, butthese studies provide inconsistent estimates of M (cid:63) ranging from 5 to × M (cid:12) , possibly dueto contributions from a starburst and AGN. For example, the AGN contribution in the opticalportion of the SED can mimic an old stellar population. We adopt the arithmetic average of thesix independent determinations ( M (cid:12) ) as our initial guess. This mass may be distributedbetween the disk and bulge components, so we center the disk prior at log(Υ disk ) = 0 . and thebulge prior at log(Υ bul ) = 0 . , assuming a standard deviation of 0.5 dex for both components.Formally, Υ bul includes the mass of the central super-massive black hole, which cannot beindependently measured with the available data. In nearby galaxies, the mass of the centralblack hole is typically 1 % (or less) of the mass of the bulge ( ). In high- z galaxies the situationmight be different, but super-massive black holes have been found to contribute < of thecentral dynamical mass ( ). Thus, the central black hole can have only a small effect on themeasured Υ bul .To constrain the halo mass ( M ), we impose the stellar mass − halo mass relation fromabundance − matching techniques, i.e., by matching the cumulative abundance of predicted DMhalos to that of observed galaxies (
65, 66 ). The M (cid:63) − M relation is very uncertain in the earlyUniverse. We adopt a simple log-linear relation: log( M ) = 0 .
73 log( M (cid:63) ) + 4 . , (S3)which we derive from data at z = 4 . − . ( ). This relation assigns a galaxy with M (cid:63) = × M (cid:12) a DM halo with M = 2 . × M (cid:12) (on average) giving a “condensed”baryon fraction of 0.019, which is about one order of magnitude lower than the cosmic value Ω b / Ω m = 0 . ( ). We assume a scatter of 0.45 dex in halo mass at fixed stellar mass ( ).The halo concentration ( C ) becomes nearly independent of halo mass at z > , thus weassume a Gaussian prior centered at C = 3 . with a standard deviation of 25 % ( ). Withoutimposing these two priors, the halo parameters would be unconstrained. The addition of a DMhalo, however, is necessary to determine whether the observed rotation curve is consistent withthe expectations of the Λ CDM cosmology.The posterior probability distributions of the parameter set are mapped using the affine-invariant ensemble sampler
EMCEE ( ). The MCMC chains are initialized with 200 walkers.Their starting positions are randomly assigned within these ranges for the fit parameters: ◦
MCMC results: the panels show the posterior probability distribution of pairsof mass model fitting parameters, and the marginalized probability distribution of each fittingparameter (histograms). Individual MCMC samples outside the 2 σ confidence region are shownwith black dots, while binned MCMC samples inside the 2 σ confidence region are shown by thegreyscale; the black contour corresponds to the 1 σ confidence region; the red squares and solidlines show median values. In the histograms, solid and dashed lines correspond to the medianand ± σ values, respectively. Υ gas , Υ disk and Υ bul are normalized masses in units of 10 M (cid:12) .28ial mass model. We use CORNER . PY ( ) and calculate 1 σ errors corresponding to the 68%confidence interval of the marginalized 1-dimensional posterior probability distributions. Asexpected, both M disk and M bul are degenerate with i . This occurs because the overall normal-ization of the rotation curve is set by / sin( i ) . For our inclination prior of ◦ ± ◦ (from the B AROLO models), the final best-fitting inclination is . ◦ ± . ◦ , giving rotation velocitiesbetween 300 and 400 km s − and pushing the baryonic masses towards the lower 1 σ boundaryof the respective priors. Figure S5 shows that M disk and M bul are only weakly degenerate witheach other, implying that the ratio M bul /M disk is well constrained and nearly independent of i . Alternative Mass Models.
We build more conservative mass models that have three gravita-tional components instead of four, excluding (in turn) the stellar bulge, the stellar disk, the gasdisk, and the DM halo. Some of these models (e.g. the one without a DM halo) are statisticallyfavored over our fiducial model because they provide a similarly accurate fit with a smaller num-ber of free parameters. They are, however, less plausible within the standard Λ CDM paradigmof galaxy formation. For mass models without the stellar disk or stellar bulge, the Bayesianpriors are revised such that the total fiducial stellar mass (10 M (cid:12) ) is preserved.The model without the bulge (Fig. S5A) is unsatisfactory because the inner rotation veloc-ities cannot be reproduced. The MCMC fitting tries to compensate for the lack of the bulgecomponent by increasing the gravitational contribution of the stellar disk, as well as decreasingthe total rotation velocities by pushing the inclination angle up to 23 ◦ . This overshoots the rota-tion curve at large radii. A massive central component, which is not traced by the dust emission,is therefore necessary to explain the observed kinematics of ALESS 073.1.The model without stellar disk (Fig. S5B) provides a better fit but does not capture theobserved shape of the rotation curve. This model assumes that the star-formation activity, tracedby the dust emission, has just started in the disk and no substantial stellar component has been29igure S6: Alternative Mass Models: the observed rotation curve (black dots) is comparedwith mass models fitted by assuming no bulge ( A ), no stellar disk ( B ), no gas ( C ), and no darkmatter halo ( D ). The best-fitting model (black, solid line) in each case is the sum of gravitationalcontributions from (when present): stellar bulge (dot-dashed, red line), stellar disk (dashed, blueline), cold gas disk (green, dotted line), and DM halo (dash-dotted, magenta line).30uilt yet. We regard this as unlikely because the estimated star-formation rate of 1000 M (cid:12) yr − can build a stellar disk of 10 M (cid:12) in 10 years. This model gives M bul /M baryon = 0 . whichis higher than M bul /M baryon = 0 . from our fiducial model because the relative contributionof the bulge stays almost the same as it is fixed by the inner rotation curve, while the relativecontributions of gas and DM increase to explain the outer rotation velocities. M bul /M baryon remains > . confirming that a prominent bulge is in place at z (cid:39) .The model without gas disk (Fig. S5C) is almost identical to the fiducial mass model be-cause the subdominant gas contribution is compensated by a slightly heavier stellar disk. Thisgives M bul /M baryon = 0 . which is equal to our fiducial value within the uncertainties (TableS2). We regard this model as unphysical because [C II ] and CO( J = 2 → ) emission lines aredetected in ALESS 073.1.The model without DM halo (Fig. S5D) is similar to our fiducial model but has a lowervalue of M bul /M baryon = 0 . . This indicates that DM within 3.5 kpc is not strictly necessary.This model is equivalent to a rotation curve in the context of Modified Newtonian Dynamics(MOND), which replaces DM with a departure from the classical laws of Newtonian dynamicsat accelerations smaller than a (cid:39) − m s − ( ). Assuming that the acceleration scale a does not vary with z and that the redshift-distance relationship in a MOND cosmology is similarto that in Λ CDM ( ), the observed accelerations in the outer regions of ALESS 073.1 are about7 times higher than a , so the galaxy is expected to behave as a classic Newtonian system withno mass discrepancy. Captions for Data
Data S1: observed [C II ] data cube after cleaning and continuum subtraction. File in FlexibleImage Transport System (FITS) format with a single layer. Data S2: observed [C II ] intensity map (moment zero). File in FITS format with a single layer31 ata S3: observed [C II ] velocity map (moment one). File in FITS format with a single layer. Data S4: model [C II ] data cube from B AROLO . File in FITS format with a single layer.
Data S5: model [C II ] intensity map from B AROLO . File in FITS format with a single layer.
Data S6: model [C II ] velocity map from B AROLO . File in FITS format with a single layer.
Data S7: continuum maps from G
ALFIT . File in FITS format with four layers: Ext[1] withthe observed continuum map, Ext[2] with the model continuum map, Ext[3] with the residualcontinuum map. The primary layer (Ext[0]) is empty.
Data S8: best-fitting mass model. File in machine-readable table (mrt) format.
References
1. S. D. M. White, M. J. Rees,
Mon. Not. R. Astron. Soc. , 341 (1978).2. A. Dekel, Y. Birnboim,
Mon. Not. R. Astron. Soc. , 2 (2006).3. M. Vogelsberger, et al. , Nature , 177 (2014).4. J. Kormendy, J. Kennicutt, Robert C.,
Annu. Rev. Astron. Astrophys. , 603 (2004).5. E. Ba˜nados, et al. , Astrophys. J. , L23 (2019).6. C. De Breuck, et al. , Astron. Astrophys. , A59 (2014).7. G. C. Jones, et al. , Astrophys. J. , 180 (2017).8. A. Pensabene, et al. , Astron. Astrophys. , A84 (2020).9. Materials and methods are available as supplementary materials.10. K. E. K. Coppin, et al. , Mon. Not. R. Astron. Soc. , 1905 (2009).11. R. Gilli, et al. , Astron. Astrophys. , A67 (2014).322. R. Gilli, et al. , Astrophys. J. , L28 (2011).13. R. Sancisi, F. Fraternali, T. Oosterloo, T. van der Hulst,
Astron. Astrophys. Rev. , 189(2008).14. S. M. Sweet, et al. , Mon. Not. R. Astron. Soc. , 5700 (2019).15. F. Fraternali, T. Oosterloo, R. Sancisi, G. van Moorsel,
Astrophys. J. , L47 (2001).16. G. de Vaucouleurs,
Annales d’Astrophysique , 247 (1948).17. E. M. Di Teodoro, F. Fraternali, Mon. Not. R. Astron. Soc. , 3021 (2015).18. S. Casertano, J. H. van Gorkom,
Astron. J. , 1231 (1991).19. E. Noordermeer, J. M. van der Hulst, R. Sancisi, R. S. Swaters, T. S. van Albada,
Mon. Not.R. Astron. Soc. , 1513 (2007).20. F. Lelli, S. S. McGaugh, J. M. Schombert,
Astron. J. , 157 (2016).21. T. A. Davis, et al. , Mon. Not. R. Astron. Soc. , 214 (2016).22. E. Wisnioski, et al. , Astrophys. J. , 209 (2015).23. F. Lelli, et al. , Mon. Not. R. Astron. Soc. , 5440 (2018).24. R. Genzel, et al. , Nature , 397 (2017).25. H. ¨Ubler, et al. , Astrophys. J. , 48 (2019).26. M. D. Lehnert, et al. , Astrophys. J. , 1660 (2009).27. E. M. Di Teodoro, F. Fraternali, S. H. Miller,
Astron. Astrophys. , A77 (2016).28. B. P. Venemans, et al. , Astrophys. J. , L30 (2019).339. W. Ishibashi, A. C. Fabian,
Mon. Not. R. Astron. Soc. , 1474 (2014).30. R. Maiolino, et al. , Nature , 202 (2017).31. B. Gullberg, et al. , Astrophys. J. , 12 (2018).32. J. P. McMullin, B. Waters, D. Schiebel, W. Young, K. Golap,
Astronomical Data Anal-ysis Software and Systems XVI , R. A. Shaw, F. Hill, D. J. Bell, eds. (2007), vol. 376 of
Astronomical Society of the Pacific Conference Series , p. 127.33. D. S. Briggs,
American Astronomical Society Meeting Abstracts (1995), vol. 187 of
Amer-ican Astronomical Society Meeting Abstracts , p. 112.02.34. M. G. R. Vogelaar, J. P. Terlouw,
Astronomical Data Analysis Software and Systems X ,F. R. Harnden, Jr., F. A. Primini, H. E. Payne, eds. (2001), vol. 238 of
Astronomical Societyof the Pacific Conference Series , p. 358.35. Planck Collaboration, et al. , Astron. Astrophys. , A6 (2020).36. F. Lelli, M. Verheijen, F. Fraternali,
Mon. Not. R. Astron. Soc. , 1694 (2014).37. F. Lelli, M. Verheijen, F. Fraternali, R. Sancisi,
Astron. Astrophys. , A72 (2012).38. F. Lelli, M. Verheijen, F. Fraternali, R. Sancisi,
Astron. Astrophys. , A145 (2012).39. T. A. Davis, M. Bureau, M. Cappellari, M. Sarzi, L. Blitz,
Nature , 328 (2013).40. C. Trachternach, W. J. G. de Blok, F. Walter, E. Brinks, J. Kennicutt, R. C.,
Astron. J. ,2720 (2008).41. C. Y. Peng, L. C. Ho, C. D. Impey, H.-W. Rix,
Astron. J. , 2097 (2010).42. J. A. Hodge, et al. , Astrophys. J. , 130 (2019).343. E. Athanassoula,
Astron. Astrophys. , 395 (1978).44. M. Thomasson, et al. , A&A , 25 (1989).45. S. Casertano,
Mon. Not. R. Astron. Soc. , 735 (1983).46. J. A. Sellwood,
Galaxy Dynamics - A Rutgers Symposium , D. R. Merritt, M. Valluri, J. A.Sellwood, eds. (1999), vol. 182 of
Astronomical Society of the Pacific Conference Series ,p. 351.47. P. C. van der Kruit, K. C. Freeman,
Annu. Rev. Astron. Astrophys. , 301 (2011).48. E. Noordermeer, Mon. Not. R. Astron. Soc. , 1359 (2008).49. C. Gall, J. Hjorth, A. C. Andersen,
Astron. Astrophys. Rev. , 43 (2011).50. G. Pezzulli, F. Fraternali, S. Boissier, J. C. Mu˜noz-Mateos, Mon. Not. R. Astron. Soc. ,2324 (2015).51. S. M. Kent,
Astron. J. , 1301 (1986).52. N. A. Grogin, et al. , Astrophys. J. Suppl. Ser. , 35 (2011).53. J. A. Hodge, et al. , Astrophys. J. , 103 (2016).54. P. Lang, et al. , Astrophys. J. , 54 (2019).55. J. F. Navarro, C. S. Frenk, S. D. M. White,
Astrophys. J. , 563 (1996).56. H. Tak, S. K. Ghosh, J. A. Ellis,
Mon. Not. R. Astron. Soc. , 277 (2018).57. K. E. K. Coppin, et al. , Mon. Not. R. Astron. Soc. , L103 (2010).58. L. J. Tacconi, et al. , Astrophys. J. , 246 (2008).359. C. De Breuck, et al. , Astron. Astrophys. , L8 (2011).60. D. P. Stark, A. J. Bunker, R. S. Ellis, L. P. Eyles, M. Lacy,
Astrophys. J. , 84 (2007).61. J. L. Wardlow, et al. , Mon. Not. R. Astron. Soc. , 1479 (2011).62. J. M. Simpson, et al. , Astrophys. J. , 125 (2014).63. T. Wiklind, et al. , Astrophys. J. , 111 (2014).64. J. Kormendy, L. C. Ho,
Annu. Rev. Astron. Astrophys. , 511 (2013).65. A. Vale, J. P. Ostriker, Mon. Not. R. Astron. Soc. , 189 (2004).66. Q. Guo, S. White, C. Li, M. Boylan-Kolchin,
Mon. Not. R. Astron. Soc. , 1111 (2010).67. L. Legrand, et al. , Mon. Not. R. Astron. Soc. , 5468 (2019).68. A. A. Dutton, A. V. Macci`o,
Mon. Not. R. Astron. Soc. , 3359 (2014).69. D. Foreman-Mackey, D. W. Hogg, D. Lang, J. Goodman,
PASP , 306 (2013).70. D. Foreman-Mackey,
The Journal of Open Source Software , 24 (2016).71. M. Milgrom, Astrophys. J. , 365 (1983).72. C. Skordis, T. Złosnik, arXiv e-printsarXiv e-prints