Simulation of Head-on Collisions Between Filamentary Molecular Clouds Threaded by a Lateral Magnetic Field and Subsequent Evolution

In this study, we investigated a head-on collision of two identical filaments penetrated by a magnetic field such that they share the same magnetic flux. To simulate the filament–filament collision, we solve the following ideal magnetohydrodynamic equations, which can be written in a conservative form as follows:

where $\rho$ is the density, $\bm{v}$ is velocity, $\bm{B}$ is the magnetic field, $P_{\rm tot}$ is the total pressure, which is the sum of the gas pressure and the magnetic pressure and is written as,

, and $\psi$ is the gravitational potential obtained by solving Poisson’s equation expressed in the following equation:

where $G$ is the gravitational constant. We assumed that the isothermal filaments are confined within an isothermal hot ambient medium, leading to the presence of two gases with different temperature. To reproduce this, we implemented a scalar field, in which the equation of state is given as,

where $S$ is the scalar field proportional to the temperature, and we assumed that $S$ evolves according to the advection of gas as follows:

Further details regarding the scalar field is described in the next section (see section 3.2).

To solve the basic equations, we used Athena++ (Stone et al., 2020) with the following setting: the time evolution of these magnetohydrodynamic equations was solved using the third-order Runge–Kutta time integrator scheme, which was proposed by Gottlieb et al. (2009), and the Riemann problem was solved using the HLLD method (Miyoshi & Kusano, 2005). For spatial construction, the piecewise linear method was applied to the characteristic variables. The constrained transport method was applied for maintaining the initial $\nabla\cdot\bm{B}=0$ condition (Evans & Hawley, 1988; Gardiner & Stone, 2008). The multi-grid algorithm was used to solve Poisson’s equation (Equation (9)) and was implemented in Athena++ by Tomida & Stone (2023).

In addition, we modified the multi-grid module of our two-dimensional calculations because it was originally designed for three-dimensional simulations. The Jeans criterion (Truelove et al., 1997) was used to avoid unexpected numerical fragmentation. We stopped the calculations when the relation between the Jeans length $L_{\rm J}=\sqrt{\pi c_{s}^{2}/(G\rho)}$ and the grid size $\Delta x$ breaks the condition of $4\Delta x\leq L_{\rm J}$.

We assume that the filaments are initially in magnetohydrostatic equilibrium threaded by the global magnetic field, which is perpendicular to their long axis. This magnetohydrostatic equilibrium solution is obtained by numerical calculations (Tomisaka, 2014; Kashiwagi & Tomisaka, 2021). We then consider the initial condition of two colliding filaments, which are formed through fragmentation from the sheet, in which they shared the same magnetic flux (Nagai et al., 1998). The filaments are assumed to extend infinitely and uniformly in the $z$-axis direction in the Cartesian coordinate system $(x,y,z)$. In the $x-y$ plane, we assume a computational domain of $-2.5R_{0}\leq x,y\leq 2.5R_{0}$, where $R_{0}$ is the radius of the hypothetical parent cloud of the filament in magnetohydrostatic equilibrium (see Equation (4)). The magnetic field outside the filament is considered to be force-free and converged to the uniform magnetic field of $B_{0}$ with an increase in the distance from the filament $r\gg R_{0}$, following Tomisaka (2014). Figure 1 shows one of the initial states assumed in this study. The two identical filaments are adjacent to each other and moving in $\pm x$-directions. The global magnetic field is parallel to the $x$-axis.

In our simulation, we assumed that the gas temperature in the molecular cloud is constant, $T\simeq 5-10\,\mathrm{K}$; thus, the equation of state is isothermal. However, we consider the possibility that the temperature of the gas may depend on its origin (filament or external medium). Therefore, we use the scalar field (Equation (10)) in the equation of state. Initially, we assume that the filament and the external medium achieved pressure balance at the surface as follows:

where $p_{s}$ is the pressure at the filament surface; $\rho_{s}$ and $\rho_{\mathrm{ext}}$ represent the density at the surface and that of the external medium, respectively; and $S_{\mathrm{fil}}$ and $S_{\mathrm{ext}}$ denote values of scalar $S$ inside and outside of the filament, respectively. We set the density contrast of the external medium and the filament surface as $\rho_{\mathrm{ext}}/\rho_{s}=0.01$. Therefore, the scalar $S$ of the external medium is $S_{\mathrm{ext}}=100c_{s}^{2}$ and that of inside the filament is $S_{\mathrm{fil}}=c_{s}^{2}$, where $c_{s}$ is the sound speed inside the filament. Finally, in our simulation, we maintain the value of the sound speed of the external medium at $S_{\mathrm{ext}}^{1/2}=10c_{s}$. In other words, as the density of the external medium is lower than that of the filament, we assume that the temperature of the external medium is higher than the filament, and the initial temperature is maintained after the simulations are started.

The boundary conditions of the hydrodynamic variables were imposed to be periodic in the $y$-direction and outflow in the $x$-direction. The boundary condition of the gravitational potential was assumed to be the Dirichlet condition, and the potential is calculated using multipole expansion.

In this study, physical variables are normalized using the following quantities: external pressure ($p_{\mathrm{ext}}$), density of the filament surface ($\rho_{s}$), and sound velocity in the filament ($c_{s}$). The scale length $L$ is defined by the free-fall time $t_{\mathrm{ff}}$ at the filament’s surface density $\rho_{s}$ and the isothermal sound speed inside the filament as $L=c_{s}t_{\mathrm{ff}}=c_{s}/(4\pi G\rho_{s})^{1/2}$. The physical scales characterizing the system are given in Table 1. We define the normalized variables as follows: $\bm{v}^{\prime}\equiv\bm{v}/c_{s}$, $p^{\prime}\equiv p_{\mathrm{gas}}/p_{\mathrm{ext}}$, $\rho^{\prime}\equiv\rho/\rho_{s}$, $\Phi^{\prime}\equiv\Phi/(B_{u}L)$, $\bm{B}^{\prime}\equiv\bm{B}/B_{u}$, $\lambda^{\prime}\equiv\lambda/(\rho_{s}L^{2})$, $x^{\prime}\equiv x/L$, and $y^{\prime}\equiv y/L$, where the prime represents the normalized variables. The scalar field is normalized using the isothermal sound speed in the filament as $S^{\prime}\equiv S/c_{s}^{2}$.

In this study, we consider three parameters: the field strength $B_{0}$, the total line mass $\lambda_{\mathrm{tot}}$, and the collision speed $V_{\mathrm{int}}$. The former two characterize the magnetohydrostatic state of each of the filaments, while $V_{\mathrm{int}}$ specifies the condition of the collision.

The nondimensional form of the field strength $B_{0}$ is given by $B^{\prime}_{0}=B_{0}/(4\pi p_{\mathrm{ext}})^{1/2}=\sqrt{2}\beta_{0}^{-1/2}$, where $\beta_{0}=8\pi p_{\mathrm{ext}}/B_{0}^{2}$ is the plasma beta of the external gas far from the filaments. In this study, $\beta_{0}$ is used instead of $B_{0}$ to specify the field strength. We examine three plasma beta values: $\beta_{0}=1$ (weak), $\beta_{0}=0.1$ (fiducial), and $\beta_{0}=0.01$ (strong).

The total line masses are determined based on the magnetic critical line mass, whose nondimensional form is given by

where the nondimensional magnetic flux $\Phi^{\prime}=R_{0}^{\prime}\sqrt{2/\beta_{0}}$ is a function of $R_{0}^{\prime}$ and $\beta_{0}$, that is, $\lambda_{\mathrm{tot}}^{\prime}/2\leq\lambda_{\mathrm{crit,B}}^{\prime}$. For simplicity, the radius of the parental cylinder $R_{0}^{\prime}$ is set to $2$. The magnetic critical line masses for the three $\beta_{0}$ are $\lambda^{\prime}_{\rm crit,B}(\beta_{0}=1)=1.17\lambda^{\prime}_{\rm crit}$, $\lambda^{\prime}_{\rm crit,B}(\beta_{0}=0.1)=1.91\lambda^{\prime}_{\rm crit}$, and $\lambda^{\prime}_{\rm crit,B}(\beta_{0}=0.01)=4.25\lambda^{\prime}_{\rm crit}$, where $\lambda^{\prime}_{\mathrm{crit}}$ is the nondimensional form of the critical line mass of Equation (2) as $\lambda^{\prime}_{\mathrm{crit}}=\lambda_{\mathrm{crit}}/(\rho_{s}L^{2})=8\pi$. We examine three types of line masses for the respective magnetic strengths, in which the filament in models L (massive), M (intermediate), and S (less massive) has a line mass of $0.9$, $0.75$, and $0.5$ times as large as the magnetically critical line mass, $\lambda_{\mathrm{crit,B}}^{\prime}(\beta_{0})$, respectively. We also calculated additional models of S’ (0.6$\lambda_{\mathrm{crit,B}}^{\prime}(\beta_{0})$) and M’ (0.7$\lambda_{\mathrm{crit,B}}^{\prime}(\beta_{0})$), if necessary.

For calculating the collision speed $V_{\mathrm{int}}$, the cases with $V_{\mathrm{int}}=c_{\mathrm{s}}$ and $10c_{\mathrm{s}}$ are considered.

In short, the head-on collision was reproduced using the plasma beta value ($\beta_{0}$), total line mass ($\lambda^{\prime}_{\mathrm{tot}}$), and initial relative velocity ($V^{\prime}_{\rm int}$) as the parameters. In this study, we simulate 16 different models, shown in Table 2. Hereafter, we omit ^′ with normalized variables, unless the quantities are misunderstood as dimensional variables.

In this section, we describe the typical evolution of the radial collapse model. The model with $\beta_{0}=0.1$, $\lambda_{\mathrm{tot}}=2.87\lambda_{\rm crit}$, and $V_{\rm int}=1$ is selected as the fiducial model and corresponded to $\beta 01$M in Table 2.

Figure 2 shows the images of the time evolution of the filament–filament collision. In Figure 2(a), at $t=0.350$, the collision of two filaments resulted in the formation of a sheet-like dense structure, which is elongated in the $y$-direction and with a thickness of $\sim 0.1$ measured on the $x$-axis at the center of the merged filament, referred to as the shocked region.

Figures 3(a) and (b) show the density and velocity profiles on the $x$- and $y$-axes, respectively. In the upper panel of Figure 3(a), a shock front is visible at $|x|\simeq 0.06$ in the density profile. This is visualized as the edge of the shocked region, which extends in the $y$-direction in Figure 2(a). Along the $y$-axis, shock fronts have not yet formed at $t=0.350$ (see the upper panel of Figure 3(b)). By comparing the density profiles on the $x$- and $y$-axes, we can observe that the density distribution of the merged filament in the $y$-direction ($|y|\simeq 1.5$) is wider than that in the $x$-direction ($|x|\simeq 0.06$). In the lower panel of Figure 3(a), the velocity distribution ($v_{x}$) at $t=0.350$ shows that a global inflow to the shocked region occurs in the $x$-direction outside the shocked region ($|x|\gtrsim 0.06$). On the $y$-axis, the velocity distribution ($v_{y}$) does not show such a global inflow; however, an accretion flow is newly formed by the self-gravity near $|y|\simeq 1$ (see the lower panel of Figure 3(b)).

Figure 2(b) shows that the shocked region becomes denser at $t=0.670$ and is observed more clearly with time. In addition, the magnetic field lines are dragged toward the center of the shocked region. Magnetic field lines bend at the shock front reaching the front in the post-shock (inner) region, indicating that the shock wave created by the collision is a fast shock.

In the upper panel of Figure 3(a), By comparing the density profiles at $t=0.350$ (red line) and $t=0.670$ (yellow line), We show that the overall density of the shocked region increases approximately by a factor of $100$.

The position of the shock front moved from $|x|\simeq 0.06$ to $|x|\simeq 0.13$, indicating that the shock wave is expanding. The shock front, which is at $|x|\simeq 0.06$ at $t=0.350$, expands to $|x|\simeq 0.13$ at $t=0.670$.

Furthermore, another shock is observed at $|x|\simeq 0.04$, which continues to expand to the edge of the filament $|y|\lesssim 1.5$, as observed in Figure 3(b). The panel shows that the magnetic field lines bend at the front, departing from the front in the post-shock region $|x|\lesssim 0.04$. In addition, the lower panel of Figure 3(c) shows that the plasma beta attains $\beta\simeq 0.13$ in the pre-shock region $|x|\gtrsim 0.1$, decreases with a jump passing the fast shock front at $|x|\simeq 0.13$, and increases with a positive jump at the accretion shock region $|x|\simeq 0.04$, which indicates that the second shock is a slow shock. In other words, as a consequence of the supersonic collision of the filaments, a pair of outward-facing shock fronts are formed (outer fast and inner slow shocks).

On the contrary, along the $y$-axis, although the overall density of the shocked region increases with time, there is no clear evidence of the presence of a shock wave (see the upper panel of Figure 3(b)). In the lower panel of Figure 3(b), the velocity distribution ($v_{y}$) at $t=0.670$ shows that the central directional velocity component ranges from $|y|\simeq 10^{-2}$ to $|y|\simeq 2$. The lower panel of Figure 3(b) shows that the gas in $10^{-2}\lesssim|y|\lesssim 2$ moves toward the center. By comparing the velocity structures at $t=0.350$ and $t=0.670$, we show that the inflow is accelerated with time by the effect of self-gravity.

At $t=0.710$, Figure 2(c) shows that the fast shock front, which was formed by the collision, expands to $|x|\simeq 0.2$. Another slow shock is almost stalled around $|x|\simeq 0.04$. However, other differences from the previous epoch are not easily distinguished from the two-dimentional density profiles. Therefore, to more effectively analyze the evolution after the collision, we primarily focused on the radial profiles. In the upper panel of Figures 3(a) and (b), the density profiles at $t=0.710$ (green lines) showed that the overall density of the shocked region continues to increase compared with the last epoch with $t=0.670$ (yellow). Although the slow shock facing toward the $x$-axis was almost stagnant near $|x|\simeq 0.04$, the density profile illustrated that the accretion shock is newly formed at $|y|\simeq 0.03$ along the $y$-axis. This accretion shock is formed by the infalling gas as it reaches the high-density region, and the jump in the density and pressure at this location indicates the formation of an accretion shock. This indicates that the gas is continuously accreted by the self-gravity of the shocked region, leading to an increase in the density and pressure of the region. The velocity profiles clearly show that the high-density shocked gas, even inside the accretion shock, contracts as a whole, primarily by the inflow in the $y$-direction (see the lower panel of Figures 3(a) and (b)).

At $t=0.735$, the 1D distribution of the merged filament near its final state is represented by the blue solid lines in Figures 3(a) and (b). The shock front contracts as a whole and reaches $|x|\simeq|y|\simeq 0.025$. Inside $r\lesssim 5\times 10^{-3}$, the so-called “Hubble-like inflow” was observed as $v_{x}/x\simeq v_{y}/y\simeq\mathrm{const}$, and a uniform density core is contractiong in this region.

Next, we focus on the evolution of the magnetic field driven by the collision. Figures 3(c) and (d) show the time evolution of the gas ($p_{\rm gas}$), magnetic ($p_{\mathrm{B}}$) pressures, and plasma beta ($\beta\equiv p_{\rm gas}/p_{\mathrm{B}}$) along the $x$- and $y$-axes, respectively. The magnetic pressure is defined as $p_{\rm B}\equiv B^{2}/2$ in our nondimensional variables. In the upper panels of Figures 3(c) and (d), we showed that the magnetic pressure within the shocked region is consistently smaller than the gas pressure throughout the evolution. As an example, at $t=0.735$, the plasma beta in the central part of the shocked region is $\beta=2.83$, meaning that, even in the presence of a magnetic field, thermal pressure dominates near the center of the shocked region.

The collision process leads to the amplification of the magnetic field strength, as observed by comparing the magnetic pressure of the central part at $t=0.350$ (red dashed line) and $t=0.735$ (blue dashed line) in the upper panel of Figures 3(c) and (d), which indicate that the magnetic pressure at the central part increases from $p_{\rm B}\simeq 30$ to $p_{\rm B}\simeq 2.6\times 10^{5}$, and the amplification reaches a factor $\sim 10^{4}$.

This suggests that the magnetic field strength is amplified by a factor of $\sim 10^{2}$ owing to the contraction induced by the collision.

In Figure 4, we plot the evolution of maximum density $\rho_{\mathrm{max}}$ for the collapse ($\beta$01M) and the stable models ($\beta$01S). Except for a special circumstance in the early phase of collision, maximum density is attained at the center of the merged filament as $\rho_{\mathrm{max}}=\rho_{c}\equiv\rho(0,0)$. Figure 4 clearly shows that the central density increases monotonically with time, and the contraction never stops in the collapse model.

For stable runs, we develop the model $\beta 01$S ($\beta_{0}=0.1$, $\lambda_{\mathrm{tot}}=1.91\lambda_{\mathrm{crit}}$, and $V_{\rm int}=1$).

The dynamical evolution of the stable models is similar to that of the radial collapse model during the stage at which the shock fronts are sweeping two colliding filaments. Figure 5 shows four images of the time evolution after the collision of two filaments. In Figure 5(a), at $t=0.900$, we can observe a sheet-like structure enclosed with two shock fronts facing outwardly, which is the same as that observed in the radial collapse model (Figure 2(a)). Figures 6(a) and (b) represent the density and velocity profiles along the $x$- and $y$-axes, respectively. The red lines correspond to the same epoch, $t=0.900$, as observed in Figure 5(a), which reveals that the shock wave is formed only along the $x$-axis ($|x|\simeq 0.1$). In the lower panel, velocity distributions (red lines) confirm this finding by showing that the global inflow to the shocked region is only observed along the $x$-axis.

At $t=1.300$, Figure 5(b) shows that the central part of the shocked region undergoes contraction and becomes denser, with the magnetic field lines also being dragged towards the center. In Figure 6, the yellow lines indicate the density and velocity profiles at $t=1.300$. By comparing with the density profile of the previous epoch ($t=0.900$), we show that the central density of the shocked region increased by a factor of $\simeq 6$. In addition, the density profile along the $y$-axis shows that a new accretion shock is formed at $|y|\simeq 0.3$. In the lower panel of Figure 6(b), the velocity distribution further supports this finding as $v_{y}$ indicates the presence of the accretion shock at $|y|\simeq 0.3$, resulting from the global inflow occurring outside the new accretion shock ($|y|\gtrsim 0.3$). Early phase evolution is qualitatively the same as the collapsing model; for example, the formation of the shocked region initially occurred because of the $x$-axis global inflow (stage of $t=0.900$) and subsequently because of the $y$-axis accretion flow (stage of $t=1.300$). However, in the stable model, the velocity profile on the $x$-axis demonstrated that the shock fronts observed in Figure 5(a) have swept across the merged filament, the boundary of which is observed at $|x|\simeq 0.50$, and are traveling in an external low-density medium (see the yellow line in the lower panel of Figure 6(a)).

In Figure 4, we plotted the evolution of maximum density $\rho_{\mathrm{max}}(t)$ of this stable model in blue. Figures 5(c) and (d) depict the structures when $\rho_{\mathrm{max}}$ takes a maximum ($t\simeq 1.9$) and a minimum ($t\simeq 3.5$), respectively. In Figure 6, a comparison between these density profiles at $t=1.9$ and $3.5$ shows that the density near the center decreases from $t=1.900$ to $t=3.500$, more than one order of magnitude; that is, by passing the maximum at $t=1.900$, the shocked region does not collapse but instead expands, as shown in Figure 5(d). Furthermore, the magnetic field lines are dragged outward. After the expansion of the shocked region, the central density of the shocked region begins to increase again around $t\simeq 3.7$ and subsequently decreases again around $t\simeq 5.0$, as shown in Figure 4. Therefore, this stable model does not exhibit a global collapse but the shocked region continues to oscillate after the collision.

We then focus on the structure of the shocked region of the stable model. Figure 7 shows the density profile after oscillating for two cycles ($t=6.500$), in which the density profiles are plotted on the $x$- (magenta solid line) and $y$-axes (magenta dashed line). In this figure, we also plot the magnetohydrostatic equilibrium solution of Tomisaka (2014), which has the same line mass and magnetic flux as the merged filament, represented in cyan. Figure 7 clearly indicates that the density profiles of the shocked region are quite similar to those of the magnetohydrostatic equilibrium state.

In the previous sections 4.1 and 4.2, we show the models with a plasma beta value of $\beta_{0}=0.1$. Here, we compared the models with plasma beta values of $\beta_{0}=1$ (models $\beta 1$S and $\beta 1$M) and $0.01$ (models $\beta 001$S and $\beta 001$M).

Figure 8 shows the maximum density evolution of the radial collapse (a) and stable (b) models. For the radial collapse model (a), the strength of the initial magnetic field does not affect the outcome of radial collapse. However, the stable model (b) maintains a stable state while undergoing oscillations, regardless of the initial plasma beta value. We assumed that the line masses of intermediate (M) and less-massive (S) filament are 0.75 and 0.5 times as large as the magnetically critical line mass of Equation (13), respectively. As the magnetically critical line mass increases with a decrease in plasma beta $\beta_{0}$, line masses of the respective models satisfy the following relations: $\lambda(\beta\mathrm{1M})<\lambda(\beta\mathrm{01M})<\lambda(\beta\mathrm{001M})$, and $\lambda(\beta\mathrm{1S})<\lambda(\beta\mathrm{01S})<\lambda(\beta\mathrm{001S})$. The line mass increases as the plasma beta value decreases, leading to a stronger gravitational force in the merged filament. Consequently, the radial collapse model shows a shorter timescale, and the stable model exhibits a shorter oscillation period. We defined a time for the maximum density reaching $\rho_{\mathrm{lim}}$ as $t_{\mathrm{lim}}$. In the radial collapse model, $t_{\mathrm{lim}}$ changes from $t_{\mathrm{lim}}=1.103$($\beta 1$M) to $t_{\mathrm{lim}}=0.443$($\beta 001$M). In the stable model, the oscillation period $T_{\mathrm{osi}}$ becomes shorter from $T_{\mathrm{osi}}\simeq 4.0$($\beta 1$S) to $T_{\mathrm{osi}}\simeq 1.8$($\beta 001$S).

Next, we study the density structure of both models. Figure 9 shows the comparison of the density profiles for the radial collapse models with varying plasma beta values ($\beta_{0}$). The profiles are shown for the final epochs, which have the same central density of $\rho_{c}=\rho_{\mathrm{lim}}\simeq 1.66\times 10^{6}$. In Figure 9, two analytic solutions are shown, that of the self-gravitating isothermal plane-parallel gas disk (solid line; $\rho=\rho_{c}\mathrm{sech}^{2}\left(r\sqrt{\rho_{c}/2}\right)$) and that of the self-gravitating isothermal cylinder (dashed line; Equation (1)). By considering the central high-density part $\rho\gtrsim 10^{4}$ in Figure 9(a), we realize that the density profiles on the $x$-axis are more similar to the hydrostatic solution of the gas disk (solid line) than the solution of the gas cylinder (dashed line), regardless of $\beta_{0}$. This is because the Lorentz force is weak along the $x$-axis, leading to the thermal pressure playing a dominant role against the self-gravity of the gas disk. On the contrary, in the central high-density part on the $y$-axis (b), the profiles seem to be more expanded than that of the isothermal gas cylinder as the plasma beta value decreased from $\beta_{0}=1$ (red) to 0.01 (yellow) owing to the strong Lorentz force.

Figure 10 shows the density distribution of models $\beta 1$S (a), $\beta 01$S (b), and $\beta$001S (c) achieved at the respective final epochs as $t=8.100$, 6.600, and 4.800. Figure 10 shows the comparison of these stable models with the magnetohydrostatic equilibrium states described in Tomisaka (2014). Of note, the equilibrium state is selected to have the same magnetic flux and the same central density as the final state of the corresponding stable model. The density structures in all the models resemble the equilibrium states; in particular, in the models with a strong magnetic field (b) and (c), two profiles are well-matched. Thus, in general, less massive filaments lead to a merged filament stable with oscillation, and its density profile is highly similar to that of the equilibrium state.

In this section, we describe the evaluation of the effect of the initial velocity on the evolution of filament collision. We compare four models: $\beta$01M, $\beta$01M_highV, $\beta$01S, and $\beta$01S_highV. In model $\beta$01M_highV, we assume the same filament (line mass and the magnetic flux) as model $\beta$01M but with a larger collision velocity, $V_{\mathrm{int}}=10$. Furthermore, model $\beta$01S_highV is a high-speed collision model of $\beta$01S. Figure 11 shows the evolution of the maximum density for these four models. Models $\beta$01M_highV and $\beta$01M exhibit runaway collapse, and models $\beta$01S_highV and $\beta$01S result in oscillation in $\rho_{\mathrm{max}}$.

Panel (a) shows the effect of the initial velocity on the time scale of reaching the radial collapse. We defined a time for the maximum density reaching $\rho_{\mathrm{lim}}$ as $t_{\mathrm{lim}}$. The time is approximately equal to $t_{\mathrm{lim}}=0.735$ ($V_{\mathrm{int}}=1$) and $0.600$ ($V_{\mathrm{int}}=10$), respectively. Thus, we show that the time scale required to reach radial collapse $t_{\mathrm{lim}}$ decreases as the initial velocity increases. This is because a higher initial velocity leads to a more efficient inflow of gas toward the shocked region. In other words, the time required for the shock wave to sweep the filament becomes shorter compared with a low-speed collision. However, if the initial velocity is too large (e.g., $V_{\mathrm{int}}\gg 10$), it is likely for the merged filament to expand more before undergoing radial collapse due to its momentum. We confirmed this finding using a model colliding with an extreme collision speed $V_{\mathrm{int}}=15$ (not shown in the paper). Such high-speed collisions result in a long time to collapse compared with the case with $V_{\mathrm{int}}=1$. By contrast, panel (b) shows that, even by considering the supersonic collision with a high Mach number ($V_{\mathrm{int}}=10$), less massive filaments of $\beta 01$S_highV model evolve into a stable state and form a static filament.

Next, we focus on the density structure of the radial collapse and stable models with a high Mach number ($V_{\mathrm{int}}=10$). Figure 12 shows the density profile on the $x$- and $y$-axes of the intermediate-mass model, $\beta 01$M_highV. By comparing the models with $V_{\rm int}=1$ (red line) and that with $V_{\rm int}=10$ (blue line), we observed that, despite the high initial velocity, the filament width of the shocked region remains similar to that of the result for $V_{\rm int}=1$. In other words, the structure of the dense part ($\rho\gtrsim 10^{4}$) of the shocked region remains almost similar, regardless of the initial velocity.

In addition, Figure 13 shows the density profile on the $x$- and $y$-axes of the less-massive model ($\beta 01$S_highV). The structure of the shocked region overlaps with the magnetohydrostatic equilibrium state, similar to that of $V_{\rm int}=1$ (see Figure 7).

In this study, we attempt to determine the conditions under which the merged filament becomes unstable in the radial direction by filament collision. Figure 14 summarizes the results of the collision, indicating whether the merged filament underwent radial collapse or not. This shows that models with larger line masses and/or models with a weaker magnetic field (larger $\beta_{0}$) are likely to experience radial collapse. The outcome does not appear to depend on the collision velocity, $V_{\mathrm{int}}$. In Figure 14, we plotted the critical line mass of magnetized filaments, $\lambda_{\mathrm{tot}}=\lambda_{\mathrm{crit,B}}(\beta_{0})$. The results can be separated by comparison with the magnetically critical line mass (Equation (13)). Merged filaments with $\lambda_{\mathrm{tot}}\gtrsim\lambda_{\mathrm{crit,B}}$ undergo radial collapse, whereas those with $\lambda_{\mathrm{tot}}\lesssim\lambda_{\mathrm{crit,B}}$ exhibit stable oscillation. From this finding, we conclude that the necessary condition for the occurrence of radial collapse in the filament collision is that the total line mass must exceed the magnetically supported critical line mass. In other words, the merged filament should be a magnetically supercritical filament. Even if the magnetic field is extremely strong, the magnetically supercritical filament is likely to undergo radial collapse. If the collision velocity is too fast, as explained in Section 4.4, collapse may be delayed. Therefore, the collapse appears not to occur within a finite time for the collision with an ultimately large initial velocity. However, in the case of a filament collision originated from typical velocity dispersion ($\Delta v\lesssim\mathrm{a\,few\,km\,s^{-1}}$) in molecular clouds, as considering in this study, this delayed collapse is not realized.

It is evident that the critical line mass is a dominant factor in determining the radial instability of the shocked region. This conclusion may be restricted to a specific scenario of a head-on collision described here. In such a collision, a merged filamentary structure elongated in the $z$ direction is retained as the shocked region because there is no inclination between the long axes of two filaments. As a result, we expect the radial instability of the shocked region to be described by the magnetically critical line mass (Equation (3)). This finding can help in considering the evolution of a filament–filament collision in a more general configuration. An elongated merger of a filament–filament collision is made even by a collision with some inclinations, although the filaments far from the collision point pass by freely. When such colliding filaments share the same magnetic flux, the necessary condition of the gravitational collapse given as the merged filament (or a hub) is magnetically supercritical.

In this section, we describe the examination of self-similar solutions of a collapsing filament. This solution has been previously studied to describe the properties of the collapsing filament in (Kawachi & Hanawa, 1998). In Kawachi & Hanawa (1998), for similarity coordinates, the defined zooming coordinate $\bm{\chi}$ and the density in the zooming coordinate $\varrho$ as,

where $t_{0}$ represents the time when the central density becomes infinity.³³3 To determine $t_{0}$, we first examine the time evolution of $1/\sqrt{\rho_{\mathrm{max}}}$. We identify the part that follows a linear relation similar to$1/\sqrt{\rho_{\mathrm{max}}}=at+b$. The time at which $1/\sqrt{\rho_{\mathrm{max}}}=0$ corresponds to $t_{0}=-b/a$.

Figure 15 shows the time evolution of the density profiles written in the similarity coordinate $(\bm{\chi},\varrho(\bm{\bm{\chi}}))$ for model $\beta$1M. The color gradation indicates a time sequence as the central density of the merged filament increases from $\rho_{c}=1460$ (cyan) to $\rho_{c}=71484$ (magenta). Panel (a) shows the density profiles on the $x$-axis. Although $\rho_{c}$ varies by a factor of $\sim 50$, the density profiles appear not to depend on time significantly when the similarity coordinates are used. Panel (b) shows the density profiles on the $y$-axis. Compared with the density profiles on the $x$-axis, the density profiles on the $y$-axis demonstrated an increase in their width as the central density increased. Thus, although the convergence to the self-similarity depends on the direction, we expect that the contraction proceeds in a self-similar manner, especially for the collision direction.

In this section, we describe the examination of the initial line mass dependence of the shocked region. Here, we assume that the line mass of the shocked region $\lambda_{\rho_{\mathrm{max}}/100}$ is estimated by integrating the density larger than $\rho_{c}/100$ as $\lambda_{\rho_{\mathrm{max}}/100}\equiv\iint_{\rho\geq\rho_{\mathrm{max}}/100}\rho dxdy$. This assumption aims to follow the evolution of the dense part in the shocked region. We compare the models with the same plasma beta value $\beta_{0}$ but different initial line masses $\lambda_{\mathrm{tot}}$ for three pairs of models: ($\beta 1$M, $\beta 1$L), ($\beta 01$M, $\beta 01$L), and ($\beta 001$M, $\beta 001$L).

In Figure 16, the line mass of the shocked region $\lambda_{\rho_{\mathrm{max}}/100}$ is plotted against the maximum density $\rho_{\mathrm{max}}$. As the maximum density rarely decreases with time, $\rho_{\mathrm{max}}$ substitutes for the elapsed time $t$. The initial relative velocity is set to $V_{\rm int}=1$. Although the line mass of the dense part $\lambda_{\rho_{\mathrm{max}}/100}$ differs initially (left end of each line), both line masses $\lambda_{\rho_{c}/100}$ come close to almost the same value during the time evolution for all pairs of different plasma beta values (right end of each line). In particular, in the models with $\beta_{0}=1$ (red line), the line mass maintains a constant value $\lambda_{\rho_{\mathrm{max}}/100}\simeq\lambda_{\mathrm{crit}}$ after the maximum density exceeds $\rho_{\mathrm{max}}\gtrsim 10^{5}$, regardless of the initial line mass.

In the model with a small plasma beta value of $\beta_{0}=0.01$, the final $\lambda_{\rho_{\mathrm{max}}/100}$ of the two models agrees with each other. However, it is unclear whether $\lambda_{\rho_{\mathrm{max}}/100}$ maintains the final value for further contraction. Evolution may reach only the extremely early stage with small beta models; thus, we require further evaluation to confirm the true convergence in the models with $\beta_{0}=0.01$. In conclusion, although the degree of convergence depends on the magnetic field strength, the typical line mass of the high-density portion of the contracting filament $\lambda_{\rho_{\mathrm{max}}/100}$ does not depend on the initial line mass of the filament.

Figure 17 shows a comparison of the density profiles realized in the final state. The density distribution is shown to be normalized by the maximum density at the final state such that the vertical axis represents $\rho/\rho_{\mathrm{max}}$. In the horizontal axis, we consider the distance to be normalized by the scale-length given at the center as $x^{\prime}=x/\tilde{L}$, where $\tilde{L}=c_{s}/(4\pi G\rho_{\mathrm{max}})^{1/2}$. Figure 17(a) shows that the dense part of $\rho/\rho_{\mathrm{max}}\gtrsim 10^{-2}$ has an extremely similar density distribution along the $x$-axis, which is shown in Figures 9(a) and 12(a). By contrast, density distribution in the $y$-direction has a relatively poor similarity. Models $\beta$1M and $\beta$1L have similar distributions in the $y$-direction. This is valid for models $\beta$01M and $\beta$01L. However, model $\beta$001M has a wider distribution than $\beta$001L. This implies that the time evolution is not sufficient for models with $\beta_{0}=0.01$ compared with the other two models with $\beta_{0}=1$ and 0.1 (Figure 16).

The convergence of $\lambda_{\rho_{\mathrm{max}}/100}$ (Figure 16) and the density distributions resembling each other (Figure 17) indicate that the collapse is expressed using the same solution, independent of the initial line mass after $\rho_{c}$ increases sufficiently. Figure 16 shows that the difference in $\lambda_{\rho_{\mathrm{max}}/100}$ owing to different $\beta_{0}$ also decreases as the collapse proceeds. However, it is not clear whether the collapse is also expressed with a unique solution, irrespective of $\beta_{0}$, similar to the axisymmetric three-dimensional cloud (Nakamura et al., 1999).

In this section, we describe the investigation of the fragmentation of merged filaments resulting from head-on collisions. Although our investigation primarily focused on the radial instability caused by such collisions, the fragmentation process is essential for the formation of molecular cloud cores from the filaments. Studies on the fragmentation of magnetized filaments have been conducted by Stodółkiewicz (1963); Nagasawa (1987); Hanawa et al. (1993); Nakamura et al. (1993); Fiege & Pudritz (2000), with the assumption of a magnetic field that is either parallel or helical to the filament axis. By contrast, Hanawa et al. (2017, 2019) assumed a magnetic field that is perpendicular to the filament axis.

For the nonmagnetized filament, the maximum growth rate of the gravitational instability was found to be $|\omega_{\mathrm{max}}|=0.339(4\pi G\rho_{c})^{1/2}$ by Nagasawa (1987), resulting in a timescale of the maximum growth as $t_{\mathrm{min}}\equiv 1/|\omega_{\mathrm{max}}|=2.94(4\pi G\rho_{c})^{-1/2}$. If the dynamical timescale is slower than the timescale of the gravitational instability and if this continues for a sufficient duration, fragmentation occurs in the filament. For a filament penetrated with a uniform magnetic field perpendicular to the axis, ⁴⁴4Notably, although the model of a uniform magnetic field does not influence the cloud’s radial stability, it is not in the magnetohydrostatic equilibrium state. the wave number ($k$) and growth time showing the maximum growth are expressed approximately as $k\simeq 0.15\sqrt{4\pi G\rho_{c}}/c_{s}$, and $t_{\mathrm{min}}\simeq 5.78(4\pi G\rho_{c})^{-1/2}$ according to Hanawa et al. (2017). As the Lorentz force suppresses the fragmentation process, the growth rate becomes small, and the wavelength of the fragmentation becomes longer compared with the nonmagnetized ones.

We assume that density is perturbed along the $z$-direction and investigate whether it can grow and form fragmentations in the merged filament. For the radial collapse models, the central (maximum) density evolves almost monotonically, making fragmentation difficult (see Figures 4 and 8). However, we find that the contraction stops temporally for a certain period, even in the radial collapse models such as $\beta$001M in Figure 8(a) and $\beta$01M_highV in Figure 11(a). Here, we studied the possibility of whether gravitational instability developing during contraction is delayed. For instance, model $\beta 01$M_highV expands after the collision at approximately $t=0.3$ (see Figure 11). We found that the central density is maintained at $\rho^{\prime}_{c}\simeq 250$ during $t^{\prime}_{\mathrm{stall}}\simeq 0.08$, and the growth time of gravitational instability $t^{\prime}_{\mathrm{min}}\simeq 0.36(\rho^{\prime}_{c}/250)^{-1/2}$ is longer than the lag time of contraction ($t_{\mathrm{stall}}<t_{\mathrm{min}}$), and, even in the phase of maxima at $t\sim 0.16$, $t_{\mathrm{stall}}$ is shorter than $t_{\mathrm{min}}$. Another radial collapse model ($\beta 001M$) with a bump is also $t_{\mathrm{stall}}<t_{\mathrm{min}}$ for the entire range of the simulation.

contrastingly, the merged filaments in the stable models can undergo fragmentation around the maxima or minima during oscillation, regardless of the plasma beta value. For example, in the model of $\beta 01$S, the maximum density is maintained at $\rho^{\prime}_{c}\simeq 226$ for $t^{\prime}_{\mathrm{stall}}\simeq 0.7$ (Figure 8(b)); thus, $t^{\prime}_{\mathrm{stall}}>t^{\prime}_{\mathrm{min}}\simeq 0.38(\rho^{\prime}_{c}/226)^{-1/2}$.

In this section, we investigate the case where the collision direction of the filamentary molecular clouds is perpendicular to the direction of magnetic field lines. We studied the models in which the magnetic field lines ran parallel to the colliding direction. In reality, this happens in a case where the magnetic field lines run perpendicular to the colliding direction.

Figure 18 shows the time evolution of models $\beta 1$M’$\perp$ (upper) and $\beta 1$S’$\perp$ (lower), in which the total line mass is considered as $\lambda_{\mathrm{tot}}=1.64\lambda_{\mathrm{crit}}$ and $1.40\lambda_{\mathrm{crit}}$. The other parameters are set as $\beta_{0}=1$ and $V_{\rm int}=1$.

In the model of $\beta$1M’$\perp$, the shocked region begins to collapse, whereas, in the model of $\beta$1S’$\perp$, the shocked region does not collapse and is stable while oscillating. As the filament with $\beta_{0}=1$ has a magnetically critical line mass of $\lambda_{\mathrm{crit,B}}=1.17\lambda_{\mathrm{crit}}$ (Equation (13)), the total line mass of the two models exceeds the critical line mass as $1.4$ ($\beta$1M’$\perp$) and $1.2$ ($\beta$1S’$\perp$) times $\lambda_{\mathrm{crit,B}}$. Thus, if a collision occurs in the parallel direction to the magnetic field, both models would undergo the radial collapse (sec. 5.1) as $\lambda_{\mathrm{tot}}>\lambda_{\mathrm{crit,B}}$. However, model $\beta 1$S’$\perp$ results in a stable oscillation (Figure 18 lower panels).

This difference is attributed to the collision direction. In this case, the magnetic flux of the merged filament become twice as large as that in the previous models ($\beta 1$S and $\beta 1$M), in which the colliding direction is parallel to the magnetic field lines. Thus, the expression of the threshold line mass should be modified by the magnetic field’s geometry. When we consider that the magnetic field lines run perpendicular to the collision direction, the magnetically critical line mass (Equation (3)) of the merged filament is changed as

We rewrite the above equation in the normalized form as,

By adopting this equation, the critical line mass of a merged filament is estimated as, $\lambda_{\mathrm{crit,B}}=1.51\lambda_{\mathrm{crit}}$ instead of $1.17\lambda_{\mathrm{crit}}$. Therefore, the total line mass of $\beta$1M’$\perp$ model $\lambda_{\mathrm{tot}}=1.64\lambda_{\mathrm{crit}}$ exceeds the critical value $\lambda_{\mathrm{crit,B}}=1.51\lambda_{\mathrm{crit}}$ and the collapse of the shocked region, even in the perpendicular collision. On the contrary, for the model $\beta$1S’$\perp$, the total line mass $\lambda_{\mathrm{tot}}=1.17$ is significantly smaller than the critical line mass $1.51\lambda_{\mathrm{crit}}$. Thus, the merged filament of model $\beta$1S’$\perp$ is magnetically subcritical $\lambda_{\mathrm{tot}}<\lambda_{\mathrm{crit,B}}$ and exhibits stable oscillations. Thus, we conclude that the stability of the perpendicular model can be explained by considering whether the merged filament is larger than the magnetically critical line mass if we consider the fact that the magnetic flux increases twice due to the collision.