Astrophysical Blastwaves Citation Zukin, Phillip Gregory (2006) Astrophysical Blastwaves. Senior thesis (Major), California Institute of Technology. doi:10.7907/j6yy-wr15. https://resolver.caltech.edu/CaltechTHESIS:02122018-164440685 Abstract Shock waves are known to significantly alter the medium through with they propagate. Shocks in the interstellar medium have been known to affect star formation. Likewise, shocks in the intergalactic medium are expected to affect galaxy formation and the confinement of intergalactic clouds. Because of their significant influence, it is important to have a sound understanding of the different solutions used to study shocks and see how they differ. In this paper, I will focus on the self similar solution, thoroughly developed by Ostriker and McKee, and the explosion model, developed by Tegmark, Silk, and Evrard. Item Type: Thesis (Senior thesis (Major)) Subject Keywords: Physics ; Astrophysics Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Physics Thesis Availability: Public (worldwide access) Research Advisor(s): Kamionkowski, Marc P. Group: Astronomy Department Thesis Committee: None, None Defense Date: 2006 Record Number: CaltechTHESIS:02122018-164440685 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:02122018-164440685 DOI: 10.7907/j6yy-wr15 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 10693 Collection: CaltechTHESIS Deposited By: Joy Painter Deposited On: 13 Feb 2018 17:37 Last Modified: 14 Nov 2025 19:48 Thesis Files Preview PDF - Final Version See Usage Policy. 11MB Repository Staff Only: item control page

Senior Thesis: Astrophysical Blastwaves Phillip Zukin - Caltech 2006 Mentors: Dr. Marc Kamionkowski Dr. Steve Furlanetto 2 Table of Contents Introduction 4 Background 4 Self Similar Solution 5 Static Universe 5 Expanding Universe 9 Explosion Model 12 Comparing Models 14 Static Universe 15 Expanding Universe 18 Energy Analysis 21 First Phase of the Remnant 25 Conclusion 26 Appendix A: Accuracy of the numerical analysis 27 Appendix B: Post- Shock Temperature 29 Appendix C: Code 30 References 41 3 Introduction Shock waves are known to significantly alter the medium through with they propagate. Shocks in the interstellar medium have been known to affect star formation. Likewise, shocks in the intergalactic medium are expected to affect galaxy formation and the confinement of intergalactic clouds. Because of their significant influence, it is important to have a sound understanding of the different solutions used to study shocks and see how they differ. In this paper, I will focus on the self similar solution, thoroughly developed by Ostriker and McKee, and the explosion model, developed by Tegmark, Silk, and Evrard. Background A shock resulting from an explosion undergoes four phases. In the first phase, the external medium can be neglected, as the ejected mass undergoes free expansion. Then, when the shock absorbs a mass equal to the initially ejected mass, the shock begins its adiabatic phase, where it can be assumed that radiative losses are negligible, and all of the shock's energy is conserved. At this point, the well known Taylor Sedov solution is applicable and the shock is said to be self-similar. Later, when radiative losses come to total the amount of the initially injected energy, the shock enters the snow plow phase, where momentum is conserved. Last, as the pressure in the interior of the shock drops to the pressure of the external medium, the shock merges with the outside medium. This paper focuses mainly on the second phase of the supernova remnant. I will analyze the models in both a static and expanding universe. The static case applies when the shock remains in its galaxy. The expanding universe case applies when the shock reaches the IGM. In order for the shock to reach the IGM, it must either be infused with energy from multiple supernova bursts, with a total energy input estimated at around 1053 ergs, or the galaxy must be small, which was the case for early galaxies. In addition, I will use the models to analyze upper limits on how much energy can be inputted into the system so that the shock becomes weak on time scales smaller than the Hubble time. Last, I will give the solution to a simplified model that takes into account the first correction to a blastwave during its free expansion phase. 4 Self Similar Solution A self similar flow means that the flow at any point and time looks the same as it did at a prior point and a prior time. In other words, while the flow evolves in real space, it is invariant once suitably scaled. This implies that any non-dimensional quantity constructed out of parameters describing the system must remain constant for all time. Hence, every quantity must have power law dependences on a specific parameter, the most convenient of which is the normalized shock radius. Solution in a static universe: The parameters used to describe the system in this case are the initial injected energy (E), the outside density (p 1 ) , the radius of the shock (r), and time (t). Given these quantities, we can make a non-dimensional constant ~ as shown below, and then solve for the radius of the shock. 1.1.1 I Et 2 \ ~ r sh (t ) = .; 0 I - - I ~ Pl ) where .;o = nume rically computed c onstant 1.1.2 For this case, a single supernova is injecting all of the energy. However, if multiple supernovas are firing, the energy is additive, and when looking at the solution at time scales when the self similar solution is valid, it is as if a single supernova injected the sum total of all the multiple energy inputs. Note that ~ 0 is a specific value chosen to enforce energy conservation (see equation 1.1.12) Given this power law solution for the radius of the shock, it is important to get a feel for the magnitude of quantities we are dealing with and approximately when the solution is valid. Look to the table below for reference. Clearly, the velocity of the shock cannot exceed that of the initial ejected mass. For an initial energy input of 105 1 ergs and an external density of 2 * 10-24 g cm- 3 , assuming an ejection of 1 Mo, we will have an initial velocity of 104 km s- 1• So we cannot expect the solution to be valid for times much smaller than lOOyr. Now consider how much energy is radiated away. The cooling function has a value of 10-2 1 - 10-22 erg 5 cm-3 s- 1 for a temperature of 2 * 105 K. With these values, you'd expect by a time of 10 5 yr for about 1051 ergs to be radiated away, which implies that the solution will not be valid for times much larger than 105 yr. (For a derivation of the post shock temperature, shown in the table, look to appendix B. To see how the radiated energy was found, look to the Energy Analysis section) t (yr) rsh (pc) Ush (km S- 1) T2 (kelvin) 1 .315 124,000 2.08El l 10 .791 31,000 1.3El0 100 1.99 7,820 8.3E8 1,000 4.99 1,970 5.3E7 10,000 12.5 494 3.3E6 100,0000 31.5 124 2.1E5 Aside from the shock radius, it is also important to analyze the interior structure of the shock. In order to know how the pressure, velocity, and density in the interior of the shock act, we first need to know the values of these quantities directly after the shock. These quantities, labeled with a subscript 2, are called the post-shock values. They are determined by the Rankine - Hugoniot jump conditions, which come about by integrating the continuity equations over the shock. The post shock values, in terms of the external values, in the limit of a strong shock, are shown below. ( y + 1 I P2 = 1- - I Pl 2 --ush \ )' - 1 ) 2 P2 = · -- y+ l )'+ 1 Pl Ush 2 1.1.3 Once these quantities are known, we can begin to determine the internal structure, by using the fluid conservation equations shown below (mass, momentum, and energy). op ·- + at 1 a (r2 p u ) = O r2 ar au au at ar - - + u - 1.1.4 1 oP =- - p or 1.1.5 _!!__ [P (E + ~ u 2 ) ] + ]_ -~- [r2 p u rE + -~ + ~ u 2J] = 0 at 2 r2 a r \ p 2 ) 1.1.6 6 We first assume the hydrodynamic variables have the following dependences on c; : I .;- 12 P ( r, t ) = P2 I - I \ .;-oJ p ( ,;- ) 1.1.7 Plugging in these dependences into the fluid conservation equations, we are left with 3 coupled ordinary differential equations. - .;-~ + -~ [ 3av +.;-~- (av) ]= 0 2 dv 5 ct.;- -v - --- ,;- - - + - d i;- '( + 1 4 I 2 5 ('( + 1 ) \ 1.1.8 dv\ 2 l y- 1 \ 1 I dp \ ct.;- } 5 \ '( + 1 } O'. ct.;- } - - - I v + v,;- - 2 2 d 2 - 2 ( p + O'.V ) - - - ,;- --- (p + O'.V ) + 5 di;- I 1- - 1 12 p +.;- \ 1 1.1. 9 4 d [ V ( yp + O'.V2 ) ) } = 0 - - - { 5 V ( yp + O'.V 2 ) + ,;- d,;- di;- 5 (y + 1 ) 1.1.10 Note though, as mentioned before, that only one particular value of .;-o is appropriate. In order to determine this value, we have a normalizing condition, that forces energy conservation. l r sh( t ) [ 0 P 1 y -1 2 --- + - pu 2] 2 4 lfr di r =E 1.1.11 The first term represents the thermal energy in the shock while the second represents the kinetic energy. Note that radiative losses such as ionizing the external medium and Compton cooling off the CMB are not taken into account. This equation, once the above dependences are substituted in, yields the following integral. 32 7f 2 5 (y 2- 1) ---------- J(o 0 [ p (.;- ) + O'. (.;-) v 2 4 (.;-) l .;- di.;- = 1 1.1.12 Only for an appropriate value of c; 0 will the integral be satisfied. In our case the value is approximately 1.17. It is important to notice that self similarity forces the thermal and kinetic energy to always remain in the same ratio. Numerically integrating the above 7 equations ( 1.1.8 - 1.1.10), we get the following graphs for the internal structure of a blastwave in a static medium. p p -P,- P2 11 11 o.81 I I I 0 .8 J 0. 6 :::L 1 0.4 0. 2 0.2 0. 2 - -0.= = 4 ---~---~- -- ~~__._~ 0.2 0.4 Figure 1.1.1 0.6 0.8 I _ 1 Rs __,_ r Figure 1.1 .2 u U2 1 0 .8 0.6 0.4 0.2 ---·-"" --- -- · - -- ~~ --- ~~ 0.2 0. 4 0.6 0 .8 -- ·-· _____._ 1 r R, Figure 1.1.3 Figures: 1.1.1 normalized density profile, 1.1.2 normalized pressure profile, 1.1.3 normalized velocity profile Notice that the density profile is biased towards the shock radius. Since radiative cooling has an n2 dependence, it is clear that radiative effects will be most prevalent towards the shock radius. As a result the particles at a smaller normalized radius will catch up with the mass in front of them, and as time progresses, a shell structure will form. Also note that both the velocity and density go to zero at the shock center, while the pressure remains nonzero, signifying a very hot inner core. A plot of the temperature of the blastwave, as a function of normalized shock radius, is shown below. 8 T Tz 1. x 10 6 10000 100 1 Figure 1.1.4 - normalized internal temperature. Solution in an expanding universe In an expanding universe, the density is time dependent. As a result, the time dependence of the shock wave changes. We are interested in analyzing a matter dominated universe with Q = 1, where Q is the ratio of the energy density to the critical density today. Given this type of universe, the scaling of the density is shown below. 2 whe re h = -- f or Q = 1, r is the cosmic time 3 a nd Pu is the density of the external medium when the shock begins. PO = Pu 1: -3h 1.2.1 In addition, it is possible to make the solution more general, given a new formalism for expressing different parameters in terms of power law dependences on a normalized shock radius. The formalism is shown below. x = x ( 1 ) R-kx 1.2.2 Where x is any hydrodynamic variable, x(l) represents the post shock value of the variable, and R is a normalized shock radius. This new formalism allows for the possibility of introducing energy decay and other complications. The new time dependence of the shock wave, assuming an expanding universe and no energy decay, is expressed below. Rs E d = .;- II ----I 1: (2+3h) /5 \ Pu} 1.2.3 9 Since for Q = 1, h = 213, the exponent of time is 4/5. It is intuitive that the exponent for the shock should be greater in a moving medium. Primarily, since the density of the outside medium is decreasing with time, there is less frictional force stopping the shock from propagating. In addition, the shock picks up kinetic energy from the Hubble flow . Note that because of self similarity, the ratio of the shock velocity to the Hubble velocity \ must be constant for all time. This ratio, which depends on the cosmology, is shown Y(0'"~ 5 6 2+3 h 1.2.4 Now given these expressions for the shock radius, we can look for the internal structure of the blastwave. Once again, we look to the Rankine - Hugoniot jump conditions to establish a post shock value. Because of the velocity of the external medium, the post shock values have a different form. In the limit of a strong shock, they are shown below. ( y + 1 I P2 = 1 - - I Pl \ y - 1 } 2+(y - l )VH - -- - u sh y+l 1.2.5 Given our new formalism, it is more convenient to 1) use the entropy equation instead of the energy equation and 2) express the coupled differential equations in terms of logarithmic derivatives. The definitions used in the differential equations are shown below. x = x1 ( t ) x (.A) x* = a lnx a lnr r wh ere.A = - Rs 1.2.6 dl n x dl n.A 1.2.7 Note that x once again represents any hydrodynamic variable, and x 1 represents the post shock value. The differential equations, expressed in terms of the above definitions then become: 10 p * ( 1- v) +kp-2v vv * (1- v) v * = e P* v 1.2.8 2 1 Q Yh g + - (kp - k p) + -----2 y2 _A.2 2 1.2.9 ( 1- v) p * = y ( 1- v) p* +'(kp - kp p where e = --- v 1.2.10 v P VJ..2 Note that the coupled differential equations all depend on gravity. But the gravity term is dependent on the density distribution inside the shock. Hence, it is best to assume a certain density profile, compute the internal distributions, and then recompute the profiles based on a gravity term that coincides with the previously computed density distribution. I iterated this approach until I reached a steady state solution. Plots of the internal distributions for an expanding universe are shown below. p p \~ 0.975 !" 5 -~- : 0.98 0 . 985 0.99 0.995 1. 4 I 1.2 ~ 975 21 ~-~~ --~---~-~--- -~--~---~--l0 .97 / . r ------ -~~ P2 r Rs 0.98 Figure 1.2.1 0.985 0.99 0.8 0.6 I o.4 I Figure 1.2.2 u U2 ~-~----- ~ 1.01 '~\: I l~ I r -~-- -~---· · -- . -~---------~--\l- R, 0 . 97 0. 975 0. 98 0. 985 0 . 99 0 . 995 Figure 1.2.3 Figures: 1.2. l normalized density profile, 1.2.2 normalized pressure profile, 1.2.3 normalized velocity profile 11 Notice that all the material is concentrated between A= .965 and A= 1. In addition, notice that the pressure approaches zero at A= .965. If this were not the case, then the nonzero pressure would push the material in all directions, and the shell structure would not exist. Moreover, the steadily decreasing velocity profile implies that the velocity of objects at smaller A will catch up the mass in front of it, further contributing to the shell structure. Explosion Model Since it takes a while for the self similar solution to become valid, it seems unnatural to use the energy of the initial injected energy in the formula, because of losses during the free expansion phase. As a result, Tegmark, Evrard, and Silk developed a different solution that models a blastwave from the time of the explosion to the time the shock merges into the IGM. This model has 3 assumptions about the blastwave. 1) It assumes that most of the mass in the shock is concentrated in a thin spherical shell centered at the shock radius. 2) It assumes a neutral, uniform, pressure-less IGM. 3) It assumes that the internal plasma is thin and has a constant pressure and temperature. Given these assumptions, they derive a force formula which characterizes the system. 2.1 where p is pressure, H is the Hubble constant, Qb is the density of baryonic matter today divided by the critical density, and Qd is the ratio of dark matter today divided by the critical density. The first term in the equation represents the pressure of the internal plasma driving the shock outward. The second term is a braking force, which describes the necessary force needed to increase the velocity of the outside material from the Hubble velocity to the shock velocity, as it is picked up. The last term is the gravitational force. Note that the gravitational term resulting from the baryonic matter is half that of the dark matter. The Yi factor comes from the self-gravity of a thin shell. Also, note that the dark matter does not interact with the shock. This model assumes an Q = I matter dominated universe, only composed of baryonic and dark matter. This is a valid 12 assumption for redshifts smaller than 3000, when the matter and radiation energy density were equal , till recent redshifts, when the cosmological constant has taken over. In addition to an equation for the radius, they also use an equation that describes how the internal thermal energy evolves with time. 3 Et= - 2 dEt dt Lsn pV 211pR 3 dV 2.2 2 dR L - p -- = L - 4 11pR --- dt dt 2.3 L = Lsn - Lcarp - Ltrems - Lion + Lmss 2.4 Mt 1. 2 lo -- e (t bum - t ) !vb 3rn Lmss = fd dR 2R ( dt - HR) 3 2.5 where L is the sum of all radiative losses. Since the internal plasma is assumed to have constant pressure, the internal energy is straight forward to compute. A change in internal energy occurs because of radiative losses and the work necessary to expand the shock. The only terms we focused on for the radiative effects were Lsn and Laiss· Lsn represents a sourcing term that occurs while multiple supernovas are going off. Notice that is has the form of a constant energy injection into the shock for times smaller than tbum (~ 5 * 10 7 years). Laiss occurs because we assume the collision of the shock with the external medium is inelastic. As a result, some of the excess energy, after interacting with the shell, leaks inside to reheat the interior plasma. However, the amount of energy that does this is tough to gauge so fa is introduced to account for the complicated microphysics. The parameter fa is very important, and we will vary its value in our analysis. These two equations give us two coupled differential equations. Moreover, normalizing the variables in the context of a matter dominated Q = 1 universe, we get the following equations, which can be numerically integrated. r ,, 18 11 (r) = - - T/ ( r ) Qb [ 2 - - · - 3 1- r(r) 3 _2 q ( r) T/ ( r) r ( r) ] 2 r ' ( r) r ' (r) 2 r ( r) 13 - ( 92 12ct + 91 Qb) 77 ( i:) 2 r ( i:) l (I:) r I (I:) q I ( "[; ) = - - - - - 5 2 JTr ( i: ) 3 r ( i:) 2.6 q ( "[;) 2.7 t wh e re i: = - - 1 and t* t* = hubble time . q, 1, and r are all the normalized versions of p, L, and R respectively. rt is the ratio of the Hubble constant at -c to the Hubble constant when the shock begins expanding. Note that the Hubble time t. is a function of redshift, and is the time at which the shock begins expanding. Just to get a feel for the order of magnitudes we are dealing with, below shows a plot of the shock radius as a function of time. The upper line is for fd = 1, while the lower is for fd = 0. 0.5 0. 4 0.3 0. 2 0.1 5 10 15 20 25 30 35 Figure 2.1 - radius of the shock as a function of normalized time. Comparing the Models We have chosen several possible ways to compare the two models. First is the time dependence of the solution. Given the clear time dependence in the self similar solution, we can look at how the explosion model scales with time at late times in the solution. We have also decided to test the accuracy of the explosion model's assumptions by comparing its internal structure to the internal structure of the self similar solution. Then, given these distributions, we want to compare the energy distributions of both solutions. 14 Comparisons Static Universe It would be best to first compare the models, in a simpler context, by neglecting the Hubble flow and gravity. In this case, the explosion model's equations of motion reduce dramatically and are shown below. r " (r) 4 7r -~-~- - 3 _:_~_l:_)_~ r (r) 1 (I) r (r) 3.1.1 r 1 (I ) q' ( I ) = - - -- ----- - 5 - - - - q ( I ) 2 7rr ( r ) 3 r ( r) 3.1.2 Note that the explosion model, even neglecting radiative losses, does not conserve energy because of the losses due to the work in expanding the shock. Since the self similar solution conserves energy, it is most natural to first force energy conservation in the explosion model, and hope the time dependence matches at late times. The new equation for the loss in thermal energy is shown below. E = ~ pV - qr 3 2 cans t _, q' r 3 + 3 qr 2 r ' -? q' r' -3- r 0 3.1.3 q 3.1.4 With these new equations, numerically we found that the solution scales as 2/5, which matches exactly with the self similar solution. Next we looked at the scaling of the explosion model using its original equations of motion. This time however, we were curious to find for what value of fct the solution looked like the self similar solution. Below is a chart of the exponent of time for various values of fct at late times in the solution. Note that for this simplified case, there is no appropriate time to check the time dependence of the solution, but it was clear that the time dependence of all solutions were asymptotic. As a result, we checked the time exponent at late times. 15 fct Ex of time 1 .4004 .8 .6 .3812 .3616 .4 .3410 .2 .3184 0 .2901 We can see from the chart that fct = 1 best approximates the self similar solution. Analytically, I was also able to verify the exponents for ~ = 0 and fct = 1. This was a good check that the program was working appropriately. The analysis is shown below. r I fct = 0 --7 q I = - 5 - - q --7 q- r -5 r --7 r' ' - t n-2 _ S - r-6 - t - 6n --7 n = 2 7 r r I fct = 1 --7 q I -3 - - q--? q-r 3.1.5 -3 r --?r 1 1 - t n-2 - -q - r - 4 - t -4 --? D= r 2 5 3.1.6 It was surprising that, in the case of fct = 1, a dissipative system acts like a system that conserves energy. In order to get a better understanding of how the system is behaving, we looked at the thermal and kinetic energy distributions as a function of time. They are shown below. K Th E i nit 1 0.5 I I E1r.iL 1 o.s IL~- -- - --------- 1 0.2 0 .2 0.1 0 .1 0 .05 0.05 0.02 0.02 .~~- "- · -·---~ 0 200 400 600 800 1000 Figure 3.1.l - Kinetic energy vs time. I I L_._.~· ----~·---~----- ____. __ ~ 0 200 400 600 800 I 1000 Figure 3.1.2 - Thermal energy vs time. Notice that after the first bump, the system conserves its energy. Looking at equation (3.1.6), which is valid at late times, this should not be considered a coincidence, since for this case q ~ ( 3 • Hence, the total thermal energy, which goes as qr3 , should be constant. 16 Another way of understanding the conservation of energy is that the only energy lost by the system, which results from the inelastic collisions, is regained by the system as thermal energy. In addition, note that since the self similar solution conserves energy, the time it takes for this particular case to stabilize may be a good indication for how long it takes an explosion to become self similar. Next in the analysis is verifying the assumptions the explosion model made about the internal structure of the shock. Look to figures 1.1.1 - 1.1.3 , for reference of the internal structure of the shock in the self similar solution. A plot of the explosion model's internal structure is not shown because it is not instructive. It consists of all the mass concentrated in a thin shell at the shock radius, with uniform pressure inside. Note that the self similar solution does not include certain microphysics taken into account by the explosion model. In the explosion model, it is assumed that a small amount of the ionized medium leaks from the shell into the interior, providing the internal pressure that drives the shock. Moreover, they do not consider the thermal energy of the shell itself. The self similar solution, on the other hand, does not take into account the ionized medium produ67by the shock, but does compute the pressure associated with the mass inside the shock. Hence, the only characteristic we can compare between the two internal structures is the density distribution. Looking at the density profile, for the self similar case, we see that the mass is concentrated towards the shock radius. But, it is clear that all of the mass is not concentrated in a thin spherical shell. Moreover, assuming the pressure for the self similar case roughly corresponds to the pressure of the ionized medium, we can see that it is approximately constant for r/R smaller than .6. Hence@ ems to be poor agreement between the internal structure of the Explosion Model and the internal structure of the self similar solution in the case of a shock in a static medium. Aside from comparing internal structures, we can also compare energy distributions. We could ask again for what value of fd will the energy distributions in the explosion model be equal to that of the self similar solution. Note that for a self similar solution, the shock contains 70% thermal energy and 30% kinetic energy. The table below shows the 17 computation done for several values of fd, where the first pair of energies are scaled by the initial energy inputted into the shock, and the last pair is scaled by the final energy . fd T I Einit K / Einit T I Etin K / Efin 1 .6 .4 .6 .4 .8 .26 .19 .577 .423 .6 .110 .088 .556 .444 .4 .0433 .0407 .515 .485 .2 .015 .0176 .460 .540 0 .0035 .0067 .343 .657 We can see that no value of fd will match the energy distributions in the self similar solution. Moreover, the energy in the system is only conserved for fd = 1. No matter how high fd is, the model will never be able to convert enough collisional energy into thermal energy that reheats the interior of the shock to the self similar value. This implies that there may be another mechanism for transferring heat into the interior that the explosion model has not taken into account. Or it's possible, as we saw before, that the thin shell approximation is not valid. It may be the case that the majority of the mass is not moving at the shock velocity, but some velocity slightly lower. As a result, the kinetic energy would be smaller and the energy distributions could match the self similar case. Comparisons Expanding Universe Next, we moved on to a more general problem, taking into account the expansion of the um verse. As before, we first compare time scalings for the modified explosion model, which conserves energy (see equation 3.1.4). In this case, just as before, we get a match of .8. Next, we analyzed how the time scales for the original explosion model. Below is a table showing the computed time dependences for various values of fd. Note though that for the case of the expanding universe, there is a specific time, where the solution'-s time exponent can be tested. Since we know, from before, that for a cosmological ~ave, the ratio of the shock velocity to the Hubble velocity is a specified value, we can check the exponent of time at a point in the evolution of the blastwave where its velocity is 6/5 18 of the Hubble velocity. Since it's expected that all solutions to the explosion model will eventually match this criteria, given that the velocity of the shock begins at a velocity higher than this value and decays to the Hubble velocity monotonically, I included the time at which this criteria is met in the table below. fct Exp oft Time 1 .7999 7.87 .8 .7999 7.27 .4 .7999 6.18 .6 .7999 6.71 .2 .7999 5.67 0 .7999 5.18 First we can see that the exponent of time for all values of fct matches the self similar solution exactly. This however, should be expected, because we are measuring this value at a time when the explosion model is approximately self similar. Note that the solution does not stay at this exponent, and hence does not remain self similar, for all times after this period. Plots of the scaling of the shock radius as a function of time are shown below, for multiple values of fct. We can see that the solutions begin to asymptote to the same value at late times, which is expected as the shock begins to merge with the IGM. Also, we can gauge how long the solution is approximately self similar by measuring when the scaling is between some E of .8. 0.6 0.4 0.2 ---~~--~.~~---~-~ 2 4 6 8 10 Figure 3 .2.1 - dln(R) I din T: (the exponent of time). Line above is for fd = 1. Line below is for fd = 0. What's more interesting is the time it takes for the explosion model to become approximately self similar, as a function of fct. We can see, as we expect, that it takes longer for the solution to become self similar for higher values of fct because, there is more internal pressure driving the shock at speeds too high to be self similar. Moreover, 19 the relationship between this time and fd seems linear, as the following plot shows. This is surprising, in that the physics of is fd complicated, as discussed before, and you would not expect it to have a simple linear effect. 7.5 7 6. 5 - - ·- - - -- - - - - - . - -·---·--· ·- - --· ·- 0 .2 0.4 0.6 0 .8 fct 1 5.5 Figure 3.2.2 - t: vs fd. Next, we compared the internal structure of the blastwave computed in the self similar solution to the assumptions of the explosion model. Look to figures 1.2.1 - 1.2.3 for the computed internal structure of the self similar solution. As discusses before, both models take into account thermal pressure from different sources. The explosion model has its thermal energy sourced by the ionized plasma, while the self similar solution's thermal energy comes from all material inside the shock. Because of this discrepancy, it is tough to gauge the difference in the internal structures of both models. However, we can compare density distributions, and we see that the shell approximation is valid in an expanding universe. Notice however that the shell is not uniformly distributed, and does not travel at the shock velocity, as assumed in the explosion model. Now we look to the energy distributions of each model. Below is a table listing the energy distributions for different values of fd. The first set is normalized by the initial energy input, while the second is normalized by the energy at the time of measurement. Note that once again I chose to compare the distributions when the shock velocity was 615 of the Hubble velocity. Also note that a self similar blastwave has 68.5% kinetic energy and 31.5 % thermal energy. 20 fd K / Einit T I Einit K / Efin T I Elin 1 .9329 .0338 .9650 .0350 .8 .7753 .0257 .9679 .0321 .6 .6474 .0192 .9712 .0282 .4 .5430 .0140 .9749 .0251 .2 .4571 .0097 .9792 .0208 0 .3861 .0062 .9842 .0158 We can see that no matter what value of fd, the energy distributions do no match the self similar case. Note that in an expanding medium, for the self similar solution, the kinetic energy dominates while in a static medium, the thermal energy dominates. This is to be expected because the shock absorbs the kinetic energy of the Hubble flow. Also, it's expected that the energy distributions above will not match that of the self similar solution because, just as in the static case, the model cannot account for enough thermal energy. Energy Analysis Next we wanted to find out on what time scales the shock becomes weak. A shock will be considered weak when it merges with an already ionized IGM, which is predicted to have a temperature of about 1OOOOK. Setting the mean thermal velocity equal to that of a gas of protons at this temperature, we get an approximate velocity for when the shock becomes weak. Vf = 3i1:6 cm [ -m- - 1. 9 * 10 - s 4.1 Setting this velocity equal to the velocity of a self similar shock in a static medium, we get a relation for the time at which the shock becomes negligible. Vs = 4.2 21 Note that this time is a function of redshift because the external density is a function of redshift. In this approximation though, we are not taking into account the time dependence of the external density, as the shock propagates. We are assuming that the external density at the time of the explosion is approximately correct throughout the shock's propagation on the time scales we are interested in. Adjusting for the redshift dependence, we get the following. t~ (3- _£~ ( ~-- I ~ ) \ 5 Vf \ Qb Per ) ) 5 3 1 - - i. 34 * 10 16 s (__!__I 1 + z* \ E 51 ) ~ 1 1 + z* 4.3 where z* = redshift when shock began e xpanding Per= critical density today Qb = ratio of baryonic density to critical density Next, we want to compare this time to the hubble time when the shock began. This can be easily solved for, assuming a matter dominated universe with Q = 1. tH= 2 3 Ho I 1 \ ~ - 2.90*10 n SlI ----1 d 1 1- - 1 \ 1 + z* ) \ 1 + z* ) 4.4 where Ho = hubble constant today Setting these two expressions equal to each other, we get an upper limit on how much energy can@ itially injected such that the hubble flow can be neglected. 3 E < 10 55 I 1 I2 ergs I - - - I \ 1 + z* ) Note that most supernova explosions are on the order of 10 4.5 51 ergs. Only when multiple supernovas are firing at the same time, as is assumed in the explosion model, will the injected energy reach a value of 1053 ergs. Therefore we can see from the above limit that for redshifts of about 15, when shocks become important for their influence in ionizing the IGM, we know that the shock will become weak before a Hubble time elapses, hence the Hubble flow can be neglected, as is shown in the plot below. 22 t ( sec ) I i. x 101 7 5.x10 16 1. x 10 16 5. x 10 1 5 I 1. x lO " I 5. x 10 1 4 l_______ ~ ---~~----·--~ -~1 5 2 -~--- 20 10 · -·~·- 50 z 100 Figure 4. 1: Plot of tH and tweak as a function of redshift for E = 10 53 ergs. (tH is plot above for low z). We can see that for low redshifts, tH is large compared to tweak· For high redshift past the crossing point, the Hubble flow cannot be neglected, while for lower redshifts before the crossing point, it can. Given the upper limit, it is important to verify the limit's validity. First we can ask whether the solution is approximately self similar on these time scales. To address these issues, we need to determine approximately when the solution becomes self similar, and when, because of radiative losses, the solution is not self similar. The solution becomes self similar on time scales when the shock absorbs an amount equal to the initially ejected mass. 4 3 -- lTPo Rs 3 4 r~o (E d I I 3 \ Mo ~ -- lTP o 1) t 5 \ PO } I 3 Mo ) 5 ~ t 16 11 1 i~( M li I -- -- __2___ -- -3 1 - 2 . 83 * 10 s 11:-;. JI I - I ( 3Mo \4 JTp 0 5,o 3 E5 } ' \Mo} 4.6 Given this computation, and looking to the figure above, we can see that the time scales we are interested in are well past the first phase of a SNR. Now we look at radiative losses. The energy dissipated goes as: 2 E - Ln Vt 4.7 23 8 Qb P c r where n - - - - - ( 1 + z* ) 3 m 4.8 where L is the cooling function, n is the number density, V is the volume of the shock, and tis time. The expression for n is the external baryonic density times a factor of 8. The extra factor comes from the Rankine Hugoniot jump conditions which states that, for a gas with y-= 513, the post shock density is 8 times the external density. Moreover, since most of the cooling occurs near the shock, it is appropriate to use this value. At a temperature of 1OOOOK, the cooling function is at approximately 10-24 erg cm3 s- 1. Hence we can show that radiative losses become severe at a time shown below. L *n 2 ( E \ 51 4 * t coo1 * - 7f 3 --) t coo1 - 7. 87 * 10 \ 4.3 , we see tweak ~ 6.4 * 10 \ PO } 22 ( E For a redshift of 20, with an energy of 10 14 ( £ 0 I ---- I I - 2 \ 11 I l t 5 I 3 - E J - 63 ( 1 + Z*) l l \ Es1 } 51 ergs, tcool ~ 2.1 * 10 4.9 15 s. Looking to equation s. Therefore, for times that we are interested in, the shock remains self similar, and the above limit is valid. We can also form a lower limit on the energy necessary to produce a shock wave, as a function of redshift. In order to produce a shock, the ejected material must have a speed greater than the sound speed. Note that the sound speed is dependent on the temperature of the medium, which approximately scales with redshift as (l+z). Vs= [~kT- 1 ~ --) E - - m 2 2 41 MejectV --) E > 7. 5 * 10 ( l + z *) ergs M -- M::J 4.10 This limit is easily met by most explosions at a redshift of 20. 24 First Phase of the Remnant I was also hoping to have a solution for the blastwave during its initial free expansion phase. However, instead of having a simple free expansion, I wanted to take into account a first order correction term. It seemed logical to assume a frictional term dependent on the cross section of the particles of the outside medium, as well as the density of the external medium. By dimensional analysis, we get the following force equation. mx'' = - pav2 5.1 Note however that the first phase of the supernova remnant occurs on time scales small relative to the expanding universe. As a result, I will neglect the time dependence of the external density. Integrating the above equation, we get: a pv0 v (t ) = vo ( 1 + - - ; - t )-1 5.2 The coefficient in front oft is particularly interesting: p I m is the number density of the outside medium, which implies m / ap is the mean free path length, which makes m l a pv 0 the time it takes to travel the mean free path length. The above equation can be integrated again to get the radius of the shock. apv o ) x ( t ) = -m ln ( 1 + - t m ap 5.3 Notice that in the limit of apv 0t/m << 1, the solution reduces to free expansion, as you' d expect. Previously I mentioned that the free expansion phase breaks down when the shock accumulates a mass equal to the ejected mass . I'd like to verify this claim based on the above model. 4 - 3 npx 3 (3~\ ~ = ~ --? x = I - - - I --7 t * = \ 4 np ) m ap I 3~ \ 31 I {exp ( -II - 1} \ m \ 4 np ) ) 5.4 where M ab is the absorbed mass at t.. Now, plugging in this time to the velocity, we get: 25 1 op ( 3 Mab \ 3 } I --- 1 m \ 4 JTP ) v ( t *) = vo e xp { - -- 5.5 Clearly, the free expansion assumption is not valid anymore when the exponent is of order 1, at which time the velocity has decreased by more than a factor of 2. Therefore we have: 5.6 where l is the mean free path length previously discussed. Note that in a mean free path length, the shock on average has no time to interact, let alone, absorb any mass. Hence, the quantity we see above is approximately all of the initially ejected mass. The above model therefore, correctly predicts when the solution transitions out of the free expansion phase. Conclusion The above analysis compared the self similar solution to the explosion model in both a static and expanding medium. Doing this allowed us to understand the benefits of both models. While the self similar solution is approximately exact during the second phase of a supernova remnant, it does not take into account more complicated processes, like ionizing the IGM. The explosion model, on the other hand, which accounts for different cooling effects, is useful in that it models a blastwave throughout its entire evolution. However, it is not able to accurately account for the energy distributions during the second phase of the remnant. The energy analysis was useful because it revealed that for most supernova at a redshift of 15, the time at which shocks are important for ionizing the IGM, the Hubble flow can be neglected, which greatly simplifies the problem. Last, using a simplistic model that describes the first phase of a remnant, we were able to verify the claims made that the free expansion solution breaks down when the shock absorbs a mass equal to the initially ejected mass. The above analysis is fundamental because explosions, and their resulting shocks, greatly alter the surrounding medium. Since the surrounding medium is important to most astrophysical questions, knowledge of shocks is necessary. 26 Appendix A: Accuracy of Numerical Analysis The coding was all straightforward for most of the cases described above. I typically used a second order Runge-Kutta algorithm. However, numerically finding the internal structure of the cosmological blastwave was more difficult because there were singularities involved. As a result, to show the accuracy of my code, I will analyze the solutions of the coupled differential equations at the singularity, using the Ostriker and McKee formalism, and plot the difference between my numerical results and the actual solutions near the singularity. For a hollow blastwave, which was the case for the cosmological solution, the internal radius of the shock (behind which no mass exists), is defined as Ai. Since the structure of the blastwave is self similar, Ai remains constant, which implies the velocity at Ai is AiVs. Vi =1 A.I Looking at the entropy equation (1.2.10), we see that at .Ai, p* and P* must diverge. Ostricker and McKee, in order to solve this problem, then define the following : p+ = ( 1- v) p* A.2 And noting that as .A -> Ai e = ()" - 1 ) ( 1 - y) v 2 2 (yv -1 ) A.3 vv* = p + + kp - 2 v A.4 The continuity equations then reduce to : ()" - 1) y p + 1 2 (yv - l) 2 ( 1- v) v* = -- ·· --·-- + --- ( kp - k p) A.5 A.6 Solving these equations at .A = Ai , we get: A.7 27 Pi + = '( A.8 6 + 2 kE - kp - 2 '( Vi * = '( A.9 Assuming the following solution for a hydrodynamic variable x: 1 . x - (/\ - /\ i ) xi A.10 We get the following relation for the exponent: Xi + Vi * A.11 which implies: l pi = 2 ( 3 + kE ) - (y + 1 ) kp kp + 3 '( - 2 ( 3 + kE) A.12 A.13 Given these power laws, the following plots show how close my numerically computed solutions are to the actual solutions near the singularity. p p P2 P2 ~~1 .................. ..... .• ..... / 0.6 1' 0.5 . o.4 I 0.31' ..... ..•* 0 .2 .t* 0 .1 I -·---~· -·----~-·- ~-· -- ....__ ,_ - ----·---! 0 . 9685 0. 9686 Figure A. l - Actual solution is to the left 0 . 9687 0. 9688 0.9689 r Rs Figure A.2 - Actual Solution goes to zero 28 u 1.01056 1. 01055 L,- 1.01054 .........T 1. 01053 · -·-L~. ~--·-~----=~:I 0 . 9685 0. 9686 0. 9687 0. 9688 0 . 9689 r R;. Figure A.3 - Actual Solution is above For figure A.1, we see that the actual solution has a stronger singularity than my numerically computed solution. However, both profiles go to infinity at approximately the same normalized radius. For figure A.2, the actual solution reaches zero identically, while mine stops at .4. I believe this is just a result of the numerical integrator not being able to handle the coupled singularity in the density. For figure A.3, the difference in the actual verse numerical solution always stays within .0000 I, which is a trivial difference. In addition, to make sure difference would not affect the data, I computed the energy distributions for two sets of internal structures. The first set was my numerical solutions. The second set was my numerical solutions connected to the correct analytical solutions near the singularity. The difference in energy distributions between these two sets was 1 part in a 100. Appendix B: Post-Shock Temperature The specific internal energy for a gas particle is: E = E m p - -- - where y = ratio of specific heats (Y - 1 ) p B.1 Setting E equal to the energy in terms of T, we get: Pm (Y - 1 ) P 3 -- kT 2 ~ T = P 2m (y - 1) p 3k B.2 Then, plugging in the values for the post shock values of P and p given in equation 1.1.3, we get: 29 T 4rn 3k y-1 - (y+l)2 Ush 2 3rn = - - U sh 2 16k B.3 Where in the last step I took y = 513. Note that using the post shock temperature to determine the cooling function's value, even though the temperature increases dramatically inside the shock, is valid because most of the cooling occurs near the shock radius, where the density is highest. Appendix C: Code Code used for Self Similar Solution II Omega= 1 solution for expanding universe II including gravity #include #include #include #include <ios tream> <vector> <cmath> <fstream> using namespace std; vector <doub le > den; vector <double> vel; vector <double > pres; vector <double> gr av; double den 1, vel 1, pres 1; II important parameters (fix) double g 6.672E-8; doublet 365E2 * 24 * 60 * 60; II seconds //double t 10.; doublers= 1.898 * pow(1E51 * g * pow(t, 4) , 1 /5 .); //cm //doublers= pow(1E51 I 2E-24 * pow(t, 2), 1/5.) * 1.17; double out den= 1 I 6. IM PI I g I pow(t, 2); II gm I cm3 //double out den= 2E-30; double v_s; double step; int accuracy; II # of steps away from singularity to recompute double lambda 0 = l.; double gam = 5 I 3.; double eta= 4 I 5.; int length; int temp = 0; 30 double double double double k- p 3.; 2. 5; k d v h 5 I 6.; 1.; omega double double double double lam i; l_p = gam * (3. - k_d) I (k_d + 3. * gam - 6.); l_d = (6. - (gam + 1) * k_d) I (k_d + 3. * gam - 6.); c_p, c d, c v; void Init funcs() length = (int) (lambda 0 I step) + 1; den.resize(length); vel.resize(length); pres.resize(length); grav.resize(length); lambda 0 =step* (length - l); v s = eta * r s I t; II Initial values set according to jump conditions den 1 = (gam + 1) I (gam - 1) * out_den; vel_l = (2 + (gam - 1) * v_h) I (gam + 1) * v s; pres 1 = 2 I (gam + 1) *out den* pow(v_s * (1 - v_h), 2); II normalized values den[length - l] = l; vel[length - l] = 1; pres[length - l] = 1; II intially define gravity to have a lambda distribution for (int a = O; a < length; a++) grav[a] = step * a; vector <double> Compute der(double 1, double d, double v, double p, double g) { temp++; vector <double> ret(3); double v_log, d_log, p_log; double v_der, d_der, p_der; double theta p * pres_l Id I den 1 I pow(l * v_s, 2); double v n = v * v el 1 I 1 I v_s; double v_log num, v_log_den, d_log num, d log den, p_log num, p log den; double g_mine = omega * pow(v_h I 1, 2) // do uble g theirs= omega * pow(v_h, 2) * g I v_n; I pow(l, 3) * g; 31 v_log_num = 2. *theta* (k_p - 2. * gam * v n) + (1 - v_n) k_d) * v_n - g_mine); v log den= 2. * v n * (gam *theta - pow(l - v n, 2)); v=log-= v log_num 7 v_log_den; d_log_num = 2. * theta * (k_p - gam * k_d) I (1 - v_n) 2.) * k_d + k_p * v_n - 4. * v_n * (1 - v_n) - g_mine; d l og den= 2 . * (gam *theta - pow(l - v_n, 2)); d_log = d_log_num I d_log_den; * ((k_p - (3. * v n - p log_num = -gam * v_n * k_d + k_p * (gam * v n + 2. * (1 - v_n)) 4. * gam * v_n * (1 - v_n) - gam * g mine; p log den= 2. * (gam *theta - pow(l - v_n , 2)); p=log- p log_num I p_log_den; v der d der p_der v_l og * v I l; d_log * d I l; p_log * p I l; ret[OJ ret[l) ret[2) step * d_der; step * v_der; step * p_der; //cout //cout << endl; //cout << endl; //cout << endl; //cout << l << endl; << 11 v_log 11 << v_log_num << - 11 11 << v_log_den << 11 11 << v_log << 11 d_log 11 << d_log_num << 11 11 << d_log_den << 11 11 << d log << 11 p_log 11 << p_log_num << 11 11 << p log_den << 11 11 << p_log << endl; //if (temp % 5 //cin >> temp; 0) return ret; void Compute funcs() double lambda; vector <double > kl(3); vector <double> k2(3); //vector <double> k3(3); //vector <double> k4(3); int a = length - 2; lambda= lambda 0 - step* (length - 2 - a) ; while (lambda> .9 0) II for (int a= length - 2 ; a>= O; a--) 32 //lambda= lambda 0 - step* //cout << lambda<< endl; (length - 2 - a); kl= Compute der(lambda, den[a + 1), vel[a + 1), pres[a + 1), grav[a)); k2 = Compute_der(lambda - .5 * step, den[a + lJ - .5 * kl[OJ, ve l[a + lJ - .5 * kl[lJ, pres[a + lJ - .5 * kl[2J, grav[a + lJ); //k3 =Compute der(lambda - .5 * step, den[a + lJ - .5 * k2[0J, vel[a + lJ - .5 * k2[1), pres[a + lJ - .5 * k2[2J); //k4 =Compute der(lambda - step, den[a + lJ - k3[0J, vel[a + lJ k3[1), pres[a + 1) - k3[2)); //cout << k2[0J << endl; //cout << k2[1) << endl; //cout << k2[2) << endl; //cin >> temp; den[aJ = den[a + lJ - k2 [OJ; vel[aJ = vel[a + 1) - k2 (1 J; pres[aJ = pres[a + lJ - k2 [ 2); //den(aJ den[a + lJ 1 I 6. * (kl [ 0 l + 2 * k2 [OJ + 2 * k3 [OJ + k4 (OJ); 1 I 6. * (kl (1 J + 2 * k2(1J + 2 * k3[1J + I Ivel [ a J ve 1 [a + 1 J k4 [ 1 J); //pres[a) = pres[a + 1) - 1 I 6. * (kl[2) + 2 * k2[2J + 2 * k3[2J + k4 [2 J); a--; lambda lambda 0 - step* (length - 2 - a); vo id Determine l am() double lambda; int a = (int) (length * . 96); lambda = step * a; double dif = abs(vel[aJ - lambda* 24 I 23.); while ((abs(vel[aJ - lambda* 24 I 23.) <= dif) && (a< length)) { dif = abs(vel[aJ - lambda* 24 I 23.); a++; lambda= step* a;; lam i lambda - step; cout << l am i << endl; void Recompute() { 33 II fit proportionality constants; int a i int a f (int) (lam_i I step); a _ i + accuracy; if (a f >= length) { a f length - l; pres[a_f] I pow(step * (a f - a i), l_p); den[a f] I pow(step * (a f - a i), 1 d); (1 - ;el[a_f] * 23 I 24.-1 (step* a= f)) I c_p c d c v a (step* (a f i)) ; II recompute function from a i to a f for (int a a i; a <= a_f; a++) den[a] c_d * pow(step * (a - a i), l_d); pres[a] = c_p * p o w(step * (a - a_i), l_p); vel [a] = (1 - c v * step * (a - a i)) * step * a * 24 I 23.; v oid Determine_grav() double lam; double sum O; for (int a 0; a < length; a ++) { lam = step * a; if (a < (int) (length * lam_i)) grav[a] = O; else sum+= den[a] * pow(lam, 2) * step; grav[a] = sum; for (int b = (int) (length * lam_i); b < length; b++) { lam = step * b; grav[b] = grav [b] I grav[l e ngth - l] I p ow (lam, 2); 34 void Output va l s() double lambda ; of stream out_ A ( " den . dat " ) ; ofstream out_ B ( " vel . da t") ; ofstream out_ C( " pres . dat " ); o fs t ream out D( " grav .dat " ) ; for (int a = (int) (length * lam_i); a < l eng t h ; a++ ) l ambda step * a ; //if (abs(den [a ] ) < 5) ou t A<< lambda << " " << den[a] << end l; // i f (abs(vel[a]) < 5 ) ou t B << lambda << " " << ve l[ a] << endl; / / if (abs(pres[a]) < 5) out C << lambd a << " " << p r es [a ] << endl; out D << l ambda << " " << grav [a] << end l ; /la--; out_ A.close () ; out_B . close () ; out C . close(); out D. close(); v o i d Determine_ Energy() double therm = O; doub l e kin = O; double lambda; fo r ( i nt a = length - 1 ; a > (i nt) (lengt h * l am_i ) ; a -- ) lambda a * step ; therm += pres[a] * pow(lambda , 2) *step ; k i n+= den[a ] * p o w( v el [a] * lambda, 2) * step ; cout << " Thermal: " << t h erm I (t herm + ki n) << endl; cout << " Kinetic : " << kin I (the r m + k i n) << endl ; int main () 35 int iterate; //cout << "Starting value fort (after explosion): "; //cin >> time; cout << "Step size: "; cin >> step; cout << "# of steps: "; cin >> accuracy; cout << "# of iterations: "; cin >> iterate; Init_funcs (); Compute_funcs(); for (int a = 0; a < iterate; a++) Determine_lam(); //Recompute(); Determine_grav(); //concerned about my definition Compute funcs(); Determine_lam(); Recompute(); Determine_Energy(); Output_vals (); return 0; Code used for Explosion Model #include #include #include #include <iostream> <vector> <cmath> <fstream> using namespace std; II notice f - d and M- 5 vecto r <double> r; vector <double> r_prime; vec tor <doub le> q; double omega_b = 1 /7.; double omega_d = 6/7.; double t_fin, step; int length, temp; double f d; double f_m = .l; int a ve l = O; double t_burn = 5E7; double M 5 = omega_b * 20.; double H 0 = 3 .24E-18; double h = 0.5; double t hubble; void Init funcs() length = (int) (t fin I step) + l; 36 r.resize(length); r_prime.resize(length); q.resize(length); t fin= step* (length - l); r[OJ = 1.; r_prime[OJ q[OJ = 1.; O.; II doesnt matter l ong run vector <double > Compute der(double t, double R, double Rprime, double P) vector <double> ret(3); double rDPrime, qPrime, eta, grav, l diss; eta=ll(l+t); grav = -(2 I 9. *omega d + 1 I 9. * omega_b) * pow(eta, 2) * R; rDPrime = 18 * M PI I omega b I pow(eta,2) * P I R - 3 . * (pow(Rprime, 2) - 4 I 3. *eta* R * Rprime + 4 I 9. * pow(eta * R, 2)) I R + grav; 1 I 3. * f d * omega_b * pow(eta * R, 2) * pow(Rprime - 2 I l diss 3. * eta * R, 3) ; Ill diss O; II changed to make NOT constant energy. if (t hubble * t < t_burn) l diss += 1; qPrime = l diss I 2. I M_PI I pow(R, 3) - 5 * Rprime I R * P; llcout << l diss I 2. IM PI I pow(R, 3) I (5 * Rprime I R * P) endl ; ret[O] ret[l) ret[2) << step * rDPrime; step * Rprime; step * qPrime; llcout << t << 11 11 << rDPrime << endl; return ret; void Compute funcs() double t; vector <double> k1(3); vector <double> k2(3); 37 for (int a = l; a < length; a++) { t = (a - 1) * step; kl= Compute_der(t, r[a-1], r_prime[a-1], q[a-1]); k2 = Compute_der(t + . 5 * step, r[a-1] + .5 * kl[l], .5 * kl[OJ, q[a-1] + .5 * kl(2]); r_prime[a-1] + = r [a - l ] + k2 [ l] ; r_prime[a] r_prime[a - l] + k2[0] ; q [a J = q [a - l J + k2 [ 2 J ; r [a] //if ((q[a] > q[a-1]) II q[aJ = O; 11 q[a-1] 0) I /cout << "Z: 11 << z << endl; //cout << 11 R Doulbe prime: 11 << rDPrime << endl; //cout << 11 R prime: 11 << r_prime[a - l] << endl; //cout << " Q prime: 11 << qPrime << end l; //cin >> temp; //cout << z << 11 11 << rDPrime << endl; void Determine_vel() double double double double for eta = l; diff min t; t vel; abs(2 I 3 . *eta* r[OJ I r_prime[OJ (int a = l; a < length; a++) t = a * step; eta= 1 I (l + t); if (abs(2 I 3. *eta* r[a] a vel = a; diff min= abs(2 t vel = t; - 5 I 6 .); { I r_prime[a] - 5 I 6 .) < diff_min) I 3. *eta* r[a] I r_prime[a] - 5 I 6 .); // cout << 11 a 11 << a ve l << endl; cout << 11 t vel 11 << t vel << endl; void Output_vals() 38 double t, eta, l_diss; double e_init = 2. * M PI * pow(r[O], 3) * q[O] + omega_b I 9. * pow(r[O], 3) * pow(r_prime[O], 2); double e last; doubler factor= .1 3 I h I (1 + t fin) * pow(M_S, .2); ofstream out_A("r.dat"); ofstream out_B( 11 q.dat 11 ) ; o f st r eam out_C( 11 log.dat"); ofstream out D("slope log.dat 11 ) ; ofstream out_E( 11 l_vs_w.dat 11 ) ; ofstream out F( 11 therm.dat 11 ) ; ofstream out_G( 11 kin.dat 11 ) ; ofstream out_H ( 11 rel_vel.dat"); ofstream out I( 11 temp.dat 11 ) ; o fstre am out_J( 11 slope.dat 11 ) ; double temp; double temp_l = 10000; int counter = O; for (int a = O; a < length; a++) t = a * step; //if (abs(r[a]) < 1000) out A << t << 11 11 << r[a] * r fa cto r << endl; //if (abs (q [a]) < 1000) out B << t << " " << q[a] << endl; out C << log(t + 1) <<" "<< log(r[a]) << endl; if (a != length - 1) out D << t « 11 11 << (log(r[a + l]) - log(r[a])) I step) - log(t + 1)) << endl; out J << t << 11 11 (log(t + 1 + << r_prime[a] << endl; if (a 2) length cout <<"scaling 11 << (log(r[a + l]) - log(r[a])) I step) - log(t + 1)) << endl; (log(t + 1 + } eta = 1 I (1 + t); 1 diss = 1 I 3. * f d * omega_b * pow(eta * r[a], 2) * pow(r_prime[a] - 2 I 3. *eta* r[a], 3); out E << t << 11 11 << 1 diss - 4 * M PI* q[a] * r_prime[a] pow(r[a], 2) << endl; * out F << t << " << 2. * M PI * pow(r[a], 3) * q[a] I e init << endl; << omega_b I 9. * pow (r [a], 3) * pow(r_prime[a) * out - G << t << eta, 2) I e init << endl; II II II if ((a== length - 1) 11 (a a_vel)) 39 e_last = 2. * M_PI * pow(r[a], 3) * q[a] + omega_b I 9. * pow(r[a], 3) * pow(r_prime[a] *eta , 2); cout << "Thermal: "<< 2. * M PI * pow(r[a], 3) * q[a] I e init << endl; cout <<"Kinetic: "<< omega_b I 9. * pow(r[a], 3) * pow(r_prime[a] * eta, 2) I e _init << endl ; if (a == a_vel) { cout <<"scaling"<< (log(r[a + 1]) - log(r[a])) I (log(t + 1 + step) - log(t + 1)) << endl; } cout << endl; out H << t << " "<< r_prime[a] I (2 I 3. *eta* r[a]) << endl; temp= 4.5E5 * pow(M_S, .4) * q[a] I f m I omega_b I pow(eta, 2); out I << t << " " << temp << endl; if ((counter== 0) & (temp< temp_l)) cout <<"t of temperature: " << t << endl; cout << endl; counter++; out_A . close(); out B.close(); out C . close(); out_D.c l ose (}; out_E.close (}; out_F.close(); out_ G.close (}; out_H.close() ; out I.close(); out J.close(); int main() cout << "End val of t: "; cin >> t fin; cout << "St ep size: "; c in >> step; cout << "f d: "; cin >> f_d; //cout << "T burn: "; //cin >> t_burn; cout << "I nitial z: " << pow(l + t fin, wrong 2/3.) - 1 << endl; //possible 40 cou t << endl ; t h ubb l e = 2 I 3 . I (H_ O * h ) * pow( l + t _ f i n , - 1 ); cout << " t hubb l e : " << t h ubble << e ndl; cou t << e ndl; I n i t _f uncs (); Compute_funcs() ; Dete r mine_ve l(); Output_ v a ls (); re tu r n O; References 1) Astrophysical Blastwaves, Reviews of Modem Physics, Vol. 60, No. 1, January 1988. 2) The Astrophysical Journal, 417:54-62, 1993 November 1 3) Shu, Frank H. The Physics of Astrophysics: Volume 2 Gas Dynamics. University Science Books, 1992. 41