Dynamical analysis of Maclaurin disk with velocity dispersion and its influence on bar formation

The initial condition in this study is the isolated two-dimensional Maclaurin disk whose surface density as a function of radial distance $\Sigma(r)$ is given by

where $\Sigma_{0}$ is the density at the center and $R_{0}$ is the disk radial size (see, e.g., Freeman 1966; Schulz 2009; Roshan et al. 2016 for some examples of the structure analysis). By mass normalization, $\Sigma_{0}$ relates to the disk mass $M$ and $R_{0}$ by

This surface density yields the potential energy

for $r<R_{0}$ where

in which $\Omega_{0}$ is the disk uniform angular frequency. Correspondingly, the uniform epicyclic frequency $\kappa_{0}$ can be calculated by the following definition

For the stellar disk analogue, the uniformly rotating disk of stars in dynamical equilibrium having the surface density (1) is known as the Kalnajs disk (Kalnajs, 1972). The distribution function of that disk for the rotational frequency $\Omega\leq\Omega_{0}$ as a function of the radial distance $r$, the radial velocity $v_{r}$ and the tangential velocity $v_{\theta}$ is given by

which is valid for the positive argument in the square root. Otherwise, $F=0$. In the expression (2.1), we choose $R_{0}$, $\Sigma_{0}$ and $\Omega_{0}$ to be the units of length, surface density and angular frequency, respectively. Note that the distribution function (2.1) can also be expressed as a function of the energy per unit mass $E=\frac{1}{2}(v_{r}^{2}+v_{\theta}^{2}+r^{2})$ and the orbital angular momentum $J=rv_{\theta}$. This distribution function is also known as the $\Omega$-model and it can be proved that the density profile for any $\Omega$ takes the form of Eq. (1). We limit ourselves to the prograde rotation, i.e., $0\leq\Omega\leq 1$. The case where $\Omega=0$ corresponds to the non-rotating disk while for $\Omega=1$, the disk is purely rotating with no random motion. The mean radial and tangential velocities of the $\Omega$-model read

respectively. Specific to this disk family, the velocity dispersion is isotropic thus the radial and tangential velocity dispersions are identical and they are given by

For the consideration of the stability problem, the radial velocity dispersion is typically parametrized by the Toomre’s $Q$ parameter by

The numerical value $Q=1$ is the critical value for the disk stability against the local axisymmetric perturbations at any wavelength. For simplicity, we set $Q$ to be constant. This constant-$Q$ profile yields $\Omega$ that depends on $Q$ only. The relation between $Q$ and $\Omega$ can be obtained via $\sigma_{r}$ as functions of both parameters in Eq. (8) and (9). We finally obtain the relation between $Q$ and $\Omega$ to be

which is used in constructing the $\Omega$-model disk for any $Q$. From the relation (10), the purely rotating disk yields $Q=0$ whereas the non-rotating disk has $Q=1.696$ that marks the upper limit for this disk family.

In continuity with the disk physical and kinematical details in Sec. 2.1, we will analyze the validity of the nearly-circular orbit approximation when $Q$ is present. When the Toomre’s criterion requiring $Q\geq 1$ is applied on the disk of particles, most particle orbits deviate significantly from the circular shape due to the random velocity component. On the other hand, the Lindblad analysis assumes the nearly circular orbit of particles where the frequencies are calculated. To resolve this, we first recall the axisymmetric Jeans equation that describes the interplay between the gravitational potential and the velocity moments and we express it in the following form

where $\bar{v}_{\theta}=\Omega r$ corresponds to the mean tangential velocity written in terms of $\Omega$ by Eq. (7). Note that the first term on the right-hand side stands for $\Omega_{0}^{2}$. We re-arrange Eq. (11) to be

where, specific to the Maclaurin/Kalnajs disk with uniform $Q$,

We consider $\tilde{\Phi}$ as the effective potential which corresponds to the combination of the gravitational potential and the pressure potential arising from the velocity dispersion. Eq. (12) can be interpreted by that the particle orbit remains circular with effective orbital frequency $\Omega$ in the effective potential. The random velocity component is now regarded as it oversees the pressure potential field leading to the pressure force

where the term in the derivative stands for the pressure. From the definition (14), the presence of the pressure gradient always leads to the pressure force regardless of the model of velocity dispersion profile which can be chosen to depend on $Q$ or not. Note that if $Q=0$, the disk is pressure-less and it is prone to the local gravitational instability.

From the effective potential (13), we define the effective epicyclic frequency

This expression differs from the standard definition of the epicyclic frequency which takes only the gravitational potential into account. By substituting $\tilde{\Phi}$ from Eq. (13) into this expression, we finally obtain

This effective epicyclic frequency is conceptually extended from the definition of the asymmetric drift derived from the Jeans equation (11), in which the effect of the velocity dispersion shifts the orbital frequency. Here, we make a hypothesis that the same effect can similarly shift the epicyclic frequency. The main advantage of choosing the Maclaurin/Kalnajs disk in this study is that the frequencies are uniform throughout the disk. The inspection of the deviation from the original frequencies and its variation with $Q$ is straightforward. The effective potential and the corresponding effective frequencies are our central frames of references to inspect the disk evolution in simulations for different $Q$. Note that in a more realistic disk, anisotropic velocity dispersion gives rise to the third term on the right-hand side of Eq. (11). We also note possible radial dependences of $\omega$, $\kappa$ or $Q$ profiles for a more general disk. Despite those complications, the effective potential approach can still be applied without loss of generality as the Jeans equation (11) allows any form of those functions. The effective potential might have a more complicated form given those complexities.

For the initial disk of particles in isolation, we generate it with particle positions and velocity moments prescribed by the expressions in Sec. 2.1. For the random velocity component, we draw it from the cut-off Gaussian distribution to avoid the high velocity tail. We simulate disks of particle number $N=1.024\times 10^{6}$ and $Q=0.9,0.95,1,1.1,1.2,1.35,1.5$ and $1.65$. Although the Toomre’s criterion requires $Q\geq 1$, we nevertheless inspect disks with $Q=0.9$ and $0.95$ to test the validity of the effective potential hypothesis when the local axisymmetric instability is involved.

In addition to the isolated disk which is a system of interest to test the effective potential hypothesis, the Maclaurin disk in spherical halo potential can be of equal interest for the question of the disk stability when the dark matter halo is present. We adopt the Hernquist profile to represent the dark matter halo in which the Maclaurin disk resides whose potential as a function of radius $r$ is given by

where $M_{h}$ is the halo mass and $a_{h}$ is the halo scale-radius. We choose $M_{h}=25M$ and $a_{h}=4R_{0}$. We perform the simulations of Maclaurin disks in both rigid and live halos. With the Hernquist potential, the composite disk-halo potential causes the Maclaurin disk to rotate differentially in dynamical equilibrium. The rotation curve and the radial profile of the angular frequency of the Maclaurin disk in Hernquist potential with our choices of mass and scale-radius is shown in Fig. 1. The tangential velocity is presented in units of $V_{0}\equiv R_{0}\Omega_{0}$ that corresponds to the velocity at the isolated disk cut-off radius, while the angular frequency is presented in units of $\Omega_{0}$. Testing the bar instability in this composite disk-halo system not only allows us to investigate the effect from spherical halo but we are also able to inspect the effect from the differential rotation on the bar instability. This composite system is more relevant to the real galaxies because they reside in dark matter halos and they rotate differentially. We adopt the prescription of Hernquist (1993) to construct the disk-halo system in equilibrium which bases on the Jeans equations of the composite potential. The disk rotational velocity structure is constructed from the axisymmetric Jeans equation with a chosen velocity dispersion profile given by

where $\Omega_{0}$ is now calculated from the composite disk-halo potential. Specific to the live halo scheme, the halo velocity dispersion $\sigma_{h}$ is isotropic and its radial profile is numerically determined from the spherically symmetric Jeans equation

where $\rho_{h}(r)$ is the halo radial density profile and $M_{tot}(r)$ is the total mass enclosed inside $r$. Random velocity components are drawn from cut-off Gaussian distribution. We keep the constant-$Q$ profile as employed in the case of isolated disk. In the two scenarios, disks consist of equal number of particles as in the isolated case, that is $N=1.024\times 10^{6}$, in order to keep the same magnitude of the Poisonnian fluctuations, while for the live halo simulations, the halo consists of $3.072\times 10^{6}$ particles. The live halo is truncated at $a_{h}$ which is $4$ times the disk cut-off radius. We simulate the evolutions of disk-halo systems with $Q=1.8,2.0,2.2$ and $2.25$.

Evolution of disk of particles is simulated by GADGET-2 (see Springel et al., 2001; Springel, 2005). We employ the same choices of units as in the analysis in Sec. 2.1, i.e., the disk cut-off radius $R_{0}$, the rotational frequency of the cold disk $\Omega_{0}$ and the disk surface density at the center $\Sigma_{0}$ for the units of length, frequency and surface density, respectively. We choose the rotation period of the purely rotating disk, or

as the unit of time. To give an example of the rotation period, a Maclaurin disk of Milky Way mass and radius, i.e. $M\sim 10^{12}\ M_{\odot}$ and $R_{0}\sim 30\ \text{kpc}$, has the rotation period equal to $317\ \text{Myr}$. For velocity unit, we adopt $V_{0}\equiv R_{0}\Omega_{0}$, as employed in Fig. 1. The gravitational force calculation is facilitated by the tree code. The gravitational force in the code is spline-softened with softening length $\epsilon\sim R_{0}/2000$. We adjust the opening angle $\theta$ equal to $0.5$ for simulations of isolated disks and disks in rigid halo while it is fixed to $0.7$ for the live halo simulations. The integration time-step is controlled to be not greater than $T_{0}/25000$ for the cases of isolated disks and disks in rigid halo. For live halo simulations, it is controlled to be below $T_{0}/8000$. With this simulation setting, the accuracy of the integration is such that the deviation of the total energy from the initial value at any time is less than $0.6\%$ until the end of simulation that is fixed to $7.57T_{0}$ for all cases.

The strength of bi-symmetric modes can be evaluated by the $m=2$ Fourier amplitude as a function of radius $\tilde{A}_{2}(r)$ given by

where $a_{2}$ and $b_{2}$ are the Fourier coefficients at $r$ computed from

The summation $j$ includes only particles in the annulus of radius $r$, each of which has angular position $\varphi_{j}$ and mass $\mu_{j}$, with total number $N_{r}$. In the expression (21), $A_{0}$ is the corresponding Fourier amplitude of $m=0$ modes. The bar strength $A_{2}$ is defined as the maximum $\tilde{A}_{2}$ within the bar reach, or

In addition, the radial phase change of $\tilde{A}_{2}$ is capable of capturing the presence of two-armed spiral modes. The winding degree of the two-armed modes is evaluated in terms of the logarithmic pitch angle $i$ that can be written as

where $\tilde{\phi}_{2}(r)$ is the phase of the $m=2$ modes at $r$. Our sign convention is that $\cot i$ is positive and negative for the leading and trailing spiral arms, respectively.

In the entire of Sec. 4, we focus on the evolution of isolated Maclaurin disks in the framework of the effective potential hypothesis. We start the section by inspecting the overall evolutions. Shown in Fig. 2 is the time evolution of bar strength $A_{2}$ for disks with different $Q$. The plots for each of $Q=1,1.35$ and $1.5$ are drawn from a single realization among the three different realizations. For $Q=1.65$, only one realization is simulated. In addition, disk surface density maps at different time for some bar-unstable $Q$ are depicted in Fig. 3. From the $A_{2}$ plot, it is evident that if $1\leq Q\leq 1.5$, the disks are bar-unstable: $A_{2}$ grows from the initial state before the growth is put to an end when $A_{2}$ reaches the first maximum within a few $T_{0}$ and it remains high afterwards. The evolution of $A_{2}$ is in accordance with the development of persistent barred structures in Fig. 3 for both $Q=1$ and $1.5$. As $Q$ reaches $1.65$, the bar formation is not triggered as $A_{2}$ remains at the initial noise level. This implies that the critical value of $Q$ to stabilize the disk is in the range of $1.5-1.65$. That range is close to the value estimated by the stability analysis against the two-armed (or bar-like) modes for this disk family by Kalnajs & Athanassoula-Georgala (1974) who suggested that the numerical value $\Omega<0.5072$, which is equivalent to $Q>1.46$, can suppress those modes. However, when concerning the same problem in a disk of particles, there are two factors that might affect the critical value. The first one is the force softening in simulations which acts as the stabilizer. Secondly, the random Poissonian noises embedded in the particle system might be another source of the perturbations in addition to the imposed multiple-armed modes perturbations as the original analysis took into account. Another analysis of the uniformly rotating disk was also performed by Christodoulou et al. (1995), using a diagnostic involving more disk physical and kinematics parameters. This analysis led to the same conclusion that the disk with lower rotation frequency tended to be more stable.

With a closer look on $A_{2}$ and morphological evolutions, we remark distinct evolution patterns for $Q=1$ and $1.5$. For $Q=1$, $A_{2}$ rapidly increases at first before it oscillates strongly afterwards. At $t=1.96$, or the time around the first peak, we observe the spiral pattern in opposite directions formed by the particles driven outwards by the bi-symmetric force. The growth of the spiral pattern can be explained as follows. In a bar-unstable disk, the global linearly unstable two-armed modes are the fastest growing modes, which are initiated from the Poissonian noise. These modes dominate the early stage and grow exponentially because of its linear nature (Dubinski et al., 2009; Fujii et al., 2018; Collier, 2020). The linear instability marks an end when the two-armed modes, as quantified by $A_{2}$, reach their peak. Then, the disk remains in a barred state by the non-linear bar instability. In this stage, particles are trapped by the bi-symmetric bar potential induced by the remnant of the preceding two-armed modes. This can be seen by the persisting bar in the following snapshots which lasts for time much longer than the formation time scale. This result is in agreement with the pioneering simulation of Hohl (1971). For the $Q=1.5$ disk, $A_{2}$ increases more slowly and it oscillates more weakly in the barred state. That $A_{2}$ evolves more gently can be seen in Fig. 3 where no visible spiral trace is observed at $t=3.79$ which is the time around the $A_{2}$ peak. The case where $Q=1.35$ exhibits the intermediate evolution pattern of $A_{2}$. For $Q=0.9$ and $0.95$, the bar is able to be formed amidst the local instability and the overall disk evolution is qualitatively similar to that for $Q=1$: the two-armed spiral modes grow rapidly leading to the robust bar.

We further investigate the unstable two-armed modes by measuring $\cot i$ profile where $i$ is the logarithmic pitch angle around the time of $A_{2}$ peak for different $Q$ and the results are shown in Fig. 4. From the plot, we capture the trailing spiral two-armed modes in all cases beyond the region where $\cot i\sim 0$, which corresponds to the barred region, although the spiral arms are not visible in $Q=1.5$ configuration. It turns out that the underlying processes starting from the growth of the two-armed modes to the bar instability are generic for all $Q$. Increasing $Q$ makes the perturbation evolve more slowly.

The prolonged formation time and the gentleness of the bar-forming process in higher-$Q$ disks can otherwise be understood in the dynamical point of view that regards the interplay between two competing forces. First, the bi-symmetric modes perturb the disk by means of growing the two-armed spiral structure outwards that compresses the disk environment in the process. The compressed background then exerts the counteractive pressure force, i.e., $\Delta F_{P}$, in response to the growth. As the pressure is an increasing function of $Q$, $\Delta F_{P}$ then increases with $Q$ in a similar way. This explains why the high-$Q$ disk develops the bar more slowly because the disk environment counteracts the growth with a stronger force. The fact that the bar is established for $Q\in[0.9,1.5]$ proves that the two-armed perturbative forces are stronger than the counteractive pressure forces from the disk background. Conversely, the situation for $Q=1.65$ disk implies that the perturbative force cannot surpass $\Delta F_{P}$.

In this section, we examine the $Q$-dependence of the resonance mechanism responsible for the bar formation and the resulting bar kinematics. To do this, we plot the ensemble-averaged angular frequencies of the $m=2$ modes before and after the fully formed bar for different $Q$ in the top panel of Fig. 5. These frequencies correspond to the dominant Fourier frequencies of the bar phase calculated in the time windows of widths $0.63T_{0}$ and $2.52T_{0}$ before and after the first peak of $A_{2}$, respectively. The first frequency corresponds to the angular frequency of the linearly unstable two-armed modes while the latter frequency represents the bar pattern speed. The disk effective rotational frequency as a function of $Q$ from Eq. (10) is provided for comparison. In the bottom panel of Fig. 5, we plot the ensemble-averaged fraction of number of particles that are trapped by the resonance during the bar-forming stage for different $Q$. More specifically, trapped particles are those with orbital frequencies lie within $[0.95\omega_{2},1.05\omega_{2}]$ near the $m=2$ modes axis where $\omega_{2}$ is the $m=2$ modes angular frequency. The fraction is determined in the time window of width $0.76T_{0}$ around the first peak of $A_{2}$.

When $Q\geq 1$, the two-armed modes frequency in the bar-forming stage decreases, on average, with $Q$ and it is close to the effective rotational frequency from Eq. (10). That these two frequencies are close to each other implies that the corotation resonance is responsible for the growth of the linearly unstable two-armed modes for all $Q$. After the bar is fully formed, the bar angular frequency is significantly declined from the value before in all cases. This signifies the angular momentum loss during the bar formation. The correlation between the bar angular frequency and the initial effective frequency is still retained. Although the correlation between the bar angular frequency and $Q$ has been pointed out in past studies (see, e.g., Athanassoula (2003); Hozumi (2022)), here we are able to give a more detailed explanation of that $Q$-dependence by the effective frequency hypothesis. On the contrary, the cases where $Q<1$ exhibit the anomalies from the other cases. For $Q=0.95$, the rotational frequency of the fully formed bar does not decrease much from the initial frequency. As $Q$ falls to $0.9$, the frequencies deviate from the trend by a wide margin. We also note relatively large error bars for these $Q$. The inconsistencies with results for $Q\geq 1$ disks are attributed to the unstable local perturbations originated from the random Poissonian noises that disturb the collective motion differently from realization to realization. In this regime, we conclude in the breakdown of the effective potential assumption.

Considering now the fraction of number of trapped particles, the fraction tends to decrease with $Q$ for $Q\geq 1$. This can be explained by that the velocity is more dispersed when $Q$ is higher, so there are less particles close to the $m=2$ modes angular frequency. At the other end of the plot where $Q<1$, we notice the significant drop of the fraction. This is another indication of the involvement of the unstable local fluctuations. Such fluctuations perturb the orbital motion of particles so that some of them are dragged away from the resonance and do not become part of the bar.

In past literatures, the loss of bar angular momentum has been conjectured to be caused by the dynamical friction acting on the rigid bar when it is spinning through the background particles (Weinberg, 1985). Another hypothesis is that the angular momentum loss is regulated by the transfer between the inner and outer resonances (Little & Carlberg, 1991; Athanassoula, 2003; Martinez-Valpuesta et al., 2006). In our case, the mechanism of the angular momentum loss is that the particles are trapped by the corotation resonance and some of them are transferred outwards, carrying the angular momentum away from the disk center to the outer spiral components that rotate more slowly.

In this section, we focus on the evolution of $A_{2}$ further from the formation phase. From Fig. 2, we remark the oscillation of $A_{2}$ after the bar is fully developed. The oscillation is more noticeable when $Q$ is close to $1$ whereas it is much weakened, but still observable, for $Q=1.5$. We speculate that this oscillation is originated from the radial oscillation of particles forming the bar. Therefore, the underlying mechanism of this $A_{2}$ oscillation is the epicyclic motion. We will inspect the nature of that oscillation but the difficulty is that $A_{2}$ oscillates non-steadily that makes the determination of the oscillation frequency difficult. To resolve this, we linearly adjust $A_{2}$ relative to its adjacent principal maximum and minimum so that it oscillates in $[-0.25,0.25]$ while maintaining the same frequency. An example of the pre- and post-processed $A_{2}$ for a selected case in Fig. 2 is shown in Fig. 6. The Fourier frequency of the processed $A_{2}$ is then computed in the time window of width $2.53T_{0}$ just after the first peak and the ensemble-averaged frequency is plotted in Fig. 7, presented in units of the unperturbed epicyclic frequency $\kappa_{0}$. In the plot, we also provide the averaged radial oscillation frequency of particles in nearly circular orbits in the first few $T_{0}$ for each $Q$ to represent the initial epicyclic frequency. Those particles correspond to particles with initial radial and tangential random velocities relative to their unperturbed orbital velocity below $0.02$ and $0.025$ for $Q=0.9-1.2$ and $Q=1.35-1.5$, respectively. In addition, the averaged radial oscillation frequency of particles in more eccentric orbits, i.e., those with either of the radial or tangential random velocity ratio between $0.15-0.3$, is also plotted for consideration. The effective epicyclic frequency as a function of $Q$ from Eq. (16) is provided for comparison.

For $Q\gtrsim 1$, the effective $\kappa$ explains reasonably well the radial oscillation frequency of nearly circular orbit particles in the first few $T_{0}$. It slightly underestimates the measured $\kappa$ but the decrease with $Q$ is evident. There is not only the disk angular frequency that is shifted by finite $Q$ but the epicyclic frequency is also shifted in the similar manner. The difference between the theoretical effective frequency and the simulated frequency could be explained by that the derivation takes only the initial density and velocity structures into calculation. Disks that evolve might re-adjust their internal configurations so they alter from the initial profiles. For particles in more eccentric orbits, the calculated frequencies are higher than those of particles in nearly-circular orbits but the decrease with $Q$ is still evident and these frequencies are well below $1$. It turns out that the effective epicyclic frequency provides a better description than adopting $\kappa_{0}$ for a disk with non-zero $Q$. When $Q<1$, both calculated frequencies are lower than their $Q=1$ counterparts which contradicts our assumption. These anomalies are, similar to the situation for angular motions in Fig. 5, attributed to the local instability that perturbs the particle epicyclic motion. In past literatures, it was a usual practice to derive the epicyclic frequency only from the disk gravitational potential, both in simulations (Dehnen, 1999; Rautiainen & Salo, 2000; Saha & Jog, 2014; Wu et al., 2016; Vasiliev, 2019) and in observations (Monari et al., 2019; Kawata et al., 2021; Lee et al., 2022). On the other hand, some other studies measured it directly (Athanassoula, 2003; Martinez-Valpuesta et al., 2006; Dubinski et al., 2009). We respond to those points by the fact that the radial motion of a particle in a disk of particles cannot be considered as being under the disk gravity only. It is also subject to the additional pressure force created by the velocity dispersion of particles forming a disk in equilibrium. This pressure force alters significantly and systematically the epicyclic motion as demonstrated by our simulation.

The second point captured from Fig. 7 is that the epicyclic frequency of the fully formed bar is lower than $\kappa$ in the initial phase, although the decline is modest if $Q<1$, and the bar oscillation frequency keeps the same tendency with the $\kappa(Q)$ line. This indicates that the effect from initial $Q$ does not only imprint on the bar angular frequency as demonstrated in Sec. 4.2 but the final bar epicyclic frequency also depends systematically on $Q$. About the decline of the epicyclic frequency, we will seek the theoretical explanation of the transfer of the radial action. We adopt the analysis of Athanassoula (2003) for the angular momentum transfer and we will adapt it for the radial counterpart. First of all, we recall that the distribution function of a flat disk in an equilibrium can generally be expressed as a function of the two action variables, i.e., $f(J_{1},J_{2})$. The first one is the radial action, or $J_{1}\equiv J_{r}$, that is proportional to the epicyclic frequency and the second one $J_{2}$ is the azimuthal action and it is equal to the orbital angular momentum. The action angle $w_{i}$ and frequency $\Omega_{i}$ of the corresponding $J_{i}$ can be derived from the unperturbed Hamiltonian $H_{0}$ as follows

By definition, $\Omega_{1}$ and $\Omega_{2}$ correspond to the epicyclic and orbital frequencies, respectively. Next, we write the disk potential with time-evolving perturbations as

where $\Psi_{0}$ is the unperturbed axisymmetric component, $\psi$ is the non-axisymmetric perturbation and $\varpi$ is the complex perturbation frequency. We separate $\varpi$ into the real part $\varpi_{R}$ that represents the pattern speed and the imaginary part $\varpi_{I}$ that designates the growth rate provided that it is negative. We remind that the analysis in action-angle coordinates is valid if the perturbation grows slowly. The potential perturbation $\psi$ can then be expanded in Fourier series as

where $\psi_{lm}$ is the Fourier coefficient calculated by

We follow the derivation of Athanassoula (2003) to obtain the rate of change of $J_{r}$ and it yields

For the slowly-growing perturbations, i.e., $\omega_{I}\rightarrow 0$, the change of the radial action variable is possible only if

which corresponds to the resonance condition of particles with perturbations frequencies $\varpi$. From Eq. (29), the sign of $\dot{J}_{r}$ depends on the sign of $l$ and $m$ in front of the derivatives. In general, the derivative of $f$ with respect to $J_{i}$ is negative. Therefore, the part that gains the radial action from the bar, or $\dot{J}_{r}>0$, is that with positive $l$ and $m$, which corresponds to the outer Lindblad resonance. From the derivation, the radial action transfer from the bar to the outer component that rotates more slowly is justified. That the epicyclic frequency evolves in time during the bar formation is an important information when considering the bar evolution as the epicyclic frequency is equally important as the angular frequency. In past studies, the measurement of the oscillation frequency of the bar parameter has been performed by, for instance, Miller & Smith (1979) and Hilmi et al. (2020). The latter work concluded that the oscillation was attributed to the interaction with the spiral pattern. According to our simulation, we ascertain that the epicyclic motion takes a major role in that oscillation since it correlates with $Q$ according to the effective potential assumption.

In this section, we inspect the kinematical and structural details of the disk hosting the bar. Shown in Fig. 8 is the disk intrinsic velocity fields on the face-on view, in arbitrary units, for $Q=1$ and $1.5$ at different time which represent the violent and gentle bar-forming processes, respectively. The velocity field is computed by locally averaging the velocities of particles at each position. We consider first the case where $Q=1$. At $t=2.21$, we note the significant radial velocity component on the spiral pattern in the linear instability phase. Later at $t=4.16$, no prominent radial motion is observed anymore on the spiral arm. For $Q=1.5$, the radial velocity near the bar ends is weaker than the $Q=1$ counterpart. This justifies that the radially outward motion in the linear instability phase is much suppressed if $Q$ is higher.

Last but not least, we examine the radial surface density profile of the disk in barred states. Shown in Fig. 9 is the annularly-averaged surface density as a function of radial distance $\Sigma(r)$ for disks with $Q=1$ and $1.5$ at $t=6.31$. We perform the fitting with the exponential decay function for different disk parts. In general, the surface brightness of disk galaxies in observations or the surface density of disks in simulations fall into one of the three types. The Type-I profile represents the single exponential decay while for the Type-II and the Type-III profiles, they exhibit the double exponential decays separated by break radii. The Type-II (or down-tuning) profile has the density beyond the break radius decreasing more steeply than the inner part and vice versa for the Type-III (or up-tuning) profile. From Fig. 9, it turns out that both disks manifest the Type-III profile. The obtained Type-III profile is in accordance with past literatures that also obtained this profile (Mayer & Wadsley, 2004; Debattista et al., 2006). This is another indication of the common process underlying the bar formation in low- and high-$Q$ regimes.

In the entire of Sec. 5, we inspect the bar instability in the composite disk-halo systems in many aspects. We examine first of all if the presence of halo affects the stability criterion for the disk with the same surface density and $Q$ profile comparing to the value reported in Sec. 4.1. Shown in Fig. 10 is the time evolution of $A_{2}$ for disks in rigid and live halos and for various $Q$. In rigid halo simulations, the bar formation is observed until $Q=2$. When $Q$ reaches $2.2$, $A_{2}$ evolves quietly all along the simulation, which indicates that this case is bar-stable. For the disk in live halo, it exhibits a stronger amplification as $A_{2}$ attains, on average, to higher values. The bar formation is still observed until $Q=2.2$ where $A_{2}$ is able to reach $0.2$. When $Q=2.25$, the increase of $A_{2}$ is still seen but it fails to establish the prominent barred structure as the value of $A_{2}\sim 0.1$ is too low to render the apparent barred structure according to the study of Algorry et al. (2017). It is therefore acceptable to classify this case as unbarred.

That the disk in live halo is less stable and it tends to build a stronger bar is possibly due to the particle nature of the halo. That particle halo adds up the finite-$N$ fluctuations to those of the disk, resulting in a stronger seed of the two-armed spiral modes. These additional fluctuations do not exist in the smooth rigid halo. As a consequence, the critical $Q$ for both scenarios slightly differ as the value of $Q=2.2$ proves sufficient to stabilize the disk in rigid halo while it has to be $2.25$ to suppress the amplification of $m=2$ modes in live halo case. Another difference is that the modest amplification of bar modes is still observed in live halo scheme although it is classified as bar-stable.

In continuity with the $A_{2}$ plots, it is evident that the disk is bar-unstable although $Q$ exceeds $1.65$, which was sufficient to stabilize the isolated disk (see Fig. 2). The critical $Q$ for disk-halo systems is somewhere between $2-2.25$, depending on the halo type. Because we control the disk surface density and the $Q$ profile to be the same as those in the isolated cases, the elevation of critical $Q$ is due to the only different factor which is the differential rotation. It engenders the swing amplification and, from that process, the spiral density waves are amplified by the disk shearing in combination with the synchronized epicyclic motion (Julian & Toomre, 1966; Toomre, 1981). Such mechanism enhances the non-axisymmetric spiral modes in addition to the exponential growth of the two-armed modes by the linear instability that was addressed in Sec. 4. This explains why the composite disk-halo system with $Q$ greater than $1.65$ is still prone to the bar instability because the combined procedures amplify the two-armed modes more efficiently. The swing amplification is absent in isolated disks that rotate rigidly.

We compare our stability test with the past renowned criteria of bar stability. We revisit two frameworks which are the Ostriker-Peebles (OP) criterion (Ostriker & Peebles, 1973) and the Efstathiou-Lake-Negroponte (ELN) criterion (Efstathiou, Lake, & Negroponte, 1982). The former criterion proposed the ratio of the rotational kinetic energy to the potential energy as an indicator, i.e.,

The disk is supposed to be bar-stable when $t_{OP}\lesssim 0.14$. Otherwise, the ELN indicator incorporated both the effects from the random motion and the mass concentration and it can be written as

where $v_{max}$ is the maximum tangential velocity, $M_{D}$ is the disk mass and $R_{D}$ is the characteristic disk radius, which is set to $R_{0}$ for our case. The disk is supposed to be stable if $t_{ELN}\lesssim 1.1$. We numerically calculate both indicators for our disk-halo initial states and we obtain $t_{OP}=0.186,0.164,0.139$ and $0.133$, and $t_{ELN}=1.13,1.04,1.018$ and $1.016$ for $Q=1.8,2.0,2.2$ and $2.25$, respectively. From our numerical results, it turns out that our critical $t_{OP}$ is slightly above the proposed value as it is around $0.133-0.164$. On the other hand, our critical $t_{ELN}$ is slightly below the theoretical value, being in the range of $1.016-1.04$. Note that the disk family that we adopt differs from those employed in their works. Different disk physical and kinematical structures might possibly lead to the offset of critical value.

Another factor which plays an important role in the bar evolution is the vertical effect, via the buckling instability by which the kinematics of particles forming the bar evolves vertically (Combes et al., 1990; Aumer & Binney, 2017; Li et al., 2023). Apart from the buckling instability, the effect from disk thickness alone can also affect the bar-forming dynamics and the bar kinematics in many ways. For instance, the bar formation in a thick disk can be delayed (Ghosh et al., 2023) and the established bar tends to be longer and more slowly-rotating (Klypin et al., 2009) than the situation in a thin disk. Moreover, the spiral modes are more short-lived in the thick disk because of the suppressed swing amplification (Ghosh & Jog, 2018, 2022). All of these effects are important to the fate of the entire system consisting of disk, bar and spiral arms. On the contrary, our study employs the razor-thin disk in which the effect from disk thickness is absent. This might be another cause of the discrepancy with other simulations.

In this section, we investigate the bar evolution in live halo to a greater detail comparing to past studies. There were reports of two important mechanisms underlying the bar evolution. The first one is the angular momentum transfer between disk and halo particles that causes the bar slow-down (Athanassoula, 2003; Holley-Bockelmann et al., 2005; Athanassoula, 2014). To address this, time evolutions of disk and halo angular momenta for $Q=1.8$ and $2.25$ are plotted in Fig. 11, in units of $L_{0}\equiv MR_{0}^{2}\Omega_{0}$. The former case exhibits clearly the angular momentum transfer from disk to halo in coherence with the bar development, in line with the result of Dubinski et al. (2009). We further demonstrate that, for the case of $Q=2.25$ in which the bar is not effectively formed, the angular momentum transfer does not take place.

Another important mechanism is the radial heating, or the excitation of star kinematics by non-axisymmetric forces from bar or spiral arms (Minchev & Quillen, 2006; Gustafsson et al., 2016; Ghosh et al., 2023). To verify this, we plot the radial velocity dispersion profile, or $\sigma_{r}(r)$, for $Q=1.8$ and $2.25$ at different time in Fig. 12. For $Q=1.8$, $\sigma_{r}$ profile in the bar region is elevated in accordance with the progression of $A_{2}$. On the contrary, $\sigma_{r}$ profile alters modestly for $Q=2.25$ case which is bar-stable. These results affirm the occurrence of radial heating by the bar force as spotted in past literatures, while it is not active without the bar formation.

Finally, we explore the possibility of the bar formation in supercritical-$Q$ disks. This inquiry aims to challenge recent observations in which barred galaxies have been spotted at redshift greater than $1$ (Guo et al., 2023). In that era, disk galaxies were typically gas-rich and turbulent. These indicated the hypothetically higher $Q$ than those in the local universe. The fact that bars were spotted in such galaxies implied that they were bar-unstable, even though their kinematical conditions were not in favor. Furthermore, other observational evidences and simulations revealed the barred structures that were able to emerge in thick galactic disks, which also implied that they were kinematically hot (Kasparova et al., 2016; Martig et al., 2021; Ghosh et al., 2023). From our past results, we have shown that there existed a critical $Q$ above which the circular disk was stable against the bi-symmetric perturbations in all scenarios. To resolve this puzzle, we perform further test to verify if the bar is able to be formed in supercritical $Q$ regime. We employ the rigid halo framework in which the stable disk is more quiet than the live halo counterpart. We choose $Q=2.25$ that is even above the highest $Q$ examined in Sec. 5.1, that was equal to $2.2$, and it was sufficient to stabilize the disk. In other words, the chosen disk kinematical condition here is hotter than the hottest and bar-stable disk in Sec. 5.1. This chosen $Q$ almost reaches the limiting value at which the Jeans equations yield the imaginary tangential velocity. We will make some modification of the disk center in the attempt to enhance the initial bi-symmetric perturbations to see if it can trigger the bar formation in such disk environment. To do so, we modify the positions of particles inside a core radius $R_{c}$ so that they are now confined in an ellipse of semi-major axis $R_{c}$ and flattening $f$. Finally, the velocity structure of the disk is assigned according to the new disk configuration. This modification introduces the large-scale bi-symmetric perturbations which are supposed to be stronger than those originated from the Poissonian noise alone. An example of the initial condition with modified core of $R_{c}=0.5$ and $f=0.4$ is depicted in the top panel of Fig. 13. To inspect the following evolution, time evolution of the bar strength is plotted in Fig. 14 for different $R_{c}$ and $f$. We capture the strong breakthrough bars in the two cases as $A_{2}$ is able to reach $0.4$ at the maximum and it remains high afterwards. One another case also indicates a weaker bar with $A_{2}$ almost reaching $0.2$. In those cases, the maximum $A_{2}$ clearly exceeds the value in the early stage which indicates the efficient bar formation. That breakthrough bar is further examined by inspecting visually the surface density map of the case in the top panel of Fig. 13 while in the bottom panel, it depicts that disk at $t=3.03$. We spot an evident barred structure residing in that disk.

In addition to our proposed mechanism of the bar formation in supercritical-$Q$ disks, the breakthrough bar was alternatively conjectured to be built by the close galactic encounter. That hypothesis has been initiated and tested by a number of simulations of the fly-by between a pair of galaxies and the bar formation in an initially bar stable disk has been reported, provided that the tidal interaction was sufficiently strong (Noguchi, 1987; Gerin et al., 1990; Gajda et al., 2017; Ghosh et al., 2021). In a larger frame, the analysis of evolution of interacting galaxies in the large-scale cosmological simulations also validated that possibility (Peschken & Łokas, 2019; Łokas, 2021). The fly-by scenario shares one similarity with our scenario: the hot disk can become unstable if it is triggered by sufficiently strong bi-symmetric perturbations of any origins. However, there is one major difference. The galactic fly-by scenario requires the perturbations from another encountering object. On the other hand, our scenario is the spontaneous scenario given an asymmetric seed. We do not rule out the possibility of this bar formation scenario as, according to some fly-by simulations (see, e.g., Lang et al. 2014; Łokas 2018), the emergence of unstable asymmetric seed was possible when the tidal interaction was not sufficiently strong to build up the tidal bar but it nevertheless left the disk center more eccentric than the initial state. It also hinted a sign of ongoing evolution. In some other cases with intermediate perturbations, the bar strength evolution exhibited a two-step pattern: it was boosted at first to a certain value by the encounter and remained there for a period of time. Afterwards, it evolved again to a higher value. We speculate that this was the secondary (or post-encounter) bar formation process that was self-induced by the asymmetric seed left by the past fly-by. This process is a distinct process from the direct tidal bar formation in which the tidal force forms the bar directly. This post-encounter scenario can potentially explain the origin of a number of observed hot barred galaxies aside the tidal bar hypothesis.