{'fraction_non_alphanumeric': 0.03603682136755117, 'fraction_numerical': 0.019613654590466892, 'mean_word_length': 4.505230828294462, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetohydrodynamic Models of Molecular Tornadoes

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.


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 ∼ 15 pc pitch that traces at least two rounds of a helix in 12 CO, extending to an angular size of ∼ 0.15 • × 0.15 • and corresponding to a spatial size of ∼ 20 × 20 pc 2 with a mass of (2 − 6) × 10 5 M (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/ 13 CO 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 ∼ 30 km s −1 with an estimated mass of M ∼ 1.2 × 10 6 M . 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 ∼ 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 φ (r) to the rotation law Ω(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.

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 where ρ is the density, v is the velocity, P is the pressure, B is the magnetic field, Φ is the gravitational potential, and σ is the velocity dispersion. We consider the velocity dispersion to contain both thermal and nonthermal contributions: σ = σ 2 thermal + σ 2 nonthermal . 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: and In these scaling relations, G is the universal gravitational constant, g is the gravitational field, m is the linear mass, v φ is the rotational speed, and Ω 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 σ, 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, φ, z) and rotational motion only, the velocity and magnetic field take the form and Using the assumption of equation 3, the only non-trivial terms of equation 1 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 and respectively. The flux-to-mass ratios are free functions of the magnetic flux Φ M in the context of perfect MHD due to Ferraro's law of isorotation (Ferraro 1937), but flux Φ M (r) is a function of r in cylindrical symmetry. Fiege & Pudritz (2000a) assumed that both Γ z and Γ φ were constant. However, we consider the possibility that these functions might vary as a function of Φ 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) := −∂Φ(r)/∂r, is described by Poisson's equation for gravity , which becomes Having defined g, equation 4 can be expressed as 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.

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 Our analytical solution requires that the centrifugal force balances exactly with the inward toroidal magnetic stress, so that We further assume that the ratio of pressure and total magnetic pressure is constant: We define and where β z and β z are also constant, so that equation 10 can be written as Using the definition of β φ in equation 9, and the definitions of σ from equation 1d, it is straight- Equation 13 implies that the velocity field is highly sheared at r = 0, which violates the usual boundary condition v φ (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 φ 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 P tot = P + P mag,z + P mag,φ = P 1 + 1 Furthermore, we define an effective sound speed given by where dimensional units are restored for this equation only, so that we can redefine the dimensionless gravitational potentialΦ 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 which can be integrated to find a solution for ρ so that Poisson's equation can be written as where r 0 is redefined as 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: and where 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.

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.

Angular Frequency as a Power Law
We consider a rotation law where Ω scales with the radius r according to a power law with index α: Keplerian rotation and solid body rotation are considered to be limiting cases. We show below that this assumption constrains α to the range −1 ≤ α ≤ 0, where α = 0 is the case for solid body rotation. For clarity, we note that the α = −1 index would be more accurately described as "Keplerian-like" where gas orbits a line mass in the absence of pressure gradients.

Keplerian Rotation
Poisson's equation and the divergence theorem are used to find the following relationship for g in terms of r: 2πgrL = 4π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 Thus we find that Keplerian rotation is just a special case of the radial power law (equation 21) with α = −1.

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 Φ M , so that v φ = Ω(Φ M )r. Assuming Γ z is constant, we find that ∂Φ M /∂r = 2πΓ z ρr. The angular frequency Ω is solved numerically from ∂Φ M /∂r along with other quantities mentioned later in Section 2.5. The rotational frequency is then assumed to follow a power law The lower limit on α M is determined by the proportionality relations: For the sake of interest, we also analyze α M ≥ −1 to understand how α M affects the behaviour of the radial density.

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 and The MHD equations (equation 1) yield which is the general wave equation describing transverse perturbations propagating along the flux tube at the Alfvén speed (Jackson 1975) 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 ρ 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 φ dt. As the wave propagates along the flux tube, a toroidal magnetic field B φ is generated. Consequently, Faraday's Law reduces to 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 φ = 0. Below the interface, the magnetic field has been perturbed and exhibits a toroidal component B φ , and rotates at v φ . 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 and follow through with the integration. Thus the field angle θ is perturbed according to 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 where Γ φ = Γ φ (r), in general.

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 where script characters denote quantities per-unit-length; thus K is the total kinetic energy per unit length, M is the magnetic energy per unit length, and W is the gravitational energy per unit length. The total kinetic energy contains the rotational kinetic energy K rot , surface pressure term K P , and internal turbulence K σ so that The rotational kinetic energy is given by K rot = (1/2)IΩ 2 , where 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 K P = −(3/2)P S V where V = πr 2 S is the cross-sectional area of the filament. The energy due to internal turbulence is given by K σ = (3/2)mσ 2 . Gravitational energy is described by 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 |W| yields 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 M/|W|, which determines the dominating component of the magnetic field (Fiege & Pudritz 2000a). For M/|W| < 0, the toroidal component is dominant in pinching the filament, while M/|W| > 0 implies net support by the poloidal field. The numerical results in Section 3 are examined with the help of equation 33. 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.

Guidelines for Stability
In the case where there is no B z external to the cylinder, the condition for instability is (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 and respectively. Equivalently, the stability conditions above can be expressed in terms of the rotational period P = 2π/Ω S , and an Alfvén crossing time τ A := λ/v A that describes the time it takes for an Alfvén wave to propagate across one wavelength or pitch λ. The stability conditions in equation 35 for m = 0, 1, and 2, respectively, are expressed as and 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 /λ, as expressed in equation 36. For a given v φ /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 Ω S > √ 2v A /r S , according to equation 35b. In other words, a molecular tornado observed with a kink instability must have been spinning at Ω S > √ 2v A /r S for the instability to trigger.

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 ξ defined by ξ ≡ ln(r). The radial differential operators therefore transform according to and The equations for g(r) = −∂Φ(r)/∂r, dg/dr (equation 6), and dρ/dr (equation 7) are trivially transformed to equivalent equations using ξ as the independent variable, which is what we solve numerically.
In the case of equation 30, where B φ is due to a torsional Alfvén wave, it is useful to further simplify Equation 7. If Γ z is constant, then equation 7 can be expressed explicitly in terms of the dimensionless quantities Ω, and dΩ/dξ by using the flux-to-mass ratios (equation 5) along with the expression above for Ω (equation 30). Algebraic simplification results in the expression where we have used equation 30 to write Γ φ in terms of Ω.

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 Ω 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:r S ,P S ,m S ,Ω S , and ρ 0 . The readily observable (or at least somewhat observationally constrained) quantities are σ, r S , m S , and P S . The central density ρ 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 ρ 0 from the equations. This reduces the scaling laws forr S ,P S ,m S ,Ω S , and ρ 0 to three equations that describe combinations of the dimensionless quantities on the left side, in terms of observables on the right:r 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 10 7 K cm −3 , and can possibly be as high as ∼ 10 8 K cm −3 . The velocity dispersion for molecular clouds in the CMZ are roughly σ ∼ 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.

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, namelyr S , P S ,m S ,Ω S , and σ. 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 σ = 13.2 km s −1 , and n(H 2 ) = 10 3.5±0.25 cm −3 (σ 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±0.05, we estimate via the isothermal equation of state that P S /k = (2 ± 1) × 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 S , are related by the pitch angle via r S = r S cos(θ). 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 ∼ 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 0.61 r 2 SP S 7.32, There is more uncertainty in the physical parameters of the DHN and GCT, however. In particular, σ, Ω S , and P S are unclear. Generally, molecular clouds of the CMZ exhibit σ ∼ 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 P S 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±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 ∼ −35 km s −1 and 0 km s −1 having lengths ∼ 150 pc, and masses of 0.8×10 4 M and 3.3 × 10 4 M , respectively 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 ∼ 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 2 × 10 −4 r 2 0.08 m S 0.32, and 80 Ω 2 S P S 3 × 10 5 .
Assuming σ = 15 km s −1 , and provided that the GCT has a density of ρ ∼ 2.7 × 10 −21 g cm −3 (Sofue 2007), the isothermal equation of state suggests an internal pressure of 4.4 × 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 0.11 r 2 SP S 8.71, 1.48 m S 12.04, and 1.70 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).

RESULTS
Results of numerically solving equations 7 and 38, along with dΦ/dξ and dg/dξ , and the rotation laws described in Section 2.2 are presented in this section. The density ρ 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 ρ 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, nonmagnetic, self-gravitating cylinders as a benchmark. At large radii, the Ostriker solution behaves as ρ ∼ 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, Ω 0 , Γ z , and Γ φ . 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 ≤ r ≤ r max , where the lower and upper limits to the dimensionless radius are r min = 10 −6 and r max = 10 6 , respectively. If Ω follows the power law rotation given by equation 21, then Ω 0 , Γ z , and Γ φ are free parameters. In this case, Ω 0 , Γ z , and Γ φ are each randomly chosen over the intervals 0 ≤ Ω 0 ≤ 12, (43a) and The limits of equation 43 are conservatively chosen, with the limits of Γ φ and Γ z extended beyond that of Fiege & Pudritz (2000a) to account for the environment of the Galactic Centre. Each random set of Ω 0 , Γ z , and Γ φ define a particular untruncated filament model. Otherwise if Ω follows a power law of magnetic flux (equation 23), then Γ z , and Γ φ 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 Figure 1. Density profiles for filaments obeying the rotation law Ω = Ω 0 (r/r 0 ) α and the torsional Alfvén wave constraint. Green, dashed lines are theoretical density profiles. The Ostriker solution is the solid blue line for all plots, henceforth. Notice that as the power law index becomes more negative, solutions that behave as ρ ∼ r −4 or steeper occur more frequently. In each plot, there are 50 density profiles. It should be noted that each plot is magnified for clarity and the integration is from r min ≤ r ≤ r max . Figure 2. Density profiles for filaments obeying the rotation law Ω = Ω 0 Φ α M M and the torsional Alfvén wave constraint Green, dashed lines are theoretical density profiles. Notice that as the power law index becomes more negative, solutions can become steeper. In each plot, there are 50 density profiles. It should be noted that each plot is magnified for clarity and the integration is from r min ≤ r ≤ r max . The inset has the same axis quantities as the host plot, and show the boundary condition ρ(r min )/ρ 0 = 1. The extreme density inversions shown in the inset are indeed part of the solutions, which we note are unphysical density inversions, and therefore excluded as feasible models. These extreme density inversions are due to Ω behaving as r −1 or steeper, so there is extreme rotation near the origin. The α M = −0.75 plot is included in this figure for the sake of interest only. and 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ρ/dr > 0), as recommended by Recchi et al. (2014). This condition also limits the density so that ρ(r) ≤ ρ 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   Constant parameters are set to Γ z = 1, Γ φ = π/50, and Ω 0 = π/25. Changing parameters are ordered as the dashed blue, green, red, and black lines, respectively: (Left) Γ φ = π/100, π/2, 2π, 4π; (Centre) Ω 0 = π/25, π/50, π/75, π/100; (Right) Γ z = π/2, π, 2π, 3π. It should be noted that each plot is magnified for clarity and the integration is indeed from r min ≤ r ≤ r max . Notice that as Γ φ increases, the density profile becomes more pinched, and the inversions smooth out. As Ω 0 decreases, the inversions begin at lower densities. A larger Γ z causes the inversions to increase in amplitude. Such density inversions are unrealistic, and the density profiles shown here are not truncated for the purpose of illustration only.
frequency, constant flux-to-mass ratios, and excluding the torsional Alfvén wave condition, changing Γ z changes the asymptotic density value that is reached at large radii only by a small amount. This asymptotic value of ρ at large radii can be approximated from equation 7, which reveals that 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 (ρ, P S , Ω S , Γ φ , Γ 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 (ρ, P S , Ω S , Γ φ , Γ 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 Figure 6. Sample of changing model parameters excluding self-gravity and the torsional Alfvén wave condition. All plots have parameters set such that the dashed line colours corresponding to parameters that change linearly (step size of 4) in the range Ω 0 = [0, 3], Γ z = [10 −3 , 25], and Γ φ = [0, 25] . Green lines correspond to changes in Γ z , red lines to Γ φ , and black lines to Ω 0 . If a parameter is constant then they are set to: Γ z = 1, Γ φ = π/50, and Ω 0 = π/25. For a set of coloured profiles, as Γ z is increasing, the density profiles extend radially. As Γ φ is increasing, the density profiles pinch (decrease in density overall). As Ω 0 increases, the overall density increases. It should be noted that each plot is magnified for clarity and the integration is indeed from r min ≤ r ≤ r max . The boundary condition ρ(r min )/ρ 0 = 1 is always met. In the case of α M = −0.75, ρ increases rapidly at r/r 0 << 1 due to unrealistic, extreme rotation near the origin, which is outside of the plot frame. All profiles which exceed ρ 0 are rejected from analysis and are considered unphysical models. These density profiles are not truncated for the purpose of illustration only. 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 ρ 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- Figure 7. Relative energies of observationally constrained models including self-gravity, power law rotation (Ω = Ω 0 (r/r 0 ) −0.25 ), and the torsional Alfvén wave condition. Observational constraints are based on measurements of the Pigtail. The energy terms are taken from equation 33. Notice that the pressure and magnetic stresses, in particular, play an important role in the equilibrium filament structure. The inset has the same axis labels as the host plot.
to-gravitational energy ratio is and the ratio of the surface pressure term to gravitational energy is Comparing the surface pressure term to the magnetic energy reveals that |M| > |K P | 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.

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.

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, ρ often tends to increase with radius and may exceed ρ 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 Γ φ /Ω 0 ratio tends to keep ρ < ρ 0 and dρ/dr < 0. The effect of Γ z does not seem to be as significant compared to Γ φ and Ω 0 , unless Γ 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 α M seems to have a dominating effect on the behaviour of ρ(r). There are density profiles where ρ > ρ 0 , regardless of the value of α M . A larger portion of these models exceed ρ 0 as α M becomes steeper. All of these models, except where Γ φ = 0, decrease in density as r becomes large. The Γ φ /Ω 0 ratio has a similar effect as in the models that follow equation 21. See Figure 6 for the associated density profiles.

Density Inversions
Density inversions are found in some models where ρ(r) < ρ 0 (for all r), but curiously only if selfgravity 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 Γ φ is weak.
It is easier to see the effects due to Ω 0 = Ω(r), Γ φ , and Γ z , by first assuming that they are constant. The effect of increasing Γ φ is to dampen the density inversions (see Figure 5). As Ω 0 increases, the entire density profile along with the inversions, shifts to greater densities, and the inversion amplitude decreases slightly. The effect of increasing Γ 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 Γ φ seems to have a more noticeable effect on the behaviour of ρ than Γ z . Generally, Γ φ and Γ z behave as Fiege & Pudritz (2000a) describe -Γ z tends to support/radially extend the filament compared to the Ostriker solution, while Γ φ 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 ≤ α < 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.

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 Γ z /Ω 0 ratio tends to produce steeper solutions initially, but then achieve an asymptotically constant (dρ/dr ≈ 0) value of ρ at large radii. As α decreases, asymptotic solutions can be seen with smaller Γ z /Ω 0 ratios. Non-self-gravitating filaments following the same rotation law generally have solutions that are not asymptotic. As α becomes progressively steeper, the density profiles tend to become more shallow. At α = −1, all profiles are truly asymptotically constant at ρ = ρ 0 . All profiles decrease when α > −1 .
Models that are self-gravitating and rotate according to the flux power law (equation 23) have decreasing density when α M = −0.25, but more negative values have profiles that were deemed unrealistic because ρ > ρ 0 at some point(s) (see Figure 2). All profiles begin to decrease at a sufficiently large radius, but the density inversion, where dρ/dr > 0, becomes more pronounced as α 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 α 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 Γ z /Ω 0 ratio. The ratio of Γ z /Ω 0 does not need to be as large to observe asymptotic solutions as α M decreases. Non-self-gravitating filaments following the same rotation law do not have asymptotic solutions at α M = −0.25, and the density profiles decrease with the radius. Density profiles with α M −0.50 were deemed unrealistic due to density inversions that increase the density above ρ 0 at some point(s). This α M ≈ −0.50 boundary coincides with the limit of α M that was stated in Section 2.2.3. More negative values of α M cause Ω 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 ρ 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 ρ(r) ≤ ρ 0 (or are at least asymptotically constant at ρ = ρ 0 when α = −1 for non-self-gravitating models), in both self-gravitating and non-self-gravitating cases regardless of any α < 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 ρ(r) ≤ ρ 0 in both self-gravitating and non-self-gravitating cases, as long as α 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, M/|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), W is compared to M and 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 α = −0.25. The magnetic energy ratio with gravitational energy within observational constraints is found to be −10 5 M/|W| 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 |K P |/|W| 10 3 , which further indicates that these molecular tornadoes are only weakly self-gravitating. It is also found that |M| > |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 W ≈ 0 may be justified. If the toroidal magnetic stress is dominant, equation 33 may best be rewritten as where it is easier to see that if M < 0 then (K rot + K σ ) < −M/2, and if M > 0 then (K rot + K σ ) > −M/2. It was also found that |M| > |K P |, and since models follow M < 0, this means that −(1/2)|M| < (K rot + K σ ) < −(3/2)|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 |M| and |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 Ω 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 Ω 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 |K rot |, and would then require a much stronger |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.

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.
1. 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 Γ φ and Ω (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.
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).

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).
Thus, the the confinement of magnetic fields dominates external pressure, which dominates selfgravity. We conclude that self-gravity is relatively unimportant in the equilibrium structure of molecular tornadoes. The relatively large values of M/|W| and |K P |/|W| are due to the weakness of self-gravity. Since M/|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. 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. We assume a torsional Alfvén wave of the form B φ = B φ,0 e −i(kz−ωt) and v φ = v φ,0 e −i(kz−ωt) , where B φ,0 and v φ,0 are the corresponding amplitudes, k is the wave number, and ω is the angular frequency of the perturbation. Notice that the differential operators ∂ t and ∇ become ∂ t → iω and ∇ → −ikẑ when they act on the perturbed terms. By linearizing Faraday's Law, it becomes Since the perturbations are small, it can be shown from the wave equation (equation 25) that so equation A1 can be simplified to which is identical to equation 29. Thus, the torsional Alfvén wave condition (equation 30) follows.