Formulate amplitude model

Show code cell content Hide code cell content

%config InlineBackend.figure_formats = ['svg']

import sympy as sp
from IPython.display import Math

from ampform.io import aslatex, improve_latex_rendering

improve_latex_rendering()

Generate transitions

In Generate transitions, we used generate_transitions() to create a list of allowed Transitions for a specific decay channel:

import qrules

reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["gamma", "pi0", "pi0"],
    allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"],
    allowed_interaction_types=["strong", "EM"],
    formalism="canonical-helicity",
)

Show code cell source Hide code cell source

import graphviz

dot = qrules.io.asdot(reaction, collapse_graphs=True)
graphviz.Source(dot)

Build model

We can now use the ReactionInfo to formulate an amplitude model. The type of this amplitude model is dependent on the formalism. The function get_builder() helps to get the correct amplitude builder class for this formalism:

from ampform import get_builder

model_builder = get_builder(reaction)
type(model_builder)

ampform.helicity.CanonicalAmplitudeBuilder

If we now use the HelicityAmplitudeBuilder.formulate() method of this builder, we get a HelicityModel without any dynamics:

model_no_dynamics = model_builder.formulate()

See also

Helicity versus canonical

Main expressions

A HelicityModel has a few attributes. The expression for the total intensity is given by intensity:

As can be seen in Generate transitions, the transition that this HelicityModel describes features only one decay topology. This means that the main intensity expression is relatively simple. In case there are more decay topologies, AmpForm aligns spin, which makes the main intensity expression more complicated.

model_no_dynamics.intensity

\[\displaystyle \sum_{m_{A}\in\left\{1,-1\right\}} \sum_{m_{0}\in\left\{1,-1\right\}} \sum_{m_{1}=0} \sum_{m_{2}=0}{\left|{A^{12}_{m_{A}, m_{0}, m_{1}, m_{2}}}\right|^{2}}\]

This shows that the main intensity is an incoherent sum of the amplitude for each spin projection combination of the initial and final states. The expressions for each of these amplitudes are provided with the amplitudes attribute. This is an OrderedDict, so we can inspect the first of these amplitudes as follows:

(symbol, expression), *_ = model_no_dynamics.amplitudes.items()

\[\begin{split}\displaystyle \begin{align*} A^{12}_{-1, -1, 0, 0} = & C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,-1}_{0,0,1,-1} C^{1,-1}_{1,-1,0,0} D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{-1,-1}\left(- \phi_{0},\theta_{0},0\right) \\ & + C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,-1}_{0,0,1,-1} C^{1,-1}_{1,-1,0,0} D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{-1,-1}\left(- \phi_{0},\theta_{0},0\right) \\ & + C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,-1}_{1,-1,0,0} C^{1,-1}_{2,0,1,-1} D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{-1,-1}\left(- \phi_{0},\theta_{0},0\right) \\ & + C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,-1}_{1,-1,0,0} C^{1,-1}_{2,0,1,-1} D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{-1,-1}\left(- \phi_{0},\theta_{0},0\right) \end{align*}\end{split}\]

The intensity and its amplitudes are recombined through the expression attribute. This is just a sympy.Expr, which we can pull apart by using its args (see Advanced Expression Manipulation). Here’s an example:

intensities = model_no_dynamics.expression.args
intensity_1 = intensities[0]
base, power = intensity_1.args
abs_arg = base.args[0]
amplitude_sum = abs_arg.args
some_amplitude = amplitude_sum[0]
some_amplitude

\[\displaystyle C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,1}_{0,0,1,1} C^{1,1}_{1,1,0,0} D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{1,1}\left(- \phi_{0},\theta_{0},0\right)\]

some_amplitude.doit()

\[\displaystyle C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \left(\frac{\cos{\left(\theta_{0} \right)}}{2} + \frac{1}{2}\right) e^{i \phi_{0}}\]

Parameters and kinematic variables

As can be seen, the expression contains several Symbols. Some of these represent (kinematic) variables, such as the helicity angles \(\phi_0\) and \(\theta_0\) (see get_helicity_angle_symbols() for the meaning of their subscripts). Others will later on be interpreted parameters when fitting the model to data.

The HelicityModel comes with expressions for these kinematic_variables, so that it’s possible to compute them from 4-momentum data.

kinematic_variables = set(model_no_dynamics.kinematic_variables)
kinematic_variables

{m_0, m_012, m_1, m_12, m_2, phi_0, phi_1^12, theta_0, theta_1^12}

theta = sp.Symbol("theta_0", real=True)
expr = model_no_dynamics.kinematic_variables[theta]
expr.doit()

\[\displaystyle \operatorname{acos}{\left(\frac{\left({p}_{12}\right)\left[:, 3\right]}{\sqrt{\sum_{\mathrm{axis1}}{\left({p}_{12}\right)\left[:, 1:\right]^{2}}}} \right)}\]

See also

Lorentz vectors

The remaining symbols in the HelicityModel are parameters. Each of these parameters comes with suggested parameter values (parameter_defaults), that have been extracted from the ReactionInfo object where possible:

Math(aslatex(model_no_dynamics.parameter_defaults))

\[\begin{split}\displaystyle \begin{array}{rcl} C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ \end{array}\end{split}\]

Helicity couplings

If you prefer to characterize the strength of each partial two-body decay, set use_helicity_couplings to True and formulate a model:

model_builder.config.use_helicity_couplings = True
model_with_couplings = model_builder.formulate()
Math(aslatex(model_with_couplings.parameter_defaults))

\[\begin{split}\displaystyle \begin{array}{rcl} H_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma} &=& 1+0i \\ H_{f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ H_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma} &=& 1+0i \\ H_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma} &=& 1+0i \\ H_{f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ H_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma} &=& 1+0i \\ \end{array}\end{split}\]

Scalar masses

By default, the HelicityAmplitudeBuilder creates sympy.Exprs for each kinematic variable, including all ‘invariant’ final state masses (\(m_0, m_1, \dots\)). However, it often happens that certain particles in a final state are stable. In that case, you may want to substitute these symbols with scalar values. This can be achieved by specifying which final state IDs are to be considered stable. Their corresponding mass symbols will then be considered parameters and a scalar suggested parameter value will be provided.

model_builder.config.stable_final_state_ids = [0, 1, 2]
model_stable_masses = model_builder.formulate()
Math(aslatex(model_stable_masses.parameter_defaults))

\[\begin{split}\displaystyle \begin{array}{rcl} C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ m_{0} &=& 0.0 \\ m_{1} &=& 0.1349768 \\ m_{2} &=& 0.1349768 \\ \end{array}\end{split}\]

set(model_stable_masses.kinematic_variables)

{m_012, m_12, phi_0, phi_1^12, theta_0, theta_1^12}

We can reset this option as follows:

model_builder.config.stable_final_state_ids = None

Similarly, it can happen that the initial state mass is fixed (for instance when analyzing generated data or when applying a kinematic fit while reconstructing measured data). In that case, set scalar_initial_state_mass to True.

model_builder.config.scalar_initial_state_mass = True
model_stable_masses = model_builder.formulate()
Math(aslatex(model_stable_masses.parameter_defaults))

\[\begin{split}\displaystyle \begin{array}{rcl} C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} &=& 1+0i \\ m_{012} &=& 3.0969 \\ \end{array}\end{split}\]

Extend kinematic variables

The HelicityAmplitudeBuilder by default only generates kinematic_variables (helicity angles and invariant masses) for the topologies that are available in the ReactionInfo object that it was created with. If you want to calculate more kinematic variables, you can use the method register_topology() of its helicity HelicityAmplitudeBuilder.adapter to register more topologies and generate more kinematic variables. This is especially useful when generating data later on with TensorWaves.

To make this a bit easier, there is permutate_registered_topologies(), which is a small convenience function makes it possible to generate permutations of all registered_topologies and register them as well. Note that initially, only one Topology is registered in the HelicityAmplitudeBuilder.adapter, namely the one for the decay \(J/\psi \to \gamma f_0, f_0 \to \pi^0\pi^0\):

dot = qrules.io.asdot(model_builder.adapter.registered_topologies)
graphviz.Source(dot)

We now permutate_registered_topologies() to register permutations of this Topology:

model_builder.adapter.permutate_registered_topologies()

There are now three registered_topologies―one for each permutation:

len(model_builder.adapter.registered_topologies)

3

assert len(model_builder.adapter.registered_topologies) == 3
dot = qrules.io.asdot(model_builder.adapter.registered_topologies)
graphviz.Source(dot)

And if we formulate() a new HelicityModel, we see that it has many more kinematic_variables:

set(model_builder.formulate().kinematic_variables)

{m_0,
 m_01,
 m_02,
 m_1,
 m_12,
 m_2,
 phi_0,
 phi_01,
 phi_02,
 phi_0^01,
 phi_0^02,
 phi_1^12,
 theta_0,
 theta_01,
 theta_02,
 theta_0^01,
 theta_0^02,
 theta_1^12}

Tip

To register even more topologies, use e.g. create_isobar_topologies() to generate other, non-isomorphic topologies that cannot be created with permutations. This is relevant for more than three final states.

Set dynamics

To set dynamics for specific resonances, use DynamicsSelector.assign() on the same HelicityAmplitudeBuilder.dynamics attribute. You can set the dynamics to be any kind of Expr, as long as you keep track of which Symbol names you use (see Custom dynamics).

AmpForm does provide a few common dynamics functions, which can be constructed as Expr with the correct Symbol names using DynamicsSelector.assign(). This function takes specific dynamics.builder functions and classes, such as RelativisticBreitWignerBuilder, which can create relativistic_breit_wigner() functions for specific resonances. Here’s an example for a relativistic Breit-Wigner with form factor for the intermediate resonances and use a Blatt-Weisskopf barrier factor for the production decay:

from ampform.dynamics import PhaseSpaceFactor  # optional
from ampform.dynamics.builder import RelativisticBreitWignerBuilder

bw_builder = RelativisticBreitWignerBuilder(
    energy_dependent_width=True,
    form_factor=True,
    phsp_factor=PhaseSpaceFactor,  # optional
)
for name in reaction.get_intermediate_particles().names:
    model_builder.dynamics.assign(name, bw_builder)

Note that this RelativisticBreitWignerBuilder can also be initialized with a different PhaseSpaceFactorProtocol. This allows us to insert different phase space factors (see Analytic continuation and create_analytic_breit_wigner()).

See also

Custom dynamics

And we use the reconfigured HelicityAmplitudeBuilder to generate another HelicityModel, this time with relativistic Breit-Wigner functions and form factors:

model = model_builder.formulate()

Export

There is no special export function to export an HelicityModel. However, we can just use the built-in pickle module to write the model to disk and load it back:

import pickle

with open("helicity_model.pickle", "wb") as stream:
    pickle.dump(model, stream)
with open("helicity_model.pickle", "rb") as stream:
    model = pickle.load(stream)

Cached expression ‘unfolding’

Amplitude model expressions can be extremely large. AmpForm can formulate such expressions relatively fast, but sympy has to ‘unfold’ these expressions with doit(), which can take a long time. AmpForm provides a function that can cache the ‘unfolded’ expression to disk, so that the expression unfolding runs faster upon the next run.

from ampform.sympy import perform_cached_doit

full_expression = perform_cached_doit(model.expression)
sp.count_ops(full_expression)

815

See perform_cached_doit() for some tips on how to improve performance.

Visualize

Mathematical formula

It’s possible to view the complete amplitude model as an expression. This would, however, clog the screen here, so we instead just view the formula of one of its components:

some_amplitude, *_ = model.components.values()
some_amplitude.doit()

\[\displaystyle \frac{C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \Gamma_{f_{0}(980)} m_{f_{0}(980)} \left(\frac{\cos{\left(\theta_{0} \right)}}{2} + \frac{1}{2}\right) e^{- i \phi_{0}}}{- \frac{i \Gamma_{f_{0}(980)} m_{f_{0}(980)}^{2} \sqrt{\frac{\left(m_{12}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{12}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{12}^{2}}}}{m_{12} \sqrt{\frac{\left(m_{f_{0}(980)}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{f_{0}(980)}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{f_{0}(980)}^{2}}}} - m_{12}^{2} + m_{f_{0}(980)}^{2}}\]

Note

We use doit() to evaluate the Wigner-\(D\) (Rotation.D) and Clebsch-Gordan (CG) functions in the full expression.

The HelicityModel.parameter_defaults attribute can be used to substitute all parameters with suggested values:

some_amplitude.doit().subs(model.parameter_defaults)

\[\displaystyle \frac{0.0594 \left(\frac{\cos{\left(\theta_{0} \right)}}{2} + \frac{1}{2}\right) e^{- i \phi_{0}}}{- m_{12}^{2} + 0.9801 - \frac{0.05821794 i \sqrt{\frac{\left(m_{12}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{12}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{12}^{2}}}}{m_{12} \sqrt{\left(0.9801 - \left(m_{1} - m_{2}\right)^{2}\right) \left(0.9801 - \left(m_{1} + m_{2}\right)^{2}\right)}}}\]

Tip

To view the full expression for the amplitude model without crashing Jupyter Lab, install jupyterlab-katex.

Plotting

In this case (\(J/\psi \to \gamma f_0, f_0 \to \pi^0\pi^0\)) without dynamics, the total intensity is only dependent on the \(\theta\) angle of \(\gamma\) in the center of mass frame (see get_helicity_angle_symbols()):

no_dynamics = model_no_dynamics.expression.doit()
no_dynamics_substituted = no_dynamics.subs(model.parameter_defaults)
no_dynamics_substituted

\[\displaystyle 0.8 \sqrt{10} \cos^{2}{\left(\theta_{0} \right)} + 4.4 \cos^{2}{\left(\theta_{0} \right)} + 0.8 \sqrt{10} + 4.4\]

assert len(no_dynamics_substituted.free_symbols) == 1

Show code cell source Hide code cell source

theta = next(iter(no_dynamics_substituted.free_symbols))
sp.plot(
    no_dynamics_substituted,
    (theta, 0, sp.pi),
    axis_center=(0, 0),
    ylabel="$I$",
    ylim=(0, 16),
);

For this decay channel, the amplitude model is built up of four components:

no_dynamics.subs(zip(no_dynamics.args, sp.symbols("I_{:4}")))

\[\displaystyle I_{0} + I_{1} + I_{2} + I_{3}\]

This can be nicely visualized as follows:

Show code cell source Hide code cell source

plots = []
colors = ["red", "blue", "green", "purple"]

total = 0
for i, intensity in enumerate(no_dynamics.args):
    total += intensity.subs(model.parameter_defaults).doit()
    plots.append(
        sp.plot(
            total,
            (theta, 0, sp.pi),
            axis_center=(0, 0),
            ylabel="$I$",
            ylim=(0, 16),
            line_color=colors[i],
            show=False,
            label=f"$I_{i}$",
            legend=True,
        )
    )
for i in range(1, 4):
    plots[0].extend(plots[i])
plots[0].show()

Plot the model

Tip

See Inspect model interactively for a much more didactic way to visualize the model!

In the model with dynamics, we have several free symbols, such as the mass and width of the resonances. For the fitting package these will be considered parameters that are to be optimized and (kinematic) variables that represent the data set. Examples of parameters are mass (\(m_\text{particle}\)) and width (\(\Gamma_\text{particle}\)) of the resonances and certain amplitude coefficients (\(C\)). Examples of kinematic variables are the helicity angles \(\theta\) and \(\phi\) and the invariant masses (\(m_{ij...}\)).

sorted(model.expression.free_symbols, key=lambda s: s.name)

[C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
 C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
 C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
 C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
 \Gamma_{f_{0}(1500)},
 \Gamma_{f_{0}(980)},
 d_{f_{0}(1500)},
 d_{f_{0}(980)},
 m_1,
 m_12,
 m_2,
 m_{f_{0}(1500)},
 m_{f_{0}(980)},
 phi_0,
 phi_1^12,
 theta_0,
 theta_1^12]

Let’s say we want to plot the amplitude model with respect to \(m_{3+4}\). We would have to substitute all other free symbols with some value.

First, we can use HelicityModel.parameter_defaults to substitute the parameters with suggested values:

suggested_expression = model.expression.subs(model.parameter_defaults)
suggested_expression.free_symbols

{m_1, m_12, m_2, phi_0, phi_1^12, theta_0, theta_1^12}

Ideally, we would now ‘integrate out’ the helicity angles. Here, we however just set these angles to \(0\), as computing the integral would take quite some time:

angle = 0
angle_substitutions = {
    s: angle
    for s in suggested_expression.free_symbols
    if s.name.startswith("phi") or s.name.startswith("theta")
}
evaluated_angle_intensity = suggested_expression.subs(angle_substitutions)
evaluated_angle_intensity.free_symbols

{m_1, m_12, m_2}

By now we are only left with the masses of the final state particles (\(m_1\) and \(m_2\)), since they appear as symbols in the relativistic_breit_wigner_with_ff(). Final state particles 3 and 4 are the \(\pi^0\)’s, so we can just substitute them with their masses. (Alternatively, see Scalar masses.)

from qrules import load_pdg

pi0 = load_pdg()["pi0"]
plotted_intensity = evaluated_angle_intensity.doit().subs(
    {
        sp.Symbol("m_1", nonnegative=True): pi0.mass,
        sp.Symbol("m_2", nonnegative=True): pi0.mass,
    },
    simultaneous=True,
)

Tip

Use subs() with simultaneous=True, as that avoids a bug later on when plotting with matplotlib.pyplot.

That’s it! Now we are only left with the invariant mass \(m_{3+4}\) of the two pions:

assert len(plotted_intensity.free_symbols) == 1
m = next(iter(plotted_intensity.free_symbols))
m

\[\displaystyle m_{12}\]

…and we can plot it with with sympy.plot:

sp.plot(
    plotted_intensity,
    (m, 2 * pi0.mass, 2.5),
    axis_center=(2 * pi0.mass, 0),
    xlabel=Rf"$m(\pi^{0}\pi^{0})$",
    ylabel="$I$",
    backend="matplotlib",
);

The expression itself looks like this (after some rounding of the float values in this expression using tree traversal):

Show code cell source Hide code cell source

def round_nested(expression: sp.Expr, n_decimals: int = 2) -> sp.Expr:
    for node in sp.preorder_traversal(expression):
        if node.free_symbols:
            continue
        if isinstance(node, (float, sp.Float)):
            expression = expression.subs(node, round(node, n_decimals))
        if isinstance(node, sp.Pow) and node.args[1] == 1 / 2:
            expression = expression.subs(node, round(node.n(), n_decimals))
    return expression

Show code cell source Hide code cell source

rounded = round_nested(plotted_intensity)
round_nested(rounded)

\[\displaystyle 2 \left|{\frac{0.08}{m_{12}^{2} - 0.98 + \frac{0.06 i \sqrt{m_{12}^{2} - 0.07}}{m_{12}}} + \frac{0.22}{m_{12}^{2} - 2.32 + \frac{0.17 i \sqrt{m_{12}^{2} - 0.07}}{m_{12}}}}\right|^{2}\]