Skip to content

Code Generation

oqd_trical.light_matter.compiler

codegen

ConstructHamiltonian

Bases: ConversionRule

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

Source code in src\oqd_trical\light_matter\compiler\codegen.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
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 src\oqd_trical\light_matter\compiler\utils.py 
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
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 src\oqd_trical\light_matter\compiler\utils.py 
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
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 src\oqd_trical\light_matter\compiler\utils.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
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