EMPRESS. IX. Extremely Metal-poor Galaxies are Very Gas-rich Dispersion-dominated Systems: Will the James Webb Space Telescope Witness Gaseous Turbulent High-z Primordial Galaxies?

We present kinematics of six local extremely metal-poor galaxies (EMPGs) with low metallicities (0.016–0.098 Z ⊙) and low stellar masses (104.7–107.6 M ⊙). Taking deep medium/high-resolution (R ∼ 7500) integral-field spectra with 8.2 m Subaru, we resolve the small inner velocity gradients and dispersions of the EMPGs with Hα emission. Carefully masking out substructures originating by inflow and/or outflow, we fit three-dimensional disk models to the observed Hα flux, velocity, and velocity dispersion maps. All the EMPGs show rotational velocities (v rot) of 5–23 km s−1 smaller than the velocity dispersions (σ 0) of 17–31 km s−1, indicating dispersion-dominated (v rot/σ 0 = 0.29–0.80 < 1) systems affected by inflow and/or outflow. Except for two EMPGs with large uncertainties, we find that the EMPGs have very large gas-mass fractions of f gas ≃ 0.9–1.0. Comparing our results with other Hα kinematics studies, we find that v rot/σ 0 decreases and f gas increases with decreasing metallicity, decreasing stellar mass, and increasing specific star formation rate. We also find that simulated high-z (z ∼ 7) forming galaxies have gas fractions and dynamics similar to the observed EMPGs. Our EMPG observations and the simulations suggest that primordial galaxies are gas-rich dispersion-dominated systems, which would be identified by the forthcoming James Webb Space Telescope observations at z ∼ 7.


INTRODUCTION
Primordial galaxies characterized by low metallicities and low stellar masses are important to understand galaxy formation.Cosmological and hydrodynamical simulations (e.g., Wise et al. 2012;Yajima et al. 2022) have predicted that primordial galaxies at z ∼ 10 would form in low-mass halos with halo masses of ∼ 10 8 M .Such simulated primordial galaxies at z 7 show low gas-phase metallicities of 1% of the solar abundance (Z ), specific star-formation rates (sSFRs) of ∼ 100 Gyr −1 , and low stellar masses of 10 6 M .Despite their scientific relevance, it is difficult to observe primordial galaxies due to their faintness.Isobe et al. (2022;hereafter Paper IV) have estimated an Hα flux of a primordial galaxy with M * = 10 6 M as a function of redshift, and demonstrated that even the current-best spectrographs such as Keck/LRIS or Keck/MOSFIRE cannot observe the primordial galaxy at z 0.5 without any gravitational lensing magnification (e.g., Kikuchihara et al. 2020).
Kinematics of high-z primordial galaxies can provide us a hint of what kind of mechanism (e.g., inflow/outflow) would impact on the early galaxy formation.We are still lacking a good handle on the detailed gas dynamical state, often quantified in broad terms by the relative level of rotation, via the v rot /σ 0 ratio (e.g., Förster Schreiber et al. 2009).Integral-field unit (IFU) observations (e.g., Wisnioski et al. 2015;Herrera-Camus et al. 2022) have reported that star-forming galaxies show v rot /σ 0 ratios decreasing from v rot /σ 0 ∼ 10 to ∼ 2 with increasing redshift from z ∼ 0 to 5.5, while such observations are currently limited to massive (∼ 10 10 M ) galaxies.It is important to directly determine whether low-mass ( 10 6 M ) primordial galaxies are truly dominated by dispersion.
Complementary to the high-z kinematics studies, some studies have reported v rot /σ 0 values of local galax-ies with lower stellar masses.Local galaxies are advantageous for conducting deep observations with high spectral and spatial resolutions.The SHαDE survey (Barat et al. 2020) has made a significant progress in extracting v rot and σ 0 values of local dwarf galaxies with stellar masses down to ∼ 10 6 M .These observations suggest that the ratio v rot /σ 0 decreases down to 1 with decreasing M * down to ∼ 10 6 M .However, the SHαDE galaxies have gas-phase oxygen metallicities1 higher than 12 + log(O/H) ∼ 7.69 that correspond to Z ∼ 0.1 Z (Asplund et al. 2021), which are still 1 dex higher than that of the simulated high-z primordial galaxy (Wise et al. 2012).
Recently, a project "Extremely Metal-Poor Representatives Explored by the Subaru Survey (EMPRESS)" has been launched (Paper I).EMPRESS aims to select faint EMPG photometric candidates from Subaru/Hyper Suprime-Cam (HSC; Miyazaki et al. 2018) deep optical (i lim = 26 mag; Aihara et al. 2019) images, which are ∼ 2 dex deeper than those of SDSS.Conducting follow-up spectroscopic observations of the EMPG photometric candidates, EMPRESS has identified new 12 EMPGs with low stellar masses of 10 4.2 -10 6.6 M (Papers I, IV; Nakajima et al. 2022, hereafter Paper V;Xu et al. 2022, hereafter Paper VI). Remarkably, J1631+4426 has been reported to have a metallicity of 0.016 Z , which is the lowest metallicity idenitifed so far (Paper I).
Including the 12 low-mass EMPGs found by EM-PRESS, Paper V summarizes 103 local EMPGs identified so far whose metallicities are accurately measured with the direct-temperature method (e.g., Izotov et al. 2006).The 103 EMPGs show low metallicities of 0.016-0.1 Z , low stellar masses of ∼ 10 4 -10 8 M , and high sSFRs of ∼ 1-400 Gyr −1 (Paper V).These features resemble the simulated primordial galaxy at z 7 (Wise et al. 2012), suggesting that EMPGs would be good local analogs of high-z primordial galaxies (but see also Isobe et al. 2021, hereafter Paper III).
This paper is the ninth paper of EMPRESS, reporting Hα kinematics of EMPGs observed with Subaru/Faint Object Camera and Spectrograph (FOCAS) IFU (Ozaki et al. 2020) in a series of the Subaru Intensive Program entitled EMPRESS 3D (PI: M. Ouchi).So far, EM-PRESS has released 8 papers related to EMPGs, each of which reports the survey design (Paper I), high Fe/O ratios suggestive of massive stars (Kojima et al. 2021, hereafter Paper II; Paper IV), morphology (Paper III), low-Z ends of metallicity diagnostics (Paper V), outflows (Paper VI), the shape of incident spectrum that reproduces high-ionization lines (Umeda et al. 2022, hereafter Paper VII), and the primordial He abundance (Matsumoto et al. 2022, hereafter Paper VIII).
This paper is organized as follows.Section 2 explains our observational targets.The observations are summarized in Section 3. Section 4 reports Hα flux, velocity, and velocity-dispersion maps of our targets.Our kinematic analysis is described in Section 5. Section 6 lists our results.We discuss kinematics of primordial galaxies in Section 7. Our findings are summarized in Section 8. Throughout this paper, magnitudes are in the AB system (Oke & Gunn 1983).We assume a standard ΛCDM cosmology with parameters of (Ω m , Ω Λ , H 0 ) = (0.3, 0.7, 70 km s −1 Mpc −1 ).The solar metallicity Z is defined by 12 + log(O/H) = 8.69 (Asplund et al. 2021).
We note that the 4 EMPGs other than J1631+4426 or J2115−1734 have radio and/or optical integral-field spectroscopy conducted by previous studies.IZw18 has H i observations with VLA (Lelli et al. 2012).SBS0335−052E also has VLA H i observations (Pustilnik et al. 2001) as well as VLT/MUSE Hα observations with a spectral resolution of R ∼ 3000 (Herenz et al. 2017).Our observations with FOCAS IFU are com-  1.The 6 EMPGs have low metallicities of 12 + log(O/H) = 6.90-7.68,low stellar masses of log(M * /M ) = 4.7-7.6,and high specific star-formation rates of log(sSFR/Gyr −1 ) = 0.5-2.7.These properties indicate that the 6 EMPGs are the best available analogs of primordial galaxies in the early universe.We emphasize that the 6 EMPGs include J1631+4426 having 12 + log(O/H) = 6.90 (Paper I), which is the lowest gas-phase metallicity among galaxies identified so far.

Observations and Data Reduction
This section reports our spectroscopic observations with FOCAS IFU (Ozaki et al. 2020).FOCAS IFU is an IFU with an image slicer installed in FOCAS (Kashikawa et al. 2002) mounted on a Cassegrain focus of the Subaru 8.2-m telescope.The large diameter of the Subaru Telescope allows FOCAS IFU to perform deep integral-field spectroscopy with a 5σ limiting flux of ∼ 1 × 10 −17 erg s −1 cm −2 arcsec −2 (Ozaki et al. 2020) 2 .The FoV of FOCAS IFU is 13. 5 (slice length direction; hereafter X direction) × 10. 0 (slice width di-rection; hereafter Y direction).The pixel scale in a reduced data cube is 0. 215 (X direction) and 0. 435 (Y direction).
We carried out spectroscopy with FOCAS IFU for the 6 EMPGs (Section 2).The observing nights were 2021 August 13, November 24, and December 13.We set pointing positions so that the whole structures of EMPGs are covered with single FoVs.
We used the mid-high-dispersion grism of VPH680 offering the spectral resolution of R ∼ 7500.We took comparison frames of the ThAr lamp.We observed Feige67, HZ44, BD28, Feige110, and Feige34 as standard stars.All the nights were clear with seeing sizes of ∼ 0. 7.
We use a reduction pipeline software of FOCAS IFU3 based on PyRAF (Tody 1986) and Astropy (Astropy Collaboration et al. 2013).The software performs bias subtraction, flat-fielding, cosmic-ray cleaning, sky subtraction, wavelength calibration with the comparison frames, and flux calibration.We estimate flux errors containing read-out noises and photon noises of sky and object emissions.
It should be noted that there remain systematic velocity differences among the slices even after the wavelength calibration with the comparison frames.We conducted additional wavelength calibration with sky emission lines around observed wavelengths of Hα.

Flux, Velocity, and Velocity Dispersion Measurement
At all the spaxels, we fit Hα lines using a single Gaussian function (+ constant) to measure Hα fluxes, (relative) velocities, and velocity dispersions.We derive a The X (Y) direction denotes the slice length (width) direction.The blue solid lines describe the position of the single slice.The light coming into the slice (blue shaded) is projected onto the CCD, where the Y direction in the sky is parallel to the wavelength (λ) direction on the CCD.The original spatial flux profile projected onto the CCD (black dashed curve) provides fluxes of (i − 1), i, and (i + 1)th slices (fi−1, fi, and fi+1, respectively).The original flux profile is assumed to be approximated by the quadratic function f (λ) = aλ 2 + bλ + c (red solid curve) within the λ range from −2d to 2d.We estimate the systematic wavelength shift originating from the slit-width effect (λ shift ) from a barycenter of the flux following f (λ) intercepted by the slice.line-spread function (LSF) by measuring widths of sky lines.A typical sky line has a full-width half maximum (FWHM) of ∼ 0.8-1 Å.Assuming that both the Hα lines from the EMPGs and the sky lines agree well with the Gaussian function, we obtain the intrinsic velocity dispersions by subtracting the LSFs from the observed velocity dispersions quadratically.We confirm that the assumption is reasonable for most of the Hα lines and the sky lines other than some turbulent regions with multiple velocity components (Section 4).
We estimate the errors of the velocities and velocity dispersions by running the Monte Carlo simulations similar to the procedure of Herenz et al. (2016).We make 100 mock data cubes from our data cubes perturbed by the noise data cubes, and measure the velocities and velocity dispersions from each mock data cube.We regard the standard deviations of the derived velocities and velocity dispersions as the errors.The errors include the uncertainties of the LSFs, the additional wavelength calibration with the sky lines (Section 3.1), and the slitwidth effect correction (see Section 3.3), because we perform these corrections using each mock data cube.Typical uncertainties of the velocity and velocity dispersion in each spaxel are ∼ 1.7 and ∼ 1.6 km s −1 , respectively.We self-consistently obtain the errors of the kinematic parameters in Sections 5.2 and 5.3 based on the 100 mock data cubes, in the same manner as we measure the errors of the velocities and velocity dispersions.

Slit-width Effect Correction
It should be noted that the observed velocities suffer from systematic wavelength shifts known as the slitwidth effect (e.g., Bacon et al. 1995).Figure 1 illustrates the mechanism of the slit-width effect.The artificial wavelength shift is caused by a flux gradient in the Y direction (in each slice) that follows the dispersion (wavelength) axis (λ direction thereafter) on the CCD.
Assuming a 2D Gaussian function for the (spatial) flux profile (i.e., assuming a point source), Bacon et al. (1995) have derived the wavelength shifts originating from the slit-width effect.Considering that all our targets are extended and complex sources, we expanded the method to more flexible flux profiles.The wavelength shift λ shift can be estimated from a barycenter of the flux intercepted by the slice as follows: where w and f (λ) are a slice width projected onto the CCD and a flux profile in the Y direction, respectively.Because the pseudo slit width is sampled by 4 pixels for each slice, w in the unit of Å is given by w = 4d, where d is a dispersion in the unit of Å pixel −1 .We approximate f (λ) by a quadratic function f (λ) = aλ 2 + bλ + c so that the function form is determined by the 3 points (λ i−1 , f i−1 ), (λ i , f i ), and (λ i+1 , f i+1 ) in Figure 1, where f i−1 , f i , and f i+1 represent fluxes of (i − 1)th, ith, and (i + 1)th slices, respectively.In this case, λ shift can be calculated as follows: We derive a, b, and c from the equations below: Using Equations 2 and 3, we obtain We infer from Equation 4 that high-dispersion dispersers make the slit-width effect weak.VPH680 has only  To test our slit-width effect correction, we observed SBS0335−052E with 2 frames whose position angles are different by 180 degree.We expect that an average of the 2 velocity maps obtained from the 2 frames cancel out the slit-width effect.Panel (d) of Figure 2 shows the average map, and panel (e) represents residuals of the corrected map and the average map.The residuals are smaller than the velocity shifts caused by the slit-width effect.Figure 3 compares residuals of the observed velocity map and the average map (i.e., residuals NOT corrected for the slit-width effect; left) and residuals of the corrected velocity map and the average map (i.e., residuals corrected for the slit-width effect; right).
Figure 3 shows that the residuals get much smaller after the slit-width effect correction.The small residuals of ∼ ±5 km s −1 also indicate that the corrected map agrees well with the average map.Thus, we conclude that our slit-width correction works well.
We evaluate the uncertainty of the slit-width effect correction by performing Monte Carlo simulation based on flux errors.A typical value of the uncertainty is 0.3 km s −1 .clump and the EMPG tail are not in the same kinematic structure.The blue clump itself shows a weak velocity gradient of ∼ 10 km s −1 from east to west, while the velocity dispersion is relatively high (∼ 25 km s −1 ).#2 IZw18: IZw18 has 2 main blue clumps of IZw18 Northwest (NW) and IZw18 Southeast (SE), both of which have been confirmed to be EMPGs (Izotov & Thuan 1998).IZw18NW is blueshifted by ∼ 40 km s −1 with respect to IZw18SE (at around their flux peaks), consistent with the H i gas kinematics (Lelli et al. 2012).We thus think that IZw18NW and IZw18SE are not in the same kinematic structure.IZw18NW shows a complex morphokinematic structure.In the Hα flux map, we find that the arc-like structure indicated by the gray curve (Hα arc; Dufour & Hester 1990) has a velocity similar to that of IZw18NW, which indicates that the Hα arc and IZw18NW belong to the same kinematic structure.In the velocity map, we also identify a velocity structure that is redshifted by ∼ 20 km s −1 with respect to the flux peak of IZw18NW (red circle).The velocity structure has a relatively high velocity dispersion of ∼ 30 km s −1 , which implies that the structure is not settled.IZw18NW itself does not show a clear bulk rotation (i.e., not likely to be dominated by rotation).
#3 SBS0335−052E: The gri image shows that SBS0335−052E consists of a large blue clump (thick cyan arrow) and a western subclump (thin cyan arrow).In the velocity map, we confirm that SBS0335−052E has a redshifted (∼ +20 km s −1 ) region at northeast (red circle) and a blueshifted (∼ −30 km s −1 ) region at northwest (blue circle).The redshifted and blueshifted regions seem connected to the northeast and north filaments identified with the wide-field MUSE Hα observations (Herenz et al. 2017), respectively.Both regions have relatively high velocity dispersions of ∼ 60 km s −1 , which agree with the scenario that the 2 structures are created by outflows (Herenz et al. 2017).The southern area of SBS0335−052E generally shows a relatively low velocity dispersion of ∼ 30 km s −1 , which indicates that the southern area is relatively settled.The southern area shows a bulk velocity gradient of ∼ 20 km s −1 from northwest to southeast.
#4 HS0822+3542: HS0822+3542 consists of a blue clump (cyan arrow) and a very diffuse EMPG tail elongated from southeast to northwest (white arrow).The Hα map indicates that the major star formation occurs in the blue clump.We find a very weak velocity gradient of ∼ 5 km s −1 in the blue clump from north to south.HS0822+3542 shows a high velocity dispersion of ∼ 20 km s −1 , which is comparable to the previous observations (Moiseev et al. 2010).#5 J1044+0353: The deep HSC image shows that J1044+0353 is composed of a western giant blue clump (thick cyan arrow), an eastern small blue clump (thin cyan arrow), and 3 white clumps between the 2 blue clumps (white arrow)4 .We regard the 3 white clumps as the EMPG tail.The Hα map indicates that the major star formation occurs in the giant blue clump.We obtain similar velocity and velocity dispersion maps to Moiseev et al. (2010)'s result, while we show more clearly that the small blue clump emits Hα.We find that the 2 blue clumps are redshifted by ∼ 20 km s −1 with respect to the EMPG tail, which indicate that the 2 blue clumps and the EMPG tail are not in the same kinematic structure.We find that the giant blue clump has a weak velocity gradient of ∼ 20 km s −1 from southeast to northwest, while the velocity dispersion (∼ 30 km s −1 ) is quite high.
#6 J2115−1734: J2115−1734 consists of a blue clump (cyan arrow) and a north-east EMPG tail (white arrow).The majority of star formation is likely to occur in the blue clump.The EMPG tail is redshifted by ∼ 40 km s −1 with respect to the blue clump.Although the observed velocity field continuously changes between the blue clump and the EMPG tail, J2115−1734 shows a relatively high velocity dispersion of ∼ 40 km s −1 originating from 2 different velocity components.The blue clump itself has a weak velocity gradient of ∼ 20 km s −1 from northeast to southwest with a velocity dispersion of ∼ 30 km s −1 .
In summary, all 6 EMPGs look irregular and dominated by dispersion rather than rotation.Their rotation velocities are not likely to be significant.In Section 5, we analyze the kinematics of the EMPGs more quantitatively.

Masking
Below, we quantify the level of rotation of the EMPGs by assuming that those systems are represented by a single rotating disk.The dynamical center of the disk is thought to be located at the Hα flux peak of the main blue clump (IZw18NW for IZw18) under the assumptions of the gas-mass dominance on the galactic scale (cf.e.g., Herrera-Camus et al. 2022) and the empirical positive correlation between the gas mass surface density Σ gas and the SFR surface density Σ SFR (a.k.a.Kennicutt-Schmidt law; Kennicutt 1998).We confirm that these assumptions are reasonable by deriving the mass profile of each EMPG (Section 6.2).In the following analysis, we use only the region within the Kron radius of the main blue clump to remove the contamination from EMPG tails (IZw18SE for IZw18).We also mask spaxels with velocity dispersions higher than the flux-weighted 84th percentile of the distribution of the velocity dispersions because such regions are thought to be highly turbulent (Egorov et al. 2021).The white regions of Figure 5 show the masked regions.

Non-parametric Method
We calculate a shearing velocity v shear of each EMPG, which is a non-paramteric kinematic property widely used in the literature (e.g., Law et al. 2009;Herenz et al. 2016).We derive v shear from where v max (v min ) is the 95th (5th) percentile of the velocity distribution (Herenz et al. 2016).We note that high v shear values do not necessarily imply the existence of rotation because different velocity components in the complex velocity fields can mimic rotation patterns at small scales.In this sense, v shear can be regarded as an upper limit of the rotation velocity.The global velocity dispersion can be quantified by the flux-weighted median of the distribution of the velocity dispersions (σ med ).
We calculate the errors of the v shear and σ med values (∆ v shear and ∆ σ med , respectively) by running the Monte Carlo simulations explained in Section 3.2.We list v shear and σ med of the 6 EMPGs in

Parametric Method
We conduct detailed dynamical modeling with the 3D parametric code GalPaK 3D (Bouché et al. 2015).Constraining morphological and kinematic properties simultaneously from 3D data cubes, GalPaK 3D provides the deprojected maximum rotational velocity (v rot ) that is irrespective of the inclination.GalPaK 3D also calculates σ 0 free from the velocity dispersions driven by the self gravity (e.g., Genzel et al. 2008) or mixture of the lineof-sight velocities.Because GalPaK 3D does not support rectangular spaxels, we divide the spaxels in the Y direction based on the linear interpolation so that the pixel scale in the Y direction is the same as that in the X direction.We have 10 free parameters of XY coordinates of the disk center, systemic redshift, flux, inclination (i), position angle, effective radius of Hα (r e,Hα ), turnover radius (r t ), v rot , and intrinsic velocity dispersion σ 0 .We assume that all 6 EMPGs have thick disks with rotation curves of the arctan profiles: where r is a radius from the dynamical center.We note that i is determined by the axis ratio and the disk height, where the disk height of the thick-disk model is fixed to 0.15r e,Hα (Bouché et al. 2015).We also assume that the surface-brightness (SB) profiles follow Sérsic profiles with Sérsic indices n = 1, which are inferred from iband SB profiles (Paper III).We have confirmed that these assumptions do not change v rot and σ 0 much.We also use a point-spread function (PSF) obtained from standard stars with a Moffat profile whose FWHM and power index are 0. 7 and 2.5, respectively.Figure 6 summarizes fitting results.Table 2 lists kinematic properties extracted by the GalPaK 3D analysis.We estimate the errors of the v rot and σ 0 values by carrying out the Monte Carlo simulations (Section 3.2).We note that v rot can be regarded as an upper limit of the actual rotation velocity because the small-scale velocity differences in the complex velocity field can mimic rotation patterns (see also Section 5.2).

Mass Profile
The best-fit disk models obtained in Section 5.3 provide an estimate of the radial profiles for dynamical masses M dyn .The dynamical mass M dyn is expected to be a sum of stellar mass M * , gas mass M gas , darkmatter (DM) mass M DM , and dust mass M dust .Note that M dust is negligible in all 6 EMPGs because of their low E(B −V ) values (e.g., Paper I). Hereafter, we derive radial profiles of M dyn , M * , M gas , and M DM .We derive M dyn enclosed within the radius r from the equation

Dynamical mass profile
under the assumption of the virial equilibrium.Figure 7 presents mass profiles of all 6 EMPGs.The red curves correspond to M dyn (< r).Dynamical masses within r e are listed in Table 4.
We use the spectral energy distribution (SED) interpretation code of beagle (Chevallard & Charlot 2016).
The beagle code calculates both the stellar continuum and the nebular emission using the stellar population synthesis code (Bruzual & Charlot 2003) and the nebular emission library of Gutkin et al. (2016) that are computed with the photoionization code cloudy (Ferland et al. 2013).We fit the SED models to SDSS ugrizband photometries (DR16; Ahumada et al. 2020).We run the beagle code with 4 free parameters of maximum stellar age t max , stellar mass M * , ionization parameter U , and V -band optical depth τ V whose parametric ranges are the same as those adopted in Papers III and IV.We fix metallicities Z to be the gas-phase metallicities 12 + log(O/H) listed in Table 1.We also assume a constant star-formation history and the Chabrier (2003) IMF in the same manner as Papers III and IV.
We note that we conduct the SED fitting with almost the same setting as that of Paper I for J1631+4426 and J2115−1734, except for Z and τ V .In Paper I, the Z is treated as a free parameter, and the τ V is fixed to be 0, while it has no impact on M * of J1631+4426 and J2115−1734 because both 2 EMPGs are dust-poor (Paper I).We also confirm that the stellar masses of the 6 EMPGs are different only by 0.3 dex from those derived from i-band magnitudes based on the same massto-light ratio.We include this systematic error of 0.3 dex into the total error of the M * profile.
To obtain M * profiles, we assume that azimuthallyaveraged M * distributions of the EMPGs follow Sérsic profiles.J1631+4426 has an i-band effective radius (r e,i ) of 137 +9 −7 pc and an i-band Sérsic index (n i ) of 1.08 +0.15 −0.13 (Paper III).Because the i-band surface-brightness distribution is expected to trace the M * distribution well (Paper III), we assume that J1631+4426 has a stellar effective radius (r e, * ) and a stellar Sérsic index (n * ) equal to r e,i and n i , respectively.We also assume that the other 5 EMPGs have r e, * within the range from r e,Hα /2 to r e,Hα because r e,i of J1631+4426 is ∼ 2 times smaller than r e,Hα of J1631+4426 (see Table 3)5 .We confirm that the assumed r e, * of the 5 EMPGs are comparable to r e,i of EMPGs (Paper III).The 5 EMPGs are assumed to have n * within the range from 0.7 to 1.7 inferred from the typical n i value of EMPGs (Paper III).We include these uncertainties of r e, * and n * into the error of the M * profile.The yellow shaded regions in Figure 7 represent cumulative M * profiles with their uncertainties.Stellar masses within r e are listed in Table 4.

Gas mass profile
We obtain M gas distributions from the Hα flux distributions, using the Kennicutt-Schmidt law (Section 5.1).However, observational studies (e.g., Shi et al. 2014) have reported that some EMPGs have Σ gas values ∼ 1 dex larger than those inferred from the Kennicutt-Schmidt law, which can be interpreted as the lack of metals suppressing efficient gas cooling and succeeding star formation (e.g., Ostriker et al. 2010;Krumholz 2013).Given the uncertainty of the Kennicutt-Schmidt law at the low-metallicity end, we add this 1 dex upper error to the original scatter of the Kennicutt-Schmidt law (∼ 0.3 dex; Kennicutt 1998).The cyan shaded regions in Figure 7 indicate cumulative M gas profiles with their uncertainties.We note that the H i observations of Lelli et al. (2012) and Pustilnik et al. (2001) have reported that IZw18NW and SBS0335−052E have M gas ∼ 1 × 10 8 and ∼ 1 × 10 9 M within wide scales of ∼ 0.2 and ∼ 3 kpc, respectively, which are consistent with the extrapolations of the M gas profiles that we derive.Gas masses within r e are listed in Table 4.We obtain gas mass fractions at r e,Hα from the following equation f gas = M gas (< r e,Hα )/[M gas (< r e,Hα ) + M * (< r e,Hα )].

Dark-matter mass profile
We estimate M DM profiles under the assumption that the M DM (density) profile follows the NFW profile   ( Navarro et al. 1996).At the virial radius r 200 within which the spherically-averaged mass density is 200 times the critical density ρ c , the NFW profile can be described by (10) We derive M 200 from an empirical stellar-to-halo mass ratio of dwarf galaxies (Brook et al. 2014) of We note that the scatter of the observed stellar-to-halo mass ratios of dwarf galaxies around Equation 11 is ∼ 0.3 dex (Prole et al. 2019).We include this uncertainty into the error of the M DM profile.We also obtain c 200 from M 200 using a halo mass-concentration relation for the Planck cosmology (Dutton & Macciò 2014): where h = H 0 /(100 km s −1 Mpc −1 )6 .The gray shaded regions in Figure 7 denote M DM profiles with their uncertainties.DM masses within r e are listed in Table 4.
The top left panel of Figure 8 shows v rot /σ 0 values of the 6 EMPGs (red circle) as a function of metallicity.We compare our results with other surveys of Hα kinematics of star-forming dwarf galaxies (SHαDE; Barat et al. 2020) and star-forming galaxies (SAMI; Barat et al. 2019;DYNAMO;Green et al. 2014), whose metallicities are drawn from the SDSS MPA-JHU catalog (Tremonti et al. 2004;Brinchmann et al. 2004).The SHαDE, SAMI, and DYNAMO galaxies have 12+log(O/H) ∼ 8-9.We find that v rot /σ 0 decreases with decreasing 12 + log(O/H).The top middle and top right panels of Figure 8 show v rot /σ 0 as a function of M * and sSFR, respectively.The SFRs are derived from Hα luminosity.The M * and sSFR potentially have uncertainties of ∼ 0.3 dex under different assumptions such as initial mass functions (cf.Section 5.4.2).Although it should be noted that the EMPGs are biased toward lower metallicities, we also find that v rot /σ 0 decreases as M * decreases and sSFR increases.These results suggest that galaxies in earlier stages of the formation phase may have lower v rot /σ 0 .

Mass Fraction
Figure 7 summarizes radial profiles of M dyn (red), M * (yellow), M gas (cyan), and M DM (black).We find that the 4 EMPGs other than IZw18NW or SBS0335−052E have the M * profiles ∼ 2 dex below the M dyn profiles within radii up to several times r e,Hα , which means that the 4 EMPGs are dominated by M gas or M DM on galactic scales.On the other hand, IZw18NW and SBS0335−052E have the M * profiles comparable to the M dyn profiles.We also find that the M dyn profiles of all 6 EMPGs can be explained by the M gas profiles within the uncertainties (see Section 5.4.3).We confirm that the M DM profiles of all 6 EMPGs are ∼ 1 dex below the M dyn profiles within radii up to several times r e,Hα .We thus conclude that the masses of the 4 EMPGs except for IZw18NW and SBS0335−052E are dominated by M gas on galactic scales.We note that IZw18NW and SBS0335−052E indeed have large M gas values of ∼ 1 × 10 8 and ∼ 1 × 10 9 M inferred from the H i observations within ∼ 0.2 and ∼ 3 kpc, respectively (Section 5.4.3).Within these larger scales, we can say that both IZw18NW and SBS0335−052E are gas-rich (i.e., f gas ∼ 1).This conclusion is consistent with those in Pustilnik et al. (2020aPustilnik et al. ( ,b, 2021) ) based on H i observations.It should also be noted that Equations 10 and 11 suggest that EMPGs with ∼ 10 6 M have M 200 ∼ 7 × 10 9 M at r 200 ∼ 30 kpc, whereas we can only observe the area within at most ∼ 1 kpc.Therefore, it is natural that M DM has negligible effects on the mass profile in this study.
The middle left panel of Figure 8 shows f gas of the 6 EMPGs as a function of metallicity.Except for IZw18NW and SBS0335−052E with large f gas uncertainties, we find that the EMPGs are gas-rich with high f gas values of ∼ 0.9-1.0.Comparing with the literature, we find that galaxies with lower metallicities tend to have higher f gas values.The f gas values also increase at lower M * and higher sSFR.These results indicate that galaxies in the earlier formation phase would have higher f gas , which are consistent with previous observations (e.g., Maseda et al. 2014) as well as models of galaxy formation based on the ΛCDM model (e.g., Geach et al. 2011).It should be noted that M gas of the SHαDE and DYNAMO galaxies are estimated from  4), we exclude the data point of this EMPG from the bottom panels.
Kennicutt-Schmidt law (i.e., correlation between Σ SFR and Σ gas ) in the similar way as our analysis (cf.Section 5.4.3).Thus, it is natural that we find a tight correlation between sSFR and f gas .In Section 6.1, we report that all 6 EMPGs have low v rot /σ 0 < 1.We also find that v rot /σ 0 decreases with decreasing 12+log(O/H), decreasing M * , and increasing sSFR (Figure 8).Below, we investigate well-discussed three contributors (e.g., Glazebrook 2013;Barat et al. 2020) to such a low v rot /σ 0 < 1.

Thermal expansion
The first possible contributor to the low v rot /σ 0 is the thermal expansion of H ii regions (e.g., Krumholz et al. 2018;Fukushima & Yajima 2021, 2022).We estimate a velocity dispersion of the thermal expansion (σ th ) from the line-of-sight component of the Maxwellian velocity distribution (σ th = kT e /m; e.g., Chávez et al. 2014;Pillepich et al. 2019), where k, T e , and m represent the Boltzmann constant, the electron temperature, and the hydrogen mass, respectively.We obtain σ th = 9.1 km s −1 under the assumption of T e = 10000 K, which is consistent with the typical electron temperature of O ii (i.e., H ii) regions of the EMPGs (e.g., Paper I). Subtracting σ th from σ 0 quadratically, we obtain velocity dispersions being free from the thermal line broadening (σ no therm = 14-30 km s −1 ).We confirm that v rot /σ no therm values are still lower than unity (0.31-0.84) for all 6 EMPGs.We thus conclude that the thermal expansion cannot explain v rot /σ 0 < 1.

Merger/inflow
The second possible contributor is merger/inflow events (e.g., Glazebrook 2013).The merger (inflow) can raise velocity dispersions by tidal heating (releasing potential energies of infalling gas).This scenario can explain all the trends seen in Figure 8 because we can ex-pect that both gas-rich minor merger and inflow would supply metal-poor gas and trigger succeeding starbursts.Especially for J1631+4426, IZw18NW, SBS0335−052E, and J2115−1734, we find the velocity differences from the EMPG tails (IZw18SE for IZw18NW) suggestive of merger (Section 4).

Stellar feedback
The third possible contributor is the stellar feedback (e.g., Lehnert et al. 2009).This includes the supernova (SN) feedback (Dib et al. 2006) as well as stellar winds and the radiative pressure from young massive stars (Mac Low & Klessen 2004).Outflowing gas from SNe and/or young massive stars would raise velocity dispersions.This stellar feedback scenario can directly explain the decreasing trend between v rot /σ 0 and sSFR (the top right panel of Figure 8).Given the expectation that young galaxies (i.e., with high sSFRs) would have low 12 + log(O/H) and M * , this scenario can indirectly reproduce the trend that v rot /σ 0 decreases with decreasing 12 + log(O/H) and M * .IZw18NW, SBS0335−052E, and J1044+0353 especially show outflow signatures in flux, velocity, and velocity-dispersion maps (see Section 4), which imply the dominance of the stellar feedback.We thus conclude that the stellar feedback can also be one of the main contributors of v rot /σ 0 < 1 at the lowmetallicity, low-M * , and high-sSFR ends.Cosmological zoom-in simulations will provide a hint of what kind of feedback is the main contributor.

Toomre Q Parameter
To compare with other kinematic studies, we derive the Toomre Q parameter (Toomre 1964) of the EMPGs.In general, if the Q value of a rotating disk is greater than unity (i.e., Q > 1), the disk is thought to be gravitationally stable.On the other hand, the disk is gravitationally unstable if Q < 1.However, note that it is unclear if this criterion is applicable for the EMPGs because they may not have rotating disks (Section 4).An average of Q within a disk so called the global Q (e.g., Aumer et al. 2010) is calculated from the equation where a is a parameter with values ranging from 1 to 2 depending on the gas distribution (Genzel et al. 2011).
Here, we adopt a = √ 2 assuming the constant rotational velocity.
Table 4 lists the global Q of the 6 EMPGs.We find that all the 6 EMPGs show Q > 1, albeit one of the 6 EMPGs (IZw18NW) with the large Q uncertainty.The bottom panels of Figure 8 illustrate the global Q of the EMPGs as a function of metallicity (bottom left), M * (bottom center), and sSFR (bottom right), while we exclude IZw18NW due to its large Q uncertainty.We also find that the global Q increases with decreasing 12 + log(O/H), decreasing M * , and increasing sSFR.However, the large global Q values are inconsistent with the large sSFR values because star-formation activities are not likely to become aggressive in a gravitationallystable disk.
This inconsistency probably originates from the fact that the global Q parameter is not a reliable indicator for gas-rich galaxies (Romeo et al. 2010;Romeo & Agertz 2014).Instead, it can be important to investigate if the EMPGs lie on a tight relation based on observables such as H i angular momentum (Romeo 2020;Romeo et al. 2020).We need high-resolution H i observations to discuss gravitational instability of EMPGs precisely.

Connection to High-z Primordial Galaxies
The trend that v rot /σ 0 (f gas ) decreases (increases) with decreasing (increasing) 12 + log(O/H) and M * (sSFR) suggests that primordial galaxies at high redshifts would be dispersion-dominated gas-rich systems.This suggestion agrees with the decreasing trend of v rot /σ 0 with the redshift reported by both observational (e.g., Wisnioski et al. 2015) and simulation (Pillepich et al. 2019) studies.
Here, we investigate kinematics of simulated primordial galaxies at high redshifts.In this study, we choose a z = 7.3 primordial galaxy of Wise et al. (2014)'s cosmological radiation hydrodynamics simulation because the simulated galaxy has a low gas-phase metallicity (∼ 4% Z ), a low stellar mass (3.8×10 6 M ), a large f gas (∼ 1), a large sSFR (∼ 5 Gyr −1 ), and a small half-mass radius (∼ 200 pc), all of which are comparable7 to those of the EMPGs (Sections 2 and 5.4.2).Wise et al. in prep. (hereafter W23) extract Hα flux, velocity, and velocitydispersion maps from the simulated galaxy.The FoV of the extracted region is 3.71 kpc with the spatial resolution of 3.71/252 = 0.015 kpc pixel −1 .The spectral resolution of the datacube is 0.1 Å. W23 use cloudy (Ferland et al. 2013) to derive Hα fluxes from the following 5 quantities: Hydrogen number density, temperature, metallicity, incident radiation intensity, and H i column density, along with enough parametric ranges and number of datapoints (1.8 million points in total).For each cell in the data cube of the 5 quantities, W23 linearly interpolate the emissivity table in 5D.To calculate the velocity map, W23 sum the line-of-sight veloc- ities weighted by the H-alpha fluxes.To calculate the velocity dispersion map, W23 use the same method as Pillepich et al. (2019;Equation A3), using the H-alpha flux as the weights.We note that the velocity dispersions include the thermal broadening.
We coarsen the data cube to a spatial resolution of 180 pc, similar to the resolution of our observations, to make a more relevant comparison.The top left, top right, and bottom left panels of Figure 9 are Hα flux, velocity, and velocity-dispersion maps of the simulated galaxy, respectively.We find that the simulated galaxy has an irregular morphology with multiple kinematic sub-structures and localized turbulent regions.These features can be seen in the EMPGs as well (Figure 4).Masking out the kinematic sub-structures and the turbulent regions, we also derive kinematics properties of v shear and σ med in the same manner as we do for the EMPGs (Sections 5.1 and 5.2).We find that the simulated galaxy has a low v shear (8.5 km s −1 ) and a high σ med (18.2 km s −1 ), which are comparable to those of the EMPGs (v shear = 5.5-14.3km s −1 , σ med = 16.9-30.8km s −1 ; see Section 6.1).Consequently, the simulated galaxy has a low v shear /σ med = 0.47 below unity suggesting that the simulated galaxy is dominated by dispersion as well as the EMPGs.Given that the local EMPGs, analogs of high-z primordial galaxies, and the simulated high-z primordial galaxies are both dispersion-dominated systems, we expect that high-z primordial galaxies are likely to dispersion-dominated galaxies.
The forthcoming James Webb Space Telescope (JWST) can directly investigate Hα kinematics of highz primordial galaxies.Paper IV has simulated Hα fluxes of primordial galaxies with M * = 10 6 M at redshifts ranging from 0 to 10. Comparing the Hα fluxes with the limiting flux of JWST/NIRSpec, we estimate that NIR-Spec can detect Hα fluxes of primordial galaxies with M * = 10 6 M at z 1 without gravitational lensing.With a ∼ 2 dex magnification from gravitational lens-ing, NIRSpec can detect Hα fluxes of primordial galaxies with M * = 10 6 M at z ∼ 7. We infer from this estimation that NIRSpec could observe primordial galaxies with M * ∼ 10 7 M at z ∼ 7 with a realistic magnification of ∼ 1 dex.The top middle panel of Figure 8 shows that most of galaxies with M * ∼ 10 7 M already have low v rot /σ 0 < 1 and high f gas > 0.5.Lensing cluster surveys using JWST such as GLASS (PI: T. Treu) and CANUCS (PI: C. Willott) will potentially pinpoint low-mass galaxies with M * ∼ 10 7 M at z ∼ 7, and follow-up IFU observations with JWST may identify dispersion-dominated gas-rich galaxies.
Taking deep medium-high resolution (R ∼ 7500) integral-field spectra with 8.2-m Subaru (Section 3), we resolve the small inner velocity gradients and dispersions of the EMPGs with Hα emission.Carefully masking out sub-structures originated by inflow and/or outflow, we fit 3-dimensional disk models to the observed Hα flux, velocity, and velocity-dispersion maps (Sections 4 and 5).All the EMPGs show rotational velocities (v rot ) of 5-23 km s −1 smaller than the velocity dispersions (σ 0 ) of 17-31 km s −1 , indicating dispersion-dominated systems with small ratios of v rot /σ 0 = 0.29−0.80(Section 6.1) that can be explained by turbulence driven by inflow and/or outflow (Section 7.1).Except for two EMPGs with large uncertainties, we find that the EMPGs have very large gas-mass fractions of f gas 0.9 − 1.0 (Section 6.2).Comparing our results with other Hα kinematics studies, we find that v rot /σ 0 (f gas ) decreases (increases) with decreasing metallicity.We compare numerical simulations of first-galaxy formation and identify that the simulated high-z (z ∼ 7) forming galaxies have gas-fractions and dynamics similar to the observed EMPGs.Our EMPG observations and the simulations suggest that primordial galaxies are gas-rich dispersion-dominated systems, which would be identified by the forthcoming JWST observations at z ∼ 7 (Section 7.3).
We are grateful to Yoshiaki Sofue, Alessandro Romeo, and Roberto Terlevich for having useful discussions.We also thank the staff of Subaru Telescope for their help with the observations.This research is based on data collected at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan (NAOJ).We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical, and natural significance in Hawaii.The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University.The HSC instrumentation and software were developed by the NAOJ, the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan

Figure 1 .
Figure1.Explanation of the slit-width effect.The inset panel shows a schematic figure of an object on the celestial sphere.The X (Y) direction denotes the slice length (width) direction.The blue solid lines describe the position of the single slice.The light coming into the slice (blue shaded) is projected onto the CCD, where the Y direction in the sky is parallel to the wavelength (λ) direction on the CCD.The original spatial flux profile projected onto the CCD (black dashed curve) provides fluxes of (i − 1), i, and (i + 1)th slices (fi−1, fi, and fi+1, respectively).The original flux profile is assumed to be approximated by the quadratic function f (λ) = aλ 2 + bλ + c (red solid curve) within the λ range from −2d to 2d.We estimate the systematic wavelength shift originating from the slit-width effect (λ shift ) from a barycenter of the flux following f (λ) intercepted by the slice.

Figure 2 .
Figure 2. Confirmation of the slit-width effect correction using SBS0335−052E.(a) Observed velocity map.(b) Velocity shift generated by the slit-width effect.(c) Velocity map corrected for the slit-width effect.(d) Average of the 2 velocity maps whose position angles are different by 180 degree.(e) Residual of the corrected map and the average map.The gray contours illustrate SN ratios in the order of 5, 10, 20, ... from the outside.

km s - 1 Figure 3 .
Figure 3. (Left) Residual of the observed velocity map and the average map in Figure 2. (Right) Residual of the corrected map and the average map, i.e., same as panel (e) in Figure 2, but with the smaller velocity range.

Figure 4 .
Figure 4 summarizes Hα flux, velocity, and velocity dispersion maps of the 6 EMPGs as well as gri images cut out of the FoV of FOCAS IFU.We note that only 2 EMPGs of #1 and 5 (i.e., J1631+4426 and J1044+0353) have images of the HSC-Subaru Strategic Program (SSP) Public Data Release (PDR) 3 (Aihara et al. 2022).The deep HSC images of the 2 EMPGs are shown in Figure 4, while the other gri images are taken from the Pan-STARRS catalog (Flewelling et al. 2020).We report morphological and kinematic features of each EMPG below, checking the consistency with previous IFU studies.#1 J1631+4426: The HSC gri image illustrates that J1631+4426 consists of a blue clump (indicated by the cyan arrow) and a white diffuse structure elongated from north to south (white arrow).We refer to the white structure as the "EMPG tail" (Paper III).The Hα flux map shows that the Hα flux of J1631+4426 is dominated by the blue clump, which indicates that star formation in J1631+4426 should mainly occur in that region.Paper I has confirmed that the blue clump has the very low metallicity of 12 + log(O/H) = 6.90 (Section 2) based on the direct-temperature method, while the EMPG tail can have a metallicity even lower than the blue clump based on the [O iii]λ5007/Hα ratio (Kashiwagi et al. 2021).The velocity map shows that the blue clump is discontinuously redshifted by ∼ 20 km s −1 with respect to the EMPG tail, which indicates that the blue

Figure 5 .
Figure 5. Mask map.The white regions correspond to the masked regions (Section 5.1).

Figure 6 .
Figure 6.GalPaK 3D fitting of J1631+4426.The top, middle, and bottom panels show Hα flux, velocity, and velocitydispersion maps, respectively.The left, center, and right panels illustrate observed, model, and residual (= observed − model) maps, respectively.Note that the velocity maps in this figure denote relative velocities with respect to the systemic redshift obtained with GalPaK 3D (Section 5.3).

Figure 7 .
Figure7.Enclosed mass profiles of the 6 EMPGs.The red, yellow, cyan, and black curves represent dynamical, stellar, gas, and dark-matter mass profiles, respectively.The vertical dotted lines show re,Hα.The edge of the plots correspond to the outer most radii used for the kinematic analysis.
c 200 ) − c 200 /(1 + c 200 ) (9) and x = r/r 200 .The concentration parameter c 200 is defined as the ratio r 200 /r d controlling where the slope equals −2 as it changes from −3 at large radii to a central value of −1.We obtain r 200 from the virial mass M 200 using the equation r = 3M 200 4π(200ρ c ) 1/3 .

Figure 8 .
Figure 8. Diagrams of vrot/σ0, fgas, and global Q as functions of 12 + log(O/H), M * , and sSFR.The red circles represent the EMPGs, while the blue squares, the green diamonds, and the black circles indicate the SHαDE (Barat et al. 2020), SAMI(Barat et al. 2019), and DYNAMO galaxies(Green et al. 2014), respectively.Because IZw18NW has the large uncertainty of the global Q (see Table4), we exclude the data point of this EMPG from the bottom panels.
Origin of Low v rot /σ 0

Table 2 .
The medians of our v shear /∆ v shear and σ med /∆ σ med values are 22 and 73, respectively.These values are comparable to those of Herenz et al. (2016)'s observations, whose spectral resolutions and SN ratios of Hα are similar to those of our observations.Herenz et al. (2016)'s v shear /∆ v shear ∼ 22 and σ med /∆ σ med ∼ 108, respectively.

Table 2 .
Kinematic properties of the 6 EMPGs