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

flowchart TD
    A["Step 1: Initialize n qubits in |0>"] --> B["Step 2: Apply Hadamard to qubit 0"]
    B --> C["Step 3: Apply CNOT(0,1)"]
    C --> D["Step 4: Apply CNOT(0,2)"]
    D --> E["... Apply CNOT(0,k) for k = 1 to n-1"]
    E --> F["Step 5: Measure all qubits"]

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)

graph LR
    subgraph "3-qubit GHZ Circuit"
        Q0["q[0]: |0> ---[H]---*-------*---"]
        Q1["q[1]: |0> ----------|-[+]---|---"]
        Q2["q[2]: |0> ------------------[+]-"]
    end

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

运行结果：

3-qubit GHZ measurement counts:
  |000> : 5089
  |111> : 4911

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}")

运行结果：

n=3: GHZ fidelity estimate = 1.0000
n=5: GHZ fidelity estimate = 1.0000
n=8: GHZ fidelity estimate = 1.0000

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(3) Verification (20000 shots) ---
  |000> : 10015 (50.1%)
  |111> : 9985 (49.9%)
  Other outcomes : 0 (0.00%)
  PASS: State is consistent with GHZ.
--- GHZ(4) Verification (20000 shots) ---
  |0000> : 10095 (50.5%)
  |1111> : 9905 (49.5%)
  Other outcomes : 0 (0.00%)
  PASS: State is consistent with GHZ.
--- GHZ(5) Verification (20000 shots) ---
  |00000> : 9908 (49.5%)
  |11111> : 10092 (50.5%)
  Other outcomes : 0 (0.00%)
  PASS: State is consistent with GHZ.
--- GHZ(6) Verification (20000 shots) ---
  |000000> : 10034 (50.2%)
  |111111> : 9966 (49.8%)
  Other outcomes : 0 (0.00%)
  PASS: State is consistent with GHZ.

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)")

运行结果：

Exact probabilities (3-qubit GHZ):
  |000>: 0.500000
  |111>: 0.500000
Density matrix shape: (8, 8)
Purity Tr(rho^2): 1.000000
Reduced density matrix (qubit 0):
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]]
Reduced purity: 0.500000 (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)

运行结果：

--- GHZ(3) Noise Comparison (10000 shots) ---
  Ideal: {'000': 5024, '111': 4976}
  Noisy: {'000': 4761, '001': 51, '010': 197, '011': 57, '100': 41, '101': 198, '110': 52, '111': 4643}
  GHZ outcomes: 9404/10000 (94.0%)
--- GHZ(5) Noise Comparison (10000 shots) ---
  Ideal: {'00000': 4966, '11111': 5034}
  Noisy: {'00000': 4105, '00001': 38, '00010': 332, '00011': 2, '00100': 137, '00101': 1, '00110': 139, '00111': 1, '01000': 85, '01001': 2, '01010': 10, '01011': 3, '01100': 8, '01101': 5, '01110': 68, '01111': 42, '10000': 41, '10001': 77, '10010': 2, '10011': 9, '10100': 3, '10101': 11, '10110': 2, '10111': 80, '11000': 3, '11001': 139, '11010': 1, '11011': 136, '11100': 4, '11101': 341, '11110': 34, '11111': 4139}
  GHZ outcomes: 8244/10000 (82.4%)

Comparing Ideal vs. Noisy GHZ Correlations

The ideal GHZ state can be analyzed exactly via the density matrix simulator. Under noise, the CPUQVM with shots provides an approximate comparison:

from pyqpanda3 import core
import numpy as np


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

    # Ideal (exact via density matrix)
    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 (sampled via CPUQVM)
    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]
        )

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

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

    # Print comparison
    dim = 2 ** n
    print(f"\n--- GHZ({n}) Probability Comparison (ideal exact vs noisy {shots} shots) ---")
    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')
        noisy_prob = counts.get(bs, 0) / shots
        if probs_ideal[i] > 1e-6 or noisy_prob > 0:
            d = abs(probs_ideal[i] - noisy_prob)
            print(f"|{bs}> {probs_ideal[i]:10.6f} {noisy_prob:10.6f} {d:10.6f}")

    print(f"\nIdeal purity:  {purity_ideal:.6f}")


compare_ideal_noisy_ghz(3)
compare_ideal_noisy_ghz(4)

运行结果：

--- GHZ(3) Probability Comparison (ideal exact vs noisy 50000 shots) ---
Basis      Ideal      Noisy      Delta
--------------------------------------
|000>   0.500000   0.471260   0.028740
|001>   0.000000   0.004740   0.004740
|010>   0.000000   0.018980   0.018980
|011>   0.000000   0.004600   0.004600
|100>   0.000000   0.004420   0.004420
|101>   0.000000   0.019840   0.019840
|110>   0.000000   0.004620   0.004620
|111>   0.500000   0.471540   0.028460

Ideal purity:  1.000000

--- GHZ(4) Probability Comparison (ideal exact vs noisy 50000 shots) ---
Basis       Ideal      Noisy      Delta
---------------------------------------
|0000>   0.500000   0.440620   0.059380
|0001>   0.000000   0.004520   0.004520
|0010>   0.000000   0.025600   0.025600
|0011>   0.000000   0.000280   0.000280
|0100>   0.000000   0.009860   0.009860
|0101>   0.000000   0.000520   0.000520
|0110>   0.000000   0.009480   0.009480
|0111>   0.000000   0.004720   0.004720
|1000>   0.000000   0.004580   0.004580
|1001>   0.000000   0.009480   0.009480
|1010>   0.000000   0.000260   0.000260
|1011>   0.000000   0.010140   0.010140
|1100>   0.000000   0.000400   0.000400
|1101>   0.000000   0.027920   0.027920
|1110>   0.000000   0.005100   0.005100
|1111>   0.500000   0.446520   0.053480

Ideal purity:  1.000000

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_sampled(n: int, noise: core.NoiseModel, shots: int = 20000) -> float:
    """Estimate GHZ fidelity by sampling |0...0> and |1...1> outcomes."""
    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


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_sampled(n, noise)
    print(f"{p:>12.3f}  {fid:>10.6f}")

运行结果：

GHZ(4) fidelity vs. depolarizing error strength:
  Error rate    Fidelity
-------------------------
       0.000    1.000000
       0.005    0.968300
       0.010    0.942250
       0.020    0.886950
       0.050    0.748150
       0.100    0.566000

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: 1.0})
        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()

运行结果：

  <X0 X1 X2> = 1.000000, contribution = 1.000000
  <X0 Y1 Y2> = -1.000000, contribution = 1.000000
  <Y0 X1 Y2> = -1.000000, contribution = 1.000000
  <Y0 Y1 X2> = -1.000000, contribution = 1.000000

Mermin operator <M> = 4.000000
Classical bound: |<M>| <= 2
Quantum prediction: <M> = 4
RESULT: Mermin inequality VIOLATED!

---

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.

sequenceDiagram
    participant State as Quantum State
    participant H as Hadamard(0)
    participant CX1 as CNOT(0,1)
    participant CX2 as CNOT(0,2)

    Note over State: |000> = [1, 0, 0, 0, 0, 0, 0, 0]
    State->>H: Apply H on qubit 0
    Note over H: H|0> = (|0>+|1>)/sqrt(2)
    H->>CX1: (|000>+|100>)/sqrt(2)
    Note over CX1: |000> stays |000>, |100> becomes |110>
    CX1->>CX2: (|000>+|110>)/sqrt(2)
    Note over CX2: |000> stays |000>, |110> becomes |111>
    CX2->>State: (|000>+|111>)/sqrt(2) = |GHZ_3>
    Note over State: Genuine 3-partite entanglement

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:

graph TB
    subgraph "Bell State (2 qubits)"
        B0["q[0]"] ---|"maximally entangled"| B1["q[1]"]
    end
    subgraph "GHZ State (n qubits)"
        G0["q[0]"] --- G1["q[1]"]
        G0 --- G2["q[2]"]
        G0 -.-> Gk["q[n-1]"]
    end

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.

graph LR
    subgraph "Decoherence Impact on GHZ State"
        A["Ideal GHZ: Only |0...0> and |1...1>, 50/50"] -->|"Phase noise"| B["Noisy GHZ: Mixed outcomes, Purity < 1"]
        B -->|"Strong noise"| C["Fully mixed: No entanglement"]
    end

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)

运行结果：

==================================================
Full GHZ(3) Verification (20000 shots)
==================================================
  Population test:
    |000>: 0.4960  (expected 0.5000)
    |111>: 0.5040  (expected 0.5000)
    Other:        0.0000  (expected 0.0000)
    Result: PASS
  Parity test (adjacent qubit XOR):
    q[0] XOR q[1] = 0: 1.0000  (expected 1.0000)
    q[1] XOR q[2] = 0: 1.0000  (expected 1.0000)
    Result: PASS (0 violations)
  Coherence test (X-basis even parity):
    Even parity fraction: 1.0000  (expected 1.0000)
    Result: PASS

  Overall: ALL TESTS PASSED

==================================================
Full GHZ(4) Verification (20000 shots)
==================================================
  Population test:
    |0000>: 0.5019  (expected 0.5000)
    |1111>: 0.4981  (expected 0.5000)
    Other:        0.0000  (expected 0.0000)
    Result: PASS
  Parity test (adjacent qubit XOR):
    q[0] XOR q[1] = 0: 1.0000  (expected 1.0000)
    q[1] XOR q[2] = 0: 1.0000  (expected 1.0000)
    q[2] XOR q[3] = 0: 1.0000  (expected 1.0000)
    Result: PASS (0 violations)
  Coherence test (X-basis even parity):
    Even parity fraction: 1.0000  (expected 1.0000)
    Result: PASS

  Overall: ALL TESTS PASSED

==================================================
Full GHZ(5) Verification (20000 shots)
==================================================
  Population test:
    |00000>: 0.5037  (expected 0.5000)
    |11111>: 0.4964  (expected 0.5000)
    Other:        -0.0000  (expected 0.0000)
    Result: PASS
  Parity test (adjacent qubit XOR):
    q[0] XOR q[1] = 0: 1.0000  (expected 1.0000)
    q[1] XOR q[2] = 0: 1.0000  (expected 1.0000)
    q[2] XOR q[3] = 0: 1.0000  (expected 1.0000)
    q[3] XOR q[4] = 0: 1.0000  (expected 1.0000)
    Result: PASS (0 violations)
  Coherence test (X-basis even parity):
    Even parity fraction: 1.0000  (expected 1.0000)
    Result: PASS

  Overall: ALL TESTS PASSED

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


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).dagger()  # S^\dagger for correct Y-basis rotation
            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)

运行结果：

==================================================
IDEAL (noiseless) Mermin test
==================================================
Mermin Inequality Test (3-qubit GHZ, 50000 shots)
    Term   Expectation   Contribution
----------------------------------------
     XXX      1.000000       1.000000
     XYY     -1.000000       1.000000
     YXY     -1.000000       1.000000
     YYX     -1.000000       1.000000

Mermin operator <M> = 4.000000
Classical bound:     |<M>| <= 2
Quantum prediction:  <M> = 4
Violation ratio:     2.0000
RESULT: Mermin inequality VIOLATED (|<M>| = 4.0000 > 2)

==================================================
NOISY Mermin test (1% depolarizing)
==================================================
Mermin Inequality Test (3-qubit GHZ, 50000 shots)
    Term   Expectation   Contribution
----------------------------------------
     XXX      0.860200       0.860200
     XYY     -0.860320       0.860320
     YXY     -0.858320       0.858320
     YYX     -0.858160       0.858160

Mermin operator <M> = 3.437000
Classical bound:     |<M>| <= 2
Quantum prediction:  <M> = 4
Violation ratio:     1.7185
RESULT: Mermin inequality VIOLATED (|<M>| = 3.4370 > 2)

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


# 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}  {'Sampled F':>10s}  {'Theoretical':>12s}")
print("-" * 45)

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]
        )

    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}  {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}  {'Fidelity':>10s}")
print("-" * 27)

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]
            )

    fid = ghz_fidelity(n, noise)
    print(f"{p:>12.3f}  {fid:>10.6f}")

运行结果：

GHZ Fidelity vs. Qubit Count (depolarizing error = 2%)
 Qubits   CNOTs   Sampled F   Theoretical
---------------------------------------------
      2       1      0.9811      0.947378
      3       2      0.9420      0.922114
      4       3      0.8877      0.897525
      5       4      0.8200      0.873591
      6       5      0.7442      0.850295
      7       6      0.6631      0.827620
      8       7      0.5776      0.805551


GHZ(4) Fidelity vs. Noise Strength
  Error rate    Fidelity
---------------------------
       0.000    1.000000
       0.005    0.975200
       0.010    0.941500
       0.020    0.885800
       0.030    0.836000
       0.050    0.740700
       0.080    0.631800
       0.100    0.564700
       0.150    0.435100
       0.200    0.355200

Explanation

The sweep reveals several important patterns:

Qubits	CNOTs	Fidelity at 2% error
2	1	~0.98
3	2	~0.94
4	3	~0.89
5	4	~0.82
6	5	~0.74
7	6	~0.66
8	7	~0.58

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 98% to roughly 58%.

2. Noise accumulation beyond theory: The sampled fidelity falls faster than the simple theoretical estimate $(1-4p/3{)}^n$ for larger $n$, because CNOT gates couple noise between qubits in ways that single-qubit approximations do not capture.

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