Quantum Spintronics

2.1 Spin Qubit Fundamentals

A spin qubit uses the quantum spin state of an electron (or nucleus) as the basis for quantum information processing. The two-level nature of spin-1/2 systems makes them natural candidates for quantum bits.

Bloch Sphere Representation

Any pure qubit state can be visualized as a point on the Bloch sphere:

where $\theta \in [0, \pi]$ is the polar angle and $\phi \in [0, 2\pi)$ is the azimuthal angle. The Bloch vector $\mathbf{r} = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$ points to the state.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def bloch_sphere_visualization(states, labels=None):
    """
    Visualize qubit states on the Bloch sphere

    Parameters:
    -----------
    states : list of tuples, [(theta1, phi1), (theta2, phi2), ...]
    labels : list of strings, state labels
    """
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')

    # Draw sphere wireframe
    u = np.linspace(0, 2*np.pi, 50)
    v = np.linspace(0, np.pi, 50)
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones(np.size(u)), np.cos(v))
    ax.plot_wireframe(x, y, z, alpha=0.1, color='gray')

    # Draw axes
    ax.quiver(0, 0, 0, 1.3, 0, 0, color='r', arrow_length_ratio=0.1, alpha=0.5)
    ax.quiver(0, 0, 0, 0, 1.3, 0, color='g', arrow_length_ratio=0.1, alpha=0.5)
    ax.quiver(0, 0, 0, 0, 0, 1.3, color='b', arrow_length_ratio=0.1, alpha=0.5)
    ax.text(1.4, 0, 0, 'X', fontsize=12)
    ax.text(0, 1.4, 0, 'Y', fontsize=12)
    ax.text(0, 0, 1.4, 'Z (|0⟩)', fontsize=12)
    ax.text(0, 0, -1.4, '|1⟩', fontsize=12)

    # Plot states
    colors = plt.cm.viridis(np.linspace(0, 1, len(states)))
    for i, (theta, phi) in enumerate(states):
        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)
        ax.quiver(0, 0, 0, x, y, z, color=colors[i], arrow_length_ratio=0.1)
        ax.scatter([x], [y], [z], s=100, color=colors[i])
        if labels:
            ax.text(x*1.1, y*1.1, z*1.1, labels[i], fontsize=10)

    ax.set_xlim([-1.5, 1.5])
    ax.set_ylim([-1.5, 1.5])
    ax.set_zlim([-1.5, 1.5])
    ax.set_box_aspect([1, 1, 1])
    ax.set_title('Bloch Sphere Representation')
    plt.show()

# Example states
states = [
    (0, 0),              # |0⟩
    (np.pi, 0),          # |1⟩
    (np.pi/2, 0),        # |+⟩ = (|0⟩+|1⟩)/√2
    (np.pi/2, np.pi),    # |-⟩ = (|0⟩-|1⟩)/√2
    (np.pi/2, np.pi/2),  # |+i⟩ = (|0⟩+i|1⟩)/√2
]
labels = ['|0⟩', '|1⟩', '|+⟩', '|-⟩', '|+i⟩']

bloch_sphere_visualization(states, labels)

# Calculate state vector from Bloch angles
def bloch_to_state(theta, phi):
    alpha = np.cos(theta/2)
    beta = np.exp(1j*phi) * np.sin(theta/2)
    return np.array([alpha, beta])

print("State vectors:")
for (theta, phi), label in zip(states, labels):
    state = bloch_to_state(theta, phi)
    print(f"{label}: α={state[0]:.3f}, β={state[1]:.3f}")

Single Qubit Gates

Quantum gates are represented by unitary operators. The key single-qubit gates are:

import numpy as np

class SpinQubit:
    """Single spin qubit with gate operations"""

    def __init__(self, state=None):
        """Initialize qubit, default to |0⟩"""
        if state is None:
            self.state = np.array([1, 0], dtype=complex)
        else:
            self.state = np.array(state, dtype=complex)
            self.state /= np.linalg.norm(self.state)

    # Pauli matrices
    X = np.array([[0, 1], [1, 0]], dtype=complex)
    Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
    Z = np.array([[1, 0], [0, -1]], dtype=complex)

    # Hadamard gate
    H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)

    def apply_gate(self, gate):
        """Apply a single qubit gate"""
        self.state = gate @ self.state
        return self

    def rotation(self, axis, angle):
        """Apply rotation around axis by angle"""
        if axis == 'x':
            gate = np.cos(angle/2)*np.eye(2) - 1j*np.sin(angle/2)*self.X
        elif axis == 'y':
            gate = np.cos(angle/2)*np.eye(2) - 1j*np.sin(angle/2)*self.Y
        elif axis == 'z':
            gate = np.cos(angle/2)*np.eye(2) - 1j*np.sin(angle/2)*self.Z
        self.state = gate @ self.state
        return self

    def measure(self):
        """Measure in computational basis, returns 0 or 1"""
        prob_0 = np.abs(self.state[0])**2
        result = 0 if np.random.random() < prob_0 else 1
        # Collapse state
        self.state = np.array([1, 0] if result == 0 else [0, 1], dtype=complex)
        return result

    def expectation(self, observable):
        """Calculate expectation value of observable"""
        return np.real(np.conj(self.state) @ observable @ self.state)

    def __repr__(self):
        return f"|ψ⟩ = {self.state[0]:.3f}|0⟩ + {self.state[1]:.3f}|1⟩"

# Demonstration
qubit = SpinQubit()
print(f"Initial state: {qubit}")

qubit.apply_gate(SpinQubit.H)
print(f"After Hadamard: {qubit}")

qubit.rotation('z', np.pi/4)
print(f"After Z rotation (π/4): {qubit}")

# Expectation values
print(f"\n⟨X⟩ = {qubit.expectation(SpinQubit.X):.3f}")
print(f"⟨Y⟩ = {qubit.expectation(SpinQubit.Y):.3f}")
print(f"⟨Z⟩ = {qubit.expectation(SpinQubit.Z):.3f}")

2.2 Quantum Coherence and Relaxation

Quantum coherence is essential for quantum information processing. Understanding and mitigating decoherence mechanisms is crucial for practical spin qubits.

T1 and T2 Relaxation Times

Bloch Equations

The dynamics of a spin qubit in a magnetic field are described by the Bloch equations:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def bloch_equations(M, t, gamma, B, T1, T2, M0):
    """
    Bloch equations for spin relaxation

    Parameters:
    -----------
    M : array, magnetization vector [Mx, My, Mz]
    t : float, time
    gamma : float, gyromagnetic ratio
    B : array, magnetic field [Bx, By, Bz]
    T1, T2 : float, relaxation times
    M0 : float, equilibrium magnetization
    """
    Mx, My, Mz = M
    Bx, By, Bz = B

    # Precession term (M × B)
    dMx = gamma * (My * Bz - Mz * By) - Mx / T2
    dMy = gamma * (Mz * Bx - Mx * Bz) - My / T2
    dMz = gamma * (Mx * By - My * Bx) - (Mz - M0) / T1

    return [dMx, dMy, dMz]

# Parameters
gamma = 1.76e11  # Electron gyromagnetic ratio (rad/s/T)
B0 = 0.1  # Static field (T)
B = [0, 0, B0]
T1 = 1e-6  # 1 μs
T2 = 0.5e-6  # 0.5 μs
M0 = 1.0

# Initial state: spin in x-y plane
M_init = [1, 0, 0]

# Time evolution
t = np.linspace(0, 3e-6, 1000)  # 3 μs
solution = odeint(bloch_equations, M_init, t, args=(gamma, B, T1, T2, M0))

Mx, My, Mz = solution.T

# Plot
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Time evolution
axes[0, 0].plot(t*1e6, Mx, 'b-', label='Mx')
axes[0, 0].plot(t*1e6, My, 'g-', label='My')
axes[0, 0].plot(t*1e6, Mz, 'r-', label='Mz')
axes[0, 0].set_xlabel('Time (μs)')
axes[0, 0].set_ylabel('Magnetization')
axes[0, 0].legend()
axes[0, 0].set_title('Bloch Equation Dynamics')
axes[0, 0].grid(True, alpha=0.3)

# Transverse magnetization decay
Mxy = np.sqrt(Mx**2 + My**2)
axes[0, 1].plot(t*1e6, Mxy, 'b-', label='|M⊥|')
axes[0, 1].plot(t*1e6, np.exp(-t/T2), 'r--', label=f'exp(-t/T2), T2={T2*1e6:.1f}μs')
axes[0, 1].set_xlabel('Time (μs)')
axes[0, 1].set_ylabel('Transverse Magnetization')
axes[0, 1].legend()
axes[0, 1].set_title('T2 Decay')
axes[0, 1].grid(True, alpha=0.3)

# Longitudinal recovery
axes[1, 0].plot(t*1e6, Mz, 'r-', label='Mz')
axes[1, 0].plot(t*1e6, M0*(1-np.exp(-t/T1)), 'b--', label=f'M0(1-exp(-t/T1)), T1={T1*1e6:.1f}μs')
axes[1, 0].set_xlabel('Time (μs)')
axes[1, 0].set_ylabel('Longitudinal Magnetization')
axes[1, 0].legend()
axes[1, 0].set_title('T1 Recovery')
axes[1, 0].grid(True, alpha=0.3)

# 3D trajectory
ax3d = fig.add_subplot(2, 2, 4, projection='3d')
ax3d.plot(Mx, My, Mz, 'b-', linewidth=0.5)
ax3d.scatter([Mx[0]], [My[0]], [Mz[0]], c='g', s=100, label='Start')
ax3d.scatter([Mx[-1]], [My[-1]], [Mz[-1]], c='r', s=100, label='End')
ax3d.set_xlabel('Mx')
ax3d.set_ylabel('My')
ax3d.set_zlabel('Mz')
ax3d.set_title('Bloch Vector Trajectory')
ax3d.legend()

plt.tight_layout()
plt.show()

Decoherence Mechanisms

2.3 Spin Qubit Platforms

Semiconductor Quantum Dots

Electrons confined in semiconductor quantum dots form highly controllable spin qubits. Gate electrodes define the potential well and control the qubit.

NV Centers in Diamond

Nitrogen-vacancy (NV) centers in diamond are optically addressable spin defects with exceptional coherence properties at room temperature.

import numpy as np
import matplotlib.pyplot as plt

def nv_center_hamiltonian(B, theta, D=2.87e9, E=0):
    """
    Calculate NV center Hamiltonian in ground state triplet

    Parameters:
    -----------
    B : float, magnetic field strength (T)
    theta : float, angle between B and NV axis (rad)
    D : float, zero-field splitting (Hz)
    E : float, strain-induced splitting (Hz)

    Returns:
    --------
    H : array, 3x3 Hamiltonian matrix
    energies : array, eigenvalues
    """
    gamma_e = 28.024e9  # Electron gyromagnetic ratio (Hz/T)

    # Spin-1 matrices
    Sz = np.diag([1, 0, -1])
    Sx = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / np.sqrt(2)
    Sy = np.array([[0, -1j, 0], [1j, 0, -1j], [0, 1j, 0]]) / np.sqrt(2)

    # Magnetic field components
    Bz = B * np.cos(theta)
    Bx = B * np.sin(theta)

    # Hamiltonian: H = D*Sz^2 + E*(Sx^2 - Sy^2) + gamma*B·S
    H = (D * Sz @ Sz +
         E * (Sx @ Sx - Sy @ Sy) +
         gamma_e * (Bz * Sz + Bx * Sx))

    energies, states = np.linalg.eigh(H)
    return H, energies, states

# Calculate energy levels vs magnetic field
B_range = np.linspace(0, 0.1, 200)  # 0 to 100 mT
energies_vs_B = []

for B in B_range:
    _, E, _ = nv_center_hamiltonian(B, theta=0)  # B parallel to NV axis
    energies_vs_B.append(E / 1e9)  # Convert to GHz

energies_vs_B = np.array(energies_vs_B)

# Plot
plt.figure(figsize=(10, 6))
labels = ['|ms=-1⟩', '|ms=0⟩', '|ms=+1⟩']
for i in range(3):
    plt.plot(B_range*1000, energies_vs_B[:, i], label=labels[i])

plt.xlabel('Magnetic Field (mT)')
plt.ylabel('Energy (GHz)')
plt.title('NV Center Energy Levels vs Magnetic Field')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Zero-field splitting
print(f"Zero-field splitting D = 2.87 GHz")
print(f"ms=0 to ms=±1 transition at B=0: {2.87:.2f} GHz")

# ODMR frequencies at 10 mT
B = 0.01  # 10 mT
_, E, _ = nv_center_hamiltonian(B, theta=0)
f1 = (E[2] - E[1]) / 1e9
f2 = (E[1] - E[0]) / 1e9
print(f"\nODMR frequencies at B=10mT:")
print(f"  f+ = {f1:.3f} GHz")
print(f"  f- = {f2:.3f} GHz")

2.4 Entanglement and Two-Qubit Operations

Quantum entanglement is essential for universal quantum computation. Two-qubit gates create entangled states that cannot be described as products of individual qubit states.

Bell States

The four maximally entangled two-qubit states are:

import numpy as np

class TwoQubitSystem:
    """Two-qubit system for entanglement studies"""

    def __init__(self):
        """Initialize to |00⟩ state"""
        self.state = np.zeros(4, dtype=complex)
        self.state[0] = 1  # |00⟩

    # Single qubit gates (tensor products)
    I = np.eye(2, dtype=complex)
    X = np.array([[0, 1], [1, 0]], dtype=complex)
    Z = np.array([[1, 0], [0, -1]], dtype=complex)
    H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)

    # Two-qubit gates
    CNOT = np.array([[1,0,0,0], [0,1,0,0], [0,0,0,1], [0,0,1,0]], dtype=complex)
    SWAP = np.array([[1,0,0,0], [0,0,1,0], [0,1,0,0], [0,0,0,1]], dtype=complex)

    def apply_single(self, gate, qubit):
        """Apply single qubit gate to specified qubit (0 or 1)"""
        if qubit == 0:
            full_gate = np.kron(gate, self.I)
        else:
            full_gate = np.kron(self.I, gate)
        self.state = full_gate @ self.state
        return self

    def apply_two_qubit(self, gate):
        """Apply two-qubit gate"""
        self.state = gate @ self.state
        return self

    def create_bell_state(self, which='phi+'):
        """Create Bell state"""
        self.state = np.zeros(4, dtype=complex)
        self.state[0] = 1  # Reset to |00⟩

        # Apply H to first qubit, then CNOT
        self.apply_single(self.H, 0)
        self.apply_two_qubit(self.CNOT)

        # Modify for different Bell states
        if which == 'phi-':
            self.apply_single(self.Z, 0)
        elif which == 'psi+':
            self.apply_single(self.X, 1)
        elif which == 'psi-':
            self.apply_single(self.Z, 0)
            self.apply_single(self.X, 1)

        return self

    def concurrence(self):
        """Calculate concurrence (entanglement measure)"""
        # Reshape to 2x2 matrix
        rho = np.outer(self.state, np.conj(self.state))

        # Spin-flip transformation
        sigma_y = np.array([[0, -1j], [1j, 0]])
        sigma_yy = np.kron(sigma_y, sigma_y)
        rho_tilde = sigma_yy @ np.conj(rho) @ sigma_yy

        # Calculate concurrence
        product = rho @ rho_tilde
        eigenvalues = np.sqrt(np.maximum(0, np.linalg.eigvalsh(product)))
        eigenvalues = np.sort(eigenvalues)[::-1]

        return max(0, eigenvalues[0] - eigenvalues[1] - eigenvalues[2] - eigenvalues[3])

    def chsh_test(self, n_measurements=10000):
        """
        CHSH Bell inequality test
        Classical bound: |S| <= 2
        Quantum maximum: |S| <= 2√2 ≈ 2.83
        """
        # Measurement angles
        a1, a2 = 0, np.pi/2  # Alice's settings
        b1, b2 = np.pi/4, -np.pi/4  # Bob's settings

        def measure_correlation(theta_a, theta_b):
            """Measure correlation E(a,b)"""
            # Measurement operators
            A = np.cos(theta_a)*self.Z + np.sin(theta_a)*self.X
            B = np.cos(theta_b)*self.Z + np.sin(theta_b)*self.X
            AB = np.kron(A, B)

            # Expectation value
            return np.real(np.conj(self.state) @ AB @ self.state)

        # Calculate CHSH parameter S
        E11 = measure_correlation(a1, b1)
        E12 = measure_correlation(a1, b2)
        E21 = measure_correlation(a2, b1)
        E22 = measure_correlation(a2, b2)

        S = E11 - E12 + E21 + E22
        return S

    def __repr__(self):
        labels = ['|00⟩', '|01⟩', '|10⟩', '|11⟩']
        terms = []
        for i, amp in enumerate(self.state):
            if np.abs(amp) > 1e-10:
                terms.append(f"{amp:.3f}{labels[i]}")
        return "|ψ⟩ = " + " + ".join(terms)

# Demonstration
system = TwoQubitSystem()

# Create and analyze Bell states
bell_states = ['phi+', 'phi-', 'psi+', 'psi-']
print("Bell States Analysis:")
print("-" * 50)

for bell in bell_states:
    system.create_bell_state(bell)
    C = system.concurrence()
    S = system.chsh_test()
    print(f"|{bell}⟩: {system}")
    print(f"   Concurrence: {C:.3f}")
    print(f"   CHSH S value: {S:.3f} (quantum violation: |S| > 2)")
    print()

# Compare with product state
system.state = np.array([1, 0, 0, 0], dtype=complex)  # |00⟩
print("Product state |00⟩:")
print(f"   Concurrence: {system.concurrence():.3f}")
print(f"   CHSH S value: {system.chsh_test():.3f}")

2.5 Spin-Photon Coupling

Coupling spins to photons enables long-range quantum information transfer and provides a key interface between spin qubits and quantum communication networks.

Jaynes-Cummings Model

The fundamental model for a single spin coupled to a single cavity mode is:

where $\omega_c$ is the cavity frequency, $\omega_s$ is the spin transition frequency, $g$ is the coupling strength, and $\hat{a}^\dagger, \hat{a}$ are photon creation/annihilation operators.

import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt

def jaynes_cummings_hamiltonian(omega_c, omega_s, g, n_photons=5):
    """
    Build Jaynes-Cummings Hamiltonian

    Parameters:
    -----------
    omega_c : float, cavity frequency
    omega_s : float, spin frequency
    g : float, coupling strength
    n_photons : int, maximum photon number

    Returns:
    --------
    H : array, Hamiltonian matrix
    """
    dim = 2 * (n_photons + 1)

    # Photon operators
    a = np.zeros((n_photons+1, n_photons+1))
    for n in range(n_photons):
        a[n, n+1] = np.sqrt(n+1)
    a_dag = a.T

    # Spin operators
    sigma_z = np.array([[1, 0], [0, -1]])
    sigma_plus = np.array([[0, 1], [0, 0]])
    sigma_minus = np.array([[0, 0], [1, 0]])

    # Tensor products
    I_photon = np.eye(n_photons + 1)
    I_spin = np.eye(2)

    H_cavity = omega_c * np.kron(a_dag @ a, I_spin)
    H_spin = 0.5 * omega_s * np.kron(I_photon, sigma_z)
    H_int = g * (np.kron(a_dag, sigma_minus) + np.kron(a, sigma_plus))

    H = H_cavity + H_spin + H_int
    return H

def vacuum_rabi_oscillation(omega_c, omega_s, g, t_max, n_points=200):
    """Simulate vacuum Rabi oscillations"""
    n_photons = 5
    H = jaynes_cummings_hamiltonian(omega_c, omega_s, g, n_photons)

    # Initial state: spin excited, cavity empty (|e,0⟩)
    dim = 2 * (n_photons + 1)
    psi0 = np.zeros(dim, dtype=complex)
    psi0[1] = 1  # |0⟩_photon ⊗ |e⟩_spin

    # Time evolution
    t = np.linspace(0, t_max, n_points)
    excitation_prob = []
    photon_number = []

    for ti in t:
        U = expm(-1j * H * ti)
        psi_t = U @ psi0

        # Spin excitation probability
        p_excited = 0
        for n in range(n_photons + 1):
            idx = 2*n + 1  # |n⟩_photon ⊗ |e⟩_spin
            p_excited += np.abs(psi_t[idx])**2
        excitation_prob.append(p_excited)

        # Average photon number
        n_avg = 0
        for n in range(n_photons + 1):
            p_n = np.abs(psi_t[2*n])**2 + np.abs(psi_t[2*n+1])**2
            n_avg += n * p_n
        photon_number.append(n_avg)

    return t, np.array(excitation_prob), np.array(photon_number)

# Parameters
omega_c = 1.0  # Cavity frequency (normalized)
g = 0.1  # Coupling strength

# Resonant case
omega_s_res = omega_c
t, P_e_res, n_ph_res = vacuum_rabi_oscillation(omega_c, omega_s_res, g, 100)

# Detuned case
omega_s_det = omega_c + 0.2
t, P_e_det, n_ph_det = vacuum_rabi_oscillation(omega_c, omega_s_det, g, 100)

# Plot
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].plot(t*g, P_e_res, 'b-', label='Resonant')
axes[0].plot(t*g, P_e_det, 'r--', label='Detuned (δ=0.2)')
axes[0].set_xlabel('Time (gt)')
axes[0].set_ylabel('Spin Excitation Probability')
axes[0].set_title('Vacuum Rabi Oscillations')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t*g, n_ph_res, 'b-', label='Resonant')
axes[1].plot(t*g, n_ph_det, 'r--', label='Detuned')
axes[1].set_xlabel('Time (gt)')
axes[1].set_ylabel('Average Photon Number')
axes[1].set_title('Cavity Photon Number')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Rabi frequency
Omega_R = 2 * g
print(f"Vacuum Rabi frequency: Ω_R = 2g = {Omega_R:.2f}")
print(f"Rabi period: T_R = π/g = {np.pi/g:.2f}")

Strong Coupling Regime

Strong coupling is achieved when $g > \kappa, \gamma$ (cavity decay and spin relaxation rates). This enables:

• Reversible energy exchange between spin and photon
• Formation of polariton states
• Quantum state transfer between spin and photon

2.6 Topological Quantum States

Topological protection of quantum states offers robustness against local perturbations. Majorana fermions in topological superconductors are promising candidates for fault-tolerant qubits.

Kitaev Chain Model

The Kitaev chain is a 1D model of spinless fermions with p-wave superconducting pairing:

where $\mu$ is the chemical potential, $t$ is the hopping amplitude, and $\Delta$ is the pairing amplitude.

import numpy as np
import matplotlib.pyplot as plt

def kitaev_chain(N, mu, t, Delta):
    """
    Build Kitaev chain Hamiltonian in BdG form

    Parameters:
    -----------
    N : int, number of sites
    mu : float, chemical potential
    t : float, hopping amplitude
    Delta : float, pairing amplitude

    Returns:
    --------
    H : array, BdG Hamiltonian (2N x 2N)
    """
    H = np.zeros((2*N, 2*N), dtype=complex)

    for i in range(N):
        # On-site terms
        H[i, i] = -mu/2
        H[N+i, N+i] = mu/2

        # Hopping terms
        if i < N-1:
            H[i, i+1] = -t
            H[i+1, i] = -t
            H[N+i, N+i+1] = t
            H[N+i+1, N+i] = t

            # Pairing terms
            H[i, N+i+1] = Delta
            H[N+i+1, i] = Delta
            H[i+1, N+i] = -Delta
            H[N+i, i+1] = -Delta

    return H

# Parameters for topological phase
N = 50
t = 1.0
Delta = 1.0

# Scan chemical potential
mu_range = np.linspace(-4, 4, 100)
energies = []
gaps = []

for mu in mu_range:
    H = kitaev_chain(N, mu, t, Delta)
    E = np.linalg.eigvalsh(H)
    energies.append(E)
    # Energy gap (smallest positive energy)
    gaps.append(np.min(np.abs(E)))

energies = np.array(energies)

# Plot energy spectrum
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
for i in range(2*N):
    plt.plot(mu_range, energies[:, i], 'b-', alpha=0.3, linewidth=0.5)
plt.axhline(0, color='r', linestyle='--', alpha=0.5)
plt.axvline(-2*t, color='g', linestyle='--', alpha=0.5, label='Phase boundary')
plt.axvline(2*t, color='g', linestyle='--', alpha=0.5)
plt.xlabel('Chemical Potential μ/t')
plt.ylabel('Energy E/t')
plt.title('Kitaev Chain Energy Spectrum')
plt.legend()
plt.ylim(-3, 3)

plt.subplot(1, 2, 2)
plt.semilogy(mu_range, gaps)
plt.axvline(-2*t, color='g', linestyle='--', alpha=0.5, label='Phase boundary')
plt.axvline(2*t, color='g', linestyle='--', alpha=0.5)
plt.xlabel('Chemical Potential μ/t')
plt.ylabel('Energy Gap')
plt.title('Gap Closing at Topological Phase Transition')
plt.legend()

plt.tight_layout()
plt.show()

# Analyze zero modes in topological phase
mu_topo = 0  # In topological phase
H_topo = kitaev_chain(N, mu_topo, t, Delta)
E_topo, V_topo = np.linalg.eigh(H_topo)

# Find zero modes
zero_mode_idx = np.where(np.abs(E_topo) < 0.01)[0]
print(f"Number of near-zero modes: {len(zero_mode_idx)}")

if len(zero_mode_idx) >= 2:
    # Plot wavefunction localization
    plt.figure(figsize=(10, 4))
    for idx in zero_mode_idx[:2]:
        psi = V_topo[:, idx]
        prob = np.abs(psi[:N])**2 + np.abs(psi[N:])**2
        plt.plot(range(N), prob, label=f'E = {E_topo[idx]:.4f}')

    plt.xlabel('Site')
    plt.ylabel('Probability Density')
    plt.title('Majorana Zero Mode Localization')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

print("\nTopological phase: |μ| < 2t")
print("Trivial phase: |μ| > 2t")