Circuit including phase factor in $XY(\beta, \theta)$ gate

Question

In Implementation of the XY interaction family with calibration of a single pulse, the $XY(\beta, \theta)$ gate is defined as

$$ XY(\beta, \theta) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos(\theta/2) & i\sin(\theta/2)e^{i\beta} & 0 \\ 0 & i\sin(\theta/2)e^{-i\beta} & \cos(\theta/2) & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} $$

Later on, they simplify it so that $\beta$ defaults to $0$ and therefore the following decomposition can be used to build the gate:

from qiskit import QuantumCircuit
import numpy as np

qc = QuantumCircuit(2)
qc.rz(-np.pi/2, 0)
qc.sx(0) 
qc.rz(np.pi/2, 0)
qc.s(1)
qc.cx(0, 1)
qc.ry(theta/2, [0, 1])
qc.cx(0, 1)
qc.sdg(1)
qc.rz(-np.pi/2, 0)
qc.sxdg(0) 
qc.rz(np.pi/2, 0)

However, I'm interested in also being able to control the phase factor in the two off diagonal center elements. I tried doing controlled-phase gates surrounded with $X$ gates to add phase to $|01\rangle$ and $|10\rangle$ but adds a phase to everything in the central four elements, plus it's costly.

Does anyone have a reference or know a way to implement the two-parameter $XY(\beta, \theta)$ gate that, hopefully, only requires the two $CX$ gates from the $\beta=0$ case or at least doesn't add a lot more $CX$ gates to the original circuit?

Answer 1

There's probably some simple direct way to have operations that depend on Beta, but as a fallback you can always use the KAK decomposition to turn the 4x4 matrix into the optimal number of two qubit gates.

For example, cirq can do the kak decomposition and also has functionality that uses it to spit out CZ-optimal circuits.

Example:

import cirq
import numpy as np


def compile(theta, beta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    matrix = np.array([
        [1, 0, 0, 0],
        [0, c, 1j * s * np.exp(1j * beta), 0],
        [0, 1j * s * np.exp(-1j * beta), c, 0],
        [0, 0, 0, 1],
    ])
    print(cirq.kak_decomposition(matrix))

    circuit = cirq.Circuit(cirq.two_qubit_matrix_to_operations(
        cirq.LineQubit(0),
        cirq.LineQubit(1),
        matrix,
        allow_partial_czs=False,
    ))
    cirq.MergeSingleQubitGates(synthesizer=lambda q, arr: cirq.MatrixGate(arr).on(q)).optimize_circuit(circuit)
    print(circuit.to_text_diagram(use_unicode_characters=False))

Example output:

# >>> compile(theta=np.pi / 3, beta=np.pi / 16 * 7)

KAK {
    xyz*(4/π): 0.333, 0.333, 3.53e-17
    before: (-0.0312*π around Z) ⊗ (-0.469*π around Z)
    after: (0.0312*π around Z) ⊗ (0.469*π around Z)
}

0: ---[[0.5  +0.5j   0.449-0.547j]----@---[[ 0.933+0.25j  -0.042+0.255j]----@---[[ 0.5  +0.5j   -0.449+0.547j]----
       [0.547-0.449j 0.5  +0.5j  ]]   |    [-0.091+0.242j  0.933+0.25j ]]   |    [-0.547+0.449j  0.5  +0.5j  ]]
                                      |                                     |
1: ---[[0.5  +0.5j   0.547-0.449j]----@---[[0.933+0.25j  0.091-0.242j]------@---[[ 0.5  +0.5j   -0.547+0.449j]----
       [0.449-0.547j 0.5  +0.5j  ]]        [0.042-0.255j 0.933+0.25j ]]          [-0.449+0.547j  0.5  +0.5j  ]]

Answer 2

After working on this for a while, I came up with the following decomposition:

from qiskit.circuit.quantumcircuit import QuantumCircuit

qc = QuantumCircuit(2)
    
qc.rz(beta, 1)
qc.rz(-np.pi / 2, 0)
qc.sx(0) 
qc.rz(np.pi/2, 0)
qc.s(1)
qc.cx(0, 1)
qc.ry(theta / 2, [0, 1])
qc.cx(0, 1)
qc.sdg(1)
qc.rz(-np.pi/2, 0)
qc.sxdg(0) 
qc.rz(np.pi/2, 0)
qc.rz(-beta, 1)

Which can be tested to be correct with the following snippet:

from qiskit.quantum_info import Operator
import numpy as np

xy = lambda theta, beta: np.array(
    [
        [1, 0, 0, 0],
        [0, np.cos(theta / 2), 1j * np.sin(theta / 2) * np.exp(1j * beta), 0],
        [0, 1j * np.sin(theta / 2) * np.exp(-1j * beta), np.cos(theta / 2), 0],
        [0, 0, 0, 1],
    ]
)

for i in range(100):
    theta = np.random.uniform(0, 2 * np.pi)
    beta = np.random.uniform(0, 2 * np.pi)

    qc = QuantumCircuit(2)
    qc.rz(beta, 1)
    qc.rz(-np.pi / 2, 0)
    qc.sx(0) 
    qc.rz(np.pi/2, 0)
    qc.s(1)
    qc.cx(0, 1)
    qc.ry(theta / 2, [0, 1])
    qc.cx(0, 1)
    qc.sdg(1)
    qc.rz(-np.pi/2, 0)
    qc.sxdg(0) 
    qc.rz(np.pi/2, 0)
    qc.rz(-beta, 1)

    assert np.allclose(Operator(qc).data, xy(theta, 0)), f"Not correct for {theta}, {beta}"

As you can see, this decomposition only adds two more $R_z$ gates, so the $\rm CX$ count remains the same as I wanted on my original question.