Modelling HII regions. RT with RADMC-3D

With the Model module you can create electronic density distributions and use the Radmc3d and/or Radmc3dRT classes to obtain formatted data files that may be used later to predict the emission of your model with RADMC-3D.

Example 1

Source codes and figures on GitHub: ionized_constant

Note

Uniform spherical HII region

Model: Uniform density and constant temperature

Radiative transfer: free-free emission

Useful references: Pratap+1992, Keto+2008

The preamble:

#------------------
#Import the package
#------------------
from sf3dmodels import Model, Plot_model
import sf3dmodels.utils.units as u
import sf3dmodels.rt as rt
#-----------------
#Extra libraries
#-----------------
import numpy as np
import time

a. Define the general parameters and the GRID:

Note

The radmc3d flag must be turned on at the GRID definition (see below)

#------------------
#General Parameters
#------------------
r_max = 2530 * u.au #HII sphere size
dens_e = 1.4e5 * 1e6 #Electronic numerical density, from cgs to SI
t_e = 1.e4 #K

#---------------
#GRID Definition
#---------------
sizex = sizey = sizez = 2600 * u.au
Nx = Ny = Nz = 63 #Number of divisions for each axis
GRID = Model.grid([sizex, sizey, sizez], [Nx, Ny, Nz], rt_code = 'radmc3d')
NPoints = GRID.NPoints #Final number of nodes in the grid

b. Invoke the Model module to compute the physical properties at each GRID node:

#-------------------
#PHYSICAL PROPERTIES
#-------------------
density = Model.density_Constant(r_max, GRID, envDens = dens_e)
temperature = Model.temperature_Constant(density, GRID, envTemp = t_e, backTemp=2.725)

c. Write the data into a file with the RADMC-3D format:

#----------------------
#WRITING RADMC-3D FILES
#----------------------
prop = {'dens_e': density.total,
        'dens_ion': density.total,
        'temp_gas': temperature.total}

Rad = rt.Radmc3dDefaults(GRID)
Rad.freefree(prop)

d. Plot the 3D points distribution:

#------------------------------------
#3D PLOTTING (weighting with density)
#------------------------------------
tag = 'HII'
weight = dens_e
Plot_model.scatter3D(GRID, density.total, weight, NRand = 4000,
                     colordim = density.total / 1e6 / 1e5, axisunit = u.au,
                     cmap = 'winter', marker = 'o',
                     colorlabel = r'$n_{\rm e}$ [cm$^{-3}$]',
                     output = '3Ddens_%s.png'%tag, show = True)

Plot_model.scatter3D(GRID, density.total, weight, NRand = 4000,
                     colordim = temperature.total, axisunit = u.au,
                     cmap = 'winter', marker = 'o',
                     colorlabel = r'$T_{\rm e}$ [Kelvin]',
                     output = '3Dtemp_%s.png'%tag, show = True)

Running RADMC-3D

Making SEDs: In the folder where you stored the sf3dmodels output data files (.inp’s) you should run the following command:

$ radmc3d sed

This command writes the file spectrum.out in your working directory; it has two columns: 1. flux in cgs units and 2. wavelength in microns. Let’s use that information to construct the Spectral Energy Distribution (SED) of the region, at a distance of 4 kpc:

from radmc3dPy.analyze import *
import matplotlib.pyplot as plt

tag = 'ctsphere'

s = readSpectrum(fname = 'spectrum.out') #column 0: wavelength in microns; column 1: Flux in cgs.
distance = 4000. #in pc. The spectrum.out file is still normalized to a distance of 1 pc (see radmc3d docs)
F_nu = s[:,1] * distance**-2 * 1e23 #to Jy at the set distance
nu = 3e8 * s[:,0]**-1 * 1e6 * 1e-9 #microns to GHz
plt.plot(nu, F_nu)
plt.title('%s - distance: %d pc'%(tag,distance))
plt.xlabel('Frequency [GHz]'); plt.ylabel('Flux Density [Jy]')
plt.xscale('log'); plt.yscale('log')
plt.savefig('sed_'+tag+'.png')
plt.show()

Now let’s have a look at the emission of this region at 300 GHz (or 1000 microns). The -simple- command for radmc3d would be:

$ radmc3d image lambda 1000

which writes a file named image.out. The following commands make a simple 2D plot from it:

from radmc3dPy.image import readImage, plotImage
from matplotlib import cm
a=readImage()
plotImage(a,log=True,maxlog=4,cmap=cm.hot,bunit='snu',dpc=4000,arcsec=True) #or au=True

Example 2

Source codes and figures on GitHub: ionized_powerlaw

Note

Power-law spherical HII region

Model: Power-law density distribution and constant temperature

Radiative transfer: free-free emission

Useful references: Keto+2008, Galvan-Madrid+2009

The only difference with respect to the Example 1 is the density model (density_Powerlaw_HII()):

#------------------
#General Parameters
#------------------
#from Galvan-Madrid et al. 2009, Table 3:

MStar = 34 * u.MSun
r_max = 2530 * u.au #H II sphere size
r_min = r_max / 200 #Minimum distance (!= 0 to avoid indeterminations).
r_s = r_max #Normalization distance
rho_s = 1.4e5 * 1e6 #from cgs to SI. Density at r_s
q = 1.3 #Density powerlaw
t_e = 1.e4 #K

#-------------------
#PHYSICAL PROPERTIES
#-------------------
density = Model.density_Powerlaw_HII(r_min, r_max, r_s, rho_s, q, GRID)
temperature = Model.temperature_Constant(density, GRID, envTemp = t_e, backTemp=2.725)

And the resulting plots: