Files
optical-trap-forces/README.ru.md
2023-03-28 23:25:44 +04:00

200 lines
17 KiB
Markdown
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# Силы в оптической ловушке
В данном репозитории приведен код на Python и Matlab для расчета однолучевой лазерной ловушки, также называемой оптическим пинцетом, захватывающей микрометровую диэлектрическую сферу. Это можно рассматривать как простую модель системы, описывающей лазерный захват и манипулирование живыми организмами и органелл внутри клеток. Градиентная сила и сила рассеяния определены для пучка с равномерным распределением интенсивности гауссова и бесселевых пучков (равномерный и бесселев только в коде).
Для диэлектрических сфер диаметром много больше длины волны можно использовать геометрическую оптику для расчета сил, действующих на нее со стороны пучка. Пучок света раскладывается на отдельные лучи, каждый со своим значением интенсивности, направлением и состоянием поляризации. Они распространяются по прямым линиям в однородной среде с однородным показателем преломления. Дифракционные эффекты в этом рассмотрении отсутствуют.
Простая геометрооптическая модель однолучевой градиентной ловушки показана на рис. 1.
|![](images/fig_1.png)|
|---|
|Рис. 1. (A) Однолучевая градиентная оптическая ловушка в рассмотрении геометрической оптики с фокусом $f$, расположенным на оси $Z$. (B) Геометрия падающего луча, дающего вклад в силы $F_g$ и $F_S$|
Ловушка состоит из падающего параллельного пучка лучше произвольного модового состава и поляризации, который попадает в высокоапертурный микрообъектив и фокусируется в фокусе $f$. Максимальный угол отклонения лучей на выходе микрообъектива определяется его числовой апертурой
$$NA = n_1 \cdot \sin{\sigma_A}$$
где $n_1$ показатель преломления среды;
$\sigma$ апертурный угол, рад.
Для типичного водно-иммерсионного объектива со значением $NA = 1.25$ и $n_1 = 1.33$ максимальный угол отклонения луча на выходе микрообъектива будет равен
$$\phi_{max} = \arcsin{\frac{NA}{n_1}} = 70^\circ$$
Для единичного луча мощностью $P$, падающего на диэлектрическую сферу под углом падения $\theta$ с импульсом $\frac{n_1 P}{c}$ ($c$ скорость света в вакууме), силы, которые он оказывает на сферу, определяются формулами:
$$F_s = \frac{n_1 P}{c} \left[ 1 + R \cos{2\theta} - \frac{T^2[\cos(2\theta - 2r) + R\cos{2\theta}]}{1 + R^2 + 2R\cos{2r}} \right] = \frac{n_1 P}{c} Q_s$$
$$F_g = \frac{n_1 P}{c} \left[ R \sin{2\theta} - \frac{T^2[\sin(2\theta - 2r) + R\sin{2\theta}]}{1 + R^2 + 2R\cos{2r}} \right] = \frac{n_1 P}{c} Q_g$$
где $F_s$ сила рассеяния, Н;
$F_g$ градиентная сила (возвращающая сила), Н;
$Q_s$ и $Q_g$ эффективности этих сил, соответственно;
$\theta$ и $r$ углы падения и отражения, соответственно, рад;
$R$ и $T$ коэффициенты отражения и пропускания, определяющиеся формулами Френеля:
$$R(\theta, \psi) = R_{\perp} \cos{\psi} + R_{\parallel} \sin{\psi} = \frac{\sin^2(\theta - r)}{\sin^2(\theta + r)}\sin{\psi} + \frac{\tan^2(\theta - r)}{\tan^2(\theta + r)}\cos{\psi}$$
$$T(\theta, \psi) = 1 - R(\theta, \psi)$$
где $\psi$ угол поворота плоскости поляризации, рад.
Углы падения и отражения связаны между собой законом Снеллиуса:
$$n_1 \sin{\theta} = n_2 \sin{r}$$
где $n_2$ показатель преломления материала сферы.
Сила $F_g$ является полезной силой, возвращающей частицу в положение равновесия. $F_s$ сила, выталкивающая частицу из ловушки. $F_g$ должна быть больше, чем $F_s$ для стабильного захвата.
Видно, что силы, действующие на частицу, зависят от поляризации света, т.к. от нее зависят коэффициенты Френеля.
## Эффективность луча
Построим графики для эффективностей $Q$ полученных сил. Введем относительный показатель преломления $n = \frac{n_2}{n_1}$. Для расчета будем использовать значения $NA = 1.25$, откуда $\phi_{max} = 70^\circ$;
$n_1 = 1.33$; $𝑛 = 1.2$, откуда $n_2 = 1.6$ (полистирол). Величину $f$ и радиус частицы вводить не нужно, т.к. в геометрооптическом рассмотрении от их значений результат не зависит.
|![](images/fig_2.png)|
|---|
|Рис. 2. Эффективности сил, действующих на частицу при $n = 1.2$|
На рис. 2 $Q = \sqrt{Q_s^2 + Q_g^2}$ эффективность результирующей силы, действующей на частицу. Этот график показывает вклад каждого луча, падающего на частицу под определенным углом. Видно, что максимальное значение $Q_g$ находится примерно при значении угла падения $70^\circ$, что показывает необходимость использования высоких значений $NA$. Для сравнения приведем графики для других $n$ на рис. 3 и 4.
|![](images/fig_3.png)|
|---|
|Рис. 3. Эффективности сил, действующих на частицу при $n = 1.1$|
|![](images/fig_4.png)|
|---|
|Рис. 4. Эффективности сил, действующих на частицу при $n = 1.8$|
Видно, что при $n \rightarrow 1$ необходимый максимальный угол падения становится слишком высок, а при $n \rightarrow 1.4\$ $Q_s$ равна или превышает $Q_g$ в большей части диапазона углов, что говорит о сложности получения стабильного захвата.
## Сила вдоль оси Z
Для получения результирующей силы, действующей на частицу, необходимо просуммировать вклад всех лучей, т.е. проинтегрировать. Интегрирование проводится в цилиндрической системе координат, следовательно, Якобиан равен $J = r$. Проекции сил $F_g$ и $F_s$ на ось $Z$ в данном случае определяются как
$$F_{gZ} = -F_g \sin{\phi}$$
$$F_{sZ} = F_s \cos{\phi}$$
Угол $\phi$ угол отклонения луча после прохождения микрообъектива (см. рис. 1). Угол падения $\theta$ определяется как
$$\theta = \arcsin\left(\frac{z}{a} \sin{\phi}\right)$$
Таким образом, результирующая сила, действующая на частицу вдоль оси $Z$ определяется интегралом
$$F_Z^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} (F_{gZ} + F_{sZ}) r d\beta dr$$
где $r_{max} = f \tan{\phi_{max}}$.
Аналогично, отдельные компоненты градиентной силы и силы рассеивания определяются следующими интегралами:
$$F_{gZ}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{gZ} r d\beta dr$$
$$F_{sZ}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{sZ} r d\beta dr$$
где $\frac{1}{\pi r_{max}^2}$ есть нормировочный коэффициент.
Эффективности сил, действующих вдоль оси $Z$ $Q_{tz}$, $Q_{gz}$, $Q_{sz}$ аналогично получаются путем деления сил на величину $\frac{n_1 P}{c}$.
Построим графики на основе данных, приведенных выше (рис. 5).
|![](images/fig_5.png)|
|---|
|Рис. 5. Силы, возникающие при смещении частицы вдоль оси $Z$|
Из графика видно, что при смещении частицы из положения равновесия растут силы, которые возвращают частицу обратно. Максимальные значения сила приобретает в координатах $z = 1.06\ \mu m$ и $z = 1.02\ \mu m$. От конкретных значений этих сил зависит стабильность захвата и жесткость ловушки. Также можно заметить, что положение равновесия частицы немного смещено относительно нуля и то, что благодаря симметрии в данном случае сила не зависит от поляризации света.
## Сила вдоль оси Y
Таким же образом можно получить силы, действующие в латеральной плоскости. Здесь мы рассматриваем ось $Y$, но это же рассмотрение справедливо для смещения вдоль любого направления в плоскости, перпендикулярной направлению лазерного пучка. Однако, в этом случае появляются новые углы, показанные на рис. 6.
|![](images/fig_6.png)|
|---|
|Рис. 6. (A) Геометрия ловушки с фокусом $f$, расположенным на расстоянии по оси $Y$ от центра сферы $O$. (B) Геометрия плоскости падения, показывающая направления градиентной силы $F_g$ и силы рассеивания $F_s$ для падющего луча|
Они определяются следующими формулами:
$$\alpha = 90^\circ - \arctan{\frac{r}{f}}$$
$$\gamma = \arccos\left(\cos{\alpha} \cos{\beta}\right)$$
$$\theta = \arcsin\left(\frac{y}{a} \sin{\gamma}\right)$$
$$\mu = \arccos\left(\frac{\tan{\alpha}}{\tan{\gamma}}\right)$$
В данном случае силы, действующие на частицу, будут зависеть от поляризации. Вклад каждой поляризации определяется коэффициентами
$$f_\parallel = \left(\cos{\beta} \sin{\mu - \sin{\beta} \cos{\mu}}\right)^2$$
$$f_\perp = \left(\cos{\beta} \cos{\mu - \sin{\beta} \sin{\mu}}\right)^2$$
Однако, мы будем полагать, что свет имеет круговую поляризацию, тогда
$$f_\parallel = f_\perp = 1/2$$
Здесь проекции сил на ось $Y$ определяются следующими выражениями:
$$F_{gY} = F_g \cos{\phi}$$
$$F_{sY} = F_s \sin{\gamma}$$
А соответствующие интегралы выражениями:
$$F_Y^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} (F_{gY} + F_{sY}) r d\beta dr$$
$$F_{gY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{gY} r d\beta dr$$
$$F_{sY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{sY} r d\beta dr$$
По тем же данным построим графики эффективностей сил вдоль оси $Y$ (рис. 7).
|![](images/fig_7.png)|
|---|
|Рис. 7. Силы, возникающие при смещении частицы вдоль оси $Y$|
Видно, что возвращающая сила симметрично действует при смещении частицы, и положение равновесия уже находится в нуле.
## Воздействие гауссовым пучком на сферическую частицу
В случае гауссова пучка (основной $TEM_{00}$ моды) распределение интенсивности выглядит следующим образом:
$$I(r) = \exp\left(-\frac{2r^2}{w_0}\right)$$
Графики сил, действующих вдоль осей $Z$ и $Y$ при условии $\varepsilon = 1.0\ (w_0 = r_{max})$ представлены на рис. 7 и 8.
Рассчитаем распределения для других значений $\varepsilon$ и сведем полученные данные в таблицу 1. В этой таблице величина $A = 1 - \exp\left(-\frac{2r_{max}^2}{w_0}\right)$ - доля мощности, попадающая на входной зрачок микрообъектива.
Таблица 1 Эффективности сил гауссова пучка при различных коэффициентах заполнения входного зрачка микрообъектива
|$\varepsilon$|$A$|$Q_{-z}^{max}$|$-Z_{max},\ \mu m$|$Q_y^{max}$|$Y_{max},\ \mu m$|$Q_{+z}^{max}$|$-Z_{max},\ \mu m$|$S_{E},\ \mu m$|
|---|---|---|---|---|---|---|---|---|
|1,7|0,5|-0,188|1,02|0,272|0,98|0,336|1,06|0,08|
|1,0|0,86|-0,098|1,04|0,182|0,98|0,180|1,06|0,1|
|0,727|0,98|-0,048|1,04|0,124|0,98|0,091|1,06|0,13|
|0,364|1,0|-0,005|1,16|0,044|0,98|0,014|1,3|0,31|
|0,202|1,0|-0,001|1,44|0,016|0,98|0,003|1,9|0,81|
Из таблицы 1 видно, что ловушка с $\varepsilon = 1.7$ имеет наибольшее значение силы выхода $Q_{-z}^{max} = 0.336$. Также видно, что ловушка с $\varepsilon = 1.7$ самая однородная ловушка, т.к. имеет наименьшую разницу между силами $Q_{-z}^{max}$ и $Q_{+z}^{max}$. Однако, стоит учитывать тот факт, что при таком соотношении радиуса пучка к апертуре объектива мощность, попадающая на частицу (см. величину $A$ в таблице 2), снижается. Поэтому оптимальным являются случаи $\varepsilon = 1.0$ и $\varepsilon = 0.727$, т.к. в пучке будет сосредоточена большая часть мощности.
## Как пользоваться
Все описанное выше оформлено в виде Python и Matlab файлов:
- Эффективности лучей: [trap_forces_efficiencies.py](Python/trap_forces_efficiencies.py) и [trap_forces_efficiencies.m](Matlab/trap_forces_efficiencies.m)
- Силы вдоль оси $Z$: [trap_forces_axial.py](Python/trap_forces_axial.py) и [trap_forces_axial.m](Matlab/trap_forces_axial.m)
- Силы вдоль оси $Y$: [trap_forces_transverse.py](Python/trap_forces_transverse.py) и [trap_forces_transverse.m](Matlab/trap_forces_transverse.m) (**занимает больше времени на вычисление!**)
Запустите один из этих файлов, и будут построены графики. Вы можете изменить начальные константы и провести свои эксперименты.
## Ссылки
- [A. Ashkin. Forces of a Single-Beam Gradient Laser Trap on a Dielectric Sphere in the Ray Optics Regime (1997). DOI: 10.1016/S0091-679X(08)60399-4](https://www.sciencedirect.com/science/article/abs/pii/S0091679X08603994)