Can you add support for Generalized Drude-Lorentz model?

[leSemaleon]

I have a material that's well approximated by the generalized Lorenz model, but not so well by the standard one. I'm referring to the model from the paper: "Optimizing the Drude-Lorentz model for material permittivity: Method, program, and examples for gold, silver, and copper"
Is it possible to add support for it to MEEP?

In essence, the dielectric constant is written in it using the formula:

When in meep:

[stevengj]

In the long run, I would like to see Meep simply support arbitrary IIR filters, which would give us arbitrary rational functions of $\omega$.

(Technically, an IIR filter is a rational function $e^{i\omega \Delta t}$, but in the $\Delta t \to 0$ limit this is equivalent to a rational function of $\omega$, though some additional care would be required to ensure second-order convergence with $\Delta t$.)

At first glance, the problem with an $i \omega$ term in the denominator is that it correponds to a $d\mathbf{E}/dt$ term on the right-hand side of the auxiliary ODE, which is nontrivial to implement to second-order accuracy unless there is a trick here I'm forgetting. If we're going to spend the effort on this kind of thing, I'd rather spend it on the more general case of IIR filters.

[leSemaleon]

I agree that a more general case is preferable.
I've tried making a subclass of Susceptibility with the required formula, but it doesn't work correctly.

Code

class GeneralizedLorentzianSusceptibility(mp.geom.Susceptibility):
    def __init__(self, Omega=0.0+1j*0.0, Sigma=0.0+1j*0.0, **kwargs):
        super().__init__(**kwargs)
        self.Omega = Omega
        self.Sigma = Sigma

    def eval_susceptibility(self, freq):
        sigma = np.expand_dims(
            mp.geom.Matrix(diag=self.sigma_diag, offdiag=self.sigma_offdiag), axis=0
        )
        return (
                (
                    1j*self.Sigma/(freq - self.Omega) +
                    1j*self.Sigma.conjugate()/(freq + self.Omega.conjugate())
                ) * sigma
            )

Propagation spectra of Medium with my function (custom) and meep function (legacy)

My material have increasing absorption after absorption peak, but using meep function I can't achieve this

Blue: Fit refl index (top) and absorp coef (bottom); Orange: they are from literature

May be there is a way to solve it with meep functional?

[stevengj]

eval_susceptibility is just for output and diagnostic purposes, it has no effect on the timestepping.

[leSemaleon]

Then I don't quite understand where the substance is defined. I can define a medium with a list of susceptibilities, and all these susceptibilities are classes specified in geom.py.

As far as I understand, the meep kernel runs in C. So, the same susceptibility functions are defined there? How do we pass parameters to LorentzSusceptibilities? Could you please explain?

[stevengj]

The substance is defined mainly in the update_P function, which is what timesteps the auxiliary differential equation for the polarizability.