def two_qubit_matrix_to_ion_operations(q0: ops.Qid,
                                       q1: ops.Qid,
                                       mat: np.ndarray,
                                       atol: float = 1e-8
                                       ) -> List[ops.Operation]:
    """Decomposes a two-qubit operation into MS/single-qubit rotation gates.

    Args:
        q0: The first qubit being operated on.
        q1: The other qubit being operated on.
        mat: Defines the operation to apply to the pair of qubits.
        tolerance: A limit on the amount of error introduced by the
            construction.

    Returns:
        A list of operations implementing the matrix.
    """
    kak = linalg.kak_decomposition(mat, atol=atol)
    operations = _kak_decomposition_to_operations(q0,
        q1, kak, atol)
    return _cleanup_operations(operations)