到1988年的时候，非线性光波耦合，被人为地分为连续光波或准连续光波耦合和光孤子耦合。那时候，有两个理论，一个用来分析连续或准连续光波的耦合，这就是所谓的Jensen 方程，另外一个用来分析光孤子耦合，这就是非常有名的孤子耦合理论，孤子耦合方程。 这个方程由 S. Trillo, S. Wabnitz, E. M. Wright and G. I. Stegman在1988年给出。 这几个人，因为这一工作，成为我们现在年轻人讲的所谓牛人，尤其是Stegman,大牛。在此之后，孤子耦合，顿时洛阳纸贵，发表了成千上万的论文，有好多论文发表在我们国人引以为豪的PR 和PRL上。 这些论文都是使用孤子耦合方程，并且基于孤子耦合理论或假设: 孤子可以不畸变地在非线性耦合器中来回耦合。我们之后证明，这一光孤子重要理论是不成立的。大量的论文实际上，是基于一个错误的假设或理论。 到大约2002年的时候，我们注意到，连续光波耦合理论和光孤子耦合理论之间有间隙，不具有兼容性，不能从光孤子耦合方程推到连续光波耦合方程。理论上讲，描述同样物理规律的理论之间，必须有兼容性。从一个理论，在极限条件下，要能推导出另外一个理论。比如波动光学和几何光学，在波长为0 的条件下，波动光学可以导出几何光学。比如相对论，量子力学和经典力学等等。当时，我们想应该有一个中间的光波耦合方程，将他们联系在一起。我们就开始着手用不同的归一化方式来检讨他们之间的联系。很快，我们就遇到一个问题，就是如何求解非线性偏微分耦合方程组。分步傅立叶法只适合求解非线性偏微分方程，而不能用来求解非线性偏微分耦合方程组。我们知道用差分方法可以求解非线性偏微分耦合方程组，实际上，Stegman等人就是使用的有限差分法来分析光孤子耦合的。但是，我们对差分方法不放心。做过数值计算的人都知道，差分法有其局限性：他们的解有时可能是错误的。我们也不习惯用纯数值计算方法，以前我们所有理论研究结果几乎全部是以解析的公式给出的。因此，我们决定先寻找一种非纯数值计算的计算方法。经过一年左右的尝试，功夫不负有心人，在2003年5月，我们终于成功地找到了一种简单而可靠的算法。这个算法在2004年，发表到IEEE Photonics Technologiess Letter上。在寻找算法的过程中，我们的自信以及确信一定可以找到一种方法的意志，起了很重要的作用。 完成算法发现之后，我们重新着手对描述非线性脉冲耦合的非线性薛定格耦合方程的重新归一化。因为我们已经知道我们的目标是：归一化的耦合方程可以在一定条件下退化为 Jensen 方程，所以我们只有有限的几种归一化方式去试一试。尝试了几个月之后，我们就找到了一种归一化方法，以及相应的非线性耦合方程组。当我们使用这组新的非线性耦合方程去研究光脉冲耦合和光孤子耦合时，使我们大吃一惊的是，原来的孤子耦合理论的两个重要理论或假设是错误的。 他们是：1）孤子可以不畸变地在非线性耦合器开关中来回耦合，其原因是由于非线性效应和色散平衡。 这是 S. Trillo, S. Wabnitz, E. M. Wright and G. I. Stegman在他们1988年那篇著名的孤子耦合的论文中给出的一个当时非常重要结论。他们 根据自己的一个仿真图，得到了这一结论。事实上，作者仿真的区域属于弱非线性区或线性区，在这个区域，非线性效应不能改变耦合行为，没有起到作用。它不是非线性耦合，而是线性耦合。光孤子在线性耦合器中来回耦合没什么奇怪，没人特别关注这个。采用孤子耦合方程，人们无法判断他们所研究的区是真的非线性或只是线性区。这是孤子耦合方程的一大缺陷。2）光开关耦合效率取决于光脉冲形状，孤子可以提高非线性耦合器的开关效率。光孤子耦合之所以，一出来就受到青睐，就是因为这两个假设或特点。事实上，我们后来发现，所有光脉冲都可以得到同样的开关效率，只要开关的特征参数满足一定条件。 我们发现非线性光脉冲的耦合行为由一个特征参数决定，这样，使得非线性耦合行为的理解和应用变得非常简单明了。这篇文章2004年发表在Applied physics B 上。Applied physics B的主编特意评审过我们这篇文章。在这篇文章之前，成千上万的论文对孤子耦合进行了研究，结果只是盲人摸象，管中窥豹，根本没有办法描绘一个脉冲或孤子耦合的物理图像，并且存在很多错误的理解和结论，孤子开关的条件都没法正确的给出。现在有些书上给出的孤子耦合条件，就是错误的。 因为孤子耦合只是讨论了，二阶色散。接下来，我们研究了高阶色散，模间色散，以及高阶非线性和非线性延迟效应对光脉冲耦合的影响。我们花了大概两年的时间，着手将新的归一化非线性耦合方程推广到含有高阶非线性和高阶色散。到2006年的时候，超快光脉冲的耦合和孤子耦合的新归一化方程就完善了。这个方程有一个特征参量，当特征参量远远大于1时，方程就退化为Jensen 方程。除此之外，我们还发现，一个光脉冲是否属于超快脉冲还是准连续波，也取决于这一个特征量。一个皮秒脉宽的脉冲，也可以遵从连续波方程。这样使得，早期从连续波得到的一些应用，也同样适合于皮秒脉冲。 根据这些工作，我们写了两篇论文，这两篇文章在2006年发表在IEEE/OSA Journal of lightwave technology. Agrawal教授在他新版的教科书上简单介绍我们的工作。之后，我们也将方程推广到，偏振耦合。在此之后和同时，我们有运用这一理论，发现一些新的效应，比喻脉宽控制的自开关效应。 在1988年时， S. R. Friberg, A. M. Weiner, Y.Silberber, B.G.Sfez and P.S.Smith在非线性耦合器中，用脉宽为100fs 到200fs 光脉冲成功地实现了全光开关。实验结果与连续波光开关理论吻合，而与孤子耦合理论相去甚远。无论从孤子理论的角度来讲，还是从连续波非线性耦合理论的角度来讲，这个结果都不合逻辑，没法解释。但是，按照我们最近这几年发展的非线性光脉冲耦合理论讲，理论与实验相符合。按照我们的理论，当光开关的特征参数Kd*k1时，光脉冲非线性耦合方程退化为Jensen 方程，光脉冲耦合行为与连续波非线性耦合行为一致。我们不难发现，在Friberg等人的实验中，光开关的特征参量Kd*k=601.所以实验结果应该符合连续波非线性耦合的结果。
这个程序用于仿真孤子的演变波形，这里Lcd 可以取不同的值。 % This Matlab script file solves the nonlinear Schrodinger equation % simulate Fig.2 and Fig.3in paper: Y. Wang and W.Wang, “Study of ultra-fast optical pulse propagation in a % nonlinear directional coupler”, Appl. Phys. B79, p51-55 (2004) %Here the coupled soliton equationsare devided by Lcd,so that the length becomes physical length % not normalized length. N = 1024; % Number of Fourier modes (Time domain sampling points) dz =0.5*0.0031415926536/4;% space step, make sure nonlinear0.05 M3=8000; J =100; % Steps between output of space T =30; % length of time:T*T0 lcd=1;% LCD is CL_D in the paper b=1.; MN3=0; dt = T/N; % time step n = '; % Indices t = n.*dt; ww = 4*n.*n*pi*pi/T/T; % Square of frequency. Note i^2=-1. w=2*pi*n./T; g1=-i*ww./(2*lcd); g2=-i*ww./(2*lcd); % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T %s1=1.5*exp(-1.5*1.5*1.5*1.5*t.*t.*t.*t./((8.25/3.14159)*(8.25/3.14159))/0.36); s1=b*sech(b*t); % input s2=s1.*0.; s10=s1; s20=s10*0.; S1=s1; S2=s2;P=1; p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s1(1,1)*s1(1,1))));%energy in waveguide 1 p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1))));%energy in waveguide 2 P1=p10/p10; P2=p20/p10; for m3 = 1:1:M3 % Start space evolution s1 = exp(dz*i*(abs(s1).*abs(s1))/lcd).*s1; % Solve nonlinear part of NLS s2 = exp(dz*i*(abs(s2).*abs(s2))/lcd).*s2; sca1 = fftshift(fft(s1)); % Take Fourier transform sca2 = fftshift(fft(s2)); sc2=exp(g2.*dz).*(sca2+i*sca1.*dz); sc1=exp(g1.*dz).*(sca1+i*sca2.*dz); % frequency domain phase shift % Advance in Fourier space s2 = ifft(fftshift(sc2)); % Return to physical space s1 = ifft(fftshift(sc1)); p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s10(1,1)*s10(1,1)))); p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s20(1,1)*s20(1,1)))); if rem(m3,J) == 0 % Save output every J steps. S1 = ; % put solutions in U array S2= ; P1= ; P2= ; P= ; MN3= ; z3=dz*MN3'; % output location end end figure(1) waterfall(t',z3',abs(S1').*abs(S1')) %t' is 1xn, z' is 1xm, and U1' is mxn figure(2) waterfall(t',z3',abs(S2').*abs(S2')) %t' is 1xn, z' is 1xm, and U1' is mxn figure(3) plot(P,P1,P,P2); figure(4) plot(t, abs(s10').*abs(s10'),'.', t, abs(s2').*abs(s2'),'--', t, abs(s1').*abs(s1'))
基于求解非线性薛定格方程的分步傅立叶方法和参变量积分法，在该论文中，我们提出了一种求解光脉冲在N芯光纤耦合器中传输信号的耦合模方程组的数值方法，同时用此方法对光脉冲在双芯和三芯非线性耦合器中的传输演化进行了仿真。 原文请参见； Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear pulse propagation in N-core optical couplers”, IEEE Photo. Tech. lett. Vol.16, No.4, p1077-1079(2004) 该文可以在这里下载； http://www.qiji.cn/eprint/abs/1022.html % This Matlab script file solves the coupled nonlinear Schrodinger equations of % soliton in 2 cores coupler. The output pulse evolution plot is shown in Fig.1 of % Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear % pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 %fid=fopen('e21.dat','w'); N = 128; % Number of Fourier modes (Time domain sampling points) M1 =3000; % Total number of space steps J =100; % Steps between output of space T =10; % length of time windows:T*T0 T0=0.1; % input pulse width MN1=0; % initial value for the space output location dt = T/N; % time step n = '; % Index t = n.*dt; u10=1.*sech(1*t); % input to waveguide1 amplitude: power=u10*u10 u20=u10.*0.0; % input to waveguide 2 u1=u10; u2=u20; U1 = u1; U2 = u2; % Compute initial condition; save it in U ww = 4*n.*n*pi*pi/T/T; % Square of frequency. Note i^2=-1. w=2*pi*n./T; g=-i*ww./2; % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T L=4; % length of evoluation to compare with S. Trillo's paper dz=L/M1; % space step, make sure nonlinear0.05 for m1 = 1:1:M1 % Start space evolution u1 = exp(dz*i*(abs(u1).*abs(u1))).*u1; % 1st sSolve nonlinear part of NLS u2 = exp(dz*i*(abs(u2).*abs(u2))).*u2; ca1 = fftshift(fft(u1)); % Take Fourier transform ca2 = fftshift(fft(u2)); c2=exp(g.*dz).*(ca2+i*1*ca1.*dz); % approximation c1=exp(g.*dz).*(ca1+i*1*ca2.*dz); % frequency domain phase shift u2 = ifft(fftshift(c2)); % Return to physical space u1 = ifft(fftshift(c1)); if rem(m1,J) == 0 % Save output every J steps. U1 = ; % put solutions in U array U2= ; MN1= ; z1=dz*MN1'; % output location end end hg=abs(U1').*abs(U1'); % for data write to excel ha= ; % for data write to excel t1= ; hh= ; % for data write to excel file %dlmwrite('aa',hh,'\t'); % save data in the excel format figure(1) waterfall(t',z1',abs(U1').*abs(U1')) % t' is 1xn, z' is 1xm, and U1' is mxn figure(2) waterfall(t',z1',abs(U2').*abs(U2')) % t' is 1xn, z' is 1xm, and
标签: 孤子仿真
