Я получаю отрицательные элементы матрицы при эволюции по уравнению ГКСЛ. ⇐ Python

[Anonymous]

В этом коде я получаю отрицательные записи в матрице плотности. Кто-нибудь знает какой-нибудь метод, который мы можем использовать? Модуль собственных значений не меняет того факта, что я получаю отрицательный квантовый дискорд, а реализация других подходов, позволяющих убедиться, что матрица всегда является положительной полуопределенной, не предотвращает превращение квантового дискорда в отрицательный.
Вот код. Любые комментарии приветствуются!
import numpy as np
import matplotlib.pyplot as plt

# Define the mixed states
p = 0.6
q = 0.4

# Define the pure states
psi3 = (1/np.sqrt(2)) * np.array([1, 0, 0, 1]) # |psi3> = (|00> + |01>)/sqrt(2)
psi4 = np.array([0, 0, 1, 0]) # |psi4> = |10>
psi5 = np.array([0, 0, 0, 1]) # |psi5> = |11>
psi6 = (np.sqrt(3)/2) * np.array([1, 0, 0, 0]) + (1/2) * np.array([0, 0, 1, 0]) # |psi6> = (sqrt(3)/2)|00> + (1/2)|10>

# Mixed states
rho1 = p * np.outer(psi3, np.conj(psi3)) + (1 - p) * np.outer(psi4, np.conj(psi4))
rho2 = q * np.outer(psi5, np.conj(psi5)) + (1 - q) * np.outer(psi6, np.conj(psi6))

# Define the Hamiltonian and the coupling constant (2 qubits)
H = np.kron(np.array([[0, 1], [1, 0]]), np.array([[0, 1], [1, 0]])) # Pauli-X for two qubits

# Define collapse operator C
C = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2)) # Raise first qubit, identity on second qubit

gamma = 1 # Coupling constant

# Function to compute the entropy of a density matrix
def entropy(rho):
eigenvalues = np.linalg.eigvals(rho)
eigenvalues = np.maximum(np.real(eigenvalues), 1e-10) # Avoid log(0) and negative values
return -np.sum(eigenvalues * np.log(eigenvalues))

# Regularization function to ensure positive semi-definiteness of the density matrix
def regularize_matrix(rho, epsilon=1e-10):
eigvals, eigvecs = np.linalg.eigh(rho)
eigvals = np.maximum(eigvals, epsilon) # Set any negative eigenvalues to epsilon
return np.dot(eigvecs, np.dot(np.diag(eigvals), eigvecs.T))

# Function to compute the partial trace
def partial_trace(rho, dims, subsystem):
"""
Compute the partial trace of a density matrix `rho`.

Args:
rho (ndarray): Density matrix (shape (d1*d2, d1*d2)).
dims (tuple): Dimensions of the subsystems (d1, d2).
subsystem (int): Subsystem to trace out (0 for the first subsystem, 1 for the second).

Returns:
ndarray: Reduced density matrix.
"""
d1, d2 = dims
rho = rho.reshape((d1, d2, d1, d2)) # Reshape into a 4D tensor
if subsystem == 0:
return np.trace(rho, axis1=0, axis2=2) # Trace out the first subsystem
elif subsystem == 1:
return np.trace(rho, axis1=1, axis2=3) # Trace out the second subsystem
else:
raise ValueError("Subsystem must be 0 or 1.")

# Lindblad evolution step (discrete time step)
def lindblad_step(rho, dt, gamma, H, C):
commutator = -1j * (np.dot(H, rho) - np.dot(rho, H))
dissipative = gamma * (np.dot(C, np.dot(rho, C.conj().T)) - 0.5 * (np.dot(np.dot(C.conj().T, C), rho) + np.dot(rho, np.dot(C.conj().T, C))))
rho_new = rho + dt * (commutator + dissipative)
rho_new = regularize_matrix(rho_new)
return rho_new

# Calculate the quantum discord
def quantum_discord(rho):
"""
Compute an approximation of the quantum discord for a 2-qubit system.

Args:
rho (ndarray): Full density matrix of the 2-qubit system.

Returns:
float: Quantum discord approximation.
"""
dims = (2, 2) # Dimensions of the two qubits

# Partial traces to get reduced density matrices
rho_A = partial_trace(rho, dims, 1) # Trace out subsystem B
rho_B = partial_trace(rho, dims, 0) # Trace out subsystem A

# Compute entropies
S_A = entropy(rho_A) # Entropy of subsystem A
S_B = entropy(rho_B) # Entropy of subsystem B
S_AB = entropy(rho) # Entropy of the total system

# Total correlation: mutual information
I_A_B = S_A + S_B - S_AB

# Approximation of classical correlation (simplified but improved)
classical_corr = max(S_A, S_B) # Classical correlation: use max(S_A, S_B) as a more reasonable approximation

return I_A_B - classical_corr

# Time evolution
time_steps = 1000
dt = 0.01
times = np.linspace(0, time_steps * dt, time_steps)

# Store discord values for both mixed states
discord_values_rho1 = []
discord_values_rho2 = []

# Evolve the systems and compute discord at each time step
for t in times:
rho1 = lindblad_step(rho1, dt, gamma, H, C)
discord_values_rho1.append(quantum_discord(rho1))

rho2 = lindblad_step(rho2, dt, gamma, H, C)
discord_values_rho2.append(quantum_discord(rho2))

# Plot the discord over time for the mixed states
plt.plot(times, discord_values_rho1, label=r'$\rho_1$', color='tab:red')
plt.plot(times, discord_values_rho2, label=r'$\rho_2$', color='tab:blue')
plt.xlabel('Time')
plt.ylabel('Quantum Discord')
plt.title(r'Quantum Discord of the Mixed States $\rho_1$ and $\rho_2$ Over Time')
plt.grid(True)
plt.legend()
plt.show()


Подробнее здесь: https://stackoverflow.com/questions/792 ... l-equation