How can I convolve a template spectrum with a photometric filter response spectrum?

  • Suppose I have a template stellar population spectrum (say, from Bruzual & Charlot 2003) which runs from like 1000 Angstroms to 160,000 Angstroms and which has x-axis wavelength units of Angstroms and y-axis flux density units of $L_{\odot}$ per Angstrom. Now, let's say I want to get the integrated luminosity within a bandpass that runs from $\lambda_1$ to $\lambda_2$ Angstroms, and I have another "spectrum" which gives me the filter response function for this photometric bandpass.



    I know that I'm supposed to "convolve the stellar template spectrum with the bandpass response function" in order to get the integrated stellar luminosity within my bandpass, but how do I do this exactly (preferably in python)? The stellar template and filter response spectra don't have the same number of data points and are not defined on the same wavelength grid (the stellar template spectrum obviously extends much further than the filter response function).



    I feel like I could interpolate the filter response spectrum onto the wavelength array of the stellar template spectrum, and assign a value of 0 to the interpolated response spectrum for wavelengths outside $\lambda_1$ to $\lambda_2$. Then I would be able to do a simple element-wise multiplication, but this does not seem like the correct approach using traditional "convolution."


  • pela

    pela Correct answer

    5 years ago

    Your feeling is right: You shouldn't convolve the spectrum and the filter, you should only multiply so that flux outside the bandpass is suppressed. Subsequently you integrate the resulting function over wavelength, so that flux density (in energy/time/area/wavelength) becomes flux (in energy/time/area).



    Simply setting the flux to 0 outside $\lambda_1$ and $\lambda_2$ (or, equivalently, just integrating from $\lambda_1$ to $\lambda_2$) corresponds to a "top-hat filter". Most realistic filters are more smooth.



    So, something like (untested):



    import numpy as np
    from scipy.integrate import simps

    lamS,spec = np.loadtxt('spectrum.dat',unpack=True) #Two columns with wavelength and flux density
    lamF,filt = np.loadtxt('filter.dat' ,unpack=True) #Two columns with wavelength and response in the range [0,1]
    filt_int = np.interp(lamS,lamF,filt) #Interpolate to common wavelength axis
    filtSpec = filt_int * spec #Calculate throughput
    flux = simps(filtSpec,lamS) #Integrate over wavelength
    print 'Total flux is {0:8.3e}'.format(flux)


    Magnitude

    I'm wondering, however, if "the total luminosity within the bandpass" is really what you're interested in. This quantity doesn't really bear any physical significance for the source, as it depends on the particular filter. Usually, one would normalize to the filter, thus getting the magnitude, which is independent of the filter shape.



    In the AB magnitude system, the magnitude $m_\mathrm{AB}$ of a source is (Oke & Gunn 1983; note that they make a sign error in their own definition):
    $$
    m_\mathrm{AB} = -2.5 \log f_\nu -48.6,
    $$
    where the average $f_\nu$ in $\mathrm{erg} \,\mathrm{s}^{-1} \,\mathrm{cm}^{-2} \,\mathrm{Hz}^{-1}$ is given by
    $$
    f_\nu = \frac{1}{c}\frac{\int ST\lambda\,d\lambda}{\int T / \lambda\,d\lambda},
    $$
    where $S$ is the spectrum with units $\mathrm{erg} \,\mathrm{s}^{-1} \,\mathrm{cm}^{-2}$ Å$^{-1}$, $T$ is the filter curve (un-normalized as above), $\lambda$ is the wavelength in Ångström, and $c$ is the speed of light in Å $\mathrm{s}^{-1}$. That is, a flat spectrum (in $\nu$) of height $f_\nu$ would have the same integrated flux over $T$ as $S$ would.



    So:



    import numpy as np
    from scipy.integrate import simps
    c_AAs = 2.99792458e18 # Speed of light in Angstrom/s
    lamS,spec = np.loadtxt('spectrum.dat',unpack=True) #Two columns with wavelength (in Angstrom) and flux density (in erg/s/cm2/AA)
    lamF,filt = np.loadtxt('filter.dat' ,unpack=True) #Two columns with wavelength and response in the range [0,1]
    filt_int = np.interp(lamS,lamF,filt) #Interpolate to common wavelength axis
    I1 = simps(S*T*lam,lam) #Denominator
    I2 = simps( T/lam,lam) #Numerator
    fnu = I1/I2 / c_AAs #Average flux density
    mAB = -2.5*np.log10(fnu) - 48.6 #AB magnitude

    Thanks so much @pela. One quick question: although my filter response function is normalized so that the minimum and maximum values are 0 and 1 respectively, the area under the curve is not itself 1.0 (it's ~0.26). So do I need to divide your filtSpec variable above by the area of the filter response function (i.e., normalize by ~0.26)? For reference, my filter is the 2MASS Ks-band filter provided at the bottom of http://www.ipac.caltech.edu/2mass/releases/allsky/doc/sec6_4a.html

    @quantumflash: The filter shouldn't necessarily go all the way to 1, it just shouldn't go above. A value of 1 (= 100%) means that it transmits _all_ light at that wavelength, which is not necessarily the case. Likewise, the total area under the curve is not of any particular importance. So, if your filter has values all the way from 0 to 1, you should make sure that it is hasn't just been normalized like that, because if so, you have lost the information of the actual throughput.

    @quantumflash: I'm wondering if "the total luminosity within the bandpass" is really what you're interested in. This quantity doesn't really bear any physical significance for the source, as it depends on the particular filter. Usually, one would normalize to the filter, thus getting the _magnitude_, which is independent of the filter shape. If this is actually what you want, I can edit my answer.

    I have a stellar population template spectrum from Bruzual & Charlot 2003 which has y-axis units of $L_{\odot}$ per Angstrom, and the total spectrum is normalized to correspond to luminosity per wavelength output by 1 $M_{\odot}$. If I can integrate this spectrum in the 2MASS Ks-band given the filter from the 2MASS website above, then I can compute the mass-to-light ratio in the K-band for my synthetic (old) stellar population, which is what I want. I don't think I need the magnitude since the y-axis of the spectrum is already in $L_{\odot}$ per Angstrom. Thank you for your help!

    @pela, it's not quite true that the magnitude is independent of the filter right? Rather it's defined with respect to the filter.

    @John: You're right that, in general, it's dependent on the filter in the sense that in one wavelength range an object may be one magnitude, while in another range it may be another. Various magnitude systems can then be defined so that the object _does_ have the same magnitude (e.g. in the AB system, a flat spectrum (in $\nu$) is the same, while in the Vega system, Vega has magnitude 0 in all bands). [cont'd below]

    What I meant was that the exact filter _shape_ does not influence the result, since this is taken care of in the integral. Measure the magnitude in a given range with a tophat filter, a triangular filter, and a Gaussian filter, and you get the same result.

    @Pela: That's only true for a flat spectrum though isn't it? Otherwise, even though two filters may have the same wavelength range, their response curves may give different results I reckon.

    @John: You're right, it's only completely true for a truly flat spectrum. However, most realistic spectra are flat enough over the range of a filter, and most realistic filters enough alike, that in reality the exact filter shape doesn't matter. Astronomers are not mathematicians, so they're generally happy with this. But you're absolutely right that it's not exactly true.

    @pela I'm returning to this and I was wondering if you can help clear some confusion (you can skip the stuff above). Theoretical templates like Bruzual & Charlot output spectra with luminosity density units (Lsun/AA) vs AA. But to compute the AB mag within say the SDSS r-band filter, I need the y-axis units to be in natural flux density units (erg/s/cm**2/AA). What do people usually do for this luminosity density to flux density conversion for theoretical templates? I'm wondering if I can use the normal F=L/(4piD^2) formula with D=10 pc, since abs mags and luminosities assume D=10 pc? Thanks!

    @quantumflash You're right wrt. the formula, but you should only use d = 10 pc if the distance to the observed object is 10 pc. I guess you're thinking of converting between apparent and absolute magnitudes. The purpose of a theoretical template spectrum is to be able to compare it to an observed one, and the observed object will be at some distance, so you put the theoretical one at the same distance.

    Thanks @pela -- I'm just doing "mock tests" with simple stellar populations. Often these are given in units of Lsun/AA vs AA and normalized to 1 Msun of stars formed (so that you can easily compute things like stellar mass-to-light ratio in some bandpass). But the problem now is that I want to compare two different model spectra and compute the magnitude difference in some bandpasses. I guess in principle if I only care about $\Delta$ mag difference, I can use any $D$ as long as its the same for the two individual mags since ultimately it'll cancel when I subtract $m_1-m_2$ or (abs) $M_1-M_2$

    @pela yes I think I got this worked out, thanks a lot again!

    @quantumflash Yes exactly, to compare two different templates, you may assume any distance (and 10 pc is a good choice).

License under CC-BY-SA with attribution


Content dated before 7/24/2021 11:53 AM