Распространение света в среде
Лекция 1. Уравнение переноса
ОИЯИ
Уравнение переноса света
При упругом рассеянии длина волны не меняется. Рассматриваем монохроматический свет; оптические коэффициенты задаём для выбранной длины волны и дальше её не обозначаем.
Задаём источник \(q\) и оптические свойства среды. Ищем интенсивность \(I(\mathbf x,\hat{\mathbf s},t)\):
\[ \boxed{ \begin{aligned} \frac1{v_g}\frac{\partial I}{\partial t} +\hat{\mathbf s}\cdot\nabla I ={}&-(\mu_a+\mu_s)I\\ &+\mu_s\int_{S^2} p(\hat{\mathbf s},\hat{\mathbf s}') I(\mathbf x,\hat{\mathbf s}',t)\,d\Omega' +q. \end{aligned} } \]
Слева — изменение поля при движении фотонов. Справа — поглощение, рассеяние и создание фотонов.
| Обозначение | Значение | Единица |
|---|---|---|
| \(\mathbf x\) | Положение | м |
| \(\hat{\mathbf s}\) | Направление движения, \(|\hat{\mathbf s}|=1\) | Безразмерная |
| \(t\) | Время | нс |
| \(I\) | Угловая интенсивность числа фотонов | фотон/(м\(^2\)·нс·ср) |
| \(d\Omega\) | Элемент телесного угла на сфере направлений \(S^2\) | ср |
В интеграле \(\hat{\mathbf s}'\) — входящее направление, \(\hat{\mathbf s}\) — направление после рассеяния.
| Обозначение | Значение | Единица |
|---|---|---|
| \(v_g=c_0/n_g\) | Групповая скорость света | м/нс |
| \(\mu_a\) | Коэффициент поглощения | м\(^{-1}\) |
| \(\mu_s\) | Коэффициент рассеяния | м\(^{-1}\) |
| \(p(\hat{\mathbf s},\hat{\mathbf s}')\) | Плотность вероятности исходящего направления при рассеянии | ср\(^{-1}\) |
| \(q\) | Скорость создания фотонов на единицу объёма и телесного угла | фотон/(м\(^3\)·нс·ср) |
\(n_g\) — групповой показатель преломления, \(c_0=0.299792458\) м/нс — скорость света в вакууме.
| Величина | Определение | Единица |
|---|---|---|
| Плотность числа фотонов \(\rho\) | \(\rho=\int_{S^2}f\,d\Omega,\quad f=I/v_g\) | фотон/м\(^3\) |
| Скорость флюенса \(\Phi\) (скалярный поток) | \(\Phi=\int_{S^2}I\,d\Omega=v_g\rho\) | фотон/(м\(^2\)·нс) |
| Вектор потока \(\mathbf J\) | \(\mathbf J=\int_{S^2}\hat{\mathbf s}\,I\,d\Omega\) | фотон/(м\(^2\)·нс) |
| Флюенс \(F\) | \(F=\int\Phi\,dt\) | фотон/м\(^2\) |
\(\rho\) и \(\Phi\) не зависят от направления: направления сложены. \(\mathbf J\cdot\hat{\mathbf n}\) — чистый поток через площадку с нормалью \(\hat{\mathbf n}\) (вперёд минус назад); в изотропном поле \(\mathbf J=0\). Флюенс — число фотонов, попавших на малую сферу, на единицу площади её поперечного сечения.
| Член уравнения | Что он описывает |
|---|---|
| \(v_g^{-1}\partial_t I\) | Изменение числа фотонов в данной точке |
| \(\hat{\mathbf s}\cdot\nabla I\) | Перенос вдоль направления движения |
| \(-\mu_a I\) | Поглощение фотонов |
| \(-\mu_s I\) | Убыль из выбранного направления при рассеянии |
| \(+\mu_s\int p I(\hat{\mathbf s}')\,d\Omega'\) | Приход в выбранное направление после рассеяния |
| \(+q\) | Излучение источника |
Теперь определим нормировку \(I\) и выведем каждый из этих членов.
Что измеряет интенсивность
Монохроматическая интенсивность \(I\) учитывает все фотоны выбранного света. Её нормировку задаём ожидаемым числом фотонов в элементе объёма в момент \(t\):
\[ \boxed{ d\langle N\rangle =\frac{I(\mathbf x,\hat{\mathbf s},t)}{v_g} \,d^3x\,d\Omega. } \]
Угловая плотность числа фотонов равна
\[ f=\frac{I}{v_g},\qquad [f]=\frac{\text{фотон}}{\text{м}^3\,\text{ср}}. \]
\(d^3x\) — пространственный объём, \(\int_{S^2}d\Omega=4\pi\).
Площадка \(dA\) имеет внешнюю нормаль \(\hat{\mathbf n}\). Для \(\hat{\mathbf s}\cdot\hat{\mathbf n}>0\) за время \(dt\) её пересекут фотоны из объёма
\[ dV=v_g\,dt\,dA\,(\hat{\mathbf s}\cdot\hat{\mathbf n}). \]
Подставляем его в определение плотности:
\[ \begin{aligned} d\langle N_{\rm cross}\rangle &=f\,dV\,d\Omega\\ &=I(\hat{\mathbf s}\cdot\hat{\mathbf n}) \,dA\,dt\,d\Omega. \end{aligned} \]
Интенсивность задаёт поток фотонов через площадку.
Условия модели и нормировки
Среда однородна и стационарна. Рассеяние упругое. Поляризацию в этой модели не учитываем.
На коротком участке пути \(d\ell\):
\[ \begin{aligned} dP(\text{поглощение})&=\mu_a\,d\ell,\\ dP(\text{рассеяние})&=\mu_s\,d\ell. \end{aligned} \]
В вероятностях за короткое время используем \(d\ell=v_g\,dt\).
Для фотона, движущегося вдоль \(\hat{\mathbf s}'\), вероятность выйти в угол \(d\Omega\) около \(\hat{\mathbf s}\) равна
\[ dP(\hat{\mathbf s}\mid\hat{\mathbf s}') =p(\hat{\mathbf s},\hat{\mathbf s}')\,d\Omega. \]
Первый аргумент — исходящее направление, второй — входящее.
\[ \boxed{ p\ge0,\qquad \int_{S^2}p(\hat{\mathbf s},\hat{\mathbf s}')\,d\Omega=1. } \]
Нормировка выражает сохранение фотона при рассеянии.
\(d\langle N_{\rm created}\rangle=q\,d^3x\,d\Omega\,dt\). Сравним угловые распределения источника:
\(q=q_0h(\hat{\mathbf s})\), \(\int_{S^2}h\,d\Omega=1\).
\[ h=\delta_{S^2}(\hat{\mathbf s}-\hat{\mathbf s}_0). \]
\[ h=\frac1{4\pi}. \]
\(\hat{\mathbf u}\) — ось сегмента, \(\beta c_0\) — скорость частицы, \(n_{\rm ph}\) — фазовый показатель преломления. Черенковское излучение возникает при \(\beta n_{\rm ph}>1\).
\[ \begin{gathered} h_C=\frac{\delta(\hat{\mathbf s}\cdot\hat{\mathbf u}-\cos\theta_C)}{2\pi},\\ \cos\theta_C=\frac1{\beta n_{\rm ph}}. \end{gathered} \]
Вывод: баланс фотонов
Возьмём неподвижный объём \(V\) и выбранное направление \(\hat{\mathbf s}\) с телесным углом \(d\Omega\). Число фотонов в нём равно \(d\Omega\int_Vf\,d^3x\). За время \(dt\) оно меняется из-за пяти процессов. Вклады даны на единицу \(d\Omega\,dt\):
| Процесс | Вклад |
|---|---|
| Вылет минус влёт через границу | \(-\oint_{\partial V}I\,(\hat{\mathbf s}\cdot\hat{\mathbf n})\,dA\) |
| Поглощение | \(-\int_V\mu_aI\,d^3x\) |
| Уход из направления при рассеянии | \(-\int_V\mu_sI\,d^3x\) |
| Приход в направление при рассеянии | \(+\int_V\mu_s\int_{S^2}p(\hat{\mathbf s},\hat{\mathbf s}')I(\hat{\mathbf s}')\,d\Omega'\,d^3x\) |
| Источник | \(+\int_Vq\,d^3x\) |
Скорость изменения \(\partial_t\int_Vf\,d^3x\) равна сумме этих вкладов. Выведем их по очереди: сначала поток через границу.
Через площадку \(dA\) с внешней нормалью \(\hat{\mathbf n}\) за время \(dt\) проходит \(I(\hat{\mathbf s}\cdot\hat{\mathbf n})\,dA\,dt\,d\Omega\) фотонов (таб «Фотоны через площадку»). Там, где \(\hat{\mathbf s}\cdot\hat{\mathbf n}<0\), фотоны входят, и вклад отрицателен.
Вылет минус влёт — один интеграл по замкнутой границе. Теорема Гаусса при постоянном \(\hat{\mathbf s}\):
\[ \begin{aligned} \oint_{\partial V}I(\hat{\mathbf s}\cdot\hat{\mathbf n})\,dA &=\int_V\nabla\cdot(I\hat{\mathbf s})\,d^3x\\ &=\int_V\hat{\mathbf s}\cdot\nabla I\,d^3x. \end{aligned} \]
Вывод: поглощение и рассеяние
За время \(dt\) поглощается доля фотонов \(dP_a=\mu_a v_g\,dt\). В элементе \(d^3x\,d\Omega\) их \(f\,d^3x\,d\Omega\). Убыль равна
\[ -dP_a\,f\,d^3x\,d\Omega =-\mu_a v_gf\,d^3x\,d\Omega\,dt. \]
Вклад на единицу \(d^3x\,d\Omega\,dt\), с учётом \(v_gf=I\):
\[ \boxed{-\mu_aI}. \]
За время \(dt\) рассеивается доля \(dP_s=\mu_s v_g\,dt\) фотонов. Убыль из выбранного направления \(\hat{\mathbf s}\):
\[ -dP_s\,f\,d^3x\,d\Omega =-\mu_s v_gf\,d^3x\,d\Omega\,dt. \]
Вклад на единицу \(d^3x\,d\Omega\,dt\):
\[ \boxed{-\mu_sI}. \]
Фотон не пропадает: после рассеяния он движется в другом направлении. Его вклад — в приходе.
Из входящего угла \(d\Omega'\) за время \(dt\) рассеивается \(dP_s\,f(\hat{\mathbf s}')\,d^3x\,d\Omega'\) фотонов. В выбранный угол \(d\Omega\) из них попадает доля \(dP(\hat{\mathbf s}\mid\hat{\mathbf s}')=p(\hat{\mathbf s},\hat{\mathbf s}')\,d\Omega\).
Суммируем по \(d\Omega'\) и делим на \(d^3x\,d\Omega\,dt\):
\[ \boxed{\mu_s\int_{S^2} p(\hat{\mathbf s},\hat{\mathbf s}')\, I(\hat{\mathbf s}')\,d\Omega'}. \]
Складываем полученные члены
Подставим вклады в баланс. Так как \(\partial_t\int_Vf\,d^3x=\int_V v_g^{-1}\partial_tI\,d^3x\),
\[ \begin{aligned} \int_V\Bigl[\frac1{v_g}\partial_tI +\hat{\mathbf s}\cdot\nabla I &+(\mu_a+\mu_s)I\\ &-\mu_s\int_{S^2}p(\hat{\mathbf s},\hat{\mathbf s}') I(\hat{\mathbf s}')\,d\Omega'-q\Bigr]d^3x=0. \end{aligned} \]
Объём \(V\) произволен, поэтому равно нулю подынтегральное выражение. Это RTE, с которого начали лекцию.
Подставляем \(I=v_g f\). В однородной стационарной среде:
\[ \boxed{ \begin{aligned} \partial_t f+v_g\hat{\mathbf s}\cdot\nabla f ={}&-v_g(\mu_a+\mu_s)f\\ &+v_g\mu_s\int_{S^2} p(\hat{\mathbf s},\hat{\mathbf s}') f(\hat{\mathbf s}')\,d\Omega'+q. \end{aligned} } \]
Все члены имеют размерность \(\text{фотон}/(\text{м}^3\,\text{нс}\,\text{ср})\).
Проверка сохранения числа фотонов
Интегрируем приход по исходящим направлениям:
\[ \begin{aligned} &\mu_s\int_{S^2}d\Omega\int_{S^2}d\Omega'\, p(\hat{\mathbf s},\hat{\mathbf s}')I(\hat{\mathbf s}')\\ &\quad=\mu_s\int_{S^2}I(\hat{\mathbf s}') \underbrace{\int_{S^2}p(\hat{\mathbf s},\hat{\mathbf s}')\,d\Omega}_{1} \,d\Omega'\\ &\quad=\mu_s\int_{S^2}I(\hat{\mathbf s}')\,d\Omega'. \end{aligned} \]
Этот приход сокращается с убылью \(-\mu_s\int I\,d\Omega\). Рассеяние сохраняет суммарное число фотонов.
Введём плотность числа и плотность потока:
\[ \rho=\int_{S^2}f\,d\Omega,\qquad \mathbf J=\int_{S^2}\hat{\mathbf s}I\,d\Omega. \]
Интегрирование RTE по направлениям даёт
\[ \boxed{ \partial_t\rho+\nabla\cdot\mathbf J =-\mu_a v_g\rho+q_0. } \]
Поглощение уменьшает число фотонов, источник увеличивает его. Перенос меняет число внутри конечного объёма через поток на границе.
Пусть источник мгновенный и действует при \(t=0\). Подставим в RTE \(I=e^{-\alpha t}I_0\), где \(\alpha=v_g\mu_a\):
\[ \frac1{v_g}\partial_t\left(e^{-\alpha t}I_0\right) =e^{-\alpha t}\left(\frac1{v_g}\partial_tI_0-\mu_aI_0\right). \]
Член \(-\mu_aI\) сокращается: \(I_0\) решает RTE без поглощения с тем же начальным условием. Поэтому точно
\[ \boxed{I(\mathbf x,\hat{\mathbf s},t)=e^{-\alpha t}I_0(\mathbf x,\hat{\mathbf s},t)}. \]
Однородное поглощение меняет число фотонов, но не форму поля.
Направленная вспышка: центр облака
В начале координат при \(t=0\) испущено \(N\) фотонов вдоль \(\hat{\mathbf n}\).
Средний косинус угла рассеяния \(g\) задаётся равенством \[ \int_{S^2}\hat{\mathbf s}\,p\,d\Omega=g\,\hat{\mathbf s}'. \]
Определим: \[ \mu_s'=\mu_s(1-g), \qquad\tau_*=1/(v_g\mu_s'),\qquad\ell_*=v_g\tau_*. \]
Усредняем по всему пространству и всем направлениям среди сохранившихся фотонов:
\[ \pi=\frac{f}{Ne^{-\alpha t}},\qquad \langle h\rangle=\int h\pi\,d^3x\,d\Omega, \qquad\int\pi\,d^3x\,d\Omega=1. \]
Однородное поглощение сокращается из этих нормированных средних.
Интегрирование RTE даёт точные уравнения:
\[ \frac{d\langle\hat{\mathbf s}\rangle}{dt} =-\frac{\langle\hat{\mathbf s}\rangle}{\tau_*}, \qquad \frac{d\langle\mathbf x\rangle}{dt}=v_g\langle\hat{\mathbf s}\rangle. \]
При \(\langle\hat{\mathbf s}(0)\rangle=\hat{\mathbf n}\), \(\langle\mathbf x(0)\rangle=0\):
\[ \boxed{\langle\hat{\mathbf s}(t)\rangle=e^{-t/\tau_*}\hat{\mathbf n}}, \qquad \boxed{\mathbf m(t)\equiv\langle\mathbf x(t)\rangle =\ell_*(1-e^{-t/\tau_*})\hat{\mathbf n}}. \]
При \(t\ll\tau_*\) центр движется как \(v_gt\hat{\mathbf n}\); при \(t\gg\tau_*\) его смещение стремится к \(\ell_*\hat{\mathbf n}\), и \(|\mathbf m(t)|<\ell_*\) при любом \(t\).
Центр останавливается, но фотоны продолжают лететь: направление забыто, и облако расширяется диффузией вокруг \(\ell_*\hat{\mathbf n}\) (ширина — на следующем слайде).
Направленная вспышка: второй порядок
Обозначения \(S\), \(C\), \(X\) — матрицы \(3\times3\) с элементами \(S_{ij}\), \(C_{ij}\), \(X_{ij}\). \(\hat{\mathbf n}\hat{\mathbf n}^T\) — матрица с элементами \(n_in_j\); \(\mathbb I\) — единичная матрица. \(X\) — не вектор \(\mathbf m=\langle\mathbf x\rangle\), а матрица \(\langle x_ix_j\rangle\).
Определим глобальные совместные моменты:
\[ S_{ij}=\langle\hat s_i\hat s_j\rangle,\qquad C_{ij}=\langle x_i\hat s_j\rangle,\qquad X_{ij}=\langle x_ix_j\rangle. \]
Нужно знать, как быстро рассеяние стирает угловую зависимость.
Пусть фазовая функция зависит только от \(\hat{\mathbf s}\cdot\hat{\mathbf s}'\). Тогда сферическая гармоника степени \(l\) при рассеянии умножается на момент фазовой функции
\[ g_l=\int_{S^2}P_l(\hat{\mathbf s}\cdot\hat{\mathbf s}')\, p(\hat{\mathbf s},\hat{\mathbf s}')\,d\Omega, \qquad g_1=g, \]
и затухает со скоростью
\[ \gamma_l=v_g\mu_s(1-g_l). \]
Здесь \(P_l\) — полином Лежандра, \(p\) — фазовая функция. Для \(l=1\) получаем \(\gamma_1=1/\tau_*\). Для фазовой функции Хеньи–Гринстейна \(g_l=g^l\), в частности \(g_2=g^2\). Для произвольной фазовой функции \(g_2\) нужно вычислить отдельно; это среднее \(P_2\), а не квадрат среднего косинуса.
Обозначения \(S\), \(C\), \(X\) — матрицы \(3\times3\) с элементами \(S_{ij}\), \(C_{ij}\), \(X_{ij}\). \(\hat{\mathbf n}\hat{\mathbf n}^T\) — матрица с элементами \(n_in_j\); \(\mathbb I\) — единичная матрица. \(X\) — не вектор \(\mathbf m=\langle\mathbf x\rangle\), а матрица \(\langle x_ix_j\rangle\).
Квадрупольная часть углового распределения затухает со скоростью \(\gamma_2\):
\[ \boxed{S(t)=\frac13\,\mathbb I +e^{-\gamma_2t}\left(\hat{\mathbf n}\hat{\mathbf n}^{T}-\frac13\,\mathbb I\right)}. \]
При \(t=0\) все фотоны летят вдоль \(\hat{\mathbf n}\): \(S=\hat{\mathbf n}\hat{\mathbf n}^T\). При \(t\to\infty\) поле изотропно: \(S=\mathbb I/3\).
Обозначения \(S\), \(C\), \(X\) — матрицы \(3\times3\) с элементами \(S_{ij}\), \(C_{ij}\), \(X_{ij}\). \(\hat{\mathbf n}\hat{\mathbf n}^T\) — матрица с элементами \(n_in_j\); \(\mathbb I\) — единичная матрица. \(X\) — не вектор \(\mathbf m=\langle\mathbf x\rangle\), а матрица \(\langle x_ix_j\rangle\).
Вдоль \(\hat{\mathbf n}\) и поперёк (\(\mathbf e_\perp\perp\hat{\mathbf n}\), \(|\mathbf e_\perp|=1\)):
\[ \hat{\mathbf n}^{T}S\,\hat{\mathbf n}=\frac13+\frac23e^{-\gamma_2t}, \qquad \mathbf e_\perp^{T}S\,\mathbf e_\perp=\frac13\left(1-e^{-\gamma_2t}\right). \]
Продольная доля падает от 1 до \(1/3\), поперечная на каждую ось растёт от 0 до \(1/3\). Для HG \(\gamma_2=(1+g)\gamma_1\): квадруполь стирается быстрее, чем направленность.
Обозначения \(S\), \(C\), \(X\) — матрицы \(3\times3\) с элементами \(S_{ij}\), \(C_{ij}\), \(X_{ij}\). \(\hat{\mathbf n}\hat{\mathbf n}^T\) — матрица с элементами \(n_in_j\); \(\mathbb I\) — единичная матрица. \(X\) — не вектор \(\mathbf m=\langle\mathbf x\rangle\), а матрица \(\langle x_ix_j\rangle\).
При движении \(\dot{\mathbf x}=v_g\hat{\mathbf s}\), поэтому
\[ \frac{d}{dt}\langle x_ix_j\rangle =v_g\bigl(\langle\hat s_ix_j\rangle+\langle x_i\hat s_j\rangle\bigr). \]
Рассеяние меняет направление, но не положение в момент столкновения. В \(C_{ij}=\langle x_i\hat s_j\rangle\) направление входит как гармоника степени 1 и затухает со скоростью \(\gamma_1\). Отсюда точно:
\[ \dot C=v_gS-\gamma_1C,\qquad \dot X=v_g(C+C^T),\qquad C(0)=X(0)=0. \]
Первое уравнение решаем интегрирующим множителем, затем интегрируем второе:
\[ C(t)=v_g\int_0^t e^{-\gamma_1(t-u)}S(u)\,du, \qquad X(t)=v_g\int_0^t[C(u)+C(u)^T]\,du. \]
Обозначения \(S\), \(C\), \(X\) — матрицы \(3\times3\) с элементами \(S_{ij}\), \(C_{ij}\), \(X_{ij}\). \(\hat{\mathbf n}\hat{\mathbf n}^T\) — матрица с элементами \(n_in_j\); \(\mathbb I\) — единичная матрица. \(X\) — не вектор \(\mathbf m=\langle\mathbf x\rangle\), а матрица \(\langle x_ix_j\rangle\).
Ширина облака определяется центральным вторым моментом:
\[ \boxed{\Sigma(t)=X(t)-\mathbf m(t)\mathbf m(t)^T}. \]
В силу осевой симметрии:
\[ \Sigma=\sigma_\perp^2(\mathbb I-\hat{\mathbf n}\hat{\mathbf n}^T) +\sigma_\parallel^2\hat{\mathbf n}\hat{\mathbf n}^T. \]
Эти моменты не являются локальным тензором \(\Pi(\mathbf x,t)\) из вывода диффузии (лекция 2).
Те же моменты суммированием по порядкам рассеяния для фазовой функции Хеньи–Гринстейна (§4, формулы (40)–(54)): В. Аллахвердян, Д. В. Наумов, Infinite series solution of the time-dependent radiative transfer equation in anisotropically scattering media, J. Quant. Spectrosc. Radiat. Transf. 326, 109126 (2024), doi:10.1016/j.jqsrt.2024.109126.
Практика 1. Установить LightHit
- LightHit — пакет DN для ожидаемого отклика оптических модулей.
- Нужен Python 3.11 или новее; в курсе фиксируем опубликованную версию 0.3.0a4.
- Метры и наносекунды; скорость получается \(0.218615\) м/нс.
- Первые два аргумента — \(\mu_a,\mu_s\) в м\(^{-1}\), затем \(g\) и \(n_g\).
- Здесь 450 нм, образец BaikalWater; HG \(g=0.9\) задан отдельно.
- Распакуйте примеры и откройте папку
examples.
summary.json: коэффициенты, скорость и оптические длины.moments.csv: выживание и точный центр направленной вспышки.centre.svg: центр облака в транспортных единицах.- Разбор примера · Задачи 1.1–1.6.
