The Milky Way as modeled by Percolation and Superbubbles
PPercolation and Superbubbles 1
Chapter 1 T HE M ILKY W AY AS MODELED BY P ERCOLATIONAND S UPERBUBBLES
Lorenzo Zaninetti ∗ Physics Department, via P.Giuria 1,I-10125 Turin,Italy
PACS
Keywords:
Characteristics and properties of the Milky Waygalaxy; Supernova remnants
Abstract
The spiral structure of the Milky Way can be simulated by adopting percolation theory,where the active zones are produced by the evolution of many supernova (SN). Here weassume conversely that the percolative process is triggered by superbubbles (SB), the resultof multiple SNs. A first thermal model takes into account a bursting phase which evolvesin a medium with constant density, and a subsequent adiabatic expansion which evolves ina medium with decreasing density along the galactic height. A second cold model followsthe evolution of an SB in an auto-gravitating medium in the framework of the momentumconservation in a thin layer. Both the thermal and cold models are compared with the resultsof numerical hydro-dynamics. A simulation of GW 46.4+5.5, the Gould Belt, and the ∗ Email address: [email protected] a r X i v : . [ a s t r o - ph . GA ] O c t ZaninettiGalactic Plane is reported. An elementary theory of the image, which allows reproducingthe hole visible at the center of the observed SB, is provided.
1. Introduction
The term supershell was observationally defined by [1], where eleven H I objects were ex-amined. Supershells have been observed as expanding shells, or holes, in the H I-columndensity distribution of our galaxy, in the Magellanic Clouds, in the dwarf irregular galaxyHo-II [2], and in the nearby dwarf galaxy IC 2574 [3]. The dimensions of these objects spanfrom 100 pc to 1700 pc and often present elliptical shapes or elongated features, which aredifficult to explain based on an expansion in a uniform medium. These structures are com-monly explained through introducing theoretical objects named bubbles or superbubbles(SB); these are created by mechanical energy input from stars (see for example [4]; [5]).Thus, the origin of a supershell is not necessarily an SB, and other possible origins includecollisions of high-velocity clouds (see for example [6]; [7]). A worm is another observedfeature that may, or may not, be associated with a wall of an SB. Galactic worms were firstidentified as irregular, vertical columns of atomic gas stretching from the galactic plane;now, similar structures are found in radio continuum and infrared maps (see for example[8]). SBs are also observed in other galaxies, we cite two examples among many: a kine-matical and photometric catalog was compiled for 210 H II regions in the Large MagellanicCloud [9] and 12 SBs were classified in the irregular galaxy NGC 1569 [10].The models that explain SBs as being due to the combined explosions of super-nova in a cluster of massive stars are now briefly reviewed. Semi-analytical and hydro-dynamical calculations are generally adopted. In semi-analytical calculations, the thin-shell approximation can be the key to obtaining the expansion of the SB; see, for example,[11, 12, 13, 14, 15].The thin-shell approximation has already been used in a variety of different problems,and with both analytical and numerical approaches (besides the review by [16], see [17]and [18]). Thus the validity and limitations of the method are well known. For instance,modeling is fast and simple because the shell dynamics are only included in an approximateway, while the fluid variables are not included at all. The price that one has to pay is a lackof knowledge of the density and velocity profiles and, obviously, neither the onset anddevelopment of turbulence, instabilities, or mixing can be followed. These facts limit theapplications of this method to deriving only the general shapes, approximate expansionrates, and gross features of the surface mass-density distributions of the shells.The hydro-dynamical approximation was used by [19]. As for the effect of magneticfields, a semi-analytical method was introduced by [20] and a magneto-hydrodynamic codehas been adopted by various authors [21, 22, 23]. The often adopted exponential and verticalprofiles in density do not correspond to a physical process of equilibrium. The case ofan isothermal self-gravitating disk (ISD) is an equilibrium vertical profile which can becoupled with momentum conservation. The models cited leave some questions unansweredor only partially answered: • Is it possible to calibrate the vertical profile of an ISD?ercolation and Superbubbles 3 • Is it possible to deduce an analytical formula for the temporal evolution of an SB inthe presence of a vertical profile density as given by an ISD? • Is it possible to deduce numerical results for an SB when the expansion starts at agiven galactic height? • What is the influence of galactic rotation on the temporal evolution of an SB? • Can we explain the worms as a particular effect using image theory applied to SBs?In order to answer these questions, Section 2. reviews the percolation theory whichallows modeling a spiral galaxy. Section 3. introduces a numerical code which solves themomentum equation coupled with the variation of pressure in the presence of the injectionof mechanical luminosity and the adiabatic losses, the so called thermal model . Section 4.models an aspherical expansion, in a vertical profile in the number of particles as given byan ISD, the so called cold model . Section 5. reports a comparison between the numericalhydro-dynamical calculations and the two codes here presented. Section 6. models theSBs associated with GW 46.4+5.5, the Gould Belt, and the Galactic Plane. Section 7.contains detailed information on how to build an image of an SB both in the symmetric andasymmetric case.
2. The Percolation
The appearance of arms can be simulated through percolation theory [24, 25, 26, 27, 28].The fundamental hypotheses and the parameters adopted in the simulation will now bereviewed.1. The motion of a gas on the galactic plane has a constant rotational velocity, denotedby V G (in the case of spiral type Sb 218 km s − ). Here the velocity, V G , is expressedin units of 200 km s − and will therefore be V G =1.09.2. The polar simulation array made by rings and cells has a radius R G = 12 kpc. Thenumber of rings, (59), by the multiplication of R G with the number of rings foreach kpc, denoted by nring/kpc, which in our case is 5. Every ring is then made up ofmany cells , each one with a size on the order of the galactic thickness, ≈ p sp (for example 0.01), allows the process to start (with the previ-ous parameters, 111 new clusters were generated). Each one of these sources has sixnew surroundings that are labeled for each ring.4. In order to better simulate the decrease of the gas density along the radius, a stimu-lated probability of forming new clusters with a linear dependence by the radius , p s t = a + bR , was chosen. The values a and b are found by fixing prmax ( for exam-ple 0.18), the stimulated probability at the outer ring, and prmin (for example 0.24)in the inner ring; of course, prmin ≥ prmax . This approach is surprisingly similar to ZaninettiFigure 1. A typical simulation of the galaxy’s arms. Here ‘active’ is the final number of ac-tive cells, ‘random’ the initial number of randomly distributed sources, ‘pc’ the consideredprobability, ‘cells’ the total number of cells, R the radius in Kpc, ‘ring/kpc’ the numberof rings considered per kpc, the size of the cells is in kpc, v the velocity in units of 200Km/sec, and the time is expressed in units of millions of years, see the numerical values inTable 1the introduction of an anisotropic probability distribution in order to better simulatecertain classes of spirals [29].5. Now, new sources are selected in each ring based on the hypothesis of different stim-ulated probabilities. A rotation curve is imposed so that the array rotates in the samemanner as the galaxy. The procedure repeats itself n times (100); we denote by t G theage of the simulation, viz., · yr, where yr is the astrophysical counterpartof one time step.6. In order to prevent catastrophic growth, the process is stopped when the number ofsurroundings is greater than max (1000) and restarts by spontaneous probability.7. The final number of active cells (3824) is plotted with the size, which decreaseslinearly with the cluster age. In other words, the young clusters are bigger than theold ones. Only ten cluster ages are shown; only cells with an age of less than life (inour case · yr) are selected. A more deterministic approach has been followedby [30], where a model that connects the energy injection by star formation with theresulting interstellar structures in a differentially rotating disc has been introduced.A typical run is shown in Figure 2. and in Table 1 we summarize the parameters usedin this simulation.A real spiral galaxy as observed in the infrared is reported in Figure 2 as observedby the Spitzer Space Telescope, more details are available at ercolation and Superbubbles 5Table 1. Parameters of the model galactic radius kpcrings/kpc stochastic probability . timestep yearscells active Figure 2. Composite image in the infrared of M100. caltech.edu/ .
3. The thermal model for SBs
The starting equation for the evolution of the SB [11, 12, 31] is momentum conservationapplied to a pyramidal section, characterized by a solid angle ∆Ω j : ddt (cid:16) ∆ M j ˙ R j (cid:17) = pR j ∆Ω j , (1)where the pressure of the surrounding medium is assumed to be negligible and the mass isconfined in a thin shell with mass ∆ M j . The subscript j was added here in order to notethat this is not a spherically symmetric system. Due to the fact that p is uniform in the cavityof the SB, a summation to obtain the total volume, such as V = (cid:88) ∆Ω j R j / , (2) Zaninettiis necessary to determine its value. The mass conservation equation for a thin shell is ∆ M j = 13 R j ¯ ρ j ∆Ω j . (3)The pressure is enclosed in the energy conservation equation, γ − d ( pV ) dt = L − p dVdt , (4)where L = E R SN / π denotes the mechanical luminosity injected into a unit solid angleand γ = 5 / . Equation (4) can be expanded, obtaining dpdt = L ( γ − V − γ pV dVdt . (5)Formulae (1) and (5) will be our basic equations to numerically integrate in the burstingphase. An approximation concerning the pressure can be obtained by ignoring the secondterm of the rhs in (4); this leads to p = 2 E R SN t V , (6)where t is the considered time, R SN the rate of supernova explosions, E the energy of eachsupernova, and V is computed as in Equation (2). On continuing to consider the case ofconstant density, the volume becomes V = 4 π R . (7)Equations (1) and (6) lead to ddt (cid:16) ¯ ρR ˙ R (cid:17) = 3 E R SN π tR . (8)In order to integrate the above equation, R ˙ R = AR α is imposed. After adopting the initialcondition of R = 0 at t = 0 and assuming ¯ ρ is constant irrespective of R or t , the equationfor the expansion speed is obtained: ˙ R = α + 1 α E R SN π ¯ ρR t . (9)By integrating the previous equation, “the expansion law” is obtained: R = (cid:20) α + 1)4 πα (cid:21) / (cid:18) E R SN ¯ ρ (cid:19) / t / , (10)which is identical to equation (10.34) of [11], since α = 7 / .ercolation and Superbubbles 7 It is clear that an upper limit should be inserted into the basic equation (6); this is the timeafter which the bursting phenomenon stops, t burst .Since there is a pdV term in the the First Law of Thermodynamics ( dQ = 0 = dU + pdV ), the total thermal energy decreases with time. The pressure of the internalgas decreases according to the adiabatic law, p = p ( V V ) / , (11)where p = 23 E R SN t burst V , (12)and the equation for the conservation of momentum becomes ddt (cid:16) ¯ ρ j R j ˙ R j (cid:17) = 3 p V / V / R j . (13)It is important to remember that this phase occurs after the SN burst stops. At that time( t = t burst ), the volume of the bubble is computed by using R j ( t burst ) , and the expansionspeed is equal to ˙ R j ( t burst ) , both of which are obtained from the numerical solution ofEquation (8). On assuming (also here) that ¯ ρ is constant irrespective of R or t , Equation (12)is now inserted into Equation (1) and the expansion law is obtained for the after-burst phase, R = (cid:18) π (cid:19) / (cid:32) E R SN t burst R ¯ ρ (cid:33) / t / , (14)which is identical to equation (10.33) of [11]. The differential equations need to be solved by sectors, with each sector being treated asindependent from the others, except for the coupling in computing the volume of the bubble,see Equation (2). The range of the polar angle θ ( ◦ ) will be divided into n θ steps, andthe range of the azimuthal angle φ ( ◦ ) into n φ steps. This will yield ( n θ +1) ( n φ +1)directions of motion that can also be identified with the number of vertices of the polyhedronrepresenting the volume occupied by the expansion; this polyhedron varies from a sphere tovarious morphological shapes based on the swept up material in each direction. In 3D plotsshowing the expansion surface of the explosion, the number of vertices is n v = ( n θ +1) · ( n φ +1) and the number of faces is n θ · n φ , typically we have n θ =50 and n φ =50 (in this case thesubscript varies between 1 and 2601). However, all calculated models are axisymmetric,and the essential number of points to draw such figures is only n θ + 1 =51. R up , R eq and R down are now introduced, which represent the distances from the position of the the OBassociations (denoted by z OB ) to the top, to the left and to the bottom of the bubble. Zaninetti At each time step, ∆ t , the volume V of the expanding bubble is computed, see Equation(2). In other words, the volume V swept up from the explosions is no longer a sphere,but becomes an egg or an hourglass. The pressure in the first phase, see Equation (5), iscomputed through the following finite-difference approximation: p k = p k − + (cid:34) L γ − V k − γ p k − V k V k − V k − ∆ t (cid:35) ∆ t , (15)where k is the number of steps employed.Equation (1) now leads to ddt (cid:16) ¯ ρ j R j ˙ R j (cid:17) = 3 pR j , (16)which may rewritten as R j ˙ R j + R j ¨ R j = 3 p R j ¯ ρ j − ˙¯ ρ j ¯ ρ j R j ˙ R j . (17)The first term represents the ram pressure of the stratified ISM on the expanding surfaceand the second represents the inertia of the bubble. The average density ¯ ρ j is numericallycomputed according to the algorithm outlined in subsection 3.5.; the time derivative of ¯ ρ j at each time step ∆ t is computed according to the finite difference-approximation, ˙¯ ρ j = ¯ ρ jk − ¯ ρ jk − ∆ t , (18)where k is the number of steps considered.Equation (17) can be re-expressed as two differential equations (along each direction j )of the first order suitable to be integrated: dy ,j dt = y ,j , (19) dy ,j dt = 3 p ρ j y ,j − y ,j y ,j − ˙¯ ρ j ¯ ρ j y ,j . (20)In the new regime ( t ≥ t burst ), Equation (13) now becomes R j ˙ R j + R j ¨ R j = 3 p V / V / ¯ ρ j R j − ˙¯ ρ j ¯ ρ j R j ˙ R j . (21)Equation (21) can be re-expressed as two differential equations (along each direction j ) ofthe first order suitable to be integrated: dy ,j dt = y ,j , (22) dy ,j dt = 3 p V / V / ¯ ρ j y ,j − y ,j y ,j − ˙¯ ρ j ¯ ρ j y ,j . (23)ercolation and Superbubbles 9Figure 3. The average structure of the gaseous disk in the z -direction; z is allowed to varybetween 0 pc and 400 pc.The integrating scheme used is the Runge–Kutta method, and in particular the sub-routine RK4 ([32]); the time derivative of the density along a certain direction j , ˙¯ ρ j , iscomputed at each time step according to formula (18).The integration time, t age , and the time steps, are always indicated in the captions of thevarious diagrams. The pressure across the bursting time is continuous. This is because thefirst value of the pressure after the bursting time is that of the last step in the bursting phase,modified according to the adiabatic law modelled by formula (11).The upper limit chosen to integrate the differential equations is t age = 2 . · yr ; thisis the approximate age of GSH 238. The vertical density distribution of galactic H I is well-known; specifically, it has the fol-lowing three component behavior as a function of z , the distance from the galactic plane inpc: n ( z ) = n e − z /H + n e − z /H + n e −| z | /H . (24)We took [33, 34, 35] 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. This distribution of galactic H I isvalid in the range 0.4 ≤ R ≤ R , where R = 8.5 kpc and R is the distance from the centerof the galaxy. A plot showing such a dependence of the ISM density on z is shown in Figure3. The ISM density is not constant, but varies in the z -direction according to Equation (24).The swept mass can therefore be computed in a certain direction j using the followingalgorithm:1. The pyramidal sector is divided into layers (for example 1000) whose radii rangefrom R L − / to R L +1 / .0 Zaninetti2. In each layer the volume ∆ V L = ( R L +1 / − R L − / ) is computed as well as thecorresponding mass, ∆ M L = ∆ V L ρ j ( z = R L cos θ j + z OB ) , where θ j represents anangle between the z -axis and the j th radial path.3. The various contributions, ∆ M L , are added in order to obtain the total swept mass inthe considered sector.4. The average density along a sector j as ¯ ρ j is calculated using ∆ M L and ∆ V L as ¯ ρ j = (cid:80) L ∆ M L / (cid:80) L ∆ V L . Our basic units are: time ( t ), which is expressed in units of yr; E , the energy in erg; n , the density expressed in particles cm − (density ρ = n m, where m=1.4 m H );and N ∗ , which is the number of SN explosions in . · yr.By using the previously defined units, formula (10) concerning “the expansion law” inthe bursting phase becomes R = 111 .
56 pc( E t N ∗ n ) . (25)Conversely, Equation (14) concerning “the expansion law” in the adiabatic phase is R = 171 .
85 pc( E N ∗ n ) ( t burst7 ) ( t ) . (26)This is the approximate radius derived from a spherical model. A perfect coincidence is notexpected upon making a comparison with that expected in a stratified medium obtained withthe multidimensional thin shell approximation. Another useful formula is the luminosity ofthe bubble, see [11], L SN = E R SN = 0 . erg s − E N ∗ . (27)This formula can be useful in order to derive the parameter N ∗ from observations; forexample, a luminosity of . erg s − corresponds to N ∗ = 248 . The total depositedenergy, E tot, is E tot = E N ∗ t age · yr erg , (28)when t burst = t age and E tot = E N ∗ t burst · yr erg , (29)when t burst < t age .The spectrum of the radiation emitted depends on the temperature behind the shockfront, see for example formula 9.14 in [36], T = 316 µk v s K ◦ , (30)ercolation and Superbubbles 11where µ is the mean mass per particle, k the Boltzmann constant and v s the shock velocityexpressed in cm sec − . A formula which is useful for the implementation in the code iseasily derived, T = 31 . v sk K ◦ , (31)when v sk is expressed in km sec − .
4. The Cold model for SBs
This Section reviews the standard thin layer approximation and then introduces two recur-sive equations for the dynamical evolution of a SB in an auto-gravitating medium, see [37].
The thin layer approximation assumes that all the swept-up gas accumulates infinitely in athin shell just after the shock front. The conservation of radial momentum requires that πR ρ ˙ R = M , (32)where R and ˙ R are the radius and the velocity of the advancing shock, ρ is the density ofthe ambient medium, M is the momentum evaluated at t = t , R is the initial radius, and ˙ R the initial velocity, see [38, 39]. The law of motion is R = R (cid:32) R R ( t − t ) (cid:33) , (33)and the velocity ˙ R = ˙ R (cid:32) R R ( t − t ) (cid:33) − . (34) Given the Cartesian coordinate system ( x, y, z ) , the plane z = 0 will be called the equatorialplane, z = R sin( θ ) , where θ is the latitude angle which has range [ − ◦ ↔ +90 ◦ ] ,and R is the distance from the origin. The latitude angle is often used in astrophysics tomodel asymmetries in the polar lobes, see the example of the nebula around η -Carinae(Homunculus) shown in Table 1 in [40]. In our framework, the polar angle of the sphericalcoordinate system is − θ . The vertical number density distribution of galactic H I isusually modeled by the three component function as given by eqn. (24). Here, conversely,we adopt the density profile of a thin self-gravitating disk of gas which is characterized bya Maxwellian distribution in velocity and distribution which varies only in the z -direction(ISD). The number density distribution is n ( z ) = n sech ( z h ) , (35)2 ZaninettiFigure 4. Profiles of density versus scale height z : the self-gravitating disk as given byEq. (35) when h = 90 pc (dashed) and the three-component exponential distribution asgiven by Eq. (24) (full line).where n is the density at z = 0 , h is a scaling parameter, and sech is the hyperbolic secant([41, 42, 43, 44]).Figure (4) compares the empirical function sum of three exponential disks and thetheoretical function as given by Eq. (35). Assuming that the expansion starts at z = 0 , wecan write z = R sin( θ ) , and therefore n ( R, θ ) = n sech ( R sin( θ )2 h ) , (36)where R is the radius of the advancing shell.The 3D expansion that starts at the origin of the coordinates will be characterized bythe following properties. • The dependence of the momentary radius of the shell on the latitude angle θ over therange [ − ◦ ↔ +90 ◦ ] . • The independence of the momentary radius of the shell from φ , the azimuthal anglein the x - y plane, which has a range [0 ◦ ↔ ◦ ] .The mass swept, M , along the solid angle ∆ Ω between 0 and R is M ( R, θ ) = ∆ Ω3 m H n I m ( R ) + 43 πR n m H , (37)where I m ( R ) = (cid:90) RR r sech ( r sin( θ )2 h ) dr , (38)where R is the initial radius and m H the mass of hydrogen. The integral is I m ( R ) = − hr (sin( θ )) − (1 + e r sin( θ ) h ) − + 4 hr sin( θ ) − h r ln(1 + e r sin( θ ) h )(sin( θ )) − − h Li ( − e r sin( θ ) h )(sin( θ )) − , (39)ercolation and Superbubbles 13where Li ( z ) = ∞ (cid:88) n =1 z n n , (40)is the dilogarithm, see [45, 46, 47].The conservation of momentum along the solid angle ∆ Ω gives M ( R, θ ) ˙ R ( θ ) = M ( R ) ˙ R , (41)where ˙ R ( θ ) is the velocity at R and ˙ R is the initial velocity at R = R . Using the previousequation, an analytical expression for ˙ R ( θ ) along the solid angle can be found, but it iscomplicated, and therefore we omit it. In this differential equation of the first order in R ,the variables can be separated and integrating term by term gives (cid:90) RR M ( r, θ ) dr = M ( R ) ˙ R ( t − t ) , (42)where t is the time and t the time at R . We therefore have an equation of the type F ( R, R , h ) NL = 13 R ˙ R ( t − t ) , (43)where F ( R, R , h ) NL has an analytical but complicated form. The case of expansion thatstarts from a given galactic height z , denoted by z OB , which represent the OB associations,cannot be solved by Eq. (43), which is derived for a symmetrical expansion that startsat z = 0 . It is not possible to find R analytically and a numerical method should beimplemented. In our case, in order to find the root of the nonlinear Eq. (43), the FORTRANsubroutine ZRIDDR from [32] has been used.The following two recursive equations are found when momentum conservation is ap-plied: R n +1 = R n + V n ∆ tV n +1 = V n (cid:16) M n ( r n ) M n +1 ( R n +1 ) (cid:17) , (44)where R n , V n , M n are the temporary radius, the velocity, and the total mass, respectively, ∆ t is the time step, and n is the index. The advancing expansion is computed in a 3D Carte-sian coordinate system ( x, y, z ) with the center of the explosion at (0,0,0). The explosion isbetter visualized in a 3D Cartesian coordinate system ( X, Y, Z ) in which the galactic planeis given by Z = 0 . The following translation, T OB , relates the two Cartesian coordinatesystems. T OB X = xY = yZ = z + z OB , (45)where z OB is the distance in parsecs of the OB associations from the galactic plane.The physical units have not yet been specified: parsecs for length and yr for time,see also subsection 3.6., are an acceptable astrophysical choice. With these units, the initial4 Zaninettivelocity V = ˙ R is expressed in units of pc/( yr) and should be converted into km/s;this means that V = 10 . V where V is the initial velocity expressed in km/s.Analytical results can also be obtained solving the Kompaneets equation, see [48], forthe motion of a shock wave in different plane-parallel stratified media such as exponential,power-law type, and a quadratic hyperbolic-secant, see synoptic Table 4 in [49].
5. A Comparison with Hydrosimulations
This Section reports the comparison of the thermal and cold models with numerical hydro-dynamics calculations.
The shape of the SB in the thermal model, see Section 3., depends strongly on the timeelapsed since the first explosion, t age , the duration of SN burst, t burst , on the number of SNexplosions in the bursting phase, and on the adopted density profile. Ordinarily, the dynam-ics of an SB is studied with hydro simulations (partial differential equations). Conversely,here under the hypothesis called “thin-shell approximation”, the dynamics is studied bysolving ordinary differential equations. This choice allows the calculation of a much largernumber of models compared with hydro-dynamics calculations.The level of confidence in our results can be given by a comparison with numericalhydro-dynamics calculations see, for example, [19]. The vertical density distribution theyadopted [ see equation(1) from [19] and equation(5) from [50] ] has the following it zdependence, the distance from the galactic plane in pc: n hydro = n d (cid:40) Θexp[ − V ( z ) σ ] + (1 − Θ)exp (cid:34) − V ( z ) σ (cid:35)(cid:41) , (46)with the gravitational potential as V ( z ) = 68 . (cid:20) . (cid:18) . zz (cid:19)(cid:21) (kms − ) . (47)Here n d = 1 particles cm − , Θ = 0 . , σ IC = 14 . − , σ C = 7 . − and z = 124 pc . Table 2 reports the results of ZEUS (see [19]), a two-dimensional hydro-dynamic code, when t age = 0 . · yr . The supernova luminosity is . · erg s − , z OB =100 pc and the density distribution is given by formula (46). The ZEUS code wasoriginally described by [51].In order to make a comparison our code was run with the parameters of the hydro-code(see Figure 5); a density profile as given by formula (46) was adopted. The percentage ofreliability of our code can also be introduced, (cid:15) = (1 − | ( R hydro − R num ) | R hydro ) · , (48)where R hydro is the radius, as given by the hydro-dynamics, and R num the radius obtainedfrom our simulation. In the already cited table 2, our numerical radii can also be found in theercolation and Superbubbles 15Figure 5. Section of the on the x-z plane when the explosion starts at z OB = 100 pc. Thedensity law is given by Equation 46. The code parameters are t age = .
45 10 yr , ∆ t =0 . · yr , t burst7 =0.45, and N ∗ =248. The points represented by the small crossesindicate the inner section from Figure 3a of [19].Table 2. Code reliability. R up (pc) R down (pc) R eq (pc) R hydro ( ZEU S ) 330 176 198 R num ( our code ) 395 207 237 efficiency (%) 80 81 80 upward, downward and equatorial directions, and the efficiency as given by formula (48).The value of the radii are comparable in all of the three chosen directions and Figure 5 alsoreports on the data from Figure 3a of [19]. In the framework of the cold model developed in Section 4. Figure 6 compares the hydronumber density as given by Eq. 46 and the theoretical function as given by Eq. (35). Thedifference in the density profiles from hydrosimulations and from the cold model adoptedare due to the fact that the density at z = 0 is assumed to be n hydro = 1 particles cm − in the hydrosimulations, see [19]. In our cold model conversely at z = 0 we have n =0 . particles cm − as in [35].A typical run of ZEUS (see [19]), a two-dimensional hydrodynamic code, is done for t = 0 . and supernova luminosity of . · erg s − when z OB =100 pc. In order to make6 ZaninettiFigure 6. Profiles of density versus scale height z : the self-gravitating disk as in Eq. (35)when h = 90 pc (dashed line) and the hydro number density as given by Eq. 46 (full line).a comparison with our cold code, we adopt the same time and we search the parameterswhich produce similar results, see Figure 7.
6. Astrophysical SBs
This section evaluates the significance of the galactic rotation for the shape of the SBs,models the thermal GW 46.4+5.5, the thermal Gould Belt, the cold GW 46.4+5.5 and thethermal Galactic Plane.
The influence of the Galactic rotation on the results can be be obtained by introducing thelaw of the Galactic rotation as given by [52], V R ( R ) = 220( R [pc]8500 ) . km sec − , (49)where R is the radial distance from the center of the Galaxy expressed in pc. The transla-tion of the previous formula to the astrophysical units adopted gives V R ( R ) = 2244 (cid:18) R [pc]8500 (cid:19) . pc10 yr , (50) Ω( R ) = 2244 ( R [pc]8500 ) . R [pc] rad10 yr . (51)Here, Ω( R ) is the differential angular velocity and φ = 2244 ( R [pc]8500 ) . R [pc] t rad . (52)ercolation and Superbubbles 17Figure 7. Section of the SB in the X-Z plane when the explosion starts at z OB = 100 pc(empty stars). The cold code parameters are h = 90 pc, t = 0 . , t , = 0.00045, r = 24.43, V = 3191 km s − , N SN = 180 and N ∗ =2000000. The points represented bythe small crosses indicate the inner section from Figure 3a of MacLow et al. 1989. Theexplosion site is represented by a big cross.8 ZaninettiHere, φ , is the angle made on the circle and t the time expressed in units of yr . Uponconsidering only a single object, an expression for the angle can be found once R = R + y is introduced, φ ( y ) = V R ( R ) R + y t. (53)The shift in the angle due to differential rotation can now be introduced, ∆ φ = φ ( y ) − φ (0) , (54)where x and y are the SB coordinates in the inertial frame of the explosion, denoted by x =0, y =0 and z = z OB . The great distance from the center allows us to say that ∆ xR = ∆ φ, (55)where ∆ x is the shift due to the differential rotation in the x coordinate. The shift can befound from (55) and (54) once a Taylor expansion is performed, ∆ x ≈ − V R ( R ) yR t. (56)Upon inserting (50) in (56), the following transformation, T r , due to the rotation is obtainedfor a single object in the solar surroundings: T r x (cid:48) = x + 0 . y ty (cid:48) = yz (cid:48) = z, (57)where y is expressed in pc and t in units of yr . A careful study of the worms 46.4+5.5 and 39.7+5.7 [53] has led to the conclusion that theybelong to a single super-shell. Further on, the dynamical properties of this H I supershellcan be deduced by coupling the observations with theoretical arguments, [54]; the derivedmodel parameters to fit the observations are reported in Table 3, where the altitudes ofKK 99 3 and KK 99 4 (clouds that are CO emitters) have been identified with z OB by theauthor.These parameters are the input for our thermal computer code (see the captions of Fig-ure 8). The problem of assigning a value to z OB now arises, and the following two equationsare set up: R up + z OB = 540 ,R up + z OB R down − z OB = D up D down = 15 ◦ ◦ = 3 . (58)The algebraic system (58) consists of two equations and three variables. The value chosenfor the minimum and maximum latitudes ( ◦ and − ◦ ) are in rough agreement with theercolation and Superbubbles 19Table 3. Data of the supershell associated with GW 46.4+5.5.Size (pc ) 345 · Expansion velocity (km s − ) 15 Age (10 yr ) 0 . z OB (pc) Total energy ( erg) 15 Table 4. Radii concerning GW 46.4+5.5.
Direction R num (pc) R num (pc) with the Euler process R obs (pc)Equatorial 238 233 172 . position of the center at +5 ◦ of the galactic latitude (see [55]). One way of solving 58 is toset, for example, z OB =100 pc. The other two variables are easily found to be as follows: z OB = 100 pc ,R down = 235 pc , (59) R up = 305 pc . For this value of z OB , we have the case where a transition from an egg-shape to a V-shape isgoing on, and the simulation gives the exact shape. In order to obtain E tot = 15 . · erg (see formula (29)) and t burst = 0 . · yr , we have inserted N ∗ = 150.In order to test our simulation, an observational percentage of reliability is introducedthat uses both the size and the shape, (cid:15) obs = 100(1 − (cid:80) j | R obs − R num | j (cid:80) j R obs j ) , (60)where R obs is the observed radius, as deduced by using the following algorithm. The radiusat regular angles from the vertical ( ◦ , ◦ , ◦ ) is extracted from Table 3, giving theseries (305 pc, 172.5 pc, 235 pc); the theory of cubic splines [32] is then applied to computethe various radii at progressive angles ( ◦ , . ◦ , . ◦ , ... ◦ ), which is a series computedby adding regular steps of ◦ /n θ . The data are extracted from the dotted ellipse visiblein Figure 7 of [55]; this ellipse represents the super-shell at v LSR = 18 . km sec − .We can now compute the efficiency over 50+1 directions [ formula (60) ] of a section y-z when x =0 which turns out to be (cid:15) obs = 68 . ; the observed and numerical radii along thethree typical directions are reported in Table 4. The physical parameters adopted from [53]turn out to be consistent with our thermal numerical code. The results of the simulation canbe visualized in Figure 8, or through a section on the x-z plane (see Figure 9). Another0 ZaninettiFigure 8. Model of GW 46.4+5.5. The parameters are t age = . · yr , ∆ t = 0 . · yr , t burst7 = 0.5, N ∗ = 150, z OB =100 pc, and E =1. The three Eulerian angles characterizingthe point of view of the observer are Φ = 0 ◦ , Θ = 90 ◦ and Ψ = 0 ◦ .Figure 9. Section of the SB on the x-z plane when the physical parameters are the same asin Figure 8.ercolation and Superbubbles 21Figure 10. Map of the expansion velocity relative to a simulation of GW 46.4+5.5 when190000 random points are selected on the surface. The physical parameters are the same asin Figure 8 and the three Eulerian angles characterizing the point of view of the observerare Φ = 0 ◦ , Θ = 90 ◦ , and Ψ =0 ◦ .important observational parameter is the H I 21 cm line emission: the observation of H I gasassociated with GW 46.4+5 [56, 53] reveals that the super-shell has an expansion velocityof V exp ≈
15 km s − [53]. In order to see how our model matches the observations, theinstantaneous radial velocities are computed in each direction. A certain number of randompoints are then generated on the theoretical surface of expansion. The relative velocity ofeach point is computed by using the method of bilinear interpolation on the four grid pointsthat surround the selected latitude and longitude [32].Our model gives radial velocities of V theo , 31 km s − ≤ V theo ≤
71 km s − (27 km s − ≤ V theo ≤ km s − with the Euler process) and a map of the expansionvelocity is given in Figure 10, from which it is possible to visualize the differences in theexpansion velocities between the various regions.Perhaps it is useful to map the velocity ( V ptheo ) in the y -direction, V ptheo = v ( θ, φ ) · sin( θ ) · (sin φ ) . (61)This is the velocity measured along the line of sight when an observer stands in front ofthe super-bubble; θ and φ were defined in subsection 3.2.. The structure of the projectedvelocity, V ptheo , is mapped in Figure 11 by using different colors; the range is 0 km s − V ptheo , relative to a simulationof GW 46.4+5.5 when 190000 random points are selected on the surface. The physicalparameters are the same as in Figure 11. ≤ V ptheo ≤
36 km s − (0 km s − ≤ V ptheo ≤
32 km s − with the Euler process) andthe averaged projected velocity is ≈ km s − ( ≈ . km s − with the Euler process), avalue that is greater by ≈ km s − ( ≈ . km s − with the Euler process) than the alreadymentioned observed expansion velocity, V exp = 15 km s − .As is evident from the map in Figure 11, the projected velocity of expansion is not uni-form over all of the shell’s surface, but is greater in the central region than in the externalregion. In this particular case of an egg-shape, we observe a nearly circular region con-nected with the maximum velocities in the upper part of the shell. It is therefore possibleto speak of egg-shaped appearances in the Cartesian physical coordinates and spherical-appearances in the projected maximum velocity. The physical parameters concerning the Gould Belt as deduced in [57], are reported inTable 5. The total energy is such as to produce results comparable with the observationsand the kinematic age is the same as in [57, 58]. In order to obtain E tot = 6 . erg with t burst = 0 .
015 10 yr , we have inserted N ∗ = 2000. The time necessary to cross the Earth’sorbit, that lies 104 pc away from the Belt center, turns out to be .
078 10 yr , which means .
52 10 yr from now. The 2D cut at z =0 of the SB can be visualized in Figure 12. Ourercolation and Superbubbles 23Table 5. Data of the super-shell associated with the Gould BeltSize (pc ) 466 · at b=0Expansion velocity (km s − ) 17 Age (10 yr) 2 . Total energy (10 erg) 6 Figure 12. Rhombi represent the circular section, the stars the rotation-distorted sectionand the big star the Sun. The Galaxy’s direction of rotation is also shown. The parametersare t age = . · yr , ∆ t = 0 . · yr , t burst7 = 0.015, N ∗ = 2000, z OB =0 pc, and E =1.4 ZaninettiFigure 13. The stars represent the rotation-distorted section of the Gould Belt and the bigstar the Sun. The velocity field of the expansion modified by the shear velocity is mapped.The Galaxy’s direction of rotation is also shown.Table 6. Simulated data of the SB associated with GW 46.4+5.5.Size (pc ) 454 · Averaged expansion velocity (km s − ) 19 . Age (10 yr ) 0 . z OB (pc) thermal model gives a radial velocity at z =0 of V theo =3.67 km s − . The influence of theGalactic rotation on the direction and modulus of the field of the radial velocity is obtainedby an application of transformation (50), see Figure 13. A comparison should be made withfigure 5 and figure 9 in [57]. The analytical cold model as given by the solution of the nonlinear Eq. (43) can be used onlyin the case z OB =0. As an example, we give a model of the SB associated with GW 46.4+5.5,see Figure 14. The numerical solution as given by the recursive relation (44) can be foundadopting the same input data as the analytical solution, see the crosses in Figure (14); thenumerical solution agrees with the analytical solution within 0.65%. We are now ready topresent the numerical evolution of the SB associated with GW 46.4+5 when z OB = 100 pc,see Figure 15 and Table 6.ercolation and Superbubbles 25Figure 14. Section of the SB GW 46.4+5 in the Y-Z;X=0 plane when the explosion startsat z OB = 0 pc. The analytical results are shown as a full line and the numerical resultsas crosses. The cold code parameters of the solution of Eq. (43) as well of the numericalcouple (44) are h = 90 pc, t = 0 . , t , = 0.0045, r = 49.12, V = 641 . − , N SN = 93 and N ∗ =103000.Figure 15. Section of the SB GW 46.4+5 in the Y-Z;X=0 plane when the explosion startsat z OB = 100 pc. The cold code parameters for the numerical couple (44) are h = 90 pc , t = 0 . , t , = 0.00045, r = 24.43, V = 3191km s − , N SN = 180 and N ∗ =2000000.The explosion site is represented by a cross.6 ZaninettiFigure 16. Structure of the galactic plane in the Hammer–Aitoff projection, as resultingfrom the SB/percolation network. The value of pselect is 0.08, corresponding to 260 se-lected clusters. The parameters were ∆ t = 0 . · yr, t burst =0.5, N ∗ = 100 and each SBhad 200 random points on its surface. In order to simulate the structure of the galactic plane, the same basic parameters as inFigure 1 can be chosen, but now the SBs are drawn on an equal-area Aitoff projection. Inparticular, a certain number of clusters will be selected through a random process accordingto the following formula: selected clusters = pselect · number of clusters f rom percolation, (62)where pselect has a probability lower than one. The final result of the simulation for thethermal case is reported in Figure 16.
7. Theory of the Image
The transfer equation in the presence of emission only, see for example [59] or [60], is dI ν ds = − k ν ζI ν + j ν ζ , (63)where I ν is the specific intensity, s is the line of sight, j ν the emission coefficient, k ν a massabsorption coefficient, ζ the mass density at position s , and the index ν denotes the relevantfrequency of emission. The solution to Eq. (63) is I ν ( τ ν ) = j ν k ν (1 − e − τ ν ( s ) ) , (64)where τ ν is the optical depth at frequency νdτ ν = k ν ζds . (65)We now continue analysing the case of an optically thin layer in which τ ν is very small (or k ν very small ) and the density ζ is replaced with our number density C ( s ) of particles.One case is taken into account: the emissivity is proportional to the number density j ν ζ = KC ( s ) , (66)ercolation and Superbubbles 27where K is a constant. This can be the case for synchrotron emission. We select as anexample the [S II] continuum of the synchrotron SB in the irregular galaxy IC10, see [61],and the X-ray emission below 2 keV around the OB association LH9 in the H II complexN11 in the Large Magellanic Cloud, see [62]. The intensity at a given frequency is I ( ν ) ∝ lν β , (67)where l is the length of the radiating region along the line of sight. The synchrotron lumi-nosity is assumed to be proportional to the flow of kinetic energy L m , L m = 12 ρAV , (68)where A is the considered area, see formula (A28) in [63]. In our case, A = 4 πR , whichmeans L m = 12 ρ πR V , (69)where R is the instantaneous radius of the SB and ρ is the density in the advancing layerin which the synchrotron emission takes place. The astrophysical version of the the rate ofkinetic energy is L ma = 1 . × n R V ergss , (70)where n is the number density expressed in units of particlecm , R is the radius in parsecs,and V is the velocity in km/s. The spectral luminosity, L ν , at a given frequency ν is L ν = 4 πD S ν , (71)with S ν = S ( νν ) β , (72)where S is the flux observed at the frequency ν and D is the distance. The total observedsynchrotron luminosity, L tot , is L tot = (cid:90) ν max ν min L ν dν , (73)where ν min and ν max are the minimum and maximum frequencies observed. The totalobserved luminosity can be expressed as L tot = (cid:15)L ma , (74)where (cid:15) is a constant of conversion from the mechanical luminosity to the total observedluminosity in synchrotron emission.The fraction of the total luminosity deposited in a band f c is f c = ν c , min β +1 − ν c , max β +1 ν min β +1 − ν max β +1 , (75)where ν c,min and ν c,max are the minimum and maximum frequencies of the band. Table7. shows some values of f c for the most important optical bands. An analytical solution8 ZaninettiTable 7. Table of the values of f c when ν min = 10 Hz , ν max = 10 Hz and β = − . .band λ ( ˚ A ) FWHM ( ˚ A ) f c U 3650 700 6.86 × − B 4400 1000 7.70 × − V 5500 900 5.17 × − Hα × − [ SII ] continuum 7040 210 0.92 × − Figure 17. The two circles (sections of spheres) which include the region with constantdensity are represented by a full line. The observer is situated along the x direction, threelines of sight are indicated, and the relativistic electrons have radius r in the region a < r
8. Conclusions
Percolation and Spiral Galaxies
The active regions of the percolation theory for spiralsare usually connected with the SNRs but also the SBs can trigger the formation of new stars.
Thermal law of motion
The expansion of a super-bubble in the ISM belonging to ourgalaxy can be simulated by applying Newton’s second law to different pyramidal sectors.The following results were achieved:1. Single objects, like the super-shells associated with GW 46.4+5 and GSH 238, weresimulated with an efficiency of (cid:15) obs ≈ and (cid:15) obs ≈ , respectively.2. The network of many explosions that originate from the galactic plane could be ten-tatively simulated. Cold and asymmetrical law of motion
The temporal evolution of an SB in a mediumwith constant density is characterized by a spherical symmetry. The presence of a thin self-gravitating disk of gas which is characterized by a Maxwellian distribution in velocity anda distribution of density which varies only in the z -direction produces an axial symmetry inthe temporal evolution of an SB. The resulting Eq. (43) has an analytical form which canbe solved numerically when z OB =0. The case of z OB (cid:54) = 0 can be attached solving tworecursive equations, see (44). Astrophysical Applications
The thermal model was applied to GW 46.4+5.5, with an efficiency of (cid:15) obs = 68 . ,see Section 6.2., and to the Gould Belt, for which a map of the velocities as modified by thegalactic rotation is reported, see Figure 13. The cold model was applied to GW 46.4+5.5. Images
The combination of different processes such as the complex shape of SBs, the thin layerapproximation (which means an advancing layer having a thickness ≈ / of the momen-tary radius), the conversion of the rate of kinetic energy into luminosity and the observer’spoint of view allows simulating the map of the theoretical intensity of GW 46.4+5.5, seeFigure 19. The grand design of the SBs in a spiral network is also simulated, see Figure 22. Acknowledgments
The credit for Figure 2 is given to JPL/Caltech. At the moment of writing some animationsconcerning the spiral galaxies as modeled by percolation theory are available at http://personalpages.to.infn.it/˜zaninett/animations.html . References [1] Heiles, C., H I shells and supershells, ApJ 229 (1979) 533.[2] Puche, D., Westpfahl, D., Brinks, E., & Roy, J.-R., Holmberg II—A laboratory forstudying the violent interstellar medium, AJ 103 (1992) 1841.[3] Walter, F., Kerp, J., Duric, N., Brinks, E., & Klein, U., X-Ray Emission from anExpanding Supergiant Shell in IC 2574, ApJ 502 (1998) L143.4 Zaninetti[4] Pikel’Ner, S. B., Interaction of Stellar Wind with Diffuse Nebulae, Astrophys. Lett. 2(1968) 97.[5] Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R., Interstellar bubbles.II—Structure and evolution, ApJ 218 (1977) 377.[6] Tenorio-Tagle, G. & Bodenheimer, P., Large-scale expanding superstructures in galax-ies, ARA&A 26 (1988) 145.[7] Santill´an, A., Franco, J., Martos, M., & Kim, J., The Collisions of High-VelocityClouds with a Magnetized Gaseous Galactic Disk, ApJ 515 (1999) 657.[8] Koo, B.-C., Heiles, C., & Reach, W. T., Galactic worms. I—Catalog of worm candi-dates, ApJ 390 (1992) 108.[9] Ambrocio-Cruz, P., Le Coarer, E., Rosado, M., et al., The kinematical properties ofsuperbubbles and H II regions of the Large Magellanic Cloud derived from the 3D H α Survey, MNRAS 457 (2016) 2048.[10] S´anchez-Cruces, M., Rosado, M., Rodr´ıguez-Gonz´alez, A., & Reyes-Iturbide, J.,Kinematics of Superbubbles and Supershells in the Irregular Galaxy, NGC 1569, ApJ799 (2015) 231.[11] McCray, R. A., Coronal interstellar gas and supernova remnants, in: A. Dalgarno &D. Layzer (Ed.),
Spectroscopy of Astrophysical Plasmas , Cambridge University Press,Cambridge, 1987, pp. 255–278.[12] McCray, R. & Kafatos, M., Supershells and propagating star formation, ApJ 317(1987) 190.[13] Mac Low, M.-M. & McCray, R., Superbubbles in disk galaxies, ApJ 324 (1988) 776.[14] Igumenshchev, I. V., Shustov, B. M., & Tutukov, A. V., Dynamics of supershells—Blow-out, A&A 234 (1990) 396.[15] Basu, S., Johnstone, D., & Martin, P. G., Dynamical Evolution and Ionization Struc-ture of an Expanding Superbubble: Application to W4, ApJ 516 (1999) 843.[16] Bisnovatyi-Kogan, G. S. & Silich, S. A., Shock-wave propagation in the nonuniforminterstellar medium, Reviews of Modern Physics 67 (1995) 661.[17] Begelman, M. C. & Li, Z.-Y., An axisymmetric magnetohydrodynamic model for theCrab pulsar wind bubble, ApJ 397 (1992) 187.[18] Moreno, E., Alfaro, E. J., & Franco, J., The Kinematics of Stars Emerging from Ex-panding Shells: An Analysis of the Gould Belt, ApJ 522 (1999) 276.[19] Mac Low, M.-M., McCray, R., & Norman, M. L., Superbubble blowout dynamics,ApJ 337 (1989) 141.ercolation and Superbubbles 35[20] Ferriere, K. M., Mac Low, M.-M., & Zweibel, E. G., Expansion of a superbubble in auniform magnetic field, ApJ 375 (1991) 239.[21] Tomisaka, K., The evolution of a magnetized superbubble, PASJ 44 (1992) 177.[22] Tomisaka, K., Superbubbles in magnetized interstellar media: Blowout or confine-ment?, MNRAS 298 (1998) 797.[23] Kamaya, H., Final Size of a Magnetized Superbubble, ApJ 493 (1998) L95.[24] Seiden, P. E. & Gerola, H., Properties of spiral galaxies from a stochastic star forma-tion model, ApJ 233 (1979) 56.[25] Seiden, P. E., The role of the gas in propagating star formation, ApJ 266 (1983) 555.[26] Schulman, L. S. & Seiden, P. E., Percolation and galaxies, Science 233 (1986) 425.[27] Zaninetti, L., Percolation and synchrotron emission. I—The case of spiral galaxies,A&A 190 (1988) 17.[28] Seiden, P. E. & Schulman, L. S., Percolation model of galactic structure, Advances inPhysics 39 (1990) 1.[29] Jungwiert, B. & Palous, J., Stochastic self-propagating star formation with anisotropicprobability distribution, A&A 287 (1994) 55.[30] Palous, J., Tenorio-Tagle, G., & Franco, J., Star formation in differentially rotatinggalactic disks: The physics of self-propagation, MNRAS 270.[31] Zaninetti, L., On the Shape of Superbubbles Evolving in the Galactic Plane, PASJ 56(2004) 1067.[32] Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P.,
Numerical Recipesin FORTRAN. The Art of Scientific Computing , Cambridge University Press, Cam-bridge, 1992.[33] Bisnovatyi-Kogan, G. S. & Silich, S. A., Shock-wave propagation in the nonuniforminterstellar medium, Rev. Mod. Phys. 67 (1995) 661.[34] Dickey, J. M. & Lockman, F. J., H I in the Galaxy, ARA&A 28 (1990) 215.[35] Lockman, F. J., The H I halo in the inner galaxy, ApJ 283 (1984) 90.[36] McKee, C. F., Astrophysical shocks in diffuse gas, in: Dalgarno, A. & Layzer, D.(Eds.),
Spectroscopy of Astrophysical Plasmas , 1987, pp. 226–254.[37] Zaninetti, L., Evolution of superbubbles in a self-gravitating disc, Monthly Notices ofthe Royal Astronomical Society 425 (2012) 2343.[38] Dyson, J. E. and Williams, D. A.,
The Physics of the Interstellar Medium , Institute ofPhysics Publishing, Bristol, UK, 1997.6 Zaninetti[39] Padmanabhan, P.,
Theoretical Astrophysics. Vol. II: Stars and Stellar Systems , Cam-bridge University Press, Cambridge, 2001.[40] Smith , N., Dissecting the Homunculus nebula around Eta Carinae with spatially re-solved near-infrared spectroscopy, MNRAS 337 (2002) 1252.[41] Spitzer, Jr., L., The Dynamics of the Interstellar Medium. III. Galactic Distribution,ApJ 95 (1942) 329.[42] Rohlfs, K. (Ed.),
Lectures on Density Wave Theory , Vol. 69 of Lecture Notes inPhysics, Springer-Verlag, Berlin, 1977.[43] Bertin, G., Dynamics of Galaxies, Cambridge University Press, Cambridge, 2000.[44] Padmanabhan, P.,
Theoretical Astrophysics. Vol. III: Galaxies and Cosmology , Cam-bridge University Press, Cambridge, 2002.[45] Hill, C. J., Ueber die Integration logarithmisch-rationaler Differentiale, J. ReineAngew. Math. 3 (1828) 101.[46] Lewin, L.,
Polylogarithms and Associated Functions , North-Holland, New York,1981.[47] Olver, F. W. J., Lozier, D. W., Boisvert, R. F., & Clark, C. W., NIST Handbook ofMathematical Functions, Cambridge University Press, Cambridge, 2010.[48] Kompaneyets, A. S., A Point Explosion in an Inhomogeneous Atmosphere, SovietPhys. Dokl. 5 (1960) 46.[49] Olano, C. A., The propagation of the shock wave from a strong explosion in a plane-parallel stratified medium: The Kompaneets approximation, A&A 506 (2009) 1215.[50] Tomisaka, K. & Ikeuchi, S., Evolution of superbubble driven by sequential supernovaexplosions in a plane-stratified gas distribution, PASJ 38 (1986) 697.[51] Stone, J. M. & Norman, M. L., ZEUS-2D: A radiation magnetohydrodynamics codefor astrophysical flows in two space dimensions. I—The hydrodynamic algorithmsand tests., ApJS 80 (1992) 753.[52] Wouterloot, J. G. A., Brand, J., Burton, W. B., & Kwee, K. K., IRAS sources beyondthe solar circle. II—Distribution in the Galactic warp, A&A 230 (1990) 21.[53] Kim, K.-T. & Koo, B.-C., An Infrared and Radio Study of the Galactic Worm GW46.4+5.5, ApJ 529 (2000) 229.[54] Igumenshchev, I. V., Tutukov, A. V., & Shustov, B. M., Shapes of Supernova Rem-nants, Soviet Astronomy36 (1992) 241.[55] Koo, D. C. & et al., Deep Pencil-Beam Redshift Surveys as Probes of Large ScaleStructures, in:
ASP Conf. Ser. 51: Observational Cosmology , 1993, 112–+.ercolation and Superbubbles 37[56] Hartmann, D. & Burton, W. B.,
Atlas of Galactic Neutral Hydrogen , Cambridge Uni-versity Press, Cambridge, 1997.[57] Perrot, C. A. & Grenier, I. A., 3D dynamical evolution of the interstellar gas in theGould Belt, A&A 404 (2003) 519.[58] Zaninetti, L., Models of Diffusion of Galactic Cosmic Rays from Superbubbles, Inter-national Journal of Modern Physics A 22 (2007) 995.[59] Rybicki, G. & Lightman, A.,
Radiative Processes in Astrophysics , Wiley-Interscience,New York, 1991.[60] Hjellming, R. M.,
Radio Stars in Galactic and Extragalactic Radio Astronomy ,Springer-Verlag, New York, 1988.[61] Lozinskaya, T. A., Moiseev, A. V., Podorvanyuk, N. Y., & Burenkov, A. N., Syn-chrotron superbubble in the galaxy IC 10: The ionized gas structure, kinematics, andemission spectrum, Astronomy Letters 34 (2008) 217.[62] Maddox, L. A., Williams, R. M., Dunne, B. C., & Chu, Y.-H., Nonthermal X-rayEmission in the N11 Superbubble in the Large Magellanic Cloud, ApJ 699 (2009)911.[63] De Young, D. S.,
The Physics of Extragalactic Radio Sources , University of ChicagoPress, Chicago, 2002.[64] Zaninetti, L., Scaling for the intensity of radiation in spherical and aspherical planetarynebulae, MNRAS 395 (2009) 667.[65] Gray, M. D., Matsuura, M., & Zijlstra, A. A., Radiation transfer in the cavity and shellof a planetary nebula, MNRAS 422 (2012) 955.[66] Opsenica, S. & Arbutina, B., Numerical Code for Fitting Radial Emission Profile of aShell Supernova Remnant. Application, Serbian Astronomical Journal 183 (2011) 75.[67] Goldstein, H., Poole, C., & Safko, J.,