Add a TrappedIon ModelProcessor

[ajrazander]

Trapped ions need more detailed gate analysis (at the laser pulse level) in open source packages. Adding a TrappedIon ModelProcessor class would make qutip a great place to do detailed noise analysis on quantum gates in trapped ions--maybe the only place that isn't behind a pay wall.

There is already a nice list of ModelProcessors, but trapped ions don't quite fit any on the list. A SpinChain is not detailed enough. For example, it ignores the phonon modes of the ions which mediate the multi-qubit Molmer Sorensen gate used in most experiments. These modes are important to keep track of in error analysis. The CavityQED is actually bit closer with the inclusion of a phonon mode, but there is only one phonon mode while for N ions there are N phonon modes.

I propose that a TrappedIon ModelProcessor should be added.

I'd be happy to add it. This would be my first contribution to qutip-qip, so I'd appreciate some guidance.

[BoxiLi]

Hi @ajrazander, thanks for your interest!

It is also what we are quite interested in. A model for ion would fit nicely into the collection of different quantum devices.

I guess the starting point would be formulating the implementation of native gates for trapped ion in the qutip (or, more generally, quantum control) language:
H(t) = H0 + c_1(t) * H1 + c_2(t) * H2 + …
where c_j(t) are time-dependent complex numbers representing the drive pulse and H_j are constant Qobj. I saw in your profile that you are an expert in trapped ions so you likely know this better than me.
Then one can add a TrappedIonModel/TrappedIonProcessor and a TrappedIonCompiler or something like that to /device and /compiler.
Lastly, there will need to be some tests that verify if everything works properly.

I’m happy to discuss this in more detail or explain the code structure a bit more through a call if needed.

[ajrazander]

Great! This sounds promising!

I'll share an update and probably ask for some feedback once I take a first pass at writing a TrappedIonModel/TrappedIonProcessor

[ajrazander]

I'm going to need to pad operators with identities to get spin and phonon operators to work together. Is there something like ikron from quimb in qutip? https://quimb.readthedocs.io/en/latest/autoapi/quimb/core/index.html#quimb.core.ikron

[BoxiLi]

Is this what you are looking for?
https://qutip.org/docs/latest/guide/guide-tensor.html#tensor-products

[ajrazander]

I'm looking for a nice way to handle tensoring all the identity operators. For example, if I have 3 qubits, and I want to define a $\sigma_x$ on each qubit for the control Hamiltonian, I'll need the operators $I \otimes I \otimes \sigma_x$, $I \otimes \sigma_x \otimes I$, and $\sigma_x \otimes I \otimes I$. The nice part about ikron is you can pick where you want your operator ( $\sigma_x$ in this example), specify what kind of dimension it's supposed to operate in ([2,2,2]), where you want the operator, and it tensors on all the identities for you. Is there a best practice for this in qutip? Or is a for loop carefully working with tensor the standard way to go?

[hodgestar]

Are you perhaps looking for qutip.expand_operator?

[ajrazander]

Thank you @hodgestar! expand_operator is exactly what I was looking for. It's also in qutip_qip.operations.gates docs

@BoxiLi I am working on the compiler for the Molmer Sorensen gate, and I'm a bit confused about what label to give the two-qubit drive. I'm learning from this notebook. In input cell 19, the pulse_info lists two tuples. My understanding is those are the separate drives on each qubit? Let's assume for simplicity that's how the MS gate is physically driven, I'm confused about how to label each pulse though. I'd like to label them both as "ms_01" as that's the label I gave the drive in the control Hamiltonian, but a duplicate error is thrown if I can give the same label to each qubit in the pulse_info list. In the example notebook the two drives are labeled as "sx_[index]" and "sx_[index+1]" but that references the single qubit drive from the Hamiltonian at the beginning of the notebook. It should reference "g0" from the original Hamiltonian, yes? Would you explain more about these pulse labels and what's going on behind the scenes here?

[BoxiLi]

Each pair in pulse_info corresponds to one time-dependent coefficient c_j(t) of the Hamiltonian c_j(t) * H_j.

In that particular cell, I was just showing how you can drive multiple Hamiltonian terms to implement a gate. The example does not have any physical meaning because it is just implementing simultaneous single-qubit X gates on both qubits. Which can be done separately.

With this explained, I'm wondering why you want to label them both as ms_01? Here the label should match the label of the Hamiltonian. It doesn't really make much sense to have two different pulses c_j(t), addressing the same Hamiltonian, right?

[ajrazander]

Ah, thank you for the explanation! Now, I'm labeling the one pulse as 'ms_01', and I'd like it to show up on the .plot_pulses()--but nothing shows up. There's a space in the timeline for the MS gate, but there's no pulse shape.

Wait, I just realized I need to use the use_control_latex=False option. Nevermind! The MS pulse is showing :)

[BoxiLi]

Great! Although I think even with use_control_latex=True, it should still show the MS gate. There might indeed be a mistake in the code somewhere. I'll double check.

[ajrazander]

I'm seeing an error in the calculation of the pulse area. Specifically line 416 of gatecompiler.py. The time axis (tlist) of the pulse is being scaled to get the area of the pulse to match the input area, but the current calculation tlist *= abs(area) / np.abs(maximum) is not general enough. I believe this only works for square pulses. For a generic pulse shape, the normalization factor has to be calculated from the integral of the pulse shape. Something like this current_pulse_area = (tlist[1] - tlist[0]) * sum(coeff * tlist) and then normalize tlist tlist /= np.sqrt(current_pulse_area).

In summary replace:
tlist *= abs(area) / np.abs(maximum)

with

dt = tlist[1] - tlist[0]
pulse_area = dt * sum(tlist * coeff)  # compute pulse area
tlist *= np.sqrt(abs(area) / pulse_area)  # normalize tlist so the pulses's area = area

Should I write this up as a separate issue?

[BoxiLi]

I guess the documentation might be a bit unclear here

qutip-qip/src/qutip_qip/compiler/gatecompiler.py

Lines 413 to 417 in 5790ddc

	coeff, tlist = _normalized_window(shape, num_samples)
	sign = np.sign(area)
	coeff *= np.abs(maximum) * sign
	tlist *= abs(area) / np.abs(maximum)
	return coeff, tlist

The _normalized_window returns coeff and tlist that integrates to 1. These are precalculated to avoid recomputing this at each compilation run.
The rest lines stretch the pulse so that the maximum matches the input maximum. If you ignore the sign, you see that coeff is multiplied by it and tlist is divided by it. In addition to that, tlist is also multiplied by area. So the resulting integral should be area and the maximal coeff equals to maximum.

[ajrazander]

Thank you for the thorough explanation! I see now that I was miscalculating the pulse area. My mistake! I did notice some odd behavior with the boxcar window. The tlist and coeff output looks right and it plots as a rectangle window (as expected), but

processor.plot_pulses(use_control_latex=False, show_axis=True)
plt.show()

Shows a right triangle

The state dynamics look correct though.

[ajrazander]

@BoxiLi I have a more general question for gates. I'd like to define a custom gate such as a general rotation operator $R(\theta, \phi) = e^{i \theta/2 ( \cos(\phi) X + \sin(\phi) Y)}$. The circuit definition could look like this

def R_phi(angles):
    theta, phi = angles
    # (cos(phi) X + sin(phi) Y) gate
    R_phi_op = Qobj([[np.cos(theta / 2), np.exp(-1j * phi) * np.sin(theta/2)],
          [np.exp(1j * phi) * np.sin(theta/2), np.cos(theta / 2)]])
    return R_phi_op

I tried integrate this gate with the TrappedIonModel, but I'm not sure the best way to go about it. The users sets the angle phi when they make add R_phi((theta, phi)) to a quantum circuit, but the control Hamiltonian would need to be adjusted on a per gate basis according the phi. Is that possible? Are phase terms adjustable in the control Hamiltonian at all, or is it just global coefficients? Using this kind of rotation gate is also mentioned in the second after of this comment, but there was no follow up.

[BoxiLi]

The tlist and coeff output looks right and it plots as a rectangle window (as expected), but ...

That figure could indeed be a mistake. Usually if one wants to use a discrete pulse like boxcar. One can directly set processor.pulse_mode="discrete". boxcard was not much tested. If you still have it, could you post the coeff and list when you got this plot?

For the second question, yes, it can be done. At the Hamiltonian level, the key problem is how to implement
exp(1.j*theta/2) * (cos(phi) X + sin(phi) Y)
with a different theta and phi for each gate. I would suggest the following:

• Define two independent control Hamiltonian term X and Y.
• Extract the gate parameter theta and phi in the compiler. And compute the coefficient exp(1.j*theta/2) * cos(phi) for X and exp(1.j*theta/2) * sin(phi) for Y.

Here is a simplified version of the code, where the phase is hard coded. In your case you need to extract it from the Gate parameters:

qutip-qip/doc/pulse-paper/customize.py

Lines 70 to 107 in 5790ddc

	def generate_pulse(self, gate, tlist, coeff, phase=0.0):
	"""Generates the pulses.
	Args:
	gate (qutip_qip.circuit.Gate): A qutip Gate object.
	tlist (array): A list of times for the evolution.
	coeff (array): An array of coefficients for the gate pulses
	phase (float): The value of the phase for the gate.
	Returns:
	Instruction (qutip_qip.compiler.instruction.Instruction): An instruction
	to implement a gate containing the control pulses.
	"""
	pulse_info = [
	# (control label, coeff)
	("sx" + str(gate.targets[0]), np.cos(phase) * coeff),
	("sy" + str(gate.targets[0]), np.sin(phase) * coeff),
	]
	return [Instruction(gate, tlist=tlist, pulse_info=pulse_info)]
	def single_qubit_gate_compiler(self, gate, args):
	"""Compiles single qubit gates to pulses.
	Args:
	gate (qutip_qip.circuit.Gate): A qutip Gate object.
	Returns:
	Instruction (qutip_qip.compiler.instruction.Instruction): An instruction
	to implement a gate containing the control pulses.
	"""
	# gate.arg_value is the rotation angle
	tlist = np.abs(gate.arg_value) / self.params["pulse_amplitude"]
	coeff = self.params["pulse_amplitude"] * np.sign(gate.arg_value)
	coeff /= 2 * np.pi
	if gate.name == "RX":
	return self.generate_pulse(gate, tlist, coeff, phase=0.0)
	elif gate.name == "RY":
	return self.generate_pulse(gate, tlist, coeff, phase=np.pi / 2)

Also, the circuit QED model also implemented this. The DRAG pulse adds a Y correction pulse to an X drive:

qutip-qip/src/qutip_qip/compiler/circuitqedcompiler.py

Lines 103 to 159 in 5790ddc

	def _rotation_compiler(self, gate, op_label, param_label, args):
	"""
	Single qubit rotation compiler.
	Parameters
	----------
	gate : :obj:`~.operations.Gate`:
	The quantum gate to be compiled.
	op_label : str
	Label of the corresponding control Hamiltonian.
	param_label : str
	Label of the hardware parameters saved in
	:obj:`GateCompiler.params`.
	args : dict
	The compilation configuration defined in the attributes
	:obj:`.GateCompiler.args` or given as a parameter in
	:obj:`.GateCompiler.compile`.
	Returns
	-------
	A list of :obj:`.Instruction`, including the compiled pulse
	information for this gate.
	"""
	targets = gate.targets
	coeff, tlist = self.generate_pulse_shape(
	args["shape"],
	args["num_samples"],
	maximum=self.params[param_label][targets[0]],
	area=gate.arg_value / 2.0 / np.pi,
	)
	if args["DRAG"]:
	pulse_info = self._drag_pulse(op_label, coeff, tlist, targets[0])
	else:
	pulse_info = [(op_label + str(targets[0]), coeff)]
	return [Instruction(gate, tlist, pulse_info)]
	def _drag_pulse(self, op_label, coeff, tlist, target):
	dt_coeff = np.gradient(coeff, tlist[1] - tlist[0]) / 2 / np.pi
	# Y-DRAG
	alpha = self.params["alpha"][target]
	y_drag = -dt_coeff / alpha
	# Z-DRAG
	z_drag = -(coeff**2) / alpha + (np.sqrt(2) ** 2 * coeff**2) / (
	4 * alpha
	)
	# X-DRAG
	coeff += -(coeff**3 / (4 * alpha**2))
	pulse_info = [
	(op_label + str(target), coeff),
	("sz" + str(target), z_drag),
	]
	if op_label == "sx":
	pulse_info.append(("sy" + str(target), y_drag))
	elif op_label == "sy":
	pulse_info.append(("sx" + str(target), -y_drag))
	return pulse_info

In your case, you need additionally

theta, phi = gate.args

And then multiply the coeff with the corresponding phase.

[ajrazander]

Ok, excellent! I've attached a test file that defines the classes and tests them out on a circuit. Based on my checks, everything is passing at this point. If this looks like it's heading in the right direction, I can put the changes in a pull request.
trappedion_test.zip

[BoxiLi]

Looks great! The fidelity is a bit low though. But I guess one needs to fine-tune the Hamiltonian to really generate a good one. You could move some inline comments to the documentation of the class. Also, it would be nice to move those tests to a separate file e.g. tests/test_trappedion, which will then be executed by pytest.

[BoxiLi]

I realized that there is actually a MS gate and a QROT gate defined in the master branch. But for some reason, the QROT gate is not correctly exposed to the circuit interface. I'll fix that.

[ajrazander]

The pulse fidelity is my fault! I picked the wrong final state haha I'll change that!

Your comment about the MS and QROT gates already in qutip got me thinking... The MS gate can rotate around XX or YY "axes" depending on a phase. Could we add an Rxx and Ryy gate to qutip? They're effectively two-body MS gates with different phases. That way we won't need to define any custom gates.

[BoxiLi]

Hey sorry about the delay. I just merged the PR for adding the R gate.

Could we add an Rxx and Ryy gate to qutip?

This sounds like a great idea! Especially if it is going to be useful for the Trapped Ion. Would you like to open a PR for this first, since you might have the formula and everything with you? This will be a simpler PR compared to adding the Trapped ions.

[ajrazander]

Busy week, but I just added the PR for the Rxx and Ryy gates here

[BoxiLi]

Hey @ajrazander

Is there any future plan on adding the PR for the full ion processor?

[ajrazander]

@BoxiLi Thanks for following up. Yes, the plan is to incorporate the new MS gate, check things are working as expected, then move the example ion processor into qutip. It might take a week or two.