Detecting a disc bending wave in a barred-spiral galaxy at redshift 4.4

ABSTRACT

The recent discovery of barred spiral galaxies in the early Universe (z > 2) poses questions of how these structures form and how they influence galaxy evolution in the early Universe. In this study, we investigate the morphology and kinematics of the far-infrared (FIR) continuum and [C ii] emission in BRI1335-0417 at z ≈ 4.4 from ALMA observations. The variations in position angle and ellipticity of the isophotes show the characteristic signature of a barred galaxy. The bar, |$3.3^{+0.2}_{-0.2}$| kpc long in radius and bridging the previously identified two-armed spiral, is evident in both [C ii] and FIR images, driving the galaxy’s rapid evolution by channelling gas towards the nucleus. Fourier analysis of the [C ii] velocity field reveals an unambiguous kinematic m = 2 mode with a line-of-sight velocity amplitude of up to ∼30–40 km s⁻¹; a plausible explanation is the disc’s vertical bending mode triggered by external perturbation, which presumably induced the high star formation rate and the bar/spiral structure. The bar identified in [C ii] and FIR images of the gas-rich disc galaxy (≳ 70 per cent of the total mass within radius R ≈ 2.2 disc scale lengths) suggests a new perspective of early bar formation in high redshift gas-rich galaxies – a gravitationally unstable gas-rich disc creating a star-forming gaseous bar, rather than a stellar bar emerging from a pre-existing stellar disc. This may explain the prevalent bar-like structures seen in FIR images of high-redshift submillimeter galaxies.

galaxies: bar, galaxies: disc, galaxies: kinematics and dynamics, galaxies: spiral

1 INTRODUCTION

Bar structure plays a crucial role in driving galaxy evolution and shaping disc structure. In galaxies an axisymmetric stellar bar exerts gravitational torque on the gas, driving it towards the galactic centre and forming a centralized stellar structure such as a bulge and nuclear disc (Athanassoula 1992; Wada & Habe 1992). This process may also promote gas accretion onto the black hole observed as active galactic nuclei (AGN; Hopkins & Quataert 2010; Emsellem et al. 2015). Bars can also drive radial migration of gas and stars, which is essential for explaining the observed stellar kinematics in Milky Way galaxies (e.g. Kawata et al. 2021).

Numerical simulation suggests that stellar bar formation in galaxies leads to an immediate gas inflow into the central region and the formation of a nuclear disc, making the stellar age of the nuclear disc a good indicator of the epoch when the bar first formed (Baba & Kawata 2020). Rather earlier however, the idea was applied to observations for estimating bar formation age in Gadotti et al. (2015).

Recent observations using integral field spectroscopy (IFS) have provided insight into the stellar populations of nuclear discs in barred galaxies (Gadotti et al. 2015; Bittner et al. 2020; de Sá-Freitas et al. 2023), suggesting that the oldest nuclear discs in barred galaxies are at least 10 Gyr old (z > 2; de Sá-Freitas, Gadotti & Fragkoudi 2023). The findings, pointing to early bar formations, align with the latest high redshift observations. For instance, recent Atacama Large Millimeter/submillimeter Array (ALMA) observations show gas discs already formed at z > 4 (e.g. Neeleman et al. 2020; Rizzo et al. 2020; Lelli et al. 2021; Tsukui & Iguchi 2021), while James Webb Space Telescope (JWST) observations recently discovered early barred galaxies at redshift 1 < z < 3 (Guo et al. 2023; Le Conte et al. 2023) and at z = 4.2 (Smail et al. 2023). In addition, numerous disc-like systems (⁠|$\sim 50~{{\ \rm per \, cent}}$|⁠) have been found at 3 < z < 6 by JWST (Ferreira et al. 2022; Nelson et al. 2023). While these recent discoveries are consistent with some pre-JWST and ALMA studies showing disc prevalence already at z ∼ 2.6 (Wuyts et al. 2011; Wisnioski et al. 2015, but see Conselice 2014), they show surprisingly numerous discs up to z ∼ 6 and earlier spiral and bars than expected (Elmegreen & Elmegreen 2014).

Numerical simulations motivated by these discoveries show that baryon-dominated disc galaxies promptly form a bar (Fujii et al. 2018; Bland-Hawthorn et al. 2023). It remains an open question if this rule applies to galaxies with high gas fractions (more than ∼ 50 per cent of total baryonic mass at z > 3; Carilli & Walter 2013), as some simulations suggest molecular gas can suppress bar formation or result in a weaker bar (Athanassoula, Machado & Rodionov 2013; Łokas 2020). Conversely, ALMA observations reveal prevalent bar morphologies in dust continuum images of gas-rich submillimeter galaxies (Gullberg et al. 2019; Hodge et al. 2019; Smail et al. 2023). Tidal interaction is another promising avenue to form bars (e.g. Noguchi 1996; Łokas et al. 2014, 2016), even for gas-rich systems (Gajda, Łokas & Athanassoula 2018). Recent cosmological simulations, including a realistic high gas fraction of high-redshift galaxies and external effects such as mergers, suggest a high bar fraction out to z ∼ 4 (Rosas-Guevara et al. 2022). Different sub-grid physics implementations in simulations can result in different galaxy properties and bar dynamics (Fragkoudi et al. 2021).

Recently, Tsukui & Iguchi 2021 revealed a spiral morphology in BRI 1335–0417 at z = 4.4074 (Guilloteau et al. 1997), hosting an optical quasar initially identified in Automatic Plate Measuring (APM) survey (Storrie-Lombardi et al. 1996). The galaxy exhibits a high star formation rate (SFR) ∼1700 M_⊙ yr⁻¹ estimated from the spectral energy distribution modelling with the AGN contribution being corrected using spatially resolved information (point source and extended source separation; Tsukui et al. 2023b), making it one of the brightest unlensed submillimeter source at z > 4 (Jones et al. 2016). The [C ii] and dust continuum maps of the galaxy revealed a two-armed structure with a pitch angle of |${26.7^{+4.1}_{-1.6}}^\circ$| (Tsukui & Iguchi 2021). These arms extend from 2 to 5 kpc and appear to start at the end of an elongated bar-like structure that bridges them.

Despite new observational and theoretical results, it remains unclear what was the dominant cause for early bar and spiral formation – internal or external processes. As the brightest and earliest barred spiral example, BRI 1335–0417 allows us to study the detailed [C ii] line and FIR continuum morphology in unprecedented detail and without uncertainties due to lens model reconstruction. The spatially resolved [C ii] line kinematics together with numerical simulations provide an excellent laboratory to provide new insights into early bar formation.

Gas disc kinematics tell us not only the dynamics of the galaxy, such as disc stability and underlying mass distribution of the galaxy. Subtracting the overall rotation and examining the more subtle residual velocity field allows us to explore further: gas inflows (e.g. Di Teodoro & Peek 2021; Genzel et al. 2023), bar and dynamical effect of the bar/spiral structure (e.g. Grand et al. 2016; Monari et al. 2016; Gómez et al. 2021), and even bending waves of the disc – seismic ripple propagating through the disc due to perturbation by a recent interaction with satellites (Bland-Hawthorn & Tepper-García 2021; Tepper-García, Bland-Hawthorn & Freeman 2022; Urrejola-Mora et al. 2022) or misaligned gas accretion (Khachaturyants et al. 2022).

In this paper, we report a new analysis of the bar structure in the quasar host galaxy BRI 1335–0417 using ALMA Band 7 data of the far-infrared (FIR) continuum (observed frame |$\sim 869\, \mu \mathrm{m}$| or rest-frame |$\sim {160}\, \mu \mathrm{m}$|⁠) and ionized carbon [C ii]. We also demonstrate that the rotation-subtracted residual velocity field is consistent with the dynamical imprint of a recent interaction by an external perturber, which likely induced the bar and spiral density wave in the gas disc.

Throughout this paper, we assume a flat Λ-dominated cold dark matter cosmology with a present-day Hubble constant of H₀ = 70 km s⁻¹ Mpc⁻¹, and a density parameter of pressureless matter Ω_M = 0.3.

2 OBSERVATION AND DATA REDUCTION

This study uses rest-frame |$160\, \mu \mathrm{m}$| FIR continuum image and [C ii] line cube from the observation program #2017.1.00394.S (PI = González López, Jorge), which was carried out on 2018 January 21. The data calibration and reduction details were described in (Tsukui & Iguchi 2021; Tsukui et al. 2023b). In order to accurately quantify the elongated bar structure, both [C ii] and FIR intensity maps were convolved to have a circular-shaped beam (point spread function; PSF) with the same resolution [full width half-maximum (FWHM) = 0.195 arcsec = 1.3 kpc]. The [C ii] line data are binned in the spectral axis to have a spectral resolution of 20 km s⁻¹.

3 RESULTS AND DISCUSSION

In this section, we present the results of structural investigations in the first two subsections and then investigations from the dynamical point of view using first-order and higher-order [C ii] kinematics in the subsequent subsections.

3.1 Bar identification by the ellipse fitting

We employ the commonly used ellipse-fitting method to examine the photometric structure of the disc (e.g. Jedrzejewski 1987; Wozniak et al. 1995; Erwin 2005; Gadotti et al. 2007). We fit ellipses to isophotes of the [C ii] and FIR intensity images, with position angle (PA) and ellipticity (ϵ) allowed to vary for each ellipse and the central position fixed to the best fit of the smallest ellipse. To estimate the statistical uncertainty in our measurements, we performed the fitting procedure 300 times, each time adding realistic correlated noise to the original images. The noise properties are measured using the noise auto-correlation function (Tsukui et al. 2023a).

The left column of Fig. 1 shows the ellipticity and PA of the best-fitting ellipses as a function of radius. The right column presents the best-fitting isophote ellipses overlaid on [C ii] and FIR images. The ellipticity profile of the best-fitting ellipses exhibits a characteristic profile common to barred galaxies (Wozniak et al. 1995). At smaller radii, the ellipticity is small as the ellipses trace the centrally concentrated light. Then the ellipticity reaches a maximum as the ellipse aligns with the elongated bar shape. Subsequently, the ellipticity decreases as the ellipse traces the disc more circular than the bar. The PA changes slowly as the ellipse traces the bar and rapidly changes at the end of the bar as the ellipse starts tracing the disc because the disc PA is offset from the bar in most cases.¹

Figure 1.

Left column: Ellipticity (top) and PA (bottom) of isophote ellipses as a function of the ellipse radius (semimajor axis) of [C ii] (blue) and FIR (orange) images of BRI 1335–0417. The vertical lines indicate the radii of the maximum ellipticity for [C ii] (blue) and FIR (orange). For each measurement, 1σ uncertainty is shown with a shaded region. The grey-shaded region indicates the PSF FWHM, where accurate ellipticity and PA estimates are not possible. Right column: [C ii] and FIR images of BRI 1335–0417 overlaid with the best-fitting isophote ellipses (black ellipses), the bar ellipse (a blue ellipse) and its PA in the solid blue line. The disc kinematic PA (Tsukui & Iguchi 2021) is also shown in the dashed blue line. The FWHM of PSF is shown in the bottom left corner of both images.

In both [C ii] and FIR continuum, the ellipticity reaches a maximum at the same radius within some uncertainty. Therefore, we adopt the radius of the maximum ellipticity of [C ii] as a fiducial sky-projected bar length, |$R^{\mathrm{sky}}_{\textrm {bar}, \epsilon _\textrm {max}} = 3.1^{+0.2}_{-0.1}$| kpc. This bar length is larger than the 2×FWHM of the beam, which is the required minimum criterion for detecting the bar (Erwin 2018). Another way to define the bar length from the ellipse fitting is the radius where the PA changes by 5° relative to its value at the radius of maximum ellipticity. This alternative definition provides slightly larger values, |$R^{\mathrm{sky}}_{\textrm {bar}, \Delta \textrm {PA} = 5\textrm {deg}} = 3.5^{+0.3}_{-0.1}$| kpc yet similar to our adopted fiducial value, confirming that the choice of the definition does not affect the conclusion.

We find the intrinsic bar length |$R^{\mathrm{int}}_{\textrm {bar}, \epsilon _\textrm {max}} = 3.3^{+0.2}_{-0.2}$| kpc by the analytical deprojection assuming a planer bar ellipse (Gadotti et al. 2007) and the disc inclination (37.3°) and PA (7.6°) estimated by Tsukui & Iguchi (2021). The bar length relative to the [C ii] disc scale length R_d = 1.83 ± 0.04 kpc (Tsukui & Iguchi 2021) is larger than stellar bars seen in nearby galaxies (Erwin 2019; their derived scaling relation predicts a bar length of 1.7 ± 0.1 kpc given the same disc scale length).

We note that the PA of the outermost isophote is misaligned with the kinematic disc PA. The outer ellipses are influenced not only by noise, as evidenced by the increased uncertainty shown in the bottom left panel of Fig. 1, but also by the spiral structure and faint tidal tail-like structure extending from north-east to south-west. Therefore, in later analysis, we will use the kinematic disc PA rather than the photometric isophote PA as the disc’s PA.

3.2 Interpretation of [C ii] and FIR continuum bars

The [C ii] line traces not only overall star formation but also multiphase gas distribution, from neutral gas (Herrera-Camus et al. 2018) to molecular and ionized gas (Pineda et al. 2013). On the other hand, the FIR continuum represents thermal emission from dust heated by ultraviolet emission produced by young massive stars. Although the FIR continuum is commonly used as the tracer for dust mass and SFR (young massive stars heating the dust), the relationships of the single band FIR continuum to these two quantities are complex depending on the spatial temperature and opacity variation (da Cunha et al. 2021; Tsukui et al. 2023b). The FIR continuum is proportional to the dust mass and SFR² specifically in a limiting case of dust being optically thin with a spatially constant temperature.

The ellipse fitting result confirms both [C ii] and FIR continuum images have elliptical elongated bar shapes, which are shown in Fig. 2 along with the two-armed spiral confirmed in Tsukui & Iguchi 2021 (projected on the sky). The configuration of the bar and spiral is similar to the nearby grand-design spiral galaxies observed (e.g. Stuber et al. 2023) and found in simulations (e.g. Baba 2015).

Figure 2.

[C ii] and FIR images of BRI 1335–0417 overlaid with our identified bar ellipses (blue ellipse) and two-armed spiral structure identified in Tsukui & Iguchi (2021, black solid line). Contours start at 2σ in both maps but are linearly spaced by 2σ for [C ii] map and logarithmically spaced in powers of 2 for FIR map (2σ, 4σ, 8σ,...).

In the nearby Universe, not all barred galaxies show star formation and gas reservoirs in the bar region. Star formation or gas lanes along the bar region are preferentially found in gas-rich late-type barred galaxies (Díaz-García et al. 2020; Fraser-McKelvie et al. 2020) and more massive gas-rich barred systems (Stuber et al. 2023). FIR continuum bars are often seen in high redshift (z ∼ 2–4) submillimeter galaxies (Gullberg et al. 2019; Hodge et al. 2019), implying that the star-forming gasious bar may be common in the gas-rich early galaxies (Carilli & Walter 2013).

Interestingly, there is a slight yet statistically significant offset in PA of the bar ellipses between [C ii] and FIR continuum (PA|$_{\rm [C\, {\small II}]} = -15.9^{+4.0}_{-3.1}$|⁠, PA|$_{\mathrm{FIR}} = -6.9^{+4.2}_{-4.7}$|⁠). Also, the bar ellipse of [C ii] is more elongated than that of the FIR continuum (⁠|$\epsilon _{\rm [C\, {\small II}]} = 0.27^{+0.02}_{-0.02}$|⁠, |$\epsilon _{\mathrm{FIR}} = 0.39^{+0.03}_{-0.02}$|⁠). Although the rounder shape of the FIR bar may be attributed to the significant AGN contribution to the FIR image suggested by Tsukui et al. (2023b) these morphological differences, if further confirmed by higher resolution data, may tell the dynamics of the bar.

These differences may be consistent with the bar characteristics commonly confirmed in local observations (Fraser-McKelvie et al. 2020) and simulation (Emsellem et al. 2015) – the compressed gas flow (traced by [C ii]) which usually forms a straight line into the centre of the galaxy along the leading edge of the bar made up by the stellar components (traced by FIR emission), making the [C ii] bar thinner than the FIR bar. The [C ii] bar is also inclined towards the leading edge of the FIR bar.

This interpretation requires two assumptions. (1) The spiral arm in this galaxy is trailing as the majority of spiral galaxies are trailing (e.g. Iye, Tadaki & Fukumoto 2019). If so, the right (west) side in Fig. 2 is the far side of the disc, and the disc rotates clockwise, given that the north/south is the receding/approaching side of the disc (see Fig. 3a). The bar ellipse of the [C ii] line is displaced towards the leading edge of the bar (on the downstream side of the gas motion). (2) The FIR emission may trace more towards young stellar components than interstellar medium (ISM) traced by [C ii], as the FIR emission can depend on the distribution of the dust (ISM) and young massive stars which heat the dust. Especially in the early universe, ∼ 1 Gyr, the main producers of dust grains are asymptotic giant branch (AGB) stars or core-collapse supernovae (Gall, Hjorth & Andersen 2011). AGB stars start contributing at 10–100 Myr after the stellar evolution of stars with the mass of 3–8 M_⊙ while core-collapse supernovae contribute after the lifetime of 10 Myr. Within the relatively longer dynamical time-scale of the disc, 120 Myr, the stellar distribution producing the dust would not have displaced relative to the dust distribution.

Figure 3.

(a) [C ii] velocity field derived fitting a Gaussian at each spatial pixel of the data cube. (b) Best-fitting kinemetry model using odd Fourier series (m = 1, 3, 5) to expand the azimuthal velocity profile. (c) Best-fitting kinemetry model using odd + even Fourier series (m = 1, 2, 3, 4, 5). (d) Residual map subtracting the model with odd terms only (b) from the data velocity field (a). (e) Residual map subtracting the model with both odd and even terms (C) from the data velocity field (a). In (a), (d), and (d), [C ii] contours, bar ellipse, and two armed-spiral are overplotted in the same way as in Fig. 2.

As discussed in the next section, the stellar content is quite small relative to gas, so the bar dynamics and formation may differ from the nearby stellar-dominated bar paradigm.

3.3 Dynamical constraint from [C ii] rotation curve

The detailed mass distribution within galaxies is an essential factor for bar formation (e.g. Efstathiou, Lake & Negroponte 1982; Fujii et al. 2018; Bland-Hawthorn et al. 2023; Romeo, Agertz & Renaud 2023). A disc-dominated system is prone to instability and bar formation, while spherical structures such as dark matter and bulges as well as the presence of gas can suppress the bar formation. In this section, we derive a lower boundary on the disc fraction and gas fraction using (1) the dynamical mass distribution estimated by Tsukui & Iguchi (2021) whose rotation curve is consistent with the independent modelling by Roman-Oliveira, Fraternali & Rizzo (2023) and (2) the CO(2-1) line luminosity which traces molecular gas mass (Jones et al. 2016). Then, we discuss the implication for the bar formation in this galaxy.

Tsukui & Iguchi (2021) decomposed the rotation curve into contributions from the bulge and disc, assuming a de Vaucouleurs mass distribution for the bulge and an exponential for the disc (with the scale radius from the [C ii] emission). The derived disc mass, |$M_\mathrm{disc} = 4.9^{+1.7}_{-2.5}\times 10^{10}\, {\rm M}_{\odot }$|⁠, is consistent with the lower limit of the molecular gas mass, |$M_\mathrm{gas} = 5.1\times 10^{10}\, {\rm M}_{\odot }$| estimated by Jones et al. (2016) assuming solar metalicity and a typical CO line ratio r₂₁ = L_CO(2-1)/L_CO(1-0) = 0.85 for submillimeter galaxies (Carilli & Walter 2013), suggesting that the galaxy is baryon-dominated and gas-rich.

As Tsukui et al. (2023b) revealed the optical emission of BRI 1335–0417 is dominated by a quasar, using a typical r₂₁ = 0.99 for quasars (Carilli & Walter 2013) further reduces the estimated molecular gas mass by 10 per cent to |$M_\mathrm{gas} = 4.6\times 10^{10}\, {\rm M}_{\odot }$|⁠. Using this lower limit on the molecular gas estimate, we derive a lower limit for the disc mass fraction within the radius at which disc dynamics dominates, 2.2R_d = 4.0 kpc,

$$\begin{eqnarray} f_{\mathrm{disc}} = \left(\frac{v_{\mathrm{disc}}(R)}{v_{\mathrm{total}}(R)}\right)^2_{R = 2.2R_{\mathrm{disc}}}\gt \left(\frac{v_{\mathrm{gas}}(R)}{v_{\mathrm{total}}(R)}\right)^2_{R = 2.2R_{\mathrm{disc}}}\gt 0.73\pm 0.07 \\ \end{eqnarray}$$

(1)

where v_total(2.2R_d) = 200 ± 10 km s⁻¹ is the total circular velocity at 2.2R_d (Tsukui & Iguchi 2021), and v_gas(2.2R_d) = 179 km s⁻¹ is the circular velocity of the exponential gas disc at 2.2R_d with the lower limit mass.³ This lower limit on the disc mass fraction can also be interpreted as a lower limit on the gas mass fraction and baryon fraction in the disc, suggesting that the gravitational potential of BRI 1335–0417 is dominated by gas rather than stars and dark halo (see also Carilli et al. 2002; Riechers et al. 2008; suggesting gas dominance in this galaxy).

Recent numerical studies have primarily focused on forming stellar bars out of stellar discs alternating the orbits of disc stars (Bland-Hawthorn et al. 2023) in which the gas delay the process. However, a disc with such a high gas fraction as in BRI 1335–0417 may behave differently. The dominant gas disc could potentially lead to the formation of a gas bar (Barnes & Tohline 2001), as opposed to a dominant stellar disc forming a stellar bar, which then influences the gas kinematics. Theoretically, an axisymmetric 100 per cent gas disc is found to be able to form a self-gravitating stable gaseous bar structure (Cazes & Tohline 2000). This scenario seems plausible explaining the prevalent high-redshift bar structures seen in far-infrared emissions (Gullberg et al. 2019; Hodge et al. 2019) and another high-z barred galaxy at z = 4.3 with an extreme gas fraction (Smail et al. 2023). However, whether the stars forming out of a gaseous bar can lead to a stellar bar as we see in later epochs (Guo et al. 2023; Le Conte et al. 2023) remains an open question, as only idealised simulation experiments have been conducted thus far.

3.4 Second order disc kinematics

We analyzed the [C ii] gas velocity field of BRI 1335–0417 (Fig. 3a) using kinemetry (Krajnović et al. 2006)⁴ to determine if there is a dynamical imprint from a bar, recent interaction, and or inflow/outflow in the gas kinematics.kinemetry expands the velocity field profile v(a, θ) along ellipses into Fourier series,

$$\begin{eqnarray} v(a, \theta) = A_0(a)+\sum _{m = 1}^{N} A_m(a)\sin (m\theta)+B_m(a)\cos (m\theta), \end{eqnarray}$$

(2)

where the ellipses are defined by the semimajor axis a, axis ratio q, and PA from which azimuthal angle θ is measured. As the velocity field of a rotating disc can be expressed by v(a, θ) = A₀ + B₁(a)cos (θ) (Krajnović et al. 2006), the deviation from pure rotational motion manifests as higher order coefficients (A₁, A₂, B₂,...). By this approach, we can investigate non-circular velocities in a non-parametric manner, in contrast to the dynamical modelling method used in Tsukui & Iguchi (2021), where they assumed a pure circular rotation and bulge-disc mass profile.

We performed the kinemetry expansion of the azimuthal velocity profile up to the fifth-order Fourier terms (m = 5) along concentric ellipses. We used ellipses with a PA 7.6° and an axis ratio q = 0.79, which are determined as the global kinematic PA using kinemetry and the photometric axis ratio of the dust continuum distribution in Tsukui & Iguchi (2021), consistent with the independent measurements by Roman-Oliveira, Fraternali & Rizzo (2023). The centre of each ellipse is fixed at the peak position of the continuum image, which matches with the centre of rotation [see Fig. 3A Tsukui & Iguchi (2021)] and the optical quasar position (Tsukui et al. 2023b).

Fig. 3(b) shows the best-fitting expansion with odd-term harmonics m = 1, 3, 5. The data-model residual (Fig. 3d) reveals the significant even (m = 2) component, with two redshifted regions and two blueshifted regions in symmetric positions. The redshifted residuals spatially coincide with the two-armed spiral morphology in the [C ii] and FIR continuum maps. The features are already visible in the velocity field before subtraction in Fig. 3(a); the redshifted/blueshifted high velocity/low velocity are aligned with spiral arms.

By definition, the m = 2 mode cannot be compensated by any odd modes which are the orthogonal basis. This is illustrated by Figs A1–A3. Fig. A1 shows the ellipses along which the azimuthal velocity profiles are extracted by kinemetry. Figs A2 and A3 show (top) the azimuthal velocity profiles v(θ) for ellipses at radii 0.54 and 0.68 arcsec with the best-fitting circular velocity A₀ + B₁cos θ and the circular velocity residual v(θ) − A₀ + B₁cos θ with higher order expansion with only odds terms and odds plus even terms (bottom), clearly showing that m = 2 mode is required to reproduce the data profile.

We also confirm that the m = 2 mode, with the amplitude of up to ∼ 30–40 km s⁻¹, cannot be diminished by changing the PA and axis ratio of the kinematic ellipses by ±20° and ±0.1, respectively (Fig. A4). The addition of the even components fully characterizes the velocity field of the BRI 1335–0417 [see the best-fitting harmonic expansion in Fig. 3(c)] leaving only a small data-model residual over the disc shown in Fig. 3(e). Recently, Bagge et al. (2023) also explored the importance of the even components in spatially resolved gas kinematics for galaxies at intermediate redshift z ∼ 0.3.

3.5 Interpretation of the m = 2 mode in velocity field

We conclude that the observed m = 2 kinematic mode in the velocity residual is less likely due to inflow or outflow, as symmetric flow would show line of sight velocity residuals with the opposite sign in symmetric positions (Di Teodoro & Peek 2021; Genzel et al. 2023), leaving the peculiar possibility where the residual on one side is due to outflow while the opposite side is due to inflow. Similarly, it is unlikely due to in-plane motion induced by the dynamical influence of spiral and bar structures (e.g. Grand et al. 2016; Monari et al. 2016; Gómez et al. 2021). The analytical model suggests a m-fold symmetric mass density structure would only induce m − 1 or m + 1 velocity field distortions (Canzian 1993; Schoenmakers, Franx & de Zeeuw 1997), inconsistent with the observed m = 2 mode velocity residual for a m = 2 spiral and bar density structures in the galaxy. Also, the m = 2 mode velocity residual cannot be attributed to a PA change of the disc due to disc tilting.

Using an analytic model for a galaxy with similar rotational velocity and inclination to BRI 1335–0417, Gómez et al. (2021) show that the radial velocity induced by the spiral density structure is small (≪10 km s⁻¹) compared to the disc vertical velocity caused by an interaction. Assuming the line of sight velocity residual of the m = 2 mode, ∼ 35 km s⁻¹ (Fig. 3d), is purely due to in-plane radial motion, the induced radial velocity would be 57 km s⁻¹ after inclination correction. With and without inclination correction the magnitude of the velocity residual seems too large to attribute for the radial motion caused by the spiral structure (e.g. Grand et al. 2016; Monari et al. 2016).

A more plausible explanation for the observed m = 2 mode signature is vertical motion relative to the disc (disc bending mode) induced by external forces, such as recent interactions with satellite galaxies (e.g. Gómez et al. 2017; Bland-Hawthorn & Tepper-García 2021) or misaligned gas accretions (Khachaturyants et al. 2022). Strikingly, the kinematic m = 2 mode in BRI 1335–0417 spatially coincides with the spiral structure in the intensity map, as shown in Fig. 3(d). This overlap is consistent with the vertical wave of stellar and gas discs induced by satellite interactions in some simulations (Tepper-García, Bland-Hawthorn & Freeman 2022). The simulations reveal the m = 1 mode in the disc vertical velocity develops ∼100 Myr after the Milky-Way-like galaxy is perturbed by the high-speed encounter of the Sagittarius dwarf system, followed by the development of the m = 2 mode vertical velocity as soon as 200–400 Myr (Bland-Hawthorn & Tepper-García 2021; Tepper-García, Bland-Hawthorn & Freeman 2022). The vertical velocity pattern initially aligns with the spiral pattern in density and then decouples.

The observed co-spatial spiral arm patterns in velocity residual and intensity (Fig. 3d) are consistent with the simulation of bending waves induced by interactions, or at least imply a common cause for the two. Given that BRI 1335–0417 has a gas-rich, dynamically hot disc with velocity dispersion of ∼70 km s⁻¹ (Tsukui & Iguchi 2021), conditions typically unfavourable for spontaneous spiral arm formation (Elmegreen & Elmegreen 2014), it is reasonable to interpret that the external perturbation responsible for the m = 2 bending mode induced the observed spiral arms even in a dynamically hot disc (Law et al. 2012) and potentially formed bar (Łokas et al. 2014, 2016). The disc is turbulent but still baryon-dominated and gravitationally unstable (Tsukui & Iguchi 2021), making it prone to forming such substructures through perturbations (Law et al. 2012).

The gas velocity dispersion, σ, of BRI 1335–0417 is estimated to be |$71^{+14}_{-11}$| km s⁻¹ (Tsukui & Iguchi 2021). This value is relatively high compared to average values at lower redshifts (45 km s⁻¹ at z = 2.3, 30 km s⁻¹ at z = 0.9; Übler et al. 2019) but is in agreement with the scatter seen around cosmic noon (Kassin et al. 2012; Wisnioski et al. 2015; Übler et al. 2019; Wisnioski et al. 2019). Using the EAGLE simulations, Jiménez et al. (2023) highlights the importance of the gas accretion rate and its relative orientation to the disc in driving the redshift evolution of disc velocity dispersion. The gas accretion (and accompanying satellites) with a large angle relative to the disc plane effectively increases the velocity dispersion of the disc. In general bending mode (buckling instability) is suppressed for a disc with isotropic velocity dispersion (Toomre 1964; Araki 1985; Merritt & Sellwood 1994, see for review Sellwood 2013) which may hold for high-redshift galaxies (Genzel et al. 2017; Genzel et al. 2023). The bending wave induced by the vertical perturbation of the high-angle accreters may be immediately damped and contribute to the kinetic energy of the disc velocity dispersion.

4 CONCLUSION

We identify an elongated bar-like structure in the z = 4.4 spiral galaxy BRI 1335–0417 using both [C ii] and rest-frame 160 |$\, \mu \mathrm{m}$| FIR images by fitting ellipses to the isophotes. The variation of ellipticity and PA shows the characteristic profile of barred galaxies. The identified bar, |$3.3^{+0.2}_{-0.2}$| kpc long in radius, appears to bridge the two arm spirals identified in the previous study (Tsukui & Iguchi 2021). When compared to the disc scale length, the bar length is larger than stellar bars in nearby galaxies (Erwin 2019).

The [C ii] bar is more elongated than the FIR bar and displaced towards the leading edge of the FIR bar (at the downstream side of the gas motion), assuming the spiral arms of BRI 1335–0417 are trailing spiral arms. If we attribute the structural difference of FIR and [C ii] to the fact that FIR emission can depend on the distribution of both dust (ISM roughly traced by [C ii]) and young massive stars heating the dust, the observations seem to align with the established picture from both observations and simulations – the compressed gas flow (traced by [C ii]) forms a straight line into the centre of the galaxy along the leading edge of the young stellar bar (more traced by FIR emission).

The galaxy is shown to have a significant gas fraction (>73 per cent), hosting the star-forming gas bar, which prevents us from assuming the presence of a stellar bar dominating the gravitational potential and influences the gas bar kinematics. Along with this object, the abundant bar-like structures in FIR continuum observed in bright sub-millimetre galaxies (Gullberg et al. 2019; Hodge et al. 2019) and another example high-redshift barred galaxy at z = 4.3 with an extremely high gas fraction (Smail et al. 2023) may prompt a different perspective on the bar formation and dynamics rather than simulations focused on how disc stars form stellar bars. The dominant gas disc, by itself gravitationally unstable, potentially leads to the formation of a gas bar (Barnes & Tohline 2001), as an axisymmetric 100 per cent gas disc is theoretically proven to form a self-gravitating stable gaseous bar structure (Cazes & Tohline 2000).

By applying a harmonic expansion of azimuthal profiles of the [C ii] velocity field with kinemetry we reveal a dominant m = 2 mode with the amplitude of up to ∼30–40 km s⁻¹ in the velocity residual. The m = 2 mode cannot be explained by large-scale inflow/outflow or the non-circular motions caused by bar/spirals, which generally produce odd-numbered modes (e.g. Canzian 1993; Schoenmakers, Franx & de Zeeuw 1997). Therefore we interpret the m = 2 mode as vertical motion relative to the disc (disc bending wave). The redshifted velocity regions contributing to the m = 2 mode coincide spatially with the spiral arms, consistent with the scenario where the galactic disc is perturbed by an interaction creating the initially co-spatial spiral density wave and vertical bending mode in the disc (e.g. Bland-Hawthorn & Tepper-García 2021; Tepper-García, Bland-Hawthorn & Freeman 2022).

As it is assumed that a gas disc is stable for the bending mode (buckling instability), the m = 2 mode velocity perturbation and spiral arms are likely to be induced by a recent interaction and/or gas accretion. It is natural to suppose that such an interaction would also activate the high star formation activity.

For the first time, this study detects an unambiguous m = 2 mode in a disc at high redshift (z > 4), lending support to the use of disc seismic ripple as evidence of a recent strong perturbation. With further numerical simulations tuned for galaxies with realistic gas, halo, and stellar masses and gas accretion histories, the m = 2 mode may provide a useful constraint on the exact timing and primary origin of external perturbations.

ACKNOWLEDGEMENTS

TT is grateful to the conference organizers of Galactic Bars: driving and decoding galaxy evolution held in Granada, Spain in 2023 July, which greatly helped in writing this paper and Trevor Mendel for encouraging me to attend. TT also thanks Camila de Sá-Freitas, Karin Menedez Delmestre, Ewa Luiza Lokas, and Andreas Burkert for insightful discussions. This research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. Data analysis was partly carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00394.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC 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. Finally, TT thanks Takaho Masai for his kind support at NAOJ.

DATA AVAILABILITY

The ALMA data we use in this work are publicly available in https://almascience.nrao.edu/aq/.

Footnotes

1

The elongated bar structure has a random orientation relative to the disc PA.

2

With several assumptions such that the dust is solely heated by the young massive stars and archived single equilibrium temperature. Stellar initial mass function and dust geometry are known.

3

We assumed the dispersion supported disc with finite thickness and the scale radius R_d of [C ii] (Tsukui & Iguchi 2021).

4

See the derivation of [C ii] velocity field (Tsukui et al. 2023b), Briefly, we derived the velocity field by fitting a Gaussian function (including Hermite parameters h3 and h4; Cappellari 2017) to the [C ii] spectrum at each pixel, where the three channels are available with the emission detected more than 3σ.

APPENDIX A: SUPPLEMENTAL INFORMATION OF kinemetry ANALYSIS

In this section, we provide the detailed results of the kinemetry analysis to verify the presence of m = 2 mode. Fig. A1 shows the ellipses to extract the azimuthal profile of the velocity field. Figs A2 and A3 illustrate the harmonic expansion of the azimuthal velocity profiles at two specific radii: 0.54 and 0.68 arcsec. For the sake of visualization, the pure rotational motion (top: blue line) is subtracted from the original velocity profile (top: black points) and we have shown a higher-order expansion for the residual (bottom: black points). The expansion involves odds terms with or without even terms up to m = 5 modes (shown as red solid and dashed lines respectively). Even terms, especially m = 2 mode are required to reproduce the velocity fields. The dominance of the m = 2 relative to the other modes obtained by expansion with full terms (odds plus even) is shown in Fig. A4; the m = 2 is dominant after 0.4 arcsec and, opposed to other odds modes, is not diminished by the vast change of PA and ellipticity of sampling ellipses.

Figure A1.

Blue points: Ellipses and azimuthal data sampling of each ellipse used for kinemetry analysis overlain on the BRI 1335–0417 velocity field. Black solid lines show ellipses with a radius of 0.54 arcsec and 0.68 arcsec whose velocity profiles are shown in Figs A2 and A3.

Figure A2.

Top: Azimuthal velocity profile v(θ) at a radius of 0.54 arcsec (black points) with the best-fitting circular velocity A₀ + B₁cos θ (blue solid line). Bottom: The residual velocity after subtracting the circular velocity v(θ) − A₀ + B₁cos θ (black points) and the best-fitting expansion with higher order harmonics. The red dashed line is the best fitted with only odd terms (m = 1, 3, 5) and the solid line is the best fit including even terms (m = 1, 2, 3, 4, 5).

Figure A3.

Same as Fig. A2 but for a radius of 0.68 arcsec.

Figure A4.

Solid lines; Fourier coefficients (m = 2, 3, 4, 5) of azimuthal velocity profile as a function of the semimajor axis of the ellipse obtained by kinemetry, using concentric ellipse with an axial ratio q = 0.79 and a PA of 7.6°. Shaded regions; the maximal variation allowed when varying the axial ratio (±0.1) and PA (±20°) of the ellipse used sampling the azimuthal velocity profile.