这篇文章将介绍微腔系统中的两个光学模式和同一个布里渊力学模式之间的耦合

基本方程

基本示意图
$$
\begin{aligned}
H&=\hbar\omega_1 a^\dagger_1 a_1+\hbar\omega_2 a^\dagger_2 a_2+\hbar\Omega_m b^\dagger b+\hbar g(a^\dagger_1 b^\dagger a_2+a_1 b a^\dagger_2)\\
&+i\hbar\sqrt{\kappa_{ex,1}}s_{in,1}(a^\dagger_1 e^{-i\omega_{c1} t}-a_1 e^{i\omega_{c1} t})\\
&+i\hbar\sqrt{\kappa_{ex,2}}s_{in,2}(a^\dagger_2 e^{-i\omega_{c2} t}-a_2 e^{i\omega_{c2} t})
\end{aligned}
$$
采取变换$U=\exp{(i\omega_{c1}a_1^\dagger a_1 t+i\omega_{c2}a_2^\dagger a_2 t+i(\omega_{c2}-\omega_{c1})b^\dagger b t})$
最后得到的新哈密顿量为
$$\begin{aligned}
H&=-\hbar\Delta_1a_1^\dagger a_1-\hbar \Delta_2 a_2^\dagger a_2-\hbar (\omega_{c2}-\omega_{c1}-\Omega_m)b^\dagger b\\
&+\hbar g(a_1^\dagger b^\dagger a_2+a_1 b a_2^\dagger)\\
&+i\hbar\sqrt{\kappa_{ex,1}}s_{in,1}(a_1^\dagger-a_1)\\
&+i\hbar \sqrt{\kappa_{ex,2}}s_{in,2}(a_2^\dagger-a_2)
\end{aligned}$$
设$\Delta_1=\omega_{c1}-\omega_1,\Delta_2=\omega_{c2}-\omega_2,\Delta_m=\omega_{c2}-\omega_{c1}-\Omega_m$

代入运动学方程
$$
\begin{aligned}
\frac{da_1}{dt}&=(i\Delta_1-\kappa_1/2) a_1-igb^\dagger a_2+\sqrt{\kappa_{ex,1}}s_{in,1}\\
\frac{da_2}{dt}&=(i\Delta_2-\kappa_2/2) a_2-iga_1 b+\sqrt{\kappa_{ex,2}}s_{in,2}\\
\frac{db}{dt}&=(i\Delta_m -\Gamma_m)b-iga_1^\dagger a_2\\
\end{aligned}$$
将变量分成实部和虚部，得到
$$
\begin{aligned}
\frac{da_1^R}{dt}&=-\Delta_1 a_1^I-\kappa_1 a_1^R/2+g(b^R a_2^I-b^I a_2^R)+\sqrt{\kappa_{ex,1}}s_{in,1}\\
\frac{da_1^I}{dt}&=\Delta_1a_1^R -\kappa_1a_1^I/2-g(b^Ra_2^R+b^Ia_2^I)\\
\frac{da_2^R}{dt}&=-\Delta_2a_2^I-\kappa_2a_2^R/2+g(a_1^Rb^I+a_1^Ib^R)+\sqrt{\kappa_{ex,2}}s_{in,2}\\
\frac{da_2^I}{dt}&=\Delta_2a_2^R-\kappa_2a_2^I/2-g(a_1^Rb^R-a_1^Ib^I)\\
\frac{db^R}{dt}&=-\Delta_mb^I-\Gamma_mb^R+g(a_1^Ra_2^I-a_1^Ia_2^R)\\
\frac{db^I}{dt}&=\Delta_mb^R-\Gamma_mb^I-g(a_1^Ra_2^R+a_1^Ia_2^I)\\
\end{aligned}
$$

考虑透射谱

$a_1$的透射谱

我们有输入输出关系$s_{out,1}=s_{in,1}-\sqrt{\kappa_{ex,1}}a_1$，由于我们要模拟脉冲输入的情况，所以输入会出现零的情况，那么我们就不直接计算透射系数，而计算$s_{out,1}$,于是，我们将前面的运动方程进行一定的变形
$$
\begin{aligned}
\frac{dA_1}{dT}&=(i\bar\Delta_1-\kappa_1’/2)A_1-iGB^\dagger A_2/r+s_{in,1}\\
\frac{dA_2}{dT}&=(i\bar\Delta_2-\kappa_2’/2)A_2-iGBA_1 r+r^2 s_{in,2}\\
\frac{dB}{dT}&=(i\bar\Delta_m-\Gamma_m’)B-iGA_1^\dagger A_2/r
\end{aligned}$$
其中
$$\begin{aligned}
T&=\kappa_{ex,1}t,&r=\sqrt{\frac{\kappa_{ex,2}}{\kappa_{ex,1}}},\quad &G=g/\kappa_{ex,1}^{3/2}\\
A_1&=\sqrt{\kappa_{ex,1}}a_1,&A_2=\sqrt{\kappa_{ex,2}}a_2,\quad &B=\sqrt{\kappa_{ex,1}}b\\
\bar\Delta_1&=\Delta_1/\kappa_{ex,1},&\bar\Delta_2=\Delta_2/\kappa_{ex,1},\quad &\bar\Delta_m=\Delta_m/\kappa_{ex,1}\\
\kappa_1’&=\kappa_1/\kappa_{ex,1},&\kappa_2’=\kappa_2/\kappa_{ex,1},\quad &\Gamma_m’=\Gamma_m/\kappa_{ex,1}
\end{aligned}$$

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
from pylab import *
import numpy as np
from scipy.integrate import ode
from manipulate import mplot
def gaussian(x,mu,sig):
    # try:
    #     temp=copy(x)
    #     temp[temp>mu]=mu
    #     return np.exp(-np.power(temp-mu,2.)/(2*np.power(sig,2.)))
    # except:
    #     if x>=mu:
    #         return 1
    #     else:
    #         return np.exp(-np.power(x-mu,2.)/(2*np.power(sig,2.)))
    return np.exp(-np.power(x-mu,2.)/(2*np.power(sig,2.)))
def fun(t,y,arg):
    d1,d2,dm,k1,k2,gm,g,r,s1,s2,mu1,sig1,mu2,sig2=arg
    s1=s1*gaussian(t,mu1,sig1)
    s2=s2*gaussian(t,mu2,sig2)
    return [(1.0j*d1-k1/2)*y[0]-1.0j*g*conj(y[2])*y[1]/r+s1,
            (1.0j*d2-k2/2)*y[1]-1.0j*g*y[2]*y[0]*r+r**2*s2,
            (1.0j*dm-gm)*y[2]-1.0j*g*conj(y[0])*y[1]/r]
def jac(t,y,arg):
    d1,d2,dm,k1,k2,gm,g,r,s1,s2=arg
    return [[1.0j*d1-k1/2,-1.0j*g*conj(y[2])/r,-1.0j*g*y[1]/r],
            [-1.0j*g*y[2]*r,1.0j*d2-k2/2,-1.0j*g*y[0]*r],
            [-1.0j*g*y[1]/r,-1.0j*g*conj(y[0])/r],1.0j*dm-gm]
def time_domain(t,dd=0,d1=0,k1=2,gm=10,g=1,r=1,s1=2.2,s2=0.5,mu1=2.5,sig1=1,mu2=4,sig2=1.6):
    para=[d1,dd,dd,k1,k1,gm,g,r,s1,s2,mu1,sig1,mu2,sig2]
    tran=ode(fun).set_integrator('zvode')
    y0,t0=[0,0,0],0
    tran.set_initial_value(y0,t0).set_f_params(para)
    out=[]
    dt=t[1]-t[0]
    for i in t:
        tran.integrate(tran.t+dt)
        out.append(tran.y)
    out=array(out)
    s=ones(out.shape)
    s[:,0]=s1*gaussian(t,mu1,sig1)
    s[:,1]=s2*gaussian(t,mu2,sig2)
    # s[t<t_stop,1]=s2
    ou=zeros([len(t),4])
    ou[:,0:2]=(abs(s-out)**2)[:,0:2]
    ou[:,2:]=s[:,0:2]
    return ou
mplot(time_domain,linspace(0,20,400),dd=(-75,75),gm=(0.1,10),g=(0.1,10),s1=(0,5),s2=(0,5),mu1=(1,15),mu2=(1,15))
show()

manipulate.py文件