关于3维球对称势阱中的拉比震荡

计算 Rabi 振荡的振荡频率需要知道偶极跃迁矩阵元，即形如

的矩阵元，分别对应  的情形。而这可以用CG系数表示为（樱井纯的课后题）

其中 R 代表径向波函数。如果引入等效的一维本征波函数

那么就是

可以方便地用 matlab 数值计算第一个括号内的积分，

后面部分查CG系数（）的表再带入就行了

至于含时势能情形的Schrodinger方程数值求解，只用在之前的代码上稍作修改：

1. 修改 UMmt

UMmt = exp(-1i*dT*(Kx.^2+Ky.^2+Kz.^2)/4); % exp(-i*(-0.5*Laplacian)*dT/2), 处于动量空间

2. 球对称势能 Up 固定，总势能修改为

Zr = R.*cos(Tt);

Xr = R.*sin(Tt).*cos(P);

Yr = R.*sin(Tt).*sin(P);

F0 = 1.5;%按需修改

Om = 6.25; %按需修改，应当是能级间距+失谐

Utm = @(z,t) F0*sin(Om*t)*z; % q=0

% Utm = @(x,y,t) F0*(sin(Om*t)*x+cos(Om*t)*y)/sqrt(2); % |q|=1

Unow = Up+ Utm(Zr,0);

% Unow = Up+ Utm(Xr,Yr,0);

3. for 循环中时间演化的部分修改为

Unow = Up+Utm(Zr,t*dT);

% Unow = Up+Utm(Xr,Yr,t*dT);

UVh = exp(-1i*Unow*dT);

%exp(-i*(-0.5*Laplacian)*dT/2)

psik = UMmt.*fftshift(fftn(psi3));

psi3 = ifftn(ifftshift(psik));

% exp(-i*U*dT)

psi3 = UVh.*psi3;

%exp(-i*(-0.5*Laplacian)*dT/2)

psik = UMmt.*fftshift(fftn(psi3));

psi3 = ifftn(ifftshift(psik));

至于Rabi振荡的角频率，零失谐时就是

M 就是之前提到的矩阵元。如果还有非零的失谐，那就是 计算失谐与上式的平方和，再开根号。

矩阵元 M 中积分计算的一个例子：

N1 = 4000;

x1 = linspace(1e-3,L,N1)';

l = 0; l2 = 1; 

m = 0; m2 = 0; 

nl1 = 1;nl2 = 1;

[En1,Vl] = RadEig(x1,SphP(x1),l);

[En2,Vl2] = RadEig(x1,SphP(x1),l2);

I=Vl2(:,nl2)'*((x1).*Vl(:,nl1))

I 就是想要的积分（的近似值）了。

自定义函数 RadEig 见上一个文章