參考資料：Notes on Perfectly Matched Layers (PMLs)

波動方程

以無源標量波動方程為例 nabla cdot(a nabla u)=frac{1}{b} frac{partial^2 u}{partial t^2}=frac{ddot{u}}{b}tag{1} 其中 u(mathbf{x},t) 是標量波的幅度； c=sqrt{ ab } 是波關於介質（可能不均勻）參數 a(mathbf{x}) 和 b(mathbf{x}) 的相速度；對於無損的 propagating waves， a 和 b 應該是實數和正的。 為瞭計算方便（使用 staggered-grid leapfrog discretization）和分析，通過引入輔助場 mathbf{v}(mathbf{x},t) 將二階方程方程拆分為兩個耦合的一階方程 begin{gathered} frac{partial u}{partial t}=b nabla cdot mathbf{v} \ frac{partial mathbf{v}}{partial t}=a nabla u end{gathered}tag{2} 上述方程可以更概括性的寫為 frac{partial mathbf{w}}{partial t}=frac{partial}{partial t}left(begin{array}{c}u \ mathbf{v}end{array}right)=left(begin{array}{cc} & b nabla cdot \ a nabla & end{array}right)left(begin{array}{l} u \ mathbf{v} end{array}right)=hat{D} mathbf{w}tag{3} 其中 hat{D} 是一個 4times4 的線性算子； mathbf{w}=(u;mathbf{v}) 是一個 4 分量矢量（ mathbf{v} 是在三維空間）。隻要選擇瞭恰當內積的 anti-Hermitian 算子 hat{D} ，所有常見的波動方程就都可以被寫成式 (3) 的形式。

復數坐標拉伸

下圖是一個典型的波動方程問題的示意圖，我們希望通過吸收輻射波的方式來截斷無限大的空間，隻保留感興趣（Region of Interest, ROI）的計算區域。

image.png

這種截斷方法有三個概念性的步驟，總結如下： 1. 在無限空間中，將解和方程解析延拓（analytically continue）到一個復數 x 輪廓，這將使振蕩波在感興趣區域之外變成指數衰減的波，而且沒有反射。 2. 在無限空間中，進行坐標變換（coordinate transformation），將復數 x 表示為實數坐標的函數。在新的坐標系中，我們有 real coordinates 和 complex materials。3. 將這個新實數坐標的定義域在復數材料區域內截斷：由於解在那裡衰減，隻要我們在足夠遠的距離之後截斷它（指數尾巴很小的地方），使用哪種邊界條件都無關緊要（硬截斷也是可以的）。

為瞭便於簡化，假設遠離 ROI 的空間是： – 均勻 – 線性時不變

在這些假設下，無限區域的輻射解必須采用平面波（planewaves）疊加的形式： mathbf{w}(mathbf{x},t)=sum_{mathbf{k},omega}mathbf{W}_{mathbf{k},omega} e^{i(mathbf{k}cdot mathbf{x}-omega t)}tag{4} 其中 mathbf{W}_{mathbf{k},omega} 是恒定的幅度； omega 是（角）頻率； mathbf{k} 是波矢量。 輻射解可以分解成下面函數的形式 begin{equation} mathbf{W}(y,z)e^{i(kx-omega t)} tag{5} end{equation} 其中 omega /k 是相速度，它不同於群速度 domega /dk (無損介質中的能量傳輸速度)。對於沿著 +x 方向傳播的波， k 是正的。

解析延拓

因為式 (5) 是 x 的解析函數，所以我們可以自由地通過解析延拓，在復數值的 x 上進行求解。 例如：下圖的 Top 部分為原始波動方程對應實數軸 x 方向上的解 e^{ikx} 。Bottom 部分則將輪廓從實數軸 x 延拓到瞭復數（在 mathrm{Re}:x>5 時，添加一個線性增長的虛部），此時解的形式為 e^{ik (mathrm{Re}:x+i:mathrm{Im}:x)}=e^{ik:mathrm{Re}: x}e^{-k:mathrm{Im}:x} ，在 k>0 的時候將呈現指數衰減。 不同於吸收材料，解在未變形的區域保持不變，所以變形後在復平面上的輪廓區域就對應於無反射的吸收區域，即 PML。

image.png

註意，解析延拓的前提是微分方程在這個區域 x-invariant，所以 x 隻出現在像 partial / partial x 這樣的導數中，且解析函數的導數在復平面的任何 dx 方向上都是相同的。

坐標變換回實數$x$

將復數 x 輪廓表示為 widetilde{x} ，實部為 x ，因此滿足 widetilde x(x)=x+if(x) tag{6} 其中 f(x) 是如何沿著虛軸改變輪廓的函數。因為復數坐標 widetilde x 不方便，所以隻需要改變變量，以便用實部 x 表示方程。所以我們有 partial widetilde x=left ( 1+i frac{df}{dx} right)partial x 。方便起見，記 frac{df}{dx}=frac{sigma_{x}(x)}{omega} 。例如在解析延拓部分的示例中，Bottom 處 sigma_{x}(x) 是一個階梯函數。 因為我們假設原始的波動方程是 x-invariant 的 (至少在 fneq 0 的 large-x 區域），所以我們不需要做其他替換。 PML 的全部處理過程可以在概念上總結為對原始微分方程的單一變換： frac{partial}{partial x}to frac{partial}{partial widetilde x}=frac{1}{1+i frac{sigma_{x}(x)}{omega}} frac{partial}{partial x} tag{7} 之所以選擇 frac{df}{dx}=frac{sigma_{x}(x)}{omega} 而不是 sigma_{x}(x) ，是因為在新坐標系下解如下所示： e^{ikx}to e^{ikwidetilde x}=e^{i k x} e^{-frac{k}{omega} int^x sigma_xleft(x^{prime}right) d x^{prime}}tag{8} 因子 k/omega=1/c_{x} ，其中 c_{x} 為 x 方向的相速度。其中在無色散材料中（例如：在真空中傳播的光），對於固定角度， c_{x} 是與速度無關的常數。 所以在這種情況下，PML 中的衰減率與頻率 omega 無關，即所有波長以同樣的速率衰減。如果忽略 1 / omega ，短波長將衰減的比長波長更快。此外，衰減率並不獨立於光的角度，這個困難將留到後面討論。

計算區域的截斷

對波動方程進行瞭 PML 轉換，如式 (7) 所示，解在我們感興趣的區域（small-x）內保持不變，在外部區域（large-x）則呈指數衰減，這意味著我們可以在某個足夠大的 x 處截斷計算區域，例如放一個 Dirichlet 邊界條件。 原則上來說，我們可以通過讓 sigma_{x} 足夠大（指數衰減速率快），來使得 PML 區域足夠薄。然而在離散化的時候，一個非常大的 sigma_{x} 會造成數值反射；相反，在長度約為半個波長的區域內將 sigma_{x}(x) 二次或三次地從零逐漸增加，實際上反射會非常小。

其他方向上的 PML 邊界

前述理論基於 +x 方向，對於其他方向，如 -x 方向，仍然是在負 large-x 區域處應用 sigma_{x}>0 的 PML 變換，如式 (7) 所示。因為對於 x<0 的區域，沿著 -x 方向傳播的輻射波的波矢量 k<0 ，所以 PML 的解對於相同的正 sigma_{x} 來說，在 xto- infty 上指數衰減。在 pm y,pm z 方向上同理即可，此外在計算單元的角落，將同時具有沿著兩個或三個方向的 PML 區域（即兩個或者三個 sigma 非零）。

坐標變換和材料

在標量波動方程的背景下，PML 坐標變換的 1+isigma /omega 項在波動方程中表現為一種人工各向異性吸收材料（將 a 和 b 轉變為復數，其中 a 是一個張量），具體實例見後文。 以麥克斯韋方程組為例，麥克斯韋方程在任何坐標變換下都可以表達為在笛卡爾坐標系中具有變換材料的麥克斯韋方程。即坐標變換被“吸收”到介電常數 varepsilon 和磁導率 mu 的變化中（通常是各向異性張量）。這就是為什麼構造無反射各向異性吸收器的 UPML 等效於復坐標拉伸的原因：它隻是將坐標拉伸吸收到材料張量中。

頻域和時域中的 PML 實例

在頻域中，當求解具有時間依賴性 e^{-iomega t} 的解時，我們隻需要應用式 (7) 即可。 在時域中，為瞭使衰減速率頻率無關，復數拉伸因子 1+isigma /omega 是頻率相關的。而時域中沒有 omega (例如，時域波函數可以同時疊加多個頻率)，所以我們通過輔助微分方程（auxiliary differential equation, ADE）的方法在時域中直接實現理想的 1/omega 依賴關系

一個 2d 標量波的 PML 實例

對於 2d 標量波方程，如式 (2) 所示： begin{aligned} frac{partial u}{partial t}=b nabla cdot v & =b frac{partial v_x}{partial x}+b frac{partial v_y}{partial y}=-i omega u \ frac{partial v_x}{partial t} & =a frac{partial u}{partial x}=-i omega v_x \ frac{partial v_y}{partial t} & =a frac{partial u}{partial y}=-i omega v_y . end{aligned}tag{9} 再對前兩個方程中的 partial / partial x 執行式 (7) 的 PML 變換，可得 begin{gathered} b frac{partial v_x}{partial x}+b frac{partial v_y}{partial y}left(1+i frac{sigma_x}{omega}right)=-i omega u+sigma_x u \ a frac{partial u}{partial x}=-i omega v_x+sigma_x v_x . end{gathered}tag{10} 第二個方程可以很容易轉為時域，但第一個方程中的 frac{i b sigma_x}{omega} frac{partial v_y}{partial y} 項中含有顯式的 1/omega 因子。在傅裡葉變換中， -i omega 對應對應微分， i / omega 對應積分，所以引入一個新的輔助場變量 psi 滿足 -i omega psi=b sigma_x frac{partial v_y}{partial y}, 所以式 (10) 的第一個方程為 b frac{partial v_x}{partial x}+b frac{partial v_y}{partial y}+psi=-i omega u+sigma_x u . 現在再將所有內容傅裡葉變換回時域，可以得到一組帶有 PML 吸收邊界的沿x方向的四個時域方程

begin{gathered} frac{partial u}{partial t}=b nabla cdot v-sigma_x u+psi \ frac{partial v_x}{partial t}=a frac{partial u}{partial x}-sigma_x v_x \ frac{partial v_y}{partial t}=a frac{partial u}{partial y} \ frac{partial psi}{partial t}=b sigma_x frac{partial v_y}{partial y}, end{gathered}

其中最後一個關於 psi 的方程是我們的輔助方程（初始條件 psi=0 ）。可以註意到在 u 和 v_x 的方程中有 sigma_x 吸收項，而 v_y 沒有：PML 對應於對應於一個各向異性吸收體，就好像將 a 替換為 2 times 2 復張量： left(begin{array}{cc} frac{1}{a}+i frac{sigma_x}{omega a} & \ & frac{1}{a} end{array}right)^{-1}

不均勻介質中的 PML

上述推導並不依賴假設：對於 x 的 PML 層，介質在 (y,z) 平面上需要是均勻的。我們隻需要假設介質在足夠大的 x 方向上是不變的即可。 要點是：隻要解/方程在 x 方向仍然是解析的，PML 就是無反射的。

漸逝波的 PML

上述討論中假設瞭波矢量 k 是正實數，但這並非必要。在二維或者更高維度上，波動方程可以有衰減解（evanescent solutions），其中 k 是復數。 例如，考慮在相速度為 c 的的均勻二維介質中的平面波 e^{i(mathbf{k} cdot mathbf{x}-omega t)} ，即 omega=c|mathbf{k}|=c sqrt{k_x^2+k_y^2} 。在這種情況下： k=k_x=sqrt{frac{omega^2}{c^2}-k_y^2} . 對於足夠大的 k_y （例如 y 方向上的高頻傅裡葉成分）， k 是純虛數，所以當 x rightarrow infty 時，隻有 operatorname{Im} k>0 才能使 e^{i k x} 指數衰減。 在 PML 的媒質中這樣一個衰減的、虛數 k 的漸逝波（evanescent wave）會發生什麼？令 k=i kappa ，在 PML 中：也就是說，PML 在漸逝波中增加瞭振蕩，但並沒有增加衰減率，PML 無反射，但沒有任何幫助。如果我們想加快衰減，隻需要令 mathrm{Im}:sigma_{x}<0 （對應於實坐標拉伸），就能加快漸逝波的衰減。但是增添 sigma_{x} 的虛部也有代價，即讓傳播（ k 是實數）的波振蕩的更快，這加劇瞭數值反射。 簡而言之，一切要適度。

PML 的限制

離散化和數值反射

隻有在求解精確的波動方程時，PML 才是無反射的。一旦對問題進行離散化處理（無論是有限差分還是有限元），就變成瞭一個近似問題，PML 的分析完美性失效瞭。雖然 PML 依舊是一個吸收材料，但是並非無反射，而是會反射很小一部分回來。 怎樣才能讓反射足夠小？一個關鍵事實是：根據絕熱定理（adiabatic theorem），即使沒有 PML，隻要介質變化緩慢，反射可以無限小。 使用非 PML 吸收材料可能需要非常緩慢的進行（非常厚的吸收層）才能得到可接受的反射，然而使用 PML 時，常數因子是一個非常好的開始。經驗表明通過簡單的二次或者三次的 PML 吸收開啟可以得到可忽視的反射（在半個波長或者更薄的 PML 層內）。

角度相關的吸收

PML 的吸收與角度有關。以式 (8) 為例，衰減率是正比於 k /omega ，但是 k 在這裡隻是 k_{x} ，即波矢量 mathbf{k} 在 x 方向上的成分，所以衰減率實際上是正比於 |mathbf{k}|cos theta ，其中 theta 是輻射波與 x 軸的夾角。 當輻射波掠入射（glancing incidence）時，衰減率趨於 0，對於任何固定的 PML 厚度，波都會在 PML 中發生大量“往返”反射。但在大多數情況下，輻射起源於 ROI 的中心，傳播到 PML 的輻射將限制在 theta<55degree approx cos^{-1}(1 / sqrt{ 3 }) 。所以，如果邊界足夠遠，就能保證一個最大角度，並使 PML 足夠厚以便吸收所有角度錐內的波。

令 PML 失效的不均勻介質

PML 在介質不具有 x 不變性（對於一個 x 的邊界）的時候完全失效。 PML 的中心思想是波動方程和解在與邊界垂直的方向上是解析函數，從而能夠解析延拓到復數坐標平面。而介質如果變化則很可能是不連續變化的，從而使得解析延拓的概念不適用瞭。 這種情況下，傳統的吸收邊界條件也失效瞭，隻能采用上面提到的絕熱定理實現無反射。