There are effectively two main issues that explain why your CHSH value sometimes reaches 4 (in the case without Eve) or goes below/above the expected value (in the case with Eve):
You are not implementing the ‘classic’ CHSH measurement settings
In the standard CHSH setup, one (entangled) Bell pair is used, and Alice has two possible measurement bases $\{A_0, A_1\}$, while Bob has two possible measurement bases $\{B_0, B_1\}$. Concretely, for a maximal quantum mechanical CHSH violation (Tsirelson bound $(2\sqrt{2})$, the angles typically chosen (for example) are:
- $A_0 = 0$, $A_1 = \pi/4$
- $B_0 = \pi/8$, $B_1 = -\pi/8$
In your code, however, you choose the bases (looking at the angles array) from the sets
and then completely randomly $((\mathrm{randint}(3)))$. This is not the well‐known CHSH protocol setup, where typically only two angles per party are needed, and certainly not the pair of ‘optimal’ angles that yield $2\sqrt{2}$.
- If Alice and Bob happen by chance to choose the same angle (or angles that produce perfect correlation for the Bell or GHZ state you prepared), the result in the formula for $S$ can be very strongly correlated. In the extreme, this can lead to $|S|=4$ (this is the mathematically maximum value, if all measurements “perfectly” reinforce each other).
- In other runs, especially when a third qubit (Eve) is also involved, that extra measurement disturbs the total quantum superposition, leading to strongly varying outcomes for $S$. This can easily fall below the classical limit of 2 or slightly above it, depending on chance and the randomly assigned basis label.
You are using (with ‘WITH_EVE=False’) a non‐standard state and measurement sequence
The code does create a Bell‐ or GHZ‐like state:
circuits[i].h(qr[A])
circuits[i].cx(qr[A], qr[B])
if WITH_EVE:
circuits[i].cx(qr[B], qr[E])
But immediately afterward, you apply $\mathrm{RY}(\theta)$ rotations to all qubits (including Alice and Bob) and then measure in the computational basis. For a standard CHSH test, you would normally want:
- Exactly two measurement settings for Alice (A0, A1) and two for Bob (B0, B1),
- A singlet/fully entangled state (e.g., $\tfrac{1}{\sqrt{2}}(\lvert 01\rangle - \lvert 10\rangle)$),
- Randomly choose whether you use $\{A_0\}$ or $\{A_1\}$ (and for Bob $\{B_0\}$ or $\{B_1\}$) on each measurement,
- Then compute the correlations specifically for the four combinations $(A_0, B_0)$, $(A_0, B_1)$, $(A_1, B_0)$, $(A_1, B_1)$.
In your code, for each ‘shot’ the angles $\theta$ and the corresponding basis label $\{0,1,2\}$ are chosen completely randomly, and then you derive a CHSH‐type correlation from that. This does not produce a standard CHSH measurement, but rather a mix of many angles. Hence you get a very erratic result, which in some cases can go up to 4 (because you happened to choose angles that perfectly correlate the measurements), or in the presence of Eve can go off in any direction.
How can you get $2\sqrt{2}$?
Create a pure Bell state
The simplest is $\frac{1}{\sqrt{2}}(\lvert 00\rangle + \lvert 11\rangle)$ or the singlet $\frac{1}{\sqrt{2}}(\lvert 01\rangle - \lvert 10\rangle)$. In Qiskit, for example, you can do:
qc.h(qr[A])
qc.cx(qr[A], qr[B])
to create $\frac{1}{\sqrt{2}}(\lvert 00\rangle + \lvert 11\rangle)$ (if both qubits start in $\lvert 0\rangle$).
Apply the CHSH basis rotations
Instead of three possible $\mathrm{RY}$ angles, pick for each measurement on Alice either $\theta = 0$ (or $\pi/2$, depending on the convention) or $\theta = \pi/4$. For Bob, pick either $\theta = \pi/8$ or $\theta = -\pi/8$.
Why do you often get < 2 (or sometimes just above) with Eve?
As soon as an additional party (Eve) measures the same entangled qubits, the correlations between Alice and Bob become disturbed. In a genuine CHSH setup, there are only two parties. If Eve also carries out (partial) measurements, the quantum state of the system changes (collapse, decoherence). As a result, the nice two‐party correlation disappears or is disturbed in unpredictable ways. Hence, by choosing random angles and having a third qubit measured, you get strong fluctuations in $S$.
With the modified code S reached 2.47, so not the maximum value, could be due to noise.
import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_aer import AerSimulator
First we define the angles:
Alice: A0 = 0, A1 = pi/4
Bob: B0 = pi/8, B1 = -pi/8
A0 = 0.0
A1 = np.pi/4
B0 = np.pi/8
B1 = -np.pi/8
Eve chooses to measure Bob's qubit in one of Bob's two angles:
This means that Eve uses the same measurement angles as Bob
E0 = B0
E1 = B1
(Optionally, you can force Eve to use only one angle or random. Here we choose random.)
If EVE_ACTIVE = True, then Eve is active.
EVE_ACTIVE = False
We perform 1000 runs and gather data
N_RUNS = 1000
Measurement of qubit output maps |0> -> +1, |1> -> -1 in correlation calculations
def bit_to_pm1(bit_value: int) -> int:
return +1 if bit_value == 0 else -1
Data structure for storing results:
results[(a_setting, b_setting)] = list of (alice_out, bob_out)
results = {
('A0', 'B0'): [],
('A0', 'B1'): [],
('A1', 'B0'): [],
('A1', 'B1'): []
}
backend = AerSimulator()
for _ in range(N_RUNS):
# Randomly choose Alice's measurement basis: A0 or A1
alice_choice = np.random.choice(['A0','A1'])
# Randomly choose Bob's measurement basis: B0 or B1
bob_choice = np.random.choice(['B0','B1'])
# Build a circuit
# - qr[0] = Alice's qubit
# - qr[1] = Bob's qubit
qr = QuantumRegister(2, 'q')
cr = ClassicalRegister(2, 'c')
qc = QuantumCircuit(qr, cr)
# 1. Create Bell state (|Phi+>) = (|00> + |11>)/sqrt(2)
qc.h(qr[0])
qc.cx(qr[0], qr[1])
# 2. Eve intercepts Bob's qubit (optional):
if EVE_ACTIVE:
# Randomly choose Eve's measurement basis: E0 or E1
eve_choice = np.random.choice(['E0','E1'])
eve_angle = E0 if eve_choice == 'E0' else E1
# Eve's rotation and measurement:
# - Rotate Bob's qubit into Eve's basis
qc.ry(-eve_angle, qr[1]) # rotating by -theta => measure in "theta" basis
qc.measure(qr[1], cr[1]) # measure Bob's qubit in Eve's basis
# reset Bob's qubit and set it back to the post-measurement state
qc.reset(qr[1])
# c_if (conditional gates) is one approach, but here we use transpile tricks for simplicity.
# In standard Qiskit, it is a bit more complicated to precisely restore the measured state.
# Simpler is to ignore the measured outcome and leave Bob's qubit in, for example, |0>.
#
# => This, however, reveals some info about the realism of Eve's action...
# In a real eavesdropping, Eve would theoretically know the measured outcome
# and then forward exactly that state.
#
# For the core of the CHSH-protocol, this difference is not very important:
# the fact is that Eve's measurement 'disturbs' the entanglement.
#
# Here we set Bob's qubit to |0> after Eve's measurement and reset.
# 3. Alice measurement in A0 = 0 or A1 = pi/4:
if alice_choice == 'A0':
angle_a = A0
else:
angle_a = A1
qc.ry(-angle_a, qr[0]) # rotate qubit 0 from Alice's basis to Z-basis
# 4. Bob measurement in B0 = pi/8 or B1 = -pi/8:
if bob_choice == 'B0':
angle_b = B0
else:
angle_b = B1
qc.ry(-angle_b, qr[1]) # rotate qubit 1 from Bob's basis to Z-basis
# 5. Measure in the computational basis
qc.measure(qr[0], cr[0]) # Alice's outcome in cr[0]
qc.measure(qr[1], cr[1]) # Bob's outcome in cr[1]
# Simulate
transpiled_qc = transpile(qc, backend)
result = backend.run(transpiled_qc, shots=1).result()
counts = result.get_counts()
# There is only one outcome (shots=1), so take the only key
measured_state = list(counts.keys())[0]
# measured_state is, for example, '10' => cr[1]cr[0], so reversed indexing
# Note: the classical register index c[1] is Bob, c[0] is Alice
bob_bit = int(measured_state[0]) # MSB
alice_bit = int(measured_state[1]) # LSB
# In CHSH correlation: 0 -> +1, 1 -> -1
a_val = bit_to_pm1(alice_bit)
b_val = bit_to_pm1(bob_bit)
# Save in the results dict
results[(alice_choice, bob_choice)].append((a_val, b_val))
Now the results contain the outcomes per (A_i, B_j)
Function to compute E(A,B) = <a*b> from the measured data
def correlation(list_of_pairs):
# list_of_pairs contains tuples (a_val, b_val) where each is in {+1, -1}
if len(list_of_pairs) == 0:
return 0
return sum(a*b for (a,b) in list_of_pairs) / len(list_of_pairs)
Compute E(A0,B0), E(A0,B1), E(A1,B0), E(A1,B1)
E_A0B0 = correlation(results[('A0','B0')])
E_A0B1 = correlation(results[('A0','B1')])
E_A1B0 = correlation(results[('A1','B0')])
E_A1B1 = correlation(results[('A1','B1')])
CHSH value:
Note, the 'classical' form is S = E(A0,B0) - E(A0,B1) + E(A1,B0) + E(A1,B1)
or small variants thereof, depending on the sign convention.
Note that in some sources -E(A1,B1) is used. There are different definitions,
so you can check what you want.
Here we use a common convention:
S = (E_A0B0 + E_A0B1 + E_A1B0 - E_A1B1)
Print results
print("\n=== CHSH RESULTS ===")
print(f"E(A0, B0) = {E_A0B0:.3f}")
print(f"E(A0, B1) = {E_A0B1:.3f}")
print(f"E(A1, B0) = {E_A1B0:.3f}")
print(f"E(A1, B1) = {E_A1B1:.3f}")
print(f"S = {S:.3f}")
print("\nEVE_ACTIVE =", EVE_ACTIVE)
print("Number of runs:", N_RUNS)
Ouput
=== CHSH RESULTS ===
E(A0, B0) = 0.919
E(A0, B1) = 0.941
E(A1, B0) = 0.949
E(A1, B1) = 0.333
S = 2.476
EVE_ACTIVE = False
Number of runs: 1000
```
