# Possible Systematic Rotation in the Mature Stellar Population of a $z=9.1$
Galaxy

Tsuyoshi Tokuoka  Department of Pure and Applied Physics, Graduate School of
Advanced Science and Engineering,
Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku,
Tokyo 169-8555, Japan  Akio K. Inoue Department of Pure and Applied Physics,
Graduate School of Advanced Science and Engineering,
Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku,
Tokyo 169-8555, Japan  Waseda Research Institute for Science and Engineering,
Faculty of Science and Engineering, Waseda University,
3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan<EMAIL_ADDRESS>Takuya
Hashimoto Tomonaga Center for the History of the Universe (TCHoU), Faculty of
Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571,
Japan  Richard S. Ellis  Department of Physics & Astronomy, University College
London, Gower Street, London WC1E 6BT, United Kingdom  Nicolas Laporte Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge CB3 0HE, UK  Yuma Sugahara Waseda Research Institute for
Science and Engineering, Faculty of Science and Engineering, Waseda
University,
3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan  National Astronomical
Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan  Hiroshi
Matsuo National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka,
Tokyo 181-8588, Japan  Yoichi Tamura Division of Particle and Astrophysical
Science, Graduate School of Science, Nagoya University, Aichi 4648602, Japan
Yoshinobu Fudamoto  Waseda Research Institute for Science and Engineering,
Faculty of Science and Engineering, Waseda University,
3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan  National Astronomical
Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan  Kana
Moriwaki Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo,
Tokyo 113-0033, Japan  Guido Roberts-Borsani Department of Physics and
Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los
Angeles, CA 90095, USA  Ikkoh Shimizu  Department of Literature, Shikoku
Gakuin University, 3-2-1 Bunkyocho, Zentsuji, Kagawa 765-8505, Japan  Satoshi
Yamanaka General Education Department, National Institute of Technology, Toba
College, 1-1, Ikegami-cho, Toba, Mie 517-8501, Japan  Naoki Yoshida Department
of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033,
Japan  Kavli Institute for the Physics and Mathematics of the Universe (WPI),
UT Institutes for Advanced Study, The University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8583, Japan  Research Center for the Early Universe, School
of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033,
Japan  Institute for Physics of Intelligence, School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Erik
Zackrisson Observational Astrophysics, Department of Physics and Astronomy,
Uppsala University, Box 516, SE-751 20 Uppsala, Sweden  Wei Zheng Department
of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

We present new observations with the Atacama Large Millimeter/submillimeter
Array for a gravitationally-lensed galaxy at $z=9.1$, MACS1149-JD1. [O iii]
88-$\mu$m emission is detected at 10$\sigma$ with a spatial resolution of
$\sim 0.3$ kpc in the source plane, enabling the most distant morpho-kinematic
study of a galaxy. The [O iii] emission is distributed smoothly without any
resolved clumps and shows a clear velocity gradient with $\Delta V_{\rm
obs}/2\sigma_{\rm tot}=0.84\pm 0.23$, where $\Delta V_{\rm obs}$ is the
observed maximum velocity difference and $\sigma_{\rm tot}$ is the velocity
dispersion measured in the spatially-integrated line profile, suggesting a
rotating system. Assuming a geometrically thin self-gravitating rotation disk
model, we obtain $V_{\rm rot}/\sigma_{V}=0.67_{-0.26}^{+0.73}$, where $V_{\rm
rot}$ and $\sigma_{V}$ are the rotation velocity and velocity dispersion,
respectively, still consistent with rotation. The resulting disk mass of
$0.65_{-0.40}^{+1.37}\times 10^{9}$ M⊙ is consistent with being associated
with the stellar mass identified with a 300 Myr-old stellar population
independently indicated by a Balmer break in the spectral energy distribution.
We conclude that the most of the dynamical mass is associated with the
previously-identified mature stellar population that formed at $z\sim 15$.

Galaxy dynamics (591) — Galaxy evolution (594) — Galaxy formation (595) —
High-redshift galaxies (734)

##  1 Introduction

The Atacama Large Millimeter/submillimeter Array (ALMA) has revolutionized
high-redshift galaxy observations, allowing galaxies to be characterized well
into the epoch of reionization. For example, dust continuum as well as the [O
iii] 88- and [C ii] 158-$\mu$m emission lines of galaxies at $z>7$ have been
successfully observed (e.g., Watson et al. 2015; Inoue et al. 2016; Hashimoto
et al. 2019). In particular, MACS1149-JD1 (Zheng et al., 2012; Bouwens et al.,
2014; Huang et al., 2016; Zheng et al., 2017; Hoag et al., 2018) is a
gravitationally-lensed galaxy emitting the [O iii] line at $z=9.1$, one of the
most distant objects spectroscopically confirmed (Hashimoto et al. 2018,
hereafter H18). This galaxy also shows a Balmer break consistent with a
stellar population of a few hundred Myrs old, suggesting its formation epoch
is $z\sim 15$ (H18, Binggeli et al. 2019; Roberts-Borsani et al. 2020; Laporte
et al. 2021).

Beyond finding high-redshift galaxies, studying their dynamics based on the
kinematics of their interstellar medium provides further motivation for
probing the early physics of galaxy formation. Such studies have been mostly
conducted with three-dimensional (3D) near-infrared spectroscopy for galaxies
at $z<4$ (e.g., Förster Schreiber et al. 2009; Jones et al. 2010; Wisnioski et
al. 2015). However, the high sensitivity and high spatial and frequency
resolution of ALMA also make it possible to analyze morpho-kinematics of
galaxies at $4<z<6$ (Rizzo et al., 2020, 2021; Lelli et al., 2021) and even at
$z\sim 7$ (Smit et al., 2018). In this Letter, we present the most distant
example of a morpho-kinematic analysis of the [O iii] emission in MACS1149-JD1
at $z=9.1$ and discuss when the rotational motion in galaxies first appears.

Throughout this Letter, we use a flat $\Lambda$CDM cosmology with a parameter
set of $(h,\Omega_{m},\Omega_{\Lambda})=(0.704,0.272,0.728)$ (Komatsu et al.,
2011).

##  2 Observational data

New observations for [O iii] 88-$\mu$m emission from MACS1149-JD1 were
performed in Band 7 during ALMA Cycle 6 (2018.1.00616.S, PI: T. Hashimoto) to
improve the spatial resolution of the emission line. The antenna
configurations were C43-4, -5, and -6 (minimum baseline = 15.1 m and maximum
baseline = 783.5–2516.9 m). We set a spectral window (SPW) centered at a
335.625-GHz emission line frequency (H18) with a bandwidth of 1.875 GHz and
240 channels, corresponding to a velocity resolution of 7.0 km s-1. Three
other SPWs were set at the central frequencies of 337.375, 347.417, and
349.176 GHz to obtain the continuum emission, with bandwidths of 2.000 GHz and
128 channels. The observations were executed in a series of 13 sets between
October 18 and December 15, 2018. The precipitable water vapor during the
observations spread over 0.3–1.0 mm, and the mean was 0.5 mm. The total on-
source exposure time was 9.6 h, compared to 2.0 h in the previous observations
(H18). Raw data were processed using Common Astronomy Software Applications
(CASA; McMullin et al. 2007) version 5.4.0-68, Pipeline version 42030M.

##  3 Data analysis

###  3.1 Overview of new results

First, we created a dust continuum image from all data using CASA task tclean,
which resulted in a null-detection. A new $3\sigma$ upper limit on the dust
continuum is 19 $\mu$Jy beam-1, a factor of 3 improvement over the previous
limit (H18). Next, we created a data cube from the Cycle 6 SPWs supposed to
contain the [O iii] line using tclean with natural weighting and 50-km s-1
velocity binning. We successfully confirmed the [O iii] line in this
independent dataset. The redshift is $z=9.1111\pm 0.006$, which is consistent
with $z=9.1096\pm 0.006$ found by H18. Combining all SPWs of H18 and Cycle 6
data, which contain the [O iii] line, we also created a data cube using tclean
with natural weighting and 50-km s-1 velocity binning. We use this dirty
imaging data cube throughout this Letter.111 We use the dirty beam to convolve
the dynamical model in §4, avoiding any biases induced by modeling of the
clean beam for the clean data cube, which is not straightforward for the
combined datasets of different array configurations.

We created a velocity-integrated intensity map, i.e. moment0 map, of the [O
iii] line (Figure 1, top left) using CASA task immoments from the data cube
within a velocity range of $-150$ to $+200$ km s-1 relative to the line
redshift $z=9.1096$ (H18). We also extracted the total line spectrum (Figure
1, top right) from the data cube integrated over the area where the line
emission was detected at $>3\sigma$ in the moment0 map. The velocity
dispersion of the total line profile was measured at $\sigma_{\rm tot}=72.7\pm
8.1$ km s-1, also consistent with $65.4\pm 16.6$ km s-1 in H18.

In the moment0 map, the peak signal-to-noise ratio increased from $7.4\sigma$
(H18) to $10\sigma$, and the beam full width at the half maximum (FWHM)
improved from $0.^{\prime\prime}62\times 0.^{\prime\prime}52$ (H18) to
$0.^{\prime\prime}39\times 0.^{\prime\prime}33$. The deconvolved FWHM of the
emission was measured at $0.^{\prime\prime}81\times 0.^{\prime\prime}47$,
consistent with H18. Therefore, the ionized gas is distributed smoothly
without clumps, even in the improved resolution.

Figure 1: Top left: ALMA [O iii] 88-$\mu$m emission moment0 map contours on
the HST/WFC3 F160W image of MACS1149-JD1 at $z=9.1$. The contours show
$+3,4,5$, and $6\sigma$ with $\sigma=9.4$ mJy km s-1 beam-1. The synthesized
beam ellipse is shown in the bottom left corner. The circle and inverse-
triangle indicate the centers of [O iii] emission dynamical disk (§4) and
ultraviolet (UV) emission, respectively. Top right: [O iii] line spectrum
integrated over the area where the line was detected at $>3\sigma$. The gray
shaded region indicates the $\pm 1\sigma$ noise level. The solid curve is the
best-fit Gaussian profile. Bottom left: [O iii] line velocity field overlaid
on the line moment0 contours. The velocity field is depicted only in the area
where the Gaussian line profile fitting is favored with confidence $>5\sigma$
(see §3). Bottom right: [O iii] line velocity dispersion map overlaid on the
line moment0 contours. The depicted area is the same as the velocity field.

###  3.2 Velocity field

To examine the velocity structure of MACS1149-JD1, we adopted a method of Smit
et al. (2018): a Gaussian line profile fit for the spectrum at each spatial
pixel in the dirty data cube. We selected spatial pixels, where a $>5\sigma$
significance of the Gaussian fit was obtained. Namely,
$\Delta\chi^{2}=\chi^{2}_{\rm no~{}line}-\chi^{2}_{\rm Gauss}\geq 25$, where
$\chi^{2}_{\rm no~{}line}=\sum_{i}\frac{|I_{\rm obs}(V_{i})-I_{\rm
cont}|^{2}}{{\delta_{I_{\rm obs}(V_{i})}}^{2}}$ and $\chi^{2}_{\rm
Gauss}=\sum_{i}\frac{|I_{\rm obs}(V_{i})-\\{I_{\rm cont}+I_{\rm
Gauss}(V_{i})\\}|^{2}}{{\delta_{I_{\rm obs}(V_{i})}}^{2}}$. $I_{\rm
obs}(V_{i})$ and $I_{\rm Gauss}(V_{i})$ are, respectively, the observed
intensity and Gaussian line profile at velocity $V_{i}$. $I_{\rm cont}$ is the
continuum level, and we set it to zero because of its null-detection.
$\delta_{I_{\rm obs}(V_{i})}$ is the observed uncertainty of the intensity at
velocity $V_{i}$, which was measured as the root-mean-square (RMS) in the
velocity channel of the dirty cube. The resultant maps of the line velocity
and the velocity dispersion are shown in the bottom left and right panels in
Figure 1, respectively. There is a clear velocity gradient along the
north–south direction. The maximum velocity difference was measured at $\Delta
V_{\rm obs}=122\pm 30$ km s-1. The uncertainty was calculated from those of
the reddest and bluest velocities. The velocity dispersion ranges from a few
tens to a hundred km s-1, and its average is $\sim 70$ km s-1.

We have obtained a kinematic ratio of $\Delta V_{\rm obs}/2\sigma_{\rm
tot}=0.84\pm 0.23$, which is a factor of 1.5 greater than those observed in
two rotation-dominated galaxies at $z\sim 7$ reported by Smit et al. (2018),
as shown in the top panel of Figure 2. The ratio is also $\sim 2\sigma$ above
a criterion for determining whether a galaxy is rotation- or dispersion-
dominated, $(\Delta V_{\rm obs}/2\sigma_{\rm tot})_{\rm crit}=0.4$,
empirically derived from a set of simulations and H$\alpha$ line observations
of galaxies (e.g., Förster Schreiber et al. 2009; Smit et al. 2018).

Figure 2:  (Top) Observed kinematic ratio, $\Delta V_{\rm obs}/2\sigma_{\rm
tot}$, as a function of the stellar mass. Our measurement of MACS1149-JD1
shown by the five-pointed star is compared with those of two $z\sim 6$
galaxies (Smit et al., 2018) shown by the circles and the distribution of
measurements of $z\sim 2$ star-forming galaxies (Förster Schreiber et al.,
2009) shown by the shaded area. The horizontal line at $\Delta V_{\rm
obs}/\sigma_{\rm tot}=0.4$ is an empirical boundary between the rotation-
dominated and dispersion-dominated systems (Förster Schreiber et al., 2009).
(Bottom) Kinematic ratio, $V_{\rm rot}/\sigma_{V}$, based on dynamical
modeling (§4) as a function of redshift. The estimation for MACS1149-JD1 is
shown by the five-pointed star. The box plots at $z<4$ are taken from
Wisnioski et al. (2015) and show the median (middle bar) as well as the
central 50-percentile (box) and 90-percentile (vertical lines) of the sample
distribution in each survey (see Wisnioski et al. 2015 for the references).
The box plot at $z\sim 6$ is the result of quasar host galaxies (Neeleman et
al., 2021). The circles at $z\sim 4.5$ are results of massive rotational
galaxies (Lelli et al., 2021; Rizzo et al., 2021). The blue and magenta solid
lines show semi-empirical models based on Toomre’s disk instability parameter,
$Q_{\rm crit}=1.0$ (quasistable thin gas disk) (Wisnioski et al., 2015) with a
stellar mass of $10^{10.5}$ M⊙ and $10^{9.4}$ M⊙, respectively. The shaded
areas around the lines indicate the range of $Q_{\rm crit}=0.67$ (thick gas
disk) and $Q_{\rm crit}=2.0$ (star and gas composite disk) cases (Wisnioski et
al., 2015).

##  4 Dynamical modeling

###  4.1 Procedure

Motivated by the large kinematic ratio indicative of a rotation-dominated
system (§3), we performed dynamical modeling of MACS1149-JD1, assuming a
geometrically thin rotating disk. The model constraint was obtained by fitting
the 3D dirty data cube, rather than the two-dimensional (2D) velocity field,
to avoid rotational velocity underestimation and velocity dispersion
overestimation due to the beam smearing effect (Di Teodoro & Fraternali,
2015). The fitting procedure comprises four steps: (1) construction of a model
of the 3D emission line data cube in the source plane, (2) coordinate mapping
from the source plane to the image plane using a gravitational lensing model,
(3) convolution with the dirty beam profile in the image plane, and (4)
optimization. Before explaining each step below, we summarize the 9 fitting
parameters in the modeling: $M_{\rm disk}$ (mass), $r_{\rm disk}$ (scale-
length), $i_{\rm disk}$ (inclination), $PA_{\rm disk}$ (position angle),
$x_{0},~{}y_{0}$ (central position) $A_{0}$ ([O iii] intensity at the disk
center), $\sigma_{V}$ (velocity dispersion), and $\Delta V_{\rm sys}$
(velocity offset from the systemic redshift).

The construction of the line data cube was made by adopting Freeman’s formula
for the velocity field along the radial coordinate (Freeman, 1970). Namely, we
assumed a geometrically thin, self-gravitating rotation disk with an
exponential surface mass distribution along the radial distance, which is
described by $M_{\rm disk}$, $r_{\rm disk}$, $i_{\rm disk}$, and $PA_{\rm
disk}$. The maximum rotation velocity of the disk, $V_{\rm rot}$, is
calculated from $M_{\rm disk}$ and $r_{\rm disk}$ by $V_{\rm
rot}=0.88\sqrt{GM_{\rm disk}/2r_{\rm disk}}$ (Freeman, 1970), where $G$ is the
gravitational constant. We allowed a spatial offset of the disk center,
$(x_{0},y_{0})$, from a reference point in the source plane. The [O iii] line
velocity field was given by this disk model, and the line profile was assumed
to be a Gaussian function with a constant velocity dispersion, $\sigma_{V}$,
throughout the system. We allowed a constant velocity offset, $\Delta V_{\rm
sys}$, from the line redshift of $z=9.1096$ (H18). The velocity-integrated [O
iii] line intensity distribution is assumed to follow the surface mass
distribution of the dynamical disk with the last parameter describing the
central line intensity, $A_{0}$.

Coordinate mapping was performed by adopting lensing models (convergence,
$\kappa$, and shear, $\gamma_{1}$ and $\gamma_{2}$) of the MACS1149 cluster
released by the Hubble Frontier Field (HFF) project (Lotz et al.,
2017).222https://archive.stsci.e.,du/pub/hlsp/frontier/macs1149/models/ We
confirmed that all six lensing models in the HFF gave qualitatively same
results. In this Letter, we present the results with the model of Kawamata et
al. (2015) based on glafic (Oguri, 2010).

Convolution with the dirty beam in the image plane was performed to the
constructed line data cube model. We adopted a synthesized dirty beam model at
a velocity of 0 km s-1 in the observed dirty data cube.

Optimization was performed using a least-square method based on the
Levenberg–Marquardt algorithm in a package of
scipy.optimize.least_squares.333https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html
The model fitting was done in the image plane by comparing the modeled dirty
cube with the observed cube in the velocity range of $-150$ to $+200$ km s-1
with a 50-km s-1 binning. The spatial area used in the fitting was the region
where the [O iii] line was detected at $>3\sigma$ (i.e. the region enclosed by
the outer-most contour in the top left panel of Figure 1). The chi-square was
defined by $\sum_{i,j,k}|I^{\rm obs}_{i,j}(V_{k})-I^{\rm
model}_{i,j}(V_{k})|^{2}/{\delta_{I^{\rm obs}_{i,j}(V_{k})}}^{2}$, where
$I^{\rm model}_{i,j}(V_{k})$, $I^{\rm obs}_{i,j}(V_{k})$, and $\delta_{I^{\rm
obs}_{i,j}(V_{k})}$ are, respectively, the model intensity, the observed
intensity, and the observed RMS at the spatial pixel $i,j$ and the velocity
$V_{k}$.

The initial parameter set is crucial to obtain a converged solution. We
adopted an iterative method to ensure convergence, repeating the 3D fitting
where the best-fit parameters in the previous cycle were injected as the
initial parameter set. The initial guess for the first cycle was obtained from
a preparatory 2D fitting for the line moment0 map and velocity field. The
best-fit parameters and their uncertainties converged within eight cycles.

Uncertainties of the fitting parameters were obtained by a Monte Carlo method.
We repeated 3D fittings for mock observed data cubes with randomly-produced
image-plane noise maps. The noise was spatially correlated on the appropriate
beam-scale and its RMS was scaled to the observed one. We adopted the best-fit
values obtained from the real observed data cube as the most likely solution
and the values in the central 68-percentile in the distribution obtained by
the Monte Carlo method as their uncertainties.

In addition, we separately performed an exponential model fitting to the UV
continuum image of the Hubble Space Telescope (HST) Wide Field Camera 3
(WFC3)/F160W filter. The fitting parameters are the central brightness,
$A_{\rm UV}$, scale-length, $r_{\rm UV}$, inclination, $i_{\rm UV}$, and
position angle, $PA_{\rm UV}$. We also allowed a spatial offset of the UV
central position $(x_{\rm UV,0},y_{\rm UV,0})$ from the reference point in the
source plane. A Gaussian point-spread-function (PSF) with an FWHM of
$0.^{\prime\prime}15$ of the HST/WFC3 IR channel (Skelton et al., 2014) was
adopted for the convolution in the fitting procedure in the image plane. In
this particular case, we did not perform a Monte Carlo estimation of the
uncertainties.

###  4.2 Result

Figure 3 shows a comparison between the observed and best-fit [O iii] line
cubes in the image plane. Figure 4 shows the source-plane reconstruction of
the 2D [O iii] line moment0 map and velocity field from the best-fit line cube
as well as the UV intensity map. Low residuals shown in these figures
demonstrate a reasonably good fit. Table 1 presents a summary of the observed
and derived properties of MACS1149-JD1. The inclinations, PAs, and central
positions of the [O iii] disk and the UV disk coincide with each other, while
the scale-lengths are different.

We have found $V_{\rm rot}/\sigma_{V}=0.67_{-0.26}^{+0.73}$, whose range still
permits a value greater than 1 consistent with a rotation-dominated system.
This ratio is compared with similarly estimated values at lower redshift from
the literature in the bottom panel of Figure 2. $V_{\rm rot}/\sigma_{V}$ of
star-forming galaxies at $z<4$ and $z\sim 6$ quasar host galaxies are well-
explained by a semi-empirical model with a range of Toomre’s disk instability
parameter of $0.67<Q_{\rm crit}<2$ (Wisnioski et al., 2015). An extrapolation
of the model indicates that MACS1149-JD1 at $z=9.1$ has a quasi-stable disk
composed of stars and gas: $Q_{\rm crit}\sim 2$ (Wisnioski et al., 2015).

Figure 3: Image plane channelmap from $-150$ km s-1 (left) to $+200$ km s-1
(right) with a 50-km s-1 step. Top: the observed cube; middle: the model cube;
bottom: the residual.

Figure 4: Top row: the [O iii] 88 $\mu$m line moment0 map, middle row: the [O
iii] line velocity field, and bottom row: the UV continuum image. Left: the
observations, middle: the models, and right: the residuals. All images are
shown in the source plane corrected for the gravitational lensing effect.
Table 1: A summary of properties of MACS1149-JD1.

Property | Value | Remark
---|---|---
Observed values
R.A. (ICRS) | $11^{h}49^{m}33.58s$ | H18
Dec. (ICRS) | $+22^{\circ}24^{\prime}45.7^{\prime\prime}$ | H18
$z_{\rm spec}$ | $9.1096\pm 0.0006$ | H18
$M_{*}$ [$10^{9}$ M⊙] | $1.08^{+0.53}_{-0.18}\times(10/\mu)$ | H18
$SFR$ [M⊙ yr-1] | $4.2^{+0.8}_{-1.1}\times(10/\mu)$ | H18
$\Delta V_{\rm obs}/2\sigma_{\rm tot}$ | $0.84\pm 0.23$ | This work
Assumed value: Lensing magnification
$\mu$ | 4.7 |  glafic model
Best-fit values and their 68-percentile of the dynamical model.
$\Delta V_{\rm sys}$ [km s-1] |  $10.3$ ($-0.5$–$20.2$) |
$V_{\rm rot}$ [km s-1] |  $46.7$ ($28.3$–$87.1$) | ($\dagger$)
$\sigma_{V}$ [km s-1] |  $69.9$ ($55.6$–$75.9$) |
$V_{\rm rot}/\sigma_{V}$ |  $0.67$ ($0.41$–$1.40$) | ($\ddagger$)
$M_{\rm disk}$ [$10^{9}$ M⊙] |  $0.65$ ($0.25$–$2.02$) | (s)
$r_{\rm disk}$ [kpc] |  $0.50$ ($0.40$–$0.61$) | (s)
$i_{\rm disk}$ [∘] |  $56.8$ ($44.9$–$66.3$) |
$PA$(disk) [∘] |  $157.5$ ($143.2$–$168.4$) |
Reduced $\chi^{2}$ (Data cube) | 1.18 |
Best-fit values of the UV emission model.
$r_{\rm UV}$ [kpc] | $0.31$ | (s)
$i_{\rm UV}$ [∘] | $52$ |
$PA$(UV) [∘] | $138$ |
$\Delta$(UV–disk) [kpc] | $0.13\pm 0.2$ | (s) (*)
Reduced $\chi^{2}$ (UV) | 1.08 |

Note. — H18: Hashimoto et al. (2018). glafic model: Kawamata et al. (2015);
Oguri (2010). ($\dagger$): Derived from $M_{\rm disk}$ and $r_{\rm disk}$, not
a fitting parameter. ($\ddagger$): Derived from $V_{\rm rot}$ and
$\sigma_{V}$, not a fitting parameter. (s): Values in the source plane. (*)
Systematic uncertainty in astrometory between the HST UV image and ALMA data
cube is $0.^{\prime\prime}1$ (H18), corresponding to $\sim 0.2$ kpc in the
source plane.

##  5 Discussion

MACS1149-JD1 satisfies an empirical criterion for a rotation-dominated system,
$\Delta V_{\rm obs}/2\sigma_{\rm tot}>0.4$ (Förster Schreiber et al., 2009).
Among the five criteria for a rotation disk (Wisnioski et al., 2015), the
galaxy also satisfies (i) a continuous velocity gradient along a single axis,
(ii) $V_{\rm rot}/\sigma_{V}>1$ (but marginally), (iv) an agreement between
photometric (UV) and kinematic ([O iii]) axes ($<30^{\circ}$), and (v) a
positional agreement of the dynamical center ([O iii]) and continuum centroid
(UV). However, it does not satisfy criterion (iii) a positional agreement
between the dynamical center and the velocity dispersion peak.

For the important criterion (ii), our dynamical model gives $V_{\rm
rot}/\sigma_{V}=0.67$ ($0.41$–$1.40$, 68% range), still indicative of a
rotation, taking into account the possible range. We also examined a case
where the [O iii] emission and the dynamical disk are two different components
in the light of the two-component stellar populations discussed in H18.
Namely, we considered 5 additional parameters of the central positions
($x,y$), inclination, PA, and scale-length for the [O iii] emission disk
“decoupled” from the dynamical disk introduced in §4.1. As a result, we
obtained a higher kinematic ratio of $V_{\rm rot}/\sigma_{V}=1.4$. The
dynamical disk center was spatially offset from the centers of [O iii]
emission and UV continuum arisen by the star-forming population, as found in
cosmological simulations (e.g., Moriwaki et al. 2018). Although the fitting
uncertainty was large due to the larger number of free-parameters, this
decoupled disk scenario would be very interesting to examine with James Webb
Space Telescope (JWST).

The criterion (iii) should also be discussed. The velocity dispersion peak is
located $\sim 1$ kpc away from the [O iii] disk center, possibly indicating a
merger or an outflow. A possible small kink in UV emission around the
dispersion peak (Fig. 1) might be a sign of a merger. The weaker Ly$\alpha$
line showing blueshift from the [O iii] line (H18) may be another sign of a
different component in the galaxy. However, there is no distinct structure in
[O iii] emission. We could not examine this possibility further with the
current dataset. The higher spatial resolution and sensitivity offered by JWST
is required to address the merger possibility. In the following, we consider
that MACS1149-JD1 is a rotating disk galaxy.

The stellar mass of MACS1149-JD1 corrected by the glafic model ($\mu=4.7$) is
$2.3^{+1.1}_{-0.4}\times 10^{9}$ M⊙. The bulk of the mass is associated with
the $\sim 300$ Myr-old mature stellar population, whereas $\sim 3$ Myr-old
star-forming population is a minor component (H18). Crucially, the dynamical
disk mass of $0.65_{-0.40}^{+1.37}\times 10^{9}$ M⊙ is consistent with that
independently-determined for the mature stellar population. If we consider the
effect of the velocity dispersion in the case of $V_{\rm rot}/\sigma_{V}\sim
1$ (e.g., Burkert et al. 2010), the dynamical mass can be $\sim 5$ times
larger,444From equation (11) in Burkert et al. (2016), one can derive the
dynamical mass corrected for the turbulent pressure as $M_{\rm dyn}=M_{\rm
disk}\times\\{1+4.4(\sigma_{V}/V_{\rm rot})^{2}\\}$ for $M_{\rm disk}$ and
$V_{\rm rot}$ evaluated at $r=2.2r_{\rm disk}$ where the Freeman disk has the
maximum rotation velocity (Freeman, 1970). yielding $\sim 3\times 10^{9}$ M⊙,
that is also consistent with the stellar mass of $\sim 2\times 10^{9}$ M⊙
(H18). Therefore, we conclude that the dynamical mass is attributable to the
mature stellar population that formed at $z\sim 15$.

A cold gas component should exist in the galaxy because of the ongoing star
formation. The gas mass fraction can be expressed as $f_{\rm gas}=(a/Q_{\rm
crit})(\sigma_{V}/V_{\rm rot})$ with $a=1$–2, depending on the velocity radial
profile (Genzel et al., 2011; Wisnioski et al., 2015). Because $Q_{\rm
crit}=1$–2 for a stellar-plus-gas disk (e.g., Kim & Ostriker 2007; Genzel et
al. 2011), the obtained $V_{\rm rot}/\sigma_{V}$ suggests $f_{\rm gas}\gtrsim
0.3$. Considering uncertainties in the estimated masses, it is possible that
$\sim 1\times 10^{9}$ M⊙ cold gas component and $\sim 2\times 10^{9}$ M⊙
mature stellar population coexist, yielding $f_{\rm gas}\sim 0.3$ in a total
mass of $\sim 3\times 10^{9}$ M⊙. The star formation rate of the galaxy is
$8.9^{+1.7}_{-2.3}$ M⊙ yr-1 for the magnification $\mu=4.7$ (H18). Hence, we
obtain a gas depletion time of $t_{\rm dep}\sim 100$ Myr when the gas mass of
$1\times 10^{9}$ M⊙, which agrees with an extrapolation of a mean relation of
$t_{\rm dep}\sim 1.5/(1+z)$ Gyr (Wisnioski et al., 2015) to $z=9$.

The scale-length in the source plane of [O iii] disk is 1.6-times larger than
that of UV continuum. Notably, both scale-lengths should be readily resolved
by the ALMA beam and HST PSF sizes in the source plane ($\sim 0.3$ and $\sim
0.1$ kpc in radius, respectively). The [O iii] emitting ionized gas should be
powered by the young stellar population traced by UV. The extended
distribution of the ionized gas compared to the young star cluster possibly
suggests a significant escape of ionizing photons to a larger-scale, and the
galaxy may contribute to cosmic reionization. This may also explain the
blueshift of the Ly$\alpha$ line (H18).

We have considered only a single mass component in the dynamical disk
modeling. Because we are considering the very central part ($r\lesssim 1$ kpc)
of the galaxy, the dark matter contribution is generally negligible (e.g., van
Albada et al. 1985). According to Behroozi et al. (2019), the dark matter halo
mass for a stellar mass of $\sim 10^{9}$ M⊙ galaxy at $z=9$ is $\sim 10^{11}$
M⊙. The virial radius of such a halo would be $\sim 7$ kpc in the proper
coordinate (Dayal & Ferrara, 2018). Therefore, our observations do not reach
the halo scale yet; hence, neglecting its contribution is reasonable.

In conclusion, MACS1149-JD1 at $z=9.1$, is the most distant galaxy with a
signature of rotation. This is not contradictory to the concordance
cosmological structure formation. Some theoretical studies predicted such a
rotational disk in the earliest universe (e.g., Robertson et al. 2006; Katz et
al. 2019). It is interesting to understand the role of two stellar populations
proposed by H18 in the rotation disk formation. We suggest that the mass of
the rotational disk is dominated by the mature stellar population formed in
the first major star formation episode in this galaxy at $z\sim 15$ (H18).
JWST’s Guaranteed Time Observation programs targeting this galaxy will resolve
the different spatial distributions of the young and mature stellar
populations and confirm (or revise) the scenario.

We thank Renske Smit for a discussion about the method for examining the
velocity structure, Masamune Oguri for a discussion about gravitational
lensing models, Yuxing Zhong for a discussion about 3D Baroro. AKI, YS, and YF
are supported by NAOJ ALMA Scientific Research Grant Numbers 2020-16B. TH was
supported by Leading Initiative for Excellent Young Researchers, MEXT, Japan
(HJH02007) and by JSPS KAKENHI Grant Number (20K22358 and 22H01258) RSE
acknowledges funding from the European Research Council under the European
Union Horizon 2020 research and innovation programme (grant agreement No.
669253). NL acknowledges the Kavli foundation. This paper makes use of the
following ALMA data: ADS/JAO.ALMA#2015.1.00428.S and
ADS/JAO.ALMA#2018.1.00616.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ.

## References

  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Binggeli et al. (2019) Binggeli, C., Zackrisson, E., Ma, X., et al. 2019, MNRAS, 489, 3827, doi: 10.1093/mnras/stz2387
  * Bouwens et al. (2014) Bouwens, R. J., Bradley, L., Zitrin, A., et al. 2014, ApJ, 795, 126, doi: 10.1088/0004-637X/795/2/126
  * Burkert et al. (2010) Burkert, A., Genzel, R., Bouché, N., et al. 2010, ApJ, 725, 2324, doi: 10.1088/0004-637X/725/2/2324
  * Burkert et al. (2016) Burkert, A., Förster Schreiber, N. M., Genzel, R., et al. 2016, ApJ, 826, 214, doi: 10.3847/0004-637X/826/2/214
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021, doi: 10.1093/mnras/stv1213
  * Förster Schreiber et al. (2009) Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364, doi: 10.1088/0004-637X/706/2/1364
  * Freeman (1970) Freeman, K. C. 1970, ApJ, 160, 811, doi: 10.1086/150474
  * Genzel et al. (2011) Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101, doi: 10.1088/0004-637X/733/2/101
  * Hashimoto et al. (2018) Hashimoto, T., Laporte, N., Mawatari, K., et al. 2018, Nature, 557, 392, doi: 10.1038/s41586-018-0117-z
  * Hashimoto et al. (2019) Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71, doi: 10.1093/pasj/psz049
  * Hoag et al. (2018) Hoag, A., Bradač, M., Brammer, G., et al. 2018, ApJ, 854, 39, doi: 10.3847/1538-4357/aaa9c2
  * Huang et al. (2016) Huang, K.-H., Bradač, M., Lemaux, B. C., et al. 2016, ApJ, 817, 11, doi: 10.3847/0004-637X/817/1/11
  * Inoue et al. (2016) Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559, doi: 10.1126/science.aaf0714
  * Jones et al. (2010) Jones, T. A., Swinbank, A. M., Ellis, R. S., Richard, J., & Stark, D. P. 2010, MNRAS, 404, 1247, doi: 10.1111/j.1365-2966.2010.16378.x
  * Katz et al. (2019) Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902, doi: 10.1093/mnras/stz1672
  * Kawamata et al. (2015) Kawamata, R., Ishigaki, M., Shimasaku, K., Oguri, M., & Ouchi, M. 2015, ApJ, 804, 103, doi: 10.1088/0004-637X/804/2/103
  * Kim & Ostriker (2007) Kim, W.-T., & Ostriker, E. C. 2007, ApJ, 660, 1232, doi: 10.1086/513176
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18
  * Laporte et al. (2021) Laporte, N., Meyer, R. A., Ellis, R. S., et al. 2021, MNRAS, 505, 3336, doi: 10.1093/mnras/stab1239
  * Lelli et al. (2021) Lelli, F., Di Teodoro, E. M., Fraternali, F., et al. 2021, Science, 371, 713, doi: 10.1126/science.abc1893
  * Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97, doi: 10.3847/1538-4357/837/1/97
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Moriwaki et al. (2018) Moriwaki, K., Yoshida, N., Shimizu, I., et al. 2018, MNRAS, 481, L84, doi: 10.1093/mnrasl/sly167
  * Neeleman et al. (2021) Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141, doi: 10.3847/1538-4357/abe70f
  * Oguri (2010) Oguri, M. 2010, PASJ, 62, 1017, doi: 10.1093/pasj/62.4.1017
  * Rizzo et al. (2021) Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952, doi: 10.1093/mnras/stab2295
  * Rizzo et al. (2020) Rizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201, doi: 10.1038/s41586-020-2572-6
  * Roberts-Borsani et al. (2020) Roberts-Borsani, G. W., Ellis, R. S., & Laporte, N. 2020, MNRAS, 497, 3440, doi: 10.1093/mnras/staa2085
  * Robertson et al. (2006) Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, ApJ, 645, 986, doi: 10.1086/504412
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24, doi: 10.1088/0067-0049/214/2/24
  * Smit et al. (2018) Smit, R., Bouwens, R. J., Carniani, S., et al. 2018, Nature, 553, 178, doi: 10.1038/nature24631
  * van Albada et al. (1985) van Albada, T. S., Bahcall, J. N., Begeman, K., & Sancisi, R. 1985, ApJ, 295, 305, doi: 10.1086/163375
  * Watson et al. (2015) Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327, doi: 10.1038/nature14164
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209, doi: 10.1088/0004-637X/799/2/209
  * Zheng et al. (2012) Zheng, W., Postman, M., Zitrin, A., et al. 2012, Nature, 489, 406, doi: 10.1038/nature11446
  * Zheng et al. (2017) Zheng, W., Zitrin, A., Infante, L., et al. 2017, ApJ, 836, 210, doi: 10.3847/1538-4357/aa5d55