# Coherence Protection of Electron Spin in Earth-field Range by All-optical
Dynamic Decoupling

Peiyu Yang 1 Guzhi Bao 2<EMAIL_ADDRESS>L. Q. Chen 1,4
<EMAIL_ADDRESS>Weiping Zhang 2,3,4<EMAIL_ADDRESS>1 State Key
Laboratory of Precision Spectroscopy, Quantum Insitute for Light and Atoms,
School of Physics and Electronic Science, East China Normal University,
Shanghai 200241, China.
2 School of Physics and Astronomy, and Tsung-Dao Lee Institute, Shanghai Jiao
Tong University, Shanghai 200240, China.
3 Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan, Shanxi 030006, China.
4 Shanghai Research Center for Quantum Sciences, Shanghai 201315, China.

###### Abstract

In recent years, unshielded atomic systems have been attracting researchers’
attention, in which decoherence is one of the major problems, especially for
high precision measurements. The nonlinear Zeeman effect and magnetic field
gradient are the main decoherence sources of atomic electron spin in Earth-
field range. Here, we propose a method to cancel out the two dominant
broadening effects simultaneously by an all-optical dynamic decoupling
approach based on Raman scattering in the ${}^{87}{\rm Rb}$ Zeeman sublevels.
By adjusting the parameters of the Raman lasers, we realize spin control along
an arbitrary direction. We analyze the state evolution of atomic spin under
the Raman light control sequence in detail. The results show that both the
nonlinear Zeeman effect and magnetic field gradient can be significantly
suppressed.

##  I Introduction

Atomic electron spin is a basic information carrier in quantum science
kominis2003subfemtotesla ; fixler2007atom ; dunning2014composite ;
chen2015atom ; bao2018suppression ; bao2018all ; fu2020sensitive ;
saywell2020optimal ; jiang2020interference ; dudin2013light ;
chen2010observation ; rui2015operating ; guo2019high ; brennen2005criteria ;
wu2016realization . Coherence protection is of major importance for precision
measurements and for scalable quantum information networks. Short coherence
time limits the sensitivity of atomic sensorskominis2003subfemtotesla ;
fixler2007atom ; dunning2014composite ; chen2015atom ; bao2018suppression ;
bao2018all ; fu2020sensitive ; saywell2020optimal  and quantum memory
timedudin2013light ; rui2015operating ; guo2019high . Such as, the sensitivity
of atomic magnetometer is described as $\delta
B=1/(\gamma\sqrt{nT_{2}Vt})$kominis2003subfemtotesla , in which $T_{2}$
represents the coherence time. Longer coherence time indicates better
sensitivity. Usually, high-performance atomic systems are protected by
magnetic shields kominis2003subfemtotesla ; chen2010observation ; chen2015atom
; guo2019high ; dudin2013light . However, practical application often requires
that atomic systems work in challenging conditions where background magnetic
fields with large amplitudes or large gradients exist fu2020sensitive . When
an atomic sample operates in the Earth’s magnetic field, the nonlinear Zeeman
(NLZ) effect and magnetic field gradient are major decoherence sources of
atomic electron spins acosta2006nonlinear ; seltzer2007synchronous ;
jensen2009cancellation ; fu2020sensitive ; bao2018suppression ; bao2018all .
These decoherence effects broaden the linewidth of magnetic resonance and
decrease the signal intensity, which leads to the drops of the sensitivity
quadratically to the coherence time. Furthermore, the asymmetries of magnetic-
resonance signal due to the NLZ splitting is particularly troublesome in
airborne and marine systems alexandrov2003recent . Several ways have been
reported to cancel the NLZ effect, e.g., double-modulated synchronous optical
pumping seltzer2007synchronous , high-order polarization moments
acosta2006nonlinear , tensor light shift effects jensen2009cancellation  and
spin locking bao2018suppression . Among them, spin locking is robust against
orientation changes bao2018suppression . However, spin locking brings power
broadening to the magnetic resonance, limiting the sensitivity bao2018all . To
achieve a narrow linewidth, the power of spin locking needs to be finely
controlled to reach a balance between spin locking and power broadening. With
the exisitence of a tilted angle between the sensor and the leading field, the
microwave fields and the atomic spins are no longer rotating in the same plane
which decreases the effect of spin locking.

Another choice is dynamic decoupling (DD) dudin2013light ; rui2015operating ;
antonijevic2003refocussing ; de2010universal ; naydenov2011dynamical ;
shaniv2019quadrupole ; pedernales2020motional ; Miaoeabc5186 ; hahn1950spin ;
barthel2010interlaced ; ahmed2013robustness ; pokharel2018demonstration ;
upadhyay2020ultralong , which is a quantum control technique used to extend
the atomic or nuclear coherence time by periodic sequences of instantaneous
modulation pulses and has been widely used to cancel quadrupolar energy shifts
antonijevic2003refocussing ; shaniv2019quadrupole  and inhomogeneous
broadenings hahn1950spin ; de2010universal ; rui2015operating . Conventional
DD is realized by pulsed microwave fields dudin2013light ;
antonijevic2003refocussing ; de2010universal ; naydenov2011dynamical ;
shaniv2019quadrupole ; pedernales2020motional ; Miaoeabc5186 . However, global
applied microwave fields might lead to crosstalk between closely located
samples bao2018all ; bison2009room ; lembke2014optical ; borna201720 .
Additionally, in an atom interferometer, the microwave fields cannot generate
a recoil velocity to separate the atomic energy levels fixler2007atom ;
dunning2014composite ; saywell2020optimal . Moreover, the microwave fields
cannot be applied to a faraway atomic sample directly, which limits the
applicability of this technique to remote control bao2018all ;
bustos2018remote ; patton2012remotely . An alternative solution is needed to
overcome these drawbacks above.

We construct an all-optical pulsed control field to achieve spin rotation by
Raman scattering. Raman scattering is a simple way to couple atomic spins with
only two optical fields fixler2007atom ; chen2010observation ;
dunning2014composite ; rui2015operating ; chen2015atom ; guo2019high ;
saywell2020optimal . Here, we manage to rotate the spin state by an arbitrary
angle along any direction by adjusting the Rabi frequencies, detunings,
propagating directions and pulse lengths of Raman lasers. By inserting Hahn
echo $\tau-180^{\circ}y-\tau$ into quadrupolar echo $\tau-90^{\circ}y-\tau$,
we design an all-optical sequence to cancel the NLZ effect and magnetic field
gradient simultaneously in the Earth-field range. Since all the control fields
are along $\hat{y}$ axis, even with a large tilted angle, DD works
successfully. To evaluate the performance of this all-optical scheme, we
exhibit the state evolution process and analyze the fidelity of the atomic
spin state with and without DD. The results show that the fidelity up to
$99.99\%$ with DD is much better than the free evolution case, indicating that
the coherence of the atomic spin state is well maintained by DD.

In this work, we propose and design an all-optical DD sequence by Raman
process to prolong the coherence time of atomic system in the Earth-field
range. The article is organised as follows. In Sec. II we study the underlying
dephasing mechanism limiting the coherence time. In Sec. III we give an
effective rotation Hamiltonian by Raman process to manipulate the atomic spin
in an all-optical way. In Sec. IV we design an all-optical DD sequence and
analyze the DD performance. In the last section, we draw some conclusions and
discussions.

##  II Dephasing mechanism of the electron spin

The dephasing effects of the electron spin contain at least six components:
spin-exchange relaxation, spin-destruction relaxation, wall-collision
relaxation, optical power broadening, the NLZ effect and magnetic field
gradient broadening. Among them, the NLZ effect and magnetic field gradient
are dominant in most cases.

The eigenvalues for the ground state $5S_{1/2}$ of ${}^{87}\mathop{\rm Rb}$ in
an external magnetic field B are given by the Breit-Rabi formula
breit1931measurement ; auzinsh2010optically  and illustrated in Fig. 1. The
energy splittings between adjacent Zeeman sublevels are

$E_{m_{F}+1}-E_{m_{F}}=\pm\frac{\mu_{B}B}{2}\mp\hbar\omega(2m_{F}+1),$ (1)

where $m_{F}$ is the magnetic quantum number, $\mu_{B}$ is the Bohr magneton,
B is the magnetic field intensity, $\omega=(\mu_{B}B)^{2}/(4\hbar\delta)$ is
the quantum-beat revival frequency caused by the NLZ effect, and $\delta$ is
the hyperfine-structure energy splitting; $\pm$ refer to the $F=I\pm 1/2$
hyperfine components, and $I$ is the nuclear spin. In the Earth-field range
(50 $\mu$T), the revival frequency for ${}^{87}\mathop{\rm Rb}$ is
$\omega=2\pi\times 17.9$ Hz, which is comparable to the width of magnetic
resonance. The right side of Fig. 1 gives the split magnetic-resonance signal
in the presence of the NLZ effect.

Figure 1:  Hyperfine structure of the $5S_{1/2}\,F=1$ and $F=2$ state
manifolds of the ${}^{87}\mathop{\rm Rb}$ atom in an external magnetic field.
The right side shows the split magnetic-resonance signal for the $F=2$ states
(up) and $F=1$ states (down). The signals are fitted with Lorentzian peaks
arising due to the NLZ effect. $\omega_{1}=\Omega_{L}-2\omega$,
$\omega_{2}=\Omega_{L}-\omega$, $\omega_{3}=\Omega_{L}+\omega$ and
$\omega_{4}=\Omega_{L}+2\omega$ are the quantum beating frequencies of $F=2$.
$\omega_{5}=\Omega_{L}+\omega$ and $\omega_{6}=\Omega_{L}-\omega$ are the
quantum beating frequencies of $F=1$.

When atomic systems work in an open environment, without the protection of
magnetic shields, background magnetic fields have large gradients. Due to the
inhomogeneous distribution of the external magnetic field, atoms located in
different positions have different Larmor precession frequencies, which leads
to dephasing of the atomic electron spin.

With consideration of the two main dephasing effects (NLZ effect and magnetic
field gradient), the total Hamiltonian of our system on the basis of the
Zeeman sublevels with a magnetic field along the $\hat{z}$ axis is

$H=\hbar(\Omega_{L}J_{z}+\Delta\Omega_{L}J_{z}+\omega J_{z}^{2}),$ (2)

where $\hbar$ is Plank’s constant, $J_{z}$ is the angular momentum along the
$\hat{z}$ axis, $\Omega_{L}$ is the Larmor precession frequency,
$\Delta\Omega_{L}$ is the inhomogeneous frequency broadening caused by the
magnetic field gradient, and $\omega$ is the revival frequency due to the NLZ
effect. Our goal is to cancel out the NLZ effect and magnetic field gradient
simultaneously.

##  III Spin control by Raman processes

A dynamic decoupling experiment often requires fast spin rotations by certain
degrees in the $\hat{x}$ or $\hat{y}$ direction, denoted as $\Phi_{x}$ or
$\Phi_{y}$, where $\Phi$ is the angle of rotation. Raman scattering is a
typical method for manipulating the atomic electron spin. For the all-optical
scheme, we construct an equivalent $J_{x}$ or $J_{y}$ Hamiltonian by quantum
interference of Raman channels.

Figure 2: Transition channels of the ${}^{87}\mathop{\rm Rb}$ D1 line from the
ground state $5S_{1/2},F=1$ Zeeman sublevels to excited energy levels
$5P_{1/2},F^{\prime}=1$ and $F^{\prime}=2$. $\delta$ is the hyperfine energy
level splitting. $\Delta$ denotes Raman detuning. $m_{F}$ is the magnetic
quantum number. Solid lines represent $\sigma_{+}$ transitions. Dashed lines
represent $\sigma_{-}$ transitions. Dot-dashed lines represent $\pi$
transitions.

Here, we consider the D1 transition of ${}^{87}\mathop{\rm Rb}$ from
$5S_{1/2},F=1$ to $5P_{1/2},F^{\prime}=1$ and $F^{\prime}=2$. We apply laser
lights with three polarizations, i.e., $\sigma_{+}$, $\pi$ and $\sigma_{-}$,
the Rabi frequencies of which are $\Omega_{1},\Omega_{2}$ and $\Omega_{3}$,
respectively. The propagating directions of the Raman lasers are all along the
$\hat{x}$ axis. Possible transition channels are depicted in Fig. 2. The
hyperfine energy level splitting of the excited state $\delta$ is 800 MHz and
the natural linewidth of the excited state is around 6 MHz. When the Raman
detuning $\Delta$ is far enough, the effective Hamiltonian $H_{\rm eff}$ on
the basis of the Zeeman sublevels can be derived from the method in the
Appendix A. To achieve a $J_{x}$-rotation Hamiltonian, we assume
$\Omega_{3}=\Omega_{1}$ and that $\pi$ polarized light has a $\pi/2$ phase
difference with circularly polarized lights $\sigma_{+}$, $\sigma_{-}$. We can
finally write $H_{\rm eff}$ as

$\displaystyle H_{\rm eff}$
$\displaystyle=\alpha_{1}\left(\begin{array}[]{ccc}0&1&0\\\ 1&0&1\\\ 0&1&0\\\
\end{array}\right)+\alpha_{2}\left(\begin{array}[]{ccc}0&0&1\\\ 0&0&0\\\
1&0&0\\\ \end{array}\right)$ (3)
$\displaystyle+\beta_{1}\left(\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\ 0&0&1\\\
\end{array}\right)+\beta_{2}\left(\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&1\\\ \end{array}\right),$

where $\alpha_{1},\alpha_{2}$, $\beta_{1}$ and $\beta_{2}$ are the coupling
coefficients, with

$\displaystyle\alpha_{1}$
$\displaystyle=\frac{(7\delta+4\Delta)\Omega_{1}\Omega_{2}}{6(3\delta-4\Delta)(5\delta+4\Delta)},$
(4a) $\displaystyle\alpha_{2}$
$\displaystyle=-\frac{\delta\Omega_{1}^{2}}{3(3\delta-4\Delta)(5\delta+4\Delta)},$
(4b) $\displaystyle\beta_{1}$
$\displaystyle=-\frac{2(3\delta+4\Delta)\Omega_{1}^{2}+(5\delta+4\Delta)\Omega_{2}^{2}}{6(3\delta-4\Delta)(5\delta+4\Delta)},$
(4c) $\displaystyle\beta_{2}$
$\displaystyle=-\frac{\delta_{1}(\Omega_{1}-\Omega_{2})(\Omega_{1}+\Omega_{2})}{3(3\delta-4\Delta)(5\delta+4\Delta)},$
(4d)

where $\delta$ is the hyperfine energy level splitting and $\Delta$ is the
Raman detuning.

The first term in Eq. (3) containing $\alpha_{1}$ denotes the $\Delta m_{F}=1$
transitions, which is a typical $J_{x}$ form for the $F=1$ system. The second
term containing $\alpha_{2}$ represents transitions between the two outermost
Zeeman energy levels $m_{F}=-1$ and $m_{F}=1$. The third term containing
$\beta_{1}$ gives a scalar shift independent of the Zeeman sublevels. The
scalar shift results in a global phase factor relative to the wavefunction,
which has no influence on the density matrix. The last term proportional to
$\beta_{2}$ is a tensor shift similar to the NLZ shift. To obtain a clean
$J_{x}$ Hamiltonian, we need to adjust the detunings and Rabi frequencies of
the Raman fields to eliminate the second term and the last term.

Figure 3: (a) The dependence of the fidelity on $\Omega_{1}/\Omega_{2}$ and
$\Delta/\delta$ in the basis of $5S_{1/2},F=1$. (b) The dependence of the
fidelity on $\Omega_{1}/\Omega_{2}$ and $\Delta/\delta$ on the basis of
$5S_{1/2},F=2$. The red area where the fidelity is close to one indicates the
optimized conditions to realize the $J_{x}$ Hamiltonian. The black dashed line
is a contour line where the fidelity is 0.99.

Since we want to obtain a pure $J_{x}$ Hamiltonian by Raman processes, we
simulate how $J_{z}$ evolves with the application of $H_{\rm eff}$ in Eq. (3).
$J_{z}$ is turned into $-J_{y}$ after a $90^{\circ}x$ pulse. To find the
optimal detunings and Rabi frequencies of the Raman fields for a $J_{x}$
Hamiltonian, we calculate the fidelity of $-J_{y}$ and the state after a
$90^{\circ}$ rotation operation by $H_{\rm eff}$. The fidelity formula is
$\mathcal{F}={\rm Tr}|\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}|^{2}$
jeong2007quantum . $\rho_{1}$ is $-J_{y}$ and $\rho_{2}$ is the state after a
$90^{\circ}$ rotation operation by $H_{\rm eff}$. If the fidelity is close to
one, we can say that $H_{\rm eff}$ is a perfect $J_{x}$ Hamiltonian.

Figure 3a shows the dependence of the fidelity on the dimensionless parameters
$\Omega_{1}/\Omega_{2}$ and $\Delta/\delta$ for $F=1$. When the Raman detuning
is large enough ($\Delta\gg\delta$), the fidelity is close to one. Then, we
can write the effective Hamiltonian as

$\displaystyle H_{\rm
eff}=\frac{\Omega_{1}\Omega_{2}}{24\Delta}\left(\begin{array}[]{ccc}0&1&0\\\
1&0&1\\\ 0&1&0\\\ \end{array}\right).$ (5)

The coefficient $\Omega_{1}\Omega_{2}/(24\Delta)$ is the effective Rabi
frequency.

Figure 4:  The density matrix evolution under the quadrupolar echo
$\tau-90^{\circ}y-\tau$ in the $J_{y}$ representation. (a-d) are real parts of
the density matrix. (e-h) are imaginary parts of the density matrix.
$\rho_{1}$ is the state $-J_{y}$. $\rho_{2}$ is the state after free evolution
with the NLZ dephasing effect. $\rho_{3}$ is the state after the $90^{\circ}y$
pulse. $\rho_{4}$ is the final state after DD for a whole Larmor precession
cycle.

Note that in Eq. (5), the coherence between $m_{F}=1$ and $m_{F}=-1$
disappears. This is caused by quantum interference between the two Raman
transition channels from $5S_{1/2},F=1$ to $5P_{1/2},F^{\prime}=1$ and
$F^{\prime}=2$. Under appropriate detunings and Rabi frequencies of the Raman
lights, the two transition amplitudes form destructive interference.

Following the same steps, we can construct an efficient $J_{x}$ Hamiltonian
for the energy level $5S_{1/2},F=2$. Similarly, we can plot the fidelity of
$-J_{y}$ and the state after a $90^{\circ}$ rotation operation by the $H_{\rm
eff}$ on the basis of $5S_{1/2},F=2$, as shown in Fig. 3b. Compared with the
$F=1$ case, the optimal condition for $F=2$ is more rigorous. We need a larger
Raman detuning to achieve the same fidelity (dashed line in Fig. 3).

In total, there are three necessary conditions for Raman lights to construct a
perfect $J_{x}$ rotation Hamiltonian: (i) The Rabi frequencies and phases of
two circularly polarized lights are equal. (ii) The $\pi$ polarized light has
a $\pi/2$ phase difference with the two circularly polarized lights. (iii) The
Raman detuning $\Delta$ is far off resonance.

For both the $F=1$ and $F=2$ systems, an effective $J_{y}$ Hamiltonian can be
realized in the same way by changing the propagating directions of the Raman
lasers from the $\hat{x}$ direction to the $\hat{y}$ direction.

##  IV Dynamic decoupling

In Section III, we give an effective rotation Hamiltonian by Raman process.
Next, we utilize the effective Hamiltonian to construct a DD sequence. Through
the dynamic evolution of the spin state, we give the optimal DD schemes for
different atomic systems.

Consider an atomic system with a total angular momentum $F$ (for the ground
state of ${}^{87}\mathop{\rm Rb}$ $5S_{1/2}$, $F=1$ or $2$) interacting with a
leading magnetic field along $\hat{z}$ in the Earth-field range, we assume
that the atomic spins are prepared in the $m_{F}=F$ state along the $\hat{z}$
direction ($J_{z}$) by an optical pumping (OP) field. The spin state is
rotated to $-J_{y}$ by a $90^{\circ}x$ pulse to generate the coherence between
Zeeman sublevels as quantum sensor or carrier of information. After that, the
spin state begins rotating freely in the $\hat{x}-\hat{y}$ plane because of
Larmor precession. In presence of the NLZ effect and magnetic field gradient,
the spin state begin diffusing, which broadens the linewidth of magnetic
resonance and decreases the coherence time.

###  IV.1 NLZ cancellation

To eliminate the NLZ effect by DD, we choose the quadrupolar echo
$\tau-90^{\circ}y-\tau$. With the $H_{\rm eff}$ obtained in Section III, the
spin rotation operations $90^{\circ}x$ and $90^{\circ}y$ can be realized by
controlling the Rabi frequencies, detunings, propagating directions and pulse
lengths of the Raman lasers.

For the $F=1$ system, we calculate the density matrix evolution under the
quadrupolar echo $\tau-90^{\circ}y-\tau$ in detail in the Appendix B. In terms
of the expression of $\rho(2\tau)$, the NLZ effect still exists. If we choose
an appropriate free evolution time $\tau$ such that $\Omega_{L}\tau=\pi$,
$\rho(2\tau)$ turns into

$\displaystyle\rho(2\tau)=i\frac{1}{\sqrt{2}}\left(\begin{array}[]{ccc}0&1&0\\\
-1&0&1\\\ 0&-1&0\\\ \end{array}\right).$ (6)

To realize $\Omega_{L}\tau=\pi$, atomic spin can be drived by turning on the
DD sequence with an optical rotation signal as the trigger (self-oscillating
mode) higbie2006robust . The final state is refocused to $-J_{y}$. An
alternative approach is to carry out quadrature detection, in which only the
components of $\rho(2\tau)$ irrelevant with respect to $\omega$ are detected
antonijevic2003refocussing ; levitt2013spin . Thus, the NLZ effect is removed
by $90^{\circ}y$ pulses with arbitrary free evolution time $\tau$ as long as
$\tau$ is much smaller than the NLZ dephasing time and much bigger than the
pulse length.

In the $F=2$ system, due to more NLZ splitting, the spin state dynamics under
the quadrupolar echo is more complicated. To see the DD performance clearly,
we present the density matrix evolution of the electron spin in the $J_{y}$
representation, as plotted in Fig 4. Figs. 4 (a-d) show the real parts of the
density matrix, and Figs. 4 (e-h) show the imaginary parts of the density
matrix. In the $J_{y}$ representation, $-J_{y}$ means that atoms are fully
populated in the $m_{F}=-F$ state (Fig. 4a, 4e). During the free precession of
state $-J_{y}$ in half of the Larmor precession period $\tau$ under the
influence of the Larmor precession Hamiltonian $\Omega_{L}J_{z}$, the spin
state rotates from $-J_{y}$ to $J_{y}$, as shown in Fig. 4b. At the same time,
the nondiagonal terms $\rho_{F,F-2}$ and $\rho_{F-2,F}$ occur because the NLZ
effect Hamiltonian $\omega J_{z}^{2}$ brings $\Delta m_{F}=2$ transitions, as
shown in Fig. 4f. The produced coherence terms lead to a reduction in the
atomic polarization. To eliminate the atomic decoherence effect, we apply a
$90^{\circ}y$ pulse, which adds a $\pi$ phase to the terms $\rho_{F,F-2}$ and
$\rho_{F-2,F}$ (Fig. 4g) and keeps the target term $\rho_{F,F}$ unchanged
(Fig. 4c) (see Appendix C). Finally, the dephasing effect in the second free
evolution $\tau$ counteracts the dephasing effect from the first free
evolution $\tau$, which refocuses the final spin state to $-J_{y}$ (Fig. 4d,
4h).

Figure 5: The state fidelity comparison between the cases free evolution and
with DD applied after five cycles. $\eta$ is the ratio of the revival
frequency $\omega$ and Larmor frequency $\Omega_{L}$. Insert is the zoom in
where real experimental parameters locate. Dashed line corresponds to the case
of ${}^{87}\mathop{\rm Rb}$ atoms in the Earth-field range.

Different from the $F=1$ system, the density matrix of the $F=2$ system
includes more nonzero terms induced by the NLZ effect. The $90^{\circ}y$ pulse
cannot flip them all. In the $J_{y}$ representation, the $90^{\circ}y$
operation turns density matrix $\rho_{mn}$ to $\rho^{\prime}_{mn}$ with
$\rho^{\prime}_{mn}=\exp[i(m+n)\pi/2]\rho_{mn}$, where m and n denote the
magnetic quantum number (see Appendix C). There are 25 elements in the density
matrix in total in the $F=2$ system, of which only $\rho_{2,0}$, $\rho_{0,2}$,
$\rho_{-2,0}$ and $\rho_{0,-2}$ can be inverted by the $90^{\circ}y$ pulse,
meaning that the quadrupolar echo cannot cancel out the NLZ effect completely.
The cancellation performance depends on the ratio $\eta$ of the revival
frequency $\omega$ and Larmor frequency $\Omega_{L}$. We define
$\eta=\omega/\Omega_{L}$. The fidelity is introduced to evaluate the DD
performance pokharel2018demonstration ; F2  and formulated as
$\mathcal{F}={\rm
Tr}|\sqrt{\sqrt{\rho_{1}^{\prime}}\rho_{2}^{\prime}\sqrt{\rho_{1}^{\prime}}}|^{2}$
jeong2007quantum . Here, $\rho_{1}^{\prime}$ is $-J_{y}$ and
$\rho_{2}^{\prime}$ is the state evolving for five DD cycles. Figure 5
compares the fidelity against $\eta$ in the cases with and without DD. From
Fig. 5, we find that the fidelity is maintained as high as $99.99\%$ in a wide
range $\eta<1\times 10^{-3}$, where most atomic systems are located
seltzer2008developments . For example, when ${}^{87}\mathop{\rm Rb}$ atoms
operate in the Earth-field range, the corresponding ratio $\eta$ is $5.1\times
10^{-5}$ (see Fig. 5 inset).

###  IV.2 Magnetic field gradient cancellation

A magnetic field gradient is inevitable in a natural environment. The phase
perturbation due to an inhomogeneous magnetic field accumulates over time and
causes a small spin orientation broadening along the target spin state,
formulated as $\phi=(\Omega_{L}+\Delta\Omega_{L})\tau$. First, this dephasing
effect reduces the coherence time of the atomic electron spin. Second, for a
long precession time and large phase perturbation, the condition
$\Omega_{L}\tau=\pi$ cannot be achieved for all atoms. Hence, the quadrupolar
echo cannot refocus the spin state perfectly, which compromises the
performance of NLZ effect cancellation.

Dephasing due to an inhomogeneous magnetic field can be largely alleviated by
using the Hahn echo rui2015operating ; hahn1950spin , where a $180^{\circ}y$
pulse is applied during the free evolution interval. In our scheme, the
$180^{\circ}y$ pulse is realized by Raman scattering.

Figure 6: DD control sequence by inserting Hahn echo $\tau-180^{\circ}y-\tau$
into quadrupolar echo $\tau-90^{\circ}y-\tau$ for (a) $F=1$ system. (b) $F=2$
system. The curved line represents projection of the evolved spin state to
$J_{y}$ in the Larmor precession process. The Raman pulses are all square
waves. Larmor period $T$ is 28.5 $\mu$s. The free evolution time $\tau_{1}$ is
adjustable from about 14.3 $\mu$s to about 5.59 ms and $\tau_{2}$ is 14.3
$\mu$s. $90^{\circ}y$ pulse length $t_{1}$is 2.85 $\mu$s. $180^{\circ}y$ pulse
length $t_{2}$ is 5.7 $\mu$s.

To remove the magnetic field gradient altogether, we apply a $180^{\circ}y$
pulse at the end of each Larmor cycle, as shown in Fig.6. Note that the
$180^{\circ}y$ pulses have no influence on the NLZ effect cancellation (see
Appendix C). Figure 6 gives the specific DD control schemes. The Larmor
frequency for ${}^{87}\mathop{\rm Rb}$ in the Earth’s magnetic field is
$\Omega_{L}=2\pi\times 350$ kHz. The corresponding Larmor period $T$ is 28.5
$\mu$s. For $F=1$ system, the free evolution time $\tau_{1}$ is arbitrary as
long as it is smaller than the NLZ dephasing time and much bigger than the
pulse length using quadrature detection. Based on the revival frequency
$\omega=2\pi\times 17.9$ kHz, $\tau_{1}$ is adjustable from about 14.3 $\mu$s
to about 5.59 ms. For $F=2$ system, the free evolution time $\tau_{2}$ is half
of the Larmor period 14.3 $\mu$s. To get a fast rotation operation, the length
of square Raman pulse can be set as $2.85\,\mu$s for $90^{\circ}$ rotation and
$5.7\,\mu$s for $180^{\circ}$ rotation. Thus, the effective Rabi frequency
$\Omega_{1}\Omega_{2}/(24\Delta)$ is $2\pi\times 0.875$ MHz. For a Raman
detuning 4 GHz, Rabi frequencies of the three polarized Raman lights can be
set as $2\pi\times$290 MHz equally.

##  V Conclusion and discussion

We demonstrate an all-optical DD scheme to suppress the NLZ effect and
magnetic field gradient at the same time. We design an all-optical control
sequence by inserting Hahn echo $\tau-180^{\circ}y-\tau$ into quadrupolar echo
$\tau-90^{\circ}y-\tau$ via Raman process. The NLZ effect can be canceled out
by the quadrupolar echo. There are two ways to set the free evolution time
$\tau$. One is to set $\tau$ to half of the Larmor precession period, which
can be realized by self-oscillating. The other is to use arbitrary $\tau$ with
quadrature detection as long as $\tau$ is much smaller than the NLZ dephasing
time and much bigger than the pulse length. For the $F=1$ system, the NLZ
effect can be completely removed by both methods. For the $F=2$ system, the
$90^{\circ}y$ pulses are added at every half time of the Larmor precession
period. Due to the greater amount of NLZ splitting, the NLZ effect cannot be
completely canceled out, and the decoupling results rely on the ratio $\eta$
of the revival frequency to the Larmor precession frequency. By analyzing the
fidelity of the atomic spin state with and without DD, we find that the
fidelity is maintained as high as $99.99\%$ in a wide range $\eta<1\times
10^{-3}$ after five DD cycles. The magnetic field gradient can be completely
eliminated by Hahn echo.

For ${}^{87}\mathop{\rm Rb}$ in the Earth-field range, the signal splitting
due to the NLZ effect is $2\pi\times 35.8$ Hz for $F=1$ system and $2\pi\times
71.6$ Hz for $F=2$ system. The typical relative inhomogeneity of magnetic
field is about $10^{-4}$ magnetic , leading to the signal broadening of
$2\pi\times 35$ Hz. As a consequence, the coherence time is 27.9 ms for $F=1$
system and 14 ms for $F=2$ system. After removing the NLZ effect and magnetic
field gradient by DD, the coherence time is only limited by the electron spin-
destruction collisions with the cell walls, residual relaxation due to
collisions with the reservoir, and alkali-metal–alkali-metal spin-exchange
collisions, which reaches 60 s 77s . The sensitivity of atomic magnetometer
can be improved by 46.4 times for $F=1$ system and 65.5 times for $F=2$ with
DD applied.

The Larmor frequency for ${}^{87}\mathop{\rm Rb}$ in the Earth’s magnetic
field is $\Omega_{L}=2\pi\times 350$ kHz. To get a fast rotation operation,
the length of square Raman pulse can be set as $2.85\,\mu$s for $90^{\circ}$
rotation and $5.7\,\mu$s for $180^{\circ}$ rotation. The effective Rabi
frequency is $2\pi\times 0.875$ MHz. For a Raman detuning 4 GHz, Rabi
frequencies of the three polarized Raman lights can be set as $2\pi\times$290
MHz equally. Our pulsed DD scheme suppresses the power broadening
significantly compared with the continuous spin locking approach. Furthermore,
this all-optical DD does not bring crosstalk between adjacent samples, which
is essential for multisensor detection. It has promising prospects for
practical application in quantum sensing and quantum information.

## Acknowledgments

The authors would like to thank Min Jiang from the University of Science and
Technology of China for useful discussions. The authors acknowledge financial
support from the National Key Research and Development Program of China
(2016YFA0302001); the National Natural Science Foundation of China (11874152,
11654005); the Fundamental Research Funds for the Central Universities; the
Shanghai Municipal Science and Technology Major Project (2019SHZDZX01); and
the Fellowship of China Postdoctoral Science Foundation (2020TQ0193); W. Z.
also acknowledges additional support from the Shanghai talent program.

##  Appendix A Adiabatic Approximation

When multilevel atoms interact with optical fields, the block matrix form of
the interaction Hamiltonian can be written as

$\displaystyle H=\left(\begin{array}[]{cc}H_{gg}&H_{ge}\\\ H_{eg}&H_{ee}\\\
\end{array}\right).$ (7)

The corresponding block density matrix is

$\displaystyle\rho=\left(\begin{array}[]{cc}\rho_{gg}&\rho_{ge}\\\
\rho_{eg}&\rho_{ee}\\\ \end{array}\right).$ (8)

Substituting these two expressions into the Liouville equation
$id\rho/dt=[H,\rho]$, we have

$\displaystyle i\frac{d\rho_{gg}}{dt}$
$\displaystyle=H_{gg}\rho_{gg}+H_{ge}\rho_{eg}-\rho_{gg}H_{gg}-\rho_{ge}H_{eg},$
(9a) $\displaystyle i\frac{d\rho_{ge}}{dt}$
$\displaystyle=H_{gg}\rho_{ge}+H_{ge}\rho_{ee}-\rho_{gg}H_{ge}-\rho_{ge}H_{ee},$
(9b) $\displaystyle i\frac{d\rho_{eg}}{dt}$
$\displaystyle=H_{eg}\rho_{gg}+H_{ee}\rho_{eg}-\rho_{eg}H_{gg}-\rho_{ee}H_{eg},$
(9c) $\displaystyle i\frac{d\rho_{ee}}{dt}$
$\displaystyle=H_{eg}\rho_{ge}+H_{ee}\rho_{ee}-\rho_{eg}H_{ge}-\rho_{ee}H_{ee}.$
(9d)

Usually, the ground-state Zeeman splitting is hundreds of kHz, while the Raman
detuning $\Delta$ is several GHz. Thus, compared with $H_{ee}$, $H_{gg}$ is
negligible. In addition, the Raman detuning $\Delta$ is much larger than the
natural linewidth of the excited energy level $5P_{1/2}$. Therefore, we can
carry out an adiabatic approximation for the excited state, that is,
$d\rho_{ge}/dt=0$, $d\rho_{eg}/dt=0$. Additionally, in the Raman scattering
process, the number of excited atoms is far less than the number of unexcited
atoms. Thus, we can simply set $\rho_{ee}=0$. With all the assumptions above,
Eq. (A3) is reduced to

$\displaystyle 0$ $\displaystyle=-\rho_{gg}H_{ge}-\rho_{ge}H_{ee},$ (10)
$\displaystyle 0$ $\displaystyle=H_{eg}\rho_{gg}+H_{ee}\rho_{eg}.$

Then, we have

$\displaystyle\rho_{ge}$ $\displaystyle=-\rho_{gg}H_{ge}H_{ee}^{-1},$ (11)
$\displaystyle\rho_{eg}$ $\displaystyle=-H_{ee}^{-1}H_{eg}\rho_{gg}.$

Substituting Eq. (A5) into Eq. (A3), we end up with

$\displaystyle i\frac{d\rho_{gg}}{dt}=$ $\displaystyle
H_{gg}\rho_{gg}-H_{ge}H_{ee}^{-1}H_{eg}\rho_{gg}-$ (12)
$\displaystyle\rho_{gg}H_{gg}+\rho_{gg}H_{ge}H_{ee}^{-1}H_{eg}$
$\displaystyle=$ $\displaystyle[H_{gg}+H_{eff},\rho_{gg}].$

Furthermore, we obtain

$H_{eff}=-H_{ge}H_{ee}^{-1}H_{eg}.$ (13)

With the help of the ADM package in Mathematica ADM , we can write out the
specific efficient Hamiltonian on the basis of the Zeeman sublevels.

##  Appendix B State evolution with the quadrupolar echo

For the $F=1$ system, we calculate the state evolution under the quadrupolar
echo $\tau-90^{\circ}y-\tau$. Initially, the atoms are populated in the
$J_{z}$ state. The density matrix after a $90^{\circ}x$ pulse evolves into

$\displaystyle J_{z}\stackrel{{\scriptstyle
90^{\circ}x}}{{\longrightarrow}}-J_{y}=\rho(0).$ (14)

Free precession of this state for a period $\tau$ under the influence of the
Larmor precession Hamiltonian $\Omega_{L}J_{z}$ and the NLZ effect Hamiltonian
$\omega J_{z}^{2}$ produces the following transformations:

$\displaystyle\stackrel{{\scriptstyle\Omega_{L}J_{z}\tau}}{{\longrightarrow}}-$
$\displaystyle[J_{y}\cos(\Omega_{L}\tau)-J_{x}\sin(\Omega_{L}\tau)],$ (15)
$\displaystyle\stackrel{{\scriptstyle\omega
J_{z}^{2}\tau}}{{\longrightarrow}}-$
$\displaystyle\left\\{\left[J_{y}\cos(\Omega_{L}\tau)-iC_{y}\sin(\omega\tau)\right]\cos(\Omega_{L}\tau)-\right.$
$\displaystyle\left.[J_{x}\cos(\Omega_{L}\tau)-iC_{x}\sin(\omega\tau)]\sin(\Omega_{L}\tau)\right\\}.$

Since $\Omega_{L}J_{z}$ commutates with $\omega J_{z}^{2}$, the order in which
we perform the transformations in Eq. (B2) is irrelevant.

Then, a $90^{\circ}y$ pulse is applied to the final state in Eq. (B2),
yielding

$\displaystyle\stackrel{{\scriptstyle 90^{\circ}y}}{{\longrightarrow}}-$
$\displaystyle\left\\{[J_{y}\cos(\omega\tau)+iC_{y}\sin(\omega\tau)]\cos(\Omega_{L}\tau)-\right.$
(16)
$\displaystyle\left.[-J_{z}\cos(\omega\tau)-iC_{q}\sin(\omega\tau)]\sin(\Omega_{L}\tau)\right\\}.$

After the $90^{\circ}y$ pulse, an echo is formed following a second period
$\tau$ of free precession:

$\displaystyle\stackrel{{\scriptstyle\Omega_{L}J_{z}\tau}}{{\longrightarrow}}-$
$\displaystyle\left\\{\left\\{[J_{y}\cos(\Omega_{L}\tau)-J_{x}\sin(\Omega_{L}\tau)]\cos(\omega\tau)\right.\right.$
(17)
$\displaystyle+\frac{1}{\sqrt{2}}[C_{p}\cos(\Omega_{L}\tau)-iC_{a}\sin(\Omega_{L}\tau)]$
$\displaystyle\left.\sin(\omega\tau)\right\\}\cos(\Omega_{L}\tau)\left\\{-J_{z}\cos(\omega\tau)-i[C_{q}\right.$
$\displaystyle\left.\cos(2\Omega_{L}\tau)-iC_{b}\sin(2\Omega_{L}\tau)]\sin(\omega\tau)\right\\}$
$\displaystyle\left.\sin(\Omega_{L}\tau)\right\\},$
$\displaystyle\stackrel{{\scriptstyle\omega
J_{z}^{2}\tau}}{{\longrightarrow}}-$
$\displaystyle\left\\{\left\\{[J_{y}\cos(\omega\tau)-iC_{y}\sin(\omega\tau)]\cos(\Omega_{L}\tau)\right.\right.$
(18)
$\displaystyle-[J_{x}\cos(\omega\tau)-iC_{x}\sin(\omega\tau)]\cos(\omega\tau)$
$\displaystyle+\frac{1}{\sqrt{2}}\left\\{[C_{p}\cos(\omega\tau)-iC_{i}\sin(\omega\tau)]\right.$
$\displaystyle\cos(\Omega_{L}\tau)-i[C_{a}\cos(\omega\tau)-i\sqrt{2}J_{x}$
$\displaystyle\left.\left.\sin(\omega\tau)]\sin(\Omega_{L}\tau)\right\\}\sin(\omega\tau)\right\\}\cos(\Omega_{L}\tau)-$
$\displaystyle\left\\{\left\\{-J_{z}\cos(\omega\tau)-i[C_{q}\cos(2\Omega_{L}\tau)-iC_{b}\right.\right.$
$\displaystyle\left.\left.\left.\sin(2\Omega_{L}\tau)]\sin(\omega\tau)\right\\}\right\\}\sin(\Omega_{L}\tau)\right\\}=\rho(2\tau),$

where

$C_{x}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{ccc}0&1&0\\\ -1&0&-1\\\
0&1&0\\\ \end{array}\right)$,
$C_{y}=i\frac{1}{\sqrt{2}}\left(\begin{array}[]{ccc}0&-1&0\\\ -1&0&1\\\
0&1&0\\\ \end{array}\right)$,

$C_{n}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{ccc}0&0&1\\\ 0&0&0\\\ -1&0&0\\\
\end{array}\right)$, $C_{i}=\left(\begin{array}[]{ccc}0&1&0\\\ -1&0&1\\\
0&-1&0\\\ \end{array}\right)$,

$C_{a}=\left(\begin{array}[]{ccc}0&1&0\\\ -1&0&-1\\\ 0&1&0\\\
\end{array}\right)$, $C_{p}=i\sqrt{2}C_{y},C_{q}=\sqrt{2}C_{n}$.

With these substitutions, the matrix form of $\rho(2\tau)$ at the end of the
second free precession period is

$\begin{split}\rho(2\tau)=\left(\begin{array}[]{ccc}-\cos(w\tau)\sin(\Omega_{L}\tau)&\frac{\cos(\Omega_{L}\tau)(i\cos(\Omega_{L}\tau)+\sin(\Omega\tau))}{\sqrt{2}}&-ie^{-2i\Omega_{L}\tau}\sin(w\tau)\sin(\Omega_{L}\tau)\\\
\frac{-i(1+e^{-2i\Omega_{L}\tau})}{2\sqrt{2}}&0&\frac{\cos(\Omega_{L}\tau)(i\cos(\Omega_{L}\tau)+\sin(\Omega_{L}\tau))}{\sqrt{2}}\\\
ie^{-2i\Omega\tau}\sin(w\tau)\sin(\Omega_{L}\tau)&\frac{i(1+e^{-2i\Omega_{L}\tau})}{2\sqrt{2}}&\cos(w\tau)\sin(\Omega_{L}\tau)\\\
\end{array}\right).\end{split}$ (19)

##  Appendix C Spin rotation in the $J_{y}$ representation

In the $J_{y}$ representation, the matrix form of the $90^{\circ}y$ operation
on the basis of the $F=2$ Zeeman sublevels is

$\displaystyle{\rm
rot}_{90^{\circ}y}=\left(\begin{array}[]{ccccc}-1&0&0&0&0\\\ 0&i&0&0&0\\\
0&0&1&0&0\\\ 0&0&0&i&0\\\ 0&0&0&0&-1\end{array}\right).$ (20)

With the $90^{\circ}y$ pulse applied, an arbitrary density matrix $\rho$

$\displaystyle\rho=\left(\begin{array}[]{ccccc}\rho_{2,2}&\rho_{2,1}&\rho_{2,0}&\rho_{2,-1}&\rho_{2,-2}\\\
\rho_{1,2}&\rho_{1,1}&\rho_{1,0}&\rho_{1,-1}&\rho_{1,-2}\\\
\rho_{0,2}&\rho_{0,1}&\rho_{0,0}&\rho_{0,-1}&\rho_{0,-2}\\\
\rho_{-1,2}&\rho_{-1,1}&\rho_{-1,0}&\rho_{-1,-1}&\rho_{-1,-2}\\\
\rho_{-2,2}&\rho_{-2,1}&\rho_{-2,0}&\rho_{-2,-1}&\rho_{-2,-2}\\\
\end{array}\right)$ (21)

is transformed into

$\displaystyle\rho^{\prime}=\left(\begin{array}[]{ccccc}\rho_{2,2}&i\rho_{2,1}&-\rho_{2,0}&-i\rho_{2,-1}&\rho_{2,-2}\\\
-i\rho_{1,2}&\rho_{1,1}&i\rho_{1,0}&-\rho_{1,-1}&-i\rho_{1,-2}\\\
-\rho_{0,2}&-i\rho_{0,1}&\rho_{0,0}&i\rho_{0,-1}&-\rho_{0,-2}\\\
i\rho_{-1,2}&-\rho_{-1,1}&-i\rho_{-1,0}&\rho_{-1,-1}&i\rho_{-1,-2}\\\
\rho_{-2,2}&i\rho_{-2,1}&-\rho_{-2,0}&-i\rho_{-2,-1}&\rho_{-2,-2}\\\
\end{array}\right),$ (22)

with the unitary transform $\rho^{\prime}={\rm rot}_{90^{\circ}y}\rho{\rm
rot}_{90^{\circ}y}^{\dagger}$.

A uniform expression for the density matrix elements is

$\displaystyle\rho^{\prime}_{mn}=\exp(i\frac{m+n}{2}\pi)\rho_{mn},$ (23)

where $m$ and $n$ are magnetic quantum numbers.

From Eq. (C4), we can see that the density matrix obtains a $\pi$ phase for
the elements $\rho_{2,0}$, $\rho_{0,2}$, $\rho_{-2,0}$ and $\rho_{0,-2}$ while
keeping the diagonal elements unchanged.

The $180^{\circ}y$ operation is equivalent to two $90^{\circ}y$ operations. A
uniform expression for density matrix elements under the $180^{\circ}y$
operation is

$\displaystyle\rho^{\prime}_{mn}=\exp[i(m+n)\pi]\rho_{mn}.$ (24)

## References

  * (1) I. Kominis, T. Kornack, J. Allred, and M. V. Romalis, A subfemtotesla multichannel atomic magnetometer, Nature 422, 596 (2003).
  * (2) J. B. Fixler, G. Foster, J. McGuirk, and M. Kasevich, Atom interferometer measurement of the newtonian constant of gravity, Science 315, 74 (2007).
  * (3) A. Dunning, R. Gregory, J. Bateman, N. Cooper, M. Himsworth, J. A. Jones, and T. Freegarde, Composite pulses for interferometry in a thermal cold atom cloud, Phys. Rev. A 90, 033608 (2014).
  * (4) B. Chen, C. Qiu, S. Chen, J. Guo, L. Chen, Z. Ou, and W. Zhang, Atom-light hybrid interferometer, Phys. Rev. Lett. 115, 043602 (2015).
  * (5) G. Bao, A. Wickenbrock, S. Rochester, W. Zhang, and D. Budker, Suppression of the nonlinear zeeman effect and heading error in earthfield-range alkali-vapor magnetometers, Phys. Rev. Lett. 120, 033202 (2018).
  * (6) G. Bao, D. Kanta, D. Antypas, S. Rochester, K. Jensen, W. Zhang, A. Wickenbrock, and D. Budker, All-optical spin locking in alkali-vapor magnetometers, arXiv preprint arXiv:1809.10906 (2018).
  * (7) K.-M. C. Fu, G. Z. Iwata, A. Wickenbrock, and D. Budker, Sensitive magnetometry in challenging environments, arXiv preprint arXiv:2008.00082 (2020).
  * (8) J. Saywell, M. Carey, M. Belal, I. Kuprov, and T. Freegarde, Optimal control of raman pulse sequences for atom interferometry, J. Phys. B: At. Mol. Opt. Phys. 53, 085006 (2020).
  * (9) M. Jiang, W. Xu, Q. Li, Z. Wu, D. Suter, and X. Peng, Interference in atomic magnetometry, Adv. Quantum Technol. p. 2000078 (2020).
  * (10) Y. Dudin, L. Li, and A. Kuzmich, Light storage on the time scale of a minute, Phys. Rev. A 87, 031801 (2013).
  * (11) L. Chen, G.-W. Zhang, C.-L. Bian, C.-H. Yuan, Z. Ou, and W. Zhang, Observation of the rabi oscillation of light driven by an atomic spin wave, Phys. Rev. Lett. 105, 133603 (2010).
  * (12) J. Rui, Y. Jiang, S.-J. Yang, B. Zhao, X.-H. Bao, and J.-W. Pan, Operating spin echo in the quantum regime for an atomic-ensemble quantum memory, Phys. Rev. Lett. 115, 133002 (2015).
  * (13) J. Guo, X. Feng, P. Yang, Z. Yu, L. Chen, C.-H. Yuan, and W. Zhang, High-performance raman quantum memory with optimal control in room temperature atoms, Nat. Commun. 10, 1 (2019).
  * (14) G. K. Brennen, D. P. OLeary, and S. S. Bullock, Criteria for exact qudit universality, Phys. Rev. A 71, 052318 (2005).
  * (15) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Realization of two-dimensional spin-orbit coupling for bose-einstein condensates, Science 354, 83 (2016).
  * (16) V. Acosta, M. Ledbetter, S. Rochester, D. Budker, D. J. Kimball, D. Hovde, W. Gawlik, S. Pustelny, J. Zachorowski, and V. Yashchuk, Nonlinear magneto-optical rotation with frequency-modulated light in the geophysical field range, Phys. Rev. A 73, 053404 (2006).
  * (17) S. Seltzer, P. Meares, and M. Romalis, Synchronous optical pumping of quantum revival beats for atomic magnetometry, Phys. Rev. A 75, 051407 (2007).
  * (18) K. Jensen, V. Acosta, J. Higbie, M. Ledbetter, S. Rochester, and D. Budker, Cancellation of nonlinear zeeman shifts with light shifts, Phys. Rev. A 79, 023406 (2009).
  * (19) E. Alexandrov, Recent progress in optically pumped magnetometers, Phys. Scr. 2003, 27 (2003).
  * (20) S. Antonijevic and S. Wimperis, Refocussing of chemical and paramagnetic shift anisotropies in 2h nmr using the quadrupolar-echo experiment, J. Magn. Reson. 164, 343 (2003).
  * (21) G. De Lange, Z. Wang, D. Riste, V. Dobrovitski, and R. Hanson, Universal dynamical decoupling of a single solid-state spin from a spin bath, Science 330, 60 (2010).
  * (22) B. Naydenov, F. Dolde, L. T. Hall, C. Shin, H. Fedder, L. C. Hollenberg, F. Jelezko, and J. Wrachtrup, Dynamical decoupling of a single electron spin at room temperature, Phys. Rev. B 83, 081201 (2011).
  * (23) R. Shaniv, N. Akerman, T. Manovitz, Y. Shapira, and R. Ozeri, Quadrupole shift cancellation using dynamic decoupling, Phys. Rev. Lett. 122, 223204 (2019).
  * (24) J. S. Pedernales, G. W. Morley, and M. B. Plenio, Motional dynamical decoupling for interferometry with macroscopic particles, Phys. Rev. Lett. 125, 023602 (2020).
  * (25) K. C. Miao, J. P. Blanton, C. P. Anderson, A. Bourassa, A. L. Crook, G. Wolfowicz, H. Abe, T. Ohshima, and D. D. Awschalom, Universal coherence protection in a solid-state spin qubit, Science 10, 1126 (2020).
  * (26) E. L. Hahn, Spin echoes, Phys. Rev. 80, 580 (1950).
  * (27) C. Barthel, J. Medford, C. Marcus, M. Hanson, and A. Gossard, Interlaced dynamical decoupling and coherent operation of a singlet-triplet qubit, Phys. Rev. Lett. 105, 266808 (2010).
  * (28) M. A. A. Ahmed, G. A. Álvarez, and D. Suter, Robustness of dynamical decoupling sequences, Phys. Rev. A 87, 042309 (2013).
  * (29) B. Pokharel, N. Anand, B. Fortman, and D. A. Lidar, Demonstration of fidelity improvement using dynamical decoupling with superconducting qubits, Phys. Rev. Lett. 121, 220502 (2018).
  * (30) S. Upadhyay, U. Dargyte, D. Patterson, and J. D. Weinstein, Ultralong spin-coherence times for rubidium atoms in solid parahydrogen via dynamical decoupling, Phys. Rev. Lett. 125, 043601 (2020).
  * (31) G. Bison, N. Castagna, A. Hofer, P. Knowles, J.-L. Schenker, M. Kasprzak, H. Saudan, and A. Weis, A room temperature 19- channel magnetic field mapping device for cardiac signals, Appl. Phys. Lett. 95, 173701 (2009).
  * (32) G. Lembke, S. Erné, H. Nowak, B. Menhorn, A. Pasquarelli, and G. Bison, Optical multichannel room temperature magnetic field imaging system for clinical application, Biomed. Opt. Express 5, 876 (2014).
  * (33) A. Borna, T. R. Carter, J. D. Goldberg, A. P. Colombo, Y.-Y. Jau, C. Berry, J. McKay, J. Stephen, M. Weisend, and P. D. Schwindt, A 20- channel magnetoencephalography system based on optically pumped magnetometers, Phys. Medicine & Biol. 62, 8909 (2017).
  * (34) F. P. Bustos, D. B. Calia, D. Budker, M. Centrone, J. Hellemeier, P. Hickson, R. Holzlöhner, and S. Rochester, Remote sensing of geomagnetic fields and atomic collisions in the mesosphere, Nat. Commun. 9, 1 (2018).
  * (35) B. Patton, O. Versolato, D. C. Hovde, E. Corsini, J. M. Higbie, and D. Budker, A remotely interrogated all-optical 87 Rb magnetometer, Appl. Phys. Lett. 101, 083502 (2012).
  * (36) G. Breit and I. Rabi, Measurement of nuclear spin, Phys. Rev. 38, 2082 (1931).
  * (37) D. Budker and D. F. J. Kimball, _Optically polarized atoms: understanding light-atom interactions_ (Oxford University Press, 2010).
  * (38) H. Jeong, T. C. Ralph, and W. P. Bowen, Quantum and classical fidelities for gaussian states, Josa B 24, 355 (2007).
  * (39) J. Higbie, E. Corsini, and D. Budker, Robust, high-speed, all-optical atomic magnetometer, Rev. Sci. Instrum. 77, 113106 (2006).
  * (40) M. H. Levitt, _Spin dynamics: basics of nuclear magnetic resonance_ (John Wiley & Sons, 2013).
  * (41) X. Peng, D. Suter and D. A. Lidar, High fidelity quantum memory via dynamical decoupling: theory and experiment, J. Phys. B: At. Mol. Opt. Phys. 44, 154003 (2011).
  * (42) S. J. Seltzer, Developments in alkali-metal atomic magnetometry (Princeton University, 2008).
  * (43) Y. Zhang , Y. Li, Q. Jiang, Z. Wang, T. Xia, and H. Luo, Analytical design of axial magnetic coils with systematically improved uniformity for miniature quantum devices, Rev. Sci. Instrum. 90, 114706 (2019).
  * (44) M. V. Balabas, T. Karaulanov, M. P. Ledbetter, and D. Budker, Polarized Alkali-metal vapor with minute-long transverse spin-relaxation time, Phys. Rev. Lett. 105, 070801 (2010).
  * (45) http://rochesterscientific.com/ADM/.


Improvised Aerial Object Detection approach for YOLOv3 Using Weighted
Luminance

Sai Ganesh CSa (saiganeshcs@ieee.org), Aouthithiye Barathwaj SR Ya
(aouthithiyebarathwaj@ieee.org),R. Swethaa Sb (swethaa@ieee.org), R.
Azhagumurugana<EMAIL_ADDRESS>
a Department of Electrical and Electronics Engineering, Sri Sai Ram
Engineering College, Sai Leo Nagar, West Tambaram, Chennai, India
b Department of Electronics and Communication Engineering, Sri Sai Ram
Engineering College, Sai Leo Nagar, West Tambaram, Chennai, India

Corresponding Author:
Sai Ganesh CS
Department of Electrical and Electronics Engineering, Sri Sai Ram Engineering
College, Sai Leo Nagar, West Tambaram, Chennai, India
Tel: +91 8754444820
Email<EMAIL_ADDRESS>
###### Abstract

Remote sensing is the image acquisition of a target without having physical
contact with it. Now-a-days remote sensing data is widely preferred due to its
reduced image acquisition period. The remote sensing of ground targets is more
challenging because of the various factors that affect the propagation of
light through different mediums from a satellite acquisition. Several
Convolutional Neural Network-based algorithms are being implemented in the
field of remote sensing. Supervised learning is a machine learning technique
where the data is labeled according to their classes prior to the training. In
order to detect and classify the targets more accurately, YOLOv3, an algorithm
based on bounding and anchor boxes is adopted. In order to handle the various
effects of light traveling through the atmosphere, Grayscale based YOLOv3
configuration is introduced. For better prediction and for solving the
Rayleigh scattering effect, RGB based grayscale algorithms are proposed. The
acquired images are analysed and trained with the grayscale based YOLO3
algorithm for target detection. The results show that the grayscale-based
method can sense the target more accurately and effectively than the
traditional YOLOv3 approach.

###### keywords:

Remote Sensing , Computer Vision , Convolutional Neural Networks , Rayleigh
Scattering , Grayscale , YOLOv3

††journal: Expert Systems with Applications

##  1 Introduction

Remote sensing is the process of acquiring data from the earth’s surface
without being in actual physical contact with the target. It employs
electromagnetic radiation in the form of light to obtain the information about
the target on the earth’s surface. It is being spoken of with high values in
terms of data acquisition especially where human influence is very low. Remote
sensing in the field of defence has been in practice over a long period of
time due to its explicitly easier way of communication and data collection
over a larger area and we can easily and rapidly process and transmit the
data. Remote sensing is considered as the next technological extension to
photogrammetry. This enables humans to fetch data beyond normal human eye i.e.
infrared, thermal etc. This manages to not only observe objects physically but
it also enables them to detect, observe and also record the energy of the
target observed without being in a physical contact. Computer vision in the
field of remote sensing is gaining a huge popularity as it is used in various
applications such as geo-spatial data analysis and it involves using an
exponentially large dataset, which is time consuming and involves complex
algorithms. The panchromatic camera is the most commonly used device for the
purpose of surveying, obtaining images, etc. It is a single lens camera and
both CCD (Charge Coupled Devices) arrays and photographic films can be used in
the panchromatic camera. It has the ability to sense beyond the visible
wavelength and so is named accordingly. When the light is focused on the CCD,
it turns the incident photons into voltage and which is then quantized to get
the picture in the digital form. The major hurdles to the remote sensing
process are the scattering and absorption effects of light which affects the
quality of the image during the image acquisition.

YOLOv3 (You Only Look Once v3) algorithm proposed by Joseph Redmon Redmon &
Farhadi [2018], Redmon et al. [2016], Redmon & Farhadi [2017], a popular
computer vision algorithm is adopted in this article. This algorithm is based
on the darknet53 deep learning framework. The darknet53 network comprises of
53 convolutional layers along with batch normalizing layers and Leaky-ReLU
activation layers. This network acts as a backbone to the fully-connected
YOLOv3 network. Prediction of the coordinates of the bounding boxes is
achieved through the method of regression. Dimensional clusters are applied as
anchor boxes in order to obtain the required bounding boxes. A special feature
of this algorithm is that a single neural network is applied to the full
image. Images are divided into multiple regions and coordinates of bounding
boxes are obtained. The probabilities for each regions are procured and
weighed and the high scoring regions are further processed. In this article, a
single channelled YOLOv3, which supports Grayscale images is proposed. This
also eliminates several man-made optical illusions, colour illusions and other
misinterpretationsBui et al. [2016]. The bottom surface of the Earth’s
atmosphere acts as a reflecting base for the surface and radiations inbound
Ma, and Liu, and Ren, and Luo, Yu-Feng [2019], Fasola & Veloso [2006],
Hadjimitsis et al. [2010a], Bui et al. [2017], Boyat & Joshi [2015], Abd Al-
salam Selami & Fadhil [2016], Mikołajczyk & Grochowski [2018], Hadjimitsis et
al. [2010b] This affects the quality of the image that is being collected.
Moreover, Earth’s atmosphere consists of various particles such as hydrogen,
oxygen, water vapour and many such particles which attenuates light. In remote
sensing, image data acquisition quality is affected by absorption of light
from particles such as ozone, water vapour and other atmospheric contents. In
addition to the atmospheric effects, Rayleigh and Mie scattering effects also
play a negative role in the field of remote sensing. Shorter wavelength gets
scattered more compared to that of larger wavelength. This causes a haze due
to the excessively scattered blue light in the captured aerial images hailing
to its shorter wavelength. In aerial photography various filters are prevalent
to remove these effects but they do not prove to be of much efficiency and are
of higher cost. Also, non-selective scattering of light happens due to clouds
and fog, thus making it difficult to interpret the image obtained. These
extensive parameters possess an immense challenge to the field of remote
sensing especially in the field of image acquisition..

##  2 Research Methodology

The sky would be black without the presence of atmospheric particles such as
NO2, O2 etc. and more importantly the scattering phenomena. During sunrise and
sunset, sunlight travels a comparatively longer distance in order to reach the
Earth’s atmosphere which tends to give a reddish hue to the color of the sky
since shorter wavelengths gets scattered easily Kusetogullari & Yavariabdi
[2018]. The images acquired at a higher altitude experience Rayleigh’s
scattering effect i.e. the shorter wavelengths appear pre-dominant over the
others. Thus, the targeted objects in the image dataset of RGB tend to get
affected by these factors. Also, some of the images tend to get camouflaged
with the surroundings making it completely diminished in size, which pose a
serious threat to the areas which border countries and other similar hotspots.
Fig. 1 represents the pictorial view of the Rayleigh scattering effect.

Fig. 1 (a) This represent the excessive presence of blue wavelength during
daytime and Fig. 1(b) represents the sky at sunset and sunrise where red and
yellow wavelengths presence largely. Image pre-processing in this article
varies with image acquisition time. This means that Black and white with red
filters is used to compensate for the pre-dominant blue color during the day
time. Black and white with blue filters are used during the time of sunrise
and sunset. Default grayscale images are used during the night time. Thus, all
the acquired images are pre-processed prior to the training using the above-
mentioned filters. Since the pre-processed images are produced in grayscale,
the configuration of the YOLOv3 network is also modified for the suitability
of the training. Hence the number of channels in the configuration of YOLOv3
is set to one. The default computing channel in the python interpreter is the
BGR (Blue Green Red) channel. So, in this case, the three-channeled BGR images
are converted to single-channeled grayscale images which are pre-programmed in
the YOLOv3 supporting files Kanan & Cottrell [2012], Prasad et al. [2015],
Welsh et al. [2002], Ghani et al. [2017], Wallace [1983], Xia et al. [2018],
Lu & Plataniotis [2009], Chandran [2010].

Figure 1: Pictorial representation of Rayleigh scattering effect, 1 (a)
describes the scattering of blue light in the Earth’s atmosphere, which makes
sky appear blue in color during the day time. 1 (b) describes the sunlight
traveling a comparatively longer distance to reach the Earth’s surface, making
the longer wavelength lights to dominate.

$Grayscale(x,y)=0.3\times R(x,y)+0.1\times G(x,y)+0.5\times B(x,y)$ (1)

The traditional conversion of RGB (Red, Green, Blue) channel to grayscale is
defined in the Eq (1).

During the day-time, blue light is pre-dominant making the longer wavelengths
of light less visible. In order to solve this, grayscale based red filter is
adopted to compensate the attenuation of longer wavelengths Levin et al.
[2004], Karthikeyani et al. [2007], Bala & Braun [2004], Kumar & Verma [2010].
Thus, increase in the intensity of red light is essential. Since green light
is not largely affected by the Rayleigh effect, intensity of green is reduced
to compensate the attenuation of red light. The conversion of RGB channeled
image to grayscale based red filter is defined in Eq (2).

$Grayscale_{weighted}(x,y)=W_{r}\times R(x,y)+W_{g}\times G(x,y)+W_{b}\times
B(x,y)$ (2)

The weighing factors for each channel is determined individually as given in
the equation (2). Where Wr, Wg, Wb are the weighing factors of red, green,
blue channels respectively. The summation of the above three weighing factors
is unity. The RGB histogram of the image is hence obtained from which the mean
and standard deviation is determined. The total number of image pixel as
percentage of it’s frequency is defined as the Perc Tempfli et al. [2009]. The
Perc is calculated with respect to the Table 1.

$\begin{split}w_{b}=Perc\times mean\\\ w_{g}=(1-W_{b})-(mean+StdDev)\\\
w_{r}=1-(w_{b}+w_{g})\end{split}$ (3)

During the sunset or sunrise, the attenuation of shorter wavelengths occurs,
resulting in a reddish tint in the acquired images. Hence, the grayscale based
blue filter is adopted in this scenario. Here the intensity of blue is
increased and the intensity of red is decreased. The conversion of RGB
channeled image into grayscale based blue filter is defined in equation (3).
The histogram of sunset/sunrise acquisition image is given in Figure 2 (d) and
(g).

$\begin{split}w_{b}=Perc\times mean\times(2\times stdDev)\\\
w_{g}=(1-W_{b})-(avg)\\\ w_{r}=1-(w_{b}+w_{g})\end{split}$ (4)

The conversion of RGB channeled image into grayscale based red filter is
defined in equation (4). The R (x, y) is the pixel value corresponding to Red,
G (x, y) is the pixel value corresponding to Green and B (x, y) is the pixel
value pertaining to Blue, where x and y are pixel locations. The histograms of
the original image and the proposed filters are given in Fig. 2.

DN | Npix | Perc | CumNpix | CumPerc
---|---|---|---|---
0 | 0 | 0.00 | 0 | 0.00
13 | 0 | 0.00 | 0 | 0.00
14 | 1 | 0.00 | 1 | 0.00
15 | 3 | 0.00 | 4 | 0.01
16 | 2 | 0.00 | 6 | 0.01
51 | 55 | 0.08 | 627 | 0.86
52 | 59 | 0.08 | 686 | 0.94
53 | 94 | 0.13 | 780 | 1.07
54 | 138 | 0.19 | 918 | 1.26
102 | 1392 | 1.90 | 25118 | 34.36
103 | 1719 | 2.35 | 26837 | 36.71
104 | 1162 | 1.59 | 27999 | 38.30
105 | 1332 | 1.82 | 29331 | 40.12
106 | 1491 | 2.04 | 30822 | 42.16
107 | 1685 | 2.31 | 32507 | 44.47
108 | 1399 | 1.91 | 33906 | 46.38
109 | 1199 | 1.64 | 35105 | 48.02
110 | 1488 | 2.04 | 36593 | 50.06
111 | 1460 | 2.00 | 38053 | 52.06
163 | 720 | 0.98 | 71461 | 97.76
164 | 597 | 0.82 | 72058 | 98.57
165 | 416 | 0.57 | 72474 | 99.14
166 | 274 | 0.37 | 72748 | 99.52
173 | 3 | 0.00 | 73100 | 100.00
174 | 0 | 0.00 | 73100 | 100.00
255 | 0 | 0.00 | 73100 | 100.00
Table 1: Histogram Table

The model starts training after the preprocessing of the acquired images. Fig.
2. represents the block diagram of the proposed methodology. The image dataset
is obtained from LANDSAT 8 via google earth engine which consists of numerous
images pertaining to 5 different classes namely, armed vehicles, aircraft,
helicopters, tents, and ships/boats. The LANDSAT 8 orbits around the earth at
an altitude of 705km. It is inclined at an angle of 98.2 degrees and it takes
around 99 minutes to complete one revolution. The images are captured at an
altitude of about 20-30 meters above ground level. With the help of LABELIMG
the images are annotated in the YOLO format. Thus, supervised learning method
is employed. The bounding boxes are drawn along with their respective classes.
The training and testing dataset are split as 80% and 20% respectively.

Dataset | Aircraft | Helicopters | Trucks | Ships | Tents
---|---|---|---|---|---
Training set | 220 | 220 | 220 | 220 | 220
Test set | 25 | 25 | 25 | 25 | 25
Table 2:  Dataset size and classification Image | Original Image | Processed Image
---|---|---
Mean | Std. Deviation |  Mean | Std. Deviation
$a$ | 128.73 | 50.93 | 98.44 | 51.54
$b$ | 116.92 | 57.98 | 87.60 | 55.57
$c$ | 102.46 | 56.06 | 75.04 | 52.28
$d$ | 87.76 | 31.92 | 92.38 | 32.56
$e$ | 126.58 | 28.91 | 139.36 | 26.87
$f$ | 94.01 | 35.61 | 97.86 | 36.42
Table 3:  Mean and standard deviation of the test case images for both
original and processed images

A test case of six images is taken and is processed with the proposed
algorithms. The images are obtained from the LandViewer database and Google
Earth Engine. The images from NIOP satellite are acquired during sunrise /
sunset as the angle of sun elevation is less than 10 degrees. The images
obtained from LANDSAT 8 are acquired with the angle of sun elevation greater
than 30 degrees. Hence grayscale blue filter is implemented in the NIOP images
and grayscale red filter is implemented in LANDSAT 8 images. The images are
given in the Figure 2 and the distribution of mean and standard deviation of
the histogram of RGB channel is given in Table 3.

Figure 2: Images 2(a)(i), 2(a)(ii) and 2(a)(iii) are the orginal images
obtained from NIOP satellite image database. Images 2(a)(iv), 2(a)(v) and
2(a)(vi) are the orginal images obtained from LANDSAT 8 satellite image
database 2(b)(i-vi) images are the processed images with the proposed
algorithm.

##  3 Experiment and Analysis

###  3.1 Training

The model is trained in the Google Colaboratory environment using Nvidia P4
16GB GPU, with Core i5-8259U 2.30GHz Quad core CPU with a memory of 8 GB. The
learning rate is set to 0.001.The momentum is set to 0.9, with a maximum batch
size of 2500. For every single training step, the batch size is set to 64 and
16 subdivisions i.e. 64 images for each step. The training steps were of batch
size 2000, 2250 is set at 80% and 90%.

###  3.2 Results

The Mean average precision (mAP) and F1-Score are calculated after every 100
iterations. The average loss is calculated after each iteration. The
experiment results show that the proposed method is more robust than the
YOLOv3 particularly in the field of remote sensing. The mAP results of YOLOv3
started at 67% and gradually increased to 75.5% (77% peak value) over the
iterations. In the proposed method, the mAP results started at 30%, and
gradually increased to 84.5%. The proposed model yields a maximum of 88% of
mAP result. Fig. 3 provides a graph to prove that the proposed method yields
better loss and mAP than the traditional YOLOv3.

Figure 3: . Graphical representation of the mAP and Loss calculations. X-axis represents the number of iterations and Y-axis represents the values for loss and mAP curves. The blue line represents the loss curve and the red line represents mAP. Fig 3 (a) represents the graph for training executed using YOLOv3 and Fig 3 (b) represents the graph of the training executed using the proposed method Method | Minimum Loss | mAP | Peak mAP
---|---|---|---
YOLOv3 | 0.153 | 75.5% | 77%
Proposed Method | 0.314 | 84.5% | 88%
Table 4: Experiment results of the implementation of YOLOv3 and the proposed
methodologies computing the minimum training loss and Mean Average Precision

###  3.3 Deployment

The deployment of the model is executed in PyCharm IDE with the trained
configuration and weights. With the help of time module library from Python,
we apply the respective filters corresponding to its time of image
acquisition. The confidence score for the object detection is set to 0.30.
Fig. 4 compares the images predicted by the YOLOv3 based approach and the
proposed approach.

Figure 4: List of Deployment images with different filters. 4 (a), (d), (g),
(j) are original acquired images. 4 (b), (e), (h), (k) are images with
grayscale based red filter. 4 (c), (f), (i), (l) are images with grayscale
based blue filter.

##  4 Conclusion

In this article, grayscale based algorithms are specifically implemented with
respect to the sun’s angle of elevation. The images are obtained from NIOP and
LANDSAT 8 satellite imagery using LandViewer database. All the acquired images
are analyzed and implemented with respect to its adaptable algorithm. After
the preprocessing of images, the images are trained with deep learning based
YOLOv3 algorithm for target detection. The training results are discussed and
the model is successfully deployed using PyCharm IDE. the future scope of this
project aims to program a machine that alarms any trespasses detection and
sends a warning through email and mobile applications.

## Acknowledgements

The authors would like to thank the management of Sri Sai Ram Engineering
College who gave technical support to this research work. The authors would
like to show their gratitude towards Dr. T. Porselvi and Dr. B. Meenakshi for
their continuous support and guidance.

## References

  * Bala & Braun [2004] Bala, R., & Braun, K. (2004).  Color-to-grayscale conversion to maintain discriminability.  Proc. SPIE, 5293. doi:10.1117/12.532192.
  * Boyat & Joshi [2015] Boyat, A., & Joshi, B. (2015).  A review paper: Noise models in digital image processing.  Signal & Image Processing : An International Journal, 6. doi:10.5121/sipij.2015.6206.
  * Bui et al. [2016] Bui, H., Lech, M., Cheng, E., Neville, K., & Burnett, I. (2016).  Using grayscale images for object recognition with convolutional-recursive neural network, .  (pp. 321–325). doi:10.1109/CCE.2016.7562656.
  * Bui et al. [2017] Bui, H., Lech, M., Cheng, E., Neville, K., & Burnett, I. (2017).  Object recognition using deep convolutional features transformed by a recursive network structure.  IEEE Access, PP, 1–1. doi:10.1109/ACCESS.2016.2639543.
  * Chandran [2010] Chandran, S. (2010).  Color image to grayscale image conversion, .  (pp. 196 – 199). doi:10.1109/ICCEA.2010.192.
  * Fasola & Veloso [2006] Fasola, J., & Veloso, M. (2006).  Real-time object detection using segmented and grayscale images.  Proceedings - IEEE International Conference on Robotics and Automation, (pp. 4088–4093). doi:10.1109/ROBOT.2006.1642330.
  * Ghani et al. [2017] Ghani, U., Gotoh, Y., & Nida, N. (2017).  Medical image colorization for better visualization and segmentation, .
  * Hadjimitsis et al. [2010a] Hadjimitsis, D., Papadavid, G., Agapiou, A., Themistocleous, K., Hadjimitsis, M., Retalis, A., Michaelides, S., Chrysoulakis, N., Toulios, L., & Clayton, C. (2010a).  Atmospheric correction for satellite remotely sensed data intended for agricultural applications: Impact on vegetation indices.  Natural Hazards and Earth System Sciences, 10. doi:10.5194/nhess-10-89-2010.
  * Hadjimitsis et al. [2010b] Hadjimitsis, D., Papadavid, G., Agapiou, A., Themistocleous, K., Hadjimitsis, M., Retalis, A., Michaelides, S., Chrysoulakis, N., Toulios, L., & Clayton, C. (2010b).  Atmospheric correction for satellite remotely sensed data intended for agricultural applications: Impact on vegetation indices.  Natural Hazards and Earth System Sciences, 10. doi:10.5194/nhess-10-89-2010.
  * Kanan & Cottrell [2012] Kanan, C., & Cottrell, G. (2012).  Color-to-grayscale: Does the method matter in image recognition?  PloS one, 7, e29740. doi:10.1371/journal.pone.0029740.
  * Karthikeyani et al. [2007] Karthikeyani, V., Duraiswamy, K., & Rangasamy, K. (2007).  Conversion of gray-scale image to color image with and without texture synthesis, .
  * Kumar & Verma [2010] Kumar, T., & Verma, K. (2010).  A theory based on conversion of rgb image to gray image.  International Journal of Computer Applications, 7. doi:10.5120/1140-1493.
  * Kusetogullari & Yavariabdi [2018] Kusetogullari, H., & Yavariabdi, A. (2018).  Unsupervised change detection in landsat images with atmospheric artifacts: A fuzzy multiobjective approach.  Mathematical Problems in Engineering, 2018, 1–16. doi:10.1155/2018/7274141.
  * Levin et al. [2004] Levin, A., Lischinski, D., & Weiss, Y. (2004).  Colorization using optimization.  ACM Trans. Graph., 23, 689–694. URL: https://doi.org/10.1145/1015706.1015780. doi:10.1145/1015706.1015780.
  * Lu & Plataniotis [2009] Lu, J., & Plataniotis, K. (2009).  On conversion from color to gray-scale images for face detection.  2012 IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops, 0, 114–119. doi:10.1109/CVPRW.2009.5204297.
  * Ma, and Liu, and Ren, and Luo, Yu-Feng [2019] Ma, and Liu, and Ren, and Luo, Yu-Feng (2019).  Detection of collapsed buildings in post-earthquake remote sensing images based on the improved yolov3.  Remote Sensing, 12, 44. doi:10.3390/rs12010044.
  * Mikołajczyk & Grochowski [2018] Mikołajczyk, A., & Grochowski, M. (2018).  Data augmentation for improving deep learning in image classification problem, .  (pp. 117–122). doi:10.1109/IIPHDW.2018.8388338.
  * Prasad et al. [2015] Prasad, S., Kumar, P., & Sinha, K. (2015).  Grayscale to color map transformation for efficient image analysis on low processing devices.  Advances in Intelligent Systems and Computing, 320, 9–18. doi:10.1007/978-3-319-11218-32.
  * Redmon et al. [2016] Redmon, J., Divvala, S., Girshick, R., & Farhadi, A. (2016).  You only look once: Unified, real-time object detection, .  (pp. 779–788). doi:10.1109/CVPR.2016.91.
  * Redmon & Farhadi [2017] Redmon, J., & Farhadi, A. (2017).  Yolo9000: Better, faster, stronger, .  (pp. 6517–6525). doi:10.1109/CVPR.2017.690.
  * Redmon & Farhadi [2018] Redmon, J., & Farhadi, A. (2018).  Yolov3: An incremental improvement.  CoRR, abs/1804.02767. URL: http://arxiv.org/abs/1804.02767. arXiv:1804.02767.
  * Abd Al-salam Selami & Fadhil [2016] Abd Al-salam Selami, A., & Fadhil, A. (2016).  A study of the effects of gaussian noise on image features.  Kirkuk University Journal / Scientific Studies (1992-0849), 11, 152 -- 169. doi:10.32894/kujss.2016.124648.
  * Tempfli et al. [2009] Tempfli, K., Huurneman, G., Bakker, W., & Janssen, L. (2009).  Principles of remote sensing : an introductory textbook.
  * Wallace [1983] Wallace, A. M. (1983).  Greyscale image processing for industrial applications.  Image and Vision Computing, 1, 178 -- 188. URL: http://www.sciencedirect.com/science/article/pii/0262885683900173. doi:https://doi.org/10.1016/0262-8856(83)90017-3.
  * Welsh et al. [2002] Welsh, T., Ashikhmin, M., & Mueller, K. (2002).  Transferring color to greyscale images.  ACM Trans. Graph., 21, 277--280. doi:10.1145/566570.566576.
  * Xia et al. [2018] Xia, M., Liu, X., & Wong, T. (2018).  Invertible grayscale.  ACM Transactions on Graphics (TOG), 37, 1 -- 10.