Extend model of steady state globule flows inside a Stromgren sphere to include the case with an inner hot bubble

[will-henney]

This will be a slight modification of the model given in section 3 of Henney (2013)

https://ui.adsabs.harvard.edu/abs/2003RMxAC..15..175H/abstract

• [x] First task is for @RobeReyes to read this section and understand it
• [ ] Then we want to extend it to include the interior hot bubble
• [ ] And then apply it to M 17 and possibly other regions

[will-henney]

Roberto has read the earlier sections of this paper, which presents an approximate model for the photoevaporation flow from a plane ionization front. In a later paper Henney et al (2005) we improved on this model in the Appendix A, so he should read that instead.

But, @RobeReyes still need to study section 3 of the 2003 paper since we are going to apply a similar model to the pressure balance of the photoevaporation flow with the hot gas (see #6)

[will-henney]

Notes from blackboard about calculating the pressure balance

[will-henney]

@RobeReyes will need to use equation (14) from that paper, together with the continuity equation and the equation for the total pressure (static thermal pressure plus ram pressure), which are all given on the blackboard image. We want to find the radius where the total pressure has fallen by a factor of four, since that is where we will have equilibrium with the hot gas.

[will-henney]

Roberto still need to finish this calculation:

This is the total dynamic pressure in the flow

We want it to have fallen to 0.25 times its initial value, so it is in equilibrium with the hot gas. This will be at the point "1", whereas the ionization front is point "0"

The definition of the total pressure gives us 1 equation

The second equation comes from mass conservation

And a third equation from the isothermal Bernoulli principle

We combine these equations to get a single equation for the Mach number at point 1

@RobeReyes should solve this equation by writing a python program. We can use scipy.optimize.fsolve to find $M_1$

The photoevaporation flow will form a shock, which we will assume is isothermal. We define Point 2 as being just outside this shock in a dense shell. We also want to find the density there.

Then later we can calculate the surface brightness of the shell in terms of the surface brightness of the base of the photoevaporation flow. For this we will also need the shell thickness

[RobeReyes]

To compute the values of $M_1$ and $r_1$ on equilibrium point, first we calculate the Mach number $$\frac{P}{P_0}=\frac{1+M^2}{2}e^{\frac{1-M^2}{2}}$$ and then we calculate $r_1$ $$\frac{r}{r_0}=M^{-1/alpha}e^{\frac{M^2-1}{2\alpha}}.$$

We considerate that $\frac{P}{P_0}=0.25$, but in the Table we considerate others ratios too.

	P/P_0	M_1	R/R_0
0.500000	[2.0872695034502167]	[1.6019876154558883]
0.333333	[2.3618981411410918]	[2.043998443517419]
0.250000	[2.5269089927774195]	[2.4177052766816747]
0.200000	[2.643599193146395]	[2.7485775238564956]
0.166667	[2.733198469715219]	[3.049132752352427]
0.142857	[2.8055534914820472]	[3.3267045090417726]
0.125000	[2.8660158884661424]	[3.586046462195688]
0.111111	[2.917807287745739]	[3.8304551698598424]
0.100000	[2.9630120384053216]	[4.062330887863703]
0.090909	[3.0030532411044946]	[4.2834871526475595]
0.083333	[3.0389440079280536]	[4.495334702416166]
0.076923	[3.0714301227600624]	[4.698997105639484]
0.071429	[3.1010758409032033]	[4.895386810385624]
0.066667	[3.128317957368406]	[5.085257042831157]

[RobeReyes]

Now for the density on the point 2 we have that $M_1 M_2=1 \Rightarrow M_2=M_1^{-1}$ and from $\rho_2 M_2=\rho_1 M_1\Rightarrow \frac{\rho_2}{\rho_1}=M_1^2$.

For the shock shell thickness we use $\frac{h}{r_1}=\frac{3}{4 M_1^2}$. .

	P/P_0	M_1	\rho_2/\rho_1	h/r_1
0.500000	[2.0872695034502167]	[4.356693980033314]	[0.17214888248686797]
0.333333	[2.3618981411410918]	[5.578562829125745]	[0.13444322901307856]
0.250000	[2.5269089927774195]	[6.385269057779393]	[0.117457853884207]
0.200000	[2.643599193146395]	[6.98861669400427]	[0.10731737521725093]
0.166667	[2.733198469715219]	[7.470373874853615]	[0.10039658155860326]
0.142857	[2.8055534914820472]	[7.871130393567106]	[0.09528491620631234]
0.125000	[2.8660158884661424]	[8.214047072940371]	[0.09130700047613965]
0.111111	[2.917807287745739]	[8.513599368422145]	[0.08809434970381981]
0.100000	[2.9630120384053216]	[8.779440339734858]	[0.08542685763300605]
0.090909	[3.0030532411044946]	[9.01832876890821]	[0.08316396742883411]
0.083333	[3.0389440079280536]	[9.235180683321822]	[0.08121118857527657]
0.076923	[3.0714301227600624]	[9.433682998997892]	[0.07950235343711146]
0.071429	[3.1010758409032033]	[9.61667137103351]	[0.07798956323485108]
0.066667	[3.128317957368406]	[9.786373242393635]	[0.0766371751233717]

[will-henney]

Next steps:

1. Include higher values of P/P0, starting from 1
2. Calculate the shell density in terms of the initial density: $\rho_2 / \rho_1$
3. Compare the surface brightness of the shell with the surface brightness of the ionization front

[will-henney]

The surface brightness is proportional to the emission measure, which is the integral of the density along the line of sight (assuming fully ionized hydrogen).

For the shell, we will assume constant density $n_2$ so that the EM is $n_2^2 \ell$, where $\ell$ is the path length along the line of sight through the shell.

The maximum value of $\ell$ can be found from simple geometry:

We can use the same technique for the EM from the base of the flow, but using the effective thickness of the photoevaporation flow.

@RobeReyes should do a numerical calculation to check the value of $h_0 / r_0$. I think the answer should be about 0.12

[RobeReyes]

To calculate $\frac{h}{r_o}$ and $\frac{\rho_2}{\rho_0}$ we use the previous equations,

$$\frac{h}{r_0}=\frac{h}{r_1}\frac{r_1}{r_0}=\frac{3 M^{-1/\alpha}}{4 M^2}e^{\frac{M^2 -1}{2\alpha}}$$

and

$$\frac{\rho_2}{\rho_0}=\frac{\rho_2}{\rho_1}\frac{\rho_1}{\rho_0}=M^2 e^{\frac{M^2 -1}{2}}$$

and we have

P_P_0_1 = [i+2 for i in range(14)]
Extr=[1/0.9,1/0.8,1/0.7,1/0.6]
P_P_0=Extr+P_P_0_1

Roots = []

for i in P_P_0:
    def func(x):

        return ((1+x**2)/2)*(e**((1-x**2)/2))-1/(i)

    root = fsolve(func, 1.5)

    Roots.append(root)

alpha=2

h_R_0 = [(i**(-1/alpha))*(3/(4*i**2))*e**((i**2-1)/(2*alpha)) for i in Roots]
rho2_rho0 = [(i**2)*e**((1-i**2)/2) for i in Roots]

T3 = pd.DataFrame(data = [[1/P_P_0[i],Roots[i],rho2_rho0[i],h_R_0[i]] for i in range(len(P_P_0))],
                columns = ['P/P_0','M_1', '\rho_2/\rho_0','h/r_0'])

This time we take values of $P/P_0\approx 1$.

	P/P_0	M_1	\rho_2/\rho_0	h/r_0
0	0.900000	[1.43653166229611]	[1.212460386727212]	[0.3955985401577511]
1	0.800000	[1.62750625745831]	[1.161496861145731]	[0.335172394947937]
2	0.700000	[1.7873719314700518]	[1.066245374672153]	[0.30395569363387]
3	0.600000	[1.9372255119437625]	[0.9475195204739406]	[0.28575702167852984]
4	0.500000	[2.0872695034502167]	[0.8133176911491646]	[0.2757803777585336]
5	0.333333	[2.3618981411410918]	[0.565327410056532]	[0.2748017508441885]
6	0.250000	[2.5269089927774195]	[0.43229765955875116]	[0.2839784731235524]
7	0.200000	[2.643599193146395]	[0.3499287529591649]	[0.29497012544141]
8	0.166667	[2.733198469715219]	[0.2939804856009608]	[0.3061225050545589]
9	0.142857	[2.8055534914820472]	[0.25350708403438516]	[0.31698476038720674]
10	0.125000	[2.8660158884661424]	[0.2228675143483692]	[0.3274311460311606]
11	0.111111	[2.917807287745739]	[0.1988638471617976]	[0.3374414572584375]
12	0.100000	[2.9630120384053216]	[0.17954893193761]	[0.34703216241569557]
13	0.090909	[3.0030532411044946]	[0.16366962770180535]	[0.356231786044611]
14	0.083333	[3.0389440079280536]	[0.1503829612956301]	[0.36507147422690406]
15	0.076923	[3.0714301227600624]	[0.13910101027020574]	[0.37358132869251404]
16	0.071429	[3.1010758409032033]	[0.12940121699633922]	[0.38178907920762556]
17	0.066667	[3.128317957368406]	[0.12097205764438175]	[0.3897197345388107]

[will-henney]

We now have the density and thickness of the shell as a function of the pressure ratio.

• We find that the thickness of the shell is roughly constant, equal to 0.3 to 0.4 times r0.

• It is interesting that for pressure ratios close to unity that the shell is denser than the photoevaporation flow.

• [ ] However, this means that the recombinations in the shocked shell will be important. So, that is another thing that we should calculate.

  • The recombinations in the shell are given by $\alpha_B h_2 n_2^2$, which needs to be compared with the recombinations in the photoevaporation flow, $\alpha_B h_0 n_0^2$

• [ ] Also, in order to calculate the surface brightness of the shell, we need the line of sight depth of the shell: $\ell_2 = 2(r_2 h_2)^{1/2}$

  • This can be compared with the surface brightness of the photoevaporation flow. The ratio will be $n_2^2 \ell_2 / n_0^2 \ell_0$ where $\ell_0 = 2(r_0 h_0)^{1/2}$

• [ ] We are also still missing the calculation of $h_0 / r_0$

@RobeReyes should calculate all these quantities

[will-henney]

Note that for the third item in the previous comment, the definition of $h_0$ is as follows:

[RobeReyes]

To calculate $h_0$ we use the equation for r and $rho$ $$\frac{\rho}{\rho_0}=e^{\frac{1-M^2}{2}}$$ and $$\frac{r}{r_0}=M^{-1/\alpha}e^{\frac{M^2-1}{2 \alpha}}$$ hence $$\frac{dr}{dM}=r_0\Big( e^{\frac{M^2-1}{2\alpha}}\Big(-\frac{M^{-1/\alpha-1}}{\alpha}+\frac{M}{\alpha M^{1/\alpha}}\Big)\Big)=\frac{r_0}{\alpha}e^{\frac{M^2-1}{2\alpha}}\Big(M^{\frac{\alpha-1}{\alpha}}-M^{-\frac{\alpha+1}{\alpha}}\Big)$$ therefore $$h_0=\int_0^\infty \Big(\frac{n(r)}{n_0}\Big)^2 dr=\int_1^\infty \frac{r_0}{\alpha} (exp(1-M^2))^{1-\frac{1}{2\alpha}}\Big(M^{\frac{\alpha-1}{\alpha}}-M^{-\frac{\alpha+1}{\alpha}}\Big) dM.$$

For $\alpha=2$ we have $$h_0=r_0\int_1^\infty \frac{exp(\frac{3}{4}(1-M^2))}{2}(M^{1/2}-M^{-3/4})dM$$

[RobeReyes]

for the recombinations we have $$\frac{\alpha_\beta h_2 n2^2}{\alpha\beta h_0 n_0^2}=\frac{(3/4)M^{-1/\alpha-2}exp(\frac{M^2-1}{2\alpha})r_0}{0.12 r_0}\Big(M^2 exp(\frac{M^2-1}{2})\Big)^2=\frac{3}{4*0.12}M^{2-1/\alpha}exp(\frac{M^2-1}{2\alpha})exp(M^2-1).$$

and for the surface brightness of the shell we have $$\frac{n_2^2 l_2}{n_0^2 l_0}=\Big(\frac{n_2}{n_0}\Big)^2\Big(\frac{r_2 h_2}{r_0 h_0}\Big)^{1/2}=\Big(M^2 exp(\frac{M^2-1}{2})\Big)^2\Big(\frac{3}{4*0.12}M^{-2-1/\alpha} exp(\frac{M^2-1}{2\alpha})\frac{r_2}{r_0}\Big)^{1/2}$$

[will-henney]

We have revised the notebook and found a cleaner way of doing the calculations without loops. @RobeReyes should finish this and include the relevant equations from above, and a logical description of what we have done.

[will-henney]

@RobeReyes has made a figure of the emission measures, but we want to expand this to include the other columns in the table. Also, possibly use a log scale on the vertical axis.

• [x] And we could even extend the range of pressure contrasts to be from 0.01 to 1 (and use a finer grid of points).

[will-henney]

• [ ] Also, remove the empirical calculations for LDN 1616 and put them in a separate notebook. For the densities, we should first give the empirical determination of the density from the Ha surface brightness.

• [ ] Then, for the theoretical version, we should invert the calculation. Given the density, we can calculate the required flux of ionizing photons from the three contributions: new ionizations, recombinations in the photo evaporation flow, recombinations in the shocked shell

[will-henney]

@RobeReyes has made some new figures. For example

But we still want to use a larger range and a finer grid. I suggest 100 points in P/P0 entre 0.01 y 0.99

[will-henney]

• [x] We also want to add a new column: $L_2 / L_0$, which is the ratio of the Ha luminosity of the shell to the photo evaporation flow.

This is illustrated in the above blackboard

@RobeReyes will have to evaluate the integral numerically to find $\xi_0$, which is similar to the previous calculation that we did to find $h_0$, except this time we have an extra factor of $x^2$ inside the integral