Magnetohydrodynamic Models of Molecular Tornadoes

TL;DR

This paper develops magnetohydrodynamic models to understand the structure of molecular tornadoes near the Galactic Centre, highlighting the roles of magnetic fields and external pressure in their equilibrium.

Contribution

It introduces an analytical MHD model for rotating, helically magnetized filaments and links magnetic properties to observed tornado structures, improving understanding of their equilibrium states.

Findings

Magnetic fields and external pressure dominate tornado confinement.

Analytical solutions avoid unphysical density inversions.

Models match observed properties of specific molecular tornadoes.

Abstract

Recent observations near the Galactic Centre have found several molecular filaments displaying striking helically-wound morphology, which are collectively known as "molecular tornadoes." We investigate the equilibrium structure of these molecular tornadoes by formulating a magnetohydrodynamic model of a rotating, helically magnetized filament. A special analytical solution is derived where centrifugal forces balance exactly with toroidal magnetic stress. From the physics of torsional Alfv\'{e}n waves, we derive a constraint that links the toroidal flux-to-mass ratio and the pitch angle of the helical field to the rotation laws, which we find to be an important component in describing molecular tornado structure. The models are compared to the Ostriker solution for isothermal, non-magnetic, non-rotating filaments. We find that neither the analytic model nor the Alfv\'{e}n wave model suffer from unphysical density inversions noted by other authors. A Monte Carlo exploration of our parameter space is constrained by observational measurements of the Pigtail Molecular Cloud (Pigtail), Double Helix Nebula (DHN), and Galactic Centre molecular Tornado (GCT). Observable properties such as the velocity dispersion, filament radius, linear mass, and surface pressure can be used to derive three dimensionless constraints for our dimensionless models of these three objects. A virial analysis of these constrained models is studied for these three molecular tornadoes. We find that self-gravity is relatively unimportant, whereas magnetic fields, and external pressure play a dominant role in the confinement and equilibrium radial structure of these objects.

Equations

$$\nabla \cdot (\rho \mathbf{v})=0,$$

$$\rho (\mathbf{v}\cdot \nabla )\mathbf{v}=-\nabla P-\rho \nabla \Phi +\frac{(\nabla \times \mathbf{B})\times \mathbf{B}}{4\pi },$$

$$\nabla \times (\mathbf{v}\times \mathbf{B})=0,\text{ and}$$

$$P={\sigma }^2\rho ,$$

$$\rho ={\rho }_0\tilde{\rho },$$

$$r=\frac{\sigma }{\sqrt{4\pi G{\rho }_0}}\tilde{r},$$

$$P={\sigma }^2{\rho }_0\tilde{P},$$

$$B=\sigma \sqrt{{\rho }_0}\tilde{B},$$

$$\Phi ={\sigma }^2\tilde{\Phi },$$

$$g=\sqrt{4\pi G{\rho }_0}\sigma \tilde{g},$$

$$m=\frac{{\sigma }^2}{4\pi G}\tilde{m},$$

$$v_{\varphi }=\sigma {\tilde{v}}_{\varphi },$$

$$\Omega ={\Omega }_0\tilde{\Omega }=\sqrt{4\pi G{\rho }_0}\tilde{\Omega }.$$

$$\mathbf{v}=v_{\varphi }(r)\hat{\varphi },$$

$$\mathbf{B}=B_{\varphi }(r)\hat{\varphi }+B_z(r)\hat{z}.$$

$$0=-\frac{dP}{dr}-\rho\frac{d\Phi}{dr}+\rho\frac{v_{\phi}^{2}}{r}-\frac{B_{\phi}^{2}}{4\pi r}-\frac{\partial}{\partial r}\bigg{(}\frac{B_{\phi}^{2}+B_{z}^{2}}{8\pi}\bigg{)}.$$

$${\Gamma }_{\varphi }=\frac{B_{\varphi }}{r\rho },$$

$${\Gamma }_z=\frac{B_z}{\rho },$$

$$\frac{dg}{dr}=-\frac{g}{r}-\rho .$$

$$\frac{d\rho }{dr}=\frac{\rho g+\rho {\Omega }^{2}r-\frac{1}{2\pi }r{\rho }^2{\Gamma }_{\varphi }^2-\frac{1}{8\pi }{r}^2{\rho }^2\frac{\partial }{\partial r}{\Gamma }_{\varphi }^2-\frac{1}{8\pi }{\rho }^2\frac{\partial }{\partial r}{\Gamma }_z^2}{1+\frac{1}{4\pi }{r}^2\rho {\Gamma }_{\varphi }^2+\frac{1}{4\pi }\rho {\Gamma }_z^2}.$$

$$P_{tot}=P+\frac{B_{\varphi }^2+B_z^2}{8\pi }.$$

$$\rho v_{\varphi }^2-\frac{B_{\varphi }^2}{4\pi }=0.$$

$$\beta=P\bigg{(}\frac{B_{\phi}^{2}}{8\pi}+\frac{B_{z}^{2}}{8\pi}\bigg{)}^{-1}=constant.$$

$${\beta }_z:=\frac{P}{P_{mag,z}},$$

$${\beta }_{\varphi }:=\frac{P}{P_{mag,\varphi }},$$

$$\beta=\bigg{(}\frac{1}{\beta_{z}}+\frac{1}{\beta_{\phi}}\bigg{)}^{-1}.$$

$$v_{\varphi }=\sqrt{\frac{2{\sigma }^2}{{\beta }_{\varphi }}}=constant.$$

$$P_{tot}=P+P_{mag,z}+P_{mag,\phi}=P\bigg{(}1+\frac{1}{\beta_{z}}+\frac{1}{\beta_{\phi}}\bigg{)}.$$

$$\sigma_{eff}^{2}:=\bigg{(}1+\frac{1}{\beta_{z}}+\frac{1}{\beta_{\phi}}\bigg{)}\sigma^{2},$$

$$\tilde{\Phi }:=\frac{\Phi }{{\sigma }_{eff}^2},$$

$$0=\frac{\partial }{\partial r}\log \rho +\frac{\partial }{\partial r}\Phi ,$$

$${\nabla }^2\Phi =e^{-\Phi },$$

$$r_0=\frac{{\sigma }_{eff}}{\sqrt{4\pi G}{\rho }_0}.$$

$$\rho=\frac{\rho_{0}}{\Big{(}1+\frac{r^{2}}{8r_{0}^{2}}\Big{)}^{2}},$$

$$v_{\varphi }=\sqrt{\frac{2{\sigma }^2}{{\beta }_{\varphi }}},$$

$$r_0^2=\frac{{\sigma }^2(1+{\beta }_z^{-1}+{\beta }_{\varphi }^{-1})}{4\pi G{\rho }_0}.$$

$$\Omega=\Omega_{0}\bigg{(}\frac{r}{r_{0}}\bigg{)}^{\alpha}.$$

$$\Omega =\frac{v_{\varphi }}{r}=\frac{\sqrt{2Gm(r)}}{r}.$$

Full text

ISM: clouds — magnetohydrodynamics (MHD) — methods: numerical

††journal: ApJ

\AuthorCollaborationLimit

=3

1 Introduction

The Pigtail Molecular Cloud (Pigtail) is an intriguing molecular filament that resides close to the Galactic midplane within the central molecular zone (CMZ). It exhibits a helical structure that was first noticed by Matsumura et al. (2012) from CO $J=1-0$ data observed by the Nobeyama Radio Observatory 45m radio telescope. The Pigtail helix has a $\sim\,15$ pc pitch that traces at least two rounds of a helix in 12CO, extending to an angular size of $\sim\,0.15^{\circ}\times 0.15^{\circ}$ and corresponding to a spatial size of $\sim\,20\times 20$ pc2 with a mass of $(2-6)\times 10^{5}\,M_{\odot}$ (Matsumura et al., 2012).

The proposed mechanism that triggers the helical morphology of the Pigtail is similar to that of two other molecular tornadoes known as the “Double Helix Nebula” (DHN) and the “Galactic Center molecular Tornado” (GCT) (Matsumura et al., 2012). A vertical magnetic tube that extends into the Galactic halo is twisted into its helical form. The driver of this twisting motion is hypothesized to be due to the shearing of inner ($x_{1}$) and outer ($x_{2}$) cloud orbits. In this scenario, the vertical magnetic tube is caught between two cloud orbits with different velocities, and the velocity shear twists the magnetic tube. With enough twisting, this torsional Alfvén wave causes the magnetic tube to undergo a kink ($m=1$) instability (Matsumura et al., 2012; Jackson, 1975), which winds it into its observed corkscrew-like morphology. Matsumura et al. (2012) estimate the Pigtail’s magnetic field to be on the order of 1 mG. Matsumura et al. (2012) support this formation scenario due to observations of the SiO/13CO ratio that are indicative of the shock between the cloud orbit interaction and of other molecular gas.

Another molecular tornado was found approximately $100$ pc from the Galactic Centre (GC) in the infrared by the Spitzer Space Telescope by Morris et al. (2006). They called it the “Double Helix Nebula” (DHN) because of its apparently intertwined double helix structure. The DHN is observed to wind at least $1.25$ complete rounds that stretch approximately 25 pc in length with its long axis parallel the Galaxy’s rotation axis. Morris et al. (2006) propose that the DHN is a magnetohydrodynamic (MHD) torsional Alfvén wave propagating away from the Galactic disk, driven by the circumnuclear disk (CND). Morris et al. (2006) estimated the magnetic field strength to be $B=0.1-1$ mG, depending on the assumed proton density in the region. Unlike the Pigtail Molecular Cloud, the presence of two intertwining strands of the DHN suggest that the twisted field due to the torsional wave has triggered an $m=2$ instability (Morris et al., 2006; Jackson, 1975).

Sofue (2007) first reported on the GCT, which was also observed in CO with the Nobeyama Radio Observatory. The GCT extends $170$ pc vertically from the Galactic plane and was measured to rotate at a velocity of $\sim\,30$ km s*-1* with an estimated mass of $M\,\sim\,1.2\times 10^{6}\,M_{\odot}$. Sofue (2007) note that the GCT, or at least part of it, consists of two helical strands like the DHN. Interestingly, the authors suggest that the GCT is not gravitationally bound, but that there must be an external force or pressure to keep the filament bound. They suggest that the most plausible mechanism confining the GCT against centrifugal force is the magnetic tension due to a toroidally wrapped field with $B\,\sim\,0.4$ mG. The GCT was proposed to arise from the same mechanism as that of the DHN – that a torsional Alfvén wave, driven by the epicyclic rotation of a cloud orbiting the Galactic Centre, twists a magnetic flux tube, which then undergoes an $m=2$ MHD instability (Sofue, 2007).

There has not yet been a detailed theoretical/numerical study on the structure of molecular tornadoes. In this paper, we consider the equilibrium structure of uniformly and differentially rotating, isothermal, self-gravitating and non-self-gravitating, pressure truncated cylinders threaded by somewhat general helical magnetic fields, analogous to the stationary filaments studied by Fiege & Pudritz (2000a, b). In our equilibrium models, we also take into account the propagation of torsional Alfvén waves. A simple analytical model is developed, followed by more general analytical models, and the latter are constrained using the limited observational constraints that are available (Matsumura et al., 2012; Morris et al., 2006; Sofue, 2007). Our models provide the equilibrium structure needed to study the instabilities described by Matsumura et al. (2012), Sofue (2007), and Morris et al. (2006) that are believed to trigger the coiling of a twisted magnetic tube, which we will explore in a future paper. An especially interesting question that we will address in future work is that of why some molecular tornadoes appear to exhibit an $m=1$ mode of instability, while others appear to be dominated by the $m=2$ mode. Additionally, the equilibrium models developed here, and in our forthcoming stability calculations, will enable us to predict sub-millimetre polarization maps, to hopefully further constrain our models in the future, when such maps become available.

Hansen et al. (1976) provided a theoretical study of the equilibrium and stability of uniformly rotating, isothermal, and infinitely long gas cylinders of finite radius. They noted an interesting feature, in which the density fluctuates non-monotonically as a function of radius. These fluctuations can be seen in Figure 1 of Hansen et al. (1976), where they were referred to as “density inversions.” Density inversions were also found by Recchi et al. (2014), who studied the equilibrium of uniformly and differentially rotating, non-magnetic, pressure truncated filaments. Hansen et al. (1976) suggested that density inversions are due to a battle between gravitational, centrifugal, and pressure gradient forces in maintaining the filament in equilibrium, and Recchi et al. (2014) proposed that such density profiles must be truncated by an external pressure before any density inversion occurs, because they would otherwise not be physically representative of observed filaments. We adopt this viewpoint in our helically magnetized models, which we truncate whenever the density is seen to increase with radius.

Kaur et al. (2006) presented their findings on differentially rotating, self-gravitating filaments, and have incorporated the effects of a helical magnetic field. They took into account the isothermal and logatropic equations of state, combined with constant velocity, and constant angular momentum rotation laws. The magnetic field in their model is proportional to the square root of the density of the filament. We also examine isothermal models, but with rotation laws that follow a power law in radius, as well as rotation laws that follow a power law in magnetic flux. Physics of torsional Alfvén waves are taken into account in our models, which has not been incorporated in the models of Kaur et al. (2006). The assumption of a torsional Alfvén wave couples the toroidal magnetic field $B_{\phi}(r)$ to the rotation law $\Omega(r)$ (see Section 2.2.4). Thus, the physics of torsional Alfvén waves provides a realistic constraint on the model, and results in interesting insights into the physics of molecular tornadoes.

Our model of molecular tornadoes begins with the equations of steady state MHD (Section 2). A special analytical solution is presented in Section 2.1. This closed form solution is not general, but it is interesting and does provide some insight into the physics of molecular tornadoes. Various differential rotation laws are described in Section 2.2 including a constraint due to torsional Alfvén waves. We conduct a virial analysis introduced in Section 2.3. Our numerical method is presented in Section 2.5, and the method for constraining our models by observed quantities is included in Section 2.6. Results are shown in Section 3, discussed in Section 4, and summarized in Section 5.

2 Theory

We assume ideal MHD, as described by the four ideal MHD equations in steady state – the continuity equation, Cauchy momentum equation, Faraday’s Law, and the isothermal equation of state, respectively given by

[TABLE]

where $\rho$ is the density, $\bm{v}$ is the velocity, $P$ is the pressure, $\bm{B}$ is the magnetic field, $\Phi$ is the gravitational potential, and $\sigma$ is the velocity dispersion. We consider the velocity dispersion to contain both thermal and nonthermal contributions: $\sigma=\sqrt{\sigma_{thermal}^{2}+\sigma_{nonthermal}^{2}}$.

For convenience, we define a set of dimensionless quantities, written here with an overscript tilde, where the physical quantities are scaled according to the following scaling laws:

[TABLE]

In these scaling relations, $G$ is the universal gravitational constant, $g$ is the gravitational field, $m$ is the linear mass, $v_{\phi}$ is the rotational speed, and $\Omega$ is the rotational frequency. The subscript 0 denotes the quantities at the filament’s radial core ($r=0$). It should be noted that all physical quantities, except $\sigma$, are considered functions of $r$. Hereafter, the overscript tilde notation is dropped, except where otherwise noted, since model quantities are presented in dimensionless form for most of the paper.

By assuming cylindrical symmetry ($r,\phi,z$) and rotational motion only, the velocity and magnetic field take the form

[TABLE]

Using the assumption of equation 3, the only non-trivial terms of equation 1 are in the $\hat{r}$ direction. Both the continuity equation and Faraday’s law are satisfied automatically under the assumption of equation 3, while the Cauchy momentum Equation simplifies to

[TABLE]

We introduce the magnetic flux-to-mass ratios formulated by Fiege & Pudritz (2000a). The toroidal and poloidal magnetic flux-to-mass ratios are defined by

[TABLE]

respectively. The flux-to-mass ratios are free functions of the magnetic flux $\Phi_{M}$ in the context of perfect MHD due to Ferraro’s law of isorotation (Ferraro, 1937), but flux $\Phi_{M}(r)$ is a function of $r$ in cylindrical symmetry. Fiege & Pudritz (2000a) assumed that both $\Gamma_{z}$ and $\Gamma_{\phi}$ were constant. However, we consider the possibility that these functions might vary as a function of $\Phi_{M}$, or equivalently $r$.

In addition to the equations of MHD, the effects of gravity also need to be taken into account. The gravitational field, $g(r):=-\partial\Phi(r)/\partial r$, is described by Poisson’s equation for gravity , which becomes

[TABLE]

Having defined $g$, equation 4 can be expressed as

[TABLE]

It is not generally possible to solve this equation analytically except under special assumptions. For most of the work presented in the paper, we apply numerical methods to solve equation 7, taking into account the physics of torsional Alfvén waves (see Section 2.2.4), which leads to the numerical solutions presented in Section 3. However, we first present one special solution that can be solved analytically.

2.1 Analytic Solution: Rotational Balance with Toroidal Magnetic Stress

In this section, we derive a special analytical solution for the case where centrifugal forces exactly balance magnetic stresses from the toroidal field. We acknowledge that this solution is very specific, and therefore lacks the generality of the numerical solutions presented later in this paper. All quantities are presented in dimensionless form, where we have dropped the overscript tilde notation.

Various simplifying terms and assumptions are introduced so that the momentum equation (equation 4) is simplified, which leads to an interesting analytical solution. We first introduce a “total” pressure is defined by

[TABLE]

Our analytical solution requires that the centrifugal force balances exactly with the inward toroidal magnetic stress, so that

[TABLE]

We further assume that the ratio of pressure and total magnetic pressure is constant:

[TABLE]

We define

[TABLE]

where $\beta_{z}$ and $\beta_{z}$ are also constant, so that equation 10 can be written as

[TABLE]

Using the definition of $\beta_{\phi}$ in equation 9, and the definitions of $\sigma$ from equation 1d, it is straightforward to show that

[TABLE]

Equation 13 implies that the velocity field is highly sheared at $r=0$, which violates the usual boundary condition $v_{\phi}(0)=0$. Thus our special solution may not be physical very close to the axis $r=0$. Interestingly, Kaur et al. (2006) have also considered constant $v_{\phi}$ models.

The total pressure is then given by the sum of the gas pressure (possibly including non-thermal contributions) and the magnetic partial pressures according to

[TABLE]

Furthermore, we define an effective sound speed given by

[TABLE]

where dimensional units are restored for this equation only, so that we can redefine the dimensionless gravitational potential

[TABLE]

which is used to express Poisson’s equation for gravity in dimensionless form. Once again dropping the tilde notation, by assuming that centrifugal force and toroidal magnetic stresses exactly balance (equation 9), and incorporating equations 14 and 16, equation 4 becomes

[TABLE]

which can be integrated to find a solution for $\rho$ so that Poisson’s equation can be written as

[TABLE]

where $r_{0}$ is redefined as

[TABLE]

Equation 18 was solved in cylindrical symmetry by Ostriker (1964) in his study of non-rotating, non-magnetic, isothermal filaments. Restoring dimensions, we therefore obtain an analogous solution for a rotating, magnetic filament:

[TABLE]

Thus, our analytic solution is a rescaling of the Ostriker solution, where the core radius $r_{0}$ is modified by the magnetic field. The filament rotates with constant velocity. The density inversions noted by Hansen et al. (1976) and Recchi et al. (2014) are absent from this special solution.

Nagasawa (1987) studied the stability of non-rotating isothermal cylinders threaded by a constant, trivial axial magnetic field. The lack of magnetic stresses from such a constant magnetic field do not modify the equilibrium structure from the Ostriker solution. In contrast, the analytical solution presented above has a radial scale modified by both axial and toroidal magnetic stresses. The solutions we explore throughout the rest of this work include non-trivial magnetic stresses.

2.2 Differential Rotation

The special solution presented in Section 2.1 represents a finely tuned special case, which we now generalize. Rotation plays an important role in the morphology of molecular tornadoes as seen in observations (Matsumura et al., 2012; Sofue, 2007; Morris et al., 2006) and past work on rotating filaments (Hansen et al., 1976; Recchi et al., 2014). It is useful to examine several possible rotation laws in order to accurately describe and understand the physics of molecular tornadoes.

2.2.1 Angular Frequency as a Power Law

We consider a rotation law where $\Omega$ scales with the radius $r$ according to a power law with index $\alpha$:

[TABLE]

Keplerian rotation and solid body rotation are considered to be limiting cases. We show below that this assumption constrains $\alpha$ to the range $-1\leq\alpha\leq 0$, where $\alpha=0$ is the case for solid body rotation. For clarity, we note that the $\alpha=-1$ index would be more accurately described as “Keplerian-like” where gas orbits a line mass in the absence of pressure gradients.

2.2.2 Keplerian Rotation

Poisson’s equation and the divergence theorem are used to find the following relationship for $g$ in terms of $r$: $2\pi grL=4\pi GM(r)$, where $M$ is the filament mass and $L$ is its length. By solving for $g$ and introducing the linear mass $m=M/L$, equating the angular acceleration and gravitational acceleration implies that

[TABLE]

Thus we find that Keplerian rotation is just a special case of the radial power law (equation 21) with $\alpha=-1$.

2.2.3 Rotational Frequency as a Free Function of Magnetic Field Lines

We follow the methods introduced by Mouschovias (1976) and Tomisaka et al. (1988) to describe the rotational frequency as a function the magnetic flux $\Phi_{M}$, so that $v_{\phi}=\Omega(\Phi_{M})r$. Assuming $\Gamma_{z}$ is constant, we find that $\partial\Phi_{M}/\partial r=2\pi\Gamma_{z}\rho r$. The angular frequency $\Omega$ is solved numerically from $\partial\Phi_{M}/\partial r$ along with other quantities mentioned later in Section 2.5. The rotational frequency is then assumed to follow a power law

[TABLE]

The lower limit on $\alpha_{M}$ is determined by the proportionality relations: $v_{\phi}\propto\Omega r\propto\Phi_{M}^{\alpha_{M}}r$ and $\Phi_{M}\propto B_{z}r^{2}\propto r^{2}$, so that $v_{\phi}\propto r^{2\alpha_{M}+1}$. For $2\alpha_{M}+1>0$, this indicates that $\alpha_{M}>-1/2$. For the sake of interest, we also analyze $\alpha_{M}\geq-1$ to understand how $\alpha_{M}$ affects the behaviour of the radial density.

2.2.4 Torsional Alfvén Waves

Appendix A derives equations for small amplitude oscillatory torsional Alfvén waves. Here, we present an analogous derivation for the wave front of a large amplitude torsional Alfvén wave, which is important for the numerical solutions presented in Section 3. In both Appendix A and this large amplitude scenario, we assume that the unperturbed rotation velocity and magnetic fields are described by

[TABLE]

The MHD equations (equation 1) yield

[TABLE]

which is the general wave equation describing transverse perturbations propagating along the flux tube at the Alfvén speed (Jackson, 1975)

[TABLE]

Consider an Alfvén wave propagating locally along the flux tube, where the rotation velocity and magnetic fields are finite amplitude perturbations. Furthermore, $B_{z}$ and $\rho$ are constant in $z$ and $t$. If an initially stationary flux tube starts to rotate at its base at time $t=0$, the perturbation travels a distance $v_{A}dt$ in time $dt$, displacing the flux tube by distance $v_{\phi}dt$. As the wave propagates along the flux tube, a toroidal magnetic field $B_{\phi}$ is generated. Consequently, Faraday’s Law reduces to

[TABLE]

The interface ($z=0$) of the propagating wavefront separates the conditions above ($z>0$) and below ($z<0$) the wavefront. Above the wavefront interface, $B_{z}$ is constant, and $B_{\phi}=0$. Below the interface, the magnetic field has been perturbed and exhibits a toroidal component $B_{\phi}$, and rotates at $v_{\phi}$. A torsional Alfvén wave propagates along a magnetic flux tube at the Alfvén speed $v_{A}=dz/dt$. By integrating Faraday’s Law (equation 27), substituting $dt$ from $v_{A}$, and introducing the aforementioned bounds, we obtain

[TABLE]

and follow through with the integration. Thus the field angle $\theta$ is perturbed according to

[TABLE]

We use the flux-to-mass ratios (equation 5a and equation 5b) in equation 29, as well as the Alfvén speed ($v_{A}$) to find

[TABLE]

where $\Gamma_{\phi}=\Gamma_{\phi}(r)$, in general.

2.3 Virial Analysis

The work by Sofue (2007) suggests that the GCT may not be bound by self-gravity. By extension, it may be possible that other molecular tornadoes may not require significant self-gravity to hold together. It is helpful to conduct a virial analysis in order to investigate the relative energies associated with various physical properties such as self-gravity, external pressures, magnetic fields, turbulence, and rotation.

Following Fiege & Pudritz (2000a), the virial equation that describes our models is given by

[TABLE]

where script characters denote quantities per-unit-length; thus $\mathcal{K}$ is the total kinetic energy per unit length, $\mathcal{M}$ is the magnetic energy per unit length, and $\mathcal{W}$ is the gravitational energy per unit length. The total kinetic energy contains the rotational kinetic energy $\mathcal{K}_{rot}$, surface pressure term $\mathcal{K}_{P}$, and internal turbulence $\mathcal{K}_{\sigma}$ so that

[TABLE]

The rotational kinetic energy is given by $\mathcal{K}_{rot}=(1/2)\mathcal{I}\Omega^{2}$, where $\mathcal{I}$ is the rotational inertia per unit length of a cylindrical filament spinning along its vertical axis. Energy associated with external pressure is given by $\mathcal{K}_{P}=-(3/2)P_{S}\mathcal{V}$ where $\mathcal{V}=\pi r_{S}^{2}$ is the cross-sectional area of the filament. The energy due to internal turbulence is given by $\mathcal{K}_{\sigma}=(3/2)m\sigma^{2}$. Gravitational energy is described by $\mathcal{W}=-m^{2}G$, which is valid for any equation of state, magnetic field, and internal structure (Fiege & Pudritz, 2000a). Rearranging equation 31 and normalizing by $|\mathcal{W}|$ yields

[TABLE]

which is useful to find the relative magnetic energy components scaled by the gravitational energy. Determining the relative energies of these properties allows us to gain insight into which ones are particularly important for molecular tornadoes. Equation 33 makes it clear how terms contribute to the sign of $\mathcal{M}/|\mathcal{W}|$, which determines the dominating component of the magnetic field (Fiege & Pudritz, 2000a). For $\mathcal{M}/|\mathcal{W}|<0$, the toroidal component is dominant in pinching the filament, while $\mathcal{M}/|\mathcal{W}|>0$ implies net support by the poloidal field. The numerical results in Section 3 are examined with the help of equation 33.

2.4 Guidelines for Stability

Shafranov (1956) investigated the stability criterion for MHD instabilities. Unlike our equilibrium model for molecular tornadoes, the filaments discussed by Shafranov (1956) do not include rotation, self-gravity, nor are they truncated by external pressure. Nevertheless, Shafranov (1956) provides a guideline for the stability analysis to be performed on our equilibrium model of molecular tornadoes, which will be presented in a future paper.

In the case where there is no $B_{z}$ external to the cylinder, the condition for instability is

[TABLE]

(Shafranov, 1956), where the magnetic field pitch angle relationship (equation 29) at the cylinder surface is used, and $I_{m}$ and $K_{m}$ are modified Bessel functions of the first and second kind, respectively. In the case where there is no $B_{z}$ external to the cylinder, the $m=0$ (sausage), $m=1$ (kink), and $m=2$ modes are unstable when

[TABLE]

Equivalently, the stability conditions above can be expressed in terms of the rotational period $P=2\pi/\Omega_{S}$, and an Alfvén crossing time $\tau_{A}:=\lambda/v_{A}$ that describes the time it takes for an Alfvén wave to propagate across one wavelength or pitch $\lambda$. The stability conditions in equation 35 for $m=0,1,\text{and }2$, respectively, are expressed as

[TABLE]

Equation 35 suggests that filaments that are stable against $m=1$ instabilities are also stable against $m=2$ instabilities. In other words, filaments that spin faster may still be stable against $m=2$, but not $m=1$ instabilities. It may also help to understand the stability conditions based on the axis ratio $r_{S}/\lambda$, as expressed in equation 36. For a given $v_{\phi}/v_{A}$, increasing the axis ratio will first stabilize against $m=2$ instabilities before stabilizing against $m=1$ instabilities. It will be interesting to see how these stability conditions described by Shafranov (1956) will compare to our molecular tornado models. Although Shafranov (1956)’s models are not exactly like those described in this paper, these stability relations are a useful guideline. In particular, equation 35 is used to estimate the angular frequency of our equilibrium models in Section 2.6 as it sets a lower bound for the angular frequency required for the observed instability to occur. For example, the minimum rotation rate required for a filament to undergo a kink instability would be $\Omega_{S}>\sqrt{2}v_{A}/r_{S}$, according to equation 35b. In other words, a molecular tornado observed with a kink instability must have been spinning at $\Omega_{S}>\sqrt{2}v_{A}/r_{S}$ for the instability to trigger.

2.5 Squashing Transformation

To efficiently calculate a wide range of $r$, and to better study the large $r$ asymptotic behaviour of our equations, we employ a convenient “squashing” transformation to more easily numerically integrate out to large $r$. The squashing transformation transforms the dimensionless radius $r$ into a dimensionless quantity $\xi$ defined by $\xi\equiv\ln(r)$. The radial differential operators therefore transform according to

[TABLE]

The equations for $g(r)=-\partial{\Phi(r)}/\partial r$, $dg/dr$ (equation 6), and $d\rho/dr$ (equation 7) are trivially transformed to equivalent equations using $\xi$ as the independent variable, which is what we solve numerically.

In the case of equation 30, where $B_{\phi}$ is due to a torsional Alfvén wave, it is useful to further simplify Equation 7. If $\Gamma_{z}$ is constant, then equation 7 can be expressed explicitly in terms of the dimensionless quantities $\Omega$, and $d\Omega/d\xi$ by using the flux-to-mass ratios (equation 5) along with the expression above for $\Omega$ (equation 30). Algebraic simplification results in the expression

[TABLE]

where we have used equation 30 to write $\Gamma_{\phi}$ in terms of $\Omega$.

2.6 Dimensionless Quantities: Comparing Observations with Numerical Models

Observations of molecular tornadoes measure some useful global quantities that can be used to constrain our solutions to equation 38. We denote quantities at the surface of the filament with the subscript $S$, the radius from the filament axis to its surface $r_{S}$, surface pressure $P_{S}$, linear mass $m_{S}$, and rotational frequency $\Omega_{S}$. These observables are written in terms of the dimensionless quantities of our model in the same manner as the scaling laws of equation 2. Thus, there are five unknowns in the scaling laws: $\tilde{r}_{S}$, $\tilde{P}_{S}$, $\tilde{m}_{S}$, $\tilde{\Omega}_{S}$, and $\rho_{0}$. The readily observable (or at least somewhat observationally constrained) quantities are $\sigma$, $r_{S}$, $m_{S}$, and $P_{S}$. The central density $\rho_{0}$ cannot be easily obtained observationally, which also makes $r_{0}$ difficult to obtain directly. Thus, we use the definition of $r_{0}$ to eliminate $\rho_{0}$ from the equations. This reduces the scaling laws for $\tilde{r}_{S}$, $\tilde{P}_{S}$, $\tilde{m}_{S}$, $\tilde{\Omega}_{S}$, and $\rho_{0}$ to three equations that describe combinations of the dimensionless quantities on the left side, in terms of observables on the right:

[TABLE]

The surface pressure can be reasonably constrained from estimates of the isothermal equation of state, as well as from the dynamics and surface density of stars (Rathborne et al., 2014; Swinbank et al., 2011). Rathborne et al. (2014) and Swinbank et al. (2011) suggest that $P_{S}/k\gtrsim 10^{7}$ K cm*-3*, and can possibly be as high as $\sim\,10^{8}$ K cm*-3*. The velocity dispersion for molecular clouds in the CMZ are roughly $\sigma\,\sim\,15$ km s*-1* (Rathborne et al., 2014). The angular frequency of the filament in equilbrium is estimated from the Shafranov stability relations of equation 35 as well as from the rotation mechanism, where possible.

2.6.1 Pigtail, DHN, and GCT

Our exploration focusses on the Pigtail to constrain our models. Equation 39 shows that there are five parameters that are necessary to constrain our dimensionless models by observations, namely $\tilde{r}_{S}$, $\tilde{P}_{S}$, $\tilde{m}_{S}$, $\tilde{\Omega}_{S}$, and $\sigma$. Most of these parameters were derived by observations of the Pigtail. The only exception was for rotation, which was not detected due to resolution limitations (Matsumura et al., 2012). The Pigtail’s velocity dispersion is $\sigma=13.2$ km s*-1*, and $n(H_{2})=10^{3.5\pm 0.25}$ cm*-3* ($\sigma$ and $n(H_{2})$ provided via private communications with Dr. Tomoharu Oka, Keio University; see also Matsumura et al. (2012)). Assuming hydrogen abundance $X=0.7\pm 0.05$, we estimate via the isothermal equation of state that $P_{S}/k=(2\pm 1)\times 10^{8}$ K cm*-3*. A direct measurement of the rope of the Pigtail from the figures of Matsumura et al. (2012) suggests that its radius is between $2-4$ pc. The pressure and velocity dispersion are assumed to be constant, and so it may be reasonable that the volume of the Pigtail rope between equilibrium and its present form is conserved. Thus, the radius of the filament in equlibrium $r_{S}$ and the radius of the filament when the Pigtail has undergone instability $r^{\prime}_{S}$, are related by the pitch angle via $r_{S}=r^{\prime}_{S}\cos(\theta)$. The rotation speed of the equilibrium filament can be estimated either by considering that the flux tube is spun up by the shearing of two cloud orbits (Matsumura et al., 2012) or by the stability calculations of Shafranov (1956) mentioned in Section 2.4. In this case, we opt for the shearing mechanism described by (Matsumura et al., 2012) to derive the rotation frequency as the Shafranov stability conditions have underlying assumptions that differ from our models. Although, using the $m=1$ stability criterion (equation 35b) with the equilibrium radius and Alfvén speed $v_{A}\sim 35$ km s*-1* (Matsumura et al., 2012) produces a very close estimate (well within one order of magnitude) of the angular frequency to that derived by considering the shearing mechanism. These parameters result in dimensionless constraints

[TABLE]

There is more uncertainty in the physical parameters of the DHN and GCT, however. In particular, $\sigma$, $\Omega_{S}$, and $P_{S}$ are unclear. Generally, molecular clouds of the CMZ exhibit $\sigma\,\sim\,15$ km s*-1* (Rathborne et al., 2014) which is similar to the Pigtail’s velocity dispersion. For the DHN and GCT, we consider $10^{7}$ K cm*-3* $\lesssim P_{S}\lesssim 10^{8}$ K cm*-3* (Rathborne et al., 2014) as a conservative estimate of the surrounding pressure. This range is slightly lower than the value that we have calculated for the Pigtail, and this discrepancy may be due to the density of the Pigtail, which could be lower than $n(H_{2})=10^{3.5\pm 0.25}$ cm*-3* (Matsumura et al., 2012) on its surface, and/or a hydrogen abundance larger than $X=0.7$. Our analysis focusses primarily on the Pigtail because observational measurements for the Pigtail have been more conclusive. Still, we attempt to carry out our analysis with the DHN and GCT, but the results should be taken with more caution.

CO observations have revealed two molecular counterparts that are likely associated with the DHN at radial velocities of $\sim-35$ km s*-1* and [math] km s*-1* having lengths $\sim\,150$ pc, and masses of $0.8\times 10^{4}\,M_{\odot}$ and $3.3\times 10^{4}\,M_{\odot}$, respectively (Torii et al., 2014; Enokiya et al., 2014). The equilibrium radius of the DHN filament is conservatively estimated to fall between the narrowest strand width and the overall structure width $\sim 3.5$ pc (Morris et al., 2006). The CND rotation is the proposed mechanism which drives the DHN. We assume that the CND rotation lies in the range $70-110$ km s*-1* (Morris et al., 2006; Enokiya et al., 2014). Alternatively, with $r_{S}$ and $v_{A}=10^{3}$ km s*-1* (Morris et al., 2006), an estimate of the angular frequency can be made via the $m=2$ Shafranov stability relation of equation 35c. We opt to use the angular frequency determined from the CND rotation for the analysis however, since there are assumptions made in equation 35c that differ from our models (see Section 2.4). Regardless, the general qualitative conclusions drawn from either method of determining the rotation are similar (see Section 4.3). These parameters along with the pressure estimate from Swinbank et al. (2011) and Rathborne et al. (2014), correspond to

[TABLE]

Assuming $\sigma=15$ km s*-1*, and provided that the GCT has a density of $\rho\,\sim\,2.7\times 10^{-21}$ g cm*-3* (Sofue, 2007), the isothermal equation of state suggests an internal pressure of $4.4\times 10^{7}$ K cm*-3*. As a rough guideline, this is within the bounds of the external pressure that we consider. Sofue (2007) draws the similarity between the DHN and GCT that they both exhibit two strands wound in a double helix configuration, and postulates that they arise from a similar mechanism. Thus, we assume that the radius of the GCT is approximately $7-20$ pc wide, as estimated from the current radius of the GCT, and the size of the molecular cloud it may have originated from (Sofue, 2007). The GCT’s angular frequency is estimated from the Shafranov (1956) instability condition (equation 35c). We acknowledge that the Shafranov conditions of equation 35 describe minimum bounds for instability, and that the rotation does not have an upper bound. Realistically however, there must be an upper rotation limit otherwise the filament would be unstable and the filament would not remain intact. Even if the GCT were rotating faster than the rate described by equation 35c, it should not be much faster before undergoing instability. Thus, we conservatively estimate limits on the angular frequency resulting from equation 35c by multiplying by a factor of 0.5 to 2. These parameters correspond to

[TABLE]

Most of the dimensionless constraints of equations 40, 41, and 42 are based on absolute minimum and maximum bounds of the observed measurements, except where otherwise noted. These quantities are chosen to conservatively estimate the bounds of equation 39 for each molecular tornado, in an attempt to account for the uncertainty in the observations (particularly with uncertainties associated with the DHN, and GCT).

3 Results

Results of numerically solving equations 7 and 38, along with $d\Phi/d\xi$ and $dg/d\xi$ , and the rotation laws described in Section 2.2 are presented in this section. The density $\rho$ is calculated from the general form of the Cauchy momentum equation (equation 7) for models excluding the Alfvén wave constraint (equation 30). Equation 38 was derived under the constraint of equation 30. The radial profile of $\rho$ changes due to different rotation laws (equations 21, and 23), magnetic fields, and self-gravity. These profiles are compared to the Ostriker solution, equation 20, for isothermal, non-rotating, non-magnetic, self-gravitating cylinders as a benchmark. At large radii, the Ostriker solution behaves as $\rho\,\sim\,r^{-4}$.

In our exploration, we are free to explore many different filament models. We have the flexibility to include or exclude self-gravity, choose between various rotation laws and their respective power law indices (see Section 2.2), and define the allowed range of $r$, $\Omega_{0}$, $\Gamma_{z}$, and $\Gamma_{\phi}$. While all of these models are interesting to explore, the most physically relevant models make use of the magnetic field angle constraint from torsional Alfvén waves, presented in Section 2.2.4. However, we also explore models that exclude this magnetic field constraint for comparison.

Numerical integration was performed with MATLAB’s ode45, which is a Runge-Kutta 4-5 integrator. In a few cases where stiffness was detected, we used MATLAB’s ode23s, which is a low order method based on the Rosenbrock formula, and efficiently solves rapidly changing problems. We integrate over the range $r_{min}\leq r\leq r_{max},$ where the lower and upper limits to the dimensionless radius are $r_{min}=10^{-6}$ and $r_{max}=10^{6}$, respectively. If $\Omega$ follows the power law rotation given by equation 21, then $\Omega_{0}$, $\Gamma_{z}$, and $\Gamma_{\phi}$ are free parameters. In this case, $\Omega_{0}$, $\Gamma_{z}$, and $\Gamma_{\phi}$ are each randomly chosen over the intervals

[TABLE]

The limits of equation 43 are conservatively chosen, with the limits of $\Gamma_{\phi}$ and $\Gamma_{z}$ extended beyond that of Fiege & Pudritz (2000a) to account for the environment of the Galactic Centre. Each random set of $\Omega_{0}$, $\Gamma_{z}$, and $\Gamma_{\phi}$ define a particular untruncated filament model. Otherwise if $\Omega$ follows a power law of magnetic flux (equation 23), then $\Gamma_{z}$, and $\Gamma_{\phi}$ are free, randomly chosen over the same interval is in equation 43, and define a unique density profile. The initial conditions at $r_{min}$ are

[TABLE]

The reader is reminded that all quantities in our numerical calculations are dimensionless.

Theoretical density profiles from solving equation 7 (and 38) are shown in Figures 1 and 2. We halt the integration if a density inversion is encountered ($d\rho/dr>0$), as recommended by Recchi et al. (2014). This condition also limits the density so that $\rho(r)\leq\rho_{0}$ everywhere. Some density profiles may be asymptotically constant at large radii, so these profiles were noted accordingly, but were allowed to integrate over the full range of $r$. For models including self-gravity, constant angular frequency, constant flux-to-mass ratios, and excluding the torsional Alfvén wave condition, changing $\Gamma_{z}$ changes the asymptotic density value that is reached at large radii only by a small amount. This asymptotic value of $\rho$ at large radii can be approximated from equation 7, which reveals that

[TABLE]

Truncation pressures were chosen randomly, with half chosen from a uniform distribution, and half chosen from a logarithmic distribution. This choice allows us to sample a wide range of possible external pressures that a filament may reside in. The logarithmic distribution allows the truncation values to span the entire range of $r$, which is particularly useful for probing the low density tail that most density profiles exhibit at large $r$. The uniformly distributed truncation values better sample where observationally constrained models are more likely to reside (discussed further below; see Figure 3 as well). Solutions that exhibit density inversions are halted at the first density inversion, then further truncation pressures are chosen along the truncated density profile with the aforementioned distribution. At each point, a number of important parameters ($\rho$, $P_{S}$, $\Omega_{S}$, $\Gamma_{\phi}$, $\Gamma_{z}$, etc.) are evaluated. In total, we sample approximately $10^{4}$ random combinations of the model parameters (equation 43), each describing a unique density profile for a particular rotation law. Additionally, each external truncation pressure characterizes a unique filament model. Since there are approximately $10^{4}$ external pressures, we explore approximately $10^{8}$ models for a given rotation law.

The parameters evaluated at each truncation pressure ($\rho$, $P_{S}$, $\Omega_{S}$, $\Gamma_{\phi}$, $\Gamma_{z}$, etc.) are used to constrain the parameter space where observed molecular filaments reside, according to the constraints of equation 39. By converting the observed parameters into dimensionless quantities via equation 39 (shown explicitly in Section 2.6.1), we plot where real filaments reside in the dimensionless parameter space (see Figures 3 and 4). The density profiles are also categorized based on their general behaviour or special characteristics. For example, some theoretical density profiles feature inversions (see Figure 5), some reach a constant density (at large radii), and some increase above $\rho_{0}$ (see Figures 2 and 6).

Having constrained the models within observational bounds, they can then be used to analyze the virial equation (equation 33). The result of this is shown in Figure 7 for self-gravitating models following the radial power law (equation 21) with the torsional Alfvén wave condition (equation 30), constrained by the Pigtail. Within observational constraints, our analysis shows that the magnetic-to-gravitational energy ratio is

[TABLE]

and the ratio of the surface pressure term to gravitational energy is

[TABLE]

Comparing the surface pressure term to the magnetic energy reveals that

[TABLE]

for most self-gravitating models with equation 30 within observational constraints of the Pigtail as shown in Figure 8. These results are discussed further in Section 4.3.

4 Discussion

In this section, we discuss the results presented in Section 3. We compare the behaviour of our models when varying different parameters, to the Ostriker solution. Models that exclude torsional Alfvén wave physics (equation 30) in Section 4.1 are first discussed, and in Section 4.2 we discuss models that include the constraints of torsional Alfvén waves. The results of our virial analysis are discussed in Section 4.3.

4.1 Excluding the Torsional Alfvén Wave Condition

Not surprisingly, the effect of self-gravity is to keep the filament together against the centrifugal force. In the absence of significant self-gravity, $\rho$ often tends to increase with radius and may exceed $\rho_{0}$ due to rotation, which is not realistic for a filament that eventually merges with the ISM. Several solutions of this type can be seen in Figures 2 and 6. Typically, this behaviour is associated with strong rotation and weak toroidal magnetic fields. However, if magnetic stresses are strong enough, and/or rotation slows within a reasonable radius, then the density profile may decrease. Empirically, if a power law in radius (equation 21) is followed (see Figure 6), then a greater $\Gamma_{\phi}/\Omega_{0}$ ratio tends to keep $\rho<\rho_{0}$ and $d\rho/dr<0$. The effect of $\Gamma_{z}$ does not seem to be as significant compared to $\Gamma_{\phi}$ and $\Omega_{0}$, unless $\Gamma_{z}$ is large. These features can be seen in Figure 6.

If the flux power law of rotation (equation 23) is followed (see Figure 6) with self-gravity neglected, then the value of $\alpha_{M}$ seems to have a dominating effect on the behaviour of $\rho(r)$. There are density profiles where $\rho>\rho_{0}$, regardless of the value of $\alpha_{M}$. A larger portion of these models exceed $\rho_{0}$ as $\alpha_{M}$ becomes steeper. All of these models, except where $\Gamma_{\phi}=0$, decrease in density as $r$ becomes large. The $\Gamma_{\phi}/\Omega_{0}$ ratio has a similar effect as in the models that follow equation 21. See Figure 6 for the associated density profiles.

4.1.1 Density Inversions

Density inversions are found in some models where $\rho(r)<\rho_{0}$ (for all $r$), but curiously only if self-gravity is included, and the field angle constraint in equation 30 is excluded. An important difference between our models, and the models of Hansen et al. (1976) and Recchi et al. (2014) is the inclusion of magnetic fields. Even without the constraint of torsional Alfvén waves (equation 30), if there are constant poloidal and toroidal magnetic fields, then density inversions can still be observed as long as $\Gamma_{\phi}$ is weak.

It is easier to see the effects due to $\Omega_{0}=\Omega(r)$, $\Gamma_{\phi}$, and $\Gamma_{z}$, by first assuming that they are constant. The effect of increasing $\Gamma_{\phi}$ is to dampen the density inversions (see Figure 5). As $\Omega_{0}$ increases, the entire density profile along with the inversions, shifts to greater densities, and the inversion amplitude decreases slightly. The effect of increasing $\Gamma_{z}$ is that the density inversions tend to occur at larger radii as if each inversion was shifted over, and with larger amplitude. The toroidal flux-to-mass ratio $\Gamma_{\phi}$ seems to have a more noticeable effect on the behaviour of $\rho$ than $\Gamma_{z}$. Generally, $\Gamma_{\phi}$ and $\Gamma_{z}$ behave as Fiege & Pudritz (2000a) describe – $\Gamma_{z}$ tends to support/radially extend the filament compared to the Ostriker solution, while $\Gamma_{\phi}$ pinches the filament. These effects can still be seen even with density inversions occurring (see Figure 5).

Density inversions exist in solutions where the rotation follows the radial power law of equation 21 (see Figure 5), and for those that follow the flux power law of equation 23. At large radii, the density decreases like a power law for $-1\leq\alpha<0$. The general behaviour of the density inversions for both rotation laws is similar. There is a subtle difference in that the inversions following the flux power law seem to exhibit a smaller frequency and features are stretched over a slightly larger radius than those that follow equation 21. While a particular model may exhibit density inversions, it is only realistic if it is pressure truncated before any inversions are observed.

4.2 Including the Torsional Alfvén Wave Condition

Self-gravitating filaments following the radial power law of equation 21 all have density profiles that decrease with radius (see Figure 1). A steeper power law index tends to cause a steeper decline in density, as seen in Figure 1. A large $\Gamma_{z}/\Omega_{0}$ ratio tends to produce steeper solutions initially, but then achieve an asymptotically constant ($d\rho/dr\approx 0$) value of $\rho$ at large radii. As $\alpha$ decreases, asymptotic solutions can be seen with smaller $\Gamma_{z}/\Omega_{0}$ ratios. Non-self-gravitating filaments following the same rotation law generally have solutions that are not asymptotic. As $\alpha$ becomes progressively steeper, the density profiles tend to become more shallow. At $\alpha=-1$, all profiles are truly asymptotically constant at $\rho=\rho_{0}$. All profiles decrease when $\alpha>-1$ .

Models that are self-gravitating and rotate according to the flux power law (equation 23) have decreasing density when $\alpha_{M}=-0.25$, but more negative values have profiles that were deemed unrealistic because $\rho>\rho_{0}$ at some point(s) (see Figure 2). All profiles begin to decrease at a sufficiently large radius, but the density inversion, where $d\rho/dr>0$, becomes more pronounced as $\alpha_{M}$ becomes more negative. The existence of this extreme density inversion persists even after numerous tests of the numerical integration. Thus, we deem these models to be physically unrealistic and exclude them. As $\alpha_{M}$ becomes more negative, the profiles following the flux power law for rotation (equation 23) have a greater frequency of asymptotic solutions where there is a relatively large $\Gamma_{z}/\Omega_{0}$ ratio. The ratio of $\Gamma_{z}/\Omega_{0}$ does not need to be as large to observe asymptotic solutions as $\alpha_{M}$ decreases. Non-self-gravitating filaments following the same rotation law do not have asymptotic solutions at $\alpha_{M}=-0.25$, and the density profiles decrease with the radius. Density profiles with $\alpha_{M}\lesssim-0.50$ were deemed unrealistic due to density inversions that increase the density above $\rho_{0}$ at some point(s). This $\alpha_{M}\approx-0.50$ boundary coincides with the limit of $\alpha_{M}$ that was stated in Section 2.2.3. More negative values of $\alpha_{M}$ cause $\Omega$ to behave as $r^{-1}$ or steeper, and the filament rotates extremely rapidly near the origin. These models are deemed unrealistic.

Generally, by including the physics of Alfvén waves through the description in Section 2.2.4, we see that the density profiles are better behaved in the sense that the density profiles usually do not increase above $\rho_{0}$, the density falls off with the radius, and there are no density inversions. Interestingly, every model that follows the radial power law rotation of equation 21 exhibit $\rho(r)\leq\rho_{0}$ (or are at least asymptotically constant at $\rho=\rho_{0}$ when $\alpha=-1$ for non-self-gravitating models), in both self-gravitating and non-self-gravitating cases regardless of any $\alpha<0$. Furthermore, by including the Alfvén wave constraint (equation 30), every model that follows the flux power law rotation of equation 23 also exhibit $\rho(r)\leq\rho_{0}$ in both self-gravitating and non-self-gravitating cases, as long as $\alpha_{M}>-1/2$. This suggests that the inclusion of Alfvén wave physics is an important factor in the description of molecular tornadoes.

4.3 Is Self-Gravity Important?

Self-gravity, magnetic fields, and external pressure have a significant influence on constricting the density profile in our models. On the other hand, rotation, and turbulence tend to increase the density with radius. Interestingly, the observationally constrained models (see Figure 8) indicate via equation 33 that, for all models, $\mathcal{M}/|\mathcal{W}|<0$. This is not surprising since we have noticed the pinching effect that magnetic fields have on density profiles (see Section 4.1.1, Figure 6), and the magnetic-squeezing/twisting mechanism that is proposed by Matsumura et al. (2012); Morris et al. (2006); Sofue (2007).

Having found the relative energies via our virial analysis (see Section 3), $\mathcal{W}$ is compared to $\mathcal{M}$ and $\mathcal{K}_{P}$, which are the quantities that can help keep the filament structure bound, to see which terms are particularly important. An example of these comparisons is shown in Figure 7 for the self-gravitating model with the torsional Alfvén wave constraint of equation 30, and radial power law of equation 21 with $\alpha=-0.25$. The magnetic energy ratio with gravitational energy within observational constraints is found to be $-10^{5}\lesssim\mathcal{M}/|\mathcal{W}|\lesssim 0$, which suggests that the magnetic field may be much more important than self-gravity. Furthermore, the ratio of the surface pressure to gravitational energy across all self-gravitating models within observational constraints is $10^{0}\lesssim|\mathcal{K}_{P}|/|\mathcal{W}|\lesssim 10^{3}$, which further indicates that these molecular tornadoes are only weakly self-gravitating. It is also found that $|\mathcal{M}|>|\mathcal{K}_{P}|$ for most models within observational constraints, as shown in Figure 8. Thus, the toroidal magnetic field confinement is more important than the surface pressure, and both are more important than self-gravity. This result is similar to the interpretation by Sofue (2007), who considered the magnetic tension to be a significant component, and did not focus on the external pressure. Knowing this, the approximation $\mathcal{W}\approx 0$ may be justified. If the toroidal magnetic stress is dominant, equation 33 may best be rewritten as

[TABLE]

where it is easier to see that if $\mathcal{M}<0$ then $(\mathcal{K}_{rot}+\mathcal{K}_{\sigma})<-\mathcal{M}/2$, and if $\mathcal{M}>0$ then $(\mathcal{K}_{rot}+\mathcal{K}_{\sigma})>-\mathcal{M}/2$. It was also found that $|\mathcal{M}|>|\mathcal{K}_{P}|$, and since models follow $\mathcal{M}<0$, this means that $-(1/2)|\mathcal{M}|<(\mathcal{K}_{rot}+\mathcal{K}_{\sigma})<-(3/2)|\mathcal{M}|$.

A similar analysis of the DHN and GCT is shown in Figure 9. All observationally constrained models of the GCT are dominated by the toroidal magnetic stress component, and the vast majority of models are magnetically dominated over the pressure. This lends support for the magnetic squeezing model by Sofue (2007). Analysis suggests that the dominant component of the DHN is not clear – pressure or magnetic stresses may dominate. Comparison between $|\mathcal{M}|$ and $|\mathcal{K}_{P}|$ for the DHN reveals that more of the constrained models are magnetically dominated, though the dominating component cannot be conclusively determined. This uncertainty is likely due to estimating $\Omega_{S}$ By the torsional Alfvén wave condition (equation 30), and as seen in Figures 7 and 9, rotation and the toroidal magnetic stress are directly proportional. A more accurate method of estimating the equilibrium filament radius would help constrain $\Omega_{S}$, and therefore, determine whether magnetic or pressure terms dominate. If the rotation had instead been estimated by the Shafranov condition (equation 35c), then a much faster rotation would have been resulted. This would increase $|\mathcal{K}_{rot}|$, and would then require a much stronger $|\mathcal{M}|$ to satisfy the virial relation. Thus, the DHN would convincingly be magnetically dominated. It is highly suggestive, however, that both the DHN and GCT are weakly self-gravitating, like the Pigtail Molecular Cloud. We do not analyze the results of the DHN and GCT in detail as the observational constraints carry significant uncertainty associated with observational measurements.

5 Summary

Molecular tornadoes are fascinating objects that reside in the extreme environment of the CMZ. We study equilibrium MHD models of molecular tornadoes and introduce different rotation laws (see Section 2.2), reasonably general helical magnetic fields (equations 5 and 7), external pressures, and consider an isothermal equation of state.

The physics of torsional Alfvén waves was introduced in Section 2.2.4, which led to a constraint on the magnetic field components (equation 29), and between $\Gamma_{\phi}$ and $\Omega$ (equation 30). This constraint leads to solutions that are generally more well-behaved than solutions excluding the torsional Alfvén wave condition – in the sense that their density profiles do not increase indefinitely, and there are no density inversions. We conclude that including torsional Alfvén wave physics is an important component in a realistic model of molecular tornado structure. 2. 2.

A special analytical solution, where centrifugal forces balanced exactly with toroidal magnetic stresses, was explored. The solution also assumes that the ratio between the pressure (from the isothermal equation of state) and the magnetic stresses, is constant (equation 10, and 11). These assumptions leads to solutions that are rescaled versions of the Ostriker solution, with constant rotational velocity of the filament (equation 20). 3. 3.

A Monte Carlo analysis of our models was conducted to explore the associated parameter space. Our models are constrained by limited observable constraints via equation 39. We focus our analysis on the Pigtail Molecular Cloud whose observational properties are more conclusively measured. Our analysis suggests that external pressures are important in its equilibrium structure (see Figure 8). 4. 4.

We performed a study of the virial theorem and found that, within observational constraints for the Pigtail,

[TABLE]

Thus, the the confinement of magnetic fields dominates external pressure, which dominates self-gravity. We conclude that self-gravity is relatively unimportant in the equilibrium structure of molecular tornadoes. The relatively large values of $\mathcal{M/|W|}$ and $|\mathcal{K}_{P}|/|\mathcal{W}|$ are due to the weakness of self-gravity. Since $\mathcal{M}/|\mathcal{W}|<0$, the magnetic field is dominated by the toroidal component, which lends support for the proposed magnetic-squeezing/twisting mechanism of molecular tornadoes. 5. 5.

We performed the same analysis for the DHN and GCT as we did for the Pigtail, but the conclusions we draw should be taken with more caution for these objects. Our analysis seems to suggest, however, that the GCT is dominated by the toroidal magnetic stress, whereas the DHN may be dominated by the pressure and/or magnetic stress. Self-gravity is relatively unimportant in these molecular tornadoes as well.

6 Acknowledgments

K.A. would like to acknowledge Tomoharu Oka, Mark Morris, and Yoshiaki Sofue for their insights and interpretations on the Pigtail Molecular Cloud, DHN, and GCT, respectively. J.D.F. acknowledges the support of a Discovery Grant from the National Sciences and Engineering Research Council of Canada.

Appendix A Torsional Alfvén Wave Condition via Small Perturbations

We assume a torsional Alfvén wave of the form ${B}_{\phi}=B_{\phi,0}e^{-i(kz-\omega t)}$ and ${v}_{\phi}=v_{\phi,0}e^{-i(kz-\omega t)}$, where $B_{\phi,0}$ and $v_{\phi,0}$ are the corresponding amplitudes, $k$ is the wave number, and $\omega$ is the angular frequency of the perturbation. Notice that the differential operators $\partial_{t}$ and $\bm{\nabla}$ become $\partial_{t}\rightarrow i\omega$ and $\bm{\nabla}\rightarrow-ik\hat{z}$ when they act on the perturbed terms. By linearizing Faraday’s Law, it becomes

[TABLE]

Since the perturbations are small, it can be shown from the wave equation (equation 25) that

[TABLE]

so equation A1 can be simplified to

[TABLE]

which is identical to equation 29. Thus, the torsional Alfvén wave condition (equation 30) follows.