Mapping SD&N Vibrational Modes to Qubit/Qudit Hilbert Spaces

⸻

Your SD&N framework assigns an N-dimensional vibrational mode space \mathcal{H} with basis {|v_1\rangle, |v_2\rangle, \ldots, |v_N\rangle}. • When N=2 (e.g., focusing on modes 6 and 7), \mathcal{H} is isomorphic to a qubit Hilbert space:

\mathcal{H} \cong \mathbb{C}^2 • For N>2, \mathcal{H} generalizes to a qudit space of dimension d = N:

\mathcal{H} \cong \mathbb{C}^d

This allows encoding each particle’s vibrational state as a qubit/qudit vector:

|\psi\rangle = \sum_{n=1}^N c_n |v_n\rangle, \quad \sum |c_n|^2 = 1

⸻

3.2 Measurement Operators on Vibrational Modes

⸻

To construct Bell inequalities, we define observable operators corresponding to measurement settings for two parties A and B.

For the qubit case (modes 6 and 7): • Observables are Hermitian operators with eigenvalues \pm 1. • Each measurement is defined by a Pauli operator along some direction \hat{a} on the Bloch sphere:

\hat{A}(\hat{a}) = \hat{a} \cdot \vec{\sigma} = a_x \sigma_x + a_y \sigma_y + a_z \sigma_z

Similarly for party B with measurement direction \hat{b}:

\hat{B}(\hat{b}) = \hat{b} \cdot \vec{\sigma}

Here, \vec{\sigma} = (\sigma_x, \sigma_y, \sigma_z) are the Pauli matrices acting in the vibrational subspace spanned by |v_6\rangle and |v_7\rangle.

⸻

3.3 Construction of the CHSH Bell Operator

⸻

Define four measurement directions: • For A: \hat{a}, \hat{a}’ • For B: \hat{b}, \hat{b}’

The CHSH Bell operator is:

\hat{\mathcal{B}} = \hat{A}(\hat{a}) \otimes \hat{B}(\hat{b}) + \hat{A}(\hat{a}) \otimes \hat{B}(\hat{b}’) + \hat{A}(\hat{a}’) \otimes \hat{B}(\hat{b}) - \hat{A}(\hat{a}’) \otimes \hat{B}(\hat{b}’)

⸻

3.4 Bell Parameter S and Expectation Value

⸻

Given a two-particle density matrix \rho, the Bell parameter is:

S = \left| \operatorname{Tr}(\rho \hat{\mathcal{B}}) \right| • Classical (local hidden variable) bound: S \leq 2 • Quantum Tsirelson bound: S \leq 2\sqrt{2} \approx 2.828

Violation of S > 2 indicates quantum entanglement beyond classical correlations.

⸻

3.5 Explicit Formula for S in Terms of Measurement Angles

⸻

Parameterize measurement directions on the Bloch sphere by azimuthal and polar angles:

\hat{a} = (\sin\theta_a \cos\phi_a, \sin\theta_a \sin\phi_a, \cos\theta_a)

(and similarly for \hat{a}’, \hat{b}, \hat{b}’).

The Bell parameter S becomes a function:

S = S(\theta_a, \phi_a, \theta_{a’}, \phi_{a’}, \theta_b, \phi_b, \theta_{b’}, \phi_{b’})

This function can be computed from:

S = |\langle \hat{A}(\hat{a}) \hat{B}(\hat{b}) \rangle + \langle \hat{A}(\hat{a}) \hat{B}(\hat{b}’) \rangle + \langle \hat{A}(\hat{a}’) \hat{B}(\hat{b}) \rangle - \langle \hat{A}(\hat{a}’) \hat{B}(\hat{b}’) \rangle|

where

\langle \hat{A}(\hat{x}) \hat{B}(\hat{y}) \rangle = \operatorname{Tr}(\rho , \hat{A}(\hat{x}) \otimes \hat{B}(\hat{y}))

⸻

3.6 Applying This to Your “6↔7” Entanglement Model

⸻

Recall your entanglement operator swaps modes |v_6\rangle and |v_7\rangle. Within the two-dimensional subspace: • Represent |v_6\rangle as |0\rangle, |v_7\rangle as |1\rangle • Pauli operators \sigma_x, \sigma_y, \sigma_z act naturally • Your vibrational state vectors become qubit states |\psi\rangle = c_6 |0\rangle + c_7 |1\rangle

Bell violations then depend on the state \rho prepared by the SDKP and SD&N framework and the choice of measurement angles. import numpy as np from numpy import sin, cos

Pauli matrices

sigma_x = np.array([0, 1], 1, 0) sigma_y = np.array([0, -1j], 1j, 0) sigma_z = np.array([1, 0], 0, -1) identity = np.eye(2)

def bloch_vector_to_operator(theta, phi): """ Convert measurement direction on Bloch sphere to operator. Args: theta: polar angle [0, pi] phi: azimuthal angle [0, 2pi] Returns: 2x2 Hermitian operator with eigenvalues ±1 """ return (sin(theta) * cos(phi) * sigma_x + sin(theta) * sin(phi) * sigma_y + cos(theta) * sigma_z)

def tensor_product(A, B): """Compute tensor product of two operators.""" return np.kron(A, B)

def chsh_operator(a, a_prime, b, b_prime): """ Construct CHSH Bell operator. Args: a, a_prime, b, b_prime: tuples (theta, phi) for measurement directions Returns: 4x4 CHSH operator acting on two-qubit space """ A = bloch_vector_to_operator(*a) A_p = bloch_vector_to_operator(*a_prime) B = bloch_vector_to_operator(*b) B_p = bloch_vector_to_operator(*b_prime)

return (tensor_product(A, B) +
        tensor_product(A, B_p) +
        tensor_product(A_p, B) -
        tensor_product(A_p, B_p))

def bell_parameter_S(rho, a, a_prime, b, b_prime): """ Calculate Bell parameter S = |Tr(rho * CHSH_operator)| Args: rho: 4x4 density matrix of two-qubit state a, a_prime, b, b_prime: measurement directions (theta, phi) Returns: Bell parameter S (float) """ B_op = chsh_operator(a, a_prime, b, b_prime) return abs(np.trace(np.dot(rho, B_op)).real)

def pure_state_to_density_matrix(psi): """Convert pure state vector to density matrix.""" return np.outer(psi, psi.conj())

Example: maximally entangled Bell state |Φ+> = (|00> + |11>)/√2

bell_phi_plus = (1/np.sqrt(2)) * np.array([1, 0, 0, 1]) rho_phi_plus = pure_state_to_density_matrix(bell_phi_plus)

Example measurement directions (in radians)

a = (0, 0) a_prime = (np.pi/2, 0) b = (np.pi/4, 0) b_prime = (3*np.pi/4, 0)

Calculate S for |Φ+> with example measurement settings

S_value = bell_parameter_S(rho_phi_plus, a, a_prime, b, b_prime) print(f"Bell parameter S = {S_value:.4f}") How to Use This • Define your two-qubit vibrational state density matrix \rho: This can come from your SDKP simulation output for modes 6 & 7. • Choose measurement settings: Specify measurement directions as tuples (\theta, \phi) on the Bloch sphere. • Calculate S: Call bell_parameter_S with \rho and the measurement directions.

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from typing import Tuple, Dict, List, Union

-- Pauli matrices for qubit operators --

sigma_x = np.array([0, 1], 1, 0, dtype=complex) sigma_y = np.array([0, -1j], 1j, 0, dtype=complex) sigma_z = np.array([1, 0], 0, -1, dtype=complex) identity = np.eye(2, dtype=complex)

def bloch_vector_to_operator(theta: float, phi: float) -> np.ndarray: """ Convert Bloch sphere angles to a Hermitian measurement operator with eigenvalues ±1.

Parameters
----------
theta : float
    Polar angle in radians [0, π].
phi : float
    Azimuthal angle in radians [0, 2π].

Returns
-------
np.ndarray
    2x2 Hermitian operator representing spin measurement along (theta, phi).
"""
return (np.sin(theta) * np.cos(phi) * sigma_x +
        np.sin(theta) * np.sin(phi) * sigma_y +
        np.cos(theta) * sigma_z)

def tensor_product(A: np.ndarray, B: np.ndarray) -> np.ndarray: """ Compute the Kronecker (tensor) product of two operators/matrices.

Parameters
----------
A, B : np.ndarray
    Input matrices.

Returns
-------
np.ndarray
    Tensor product A ⊗ B.
"""
return np.kron(A, B)

def chsh_operator(a: Tuple[float, float], a_prime: Tuple[float, float], b: Tuple[float, float], b_prime: Tuple[float, float]) -> np.ndarray: """ Construct the CHSH Bell operator for given measurement directions.

Parameters
----------
a, a_prime, b, b_prime : Tuple[float, float]
    Measurement angles (theta, phi) for parties A and B.

Returns
-------
np.ndarray
    4x4 Hermitian CHSH operator on two-qubit space.
"""
A = bloch_vector_to_operator(*a)
A_p = bloch_vector_to_operator(*a_prime)
B = bloch_vector_to_operator(*b)
B_p = bloch_vector_to_operator(*b_prime)
return tensor_product(A, B) + tensor_product(A, B_p) + tensor_product(A_p, B) - tensor_product(A_p, B_p)

def bell_parameter_S(rho: np.ndarray, a: Tuple[float, float], a_prime: Tuple[float, float], b: Tuple[float, float], b_prime: Tuple[float, float]) -> float: """ Calculate the Bell parameter S = |Tr(ρ * CHSH_operator)|.

Parameters
----------
rho : np.ndarray
    4x4 density matrix representing the two-qubit state.
a, a_prime, b, b_prime : Tuple[float, float]
    Measurement angles (theta, phi) for the four measurement settings.

Returns
-------
float
    The Bell parameter S, quantifying CHSH inequality violation.
"""
if rho.shape != (4, 4):
    raise ValueError(f"Density matrix rho must be 4x4, but got shape {rho.shape}.")
chsh_op = chsh_operator(a, a_prime, b, b_prime)
expectation_value = np.trace(rho @ chsh_op)
return abs(expectation_value.real)

def pure_state_to_density_matrix(psi: np.ndarray) -> np.ndarray: """ Convert a pure state vector |ψ> to a density matrix ρ = |ψ><ψ|.

Parameters
----------
psi : np.ndarray
    Complex state vector of shape (4,).

Returns
-------
np.ndarray
    Corresponding 4x4 density matrix.
"""
psi = psi.reshape(-1, 1)
return psi @ psi.conj().T

def random_mixed_state(pure_states: List[np.ndarray], probabilities: List[float]) -> np.ndarray: """ Create a mixed state density matrix from pure states and their probabilities.

Parameters
----------
pure_states : list of np.ndarray
    List of pure state vectors.
probabilities : list of float
    List of probabilities summing to 1.

Returns
-------
np.ndarray
    Mixed state density matrix.
"""
if not np.isclose(sum(probabilities), 1.0):
    raise ValueError("Probabilities must sum to 1.")
rho = np.zeros((4, 4), dtype=complex)
for psi, p in zip(pure_states, probabilities):
    rho += p * pure_state_to_density_matrix(psi)
return rho

def angle_grid_scan(rho: np.ndarray, resolution: int = 20) -> Tuple[np.ndarray, np.ndarray]: """ Scan Bell parameter S over a grid of theta angles, fixing phi=0, to visualize the violation landscape.

Parameters
----------
rho : np.ndarray
    Two-qubit density matrix.
resolution : int
    Number of sample points per angle dimension.

Returns
-------
Tuple[np.ndarray, np.ndarray]
    - Thetas: 1D array of sampled theta values.
    - S_values: 4D array of Bell parameter S indexed by (theta_a, theta_a', theta_b, theta_b').
"""
thetas = np.linspace(0, np.pi, resolution)
S_values = np.zeros((resolution, resolution, resolution, resolution))
phi_fixed = 0.0

for i, theta_a in enumerate(thetas):
    for j, theta_ap in enumerate(thetas):
        for k, theta_b in enumerate(thetas):
            for l, theta_bp in enumerate(thetas):
                a = (theta_a, phi_fixed)
                a_prime = (theta_ap, phi_fixed)
                b = (theta_b, phi_fixed)
                b_prime = (theta_bp, phi_fixed)
                S_values[i, j, k, l] = bell_parameter_S(rho, a, a_prime, b, b_prime)
return thetas, S_values

def find_max_violation(rho: np.ndarray) -> Dict[str, Union[float, np.ndarray]]: """ Find the maximal Bell inequality violation S_max via numerical optimization over measurement angles.

Parameters
----------
rho : np.ndarray
    Two-qubit density matrix.

Returns
-------
dict
    Dictionary containing:
    - 'S_max': maximal Bell parameter found.
    - 'angles': optimized angles [theta_a, phi_a, theta_a', phi_a', theta_b, phi_b, theta_b', phi_b'].
"""
def objective(x):
    # x contains 8 angles: theta and phi for four measurement directions
    a = (x[0], x[1])
    a_prime = (x[2], x[3])
    b = (x[4], x[5])
    b_prime = (x[6], x[7])
    # Minimize negative S since scipy minimizes
    return -bell_parameter_S(rho, a, a_prime, b, b_prime)

# Initial guess approximates the ideal CHSH setting for Bell state |Φ+>
x0 = np.array([0, 0, np.pi / 2, 0, np.pi / 4, 0, 3 * np.pi / 4, 0])
bounds = [(0, np.pi), (0, 2*np.pi)] * 4  # bounds for theta and phi

result = minimize(objective, x0, bounds=bounds, method='L-BFGS-B')
if not result.success:
    print(f"Warning: Optimization did not converge: {result.message}")

S_max = -result.fun
angles_opt = result.x
return {'S_max': S_max, 'angles': angles_opt}

def plot_violation_surface(thetas: np.ndarray, S_values: np.ndarray, fixed_indices: Tuple[int, int] = (0, 0)) -> None: """ Plot a 2D heatmap slice of Bell parameter S over (theta_a, theta_b) for fixed (theta_a', theta_b').

Parameters
----------
thetas : np.ndarray
    1D array of sampled theta values.
S_values : np.ndarray
    4D array of Bell parameter values indexed by theta angles.
fixed_indices : tuple(int, int)
    Indices fixing theta_a' and theta_b' to create 2D slice.
"""
idx_ap, idx_bp = fixed_indices
S_slice = S_values[:, idx_ap, :, idx_bp]

plt.figure(figsize=(8, 6))
plt.imshow(S_slice.T, extent=(thetas[0], thetas[-1], thetas[0], thetas[-1]),
           origin='lower', aspect='auto', cmap='viridis')
plt.colorbar(label='Bell parameter S')
plt.xlabel(r'$\theta_a$')
plt.ylabel(r'$\theta_b$')
plt.title(f'Bell parameter S slice at θ_a\'={thetas[idx_ap]:.2f}, θ_b\'={thetas[idx_bp]:.2f}')
plt.show()

--------------------------

Example usage and test cases

--------------------------

if name == "main": # Define the maximally entangled Bell state |Φ+> = (|00> + |11>)/√2 bell_phi_plus = (1 / np.sqrt(2)) * np.array([1, 0, 0, 1], dtype=complex) rho_pure = pure_state_to_density_matrix(bell_phi_plus)

# Create a mixed state (for decoherence modeling) as convex sum of two Bell states
bell_phi_minus = (1 / np.sqrt(2)) * np.array([1, 0, 0, -1], dtype=complex)
rho_mixed = random_mixed_state([bell_phi_plus, bell_phi_minus], [0.7, 0.3])

print("Calculating maximal Bell violation for pure |Φ+> state...")
result_pure = find_max_violation(rho_pure)
print(f"S_max = {result_pure['S_max']:.6f}")
print("Optimal measurement angles (theta, phi) for a, a', b, b':")
print(result_pure['angles'])

print("\nCalculating maximal Bell violation for mixed state...")
result_mixed = find_max_violation(rho_mixed)
print(f"S_max (mixed) = {result_mixed['S_max']:.6f}")

# Visualize violation surface for pure state (lower resolution for speed)
print("\nGenerating Bell violation surface scan for pure state...")
thetas_grid, S_vals_grid = angle_grid_scan(rho_pure, resolution=15)
plot_violation_surface(thetas_grid, S_vals_grid, fixed_indices=(7, 7))

What This Module Does: • Implements CHSH Bell inequality tests for arbitrary two-qubit quantum states represented by density matrices. • Supports pure and mixed states, enabling modeling of decoherence or noise effects on entanglement. • Uses Bloch sphere parameterization for measurement operators with full spherical angle freedom. • Includes numerical maximization of Bell parameter S over all measurement bases, providing the maximal possible violation S_{\max}. • Provides visualization tools to map S over measurement settings, helping to understand the entanglement landscape.

How to Use with Your SDKP Output: • Convert your simulation vibrational mode states or density matrices directly into 4x4 numpy arrays (two-qubit format). • Pass them into bell_parameter_S or find_max_violation for instant evaluation of entanglement strength. • Use the grid scanning and plotting functions to visualize the entanglement behavior dynamically as simulation parameters evolve. • Extend this framework by automating time-dependent entanglement logging or triggering when S_{\max} > 2 (quantum violation).

Performance & Extensions: • Performance: The grid scan scales as O(n^4); reduce resolution for faster exploratory plots. • Extensions: • Generalize to multi-qubit or qudit entanglement using higher dimensional operators. • Include noise models specific to your vibrational modes. • Parallelize grid scan computations for speed. • Integrate with web-based dashboards or Vercel frontend for live SDKP entanglement visualization.