# Stromgren Sphere Examples¶

In this section we will go thourgh some concrete examples using SphereStromgren. This class models a spherically symmetric gas distribution with a point source at the center. In particular we will examine Tests 1 and 2 in the Cosmological Radiative Transfer Comparison Project (CRTCP) [Iliev06], the work of [Friedrich12] which includes helium, and the work of [Raicevic14] which includes radiative transfer of recombination radiation.

## Preparation¶

### Import Packages¶

The following examples require three python packages to be imported,

import numpy as np
import pylab as plt
import rabacus as ra


### Plotting Functions¶

We will be making several plots during the course of these examples. In order to facilitate this, we will define a few convenience functions here,

def plot_sphere_x( s, fname, style, plot_H=True, plot_He=True ):
""" put plot of ionization fractions from sphere s into fname """

plt.figure( figsize=(7,7) )

s.r_c.units = 'kpc'
s.Edges.units = 'kpc'
if style == 'I':
xx = s.r_c / s.Edges[-1]
else:
xx = s.r_c

if plot_He:
plt.plot( xx, np.log10( s.xHe1 ),
color='green', ls='-', label = r'$x_{\rm HeI}$' )
plt.plot( xx, np.log10( s.xHe2 ),
color='green', ls='--', label = r'$x_{\rm HeII}$' )
plt.plot( xx, np.log10( s.xHe3 ),
color='green', ls=':', label = r'$x_{\rm HeIII}$' )

if plot_H:
plt.plot( xx, np.log10( s.xH1 ),
color='red', ls='-', label = r'$x_{\rm HI}$' )
plt.plot( xx, np.log10( s.xH2 ),
color='red', ls='--', label = r'$x_{\rm HII}$' )

if style == 'I':
plt.xlim( -0.05, 1.05 )
plt.xlabel( 'r_c / Rsphere' )
else:
plt.xlim( -0.50, 15.2 )
plt.xlabel( 'r_c [kpc]' )

if style == 'I':
plt.ylim( -6.0, 0.5 )
else:
plt.ylim( -3.1, 0.1 )
plt.ylabel( 'log 10 ( x )' )

plt.grid()
plt.legend(loc='lower center', ncol=2)
plt.tight_layout()
plt.savefig( 'doc/img/x' + fname )

def plot_sphere_T( s, fname, style ):
""" put plot of temperature from sphere s into fname """

plt.figure( figsize=(7,7) )

s.r_c.units = 'kpc'
s.Edges.units = 'kpc'
if style == 'I':
xx = s.r_c / s.Edges[-1]
else:
xx = s.r_c

plt.plot( xx, s.T,
color='black', ls='-', label = r'$T$' )

if style == 'I':
plt.xlim( -0.05, 1.05 )
plt.xlabel( 'r_c / Rsphere' )
else:
plt.xlim( -0.50, 15.2 )
plt.xlabel( 'r_c [kpc]' )

plt.ylim( 0.0, 5.0e4 )
plt.ylabel( 'T [K]' )

plt.grid()
plt.legend(loc='best')
plt.tight_layout()
plt.savefig( 'doc/img/T_' + fname )


### Create Sources¶

Next we will setup three point sources, one monochromatic, one with a thermal spectrum, and one with a powerlaw spectrum. We will normalize all spectra such that they emit the same number of photons per second,

Ln = 5.0e48 / ra.u.s  # set photon luminosity

q_mono = 1.0
q_min = q_mono
q_max = q_mono
src_mono = ra.PointSource( q_min, q_max, 'monochromatic' )
src_mono.normalize_Ln( Ln )

q_min = 1.0
q_max = 10.0
T_eff = 1.0e5 * ra.u.K
src_thrm = ra.PointSource( q_min, q_max, 'thermal', T_eff=T_eff )
src_thrm.normalize_Ln( Ln )

q_min = 1.0
q_max = 10.0
alpha = -1.0
src_pwr1 = ra.PointSource( q_min, q_max, 'powerlaw', alpha=alpha )
src_pwr1.normalize_Ln( Ln )


### Solution Container¶

We also initialize a dictionary to hold our results,

spheres = {}


## Iliev06 Examples¶

### Setup¶

To begin, we define a sphere as descibed in [Iliev06]. Note that we take the helium density from [Friedrich12] and also define a null helium densiy which is much lower in order to approximate a zero helium environment.

Nl = 512
T = np.ones(Nl) * 1.0e4 * ra.u.K

Rsphere = 6.6 * ra.u.kpc
Edges = np.linspace( 0.0 * ra.u.kpc, Rsphere, Nl+1 )
nH = np.ones(Nl) * 1.0e-3 / ra.u.cm**3
nHe = np.ones(Nl) * 8.7e-5 / ra.u.cm**3
nHe_null = np.ones(Nl) * 1.0e-15 / ra.u.cm**3


### Optically Thin¶

We begin with a very basic scenario, a pure hydrogen sphere with constant density and temperature. We place a monochromatic source at the center but allow only geometric dillution of the radiation by setting the keyword thin to True. By default, case A recombination rates are used. We use case B recombination rates by setting the keyword fixed_fcA to 0.0.

key = 'thin_caseB_mono_fixT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, fixed_fcA=0.0, thin=True )

plot_sphere_x(
spheres[key], 'strm_sphere_' + key + '.png', 'I', plot_He=False )


Stromgren Sphere - Thin - Case B - Mono - Fixed T

In this case, the gas remains ionized at all radii.

### Test 1¶

Test 1 of [Iliev06] also involves a pure hydrogen sphere with constant density and temperature and a central monochromatic source. We call the solver with the same input except this time we leave out the keyword thin so the radiation will be attenuated by the gas. The resulting figure can be compared directly to the right panel of Fig. 8 in [Iliev06] (reproduced below).

key = 'rt_caseB_mono_fixT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, fixed_fcA=0.0 )

plot_sphere_x(
spheres[key], 'strm_sphere_' + key + '.png', 'I', plot_He=False )


Stromgren Sphere - RT - Case B - Mono - Fixed T

[Iliev06] - Fig. 8

### Test 2¶

Test 2 of [Iliev06] uses a 1.0e5 K thermal spectrum and allows the gas temperature to vary. We model this situation by setting the find_Teq keyword to True and passing in src_thrm instead of src_mono. Note that when find_Teq is True we also have to pass in a redshift using the keyword z so that Compton cooling can be accounted for. These figures can be compared directly to the rightmost panels of Figs. 16 and 17 in [Iliev06] (reproduced below).

key = 'rt_caseB_thrm_evoT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_thrm, fixed_fcA=0.0, find_Teq=True, z=0.0 )

plot_sphere_x(
spheres[key], 'strm_sphere_' + key + '.png', 'I', plot_He=False )
plot_sphere_T(
spheres[key], 'strm_sphere_' + key + '.png', 'I' )


Stromgren Sphere - RT - Case B - Thermal - Equilibrium T

[Iliev06] - Fig. 16

Iliev Sphere - RT - Case B - Thermal - Equilibrium T

[Iliev06] - Fig. 17

## Friedrich12 Examples¶

We now focus on the tests presented in [Friedrich12]. They are based on those described in [Iliev06] except they include helium. Because of the longer mean free path of helium ionizing photons, the radius of the sphere is increased to 15 kpc.

Rsphere = 15.0 * ra.u.kpc
Edges = np.linspace( 0.0 * ra.u.kpc, Rsphere, Nl+1 )


### Test 1 A¶

Test 1 A, examines a fixed temperature sphere, uses a thermal spectrum, and case A recombination rates. The Rabacus figure can be compared to the panels in Fig. 4 of [Friedrich12] (reproduced below). Because Rabacus produces equilibrium solutions, the lines in this figure should technically be compared to the CLOUDY results in [Friedrich12].

key = 'rt_caseA_thrm_fixT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe, src_thrm, fixed_fcA=1.0 )

plot_sphere_x( spheres[key], 'strm_sphere_' + key + '.png', 'F' )


Stromgren Sphere - RT - Case A - Thermal - Fixed T

[Friedrich12] - Fig. 4

Note

The maximum energy of photons considered when constructing the spectra in [Friedrich12] is unclear. We have used q_max = 10, however this choice can have order of magnitude effects on the ionization fractions. For example, try the powerlaw example with q_max = 60 instead of 10 and plot the results.

### Test 1 B¶

Test 1 B, examines a fixed temperature sphere, uses a thermal spectrum, and case B recombination rates (i.e. the on-the-spot approximation). We note that using a case B rate for all ionic species is equivalent to what [Friedrich12] term the U-OTS or the uncoupled on-the-spot approximation. This figure can be compared to the panels in Fig. 5 of [Friedrich12] (reproduced below). Again our equilibrium solutions should technically be compared to the CLOUDY OTS solutions in that plot.

key = 'rt_caseB_thrm_fixT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe, src_thrm, fixed_fcA=0.0 )

plot_sphere_x( spheres[key], 'strm_sphere_' + key + '.png', 'F' )


Stromgren Sphere - RT - Case B - Thermal - Fixed T

[Friedrich12] - Fig. 5

### Test 2¶

Test 2 in [Friedrich12] allows the temperature to vary and examines ionization profiles in the case of powerlaw sources. They present results for spectra with powerlaw indices of and . Here we will only reproduce the solution. This figure can be compared to the panels in Fig. 6 of [Friedrich12] (reproduced below). Our equilibrium solutions should technically be compared to the CLOUDY OTS solutions in that plot. The dependence on the maximum photon energy included in the spectrum q_max is stronger for the powerlaw source than for the thermal source.

key = 'rt_caseB_pwr1_evoT'

spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe, src_pwr1, fixed_fcA=0.0, find_Teq=True, z=0.0 )

plot_sphere_x( spheres[key], 'strm_sphere_' + key + '.png', 'F' )
plot_sphere_T( spheres[key], 'strm_sphere_' + key + '.png', 'F' )


Stromgren Sphere - RT - Case B - Powerlaw - Equilibrium T

[Friedrich12] - Fig. 6

Stromgren Sphere - RT - Case B - Powerlaw - Equilibrium T

## Raicevic14 Examples¶

Now we focus on the treatment of recombination radiation. [Raicevic14] present the results of treating this radiation explicitly using a 3-D radiative transfer code. We will solve the problem by transporting this radiation in our 1-D geometries. First we resize the sphere back to the original CRTCP size,

Rsphere = 6.6 * ra.u.kpc
Edges = np.linspace( 0.0 * ra.u.kpc, Rsphere, Nl+1 )


The first test involves a simple comparison between using fixed case A rates, fixed case B rates, or following the recombination radiation explicitly. We do not have a fixed case A sphere yet so we solve that first,

key = 'rt_caseA_mono_fixT'
spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, fixed_fcA=1.0 )


Rabacus has three options for the treatment of recombination radiation. In the first, all recombination radiation is assumed to be transported radially outward. This makes the assumption that the ionized center of the sphere produces zero optical depth. This option can be activated by setting the keyword rec_meth to outward.

key = 'rt_outward_mono_fixT'
spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, rec_meth='outward' )


The second option transports recombination radiation along radial lines as well, but accounts for the optical depth encounterd by inward travelling radiation. This option can be activated by setting the keyword rec_meth to radial.

key = 'rt_radial_mono_fixT'
spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, rec_meth='radial' )


The third and most realistic option transports recombination radiation isotropically from each layer. This option can be activated by setting the keyword rec_meth to isotropic.

key = 'rt_isotropic_mono_fixT'
spheres[key] = ra.SphereStromgren(
Edges, T, nH, nHe_null, src_mono, rec_meth='isotropic' )


Below we produce a plot that can be compared to the upper right panel in Fig. 2 of [Raicevic14] (reproduced below).

plt.figure( figsize=(7,7) )

key = 'rt_caseA_mono_fixT'
s = spheres[key]
s.r_c.units = 'kpc'
s.Edges.units = 'kpc'
xx = s.r_c / s.Edges[-1]

key = 'rt_caseA_mono_fixT'
s = spheres[key]
plt.plot( xx, np.log10( s.xH1 ),
color='blue', ls=':', label = r'case A' )

key = 'rt_caseB_mono_fixT'
s = spheres[key]
plt.plot( xx, np.log10( s.xH1 ),
color='black', ls='--', label = r'case B' )

key = 'rt_outward_mono_fixT'
s = spheres[key]
plt.plot( xx, np.log10( s.xH1 ),
color='cyan', ls='-', label = r'outward' )

s = spheres[key]
plt.plot( xx, np.log10( s.xH1 ),
color='green', ls='-', label = r'radial' )

key = 'rt_isotropic_mono_fixT'
s = spheres[key]
plt.plot( xx, np.log10( s.xH1 ),
color='red', ls='-', label = r'isotropic' )

plt.xlim( -0.05, 1.05 )
plt.xlabel( 'r_c / Rsphere' )

plt.ylim( -6.0, 0.5 )
plt.ylabel( 'log 10 ( x )' )

plt.grid()
plt.legend(loc='best', ncol=2)
plt.tight_layout()
plt.savefig( 'doc/img/x_raicevic.png' )


Stromgren Sphere - RT - Case A vs Case B vs Transport

[Raicevic14] - Fig. 2

The outward and radial options produce similar results and isotropic is the same as [Raicevic14].

### Diffuse vs. Source¶

In the next Figure we show the radiation intensity from the central point source (source) compared to that from recombinations (diffuse). The photoionization rates due to the central point source and the diffuse recombinations are both stored in the returned solved object. The following plots can be compared to the lower right panel in Fig. 2 of [Raicevic14] (reproduced below). Here we show the results of using the isotropic keyword. Note that the x-axis is scaled by the stromgren radius = 5.4 kpc and not the box length.

plt.figure( figsize=(7,7) )

key = 'rt_isotropic_mono_fixT'
s = spheres[key]
s.r_c.units = 'kpc'
s.Edges.units = 'kpc'
xx = s.r_c / (5.4*ra.u.kpc) # s.Edges[-1]

geo = 4.0 * np.pi * s.r_c**2
I0 = s.H1i_src[0] * geo[0]

# source I/I0
#------------------------------------
Is = s.H1i_src * geo
yy = Is / I0
plt.plot( xx, yy,
color='red', ls='-', label='source' )

# recomb I/I0
#------------------------------------
Id = s.H1i_rec * geo
yy = Id / I0
plt.plot( xx, yy,
color='red', ls='--', label='diffuse' )

# ratio diffuse/source
#------------------------------------
yy = Id / Is
plt.plot( xx, yy,
color='red', ls=':', label='ratio' )

plt.ylabel( r'$I/I_0$', fontsize=25 )
plt.ylim( 0.0, 1.05 )

plt.xlabel( 'r_c / R_strom' )
plt.xlim( -0.05, 1.05 )

plt.grid()
plt.legend(loc='best', ncol=1)
plt.tight_layout()
plt.savefig( 'doc/img/I_raicevic_isotropic.png' )


Stromgren Sphere - Source vs Diffuse photons - Isotropic

[Raicevic14] - Fig. 2

Making this plot using the sphere solved with the radial model, exposes the differences between isotropic (and also outward) and radial. While the results are similar for all three models at large radii, the diffuse radiation makes a larger contribution at smaller radii in radial. This is due to the fact that the radial method considers the opacity encountered by inward travelling recombination radiation.

plt.figure( figsize=(7,7) )

s = spheres[key]
s.r_c.units = 'kpc'
s.Edges.units = 'kpc'
xx = s.r_c / (5.4*ra.u.kpc) # s.Edges[-1]

geo = 4.0 * np.pi * s.r_c**2
I0 = s.H1i_src[0] * geo[0]

# source I/I0
#------------------------------------
Is = s.H1i_src * geo
yy = Is / I0
plt.plot( xx, yy,
color='green', ls='-', label='source' )

# recomb I/I0
#------------------------------------
Id = s.H1i_rec * geo
yy = Id / I0
plt.plot( xx, yy,
color='green', ls='--', label='diffuse' )

# ratio diffuse/source
#------------------------------------
yy = Id / Is
plt.plot( xx, yy,
color='green', ls=':', label='ratio' )

plt.ylabel( r'$I/I_0$', fontsize=25 )
plt.ylim( 0.0, 1.05 )

plt.xlabel( 'r_c / R_strom' )
plt.xlim( -0.05, 1.05 )

plt.grid()
plt.legend(loc='best', ncol=1)
plt.tight_layout()