Skip to content

Compiler

oqd_trical.light_matter.compiler

FirstOrderLambDickeApprox

Bases: RewriteRule

Applies the Lamb-Dicke approximation to first order.

Attributes:

Name Type Description
cutoff float

Lamb-Dicke parameter cutoff below which approximation is applied.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class FirstOrderLambDickeApprox(RewriteRule):
    """
    Applies the Lamb-Dicke approximation to first order.

    Attributes:
        cutoff (float): Lamb-Dicke parameter cutoff below which approximation is applied.
    """

    def __init__(self, cutoff=1):
        super().__init__()
        self.cutoff = cutoff

        self.approximated_operators = []

    def map_Displacement(self, model):
        if isinstance(model.alpha.amplitude, MathNum):
            if np.abs(model.alpha.amplitude.value) < self.cutoff:
                self.approximated_operators.append(model)

                alpha_conj = model.alpha.conj()
                return Identity(subsystem=model.subsystem) + (
                    model.alpha * Creation(subsystem=model.subsystem)
                    - alpha_conj * Annihilation(subsystem=model.subsystem)
                )

SecondOrderLambDickeApprox

Bases: RewriteRule

Applies the Lamb-Dicke approximation to second order.

Attributes:

Name Type Description
cutoff float

Lamb-Dicke parameter cutoff below which approximation is applied.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class SecondOrderLambDickeApprox(RewriteRule):
    """
    Applies the Lamb-Dicke approximation to second order.

    Attributes:
        cutoff (float): Lamb-Dicke parameter cutoff below which approximation is applied.
    """

    def __init__(self, cutoff=1):
        super().__init__()

        self.cutoff = cutoff

        self.approximated_operators = []

    def map_Displacement(self, model):
        if isinstance(model.alpha.amplitude, MathNum):
            if np.abs(model.alpha.amplitude.value) < self.cutoff:
                self.approximated_operators.append(model)

                alpha_conj = model.alpha.conj()
                return (
                    Identity(subsystem=model.subsystem)
                    + (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                    + ConstantCoefficient(value=1 / 2)
                    * (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                    * (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                )

ConstructHamiltonian

Bases: ConversionRule

Maps an AtomicCircuit to an AtomicEmulatorCircuit replaces laser descriptions of operations with Hamiltonian description of operations

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/codegen.py
class ConstructHamiltonian(ConversionRule):
    """Maps an AtomicCircuit to an AtomicEmulatorCircuit replaces laser descriptions of operations with Hamiltonian description of operations"""

    def map_AtomicCircuit(self, model, operands):
        gates = (
            [operands["protocol"]]
            if isinstance(operands["protocol"], AtomicEmulatorGate)
            else operands["protocol"]
        )
        for gate in gates:
            gate.hamiltonian = gate.hamiltonian + operands["system"]

        return AtomicEmulatorCircuit(sequence=gates)

    def map_System(self, model, operands):
        self.N = len(model.ions)
        self.M = len(model.modes)

        self.ions = model.ions
        self.modes = model.modes

        ops = []
        for n, ion in enumerate(model.ions):
            ops.append(
                reduce(
                    lambda x, y: x @ y,
                    [
                        (
                            self._map_Ion(ion, n)
                            if i == n
                            else Identity(
                                subsystem=f"E{i}" if i < self.N else f"P{i - self.N}"
                            )
                        )
                        for i in range(self.N + self.M)
                    ],
                )
            )

        for m, mode in enumerate(model.modes):
            ops.append(
                reduce(
                    lambda x, y: x @ y,
                    [
                        (
                            self._map_Phonon(mode, m)
                            if i == self.N + m
                            else Identity(
                                subsystem=f"E{i}" if i < self.N else f"P{i - self.N}"
                            )
                        )
                        for i in range(self.N + self.M)
                    ],
                )
            )

        op = reduce(lambda x, y: x + y, ops)
        return op

    def _map_Ion(self, model, index):
        ops = [
            WaveCoefficient(amplitude=level.energy, frequency=0, phase=0)
            * KetBra(ket=n, bra=n, subsystem=f"E{index}")
            for n, level in enumerate(model.levels)
        ]

        op = reduce(lambda x, y: x + y, ops)
        return op

    def _map_Phonon(self, model, index):
        return WaveCoefficient(amplitude=model.energy, frequency=0, phase=0) * (
            Creation(subsystem=f"P{index}") * Annihilation(subsystem=f"P{index}")
        )

    def map_Beam(self, model, operands):
        I = intensity_from_laser(model)  # noqa: E741

        angular_frequency = (
            abs(model.transition.level2.energy - model.transition.level1.energy)
            + model.detuning
        )
        wavevector = angular_frequency * np.array(model.wavevector) / cst.c

        ops = []
        if self.modes:
            displacement_plus = []
            displacement_minus = []
            for m, mode in enumerate(self.modes):
                eta = np.dot(
                    wavevector,
                    mode.eigenvector[model.target * 3 : model.target * 3 + 3],
                ) * np.sqrt(
                    cst.hbar
                    / (2 * self.ions[model.target].mass * cst.m_u * mode.energy)
                )

                displacement_plus.append(
                    Displacement(
                        alpha=WaveCoefficient(
                            amplitude=eta, frequency=0, phase=np.pi / 2
                        ),
                        subsystem=f"P{m}",
                    )
                )
                displacement_minus.append(
                    Displacement(
                        alpha=WaveCoefficient(
                            amplitude=eta, frequency=0, phase=-np.pi / 2
                        ),
                        subsystem=f"P{m}",
                    )
                )

            displacement_plus = reduce(lambda x, y: x @ y, displacement_plus)
            displacement_minus = reduce(lambda x, y: x @ y, displacement_minus)

            for transition in self.ions[model.target].transitions:
                rabi = rabi_from_intensity(model, transition, I)

                ops.append(
                    (
                        reduce(
                            lambda x, y: x @ y,
                            [
                                (
                                    WaveCoefficient(
                                        amplitude=rabi / 2,
                                        frequency=-angular_frequency,
                                        phase=model.phase,
                                    )
                                    * (
                                        KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                        + KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                    )
                                    if i == model.target
                                    else Identity(subsystem=f"E{i}")
                                )
                                for i in range(self.N)
                            ],
                        )
                        @ displacement_plus
                    )
                )

                ops.append(
                    (
                        reduce(
                            lambda x, y: x @ y,
                            [
                                (
                                    WaveCoefficient(
                                        amplitude=rabi / 2,
                                        frequency=angular_frequency,
                                        phase=-model.phase,
                                    )
                                    * (
                                        KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                        + KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                    )
                                    if i == model.target
                                    else Identity(subsystem=f"E{i}")
                                )
                                for i in range(self.N)
                            ],
                        )
                        @ displacement_minus
                    )
                )

        else:
            for transition in self.ions[model.target].transitions:
                rabi = rabi_from_intensity(model, transition, I)

                ops.append(
                    reduce(
                        lambda x, y: x @ y,
                        [
                            (
                                WaveCoefficient(
                                    amplitude=rabi / 2,
                                    frequency=angular_frequency,
                                    phase=model.phase,
                                )
                                * (
                                    KetBra(
                                        ket=self.ions[model.target].levels.index(
                                            transition.level1
                                        ),
                                        bra=self.ions[model.target].levels.index(
                                            transition.level2
                                        ),
                                        subsystem=f"E{model.target}",
                                    )
                                    + KetBra(
                                        ket=self.ions[model.target].levels.index(
                                            transition.level2
                                        ),
                                        bra=self.ions[model.target].levels.index(
                                            transition.level1
                                        ),
                                        subsystem=f"E{model.target}",
                                    )
                                )
                                if i == model.target
                                else Identity(subsystem=f"E{i}")
                            )
                            for i in range(self.N)
                        ],
                    )
                )

        op = reduce(lambda x, y: x + y, ops)
        return op

    def map_Pulse(self, model, operands):
        return AtomicEmulatorGate(
            hamiltonian=operands["beam"],
            duration=model.duration,
        )

    def map_ParallelProtocol(self, model, operands):
        duration = operands["sequence"][0].duration

        for op in operands["sequence"]:
            assert op.duration == duration

        op = reduce(
            lambda x, y: x + y, map(lambda x: x.hamiltonian, operands["sequence"])
        )
        return AtomicEmulatorGate(hamiltonian=op, duration=duration)

    def map_SequentialProtocol(self, model, operands):
        return operands["sequence"]

CoefficientPrinter

Bases: ConversionRule

Prints Coefficients in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class CoefficientPrinter(ConversionRule):
    """Prints Coefficients in a pretty manner"""

    def map_MathExpr(self, model, operands):
        return Post(PrintMathExpr())(model)

    def map_MathAdd(self, model, operands):
        return f"({Post(PrintMathExpr())(model)})"

    def map_WaveCoefficient(self, model, operands):
        frequency_term = (
            f"{operands['frequency']} * t"
            if model.frequency != MathNum(value=0)
            else ""
        )
        phase_term = f"{operands['phase']}" if model.phase != MathNum(value=0) else ""

        wave_term = (
            f" * exp(1j * ({frequency_term}{' + ' if frequency_term and phase_term else ''}{phase_term}))"
            if phase_term or frequency_term
            else ""
        )

        return f"{operands['amplitude']}{wave_term}"

    def map_CoefficientAdd(self, model, operands):
        return f"{'(' + operands['coeff1'] + ')' if isinstance(model.coeff1, CoefficientAdd) else operands['coeff1']} + {'(' + operands['coeff2'] + ')' if isinstance(model.coeff2, CoefficientAdd) else operands['coeff2']}"

    def map_CoefficientMul(self, model, operands):
        return f"{'(' + operands['coeff1'] + ')' if isinstance(model.coeff1, CoefficientAdd) else operands['coeff1']} * {'(' + operands['coeff2'] + ')' if isinstance(model.coeff2, CoefficientAdd) else operands['coeff2']}"

CondensedOperatorPrettyPrint

Bases: PrettyPrint

Prints An AtomicEmulatorCircuit in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class CondensedOperatorPrettyPrint(PrettyPrint):
    """Prints An AtomicEmulatorCircuit in a pretty manner"""

    def map_Operator(self, model, operands):
        return f"Operator({Post(OperatorPrinter())(model)})"

OperatorPrinter

Bases: ConversionRule

Prints Operators in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class OperatorPrinter(ConversionRule):
    """Prints Operators in a pretty manner"""

    def map_KetBra(self, model, operands):
        return f"|{model.ket}><{model.bra}|_{model.subsystem}"

    def map_Annihilation(self, model, operands):
        return f"A_{model.subsystem}"

    def map_Creation(self, model, operands):
        return f"C_{model.subsystem}"

    def map_Identity(self, model, operands):
        return f"I_{model.subsystem}"

    def map_PrunedOperator(self, model, operands):
        return "PrunedOperator"

    def map_Displacement(self, model, operands):
        return f"D({operands['alpha']})_{model.subsystem}"

    def map_OperatorAdd(self, model, operands):
        return f"{operands['op1']} + {operands['op2']}"

    def map_OperatorMul(self, model, operands):
        return f"{'(' + operands['op1'] + ')' if isinstance(model.op1, OperatorAdd) else operands['op1']} * {'(' + operands['op2'] + ')' if isinstance(model.op2, OperatorAdd) else operands['op2']}"

    def map_OperatorKron(self, model, operands):
        return f"{'(' + operands['op1'] + ')' if isinstance(model.op1, OperatorAdd) else operands['op1']} @ {'(' + operands['op2'] + ')' if isinstance(model.op2, OperatorAdd) else operands['op2']}"

    def map_OperatorScalarMul(self, model, operands):
        return f"{'(' + operands['coeff'] + ')' if isinstance(model.coeff, CoefficientAdd) else operands['coeff']} * {'(' + operands['op'] + ')' if isinstance(model.op, OperatorAdd) else operands['op']}"

    def map_Coefficient(self, model, operands):
        return Post(CoefficientPrinter())(model)

canonicalize_emulator_circuit_factory()

Creates a new instance of the canonicalization pass for AtomicEmulatorCircuit

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
def canonicalize_emulator_circuit_factory():
    """Creates a new instance of the canonicalization pass for AtomicEmulatorCircuit"""
    return Chain(
        canonicalize_operator_factory(),
        canonicalize_coefficient_factory(),
        canonicalize_math_factory(),
        Post(PruneOperator()),
        Pre(ScaleTerms()),
        Post(CombineTerms()),
        canonicalize_coefficient_factory(),
        canonicalize_math_factory(),
        Post(PruneOperator()),
    )

compute_matrix_element(laser, transition)

Computes matrix element corresponding to a laser interacting with a particular transition

Parameters:

Name Type Description Default
laser Beam

laser to compute matrix element of

required
transition Transition

transition to compute matrix element of

required

Returns:

Name Type Description
matrix_element float

Multipole matrix elements corresponding to the interaction between the laser and the transition

Warning

Currently implemented for only E1 and E2 transitions.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def compute_matrix_element(laser, transition):
    """Computes matrix element corresponding to a laser interacting with a particular transition

    Args:
        laser (Beam): laser to compute matrix element of
        transition (Transition): transition to compute matrix element of

    Returns:
        matrix_element (float): Multipole matrix elements corresponding to the interaction between the laser and the transition

    Warning:
        Currently implemented for only `E1` and `E2` transitions.
    """

    # If this happens there's probably an error with the ion species card
    if transition.level1.nuclear != transition.level2.nuclear:
        raise ValueError(
            "Different nuclear spins between two levels in transition:", transition
        )

    # Just by convention; these orderings are set upon instantiation of an Ion object
    if transition.level1.energy > transition.level2.energy:
        raise ValueError("Expected energy of level2 > energy of level1")

    polarization = np.array(laser.polarization).T  # make it a column vector
    wavevector = laser.wavevector

    J1, J2, F1, F2, M1, M2, E1, E2, I, A = (  # noqa: E741
        transition.level1.spin_orbital,
        transition.level2.spin_orbital,
        transition.level1.spin_orbital_nuclear,
        transition.level2.spin_orbital_nuclear,
        transition.level1.spin_orbital_nuclear_magnetization,
        transition.level2.spin_orbital_nuclear_magnetization,
        transition.level1.energy,
        transition.level2.energy,
        transition.level1.nuclear,
        transition.einsteinA,
    )

    q = M2 - M1

    omega = E2 - E1

    # --- M1 transition multipole ---
    if transition.multipole == "M1":
        # The Einstein A coefficient for an M1 transition is given by:
        #   A(M1) = (16 * pi^3)/(3 * hbar) * (omega^3/c^3) * (|<J2||μ||J1>|^2/(2J2+1))

        # Solving for the reduced matrix element (expressed in units of the Bohr magneton, μ_B) :
        #   |<J2||μ||J1>|/μ_B = sqrt((3 * hbar * c^3 * A * (2J2+1))/(16 * pi^3 * omega^3)) / μ_B

        # Reference: I. I. Sobelman, "Atomic Spectra and Radiative Transitions", 2nd ed., Springer (1992).

        # A unit term is definied so that when multiplied by standard angular momentum factors,  the full matrix is reproduced

        units_term = np.sqrt(
            (3 * cst.hbar * cst.epsilon_0 * cst.c**5 * A) / (16 * np.pi**3 * omega**3)
        )

        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 1, F2, F1, I
        )

        # For M1 transitions the operator is a vector operator coupling to the magnetic field
        # Plane wave: magnetic field is proportional to cross product of wavevector and electric field polarization
        B_field = np.cross(wavevector, polarization)

        # Define a spherical basis for a vector (identical to the one used for E1):
        polarization_map = {
            -1: 1 / np.sqrt(2) * np.array([1, 1j, 0]),
            0: np.array([0, 0, 1]),
            1: 1 / np.sqrt(2) * np.array([1, -1j, 0]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * polarization_map[q].dot(B_field)
            * wigner_3j(F2, 1, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    # --- E1 transition multipole ---
    if transition.multipole == "E1":
        units_term = np.sqrt(
            (3 * np.pi * cst.epsilon_0 * cst.hbar * cst.c**3 * A)
            / (omega**3 * cst.e**2)
        )
        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 1, F2, F1, I
        )

        # q -> polarization
        polarization_map = {
            -1: 1 / np.sqrt(2) * np.array([1, 1j, 0]),
            0: np.array([0, 0, 1]),
            1: 1 / np.sqrt(2) * np.array([1, -1j, 0]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * polarization_map[q].dot(polarization)
            * wigner_3j(F2, 1, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    # --- E2 transition multipole ---
    elif transition.multipole == "E2":
        units_term = np.sqrt(
            (15 * np.pi * cst.epsilon_0 * cst.hbar * cst.c**3 * A)
            / (omega**3 * cst.e**2)
        )  # <- anomalous constants I needed to add... hmm
        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 2, F2, F1, I
        )

        # q -> polarization
        polarization_map = {
            -2: 1 / np.sqrt(6) * np.array([[1, 1j, 0], [1j, -1, 0], [0, 0, 0]]),
            -1: 1 / np.sqrt(6) * np.array([[0, 0, 1], [0, 0, 1j], [1, 1j, 0]]),
            0: 1 / 3 * np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 2]]),
            1: 1 / np.sqrt(6) * np.array([[0, 0, -1], [0, 0, 1j], [-1, 1j, 0]]),
            2: 1 / np.sqrt(6) * np.array([[1, -1j, 0], [-1j, -1, 0], [0, 0, 0]]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * wavevector.dot(np.matmul(polarization_map[q], polarization))
            * wigner_3j(F2, 2, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    else:
        raise ValueError(
            "Currently only support dipole and quadrupole allowed transitions"
        )

intensity_from_laser(laser)

Computes the intensity from a laser

Parameters:

Name Type Description Default
laser Beam

laser to compute intensity of.

required

Returns:

Name Type Description
intensity float

intensity of the laser

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def intensity_from_laser(laser):
    """Computes the intensity from a laser

    Args:
        laser (Beam): laser to compute intensity of.

    Returns:
        intensity (float): intensity of the laser
    """
    matrix_elem = compute_matrix_element(laser, laser.transition)

    if laser.transition.multipole[0] == "E":
        return (
            cst.c
            * cst.epsilon_0
            / 2
            * (cst.hbar * laser.rabi / (matrix_elem * cst.e)) ** 2
        )  # noqa: E741

    if laser.transition.multipole[0] == "M":
        return cst.c**3 * cst.epsilon_0 / 2 * (cst.hbar * laser.rabi / matrix_elem) ** 2  # noqa: E741

rabi_from_intensity(laser, transition, intensity)

Computes a transition's resonant rabi frequency when addressed by a laser with intensity

Parameters:

Name Type Description Default
laser Beam

laser to compute resonant rabi frequency of

required
transition Transition

transition to compute resonant rabi frequency of

required
intensity float

intensity of laser

required

Returns:

Name Type Description
rabi_frequency float

resonant rabi frequency corresponding to the interaction between the laser and transition

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def rabi_from_intensity(laser, transition, intensity):
    """Computes a transition's resonant rabi frequency when addressed by a laser with intensity

    Args:
        laser (Beam): laser to compute resonant rabi frequency of
        transition (Transition): transition to compute resonant rabi frequency of
        intensity (float): intensity of laser

    Returns:
        rabi_frequency (float): resonant rabi frequency corresponding to the interaction between the laser and transition
    """

    matrix_elem = compute_matrix_element(laser, transition)

    if transition.multipole[0] == "E":
        E = (2 * intensity / (cst.epsilon_0 * cst.c)) ** 0.5
        return matrix_elem * E * cst.e / cst.hbar
    if transition.multipole[0] == "M":
        B = (2 * intensity / (cst.epsilon_0 * cst.c**3)) ** 0.5
        return matrix_elem * B / cst.hbar

approximate

FirstOrderLambDickeApprox

Bases: RewriteRule

Applies the Lamb-Dicke approximation to first order.

Attributes:

Name Type Description
cutoff float

Lamb-Dicke parameter cutoff below which approximation is applied.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class FirstOrderLambDickeApprox(RewriteRule):
    """
    Applies the Lamb-Dicke approximation to first order.

    Attributes:
        cutoff (float): Lamb-Dicke parameter cutoff below which approximation is applied.
    """

    def __init__(self, cutoff=1):
        super().__init__()
        self.cutoff = cutoff

        self.approximated_operators = []

    def map_Displacement(self, model):
        if isinstance(model.alpha.amplitude, MathNum):
            if np.abs(model.alpha.amplitude.value) < self.cutoff:
                self.approximated_operators.append(model)

                alpha_conj = model.alpha.conj()
                return Identity(subsystem=model.subsystem) + (
                    model.alpha * Creation(subsystem=model.subsystem)
                    - alpha_conj * Annihilation(subsystem=model.subsystem)
                )

SecondOrderLambDickeApprox

Bases: RewriteRule

Applies the Lamb-Dicke approximation to second order.

Attributes:

Name Type Description
cutoff float

Lamb-Dicke parameter cutoff below which approximation is applied.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class SecondOrderLambDickeApprox(RewriteRule):
    """
    Applies the Lamb-Dicke approximation to second order.

    Attributes:
        cutoff (float): Lamb-Dicke parameter cutoff below which approximation is applied.
    """

    def __init__(self, cutoff=1):
        super().__init__()

        self.cutoff = cutoff

        self.approximated_operators = []

    def map_Displacement(self, model):
        if isinstance(model.alpha.amplitude, MathNum):
            if np.abs(model.alpha.amplitude.value) < self.cutoff:
                self.approximated_operators.append(model)

                alpha_conj = model.alpha.conj()
                return (
                    Identity(subsystem=model.subsystem)
                    + (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                    + ConstantCoefficient(value=1 / 2)
                    * (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                    * (
                        model.alpha * Creation(subsystem=model.subsystem)
                        - alpha_conj * Annihilation(subsystem=model.subsystem)
                    )
                )

RotatingWaveApprox

Bases: RewriteRule

Applies the rotating wave approximation.

Attributes:

Name Type Description
cutoff float

Frequency cutoff above which approximation is applied.

Warning

Currently not implmented!

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class RotatingWaveApprox(RewriteRule):
    """
    Applies the rotating wave approximation.

    Attributes:
        cutoff (float): Frequency cutoff above which approximation is applied.

    Warning:
        Currently not implmented!
    """

    def __init__(self, cutoff):
        super().__init__()

        self.cutoff = cutoff
        self.current_time = 0

    def map_AtomicEmulatorGate(self, model):
        hamiltonian = Post(
            _RotatingWaveApproxHelper(
                cutoff=self.cutoff,
                start_time=self.current_time,
                duration=model.duration,
            )
        )(model.hamiltonian)

        self.current_time += model.duration
        return model.__class__(hamiltonian=hamiltonian, duration=model.duration)

RotatingReferenceFrame

Bases: RewriteRule

Moves to an interaction picture with a rotating frame of reference.

Attributes:

Name Type Description
frame Operator

Operator that defines the rotating frame of reference.

Warning

Currently not implmented!

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/approximate.py
class RotatingReferenceFrame(RewriteRule):
    """
    Moves to an interaction picture with a rotating frame of reference.

    Attributes:
        frame (Operator): [`Operator`][oqd_trical.light_matter.interface.operator.Operator] that defines the rotating frame of reference.

    Warning:
        Currently not implmented!
    """

    def __init__(self, frame_specs):
        super().__init__()

        self.frame_specs = frame_specs

    @cached_property
    def system(self):
        return list(self.frame_specs.keys())

    def _complete_operator(self, op, subsystem):
        return reduce(
            lambda x, y: x @ y,
            [op if subsystem == s else Identity(subsystem=s) for s in self.system],
        )

    @cached_property
    def frame(self):
        ops = []
        for subsystem, energy in self.frame_specs.items():
            if subsystem[0] == "E":
                ops.extend(
                    [
                        ConstantCoefficient(value=e)
                        * self._complete_operator(
                            KetBra(ket=n, bra=n, subsystem=subsystem), subsystem
                        )
                        for n, e in enumerate(energy)
                    ]
                )
            else:
                ops.append(
                    ConstantCoefficient(value=energy)
                    * self._complete_operator(
                        Creation(subsystem=subsystem)
                        * Annihilation(subsystem=subsystem),
                        subsystem,
                    )
                )

        return reduce(lambda x, y: x + y, ops)

    def map_AtomicEmulatorCircuit(self, model):
        return model.__class__(frame=self.frame, sequence=model.sequence)

    def map_AtomicEmulatorGate(self, model):
        return model.__class__(
            hamiltonian=model.hamiltonian - self.frame, duration=model.duration
        )

    def map_KetBra(self, model):
        return (
            WaveCoefficient(
                amplitude=1,
                frequency=self.frame_specs[model.subsystem][model.ket]
                - self.frame_specs[model.subsystem][model.bra],
                phase=0,
            )
            * model
        )

    def map_Displacement(self, model):
        alpha = model.alpha
        return model.__class__(
            alpha=alpha.__class__(
                amplitude=alpha.amplitude,
                frequency=alpha.frequency + self.frame_specs[model.subsystem],
                phase=alpha.phase,
            ),
            subsystem=model.subsystem,
        )

    def map_Annihilation(self, model):
        return (
            WaveCoefficient(
                amplitude=1, frequency=self.frame_specs[model.subsystem], phase=0
            )
            * model
        )

    def map_Creation(self, model):
        return (
            WaveCoefficient(
                amplitude=1, frequency=-self.frame_specs[model.subsystem], phase=0
            )
            * model
        )

canonicalize

PruneOperator

Bases: RewriteRule

Prunes an Operator AST by removing zeros

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class PruneOperator(RewriteRule):
    """Prunes an Operator AST by removing zeros"""

    def map_OperatorAdd(self, model):
        if isinstance(model.op1, PrunedOperator):
            return model.op2
        if isinstance(model.op2, PrunedOperator):
            return model.op1

    def map_OperatorMul(self, model):
        if isinstance(model.op1, PrunedOperator) or isinstance(
            model.op2, PrunedOperator
        ):
            return PrunedOperator()

    def map_OperatorScalarMul(self, model):
        if isinstance(
            model.coeff, WaveCoefficient
        ) and model.coeff.amplitude == MathNum(value=0):
            return PrunedOperator()
        if isinstance(model.op, PrunedOperator):
            return PrunedOperator()

    def map_OperatorKron(self, model):
        if isinstance(model.op1, PrunedOperator) or isinstance(
            model.op2, PrunedOperator
        ):
            return PrunedOperator()

    def map_Displacement(self, model):
        if isinstance(
            model.alpha, WaveCoefficient
        ) and model.alpha.amplitude == MathNum(value=0):
            return Identity(subsystem=model.subsystem)

PruneZeroPowers

Bases: RewriteRule

Prunes a MathExpr AST by MathPow when base is zero

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class PruneZeroPowers(RewriteRule):
    """Prunes a MathExpr AST by MathPow when base is zero"""

    def map_MathPow(self, model):
        if model.expr1 == MathNum(value=0):
            return MathNum(value=0)

OperatorDistributivity

Bases: RewriteRule

Implements distributivity of addition over multiplication, kronecker product and scalar multiplication

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class OperatorDistributivity(RewriteRule):
    """Implements distributivity of addition over multiplication, kronecker product and scalar multiplication"""

    def map_OperatorMul(self, model):
        if isinstance(model.op1, (OperatorAdd)):
            return model.op1.__class__(
                op1=OperatorMul(op1=model.op1.op1, op2=model.op2),
                op2=OperatorMul(op1=model.op1.op2, op2=model.op2),
            )
        if isinstance(model.op2, (OperatorAdd)):
            return model.op2.__class__(
                op1=OperatorMul(op1=model.op1, op2=model.op2.op1),
                op2=OperatorMul(op1=model.op1, op2=model.op2.op2),
            )
        if isinstance(model.op1, (OperatorKron)) and isinstance(
            model.op2, (OperatorKron)
        ):
            return OperatorKron(
                op1=OperatorMul(op1=model.op1.op1, op2=model.op2.op1),
                op2=OperatorMul(op1=model.op1.op2, op2=model.op2.op2),
            )
        return None

    def map_OperatorKron(self, model):
        if isinstance(model.op1, (OperatorAdd)):
            return model.op1.__class__(
                op1=OperatorKron(op1=model.op1.op1, op2=model.op2),
                op2=OperatorKron(op1=model.op1.op2, op2=model.op2),
            )
        if isinstance(model.op2, (OperatorAdd)):
            return model.op2.__class__(
                op1=OperatorKron(op1=model.op1, op2=model.op2.op1),
                op2=OperatorKron(op1=model.op1, op2=model.op2.op2),
            )
        return None

    def map_OperatorScalarMul(self, model):
        if isinstance(model.op, (OperatorAdd)):
            return model.op.__class__(
                op1=OperatorScalarMul(op=model.op.op1, coeff=model.coeff),
                op2=OperatorScalarMul(op=model.op.op2, coeff=model.coeff),
            )
        return None

    pass

OperatorAssociativity

Bases: RewriteRule

Implements associativity of addition, multiplication and kronecker product

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class OperatorAssociativity(RewriteRule):
    """Implements associativity of addition, multiplication and kronecker product"""

    def map_OperatorAdd(self, model):
        return self._map_addmullkron(model=model)

    def map_OperatorMul(self, model):
        return self._map_addmullkron(model=model)

    def map_OperatorKron(self, model):
        return self._map_addmullkron(model=model)

    def _map_addmullkron(self, model):
        if isinstance(model.op2, model.__class__):
            return model.__class__(
                op1=model.__class__(op1=model.op1, op2=model.op2.op1),
                op2=model.op2.op2,
            )
        return model.__class__(op1=model.op1, op2=model.op2)

GatherCoefficient

Bases: RewriteRule

Gathers the coefficients of an operator into a single scalar multiplication for each term

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class GatherCoefficient(RewriteRule):
    """Gathers the coefficients of an operator into a single scalar multiplication for each term"""

    def map_OperatorScalarMul(self, model):
        if isinstance(model.op, OperatorScalarMul):
            return (model.coeff * model.op.coeff) * model.op.op

    def map_OperatorMul(self, model):
        return self._map_mulkron(model)

    def map_OperatorKron(self, model):
        return self._map_mulkron(model)

    def _map_mulkron(self, model):
        if isinstance(model.op1, OperatorScalarMul) and isinstance(
            model.op2, OperatorScalarMul
        ):
            return (
                model.op1.coeff
                * model.op2.coeff
                * model.__class__(op1=model.op1.op, op2=model.op2.op)
            )
        if isinstance(model.op1, OperatorScalarMul):
            return model.op1.coeff * model.__class__(op1=model.op1.op, op2=model.op2)

        if isinstance(model.op2, OperatorScalarMul):
            return model.op2.coeff * model.__class__(op1=model.op1, op2=model.op2.op)

CoefficientDistributivity

Bases: RewriteRule

Implements distributivity of addition over multiplication on coefficient

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class CoefficientDistributivity(RewriteRule):
    """Implements distributivity of addition over multiplication on coefficient"""

    def map_CoefficientMul(self, model):
        if isinstance(model.coeff1, (CoefficientAdd)):
            return model.coeff1.__class__(
                coeff1=CoefficientMul(coeff1=model.coeff1.coeff1, coeff2=model.coeff2),
                coeff2=CoefficientMul(coeff1=model.coeff1.coeff2, coeff2=model.coeff2),
            )
        if isinstance(model.coeff2, (CoefficientAdd)):
            return model.coeff2.__class__(
                coeff1=CoefficientMul(coeff1=model.coeff1, coeff2=model.coeff2.coeff1),
                coeff2=CoefficientMul(coeff1=model.coeff1, coeff2=model.coeff2.coeff2),
            )

CoefficientAssociativity

Bases: RewriteRule

Implements associativity of addition, multiplication on coefficient

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class CoefficientAssociativity(RewriteRule):
    """Implements associativity of addition, multiplication on coefficient"""

    def map_CoefficientAdd(self, model):
        return self._map_addmul(model=model)

    def map_CoefficientMul(self, model):
        return self._map_addmul(model=model)

    def _map_addmul(self, model):
        if isinstance(model.coeff2, model.__class__):
            return model.__class__(
                coeff1=model.__class__(coeff1=model.coeff1, coeff2=model.coeff2.coeff1),
                coeff2=model.coeff2.coeff2,
            )
        return model.__class__(coeff1=model.coeff1, coeff2=model.coeff2)

ScaleTerms

Bases: RewriteRule

Adds a scalar multiplication for each term if it does not exist

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class ScaleTerms(RewriteRule):
    """Adds a scalar multiplication for each term if it does not exist"""

    def __init__(self):
        super().__init__()
        self.op_add_root = False

    def map_AtomicEmulatorGate(self, model):
        self.op_add_root = False

    def map_OperatorAdd(self, model):
        self.op_add_root = True
        if isinstance(model.op1, Union[OperatorAdd, OperatorScalarMul]):
            op1 = model.op1
        else:
            op1 = ConstantCoefficient(value=1) * model.op1
        if isinstance(model.op2, OperatorScalarMul):
            op2 = model.op2
        else:
            op2 = ConstantCoefficient(value=1) * model.op2
        return op1 + op2

    def map_OperatorScalarMul(self, model):
        self.op_add_root = True
        pass

    def map_Operator(self, model):
        if not self.op_add_root:
            self.op_add_root = True
            return ConstantCoefficient(value=1) * model

_CombineTermsHelper

Bases: RewriteRule

Helper for combining terms of the same operator by combining their coefficients

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class _CombineTermsHelper(RewriteRule):
    """Helper for combining terms of the same operator by combining their coefficients"""

    def __init__(self):
        super().__init__()

        self.operators = []

    @property
    def coefficients(self):
        return [o[0] for o in self.operators]

    @property
    def terms(self):
        return [o[1] for o in self.operators]

    def emit(self):
        if self.operators == []:
            return PrunedOperator()
        return reduce(
            lambda op1, op2: op1 + op2,
            [o[0] * o[1] for o in self.operators],
        )

    def map_OperatorScalarMul(self, model):
        if model.op in self.terms:
            i = self.terms.index(model.op)
            self.operators[i] = (model.coeff + self.coefficients[i], model.op)
        else:
            self.operators.append((model.coeff, model.op))

CombineTerms

Bases: RewriteRule

Combines terms of the same operator by combining their coefficients

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class CombineTerms(RewriteRule):
    """Combines terms of the same operator by combining their coefficients"""

    def map_AtomicEmulatorCircuit(self, model):
        return model.__class__(frame=model.frame, sequence=model.sequence)

    def map_AtomicEmulatorGate(self, model):
        combiner = _CombineTermsHelper()
        Pre(combiner)(model)

        return model.__class__(hamiltonian=combiner.emit(), duration=model.duration)

PartitionMathExpr

Bases: RewriteRule

This separates real and complex portions of MathExpr objects.

Parameters:

Name Type Description Default
model MathExpr

The rule only acts on MathExpr objects.

required

Returns:

Name Type Description
model MathExpr
Assumptions

DistributeMathExpr, ProperOrderMathExpr

Example
  • MathStr(string = '1 + 1j + 2') => MathStr(string = '1 + 2 + 1j')
  • MathStr(string = '1 * 1j * 2') => MathStr(string = '1j * 1 * 2')
Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
class PartitionMathExpr(RewriteRule):
    """
    This separates real and complex portions of [`MathExpr`][oqd_core.interface.math.MathExpr] objects.

    Args:
        model (MathExpr): The rule only acts on [`MathExpr`][oqd_core.interface.math.MathExpr] objects.

    Returns:
        model (MathExpr):

    Assumptions:
        [`DistributeMathExpr`][oqd_core.compiler.math.rules.DistributeMathExpr],
        [`ProperOrderMathExpr`][oqd_core.compiler.math.rules.ProperOrderMathExpr]

    Example:
        - MathStr(string = '1 + 1j + 2') => MathStr(string = '1 + 2 + 1j')
        - MathStr(string = '1 * 1j * 2') => MathStr(string = '1j * 1 * 2')
    """

    def map_MathAdd(self, model):
        priority = dict(
            MathImag=5, MathNum=4, MathVar=3, MathFunc=2, MathPow=1, MathMul=0
        )

        if isinstance(
            model.expr2, (MathImag, MathNum, MathVar, MathFunc, MathPow, MathMul)
        ):
            if isinstance(model.expr1, MathAdd):
                if (
                    priority[model.expr2.__class__.__name__]
                    > priority[model.expr1.expr2.__class__.__name__]
                ):
                    return MathAdd(
                        expr1=MathAdd(expr1=model.expr1.expr1, expr2=model.expr2),
                        expr2=model.expr1.expr2,
                    )
            else:
                if (
                    priority[model.expr2.__class__.__name__]
                    > priority[model.expr1.__class__.__name__]
                ):
                    return MathAdd(
                        expr1=model.expr2,
                        expr2=model.expr1,
                    )

    def map_MathMul(self, model: MathMul):
        priority = dict(MathImag=4, MathNum=3, MathVar=2, MathFunc=1, MathPow=0)

        if isinstance(model.expr2, (MathImag, MathNum, MathVar, MathFunc, MathPow)):
            if isinstance(model.expr1, MathMul):
                if (
                    priority[model.expr2.__class__.__name__]
                    > priority[model.expr1.expr2.__class__.__name__]
                ):
                    return MathMul(
                        expr1=MathMul(expr1=model.expr1.expr1, expr2=model.expr2),
                        expr2=model.expr1.expr2,
                    )
            else:
                if (
                    priority[model.expr2.__class__.__name__]
                    > priority[model.expr1.__class__.__name__]
                ):
                    return MathMul(
                        expr1=model.expr2,
                        expr2=model.expr1,
                    )

canonicalize_math_factory()

Creates a new instance of the canonicalization pass for math expressions

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
def canonicalize_math_factory():
    """Creates a new instance of the canonicalization pass for math expressions"""
    return Chain(
        FixedPoint(
            Post(
                Chain(
                    PruneMathExpr(),
                    PruneZeroPowers(),
                    SimplifyMathExpr(),
                    DistributeMathExpr(),
                    ProperOrderMathExpr(),
                )
            )
        ),
        FixedPoint(Post(PartitionMathExpr())),
    )

canonicalize_coefficient_factory()

Creates a new instance of the canonicalization pass for coefficients

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
def canonicalize_coefficient_factory():
    """Creates a new instance of the canonicalization pass for coefficients"""
    return Chain(
        FixedPoint(
            Post(
                Chain(
                    PruneCoefficient(),
                    CoefficientDistributivity(),
                    CombineCoefficient(),
                    CoefficientAssociativity(),
                )
            )
        ),
        FixedPoint(Post(Chain(SortCoefficient(), CombineCoefficient()))),
    )

canonicalize_operator_factory()

Creates a new instance of the canonicalization pass for operators

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
def canonicalize_operator_factory():
    """Creates a new instance of the canonicalization pass for operators"""
    return FixedPoint(
        Post(
            Chain(
                OperatorDistributivity(),
                GatherCoefficient(),
                OperatorAssociativity(),
            )
        )
    )

canonicalize_emulator_circuit_factory()

Creates a new instance of the canonicalization pass for AtomicEmulatorCircuit

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/canonicalize.py
def canonicalize_emulator_circuit_factory():
    """Creates a new instance of the canonicalization pass for AtomicEmulatorCircuit"""
    return Chain(
        canonicalize_operator_factory(),
        canonicalize_coefficient_factory(),
        canonicalize_math_factory(),
        Post(PruneOperator()),
        Pre(ScaleTerms()),
        Post(CombineTerms()),
        canonicalize_coefficient_factory(),
        canonicalize_math_factory(),
        Post(PruneOperator()),
    )

codegen

ConstructHamiltonian

Bases: ConversionRule

Maps an AtomicCircuit to an AtomicEmulatorCircuit replaces laser descriptions of operations with Hamiltonian description of operations

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/codegen.py
class ConstructHamiltonian(ConversionRule):
    """Maps an AtomicCircuit to an AtomicEmulatorCircuit replaces laser descriptions of operations with Hamiltonian description of operations"""

    def map_AtomicCircuit(self, model, operands):
        gates = (
            [operands["protocol"]]
            if isinstance(operands["protocol"], AtomicEmulatorGate)
            else operands["protocol"]
        )
        for gate in gates:
            gate.hamiltonian = gate.hamiltonian + operands["system"]

        return AtomicEmulatorCircuit(sequence=gates)

    def map_System(self, model, operands):
        self.N = len(model.ions)
        self.M = len(model.modes)

        self.ions = model.ions
        self.modes = model.modes

        ops = []
        for n, ion in enumerate(model.ions):
            ops.append(
                reduce(
                    lambda x, y: x @ y,
                    [
                        (
                            self._map_Ion(ion, n)
                            if i == n
                            else Identity(
                                subsystem=f"E{i}" if i < self.N else f"P{i - self.N}"
                            )
                        )
                        for i in range(self.N + self.M)
                    ],
                )
            )

        for m, mode in enumerate(model.modes):
            ops.append(
                reduce(
                    lambda x, y: x @ y,
                    [
                        (
                            self._map_Phonon(mode, m)
                            if i == self.N + m
                            else Identity(
                                subsystem=f"E{i}" if i < self.N else f"P{i - self.N}"
                            )
                        )
                        for i in range(self.N + self.M)
                    ],
                )
            )

        op = reduce(lambda x, y: x + y, ops)
        return op

    def _map_Ion(self, model, index):
        ops = [
            WaveCoefficient(amplitude=level.energy, frequency=0, phase=0)
            * KetBra(ket=n, bra=n, subsystem=f"E{index}")
            for n, level in enumerate(model.levels)
        ]

        op = reduce(lambda x, y: x + y, ops)
        return op

    def _map_Phonon(self, model, index):
        return WaveCoefficient(amplitude=model.energy, frequency=0, phase=0) * (
            Creation(subsystem=f"P{index}") * Annihilation(subsystem=f"P{index}")
        )

    def map_Beam(self, model, operands):
        I = intensity_from_laser(model)  # noqa: E741

        angular_frequency = (
            abs(model.transition.level2.energy - model.transition.level1.energy)
            + model.detuning
        )
        wavevector = angular_frequency * np.array(model.wavevector) / cst.c

        ops = []
        if self.modes:
            displacement_plus = []
            displacement_minus = []
            for m, mode in enumerate(self.modes):
                eta = np.dot(
                    wavevector,
                    mode.eigenvector[model.target * 3 : model.target * 3 + 3],
                ) * np.sqrt(
                    cst.hbar
                    / (2 * self.ions[model.target].mass * cst.m_u * mode.energy)
                )

                displacement_plus.append(
                    Displacement(
                        alpha=WaveCoefficient(
                            amplitude=eta, frequency=0, phase=np.pi / 2
                        ),
                        subsystem=f"P{m}",
                    )
                )
                displacement_minus.append(
                    Displacement(
                        alpha=WaveCoefficient(
                            amplitude=eta, frequency=0, phase=-np.pi / 2
                        ),
                        subsystem=f"P{m}",
                    )
                )

            displacement_plus = reduce(lambda x, y: x @ y, displacement_plus)
            displacement_minus = reduce(lambda x, y: x @ y, displacement_minus)

            for transition in self.ions[model.target].transitions:
                rabi = rabi_from_intensity(model, transition, I)

                ops.append(
                    (
                        reduce(
                            lambda x, y: x @ y,
                            [
                                (
                                    WaveCoefficient(
                                        amplitude=rabi / 2,
                                        frequency=-angular_frequency,
                                        phase=model.phase,
                                    )
                                    * (
                                        KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                        + KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                    )
                                    if i == model.target
                                    else Identity(subsystem=f"E{i}")
                                )
                                for i in range(self.N)
                            ],
                        )
                        @ displacement_plus
                    )
                )

                ops.append(
                    (
                        reduce(
                            lambda x, y: x @ y,
                            [
                                (
                                    WaveCoefficient(
                                        amplitude=rabi / 2,
                                        frequency=angular_frequency,
                                        phase=-model.phase,
                                    )
                                    * (
                                        KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                        + KetBra(
                                            ket=self.ions[model.target].levels.index(
                                                transition.level2
                                            ),
                                            bra=self.ions[model.target].levels.index(
                                                transition.level1
                                            ),
                                            subsystem=f"E{model.target}",
                                        )
                                    )
                                    if i == model.target
                                    else Identity(subsystem=f"E{i}")
                                )
                                for i in range(self.N)
                            ],
                        )
                        @ displacement_minus
                    )
                )

        else:
            for transition in self.ions[model.target].transitions:
                rabi = rabi_from_intensity(model, transition, I)

                ops.append(
                    reduce(
                        lambda x, y: x @ y,
                        [
                            (
                                WaveCoefficient(
                                    amplitude=rabi / 2,
                                    frequency=angular_frequency,
                                    phase=model.phase,
                                )
                                * (
                                    KetBra(
                                        ket=self.ions[model.target].levels.index(
                                            transition.level1
                                        ),
                                        bra=self.ions[model.target].levels.index(
                                            transition.level2
                                        ),
                                        subsystem=f"E{model.target}",
                                    )
                                    + KetBra(
                                        ket=self.ions[model.target].levels.index(
                                            transition.level2
                                        ),
                                        bra=self.ions[model.target].levels.index(
                                            transition.level1
                                        ),
                                        subsystem=f"E{model.target}",
                                    )
                                )
                                if i == model.target
                                else Identity(subsystem=f"E{i}")
                            )
                            for i in range(self.N)
                        ],
                    )
                )

        op = reduce(lambda x, y: x + y, ops)
        return op

    def map_Pulse(self, model, operands):
        return AtomicEmulatorGate(
            hamiltonian=operands["beam"],
            duration=model.duration,
        )

    def map_ParallelProtocol(self, model, operands):
        duration = operands["sequence"][0].duration

        for op in operands["sequence"]:
            assert op.duration == duration

        op = reduce(
            lambda x, y: x + y, map(lambda x: x.hamiltonian, operands["sequence"])
        )
        return AtomicEmulatorGate(hamiltonian=op, duration=duration)

    def map_SequentialProtocol(self, model, operands):
        return operands["sequence"]

utils

compute_matrix_element(laser, transition)

Computes matrix element corresponding to a laser interacting with a particular transition

Parameters:

Name Type Description Default
laser Beam

laser to compute matrix element of

required
transition Transition

transition to compute matrix element of

required

Returns:

Name Type Description
matrix_element float

Multipole matrix elements corresponding to the interaction between the laser and the transition

Warning

Currently implemented for only E1 and E2 transitions.

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def compute_matrix_element(laser, transition):
    """Computes matrix element corresponding to a laser interacting with a particular transition

    Args:
        laser (Beam): laser to compute matrix element of
        transition (Transition): transition to compute matrix element of

    Returns:
        matrix_element (float): Multipole matrix elements corresponding to the interaction between the laser and the transition

    Warning:
        Currently implemented for only `E1` and `E2` transitions.
    """

    # If this happens there's probably an error with the ion species card
    if transition.level1.nuclear != transition.level2.nuclear:
        raise ValueError(
            "Different nuclear spins between two levels in transition:", transition
        )

    # Just by convention; these orderings are set upon instantiation of an Ion object
    if transition.level1.energy > transition.level2.energy:
        raise ValueError("Expected energy of level2 > energy of level1")

    polarization = np.array(laser.polarization).T  # make it a column vector
    wavevector = laser.wavevector

    J1, J2, F1, F2, M1, M2, E1, E2, I, A = (  # noqa: E741
        transition.level1.spin_orbital,
        transition.level2.spin_orbital,
        transition.level1.spin_orbital_nuclear,
        transition.level2.spin_orbital_nuclear,
        transition.level1.spin_orbital_nuclear_magnetization,
        transition.level2.spin_orbital_nuclear_magnetization,
        transition.level1.energy,
        transition.level2.energy,
        transition.level1.nuclear,
        transition.einsteinA,
    )

    q = M2 - M1

    omega = E2 - E1

    # --- M1 transition multipole ---
    if transition.multipole == "M1":
        # The Einstein A coefficient for an M1 transition is given by:
        #   A(M1) = (16 * pi^3)/(3 * hbar) * (omega^3/c^3) * (|<J2||μ||J1>|^2/(2J2+1))

        # Solving for the reduced matrix element (expressed in units of the Bohr magneton, μ_B) :
        #   |<J2||μ||J1>|/μ_B = sqrt((3 * hbar * c^3 * A * (2J2+1))/(16 * pi^3 * omega^3)) / μ_B

        # Reference: I. I. Sobelman, "Atomic Spectra and Radiative Transitions", 2nd ed., Springer (1992).

        # A unit term is definied so that when multiplied by standard angular momentum factors,  the full matrix is reproduced

        units_term = np.sqrt(
            (3 * cst.hbar * cst.epsilon_0 * cst.c**5 * A) / (16 * np.pi**3 * omega**3)
        )

        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 1, F2, F1, I
        )

        # For M1 transitions the operator is a vector operator coupling to the magnetic field
        # Plane wave: magnetic field is proportional to cross product of wavevector and electric field polarization
        B_field = np.cross(wavevector, polarization)

        # Define a spherical basis for a vector (identical to the one used for E1):
        polarization_map = {
            -1: 1 / np.sqrt(2) * np.array([1, 1j, 0]),
            0: np.array([0, 0, 1]),
            1: 1 / np.sqrt(2) * np.array([1, -1j, 0]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * polarization_map[q].dot(B_field)
            * wigner_3j(F2, 1, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    # --- E1 transition multipole ---
    if transition.multipole == "E1":
        units_term = np.sqrt(
            (3 * np.pi * cst.epsilon_0 * cst.hbar * cst.c**3 * A)
            / (omega**3 * cst.e**2)
        )
        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 1, F2, F1, I
        )

        # q -> polarization
        polarization_map = {
            -1: 1 / np.sqrt(2) * np.array([1, 1j, 0]),
            0: np.array([0, 0, 1]),
            1: 1 / np.sqrt(2) * np.array([1, -1j, 0]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * polarization_map[q].dot(polarization)
            * wigner_3j(F2, 1, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    # --- E2 transition multipole ---
    elif transition.multipole == "E2":
        units_term = np.sqrt(
            (15 * np.pi * cst.epsilon_0 * cst.hbar * cst.c**3 * A)
            / (omega**3 * cst.e**2)
        )  # <- anomalous constants I needed to add... hmm
        hyperfine_term = np.sqrt((2 * F2 + 1) * (2 * F1 + 1)) * wigner_6j(
            J1, J2, 2, F2, F1, I
        )

        # q -> polarization
        polarization_map = {
            -2: 1 / np.sqrt(6) * np.array([[1, 1j, 0], [1j, -1, 0], [0, 0, 0]]),
            -1: 1 / np.sqrt(6) * np.array([[0, 0, 1], [0, 0, 1j], [1, 1j, 0]]),
            0: 1 / 3 * np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 2]]),
            1: 1 / np.sqrt(6) * np.array([[0, 0, -1], [0, 0, 1j], [-1, 1j, 0]]),
            2: 1 / np.sqrt(6) * np.array([[1, -1j, 0], [-1j, -1, 0], [0, 0, 0]]),
        }

        geometry_term = (
            np.sqrt(2 * J2 + 1)
            * wavevector.dot(np.matmul(polarization_map[q], polarization))
            * wigner_3j(F2, 2, F1, M2, -q, -M1)
        )

        return float(
            (abs(units_term) * abs(geometry_term) * abs(hyperfine_term)).evalf()
        )

    else:
        raise ValueError(
            "Currently only support dipole and quadrupole allowed transitions"
        )

rabi_from_intensity(laser, transition, intensity)

Computes a transition's resonant rabi frequency when addressed by a laser with intensity

Parameters:

Name Type Description Default
laser Beam

laser to compute resonant rabi frequency of

required
transition Transition

transition to compute resonant rabi frequency of

required
intensity float

intensity of laser

required

Returns:

Name Type Description
rabi_frequency float

resonant rabi frequency corresponding to the interaction between the laser and transition

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def rabi_from_intensity(laser, transition, intensity):
    """Computes a transition's resonant rabi frequency when addressed by a laser with intensity

    Args:
        laser (Beam): laser to compute resonant rabi frequency of
        transition (Transition): transition to compute resonant rabi frequency of
        intensity (float): intensity of laser

    Returns:
        rabi_frequency (float): resonant rabi frequency corresponding to the interaction between the laser and transition
    """

    matrix_elem = compute_matrix_element(laser, transition)

    if transition.multipole[0] == "E":
        E = (2 * intensity / (cst.epsilon_0 * cst.c)) ** 0.5
        return matrix_elem * E * cst.e / cst.hbar
    if transition.multipole[0] == "M":
        B = (2 * intensity / (cst.epsilon_0 * cst.c**3)) ** 0.5
        return matrix_elem * B / cst.hbar

intensity_from_laser(laser)

Computes the intensity from a laser

Parameters:

Name Type Description Default
laser Beam

laser to compute intensity of.

required

Returns:

Name Type Description
intensity float

intensity of the laser

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/utils.py
def intensity_from_laser(laser):
    """Computes the intensity from a laser

    Args:
        laser (Beam): laser to compute intensity of.

    Returns:
        intensity (float): intensity of the laser
    """
    matrix_elem = compute_matrix_element(laser, laser.transition)

    if laser.transition.multipole[0] == "E":
        return (
            cst.c
            * cst.epsilon_0
            / 2
            * (cst.hbar * laser.rabi / (matrix_elem * cst.e)) ** 2
        )  # noqa: E741

    if laser.transition.multipole[0] == "M":
        return cst.c**3 * cst.epsilon_0 / 2 * (cst.hbar * laser.rabi / matrix_elem) ** 2  # noqa: E741

visualization

CoefficientPrinter

Bases: ConversionRule

Prints Coefficients in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class CoefficientPrinter(ConversionRule):
    """Prints Coefficients in a pretty manner"""

    def map_MathExpr(self, model, operands):
        return Post(PrintMathExpr())(model)

    def map_MathAdd(self, model, operands):
        return f"({Post(PrintMathExpr())(model)})"

    def map_WaveCoefficient(self, model, operands):
        frequency_term = (
            f"{operands['frequency']} * t"
            if model.frequency != MathNum(value=0)
            else ""
        )
        phase_term = f"{operands['phase']}" if model.phase != MathNum(value=0) else ""

        wave_term = (
            f" * exp(1j * ({frequency_term}{' + ' if frequency_term and phase_term else ''}{phase_term}))"
            if phase_term or frequency_term
            else ""
        )

        return f"{operands['amplitude']}{wave_term}"

    def map_CoefficientAdd(self, model, operands):
        return f"{'(' + operands['coeff1'] + ')' if isinstance(model.coeff1, CoefficientAdd) else operands['coeff1']} + {'(' + operands['coeff2'] + ')' if isinstance(model.coeff2, CoefficientAdd) else operands['coeff2']}"

    def map_CoefficientMul(self, model, operands):
        return f"{'(' + operands['coeff1'] + ')' if isinstance(model.coeff1, CoefficientAdd) else operands['coeff1']} * {'(' + operands['coeff2'] + ')' if isinstance(model.coeff2, CoefficientAdd) else operands['coeff2']}"

OperatorPrinter

Bases: ConversionRule

Prints Operators in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class OperatorPrinter(ConversionRule):
    """Prints Operators in a pretty manner"""

    def map_KetBra(self, model, operands):
        return f"|{model.ket}><{model.bra}|_{model.subsystem}"

    def map_Annihilation(self, model, operands):
        return f"A_{model.subsystem}"

    def map_Creation(self, model, operands):
        return f"C_{model.subsystem}"

    def map_Identity(self, model, operands):
        return f"I_{model.subsystem}"

    def map_PrunedOperator(self, model, operands):
        return "PrunedOperator"

    def map_Displacement(self, model, operands):
        return f"D({operands['alpha']})_{model.subsystem}"

    def map_OperatorAdd(self, model, operands):
        return f"{operands['op1']} + {operands['op2']}"

    def map_OperatorMul(self, model, operands):
        return f"{'(' + operands['op1'] + ')' if isinstance(model.op1, OperatorAdd) else operands['op1']} * {'(' + operands['op2'] + ')' if isinstance(model.op2, OperatorAdd) else operands['op2']}"

    def map_OperatorKron(self, model, operands):
        return f"{'(' + operands['op1'] + ')' if isinstance(model.op1, OperatorAdd) else operands['op1']} @ {'(' + operands['op2'] + ')' if isinstance(model.op2, OperatorAdd) else operands['op2']}"

    def map_OperatorScalarMul(self, model, operands):
        return f"{'(' + operands['coeff'] + ')' if isinstance(model.coeff, CoefficientAdd) else operands['coeff']} * {'(' + operands['op'] + ')' if isinstance(model.op, OperatorAdd) else operands['op']}"

    def map_Coefficient(self, model, operands):
        return Post(CoefficientPrinter())(model)

CondensedOperatorPrettyPrint

Bases: PrettyPrint

Prints An AtomicEmulatorCircuit in a pretty manner

Source code in oqd-trical/src/oqd_trical/light_matter/compiler/visualization.py
class CondensedOperatorPrettyPrint(PrettyPrint):
    """Prints An AtomicEmulatorCircuit in a pretty manner"""

    def map_Operator(self, model, operands):
        return f"Operator({Post(OperatorPrinter())(model)})"