延遲、極限環與相位擴散：光機械懸浮納米粒子的數學框架

想像站在光學平台前，一顆介電納米粒子被聚焦雷射懸浮在真空腔內。壓力進入 Knudsen 區域後，系統高度欠阻尼。實驗者可根據量得位置調制捕獲雷射，以冷卻粒子或操控其運動；但任何電子回饋都會引入有限時間延遲 $\tau$ 。

這就是光機械延遲困境：有限延遲把簡單的二階常微分方程變成需要完整歷史函數的延遲微分方程（DDE）。本文由無因次 DDE 開始，再加入熱雜訊形成 stochastic DDE（SDDE），並以歷史插值數值求解。最後透過 bifurcation diagram、power spectral density（PSD）與 phase diffusion，檢查模型行為是否與所假設機制一致。

模擬相圖

1. 問題與假設

模型研究非線性光學勢及延遲回饋下的一維懸浮粒子，建立於三項假設：

**偶極近似：**粒子半徑 $R$ 遠小於雷射波長 $\lambda$ ，光力可近似成強度梯度。
**Duffing 非線性：**真實光學勢為 Gaussian 形狀；在位移足夠小時，以三次項 Taylor 展開。
**Langevin 熱環境：**殘餘氣體碰撞與其他擾動以阻尼加 Gaussian white noise 表示。

目標是找出延遲引致 Hopf bifurcation 的穩定邊界，檢查非諧捕獲與延遲回饋能否形成穩定極限環，再由譜峰寬度估計相位擴散係數 $D_\phi$ 。

延遲會把目前運動與過去狀態連接，因此系統狀態不再只由 $(z,\dot z)$ 決定，而要包含一段歷史。這也是 DDE 在數學上屬無限維動力系統的原因。

2. 數學模型

2.1 無因次 delay–Duffing 方程

以 $\tilde t=\omega_0t$ 及 $z=x/x_c$ 作尺度化：

\ddot z(t)+\gamma\dot z(t)+z(t)-\alpha z^3(t) =-\kappa z(t-\tau)+\beta z^3(t-\tau).

$\alpha$ 控制保守捕獲軟化， $\kappa$ 與 $\beta$ 分別控制線性及三次延遲位置回饋。此結構把目前位置的保守非線性，與依賴過去位置的非保守能量注入結合。

「不同於標準 forced Duffing」只代表方程結構不同，並不自動代表新穎性；學術主張仍需與現有 delayed Duffing、Pyragas control 及 optomechanical feedback 文獻作系統比較。

2.2 隨機動力學

加入 white noise：

\ddot z(t)+\gamma\dot z(t)+z(t)-\alpha z^3(t) =-\kappa z(t-\tau)+\beta z^3(t-\tau) +\sigma\eta(t).

在物理座標下，平衡 Langevin 系統的阻尼與雜訊需滿足 fluctuation–dissipation relation；例如速度方程的力雜訊強度會包含質量、阻尼、 $k_BT$ 及尺度化因子。無因次式中只寫 $\sigma\propto\sqrt{2\gamma k_BT}$ 並不完整，必須由尺度化推導 $\sigma$ 的確切單位與係數。

2.3 線性穩定與 Hopf 邊界

在 $z^*=0$ 線性化：

\lambda^2+\gamma\lambda+1+\kappa e^{-\lambda\tau}=0.

令 $\lambda=i\omega$ ，分開實部與虛部：

1-\omega^2+\kappa\cos(\omega\tau)=0, \qquad \gamma\omega-\kappa\sin(\omega\tau)=0.

因此候選 $\omega_c$ 必須先滿足

(1-\omega_c^2)^2+\gamma^2\omega_c^2=\kappa^2.

其後才可由相位條件求 $\tau_c$ 。其中一種分支寫法為

\tau_c^{(n)} =\frac{1}{\omega_c} \left[ \operatorname{atan2} \!\left(\gamma\omega_c,\omega_c^2-1\right) +2n\pi \right],

並需檢查符號約定、分支與 crossing direction。原文只使用 arccos 的公式會遺漏正弦條件，不能單獨決定所有臨界延遲。

2.4 相位擴散

極限環的徑向擾動通常被吸引回軌道，沿相位方向則沒有固定參考點。弱雜訊可令相位近似 random walk：

\left\langle[\Delta\phi(t)]^2\right\rangle =2D_\phi t.

若訊號近似 $z(t)=A_0\cos(\omega_{LC}t+\phi(t))$ ，而 phase increments 為 Gaussian diffusion，相關函數按 $e^{-D_\phi|t|}$ 衰減，頻譜基本峰附近形成 Lorentzian：

S(\omega)\propto \frac{D_\phi} {D_\phi^2+(\omega-\omega_{LC})^2}.

這是弱雜訊、穩定單一極限環及純相位擴散下的近似；幅度—相位耦合、colored noise 或 switching 可產生非 Lorentzian 峰。

3. 數值方法

3.1 歷史與延遲插值

DDE 求解器要保存過去狀態。當 $\tau/\Delta t$ 不是整數，原實作以線性插值取 $z(t-\tau)$ ：

def get_z_tau(step_idx, offset_fraction=0.0):
    float_idx = (
        step_idx + offset_fraction - delay_steps_float
    )
    idx0 = int(np.floor(float_idx))
    idx1 = idx0 + 1
    weight = float_idx - idx0
    return (
        z[idx0] * (1.0 - weight)
        + z[idx1] * weight
    )

在 deterministic drift 使用 RK4、在每步 stochastic increment 使用 Euler–Maruyama，可稱為 split 或混合離散法，不能僅因 deterministic 部分是 RK4 便稱為四階強收斂 stochastic Runge–Kutta。整體 stochastic strong／weak order 需要另行分析；延遲插值、歷史初始化及步長收斂亦要測試。

3.2 譜估計

高解像模擬使用長時間序列及 Welch 方法，nperseg=65536 提高頻率網格解像度。較長 segment 會減少平均次數、增加估計方差，因此應同時報告 sampling interval、window、overlap、transient removal 及多個獨立軌跡。

3.3 Lorentzian 擬合

def lorentzian(
    w, amplitude, w0, D_phi, noise_floor
):
    return (
        amplitude * D_phi
        / ((w - w0)**2 + D_phi**2)
        + noise_floor
    )

curve_fit 可以估計 $D_\phi$ ，但相鄰 Welch bins 並非獨立同方差 Gaussian observations。擬合區間、背景模型、window broadening 與有限時間解像都會影響線寬；應以模擬 calibration 或 bootstrap 量化偏差。

3.4 評估流程

模型進行兩項一致性檢查：

**小位移近似：**比較 Gaussian 恢復力 $F_{\mathrm{true}}(z)=-z e^{-\alpha z^2}$ 與 Duffing 展開 $F_{\mathrm{Duffing}}(z)=-z+\alpha z^3$ ，檢查極限環是否仍在展開誤差小的範圍。
**零回饋平衡檢查：**令 $\kappa=\beta=0$ ，比較模擬速度方差與指定離散 Langevin 模型的理論值 $\langle v^2\rangle=\frac{\sigma^2}{2\gamma}.$

兩項檢查可以發現明顯實作錯誤，卻不能「證明」整個延遲模型物理有效。第二項只驗證零回饋、指定無因次 SDE 及數值步長下的平衡統計。

4. 實驗協議

評估三類量：

位移 $z$ ：極限環振幅及不同延遲下的局部最大值；
PSD：不同無因次頻率的訊號功率；
$D_\phi$ ：以 Lorentzian 擬合估計的 phase drift。

流程為：

在無雜訊下掃描 $\tau\in[0.1,5.0]$ ，繪製分岔結構；
在選定點 $\tau=3.5$ 、 $\sigma=0.08$ 運行長 stochastic 模擬；
對多個 $\sigma^2$ 估計 $D_\phi$ ，檢查線性尺度關係。

更嚴謹協議應在掃描前固定 transient 長度、初始歷史、每點多重初值、步長與主分岔指標，避免把有限模擬造成的振幅消失誤判為分岔。

5. 結果解讀

5.1 近似與熱平衡檢查

Duffing 與 Gaussian 力及極限環範圍

零回饋下的速度直方圖

報告極限環最大振幅約 $z=0.12$ ，位於兩種恢復力相近區域。模擬速度標準差 $0.153$ ，理論值 $0.158$ ，相差約 3%。這支持小位移與零回饋離散實作的一致性，但「3%」是否可接受仍需與 Monte Carlo uncertainty 及步長誤差比較。

5.2 相圖

τ=3.5 的位置—速度相圖

在 $\tau=3.5$ 、 $\kappa=0.3$ 的指定軌跡中，投影到 $(z,\dot z)$ 平面形成閉合曲線。要證明它是穩定 isolated limit cycle，還應由多個歷史函數出發、檢查 transverse convergence 及 Floquet multipliers；單一漂亮閉合軌跡只是支持證據。

5.3 分岔圖

局部位移最大值對延遲的分岔圖

圖中平衡在 $\tau\approx0.25$ 附近失去穩定，之後出現有限振幅運動，振幅隨延遲變化， $\tau>3.5$ 附近再下降。是否屬次級分岔、振幅 collapse、數值 transient 或 attractor switching，需以解析邊界、延長模擬及 continuation 區分。

5.4 PSD 與 Lorentzian

stochastic PSD

基本峰 Lorentzian 擬合

基本峰報告位於 $\omega_0\approx0.147$ ，擬合得到 $D_\phi\approx0.0078$ 。擬合外觀與 Lorentzian 相容，但仍需殘差、替代線型及跨 segment／seed 的區間，才能判斷模型是否足夠。

5.5 相位擴散尺度

Dφ 對 σ² 的尺度關係

在測試範圍內， $D_\phi$ 與 $\sigma^2$ 的線性回歸報告 $R^2>0.99$ 。這與弱白雜訊相位擴散預期一致，卻不是物理驗證：兩個量都由同一模擬與擬合流程產生，高 $R^2$ 亦不檢查截距、異方差、估計誤差或範圍外行為。

6. 諧波與線寬

三次非線性可產生奇次諧波 $3\omega_0,5\omega_0,\ldots$ 。若第 $n$ 次諧波相位為 $n\phi(t)$ ，純 phase diffusion 令其相位方差放大 $n^2$ ，理想 Lorentzian 的全寬半高近似

\Delta\omega_n =n^2\Delta\omega_1 =2n^2D_\phi.

觀察二次尺度可支持共用相位噪音機制，但有限解像、幅度雜訊、非線性 frequency pulling 及不同諧波訊噪比都會破壞簡單關係。

7. 限制

模型只處理一維；真實粒子三維運動與軸間耦合可產生額外分岔。
Duffing 展開只在小位移適用；真實光阱深度有限，粒子可逃逸，而截斷多項式可能錯誤地繼續束縛。
white noise 省略雷射 $1/f$ 雜訊、colored gas noise、量度雜訊與回饋電子噪音。
延遲被視為精確常數，沒有 jitter、濾波器相位或數碼取樣。
所有圖來自同一合成模型，未與實驗軌跡比較。

結語

delay–Duffing SDDE 把非線性動力學、統計力學與歷史依賴連在一起。模擬展示一條由解析 Hopf 候選、deterministic 極限環、stochastic 譜峰走到 phase diffusion 的可檢查證據鏈。

目前結果與小位移、零回饋熱平衡及弱雜訊相位擴散的預期相容，不能證明裝置或模型已獲物理驗證。下一步應使用經步長收斂的 DDE／SDDE 求解器、多初始歷史 continuation、重複 stochastic trajectories，以及獨立實驗 PSD 來檢驗 $D_\phi$ 與諧波線寬。