第2節.超短脈衝在色散介質中的傳播模擬

光束在介質中傳播，主要是考慮介質折射率以及吸收係數。由於超短脈衝具有豐富的頻率成分，而不同頻率的光感受的介質折射率會不同，因此他們的速度也不同，因此超短脈衝的時域波形也會發生變化。

同時在傳播過程中，脈衝可能會發射衍射，在這一節中，我們先不考慮脈衝的空間分布，也就是當作超短脈衝是無窮大的平面波，因此不會發生衍射，空間分布也不會發生變化。

理論基礎：

1.光在介質中的傳播：

我就不從麥克斯韋方程組開始了，反正最後都是為了得到波動方程以及本徵函數。由於我們目前考慮的不涉及非線性效應，也就是我們的系統是線性系統。

可以得到本徵函數是 。這裡的 是因為不同的頻率下，波矢是不同的，這跟材料有關。所以可以把超短脈衝分解成一系列 （這裡的分解就是用到上一節介紹的傅里葉變換），然後根據不同 計算 ，得到傳播距離 後的傅里葉變化，即 。再逆變化回來就是我們要求的結果。

2.超短脈衝的描述：實數描述和複數描述。實數描述和複數描述是我的叫法，其實我們一般直接就用複數描述了，但是其實需要一個過度，所以我加了一個叫實時描述的內容。下面我先用實數描述去寫一個模擬，然後再用複數去寫。就可以知道區別了。實數描述的過程在我看來更物理一些，但是，由於需要的計算空間比較大，所以需要用到複數描述。

因為我們這次的代碼是描述具體的物理過程，所以需要使用到物理量，那就不可避免涉及單位的換算，所以就有了下面這個文件units.m。

%% 單位% timefs = 10^-15;ps = 10^-12;ns = 10^-09;us = 10^-06;ms = 10^-03;s = 1;% frequencyTHz = 10^12;GHz = 10^9;MHz = 10^6;KHz = 10^3;Hz = 1;% lengthnm = 10^-9;um = 10^-6;mm = 10^-3;cm = 10^-2;m = 1;% energyKJ = 1e3;J=1;mJ=1e-3*J;uJ=1e-6*J;nJ=1e-9*J;% poweruW = 1e-6;mW = 1e-3;W=1;KW=1e3;MW=1e6;GW=1e9;% angledeg = pi/180;rad = 1;

我們知道，電場實際上是實數物理量，而光是震蕩的光場。那麼超短脈衝是一個由超快包絡的一個震蕩。可以這樣子描述， 。

我在下面代碼中，我 ， 。選擇的 是800nm對應的頻率。

clcclear%%units;c = 299792458*m/s;[t,ft] = domain(200*fs,2048);duration = 12*fs;f_0 = c/(800*nm);A = exp(-(t./duration).^2);E_t = A.*cos(2*pi*f_0*t);plot(fftshift(t)/fs,fftshift(E_t));

我們先來看看真空中脈衝傳播， 。傳播距離12um。可以看到，脈衝向後移動了40fs（理論計算也是40fs）。意味著脈衝在12這個位置放置一個「示波器」（目前不存在這種示波器）所測的時域波形是這樣子的。而我們這個是通過模擬模擬的，而不是計算的。

但是，我們來考慮個問題，我們是希望研究脈衝在傳播過程中的脈衝演變，而這個脈衝在不斷後移，也有可能移出窗口，從後面移出窗口的同時，又會從前面出現，因為第一節課說的，周期性延拓，或者我們叫它周期性邊界條件。在DFT看來，是窗口外周期性存在這個脈衝，那麼當前脈衝l向後退出了，那麼前一個也會進入我們的窗口。

E_f = fft(E_t);z = 12*um;k = 2*pi*ft/c;E_f = E_f.*exp(-i.*k.*z);E_t = ifft(E_f);plot(fftshift(t)/fs,fftshift(E_t));

當然，研究脈衝在傳播過程中的脈衝演變是不希望脈衝移出我們的窗口的。那麼我們就需要使用一個速度去補償這個後移的速度。也就是群速度。

這時候公式變成 。由於真空中， ,所以有 。脈衝形狀不會發生變化。

然而，在其他介質中，例如：N-BK7 (SCHOTT).(折射率公式見Refractive index)。裡面的inv_vg是 。我們可以根據我們計算的群折射率ng跟鏈接中的ng是否一樣，來驗證我們的計算是否正確。

n = @(x)sqrt(1+1.03961212./(1-0.00600069867./x.^2)+0.231792344./(1-0.0200179144./x.^2)+1.01046945./(1-103.560653./x.^2))syms w_tkz = n(2*pi*c/w_t/um)*w_t/c;dk_dwt = diff(kz,w_t);inv_vg = double(subs(dk_dwt,{w_t},{2*pi*f_0}));ng = inv_vg*c;

我們先看看不加上用群速度去補償脈衝移動，可以看到較真空情況下，脈衝來得更晚一些。

E_f = fft(E_t);z = 12*um;n_f = n(c./ft./um);n_f = real(n_f);k = 2*pi*ft/c.*n_f; E_f = E_f.*exp(-i.*k.*z);E_t = ifft(E_f);plot(fftshift(t)/fs,fftshift(E_t));

把模擬公式改成下面所示。可以看到脈衝就回到中間了。但是，好像脈衝沒有發生多大變化，這是因為我們的距離比較小，還不足以產生明顯變化。

E_f = E_f.*exp(-i.*k.*z).*exp(+i.*2*pi*ft.*inv_vg*z);

我們把距離改成5*mm。可以看到脈衝展寬了。

以上就是我說的實數描述，我們看看這種情況的傅里葉變化下的譜分布。可以看到信號主要集中在375THz在附近，其他地方就沒有信號了，為了正確描述信號，那麼採樣率要在375*2THz以上，但是，在兩個正負375THz兩個分布間，有很大的空白地帶，也就是說，很寬的頻帶沒有信息。

plot(fftshift(ft)/THz,abs(fftshift(E_f)).^2);

而複數描述就可以大大減少採樣點數目。在複數描述下， 。我在下面代碼中，我 ， 。選擇的 是800nm對應的頻率。可以看到公式裡面沒有 ，但是後期 是一個重要信息。也就是說，這個超短脈衝的描述由 和 兩個部分組成。

下面代碼可以看出，

我只是把[t,ft] = domain(200*fs,2048);變成了[t,delta_ft] = domain(200*fs,128);。採樣點變成了只需要128個點。然後按照複數描述改寫了E_t，最後加了句ft = f_0+delta_ft;

units;c = 299792458*m/s;[t,delta_ft] = domain(200*fs,128); %%改了duration = 12*fs;f_0 = c/(800*nm);A = exp(-(t./duration).^2);E_t = A; %%改了plot(fftshift(t)/fs,fftshift(E_t));ft = f_0+delta_ft; %%加了

因為是複數描述，所以，計算過程中，E_t也一直是複數描述的超短脈衝，也就是 。因此，繪圖中就需要改成下面的樣子。取模是提取 。光強需要再平方一下就可以。因為是複數，所以包含實部和虛部。當然，我們也可以提取相位。這裡就先不引申了。

plot(fftshift(t)/fs,abs(fftshift(E_t)),black);hold onplot(fftshift(t)/fs,real(fftshift(E_t)),red);plot(fftshift(t)/fs,imag(fftshift(E_t)),blue);hold off

我們來看看複數描述下的效果（顯示的是光強y，也就是 ）。

先是不用群速度補償的(z=12um)。

然後是用群速度補償的(z=12um)。

最後是5mm展寬的。

採用複數描述可以很節約計算內存，計算量少了，速度也就上去了。那為什麼這樣子可以計算正確，又達到這些優點呢。這是因為，前面說過，超短脈衝的場實際上是實數物理量，因此，他的傅里葉變換包含正頻率以及負頻率。而正頻率與負頻率中間有段沒有數據的空間。因此同時，正頻率和負頻率分布是共軛的（強度相同，相位相反 *傅里葉的書上有說）。所以，我們只需要對正頻率繼續計算就可以了。

那麼我們就可以這樣子操作：

變成 。

正頻率分布項是： 。然後我們對這個項目進行頻移，變成： ，並且記下頻移量 。並記 。

所以，前面說，這個超短脈衝的描述由 和 兩個部分組成。

因此，我們構建空間時，實際上要構建了 附近的頻率域，所以我才把ft改成了delta_ft。

[t,delta_ft] = domain(200*fs,128);

但是，我們在計算折射率以及模擬需要使用真實的頻率，因此需要使用ft = delta_ft+f0。使用真是的頻率計算折射率等參數。

所以，最後我們來看看這個複數描述的譜

plot(fftshift(ft)/THz,abs(fftshift(E_f)).^2);

可以看到，的確是不包含負頻率的。