New science pack to accommodate modified Drude profiles for PAHs/Dust.

[els1]

Implement modified (highly asymmetric) Drude profiles for fitting and train on PDRs to set constraints on parameters for science pack.

[jdtsmith]

I'm going to include a mod_asym attribute as an option for dust features. This will "turn on" the Modified Drude in the science pack. We do need to figure out what the scale of that is. Relative to FWHM perhaps? So a value of 1.5 would mean "extra asymmetry 1.5x FWHM to red side" or similar?

[karllark]

Interesting.

The size of the asymmetry is not easy to translate.

[jdtsmith]

@Ameek-Sidhu, how does the asymmetry parameterization look now? I think a "relative to the FWHM" asymmetry could be good. And if FWHM is allowed to vary, that could be "relative to the initial FWHM".

[jdtsmith]

Here's a snippet of the current idea (not for actual use, just an example):

# Dust features [attributes: wavelength, power, fwhm, mod_asym]
PAH77:
    kind: dust_feature # Drude profiles
    features:
        PAH77a:
            wavelength: 7.42
            fwhm: 0.935
        PAH77b:
            wavelength: 7.60
            fwhm: 0.3344
        PAH77c:
            wavelength: 7.85
            fwhm:
                value: 0.416
                bounds: [-5%, 15%]
            mod_asym: 0.25 # indicates a modified/enhanced-asymmetry Drude profile.

[Ameek-Sidhu]

@jdtsmith, currently the asymmetry is incorporated in the fwhm as follows fwhm = 2*fwhm_0/(1+np.exp(a*(x-x_0))) Here fwhm_0 is the one without asymmetry, a is the asymmetry parameter, and x is wavelength, x_0 is central wavelength

[karllark]

For where the modified Drude was taken from, see Eqs 4 and 5 in Gordon et al. (2021, https://ui.adsabs.harvard.edu/abs/2021ApJ...916...33G/abstract). Based on ideas given in https://www.sciencedirect.com/science/article/abs/pii/S0924203108000453?via%3Dihub. Check out this latter reference for the motivation for this functional form of the asymmetry parameter - avoids extremes.

[jdtsmith]

Right, I remember. Since a has units of inverse wavelength, I propose we alter the equation internally to read:

fwhm(x) = 2*fwhm_0/(1+np.exp((x-x_0)/(a * fwhm_0)))

to make a dimensionless. I.e. suppose a=0.5 * fwhm_0 , at the half-max point of a normal Drude, i.e. x-x_0 = fwhm_0/2, this would correspond to a fwhm = fwhm_0 * 2/(1 + e) = 0.54 fwhm_0.

Edit: Or, since I guess a is usually negative because features have extra longward asymmetry (right @Ameek-Sidhu?), we could cast as:

fwhm(x) = 2*fwhm_0/(1+np.exp(-(x-x_0)/(a * fwhm_0)))

So at the same position fwhm = fwhm_0 * 2/(1 + 1/e) = 1.46 fwhm_0.

[jdtsmith]

Question though: do we have a closed form for the integrated intensity of a ModifiedDrude? We need that both for power estimation at the end of the fit, but also for an AreaModifiedDrude1D flavor during joint fitting.

[karllark]

Question though: do we have a closed form for the integrated intensity of a ModifiedDrude? We need that both for power estimation at the end of the fit, but also for an AreaModifiedDrude1D flavor during joint fitting.

Good question. Will investigate next week when I have time to work on the Area versions. I feel this should be straightforward, but need to work out the math.

[karllark]

Right, I remember. Since a has units of inverse wavelength, I propose we alter the equation internally to read:

fwhm(x) = 2*fwhm_0/(1+np.exp((x-x_0)/(a * fwhm_0)))

to make a dimensionless. I.e. suppose a=0.5 * fwhm_0 , at the half-max point of a normal Drude, i.e. x-x_0 = fwhm_0/2, this would correspond to a fwhm = fwhm_0 * 2/(1 + e) = 0.54 fwhm_0.

Edit: Or, since I guess a is usually negative because features have extra longward asymmetry (right @Ameek-Sidhu?), we could cast as:

fwhm(x) = 2*fwhm_0/(1+np.exp(-(x-x_0)/(a * fwhm_0)))

So at the same position fwhm = fwhm_0 * 2/(1 + 1/e) = 1.46 fwhm_0.

The asymmetry parameter is already dimensionless. So can you provide a bit more motivation and details for this change? Is it just a re-formulation of the existing equation? Or is it a replacement?

[jdtsmith]

The asymmetry parameter is already dimensionless.

Can't be. In the original version above, it must have units inverse to whatever x has, i.e. inverse wavelength. This is because together they are an argument of the exponential function, which has all polynomial powers present.

The advantage of making a dimensionless is it becomes agnostic to wavelength unit, and, in this formulation, has a more tractable/interpretable value, e.g a=0.5 means "extend the FWHM to the red by ~50%" (in my 2nd version with the negative sign).

But let's see what we can figure out on the Area version, which is important. May not be worth the effort up front. It also occurred to me this morning that using ModDrude's would make it harder to update models ala DL07/D21/etc., to incorporate JWST-based tweaks to the PAH cross-sections.

[karllark]

Ahh, yes. I see this now. Guess I didn't worry about this as it is a parameter that it is not clear has a physical interpretation.

In your statement, what does extend the FWHM to the red by 50% mean? I don't have a picture in my mind of how to interpret this. Maybe you could provide plots to illustrate? Is there something about the area redward and blueward of the feature center?

[karllark]

The fwhm is continuously changing in the modDrude equation. Not yet making the connection with your comments above between the asymmetry parameter and fwhm_o.

[karllark]

Maybe something with what you talk about above, but giving the fwhm blueward and redward at the half-max point?

[karllark]

Or maybe this would be best expressed as the red and blue hwhm?

[jdtsmith]

It's just a place to put your finger and say "an extra asymmetry of a=0.5 extends the red side about 50% at the half max point".

[karllark]

Interesting. I wonder if it is then 50% less on the other side.

[jdtsmith]

Other side is just 2/(1+e)=0.54 of the normal FWHM. This is for the last "negative formulation".

[Ameek-Sidhu]

@jdtsmith When I ran the code earlier leaving asymmetry as the free parameter, a had positive values as well for some of the Drudes.

Right, I remember. Since a has units of inverse wavelength, I propose we alter the equation internally to read:

fwhm(x) = 2*fwhm_0/(1+np.exp((x-x_0)/(a * fwhm_0)))

to make a dimensionless. I.e. suppose a=0.5 * fwhm_0 , at the half-max point of a normal Drude, i.e. x-x_0 = fwhm_0/2, this would correspond to a fwhm = fwhm_0 * 2/(1 + e) = 0.54 fwhm_0.

Edit: Or, since I guess a is usually negative because features have extra longward asymmetry (right @Ameek-Sidhu?), we could cast as:

fwhm(x) = 2*fwhm_0/(1+np.exp(-(x-x_0)/(a * fwhm_0)))

So at the same position fwhm = fwhm_0 * 2/(1 + 1/e) = 1.46 fwhm_0.