## Doppler averaging

Motivation. The Doppler broadening of EIT lines are of order $\rm 100\,MHz$, which is much larger than the ~MHz EIT line width. Thus, the Doppler ensemble averaging must be considered in the numerical simulations.  

We wish to compute the Doppler averaging 
$$
\rho^D = \int_{-\infty}^\infty \rho(v)\mathcal{N}(v){\rm d}v
$$
where $\mathcal{N}(v)$ is some one-dimensional Gaussian distribution. 



Some useful conclusions: The Dawson function 
$$
D(x) = e^{-x^2}\int_0^x e^{t^2}{\rm d}t
$$
The erf function 
$$
{\rm erf}(x) := \frac{2}{\sqrt{\pi }}\int_0^x e^{-t^2}{\rm d}t
$$
The erfi function 
$$
{\rm erfi}(x) = -i {\rm erf}(ix)=\frac{2}{\sqrt{\pi}} \int_0^x e^{t^2}{\rm d}t
$$
Thus, we have ${\rm erfi}(x) = 2\pi^{-1/2}e^{x^2}D(x) $

The convolution formula for a single pole and a Gaussian distribution: 
$$
\int_{-\infty}^{\infty} \frac{1}{y-x}e^{-x^2}{\rm d}x = {2}\sqrt{\pi} D(y)
$$
Now we wish to compute 
$$
\int_{-\infty}^{\infty}\frac{1}{1+\lambda x}p(x){\rm d}x=-\frac{\sqrt{2}}{\lambda}D(-\frac{1}{\lambda\sqrt{2}})
$$

A direct numerical computation on this will do the job. 



We wish to evaluate the integral 
$$
\int_{-\infty}^{\infty} \frac{1}{z-x} p(x) {\rm d}x , \, z\in\mathbb{C}, \Im[z] > 0
$$


## The eigenspace $V_0$ 

Suppose we have two equations: 
$$
{\bf A}_0 {\bf x}_0 = {\bf 0}, \\ 
({\bf A}_0+\frac{v}{\sigma_v}{\bf A}_v){\bf x}_v = {\bf 0}.
$$
We wish to establish some relationship like this: 
$$
{\bf x}_v = \left({\bf I}_{16} + \frac{v}{\sigma_v} {\bf A}_0^+ {\bf A}_v\right)^{-1} {\bf x}_0.
$$
To establish this relationship, we need to mention the property of the operator ${\bf A}_0^+{\bf A}_0$. Let $V_\lambda({\bf A}_0)$ be the $\lambda$-eigenspace of matrix ${\bf A}_0$. Since ${\bf x}_v$ differs from ${\bf x}_0$, we consider taking its projection onto the eigenspace of ${\bf A}_0$. 

Let ${\bf A}_0={\bf V} {\bf \Lambda}{\bf V}^{-1}$ be its eigen decomposition with $\lambda_0=0$. Then ${\bf A}_0^+$ can be easily written as ${\bf A}_0^+ = {\bf V}{\rm diag}(0, \lambda_2^{-1}, \cdots, \lambda_n^{-1}){\bf V}^{-1}$. Let ${\bf x}_v$ be the solution to the second equation. Then: 
$$
({\bf I} - {\bf A}_0^+{\bf A}_0){\bf x}_v = {\bf V}{\rm diag}(1,0,\cdots,0){\bf V}^{-1}{\bf x}_v \propto {\bf x}_0
$$
This means that ${\bf x}_v$ satisfies 
$$
{\bf x}_v - {\bf A}_0^+{\bf A}_0 {\bf x}_v = c {\bf x}_0
$$
We can check that, in fact, under proper choice of $\bf V$, the coefficient $c=1$. Using the fact that ${\bf A}_0{\bf x}_v = -\frac{v}{\sigma_v}{\bf A}_v {\bf x}_v$, we get 
$$
{\bf x}_v + \frac{v}{\sigma_v}{\bf A}_0^+ {\bf A}_v{\bf x}_v = {\bf x}_0, 
$$
which yields 
$$
{\bf x}_v = \left({\bf I} + \frac{v}{\sigma_v}{\bf A}_0^+{\bf A}_v\right)^{-1}{\bf x}_0. 
$$
An eigen-decomposition on ${\bf A}_0^+{\bf A}_v = {\bf S}{\rm diag}(\mu_1, \cdots, \mu_n){\bf S}^{-1}$ yields 
$$
{\bf x}_v = {\bf S}{\rm diag}\left(\frac{1}{1+\mu_1 v/\sigma_v}, \cdots, \frac{1}{1+\mu_n v/\sigma_v}\right){\bf S}^{-1} {\bf x}_0.
$$
The Doppler average can be written as $\bar{\bf x}=\int_{\mathbb{R}} {\bf x}_{v_x}\mathcal{N}_T(v_x|0,\sigma_x^2){\rm d}v_x$, which can be computed by evaluating the integral  
$$
I(z):=\int_{-\infty}^\infty \frac{1}{1+zx}\mathcal{N}(x|0,1){\rm d}x
$$
A similar integral takes the form of a Hilbert transform of a Gaussian, which is of the shape 
$$
J(z) = \int_{-\infty}^\infty \frac{1}{z-\xi}\mathcal{N}(\xi|0,1){\rm d}\xi,
$$
where the integral $J(z)$ satisfies, for $\Im[z]\neq0$, we have 
$$
J(z) = \frac{1}{\sqrt{2\pi}} e^{-z^2/2}\left(\pi \cdot {\rm erfi}(\frac{z}{\sqrt{2}}) - \log (-1/z)-\log z\right)
$$
For the integral $I(z)$, it is easy to see that $I(z) = -z^{-1}J(-z^{-1})$. 

We recall that 
$$
\log(-1/z) + \log(z) 
$$
evaluates to $i\pi$ when $\Im[z]>0$ and $-i\pi$ when $\Im[z]<0$. Thus, we may not need to explicitly evaluate the complex logarithm function here. This change can be explained by adding a residual of the pole when shifting the P.V. integral upwards into the upper half-plane. The residual of the simple pole is -1, so the integral $J(z)$ should increase by $i\pi$ times the residual, which is $-i\pi$. 

A simple expression will look like 
$$
J(z) = \mathcal{N}(z|0,1)\left(\pi\cdot {\rm erfi}\left(z/\sqrt{2}\right) -i\pi\cdot {\rm sgn}\Im[z]\right).
$$
where the sign function ${\rm sgn}(\cdot)$ returns +1 when its argument is positive, -1 when its argument is negative, and 0 when its argument is zero. We can easily see that $J(z)$ is an odd function with $J(-z)=-J(z)$, and that $J(z)$ is analytic on the upper half-plane. 



Now we wish to evaluate the P.V. integral when $z\in\mathbb{R}$. 

We have 
$$
\sqrt{2\pi} J_\epsilon(z) = e^{-z^2/2}\int_\epsilon^\infty \frac{e^{zt}-e^{-zt}}{t}e^{-t^2/2}{\rm d}t
$$
We can see that $J_\epsilon \to J_0=J$ do converge at $\epsilon\to 0^+$, and $J(0) = 0$. 

Let 
$$
K(z) = \int_0^\infty \frac{e^{zt}-e^{-zt}}{t}e^{-t^2/2}{\rm d}t,
$$
 then we have 
$$
K'(z) = \int_0^\infty (e^{zt}+e^{-zt})e^{-t^2/2}{\rm d}t = \sqrt{2\pi}e^{z^2/2}
$$
Combining with $K(0)=0$, we have 
$$
K(z) = \sqrt{2\pi} \int_0^z e^{x^2/2}{\rm d}x = \pi {\rm erfi}(z/\sqrt{2})
$$
Thus, we have, for $z\in\mathbb{R}$, 
$$
J(z) = \frac{1}{\sqrt{2\pi}}e^{-z^2/2} \cdot \pi\cdot {\rm erfi}(z/\sqrt{2}).
$$
**Asymptotics**. We wish to evaluate the asymptotic behavior of $J(z)$ when $z$ is large. As $z\to\infty$ in the upper half-plane, we expect that $J(z)\sim 1/z$.  An asymptotic expansion requires that 
$$
J(z) \sim \frac{1}{z} + \frac{c_2}{z^2} + \frac{c_3}{z^3} + \cdots
$$
Mathematica gives me that 
$$
\frac{\pi}{2}{\rm erfi}(z/\sqrt{2}) \sim -i\sqrt{\frac{\pi}{2}} + e^{z^2/2}(z^{-1}+z^{-3} + 3z^{-5} + O(z^{-6})),\quad z\to\infty
$$
and I'm not sure whether this is correct. 





**Deprecated**. We may evaluate this with the Taylor expansion of 
$$
\frac{1}{1+zx} = \sum_{n=0}^\infty (-x)^n z^n
$$
We know that for a Gaussian variable $X\sim \mathcal{N}(0,1)$, the even moment $\mu_{2k}$ is evaluated to be 
$$
\mu_{2k} = \frac{2}{\sqrt{2\pi}}\int_0^\infty x^{2k}e^{-x^2/2}{\rm d}x = (2k-1)!! 
$$
and the odd moments all evaluates to zero. This does not work, since this will get the Taylor series to diverge. 



## Algorithm acceleration

| Method                                   | Numerical Averaging | Analytical Integration |
| ---------------------------------------- | ------------------- | ---------------------- |
| EIT-AT spectrum, 800 points, no parallel | 33s                 | 0.3s                   |
| EIT-AT spectrum, 800 points, process parallel | 5.8s                | 0.01s                  |
|                                          |                     |                        |

No significant difference observed between the analytical Doppler-average solver and the traditional numerical integration solver. 



[R1]  O. Nagib and T. G. Walker, "Fast and accurate method for Doppler averaging of Rydberg EIT signals," arXiv preprint arXiv: 2501.02141, Jan. 2025. 



## The Dynamic Case 

The $s$-domain transfer function from $H_{k\ell}$ to $\Delta\rho_{21}$ is denoted by $T_{k\ell}(s)$, which is evaluated from 
$$
{T_{k\ell}}(s) = \left[{\bf Q}\begin{bmatrix}{\bf 0} \\ {\bf I}_{15}\end{bmatrix} \left(s{\bf I}_{15}-{\bf C}_0\right)^{-1} {\bf F}_{k\ell}\bar{\bf z}_0\right]_2
$$
 where the zero-velocity steady-state $\bar{\bf z}_0$ is determined by the linear equation ${\bf C}_0\bar{\bf z}_0+(1/2){\bf w}_0={\bf 0}$. 

When considering Doppler velocity $v$, the matrices ${\bf C}_0$ is updated to be 
$$
{\bf C}_0 + \frac{v}{\sigma_v}{\bf C}_v
$$
and the vector ${\bf w}_0$ is similarly updated to be ${\bf w}_0 + (v/\sigma_v){\bf w}_v$. We can prove that ${\bf w}_v={\bf 0}$. So this part is slightly simplified. 

For the velocity-$v$ expression, we have 
$$
{T_{k\ell}}(s,v) = \left[{\bf Q}\begin{bmatrix}{\bf 0} \\ {\bf I}_{15}\end{bmatrix} \left(s{\bf I}_{15}-{\bf C}_0 - \frac{v}{\sigma_v}{\bf C}_v\right)^{-1} {\bf F}_{k\ell}\bar{\bf z}_v\right]_2,
$$
and the Doppler-averaged $T_{k\ell}$ is given by 
$$
T_{k\ell}(s) := \int_{-\infty}^\infty T_{k\ell}(s,v)\mathcal{N}(v|0, \sigma_v^2){\rm d}v
$$
where $\sigma_v^2=k_{\rm B}T/m_{\rm Cs}$.  

For the part of $(s{\bf I}_{15}-{\bf C}_0 - (v/\sigma_v){\bf C}_v)^{-1}$, for every fixed $s=i2\pi f$, this can be treated by the following transform 
$$
(s{\bf I}_{15}-{\bf C}_0 - (v/\sigma_v){\bf C}_v)^{-1} = \left({\bf I}-\frac{v}{\sigma_v}(s{\bf I}-{\bf C}_0)^{-1}{\bf C}_v \right)^{-1} (s{\bf I}-{\bf C}_0)^{-1}
$$
Following the same approach, we can assume an eigenvalue decomposition on $(s{\bf I}-{\bf C}_0)^{-1}{\bf C}_v$ as 
$$
(s{\bf I}-{\bf C}_0)^{-1}{\bf C}_v = {\bf S}(s){\rm diag}({\lambda}_1(s), \cdots, {\lambda}_n(s)){\bf S}^{-1}(s) = \sum_{m=1}^{15} \lambda_m [{\bf S}(s)]_{:,m}[{\bf S}^{-1}(s)]_{m,:}
$$
Note that the dependence of $v$ is only through $v/\sigma_v$.

and similarly, $-{\bf C}_0^{-1}{\bf C}_v$ can be decomposed as ${\bf S}(0){\rm diag}(\cdots){\bf S}^{-1}(0)$. With this formulation, the transfer function $T_{k\ell}(s)$ can be computed as 
$$
T_{k\ell}(s) = {\bf Q}\begin{bmatrix}{\bf 0} \\ {\bf I}_{15}\end{bmatrix} \sum_{m,n=1}^{15} \frac{1}{1-\lambda_m(s) \cdot v/\sigma_v}\cdot\frac{1}{1-\lambda_n(0)\cdot v/\sigma_v}{\bf s}_m(s){\bf t}_m^{\rm T}(s) (s{\bf I}-{\bf C}_0)^{-1} {\bf F}_{k\ell} {\bf s}_n(0) {\bf t}^{\rm T}_n(0)\bar{\bf z}_0
$$
where ${\bf S}(s) = [{\bf s}_1(s), \cdots, {\bf s}_{15}(s)]$, and ${\bf S}^{-1}(s) = [{\bf t}_1(s), \cdots, {\bf t}_{15}(s)]^{\rm T}$, and we notice the formula 
$$
\bar{\bf z}_v = -\frac{1}{2}({\bf C}_0 + \frac{v}{\sigma_v}{\bf C}_v)^{-1}{\bf w}_0 = ({\bf I}+ \frac{v}{\sigma_v}{\bf C}_0^{-1}{\bf C}_v)^{-1}\bar{\bf z}_0
$$
Thus, the vector $\bar{\bf z}_v$ is expressed as a function of $v$ as 
$$
\bar{\bf z}_v = {\bf S}(0) {\rm diag}\left(\frac{1}{1-\lambda_1(0)v/\sigma_v}, \cdots, \frac{1}{1-\lambda_{15}(0)v/\sigma_v}\right) {\bf S}^{-1}(0)\bar{\bf z}_0
$$


The final step is to average over a standard Gaussian distribution on $v/\sigma_v$. Recall that $I(-z) = \mathbb{E}[1/(1-zX)]$ where $X$ is standard Gaussian distributed. Thus, we have 
$$
\mathbb{E}\left[\frac{1}{(1-z_1X)(1-z_2X)}\right] = A \cdot I(-z_1) + B \cdot I(-z_2),
$$
Recall that $I(-z) = z^{-1}J(z^{-1})$. 

where $A=z_1/(z_1-z_2)$ and $B=z_2/(z_1-z_2)$. Or equivalently, 
$$
\frac{1}{z_1-z_2}(J(z_1^{-1}) - J(z_2^{-1}))
$$
For the case where $z_1\approx z_2$, we may differentiate the function $zI(z)$ to get the correct result. If any of the $z_i$ is zero, then the function will degrade to $I(z)$. 



I've succeeded to do this thing, by the defined special functions!!! 





## Response Writing 

Technical updates: 

- Doppler averaging (written to the manuscript)
- transit dephasing 
- laser noise







