How to get a noisy density matrix after partial measurements in Qiskit

Question

I'm trying to study the effect of noise on a quantum system with partial measurements using Qiskit. Suppose I have $N + M$ qubits and I run a quantum circuit and measure $M$ of them. I want to know the quantum state on the remaining $N$ qubits given the measurement outcomes on the $M$ qubits I measured. In other words, I want to post-select (condition) on some measurement outcomes.

If I'm working with a noiseless circuit, I can do this by building my usual circuit and appending the instruction circuit.snapshot(label = "psi", snapshot_type = "statevector") after my measurements. Then, I can use the AerSimulator to simulate my circuit and get out the quantum state as a snapshot after any series of measurements using:

results = simulator.run(circuit, shots = shots, memory = True).result()
states = results.data()["snapshots"]["statevector"]["psi"]

But what happens if I want to implement noise in my simulator? I know how to do this using AerSimulator, but I'm not sure how to get the state I'm looking for.

One option is to use the instruction circuit.snapshot(label = "rho", snapshot_type = "density_matrix"). Is this the recommended way to do the simulation when I want to see the effect of noise after partially measuring a quantum system?

Answer 1

You can save the simulator state using simulator instructions such as save_density_matrix. The table in the linked page shows what are the supported simulation methods for each instruction. So, you can pick a simulator that supports noisy simulation then choose one of the available instructions.

As an example, let's create a noise model and pass it to statevector simulator:

from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel
from qiskit_ibm_provider import IBMProvider

provider = IBMProvider()
backend = provider.get_backend('ibmq_manila')
noise_model = NoiseModel.from_backend(backend)

noisy_simulator = AerSimulator(method = 'statevector', noise_model = noise_model)

Now, we create a sample circuit of five qubits and measure two of them. We insert save_density_matrix instruction to save the state of the remaining qubits, and set conditional parameter to True to associate the state with the measurement result:

from qiskit.circuit import QuantumCircuit

circ = QuantumCircuit(5, 2)
circ.h(2)
circ.cx(2, 1)
circ.cx(1, 0)
circ.cx(2, 3)
circ.cx(3, 4)
circ.barrier()
circ.measure(0, 0)
circ.measure(1, 1)
circ.save_density_matrix(qubits=[2, 3, 4], label="rho", conditional=True)
circ.draw('mpl')

Now we run the experiment and display the state:

from IPython.display import display

job = noisy_simulator.run(circ)
result = job.result()

# Get statevector for each measurement:
for meas, state in result.data()['rho'].items():
    display(meas)
    display(state.draw('latex'))

The output should look like: