Energy Conservation in the Thin Layer Approximation: VI. Bubbles and Super-Bubbles

Lorenzo Zaninetti

Physics Department, via P. Giuria 1, Turin, Italy.

DOI: 10.4236/ijaa.2021.113017                

Abstract

We model the conservation of energy in the framework of the thin layer approximation for two types of interstellar medium (ISM). In particular, we analyse an ISM in the presence of self-gravity and a Gaussian ISM which produces an asymmetry in the advancing shell. The astrophysical targets to be simulated are the Fermi bubbles, the local bubble, and the W4 super-bubble. The theory of images is applied to a piriform curve, which allows deriving some analytical formulae for the observed intensity in the case of an optically thin medium.

Keywords

ISM: Bubbles, Galaxy: Disk

1. Introduction

We now summarize the first uses of some words: “super-shell” can be found in [1], where eleven H I objects are examined, “super-bubble” can be found in [2], where an X-ray region with a diameter of 450 pc connected with Cyg X-6 and Cyg X-7 is observed and “worms”, meaning gas filaments crawling away from the galactic plane in the inner Galaxy, can be found in [3]. Super-bubbles or super-shells can be defined as cavities with diameters greater than 100 pc and density of matter lower than that of the surrounding interstellar medium (ISM) [4]. Bubbles have smaller diameters, between 10 pc and 100 pc [5]. Some models which explain super-shells as being due to the combined explosions of supernova in a cluster of massive stars will now be reviewed. In semi-analytical calculations, the thin-shell approximation can be the key to obtaining the expansion of the super-bubble; for example, [5] [6] [7] [8] [9]. The Kompaneyets approximation, see [10] [11], has been used in order to model the super-bubble W4 [9] and the Orion-Eridanus super-bubble [12] [13]. The hydro-dynamical approximation, with the inclusion of interstellar density gradients, can produce a blowout into the galactic halo, see [14] [15]. Recent Planck 353-GHz polarization observations allow mapping the magnetic field, see [16] for the Orion-Eridanus super-bubble, and we recall that the expansion of super-bubbles in the presence of magnetic fields has been implemented in various magneto-hydrodynamic codes, see [17] [18]. The present paper derives the equation of motion for two different ISMs in the framework of the energy conservation for the thin layer approximation, see Section 2; compares the observed and the theoretical sections for Fermi bubbles, the local bubble, and the W4 super-bubble, see Section 3; and derives a new analytical formula for the theoretical profile in intensity using the piriform curve, see Section 4.

2. The Equations of Motion

We start with the conservation of kinetic energy in spherical coordinates in the framework of the thin layer approximation

$\frac{1}{2}{M}_0\left(r_0\right)v_0^2=\frac{1}{2}M\left(r\right)v^2\mathrm{,}$ (1)

where $M_0\left(r_0\right)$ and $M\left(r\right)$ are the swept masses at r₀ and r, while v₀ and v are the velocities of the thin layer at r₀ and r. The above equation holds for the solid angle $\Delta \Omega$, which in the following is unity. We now present two asymmetric equations of motion for bubbles and super-bubbles. The above equation is a differential equation of the first order:

$\frac{1}{2}{M}_0\left(r_0\right)v_0^2=\frac{1}{2}M\left(r\right){\left(\frac{\text{d}r}{\text{d}t}\right)}^2\mathrm{.}$ (2)

The asymmetry is due to a gradient of the number of particles with the distance or galactic height, z, which is parametrized as

$n\left(z\right)=n_1{\text{e}}^{-z^2/H_1^2}+n_2{\text{e}}^{-z^2/H_2^2}+n_3{\text{e}}^{-\left|z\right|/H_3}\mathrm{.}$ (3)

where n₁ = 0.395 particles cm⁻³, H₁ = 127 pc, n₂ = 0.107 particles cm⁻³, H₂ = 318 pc, n₃ = 0.064 particles cm⁻³, and H₃ = 403 pc [19] [20] [21]. In the framework of Cartesian coordinates, $\left(x\mathrm{,}y\mathrm{,}z\right)$, when the explosion starts at $\left(\mathrm{0,0,0}\right)$ we have an up-down symmetry, $r\left(x\mathrm{,}y\mathrm{,}-z\right)=r\left(x\mathrm{,}y\mathrm{,}z\right)$ and a right-left symmetry $r\left(x,-y,z\right)=r\left(x,y,z\right)$. Conversely, when the explosion starts at $\left(\mathrm{0,0,}{z}_{\text{OB}}\right)$, where z_OB represents the distance in pc from the position of the OB association which generate the phenomena, we have only a right-left symmetry $r\left(x,-y,z\right)=r\left(x,y,z\right)$.

2.1. Numerical Methods

In the absence of an analytical solution for the trajectory, we outline four ways which allow obtaining a numerical solution.

1) Evaluation of the numerical solution with the the Runge-Kutta method.

2) A non-linear method which obtains the trajectory by the following non-linear equation

${\displaystyle {\int }_{r_0}^r}\frac{1}{\frac{\text{d}r}{\text{d}t}}\text{d}r=t-t_0\mathrm{.}$ (4)

3) The Euler method, which solves the following recursive equations

$r_{n+1}=r_n+v_n\Delta t$ (5a)

$v_{n+1}=v_n{\left(\frac{M_n\left(r_n\right)}{M_{n+1}\left(r_{n+1}\right)}\right)}^{1/2}\mathrm{,}$ (5b)

where $r_n$, $v_n$, and $M_n$ are the temporary radius, velocity, and total mass, respectively, $\Delta t$ is the time step, and n is the index.

4) A power series solution of the form

$\begin{array}{l}r\left(t\right)\\ =a_0+a_1\left(t-t_0\right)+a_2{\left(t-t_0\right)}^2+a_3{\left(t-t_0\right)}^3+\cdots \mathrm{,}\end{array}$ (6)

see [22] [23].

The case of an expansion that starts from a given galactic height z, denoted by z_OB, which represents the OB associations, is also analysed. The advancing expansion is computed in a 3D Cartesian coordinate system $\left(x\mathrm{,}y\mathrm{,}z\right)$ with the centre of the explosion at $\left(0,0,0\right)$. The explosion is better visualized in a 3D Cartesian coordinate system $\left(X\mathrm{,}Y\mathrm{,}Z\right)$ in which the galactic plane is given by $Z=0$. The following translation, T_OB, relates the two Cartesian coordinate systems

$T_{\text{OB}}\mathrm{<mi>\; </mi>}\{\begin{array}{l}X=x\hfill \\ Y=y\hfill \\ Z=z+z_{\text{OB}}\hfill \end{array}$ (7)

where z_OB is the distance in pc of the OB associations from the galactic plane. In the case of $z_{\text{OB}}\ne 0$, the two masses which appear in Equation (5b) should be carefully evaluated.

2.2. Medium in the Presence of Self-Gravity

We assume that the number density distribution scales as

$n\left(z\right)=n_0{\text{sech}}^2\left(\frac{z}{2h}\right),$ (8)

where $n_0$ is the density at $z=0$, h is a scaling parameter, and sech is the hyperbolic secant [24] [25] [26] [27]. In order to include the boundary conditions we assume that the density of the medium around the OB associations scales with the self-gravity piece-wise dependence

$\rho \left(r;r_0\right)=\{\begin{array}{ll}{\rho }_c\hfill & \text{if}r\le r_0\hfill \\ {\rho }_c{\text{sech}}^2\left(\frac{z}{2h}\right)\hfill & \text{if}r>r_0\hfill \end{array}$ (9)

where ${\rho }_c$ is the density at $z=0$. In order to find an acceptable value of h, we make a comparison with Equation (3), after which we choose $h=90\text{pc}$, see Figure 1.

The mass M₀ swept in the interval $\left[0,r_0\right]$ is

$M_0\left({\rho }_c\mathrm{,}{r}_0\right)=\frac{4}{3}{\rho }_c\pi r_0^3\mathrm{.}$

The total mass $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}h\right)$ swept in the interval [0, r] is

$\begin{array}{l}M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}h\right)\\ =\frac{{\rho }_c{r}_0^3}{3}-4\frac{{\rho }_c{r}^{2}h}{\cos \left(\theta \right)}{\left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)}^{-1}+4\frac{{\rho }_c{r}^{2}h}{\cos \left(\theta \right)}-8\frac{{\rho }_c{h}^{2}r}{{\left(\cos \left(\theta \right)\right)}^2}\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\\ -8\frac{{\rho }_c{h}^3}{{\left(\cos \left(\theta \right)\right)}^3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)+4\frac{{\rho }_c{r}_0^{2}h}{\cos \left(\theta \right)}{\left(1+{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)}^{-1}\\ -4\frac{{\rho }_c{r}_0^{2}h}{\cos \left(\theta \right)}+8\frac{{\rho }_c{h}^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^2}\ln \left(1+{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)+8\frac{{\rho }_c{h}^3}{{\left(\cos \left(\theta \right)\right)}^3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)\end{array}$ (10)

where $\theta$ is the polar angle and the polylog operator is defined by

$polylog\left(s\mathrm{,}z\right)={\text{Li}}_s\left(z\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}\frac{z^n}{n^s}$ (11)

where ${\text{Li}}_s\left(z\right)$ is the Dirichlet series. The positive solution of Equation (1) gives the velocity as a function of the radius:

$v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}h\right)=\frac{AN}{AD}\mathrm{,}$ (12)

where

Figure 1. Profiles of density versus scale height z: the disk in presence of self-gravity as given by Equation (8) when $h=90\text{pc}$ (dashed blue line) and the three-component exponential distribution as given by Equation (3) (red full line).

$\begin{array}{c}AN=-(-24\cos \left(\theta \right)r_0\left({\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1\right)(-h^3({\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}\\ +{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1)polylog\left(\mathrm{2,}-{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\\ +h^3\left({\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1\right)polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\\ +\cos \left(\theta \right)((-h^2{r}_0\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)+h^{2}r\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\\ \begin{array}{c}{}_{}^{}\\ {}_{}^{}\end{array}-1/2\cos \left(\theta \right)\left(1/12r_0^3\cos \left(\theta \right)+h{r}^2-h{r}_0^2\right)){\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}\end{array}$

$\begin{array}{c}-h^2{r}_0\left({\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1\right)\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\\ +h^{2}r\left({\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1\right)\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\\ -1/2\cos \left(\theta \right)(\left(1/12r_0^3\cos \left(\theta \right)-h{r}_0^2\right){\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\\ +{\left(1/12r_0^3\cos \left(\theta \right)+h{r}^2\right){\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}+1/12r_0^3\cos \left(\theta \right)))))}^{1/2}\cos \left(\theta \right)v_0{r}_0\end{array}$ (13)

and

$\begin{array}{l}AD=-r_0^3{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}{\left(\cos \left(\theta \right)\right)}^3+24h^{2}r{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\cos \left(\theta \right)\\ -24h^2{r}_0{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\cos \left(\theta \right)-12h{r}^2{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}{\left(\cos \left(\theta \right)\right)}^2\\ +12h{r}_0^2{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}{\left(\cos \left(\theta \right)\right)}^2-r_0^3{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}{\left(\cos \left(\theta \right)\right)}^3-r_0^3{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}{\left(\cos \left(\theta \right)\right)}^3\\ +24h^3{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)-24h^3{\text{e}}^{\frac{\cos \left(\theta \right)\left(r+r_0\right)}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\end{array}$

$\begin{array}{c}+24h^{2}r{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\cos \left(\theta \right)+24h^{2}r{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\cos \left(\theta \right)\\ -24h^2{r}_0{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\cos \left(\theta \right)-24h^2{r}_0{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\cos \left(\theta \right)\\ -12h{r}^2{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}{\left(\cos \left(\theta \right)\right)}^2+12h{r}_0^2{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}{\left(\cos \left(\theta \right)\right)}^2-r_0^3{\left(\cos \left(\theta \right)\right)}^3\\ +24h^3{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)-24h^3{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\end{array}$

$\begin{array}{c}+24h^3{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)-24h^3{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\\ +24h^{2}r\ln \left(1+{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)\cos \left(\theta \right)-24h^2{r}_0\ln \left(1+{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\cos \left(\theta \right)\\ +24h^{3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\cos \left(\theta \right)}{h}}\right)-24h^{3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{\cos \left(\theta \right)r_0}{h}}\right)\end{array}$ (14)

The differential equation which governs the motion for the medium in the presence of self-gravity is

$\begin{array}{l}(\frac{{\rho }_c{r}_0^3}{3}-4\frac{{\rho }_c{\left(r\left(t\right)\right)}^{2}h}{\cos \left(\theta \right)}{\left(1+{\text{e}}^{\frac{r\left(t\right)\cos \left(\theta \right)}{h}}\right)}^{-1}+4\frac{{\rho }_c{\left(r\left(t\right)\right)}^{2}h}{\cos \left(\theta \right)}\\ -8\frac{{\rho }_c{h}^{2}r\left(t\right)}{{\left(\cos \left(\theta \right)\right)}^2}\ln \left(1+{\text{e}}^{\frac{r\left(t\right)\cos \left(\theta \right)}{h}}\right)-8\frac{{\rho }_c{h}^3}{{\left(\cos \left(\theta \right)\right)}^3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r\left(t\right)\cos \left(\theta \right)}{h}}\right)\\ +4\frac{{\rho }_c{r}_0^{2}h}{\cos \left(\theta \right)}{\left(1+{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)}^{-1}-4\frac{{\rho }_c{r}_0^{2}h}{\cos \left(\theta \right)}+8\frac{{\rho }_c{h}^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^2}\ln \left(1+{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)\\ +8\frac{{\rho }_c{h}^3}{{\left(\cos \left(\theta \right)\right)}^3}polylog\left(\mathrm{2,}-{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2-\frac{{\rho }_c{r}_0^3{v}_0^2}{3}=\mathrm{0,}\end{array}$ (15)

and does not have an analytical solution. Figure 2 shows the numerical solution obtained with the Runge-Kutta method.

A Taylor expansion of order 3 of Equation (15) gives

$\begin{array}{l}r\left(t\mathrm{<mo>\; ;\; </mo>}{t}_0\mathrm{,}{v}_0\mathrm{,}{r}_0\mathrm{,}h\right)\\ =r_0+v_0\left(t-t_0\right)-3\frac{v_0^2{\left(t-t_0\right)}^2}{r_0}{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}{\left({\left({\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}\right)}^2+2{\text{e}}^{\frac{r_0\cos \left(\theta \right)}{h}}+1\right)}^{-1}\mathrm{,}\end{array}$ (16)

Figure 2. Phase plane portrait and numerical solution (full yellow line) for a medium in the presence of self-gravity as given by Equation (8) when $r_0=250\text{pc}$, $h=90\text{pc}$, $t=4\times {10}^6\text{yr}$, $t_0=4\times {10}^4\text{yr}$ and $v_0=2000\text{km}\cdot {\text{s}}^{-1}$.

and Figure 3 shows the numerical solution obtained by the Runge-Kutta method and the series solution up to a time for which the percentage error is less than 10%.

2.3. Gaussian Medium

We assume that the number density distribution scales as

$n\left(z\right)=n_0{\text{e}}^{-\frac{z^2}{z_0^2}}\mathrm{,}$ (17)

where $n_0$ is the density at $z=0$ and $z_0$ is a scaling parameter. We now give the adopted piece-wise dependence for the Gaussian medium

$\rho \left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{z}_0\right)=\{\begin{array}{ll}{\rho }_c\hfill & \text{if}r\le r_0\hfill \\ {\rho }_c{\text{e}}^{-\frac{r^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\hfill & \text{if}r>r_0\hfill \end{array}$ (18)

where ${\rho }_c$ is the density at $z=0$. A comparison with Equation (3) gives $z_0=200\text{pc}$, see Figure 4. The total mass $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}{z}_0\right)$ swept in the interval [0, r] is

$M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}{z}_0\right)=\frac{BN}{12{\left(\cos \left(\theta \right)\right)}^3}\mathrm{,}$ (19)

where

$\begin{array}{c}BN=-(-4r_0^3{\left(\cos \left(\theta \right)\right)}^3+3\text{erf}\left(\frac{r_0\cos \left(\theta \right)}{z_0}\right)\sqrt{\pi }{z}_0^3-3\text{erf}\left(\frac{r\cos \left(\theta \right)}{z_0}\right)\sqrt{\pi }{z}_0^3\begin{array}{c}\\ \\ \end{array}\\ +6{\text{e}}^{-\frac{r^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r{z}_0^2-6{\text{e}}^{-\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r_0{z}_0^2){\rho }_c\mathrm{,}\end{array}$ (20)

and erf(x) [28] is the error function defined by

Figure 3. Numerical solution (red full line) and Taylor expansion of the solution (dashed blue line), parameters as in Figure 2 but $t=2.8\times {10}^5\text{yr}$.

Figure 4. Profiles of density versus scale height z: the medium is that in the presence of a Gaussian medium as given by Equation (29) when $h=90\text{pc}$ (dashed blue line) and the three-component exponential distribution as given by Equation (3) (red full line).

$\text{erf}\left(x\right)=2\frac{{\displaystyle {\int }_0^x}\exp \left(-t^2\right)\text{d}t}{\sqrt{\pi }}\mathrm{.}$ (21)

The velocity as a function of the radius is

$v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{z}_0\mathrm{,}{v}_0\right)=\frac{CN}{CD}\mathrm{,}$ (22)

where

$\begin{array}{c}CN=2(-6\cos \left(\theta \right)(-2/3r_0^3{\left(\cos \left(\theta \right)\right)}^3+{\text{e}}^{-\frac{r^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r{z}_0^2\\ -{\text{e}}^{-\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r_0{z}_0^2-1/2\text{erf}\left(\frac{\cos \left(\theta \right)r}{z_0}\right)\sqrt{\pi }{z}_0^3\\ {\begin{array}{c}\\ \\ \end{array}+1/2\text{erf}\left(\frac{r_0\cos \left(\theta \right)}{z_0}\right)\sqrt{\pi }{z}_0^3)r_0)}^{1/2}\cos \left(\theta \right)v_0{r}_0\mathrm{,}\end{array}$ (23)

and

$\begin{array}{c}CD=4r_0^3{\left(\cos \left(\theta \right)\right)}^3-6{\text{e}}^{-\frac{r^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r{z}_0^2+6{\text{e}}^{-\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\cos \left(\theta \right)r_0{z}_0^2\\ -3\text{erf}\left(\frac{r_0\cos \left(\theta \right)}{z_0}\right)\sqrt{\pi }{z}_0^3+3\text{erf}\left(\frac{\cos \left(\theta \right)r}{z_0}\right)\sqrt{\pi }{z}_0^3\mathrm{.}\end{array}$ (24)

The differential equation which governs the motion for the Gaussian medium is

$\begin{array}{l}4{\left(\cos \left(\theta \right)\right)}^3{\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{r}_0^3-4r_0^3{v}_0^2{\left(\cos \left(\theta \right)\right)}^3\\ -6r\left(t\right)\cos \left(\theta \right){\text{e}}^{-\frac{{\left(r\left(t\right)\right)}^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}{\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{z}_0^2\end{array}$

$\begin{array}{l}+6\cos \left(\theta \right){\text{e}}^{-\frac{r_0{}^2{\left(\cos \left(\theta \right)\right)}^2}{z_0{}^2}}{\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{r}_0{z}_0^2\\ -3\sqrt{\pi }\text{erf}\left(\frac{r_0\cos \left(\theta \right)}{z_0}\right){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{z}_0^3\\ +3\sqrt{\pi }\text{erf}\left(\frac{r\left(t\right)\cos \left(\theta \right)}{z_0}\right){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{z}_0^3=0.\end{array}$ (25)

Figure 5 shows the numerical solution obtained with the Runge-Kutta method.

A Taylor expansion of order 3 of Equation (25) gives

$r\left(t\mathrm{<mo>\; ;\; </mo>}{t}_0\mathrm{,}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\right)=r_0+v_0\left(t-t_0\right)-\frac{3v_0^2{\left(t-t_0\right)}^2}{4r_0}{\text{e}}^{-\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{z_0^2}}\mathrm{,}$ (26)

and Figure 6 gives the numerical solution obtained by the Runge-Kutta method

Figure 5. Phase plane portrait and numerical solution (full yellow line) in the Gaussian case when $r_0=250\text{pc}$, $z_0=90\text{pc}$, $t=7.1\times {10}^5\text{yr}$, $t_0=7.1\times {10}^3\text{yr}$ and $v_0=10000\text{km}\cdot {\text{s}}^{-1}$.

Figure 6. Numerical solution (red full line) and Taylor expansion of the solution (dashed blue line), parameters as in Figure 5 but $t=4.26\times {10}^4\text{yr}$.

and the series solution up to a time for which the percentage error is less than 9%.

3. Astrophysical Applications

In the following we will analyse the local bubble, the Fermi bubble and the super bubble W4. An observational percentage reliability, ${\epsilon }_{\text{obs}}$, is introduced over the whole range of the polar angle $\theta$,

${\epsilon }_{\text{obs}}=100\left(1-\frac{{\displaystyle {\sum }_j{\left|r_{\text{obs}}-r_{\text{num}}\right|}_j}}{{\displaystyle {\sum }_j{r}_{\text{obs}}{}_{\mathrm{,}j}}}\right)\mathrm{,}$ (27)

where $r_{\text{num}}$ is the theoretical radius of the considered bubble, $r_{\text{obs}}$ is the observed radius of the considered bubble, and the index j varies from 1 to the number of available observations. The observational percentage of reliability allows us to fix the theoretical parameters.

3.1. The Local Bubble

The local bubble (LB) has already been simulated in the framework of the conservation of momentum [29]; here we adopt the framework of the conservation of energy. The numerical solution is shown as a cut in the x-z plane: see Figure 7 for a medium in the presence of self-gravity as given by Equation (9) and Figure 8 for a Gaussian density profile as given by Equation (18).

Figure 7. Geometrical section of the LB in the x-z plane with a profile in the presence of self-gravity as given by Equation (8) (green points) and observed profile (red stars). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $h=3.5\text{pc}$, $t=8.5\times {10}^4\text{yr}$, $t_0=8.5\times {10}^2\text{yr}$ and $z_{\text{OB}}=0$. The observational percentage reliability is ${\epsilon }_{\text{obs}}=\mathrm{82.42\; <mi>\; \%\; </mi>}$.

The 3D advancing surface of the local bubble for the case of self-gravity is shown in Figure 9.

3.2. The Fermi Bubble

Fermi bubbles have already been simulated in the framework of the conservation of momentum [30]; here we apply the conservation of energy. We now test our models on the image of the Fermi bubbles available at

Figure 8. Geometrical section of the LB in the x-z plane with a Gaussian profile (green points) and observed profile (red stars). The parameters are $v_0=4000\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $z_0=9.1\text{pc}$, $t=8.5\times {10}^4\text{yr}$, $t_0=8.5\times {10}^2\text{yr}$ and $z_{\text{OB}}=0$. The observational percentage reliability is ${\epsilon }_{\text{obs}}=\mathrm{82.58\; <mi>\; \%\; </mi>}$.

Figure 9. 3D surface of the LB with parameters as in Figure 7, with a profile in presence of self-gravity as given by Equation (8). The three Euler angles are $\Theta ={90}^{\circ }$, $\Phi =0^{\circ }$ and $\Psi ={90}^{\circ }$.

https://www.nasa.gov/mission_pages/GLAST/news/new-structure.html. The numerical solution is shown as a cut in the x-z plane: see Figure 10 for a density profile in the presence of self-gravity as given by Equation (9) and Figure 11 for

Figure 10. Geometrical section of the Fermi bubbles in the x-z plane with a profile in the presence of self-gravity as given by Equation (8) (green points) and observed profile (red stars). The parameters are $v_0=2000\text{km}\cdot {\text{s}}^{-1}$, $r_0=250\text{pc}$, $z_0=12\text{pc}$, $t=4\times {10}^6\text{yr}$, $t_0=4\times {10}^4\text{yr}$ and $z_{\text{OB}}=0$. The observational percentage reliability is ${\epsilon }_{\text{obs}}=\mathrm{93\; <mi>\; \%\; </mi>}$.

Figure 11. Geometrical section of the Fermi bubbles in the x-z plane with a Gaussian profile (green points) and observed profile (red stars). The parameters are $v_0=1000\text{km}\cdot {\text{s}}^{-1}$, $r_0=250\text{pc}$, $z_0=200\text{pc}$, $t=7.1\times {10}^5\text{yr}$, $t_0=7.1\times {10}^3\text{yr}$ and $z_{\text{OB}}=0$. The observational percentage reliability is ${\epsilon }_{\text{obs}}=\mathrm{92.12\; <mi>\; \%\; </mi>}$.

a Gaussian density profile as given by Equation (18).

The 3D advancing surface of the local bubble for the Gaussian case is shown in Figure 12.

3.3. The W4 Super-Bubble

The W4 super-bubble has been analysed from the point of view of the astronomical observations [31] [32] [33], in connection with the evolution of the magnetic field [34] and from a theoretical point of view [9] [35]. The upper part of Figure 3 in [36], which combines [SII], $H\alpha$ and [OIII] images has been digitized and will be the section of reference for W4, see Figure 13.

We now simulate the egg-shape of W4 when $z_{\text{OB}}\ne 0$. The numerical solution,

Figure 12. 3D surface of the Fermi bubbles with parameters as in Figure 11, Gaussian profile. The three Euler angles are $\Theta ={90}^{\circ }$, $\Phi =0^{\circ }$ and $\Psi ={90}^{\circ }$.

Figure 13. Section of the W4 + IC 1805 complex.

which is evaluated with the Euler method, is shown as a cut in the x-z plane: see Figure 14 for a density profile in the presence of self-gravity and Figure 15 for a Gaussian profile. The two adopted profiles in density are symmetric with respect

Figure 14. Geometrical section of the W4 super-bubble in the X-Z plane with a profile in the presence of self-gravity as given by Equation (9) (green points) and observed profile (red stars). The parameters are $v_0=20000\text{km}\cdot {\text{s}}^{-1}$, $r_0=1\text{pc}$, $h=17\text{pc}$, $t=1.3\times {10}^5\text{yr}$, $t_0=10\text{yr}$ and $z_{\text{OB}}=100$.

Figure 15. Geometrical section of the W4 super-bubble in the X-Z plane with a Gaussian profile (green points) and observed profile (red stars). The parameters are $v_0=4700\text{km}\cdot {\text{s}}^{-1}$, $r_0=1\text{pc}$, $z_0=55\text{pc}$, $t=2.25\times {10}^5\text{yr}$, $t_0=10\text{yr}$ and $z_{\text{OB}}=100$.

to the galactic plane, $Z=0$, but the simulated theoretical sections do not have an up-down symmetry, due to the fact that the expansion starts at $z=z_0$. Nevertheless, we still have a right-left symmetry.

The egg shape of the W4 super-bubble is shown in Figure 16.

The curious bump visible in the upper left part of Figure 13 could be an astronomical superposition of the image of IC 1805 on W4 or an intrinsic feature in the expansion of W4. In order to reproduce this feature, we assume that the scaling factor $z_{\mathrm{0,}\theta }$ in the interval ${\theta }_{inf}<\theta <{\theta }_{sup}$ varies with the following empirical law

$z_{\mathrm{0,}\theta }=z_0+z_{0}0.0006N\left(\theta \mathrm{<mo>\; ;\; </mo>}\bar{\theta }\mathrm{,}\sigma \right)$ (28)

where

$N\left(\theta \mathrm{<mo>\; ;\; </mo>}\sigma \mathrm{,}\mu \right)=\frac{1}{\sigma {\left(2\pi \right)}^{1/2}}\exp \left(-\frac{{\left(x-\bar{\theta }\right)}^2}{2{\sigma }^2}\right)\mathrm{,}$ (29)

is the Gaussian distribution, and $\bar{\theta }=\frac{{\theta }_{inf}+{\theta }_{sup}}{2}$ and $\sigma =\bar{\theta }/9$.

Figure 17 shows an “ad hoc” simulation of the bump of W4.

4. The Theory of the Image

In the framework of an optically thin medium, we outline a new analytical model which reproduces a theoretical vertical cut in the intensity of radiation and an old numerical model which simulates the intensity of radiation as a function of the point of view of the observer.

4.1. The Piriform Model

The piriform curve, or pear-shaped quartic, in 3D Cartesian coordinates $\left(x\mathrm{,}y\mathrm{,}z\right)$ has the equation

$a^4\left(x^2+y^2\right)-b^2{z}^3\left(2a-z\right)=\mathrm{0,}$ (30)

where a and b are both positive [37], see Figure 18 where the parameters a and b

Figure 16. 3D surface of the W4 super-bubble with parameters as in Figure 15, Gaussian profile. The three Euler angles are $\Theta ={90}^{\circ }$, $\Phi =0^{\circ }$ and $\Psi ={90}^{\circ }$.

Figure 17. Geometrical section of the W4 super-bubble in the X-Z plane with a Gaussian profile (green points) and observed profile (red stars). The parameters are the same as for Figure 15.

Figure 18. 3D display of the piriform curve when $a=3\text{kpc}$ and $b=3/2\text{kpc}$.

match the Fermi bubbles.

We are interested in a section of the above curve in the x-z plane which is obtained by inserting $y=0$

$a^4{x}^2=b^2{z}^3\left(2a-z\right).$ (31)

The parametric form of the piriform curve is

$x\left(\theta \mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\right)=a\left(1+\sin \left(\theta \right)\right)$ (32a)

$z\left(\theta \mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\right)=b\cos \left(\theta \right)\left(1+\sin \left(\theta \right)\right)\mathrm{,}$ (32b)

where $-\frac{\pi }{2}\le \theta \le \frac{3\pi }{2}$ and the maximum value reached along the z axis is

$z_{max}=2a\mathrm{.}$ (33)

We assume that the emission takes place in a thin layer comprised between an internal piriform which in polar coordinates has radius

$r_{int}=\frac{z}{a^2}\sqrt{a^4+2a{b}^{2}z-b^2{z}^2}\mathrm{,}$ (34)

and an external piriform which has radius

$r_{ext}=r_{int}+c\mathrm{,}$ (35)

where c is a positive parameter, see Figure 19. We therefore assume that the number density $C_m$ is constant between the two piriforms; as an example, along the z axis the number density increases from 0 at $\left(\mathrm{0,}{z}_{max}\right)$ to a maximum value $C_m$, remains constant up to $\left(\mathrm{0,}{z}_{max}+c\right)$, and then falls again to 0. The length of sight which produces the image in the first quadrant, when the observer is situated at the infinity of the x-axis, is the locus parallel to the x-axis which crosses the position z in the Cartesian x-z plane and terminates at the external piriform. In the case of an optically thin medium, the line of sight is split into two cases

$l_I\left(z\mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\mathrm{,}c\mathrm{,}{C}_m\right)=\frac{1}{a^2}\left(C_m\sqrt{a^4{c}^2+2b^2{z}^{3}a-b^2{z}^4+2\sqrt{\left(2az-z^2\right)b^2+a^4}{a}^{2}cz}\right)$ (36)

Figure 19. Internal and external piriforms in the first quadrant when $a=3\text{kpc}$, $b=3/2\text{kpc}$ and $c=3/20$. The dotted lines represent two different lengths of sight.

when $z_{max}\le z<z_{max}+c$

$\begin{array}{l}{l}_{II}\left(z\mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\mathrm{,}c\mathrm{,}{C}_m\right)\\ =\frac{1}{a^2}\left(-C_m\left(\sqrt{z\left(2a-z\right)}zb-\sqrt{a^4{c}^2+2b^2{z}^{3}a-b^2{z}^4+2\sqrt{\left(2az-z^2\right)b^2+a^4}{a}^{2}cz}\right)\right)\end{array}$ (37)

when $0\le z<z_{max}$.

A comparison between observed and theoretical intensity can be made by replacing in the above result $C_m$ with $I_m$ and doubling the length of sight due to the contribution of the second quadrant

$\begin{array}{l}{I}_I\left(z\mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\mathrm{,}c\mathrm{,}{I}_m\right)\\ =2\times \frac{1}{a^2}\left(I_m\sqrt{a^4{c}^2+2b^2{z}^{3}a-b^2{z}^4+2\sqrt{\left(2az-z^2\right)b^2+a^4}{a}^{2}cz}\right)\end{array}$ (38)

when $z_{max}\le z<z_{max}+c$

$\begin{array}{l}{I}_{II}\left(z\mathrm{<mo>\; ;\; </mo>}a\mathrm{,}b\mathrm{,}c\mathrm{,}{I}_m\right)\\ =2\times \frac{1}{a^2}\left(-I_m\left(\sqrt{z\left(2a-z\right)}zb-\sqrt{a^4{c}^2+2b^2{z}^{3}a-b^2{z}^4+2\sqrt{\left(2az-z^2\right)b^2+a^4}{a}^{2}cz}\right)\right)\end{array}$ (39)

when $0\le z<z_{max}$.

The resulting intensity is $I_{m}2c$ at $z=0$ and increases to $I_{m}2\sqrt{c}\sqrt{4a+c}$ at $z=z_{max}$, see Figure 20 for a typical profile in intensity along the z-axis.

4.2. The Numerical Model

The source of the luminosity is assumed here to be the flux of kinetic energy, $L_m$. The observed luminosity along a given direction can be expressed as

Figure 20. The intensity profile for the piriform model along the z-axis when when $a=3\text{kpc}$, $b=3/2\text{kpc}$, $c=3/20$ and $C_m=1$.

$L=\epsilon L_m\mathrm{,}$ (40)

where $\epsilon$ is a constant of conversion from the mechanical luminosity to the observed luminosity, for more details see [30]. The image of the Fermi bubbles is shown in Figure 21 and Figure 22 shows a cut of the intensity along the z-axis.

Figure 21. Map of the theoretical intensity of the Fermi bubbles for the model in the presence of self-gravity as given by Equation (9) with parameters as in Figure 10. The three Euler angles characterizing the orientation are $\Phi =0˚$, $\Theta =90˚$ and $\Psi =90˚$.

Figure 22. The intensity profile for the Fermi bubbles evolving in a model in the presence of self-gravity as given by Equation (9) along the z-axis, parameters as in Figure 21 (green empty stars) and the intensity profile for the piriform model along the z-axis when when $a=3.3\text{kpc}$, $b=1.65\text{kpc}$, $c=0.825$ and $C_m=1$ (red full line).

Figure 22 also shows the cut of the piriform model in order to evaluate the goodness of the analytical model for complex sections.

5. Conclusions

Equations of motion We derived two equations of motion coupling the thin layer approximation with the conservation of energy. The first model implements a profile in the presence of self-gravity of density and the second a Gaussian profile of density. In the absence of analytical results for the trajectory, with the exception of a Taylor expansion, we provided a numerical solution.

Comparison with other approaches

As an example, Figure 3 in [13] models the Eridanus-Orion structure with an ellipsoid, here we introduce the mushroom shape, see Figure 10 relative to the Fermi bubble and the egg shape, see Figure 16 relative to W4. We also suggested a first model for shapes apparently impossible to be simulated, see Figure 17 for the bump of W4.

Theory of the image The introduction of the piriform curve as a model for the section of the super-bubble confirms the existence of a characteristic “U” shape which has a maximum in the internal piriform at $z=2a$ and a minimum at the centre, $z=0$, see Equation (20). The superposition of a numerical cut with the piriform’s cut, (see Figure 22), shows us that the use of the piriform curve as a model is acceptable.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Heiles, C. (1979) H I Shells and Supershells. Astrophysical Journal, 229, 533-537. https://doi.org/10.1086/156986
[2]	Cash, W., Charles, P., Bowyer, S., Walter, F., Garmire, G. and Riegler, G. (1980) The X-Ray Superbubble in Cygnus. The Astrophysical Journal, 238, L71-L76. https://doi.org/10.1086/183261
[3]	Heiles, C. (1984) HI Shells, Supershells, Shell-Like Objects, and ‘‘Worms’’. The Astrophysical Journal Supplement Series, 55, 585-595. https://doi.org/10.1086/190970
[4]	Tenorio-Tagle, G. and Bodenheimer, P. (1988) Large-Scale Expanding Superstructures in Galaxies. Annual Review of Astronomy and Astrophysics, 26, 145-197. https://doi.org/10.1146/annurev.aa.26.090188.001045
[5]	McCray, R.A. (1987) Coronal Interstellar Gas and Supernova Remnants. In: Dalgarno A. and Layzer, D., Eds., Spectroscopy of Astrophysical Plasmas, Cambridge University Press, Cambridge, 255-278. https://doi.org/10.1017/CBO9780511564659.011
[6]	McCray, R. and Kafatos, M. (1987) Supershells and Propagating Star Formation. The Astrophysical Journal, 317, 190. https://doi.org/10.1086/165267
[7]	MacLow, M.M. and McCray, R. (1988) Superbubbles in Disk Galaxies. The Astrophysical Journal, 324, 776. https://doi.org/10.1086/165936
[8]	Igumenshchev, I.V., Shustov, B.M. and Tutukov, A.V. (1990) Dynamics of Supershells—Blow-Out. Astronomy and Astrophysics, 234, 396-402.
[9]	Basu, S., Johnstone, D. and Martin, P.G. (1999) Dynamical Evolution and Ionization Structure of an Expanding Superbubble: Application to W4. The Astrophysical Journal, 516, 843. https://doi.org/10.1086/307125
[10]	Kompaneyets, A.S. (1960) A Point Explosion in an Inhomogeneous Atmosphere. Doklady Physics, 5, 46.
[11]	Olano, C.A. (2009) The Propagation of the Shock Wave from a Strong Explosion in a Plane-Parallel Stratified Medium: The Kompaneets Approximation. Astronomy and Astrophysics, 506, 1215-1228. https://doi.org/10.1051/0004-6361/200912602
[12]	Pon, A., Johnstone, D., Bally, J. and Heiles, C. (2014) Kompaneets Model Fitting of the Orion-Eridanus Superbubble. Monthly Notices of the Royal Astronomical Society, 444, 3657-3669. https://doi.org/10.1093/mnras/stu1704
[13]	Pon, A., Ochsendorf, B.B., Alves, J., Bally, J., Basu, S. and Tielens, A.G.G.M. (2016) Kompaneets Model Fitting of the Orion-Eridanus Superbubble. II. Thinking Outside of Barnard’s Loop. The Astrophysical Journal, 827, 42. https://doi.org/10.3847/0004-637X/827/1/42
[14]	MacLow, M.M., McCray, R. and Norman, M.L. (1989) Superbubble Blowout Dynamics. The Astrophysical Journal, 337, 141. https://doi.org/10.1086/167094
[15]	Melioli, C., Brighenti, F., D’Ercole, A. and de Gouveia Dal Pino, E.M. (2009) Hydrodynamical Simulations of Galactic Fountains—II. Evolution of Multiple Fountains. Monthly Notices of the Royal Astronomical Society, 399, 1089-1105. https://doi.org/10.1111/j.1365-2966.2009.14725.x
[16]	Soler, J.D., Bracco, A. and Pon, A. (2018) The Magnetic Environment of the Orion-Eridanus Superbubble as Revealed by Planck. Astrophysics of Galaxies, 609, L3. https://doi.org/10.1051/0004-6361/201732203
[17]	Tomisaka, K. (1992) The Evolution of a Magnetized Superbubble. Publications of the Astronomical Society of Japan, 44, 177-191.
[18]	Rafikov, R.R. and Kulsrud, R.M. (2000) Magnetic Flux Expulsion in Powerful Superbubble Explosions and the α-Ω Dynamo. Monthly Notices of the Royal Astronomical Society, 314, 839-848. https://doi.org/10.1046/j.1365-8711.2000.03408.x
[19]	Bisnovatyi-Kogan, G.S. and Silich, S.A. (1995) Shock-Wave Propagation in the Nonuniform Interstellar Medium. Reviews of Modern Physics, 67, 661. https://doi.org/10.1103/RevModPhys.67.661
[20]	Dickey, J.M. and Lockman, F.J. (1990) H I in the Galaxy. Annual Review of Astronomy and Astrophysics, 28, 215-259. https://doi.org/10.1146/annurev.aa.28.090190.001243
[21]	Lockman, F.J. (1984) The H I Halo in the Inner Galaxy. Astrophysical Journal, 283, 90-97. https://doi.org/10.1086/162277
[22]	Tenenbaum, M. and Pollard, H. (1963) Ordinary Differential Equations: An Elementary Textbook for Students of Mathematics, Engineering, and the Sciences. Dover, New York.
[23]	Ince, E.L. (2012) Ordinary Differential Equations. Dover, New York.
[24]	Spitzer, Jr. L. (1942) The Dynamics of the Interstellar Medium. III. Galactic Distribution. Astrophysical Journal, 95, 329. https://doi.org/10.1086/144407
[25]	Rohlfs, K. (1977) Lectures on Density Wave Theory Vol. 69 of Lecture Notes in Physics. Springer-Verlag, Berlin.
[26]	Bertin, G. (2000) Dynamics of Galaxies. Cambridge University Press, Cambridge.
[27]	Padmanabhan, P. (2002) Theoretical Astrophysics. Vol. III: Galaxies and Cosmology. Cambridge University Press, Cambridge.
[28]	Olver, F.W.J., Lozier, D.W., Boisvert, R.F. and Clark, C.W. (2010) NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge.
[29]	Zaninetti, L. (2020) On the Shape of the Local Bubble. International Journal of Astronomy and Astrophysics, 10, 11-27. https://doi.org/10.4236/ijaa.2020.101002
[30]	Zaninetti, L. (2018) The Fermi Bubbles as a Superbubble. International Journal of Astronomy and Astrophysics, 8, 200-217. https://doi.org/10.4236/ijaa.2018.82015
[31]	Normandeau, M. and Basu, S. (1999) Observations and Modeling of the Disk-Halo Interaction in Our Galaxy. In: Taylor, A.R., Landecker, T.L. and Joncas, G., Eds., New Perspectives on the Interstellar Medium Vol. 168 of Astronomical Society of the Pacific Conference Series, Astron. Soc. Pac., San Francisco, 287.
[32]	Normandeau, M. (2000) The W4 Chimney/Superbubble. In: Alloin, D., Olsen, K. and Galaz, G., Eds., Stars, Gas and Dust in Galaxies: Exploring the Links Vol. 221 of Astronomical Society of the Pacific Conference Series, Astron. Soc. Pac., San Francisco, 41.
[33]	West, J.L., English, J., Normandeau, M. and Landecker, T.L. (2007) The Fragmenting Superbubble Associated with the H II Region W4. The Astrophysical Journal, 656, 914. https://doi.org/10.1086/510609
[34]	Gao, X.Y., Reich, W., Reich, P., Han, J.L. and Kothes, R. (2015) Magnetic Fields of the W4 Superbubble. Astronomy & Astrophysics, 578, Article ID: A24.
[35]	Baumgartner, V. and Breitschwerdt, D. (2013) Superbubble Evolution in Disk Galaxies. I. Study of Blow-Out by Analytical Models Astronomy & Astrophysics, 557, Article ID: A14.
[36]	Megeath, S.T., Townsley, L.K., Oey, M.S. and Tieftrunk, A.R. (2008) Low and High Mass Star Formation in the W3, W4, and W5 Regions. In: Reipurth, B., Ed., Handbook of Star Forming Regions, Chapter 9, Astronomical Society of the Pacific, 264-295.
[37]	Lawrence, J.D. (2013) A Catalog of Special Plane Curves. Dove, New York.