1. 本記事の目的
本記事では、光圧という概念を簡単に紹介し、理解を促進するために光圧の理論式を導出します。
理論式の導出後、光圧を用いたシミュレーションを行います。
ただし、本記事では理論式の導出でかなりボリューミーになると思われるので、光圧を用いたシミュレーションは別記事にて紹介したいと思います。
2. 光圧について
この世界には、光圧と呼ばれる概念があります。
光圧とは、光と物質との相互作用により、物質が光から受ける力学的な力のことです。
光の圧力と書きますが、中二心がくすぐられる名称は素晴らしいとして、想像しがたい概念です。
光圧という概念が受け入れがたい理由の一つとして、我々の存在しているマクロな空間においては、実感することが不可能なほど小さな力だからだと思います。
どの程度小さいかですが、およそピコニュートン程度です。想像不可能ですね。
ちなみに、体重60 kgの人間同士が1 m離れたときに生じる万有引力がおよそナノニュートン程度なので、光圧はその1000倍小さい力になります。
おそらくですが、人間には感じることができません。
仮に光圧を感じれるのであれば、満員電車の中で常に周囲の人間からの万有引力を感じていることになります。
マクロな空間に存在する我々人間にとって実感のわきにくい光圧ですが、実はミクロな空間においては無視できない力となってきます。
現に、光圧を利用してミクロな微小物質を操作する技術が確立されています。
光ピンセットと呼ばれる技術です。
ピンセットで物質を捕まえて操作するかのように、光の圧力を利用してミクロな物質を捕捉して操作する、というイメージからこのように呼ばれています。
3. 光ピンセット
簡単に光ピンセットの原理について説明します。
光ピンセットはレンズを使って、レーザー光を集光します。
レーザー光を集光することで集光点付近に微小物質が捕捉されます。
以上が簡単な光ピンセットの原理になります。
直観的には、なぜ光を集光することで微小物質が捕捉されるのか、というのが疑問です。
答えは、微小物質に光圧が働いたから、です。
これだけ説明して、「なるほど微小物質に光圧が働いたから微小物質は捕捉されうるのか！」、となる人はいますか。(※いましたら、以下の記事内容を読む価値は特にないと思われます。)
このように直観に頼りずらい光圧という概念の理解を促進するために、光圧の理論式を導出していきます。
ただし、電磁気学および振動・波動学を学習していることを前提とします。
4. 光圧の理論
冒頭で述べた通り、光圧とは光と物質との相互作用により物質に働く力学的な力のことです。
ただし、本記事では物質を誘電特性が不変の剛体として取り扱います。
すなわち、物質内部での力学的な変形、誘電率の時間的な変化を無視します。
さらに、物質はもともと電荷を帯びていない中性の物質とします。
また、光は調和的であると仮定します。
すなわち、光は一定の周波数$\omega$は有しているとします。
さらに、光の波数ベクトルを$\boldsymbol{k}$とします。
以上の物理量を用いて、入射光の電場の実部$\boldsymbol{E}(\boldsymbol{r},t)$、磁束密度の実部$\boldsymbol{B}(\boldsymbol{r},t)$を次のように表します。
\boldsymbol{E}(\boldsymbol{r},t) = \text{Re}[\boldsymbol{E}(\boldsymbol{r},\boldsymbol{k})\exp(-i \omega t)]\tag{1}\\
\boldsymbol{B}(\boldsymbol{r},t) = \text{Re}[\boldsymbol{B}(\boldsymbol{r},\boldsymbol{k})\exp(-i \omega t)]\tag{2}
$\boldsymbol{E}(\boldsymbol{r},\boldsymbol{k})$および$\boldsymbol{B}(\boldsymbol{r},\boldsymbol{k})$は電場および磁束密度の複素振幅です。$\exp(-i \omega t)$の部分は光の時間振動の項を示します。
物質が光と相互作用すると、入射光の電場により物質は分極します。
ただし、この記事では簡易的に説明するために、入射光の電場に応答して生じる分極は電場に1次比例すると仮定します。
誘起分極により生じた物質内部の電荷密度を$\rho(\boldsymbol{r}, t)$、電流密度を$\boldsymbol{j}(\boldsymbol{r},t)$と置くと、物質に働くローレンツ力の総和の時間平均は次のようになります。
<\int d\boldsymbol{r} \bigl(\rho(\boldsymbol{r}, t)\boldsymbol{E}(\boldsymbol{r},t) + \boldsymbol{j}(\boldsymbol{r},t)\times\boldsymbol{B}(\boldsymbol{r},t) \bigr)>_T\tag{3}
積分範囲は物質内部の空間全範囲になります。
$<>_T$は$<>$内の物理量の時間平均を表します。
以下からの作業は、この物質に働くローレンツ力の総和の時間平均の式をひたすら変形していくことを考えます。
本記事の最後にわかりますが、光圧の理論式は(3)式を変形したものに相当します。
ただし、その変形過程が煩雑なので丁寧に式変形していきます。
戦略として、まず初めに(3)式を変形していくための下準備を行います。
電荷密度の電荷は分極により生じた分極電荷のみなので、分極の実部を$\boldsymbol{P}(\boldsymbol{r}, t)$とすると、電荷密度と分極の関係は以下の式で表されます。
\rho(\boldsymbol{r},t) = -\nabla \cdot \boldsymbol{P}(\boldsymbol{r}, t)\tag{4}
粗視化領域内に生成された電荷は誘起分極に付随した電荷に相当することを意味する式ですね。
次に物質内の電荷に関する連続の式を考えます。
物質内の粗視化領域の電荷密度の時間変化は以下のように表されます。
\frac{\partial \rho(\boldsymbol{r},t)}{\partial t} = -\nabla \cdot \boldsymbol{j}(\boldsymbol{r}, t)\tag{5}
こちらは粗視化領域内の電荷の時間変化は誘起分極に付随した電荷の移動によるものであることを意味する式ですね。
そして、(4)式と(5)式を連立すると以下の式が算出されます。
\frac{\partial \boldsymbol{P}(\boldsymbol{r}, t)}{\partial t} =\boldsymbol{j}(\boldsymbol{r}, t)\tag{6}
分極は電場に対して1次比例すると仮定したので、分極の実部$\boldsymbol{P}(\boldsymbol{r}, t)$と電場の実部$\boldsymbol{E}(\boldsymbol{r}, t)$は、
真空の誘電率$\varepsilon_0$および複素電気感受率$\boldsymbol{\chi}(\boldsymbol{r})$を用いて以下の関係で表されるとします。
\boldsymbol{P}(\boldsymbol{r}, t) = \text{Re}[\varepsilon_0 \boldsymbol{\chi}(\boldsymbol{r})\boldsymbol{E}(\boldsymbol{r},\boldsymbol{k})\exp(-i \omega t)]\tag{7}
ただし、複素電気感受率$\boldsymbol{\chi}$は2階のテンソル量です。
テンソルという言葉の響きは難しさを連想させてきますが、電気感受率に異方性のない物質を取り扱う場合は複素電気感受率$\boldsymbol{\chi}$をスカラー量ととらえても特に問題はありません。
また、複素電気感受率$\boldsymbol{\chi}$は時間的に変化しないと仮定しました。
(7)式から分極の複素振幅は以下のように表されます。
\boldsymbol{P}(\boldsymbol{r},\boldsymbol{k}) = \varepsilon_0 \boldsymbol{\chi}(\boldsymbol{r})\boldsymbol{E}(\boldsymbol{r},\boldsymbol{k})\tag{8}
さらに、マクスウェル方程式のファラデー-マクスウェルの式は電場の実部および磁束密度の実部を用いて以下のように表されます。
\frac{\partial \boldsymbol{B}(\boldsymbol{r}, t)}{\partial t} = -\nabla \times \boldsymbol{E}(\boldsymbol{r}, t)\tag{9}
ただし、(9)式の両辺にある演算子は線形演算子であるため、光の電場および磁束密度を複素表示に拡張してもファラデー-マクスウェルの式は成立します。
そこで、光の電場および磁束密度を複素表示に拡張すると以下のように表されます。
\begin{align}
\frac{\partial }{\partial t}(\boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i \omega t)) &= -\nabla \times (\boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i \omega t))\\
\Leftrightarrow(-i \omega) \boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i \omega t)) &= -\nabla \times (\boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i \omega t))\\
\Leftrightarrow i \omega \boldsymbol{B}(\boldsymbol{r},\boldsymbol{k}) &= \nabla \times \boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k})\tag{10}\\
\end{align}
以上で下準備が整いました。
次の戦略として、調達した式を(3)式のローレンツ力の総和の時間平均の式に代入して変形していきます。
(4)および(6)式を(3)式に代入すると次のように表されます。
\begin{align}
<\int d\boldsymbol{r} (\rho(\boldsymbol{r}, t)\boldsymbol{E}(\boldsymbol{r},t) + \boldsymbol{j}(\boldsymbol{r},t)\times\boldsymbol{B}(\boldsymbol{r},t) \bigr)>_T\\
= <\int d\boldsymbol{r} \Bigl(\bigl(-\nabla \cdot \boldsymbol{P}(\boldsymbol{r}, t)\bigl)\boldsymbol{E}(\boldsymbol{r}, t) + \frac{\partial \boldsymbol{P}(\boldsymbol{r}, t)}{\partial t} \times \boldsymbol{B}(\boldsymbol{r}, t)\Bigr)>_T\tag{11}\\
\end{align}
そして、(11)式の第1項および第2項を別々に変形していきます。
変形の過程で、次の関係式を頻繁に使用します。
$a$および$b$を実数として$c = a + bi$で与えられる複素数$c$の実部$a$は
a = \frac{1}{2}(c + c^*)
となりますね。
この関係式を使用して、これまで登場した複素表示の物理量とその実部を繋げながら変形していきます。
まずは第1項から
\begin{align}
<
\int d\boldsymbol{r}
\Bigl(
\bigl(
-\nabla \cdot \boldsymbol{P}(\boldsymbol{r},t)
\bigl)
\boldsymbol{E}(\boldsymbol{r}, t)
\Bigr)
>_T
&= <
\int d\boldsymbol{r}
\bigl(
-\nabla \cdot (\frac{1}{2}[
\boldsymbol{P}(\boldsymbol{r},\boldsymbol{k})\exp(-i\omega t)
+ \boldsymbol{P^*}(\boldsymbol{r},\boldsymbol{k})\exp(i\omega t)
])
\frac{1}{2}[
\boldsymbol{E}(\boldsymbol{r},\boldsymbol{k})\exp(-i\omega t)
+ \boldsymbol{E^*}(\boldsymbol{r},\boldsymbol{k})\exp(i\omega t)
]
\bigl)
>_T\\
&= -\frac{1}{2}
\int d\boldsymbol{r}
\Bigl(
\frac{1}{2}
\bigl[
(\nabla \cdot \boldsymbol{P}(\boldsymbol{r},\boldsymbol{k})) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
+(\nabla \cdot \boldsymbol{P^*}(\boldsymbol{r},\boldsymbol{k})) \boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k})
\bigr]
\Bigr)
\\
&= -\frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\bigl(
(\nabla \cdot \boldsymbol{P}(\boldsymbol{r},\boldsymbol{k})) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k}
\bigr)
\Bigr]\tag{12}
\end{align}
1番目の等式は、上述した複素表示の物理量とその実部を繋げる関係式を使用しました。
2番目の等式は、時間平均により時間振動項がある項を除去しました。これは時間平均を行う時間幅$T$が
\frac{2\pi}{\omega} << T
である場合に成り立つ近似で、本記事ではこの近似が成り立つと仮定します。
3番目の等式は、1番目の等式と同じく複素表示の物理量とその実部を繋げる関係式を使用しました。
(12)式を見ると積分の中に$\nabla \cdot \boldsymbol{P}(\boldsymbol{r},\boldsymbol{k})$の項があります。
入射光の電場により、物質から電荷が飛び出すような不安定な状況を仮定しない限り、
誘起分極は物質表面で0になるので、部分積分を利用して$\boldsymbol{P}(\boldsymbol{r} = surf,\boldsymbol{k}) = 0$を使用したいですね。
そこで、部分積分を利用して上式を変形していきます。
\begin{align}
-\frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\bigl(
(
\nabla \cdot \boldsymbol{P}(\boldsymbol{r},\boldsymbol{k})
) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
\Bigr]
&= -\frac{1}{2}
\text{Re}\Bigl[
\int dx_1dx_2dx_3
\bigl(
(
\frac{\partial}{\partial x_1} P_1(\boldsymbol{r},\boldsymbol{k})
+
\frac{\partial}{\partial x_2} P_2(\boldsymbol{r},\boldsymbol{k})
+
\frac{\partial}{\partial x_3} P_3(\boldsymbol{r},\boldsymbol{k})
)
\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k}
\bigr)
\Bigr]\\
&= -\frac{1}{2}
\text{Re}\Bigl[
\int dx_2dx_3
\bigl(
[
P_1(\boldsymbol{r},\boldsymbol{k}) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
]_{\boldsymbol{r} = surf}
- \int dx_1
P_1(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_1}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
+
\int dx_3dx_1
\bigl(
[
P_2(\boldsymbol{r},\boldsymbol{k}) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
]_{\boldsymbol{r} = surf}
- \int dx_2
P_2(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_2}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
+
\int dx_1dx_2
\bigl(
[
P_3(\boldsymbol{r},\boldsymbol{k}) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
]_{\boldsymbol{r} = surf}
- \int dx_3
P_3(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_3}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
\Bigr]\\
&= \frac{1}{2}
\text{Re}\Bigl[
\int dx_2dx_3
\bigl(
\int dx_1
P_1(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_1}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
+
\int dx_3dx_1
\bigl(
\int dx_2
P_2(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_2}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
+
\int dx_1dx_2
\bigl(
\int dx_3
P_3(\boldsymbol{r},\boldsymbol{k}) \frac{\partial}{\partial x_3}\boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
\Bigr]\\
&= \frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\bigl(
(
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \cdot \nabla
) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\bigr)
\Bigr]\tag{13}
\end{align}
ごちゃごちゃしていますが、部分積分を行い表面項を除去しようとしているだけです。
3番目の等式で表面項を除去しました。
これにて第1項の変形は終わります。
次に第2項を変形していきます。
\begin{align}
<\int d\boldsymbol{r} \Bigl(
\frac{\partial \boldsymbol{P}(\boldsymbol{r}, t)}{\partial t} \times \boldsymbol{B}(\boldsymbol{r}, t)
\Bigr)>_T
&= <
\int d\boldsymbol{r}
\Bigl(
\bigl(
\frac{\partial}{\partial t}
(
\frac{1}{2} [
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i\omega t)
+
\boldsymbol{P^*}(\boldsymbol{r}, \boldsymbol{k}) \exp(i\omega t)
]
)
\bigr)
\times
\bigl(
\frac{1}{2}[
\boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k}) \exp(-i\omega t)
+
\boldsymbol{B^*}(\boldsymbol{r}, \boldsymbol{k}) \exp(i\omega t)
]
\bigl)
\Bigr)
>_T\\
&= \frac{1}{2}
\int d\boldsymbol{r}
\Bigl(
\frac{(i\omega)}{2}
[
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \times (-\boldsymbol{B^*}(\boldsymbol{r}, \boldsymbol{k}))
+ \boldsymbol{P^*}(\boldsymbol{r}, \boldsymbol{k}) \times \boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k})
]
\Bigr)\\
&= \frac{1}{2}
\int d\boldsymbol{r}
\Bigl(
\frac{1}{2}
[
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \times (i \omega \boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k}))^*
+ \boldsymbol{P^*}(\boldsymbol{r}, \boldsymbol{k}) \times (i \omega \boldsymbol{B}(\boldsymbol{r}, \boldsymbol{k}))
]
\Bigr)\\
&= \frac{1}{2}
\int d\boldsymbol{r}
\Bigl(
\frac{1}{2}
[
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \times (\nabla \times \boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}))^*
+ \boldsymbol{P^*}(\boldsymbol{r}, \boldsymbol{k}) \times (\nabla \times \boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}))
]
\Bigr)\\
&= \frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\Bigl(
\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \times (\nabla \times \boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}))^*
\Bigr)
\Bigr]\\
&= \frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\Bigl(
(\nabla \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})) \cdot \boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k})
-(\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k}) \cdot \nabla) \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})
\Bigr)
\Bigr]\tag{14}
\end{align}
1番目の等式は複素表示の物理量とその実部を繋げる関係式を使用しました。
4番目の等式は(10)式を使用しました。
5番目の等式も同じく複素表示の物理量とその実部を繋げる関係式を使用しました。
6番目の等式はベクトル恒等式$\boldsymbol{A}\times (\nabla \times \boldsymbol{B})=(\nabla \boldsymbol{B}) \cdot \boldsymbol{A}-(\boldsymbol{A}\cdot \nabla)\boldsymbol{B}$を使用しました。
これにて第2項の変形は終わりです。
第1項および第2項の変形が終えたので、2つを統合して考えると
ローレンツ力の総和の時間平均は次のようになります。
\begin{align}
<\int d\boldsymbol{r} \Bigl(\bigl(-\nabla \cdot \boldsymbol{P}(\boldsymbol{r}, t)\bigl)\boldsymbol{E}(\boldsymbol{r}, t) + \frac{\partial \boldsymbol{P}(\boldsymbol{r}, t)}{\partial t} \times \boldsymbol{B}(\boldsymbol{r}, t)\Bigr)>_T\\
=\frac{1}{2}
\text{Re}\Bigl[
\int d\boldsymbol{r}
\Bigl(
(\nabla \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})) \cdot \boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k})
\Bigr)
\Bigr]\tag{15}
\end{align}
長い式変形になりましたが、これが光圧の一般的な表式になります。
この式から、光の電場が大きくなる空間方向に対して、物質に力学的な力が働き($=\nabla \boldsymbol{E^*}(\boldsymbol{r}, \boldsymbol{k})\cdot $の部分)、実際に物質が感受する力学的な力の度合いは、分極の大きさに比例する($=\boldsymbol{P}(\boldsymbol{r}, \boldsymbol{k})$の部分)ということが直観的に判断できます。
光ピンセットの原理の説明で、レンズを用いてレーザー光を集光していたのは空間に電場勾配を生成するためで、集光点に微小物質が捕捉されるのは、集光点に近いほど空間に生成された電場が大きいから、ということだったのです。
5. 厳密には
これまで、説明をわかりやすくするために電場の複素振幅$\boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k})$をもっぱら入射光に由来するものとして取り扱って来ましたが、厳密には入射光の電場のほかに物質内部で生じた分極により発生する電場も考慮する必要があります。
すなわち、入射光の電場の複素振幅を$\boldsymbol{E_{\text{inc}}}(\boldsymbol{r}, \boldsymbol{k})$、分極により生じた電場の複素振幅を$\boldsymbol{E_{\text{scatt}}}(\boldsymbol{r}, \boldsymbol{k})$とすると、電場の複素振幅$\boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k})$は次のように表されるべきです。
\boldsymbol{E}(\boldsymbol{r}, \boldsymbol{k}) = \boldsymbol{E_{\text{inc}}}(\boldsymbol{r}, \boldsymbol{k}) + \boldsymbol{E_{\text{scatt}}}(\boldsymbol{r}, \boldsymbol{k})
この表式を光圧の理論式に代入することで、物質に働く光圧を計算できます。
そして、実際に光圧を使用したシミュレーションにより物質が捕捉される様子を観察することは可能ですが、計算量は膨大です。
例えば、光と相互作用する物質を$N$コの粗視化領域に分けるとして、一つの粗視化領域における電場は
(入射光の電場) + (周囲の粗視化領域の分極により生じた電場)の影響を考える必要があるため、単純計算で3次元空間の場合、1つの物質に働く光圧の計算量は$O(N(3 + 3(N-1))=O(N^2)$のオーダーになります。
そして、空間に存在している捕捉対象の物質が1つではなく、複数コ存在している状況では物質間の電場の影響も計算する必要があるため計算量はさらに膨大化します。
このように物質を複数の粗視化領域に分割し、入射光の電場に加えて各々の粗視化領域の分極から生じる電場を考慮して電磁場解析を行う手法は離散双極子近似と呼ばれています。
離散双極子近似を用いて物質に働く光圧を計算し、実際に物質が電場勾配に沿って捕捉される様子を観察することは面白いですが、上述したように計算量が膨大になり、おそらく投稿主のPC環境ではスペック的に大問題です。
そこで、最終目的の光圧を用いて物質を捕捉するシミュレーションでは、本記事で導出した理論式をそのまま取り扱わず、近似して取り扱います。
近似内容については、他の記事にて紹介します。
6. さいごに
いかかでしたでしょうか？
直観的に理解しがたい光圧の正体は、光と相互作用することにより生じた振動分極が感じる力学的な力のことだったのですね。
次回以降の記事で、今回導出した光圧の理論式の近似法を紹介し、最終的には光圧を用いたシミュレーションを行いたいと思います。
7. 参考文献
石原一・芦田昌明 著, 「光圧 ―物質制御のための新しい光利用」, 朝倉書店, 2021.