Update StarkBroadenedLine line shape model to take into account Doppler broadening and Zeeman splitting

[vsnever]

Currently, StarkBroadenedLine line shape model based on Lomanowski B.A. et al, 2015, Nucl. Fusion, 55 123028, only includes a modified Lorentzian for modelling static and dynamic Stark broadening of hydrogen lines:

$$L(\lambda, \Delta\lambda_{1/2}^{(L)})=\frac{C_0(\Delta\lambda_{1/2}^{(L)})^{3/2}}{\lambda^{5/2}+(\frac{\Delta\lambda_{1/2}^{(L)}}{2})^{5/2}},$$

thus neglecting Zeeman splitting and Doppler broadening: $$G(\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(G)}),$$ where $$\Delta\lambda_{1/2}^{(G)}=2\sqrt{2ln(2)}\sigma.$$

However, in Lomanowski et al these effects are taken into account.

The combination of Doppler and Stark broadenings gives the Voigt profile:
$$V(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)})=\int_{-\infty}^{\infty}G(\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(G)})L(\lambda-\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(L)})d\lambda'$$

Calculating such a convolution for each spectral interval at each spatial point, given that thousands of paths are being traced, would take too much time. Therefore, I propose to replace this convolution with a weighted sum (so-called pseudo-Voigt):

$$V(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)})=\eta L(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(V)}) + (1-\eta)G(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(V)}),$$

with $$\Delta\lambda_{1/2}^{(V)}\equiv F(\Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)}),$$ $$\eta\equiv f(\Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)}).$$

Both function, F and f can be obtained by fitting the true Voigt profile. The F function can be fitted with 6th degree polynomials of
$$x=\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(G)}} \textrm{ for } \frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(G)}} \le 1$$ $$\textrm{ and of } x=\frac{\Delta\lambda_{1/2}^{(G)}}{\Delta\lambda_{1/2}^{(L)}}\textrm{ for } \frac{\Delta\lambda_{1/2}^{(G)}}{\Delta\lambda_{1/2}^{(L)}} \le 1:$$

The f function can be fitted with 5th degree logarithmic polynomial of $$\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(V)}}$$ in the range from 0.01 to 0.999:

The maximum relative error of the pseudo-Voigt approximation can reach 25%, but all the main characteristics of the spectrum: maximum intensity, FWHM, line wings are approximated well.

Here is the script that does the above fitting:

pseudo_voigt_fit.py

import numpy as np
from scipy.special import wofz
from scipy.optimize import minimize, brentq, minimize_scalar
from matplotlib import pyplot as plt


def gaussian(x, sigma):
    return np.exp(-0.5 * x**2 / sigma**2) / (sigma * np.sqrt(2. * np.pi))


def lorentzian(x, lambda_1_2):
    return (0.5 / np.pi) * lambda_1_2 / (x**2 + (0.5 * lambda_1_2)**2)


def lorentzian25(x, lambda_1_2):
    """
    Modified Lorentzian from B. Lomanowski, et al. Nuclear Fusion 55.12 (2015)
    123028, https://doi.org/10.1088/0029-5515/55/12/123028
    """
    return (lambda_1_2**1.5 / 7.47443) / (np.abs(x)**2.5 + (0.5 * lambda_1_2)**2.5)


def voigt(x, sigma, lambda_1_2):
    """
    Exact computation of Voigt profile (gauss * lorentzian) using the Fadeeva function.
    Useful to test the convolution integration.

    :param x: Function argument.
    :param sigma: Gausian sigma.
    :param lambda_1_2: Lorentzian FWHM.
    """

    z = (x + 1j * 0.5 * lambda_1_2) / (sigma * np.sqrt(2))

    return np.real(wofz(z)) / (sigma * np.sqrt(2 * np.pi))


def voigt_convolution(x, fwhm_gauss, fwhm_lorentz, mg=512, ml=64, modified_lorentzian=True):
    """
    Computes Voigt profile as a convolution of Gaussian and Lorentzian profiles.

    :param x: Function argument.
    :param fwhm_gauss: Gausian FWHM.
    :param fwhm_lorentz: Lorentzian FWHM.
    :param mg: Number of grid points per Gausian FWHM. Default: 512.
    :param ml: Number of grid points per Lorentzian FWHM. Default: 64.
    :param modified_lorentzian: Use modified Lorentzian function. Default: True.
    """

    x = np.asarray(x)
    if x.ndim == 0:
        x = x[None]
    res = np.zeros(len(x))
    m = 8 * int(max(ml * fwhm_gauss / fwhm_lorentz, mg))
    x1, dx1 = np.linspace(- 4 * fwhm_gauss, 4 * fwhm_gauss, m, retstep=True)
    sigma = fwhm_gauss / (2 * np.sqrt(2 * np.log(2)))
    f1 = gaussian(x1, sigma)
    for ix, x0 in enumerate(x):
        f2 = lorentzian25(x1 - x0, fwhm_lorentz) if modified_lorentzian else lorentzian(x1 - x0, fwhm_lorentz)
        res[ix] = (f1 @ f2) * dx1

    return res


def pseudo_voigt(x, fwhm_gauss, fwhm_lorentz, voigt_fwhm_func, weight_func, modified_lorentzian=True):
    """
    Computes pseudo-Voigt profile as a weighted sum of Gaussian and lorentzian profiles.
    F_voigt(x) = weight * F_lorentz + (1. - weight) * F_gauss

    :param x: Function argument.
    :param fwhm_gauss: Gausian FWHM.
    :param fwhm_lorentz: Lorentzian FWHM.
    :param voigt_fwhm_func: A function of fwhm_gauss and fwhm_lorentz that returns
                            approximated Voigt FWHM, fwhm_voigt.
    :param weight_func: A function of fwhm_lorentz/fwhm_voigt that returns
                        a weight of Lorentzian profile.
    :param modified_lorentzian: Use modified Lorentzian function. Default: True.
    """

    fwhm = voigt_fwhm_func(fwhm_gauss, fwhm_lorentz)
    weight = weight_func(fwhm_lorentz / fwhm)

    sigma_v = 0.5 * fwhm / np.sqrt(2 * np.log(2))
    lorentz = lorentzian25(x, fwhm) if modified_lorentzian else lorentzian(x, fwhm)
    gs = gaussian(x, sigma_v)

    return weight * lorentz + (1 - weight) * gs


def voigt_find_fwhm(fwhm_gauss, fwhm_lorentz):
    """
    Finds Voigt FWHM.

    :param fwhm_gauss: Gausian FWHM.
    :param fwhm_lorentz: Lorentzian FWHM.
    """
    fvmax = voigt_convolution(0, fwhm_gauss, fwhm_lorentz)[0]

    def root_func(x):
        return voigt_convolution(x, fwhm_gauss, fwhm_lorentz)[0] - 0.5 * fvmax

    return 2. * brentq(root_func, 0, fwhm_gauss + fwhm_lorentz, maxiter=1000)


def fwhm_polyfit(x, y, deg):
    """
    Constrained polynomial fit: a0 = 1, a1 >= 0.
    """

    def min_func(w):
        fit = 1.
        for i, w1 in enumerate(w):
            fit += w1 * x**(i + 1)

        return ((1. - fit / y)**2).sum()

    constraints = [{'type': 'ineq', 'fun': lambda w: w[0]}]

    return np.append(minimize(min_func, np.zeros(deg), constraints=constraints).x[::-1], [1.])


def pseudo_voigt_fwhm_fit(deg_gauss=6, deg_lorentz=6, n=1000):
    """
    Fits Voigt FWHM with two polynoms of given degrees.

    For FWHM_{G} <= FWHM_{L}:
        FWHM_{V} = 1. + Sum_{i=1}^{N_{G}} a_i * (FWHM_{G}/FWHM_{L})^i

    For FWHM_{G} > FWHM_{L}:
        FWHM_{V} = 1. + Sum_{i=1}^{N_{L}} b_i * (FWHM_{L}/FWHM_{G})^i

    :param deg_gauss: N_{G}
    :param deg_lorentz: N_{L}
    :param n: FWHM_{G}/FWHM_{L} and FWHM_{L}/FWHM_{G} grid sizes.

    :returns: a_i, b_i
    """

    fwhm_gauss = np.geomspace(0.01, 1, n)
    fwhm_voigt = np.zeros(n)

    for i, g in enumerate(fwhm_gauss):
        fwhm_voigt[i] = voigt_find_fwhm(g, 1.)

    poly_gauss = fwhm_polyfit(fwhm_gauss, fwhm_voigt, deg_gauss)

    fwhm_lorentz = np.geomspace(0.01, 1, n)
    fwhm_voigt = np.zeros(n)

    for i, l in enumerate(fwhm_lorentz):
        fwhm_voigt[i] = voigt_find_fwhm(1., l)

    poly_lorentz = fwhm_polyfit(fwhm_lorentz, fwhm_voigt, deg_lorentz)

    return poly_gauss, poly_lorentz


def pseudo_voigt_weight_fit(fwhm_gauss, fwhm_lorentz, voigt_fwhm_func, n, relative=True):
    """
    Fits the Lorentzian weight in the pseudo-Voigt profile.

    :param fwhm_gauss: Gausian FWHM.
    :param fwhm_lorentz: Lorentzian FWHM.
    :param voigt_fwhm_func: A function of fwhm_gauss and fwhm_lorentz that returns
                            approximated Voigt FWHM.
    :param n: The argument grid size for fitting the pseudo-Voigt profile.
    :param relative: Relative (shape) or absolute (value) fit. Default: True (relative).
    """

    fwhm_voigt = voigt_fwhm_func(fwhm_gauss, fwhm_lorentz)

    if relative:
        x = np.geomspace(0.01 * fwhm_voigt, 50 * fwhm_voigt, n)
    else:
        x = np.linspace(0, 50 * fwhm_voigt, n)
    fv = voigt_convolution(x, fwhm_gauss, fwhm_lorentz)

    fl = lorentzian25(x, fwhm_voigt)
    fg = gaussian(x, 0.5 * fwhm_voigt / np.sqrt(2 * np.log(2)))

    def fit_func(x):
        return np.abs(fv - x * fl - (1. - x) * fg).sum()

    def rel_fit_func(x):
        return ((1. - (x * fl + (1. - x) * fg) / fv)**2).sum()

    f = rel_fit_func if relative else fit_func

    result = minimize_scalar(f, bounds=(0, 1), method='Bounded')

    return result.x


n = 8192
nw = 256

poly_gauss, poly_lorentz = pseudo_voigt_fwhm_fit(deg_gauss=6, deg_lorentz=6, n=1000)

poly_gauss_func = np.poly1d(poly_gauss)
poly_lorentz_func = np.poly1d(poly_lorentz)


def voigt_fwhm_func(fwhm_gauss, fwhm_lorentz):
    if fwhm_gauss > fwhm_lorentz:
        return fwhm_gauss * poly_lorentz_func(fwhm_lorentz / fwhm_gauss)

    return fwhm_lorentz * poly_gauss_func(fwhm_gauss / fwhm_lorentz)


fwhm_lorentz = np.geomspace(0.01, 100, nw)
weight = np.zeros(nw)
fwhm_ratio = np.zeros(nw)
for i in range(nw):
    fwhm_ratio[i] = fwhm_lorentz[i] / voigt_fwhm_func(1., fwhm_lorentz[i])
    print('Step {} out of {}. Fitting lorentzian weight for FWHM_lorentz/FWHM_Voigt = {}.'.format(i + 1, nw, fwhm_ratio[i]))
    weight[i] = pseudo_voigt_weight_fit(1., fwhm_lorentz[i], voigt_fwhm_func, n, relative=True)


deg = 5
poly_weight = np.polyfit(np.log(fwhm_ratio), np.log(weight), deg)


def weight_func(fwhm_ratio):
    return np.exp(np.poly1d(poly_weight)(np.log(fwhm_ratio)))


fwhm_gauss = np.linspace(0.01, 1., 100)
fwhm_voigt = np.array([voigt_find_fwhm(fg, 1.) for fg in fwhm_gauss])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(fwhm_gauss, fwhm_voigt, ls='-', c='tab:blue', label='Voigt FWHM')
ax.plot(fwhm_gauss, poly_gauss_func(fwhm_gauss), ls=':', c='b', label='Polynomial fit of degree {}'.format(6))
ax.set_xlabel(r'$\frac{FWHM_{Gauss}}{FWHM_{Lorentz}}$')
ax.legend(loc=0)

fwhm_lorentz = np.linspace(0.01, 1., 100)
fwhm_voigt = np.array([voigt_find_fwhm(1., fl) for fl in fwhm_lorentz])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(fwhm_lorentz, fwhm_voigt, ls='-', c='tab:blue', label='Voigt FWHM')
ax.plot(fwhm_lorentz, poly_lorentz_func(fwhm_lorentz), ls=':', c='b', label='Polynomial fit of degree {}'.format(6))
ax.set_xlabel(r'$\frac{FWHM_{Lorentz}}{FWHM_{Gauss}}$')
ax.legend(loc=0)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(fwhm_ratio, weight, ls='-', c='tab:blue', label='Lorentzian weight coefficient')
ax.plot(fwhm_ratio, weight_func(fwhm_ratio), ls=':', c='b', label='Logarithmic polynomial fit of degree {}'.format(deg))
ax.set_xlabel(r'$\frac{FWHM_{Lorentz}}{FWHM_{Voigt}}$')
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend(loc=0)

fwhm_lorentz = [0.01, 0.1, 0.5, 1., 2., 10.]

fig = plt.figure(figsize=(10., 10.), dpi=80)

for i, fwhm_l in enumerate(fwhm_lorentz):
    fwhm_v = voigt_fwhm_func(1., fwhm_l)
    x = np.linspace(-40 * fwhm_v, 40 * fwhm_v, n)
    f_voigt = voigt_convolution(x, 1., fwhm_l)
    f_pseudo = pseudo_voigt(x, 1., fwhm_l, voigt_fwhm_func, weight_func)

    ax1 = fig.add_subplot(len(fwhm_lorentz), 3, 3 * i + 1)
    ax1.plot(x, f_voigt, ls='-', c='tab:blue', label='True Voigt profile')
    ax1.plot(x, f_pseudo, ls=':', c='b', label='Pseudo-Voigt profile')
    ax1.set_yscale('log')
    ax1.text(0.02, 0.95, r'$\frac{FWHM_{Lorentz}}{FWHM_{Gauss}} = $' + '{:g}'.format(fwhm_l), va='top', transform=ax1.transAxes, fontsize=11)

    ax2 = fig.add_subplot(len(fwhm_lorentz), 3, 3 * i + 2)
    ax2.plot(x, f_voigt, ls='-', c='tab:blue', label='True Voigt profile')
    ax2.plot(x, f_pseudo, ls=':', c='b', label='Pseudo-Voigt profile')
    ax2.set_xlim(-2 * fwhm_v, 2 * fwhm_v)

    ax3 = fig.add_subplot(len(fwhm_lorentz), 3, 3 * i + 3)
    ax3.plot(x, 1. - f_pseudo / f_voigt, ls='-', c='k', label='Residual')
    ax3.set_xlim(-2 * fwhm_v, 2 * fwhm_v)

    if i == 0:
        ax1.text(0.99, 0.95, 'Logscale', ha='right', va='top', transform=ax1.transAxes)
        ax2.text(0.99, 0.95, 'Linear scale', ha='right', va='top', transform=ax2.transAxes)
        ax3.text(0.99, 0.95, 'Relative error', ha='right', va='top', transform=ax3.transAxes)
        ax1.legend(loc=2, ncol=2, frameon=False, bbox_to_anchor=(0., 1.25))

fig.subplots_adjust(left=0.05, bottom=0.03, right=0.98, top=0.97, wspace=0.23, hspace=0.18)

plt.show()

As for the Zeeman effect, I think for now we can do the same as in Lomanowski B.A. et al, namely model Stark and Zeeman effects separately (each Zeeman component is a Stark-Doppler line). Later we can try to interpolate tabular data from J. Rosato, Y. Marandet, R. Stamm, Journal of Quantitative Spectroscopy & Radiative Transfer 187 (2017) 333, where both effects are taken into account simultaneously on quantum level.