Моделирование динамики длинного джозефсоновского перехода

Contents

Моделирование динамики длинного джозефсоновского перехода#

1. Теория и система уравнений#

1.1 Эффект Джозефсона и джозефсоновский переход#

Связь двух сверхпроводящих слоев посредством тонкого слоя несверхпроводящего барьера образует структуру, называемую джозефсоновским переходом (в честь английского ученого Брайана Джозефсона). При пропускании электрического тока через джозефсоновский переход (ДП), в зависимости от величины тока, наблюдаются стационарный и нестационарный эффекты Джозефсона.

Стационарный эффект Джозефсона. При пропускании тока \(I\) через ДП в отсутствии напряжения (\(V=0\)) течет сверхпроводящий ток \(I_{s}\), не превышающий определённого предела, называемого критическим током \(I_{c}\) (\(I<I_{c}\)). Величина тока пропорциональна синусу разности фаз \(\varphi\) параметров порядка сверхпроводящих слоев ДП: \( \begin{eqnarray} \displaystyle I_{s}=I_{c}\sin\varphi. \tag{1} \end{eqnarray} \) Нестационарный эффект Джозефсона. При увеличении значение тока выше критического \((I>I_{c})\) в ДП возникает переменное напряжение, которое пропорционально производной разности фаз по времени \( \begin{eqnarray} \displaystyle V=\frac{\hbar}{2e}\frac{d\varphi}{dt}. \tag{2} \end{eqnarray} \)

1.2 Геометрические эффекты в джозефсоновских переходах: короткие и длинные переходы#

Характерное расстояние, на которое магнитное поле может проникнуть в ДП, называется джозефсоновской глубиной проникновения \(\lambda_J\). Она определяется выражением \(\lambda_{J}=\sqrt{\hbar S /(2e\mu_{0}DI_{c})}\), где \(S\) - площадь сверхпроводящего слоя, \(e\) - заряд электрона, \(\mu_{0}\) - магнитная постоянная, \(D=2\lambda_{L}+d\) - эффективная магнитная толщина, \(\lambda_{L}\) - глубина проникновения магнитного поля в сверхпроводниках (Лондоновская глубина) и \(d\) - толщина несверхпроводящего барьера.

Если длина ДП \(L\) намного меньше джозефсоновской глубины проникновения \(\lambda_J\) (\(L<<\lambda_J\)), то такой переход называется коротким. В случае короткого ДП разность фаз полагается однородной вдоль ДП.

Если длина ДП \(L\) сравнима или превышает джозефсоновскую глубину проникновения \(\lambda_J\) (\(L\gtrapprox\lambda_J\)), то такой переход называется длинным. В случае длинного ДП пространственное распределение фазы и магнитного поля внутри перехода становится критически важным для физики процесса.

1.3 Математическая модель#

Динамика короткого джозефсоновского перехода описывается в рамках модели RCSJ [1], где переход моделируется как параллельное соединения конденсатора с емкостью \(C\), резистора с сопротивлением \(R\) и сверхпроводника.

Тогда полный ток \(I\) через переход состоит из суммы тока смещения \(I_{d}=C\frac{\partial V}{\partial t}\), квазичастичного (нормального) тока \(I_{qp}=\frac{V}{R}\) и сверхпроводящего тока \(I_{s}=I_{c}\sin\varphi\)

\( \begin{eqnarray} \tag{3} I = C\frac{\partial V}{\partial t} + \frac{V}{R} + I_{c}\sin\varphi. \end{eqnarray} \)

Напряжения \(V\) определяется соoтношением Джозефсона для нестационарного эффекта (2).

В случае длинного ДП, необходимо учитывать пространственную зависимость разности фаз. Для этого нужно модифицировать уравнения RCSJ-модели, учитывая возникающий поверхностный ток \(I_{sf}\) вдоль длины ДП. Данный ток пропорционален производной второго порядка от разности фаз по координате и определяется выражением

\( \begin{eqnarray} \tag{4} I_{sf}=I_{c}\lambda_{J}^{2}\frac{\partial^{2}\varphi}{\partial x^{2}}. \end{eqnarray} \)

Учитывая (2), (3) и (4) можно записать уравнение для описания длинного ДП

\( \begin{eqnarray} \tag{5} \frac{\hbar C}{2e}\frac{\partial ^{2}\varphi}{\partial t^{2}}-I_{c}\lambda_{J}^{2}\frac{\partial ^{2}\varphi}{\partial x^{2}}+I_{c}\sin\varphi+\frac{\hbar}{2eR}\frac{\partial \varphi}{\partial t}=I. \end{eqnarray} \)

Начальные условия для данного уравнения задаются нулевой разностью фаз и нулевым напряжением

\(\begin{eqnarray} \tag{6} \varphi(0,x)=0, \hspace{0.5 cm} V(0,x)=0, \end{eqnarray}\)

а граничные условия задаются через внешнее магнитное поле \(H_{ext}\) на краях длинного ДП, т.е.

\(\begin{eqnarray} \tag{7} \frac{\partial \varphi}{\partial x}|_{x=0}=\frac{\partial \varphi}{\partial x}|_{x=L}=\frac{2e\mu_{0}D}{\hbar} H_{ext} \end{eqnarray}\)

[1] W. C. Stewart, “Current-voltage characteristics of Josephson junctions,” J. Appl. Phys. 12, 277-280 (1968)

1.3.1 Переход к безразмерным величинам#

Для численного решения уравнении (5) необходимо записать систему уравнений в нормированных величинах

\( \begin{eqnarray} \tag{8} \frac{\partial ^{2}\varphi}{\partial t^{2}}-\frac{\partial ^{2}\varphi}{\partial x^{2}}+\sin\varphi+\beta\frac{\partial \varphi}{\partial t}=I, \end{eqnarray} \)

где \(\beta=\sqrt{\hbar/(2e I_{c}R^{2}C)}\) - параметр диссипации, время нормировано на \(\omega_{p}^{-1}\), \(\omega_{p}=\sqrt{2e I_{c}/(\hbar C)}\), координата \(x\) на \(\lambda_{J}\), внешний ток \(I\) на \(I_{c}\).

Уравнение (8) известно как уравнение синуса-Гордона с возмущением. Граничные условия в нормированных величинах записываются в виде:

\( \begin{eqnarray} \tag{9} \frac{\partial \varphi}{\partial x}|_{x=0}=\frac{\partial \varphi}{\partial x}|_{x=L}=H_{ext}. \end{eqnarray} \)

2. Задачи#

  • Расчет динамики длинного джозефсоновского перехода;

  • Расчет вольт-амперной характеристики;

  • Расчет зависимости критического тока от внешнего магнитного поля.

3. Расчет динамики длинного джозефсоновского перехода#

Расчет динамики основан на численном решении уравнения в частных производных (8) с учетом начальных и граничных условий.

3.1 Вычислительная схема#

Вычислительная схема основана на построении одномерной равномерной сетки по пространственной координате с аппроксимацией производных по пространственной координате и применением метода Рунге-Кутты четвертого порядка.

3.1.1 Построение одномерной равномерной сетки по координате#

Задается равномерная сетка по координате от \(x \in [0,L]\) c шагом \(h\). Номера узлов с координатами \(x\) и \(x\pm h\) обозначаются как \(i\) и \(i\pm 1\), соответственно: \(\omega_{h}= \{ x_{i}=ih | i=0,1,...,n_{x}, h = L / (n_{x}-1)\}\).

3.1.2 Аппроксимация производных по координате#

Аппроксимация производных по координате выполнена с помощью метода конечных разностей: для граничных точек использованы односторонние разности (вперёд и назад), для внутренних (серединных) точек — центральные разности.

Построим сетку в области изменения переменных \(x\) с шагом \(h\). Как указывалось выше, введем следующие обозначения для значений функции \(\varphi(x,t)\) в узлах сетки:

\( \begin{eqnarray} \tag{10} \varphi(x,t)=\varphi_{i}(t), \hspace{0.5 cm} \varphi(x\pm h,t)=\varphi_{i\pm 1}(t). \end{eqnarray} \)

Крайним точкам с координатами \(x=0\) и \(x=L\) соответствуют узлы с номерами \(i=0\) и \(i=n_{x}\), а внутренним точкам в интервале \(0 < x < L\) - \(i\).

Аппроксимируя производную второго порядка находим

\( \begin{eqnarray} \tag{11} \frac{\partial^{2} \varphi_{0}(t)}{\partial x^{2}} = \frac{-7\varphi_{i}(t)+ 8\varphi_{i+1}(t) - \varphi_{i+2}(t)-6hH_{ext}}{2h^{2}} + O(h^{2}), \end{eqnarray} \)

для левого края при \(i=0\),

\( \begin{eqnarray} \tag{12} \frac{\partial^{2} \varphi_{i}(t)}{\partial x^{2}} = \frac{\varphi_{i+1}(t)-2\varphi_{i}(t) + \varphi_{i-1}(t)}{h^{2}} + O(h^{2}) \end{eqnarray} \)

для внутренних точек \(0<i<n_{x}\),

\( \begin{eqnarray} \tag{12} \frac{\partial^{2} \varphi_{0}(t)}{\partial x^{2}} = \frac{-7\varphi_{i}(t)+8\varphi_{i-1}(t) - \varphi_{i-2}(t)+6hH_{ext}}{2h^{2}} + O(h^{2}) \end{eqnarray} \)

для правого края при \(i=n_{x}\). При аппроксимации мы учли граничные условия - при \(i=0\) и \(i=n_{x}\) \(\frac{\partial \varphi_{i}}{\partial x} = H_{ext}\).

После аппроксимации производных второго порядка по координате уравнение (8) преобразуется в систему \(n_x\) дифференциальных уравнений второго порядка относительно времени. Чтобы применить метод Рунге-Кутты четвертого порядка, перепишем его в виде системы \(2n_{x}\) уравнений

\( \begin{eqnarray} \tag{13} \frac{\partial V_{i}}{\partial t} &=& I + \frac{\partial ^{2}\varphi_{i}}{\partial x^{2}} + \sin\varphi_{i} + \beta\frac{\partial \varphi_{i}}{\partial t}\\ \frac{\partial \varphi_{i}}{\partial t} &=& V_{i} \end{eqnarray} \)

3.1.3 Метод Рунге-Кутты четвертого порядка#

Кратко опишем метод Рунге-Кутты четвертого порядка.

Рассмотрим задачу Коши для обыкновенного дифференциального уравнения первого порядка \( \begin{eqnarray} \tag{14} \frac{dy}{dt}=f(t,y); \hspace{0.5cm} y(t_{0})=y_{0}. \end{eqnarray} \)

Тогда приближенное значение в последующих точках вычисляется по итерационной формуле \( \begin{eqnarray} \tag{15} y_{n+1}=y_{n}+\frac{\tau}{6}(k_{n,1}+2k_{n,2}+2k_{n,3}+k_{n,4}). \end{eqnarray} \)

Вычисление нового значения осуществляется в четыре этапа: \( \begin{eqnarray} k_{n,1}&=&f(t_{n},y_{n}),\\ k_{n,2}&=&f(t_{n}+\tau/2, y_{n}+(\tau/2) k_{n,1}),\\ \tag{16} k_{n,3}&=&f(t_{n}+\tau/2, y_{n}+(\tau/2) k_{n,2}),\\ k_{n,4}&=&f(t_{n}+\tau, y_{n}+\tau k_{n,3}). \end{eqnarray} \) где \(\tau\) - величина шага сетки по \(t\).

3.2 Алгоритм численного решения (метод Рунге - Кутты четвертого порядка)#

  • Инициализация:

    • задаются значения параметров модели;

    • задается начальное условие;

    • задаются максимальное время \(T_\max\) и шаг по времени \(\Delta t\) (равномерная сетка по времени).

  • Определение правых частей уравнений в виде функции \(F(t,x)\).

  • Вычисление коэффициентов Рунге-Кутты \(k_{n,1}\), \(k_{n,2}\), \(k_{n,3}\), \(k_{n,4}\) для каждой точки сетки с многократным обращением к \(F(t,x)\) и нахождение значений функции во всех точках временной сетки.

3.3 Программная реализация численного решения#

Подключение библиотек#

# Подключение библиотек
import numpy as np
import numba as nb

import random

import matplotlib.pyplot as plt

import time

from numba import njit, prange
import numba
from numba import types
from numba.typed import Dict

import seaborn as sns
sns.set()
sns.set(style="whitegrid")

%matplotlib inline

Задание функции правых частей#

# Функция для правой стороны системы уравнений
@njit
def equations(t, y, Inoise, p):
    '''Функция для правой стороны системы уравнений
         t - значение времени
         y - массив начальных условий, который имеет 2*nx ячеек
             первый nx ячейка соответсвуют массиву ph
             следующие nx соответствуют массиву V
    Inoise - массив содержащий добавляемый шум к току
             распределенный по координате
         p - массив параметров модели
    '''
    I = p[0]
    H_ext = p[1]
    beta = p[2]
    dx = p[3]
    L = p[4]

    nx = int(L/dx) + 1

    ph = y[0:nx]
    V = y[nx:2*nx]

    dy = np.zeros(2*nx, dtype=np.float64)
    d2ph = np.empty(nx, dtype=np.float64)

    # Левая граница
    d2ph[0] = (-7*ph[0] + 8*ph[1] - ph[2] - 6*dx*H_ext) / (2*dx*dx)

    # Внутренние точки
    d2ph[1:-1] = (ph[0:-2] - 2*ph[1:-1] + ph[2:]) / (dx*dx)

    # Правая граница
    d2ph[-1] = (-7*ph[-1] + 8*ph[-2] - ph[-3] + 6*dx*H_ext) / (2*dx*dx)

    # dph
    dy[0:nx] = V

    # dV
    dy[nx:2*nx] = I + Inoise[:nx] + d2ph - np.sin(ph) - beta * V

    return dy

Определение функции для метода Рунге-Кутты четвертого порядка#

# Метод Рунге-Кутты четвертого порядка
@njit
def runge_kutta_4(t_end, deltat, y0, Inoise, p) -> (np.ndarray, np.ndarray):
    N = int(np.ceil(t_end / deltat)) + 1
    nx = int(p[4]/p[3])+1
    t_vals = np.zeros(N)
    y_vals = np.empty((N, 2*nx))
    y_vals[0] = y0

    for i in range(1, N):
        t = t_vals[i - 1]
        y = y_vals[i - 1]

        k1 = equations(t, y, Inoise, p)
        k2 = equations(t + deltat/2., y + deltat/2. * k1, Inoise, p)
        k3 = equations(t + deltat/2., y + deltat/2. * k2, Inoise, p)
        k4 = equations(t + deltat, y + deltat * k3, Inoise, p)

        y_vals[i] = y + (deltat / 6.) * (k1 + 2*k2 + 2*k3 + k4)
        t_vals[i] = t + deltat

    return t_vals, y_vals

Задание параметров модели и начальных условий#

beta = 0.2   # Параметр диссипации
H_ext = 0.0  # Внешнее магнитное поле
L = 3.0      # Длина джозефсоновского перехода

dx = 0.1
dt = dx/5.0

Tmax = 500

I = 0.45

noisemax = 10**-8

p = np.array([I, H_ext, beta, dx, L])

nx = int(L/dx) + 1
y0 = np.zeros(2*nx, dtype=np.float64)

Inoise = np.zeros(nx, dtype=np.float64)
Inoise = noisemax * np.random.random(nx)

3.4 Стационарный эффект Джозефсона#

y0[0:nx] = 0.0
y0[nx:2*nx] = 0.0

t_vals, f_spatiotemporal_vals = runge_kutta_4(Tmax, dt, y0, Inoise, p)

Построение временной зависимости при фиксированном значении пространственной координаты \(x_{fix}\)#

x_fix = 2
nx_fix = int(x_fix/dx)

ph_vals_stationar = f_spatiotemporal_vals[:, nx_fix]
V_vals_stationar = f_spatiotemporal_vals[:, nx+nx_fix]

plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(t_vals, ph_vals_stationar, label=f'$\\varphi(t,x={x_fix})$')
plt.xlabel('t')
plt.ylabel('$\\varphi$')
plt.legend()
plt.text(-0.12, 0.95, 'a)', transform=plt.gca().transAxes, fontsize=12,
         fontweight='normal', ha='right')

plt.subplot(1, 2, 2)
plt.plot(t_vals, V_vals_stationar, label=f'V(t,x={x_fix})')
plt.xlabel('t')
plt.ylabel('V')
plt.text(-0.12, 0.95, 'b)', transform=plt.gca().transAxes, fontsize=12,
         fontweight='normal', ha='right')


plt.legend()

plt.tight_layout()
plt.show()
_images/f01a4e836d935365821c575ccaf34d74faa9f93fbdb900c6e2f94c8396ffb8d7.png

Рис. 1. Временная зависимость при фиксированной координате \(x\): a) разность фаз \(\varphi\); b) напряжение \(V\).

Построения зависимости от координаты при фиксированном значении времени \(t_{fix}\)#

t_fix = 100
nt_fix = int(t_fix/dt)
x_vals = np.linspace(0, L, nx)
ph_tfix = f_spatiotemporal_vals[nt_fix, 0:nx]
V_tfix = f_spatiotemporal_vals[nt_fix, nx:2*nx]

plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(x_vals, ph_tfix, label=f'$\\varphi(t={t_fix},x)$')
plt.xlabel('x')
plt.ylabel('$\\varphi$')
plt.legend()
plt.text(-0.15, 0.95, 'a)', transform=plt.gca().transAxes,
         fontsize=12, fontweight='normal', ha='right')

plt.subplot(1, 2, 2)
plt.plot(x_vals, V_tfix, label=f'V(t={t_fix},x)')
plt.xlabel('x')
plt.ylabel('V')
plt.text(-0.12, 0.95, 'b)', transform=plt.gca().transAxes,
         fontsize=12, fontweight='normal', ha='right')

plt.legend()
plt.tight_layout()
plt.show()
_images/571e70d4592e2511b51db5d5cee4c6680b5a49c5d20509553f9da7a12c9271e5.png

Рис. 2. Пространственная зависимость при фиксированном времени \(t\): a) разности фаз \(\varphi\); b) напряжения \(V\).

3.5 Нестационарный эффект Джозефсона#

y0[0:nx] = 0.0
y0[nx:2*nx] = 2.2

t_vals, f_spatiotemporal_vals = runge_kutta_4(Tmax, dt, y0, Inoise, p)

Построение временной зависимости при фиксированном значении координаты \(x_{fix}\)#

x_fix = 2
nx_fix = int(x_fix/dx)

ph_vals_nonstationar = f_spatiotemporal_vals[:, nx_fix]
V_vals_nonstationar = f_spatiotemporal_vals[:, nx+nx_fix]

plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(t_vals, ph_vals_nonstationar, label=f'$\\varphi(t,x={x_fix})$')
plt.xlabel('t')
plt.ylabel('$\\varphi$')
plt.legend()
plt.text(-0.12, 0.95, 'a)', transform=plt.gca().transAxes,
         fontsize=12, fontweight='normal', ha='right')

plt.subplot(1, 2, 2)
plt.plot(t_vals, V_vals_nonstationar, label=f'V(t,x={x_fix})')
plt.xlabel('t')
plt.ylabel('V')
plt.text(-0.12, 0.95, 'b)', transform=plt.gca().transAxes,
         fontsize=12, fontweight='normal', ha='right')

plt.legend()
plt.tight_layout()
plt.show()
_images/a5be14196f765522ce3e0153a85903bd167dd2c7a62576ac641196190bd3365f.png

Рис. 3. Временная зависимость при фиксированной координате \(x\): a) разности фаз \(\varphi\); b) напряжения \(V\).

Построения зависимости от координаты при фиксированном значении времени \(t_{fix}\)#

t_fix = 250
nt_fix = int(t_fix/dt)
x_vals = np.linspace(0, L, nx)
ph_tfix = f_spatiotemporal_vals[nt_fix, 0:nx]
V_tfix = f_spatiotemporal_vals[nt_fix, nx:2*nx]

plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(x_vals, ph_tfix, label=f'$\\varphi(t={t_fix},x)$')
plt.xlabel('x')
plt.ylabel('$\\varphi$')
plt.legend()
plt.text(-0.15, 0.95, 'a)', transform=plt.gca().transAxes, fontsize=12,
         fontweight='normal', ha='right')

plt.subplot(1, 2, 2)
plt.plot(x_vals, V_tfix, label=f'V(t={t_fix},x)')
plt.xlabel('x')
plt.ylabel('V')
plt.text(-0.12, 0.95, 'b)', transform=plt.gca().transAxes, fontsize=12,
         fontweight='normal', ha='right')

plt.legend()
plt.tight_layout()
plt.show()
_images/0e01f43027f628f15dbcac264fa588dc733d05895bba0cdcbf210b720b45a773.png

Рис. 4. Пространственная зависимость при фиксированном времени \(t\): a) разности фаз \(\varphi\); b) напряжения \(V\).

Построения пространственно-временной зависимости напряжения V(x,t)#

from mpl_toolkits.mplot3d import Axes3D

# 1. Задаем координатную сетку
x_vals = np.linspace(0, L, nx)
X, T = np.meshgrid(x_vals, t_vals)

# 2. Определяем координату Z (поверхность)
Z = f_spatiotemporal_vals[:, nx:2*nx]

# Построение 3D графика
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

step_t = 5
step_x = 1
T_flat = T[::step_t, ::step_x].flatten()
X_flat = X[::step_t, ::step_x].flatten()
Z_flat = Z[::step_t, ::step_x].flatten()

sc = ax.scatter3D(T_flat, X_flat, Z_flat, c=Z_flat,
                  cmap='viridis', s=5, alpha=0.7)

fig.colorbar(sc, ax=ax, label='V(t,x)')

ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('V(t,x)')
ax.set_title('Пространственно-временная зависимость напряжения')

plt.tight_layout()
plt.show()
_images/61bbc3a808c046f825f8cd042e60e58dd9bb68dd054f82a6976874719a975f99.png

Рис. 5. Пространственно-временная зависимость напряжения \(V(t,x)\).

3.6 Расчет пространственно-временной зависимости магнитного поля#

Магнитное поле определяется как производная от разности фаз по координате.

Создание функции для вычисления производной разности фаз по координате для всех значений времени и координаты \(H(x,t)\)#

@njit
def derivative_by_x(f, dx):
    nt, nx = f.shape
    dfdx = np.zeros((nt, nx))

    for i in range(nt):
        # производная левого края
        dfdx[i, 0] = (-3*f[i, 0] + 4*f[i, 1] - f[i, 2]) / (2.0*dx)

        # производная внутренних точек
        for j in range(1, nx-1):
            dfdx[i, j] = (f[i, j+1] - f[i, j-1]) / 2*dx

        # производная правого края
        dfdx[i, -1] = (3*f[i, -1] - 4*f[i, -2] + f[i, -3]) / (2*dx)
    return dfdx

# Извлекается значения ph (x,t) из решения системы уравнений
# для дальнейшего расчета магнитного поля
ph_vals = f_spatiotemporal_vals[:, 0:nx]
# Расчитается производная по координате, т.е. значения магнитного поля
# для всех значений времени
dphx = derivative_by_x(ph_vals, dx)

Построение пространственно-временной зависимости магнитного поля \(H(x,t)\)#

# Построение 3D графика
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

step_t = 5
step_x = 1
T_flat = T[::step_t, ::step_x].flatten()
X_flat = X[::step_t, ::step_x].flatten()
H_flat = dphx[::step_t, ::step_x].flatten()

sc = ax.scatter3D(T_flat, X_flat, H_flat, c=H_flat,
                  cmap='viridis', s=5, alpha=0.7)

fig.colorbar(sc, ax=ax, label='H(t,x)')

ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('H(t,x)')
ax.set_title('Пространственно-временная зависимость магнитного поля')

plt.tight_layout()
plt.show()
_images/7bb6a914938fb2c0f391763c85162bd63fbe1dc0e485ee43c43e61d071621def.png

Рис. 6. Пространственно-временная зависимость магнитного поля \(H(t,x)\).

4. Расчет вольт-амперной характеристики#

4.1 Алгоритм расчета вольт-амперной характеристики#

  1. Задается значения параметров модели, параметры численного счета, начальные условия.

  2. Численно решается система уравнений (5), для фиксированного значения тока \(I\), и расчитается пространственно-временная зависимость разности фаз \(\varphi(t,x)\) и напряжения \(V(t,x)\).

  3. Полученная зависимость \(V(t,x)\) усредняется по координате

\( \begin{eqnarray} \displaystyle V_{x}(t)=\frac{1}{L}\int\limits_{0}^{L}V(t,x)dx, \end{eqnarray} \)

полученная зависимость \(V_{x}(t)\) усредняется по времени

\( \begin{eqnarray} \displaystyle <V>=\frac{1}{T_{\max}-T_{\min}}\int\limits_{T_{\min}}^{T_{\max}}V_{x}(t)dt \end{eqnarray} \)

и в результате вычисляется среднее напряжение для заданного значения тока, т.е. определяется одна точка на ВАХ.

  1. Далее меняется значение тока на \(\Delta I\) и повторяется п. 2, используя при этом \(\varphi(T_{\max},x)\) и \(V(T_{\max},x)\) из посчитанного значения тока в качестве начального условия, для полученного \(V(t,x)\) выполняется п. 3, и вычисляется новая точка ВАХ.

  2. Таким образом увеличивая и уменьшая значение тока можно расчитать ВАХ длинного ДП.

4.2 Реализация усреднения напряжения#

Функция для интегрирования по пространственной координате#

@njit
def integrate_rows_trapz(Y, dx):
    """Интегрирование методом трапеций для каждой строки двумерного массива"""
    m, n = Y.shape
    results = np.zeros(m, dtype=Y.dtype)
    for i in range(m):
        row_sum = 0.0
        for j in range(n-1):
            row_sum += (Y[i, j] + Y[i, j+1]) / 2 * dx
        results[i] = row_sum
    return results

Функция для усреднения по времени#

@njit
def averageV(Tmin, Tmax, dt, L, dx, V_vals):
    ntmin = int(Tmin/dt)
    ntmax = int(Tmax/dt) + 1
    X = np.arange(0, L + dx/2, dx)  # массив координат (включая L)
    Y = V_vals  # двумерный массив, т.е. значения V(x,t)

    # Вычисления среднего значения по координате
    V_av_x = integrate_rows_trapz(Y, dx)/L

    # Вычисления среднего значения по времени
    V_av = dt * np.sum(V_av_x[ntmin:ntmax])/(Tmax-Tmin)
    return V_av

Проверка корректности рассчета усреднения#

Tmin = 200.0
Tmax = 500.0

p[0] = 0.45
y0[0:nx] = 0
y0[nx:2*nx] = 2.2

t_vals, f_spatiotemporal_vals = runge_kutta_4(Tmax, dt, y0, Inoise, p)
ph_vals = f_spatiotemporal_vals[:, 0:nx]
V_vals = f_spatiotemporal_vals[:, nx:2*nx]

averageV(Tmin, Tmax, dt, L, dx, V_vals)

1.726447529178154

4.3 Вычисление одной точки ВАХ#

@njit
def cvc_point(Tmin, Tmax, dt, y0, Inoise, p):
    ''' Функция для вычисления одной точки ВАХ'''

    nx = int(len(y0)/2)
    dx = p[3]
    L = p[4]

    y_initial = np.zeros_like(y0)

    # Вызов метода Рунге-Кутты 4-го порядка
    t_vals, f_spatiotemporal_vals = runge_kutta_4(Tmax, dt, y0, Inoise, p)

    # Извлечения решения (временная зависимость)
    ph_vals = f_spatiotemporal_vals[:, 0:nx]
    V_vals = f_spatiotemporal_vals[:, nx:2*nx]

    # Усреднение
    Vav = averageV(Tmin, Tmax, dt, L, dx, V_vals)

    # Обновление начальных условий
    y_initial[0:nx] = ph_vals[-1, :]
    y_initial[nx:2*nx] = V_vals[-1, :] + Inoise[0:nx]

    return Vav, y_initial

p[0] = 0.45
y0[0:nx] = 0
y0[nx:2*nx] = 2.2

Vav, y0 = cvc_point(Tmin, Tmax, dt, y0, Inoise, p)
print(Vav)

1.726447529178154

4.4 Вычисление двухпетлевой ВАХ#

@njit
def comput_cvc(Tmin, Tmax, dt, y0, p, Imax1, Imin2, Imax2,
               Ibreak, delta_I, noisemax):

    nx = int(len(y0)/2)
    Inoise = np.zeros(nx+2)

    I = 0.0
    b = 1

    loop = 1
    ready_for_stop = False
    big_loop_done = False

    Vplot = []
    Iplot = []

    while True:

        for x in range(nx):
            Inoise[x] = noisemax*random.random()

        p[0] = I
        V_av, y0 = cvc_point(Tmin, Tmax, dt, y0, Inoise, p)

        Vplot.append(V_av)
        Iplot.append(I)

        if (I < (Imin2+delta_I/2.2) and big_loop_done == False and b < 0):
            big_loop_done = True
            loop = 2
            b = 1
		
		#Условие для изменения направлении тока
        if (I > Imax1-delta_I/2.0):
            b=-1
        
        if (I > (Imax2-delta_I/2.0) and big_loop_done == True and b > 0):
            ready_for_stop = True
            b = -1

        I = I + b*delta_I
        #Условие для остановки вычисления
        if ((I < Ibreak-delta_I/2.0) and (ready_for_stop == True)):
            break
            
    return np.column_stack((np.array(Iplot), np.array(Vplot)))

# Задается параметры модели (model parameters values)
L=3.0		
beta=0.2			
H_ext=0			

dx=0.1
		
Imax1 = 1.1
Imin2 = 0.0 
Imax2 = 0.0
Ibreak = 0.0

H_ext = 0.0


delta_I = 0.005

Tmax = 500.0  # Максимальное значение времени (max time)
Tmin = 200.0
dt = dx/5.0 #Шаг по времени (step size)

noisemax = 10**-8

p = np.array([I, H_ext, beta, dx, L])

y0[0:nx] = 0.0
y0[nx:2*nx] = 0.0

Imax1 = 1.1
Imin2 = 0
Imax2 = 0
Ibreak = 0

start = time.time()
res_cvc = comput_cvc(Tmin, Tmax, dt, y0, p, Imax1, Imin2, Imax2,
                     Ibreak, delta_I, noisemax)
end = time.time()

print(end-start)

27.543137073516846

Imax1 = 1.1
Imin2 = 0.4
Imax2 = 1.05
Ibreak = 1.0

start = time.time()
res_cvc_f1 = comput_cvc(Tmin, Tmax, dt, y0, p, Imax1, Imin2, Imax2,
                        Ibreak, delta_I, noisemax)
end = time.time()

print(end-start)

31.610062837600708

Imax1 = 1.1
Imin2 = 0.83
Imax2 = 1.05
Ibreak = 1.0

start = time.time()
res_cvc_f2 = comput_cvc(Tmin, Tmax, dt, y0, p, Imax1, Imin2, Imax2,
                        Ibreak, delta_I, noisemax)
end = time.time()

print(end-start)

19.821210861206055

fig = plt.figure(figsize=(12, 8))
plt.plot(res_cvc[:, 0], res_cvc[:, 1], label='Однопетлевая ВАХ', linewidth=3.0)
plt.plot(res_cvc_f1[:, 0], res_cvc_f1[:, 1],
         label='Двухпетлевая ВАХ с однофлюксонной ветвью', linewidth=3.0)
plt.plot(res_cvc_f2[:, 0], res_cvc_f2[:, 1],
         label='Двухпетлевая ВАХ с двухфлюксонной ветвью', linewidth=3.0)
plt.xlabel('I', size=12)
plt.ylabel('V', size=12)
plt.legend(loc='upper left')
plt.show()
_images/73ca45eebecb3a6440ef9e5593aaf12bb79c4ad1267348768241a0784fdf53a3.png

Рис. 7. Трехпетлевая вольт-амперная характеристика длинного джозефсоновского перехода.

5. Расчет зависимости критического тока от внешнего магнитного поля#

5.1 Алгоритм расчета зависимости критического тока от внешнего магнитного поля#

  1. Задаются значения параметров.

  2. Вычисляется ВАХ при фиксированном значении внешнего магнитного поля. В процессе вычисления сохраняется значения тока, при котором выполняется условиe \(V>\epsilon\), что и является критическим током.

  3. Увеличивается значение внешнего магнитного поля \(\Delta H\) повторяется п. 2.

  4. Пункт 3 выполняется до достижения значения \(H_{\max}\).

5.2 Программная реализация расчета зависимости критического тока от внешнего магнитного поля#

Поскольку расчет зависимости критического тока от магнитного поля при разных значениях поля не зависят друг от друга, расчет можно реализовать в параллельном режиме.#

Подключение библиотек#

from numba import config, njit, threading_layer
from numba import set_num_threads, get_num_threads
# set the threading layer before any parallel target compilation
config.THREADING_LAYER = 'threadsafe'

config.NUMBA_DEFAULT_NUM_THREADS

256

Функция для определения значении критического тока для фиксированного значения внешнего магнитного поля#

import numba as nb
import time

@njit
def find_Ic(k, Hmin, dH, Tmin, Tmax, dt, p, delta_I, noisemax, epsilon):

    H_ext = Hmin + dH * k
    p[1] = H_ext

    nx = int(p[4]/p[3]) + 1
    Inoise = np.zeros(nx)

    y0 = np.zeros(2 * nx, dtype=np.float64)

    I = 0.0
    while True:
        p[0] = I

        for x in range(nx):
            Inoise[x] = noisemax*random.random()

        # Вычисление одной точки ВАХ
        V_av, y = cvc_point(Tmin, Tmax, dt, y0, Inoise, p)
        y0 = y

        if V_av > epsilon:
            Ic = I
            break

        I += delta_I

        if I > 5.0:
            return I
    return Ic

Функция для расчета значении критического тока в параллельном режиме#

@njit(parallel=True)
def comput_Ic_on_H_parallel(H_points, Hmin, dH, Tmin, Tmax, dt, p,
                            delta_I, noisemax, epsilon):
    Ic_array = np.empty(H_points)

    for k in prange(H_points):
        # создаеся локальная копия массива p для каждого потока
        p_local = p.copy()

        # Вызов функци поиска Ic
        Ic_array[k] = find_Ic(k, Hmin, dH, Tmin, Tmax, dt, p_local,
                              delta_I, noisemax, epsilon)

    return Ic_array

delta_I = 0.001
epsilon = 0.05

H_points = 160
Hmin = 0.0
Hmax = 9.0
dH = (Hmax - Hmin) / (H_points - 1)

H_array = np.linspace(Hmin, Hmax, num=H_points)

k = 0
find_Ic(k, Hmin, dH, Tmin, Tmax, dt, p, delta_I, noisemax, epsilon)

1.0010000000000006

set_num_threads(80)

start = time.time()

Ic_array = comput_Ic_on_H_parallel(H_points, Hmin, dH, Tmin, Tmax, dt, p,
                                   delta_I, noisemax, epsilon)

end = time.time()

print("calculation time:", end - start)

calculation time: 138.09818816184998

fig = plt.figure(figsize=(8, 8))
plt.plot(H_array, Ic_array, label='Зависимость $I_c$ от $H_{ext}$',
         linewidth=3.0)
plt.xlabel('$H_{ext}$', size=12)
plt.ylabel('$I_{c}$', size=12)
plt.legend(loc='upper right')
plt.show()
_images/290dd5ba26bb9e3f1afbac88f27bd4a0f381f46e7f0fff78bd787dec70bd663b.png

Рис. 8. Зависимость критического тока от внешнего магнитного поля.

Contents