Calibrating spectrum to photometry

[tjtakys]

Hi, I am currently using Prospector to fit both photometry and spectrum simultaneously for several high-z galaxies. The spectrum I am analyzing comes from NIRSpec with a medium grating. While individual channels in the spectrum are noisy and do not show a clear continuum (left figure shows an example spectrum), summing across many channels makes the continuum detection. However, when convolving the spectrum with filter transmissions and comparing it to photometry, I observe up to a 50% discrepancy (right figure shows the relative difference for several galaxies).

To address these discrepancies, I use optimize_speccal with a second-order polynomial for spectral calibration. However, in some cases, the spectrum appears to be over-calibrated, which may strongly influence the derived parameters, such as SFH, Av, and gas logZ. This is a concern for me.

model_params = {}
# common parameters
model_params["imf_type"] = {'N':1, 'isfree':False, 'init':1} # 0:Salpeter, 1:Chabrier, 2:Kroupa
model_params["logzsol"] = {'N':1, 'isfree':True, 'init':-0.5, 'units':'log (Z/Zsun)', 'prior':priors.TopHat(mini=-2, maxi=0.2)}

# Non-parametric SFH
model_params.update(TemplateLibrary["continuity_sfh"])
nbin = 6
agelims = [0.0, 7.0] + np.linspace(start=7, stop=np.log10((cosmo.age(z=zspec)-cosmo.age(z=20)).to('yr').value), num=nbin)[1:].tolist()
agebins = np.array([agelims[:-1], agelims[1:]]).T

# This is the *total*  mass formed, as a variable
model_params["logmass"] = {"N":1, "isfree":True, "init":10, 'units':'Msun', 'prior':priors.TopHat(mini=8, maxi=12)}
# This will be the mass in each bin.  It depends on other free and fixed parameters.  Its length needs to be modified based on the number of bins
# This converts from an array of log_10(SFR_j / SFR_{j+1}) and a value of log10(\Sum_i M_i) to values of M_i.  j=0 is the most recent bin in lookback time.
model_params["mass"] = {'N':nbin, 'isfree':False, 'init':1e6, 'units':r'M$_\odot$', 'depends_on':transforms.logsfr_ratios_to_masses}
# This gives the start and stop of each age bin.  It can be adjusted and its length must match the lenth of "mass"
model_params["agebins"] = {'N':nbin, 'isfree':False, 'init':agebins, 'units':'log(yr)'}
# This controls the distribution of SFR(t) / SFR(t+dt). It has NBINS-1 components.
model_params["logsfr_ratios"] = {'N':nbin-1, 'isfree':True, 'init':np.full(nbin-1, 0.0),
                                    'prior':priors.StudentT(mean=np.full(nbin-1, 0.0),
                                                            scale=np.full(nbin-1, 0.3),
                                                            df=np.full(nbin-1, 2))}

# redshift
model_params["zred"] = {'N':1, 'isfree':False, 'init':zspec, 'units':'redshift'}

# Dust attenuation (diffuse: Kriek & Conroy 2013 + birth cloud: index=-1)
model_params["dust_type"] = {'N':1, 'isfree':False, 'init':4, 'units':'FSPS index'} # Kriek & Conroy 2013
model_params["dust_tesc"] = {'N':1, 'isfree':False, 'init':7.0, 'units':'log(yr)'} 
model_params["dust1"] = {'N':1, 'isfree':False, 'init':0.0, 'units':'optical depth towards young stars', 'depends_on':transforms.dustratio_to_dust1}
model_params["dust2"] = {'N':1, 'isfree':True, 'init':0.2, 'units':'optical depth for all stars at 5500AA', 'prior':priors.TopHat(mini=0.0, maxi=4.0)}
model_params["dust_ratio"] = {'N':1, 'isfree':True, 'init':1.0, 'units':'ratio of birth-cloud to diffuse dust', 'prior':priors.ClippedNormal(mini=0.0, maxi=2.0, mean=1.0, sigma=0.3)}
model_params["dust_index"] = {'N':1, 'isfree':False, 'init':-0.7, 'units':'power-law multiplication of Calzetti', 'prior':priors.TopHat(mini=-1.0, maxi=0.4)}
model_params["dust1_index"] = {'N':1, 'isfree':False, 'init':-1.3, 'units':'power-law multiplication of Calzetti'}

# Nebular
model_params.update(TemplateLibrary["nebular"]) 
model_params["nebemlineinspec"] = {'N':1, 'isfree':False, 'init':False} 
model_params['eline_sigma'] = {'N': 1, 'isfree':True, 'init':150.0, 'units':'km/s', 'prior':priors.TopHat(mini=50, maxi=300)}
model_params["gas_logz"] = {'N':1, 'isfree':True, 'init':0.0, 'units':'log Z/Zsun', 'prior':priors.TopHat(mini=-2, maxi=0.5)}
model_params["gas_logu"] = {'N':1, 'isfree':True, 'init':-2.0, 'units':'Q_H/N_H', 'prior':priors.TopHat(mini=-4, maxi=-1)}
    
# smoothing
model_params.update(TemplateLibrary["spectral_smoothing"])
model_params['fftsmooth'] = {'N': 1, 'isfree': False, 'init': True}
model_params['sigma_smooth'] = {'N':1, 'isfree':True, 'init':200, 'units':'km/s', 'prior':priors.TopHat(mini=50, maxi=200)} # sigma, NOT FWHM

# continuum optimize
model_params.update(TemplateLibrary["optimize_speccal"])
model_params['polyorder'] = {'N': 1, 'isfree': False, 'init': 2} # 12 is too much?
model_params['spec_norm'] = {'N': 1, 'isfree': False, 'init': 1.0, 'units': 'f_true/f_obs', 'prior':priors.Normal(mean=1.0, sigma=0.1)} 

model = PolySpecModel(model_params)

Below, I include SED plots for two galaxies. In each plot:
• Gray line: The intrinsic best-fit spectrum (without spec cal).
• Blue line: The full spectrum after spec cal.

# full sed 
obs_copy = obs.copy()
obs_copy['spectrum'] = None
obs_copy['wavelength'] = None
model_wl = sps.wavelengths * (1 + zspec) / 1e4
model_fullspec, _, _ = model.predict(result["bestfit"]["parameter"], obs_copy, sps=sps)
# nebular line
obs["wavelength"] = model._outwave
obs["spectrum"] = np.ones_like(model._outwave)
model.cache_eline_parameters(obs)
spec_neb = model.predict_eline_spec(line_indices=model._use_eline, wave=model._outwave).sum(axis=1)

ax.step(model_wl, model_fullspec+spec_neb, lw=0.5, color='gray', label="Model spec (full)")
ax.step(wspec, result["bestfit"]["spectrum"], lw=0.5, color='tab:orange', label="Best spec") # intrinsic spectrum

The first case shows minimal differences between the intrinsic and calibrated spectra, while the second case exhibits significant wavelength-dependent adjustments. Since there is no strong wavelength-dependent offset in the actual photometry-spectra comparison, I believe this large adjustment is incorrect. I am concerned that this over-calibration is skewing the derived physical parameters.

Questions:

I am considering the following approaches to address this issue:

1. Use even lower polynomial (0 or 1-order) --> But after some experimets, even the polyprder=1 case was found to cause similar problems.
2. Setting poly_regularization to a value greater than 0. However, I am unsure how to determine an appropriate value for this parameter.
3. Fit the major emission lines (e.g., Hα, Hβ, [OIII]5007) with Gaussians before performing the SED fit, subtract the continuum component, and mask the regions without emission lines. Then, input this processed spectrum into Prospector and use LineSpecModel with linespec_scaling=1 fixed during fitting. However, as I couldn’t find an example demonstration, I’m uncertain if this approach is correct.

Which approach would you recommend? Or do you have other suggestions for mitigating the effects of over-calibration while ensuring accurate physical parameter estimation?

Thank you in advance!

[bd-j]

Hello @tjtakys, sorry for the delay. I suggest beginning by trying to vary poly_regularization. This is effectively the width of a Gaussian prior (centered on zero) for the amplitude of the Chebyshev polynomial term(s). It can be a vector of length polyorder, in which case you usually want the width to get narrower for the higher order terms.

[bd-j]

Approach 3) would also work, though you only need to input the integrated continuum subtracted line fluxes from your Gaussian fit, not the full spectrum (i.e. no masking is required). And if you are concerned about spectral calibration you may still need to leave linespec_scaling free (though perhaps with a tight prior)

[tjtakys]

Hi Ben, thank you for your reply! After some trial and error, I’m now using LineSpecModel for the fitting. For the line flux inputs to LineSpecModel, should the units be erg/s/cm^2?

[bd-j]

Line fluxes are assumed erg/s/cm^2

[tjtakys]

OK, thank you very much for all of your help!