Numerical integration

Numerical integration#

%matplotlib widget

The dispersion integral of the phasespace factor cannot be solved analytically for higher angular momenta (\(\ell > 0\)). We therefore need to perform a numerical integration. However, since the integrand of the dispersion integral contains singularities, it is important to choose the numerical integration algorithm.

This notebook contains an interactive widget that shows the numerical stability and performance of different numerical integration algorithms when computing the dispersion integral along (close to) the real axis of the complex energy plane.

Tip

The conclusion is that Romberg’s method (quadax.romberg()) is the best choice for computing the dispersion integral close to the physical (real) axis. For computing the dispersion integral further away from the axis, the Gauss–Kronrod method (e.g. quadax.quadgk()) is accurate enough and faster than Romberg’s method.

Use quadax.romberg() when computing the dispersion integral for an amplitude model fit (physical axis).
Use quadax.quadgk() when determining the pole positions of the amplitude model (complex plane).

Chew–Mandelstam dispersion integral#

The general form of the Chew–Mandelstam dispersion integral \(\Sigma_\ell(s)\) is an integral over the product \(\rho(s) \, n_\ell^2(s)\) of the phase space factor \(\rho(s)\) (PhaseSpaceFactor) and the square of the form factor \(n_\ell^2(s)\) (FormFactor). The required formulas are:

\[\begin{split}\displaystyle \begin{aligned} \Sigma_{\ell}\left(s\right) \;&=\; \frac{s - \left(m_{1} + m_{2}\right)^{2}}{\pi} \int\limits_{\left(m_{1} + m_{2}\right)^{2}}^{\infty} \frac{\mathcal{F}_{\ell}\left(x, m_{1}, m_{2}\right)^{2} \rho\left(x\right)}{\left(x - \left(m_{1} + m_{2}\right)^{2}\right) \left(- i \epsilon - s + x\right)}\, dx \\ \rho\left(s\right) \;&=\; \frac{\sqrt{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt{s - \left(m_{1} + m_{2}\right)^{2}}}{s} \\ \mathcal{F}_{\ell}\left(s, m_{1}, m_{2}\right) \;&=\; \sqrt{B_{\ell}^2\left(q^2\left(s\right)\right)} \\ q\left(s\right) \;&=\; \frac{\sqrt{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt{s - \left(m_{1} + m_{2}\right)^{2}}}{2 \sqrt{s}} \\ B_{0}^2\left(z\right) \;&=\; 1 \\ B_{1}^2\left(z\right) \;&=\; \frac{2 z}{z + 1} \\ B_{2}^2\left(z\right) \;&=\; \frac{13 z^{2}}{z^{2} + 3 z + 9} \\ B_{3}^2\left(z\right) \;&=\; \frac{277 z^{3}}{z^{3} + 6 z^{2} + 45 z + 225} \\ B_{4}^2\left(z\right) \;&=\; \frac{12746 z^{4}}{z^{4} + 10 z^{3} + 135 z^{2} + 1575 z + 11025} \\ B_{5}^2\left(z\right) \;&=\; \frac{998881 z^{5}}{z^{5} + 15 z^{4} + 315 z^{3} + 6300 z^{2} + 99225 z + 893025} \\ \end{aligned}\end{split}\]

where we have used a split square root for a cleaner cut structure in the complex energy plane.

Here, \(i\epsilon\) indicates that \(\Sigma_\ell(s)\) with \(s\in\mathbb{R}\) is formulated just above the real axis in order to avoid \(s'-s=0\). However, as can be seen in the widget below, \(i\epsilon\) is not required when \(s\in\mathbb{C}\), giving us the general form:

\[\begin{split} \begin{aligned} \Sigma_\ell\left(s\right) \;&=\; \frac{s - s_\mathrm{thr}}{\pi} \int\limits_{s_\mathrm{thr}}^{\infty} \frac {\rho\!\left(s'\right) n_\ell^2\!\left(s'\right) ds'} {\left(s' - s_\mathrm{thr}\right) \left(s'- s\right)} \\ n_\ell^2(s') \;&=\; \mathcal{F}_\ell^2\!\left(s', m_1, m_2\right) \\ s_\mathrm{thr} \;&=\; (m_1 + m_2)^2 \end{aligned} \end{split}\]

The reason is that the function \(\rho(s')\,n_\ell^2(s')\) does not have a branch cut along the real axis above the \(s_\mathrm{thr}\) threshold, so there is no need to approach some contour around such a potential discontinuity. The form of the dispersion integral is explained and derived here.

When using this form of the dispersion integral to compute \(\Sigma_\ell(s)\) for \(s\in\mathbb{R}\), the caller should use the fact that \(\lim_{\epsilon\to 0^+} \Sigma_\ell(s + i\epsilon)\). In the numerical implementation, that means giving s + epsilon * 1j as input when s is a array with floats. The challenge is to find a value for \(\epsilon\): the smaller the value for \(\epsilon\), the closer we are to the physical axis, but the less accurate the numerical integration becomes. The widget below allows you to explore this trade-off for different numerical integration algorithms.

Numerical implementation#

In this notebook, we implement the formulas listed above numerical rather than lambdifying the symbolic expressions, so that we have full control over the numerical implementation and to make the implementation recognizable to pure array-oriented workflows.

In the case of \(S\)-waves (\(\ell=0\)), we can compare the result of the integration to the analytical solution to the integral:

\[\begin{split}\displaystyle \begin{aligned} \hat{\Sigma}_0\left(s\right) \;&=\; \frac{\frac{2 q^\mathrm{c}\left(s\right)}{\sqrt{s}} \log{\left(\frac{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q^\mathrm{c}\left(s\right) - s}{2 m_{1} m_{2}} \right)} - \left(m_{1}^{2} - m_{2}^{2}\right) \left(- \frac{1}{\left(m_{1} + m_{2}\right)^{2}} + \frac{1}{s}\right) \log{\left(\frac{m_{1}}{m_{2}} \right)}}{\pi} \\ \end{aligned}\end{split}\]

Interactive visualization#

Define user interface

Hide code cell source

cont = dict(continuous_update=False)
physics_sliders = dict(
    projection=w.RadioButtons(
        description="Projection",
        options=["real", "imag", "abs"],
        value="imag",
        layout=w.Layout(width="max-content"),
    ),
    m1=w.FloatSlider(value=0.13, min=0.0, max=2.0, step=0.01, description="m₁", **cont),
    m2=w.FloatSlider(value=0.98, min=0.0, max=2.0, step=0.01, description="m₂", **cont),
    ell=w.IntSlider(value=0, min=0, max=5, description="ℓ", **cont),
    y_lim=w.FloatRangeSlider(
        description="y range",
        min=-5,
        max=10,
        value=(-1, +1),
        readout_format=".1f",
        **cont,
    ),
    z_max=w.FloatLogSlider(
        value=1.0,
        min=-3,
        max=3,
        description="Color scale",
        step=0.25,
        readout_format=".3g",
        **cont,
    ),
    resolution=w.IntSlider(
        value=200,
        min=100,
        max=5000,
        description="Resolution",
        step=100,
        **cont,
    ),
    epsilon=w.FloatLogSlider(
        value=1e-4,
        min=-12,
        max=0.5,
        description="s + iϵ",
        step=0.25,
        readout_format=".0e",
        **cont,
    ),
    start_offset=w.FloatLogSlider(
        value=1e-20,
        min=-64,
        max=2,
        description="thr + ϵ",
        step=0.25,
        readout_format=".0e",
        **cont,
    ),
)
algorithm_sliders = dict(
    algorithm_name=w.RadioButtons(
        description="Integration algorithm",
        layout=w.Layout(width="150px"),
        options=Algorithm.__args__,
    ),
    order=w.RadioButtons(
        description="Integration order",
        layout=w.Layout(width="100px"),
        options=[8, 16, 32, 64, 128, 256],
        value=256,
    ),
    epsabs=w.FloatLogSlider(
        value=1e-5,
        min=-12,
        max=0,
        description="epsabs",
        step=0.5,
        readout_format=".0e",
        **cont,
    ),
    epsrel=w.FloatLogSlider(
        value=1e-5,
        min=-12,
        max=0,
        description="epsrel",
        step=0.5,
        readout_format=".0e",
        **cont,
    ),
    divmax=w.IntSlider(
        value=20,
        min=1,
        max=30,
        description="divmax",
        **cont,
    ),
    limit=w.IntSlider(
        value=50,
        min=1,
        max=100,
        description="limit",
        disabled=True,
        **cont,
    ),
    disable_limit=w.Checkbox(value=True, description="Disable limit", **cont),
)
sliders = dict(
    **physics_sliders,
    **algorithm_sliders,
)


def on_algorithm_change(change):
    algorithm_name = change["new"]
    sliders["limit"].disabled = True
    order = sliders["order"]
    divmax_slider = sliders["divmax"]
    disable_limit = sliders["disable_limit"]
    if algorithm_name == "quad_vec":
        disable_limit.disabled = False
        divmax_slider.disabled = True
        order.disabled = True
    if algorithm_name in {"quadcc", "quadgk", "quadts"}:
        disable_limit.disabled = True
        divmax_slider.disabled = True
        order.disabled = False
        if algorithm_name == "quadcc":
            _update_order(order, [8, 16, 32, 64, 128, 256])
        if algorithm_name == "quadgk":
            _update_order(order, [15, 21, 31, 41, 51, 61])
        if algorithm_name == "quadts":
            _update_order(order, [41, 61, 81, 101])
    if algorithm_name in {"romberg", "rombergts"}:
        disable_limit.disabled = True
        divmax_slider.disabled = False
        order.disabled = True


def _update_order(slider, options: list[int]) -> None:
    current_value = slider.value
    slider.options = options
    slider.value = min(options, key=lambda o: abs(o - current_value))


w.jslink((sliders["disable_limit"], "value"), (sliders["limit"], "disabled"))
sliders["algorithm_name"].observe(on_algorithm_change, names="value")
sliders["algorithm_name"].value = "romberg"  # ty:ignore[invalid-assignment]  # trigger changes
timer_box = cast("w.ValueWidget", w.HTML())
ui = w.VBox([
    tabs := w.Tab([
        w.HBox([
            physics_sliders["projection"],
            w.VBox(list(physics_sliders.values())[1:5]),
            w.VBox(list(physics_sliders.values())[5:]),
        ]),
        w.HBox([
            algorithm_sliders["algorithm_name"],
            algorithm_sliders["order"],
            w.VBox(list(algorithm_sliders.values())[2:]),
        ]),
    ]),
    timer_box,
])
tabs.titles = ["Physics", "Integration"]

Functions for interactive widget