`StateTomography` reconstructs wrong density matrix when `c_if` involved

[ARO-GZ]

Informations

• Qiskit Experiments version: 0.4.0
• Python version: 3.10.4
• Operating system: Ubuntu 20.04 (Windows10 with WSL1)

What is the current behavior?

teleport_bell_state (see code below) creates a bell state between qubits 0 and 3 and then teleports qubit 0 to qubit 2 through qubit 1. The state of qubits 2 and 3 at the end should then be the bell state initially shared between 0 and 1. However if we perform tomography on those qubits (2 and 3) then state one gets

[[ 0.36+0.j   -0.01+0.j    0.02+0.02j  0.37-0.j  ]
 [-0.01-0.j    0.14+0.j   -0.12-0.j   -0.01-0.01j]
 [ 0.02-0.02j -0.12+0.j    0.12+0.j    0.02-0.01j]
 [ 0.37+0.j   -0.01+0.01j  0.02+0.01j  0.37-0.j  ]]
Fidelity =  0.7319716283458715

I also checked that it doesn't work if you try to teleport a single qubit in a computational state and reconstruct the state with tomography.

However it does work (both two cases) if one uses the principle of deferred measurements to replace the classically conditional operations used in teleportation (teleport_bell_state_deferred).

[[0.49+0.j   0.  -0.j   0.01-0.01j 0.5 -0.j  ]
 [0.  +0.j   0.  +0.j   0.  +0.j   0.  +0.j  ]
 [0.01+0.01j 0.  -0.j   0.  +0.j   0.01+0.01j]
 [0.5 +0.j   0.  -0.j   0.01-0.01j 0.51+0.j  ]]
Fidelity =  0.9935749708197947

Steps to reproduce the problem

from qiskit_experiments.library import StateTomography
from qiskit.quantum_info import state_fidelity
from qiskit.providers.aer import AerSimulator
import qiskit

import numpy as np

def get_state(qc,qubits=[2,3]):
    backend=Aer.get_backend("qasm_simulator")

    qstexp = StateTomography(qc,measurement_qubits=qubits)
    qstdata = qstexp.run(backend, 
        seed_simulation=100).block_for_results()
    state_result = qstdata.analysis_results("state").value.data
    print(state_result.round(2))
    return state_result

def get_fidelity(state_result):
    phi = np.array([1.,0.,0.,1.])/np.sqrt(2)
    fid = state_fidelity(phi,state_result)
    print('Fidelity = ',fid)
    return fid

def teleport_bell_state():
    """Teleport qubit 0 -> 2"""
    qr = QuantumRegister(4)
    c0 = ClassicalRegister(1)
    c1 = ClassicalRegister(1)
    teleport = QuantumCircuit(qr, c0, c1)
    teleport.h(qr[0])
    teleport.cx(qr[0],qr[3])
    teleport.h(qr[1])
    teleport.cx(qr[1], qr[2])
    teleport.cx(qr[0], qr[1])
    teleport.h(qr[0])
    teleport.measure(qr[0], c0[0])
    teleport.measure(qr[1], c1[0])
    teleport.z(qr[2]).c_if(c0, 1)
    teleport.x(qr[2]).c_if(c1, 1)
    return teleport

def teleport_bell_state_deferred():
    """Teleport qubit 0 -> 2"""
    qr = QuantumRegister(4)
    c0 = ClassicalRegister(1)
    c1 = ClassicalRegister(1)
    ancillas = QuantumRegister(2,'ancilla')
    teleport = QuantumCircuit(qr, ancillas, c0, c1)
    teleport.h(qr[0])
    teleport.cx(qr[0],qr[3])
    teleport.h(qr[1])
    teleport.cx(qr[1], qr[2])
    teleport.cx(qr[0], qr[1])
    teleport.h(qr[0])
    teleport.cx(qr[0], ancillas[0])
    teleport.cx(qr[1], ancillas[1])
    teleport.cz(ancillas[0],qr[2])
    teleport.cx(ancillas[1],qr[2])
    teleport.measure(ancillas[0],c0)
    teleport.measure(ancillas[1],c1)
    return teleport

if __name__ == '__main__':
    # qc = teleport_bell_state_deferred()
    qc = teleport_bell_state()
    print(qc)
    state = get_state(qc,qubits=[2,3])
    fid = get_fidelity(state)

What is the expected behavior?

One should be able to reconstruct the correct density matrix.

Suggested solutions

It seems that there is a problem with the way the tomography experiment handles the classical registers in the case in which classically conditioned operations are used.

[chriseclectic]

Yep, it looks like something is going wrong with the classical registers. If you print the generated circuits they seem to be involved in the bell conditionals when they shouldn't be.

Original circuit:

      ┌───┐               ┌───┐┌─┐          
q0_0: ┤ H ├──■─────────■──┤ H ├┤M├──────────
      ├───┤  │       ┌─┴─┐└┬─┬┘└╥┘          
q0_1: ┤ H ├──┼────■──┤ X ├─┤M├──╫───────────
      └───┘  │  ┌─┴─┐└───┘ └╥┘  ║ ┌───┐┌───┐
q0_2: ───────┼──┤ X ├───────╫───╫─┤ Z ├┤ X ├
           ┌─┴─┐└───┘       ║   ║ └─╥─┘└─╥─┘
q0_3: ─────┤ X ├────────────╫───╫───╫────╫──
           └───┘            ║   ║   ║    ║  
  c0: ══════════════════════╬═══╩═══■════╬══
                            ║      0x1   ║  
  c1: ══════════════════════╩════════════■══
                                        0x1 

State tomography circuit

     ┌───┐               ┌───┐┌─┐                                    
q_0: ┤ H ├──■─────────■──┤ H ├┤M├────────────────────────────────────
     ├───┤  │       ┌─┴─┐└┬─┬┘└╥┘                                    
q_1: ┤ H ├──┼────■──┤ X ├─┤M├──╫─────────────────────────────────────
     └───┘  │  ┌─┴─┐└───┘ └╥┘  ║ ┌───┐┌───┐ ░ ┌────────────┐ ░ ┌─┐   
q_2: ───────┼──┤ X ├───────╫───╫─┤ Z ├┤ X ├─░─┤ PauliMeasZ ├─░─┤M├───
          ┌─┴─┐└───┘       ║   ║ └─╥─┘└─╥─┘ ░ ├────────────┤ ░ └╥┘┌─┐
q_3: ─────┤ X ├────────────╫───╫───╫────╫───░─┤ PauliMeasZ ├─░──╫─┤M├
          └───┘            ║   ║   ║    ║   ░ └────────────┘ ░  ║ └╥┘
c_0: ══════════════════════╬═══╩═══■════o═══════════════════════╬══╬═
                           ║       ║    ║                       ║  ║ 
c_1: ══════════════════════╩═══════o════■═══════════════════════╬══╬═
                                   ║    ║                       ║  ║ 
c_2: ══════════════════════════════o════o═══════════════════════╩══╬═
                                   ║    ║                          ║ 
c_3: ══════════════════════════════o════o══════════════════════════╩═

Since these measurement circuits are constructed by using QuantumCircuit.compose I expect this might be due to a bug in how that function handles flattening classical registers during composition.