Как смоделировать на питоне сверхпроводя плотность состояний суперпроводника BCS с YSR (Yu-Shiba-Rusinov) с 1D-моделью с ⇐ Python

[Anonymous]

Я пытаюсь смоделировать следующую задачу:
Вычислить плотность состояний численно у суперпроводника BCS с магнитными примеси, используя модель плотного связывания 1D -цепи с периодически граничными условиями.
На изображении есть изображение гамильтониана.

Гамильтониан и матрица для диагонализации < /p>
Чтобы решить, я хочу Чтобы написать матрицу, диагонализируйте ее и используйте решения собственных значений и собственных векторов для вычисления плотности состояний < /p>
Плотность состояний формулы < /p>
Это Код у меня есть: < /p>

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt

#parameters
N = 2000  # matrix size
delta = 1  # gap
t = 5.0  #  hopping
U = 20.0
JS = 17
eta = 0.01 # broadening

# initialize the matrix with zeros
H = np.zeros((N, N), dtype=complex)

# Fill the matrix with the gap and hoping values
for i in range(0, N - 4, 4):
H[i, i + 3] = -delta
H[i + 1, i + 2] = delta
H[i + 2, i + 1] = np.conj(delta)
H[i + 3, i] = -np.conj(delta)
H[i, i + 4] = -t
H[i + 1, i + 5] = -t
H[i + 2, i + 6] = t
H[i + 3, i + 7] = t
H[i + 4, i] = -t
H[i + 5, i + 1] = -t
H[i + 6, i + 2] = t
H[i + 7, i + 3] = t

H[N-4, 0] = -t
H[N-3, 1] = -t
H[N-2, 2] = t
H[N-1, 3] = t
H[0, N-4] = -t
H[1, N-3] = -t
H[2, N-2] = t
H[3, N-1] = t

H[0, 0] = U + JS
H[1, 1] = U - JS
H[2, 2] = -(U + JS)
H[3, 3] = -(U - JS)

# diagonalize the matrix
eigenvalues, eigenvectors = np.linalg.eigh(H)

# Compute the LDOS in one site
def compute_ldos(E_vals, E, psi, eta, site):
""" Calcula la LDOS en un sitio específico. """
a_i = psi[4 * site, :]
b_i = psi[4 * site + 1, :]
c_i = psi[4 * site + 2, :]
d_i = psi[4 * site + 3, :]
LDOS = (np.abs(a_i)**2 + np.abs(b_i)**2) * (eta / np.pi) / ((E - E_vals)**2 + eta**2) + \
(np.abs(c_i)**2 + np.abs(d_i)**2) * (eta / np.pi) / ((E + E_vals)**2 + eta**2)
return np.sum(LDOS)

#
E_range = np.linspace(-3, 3 , 2000)

# Calculate the density of states in site 250
site_index =250
LDOS_site = np.array([compute_ldos(eigenvalues, E, eigenvectors, eta, site_index) for E in E_range])

plt.figure(figsize=(8, 6))
plt.plot(E_range, LDOS_site, label=f'Densidad de estados en sitio {site_index}')
plt.xlabel('E')
plt.ylabel('LDOS')

plt.legend()
plt.grid()
plt.xlim(-3,3)

plt.show()
< /code>
Это то, что я получаю, когда пытаюсь вычислить плотность состояний на сайте 250: < /p>

Очевидно, что мой код не работает, как и ожидалось, я ожидаю, что он получит что -то Рядом с BCS LDOS. Я не понимаю, откуда поступают колебания и почему далеко от разрыва он не нормализует до одного. Кроме того, я видел это с n = 8000 

, для вычисления требуется слишком много времени, поэтому мне также интересно, есть ли более эффективный способ.

Подробнее здесь: https://stackoverflow.com/questions/794 ... cs-superco