# Emission Line Galaxy Classification#

Chris Richardson, Elon University

Funding for this tutorial made possible through the CATL Scholars Program

This tutorial is intended for students to complete in groups. A version with solutions is available to instructors upon request.

## Objectives#

• To explore classifying emission line galaxies according to the radiation source excitating gaseous clouds within each galaxy

• To understand how emission line ratios can help constrain excitation mechanisms

• To explore how models can highlight the strengths and weakness of various diagnostic diagrams

• To provide a foundation to explore your own questions about excitation mechanisms in galaxies

## Learning Outcomes#

By the end of this tutorial you should be able to:

• Use emission line ratios to explore their diagnostic potential for finding AGN

• Use diagnostic diagrams to demarcate AGN and SF boundaries

• Use photoionization models to test the validity of diagnostic diagrams

## Prerequisites#

• Computational

• Basic Python skills (e.g., for loops, conditional statements, defining functions)

• Basics of matplotlib visualization, SQL queries, and SciServer platform (e.g., see here)

• Using the astropy package for fundamental constants and unit conversion

• Astrophysics

• Introductory knowledge of the formation of spectral lines

• Basic understanding of the physical picture present in star forming galaxies and active galaxies (i.e., AGN)

• Advanced conceptual understanding of using emission lines to probe spectral energy distributions

# Import relevant libraries

import SciServer.CasJobs as CasJobs # query with CasJobs
import SciServer.SkyServer as SkyServer # look up objects and generate images
import pandas as pd # data analysis
import numpy as np # more data analysis
import seaborn as sns # plotting library
from matplotlib import pyplot as plt # another plotting library

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
1 # Import relevant libraries
----> 3 import SciServer.CasJobs as CasJobs # query with CasJobs
4 import SciServer.SkyServer as SkyServer # look up objects and generate images
5 import pandas as pd # data analysis

ModuleNotFoundError: No module named 'SciServer'

# Make notebook have a clean appearance

pd.set_option('display.max_colwidth', None)
import warnings
warnings.filterwarnings('ignore')


# Choosing a sample of galaxies#

GalSpecLine and emissionLinesPort are examples of two schema available to search for emission line galaxies. The most common emission lines are featured in both schema, but there are several lines not cross listed between the two catalogs. The methods used for measuring the emission line strengths differ between the catalogs, each with their unique merits. We will use GalSpecLine for this tutorial, however emissionLinesPort could be used in a follow up project.

While you can include as many or as little emission lines as you’d like in your query, you should always include the emission line flux errors and the galaxy redshifts. The errors on the emission line measurements enable you to improve the quality of your sample if only specific emission lines interest you. It is better to impose these restrictions after your query as it can dramatically reduce the sample size. The redshift of each galaxy is important because certain emission lines will shift into and outside the SDSS window of observation.

• At what $$z$$ will the shortest wavelength emission line shift into the observation window?

• At what $$z$$ will the longest wavelength emission line shift outside the observation window?

The SDSS aperture is fixed in size and is typically aligned with the nuclear center of each galaxy.

• How might the size of the SDSS aperture tie into your ability to discern galaxy wide properties at certain redshifts?

• What sort of bias does this introduce into your results?

# Create a query to generate a sample of galaxies

query = 'SELECT TOP 10000 s.plate, s.fiberid, s.mjd, s.z, s.zwarning, '
query += 'g.oii_3726_flux, g.oii_3729_flux, ' #-- Gather possible emission lines of interest
query += 'g.neiii_3869_flux, g.h_delta_flux, '
query += 'g.h_gamma_flux, g.oiii_4363_flux, '
query += 'g.h_beta_flux, g.oiii_4959_flux, '
query += 'g.oiii_5007_flux, g.hei_5876_flux, '
query += 'g.oi_6300_flux, g.nii_6548_flux, '
query += 'g.h_alpha_flux, g.nii_6584_flux, '
query += 'g.sii_6717_flux, g.sii_6731_flux, '
query += 'g.ariii7135_flux, '
query += 'g.oii_3726_flux_err, g.oii_3729_flux_err, ' #-- Gather emission line errors
query += 'g.neiii_3869_flux_err, g.h_delta_flux_err, '
query += 'g.h_gamma_flux_err, g.oiii_4363_flux_err, '
query += 'g.h_beta_flux_err, g.oiii_4959_flux_err, '
query += 'g.oiii_5007_flux_err, g.hei_5876_flux_err, '
query += 'g.oi_6300_flux_err, g.nii_6548_flux_err, '
query += 'g.h_alpha_flux_err, g.nii_6584_flux_err, '
query += 'g.sii_6717_flux_err, g.sii_6731_flux_err, '
query += 'g.ariii7135_flux_err, '
query += 'e.bptclass '
query += 'FROM GalSpecLine AS g, SpecObj AS s, GalSpecExtra AS e '
query += 'WHERE '
query += 's.specobjid = g.specobjid '
query += 'AND s.specobjid = e.specobjid '
query += 'AND z BETWEEN 0.04 and 0.1 '#-- Lower limit needed for global parameters w/ small aperture
query += 'AND h_alpha_flux > h_alpha_flux_err*5 ' #-- Make sure Balmer lines have adequate S/N
query += 'AND h_beta_flux > h_beta_flux_err*5 '
query += 'AND h_alpha_flux_err > 0 ' #-- Consistency check
query += 'AND h_beta_flux_err > 0 '
query += 'AND sigma_balmer < 509.55 ' #-- Exclude the broad line emitting galaxies
query += 'AND s.class = \'GALAXY\' '  #-- Looking for galaxies, not stars
query += 'AND s.zwarning = 0 '
#print (query) # useful for debugging, remove first # to uncomment

# send query to CasJobs
gals = CasJobs.executeQuery(query, "dr16")

gals

plate fiberid mjd z zwarning oii_3726_flux oii_3729_flux neiii_3869_flux h_delta_flux h_gamma_flux ... oiii_5007_flux_err hei_5876_flux_err oi_6300_flux_err nii_6548_flux_err h_alpha_flux_err nii_6584_flux_err sii_6717_flux_err sii_6731_flux_err ariii7135_flux_err bptclass
0 266 4 51630 0.064656 0 39.057380 27.096570 2.468194 6.669699 12.437470 ... 3.302969 2.745098 2.520172 1.046200 3.448143 3.155578 2.919328 2.860605 4.992754 3
1 266 120 51630 0.069341 0 -1.330543 7.919315 -0.333586 -1.529107 10.477300 ... 2.545394 1.778672 1.680469 0.703648 2.371513 2.122362 2.055846 1.841887 3.126379 2
2 266 146 51630 0.052176 0 145.169700 184.163200 26.046470 67.026570 114.044300 ... 4.271465 3.850465 3.498755 2.040700 10.994270 6.155217 4.846725 4.428813 3.586765 1
3 266 533 51630 0.085351 0 31.840760 28.041630 2.403825 12.097370 18.920390 ... 2.793131 1.910691 1.961116 0.822471 3.265500 2.480760 2.806020 2.475492 2.315306 2
4 266 630 51630 0.073687 0 12.878990 21.359510 3.436489 0.142290 5.390950 ... 1.988206 1.574178 1.350679 0.568662 2.131675 1.715214 1.609602 1.671389 1.910809 2
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
9995 1207 463 52672 0.055911 0 -0.273259 21.205140 0.332328 6.859152 4.029799 ... 2.268165 2.104548 1.901967 0.799355 2.466196 2.411038 2.136849 2.025468 2.065840 2
9996 1207 512 52672 0.040428 0 168.423800 208.261600 16.263230 33.077510 53.475750 ... 3.645575 2.059486 1.995668 0.826124 5.905511 2.491776 2.573609 2.407138 1.849287 1
9997 1208 202 52672 0.051229 0 75.635860 81.785600 4.387515 17.040500 28.414330 ... 3.219081 2.728835 2.423338 1.162954 4.954031 3.507735 3.237844 3.007984 2.462693 1
9998 1208 490 52672 0.093831 0 37.472160 44.705480 2.661830 16.227150 30.224240 ... 2.448881 1.810508 1.965327 0.843513 4.036311 2.544227 2.512585 2.451001 2.777650 1
9999 1208 593 52672 0.071304 0 71.641330 85.087980 5.179615 23.398060 36.458060 ... 2.943066 2.206493 2.083630 1.008152 4.212832 3.040816 2.661236 2.490553 2.830445 1

10000 rows × 40 columns

# Using emission lines to find active galaxies#

In the lecture-tutorial on emission line ratio diagnostics, you discovered how to find emission line ratios sensitive to the excitation source. In summary, ratios featuring emission lines with relatively different ionization potentials effectively probe the shape of the ionizing radiation field and therefore the excitation source (e.g., starlight, AGN).

As a trial run, let’s explore how well the emission line ratio [O II] $$\lambda$$3727 / [Ne III] $$\lambda$$3869 separates star forming (SF) galaxies and AGN. N.B. [O II] $$\lambda$$3727 is actually the sum of the [O II] $$\lambda$$3726 and [O II] $$\lambda$$3729 doublet.

• What are relavent ionization potentials for [O II] and [Ne III] line emission?

• Would you expect AGN to display higher or lower [O II]/[Ne III] ratios compared to SF galaxies?

# Remove possible +/- inf values
gals_flux = gals[(gals['neiii_3869_flux']>0) & (gals['oii_3726_flux']>0) & (gals['oii_3729_flux']>0)
& (gals['oi_6300_flux']>0) & (gals['oiii_5007_flux']>0)]

# Make sure relavent emission lines have well measured fluxes
quality_gals = gals_flux[(gals_flux['neiii_3869_flux']>5.0*gals_flux['neiii_3869_flux_err']) &
(gals_flux['oii_3726_flux']>5.0*gals_flux['oii_3726_flux_err']) &
(gals_flux['oii_3729_flux']>5.0*gals_flux['oii_3729_flux_err']) &
(gals_flux['oi_6300_flux']>5.0*gals_flux['oi_6300_flux_err'])]

# Split the sample into AGN and SF galaxies
gals_sf = quality_gals[quality_gals['bptclass']==1]
gals_agn = quality_gals[(quality_gals['bptclass']==4)]

oii_doublet_sf = gals_sf['oii_3726_flux']+gals_sf['oii_3729_flux']
oii_doublet_agn = gals_agn['oii_3726_flux']+gals_agn['oii_3729_flux']

oii_neiii_sf = np.log10(oii_doublet_sf/gals_sf['neiii_3869_flux'])
oii_neiii_agn = np.log10(oii_doublet_agn/gals_agn['neiii_3869_flux'])

plt.hist(oii_neiii_sf,color='blue',alpha=0.5,label='SF')
plt.hist(oii_neiii_agn,color='red',alpha=0.5,label='AGN')
plt.xlim(0.0,2.0)
#plt.title('[O II]/[Ne III] as a function of galaxy classification')
plt.xlabel('[O II]/[Ne III]')
plt.ylabel('N')
plt.legend()
plt.show()


The histogram shows the number of galaxies in each bin as a function of galaxy classification. It is useful in telling us that the $$N_{SF}$$ > $$N_{AGN}$$ and that each distribution peaks at a differnt [O II]/[Ne III]. But, the actual location of those peaks is bit uncertain because the histogram is rather blocky.

The Python package matplotlib will try to select a bin width that results in a visually appealing histogram. We can adjust the number of bins, and therefore the bin width, using the “bins” kwarg. However, histograms will always appear bumpy and depend upon the endpoints of the bins. One way around this is to use a Keneral Density Estimator (KDE) plot, which you can read about here. We can use the “seaborn” package to make publicaton quality KDE plots that show a smooth distribution of our data. The plot now shows the density of galaxies rather than the number of galaxies.

sns.kdeplot(oii_neiii_sf,color='Blue',fill=True,label='SF')
sns.kdeplot(oii_neiii_agn,color='Red',fill=True,label='AGN')
plt.xlim(0.0,2.0)
plt.xlabel('[O II]/[O III]')
plt.ylabel('Density')
plt.legend()
plt.show()


The KDE plot of [O II]/[Ne III] shows that no matter how we slice it, this ratio alone is not capable of separating SF and AGN excitation in galaxies. Let’s assess the ability of [O III] $$\lambda$$5007 /[O I] $$\lambda$$6300 to act as an effective excitation diagnostic. Try generating a KDE plot for this emission line ratio in the space below and ask yourself the following questions.

• What are relavent ionization potentials for [O III] and [O I] line emission?

• Would you expect AGN to display higher or lower [O III]/[O I] ratios compared to SF galaxies?

• Can you make a single cut to separate all galaxies with an AGN from those only experiencing SF?

# your code for KDE plot
#oiii_oi_sf =
#oiii_oi_agn =


# Diagnostic Diagrams#

We can combine multiple line ratios in a diagram to increase their utility in separating AGN and SF activity. This assesses the shape of the inicident radiation field (starlight or AGN) at more energies than possible with a single emission line ratio. Let’s combine [O II]/[Ne III] and [O III]/[O I] into a single plot.

plt.scatter(oii_neiii_sf,oiii_oi_sf,color='blue',label='SF',s=0.5,alpha=0.75)
plt.scatter(oii_neiii_agn,oiii_oi_agn,color='red',label='AGN',s=0.5,alpha=0.75)

plt.xlabel('[O II]/[Ne III]')
plt.ylabel(r'[O III]/[O I]')

plt.xlim(0.2,1.8)
plt.ylim(0.3,2.0)

plt.legend()

plt.show()

NameErrorTraceback (most recent call last)
<ipython-input-18-cf6c5fbac493> in <module>
----> 1 plt.scatter(oii_neiii_sf,oiii_oi_sf,color='blue',label='SF',s=0.5,alpha=0.75)
2 plt.scatter(oii_neiii_agn,oiii_oi_agn,color='red',label='AGN',s=0.5,alpha=0.75)
3
4 plt.xlabel('[O II]/[Ne III]')
5 plt.ylabel(r'[O III]/[O I]')

NameError: name 'oiii_oi_sf' is not defined


The [O III]/[O I] vs. [O II]/[Ne III] diagram does a reasonable job in separating AGN and SF galaxies. We can make this separation more explicit by adding a demarcation line. Since we might want to generate this plot over and over, we can make define a function that will do it for us.

def make_diagram():

o2ne3_demarc = np.linspace(0.2, 1.19)
o3o1_demarc = (0.055 / (o2ne3_demarc - 1.19)) + 1.35

plt.scatter(oii_neiii_sf,oiii_oi_sf,color='blue',label='SF',s=0.5,alpha=0.75)
plt.scatter(oii_neiii_agn,oiii_oi_agn,color='red',label='AGN',s=0.5,alpha=0.75)

plt.plot(o2ne3_demarc,o3o1_demarc,'k')

plt.xlabel('[O II]/[Ne III]')
plt.ylabel(r'[O III]/[O I]')

plt.xlim(0.2,1.8)
plt.ylim(0.3,2.0)

plt.legend()

######

make_diagram()
plt.show()


# Justifying diagnostic diagrams with models#

Up until this point, we have been using observed galaxies flagged as an AGN or star forming galaxy in the GalSpecExtra. While we would like to take this information for granted, it is important to note these tags are based on a set of assumptions that may or may not hold for all galaxies. We can use photoionization models to understand the limits of these AGN/SF distinctions and our newly created diagnostics diagram.

Photoionization models predict the observed emission line spectrum for gas clouds illuminated by an AGN or starlight in a galaxy. The set that we will use comes from this paper and a link to the simulated emission region data has been provided on Moodle. To make the visualization easier, we can restrict the model parameter space to regions that we’re interested in exploring. A brief description of the parameters is given below along with accepted values within the grid.

LOGU - dimensionaless quantity (ionization parameter) specifying the ionization level of the gas clouds
options: -4.0 to -0.5 in 0.25 increments

LOGZ - metallicity of the gas clouds relative to solar metallicity
options: [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.7, 1.0, 1.5, 2.0]

AGNFRAC - the fraction of AGN excitation as opposed to starlight excitation
options: [0.0, 0.04, 0.08, 0.16, 0.32, 0.64, 1.0]

Let’s see how the emission line predictions compare to our observed galaxies on the diagnostic diagram when we select a simulation with 100% AGN excitation and Z = 0.2 Z$$_{\odot}$$.

!wget --no-check-certificate 'https://drive.google.com/uc?export=download&id=1vQrGvdOf2w9VAu7F0lub5tEAi75cIMm4' -O model_file

#print (model_grid.columns) # useful to check out the column names in the file

# function to restrict the model parameter space
def apply_limits(grid,f_agn,U,Z):

grid = grid[(grid['AGNFRAC'] == f_agn) & (round(grid['LOGZ'],5) == round(np.log10(Z),5)) \
& (grid['MIXING'] == ' coincident') & (grid['LOGU'] >= U[0]) & (grid['LOGU'] <= U[1])]

return grid

U = [-3.75,-2.75] # range of ionization parameter
f_agn = 1.0
Z = 0.2

agn_models = apply_limits(model_grid,f_agn,U,Z)
#print (agn_models)

model_x = np.log10((agn_models[' \'O__2_372603A\'']+agn_models[' \'O__2_372881A\''])/agn_models[' \'NE_3_386876A\''])
model_y = np.log10(agn_models[' \'O__3_500684A\'']/agn_models[' \'O__1_630030A\''])

make_diagram()

plt.plot(model_x,model_y,color='Green',marker='o')

plt.show()


The AGN model assuming 0.2 Z$$_{\odot}$$ lies in the AGN part of our diagram and including a range of ionization parameter covers some of the variation in the region. It seems that the models preliminarily suggest the galaxies flagged as AGN really are AGN and the diagnostic diagram can pick them out. In the cells below, determine if this is a robust result by answering the following questions.

• Do all of the 100% AGN simulations fall within the AGN part of the diagnostic diagram?

• If so, what range of ionization parameters and metallicities do the observed AGN span? If not, what values disagree with our chosen demarcations?

• Do all of the purely SF simulations fall within the SF part of the diagnostic diagram?

• If so, what range of ionization parameters and metallicities do the observed SF galaxies span? If not, what values disagree with our chosen demarcations?

# your code to answer the questions above


# The BPT Diagram#

The diagnostic diagram that we have developed does have it’s limitations. For example, both of the [O II] and [Ne III] emission lines are unavailable for low-$$z$$ galaxies because they fall outside the SDSS observing window. This means that we’re unable to identify AGN activity for local galaxies. Another limitation is that some of the emission lines can be rather weak, making a reliable detection difficult. We can quantify this by looking at much our sample reduces down after applying our S/N cut.

print ('Fraction of initial SF galaxies included after S/N cut: ',
"{:.3}".format(len(gals_sf)/len(gals[gals['bptclass']==1])))
print ('Fraction of initial AGN galaxies included after S/N cut: ',
"{:.3}".format(len(gals_agn)/len(gals[gals['bptclass']==4])))


Therein lies the problem: most galaxies only have a handful of reliable emission lines measurements. This is one reason why the diagnostic diagram we have been analyzing isn’t a popular choice in research literature even though it does a reasonable job of distinguishing between AGN and star forming galaxies.

The so-called “gold standard” for classifying AGN with optical spectroscopy is known as the “BPT” diagram, named after the authors that developed it in this paper. The BPT diagram uses emission line ratios from [O III] $$\lambda$$5007, [N II] $$\lambda$$6584, H$$\beta$$, and H$$\alpha$$, all of which are known as “strong lines” since they’re strong enough to be reliably measured in most galaxies. Using strong lines for diagnostics drastically improves the number of galaxies we can include in a sample.

gals_flux = gals[(gals['oiii_5007_flux']>0) & (gals['nii_6584_flux']>0)]

quality_gals = gals_flux[(gals_flux['oiii_5007_flux']>5.0*gals_flux['oiii_5007_flux_err']) &
(gals_flux['nii_6584_flux']>5.0*gals_flux['nii_6584_flux_err'])]

gals_sf = quality_gals[quality_gals['bptclass']==1 | (quality_gals['bptclass']==2)]
gals_agn = quality_gals[(quality_gals['bptclass']==4)]

print ('Fraction of initial SF galaxies included after S/N cut: ',
"{:.3}".format(len(gals_sf)/len(gals[gals['bptclass']==1])))
print ('Fraction of initial AGN galaxies included after S/N cut: ',
"{:.3}".format(len(gals_agn)/len(gals[gals['bptclass']==4])))


The BPT diagram remains a hot topic in academic literature today. In fact, the GalSpecExtra schema actually uses the BPT diagram to classify galaxies as AGN or star forming. If we use the standard demarcations from literature, it is no surprise in the plot the below that all of the AGN and star forming galaxies are perfectly separated. There is even an additional class of galaxies known as “composites” that are thought to represent a transition state between purly star forming galaxies and those excited purely by an AGN.

def make_bpt_plot():

gals_comp = quality_gals[(quality_gals['bptclass']==3)]

oiii_Hb_sf = np.log10(gals_sf['oiii_5007_flux']/gals_sf['h_beta_flux'])
oiii_Hb_agn = np.log10(gals_agn['oiii_5007_flux']/gals_agn['h_beta_flux'])
oiii_Hb_comp = np.log10(gals_comp['oiii_5007_flux']/gals_comp['h_beta_flux'])

nii_Ha_sf = np.log10(gals_sf['nii_6584_flux']/gals_sf['h_alpha_flux'])
nii_Ha_agn = np.log10(gals_agn['nii_6584_flux']/gals_agn['h_alpha_flux'])
nii_Ha_comp = np.log10(gals_comp['nii_6584_flux']/gals_comp['h_alpha_flux'])

plt.scatter(nii_Ha_sf,oiii_Hb_sf,color='blue',label='SF',s=0.5,alpha=0.75)
plt.scatter(nii_Ha_agn,oiii_Hb_agn,color='red',label='AGN',s=0.5,alpha=0.75)
plt.scatter(nii_Ha_comp,oiii_Hb_comp,color='magenta',label='Comp',s=0.5,alpha=0.75)

comp_n2ha = np.linspace(-2.5,-0.1,100)
comp_o3hb = 0.61/(comp_n2ha-0.05)+1.3
plt.plot(comp_n2ha,comp_o3hb,color='k',linestyle='--',linewidth=1.0)

agn_n2ha = np.linspace(-2.5,0.3,100)
agn_o2hb = 0.61/(agn_n2ha-0.47)+1.19
plt.plot(agn_n2ha,agn_o2hb,color='k',linewidth=1.0)

plt.xlabel(r'[N II]/H$\alpha$')
plt.ylabel(r'[O III]/H$\beta$')

plt.xlim(-1.5,0.5)
plt.ylim(-1.0,1.25)

plt.legend()

######

make_bpt_plot()
plt.show()


All of the strong lines on the BPT diagram fall within the SDSS observation window for both local galaxies and those at higher redshifts. Combined with the superior statistics shown above, this makes the BPT diagram a suitable choice for tracing galaxy evolution. We can return to the photoionization models to investigate how a purely star forming model compares to the selected observations.

make_bpt_plot()

U = [-3.75,-2.75]
f_agn = 0.0
Z = 0.4

agn_models = apply_limits(model_grid,f_agn,U,Z)

#print (model_grid.columns) # useful to check out the column names in the file

model_x = np.log10(agn_models[' \'N__2_658345A\'']/agn_models[' \'H__1_656281A\''])
model_y = np.log10(agn_models[' \'O__3_500684A\'']/agn_models[' \'H__1_486133A\''])

plt.plot(model_x,model_y,color='Cyan',marker='o')

plt.show()


All diagnostic diagrams have their limitations and it is important to keep the limitations in mind when interpreting results. In the space below, use the photionization models to explore the limitations of the BPT diagram in the same way that you previously explored the limitations of the [O III]/[O I] vs. [O II]/[Ne III] diagram.

# your code to compare photoionization models to selected observations


• What are the effects of selecting models with non-binary $$f_{AGN}$$?