Energy conservation in the thin layer approximation: VI. Bubbles and super-bubbles

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.


Introduction
We now summarize the first uses of some words: 'super-shell' can be found in [1], where eleven H I objects were 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 was observed and 'worms', meaning gas filaments crawling away from the galactic plane in the inner Galaxy, can be found in [3]. Super-bubbles or supershells 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; see, 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.

The equations of motion
We start with the conservation of kinetic energy in spherical coordinates in the framework of the thin layer approximation where M 0 (r 0 ) and M (r) are the swept masses at r 0 and r, while v 0 and v are the velocities of the thin layer at r 0 and r. The above equation holds for the solid angle ∆Ω, 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: The asymmetry is due to a gradient of the number of particles with the distance or galactic height, z, which is parametrized as where n 1 =0.395 particles cm −3 , H 1 =127 pc, n 2 =0.107 particles cm −3 , H 2 =318 pc, n 3 =0.064 particles cm −3 , and H 3 =403 pc [19,20,21]. In the framework of Cartesian coordinates, (x, y, z), when the explosion starts at (0, 0, 0) we have an up-down symmetry, r(x, y, −z) = r(x, y, z) and a right-left symmetry r(x, −y, z) = r(x, y, z). Conversely, when the explosion starts at (0, 0, z OB ), 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(x, −y, z) = r(x, y, z).

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 3. The Euler method, which solves the following recursive equations where r n , v n , and M n are the temporary radius, velocity, and total mass, respectively, ∆t is the time step, and n is the index. 4. A power series solution of the form 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 (x, y, z) with the centre of the explosion at (0,0,0). The explosion is better visualized in a 3D Cartesian coordinate system (X, Y, Z) in which the galactic plane is given by Z = 0. The following translation, T OB , relates the two Cartesian coordinate systems where z OB is the distance in pc of the OB associations from the galactic plane. In the case of z OB = 0, the two masses which appear in Eq. (5b) should be carefully evaluated.

Medium in the presence of self-gravity
We assume that the number density distribution scales as 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 where ρ c is the density at z = 0. In order to find an acceptable value of h, we make a comparison with Eq. (3), after which we choose h = 90 pc, see Figure 1. The mass M 0 swept in the interval [0,r 0 ] is The total mass M (r; r 0 , ρ c , h) swept in the interval [0,r] is where θ is the polar angle and the polylog operator is defined by where Li s (z) is the Dirichlet series. The positive solution of Eq. (1) gives the velocity as a function of the radius: where The differential equation which governs the motion for the medium in the presence of self-gravity is 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 Eq. (15) gives 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%.

Gaussian medium
We assume that the number density distribution scales as 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 where ρ c is the density at z = 0. A comparison with Eq. (3) gives z 0 = 200 pc, see Figure 4. The total mass M (r; where  and erf(x) [28] is the error function defined by The velocity as a function of the radius is where and The differential equation which governs the motion for the Gaussian medium is and Figure 6 gives 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 9%.

Astrophysical applications
In the following we will analyse the local bubble, the Fermi bubble and the super bubble W4. An observational percentage reliability, obs , is introduced over the whole range of the polar angle θ, where r num is the theoretical radius of the considered bubble, r 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.

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 Eq. (9) and Figure 8 for a Gaussian density profile as given by Eq. (18). The 3D advancing surface of the local bubble for the case of self-gravity is shown in Figure 9.

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 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 Eq. (9) and Figure 11 for a Gaussian density profile as given by Eq. (18).  The 3D advancing surface of the local bubble for the Gaussian case is shown in Figure 12.

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α 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 OB = 0. The numerical solution, 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 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 0,θ in the interval θ inf < θ < θ sup varies with the following empirical law z 0,θ = z 0 + z0 0.0006N (θ; θ, σ) (28) where N (θ; σ, µ) = 1 is the Gaussian distribution, and θ = θ inf +θsup 2 and σ = θ/9. Figure 17 shows an 'ad hoc' simulation of the bump of W4.

The piriform model
The piriform curve, or pear-shaped quartic, in 3D Cartesian coordinates (x, y, z) has the equation where a and b are both positive [37], see Figure 18 where the parameters a and b 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 The parametric form of the piriform curve is where − π 2 ≤ θ ≤ 3π 2 and the maximum value reached along the z axis is We assume that the emission takes place in a thin layer comprised between an internal piriform which in polar coordinates has radius and an external piriform which has radius 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 (0, z max ) to a maximum value C m , remains constant up to (0, z max + c), 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 when 0 ≤ 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 The resulting intensity is I m 2 c at z = 0 and increases to I m 2 √ c √ 4 a + c at z = z max , see Figure 20 for a typical profile in intensity along the z-axis.

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 where 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.

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 = 2 a and a minimum at the centre, z = 0, see Eq. (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.