ЖЭТФ, 2023, том 164, вып. 3 (9), стр. 340-348

ПИК КОГЕРЕНТНОГО ОБРАТНОГО РАССЕЯНИЯ ДЛЯ

ИЗЛУЧЕНИЯ С ПОНИЖЕННОЙ ПРОСТРАНСТВЕННОЙ

КОГЕРЕНТНОСТЬЮ

В. Л. Кузьмин^a*, А. Ю. Вальков^a,b**, Ю. А. Жаворонков^b***

^a Санкт-Петербургский политехнический университет им. Петра Великого

195251, Санкт-Петербург, Россия

^b Санкт-Петербургский государственный университет

198504, Санкт-Петербург, Россия

Поступила в редакцию 25 февраля 2023 г.,

после переработки 3 апреля 2023 г.

Принята к публикации 3 апреля 2023 г.

Развита теория и выполнено численное моделирование эффекта когерентного обратного рассеяния в

сильнонеоднородной случайной среде с конечной длиной пространственной когерентности. В рамках

метода Монте-Карло показано, что ограничение числа актов рассеяния соответствует понижению коге-

рентности падающего излучения и приводит к угловому уширению пика обратного рассеяния, расши-

ряя возможности использования когерентного обратного рассеяния для биомедицинских приложений.

Впервые на основе диаграммного метода развито моделирование когерентного обратного рассеяния за

рамками лестничного приближения.

DOI: 10.31857/S0044451023090031

сравнению с такими приложениями, как NIRS и

EDN: KCCFMX

DCS, в связи с малостью ширины углового конуса

)^-1 [6-8],

КОР, определяемого величиной W₀ ∼ (kl_tr

1. ВВЕДЕНИЕ

где k волновое число, l_tr транспортная длина.

Открытие когерентных

[1-8] и корреляцион-

Известно, что конечность длины когерентности

ных

[9, 10] эффектов многократного рассеяния

объясняет [20, 21], почему отношение когерентной

излучения в случайных средах породило ши-

и некогерентной компонент обратного рассеяния

рокую область биомедицинских приложений в

не достигает теоретически предсказанного удвоения

инфракрасной и ближней инфракрасной областях,

при рассеянии строго назад. Аналогично конечность

таких как ближняя инфракрасная спектроскопия

длины когерентности используется для объяснения

(near infrared spectroscopy, NIRS) и диффузная

того, почему не выполняется соотношение Зигерта

корреляционная спектроскопия (diffuse correlation

в работах по изучению и применению временной

spectroscopy, DCS) [11-14]. Такие эффекты, как

корреляционной функции интенсивности [18, 20]. В

волны фотонной плотности [15, 16] и временные

работах [22-27] была разработана новая методика,

корреляции

[17] интенсивности, легли в основу

использующая излучение с низкой когерентностью,

инновационных методов неинвазивной диагностики

создающее более широкий угловой конус КОР, на-

биотканей, в том числе в реальной биомедицинской

блюдаемый в реальных биомедицинских условиях.

практике [13, 18, 19]. Однако применение эффекта

В работах [24, 25, 28] был рассмотрен и использован

усиления когерентного обратного рассеяния (КОР),

эффект конечности понижения пространственной

наиболее наглядно проявляющего волновую при-

когерентности падающего излучения при наблюде-

роду света в режиме многократного рассеяния, в

нии КОР. Возникающее при этом угловое уширение

биомедицинских приложениях почти не находит

конуса КОР использовалось в реальной медицин-

применения, и его развитие сильно замедлилось по

ской диагностике [25, 28].

В данной работе методом Монте-Карло выпол-

^* E-mail: kuzmin_vl@mail.ru

^** E-mail: alexvalk@mail.ru

нено численное моделирование эффекта КОР для

^*** E-mail: zhavoronkov95@gmail.com

излучения с конечной длиной когерентности. Впер-

340

ЖЭТФ, том 164, вып. 3 (9), 2023

Когерентное обратное рассеяние с пониженной пространственной ...

вые разработана теория эффектов обратного рассе-

щей среде, δε флуктуации диэлектрической про-

яния за пределами лестничного приближения с яв-

ницаемости, а произведения обозначают интеграль-

ным учетом конечности длины когерентности. Неко-

ные свертки. В (r, ω)-представлении

герентная составляющая обратного рассеяния фор-

T₀(r) = k²⁰ exp(ik₀r)/4πr,

мируется лестничными диаграммами, каждая из ко-

торых состоит из произведений двух комплексно-

k₀ =

√ε₀ω/c, ε₀ средняя диэлектрическая прони-

сопряженных полей, рассеивающихся на одной и

цаемость, c скорость света в вакууме.

той же последовательности случайных неоднород-

Усредняя по флуктуациям диэлектрической

ностей. Ряд лестничных диаграмм обеспечивает тео-

проницаемости, получаем уравнения Дайсона для

ретическую основу для диффузионного приближе-

усредненного макроскопического поля E = 〈E_rand 〉

ния. Теория когерентного обратного рассеяния ос-

и макроскопической функции Грина T = 〈T_rand 〉:

нована на суммировании максимальных перекрест-

ных (или циклических) диаграмм, топологически

E = E₀ + T₀ΠE, T = T₀ + T₀ΠT,

(2)

эквивалентных лестничным диаграммам. Низкоко-

герентное обратное рассеяние требует дополнитель-

где ядро Π

= Π(r) известно как оператор по-

ного учета интерференционных вкладов, происходя-

ляризации (ядро массового оператора/компактный

щих от корреляций диэлектрической проницаемости

блок) [30,31] и представляется в виде суммы непри-

высоких порядков рассеяния, описание которых в

водимых диаграмм. Мнимая часть его фурье-обра-

терминах диаграмм требует выхода за рамки лест-

за Π(k) определяет длину рассеяния l_s, что фак-

ничных или максимально-перекрестных диаграмм.

тически соответствует оптической теореме, связы-

Мы впервые рассматриваем такие диаграммы для

вающей полное сечение рассеяния с длиной рассея-

трехкратного рассеяния, учитывая интерференцию

ния [32]. В дальней зоне для |Π(k₀)| ≪ 1 уравнение

между парами полей, распространяющимися в слу-

Дайсона дает

чайной среде. В настоящей работе, принимая макси-

мальное число событий рассеяния n_sc в качестве па-

T (r) ≃ k²⁰ exp(ikr)/4πr,

раметра, определяющего пространственную длину

где k ≃ k₀(1 + Π(k₀)/2). Таким образом, для длины

когерентности L_c = n_sc l_s, где l_s длина рассеяния,

рассеяния получаем

мы показали, что использование низкокогерентного

излучения позволяет получить конус КОР с шири-

l_s = (Im(k₀Π(k₀)))^-1.

ной и относительной высотой, близкими к экспери-

ментальным значениям из работы [6].

Итерируя выражения (1), можно представить

В работах [24,25] высказано предположение, что

корреляционную функцию поля в виде операторно-

при понижении когерентности двукратный вклад

го ряда:

может стать доминирующим. Выполненное нами

(

)

моделирование указывает, что такое доминирование

〈δE δE^∗〉 =

Tδε + (Tδε)² + . . .

E ×

(

)

низких кратностей реализуется только при очень

T^∗δε^∗ + (T^∗δε^∗)² + . . .

E^∗ ,

(3)

больших углах обратного рассеяния, за пределами

обычно наблюдаемого конуса КОР [29]. Мы обнару-

где угловые скобки означают статистическое усред-

жили, что при больших углах обратного рассеяния

нение по случайным неоднородностям диэлектриче-

вклад членов высоких порядков приводит к интер-

ской проницаемости, δE флуктуация поля.

ференционному уменьшению обратного рассеяния.

Парный коррелятор поля может быть представ-

лен в виде ряда по порядкам рассеяния

∫

2. ДИАГРАММНОЕ ПРЕДСТАВЛЕНИЕ

∑

ПОЛЯ И ИНТЕНСИВНОСТИ

〈δE(r^′0) δE^∗(r^′0)〉 =

dr^′1dr^′′1. . . dr^′n dr^′′n ×

n=1

Пусть волновые уравнения для случайного поля

∏

E_rand и его функции Грина T_rand схематично пред-

× T(r^′0,r^′1)T^∗(r^′′0,r^′′1)

T (r^′i-1 - r^′i) T^∗(r^′′i-1 - r^′′i) ×

ставлены в виде

i=2

× K(r^′n,r^′′n,...,r^′1,r^′′1)E(r^′n)E^∗(r^′′n),

(4)

E_rand = E₀ + T₀ δε E_rand ,

(1)

T_rand = T₀ + T₀ δε T_rand,

где ядро K(r_n, r^′n, . . . , r₁, r^′1) представляет собой

где нижний индекс ¾0¿ относится к нерассеиваю-

многочастичный коррелятор флуктуаций δε.

341

В. Л. Кузьмин, А. Ю. Вальков, Ю. А. Жаворонков

ЖЭТФ, том 164, вып. 3 (9), 2023

r₂

r₁

r₂

r₁

r'₁

r'₂

r'₁

Рис. 1. Топологическая эквивалентность лестничной D^(2)L и перекрестной D^(2)C диаграмм второго порядка рассеяния

В предположении гауссова характера флуктуаций,

переходит в другую. На этом основано явление ко-

согласно теореме Вика, ядро K представляется в ви-

герентного обратного рассеяния, теоретически пред-

де суммы произведений парных корреляторов,

сказанного [1-4] и экспериментально подтвержден-

ного [5-8] в середине восьмидесятых годов прошло-

∑

∏

K(r^′n, r^′′n, . . . , r^′1, r^′′1) =

〈δε(r^′i) δε^∗(r^′′ )〉,

(5)

го века. Для обратного рассеяния строго назад мак-

j1,..., jn i=1

симально-перекрестные диаграммы равны лестнич-

ным. Обозначая явно зависимость от волновых век-

по всем перестановкам j₁, . . . , j_n чисел 1, 2, . . . , n.

торов входящей и выходящей волн, получаем

Разложение интенсивности (4) представляет со-

бой ряд по кратностям рассеяния, полученный ите-

D^(2)C(k_out, k^∗out |k_in, k^∗in) = D^(2)L(k_out, -k^∗in|k_in, -k^∗out).

рированием уравнения Бете - Солпитера. Основной

вклад для почти всех углов вносят слагаемые лест-

Четыре оставшиеся диаграммы в (7) также разби-

ничного приближения, когда пара комплексно-со-

ваются на две топологически эквивалентные пары,

пряженных полей проходит одну и ту же последо-

как показано на среднем и нижнем изображениях на

вательность флуктуаций, или неоднородностей ди-

рис. 2. Эти диаграммы сильно зависят от углов вхо-

электрической проницаемости, сохраняя при этом

да и выхода, и поэтому в широких угловых пределах

исходную разность фаз. Однако при рассеянии стро-

они могут оказаться важными. Мы будем называть

го назад когерентность сохраняется также для слу-

их случайно-интерференционными диаграммами.

чая, когда поля проходят неоднородности δε в об-

ратном порядке. В этом случае наблюдается эффект

усиления обратного рассеяния.

3. ДИАГРАММЫ 2-ГО И 3-ГО ПОРЯДКОВ

Приведем диаграммы, описывающие двукрат-

ные и трехкратные вклады рассеяния. На рис. 1 по-

При вычислении диаграмм проводится интегри-

казаны диаграммы двукратного рассеяния,

рование по пространственным координатам r^′i и r^′′i,

связывающим парную корреляцию. Удобно перей-

D⁽²⁾ = D^(2)L + D^(2)C,

(6)

ти к переменным интегрирования по координатам

а на рис. 2 трехкратного рассеяния,

¾центра масс¿ r_i = (r^′i +r^′′i)/2 и разностным коорди-

. Диапазон координат центра масс

натам r^′′′i = r^′i -r^′′i

D⁽³⁾ = D^(3)L +D^(3)C +D^(3)LC +D^(3)CL +D^(3)CLL+D^(3)LLC. (7)

r_i определяет длина рассеяния l_s, а масштаб рассто-

яния r^′′′i длина корреляций диэлектрической про-

В (7) мы опустили быстроосциллирующие диаграм-

ницаемости r_c ≪ l_s.

мы, в которых поля и сопряженные поля имеют

Функции Грина, входящие в выражение (4), для

разный порядок, так называемые грибовидные диа-

2 ≤ i ≤ n можно представить в виде

граммы.

Первые члены в (6) и (7) известны как лестнич-

k²⁰

ные диаграммы. Они описывают основной, некоге-

T (r^′i-1 - r^′i) =

T (r_i-1 - r_i) ×

рентный вклад в многократное рассеяние, посколь-

4π

(

)

ку осциллирующие фазовые коэффициенты комп-

× exp

iq_i ·

r^′′′i-1 - r^′′′i

(8)

лексно-сопряженных пар коррелированных полей в

точности компенсируются.

где q_i вектор рассеяния на i-м шаге,

Максимально-перекрестные диаграммы тополо-

гически эквивалентны своим лестничным аналогам.

T (r) = r^-1 exp ((ik₀ - µ_s/2) r)

Эта эквивалентность продемонстрирована на рис. 1

и на верхнем изображении на рис. 2 если повер-

основной пространственно-зависимый фактор

нуть нижнюю линию на 180^◦, то одна диаграмма

функции Грина, коэффициент рассеяния µ_s = l^-1s.

342

ЖЭТФ, том 164, вып. 3 (9), 2023

Когерентное обратное рассеяние с пониженной пространственной ...

r₃

r₂

r₁

r₃

r₂

r₁

r'₁

r'₂

r'₃

r'₂

r'₁

r₃

r₂

r₁

r₃

r₂

r₁

r'₁

r'₃

r'₂

r'₃

r'₁

r₃

r₂

r₁

r₃

r₂

r₁

r'₂

r'₁

r'₃

r'₁

r'₂

Рис. 2. Топологическая эквивалентность диаграмм третьего порядка: лестничной D^(3)L и максимально-перекрестной D^(3)C,

и двух пар интерференционных диаграмм третьего порядка D^(3)CLL и D^(3)LC, D^(3)LLC и D^(3)CL

В силу трансляционной инвариантности пар-

ная корреляционная функция 〈δε(r^′i) δε^∗(r^′′i)〉 зави-

сит только от разностных переменных r^′′′i:

G(r^′′′) = 〈δε(r_i + r^′′′i/2) δε^∗(r_i - r^′′′i/2)〉.

Для изотропного рассеяния радиус корреляции r_c

Рис. 3. Справа изображена пара плоских волн, падающих

много меньше длины волны λ; мы имеем, фактиче-

на систему, а слева те же волны, покидающие систему

ски, G(r^′′′) ∝ δ(r^′′′), и интегралы по r^′′′i берутся три-

в результате многократного рассеяния

виально. В общем случае эти интегралы производят

при α = δ = 1 и β = γ = 3 максимально-перекрест-

фазовую функцию p(q) с дополнительным множи-

ную диаграмму D^(3)1331 = D^(3)C, остальные случай-

телем µ_s, поскольку, согласно оптической теореме,

но-интерференционные диаграммы

∫

k⁴⁰

(3)

dr^′′′G(r^′′′) exp(iq · r^′′′) = µ_sp(q).

D₁

=D^(3)CLL, D^(3)1132 = D^(3)LC,

231

(4π)²

(3)

D₁

332

=D^(3)LLC, D^(3)1233 =D^(3)CL.

Представим вклады диаграмм второго порядка,

лестничную и перекрестную, в виде

Определим геометрию рассеяния. Пусть плоская

∫

∞

∫

волна падает нормально на слой 0 < z < T, где z

(2)

J⁽²⁾

(k_i, k_f ) = µ_s dz₁ dr₂ D_α

(k_i, k_f ).

декартова координата нормальной границы разде-

αβγδ

βγδ

ла, а T толщина слоя.

Определим элементы падающих и рас-

Индексы α, β, γ, δ обозначают номера координат

сеянных плоских волн, k_i

(0, 0, k) и

(см. рис. 3) пары входящих k_i и выходящих k_f плос-

kf = k(sinθs,0,-cosθs), где θs

угол обратного

ких волн. Набор α = β = 1, γ = δ = 2 соответ-

рассеяния, отсчитываемый от обратного направ-

ствует лестничной диаграмме D^(2)1122 = D^(2)L, набор

ления, с учетом затухания в рассматриваемой

α = δ = 1, β = γ = 2 перекрестной диаграмме

геометрии обратного рассеяния как

D^(2)1221 = D^(2)C. Для диаграмм третьего порядка

(

)

µ_s

∫

∞

∫

F^(in)αβ(k_i) = exp ik_i · (r_α - r_β) -

(z_α + z_β) ,

(3)

(

)

J⁽³⁾

(k_i, k_f ) = µ_s dz₁ dr₂ dr₃ D_α

(k_i, k_f ).

µ_s

αβγδ

βγδ

F^(out)γδ(k_f ) = exp ik_f · (r_γ - r_δ) -

(z_γ + z_δ)

Нижние индексы задают при α

= β

= 1 и

Лестничная диаграмма содержит тривиальную

γ = δ = 3 лестничную диаграмму D^(3)1133 = D^(3)L, а

зависимость, обусловленную угловой зависимостью

343

В. Л. Кузьмин, А. Ю. Вальков, Ю. А. Жаворонков

ЖЭТФ, том 164, вып. 3 (9), 2023

коэффициентов рассеяния µ_s/ cos θ при наклонном

4. МОДЕЛИРОВАНИЕ МЕТОДОМ

выходе или входе. Максимально-перекрестная диа-

МОНТЕ-КАРЛО

грамма содержит сильную угловую зависимость из-

В численном моделировании интенсивность рас-

за интерференции падающего и рассеянного полей.

сеяния представляется в виде ряда по кратностям

Другие диаграммы могут демонстрировать сильную

рассеяния:

угловую зависимость также из-за интерференции

∑

входящего и выходящего полей. Для декартовой па-

J (k_i, k_f) =

J⁽ⁿ⁾(k_i, k_f),

(10)

раметризации волновых векторов получаем

n=1

где J⁽ⁿ⁾(k_i, k_f) вклад n-го порядка рассеяния.

F^(L)in(r₁, r₁)F^(L)out(r_n, r_n) = exp(-µ_s(z₁ + z_n/cosθ_f)) ,

В рамках лестничного приближения на основе

(

F^(C)in(r₁, r_n)F^(C)out(r_n, r₁) = exp

ik(x₁ - x_n)sinθ_f +

уравнения Бете - Солпитера мы представляем член

)

+ ik(z₁ - z_n)(1 - cosθ_f) - µ_s(z₁ + z_n/cosθ_f)

n-го порядка J⁽ⁿ⁾(k_i, k_f) как среднее по выборке N_ph

падающих фотонов:

Тогда лестничные и максимально-перекрестные

диаграммы второго и третьего порядков примут вид

∑

J⁽ⁿ⁾(k_i, k_f) =

W⁽ⁱ⁾ⁿ(ki, kf) ×

D^(2)L(k_i, k_f) = µ^2s p(q₁)p(q₂)F⁽ⁱⁿ⁾¹¹(k_i)Λ(r₁₂)F^(out)22(k_f),

ph i=1

(

)

D^(2)C(k_i, k_f) = µ^2s p(q₁)p(q₂)F⁽ⁱⁿ⁾¹²(k_i)Λ(r₁₂)F^(out)21(k_f),

×p

k⁽ⁱ⁾

(11)

n n-1

kf fBLB(kf , zⁿⁱ⁾),

D^(3)L(k_i, k_f) = µ^3s p(q₁)p(q₂)p(q₃)F⁽ⁱⁿ⁾¹¹(k_i)Λ(r₁₂) ×

где

n (ki, kf) и zⁿⁱ⁾ соответственно вес и рассто-

× Λ(r₂₃)F^(out)33(k_f),

(i)

n n-го акта рассеяния.

яние до границы от точки r

D^(3)C(k_i, k_f) = µ^3s p(q₁)p(q₂)p(q₃)F⁽ⁱⁿ⁾¹³(k_i)Λ(r₁₂) ×

Экспоненциальный множитель Бугера - Ламберта -

× Λ(r₂₃)F^(out)31(k_f),

Бера f_BLB(k_f , zⁿⁱ⁾) описывает затухание рассеянно-

го излучения, распространяющегося от точки zⁿⁱ⁾

Стандартные алгоритмы моделирования много-

последнего события рассеяния до границы в при-

кратного рассеяния

[33] используют оператор

ближении Фраунгофера. Он зависит от оптических

Λ(r)

T (r)|²

пропагатор уравнения переноса

параметров среды на пути фотона, движущегося к

излучения (включая уравнение Бете - Солпитера в

границе. Для однородной среды

лестничном приближении).

В некогерентных вкладах, описываемых лест-

f_BLB(k_f , z) = exp(-µ_sz/cosθ_s).

ничными диаграммами, остаются только затуха-

ющие множители, а осциллирующие полностью

Вес

n представляет собой случайное зна-

компенсируются. В максимально-перекрестных диа-

чение многократного пространственного интегра-

граммах полная компенсация осциллирующих мно-

ла, возникшего как итерация n-го порядка урав-

жителей реализуется только при рассеянии строго

нения Бете - Солпитера. Вычисляя его, можно смо-

назад, k_f

= -k_f . Все остальные диаграммы со-

делировать стохастическую последовательность то-

держат осциллирующие множители из-за несовпа-

чек рассеяния r₁, . . . , r_n. Полную сумму лестнич-

дения как пар входящих или выходящих полей, так

ных диаграмм, практически не зависящую от уг-

и функций Грина. Например, диаграмма

ла обратного рассеяния (кроме зависимости вида

exp(-µ_sz/cosθ_s) для входящих и выходящих плос-

D₁₃₃₂(k_i, k_f) = µ^3s p(q₁)p(q₂)p(q₃)F⁽ⁱⁿ⁾¹³(k_i) ×

ких волн) в области пика обратного рассеяния, обо-

значим как

× F^(out)32(k_f )Λ(r₁₂)T (r23)T∗(r31)

(9)

∑

(n)

J_L(k_i, k_f) =

J_L

(k_i, k_f),

содержит нескомпенсированный осциллирующий

множитель. Осциллирующие множители сильно

затрудняют прямое стохастическое моделирование.

а полную сумму максимально-перекрестных диа-

Очевидно, что при больших значениях случайной

грамм как

переменной вида k(R₂₃ - R₁₃), где R_ij = |r_i - r_j |,

∑

J_C(k_i, k_f) =

J^(n)C(k_i, k_f).

или k(z₁ - z₃) вкладом таких конфигураций можно

пренебречь. Поэтому при моделировании мы выпол-

няли обрезание высокочастотных стохастических

Метод Монте-Карло в теории переноса излуче-

флуктуаций.

ния основан на хорошо известной процедуре обрат-

344

ЖЭТФ, том 164, вып. 3 (9), 2023

Когерентное обратное рассеяние с пониженной пространственной ...

ного преобразования [33,34]. Обратное преобразова-

часто используемой в биомедицинских приложени-

ние позволяет перестроить интеграл с экспоненци-

ях является фазовая функция Хеньи - Гринштейна,

альным распределением по полубесконечному ин-

важным достоинством которой является тот факт,

тервалу [0, ∞] в интеграл по случайной величине,

что обратное преобразование кумулятивной функ-

равномерно распределенной в единичном интерва-

ции (14) выполняется аналитически в явном виде.

ле [0, 1]. В рамках стандартного алгоритма затухаю-

Описание однократного рассеяния фазовой функ-

щий множитель в пропагаторе Λ(r) представляется

цией Рэлея - Ганса позволяет моделировать оптиче-

как плотность вероятности экспоненциального рас-

ские свойства ткани или биофантома на основе фи-

пределения f(r) = µ^-1s exp(-µ_sr) случайной величи-

зической модели суспензии, но требует существен-

ны r, описывающей расстояние между двумя актами

ных математических усложнений [35, 36]. В данной

рассеяния. Кумулятивная функция экспоненциаль-

работе мы ограничимся при моделировании функ-

ного распределения легко вычисляется:

цией Хеньи - Гринштейна. Отметим, что интенсив-

ность рассеяния, вычисленная с помощью выраже-

ξ(r) = 1 - exp(-µ_sr),

(12)

ния (11), может быть интерпретирована как сред-

где ξ или ξ^′ = 1 - ξ случайные величины, равно-

нее значение экспоненты exp(-µ_szⁿⁱ⁾)/cos θ_s, которая

мерно распределенные в единичном интервале [0, 1].

описывает затухание фотона, возвращающегося из

Обратное преобразование дает

объемной среды к границе после n актов рассеяния.

Поскольку среда ограничена, фотон возвращается в

r = -µ^-1slnξ^′.

(13)

среду, если он достиг границы, в силу закона отра-

В диаграммах третьего порядка и выше появля-

жения [15]. Таким образом, вес

n является про-

ются функции Грина без соответственной комплек-

изведением коэффициентов отражения Френеля. В

сно-сопряженной пары, что не позволяет описать

частном случае отсутствия отражения и адсорбции

перенос излучения в терминах пропагатора Λ(r).

он равен единице или нулю.

Для таких диаграмм мы используем метод обрат-

ного преобразования пространственных интегралов,

содержащих полевые функции Грина, а не стандарт-

5. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

ный метод, ориентированный на пропагатор Λ(r).

Основным препятствием для применения эф-

Представляя экспоненциально затухающий множи-

фекта усиления КОР в биомедицинской практике

тель функции Грина как величину, распределенную

является малая ширина пика. В работах [24, 25]

с плотностью вероятности

показано, что ширина пика значительно увеличи-

вается при понижении пространственной когерент-

f (r) = (µ_s/2)^-1 exp(-µ_sr/2),

ности падающего излучения, а его высота одно-

получим r = -(µ_s/2)^-1 ln ξ^′. Для случая, когда все

временно уменьшается. Понижение пространствен-

функции Грина можно ¾спарить¿ со своими сопря-

ной когерентности в рассматриваемой нами расчет-

женными, оба выбора функции распределения в ме-

ной схеме моделируется путем понижения макси-

тоде обратного преобразования в виде как про-

мального числа актов рассеяния n_sc . Для сравне-

пагатора переноса излучения Λ(r), так и функции

ния результатов моделирования КОР с известными

Грина, дают одинаковые результаты.

экспериментальными результатами мы взяли значе-

К интегралам по угловым переменным метод об-

ния оптических параметров из известной работы [6]:

ратного преобразования применяется, как обычно.

λ = 0.633мкм, показатель преломления n = 1.33,

Угловую переменную t = cos θ, представленную в

g = 0.93, l_s = 2.6мкм.

виде случайной величины, распределенной на осно-

На рис.4 показана индикатриса КОР, а именно,

ве заданной фазовой функции, заменяют на рав-

угловая зависимость величины

номерно распределенную случайную величину χ.

Определим кумулятивную функцию, распределен-

h(θ_s) = (J_BG + J_CBS (θ_s))/J_BG,

ную в интервале [0, 1]:

представляющей собой отношение суммарной ин-

∫^t

тенсивности обратного рассеяния к фоновой интен-

χ = F(t) = 2π p(t^′)dt^′,

(14)

сивности J_BG за пределами угловой области усиле-

-1

ния КОР. На практике J_BG это вклад лестничных

где p(t)

фазовая функция. Затем выполняем

диаграмм, J_BG = J_L. В интересующей нас срав-

обратное преобразование t

= F^-1(χ). Наиболее

нительно узкой области углов J_BG можно считать

345

ЖЭТФ, вып. 3 (9)

В. Л. Кузьмин, А. Ю. Вальков, Ю. А. Жаворонков

ЖЭТФ, том 164, вып. 3 (9), 2023

Рис. 4. Зависимости относительной интенсивности обрат-

Рис. 5. То же, что на рис. 4 для изотропного рассеяния

ного рассеяния h(θs) = (JBG + JCBS (θs))/JBG от угла

(g = 0)

обратного рассеяния θs при различных значениях мак-

симального числа актов рассеяния: сплошная линия

В частности, в работе [6] получено эксперимен-

nsc = 5 · 10⁵; пунктирная линия nsc = 340; штриховая

тальное значение h(0) = 1.64, которое совпадает

линия nsc = 200; штрихпунктирная линия nsc = 100.

с результатом нашего численного расчета при

Параметр анизотропии рассеяния g = 0.93

выборе nsc = 340. Рассчитанная при этом значе-

нии n_sc полуширина пика КОР θ_HW

= 1.50 мрад

Таблица. Уменьшение высоты и угловое уширение

также неплохо согласуется с эксперименталь-

пика КОР с понижением когерентности

ным значением θ_HW

= 1.58 мрад [6]. Добавим,

что согласно априорной оценке имеем величину

n_sc

h(0) θ_HW , мрад L_c, мкм

θ_HW ∼ W₀ = (kl_tr )^-1 ≈ 2.0 мрад, где l_tr = l_s/(1 - g),

5 · 10⁵

1.99

0.78

3.9 · 10⁵

а g параметр анизотропии среды. Диффузионное

340

1.64

1.50

884

приближение дает θ_HW

= 3.3 мрад [8], что по

200

1.55

1.73

520

порядку величины также согласуется с нашими

100

1.39

2.36

260

данными из таблицы.

Отметим, что в работе [25] эффект КОР при-

менялся в реальной биомедицинской практике при

величиной, не зависящей от θ_s. Мы использовали

nsc = 5 · 105, что соответствует практически бес-

изучении рака тканей толстой кишки человека с по-

конечной длине когерентности, а также n_sc = 340,

мощью низкокогерентного оптического излучения.

n_sc = 200 и n_sc = 100 для учета пониженной коге-

Полученная в этой работе высота пика h(0) состав-

рентности излучения. Значение n_sc = 340 соответ-

ляла около 1.07, а полуширина θ_HW ≈ 3.5 мрад. На-

ствует длине когерентности L_c ∼ n_scl_s ≈ 880 мкм,

ше моделирование методом Монте-Карло дает та-

а n_sc = 200 и n_sc = 100 соответственно длинам

кую величину h(0) при n_sc = 20.

L_c = 520 мкм и L_c = 260 мкм. Видно, что при та-

При θ_s

≳ 2мрад индикатрисы на рис.4 для

ких L_c происходит заметное уширение пика КОР,

n_sc = 5 · 10⁵, 340, 200, 100 практически совпадают,

причем ширина пика возрастает с уменьшением ко-

т. е. роль очень высоких кратностей рассеяния при

герентности излучения. В таблице приведены соот-

больших углах существенно понижается. Это про-

ветствующие высоты h(0), полуширины θ_HW пиков

исходит за счет того, что, во-первых, положитель-

КОР (HWHM полуширина на полувысоте пика,

ные вклады от лестничных диаграмм при больших

J_CBS(θ_HW ) = J_CBS(0)/2) и длины когерентности

углах рассеяния достаточно быстро убывают (по за-

L_c.

кону n

[37]), а во-вторых, потому, что в этой об-

Практически во всех экспериментальных ра-

ласти максимально-перекрестные диаграммы прак-

ботах, начиная с пионерских

[6-8, 29], значение

тически полностью компенсируют друг друга, по-

параметра усиления КОР h(0) заметно не достигает

скольку при больших углах их вклад может быть

предельного теоретического значения h(0)

= 2.

как положительным, так и отрицательным.

346

ЖЭТФ, том 164, вып. 3 (9), 2023

Когерентное обратное рассеяние с пониженной пространственной ...

терминах полевых функций Грина позволило нам

учесть вклады нелестничных диаграмм, которые мы

называем случайно-интерференционными. Модели-

рование в рамках стандартного описания КОР в тер-

минах оператора переноса Λ(r) и описания на основе

функций Грина показало, что при малых углах об-

ратного рассеяния относительный вес случайно-ин-

терференционных диаграмм составляет около 1%,

достигая 5-8% при углах рассеяния около 50 мрад,

актуальных при применении излучения с понижен-

ной когерентностью в биомедицинской практике.

Для вклада интерференционных диаграмм мы

получили оценку ширины пика обратного рассея-

ния, в несколько раз превосходящую априорную

оценку W ∼ W₀ = (kl_tr )^-1. Это свойство уширения

пика КОР наблюдается только у части перекрест-

(3)

ных диаграмм, D₁

и D^(3)1231, с явными корреляци-

Рис. 6. Вклад низших нелестничных диаграмм в h(θs) для

332

ями между входящими и выходящими полями и мо-

изотропного рассеяния в области θs ≥ 10 мрад. Длинные

штрихи сумма всех диаграмм второго и третьего поряд-

жет быть использовано в приложениях. Диаграммы

ков, короткие штрихи сумма максимально-перекрестных

D^(3)1132 и D^(3)1233 представляют собой первую поправку

диаграмм второго и третьего порядков; штрихпунктирная

к лестничному приближению уравнения Бете - Сол-

линия двукратная перекрестная диаграмма

питера. Указанные диаграммы практически не за-

висят от угла рассеяния в рассматриваемом угловом

Для сравнения на рис. 5 приводятся индикатри-

диапазоне КОР, вплоть до θ_s = 60 мрад, и представ-

сы в той же области углов, рассчитанные для изо-

ляют фактически фон пика КОР.

тропной среды. Здесь эффект уширения пика КОР

Принимая максимальное число событий рассея-

также имеет место, но для тех же длин когерентно-

ния n_sc в качестве параметра, определяющего про-

сти выражен заметно меньше.

странственную длину когерентности L_c = n_scl_s, где

Для случая изотропного рассеяния мы также

l_s

длина рассеяния, мы показали, что использова-

рассчитали угловую зависимость h(θ_s), включая

ние низкокогерентного излучения позволяет полу-

вклады нелестничных диаграмм третьего порядка

чить конус КОР с шириной и относительной высо-

(последние четыре диаграммы на рис. 2). Оказалось

той, близкими к экспериментальным значениям [6].

(рис. 6), что их вклад в области пика КОР составля-

Для сравнения с экспериментом мы взяли дан-

ет около 1% от вклада максимально-перекрестных

ные из работы [7]. Однако значения оптических па-

диаграмм, возрастая до

5-8% при углах более

раметров, представляющих в настоящее время инте-

50 мрад.

рес в биомедицинской практике, лежат в другом ин-

На рис.6 мы приводим также график угловой

тервале, где рассеяние существенно меньше. Для ти-

зависимости вклада двукратного рассеяния. Этот

пичных биологических тканей коэффициент рассея-

вклад составляет менее 10% от высоты пика КОР,

ния лежит в пределах от 200 до 500 см^-1. Для длины

что делает неоправданным предположение [24,25] об

рассеяния это дает интервал 20-50 мкм, который на

основной роли двукратного рассеяния в образова-

порядок больше типичных значений длин рассеяния

нии пика КОР для низкокогерентного излучения.

суспензий, использованных в работах [6-8]. Для ре-

ального использования КОР в биотканях требуется

добиться более широкой угловой области усиления

6. ЗАКЛЮЧЕНИЕ

обратного рассеяния, в том числе за счет различ-

ных механизмов понижения когерентности, вклю-

Мы выполнили моделирование эффекта обрат-

чая рассмотренный в настоящей работе.

ного когерентного рассеяния в широкой угловой

области. Развит алгоритм моделирования методом

Финансирование. Исследование выпол-

Монте-Карло на основе итерационного решения

нено при финансовой поддержке Россий-

уравнения Бете - Солпитера. Использование опи-

ского научного фонда, грант

№23-22-00035,

сания интенсивности многократного рассеяния в

https://rscf.ru/project/23-22-00035/.

347

В. Л. Кузьмин, А. Ю. Вальков, Ю. А. Жаворонков

ЖЭТФ, том 164, вып. 3 (9), 2023

ЛИТЕРАТУРА

19.

H. Wabnitz, J. Rodriguez, I. Yaroslavsky, A. Yaro-

slavsky, and V. V. Tuchin, Handbook of Optical

K. M. Watson, J. Math. Phys. 10, 688 (1969).

Biomedical Diagnostics. Light-Tissue Interaction,

Vol. 1, SPIE Press, Bellingham, Washington (2016).

D. А. de Wolf, IЕЕЕ Trans. Antennas Propag. 19, 254

(1971).

20.

D. A. Boas, S. Sakadzic,

J. Selb,

P. Farzam,

M. A. Franceschini, and S. A. Carp, Neurophotonics

Ю. Н. Барабаненков, Изв. вузов, Радиофизика 16,

3, 031412 (2016).

88 (1973).

21.

S. Etemad, R. Thompson, M. J. Andrejco, S. John,

А. Г. Виноградов, Ю. А. Кравцов, В. И. Татарский,

and F. C. MacKintosh, Phys. Rev. Lett.

59,

1420

Изв. вузов, Радиофизика 16, 1064 (1973).

(1987).

Y. Kuga and A. Ishimaru, J. Opt. Soc. Am. A 1, 831

22.

T. Okamoto and T. Asakura, Opt. Lett.

21,

369

(1996).

(1984).

23.

A. Wax, C. Yang, and J. A. Izatt, Opt. Lett. 28, 1230

M. P. Van Albada and A. Lagendijk, Phys. Rev. Lett.

(2003).

55, 2692 (1985).

24.

Y. L. Kim, P. Pradhan, H. Subramanian, Y. Liu,

P.-E. Wolf and G. Maret, Rev. Lett. 55, 2696 (1985).

M. H. Kim, and V. Backman, Opt. Lett. 31,

1459

(2006).

E. Akkermans, P. Wolf, R. Maynard, and G. Maret,

J. de Phys. 49, 77 (1988).

25.

Y. L. Kim, Y. Liu, V. M. Turzhitsky, H. K. Roy,

R. K. Wali, H. Subramanian, P. Pradhan, and

D. J. Pine, D.A. Weitz, P. M. Chaikin, and E. Her-

V. Backman, J. Biomed. Opt. 11, 041125 (2006).

bolzheimer, Phys. Rev. Lett. 60, 1134 (1988).

26.

H. Subramanian, P. Pradhan, Y.L. Kim, Y. Liu,

10.

P. Wolf, G. Maret, E. Akkermans, and R. Maynard,

X. Li, and V. Backman, Appl. Opt. 45, 6292 (2006).

J. de Phys. 49, 63 (1988).

27.

M. Xu, Opt. Lett. 33, 1246 (2008).

11.

T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh,

Rep. Prog. Phys. 73, 076701 (2010).

28.

J. Liu, Z. Xu, Q. Song, R. L. Konger, Y. L. Kim,

J. Biomed. Opt. 15, 037011 (2010).

12.

A. P. Tran, S. Yan, and Q. Fang, Neurophoton. 7,

29.

D. S. Wiersma, M. P. van Albada, and A. Lagendijk,

015008 (2020).

Phys. Rev. Lett. 75, 1739 (1995).

13.

A. Sabeeh and V. V. Tuchin, J. Biomed. Photon.

30.

С. М. Рытов, Ю. А. Кравцов, В. И. Татарский,

Engineer. 6, 040201 (2020).

Введение в статистическую радиофизику. Часть

2. Случайные поля, Наука, Москва (1978).

14.

В. В. Тучин, Оптика биологических тканей. Ме-

тоды рассеяния света в медицинской диагности-

31.

V. Kuzmin, V. Romanov, and L. Zubkov, Phys.Rep.

ке, IPR Media, Москва (2021).

248, 71 (1994).

15.

T. H. Pham, O. Coquoz, J. B. Fishkin, E. Anderson,

32.

В. Л. Кузьмин, ЖЭТФ 127, 1173 (2005).

and B. J. Tromberg, Rev. Sci. Instrum.

71,

2500

(2000).

33.

L. Wang, S. L. Jacques, and L. Q. Zheng, Comput.

Meth. Prog. Bio. 47, 131 (1995).

16.

R.C. Haskell, L.O. Svaasand, T.-T. Tsay, T.-C. Feng,

34.

L. Devroye, Non-Uniform Random Variate Genera-

B. J. Tromberg, and M. S. McAdams, J. Opt. Soc.

tion, Springer, New York (1986).

Am. A 11, 2727 (1994).

35.

В. Л. Кузьмин, А. Ю. Вальков, Письма в ЖЭТФ

17.

D. A. Boas, L. E. Campbell, and A. G. Yodh, Phys.

105, 261 (2017).

Rev. Lett. 75, 1855 (1995).

36.

В. Л. Кузьмин, А. Ю. Вальков, Л. А. Зубков,

18.

X. Cheng, D. Tamborini, S.A. Carp, O. Shatrovoy,

ЖЭТФ 155, 460 (2019).

B. Zimmerman, D. Tyulmankov, A. Siegel, M. Black-

well, M. A. Franceschini, and D. A. Boas, Opt. Lett.

37.

T. M. Nieuwenhuizen and J. M. Luck, Phys. Rev. E

43, 2756 (2018).

48, 569 (1993).

348