Spin alignment implementation#

Hide code cell content
import inspect
import logging
import warnings

import ampform
import graphviz
import numpy as np
import qrules
import sympy as sp
from ampform.helicity import (
    formulate_helicity_rotation_chain,
    formulate_rotation_chain,
    formulate_spin_alignment,
    formulate_wigner_d,
)
from ampform.kinematics import (
    compute_boost_chain,
    compute_wigner_angles,
    compute_wigner_rotation_matrix,
    create_four_momentum_symbols,
)
from IPython.display import Math, display
from qrules.topology import create_isobar_topologies

LOGGER = logging.getLogger()
LOGGER.setLevel(logging.ERROR)
warnings.filterwarnings("ignore")


def show_transition(transition, **kwargs):
    if "size" not in kwargs:
        kwargs["size"] = 5
    dot = qrules.io.asdot(transition, **kwargs)
    display(graphviz.Source(dot))

Helicity formalism#

Imagine we want to formulate the amplitude for the following single {external+qrules-0.9.x:class}.StateTransition:

Hide code cell source
full_reaction = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", ("Sigma+", [+0.5]), ("p~", [+0.5])],
    allowed_intermediate_particles=["Sigma(1660)~-", "N(1650)+"],
    allowed_interaction_types="strong",
    formalism="helicity",
)
graphs = full_reaction.to_graphs()
single_transition_reaction = full_reaction.from_graphs(
    [graphs[0]], formalism=full_reaction.formalism
)
transition = single_transition_reaction.transitions[0]
show_transition(transition)

The specific spin_projections for each particle only make sense given a specific reference frame. AmpForm’s HelicityAmplitudeBuilder interprets these projections as the helicity \(\lambda=\vec{S}\cdot\vec{p}\) of each particle in the rest frame of the parent particle. For example, the helicity \(\lambda_2=+\tfrac{1}{2}\) of \(\bar p\) is the helicity as measured in the rest frame of resonance \(\bar\Sigma(1660)^-\). The reason is that these helicities are needed when formulating the two-particle state for the decay node \(\bar\Sigma(1660)^- \to K^0\bar p\) (see {external+ampform-0.14.x:doc}ampform:usage/helicity/formalism).

Ignoring dynamics and coefficients, the HelicityModel for this single transition is rather simple:

Hide code cell source
builder = ampform.get_builder(single_transition_reaction)
model = builder.formulate()
model.expression.subs(model.parameter_defaults).subs(1.0, 1)
\[\displaystyle \left|{D^{\frac{1}{2}}_{- \frac{1}{2},- \frac{1}{2}}\left(- \phi^{02}_{0},\theta^{02}_{0},0\right) D^{1}_{-1,-1}\left(- \phi_{02},\theta_{02},0\right)}\right|^{2}\]

The two Wigner-\(D\) functions come from the two two-body decay nodes that appear in the {external+qrules-0.9.x:class}.StateTransition above. They were formulated as follows:

sp.Mul(
    formulate_wigner_d(transition, node_id=0),
    formulate_wigner_d(transition, node_id=1),
)
\[\displaystyle D^{\frac{1}{2}}_{- \frac{1}{2},- \frac{1}{2}}\left(- \phi^{02}_{0},\theta^{02}_{0},0\right) D^{1}_{-1,-1}\left(- \phi_{02},\theta_{02},0\right)\]

Now, as {external+ampform-0.14.x:func}~.formulate_wigner_d explains, the numbers that appear in the Wigner-\(D\) functions here are computed from the helicities of the decay products. But there’s a subtle problem: these helicities are assumed to be in the rest frame of each parent particle. For the first node, this is fine, because the parent particle rest frame matches that of the initial state in the {external+qrules-0.9.x:class}.StateTransition above. In the second node, however, we are in a different rest frame. This can result in phase differences for the different amplitudes.

If there is a single decay Topology in the ReactionInfo object for which we are formulating an amplitude model, the problem we identified here can be ignored. The reason is that the phase difference for each {external+qrules-0.9.x:class}.StateTransition (with each an identical decay Topology) is the same and does not introduce interference effects within the coherent sum. It again becomes a problem, however, when we are formulating an amplitude model with different topologies. An example would be the following reaction:

Hide code cell source
show_transition(full_reaction, collapse_graphs=True)

When formulating the amplitude model for this reaction, the HelicityAmplitudeBuilder implements the ‘standard’ helicity formalism as described in [Richman, 1984, Kutschke, 1996, Chung, 2014] and simply sums over the different amplitudes to get the full amplitude:

Hide code cell source
builder = ampform.get_builder(full_reaction)
model = builder.formulate()
latex = sp.multiline_latex(
    sp.Symbol("I"),
    model.expression.subs(model.parameter_defaults).subs(1.0, 1),
    environment="eqnarray",
)
Math(latex)

As pointed out in [Marangotto, 2020, Mikhasenko et al., 2020, Wang et al., 2021], this is wrong because of the mismatch in reference frames for the helicities.

Aligning reference frames#

In the rest of this document, we follow [Marangotto, 2020] to align all amplitudes in the different topologies back to the initial state reference frame \(A\), so that they can be correctly summed up. Specifically, we want to formulate a new, correctly aligned amplitude \(\mathcal{A}^{A\to 0,1,\dots}_{m_A,m_0,m_1,\dots}\) from the original amplitudes \(\mathcal{A}^{A\to R,S,i,...\to 0,1,\dots}_{\lambda_A,\lambda_0,\lambda_1,\dots}\) by applying Eq.(45) and Eq.(47) for generic, multi-body decays. Here, the \(\lambda\) values are the helicities in the parent rest frame of each two-body decay and the \(m\) are the canonical[1] spin projections in the rest frame of the mother particle that is the same no matter the Topology.

Just as in [Marangotto, 2020], we test the implementation with 1-to-3 body decays. We use the notation from get_boost_chain_suffix() to indicate resonances \(R,S,U\). This results in the following figure for the two alignments sums of Equations (45) and (46) in [Marangotto, 2020]:

Hide code cell source
dot1 = """
digraph {
    bgcolor=none
    rankdir=LR
    edge [arrowhead=none]
    node [shape=none, width=0]
    A
    0 [fontcolor=red]
    1 [fontcolor=green, label=<<o>1</o>>]
    2 [fontcolor=blue, label=<<o>2</o>>]
    { rank=same A }
    { rank=same 0, 1, 2 }
    N0 [label=""]
    N1 [label=""]
    A -> N0 [style=dotted]
    N0 -> N1 [label="R = 01", fontcolor=orange]
    N1 -> 0
    N0 -> 2 [style=dashed]
    N1 -> 1 [style=dashed]
}
"""
dot2 = """
digraph {
    bgcolor=none
    rankdir=LR
    edge [arrowhead=none]
    node [shape=none, width=0]
    A
    0 [label=0, fontcolor=red]
    1 [label=1, fontcolor=green, label=<<o>1</o>>]
    2 [label=2, fontcolor=blue, label=<<o>2</o>>]
    { rank=same A }
    { rank=same 0, 1, 2 }
    N0 [label=""]
    N1 [label=""]
    A -> N0 [style=dotted]
    N0 -> N1 [label="S = 02", fontcolor=violet]
    N1 -> 0
    N0 -> 1 [style=dashed]
    N1 -> 2 [style=dashed]
}
"""
display(*map(graphviz.Source, [dot1, dot2]))

The dashed edges and bars above the state IDs indicate “opposite helicity” states. The helicity of an opposite helicity state gets a minus sign in the Wigner-\(D\) function for a two-body state as formulated by {external+ampform-0.14.x:func}.formulate_wigner_d (see Helicity formalism) and therefore needs to be defined consistently. AmpForm does this with {external+ampform-0.14.x:func}.is_opposite_helicity_state.

Opposite helicity states are also of importance in the spin alignment procedure sketched by [Marangotto, 2020]. The Wigner-\(D\) functions that appear in Equations (45) and (46) from [Marangotto, 2020], operate on the spin of the final state, but the angles in the Wigner-\(D\) function are taken from the sibling state:

(1)#\[\begin{split} \begin{eqnarray} \mathcal{A}^{A \to {\color{orange}R},2 \to 0,1,2}_{m_A,m_0,m_1,m_2} &=& \sum_{\lambda_0^{01},\mu_0^{01},\nu_0^{01}} {\color{red}{D^{s_0}_{m_0,\nu_0^{01}}}}\!\left({\color{red}{\alpha_0^{01}, \beta_0^{01}, \gamma_0^{01}}}\right) {\color{red}{D^{s_0}_{\nu_0^{01},\mu_0^{01}}}}\!\left({\color{orange}{\phi_{_{01}}, \theta_{_{01}}}}, 0\right) {\color{red}{D^{s_0}_{\mu_0^{01},\lambda_0^{01}}}}\!\left({\color{red}{\phi_0^{01}, \theta_0^{01}}}\right) \\ &\times& \sum_{\lambda_1^{01},\mu_1^{01},\nu_1^{01}} {\color{green}{D^{s_1}_{m_1,\nu_1^{01}}}}\!\left({\color{green}{\alpha_1^{01}, \beta_1^{01}, \gamma_1^{01}}}\right) {\color{green}{D^{s_1}_{\nu_1^{01},\mu_1^{01}}}}\!\left({\color{orange}{\phi_{_{01}}, \theta_{_{01}}}}, 0\right) {\color{green}{D^{s_1}_{\mu_1^{01},\lambda_1^{01}}}}\!\left({\color{red}{\phi_0^{01}, \theta_0^{01}}}\right) \\ &\times& \sum_{\lambda_2^{01}} {\color{blue}{D^{s_2}_{m_2,\lambda_2^{01}}}}\!\left({\color{orange}{\phi_{_{01}}, \theta_{_{01}}}}, 0\right) \\ &\times& \mathcal{A}^{A \to {\color{orange}R},2 \to 0,1,2}_{m_A,\lambda_0^{01},\bar\lambda_1^{01},\bar\lambda_2^{01}} \end{eqnarray} \end{split}\]
(2)#\[\begin{split} \begin{eqnarray} \mathcal{A}^{A \to {\color{violet}S},1 \to 0,1,2}_{m_A,m_0,m_1,m_2} &=& \sum_{\lambda_0^{02},\mu_0^{02},\nu_0^{02}} {\color{red}{D^{s_0}_{m_0,\nu_0^{02}}}}\!\left({\color{red}{\alpha_0^{02}, \beta_0^{02}, \gamma_0^{02}}}\right) {\color{red}{D^{s_0}_{\nu_0^{02},\mu_0^{02}}}}\!\left({\color{violet}{\phi_{_{02}}, \theta_{_{02}}}}, 0\right) {\color{red}{D^{s_0}_{\mu_0^{02},\lambda_0^{02}}}}\!\left({\color{red}{\phi_0^{02}, \theta_0^{02}}}\right) \\ &\times& \sum_{\lambda_1^{02}} {\color{green}{D^{s_1}_{m_1,\lambda_1^{02}}}}\!\left({\color{violet}{\phi_{_{02}}, \theta_{_{02}}}}, 0\right) \\ &\times& \sum_{\lambda_2^{02},\mu_2^{02},\nu_2^{02}} {\color{blue}{D^{s_2}_{m_2,\nu_2^{02}}}}\!\left({\color{blue}{\alpha_2^{02}, \beta_2^{02}, \gamma_2^{02}}}\right) {\color{blue}{D^{s_2}_{\nu_2^{02},\mu_2^{02}}}}\!\left({\color{violet}{\phi_{_{02}}, \theta_{_{02}}}}, 0\right) {\color{blue}{D^{s_2}_{\mu_2^{02},\lambda_2^{02}}}}\!\left({\color{red}{\phi_0^{02}, \theta_0^{02}}}\right) \\ &\times& \mathcal{A}^{A \to {\color{violet}S},2 \to 0,1,2}_{m_A,\lambda_0^{02},\bar\lambda_1^{02},\bar\lambda_2^{02}} \end{eqnarray} \end{split}\]

This procedure also allows us to formulate the alignment summation for \(\mathcal{A}^{A \to {\color{turquoise}U},0 \to 0,1,2}_{m_A,m_0,m_1,m_2}\):

Hide code cell source
dot3 = """
digraph {
    bgcolor=none
    rankdir=LR
    edge [arrowhead=none]
    node [shape=none, width=0]
    0 [shape=none, label=0, fontcolor=red]
    1 [shape=none, label=1, fontcolor=green]
    2 [shape=none, label=2, fontcolor=blue, label=<<o>2</o>>]
    A [shape=none, label=A]
    { rank=same A }
    { rank=same 0, 1, 2 }
    N0 [label=""]
    N1 [label=""]
    A -> N0 [style=dotted]
    N0 -> N1 [label=<U =<o>12</o>>, fontcolor=turquoise, style=dashed]
    N0 -> 0
    N1 -> 1
    N1 -> 2 [style=dashed]
}
"""
graphviz.Source(dot3)

(3)#\[\begin{split} \begin{eqnarray} \mathcal{A}^{A \to {\color{turquoise}U},0 \to 0,1,2}_{m_A,m_0,m_1,m_2} &=& \sum_{\lambda_0^{12}} {\color{red}{D^{s_0}_{m_0,\lambda_0^{12}}}}\!\left({\color{red}{\phi_0, \theta_0}}, 0\right) \\ &\times& \sum_{\lambda_1^{12},\mu_1^{12},\nu_1^{12}} {\color{green}{D^{s_0}_{m_0,\nu_1^{12}}}}\!\left({\color{green}{\alpha_1^{12}, \beta_1^{12}, \gamma_1^{12}}}\right) {\color{green}{D^{s_0}_{\nu_1^{12},\mu_1^{12}}}}\!\left({\color{red}{\phi_0, \theta_0}}, 0\right) {\color{green}{D^{s_0}_{\mu_1^{12},\lambda_1^{12}}}}\!\left({\color{green}{\phi_1^{12}, \theta_1^{12}}}\right) \\ &\times& \sum_{\lambda_2^{12},\mu_2^{12},\nu_2^{12}} {\color{blue}{D^{s_2}_{m_2,\nu_2^{12}}}}\!\left({\color{blue}{\alpha_2^{12}, \beta_2^{12}, \gamma_2^{12}}}\right) {\color{blue}{D^{s_2}_{\nu_2^{12},\mu_2^{12}}}}\!\left({\color{red}{\phi_0, \theta_0}}, 0\right) {\color{blue}{D^{s_2}_{\mu_2^{12},\lambda_2^{12}}}}\!\left({\color{green}{\phi_1^{12}, \theta_1^{12}}}\right) \\ &\times& \mathcal{A}^{A \to {\color{turquoise}U},2 \to 0,1,2}_{m_A,\lambda_1^{12},\bar\lambda_1^{12},\bar\lambda_2^{12}} \end{eqnarray} \end{split}\]

Finally, the total intensity can be computed from these amplitudes by incoherently summing over the initial and final state canonical spin projections (see Equation (47) in [Marangotto, 2020]):

(4)#\[ I = \sum_{m_A,m_0,m_1,m_2}\left| \mathcal{A}^{A \to {\color{orange}R},2 \to 0,1,2}_{m_A,m_0,m_1,m_2} + \mathcal{A}^{A \to {\color{violet}S},1 \to 0,1,2}_{m_A,m_0,m_1,m_2} + \mathcal{A}^{A \to {\color{turquoise}U},0 \to 0,1,2}_{m_A,m_0,m_1,m_2} \right|^2 \]

\(J/\psi \to K^0 \Sigma^+ \bar{p}\)#

Hide code cell content
def show_all_spin_matrices(transition, functor, cleanup: bool) -> None:
    for i in transition.final_states:
        state = transition.states[i]
        particle_name = state.particle.latex
        s = sp.Rational(state.particle.spin)
        m = sp.Rational(state.spin_projection)
        display(Math(Rf"|s_{i},m_{i}\rangle=|{s},{m}\rangle \quad ({particle_name})"))
        if functor is formulate_rotation_chain:
            args = (transition, i)
        else:
            args = (transition, i, state.spin_projection)
        summation = functor(*args)
        if cleanup:
            summation = summation.cleanup()
        display(summation)

In this section, we test some of the functions from the helicity and kinematics modules to see if they reproduce Equations (1), (2), and (3). We perform this test on the channel \(J/\psi \to K^0 \Sigma^+ \bar{p}\) with resonances generated for each of the three allowed three-body topologies. The transition that corresponds to Equation (1) is shown below.

The first step is to use {external+ampform-0.14.x:func}.formulate_helicity_rotation_chain to generate the Wigner-\(D\) functions for all helicity rotations for each final state. These helicity rotations “undo” all rotations that came from each Lorentz boosts when boosting from initial state \(J/\psi\) to each final state:

Hide code cell source
transition_r = full_reaction.transitions[-1]
show_all_spin_matrices(transition_r, formulate_helicity_rotation_chain, cleanup=True)
\[\displaystyle |s_0,m_0\rangle=|0,0\rangle \quad (K^{0})\]
\[\displaystyle D^{0}_{0,0}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{0}_{\nu^{01}_{0},0}\left(\phi_{01},\theta_{01},0\right)\]
\[\displaystyle |s_1,m_1\rangle=|1/2,1/2\rangle \quad (\Sigma^{+})\]
\[\displaystyle \sum_{\lambda^{01}_{1}=-1/2}^{1/2} \sum_{\mu^{01}_{1}=-1/2}^{1/2}{D^{\frac{1}{2}}_{\mu^{01}_{1},\lambda^{01}_{1}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{\frac{1}{2}}_{\nu^{01}_{1},\mu^{01}_{1}}\left(\phi_{01},\theta_{01},0\right)}\]
\[\displaystyle |s_2,m_2\rangle=|1/2,1/2\rangle \quad (\overline{p})\]
\[\displaystyle \sum_{\lambda^{01}_{2}=-1/2}^{1/2}{D^{\frac{1}{2}}_{0.5,\lambda^{01}_{2}}\left(\phi_{01},\theta_{01},0\right)}\]
Hide code cell source
show_transition(transition_r)

The function {external+ampform-0.14.x:func}.formulate_rotation_chain goes one step further. It adds a Wigner rotation to the generated list of helicity rotation Wigner-\(D\) functions in case there are resonances in between the initial state and rotated final state. If there are no resonances in between (here, state 2, the \(\bar p\)), there is only one helicity rotation and there is no need for a Wigner rotation.

Hide code cell source
show_all_spin_matrices(transition_r, formulate_rotation_chain, cleanup=False)
\[\displaystyle |s_0,m_0\rangle=|0,0\rangle \quad (K^{0})\]
\[\displaystyle \sum_{\lambda^{01}_{0}=0} \sum_{\mu^{01}_{0}=0} \sum_{\nu^{01}_{0}=0}{D^{0}_{m_{0},\nu^{01}_{0}}\left(\alpha^{01}_{0},\beta^{01}_{0},\gamma^{01}_{0}\right) D^{0}_{\mu^{01}_{0},\lambda^{01}_{0}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{0}_{\nu^{01}_{0},\mu^{01}_{0}}\left(\phi_{01},\theta_{01},0\right)}\]
\[\displaystyle |s_1,m_1\rangle=|1/2,1/2\rangle \quad (\Sigma^{+})\]
\[\displaystyle \sum_{\lambda^{01}_{1}=-1/2}^{1/2} \sum_{\mu^{01}_{1}=-1/2}^{1/2} \sum_{\nu^{01}_{1}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{1},\nu^{01}_{1}}\left(\alpha^{01}_{1},\beta^{01}_{1},\gamma^{01}_{1}\right) D^{\frac{1}{2}}_{\mu^{01}_{1},\lambda^{01}_{1}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{\frac{1}{2}}_{\nu^{01}_{1},\mu^{01}_{1}}\left(\phi_{01},\theta_{01},0\right)}\]
\[\displaystyle |s_2,m_2\rangle=|1/2,1/2\rangle \quad (\overline{p})\]
\[\displaystyle \sum_{\lambda^{01}_{2}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{2},\lambda^{01}_{2}}\left(\phi_{01},\theta_{01},0\right)}\]

These are indeed all the terms that we see in Equation (1)!

To create all sum combinations for all final states, we can use {external+ampform-0.14.x:func}.formulate_spin_alignment. This should give the sum of Eq.(45):

alignment_summation = formulate_spin_alignment(transition_r)
alignment_summation.cleanup()
\[\displaystyle \sum_{\lambda^{01}_{1}=-1/2}^{1/2} \sum_{\mu^{01}_{1}=-1/2}^{1/2} \sum_{\nu^{01}_{1}=-1/2}^{1/2} \sum_{\lambda^{01}_{2}=-1/2}^{1/2}{D^{0}_{0,0}\left(\phi_{01},\theta_{01},0\right) D^{0}_{0,0}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{0}_{m_{0},0}\left(\alpha^{01}_{0},\beta^{01}_{0},\gamma^{01}_{0}\right) D^{\frac{1}{2}}_{m_{1},\nu^{01}_{1}}\left(\alpha^{01}_{1},\beta^{01}_{1},\gamma^{01}_{1}\right) D^{\frac{1}{2}}_{m_{2},\lambda^{01}_{2}}\left(\phi_{01},\theta_{01},0\right) D^{\frac{1}{2}}_{\mu^{01}_{1},\lambda^{01}_{1}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{\frac{1}{2}}_{\nu^{01}_{1},\mu^{01}_{1}}\left(\phi_{01},\theta_{01},0\right)}\]

Finally, here are the generated spin alignment terms for the other two decay chains. Notice that the first is indeed the same as (2):

Hide code cell source
reaction_s = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", ("Sigma+", [+0.5]), ("p~", [+0.5])],
    allowed_intermediate_particles=["N(1650)+"],
    allowed_interaction_types="strong",
    formalism="helicity",
)
transition_s = reaction_s.transitions[0]
show_all_spin_matrices(transition_s, formulate_rotation_chain, cleanup=False)
\[\displaystyle |s_0,m_0\rangle=|0,0\rangle \quad (K^{0})\]
\[\displaystyle \sum_{\lambda^{02}_{0}=0} \sum_{\mu^{02}_{0}=0} \sum_{\nu^{02}_{0}=0}{D^{0}_{m_{0},\nu^{02}_{0}}\left(\alpha^{02}_{0},\beta^{02}_{0},\gamma^{02}_{0}\right) D^{0}_{\mu^{02}_{0},\lambda^{02}_{0}}\left(\phi^{02}_{0},\theta^{02}_{0},0\right) D^{0}_{\nu^{02}_{0},\mu^{02}_{0}}\left(\phi_{02},\theta_{02},0\right)}\]
\[\displaystyle |s_1,m_1\rangle=|1/2,1/2\rangle \quad (\Sigma^{+})\]
\[\displaystyle \sum_{\lambda^{02}_{1}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{1},\lambda^{02}_{1}}\left(\phi_{02},\theta_{02},0\right)}\]
\[\displaystyle |s_2,m_2\rangle=|1/2,1/2\rangle \quad (\overline{p})\]
\[\displaystyle \sum_{\lambda^{02}_{2}=-1/2}^{1/2} \sum_{\mu^{02}_{2}=-1/2}^{1/2} \sum_{\nu^{02}_{2}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{2},\nu^{02}_{2}}\left(\alpha^{02}_{2},\beta^{02}_{2},\gamma^{02}_{2}\right) D^{\frac{1}{2}}_{\mu^{02}_{2},\lambda^{02}_{2}}\left(\phi^{02}_{0},\theta^{02}_{0},0\right) D^{\frac{1}{2}}_{\nu^{02}_{2},\mu^{02}_{2}}\left(\phi_{02},\theta_{02},0\right)}\]
Hide code cell source
show_transition(transition_s)

…and that the second matches Equation (3):

Hide code cell source
reaction_u = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", ("Sigma+", [+0.5]), ("p~", [+0.5])],
    allowed_intermediate_particles=["K*(1680)~0"],
    allowed_interaction_types="strong",
    formalism="helicity",
)
transition_u = reaction_u.transitions[0]
show_all_spin_matrices(transition_u, formulate_rotation_chain, cleanup=False)
\[\displaystyle |s_0,m_0\rangle=|0,0\rangle \quad (K^{0})\]
\[\displaystyle \sum_{\lambda^{12}_{0}=0}{D^{0}_{m_{0},\lambda^{12}_{0}}\left(\phi_{0},\theta_{0},0\right)}\]
\[\displaystyle |s_1,m_1\rangle=|1/2,1/2\rangle \quad (\Sigma^{+})\]
\[\displaystyle \sum_{\lambda^{12}_{1}=-1/2}^{1/2} \sum_{\mu^{12}_{1}=-1/2}^{1/2} \sum_{\nu^{12}_{1}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{1},\nu^{12}_{1}}\left(\alpha^{12}_{1},\beta^{12}_{1},\gamma^{12}_{1}\right) D^{\frac{1}{2}}_{\mu^{12}_{1},\lambda^{12}_{1}}\left(\phi^{12}_{1},\theta^{12}_{1},0\right) D^{\frac{1}{2}}_{\nu^{12}_{1},\mu^{12}_{1}}\left(\phi_{0},\theta_{0},0\right)}\]
\[\displaystyle |s_2,m_2\rangle=|1/2,1/2\rangle \quad (\overline{p})\]
\[\displaystyle \sum_{\lambda^{12}_{2}=-1/2}^{1/2} \sum_{\mu^{12}_{2}=-1/2}^{1/2} \sum_{\nu^{12}_{2}=-1/2}^{1/2}{D^{\frac{1}{2}}_{m_{2},\nu^{12}_{2}}\left(\alpha^{12}_{2},\beta^{12}_{2},\gamma^{12}_{2}\right) D^{\frac{1}{2}}_{\mu^{12}_{2},\lambda^{12}_{2}}\left(\phi^{12}_{1},\theta^{12}_{1},0\right) D^{\frac{1}{2}}_{\nu^{12}_{2},\mu^{12}_{2}}\left(\phi_{0},\theta_{0},0\right)}\]
Hide code cell source
show_transition(transition_u)

Compute Wigner rotation angles#

Now it’s still a matter of computing the values for the angles \(\alpha,\beta,\gamma\) in the Wigner rotation matrices. These angles represents the difference between the canonical spin frame as attained by a direct boost from the initial state versus a chain of boosts through each resonance. See Equation (36) in [Marangotto, 2020].

The kinematics module can generate an expression for the chain of Lorentz boosts from the initial state to the final state with compute_boost_chain():

Hide code cell source
dot = qrules.io.asdot(transition_u)
topology = transition_u.topology
display(graphviz.Source(dot))

momenta = create_four_momentum_symbols(topology)
for state_id in topology.outgoing_edge_ids:
    boosts = compute_boost_chain(topology, momenta, state_id)
    display(sp.Array(boosts))
\[\displaystyle \left[\begin{matrix}\boldsymbol{B}\left(p_{0}\right)\end{matrix}\right]\]
\[\displaystyle \left[\begin{matrix}\boldsymbol{B}\left({p}_{12}\right) & \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\end{matrix}\right]\]
\[\displaystyle \left[\begin{matrix}\boldsymbol{B}\left({p}_{12}\right) & \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\end{matrix}\right]\]

This chain of Lorentz boosts needs to be ‘undo’ with a direct Lorentz boost back to the initial state. A contraction of inverse Lorentz boost with the chain of Lorentz boosts can be generated with compute_wigner_rotation_matrix():

for state_id in topology.outgoing_edge_ids:
    expr = compute_wigner_rotation_matrix(topology, momenta, state_id)
    display(expr)
\[\displaystyle \boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\]
\[\displaystyle \boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\]
\[\displaystyle \boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\]

The result of this matrix product is the rotation matrix for the Wigner rotation. The function compute_wigner_angles() computes the required Euler angles from this rotation matrix by implementing Equations (B.2-3) from [Marangotto, 2020]:

angles = {}
for state_id in topology.outgoing_edge_ids:
    angle_definitions = compute_wigner_angles(topology, momenta, state_id)
    for name, expr in angle_definitions.items():
        angle_symbol = sp.Symbol(name, real=True)
        angles[angle_symbol] = expr
Hide code cell source
latex_lines = [R"\begin{eqnarray}"]
for symbol, expr in angles.items():
    latex_lines.append(Rf"{sp.latex(symbol)}&=&{sp.latex(expr)}\\")
latex_lines.append(R"\end{eqnarray}")
Math("\n".join(latex_lines))
\[\begin{split}\displaystyle \begin{eqnarray} \alpha^{12}_{0}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\left[:, 3, 2\right],\boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\left[:, 3, 1\right] \right)}\\ \beta^{12}_{0}&=&\operatorname{acos}{\left(\boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\left[:, 3, 3\right] \right)}\\ \gamma^{12}_{0}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\left[:, 2, 3\right],- \boldsymbol{B}\left(-\left(p_{0}\right)\right) \boldsymbol{B}\left(p_{0}\right)\left[:, 1, 3\right] \right)}\\ \alpha^{12}_{1}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\left[:, 3, 2\right],\boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\left[:, 3, 1\right] \right)}\\ \beta^{12}_{1}&=&\operatorname{acos}{\left(\boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\left[:, 3, 3\right] \right)}\\ \gamma^{12}_{1}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\left[:, 2, 3\right],- \boldsymbol{B}\left(-\left(p_{1}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{1}\right)\left[:, 1, 3\right] \right)}\\ \alpha^{12}_{2}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\left[:, 3, 2\right],\boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\left[:, 3, 1\right] \right)}\\ \beta^{12}_{2}&=&\operatorname{acos}{\left(\boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\left[:, 3, 3\right] \right)}\\ \gamma^{12}_{2}&=&\operatorname{atan_{2}}{\left(\boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\left[:, 2, 3\right],- \boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\left[:, 1, 3\right] \right)}\\ \end{eqnarray}\end{split}\]

Note

In the topology underlying (3), the Wigner rotation matrix with angles \(\alpha_0^{12}, \beta_0^{12}, \gamma_0^{12}\) is simply the identity matrix. This is the reason why it can be omitted in {external+ampform-0.14.x:func}.formulate_rotation_chain and we only have one helicity rotation.

Computational code#

Classes like BoostMatrix have been split up into smaller unevaluated() expression classes so that lambdification to NumPy code results in relatively small and fast code, when using cse=True in lambdify() (see NumPyPrintable).

Hide code cell source
beta = sp.Symbol("beta_1^12", real=True)
beta_expr = angles[beta]

func = sp.lambdify(momenta.values(), beta_expr.doit(), cse=True)
src = inspect.getsource(func)
n_characters = len(src)
latex = sp.latex(beta)
latex += Rf":\quad\text{{{n_characters:,} characters in generated code}}"
Math(latex)
\[\displaystyle \beta^{12}_{1}:\quad\text{2,147 characters in generated code}\]

Test on data sample#

Tip

A test with a larger data distribution is being developed in TR-013.

The following phase space mini-sample of four-momenta has been generated for the decay \(J/\psi \to K^0 \Sigma^+ \bar{p}\) with the tensorwaves.data module.

phsp = {
    "p0": np.array([
        [0.63140486, 0.13166435, -0.35734744, 0.07760603],
        [0.65169531, 0.37242432, 0.12027178, 0.15467675],
        [0.60647425, -0.22286205, -0.175258, 0.19952806],
        [0.72744323, 0.05529811, 0.30502402, -0.43064999],
        [0.76778868, -0.43557036, 0.35491651, -0.16185017],
    ]),
    "p1": np.array([
        [1.37017117, 0.173769668, 0.355893315, -0.553093198],
        [1.34556663, -5.272033e-04, -0.3074542, -0.54901747],
        [1.41660182, 0.634007973, -0.0457976658, -0.433700564],
        [1.38592340, 0.138369075, -0.258624859, 0.648189682],
        [1.37858847, 0.551405385, -0.338705615, 0.259105737],
    ]),
    "p2": np.array([
        [1.09532397, -0.30543402, 0.00145413, 0.47548716],
        [1.09963805, -0.37189712, 0.18718247, 0.39434072],
        [1.07382393, -0.41114592, 0.22105567, 0.2341725],
        [0.98353336, -0.19366719, -0.04639917, -0.21753969],
        [0.95052285, -0.11583502, -0.01621089, -0.09725557],
    ]),
}
matrix_expr = compute_wigner_rotation_matrix(topology, momenta, state_id=2)
matrix_expr
\[\displaystyle \boldsymbol{B}\left(-\left(p_{2}\right)\right) \boldsymbol{B}\left({p}_{12}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{12}\right) p_{2}\right)\]
matrix_func = sp.lambdify(momenta.values(), matrix_expr.doit(), cse=True)
matrix_array = matrix_func(*phsp.values())
np.round(matrix_array, decimals=2).real
array([[[ 1.  , -0.  , -0.  , -0.  ],
        [-0.  ,  1.  ,  0.02, -0.02],
        [-0.  , -0.02,  1.  ,  0.03],
        [ 0.  ,  0.02, -0.03,  1.  ]],

       [[ 1.  ,  0.  , -0.  ,  0.  ],
        [ 0.  ,  1.  , -0.02, -0.04],
        [ 0.  ,  0.02,  1.  , -0.  ],
        [ 0.  ,  0.04,  0.  ,  1.  ]],

       [[ 1.  ,  0.  ,  0.  ,  0.  ],
        [-0.  ,  1.  ,  0.02, -0.01],
        [ 0.  , -0.02,  1.  ,  0.02],
        [ 0.  ,  0.01, -0.02,  1.  ]],

       [[ 1.  ,  0.  , -0.  ,  0.  ],
        [ 0.  ,  1.  , -0.01,  0.02],
        [ 0.  ,  0.01,  1.  ,  0.02],
        [ 0.  , -0.02, -0.02,  1.  ]],

       [[ 1.  , -0.  ,  0.  , -0.  ],
        [-0.  ,  1.  , -0.01, -0.01],
        [ 0.  ,  0.01,  1.  ,  0.01],
        [-0.  ,  0.01, -0.01,  1.  ]]])
Hide code cell source
latex_lines = [R"\begin{eqnarray}"]
for angle_symbol, angle_expr in angles.items():
    angle_func = sp.lambdify(momenta.values(), angle_expr.doit(), cse=True)
    angle_array = angle_func(*phsp.values())
    rounded_values = np.round(angle_array, decimals=2).real
    latex_lines.append(
        Rf"{sp.latex(angle_symbol)}&=&{sp.latex(sp.Array(rounded_values))}\\"
    )
latex_lines.append(R"\end{eqnarray}")
Math("\n".join(latex_lines))
\[\begin{split}\displaystyle \begin{eqnarray} \alpha^{12}_{0}&=&\left[\begin{matrix}-1.25 & 0.0 & 0.79 & 1.33 & 2.36\end{matrix}\right]\\ \beta^{12}_{0}&=&\left[\begin{matrix}0.0 & \text{NaN} & 0.0 & 0.0 & 0.0\end{matrix}\right]\\ \gamma^{12}_{0}&=&\left[\begin{matrix}-1.89 & 3.14 & 2.36 & 1.82 & 0.79\end{matrix}\right]\\ \alpha^{12}_{1}&=&\left[\begin{matrix}2.03 & -3.04 & 1.9 & 0.74 & 2.14\end{matrix}\right]\\ \beta^{12}_{1}&=&\left[\begin{matrix}0.03 & 0.03 & 0.01 & 0.02 & 0.01\end{matrix}\right]\\ \gamma^{12}_{1}&=&\left[\begin{matrix}-2.05 & 3.06 & -1.92 & -0.73 & -2.13\end{matrix}\right]\\ \alpha^{12}_{2}&=&\left[\begin{matrix}-1.09 & 0.08 & -1.22 & -2.41 & -1.01\end{matrix}\right]\\ \beta^{12}_{2}&=&\left[\begin{matrix}0.04 & 0.04 & 0.02 & 0.03 & 0.01\end{matrix}\right]\\ \gamma^{12}_{2}&=&\left[\begin{matrix}1.11 & -0.1 & 1.25 & 2.4 & 1.0\end{matrix}\right]\\ \end{eqnarray}\end{split}\]
Hide code cell source
dot = qrules.io.asdot(transition_u, collapse_graphs=True)
graphviz.Source(dot)

Note

The nan values above come from the fact that the inverse boost on a boost results in negative values under the square root of \(\gamma=\sqrt{1-\beta^2}\). This can be ignored, because the Wigner rotation is simply omitted when formulating the chain of rotation matrices, as noted in Compute Wigner rotation angles.

Four-body decay#

The algorithm for computing Euler angles for the Wigner rotation works an arbitrary number of final states. Here, we illustrate this by formulating an expression for the Wigner rotation matrix in a four-body decay.

topology_4body = create_isobar_topologies(4)[1]
momenta_4body = create_four_momentum_symbols(topology_4body)
compute_wigner_rotation_matrix(topology_4body, momenta_4body, state_id=3)
\[\displaystyle \boldsymbol{B}\left(-\left(p_{3}\right)\right) \boldsymbol{B}\left({p}_{123}\right) \boldsymbol{B}\left(\boldsymbol{B}\left({p}_{123}\right) {p}_{23}\right) \boldsymbol{B}\left(\boldsymbol{B}\left(\boldsymbol{B}\left({p}_{123}\right) {p}_{23}\right) \boldsymbol{B}\left({p}_{123}\right) p_{3}\right)\]
Hide code cell source
dot = qrules.io.asdot(topology_4body)
graphviz.Source(dot)