Divide by a black body

Often in real-world astronomical applications, we want to see how a stellar spectrum varies in wavelength without the large and omnipresent crest of the black body curve.

In this tutorial we show how to remove the black body curve from a precomputed model spectrum.

from gollum.phoenix import PHOENIXSpectrum
from gollum.precomputed_spectrum import PrecomputedSpectrum
from astropy.modeling.physical_models import BlackBody
import astropy.units as u
import numpy as np
%config InlineBackend.figure_format='retina'

T_eff = 5_700

original_spectrum = PHOENIXSpectrum(teff=T_eff, logg=4.5, metallicity=0,
                                    wl_lo=0, wl_hi=1e10) # Get the whole spectrum

The PHOENIX spectra have units of \(\frac{\mathrm{erg}}{\mathrm{s\;cm^2 \;cm}}\)

original_spectrum.flux.unit

$\mathrm{\frac{erg}{s\,cm^{3}}}$

original_spectrum.flux.min()

$5.5922428 \times 10^{-15} \; \mathrm{\frac{erg}{s\,cm^{3}}}$

Make a black body spectrum with the same temperature

blackbody_model = BlackBody(temperature=T_eff*u.Kelvin)
blackbody_flux_per_sr = blackbody_model(original_spectrum.wavelength)

The Black Body spectra have units of \(\frac{\mathrm{erg}}{\mathrm{s\;cm^2 \;Hz \;sr}}\). To convert between the two conventions we have to do two things. First we have to pretend we are standing on the surface of the star, and multiply by \(\pi\) steradians. Second, we have to convert flux density in per frequency to per wavelength units. We can do that with astropy’s equivalencies kwarg.

blackbody_flux_per_Hz = blackbody_flux_per_sr * np.pi * u.steradian

blackbody_flux = blackbody_flux_per_Hz.to(original_spectrum.flux.unit,
                  equivalencies=u.spectral_density(original_spectrum.wavelength))

blackbody = PrecomputedSpectrum(flux=blackbody_flux,
                                spectral_axis=original_spectrum.wavelength)

ax = original_spectrum.plot()
blackbody.plot(ax=ax)
ax.set_ylim(1e0, 1e16)
ax.set_yscale('log')
ax.set_ylabel('Flux (erg/s/cm$^3$)');

The plot spans 16 orders of magnitude– a huge dynamic range! Notice that the spectra have similar—but not identical—broadband spectral shapes. They should have the identical area under the curve, by the definition of effective temperature. Let’s see if they do!

from scipy.integrate import trapezoid

original_flux = trapezoid(original_spectrum.flux, x=original_spectrum.wavelength.to(u.cm))
black_body_flux = trapezoid(blackbody.flux, x=original_spectrum.wavelength.to(u.cm))

original_flux/black_body_flux

$0.99048532 \; \mathrm{}$

The two spectral models have the same flux to within 1%, which is close enough to identical given that more spectrum resides outside the extent of the plot, and numerical artifacts limitations in the spectral modeling procedure. Let’s compute the ratio spectrum to see how flat the spectrum looks. We’ll first plot it over the same dynamic range as before to emphasize how much more compressed it is.

ratio_spectrum = original_spectrum.divide(blackbody)

The resulting spectrum is a ratio of fluxes with the same units, so it is dimensionless.

ratio_spectrum.flux.unit == u.dimensionless_unscaled

True

ax = ratio_spectrum.plot();

ax.set_ylim(1e-12, 1e4)
ax.set_yscale('log')
ax.set_ylabel('Normalized flux');

ax.axhline(1.0, linestyle='dashed', color='k');

At this dramatic zoom level, the flux looks pretty flat (except for the extreme ultraviolet portion of the spectrum). Let’s zoom in on a region of interest from \(8000 - 13000\; Å\).

ax = ratio_spectrum.plot(ylo=0, yhi=1.15);
ax.set_ylabel('Normalized flux');

ax.axhline(1.0, linestyle='dashed', color='k');
ax.set_xlim(8_000, 13_000);

OK, looks good! We have successfully used the blackbody curve to coarsely flatten the spectrum!