A Problem using Gutzwiller-RVB theory for high-temperature superconductivity

[dieon]

I'm currently working on reproducing a graph from the Gutzwiller-RVB theory of high-temperature superconductivity using Python. The graph plots χ_BBxy versus Δ, where Δ ranges from 16 to 26. I've implemented the theory's iterative parameter update scheme but encountered discrepancies in the plotted results compared to the expected graph.

Here's a simplified version of my Python code:

[CODE lang="python" title="code"]import numpy as np
import matplotlib.pyplot as plt

# Constants and parameters (initial guesses)
U = 20.0
t = 1.0
a = 1.0
N_1 = 20
N_2 = 20
N = N_1 * N_2
max_iterations = 500
tolerance = 1e-5

# Initial guesses for the parameters
m_s = 0.8
delta = 0.8
chi_AB_up = 0.8
chi_AB_down = 0.8
chi_BB_up = 0.8
chi_BB_down = 0.8
chi_BBxy_up = 0.8
chi_BBxy_down = 0.8
Delta_AB = 0.8

# Generate wave vectors
def generate_wave_vectors(a, N_1, N_2):
b1 = (2 * np.pi / a) * np.array([1, 0])
b2 = (2 * np.pi / a) * np.array([0, 1])
K = []
for n_1 in range(-N_1 // 2, N_1 // 2):
for n_2 in range(-N_2 // 2, N_2 // 2):
k = (n_1 / N_1) * b1 + (n_2 / N_2) * b2
K.append(k)
return np.array(K)

K = generate_wave_vectors(a, N_1, N_2)

def gamma_k(k):
return 2 * (np.cos(k[0]) + np.cos(k[1]))

def gamma_sc(k):
return np.cos(k[0]) - np.cos(k[1])

def gamma_k_prime(k):
return 2 * (np.cos(2*k[0]) + np.cos(2*k[1])) + 4 * (np.cos(k[0] + k[1]) + np.cos(k[0] - k[1]))

def h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U):
g1 = 1
gs = 4 / ((1 + delta)**2 - m_s**2)
h1_up_val = (U - Delta) / 2 * (1 + delta - m_s) / 2 \
- t**2 / Delta * (4 * (1 - 2 * delta) + g1 * (2 * chi_BB_up + 4 * chi_BBxy_up) + g1 * (1 - delta + m_s) / 2 * gamma_k_prime(k)) \
- 2 * t**2 * gs * m_s / (U + Delta) + 2 * t**2 / (U + Delta) * (1 - delta)
return h1_up_val

def h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U):
g1 = 1
gs = 4 / ((1 + delta)**2 - m_s**2)
h1_down_val = (U - Delta) / 2 * (1 + delta + m_s) / 2 \
- t**2 / Delta * (4 * (1 - 2 * delta) + g1 * (2 * chi_BB_down + 4 * chi_BBxy_down) + g1 * (1 - delta - m_s) / 2 * gamma_k_prime(k)) \
+ 2 * t**2 * gs * m_s / (U + Delta) + 2 * t**2 / (U + Delta) * (1 - delta)
return h1_down_val

def h2_up(k, delta, m_s, chi_AB_up, t, Delta, U):
g1 = 1
g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
gs = 4 / ((1 + delta)**2 - m_s**2)
h2_up_val = -t * g1 - t**2 / Delta * (-2 * chi_AB_up + 6 * g2 * chi_AB_up) \
- t**2 / (U + Delta) * (gs * (0.5 * chi_AB_up + chi_AB_up) + 0.5 * chi_AB_up)
return h2_up_val * gamma_k(k)

def h2_down(k, delta, m_s, chi_AB_down, t, Delta, U):
g1 = 1
g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
gs = 4 / ((1 + delta)**2 - m_s**2)
h2_down_val = -t * g1 - t**2 / Delta * (-2 * chi_AB_down + 6 * g2 * chi_AB_down) \
- t**2 / (U + Delta) * (gs * (0.5 * chi_AB_down + chi_AB_down) + 0.5 * chi_AB_down)
return h2_down_val * gamma_k(k)

def h3(k, delta, m_s, Delta_AB, t, Delta, U):
g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
gt_up = 2 * delta / (1 + delta + m_s)
gt_down = 2 * delta / (1 + delta - m_s)
gs = 4 / ((1 + delta)**2 - m_s**2)
h3_val = 4 * t**2 / Delta * (1 - g2) + 4 * t**2 / (U + Delta) * (3 * gs / 4 - 1 / 4) - 2 * t**2 / Delta * (gt_up + gt_down)
return h3_val * Delta_AB / 2 * gamma_sc(k)

def omega_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
alpha_up, beta_up, _, _ = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
return h1_up_k * (alpha_up**2 - beta_up**2) - 2 * h2_up_k * alpha_up * beta_up

def omega_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U):
_, _, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)
h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
return h1_down_k * (alpha_down**2 - beta_down**2) - 2 * h2_down_k * alpha_down * beta_down

def nu(k, delta, m_s, Delta_AB, t, Delta, U):
alpha_up, beta_up, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
return -h3(k, delta, m_s, Delta_AB, t, Delta, U) * (alpha_up * beta_down + alpha_down * beta_up)

def alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
E_up_k_val = E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
E_down_k_val = E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)

alpha_up = np.sqrt(0.5 * (1 - h1_up_k / E_up_k_val))
beta_up = np.sqrt(0.5 * (1 + h1_up_k / E_up_k_val))
alpha_down = np.sqrt(0.5 * (1 - h1_down_k / E_down_k_val))
beta_down = np.sqrt(0.5 * (1 + h1_down_k / E_down_k_val))
return alpha_up, beta_up, alpha_down, beta_down

def E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
return np.sqrt(h1_up_k**2 + h2_up_k**2)

def E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U):
h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
return np.sqrt(h1_down_k**2 + h2_down_k**2)

def u_v_squared(k, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BBxy_up, chi_BB_down, chi_BBxy_down, t, Delta, U):
omega_up_val = omega_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
omega_down_val = omega_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)
nu_k = nu(k, delta, m_s, Delta_AB, t, Delta, U)

denominator = np.sqrt((omega_up_val + omega_down_val)**2 + 4 * nu_k**2)
u1_squared = 0.5 * (1 + (omega_up_val + omega_down_val) / denominator)
u2_squared = 0.5 * (1 - (omega_up_val + omega_down_val) / denominator)
v1_squared = u2_squared
v2_squared = u1_squared
return u1_squared, u2_squared, v1_squared, v2_squared


# Function to update parameters
def update_parameters(k_points, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB, t, Delta):
sum_m_s = 0
sum_delta = 0
sum_chi_AB_up = 0
sum_chi_AB_down = 0
sum_chi_BB_up = 0
sum_chi_BB_down = 0
sum_chi_BBxy_up = 0
sum_chi_BBxy_down = 0
sum_Delta_AB = 0

for k in k_points:
N_1 = 20
N_2 = 20
N = N_1 * N_2

alpha_up, beta_up, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
E_up_k_val = E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
E_down_k_val = E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)

u1_squared, u2_squared, v1_squared, v2_squared = u_v_squared(k, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BBxy_up, chi_BB_down, chi_BBxy_down, t, Delta, U)

u1 = np.sqrt(u1_squared)
u2 = np.sqrt(u2_squared)
v1 = np.sqrt(v1_squared)
v2 = np.sqrt(v2_squared)

# Update m_s and delta
sum_m_s += ((alpha_up**2 - alpha_down**2) * v1_squared + (beta_up**2 - beta_down**2) * v2_squared)/N
sum_delta += (0.5 * (alpha_up**2 * (v1_squared - v2_squared) + beta_up**2 * (v1_squared - v2_squared) + alpha_down**2 * (v1_squared - v2_squared) + beta_down**2 * (v1_squared - v2_squared)))/N

# Update chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down
sum_chi_AB_up += 0.25 * alpha_up * beta_up * (v2_squared - v1_squared)/N
sum_chi_AB_down += 0.25 * alpha_down * beta_down * (v2_squared - v1_squared)/N
sum_chi_BB_up += (np.cos(2*k[0]) + np.cos(2*k[1])) * (alpha_up**2 * v2_squared + beta_up**2 * v1_squared)/N
sum_chi_BB_down += (np.cos(2*k[0]) + np.cos(2*k[1])) * (alpha_down**2 * v2_squared + beta_down**2 * v1_squared)/N
sum_chi_BBxy_up += 2 * np.cos(k[0]) * np.cos(k[1]) * (alpha_up**2 * v2_squared + beta_up**2 * v1_squared)/N
sum_chi_BBxy_down += 2 * np.cos(k[0]) * np.cos(k[1]) * (alpha_down**2 * v2_squared + beta_down**2 * v1_squared)/N
sum_Delta_AB += (alpha_down * beta_up * u2 * v2 - alpha_up * beta_down * u1 * v1) * gamma_sc(k)

new_m_s = sum_m_s
new_delta = sum_delta
new_chi_AB_up = sum_chi_AB_up
new_chi_AB_down = sum_chi_AB_down
new_chi_BB_up = sum_chi_BB_up
new_chi_BB_down = sum_chi_BB_down
new_chi_BBxy_up = sum_chi_BBxy_up
new_chi_BBxy_down = sum_chi_BBxy_down
new_Delta_AB = sum_Delta_AB

return new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB

# Initializing the plot data lists
Delta_values = np.linspace(16, 26, 50)
chi_BBxy_up_values = []
chi_BBxy_down_values = []

# Loop over Delta values
for Delta in Delta_values:
m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB = 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8
for iteration in range(max_iterations):
new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB = update_parameters(
K, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB, t, Delta
)

if np.abs(new_Delta_AB - Delta_AB) < tolerance:
break

m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB = new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB

chi_BBxy_up_values.append(chi_BBxy_up)
chi_BBxy_down_values.append(chi_BBxy_down)

# Plotting the results
plt.plot(Delta_values, chi_BBxy_down_values, label='χ_BBxy_down', marker='o')
plt.plot(Delta_values, chi_BBxy_up_values, label='χ_BBxy_up', marker='+')
plt.xlabel('Δ')
plt.ylabel('χ_BBxy')
plt.title('Plot of χ_BBxy vs Δ')
plt.grid(True)
plt.show()[/CODE]

Key Issues:

1. Graph Accuracy: The plotted graph of χ_BBxy against Δ does not closely match the expected results from the literature. Even small changes in initial parameters or the size of the wave vector grid (N_1 and N_2) significantly alter the output.
2. Parameter Sensitivity: I observe that varying initial guesses (m_s, delta, etc.) and the grid size (N_1, N_2) affect the final plotted values of χ_BBxy. Is this expected behavior, or could there be an error in how I'm updating parameters or generating wave vectors?

Request for Assistance:

• Clarification: Is it typical for small changes in initial parameters or the grid size to affect the output significantly in numerical simulations of the Gutzwiller-RVB theory?
• Code Review: Could you please review the provided Python code snippet to identify potential errors or improvements that might lead to more accurate results?

(I've also attached a file with the formulas I've used)

Any guidance or suggestions would be greatly appreciated!