How to inform balmer decrement Av prior to prospector?

[tjtakys]

Hi, I am using Prospector to simultaneously fit photometry and spectra in order to investigate the ratio between nebular Av and stellar continuum Av. I have observational data for the Ha/Hb ratio, which allows me to derive the nebular Av from the Balmer decrement. I would like to input this nebular Av information into Prospector to estimate the stellar continuum Av. (I often see in literatures that stellar continuum Av is derived by masking nebular lines during SED fitting, but since nebular lines are also important for constraining SFH, I have opted to include all available information within Prospector.)

Since nebular Av from the Balmer decrement includes contributions from both diffuse dust and birth cloud dust, I have defined my dust attenuation model as follows:

def dusttot_to_dust2(dust_tot=0.5, dust_mu=0.0, **extra):
    """Set the value of dust2 (diffuse) by multiplying fraction factor of the total tau_V (optical depth) calculated from the Balmer decrement measurement.(The definition of dust_mu is the same as the mu in MAGPHYS)

    Parameters
    ----------
    dust_tot : float
        The total tau_V calculated from the Balmer decrement measurement.
    mu : float (0<mu<1)
        The fraction of the total tau_V contributed by dust in the ambient (diffuse) ISM

    Returns
    -------
    dust2 : float
        The optical depth towards all stars (the FSPS ``dust2`` parameter.)
    """
    return dust_tot * dust_mu

def dusttot_to_dust1(dust_tot=0.5, dust_mu=0.0, **extra):
    return dust_tot * (1 - dust_mu)

# balmer_tauV_mean (and _sigma) is the optical depth (and its error) at V band calculated from Balmer decrement
prior_dusttot = priors.ClippedNormal(mini=0.0, maxi=10.0, mean=balmer_tauV_mean, sigma=balmer_tauV_sigma)

model_params["dust_type"] = {'N':1, 'isfree':False, 'init':4} #4: Kriek & Conroy 2013
model_params["dust_tesc"] = {'N':1, 'isfree':False, 'init':7.0} 
model_params["dust1"] = {'N':1, 'isfree':False, 'init':0.2, 'depends_on':self.dusttot_to_dust1}
model_params["dust2"] = {'N':1, 'isfree':False, 'init':0.2, 'depends_on':self.dusttot_to_dust2}
model_params["dust_tot"] = {'N':1, 'isfree':True, 'init':1.0, 'prior':prior_dusttot}
model_params["dust_mu"] = {'N':1, 'isfree':True, 'init':0.2, 'prior':priors.ClippedNormal(mini=0.0, maxi=1.0, mean=0.5, sigma=0.3)}
model_params["dust_index"] = {'N':1, 'isfree':True, 'init':-0.7, 'prior':priors.TopHat(mini=-1.0, maxi=0.4)}
model_params["dust1_index"] = {'N':1, 'isfree':False, 'init':-1.3}

However, when performing the fit, I encounter problematic solutions, such as a very high value of dust_mu, as indicated in the following trace plot. Additionally, the resulting spectra fit falls significantly below the observed flux, indicating poor fitting performance.

This behavior persists even when I specify a prior for dust_tot using priors.TopHat(mini=0.0, maxi=10.0), suggesting that the issue might not be due to incorrect Balmer decrement values, but rather a problem in how the model is configured.

When I use a typical model setup (as demonstrated in the demo files) without providing Balmer decrement information, the fitting goes well, but I am concerned about parameter degeneracy in such cases. Here is the code used in this standard model configuration:

model_params["dust_type"] = {'N':1, 'isfree':False, 'init':4} # Kriek & Conroy 2013
model_params["dust_tesc"] = {'N':1, 'isfree':False, 'init':7.0} 
model_params["dust1"] = {'N':1, 'isfree':False, 'init':0.2, 'depends_on':transforms.dustratio_to_dust1}
model_params["dust2"] = {'N':1, 'isfree':True, 'init':0.2, 'prior':priors.TopHat(mini=0.0, maxi=4.0)}
model_params["dust_ratio"] = {'N':1, 'isfree':True, 'init':1.0, 'prior':priors.ClippedNormal(mini=0.0, maxi=2.0, mean=1.0, sigma=0.3)}
model_params["dust_index"] = {'N':1, 'isfree':True, 'init':-0.7, 'prior':priors.TopHat(mini=-1.0, maxi=0.4)}
model_params["dust1_index"] = {'N':1, 'isfree':False, 'init':-1.3}

[bd-j]

Hi @tjtakys, are the Halpha and Hbeta measurements coming from the spectrum you are also fitting? If so then you should probably not include an additional prior, and let the lines inform the model. Otherwise, I think the most correct way to do this would be to include the Hbeta and Halpha line fluxes in the fit directly.

This is possible along with the photometry and spectroscpy by using the Lines observation class in the dev/v2.0 version, though this version has some other updates that may require some reworking of your script. If you don't actually have a detected continuum in the spectrum (which I'm guessing given the sigma_smooth posterior), you could also just fit the line fluxes and the photometry using prospector v1.4 with the LineSpecModel object.

Also I note the Balmer decrement to Av conversion generally assumes a dust_index value and that it corresponds to the same shape curve as the dust1_index, so the model is not quite self-consistent, especially since dust_index is going to very low values (which iirc means very steep curves)

Otherwise, I'm not sure what the issue is, I don't see anything obvious in the model specification, but it's clearly trying to put all the attenuation in the diffuse component and then make that very steep. Given the high optical depth I worry that the dust model is simply not complex enough (see e.g. Lower et al. )