The Fermi bubbles as a Superbubble

In order to model the Fermi bubbles we apply the theory of the superbubble (SB). A thermal model and a self-gravitating model are reviewed. We introduce a third model based on the momentum conservation of a thin layer which propagates in a medium with an inverse square dependence for the density. A comparison have been made between the sections of the three models and the section of an observed map of the Fermi bubbles. An analytical law for the SB expansion as function of the time and polar angle is deduced. We derive a new analytical result for the image formation of the Fermi bubbles in an elliptical framework.


Introduction
The term super-shell was observationally defined by [1] as holes in the H I-column density distribution of our Galaxy. The dimensions of these objects span from 100 pc to 1700 pc and present elliptical shapes. These structures are commonly explained through introducing theoretical objects named bubbles or Superbubbles (SB); these are created by mechanical energy input from stars (see for example [2]; [3]).
The name Fermi bubbles starts to appear in the literature with the observations of Fermi-LAT which revealed two large gamma-ray bubbles, extending above and below the Galactic center, see [4]. Detailed observations of the Fermi bubbles analyzed the all-sky radio region, see [5], the Suzaku X-ray region, see [6,7,8], the ultraviolet absorption-line spectra, see [9,10], and the very high-energy gamma-ray emission, see [11]. The existence of the Fermi bubbles suggests some theoretical processes on how they are formed. We now outline some of them: processes connected with the galactic super-massive black hole in Sagittarius A, see [12,13,14]. Other studies try to explain the non thermal radiation from Fermi bubbles in the framework of the following physical mechanisms: electron's acceleration inside the bubbles, see [15], hadronic models, see [16,17,18,19,11], and leptonic models, see [20,19]. The previous theoretical efforts allow to build a dynamic model for the Fermi bubbles for which the physics remains unknown. The layout of the paper is as follows. In Section 2 we analyze three profiles in vertical density for the Galaxy. In Section 3 we review two existing equations of motion for the Fermi Bubbles and we derive a new equation of motion for an inverse square law in density. In Section 4 we discuss the results for the three equations of motion here adopted in terms of reliability of the model. In Section 5 we derive two results for the Fermi bubbles: an analytical model for the cut in intensity in an elliptical framework and a numerical map for the intensity of radiation based on the numerical section.

The profiles in density
This section reviews the gas distribution in the galaxy. A new inverse square dependence for the gas in introduced. In the following we will use the spherical coordinates which are defined by the radial distance r, the polar angle θ, and the azimuthal angle ϕ.

Gas distribution in the galaxy
The vertical density distribution of galactic neutral atomic hydrogen (H I) is wellknown; specifically, it has the following three component behavior as a function of z, the distance from the galactic plane in pc: We took [21,22,23] 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. This distribution of galactic H I is valid in the range 0.4 ≤ r ≤ r 0 , where r 0 = 8.5 kpc and r is the distance from the center of the galaxy.
A recent evaluation for galactic H I quotes:

The thermal model
The starting equation for the evolution of the SB [30,31,32] is momentum conservation applied to a pyramidal section. The parameters of the thermal model are N * , the number of SN explosions in 5.0 · 10 7 yr, z OB , the distance of the OB associations from the galactic plane, E 51 , the energy in 10 51 erg, v 0 , the initial velocity which is fixed by the bursting phase, t 0 , the initial time in yr which is equal to the bursting time, and t the proper time of the SB. The SB evolves in a standard three component medium, see formula (1).

A recursive cold model
The 3D expansion that starts at the origin of the coordinates; velocity and radius are given by a recursive relationship, see [33]. The parameters are the same of the thermal model and the SB evolves in a self-gravitating medium as given by equation (3).

The inverse square model
In the case of an inverse square density profile for the interstellar medium ISM as given by equation (4), the differential equation which models momentum conservation is where the initial conditions are r = r 0 and v = v 0 when t = t 0 . We now briefly review that given a function f (r), the Padé approximant, after [34], is where the notation is the same of [35]. The coefficients a i and b i are found through Wynn's cross rule, see [36,37] and our choice is o = 2 and q = 1. The choice of o and q is a compromise between precision, high values for o and q, and simplicity of the expressions to manage, low values for o and q. The inverse of the velocity is where N N = (cos (θ)) 5 r 0 4 r + (cos (θ)) 4 r 0 4 z 0 + (cos (θ)) 4 r 0 3 rz 0 +3 (cos (θ)) 2 r 2 z 0 3 − 6 cos (θ) ln (r cos (θ) + z 0 ) r 0 z 0 4 Table 1. Numerical values of the parameters for the simulation in the case of the inverse square model. and The above result allows deducing a solution r 2,1 expressed through the Padè approximant with and A possible set of initial values is reported in Table 1 in which the initial value of radius and velocity are fixed by the bursting phase. The above parameters allows to obtain an approximate expansion law as function of time and polar angle with BN = 31944000 (cos (θ)) 2 + 18.8632 t cos (θ) + 5111040 cos (θ) and

Astrophysical Results
This section introduces a test for the reliability of the model, analyzes the observational details of the Fermi bubbles, reviews the results for the two models of reference and reports the results of the inverse square model.

The reliability of the model
An observational percentage reliability, ǫ obs , is introduced over the whole range of the polar angle θ, where r num is the theoretical radius, r obs is the observed radius, and the index j varies from 1 to the number of available observations.

The structure of the Fermi bubbles
The exact shape of the Fermi bubbles is a matter of research and as an example in [38] the bubbles are modeled with ellipsoids centered at 5 kpc up and below the Galactic plane with semi-major axes of 6 kpc and minor axes of 4 kpc. In order to test our models we selected the image of the Fermi bubbles available at https://www.nasa.gov/mission_pages/GLAST/news/new-structure.html which is reported in Figure 1. A digitalization of the above advancing surface is reported in Figure 2 as a 2D section. This allows to fix the observed radii to be inserted in equation (20). The actual shape of the bubbles in galactic coordinates is shown in Figure 3 and 15 of [4] and Figure 30 and Table 3 of by [39].

The two models of reference
The thermal model is outlined in Section 3.1 and Figure 3 reports the numerical solution as a cut in the x − z plane. The cold recursive model is outlined in Section 3.2 and Figure 4 reports the numerical solution as a cut in the x − z plane.

The inverse square model
The inverse square model is outlined in Section 3.3 and Figure 5 reports the numerical solution as a cut in the x − z plane.
A rotation around the z-axis of the above theoretical section allows building a 3D surface, see Figure 6. The temporal evolution of the advancing surface is reported in Figure 7 and a comparison should be done with Fig. 6 in [40].

Theory of the image
This section reviews the transfer equation and reports a new analytical result for the intensity of radiation in an elliptical framework in the non-thermal/thermal case. A numerical model for the image formation of the Fermi bubbles is reported. Section of the Fermi bubbles in the x − z plane with a thermal model (green points) and observed profile (red stars). The bursting parameters N * = 113000, z OB =0 pc, E 51 =1, t 0 = 0.036 10 7 yr when t = 90 10 7 yr give ǫ obs = 73.34%.

The transfer equation
The transfer equation in the presence of emission only in the case of optically thin layer is where K is a constant, j ν is the emission coefficient, the index ν denotes the frequency of emission and C(s) is the number density of particles, see for example [41]. As an example the synchrotron emission, as described in sec. 4 of [42], is often used in order to model the radiation from a SNR, see for example [43,44,45]. According to the above equation the increase in intensity is proportional to the number density integrated along the line of sight, which for constant density, gives where K ′ is a constant and l is the length along the line of sight interested in the emission; in the case of synchrotron emission see formula (1.175) in [46].

Analytical non thermal model
A real ellipsoid, see [47], represents a first approximation of the Fermi bubbles, see [38], and has equation in which the polar axis of the Galaxy is the z-axis. Figure 8 reports the astrophysical application of the ellipsoid in which due to the symmetry about the azimuthal angle b = d. Section of the Fermi bubbles in the x − z plane with the inverse square model (green points) and observed profile (red stars). The parameters are reported in Table 1 and the reliability is ǫ obs = 90.71%.
We are interested in the section of the ellipsoid y = 0 which is defined by the following external ellipse We assume that the emission takes place in a thin layer comprised between the external ellipse and the internal ellipse defined by see Figure 9. We therefore assume that the number density C is constant and in particular rises from 0 at (0,a) to a maximum value C m , remains constant up to (0,ac) and then falls again to 0. The length of sight, when the observer is situated at the  Table 1, inverse square profile. The three Euler angles are Θ = 40, Φ = 60 and Ψ = 60.    infinity of the x-axis, is the locus parallel to the x-axis which crosses the position z in a Cartesian x − z plane and terminates at the external ellipse. The locus length is when 0 ≤ z < (a − c) .
In the case of optically thin medium, according to equation (22), the intensity is split in two cases when (a − c) ≤ z < a where I m is a constant which allows to compare the theoretical intensity with the observed one. A typical profile in intensity along the z-axis is reported in Figure 10. The ratio, r, between the theoretical intensity at the maximum, (z = a − c), and at the minimum, (z = 0), is given by As an example the values a = 6 kpc,b = 4 kpc, c = a 12 kpc gives r = 3.19. The knowledge of the above ratio from the observations allows to deduce c once a and b are given by the observed morphology As an example in the inner regions of the northeast Fermi bubble we have r = 2, see [6], which coupled with a = 6 kpc and b = 4 kpc gives c = 1.2 kpc. The above Figure 11. The intensity profile along the z-axis for the thermal model when a = 6 kpc, b = 4 kpc and Im=1.
value is an important astrophysical result because we have found the dimension of the advancing thin layer.

Analytical thermal model
A thermal model for the image is characterized by a constant temperature in the internal region of the advancing section which is approximated by an ellipse, see equation (24). We therefore assume that the number density C is constant and in particular rises from 0 at (0,a) to a maximum value C m , remains constant up to (0,-a) and then falls again to 0. The length of sight, 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 a Cartesian x − z plane and terminates at the external ellipse in the point (0,a). The locus length is The number density C m is constant in the ellipse and therefore the intensity of radiation is A typical profile in intensity along the z-axis for the thermal model is reported in Figure 11.

Numerical model
The source of luminosity is assumed here to be the flux of kinetic energy, L m , where A is the considered area, V the velocity and ρ the density, see formula (A28) in [48]. In our case A = R 2 ∆Ω, where ∆Ω is the considered solid angle along the chosen direction. 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. A numerical algorithm which allows us to build a complex image is outlined in Section 4.1 of [49] and the orientation of the object is characterized by the Euler angles (Φ, Θ, Ψ) The threshold intensity can be parametrized to I max , the maximum value of intensity characterizing the map. The image of the Fermi bubbles is shown in Figure 12 and the introduction of a threshold intensity is visualized in Figure 13.

Conclusions
Law of motion We have compared two existing models for the temporal evolution of the Fermi bubbles, a thermal model, see Section 3.1, and an autogravitating model, see Section 3.2, with a new model which conserves the momentum in presence of an inverse square law for the density of the ISM. The best result is obtained by the inverse square model which produces a reliability of ǫ obs = 90.71% for the expanding radius in respect to a digitalized section of the Fermi bubbles. A semi-analytical law of motion as function of polar angle and time is derived for the inverse square model, see equation (17). Formation of the image An analytical cut for the intensity of radiation along the z-axis is derived in the framework of advancing surface characterized by an internal and an external ellipses. The analytical cut in theoretical intensity presents a characteristic "U" shape which has a maximum in the external ring and a minimum at the center, see equation (30). The presence of a hole in the intensity of radiation in the central region of the elliptical Fermi bubbles is also confirmed by a numerical algorithm for the image formation, see Figure 13. The theoretical prediction of a hole in the intensity map explains the decrease in intensity for the 0.3 kev plasma by = 50% toward the central region of the northeast Fermi bubble, see [6]. The intensity of radiation for the thermal model conversely presents a maximum of the intensity at the center of the elliptical Fermi bubble, see equation (33) and this theoretical prediction does not agree with the above observations.