约瑟夫森结I-V特性及非线性的数值模拟
彭加福
（江苏科技大学数理学院，应用物理，0640502112）
摘要：本文基于Matlab对约瑟夫森结（Josephson Junction）RCSJ模型的交直流I-V特性及非线性混沌现象进行数值模拟。
通过计算机数值模拟得到该模型的非线性微分方程数值解，研究了RCSJ模型中各参量对约瑟夫森结的影响，进而简要分析其I-V特性和非线性混沌现象的产生机理，绘制出约瑟夫森结的交直流I-V特性曲线、非线性微分方程的相图及因其高度非线性而引起的通过倍周期分岔和阵发性原理进入混沌状态的分岔图。
关键词：超导器件隧道效应约瑟夫森结弱耦合倍周期分岔庞加莱截面混沌
1．引言
自1911年荷兰科学家昂纳斯（H. K. Onnes）发现汞的超导现象以来，人们对超导进行了大量开拓性的研究，使超导理论]1[日趋成熟，与此同时，超导技术也在各个领域得到深入而广泛的应用]2[。
约瑟夫森效应的发现开拓了超导量子干涉仪（SQUID）在弱电方面的应用。
人们在对约瑟夫森效应进行研究的过程中发明了各种超导器件及应用电路]3[，促使超导技术应用的新领域——超导电子学逐渐发展起来。
在其中，因具有各种独特性（量子干涉、特殊的I-V 特性和高度的非线性等），约瑟夫森结得到广泛的研究和应用，并成为超导电子器件的核心部件。
实际使用中的约瑟夫森结总处于某一电路之中，因此，利用等效电路理论来研究和分析约瑟夫森结的物理行为是一种很有效的方法。
在各模型中，其物理行为均可用微分方程来描述，但这些方程大多不易直接求解析解，因而发展了很多间接解法]5[],4[。
其中，利用电路模拟（RCSJ模型和RSJ模型等等），如图1、图2所示，并用数值计算来研究约瑟夫森结的方法最直接，简易。
图1：RCSJ模型等效电路图2：RSJ模型等效电路Resistively Capacitance Shunted Junction Resistively Shunted Junction
2． 约瑟夫森效应及约瑟夫森结简介
1962年，约瑟夫森（B. D. Josephson ）提出：两块用绝缘薄层隔开且紧密地接近的超导体间，甚至在没有电势差的情况下，电子仍能够穿过绝缘薄层（隧道现象）。
在不到一年的时间内，安德森（P. W. Anderson ）和罗韦尔（J. M. Rowell ）等人从实验上观察到这一现象，证实了约瑟夫森的预言，此即约瑟夫森效应]6[（Josephson Effection ），又称隧道效应。
两块超导体通过一绝缘薄层（厚度为10埃左右）连接起来的组合称S-I-S 超导隧道结或约瑟夫森结。
绝缘层对电子来说是一势垒，一块超导体中的电子能穿过势垒进入另一超导体中，这是约瑟夫森结中特有的量子隧道效应]7[。
绝缘层太厚时，隧道效应不明显，太薄时，两块超导体实际连成一体，这两种情况下均不易观察到约瑟夫森效应，只有当绝缘层不太厚也不太薄，即弱耦合时，才能观察到显著的约瑟夫森效应。
3． 电路模型及计算机模拟
理想约瑟夫森结遵循如下方程：
ϕsin c I I = （1）
V e dt d
2=ϕ （2） 式中，ϕ是结中心量子波函数的相位差,c I 是总临界电流,与外加磁场有关（文献3中有详细介绍）。
等式（1）、（2）所描述的仅仅是电子对携带的电流，从理论上讲，结间还有位移电流、准粒子隧道电流和绝缘体漏电流，研究时应该把这些电流都考虑进去。
研究约瑟夫森结时， 一般都把一个理想的约瑟夫森结与结电阻R 和结电容C 相并联，即用电容模拟位移电流，用电阻模拟准粒子隧道电流和绝缘体漏电流。
这样的模型称为电阻、电容分路结模型，简称RCSJ 模型（Resistively Capacitance Shunted Junction ） ，如图1所示。
对于小面积隧道结、超导微桥以及点接触结，当其结电容很小时，可以不考虑电容的影响。
于是可用电阻分路结模型（Resistively Shunted Junction ）来描述：一个理想的约瑟夫森结与一个电阻相并联，简称RSJ 模型，如图2所示。
采用RCSJ 模型，并同时加上直流电流源和交流电流源后，可用如下方程描述：
dt
dV C R V I t I I c a a d ++=+)sin()sin(ϕω （3）
将（2）式代入上式并用c I 除，再令c d d I I i =，c a a I I i = ， e k 2=，t c ωτ=，C
eI c c 2=ω 即可得到无量纲化后的方程： ϕτϕβτϕωτsin )sin(2
2++=+d d d d i i c a d （4） 在上式中，c a ωωω=，22CR eI c c =β，c β为阻尼参数，ω为归一化的外加交流电流圆频率。
上述无量纲化的非线性微分方程没有确切的解析解，因此，我们将利用Matlab 求其数值解。
求解中采用Matlab 自带的四阶龙格----库塔（Runge-Kutta ）算法]8[，此算法可以在保证高精度的前提下，快速求解微分方程的数值解，其Matlab 函数为ode45()。
理论上已证明，在RCSJ 模型中接入直流电源，当d I =c I 时，约瑟夫森结间电压0≠V ，且电容的存在会使电压的产生滞后]9[。
在RCSJ 模型中接入交流电源时，会得到台阶状的I-V 特性曲线，且台阶是按一定的规律]10[分布的。
用上述的方法建立数值模型，借助Matlab 进行计算、绘图，可以精确、直观地从数值模拟实验上证明和看到上述现象。
此外，RCSJ 的数值模型还是一个很好的非线性混沌现象研究工具，因为我们不仅可以在保证精度的前提下，随意的改变电路的各个参数来研究各参量之间的关系和系统进入混沌的过程，而且可以通过绘制相图、庞加莱截面]11[、时序图和功率谱等对混沌产生机理进行分析。
4． Matlab 模拟结果及分析
a. 约瑟夫森结的直流I-V 特性
由（2）式可知约瑟夫森结平均电压为：
τϕωd d e
U c /2 = （5） 因此，计算机模拟约瑟夫森结伏安特性时可用d i —dt d /ϕ曲线代替I-V 特性曲线。
计算时用函数mean()近似地计算电压平均值，可以得到比较精确的结果，绘出的I-V 特性曲线如图3、图4所示。
图3：约瑟夫森结直流I-V 特性曲线 图4：约瑟夫森结直流I-V 特性曲线
Josephson Junction DC I-V characteristic curve Josephson Junction DC I-V characteristic curve
分析：从图中可以看出，约瑟夫森结注入直流电流后（1=c I ），当电流d I < c I 时，U=0，此为约瑟夫森结的超导态电流；当d I =c I 时是最大临界电流；当d I > c I 时，U ≠0，此时约瑟夫森结中的电流由正常态电流和超导电流组成。
电容的存在并没有改变约瑟夫森结的性质，但会使得电压产生滞后，且随着电容的增大，滞后现象会增强，如图4所示。
b. 约瑟夫森结的交流I-V 特性
图5：约瑟夫森结交流I-V 特性曲线
Josephson Junction AC I-V characteristic curve
当同时向约瑟夫森结注入直流电流d i 和交流电流a i ，交流的I-V 特性曲线上出现台阶现象，如图5所示。
数学计算的结果证明了这些台阶出现在：
f e
n U n 2 = n=0，1，2，… （6） 对于无量纲化的结果来说，交流台阶出现在：
ωτ
ϕn d d n =〉〈 n=0，1，2，… （7） c. 约瑟夫森结的高度非线性特性
从上述无量纲化过程中，可以知道归一化电压u 的表达式为：
c
I eC V d d u 2==τϕ （8） 因此，画分岔图、相图和庞加莱截面时纵轴就是归一化电压u 。
下面给出了通过Matlab 模拟计算得到的非线性混沌现象的分岔图、相图、庞加莱截面及功率谱，也给出了一些简要的分析方法。
i. 通向混沌的道路
图6：系统u 随参数d i 的分岔图 图7：系统u 随参数c β的分岔图
The bifurcation diagram of U with the parameter d i The bifurcation diagram of U with the parameter c β
图8：系统u 随参数c β的分岔图 图9：系统u 随参数a i 的分岔图 The bifurcation diagram of U with the parameter c β The bifurcation diagram of U with the parameter a i
通向混沌的道路有很多种]12[，如倍周期分岔、阵发性方式、茹厄勒——塔肯斯道路等。
它们是具有普适性的，很多动力学系统随着参数的改变会以倍周期分岔或阵发性的方式进入混沌状态，在混沌区中，由于阵发性的作用，又会出现一些具有周期性的窗口。
如图6，图7，图8，图9所示，它们分别是约瑟夫森结电压u 随参数d i 、c β、a i 的分岔图。
我们先简要分析图6和图7，然后在ii 中用相图、时序图、庞加莱截面和功率谱着重分析图9。
图6是以d i 为控制参数(0.8<d i <2.0)，其他参数为(c β=1，5.0=ω，65.0=a i )时，得到的u 随d i 变化的分岔图。
当0.8<d i <1.15时，系统处于1周期状态；当d i 等于1.15时，系统进入2周期状态，接着是4周期状态、6周期状态、4周期状态；当1.24<d i <1.43时，系统又处于2周期状态；当1.43<d i <1.6时，系统处于内部有周期窗口的混沌状态；当1.6<d i <1.68 时，系统是3周期状态；当1.68<d i <2.0时，系统处于混沌状态，且内部有周期窗口出现。
从图6我们看到了不完全级联的过程，即随直流电流的增加，系统由1周期解(单周期的运动)进入到2，4，6周期解(多周期的运动)，然后返回到4，2周期解，之后才进入混沌状态，而没有直接进入混沌，这称为不完全级联；在选择的参数范围内，系统随d i 的增加以阵发方式进入混沌状态，且内部有周期窗口。
图7，图8（局部放大图）是以c β为控制参数(0<c β<1.5)，固定其他参数(0=d i ，2.0=ω，50.1=a i )时，得到的u 随c β变化的分岔图。
从图中很容易看出：随着参数c β的减小，周期运动状态和混沌状态会不断交替出现，且混沌区中含有较狭窄的周期窗口，最后完全进入混沌状态；陆续出现周期1，周期3，周期5，周期7，…的周期状态区 ，且周期状态区域逐渐变窄，直到消失；在选择的参数范围内，系统随c β的减小以阵发性方式进入混沌状态；随着参数c β的减小，系统由周期状态进入到混沌状态的过程是级联过程。
ii. 混沌现象的分析
图9是约瑟夫森结电压u 随参数a i 的分岔图（20<<a i ）。
由图中可以看出，它与前面两种情况不同，该条件下系统是通过倍周期分岔到混沌状态，再由阵法性方式到达周期状态的。
由于单用一种方法对混沌进行研究时，在一些情况下并不能达到很好的效果。
比如相轨道很密时不易分出是几个周期；画庞加莱截面时两点靠得太近而误以为是单周期，等等。
因
此，我们用四种方法相互比较对混沌进行研究，分别是相图、庞加莱截面、功率谱和时序图。
周期状态区、混沌状态区和周期窗口在分岔图中显而易见，在此就不列举了。
我们截取有代表性的几个值进行研究，得到如下图像。
图10：系统相图（Bc=0.5,w=0.66,ia=0.8）
Phase Diagram of The System（Bc=0.5,w=0.66,ia=0.8）
图11：系统相图（Bc=0.5,w=0.66,ia=1.0636）
Phase Diagram of The System（Bc=0.5,w=0.66,ia=1.0636）
图12：系统相图（Bc=0.5,w=0.66,ia=1.0652）Phase Diagram of The System（Bc=0.5,w=0.66,ia=1.0652）
图13：系统相图（Bc=0.5,w=0.66,ia=1.08）Phase Diagram of The System（Bc=0.5,w=0.66,ia=1.08）
图14：系统相图（Bc=0.5,w=0.66,ia=1.212）Phase Diagram of The System（Bc=0.5,w=0.66,ia=1.212）
图15：系统相图（Bc=0.5,w=0.66,ia=1.32）Phase Diagram of The System（Bc=0.5,w=0.66,ia=1.32）
分析：从分岔图中可以看到，a i <1时系统处于单周期状态，在a i =1附近因阵发性而进入二周期，然后因倍周期分岔进入混沌状态，图10～13是取这一过程中的几个特殊点绘制的分析图像。
可以看出图10、图11、图12分别表示单周期、二周期和四周期，图13则表示系统进入混沌状态。
图13与图14比较，前者的庞加莱截面还是清晰的线，而后者则是分片密集的点阵；前者功率谱有局部出现噪声背景（连续普），而后者都是连续普。
这是因为a i =1.08时，倍周期分岔还在进行中，还有一点规律可循，而a i =1.212时，系统已经完全进入混沌状态。
图15表明a i =1.32时，系统又由阵发性进入周期状态，然后再次进入混沌状态，其间有微小扰动（从时序图中看出）。
5． 结束语
此工作通过计算机模拟约瑟夫森结（RCSJ 模型）的I-V 特性及非线性混沌现象，给出了：约瑟夫森结的交直流I-V 特性曲线；系统非线性混沌分岔图和与之相关的相图，庞加莱截面，功率谱及时序图；系统进入混沌状态的方式及混沌的一些分析方法。
通过这次仿真实验，我从中学到了很多物理知识和计算机仿真经验，深感浅尝辄止与深入研究的巨大差别。
6． 微分方程数值解Matlab 程序
*PART ONE:
****** Josephson_Junction.m by J. F . Peng ******
%
function dy=Josephson_Junction(t,y)
% 在程序中参数t 代表文中的τ
% Bc 为阻尼参数, w 为归一化的外加交流电流圆频率
% Bc 、ia 、id 、w 、τ均为无量纲化时引入的变量，其中ia 、id 为各电流无量纲化时引入的 % y(1)= φ 约瑟夫森结中心处波函数ψ1、ψ2的相位差
% y(2)=d φ/dt 即归一化电压u
% 以下为约瑟夫森结RCSJ 模型的无量纲化微分方程
% dφ/dt=u
% du/dt=-Bc*u-sin(φ)+id+ia*sin(w*t)
%
% ****** Author:J.F.Peng(jingyujiafu@) ******
%
global Bc ia w id;
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=id+ia*sin(w*t)-sin(y(1))-Bc*y(2);
*PART TWO:
****** plotACIV.m by J. F. Peng ******
%
% RCSJ模型的直流I-V特性
%
% ****** Author:J.F.Peng(jingyujiafu@) ******
%
clear all;
clc
tic
global Bc ia w id;
a=[0.60,1.00];
for i=1:2
Bc=a(i);
w=0.50;
ia=0.00;
v=[];
for id=-2:0.005:2
[T,Y]=ode45('Josephson_Junction',[0:0.01:80],[0 0 0]); v=[v,mean(Y(:,2))];
end
plot_id=(-2:0.005:2);
plot(v,plot_id,'.b','markersize',4)
% 将Bc=1.00和Bc=0.60的直流I-V特性曲线画一起，好比较
hold on
end
plot_id=(-2:0.005:2);
plot(v,plot_id,'.r','markersize',4)
xlabel('Volt age U~<dφ/dt>')
ylabel('Directive Current id')
desc={'Bc=0.60(blue)||1.00(red)',
'w=0.50, ia=0.00, Ic=1'};
text(-3,1.25,desc)
figure
plot_id=(-2:0.005:2);
plot(v,plot_id,'k.','markersize',4)
xlabel('Voltage U~<dφ/dt>')
ylabel('Directive Current id')
desc={'Bc=1.00, w=0.50, ia=0.00, Ic=1'};
text(-1.75,1.25,desc)
%
% RCSJ模型的交流I-V特性
%
clc
global Bc ia w id;
Bc=1.00;
w=1.00;
ia=2.00;
v=[];
for id=-4:0.005:4
[T,Y]=ode45('Josephson_Junction',[0:0.01:80],[0 0 0]); v=[v,mean(Y(:,2))];
end
plot_id=(-4:0.005:4);
figure
plot(v,plot_id,'k.','markersize',4)
xlabel('Voltage U~<dφ/dt>')
ylabel('Directive Current id')
desc={'Bc=1.00 w=1.00, ia=2.00'};
text(-3,2.5,desc)
toc
t=toc
*PART THREE:
****** plotphase.m by J. F. Peng ******
%
% 画时序图、相图、庞加莱截面及功率谱
%
% ****** Author:J.F.Peng(jingyujiafu@) ******
%
clear
clc
global Bc ia w id;
Bc=0.500;
ia=1.32;
w=0.660;
id=0.000;
%
% 解微分方程
%
[t,Y]=ode45('Josephson_Junction',[0,800],[0;0;0]);
%
% 画φ、dφ/dt、ψ（dψ/dt=w）的三维图
%
figure
plot3(Y(:,1),Y(:,2),Y(:,3))
xlabel('Phase φ');
ylabel('Voltage u=dφ/dt');
zlabel('Phase of Aternative Current ψ');
%
% 画归一化电压的时序图(Time Dependence Voltage)
%
miny=min(Y(:,2));
maxy=max(Y(:,2));
figure
subplot(2,2,4)
plot(t(700:end),Y(700:end,2))
axis([700 800 miny maxy])
xlabel('Time t');
ylabel('Voltage u=dφ/dt');
title('Time Dependence Voltage')
%
% 画相图(Phase Diagram)
%
T=2*pi/w;
[t,Y]=ode45('Josephson_Junction',[0:T/100:1000*T],[0;0;0]); n=length(t);
m=round(n/60);
subplot(2,2,1)
%plot(Y(m:n,1),Y(m:n,2),'k')
%plot(sin(Y(m:n,1)),Y(m:n,2),'k')
%
% 如果相图窗口被涂黑可选择用下个
%
plot(Y(50*m:n,1),Y(50*m:n,2),'k')
%plot(sin(Y(50*m:n,1)),Y(50*m:n,2),'k')
axis tight
%axis([560 600 -2 2])
axis([-5600 -5550 -2 1])
title('Phase Diagram')
xlabel('Phase φ')
%xlabel('sin(φ)')
ylabel('Voltage u=dφ/dt')
%
% 画庞加莱截面(Poincare Section)
%
% 通过观察Poincare截面上截点的情况可以判断是否发生混沌：
% 当Poincare截面上有且只有一个不动点或少数离散点时，运动是周期的；
% 当Poincare截面上是一封闭曲线时，运动是准周期的；
% 当Poincare截面上是一些成片的具有分形结构的密集点时，运动便是混沌。
%
Z1=[];
Z2=[];
%
% 每隔一个外加交流电周期(T=2π/ω)取一个点
%
for i=m:100:length(Y)
Z1=[Z1,Y(i,1)];
Z2=[Z2,Y(i,2)];
end
subplot(2,2,2)
plot(Z1,Z2,'b.','markersize',6)
minx=min(Y(m:n,1));
maxx=max(Y(m:n,1));
miny=min(Y(m:n,2));
maxy=max(Y(m:n,2));
axis([minx maxx miny maxy])
title('Poincare Section with（T=2π/ω）')
xlabel('Phase φ')
ylabe l('Voltage u=dφ/dt')
%
% 画功率谱(Power Spectrum)
%
subplot(2,2,3)
%
% Fs 采样频率
% CYn 自相关函数
%
Fs=1;
t=0:1/Fs:1000;
[T,Y]=ode45('Josephson_Junction',t,[0;0;0]);
n=length(t);
m=round(n/6);
Yn=Y(m:n,2);
nfft=1024;
CYn=xcorr(Yn,'unbiased');
CYk=fft(CYn,nfft);
Pyy=abs(CYk);
index=0:round(nfft/2-1);
%
% 功率谱以f=Fs/2为对称轴，左右对称，故画图频率范围为0~Fs/2 %
k=index*Fs/nfft;
plot_Pyy=1000*Pyy(index+1);
plot(k,10*log10(plot_Pyy));
title('Power Spectrum')
xlabel('Frequency f')
ylabel('power |Yk|^2')
axis tight
*PART FOUR:
****** plotfork.m by J. F. Peng ******
%
% 分岔图(Fork)
%
% ****** Author:J.F.Peng(jingyujiafu@) ******
%
clear
clc
%
% tic、toc用来计算程序运行时间，分别表示开始和结束计时
%
tic
global Bc ia w id;
id=0.00;
w=0.20;
ia=1.50;
for Bc=0.00:0.002:2
[T,Y]=ode45('Josephson_Junction',[0,600],[0;0;0]); data=Y(:,2);
n=length(data);
% 去掉一些前面的点，防止因计算机配置问题而画不出图，不影响结果 m=round(n/6);
a=data(m-2);
b=data(m-1);
for i=m:n
if b>=a&b>=data(i);
plot(Bc,b,'k');
hold on;
end
a=b;
b=data(i);
end
end
%
% 根据需要选择图像标注语句
%
% ia-u 分岔图（fork）
%
%xlabel('Aternative Current ia')
%ylabel('Voltage u=dφ/dt')
%desc={'Bc=0.50,w=0.66,id=0.00'};
%text(0.5,3,desc)
%
% id-u 分岔图 (fork)
%
%xlabel('Directive Current id')
%yl abel('Voltage u=dφ/dt')
%desc={'Bc=1.00,w=0.50,ia=0.65'};
%text(1,3,desc)
%
% Bc-u 分岔图 (fork)
%
xlabel('Damping Coefficient \betac')
ylabel('Voltage u=dφ/dt')
desc={'ia=1.50,w=0.20,id=0.00'};
text(0.75,8,desc)
toc
t=toc
参考文献：
[1] 韩汝珊. 高温超导物理[M]. 北京: 北京大学出版社,2002.10.
[2] 王惠龄,汪京荣. 超导应用低温技术[M]. 北京: 国防工业出版社,2008.1.
[3] T. Van杜塞尔,C. W. 特纳. 超导器件与电路原理[M]. 北京: 科学技术文献出版社,1986.5.
[4] 周铁戈等. 本征约瑟夫森结阵列的PSpice模型及混沌行为研究[J]. 物理学报. 2007,V ol.56,No.11:6307-6314.
[5] 徐勤军等. 高温超导Josephson结的微波响应分析[J]. 上海大学学报(自然科学版). 2000,V ol.6,No.3:237-242.
[6] 大唐搜索. 约瑟夫森效应. /Y/Y1958.htm
[7] 百度百科. 隧道效应. /view/152179.htm
[8] John H. Mathews,Kurtis D. Fink. Numerical Methods Using MA TLAB(Fourth Edition)[M]. Beijing: Publishing
House of Electronics Industry,2005.12.
[9] 张卉,胡毅飞等. 基于Matlab的约瑟夫森结I-V特性数值分析[J]. 宇航计测技术. 2007,Vol.27,No.1:15-18.
[10] 戴岭,于瑶等. 约瑟夫森结非线性特性的计算机模拟[J]. 物理实验. 2005,Vol.25,No.2:25-30.
[11] 卢新文. 自激振动的计算机模拟[J]. 甘肃联合大学学报(自然科学版). 2004, Vol.18,No.4:31-33.
[12] 黄润生,黄浩. 混沌及其应用[M]. 武昌: 武汉大学出版,2005.12.
致谢：此工作得到陈莉莉老师在非线性物理方面的指导，在此表示感谢！
江苏科技大学数理学院
2009年3月27日——2009年5月3日。
