On the Shape of the Local Bubble

Lorenzo Zaninetti

Physics Department, Turin, Italy.

DOI: 10.4236/ijaa.2020.101002                

Abstract

The shape of the local bubble is modeled in the framework of the thin layer approximation. The asymmetric shape of the local bubble is simulated by introducing axial profiles for the density of the interstellar medium, such as exponential, Gaussian, inverse square dependence and Navarro-Frenk-White. The availability of some observed asymmetric profiles for the local bubble allows us to match theory and observations via the observational percentage of reliability. The model is compatible with the presence of radioisotopes on Earth.

Keywords

ISM: Bubbles, ISM: Clouds, Galaxy: Disk

1. Introduction

The local bubble (LB) is a low-density region that surrounds the Sun. Because it is emitting in the X-rays, it is also called Local Hot Bubble (LHB), see [1] [2]. In the framework of thermal ionization equilibrium, the temperature is $kT=\left(0.097\pm 0.013\right)\text{keV}$ or $T=\left(1.1252\pm 0.15\right)\times {10}^6\text{K}$, see [3]. Recently, the following features of the LB have been discussed: the variations of the polarization degree P, see [4]; and the polarization from the interstellar medium, due to irregular dust grains aligned with the magnetic field, see [5]. The presence of ⁶⁰Fe in deep-sea measurements on Earth has triggered the study of the LB-sun interaction, see [6]. We now select some theoretical efforts that model the LB, as follows: the one-dimensional hydrocode with non-equilibrium ion evolution and dust, see [7] [8]; different tests to explain the FUSE data, see [9] [10]; the use of the parallel adaptive mesh refinement code EAF-PAMR, see [11]; hydrodynamical simulations of the LB, see [12]; and the study of the 3D structure of the magnetic field, see [13].

These models leave some questions unanswered, or only partially answered, as follows:

- Can we model the LB in the framework of the thin layer expansion of a shell in an interstellar medium (ISM) with symmetry in respect to the equatorial plane of the explosion?

- Can we compare the data of the theoretical expansion, which is a function of the latitude, with the observed profiles of expansion of the LB?

- What is the range of reliability of the Taylor expansion and Padé approximation of the theoretical expansion in the framework of the thin layer?

- Can we model the LB-Sun interaction?

To answer these questions: Section 2 analyzes four profiles of density for the interstellar medium (ISM); Section 3 derives four equations of motion for the LB; and Section 4 discusses the results for the four equations of motion in terms of reliability of the model, it also introduces the interaction of many bubbles, discusses the ⁶⁰Fe-signature and explores the interaction of many bubbles.

2. The Density Profiles

A point in Cartesian coordinates is characterized by x, y and z, and the position of the origin is the center of the LB. The same point in spherical coordinates is characterized by the radial distance $r\in \left[\mathrm{0,}\infty \right]$, the polar angle $\theta \in \left[\mathrm{0,}\pi \right]$, and the azimuthal angle $\phi \in \left[\mathrm{0,2}\pi \right]$.

The following profiles are considered: exponential, Gaussian, inverse square dependence and Navarro-Frenk-White.

2.1. An Exponential Profile

The density is assumed to have the following exponential dependence on z in Cartesian coordinates:

$\rho \left(z\mathrm{<mo>\; ;\; </mo>}b\mathrm{,}{\rho }_0\right)={\rho }_0\exp \left(-z/b\right)\mathrm{,}$ (1)

where b represents the scale. In spherical coordinates, the density has the following piecewise dependence

$\rho \left(r;r_0,b,{\rho }_0\right)=(\begin{array}{ll}{\rho }_0\hfill & \text{if}r\le r_0\hfill \\ {\rho }_0\exp \left(-\frac{r\cos \left(\theta \right)}{b}\right)\hfill & \text{if}r>r_0\hfill \end{array}$ (2)

which has a jump discontinuity at $r=r_0$ when $\theta >0$. Given a solid angle $\Delta \Omega$, the total mass swept, $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}\theta \mathrm{,}{\rho }_0\mathrm{,}\Delta \Omega \right)$, in the interval $\left[\mathrm{0,}r\right]$ is

$\begin{array}{l}M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}\theta \mathrm{,}{\rho }_0\mathrm{,}\Delta \Omega \right)\\ =(\frac{1}{3}{\rho }_0{r}_0^3-\frac{b\left(r^2{\left(\cos \left(\theta \right)\right)}^2+2rb\cos \left(\theta \right)+2b^2\right){\rho }_0}{{\left(\cos \left(\theta \right)\right)}^3}{\text{e}}^{-\frac{r\cos \left(\theta \right)}{b}}\\ +\frac{b\left(r_0^2{\left(\cos \left(\theta \right)\right)}^2+2r_{0}b\cos \left(\theta \right)+2b^2\right){\rho }_0}{{\left(\cos \left(\theta \right)\right)}^3}{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{b}})\Delta \Omega \end{array}$ (3)

2.2. A Gaussian Profile

The density is assumed to have the following Gaussian dependence on z in Cartesian coordinates:

$\rho \left(z;b,{\rho }_0\right)={\rho }_0{\text{e}}^{-\frac{1}{2}\frac{z^2}{b^2}},$ (4)

where b represents a parameter. In spherical coordinates, the density is

$\rho \left(r;r_0,b,{\rho }_0\right)=(\begin{array}{ll}{\rho }_0\hfill & \text{if}r\le r_0\hfill \\ {\rho }_0{\text{e}}^{-\frac{1}{2}\frac{z^2}{b^2}}\hfill & \text{if}r>r_0\hfill \end{array}$ (5)

and presents a jump discontinuity at $r=r_0$ when $\theta >0$. The total mass swept, $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}\theta \mathrm{,}{\rho }_0\right)$, in the interval $\left[\mathrm{0,}r\right]$ is

$\begin{array}{l}M\left(r;r_0,b,\theta ,{\rho }_0\right)\\ =(\frac{1}{3}{\rho }_0{r}_0^3+{\rho }_0\left(-\frac{r{b}^2}{{\left(\cos \left(\theta \right)\right)}^2}{\text{e}}^{-\frac{1}{2}\frac{r^2{\left(\cos \left(\theta \right)\right)}^2}{b^2}}+\frac{1}{2}\frac{b^3\sqrt{\pi }\sqrt{2}}{{\left(\cos \left(\theta \right)\right)}^3}\text{erf}\left(\frac{1}{2}\frac{\sqrt{2}\cos \left(\theta \right)r}{b}\right)\right)\\ -{\rho }_0\left(-\frac{r_0{b}^2}{{\left(\cos \left(\theta \right)\right)}^2}{\text{e}}^{-\frac{1}{2}\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{b^2}}+\frac{1}{2}\frac{b^3\sqrt{\pi }\sqrt{2}}{{\left(\cos \left(\theta \right)\right)}^3}\text{erf}\left(\frac{1}{2}\frac{\sqrt{2}\cos \left(\theta \right)r_0}{b}\right)\right))\Delta \Omega \end{array}$ (6)

where $\text{erf}\left(x\right)$ is the error function, defined by

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

2.3. The Inverse Square Dependence

The density is assumed to have the following dependence on z in Cartesian coordinates,

$\rho \left(z\mathrm{<mo>\; ;\; </mo>}{z}_0\mathrm{,}{\rho }_0\right)={\rho }_0{\left(1+\frac{z}{z_0}\right)}^{-2}\mathrm{.}$ (8)

In this paper, we will adopt the following density profile in spherical coordinates

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

where the parameter $z_0$ fixes the scale and ${\rho }_0$ is the density at $z=z_0$. The above density presents a jump discontinuity at $r=r_0$ when $\theta >0$. The mass $M_0$ swept in the interval $\left[\mathrm{0,}{r}_0\right]$ is

$M_0=\frac{1}{3}{\rho }_0{r}_0^3\Delta \Omega$ (10)

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

$M\left(r;r_0,z_0,\theta ,{\rho }_0,\Delta \Omega \right)$

$\begin{array}{l}=(\frac{1}{3}{\rho }_0{r}_0^3+\frac{{\rho }_0{b}^{2}r}{{\left(\cos \left(\theta \right)\right)}^2}-2\frac{{\rho }_0{b}^3\ln \left(r\cos \left(\theta \right)+b\right)}{{\left(\cos \left(\theta \right)\right)}^3}-\frac{{\rho }_0{b}^4}{{\left(\cos \left(\theta \right)\right)}^3\left(r\cos \left(\theta \right)+b\right)}\\ -\frac{{\rho }_0{b}^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^2}+2\frac{{\rho }_0{b}^3\ln \left(r_0\cos \left(\theta \right)+b\right)}{{\left(\cos \left(\theta \right)\right)}^3}+\frac{{\rho }_0{b}^4}{{\left(\cos \left(\theta \right)\right)}^3\left(r_0\cos \left(\theta \right)+b\right)})\Delta \Omega \end{array}$ (11)

2.4. Navarro-Frenk-White Profile

The usual Navarro-Frenk-White (NFW) distribution has a dependence on r in spherical coordinates of the type

$\rho \left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}{\rho }_0\right)=\frac{{\rho }_0{r}_0{\left(b+r_0\right)}^2}{r{\left(b+r\right)}^2}\mathrm{,}$ (12)

where b represents the scale, see [14] for more details. The NFW profile along the axis z can be obtained by substituting r with $r\cos \left(\theta \right)=z$

$\rho \left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}{\rho }_0\mathrm{,}\theta \right)=\frac{{\rho }_0{r}_0{\left(b+r_0\right)}^2}{r\cos \left(\theta \right){\left(b+r\cos \left(\theta \right)\right)}^2}\mathrm{,}$ (13)

The piece-wise density is

$\rho \left(r;r_0,b,{\rho }_0\theta \right)=\{\begin{array}{ll}{\rho }_0\hfill & \text{if}r\le r_0\hfill \\ \frac{{\rho }_0{r}_0{\left(b+r_0\right)}^2}{r\cos \left(\theta \right){\left(b+r\cos \left(\theta \right)\right)}^2}\hfill & \text{if}r>r_0\hfill \end{array}$ (14)

and has a jump discontinuity at $r=r_0$ when $\theta >0$. The total mass swept, $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}b\mathrm{,}{\rho }_0\theta \right)$, in the interval $\left[\mathrm{0,}r\right]$ is

$\begin{array}{l}M\left(r;r_0,b,\theta ,{\rho }_0,\Delta \Omega \right)\\ =(\frac{1}{3}{\rho }_0{r}_0^3+\frac{{\rho }_0\left(\left(b+r\cos \left(\theta \right)\right)\ln \left(b+r\cos \left(\theta \right)\right)+b\right){\left(b+r_0\right)}^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^3\left(b+r\cos \left(\theta \right)\right)}\text{=}\\ -\frac{{\rho }_0\left(\left(b+r_0\cos \left(\theta \right)\right)\ln \left(b+r_0\cos \left(\theta \right)\right)+b\right){\left(b+r_0\right)}^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^3\left(b+r_0\cos \left(\theta \right)\right)})\Delta \Omega \end{array}$ (15)

3. The Thin Layer Approximation

The conservation of the momentum in spherical coordinates along the solid angle $\Delta \Omega$ in the framework of the thin layer approximation states that

$M_0\left(r_0\right)v_0=M\left(r\right)v\mathrm{,}$ (16)

where $M_0\left(r_0\right)$ and $M\left(r\right)$ are the swept masses at $r_0$ and r, and $v_0$ and v are the velocities of the thin layer at $r_0$ and r. This conservation law can be expressed as a differential equation of the first order by inserting $v=\frac{\text{d}r}{\text{d}t}$ :

$M\left(r\right)\frac{\text{d}r}{\text{d}t}=M_0\left(r_0\right)v_0\mathrm{.}$ (17)

In the first phase from $r=0$ to $r=r_0$ the density is constant and the explosion is symmetrical. In the second phase the density is function of the polar angle $\theta$ and therefore the shape of the advancing expansion is asymmetrical. The equation of motion for the four profiles is now derived.

3.1. Motion with a Constant Density

In the case of constant density of the ISM, the analytical solution for the trajectory is

$r\left(t\mathrm{<mo>\; ;\; </mo>}{t}_0\mathrm{,}{r}_0\mathrm{,}{v}_0\right)=\sqrt[4]{4r_0^3{v}_0\left(t-t_0\right)+r_0^4}\mathrm{,}$ (18)

and the velocity is

$v\left(t\mathrm{<mo>\; ;\; </mo>}{t}_0\mathrm{,}{r}_0\mathrm{,}{v}_0\right)=\frac{r_0^3{v}_0}{{\left(4r_0^3{v}_0\left(t-t_0\right)+r_0^4\right)}^{3/4}}\mathrm{,}$ (19)

where $r_0$ and $v_0$ are the position and the velocity when $t=t_0$, see [15] [16].

3.2. Motion with an Exponential Profile

In the case of an exponential density profile for the ISM, as given by Equation (2), the differential equation that models momentum conservation is

$\begin{array}{l}(\frac{1}{3}{r}_0^3-\frac{b\left({\left(r\left(t\right)\right)}^2{\left(\cos \left(\theta \right)\right)}^2+2r\left(t\right)b\cos \left(\theta \right)+2b^2\right)}{{\left(\cos \left(\theta \right)\right)}^3}{\text{e}}^{-\frac{\cos \left(\theta \right)r\left(t\right)}{b}}\\ +\frac{b\left(r_0^2{\left(\cos \left(\theta \right)\right)}^2+2r_{0}b\cos \left(\theta \right)+2b^2\right)}{{\left(\cos \left(\theta \right)\right)}^3}{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{b}})\frac{\text{d}}{\text{d}t}r\left(t\right)=\frac{1}{3}{r}_0^3{v}_0.\end{array}$ (20)

There is no analytical solution.

3.3. Motion with a Gaussian Profile

In the case of a Gaussian density profile for the ISM, as given by Equation (5), the differential equation that models momentum conservation is

$\begin{array}{l}(\frac{1}{3}{\rho }_0{r}_0^3+{\rho }_0\left(-\frac{r\left(t\right)b^2}{{\left(\cos \left(\theta \right)\right)}^2}{\text{e}}^{-\frac{1}{2}\frac{{\left(r\left(t\right)\right)}^2{\left(\cos \left(\theta \right)\right)}^2}{b^2}}+\frac{1}{2}\frac{b^3\sqrt{\pi }\sqrt{2}}{{\left(\cos \left(\theta \right)\right)}^3}\text{erf}\left(\frac{1}{2}\frac{\sqrt{2}\cos \left(\theta \right)r\left(t\right)}{b}\right)\right)\\ -{\rho }_0\left(-\frac{r_0{b}^2}{{\left(\cos \left(\theta \right)\right)}^2}{\text{e}}^{-\frac{1}{2}\frac{r_0^2{\left(\cos \left(\theta \right)\right)}^2}{b^2}}+\frac{1}{2}\frac{b^3\sqrt{\pi }\sqrt{2}}{{\left(\cos \left(\theta \right)\right)}^3}\text{erf}\left(\frac{1}{2}\frac{\sqrt{2}\cos \left(\theta \right)r_0}{b}\right)\right))\frac{\text{d}}{\text{d}t}r\left(t\right)\\ =\frac{1}{3}{\rho }_0{r}_0^3{v}_0.\end{array}$ (21)

3.4. Motion with an Inverse Square Dependence

In the case of an inverse square density profile for the ISM, as given by Equation (9), the differential equation that models the momentum conservation is

$\begin{array}{l}(\frac{1}{3}{\rho }_0{r}_0^3+\frac{{\rho }_0{z}_0^{2}r\left(t\right)}{{\left(\cos \left(\theta \right)\right)}^2}-2\frac{{\rho }_0{z}_0^3\ln \left(r\left(t\right)\cos \left(\theta \right)+z_0\right)}{{\left(\cos \left(\theta \right)\right)}^3}\\ -\frac{{\rho }_0{z}_0^4}{{\left(\cos \left(\theta \right)\right)}^3\left(r\left(t\right)\cos \left(\theta \right)+z_0\right)}-\frac{{\rho }_0{z}_0^2{r}_0}{{\left(\cos \left(\theta \right)\right)}^2}+2\frac{{\rho }_0{z}_0^3\ln \left(r_0\cos \left(\theta \right)+z_0\right)}{{\left(\cos \left(\theta \right)\right)}^3}\\ +\frac{{\rho }_0{z}_0^4}{{\left(\cos \left(\theta \right)\right)}^3\left(r_0\cos \left(\theta \right)+z_0\right)})\frac{\text{d}}{\text{d}t}r\left(t\right)-\frac{1}{3}{\rho }_0{r}_0^3{v}_0=0.\end{array}$ (22)

There is not an analytical solution for this differential equation.

3.5. Motion with a Navarro-Frenk-White Profile

In the case of a Navarro-Frenk-White density profile for the ISM, as given by Equation (13), the differential equation that models momentum conservation is

$\begin{array}{l}(\frac{1}{3}{\rho }_0{r}_0^3+\frac{r_0{\rho }_0\left(\left(b+r\left(t\right)\cos \left(\theta \right)\right)\ln \left(b+r\left(t\right)\cos \left(\theta \right)\right)+b\right){\left(b+r_0\right)}^2}{{\left(\cos \left(\theta \right)\right)}^3\left(b+r\left(t\right)\cos \left(\theta \right)\right)}\\ -\frac{r_0{\rho }_0\left(\left(b+r_0\cos \left(\theta \right)\right)\ln \left(b+r_0\cos \left(\theta \right)\right)+b\right){\left(b+r_0\right)}^2}{{\left(\cos \left(\theta \right)\right)}^3\left(b+r_0\cos \left(\theta \right)\right)})\frac{\text{d}}{\text{d}t}r\left(t\right)=\frac{1}{3}{\rho }_0{r}_0^3{v}_0\mathrm{.}\end{array}$ (23)

A first approximated solution of this differential equation can be given as a series of order 4

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

Figure 1 reports a comparison between numerical and series solution.

To find a second approximate solution for this differential equation of the first order in r, we separate the variables and we integrate. The following non-linear equation is obtained

$\frac{N}{D}=t-t_0,$ (25)

where

$\begin{array}{c}N=-6{\left(b+r_0\right)}^2\left(b+r_0\cos \left(\theta \right)\right)\left(1/2r\cos \left(\theta \right)+b\right)\ln \left(b+r_0\cos \left(\theta \right)\right)\\ +6{\left(b+r_0\right)}^2\left(b+r_0\cos \left(\theta \right)\right)\left(1/2r\cos \left(\theta \right)+b\right)\ln \left(b+r\cos \left(\theta \right)\right)\\ -6\cos \left(\theta \right)\left(r-r_0\right)(-1/6r_0^3{\left(\cos \left(\theta \right)\right)}^4-1/6b{r}_0^2{\left(\cos \left(\theta \right)\right)}^3\\ +1/2r_0{\left(b+r_0\right)}^2\cos \left(\theta \right)+b{\left(b+r_0\right)}^2),\end{array}$ (26)

and

$D=v_0{r}_0^2{\left(\cos \left(\theta \right)\right)}^4\left(b+r_0\cos \left(\theta \right)\right)\mathrm{.}$ (27)

Figure 1. Numerical solution (full line) and series solution (dashed line). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $b=2.8\text{pc}$, $t=1360\text{yr}$ and $t_0=680\text{yr}$.

In this case, it is not possible to find an analytical solution for the radius, r, as a function of time. Therefore, we apply the Padé rational polynomial, see [17] [18] [19] [20]. We choose an approximation of degree 2 in the numerator and degree 1 in the denominator about the point $r=r_0$ to the left-hand side of Equation (25). The resulting equation of second degree is

$\frac{NN}{DD}=t-t_0,$ (28)

where

$\begin{array}{c}NN=-\left(r_0-r\right)(4{\left(\cos \left(\theta \right)\right)}^3{r}_0^3+2{\left(\cos \left(\theta \right)\right)}^3{r}_0^{2}r+12{\left(\cos \left(\theta \right)\right)}^2{r}_0^{2}b\\ +8r_0\cos \left(\theta \right)b^2-2\cos \left(\theta \right)b^{2}r-9r_0^3-18b{r}_0^2+9r_0^{2}r\\ {}_{{}_{}}{}^{}-9r_0{b}^2+18r_{0}br+9b^{2}r)\mathrm{,}\end{array}$ (29)

and

$DD=2\cos \left(\theta \right)v_0\left(b+r_0\cos \left(\theta \right)\right)\left(2\cos \left(\theta \right)r_0^2+r_{0}r\cos \left(\theta \right)+4b{r}_0-br\right).$ (30)

The resulting Padé approximant for the trajectory, the radius $r_{\mathrm{2,1}}$, is the second approximated solution

$r_{\mathrm{2,1}}\left(t\mathrm{<mo>\; ;\; </mo>}{t}_0\mathrm{,}{r}_0\mathrm{,}{v}_0\mathrm{,}b\right)=\frac{NNN}{DDD}\mathrm{.}$ (31)

where

$\begin{array}{c}NNN={\left(\cos \left(\theta \right)\right)}^{3}t{r}_0^2{v}_0-{\left(\cos \left(\theta \right)\right)}^3{r}_0^2{t}_0{v}_0-{\left(\cos \left(\theta \right)\right)}^3{r}_0^3\\ -6{\left(\cos \left(\theta \right)\right)}^2{r}_0^{2}b-\cos \left(\theta \right)b^{2}t{v}_0+\cos \left(\theta \right)b^2{t}_0{v}_0\\ -5r_0\cos \left(\theta \right)b^2+9r_0{b}^2+18r_0^{2}b+9r_0^3+\sqrt{A},\end{array}$ (32)

and

$DDD=2{\left(\cos \left(\theta \right)\right)}^3{r}_0^2-2b^2\cos \left(\theta \right)+9b^2+18b{r}_0+9r_0^2,$ (33)

and

$\begin{array}{c}A={\left(b+r_0\cos \left(\theta \right)\right)}^2\cos \left(\theta \right)({\left(3r_0+v_0\left(t-t_0\right)\right)}^2{r}_0^2{\left(\cos \left(\theta \right)\right)}^3\\ -2\left(3r_0+v_0\left(t-t_0\right)\right)\left(-3r_0+v_0\left(t-t_0\right)\right)r_{0}b{\left(\cos \left(\theta \right)\right)}^2\\ +{\left(-3r_0+v_0\left(t-t_0\right)\right)}^2{b}^2\cos \left(\theta \right)+54r_0{v}_0{\left(b+r_0\right)}^2\left(t-t_0\right)).\end{array}$ (34)

Figure 2 reports a comparison between the numerical and the series solution.

The two approximations that we have used here cover the range in time after which the percent error is ≈10%: 1360 yr for the Taylor series and 194,285 yr for the Padé approximant. We conclude that in our case the Padé approximant has a wider radius of convergence in respect to the Taylor series.

4. Astrophysical Results

The adopted astrophysical units are pc for length and yr for time; while the initial velocity $v_0$ is expressed in pc∙yr⁻¹. The astronomical velocities are evaluated in km∙s⁻¹ and therefore $v_0=1.02\times {10}^{-6}{v}_1$ where $v_1$ is the initial velocity expressed in km∙s⁻¹.

4.1. How to Start

The starting equations for the evolution of the SB [15] [21] are defined by the following parameters: $N^{\mathrm{<mo>\; *\; </mo>}}$, which is the number of SN explosions in 5.0 × 10⁷ yr; $Z_{\text{OB}}$, which is the distance of the OB associations from the galactic plane; $E_{51}$,

Figure 2. Numerical solution (full line) and Padé approximant (dashed line). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $b=2.8\text{pc}$, $t=194285\text{yr}$ and $t_0=680\text{yr}$.

which is the energy in 10⁵¹ erg and is usually chosen equal to one; $v_0$, which is 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, which is the proper time of the SB. The radius of the SB is

$R=111.56{\left(\frac{E_{51}{t}_7^3{N}^{\mathrm{<mo>\; *\; </mo>}}}{n_0}\right)}^{\frac{1}{5}}\text{pc}\mathrm{,}$ (35)

and its velocity

$V=6.567\frac{1}{t_7^{2/5}}\sqrt[5]{\frac{E_{51}{N}^{\mathrm{<mo>\; *\; </mo>}}}{n_0}}\text{km}/\text{s}\mathrm{.}$ (36)

In the following, we will assume that the bursting phase ends at $t=t_{7,0}$ (the bursting time is expressed in units of 10⁷ yr) when $N_{SN}$ SNs are exploded

$N_{SN}=N^{\mathrm{<mo>\; *\; </mo>}}\frac{t_{\mathrm{7,0}}\times {10}^7}{5\times {10}^7}\mathrm{.}$ (37)

The two following inverted formula allow us to derive the parameters of the initial conditions for the SB in terms of $r_0$ expressed in pc and $v_0$ expressed in km∙s⁻¹

$t_{\mathrm{7,0}}=0.05878\frac{r_0}{v_0}\mathrm{,}$ (38)

and

$N^{\mathrm{<mo>\; *\; </mo>}}=2.8289\times {10}^{-7}\frac{r_0^2{n}_0{v}_0^3}{E_{51}}\mathrm{.}$ (39)

4.2. The Astronomical Data

The LB has been recently observed in the X-ray in the 0.1 - 1.2 keV region by [3], whose Figure 7 reports six configurations of the LB along great-circle cuts through the Galactic pole and the Galactic plane. As a target of the simulation, we have chosen the cut characterized by galactic longitude, l, between 120˚ and 300˚. An observational percentage reliability, ${ϵ}_{\text{obs}}$, is introduced over the whole range of the polar angle $\theta$,

${ϵ}_{\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{,}$ (40)

where $r_{\text{num}}$ is the theoretical radius of the local LB, $r_{\text{obs}}$ is the observed radius of the local LB, 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.

4.3. The Results

The numerical solution is reported as a cut in the x-z plane: see Figure 3 for an exponential density profile as given by Equation (2); see Figure 4 for a Gaussian

Figure 3. Geometrical section of the LB in the x-z plane with an exponential profile (green points) and observed profile (red stars). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $b=3.5\text{pc}$, $t=8\times {10}^4\text{yr}$ and $t_0=80\text{yr}$. The observational percentage reliability is ${ϵ}_{\text{obs}}=\mathrm{81.93\; <mi>\; \%\; </mi>}$ and $N_{SN}=15.61$.

Figure 4. 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=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $b=5.83\text{pc}$, $t=1.35\times {10}^5\text{yr}$ and $t_0=1.35\times {10}^2\text{yr}$. The observational percentage reliability is ${ϵ}_{\text{obs}}=\mathrm{82.04\; <mi>\; \%\; </mi>}$ and $N_{SN}=15.61$.

density profile as given by Equation (5); see Figure 5 for an inverse square density profile as given by Equation (9); and see Figure 6 for a NFW density profile as given by Equation (13).

The theory of the asymmetrical expansion already developed is independent of the azimuthal angle $\phi$ and therefore the 3D advancing surface of a LB can be obtained by rotating a cut in x-z plane, see Figure 7.

4.4. ⁶⁰Fe-Signature

Some radioisotopes on Earth, such as ⁶⁰Fe (half life of 1.5 × 10⁶ yr [22]), were measured in a deep-sea ferromanganese crust: the concentration of ⁶⁰Fe increased 2.8 Myr ago, see [6]. These measurements have triggered some simulations that can explain the LB in the framework of SN explosions [12] [23] [24]. The encounter between the advancing shell of the LB and the Sun is here simulated in 2D assuming a constant density, see Equation (18). The following distances are involved:

1) $r_0$ the initial radius of the LB,

2) $r_e$ the radius of the LB when encounters the LB,

3) $r_a$ the actual radius of the LB,

4) D the distance between the sun and the LB, $D=r_a-r_e$, and they are reported in Figure 8. The times of the 2D simulation are

Figure 5. Geometrical section of the LB in the x-z plane with an exponential profile (green points) and observed profile (red stars). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $z_0=7\text{pc}$, $t=8\times {10}^4\text{yr}$ and $t_0=80\text{yr}$. The observational percentage reliability is ${ϵ}_{\text{obs}}=\mathrm{78.02\; <mi>\; \%\; </mi>}$ and $N_{SN}=15.61$.

Figure 6. Geometrical section of the LB in the x-z plane with a NFW profile (green points) and observed profile (red stars). The parameters are $v_0=3700\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $b=2.8\text{pc}$, $t=6.8\times {10}^5\text{yr}$ and $t_0=680\text{yr}$. The observational percentage reliability is ${ϵ}_{\text{obs}}=\mathrm{82.69\; <mi>\; \%\; </mi>}$ and $N_{SN}=15.61$.

Figure 7. 3D surface of the LB with parameters as in Figure 6, NFW profile. The three Euler angles are $\Theta =41$, $\Phi =-41$ and $\Psi =41$.

1) $t_0$ the time at which the radius of the LB is $r_0$,

2) $t_{{}^{60}\text{F}\text{e}}$ the time at which ⁶⁰Fe was deposited on the Earth,

3) $t_a$ the actual time of the LB,

4) $t_e$ the time of the encounter between LB and Sun, $t_e=t_a-t_{{}^{60}\text{F}\text{e}}$.

Figure 8. A sketch of the LB sun encounter. The black circle is the initial radius, the blue circle is the radius when the Sun is reached and green circle is the actual radius.

The distance LB-Sun, D, is reported in Figure 9 as function of the elapsed time.

4.5. Collective Effects

The LB is a part of other bubbles which show a Swiss-cheese structure, see Figure 10. We simulate this network with the multiple explosion of N bubbles in 2D assuming a constant density, see Equation (18). We choose $N=7$ and the time is allowed to vary in a random way in the interval $\left(t_0,t_{\max }\right)$, the position of the explosion on the two Cartesian axis is randomly generated in the interval $\left(\mathrm{0,}side\right)$, see Figure 11.

5. Conclusions

Two factors allow the comparison of different models which simulate the LB: the observational percentage reliability, see Equation (40); and acceptable observational cuts of the LB, see [3]. The best result is obtained adopting the NFW profile with a percentage reliability of ${ϵ}_{\text{obs}}=\mathrm{82.69\; <mi>\; \%\; </mi>}$. Similar results are obtained in the framework of the magnetic field model, see Figure 2 in [13]. The ⁶⁰Fe-signature is compatible with the model that we have developed here and Figure 9 reports the distance between the Sun and the LB. A simulation of the

Figure 9. Distance between Sun and the LB as function of time. The parameters of the LB are $v_0=10000\text{km}\cdot {\text{s}}^{-1}$, $r_0=5\text{pc}$, $t_{\max }=3\times {10}^6\text{yr}$, $t_{{}^{60}\text{F}\text{e}}=2.8\times {10}^6\text{yr}$ and $t_0=100\text{yr}$.

Figure 10. Map of the galactic environment of the sun with side of ≈528.12 pc.

Figure 11. Network of seven bubbles, $v_0=5000\text{km}\cdot {\text{s}}^{-1}$, $r_0=7\text{pc}$, $t_{\max }=3\times {10}^6\text{yr}$, $t_0=100\text{yr}$ and $side=400\text{pc}$.

exploding bubbles is reported in Figure 11. A more precise simulation of the exploding bubbles can be done when more detailed observations of the network, such as that reported in Figure 10, are available.

Acknowledgements

Credit for Figure 10 is given to the University of Bologna, see https://www.sslmit.unibo.it/zat/images/cartography/M-Way_2.htm.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Arnaud, M., Rothenflug, R. and Rocchia, R. (1984) The Local Hot Bubble from X-Ray Spectroscopic Measurements. Physica Scripta, 1984, 48. https://doi.org/10.1088/0031-8949/1984/T7/007
[2]	Slavin, J.D. (2016) Structures in the Interstellar Medium Caused by Supernovae: The Local Bubble. In: Alsabti, A.W. and Murdin, P., Eds., Handbook of Supernovae, Springer International Publishing, Cham, 1-13.
[3]	Liu, W., Chiao, M., Collier, M.R., et al. (2017) The Structure of the Local Hot Bubble. The Astrophysical Journal, 834, 33. https://doi.org/10.3847/1538-4357/834/1/33
[4]	Gontcharov, G.A. and Mosenkov, A.V. (2019) Interstellar Polarization and Extinction in the Local Bubble and the Gould Belt. MNRAS, 483, 299-314. https://doi.org/10.1093/mnras/sty2978
[5]	Medan, I. and Andersson, B.G. (2019) Magnetic Field Strengths and Variations in Grain Alignment in the Local Bubble Wall. The Astrophysical Journal, 873, 87. https://doi.org/10.3847/1538-4357/ab063c
[6]	Knie, K., Korschinek, G., Faestermann, T., Dorfi, E.A., Rugel, G. and Wallner, A. (2004) 60Fe Anomaly in a Deep-Sea Manganese Crust and Implications for a Nearby Supernova Source. Physical Review Letters, 93, Article ID: 171103. https://doi.org/10.1103/PhysRevLett.93.171103
[7]	Smith, R.K. and Cox, D.P. (1998) Modeling the Local Bubble Using Multiple Supernova Remnants. Proceedings of the IAU Colloquium, No. 166, Vol. 506, Garching, 133-136.
[8]	Smith, R.K. and Cox, D.P. (2001) Multiple Supernova Remnant Models of the Local Bubble and the Soft X-Ray Background. The Astrophysical Journal Supplement Series, 134, 283-309. https://doi.org/10.1086/320850
[9]	Welsh, B.Y., Sallmen, S. and Lallement, R. (2002) New Results from FUSE: A Paradigm for Testing Models of the Local Hot Bubble. American Astronomical Society Meeting, Vol. 200, 767-778.
[10]	Breitschwerdt, D., de Avillez, M.A. and Baumgartner, V. (2009) Modeling the Local Warm/Hot Bubble. American Institute of Physics Conference Series, Vol. 1156, 271-279. https://doi.org/10.1063/1.3211826
[11]	De Avillez, M.A. and Breitschwerdt, D. (2012) Non-Equilibrium Ionization Modeling of the Local Bubble. I. Tracing Civ, Nv, and Ovi Ions. Astronomy & Astrophysics, 539, L1. https://doi.org/10.1051/0004-6361/201117172
[12]	Schulreich, M., Breitschwerdt, D., Feige, J. and Dettbarn, C. (2018) A Way Out of the Bubble Trouble?—Upon Reconstructing the Origin of the Local Bubble and Loop I via Radioisotopic Signatures on Earth. Galaxies, 6, 26. https://doi.org/10.3390/galaxies6010026
[13]	Alves, M.I.R., Boulanger, F., Ferrière, K. and Montier, L. (2018) The Local Bubble: A Magnetic Veil to Our Galaxy. Astronomy & Astrophysics, 611, L5. https://doi.org/10.1051/0004-6361/201832637
[14]	Navarro, J.F., Frenk, C.S. and White, S.D.M. (1996) The Structure of Cold Dark Matter Halos. The Astrophysical Journal, 462, 563-575. https://doi.org/10.1086/177173
[15]	Dyson, J.E. and Williams, D.A. (1997) The Physics of the Interstellar Medium. Institute of Physics Publishing, Bristol. https://doi.org/10.1887/075030460X
[16]	Padmanabhan, P. (2001) Theoretical Astrophysics. Vol. II: Stars and Stellar Systems. Cambridge University Press, Cambridge.
[17]	Padé, H. (1892) Sur la représentation approchée d’une fonction par des fractions rationnelles. Annales scientifiques de l’école normale supérieure, 9, 193. https://doi.org/10.24033/asens.378
[18]	Wynn, P. (1966) Upon Systems of Recursions Which Obtain among the Quotients of the Padé Table. Numerische Mathematik, 8, 264. https://doi.org/10.1007/BF02162562
[19]	Baker, G. (1975) Essentials of Padé Approximants. Academic Press, New York.
[20]	Olver, F.W.J., Lozier, D.W., Boisvert, R.F. and Clark, C.W. (2010) NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge.
[21]	McCray, R. and Kafatos, M. (1987) Supershells and Propagating Star Formation. The Astrophysical Journal, 317, 190. https://doi.org/10.1086/165267
[22]	Rugel, G., Faestermann, T., Knie, K., Korschinek, G., Poutivtsev, M., Schumann, D., Kivel, N., Günther-Leopold, I., Weinreich, R. and Wohlmuther, M. (2009) New Measurement of the Fe60 Half-Life. Physical Review Letters, 103, Article ID: 072502. https://doi.org/10.1103/PhysRevLett.103.072502
[23]	Feige, J., Breitschwerdt, D., Wallner, A., Schulreich, M.M., Kinoshita, N., Paul, M., Dettbarn, C., Fifield, L.K., Golser, R., Honda, M., Linnemann, U., Matsuzaki, H., Merchel, S., Rugel, G., Steier, P., Tims, S.G., Winkler, S.R. and Yamagata, T. (2017) The Link between the Local Bubble and Radioisotopic Signatures on Earth. JPS Conference Proceedings, 14, Article ID: 010304. https://doi.org/10.7566/JPSCP.14.010304
[24]	Schulreich, M.M., Breitschwerdt, D., Feige, J. and Dettbarn, C. (2017) Numerical Studies on the Link between Radioisotopic Signatures on Earth and the Formation of the Local Bubble. I. 60Fe Transport to the Solar System by Turbulent Mixing of Ejecta from Nearby Supernovae into a Locally Homogeneous Interstellar Medium. Astronomy & Astrophysics, 604, A81. https://doi.org/10.1051/0004-6361/201629837