GHZ State Preparation

Create Greenberger-Horne-Zeilinger (GHZ) multi-qubit entangled states and verify their unique correlations using pyqpanda3.

---

Problem

What is a GHZ State?

A GHZ state (named after Greenberger, Horne, and Zeilinger) is a specific type of multipartite quantum entanglement that generalizes the two-qubit Bell state to three or more qubits. For a system of $n$ qubits, the GHZ state is defined as:

$$|{\text{GHZ}}_n⟩=\frac{1}{\sqrt{2}}(|0{⟩}^{\otimes n}+|1{⟩}^{\otimes n})$$$$|{\text{GHZ}}_3⟩=\frac{1}{\sqrt{2}}(|000⟩+|111⟩)$$

Unlike pairwise entanglement (Bell pairs), a GHZ state exhibits genuine multipartite entanglement -- no bipartite decomposition can fully describe the correlations.

All qubits are simultaneously in the all-zeros state and the all-ones state.

Why GHZ States Matter

GHZ states are central to quantum information science:

1. Bell inequality violations without statistics: The GHZ argument demonstrates a contradiction between quantum mechanics and local hidden variable theories using deterministic predictions, no statistical inequalities needed.
2. Quantum networking: GHZ states are a resource for quantum secret sharing, conference key agreement, and distributed quantum consensus.
3. Quantum metrology: Multipartite entanglement enables Heisenberg-limited phase estimation ($1/n$ scaling rather than $1/\sqrt{n}$).
4. Fault-tolerant quantum computing: GHZ states are building blocks for error-correcting codes and magic state distillation.
5. Foundational tests: GHZ states expose the non-local character of quantum mechanics in its starkest form.

The Mermin Inequality

For a 3-qubit GHZ state, the Mermin operator provides a Bell-type test:

$$M=X\otimes X\otimes X-X\otimes Y\otimes Y-Y\otimes X\otimes Y-Y\otimes Y\otimes X$$

Local hidden variable theories satisfy $|⟨M⟩|\le 2$, Quantum mechanics predicts $⟨M⟩=4$ -- a maximal violation.

This is the strongest known argument against local realism that does not rely on statistical bounds.

---

Solution

Step-by-step Approach

Step 1 -- Initialize all qubits to $|0⟩$

$$|{\psi }_0⟩=|0{⟩}^{\otimes n}$$

Step 2 -- Apply Hadamard to the first qubit

$$H|0⟩=\frac{1}{\sqrt{2}}(|0⟩+|1⟩)$$$$|{\psi }_1⟩=\frac{1}{\sqrt{2}}(|0{⟩}^{\otimes n}+|1⟩|0{⟩}^{\otimes (n-1)})$$

Step 3 -- Apply cascading CNOT gates

Apply $\text{CNOT}(0,k)$ for each qubit $k=1,2,\dots ,n-1$, where qubit 0 is always the control:

Each CNOT copies the control value to the target:

$$|{\psi }_3⟩=\frac{1}{\sqrt{2}}(|0{⟩}^{\otimes n}+|1{⟩}^{\otimes n})=|{\text{GHZ}}_n⟩$$

Step 4 -- Measure all qubits

All qubits yield the same outcome: either all 0 or all 1, Each with probability $1/2$.

Circuit Diagram (3-qubit GHZ)

Gate count: The circuit requires 1 Hadamard gate and $n-1$ CNOT gates, Gate depth is $n$.

---

Code

3-Qubit GHZ State

from pyqpanda3 import core

# Build a 3-qubit GHZ circuit
prog = core.QProg()
prog << core.H(0)            # Hadamard on qubit 0
prog << core.CNOT(0, 1)      # Entangle qubit 0 with qubit 1
prog << core.CNOT(0, 2)      # Entangle qubit 0 with qubit 2
prog << core.measure([0, 1, 2], [0, 1, 2])

# Run on the CPU simulator
machine = core.CPUQVM()
machine.run(prog, shots=10000)
counts = machine.result().get_counts()

print("3-qubit GHZ measurement counts:")
for outcome, count in sorted(counts.items()):
    print(f"  |{outcome}> : {count}")
# Expected: |000> ~5000, |111> ~5000, no other outcomes

General n-Qubit GHZ State Function

from pyqpanda3 import core
from typing import Dict


def build_ghz_circuit(n: int) -> core.QProg:
    """Build an n-qubit GHZ state preparation circuit."""
    if n < 2:
        raise ValueError(f"GHZ state requires n >= 2 qubits, got {n}")

    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    prog << core.measure(list(range(n)), list(range(n)))
    return prog


def run_ghz(n: int, shots: int = 10000) -> Dict[str, int]:
    """Prepare and measure an n-qubit GHZ state."""
    prog = build_ghz_circuit(n)
    machine = core.CPUQVM()
    machine.run(prog, shots)
    return machine.result().get_counts()


# --- Example usage ---
for n in [3, 5, 8]:
    counts = run_ghz(n, shots=10000)
    ideal_keys = {"0" * n, "1" * n}
    ghz_counts = sum(v for k, v in counts.items() if k in ideal_keys)
    fidelity_est = ghz_counts / sum(counts.values())
    print(f"n={n}: GHZ fidelity estimate = {fidelity_est:.4f}")

Verification via Measurement Statistics

A genuine GHZ state has two hallmark properties: only $|0{⟩}^{\otimes n}$ and $|1{⟩}^{\otimes n}$ appear, Each close to 50%.

from pyqpanda3 import core


def verify_ghz_statistics(n: int, shots: int = 20000) -> None:
    """Verify GHZ state: only |000...0> and |111...1> with ~50% each."""
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    prog << core.measure(list(range(n)), list(range(n)))

    machine = core.CPUQVM()
    machine.run(prog, shots)
    counts = machine.result().get_counts()

    all_zeros = "0" * n
    all_ones = "1" * n
    total = sum(counts.values())
    count_0 = counts.get(all_zeros, 0)
    count_1 = counts.get(all_ones, 0)
    count_other = total - count_0 - count_1

    print(f"--- GHZ({n}) Verification ({shots} shots) ---")
    print(f"  |{all_zeros}> : {count_0} ({100*count_0/total:.1f}%)")
    print(f"  |{all_ones}> : {count_1} ({100*count_1/total:.1f}%)")
    print(f"  Other outcomes : {count_other} ({100*count_other/total:.2f}%)")

    assert count_other / total < 0.01, "Too many non-GHZ outcomes"
    print("  PASS: State is consistent with GHZ.")


for n in [3, 4, 5, 6]:
    verify_ghz_statistics(n)

GHZ with DensityMatrixSimulator

The DensityMatrixSimulator gives exact probabilities without shot noise, Making it ideal for analyzing purity, reduced density matrices, and fidelity:

from pyqpanda3 import core
import numpy as np

# Build the 3-qubit GHZ circuit without measurement
prog = core.QProg()
prog << core.H(0)
prog << core.CNOT(0, 1)
prog << core.CNOT(0, 2)

# Run on the density matrix simulator
dm_sim = core.DensityMatrixSimulator()
dm_sim.run(prog)

# Get exact probabilities
probs = dm_sim.state_probs()
print("Exact probabilities (3-qubit GHZ):")
for idx, p in enumerate(probs):
    label = format(idx, '03b')
    if p > 1e-10:
        print(f"  |{label}>: {p:.6f}")

# Get the full density matrix (8x8 for 3 qubits)
dm = dm_sim.density_matrix()
print(f"Density matrix shape: {dm.shape}")

# Verify purity. For a pure state, Tr(rho^2) = 1.0.
purity = np.trace(dm @ dm).real
print(f"Purity Tr(rho^2): {purity:.6f}")
assert abs(purity - 1.0) < 1e-10, "State is not pure!"

# Compute the reduced density matrix for qubit 0.
# For the GHZ state, each subsystem is maximally mixed: I/2.
reduced_q0 = dm_sim.reduced_density_matrix([0])
print(f"Reduced density matrix (qubit 0):\n{reduced_q0}")
reduced_purity = np.trace(reduced_q0 @ reduced_q0).real
print(f"Reduced purity: {reduced_purity:.6f} (maximally mixed = 0.5)")

GHZ with Noise Simulation

Real quantum devices are noisy. This example shows how depolarizing noise degrades the GHZ state, Causing previously forbidden outcomes to appear:

from pyqpanda3 import core


def noisy_ghz(n: int, shots: int = 10000) -> None:
    """Simulate a GHZ state under depolarizing noise."""
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    prog << core.measure(list(range(n)), list(range(n)))

    # Define noise model
    noise = core.NoiseModel()
    for q in range(n):
        noise.add_quantum_error(
            core.depolarizing_error(0.01), core.GateType.H, [q]
        )
    for k in range(1, n):
        noise.add_quantum_error(
            core.depolarizing_error(0.02), core.GateType.CNOT, [0, k]
        )

    # Noisy simulation
    machine = core.CPUQVM()
    machine.run(prog, shots, model=noise)
    counts = machine.result().get_counts()

    # Ideal for comparison
    machine.run(prog, shots)
    ideal_counts = machine.result().get_counts()

    all_zeros = "0" * n
    all_ones = "1" * n
    total = sum(counts.values())
    ghz_count = counts.get(all_zeros, 0) + counts.get(all_ones, 0)

    print(f"--- GHZ({n}) Noise Comparison ({shots} shots) ---")
    print(f"  Ideal: {dict(sorted(ideal_counts.items()))}")
    print(f"  Noisy: {dict(sorted(counts.items()))}")
    print(f"  GHZ outcomes: {ghz_count}/{total} ({100*ghz_count/total:.1f}%)")


noisy_ghz(3, shots=10000)
noisy_ghz(5, shots=10000)

Comparing Ideal vs. Noisy GHZ Correlations

The DensityMatrixSimulator gives exact probabilities under noise (no shot noise), Making it ideal for precise comparison:

from pyqpanda3 import core
import numpy as np


def compare_ideal_noisy_ghz(n: int) -> None:
    """Compare ideal and noisy GHZ probability distributions exactly."""
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)

    # Ideal
    dm_ideal = core.DensityMatrixSimulator()
    dm_ideal.run(prog)
    probs_ideal = np.array(dm_ideal.state_probs())
    dm_ideal_mat = dm_ideal.density_matrix()
    purity_ideal = np.trace(dm_ideal_mat @ dm_ideal_mat).real

    # Noisy
    noise = core.NoiseModel()
    for q in range(n):
        noise.add_quantum_error(
            core.depolarizing_error(0.02), core.GateType.H, [q]
        )
    for k in range(1, n):
        noise.add_quantum_error(
            core.depolarizing_error(0.02), core.GateType.CNOT, [0, k]
        )

    dm_noisy = core.DensityMatrixSimulator()
    dm_noisy.run(prog, model=noise)
    probs_noisy = np.array(dm_noisy.state_probs())
    dm_noisy_mat = dm_noisy.density_matrix()
    purity_noisy = np.trace(dm_noisy_mat @ dm_noisy_mat).real

    # Print comparison
    dim = 2 ** n
    print(f"\n--- GHZ({n}) Exact Probability Comparison ---")
    print(f"{'Basis':<{n+2}} {'Ideal':>10} {'Noisy':>10} {'Delta':>10}")
    print("-" * (n + 35))
    for i in range(dim):
        bs = format(i, f'0{n}b')
        d = abs(probs_ideal[i] - probs_noisy[i])
        if probs_ideal[i] > 1e-6 or probs_noisy[i] > 1e-6:
            print(f"|{bs}> {probs_ideal[i]:10.6f} {probs_noisy[i]:10.6f} {d:10.6f}")

    print(f"\nIdeal purity:  {purity_ideal:.6f}")
    print(f"Noisy purity:  {purity_noisy:.6f}")
    print(f"Purity loss:   {purity_ideal - purity_noisy:.6f}")


compare_ideal_noisy_ghz(3)
compare_ideal_noisy_ghz(4)

GHZ State Fidelity Sweep

Fidelity $F=⟨\psi |\rho |\psi ⟩$ quantifies how close a noisy state is to the ideal GHZ state. This example sweeps the noise strength to observe degradation:

from pyqpanda3 import core
import numpy as np


def ghz_fidelity_dm(n: int, noise: core.NoiseModel) -> float:
    """Compute exact fidelity F = <GHZ_n|rho|GHZ_n> using density matrices."""
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)

    dm_sim = core.DensityMatrixSimulator()
    dm_sim.run(prog, model=noise)
    rho = dm_sim.density_matrix()

    dim = 2 ** n
    ghz_vec = np.zeros(dim, dtype=complex)
    ghz_vec[0] = 1.0 / np.sqrt(2)
    ghz_vec[dim - 1] = 1.0 / np.sqrt(2)

    fidelity = float(np.real(ghz_vec.conj() @ rho @ ghz_vec))
    return fidelity


n = 4
print(f"GHZ({n}) fidelity vs. depolarizing error strength:")
print(f"{'Error rate':>12s}  {'Fidelity':>10s}")
print("-" * 25)
for p in [0.0, 0.005, 0.01, 0.02, 0.05, 0.1]:
    noise = core.NoiseModel()
    if p > 0:
        for q in range(n):
            noise.add_quantum_error(
                core.depolarizing_error(p), core.GateType.H, [q]
            )
        for k in range(1, n):
            noise.add_quantum_error(
                core.depolarizing_error(p), core.GateType.CNOT, [0, k]
            )
    fid = ghz_fidelity_dm(n, noise)
    print(f"{p:>12.3f}  {fid:>10.6f}")

Mermin Inequality Test on the GHZ State

This example estimates the Mermin operator $⟨M⟩$ to demonstrate the non-local correlations of the GHZ state:

from pyqpanda3 import core
from pyqpanda3.hamiltonian import PauliOperator


def mermin_inequality_test() -> None:
    """Test the Mermin inequality on a 3-qubit GHZ state.

    M = +XXX - XYY - YXY - YYX
    Classical bound: |<M>| <= 2, Quantum prediction: <M> = 4
    """
    terms = [
        (+1.0, "X0 X1 X2"),
        (-1.0, "X0 Y1 Y2"),
        (-1.0, "Y0 X1 Y2"),
        (-1.0, "Y0 Y1 X2"),
    ]

    mermin_value = 0.0
    for coeff, pauli_str in terms:
        prog = core.QProg()
        prog << core.H(0)
        prog << core.CNOT(0, 1)
        prog << core.CNOT(0, 2)

        op = PauliOperator(pauli_str)
        exp_val = core.expval_pauli_operator(prog, op)
        contribution = coeff * exp_val
        mermin_value += contribution
        print(f"  <{pauli_str}> = {exp_val:.6f}, contribution = {contribution:.6f}")

    print(f"\nMermin operator <M> = {mermin_value:.6f}")
    print(f"Classical bound: |<M>| <= 2")
    print(f"Quantum prediction: <M> = 4")
    if abs(mermin_value) > 2:
        print("RESULT: Mermin inequality VIOLATED!")
    else:
        print("RESULT: No violation detected.")


mermin_inequality_test()

---

Explanation

Why Cascading CNOTs Create GHZ States

Starting from $|000⟩$, The Hadamard on qubit 0 creates a superposition:

$$|0⟩|00⟩\stackrel{H_0}{\to }\frac{1}{\sqrt{2}}(|000⟩+|100⟩)$$

Each CNOT copies the control (qubit 0) to a target:

$$\stackrel{{\text{CNOT}}_{01}}{\to }\frac{1}{\sqrt{2}}(|000⟩+|110⟩)\stackrel{{\text{CNOT}}_{02}}{\to }\frac{1}{\sqrt{2}}(|000⟩+|111⟩)$$

For $n$ qubits, $n-1$ CNOTs suffice because qubit 0 acts as a broadcast control that synchronizes all targets.

Entanglement Structure: GHZ vs. Bell States

The 2-qubit GHZ state $|{\text{GHZ}}_2⟩=\frac{1}{\sqrt{2}}(|00⟩+|11⟩)$ is identical to the Bell state $|{\mathrm{\Phi }}^{+}⟩$. However, for $n\ge 3$, GHZ states exhibit fundamentally different entanglement:

Property	Bell State ($n=2$)	GHZ State ($n\ge 3$)
Entanglement	Bipartite	Genuine multipartite
Tracing out one qubit	Remaining qubit is maximally mixed	Remaining qubits are separable mixed states
Robustness	Moderate	Very fragile
Bell violation	CHSH inequality ($S\le 2\sqrt{2}$)	Mermin inequality ($⟨M⟩=2^{n-1}$)

A striking property: if you trace out any single qubit, the remaining $n-1$ qubits become a completely separable mixed state -- all entanglement is destroyed. This is unlike W-states, where bipartite entanglement persists after tracing.

The Mermin Inequality and GHZ Correlations

The derivation uses $X$ swaps $|0⟩\leftrightarrow |1⟩$, while $Y|0⟩=i|1⟩$ and $Y|1⟩=-i|0⟩$:

• $X_0{X}_1{X}_2|\text{GHZ}⟩=+|\text{GHZ}⟩$ (all bits swap, superposition invariant)
• $Y_0{Y}_1{X}_2|\text{GHZ}⟩=(i^2)|111⟩+((-i{)}^2)|000⟩=-|\text{GHZ}⟩$

Yielding $⟨M⟩=(+1)-(-1)-(-1)-(-1)=4$. For the general $n$-qubit case, the MABK inequality gives:

$$|⟨M_n⟩|\le 2^{(n-1)/2}\text{(classical)},⟨M_n⟩=2^{n-1}\text{(quantum)}$$

The ratio between quantum and classical bounds grows exponentially: $2^{(n-1)/2}$.

Sensitivity to Decoherence: Why GHZ States Are Fragile

The GHZ state's coherence lives in a single off-diagonal density matrix element ${\rho }_{0,2^n-1}=1/2$. Under independent phase-damping on each qubit with rate $\gamma$, this element decays as:

$${\rho }_{0,2^n-1}(t)=\frac{1}{2}{e}^{-n\gamma t}$$

The effective decoherence rate scales linearly with $n$ -- larger GHZ states decohere $n$ times faster than a single qubit.

Qubits ($n$)	CNOT gates	Depth	Fidelity at 1% CNOT error
3	2	3	~0.96
5	4	5	~0.92
10	9	10	~0.83
20	19	20	~0.68
50	49	50	~0.40

Mitigation strategies include error detection codes, fast gates (reducing gate time relative to $T_1$, $T_2$), verification protocols using ancilla-based parity checks, and parallel CNOT architectures that reduce depth from $n$ to 2.

---

GHZ State Verification

Problem

After preparing a GHZ state, how do you verify that the circuit actually produced the correct entangled state? Simple measurement counts are not enough to confirm genuine multipartite entanglement. You need structured verification that checks both the population (correct basis states) and the coherence (correct phase relationships).

Solution: Multi-Layer Verification

A complete GHZ verification protocol checks three independent criteria:

1. Population test: Only $|0{⟩}^{\otimes n}$ and $|1{⟩}^{\otimes n}$ should appear, each near 50%.
2. Parity test: Adjacent qubit pairs must have even parity (both same). Any XOR of neighboring qubits must be 0.
3. Coherence test (optional): Rotate all qubits into the X basis and measure; all outcomes should again be correlated.

Code: Full GHZ Verification with CPUQVM

from pyqpanda3 import core
from typing import Dict, List


def ghz_population_test(counts: Dict[str, int], n: int, tol: float = 0.02) -> bool:
    """Check that only |000...0> and |111...1> appear with ~50% each."""
    total = sum(counts.values())
    all_zeros = "0" * n
    all_ones = "1" * n

    p0 = counts.get(all_zeros, 0) / total
    p1 = counts.get(all_ones, 0) / total
    p_other = 1.0 - p0 - p1

    print(f"  Population test:")
    print(f"    |{all_zeros}>: {p0:.4f}  (expected 0.5000)")
    print(f"    |{all_ones}>: {p1:.4f}  (expected 0.5000)")
    print(f"    Other:        {p_other:.4f}  (expected 0.0000)")

    pop_ok = (abs(p0 - 0.5) < tol) and (abs(p1 - 0.5) < tol) and (p_other < tol)
    print(f"    Result: {'PASS' if pop_ok else 'FAIL'}")
    return pop_ok


def ghz_parity_test(counts: Dict[str, int], n: int) -> bool:
    """Check that adjacent qubits always have even parity (XOR = 0).

    For a genuine GHZ state, every outcome has all qubits equal.
    So for any pair of adjacent qubits (i, i+1), XOR should be 0.
    """
    total = sum(counts.values())
    parity_violations = 0

    print(f"  Parity test (adjacent qubit XOR):")
    for i in range(n - 1):
        xor0_count = 0  # count of outcomes where qubit_i XOR qubit_{i+1} = 0
        for outcome, count in counts.items():
            bit_i = int(outcome[i])
            bit_j = int(outcome[i + 1])
            if bit_i == bit_j:  # XOR = 0
                xor0_count += count
        parity_frac = xor0_count / total
        print(f"    q[{i}] XOR q[{i+1}] = 0: {parity_frac:.4f}  (expected 1.0000)")
        if parity_frac < 0.98:
            parity_violations += 1

    parity_ok = parity_violations == 0
    print(f"    Result: {'PASS' if parity_ok else 'FAIL'} ({parity_violations} violations)")
    return parity_ok


def ghz_coherence_test(n: int, shots: int = 20000) -> bool:
    """Verify coherence by rotating all qubits into the X basis (H on each).

    A GHZ state rotated into the X basis produces correlated outcomes:
    all outcomes have even parity (even number of 1s).
    """
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    # Rotate into X basis: apply H to every qubit before measuring
    for k in range(n):
        prog << core.H(k)
    prog << core.measure(list(range(n)), list(range(n)))

    machine = core.CPUQVM()
    machine.run(prog, shots)
    counts = machine.result().get_counts()

    total = sum(counts.values())
    even_parity_count = 0
    for outcome, count in counts.items():
        num_ones = outcome.count("1")
        if num_ones % 2 == 0:
            even_parity_count += count

    coherence_frac = even_parity_count / total
    print(f"  Coherence test (X-basis even parity):")
    print(f"    Even parity fraction: {coherence_frac:.4f}  (expected 1.0000)")

    coherence_ok = coherence_frac > 0.98
    print(f"    Result: {'PASS' if coherence_ok else 'FAIL'}")
    return coherence_ok


def verify_ghz_full(n: int, shots: int = 20000) -> None:
    """Run all three verification layers on an n-qubit GHZ state."""
    print(f"\n{'='*50}")
    print(f"Full GHZ({n}) Verification ({shots} shots)")
    print(f"{'='*50}")

    # Prepare and measure in computational basis
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    prog << core.measure(list(range(n)), list(range(n)))

    machine = core.CPUQVM()
    machine.run(prog, shots)
    counts = machine.result().get_counts()

    results = []
    results.append(ghz_population_test(counts, n))
    results.append(ghz_parity_test(counts, n))
    results.append(ghz_coherence_test(n, shots))

    print(f"\n  Overall: {'ALL TESTS PASSED' if all(results) else 'SOME TESTS FAILED'}")


# Run verification for several qubit counts
for n in [3, 4, 5]:
    verify_ghz_full(n)

Explanation

The verification protocol works as follows:

• Population test: A perfect GHZ state only has support on $|0{⟩}^{\otimes n}$ and $|1{⟩}^{\otimes n}$. Any population in other basis states indicates errors (gate errors, decoherence, or readout errors).

• Parity test: In every correct GHZ outcome, all qubits are identical. For any adjacent pair $(i,i+1)$, the XOR of the measured bits must be 0. This test is sensitive to correlated errors that might leave the population intact but break entanglement.

• Coherence test: Applying $H^{\otimes n}$ before measurement rotates the GHZ state into a superposition of all even-parity computational basis states. For 3 qubits, the rotated state becomes $\frac{1}{2}(|000⟩+|011⟩+|101⟩+|110⟩)$. Counting even-parity outcomes verifies the relative phase between the $|0{⟩}^{\otimes n}$ and $|1{⟩}^{\otimes n}$ components is correct.

Together, these three tests catch different error modes:

• Bit-flip errors cause population and parity failures.
• Phase-flip errors leave population and parity intact but fail the coherence test.
• Depolarizing errors cause failures across all three tests.

---

Mermin Inequality Test

Problem

The Mermin inequality provides a Bell-type test for genuine multipartite non-locality. For a 3-qubit system, local hidden variable theories constrain the Mermin operator to $|⟨M⟩|\le 2$, while quantum mechanics predicts $|⟨M⟩|=4$ for the GHZ state -- a maximal violation by a factor of 2.

How do you measure this experimentally using pyqpanda3?

The Mermin Operator

For 3 qubits, the Mermin operator is:

$$M=X\otimes X\otimes X-X\otimes Y\otimes Y-Y\otimes X\otimes Y-Y\otimes Y\otimes X$$

Each term requires measuring a different combination of Pauli observables. The key insight is that:

• $X$ measurements require applying $H$ before computational basis measurement.
• $Y$ measurements require applying $S^{†}H$ (or equivalently $H{S}^{†}$) before measurement.

The expected values on the GHZ state are:

• $⟨XXX⟩=+1$
• $⟨XYY⟩=-1$
• $⟨YXY⟩=-1$
• $⟨YYX⟩=-1$

Giving $⟨M⟩=1-(-1)-(-1)-(-1)=4$.

Code: Measuring the Mermin Operator in X and Y Bases

from pyqpanda3 import core
import numpy as np


def measure_pauli_term(
    n: int,
    basis: str,
    shots: int = 50000,
    noise: core.NoiseModel = None,
) -> float:
    """Measure the expectation value of a Pauli string on the GHZ state.

    Args:
        n: number of qubits (e.g., 3).
        basis: Pauli basis string, e.g., "XXY" (one char per qubit).
               'X' = measure in X basis, 'Y' = measure in Y basis,
               'Z' = measure in Z (computational) basis.
        shots: number of measurement shots.
        noise: optional noise model for realistic simulation.

    Returns:
        Expectation value <basis> in range [-1, +1].
    """
    assert len(basis) == n, f"Basis length {len(basis)} != n qubits {n}"

    prog = core.QProg()
    # Prepare GHZ state
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)

    # Basis rotation: apply gates so that measuring in Z gives X or Y observable
    for i, pauli in enumerate(basis):
        if pauli == "X":
            prog << core.H(i)
        elif pauli == "Y":
            prog << core.S(i)       # S = phase gate (S dagger equivalent via sign)
            prog << core.H(i)
        # 'Z' needs no rotation

    prog << core.measure(list(range(n)), list(range(n)))

    machine = core.CPUQVM()
    if noise is not None:
        machine.run(prog, shots, model=noise)
    else:
        machine.run(prog, shots)
    counts = machine.result().get_counts()

    # Compute expectation: parity-weighted average
    total = sum(counts.values())
    expectation = 0.0
    for outcome, count in counts.items():
        # Count parity of measured bits
        parity = sum(int(b) for b in outcome) % 2
        sign = (-1) ** parity
        expectation += sign * count
    expectation /= total
    return expectation


def mermin_test_3qubit(shots: int = 50000, noise: core.NoiseModel = None) -> None:
    """Run the Mermin inequality test on a 3-qubit GHZ state.

    M = +XXX - XYY - YXY - YYX
    Classical bound: |<M>| <= 2
    Quantum prediction: <M> = 4
    """
    terms = [
        (+1.0, "XXX"),
        (-1.0, "XYY"),
        (-1.0, "YXY"),
        (-1.0, "YYX"),
    ]

    print(f"Mermin Inequality Test (3-qubit GHZ, {shots} shots)")
    print(f"{'Term':>8s}  {'Expectation':>12s}  {'Contribution':>13s}")
    print("-" * 40)

    mermin_value = 0.0
    for coeff, pauli_str in terms:
        exp_val = measure_pauli_term(3, pauli_str, shots=shots, noise=noise)
        contribution = coeff * exp_val
        mermin_value += contribution
        print(f"{pauli_str:>8s}  {exp_val:>12.6f}  {contribution:>13.6f}")

    print(f"\nMermin operator <M> = {mermin_value:.6f}")
    print(f"Classical bound:     |<M>| <= 2")
    print(f"Quantum prediction:  <M> = 4")
    print(f"Violation ratio:     {mermin_value / 2:.4f}")

    if abs(mermin_value) > 2:
        print(f"RESULT: Mermin inequality VIOLATED (|<M>| = {abs(mermin_value):.4f} > 2)")
    else:
        print(f"RESULT: No violation (|<M>| = {abs(mermin_value):.4f} <= 2)")


# Ideal (noiseless) test
print("=" * 50)
print("IDEAL (noiseless) Mermin test")
print("=" * 50)
mermin_test_3qubit(shots=50000)

# Test under depolarizing noise
print("\n" + "=" * 50)
print("NOISY Mermin test (1% depolarizing)")
print("=" * 50)
noise = core.NoiseModel()
noise.add_quantum_error(
    core.depolarizing_error(0.01), core.GateType.H, [0]
)
noise.add_quantum_error(
    core.depolarizing_error(0.01), core.GateType.H, [1]
)
noise.add_quantum_error(
    core.depolarizing_error(0.01), core.GateType.H, [2]
)
noise.add_quantum_error(
    core.depolarizing_error(0.01), core.GateType.CNOT, [0, 1]
)
noise.add_quantum_error(
    core.depolarizing_error(0.01), core.GateType.CNOT, [0, 2]
)
mermin_test_3qubit(shots=50000, noise=noise)

Explanation

The Mermin test works by measuring each Pauli term independently:

1. Basis rotation: To measure $⟨X⟩$ on a qubit, apply $H$ before the computational basis measurement. To measure $⟨Y⟩$, apply $S^{†}H$. This maps the eigenstates of $X$ (or $Y$) to those of $Z$.

2. Parity computation: After rotation, each shot produces a bitstring. The eigenvalue $+1$ corresponds to even parity (even number of 1s), and $-1$ to odd parity. The expectation value is the parity-weighted average.

3. Violation significance: The classical bound of 2 follows from any local hidden variable theory. Quantum mechanics achieves 4 because the GHZ state has correlations that no classical system can reproduce. Even moderate noise (1% depolarizing) reduces the violation, demonstrating why error mitigation is critical for Bell tests on real hardware.

---

GHZ State Under Noise: Fidelity Degradation

Problem

How does depolarizing noise affect GHZ state quality as the number of qubits increases? Because the GHZ circuit uses a cascade of $n-1$ CNOT gates, errors accumulate. Understanding this scaling is essential for determining the maximum feasible GHZ size on near-term hardware.

Fidelity Under Depolarizing Noise

The depolarizing channel on a single qubit replaces the state with the maximally mixed state with probability $p$, and leaves it unchanged with probability $1-p$:

$$\mathcal{E}(\rho )=(1-p)\rho +\frac{p}{3}(X\rho X+Y\rho Y+Z\rho Z)$$

For the GHZ state, fidelity decays approximately as:

$$F\approx {(1-\frac{4p}{3})}^{n-1}$$

where $n-1$ is the number of CNOT gates. This exponential decay means that even modest single-gate error rates make large GHZ states impractical without error correction.

Code: Sweep GHZ Fidelity vs. Qubit Count Under Noise

from pyqpanda3 import core
import numpy as np


def ghz_fidelity(n: int, noise: core.NoiseModel, shots: int = 10000) -> float:
    """Estimate GHZ fidelity via sampling.

    Fidelity is estimated as the fraction of shots in |000...0> or |111...1>.
    This is a lower bound on the true fidelity (does not account for
    coherent errors or phase damping that preserves populations).
    """
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)
    prog << core.measure(list(range(n)), list(range(n)))

    machine = core.CPUQVM()
    machine.run(prog, shots, model=noise)
    counts = machine.result().get_counts()

    all_zeros = "0" * n
    all_ones = "1" * n
    total = sum(counts.values())
    ghz_count = counts.get(all_zeros, 0) + counts.get(all_ones, 0)
    return ghz_count / total


def ghz_exact_fidelity(n: int, noise: core.NoiseModel) -> float:
    """Compute exact GHZ fidelity F = <GHZ_n|rho|GHZ_n> via density matrix."""
    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)

    dm_sim = core.DensityMatrixSimulator()
    dm_sim.run(prog, model=noise)
    rho = dm_sim.density_matrix()

    dim = 2 ** n
    ghz_vec = np.zeros(dim, dtype=complex)
    ghz_vec[0] = 1.0 / np.sqrt(2)
    ghz_vec[dim - 1] = 1.0 / np.sqrt(2)

    fidelity = float(np.real(ghz_vec.conj() @ rho @ ghz_vec))
    return fidelity


# Sweep: GHZ fidelity vs. qubit count at fixed noise strength
noise_rate = 0.02  # 2% depolarizing error per gate
qubit_range = range(2, 9)  # 2 to 8 qubits

print(f"GHZ Fidelity vs. Qubit Count (depolarizing error = {noise_rate:.0%})")
print(f"{'Qubits':>7s}  {'CNOTs':>6s}  {'Exact F':>10s}  {'Sampled F':>10s}  {'Theoretical':>12s}")
print("-" * 55)

for n in qubit_range:
    noise = core.NoiseModel()
    for q in range(n):
        noise.add_quantum_error(
            core.depolarizing_error(noise_rate), core.GateType.H, [q]
        )
    for k in range(1, n):
        noise.add_quantum_error(
            core.depolarizing_error(noise_rate), core.GateType.CNOT, [0, k]
        )

    exact_fid = ghz_exact_fidelity(n, noise)
    sampled_fid = ghz_fidelity(n, noise, shots=20000)
    # Theoretical approximation: (1 - 4p/3)^(n-1) * (1 - 4p/3) for the Hadamard
    theoretical = ((1 - 4 * noise_rate / 3) ** n)

    num_cnots = n - 1
    print(f"{n:>7d}  {num_cnots:>6d}  {exact_fid:>10.6f}  {sampled_fid:>10.4f}  {theoretical:>12.6f}")

# Sweep: Fidelity vs. noise strength for a fixed qubit count
print(f"\n\nGHZ(4) Fidelity vs. Noise Strength")
print(f"{'Error rate':>12s}  {'Exact F':>10s}  {'Purity':>10s}")
print("-" * 37)

n = 4
for p in [0.0, 0.005, 0.01, 0.02, 0.03, 0.05, 0.08, 0.10, 0.15, 0.20]:
    noise = core.NoiseModel()
    if p > 0:
        for q in range(n):
            noise.add_quantum_error(
                core.depolarizing_error(p), core.GateType.H, [q]
            )
        for k in range(1, n):
            noise.add_quantum_error(
                core.depolarizing_error(p), core.GateType.CNOT, [0, k]
            )

    prog = core.QProg()
    prog << core.H(0)
    for k in range(1, n):
        prog << core.CNOT(0, k)

    dm_sim = core.DensityMatrixSimulator()
    dm_sim.run(prog, model=noise)
    rho = dm_sim.density_matrix()
    fidelity = ghz_exact_fidelity(n, noise)
    purity = float(np.trace(rho @ rho).real)

    print(f"{p:>12.3f}  {fidelity:>10.6f}  {purity:>10.6f}")

Explanation

The sweep reveals several important patterns:

Qubits	CNOTs	Fidelity at 2% error
2	1	~0.97
3	2	~0.95
4	3	~0.92
5	4	~0.89
6	5	~0.86
7	6	~0.83
8	7	~0.80

Key observations:

1. Linear gate count, exponential fidelity loss: Each additional qubit adds one CNOT, and the fidelity loss compounds. Going from 2 to 8 qubits at 2% error drops fidelity from 97% to roughly 80%.

2. Sampled vs. exact fidelity: The sampled estimate (from measurement shots) is a lower bound on the true fidelity because it only detects population errors. Phase errors that destroy coherence without changing populations are invisible to simple computational basis measurement. The density matrix simulator gives the true fidelity.

3. Cross-over point: At a given error rate, there is a maximum GHZ size beyond which the state is no longer useful. For 2% depolarizing noise, the practical limit is around 8 qubits (80% fidelity). Error mitigation or error correction is needed beyond this.

4. Purity tracks fidelity: The purity $\text{Tr}({\rho }^2)$ closely tracks fidelity for depolarizing noise, since both decay as the noise mixes the pure GHZ state with the maximally mixed state.

---

Summary

Aspect	Details
State definition	$|{\text{GHZ}}_n⟩=\frac{1}{\sqrt{2}}(|0{⟩}^{\otimes n}+|1{⟩}^{\otimes n})$
Circuit	1 Hadamard + $n-1$ CNOTs, depth $n$
Measurement	All-zeros or all-ones with 50% each
Key property	Genuine multipartite entanglement; all entanglement lost on single-qubit loss
Verification	Population test, parity test, coherence test (X-basis)
Mermin violation	$⟨M⟩=4$ for 3 qubits (classical bound: 2)
Noise scaling	Fidelity decays as $(1-4p/3{)}^n$ under depolarizing noise; practical limit ~8 qubits at 2% error
Fragility	Decoherence rate scales as $n\gamma$; exponential sensitivity to qubit count
Applications	Quantum networking, metrology, secret sharing, error correction

See Also

• Bell State Preparation -- The 2-qubit special case of GHZ states
• Noise Characterization -- Detailed noise modeling techniques
• Quantum Teleportation -- Another protocol relying on entanglement