# Simulating terahertz field-induced transient ferroeletricity in quantum
paraelectric SrTiO3

Dongbin Shin<EMAIL_ADDRESS>Max Planck Institute for the Structure
and Dynamics of Matter and Center for Free Electron Laser Science, 22761
Hamburg, Germany  Simone Latini  Max Planck Institute for the Structure and
Dynamics of Matter and Center for Free Electron Laser Science, 22761 Hamburg,
Germany  Christian Schäfer  Max Planck Institute for the Structure and
Dynamics of Matter and Center for Free Electron Laser Science, 22761 Hamburg,
Germany  Shunsuke A. Sato  Max Planck Institute for the Structure and Dynamics
of Matter and Center for Free Electron Laser Science, 22761 Hamburg, Germany
Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577,
Japan  Edoardo Baldini  Department of Physics, Massachusetts Institute of
Technology, 02139 Cambridge, Massachusetts, USA  Umberto De Giovannini  Max
Planck Institute for the Structure and Dynamics of Matter and Center for Free
Electron Laser Science, 22761 Hamburg, Germany  Nano-Bio Spectroscopy Group,
Departamento de Fisica de Materiales, Universidad del País Vasco UPV/EHU-
20018 San Sebastián, Spain  Hannes Hübener  Max Planck Institute for the
Structure and Dynamics of Matter and Center for Free Electron Laser Science,
22761 Hamburg, Germany  Angel Rubio<EMAIL_ADDRESS>Max Planck
Institute for the Structure and Dynamics of Matter and Center for Free
Electron Laser Science, 22761 Hamburg, Germany  Nano-Bio Spectroscopy Group,
Departamento de Fisica de Materiales, Universidad del País Vasco UPV/EHU-
20018 San Sebastián, Spain  Center for Computational Quantum Physics (CCQ),
The Flatiron Institute, 162 Fifth avenue, New York NY 10010.

###### Abstract

Recent experiments have demonstrated that intense terahertz (THz) fields can
induce a transition from the quantum paraeletric to the ferroeletric phase of
SrTiO3. Here, we investigate this THz field-induced transient ferroeletric
phase transition by solving the time-dependent lattice Schödinger equation
based on first-principles calculations. We find that transient ferroeletricity
originates from a light-induced mixing between ground and first excited
lattice states in the quantum paraeletric phase. In agreement with the
experimental findings, our study shows that the non-oscillatory second
harmonic generation signal can be evidence of transient ferroeletricity in
SrTiO3. We reveal the microscopic details of this exotic phase transition and
highlight that this phenomenon is a unique behavior of the quantum paraeletric
phase.

Ultrafast light-induced dynamics provides attractive physical phenomena beyond
the steady state of materials. Recent studies have shown that driving
materials out of their original equilibrium state can induce complex phase
transitions or largely modify properties and functionalities of materials [1,
2, 3, 4, 5, 6]. Theoretical works have suggested that many of these THz field-
induced phenomena are closely related with nonlinear phonon interactions [7,
8, 9, 10]. For example, by applying a terahertz (THz) pulse, a high
temperature (100 K) superconducting state has been experimentally observed in
K3C60 which exhibits low critical temperature (30 K) at its ground state [1,
11, 12]. Possible topological phase transitions induced by THz pulses has also
been demonstrated in layered materials such as 1T-WTe2, ZrTe5, and graphene
[13, 14, 15]. Another paradigmatic manifestation of the THz-matter interaction
is the transformation of paraeletric SrTiO3 into ferroeletric phase [3, 4].

Unlike other transition metal perovskite materials as BaTiO3 and PbTiO3, at
low temperature SrTiO3 does not show a paraeletric-to-ferroeletric phase
transition [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. The origin of this
phenomenon lies in the potential energy surface of SrTiO3 along the
ferroeletric soft (FES) mode, which has a shallow double-well potential at the
low temperature [24, 25, 26]. Despite classically favoring a ferroelectric
phase, SrTiO3 remains in the high symmetry phase due to nuclear quantum
fluctuations [25, 26]. The corresponding low-temperature phase of SrTiO3 is
therefore a quantum paraelectric phase. Recent first-principles calculations
revealed that both quantum fluctuation and the nonlinear interaction between
the FES mode and lattice strain are required to describe the quantum
paraeletric phase in SrTiO3 [26]. The inclusion of these effects enables to
computationally reproduce temperature dependent properties, such as the energy
of the FES mode and the dielectric constant [26, 24]. Relying on such a first
principles description, an optically stabilized ferroeletric ground state of
SrTiO3 in a cavity has been theoretically demonstrated [27].

Here, we investigate the THz field-induced non-equilibrium transient
ferroeletricity of the quantum paraeletric state of SrTiO3 based on a lattice
model derived from extensive first principles calculations [26]. By solving
the time-dependent lattice Schödinger-Langevin equation, we show that a
transient ferroeletric state can be obtained in SrTiO3 by illumination with a
single cycle THz pulse. We find that this transient ferroelectricity
originates from the mixing between ground and first excited states of the
lattice wavefunction. Based on time-dependent density functional theory
(TDDFT) calculations, we obtain a non-oscillatory second harmonic generation
(SHG) signal and find that it is a signature of the induced broken symmetry,
in agreement with the conclusions drawn from experimental observations [3, 4].
Lastly, we discuss the conditions of THz pulse frequency and dissipation rate
for the transient ferroeletricity.

Figure 1:  THz field-induced transient ferroeletricity in SrTiO3. (a)
Schematic image of the THz pulse applied to SrTiO3. (b) Schematic image of the
atomic displacements involved in the FES mode ($Q_{f}$) and lattice strain
($Q_{c}$). (c) Two-dimensional potential energy surface of SrTiO3 along the
FES mode ($Q_{f}$) and lattice ($Q_{c}$) coordinates. (d) Time profile of
single-cycle THz pulse with $0.5$ THz frequency. (e) Time profiles of the
expectation value of the FES mode ($Q_{f}$) and lattice ($Q_{c}$) perturbed by
the THz pulse. (f) Variation of the time-averaged lattice density perturbed by
the THz pulse. Red dashed line in (e) is a guide to the eye for the zero
value.

By solving a time-dependent lattice Schrödinger-Langevin equation, we
investigate the THz field-induced transient ferroeletric transition from
quantum paraeletric SrTiO3 as depicted in Fig. 1(a). Recent first principles-
based study reported that the quantum paraeletric state in SrTiO3 can be
described by the two-dimensional (2D) lattice Schrödinger equation [26]. The
2D effective potential for quantum paraeletric SrTiO3 consists of FES mode
$Q_{f}$ and the lattice expansion/contraction along the FES mode direction
$Q_{c}$, as shown in Figs. 1(b) and 1(c). It is constructed from fitting the
calculated total energy as a function of geometry distortions as follows:
$\hat{V}_{\text{2D}}^{\text{FES,c}}=\sum_{i=1}^{12}k_{f,i}\hat{Q}^{2i}_{f}+\sum_{j=2}^{10}k_{c,j}\hat{Q}^{j}_{c}+\sum_{i=1}^{12}\sum_{j=1}^{10}k_{fc,i,j}\hat{Q}^{2i}_{f}\hat{Q}^{j}_{c}$
[26, 28]. For the DFT total energy calculation, we used the QUANTUM ESPRESSO
package with the projector augmented wave method and a plane-wave basis set
with 70 Ry energy cut-off [29]. The electron-electron exchange and correlation
potentials are described through the Perdew-Berke-Ernzerhof functional [30].
We considered the $\sqrt{2}\times\sqrt{2}\times 2$ tetragonal unit cell for
the low temperature phase of SrTiO3 and we sampled the Brillouin zone with a
$6\times 6\times 4$ $\mathbf{k}$-point grid. Details on the parameters used
the 2D effective potential energy surface are summarized in the Supplemental
Material [28]. To describe the effect of the THz pump pulse on quantum
paraeletric SrTiO3, we solved the time-dependent 2D lattice Schrödinger-
Langevin equation
$i\hbar\frac{d}{dt}\psi=\hat{H}_{\text{2D}}\psi+\gamma(\hat{S}-<\hat{S}>)\psi$,
with the Hamiltonian
$\hat{H}_{\text{2D}}[Q_{f},Q_{c},t]={\hat{P}^{2}_{f}}/{2M_{f}}+{\hat{P}^{2}_{c}}/{2M_{c}}+\hat{V}_{\text{2D}}^{\text{FES,c}}+E_{ext}(t)Z_{f}\hat{Q}_{f}$.
The effective mass of the FES mode is $M_{f}=1.76\times 10^{-25}$ kg. The mass
of the tetragonal cell is used as the lattice effective mass
($M_{c}=\sum_{i}M_{i}=1.22\times 10^{-24}$ kg). The $\gamma$ and
$\hat{S}[Q_{f},Q_{c},t]=arg[\psi(Q_{f},Q_{c},t)]$ are the dissipation rate and
the phase of the wavefunction respectively [31, 32]. We used a dissipation
rate of $\gamma=0.9$ THz for the FES mode as reported by recent spectroscopic
measurements [33]. We considered the interaction between the FES mode and the
THz pulse ($H_{int}=E_{ext}(t)Z_{f}\hat{Q}_{f}$) to be linear in the external
field strength and the electric dipole generated by $Q_{f}$ with its Born
effective charge ($Z_{f}=36e$) calculated under the dipole approximation.
Details on the time-dependent Schrödinger-Langevin equation are discussed in
the Supplemental Material [28].

We simulated the action of a single-cycle THz pulse on the low temperature
quantum paraeletric state of SrTiO3, choosing it to be resonant with the
frequency of the FES mode ($0.5$ THz) and with a maximum strength of $375$
kV/cm, as shown in Fig. 1(d). This THz pulse perturbs the quantum paraeletric
ground state $\psi_{0}(Q_{f},Q_{c},t=0)$ and induces a variation in the
expectation values of $Q_{f}(t)$ and $Q_{c}(t)$, as shown in Fig. 1(e). It is
notable that the asymmetric potential along $Q_{c}$ provides non-zero value
($Q_{c}(t=0)=31$ mÅ) in the ground state because of quantum fluctuations [26].
The time-profile of the FES mode exhibits a non-oscillatory component between
$7.5$ and $15$ ps time period. Because the FES mode is directly related to the
electric polarization of the system ($P=Z_{f}Q_{f}$), this plateau indicates
that the single-cycle THz pulse induces a transient ferroeletric phase. We
also found that the THz pulse transiently elongates the lattice ($Q_{c}$),
providing a deeper double well for the FES mode, as shown in Fig. 1(c) [26].
To understand the variation of the 2D lattice state, we evaluated the time-
averaged lattice densities
($\bar{n}(t_{1},t_{0})=\int_{t_{0}}^{t_{1}}|\psi(Q_{f},Q_{t},t)|^{2}dt$) for a
given period of time, as displayed in Fig. 1(f). In the ground state ($t=0$
ps), the time-averaged lattice density is delocalized over the double well
potential along the $Q_{f}$ direction, indicating a quantum paraeletric ground
state. At the time delay $t=3.5\sim 7$ ps, the applied THz pulse strongly
perturbs the lattice density. Upon application of the THz pulse ($t=7.5\sim
14$ ps), the lattice density localizes in one of the two wells signaling the
emergence of a transient ferroeletric state. Due to dissipation, the excited
state finally returns to the quantum paraeletric phase ($t=15\sim 20$ ps).

Figure 2:  Mechanism of the THz field-induced phase transition from the
quantum paraeletric state. (a) Time profiles of expectation value of FES mode
($Q_{f}$) perturbed by a THz pulse with $375$ kV/cm (top panel) and $500$
kV/cm (bottom panel) field strength. (b) Schematic image of the phase
transition from the quantum paraeletric state. (c) Ground and (d) first
excited states of the quantum paraeletric phase in SrTiO3. Time profile of the
SHG signal induced by a THz pulse with (e) $375$ kV/cm and (f) $500$ kV/cm.
The red dashed line in (a) is a guide to the eye for the zero level. The red
dashed line in (E) is a guide to the eye for range of the transient
ferroeletric stat. The inset of (f) shows the relation between the SHG signal
and the FES mode distortion evaluated by TDDFT calculations.

This THz field-induced ferroeletric transition originates from the excitation
between the ground and first excited states of the lattice wavefunction. As
shown in Fig. 1(e) and upper panel of Fig. 2(a), the THz pulse with $375$
kV/cm induces a transient ferroeletric state that is characterized by the non-
oscillatory behavior of $Q_{f}$ between $8$ and $12.5$ ps. On the other hand,
transient ferroeletricity can be barely triggered by a THz pulse with higher
strength ($500$ kV/cm). The THz-induced ferroeletricity and its strength
dependence can be explained in terms of phonon excitation by the THz pulse, as
displayed in Fig. 2(b). The THz pulse with appropriate strength allows
excitation to the first excited state ($\psi_{1}$) from the ground state
($\psi_{0}$) shown in Figs. 2(c) and 2(d). It is worth noting that the linear
combination between ground and first excited states provides two degenerate
ferroeletric states, which have opposite electric polarization
($\psi_{0}+\psi_{1}$ and $\psi_{0}-\psi_{1}$). When a single-cycle THz pulse
is applied, the symmetry of the excited state along $Q_{f}$ is broken and
ferroeletricity emerges. By contrast, a high field strength THz pulse promotes
transitions to excited states much higher in energy and produces a light-mixed
state equivalent to thermally excited paraeletric state, which spreads over
the full double-well potential. Details about the ferroeletric lattice
wavefunction in the quantum paraeletric phase are discussed in the
Supplemental Material [28].

In recent experimental results [3, 4], the SHG signal of a delayed near-
infrared probe has been interpreted as a signature of the light-induced
ferroeletricity in SrTiO3. To evaluate the SHG signal from the lattice
wavefunction, we calculated the relation between the SHG signal and $Q_{f}$
from TDDFT calculation with distorted geometry along the FES mode direction
(see the inset of Fig. 3(f)) [34]. With tetragonal SrTiO3 geometries distorted
along the FES mode direction, we evaluated the response of the current by
applying a $1.55$ eV probe pulse with $40$ fs duration. To connect the time-
profile of SHG to the 2D model quantities, the obtained SHG signal is fitted
as a function of $Q_{f}$ with polynomial parameters:
$\rm{SHG}(Q_{f})=\sum_{i}^{3}c_{i}|Q_{f}|^{i}$. Relying on this relation, we
estimated the time-profile of the SHG signal, as shown in Figs. 2(e) and 2(f).
The SHG signal induced by the $375$ kV/cm THz pulse exhibits the non-
oscillatory behavior similar to the time-profile of $Q_{f}$. On the other
hand, excitation with THz pulse of higher intensity ($500$ kV/cm) leads to a
fast decay of the SHG signal, resulting in oscillations that touch zero.
Similar to the experimental observations [3, 4], our TDDFT calculation reveals
that the non-oscillatory behavior of the SHG signal can be a key fingerprint
of the transient ferroeletric state. Details on the procedure used to evaluate
the SHG signal from TDDFT calculations are discussed in the Supplemental
Material [28].

Figure 3:  Effects of frequency and dissipation on THz-induced dynamics in
SrTiO3. (a) Time-averaged value of FES displacement and (b) time profiles of
$Q_{f}$ and $Q_{c}$ at various THZ frequencies $\omega$ and $\gamma=0.9$ THz.
(c) Time-averaged value of the FES displacement and (d) time profiles of
$Q_{f}$ and $Q_{c}$ at various dissipation rate $\gamma$ and $\omega=0.5$ THz.
In (a) and (c), the pink arrows indicate the specific points for the time-
profiles in (b) and (d).

Terahertz field-induced transient ferroeletricity can be obtained with a
single resonant THz pulse and realistic dissipation. To investigate the
condition for the THz-induced transient ferroeletricity, we evaluated the
time-averaged $Q_{f}$ value as follows:
$Q_{f}^{avg}=\frac{1}{t_{1}}\int^{t_{0}+t_{1}}_{t_{0}}Q_{f}(t)dt$, where
$t_{0}$ is the time at the end point of the pulse and $t_{1}=12$ ps. First, we
investigated the effect of the THz pulse frequency on transient
ferroeletricity, as shown in Fig. 3(a). We employed a single-cycle THz pulse
with frequency $\omega$, $E(t)=E_{0}cos(2\pi\omega
t)sin^{2}(\frac{2\pi}{10}\omega t)$, as shown in Fig. 1(d). When we apply the
resonant THz pulse ($\omega=0.5$ THz) with the dissipation rate $\gamma=0.9$
THz, we observe non-oscillatory behavior with large values of $Q_{f}^{avg}$
for a wide range of THz field strength. On the other hand, a near-resonant THz
pulse with $\omega=0.8$ THz provides suppressed $Q_{f}^{avg}$ values and an
off-resonant THz pulse with $\omega=1.0$ THz gives negligible $Q_{f}^{avg}$
values. As a consequence, while a non-oscillatory behavior of $Q_{f}(t)$ is
observed with resonant $\omega=0.5$ THz and near-resonant $\omega=0.8$ THz
frequencies, there is no plateau with off-resonant frequency ($\omega=1.0$
THz), as shown in Fig. 3(b). This result indicates that the resonant THz pulse
effectively excites the FES mode from the ground state to the first excited
state, and transient ferroelectricity is formed as a superposition of the
ground and first excited states. In contrast, off-resonant THz pulses hardly
lead to transient ferroelectricity. It is notable that single-cycle THz
pulses, which contain a broad frequency spectrum, allow for the excitation
between ground and first excited states even with near-resonant condition.

We also investigated the effect of the dissipation rate $\gamma$ on the THz-
induced transient ferroeletricity with resonant frequency as shown in Fig.
3(c). Considering that our simplified 2D lattice system does not have other
dissipation channels (e.g., electron-phonon and phonon-phonon interactions),
the experimentally observed dissipation ratio is required to describe
realistic dynamics for this lattice wavefunction. Without dissipation
($\gamma=0$ THz) and with $\gamma=0.6$ THz, SrTiO3 hardly exhibits transient
ferroeletricity but provides thermally excited behavior such as oscillating
$Q_{c}$ (see Fig. 3(d)). On the other hand, the experimentally observed
dissipation rate ($\gamma=0.9$ THz) [33] and higher rates ($\gamma=1.0$ THz)
provide a wide range of field strengths for transient ferroeletricity. These
results indicate that realistic dissipation is needed to obtain transient
ferroelectricity.

Recently, a plethora of experimental and theoretical studies on light-induced
ferroelecticity in perovskite materials have been reported [3, 4, 35, 27, 36].
Among these studies, our results well reproduce the experimental observation
of THz field-induced ferroeletricity in SrTiO3 [3]. In these experiments, a
single-cycle THz pulse was applied on SrTiO3 and resulted in a non-oscillatory
SHG signal at optical frequencies, signature of the transient ferroeletric
phase. Similar to this observation, we also obtained the non-oscillatory SHG
signal and revealed the mechanism underlying this phenomenon. In the same THz-
induced SHG traces as well as in optical depolarization signals, the
experiments also observed coherent oscillations, indicating the excitation of
the anti-ferrodistortive and FES modes in the ferroeletric state. In our
simplified lattice model Hamiltonian, the oscillation of the anti-
ferrodistortive mode is neglected. Considering the relation between the SHG
signal and the FES mode coordinate ($\rm{SHG}(Q_{f})=\rm{SHG}(-Q_{f})$), the
doubled frequency of the FES in the THz-induced SHG signal indicates
paraeletric state which shows the motion of FES mode oscillation over the zero
point ($Q_{f}=0$). We propose that this behavior of frequency doubling can be
experimentally observed under high temperature or high intensity conditions,
leading to a thermally excited paraeletric state rather than ferroeletric
state.

In conclusion, we investigated THz field-induced transient ferroelectricity in
quantum paraelectric SrTiO3. We solved the time-dependent 2D lattice
Schrödinger-Langevin equation to simulate the action of a THz pulse on the
material. We found that transient ferroeletricity is realized upon application
of a THz pulse with moderate field strength and considering a realistic
dissipation rate. The main mechanism of this transient ferroeletricity lies in
the excitation between ground and first excited states of the quantum
paraeletric phase induced by the symmetry-breaking single-cycle THz pulse. In
agreement with previous experimental reports [3, 4], we verified that the non-
oscillatory SHG signal is a signature of transient ferroelectricity. We also
found that the resonant frequency of the THz pulse and realistic dissipation
ratios are required in order to trigger this non-equilibrium phase transition.
Our study not only does reproduce the recent experimental observations, but
also provides microscopic details that emphasize the importance of the quantum
paraeletric phase in this phenomenon.

###### Acknowledgements.

We acknowledge financial support from the European Research Council
(ERC-2015-AdG-694097), Grupos Consolidados (IT1249-19), JSPS KAKENHI Grant
Number 20K14382, the Cluster of Excellence ’CUI: Advanced Imaging of Matter’
of the Deutsche Forschungsgemeinschaft (DFG) - EXC 2056 - project ID
390715994. D.S. and S.L. are supported by Alexander von Humboldt Foundation.
We also acknowledge support from the Max Planck–New York Center for Non-
Equilibrium Quantum Phenomena. The Flatiron Institute is a division of the
Simons Foundation.

## References

  * Mitrano  _et al._ [2016] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. D. Pietro, D. Pontiroli, M. Riccò, A. Subedi, S. R. Clark, D. Jaksch, and A. Cavalleri, Nature 530, 461 (2016).
  * Nova  _et al._ [2017] T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nat. Phys. 13, 132 (2017).
  * Li  _et al._ [2019] X. Li, T. Qiu, J. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, Science 364, 1079 (2019).
  * Nova  _et al._ [2019] T. F. Nova, A. S. Disa, M. Fechner, and A. Cavalleri, Science 364, 1075 (2019).
  * Kogar  _et al._ [2020] A. Kogar, A. Zong, P. E. Dolgirev, X. Shen, J. Straquadine, Y.-Q. Bie, X. Wang, T. Rohwer, I.-C. Tung, Y. Yang, R. Li, J. Yang, S. Weathersby, S. Park, M. E. Kozina, E. J. Sie, H. Wen, P. Jarillo-Herrero, I. R. Fisher, X. Wang, and N. Gedik, Nat. Phys. 16, 159 (2020).
  * Zhang  _et al._ [2019] J. Zhang, C. Lian, M. Guan, W. Ma, H. Fu, H. Guo, and S. Meng, Nano Lett. 19, 6027 (2019).
  * Knap  _et al._ [2016] M. Knap, M. Babadi, G. Refael, I. Martin, and E. Demler, Phys. Rev. B 94, 214504 (2016).
  * Kennes  _et al._ [2017] D. M. Kennes, E. Y. Wilner, D. R. Reichman, and A. J. Millis, Nat, Phys. 13, 479 (2017).
  * Guan  _et al._ [2021] M.-X. Guan, E. Wang, P.-W. You, J.-T. Sun, and S. Meng, Nat. Commun. 12, 1885 (2021).
  * Shin  _et al._ [2020] D. Shin, S. A. Sato, H. Hübener, U. D. Giovannini, N. Park, and A. Rubio, Npj Comput. Mater. 6, 182 (2020).
  * Cantaluppi  _et al._ [2018] A. Cantaluppi, M. Buzzi, G. Jotzu, D. Nicoletti, M. Mitrano, D. Pontiroli, M. Riccò, A. Perucchi, P. D. Pietro, and A. Cavalleri, Nat. Phys. 14, 837 (2018).
  * Budden  _et al._ [2021] M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, F. Schlawin, D. Jaksch, and A. Cavalleri, Nat. Phys. 17, 611 (2021).
  * Sie  _et al._ [2019] E. J. Sie, C. M. Nyby, C. D. Pemmaraju, S. J. Park, X. Shen, J. Yang, M. C. Hoffmann, B. K. Ofori-Okai, R. Li, A. H. Reid, S. Weathersby, E. Mannebach, N. Finney, D. Rhodes, D. Chenet, A. Antony, L. Balicas, J. Hone, T. P. Devereaux, T. F. Heinz, X. Wang, and A. M. Lindenberg, Nature 565, 61 (2019).
  * Vaswani  _et al._ [2020] C. Vaswani, L.-L. Wang, D. H. Mudiyanselage, Q. Li, P. M. Lozano, G. D. Gu, D. Cheng, B. Song, L. Luo, R. H. J. Kim, C. Huang, Z. Liu, M. Mootz, I. E. Perakis, Y. Yao, K. M. Ho, and J. Wang, Phys. Rev. X 10, 021013 (2020).
  * McIver  _et al._ [2020] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri, Nat. Phys. 16, 38 (2020).
  * Shirane  _et al._ [1952] G. Shirane, K. Suzuki, and A. Takeda, J. Phys. Soc. Japan 7, 12 (1952).
  * Samara [1971] G. A. Samara, Ferroelectrics 2, 277 (1971).
  * Bellaiche  _et al._ [2000] L. Bellaiche, A. García, and D. Vanderbilt, Phys. Rev. Lett. 84, 5427 (2000).
  * Guo  _et al._ [2000] R. Guo, L. E. Cross, S. E. Park, B. Noheda, D. E. Cox, and G. Shirane, Phys. Rev. Lett. 84, 5423 (2000).
  * Miyasaka  _et al._ [2003] S. Miyasaka, Y. Okimoto, M. Iwama, and Y. Tokura, Phys. Rev. B 68, 100406(R) (2003).
  * Haeni  _et al._ [2004] J. H. Haeni, P. Irvin, W. Chang, R. Uecker, P. Reiche, Y. L. Li, S. Choudhury, W. Tian, M. E. Hawley, B. Craigo, A. K. Tagantsev, X. Q. Pan, S. K. Streiffer, L. Q. Chen, S. W. Kirchoefer, J. Levy, and D. G. Schlom, Nature 430, 758 (2004).
  * Fennie and Rabe [2006] C. J. Fennie and K. M. Rabe, Phys. Rev. Lett. 97, 267602 (2006).
  * Lee and Rabe [2010] J. H. Lee and K. M. Rabe, Phys. Rev. Lett. 104, 207204 (2010).
  * Müller and Burkard [1979] K. A. Müller and H. Burkard, Phys. Rev. B 19, 3593 (1979).
  * Zhong and Vanderbilt [1995] W. Zhong and D. Vanderbilt, Phys. Rev. Lett. 74, 2587 (1995).
  * Shin  _et al._ [2021] D. Shin, S. Latini, C. Schäfer, S. A. Sato, U. De Giovannini, H. Hübener, and A. Rubio, arXiv preprint arXiv:2101.02291 (2021).
  * Latini  _et al._ [2021] S. Latini, D. Shin, S. A. Sato, C. Schäfer, U. D. Giovannini, H. Hübener, and A. Rubio, The ferroelectric photo-groundstate of srtio3: Cavity materials engineering (2021), arXiv:2101.11313 [cond-mat.mtrl-sci] .
  * [28] See supplemental material at url.
  * Giannozzi  _et al._ [2017] P. Giannozzi _et al._ , J. Condens. Matter Phys. 29, 465901 (2017).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Garashchuk  _et al._ [2013] S. Garashchuk, V. Dixit, B. Gu, and J. Mazzuca, J. Chem. Phys. 138, 054107 (2013).
  * Chou [2015] C.-C. Chou, Ann. Phys. 362, 57 (2015).
  * Kozina  _et al._ [2019] M. Kozina, M. Fechner, P. Marsik, T. v. Driel, J. M. Glownia, C. Bernhard, M. Radovic, D. Zhu, S. Bonetti, U. Staub, and M. C. Hoffmann, Nat. Phys. 15, 387 (2019).
  * Tancogne-Dejean  _et al._ [2020] N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. L. Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, U. D. Giovannini, A. Delgado, F. G. Eich, J. Flick, G. Gil, A. Gomez, N. Helbig, H. Hübener, R. Jestädt, J. Jornet-Somoza, A. H. Larsen, I. V. Lebedeva, M. Lüders, M. A. L. Marques, S. T. Ohlmann, S. Pipolo, M. Rampp, C. A. Rozzi, D. A. Strubbe, S. A. Sato, C. Schäfer, I. Theophilou, A. Welden, and A. Rubio, J. Chem. Phys. 152, 124119 (2020).
  * Ashida  _et al._ [2020] Y. Ashida, A. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, Phys. Rev. X 10, 041027 (2020).
  * Pal  _et al._ [2021] S. Pal, N. Strkalj, C.-J. Yang, M. C. Weber, M. Trassin, M. Woerner, and M. Fiebig, Phys. Rev. X 11, 021023 (2021).


constant $C:=C(\alpha,\beta,p,r,q)>0$, such that, uniformly over $t>0$,

$\|e^{t\Delta}f\|_{\mathcal{B}^{\alpha}_{p,q}}\leq
Ct^{-\frac{1}{2}\left(\frac{1}{r}-\frac{1}{p}\right)-\frac{1}{2}\left(\alpha-\beta\right)}\|f\|_{\mathcal{B}^{\beta}_{r,q}}.$
(A.10)

Secondly, if $0\leq\beta-\alpha\leq 2$ then there exists a
$C:=C(\alpha,\beta,p,q)>0$ such that for a any $t>0$,

$\|(1-e^{t\Delta})f\|_{\mathcal{B}^{\alpha}_{p,q}}\leq
Ct^{\frac{\beta-\alpha}{2}}\|f\|_{\mathcal{B}^{\beta}_{p,q}}.$ (A.11)

######  Remark A.7.

Since we have defined the Hölder-Besov spaces as the closure of
$C^{\infty}(\mathbb{T})$ under the appropriate norm, for any
$\alpha\in\mathbb{R}$, $p,\,q\in[1,\infty]$, if
$f\in\mathcal{B}^{\alpha}_{p,q}(\mathbb{T})$, we have that

$\lim_{t\rightarrow 0}\|(1-e^{t\Delta})f\|_{\mathcal{B}^{\alpha}_{p,q}}=0.$

The equivalent elliptic regularity estimate is an easy consequence of [2, Lem.
2.2] concerning Fourier multipliers.

######  Theorem A.8.

Let $f\in C^{\infty}(\mathbb{T})$ be such that $\langle f,1\rangle=0$. Suppose
$-\partial_{xx}\rho=f$ and $\langle\rho,1\rangle=0$. Then for any
$\alpha\in\mathbb{R}$, $p,q\in[1,\infty]$, there exists a constant $C\geq 0$
such that

$\|\rho\|_{\mathcal{B}^{\alpha}_{p,q}}\leq
C\|f\|_{\mathcal{B}^{\alpha-2}_{p,q}}.$ (A.12)

###### Proof.

From [2, Lem. 2.2], there exists a $C>0$ such that for any $k\geq 0$ we have
$\|\Delta_{k}\rho\|_{L^{p}}\leq C2^{-2k}\|\Delta_{k}f\|_{L^{p}}$. Since
$(\mathcal{F}f)(0)=(\mathcal{F}\rho)(0)=0$,
$\Delta_{-1}\rho(\zeta)=\Delta_{-1}f(\zeta)$. Therefore, there exists a
dimension dependent constant $C>0$, such that for any $k\geq-1$,

$2^{\alpha k}\|\Delta_{k}\rho\|_{L^{p}}\leq
C2^{(\alpha-2)k}\|\Delta_{k}f\|_{L^{p}},$

from which the claim follows. ∎

######  Remark A.9.

Applying Theorem A.8 followed by Theorem A.3 gives us the estimate, for
$f,\rho,\alpha,p,q$ as in Theorem A.8 and $r\in[1,p]$,

$\|\partial_{x}\rho\|_{\mathcal{B}^{\alpha}_{p,q}}\lesssim\|f\|_{\mathcal{B}^{\alpha-1+\left(\frac{1}{r}-\frac{1}{p}\right)}_{r,q}}.$
(A.13)

##  Appendix B Regularity of the Stochastic Heat Equation

We outline the main steps necessary to prove Theorem 2.3 using the
characterisation of Hölder–Besov spaces. We note that a similar result can be
proved using a version of the Garsia–Rodemich–Rumsey lemma, [14]. A more
detailed presentation of similar results in our context can be found in [22,
34]. We include a proof directly in our setting for the reader’s convenience.
A central tool is the following Kolmogorov regularity result which can be
found in a similar form as [22, Lem. 9 & 10]. We recall the kernels $h_{k}\in
C^{\infty}(\mathbb{T})$ from the beginning of Appendix A.

######  Lemma B.1 (Kolmogorov Regularity - [22, Lem. 9 & 10]).

Let $(t,\phi)\mapsto Z_{t}(\phi)$ be a map from $\mathbb{R}_{+}\times
L^{2}(\mathbb{T})\rightarrow L^{2}(\Omega,\mathcal{F},\mathbb{P})$ that is
linear and continuous in $\phi$ and such that $Z_{0}=0$. Assume that for some
$p>1$, $\alpha^{\prime}\in\mathbb{R}$, $\kappa^{\prime}>\frac{1}{p}$ and
$T>0$, there exists a constant $C>0$, such that for all $k\geq-1$,
$x\in\mathbb{T}$ and $s,\,t\in[0,T]$,

$\mathbb{E}\left[|Z_{t}(h_{k}(\,\cdot\,-x))-Z_{s}(h_{k}(\,\cdot\,-x))|^{p}\right]\leq
C^{p}|t-s|^{\kappa^{\prime}p}2^{-k\alpha^{\prime}p}.$ (B.1)

Then, for all $\alpha<\tilde{\alpha}^{\prime}-\frac{1}{p}$ and
$\kappa\in[0,\kappa^{\prime}-\frac{1}{p})$, there exists a modification
$\tilde{Z}$ of $Z$ such that

$\mathbb{E}\left[\|\tilde{Z}\|^{p}_{\mathcal{C}^{\kappa}_{T}\mathcal{B}^{\alpha-\kappa}_{p,p}}\right]<\infty.$

###### Proof.

It suffices to modify the proof of [22, Lem. 9 & 10] to our finite volume and
one dimensional setting. ∎

###### Proof of Theorem 2.3.

Below we show that $\mathcal{L}(\mathcal{I}_{t_{0},t})$ depends only on
$|t-t_{0}|$ so without loss of generality we set $t_{0}=0$. It is also clear,
by translation invariance of $\mathcal{H}$ in space and the law of the white
noise, that $\mathcal{I}_{t_{0},t}$ is translation invariant, i.e
$\mathcal{I}_{t_{0},t}(\phi(\,\cdot\,+x))\sim\mathcal{I}_{t_{0},t}(\phi(\,\cdot\,))$
for any $\phi\in C^{\infty}(\mathbb{T})$ and $x\in\mathbb{T}$.

For $k\geq-1$, $t\in(0,T]$ and setting $x=0$ in (B.1), we have
$v_{0,t}(h_{k})=\mathcal{I}_{0,t}(h_{k})$. So splitting the increment using
the semi-group property, Parseval’s theorem, the covariance formula for the
space-time white noise and the action of the heat kernel in Fourier space we
have for all $0\leq s<t\leq T$

$\displaystyle\mathbb{E}\left[|\mathcal{I}_{0,t}(h_{k})-\mathcal{I}_{0,s}(h_{k})|^{2}\right]$
$\displaystyle\lesssim\sum_{\begin{subarray}{c}m\in\mathrm{supp}(\mathcal{F}h_{k})\\\
m\neq 0\end{subarray}}\frac{1}{8\pi^{2}|m|^{2}}(1-e^{-(t-s)8\pi^{2}|m|^{2}})$
$\displaystyle\lesssim\sum_{\begin{subarray}{c}m\in\mathrm{supp}(\mathcal{F}h_{k})\\\
m\neq 0\end{subarray}}\frac{|t-s|^{\gamma}}{(8\pi^{2}|m|^{2})^{1-\gamma}}$
$\displaystyle\lesssim|t-s|^{\gamma}2^{-k(1-2\gamma)},$

where $\mathcal{F}h_{k}$ denotes the Fourier transform of $h_{k}$ and we may
choose any $\gamma\in[0,1]$. So then using Nelson’s estimate [24, Sec. 1.4.3],
for any $p\geq 2$, there exists a constant $C:=C(p,T)>0$ such that for all
$k\geq-1$ and $\gamma\in[0,1)$,

$\mathbb{E}\left[|\mathcal{I}_{0,t}(h_{k})-\mathcal{I}_{0,s}(h_{k})|^{p}\right]\leq
C|t-s|^{\frac{\gamma}{2}p}2^{-k\left(\frac{1}{2}-\gamma\right)p}.$

Then applying Theorem B.1 and the Besov embedding (A.3), we see that there
exists a modification of $t\mapsto v_{0,t}$ (which we do not relabel) such
that for any $p\geq 1$, $\alpha\in(0,1/2),\,\kappa\in\left[0,1/2\right)$ and
$T>0$,

$\mathbb{E}\left[\|v_{0,\,\cdot\,}\|^{p}_{\mathcal{C}^{\kappa}_{T}\mathcal{C}^{\alpha-2\kappa}}\right]<\infty.$

To see stationarity of the processes $t_{0}\mapsto v_{t_{0},t_{0}+h}$ for
fixed $h>0$, we use the fact that $v_{t_{0},t_{0}+h}(\phi)$ is Gaussian, and
has zero mean, so that its law is entirely determined by its second moment.
Therefore, letting $\phi\in C^{\infty}(\mathbb{T})$, and by similar
computations as above, we see that

$\displaystyle\mathbb{E}\left[v_{t_{0},t_{0}+h}(\phi)^{2}\right]$
$\displaystyle=\sum_{m}\frac{1}{8\pi^{2}|m|^{2}}\left(1-e^{-8\pi^{2}|m|^{2}h}\right)|\hat{\phi}_{m}|^{2},$

and hence $\mathcal{L}(v_{t_{0},t_{0}+h})$ depends only on
$h\in[t_{0},T-t_{0}]$. Finally, by the assumption that $\xi$ is spatially mean
free, for any $t_{0}\in[0,T)$ and $h\in(T-t_{0},T]$,
$\mathcal{L}(v_{t_{0},t_{0}+h})\in C^{\alpha}_{0}(\mathbb{T})$. ∎

## References

  * [1] A. Ambrosetti and G. Prodi, A Primer of Nonlinear Analysis, Cambridge University Press, 1993.
  * [2] H. Bahouri, J.-Y. Chemin, and R. Danchin, Fourier Analysis and Nonlinear Partial Differential Equations, Springer, 2011.
  * [3] V. I. Bogachev, Gaussian measures, vol. 62 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 1998.
  * [4] G. Da Prato and A. Debussche, Strong solutions to the stochastic quantization equations, Ann. Probab., 31 (2003), pp. 1900–1916.
  * [5] G. Da Prato and J. Zabczyk, Ergodicity for Infinite Dimensional Systems, London Mathematical Society Lecture Note Series, Cambridge University Press, 1996.
  * [6] G. Da Prato and J. Zabczyk, Stochastic Equations in Infinite Dimensions, Encyclopedia of Mathematics and its Applications, Cambridge University Press, 2 ed., 2014.
  * [7] R. M. Ford and D. A. Lauffenburger, Analysis of chemotactic bacterial distributions in population migration assays using a mathematical model applicable to steep or shallow attractant gradients, Bulletin of Mathematical Biology, 53 (1991), pp. 721–749.
  * [8] M. Hairer, A theory of regularity structures, Inventiones mathematicae, 198 (2014), pp. 269–504.
  * [9] T. Hillen and K. Painter, Global existence for a parabolic chemotaxis model with prevention of overcrowding, Advances in Applied Mathematics, 26 (2001), pp. 280–301.
  * [10] T. Hillen and K. J. Painter, Volume-filling and quorum-sensing in models for chemosensitive movement, Can. Appl. Math. Q., 10 (2002), pp. 501–543.
  * [11] T. Hillen and K. J. Painter, A user’s guide to PDE models for chemotaxis, J. Math. Biol., 58 (2009), pp. 183–217.
  * [12] T. Hillen and A. Potapov, The one-dimensional chemotaxis model: global existence and asymptotic profile, Mathematical Methods in the Applied Sciences, 27 (2004), pp. 1783–1801.
  * [13] D. Horstmann and M. Winkler, Boundedness vs. blow-up in a chemotaxis system, Journal of Differential Equations, 215 (2005), pp. 52–107.
  * [14] Y. Hu and K. Le, A multiparameter Garsia–Rodemich–Rumsey inequality and some applications, Stochastic Process. Appl., 123 (2013), pp. 3359–3377.
  * [15] E. F. Keller and L. A. Segel, Model for chemotaxis, Journal of Theoretical Biology, 30 (1971), pp. 225 – 234.
  * [16] , Traveling bands of chemotactic bacteria: A theoretical analysis, Journal of Theoretical Biology, 30 (1971), pp. 235–248.
  * [17] Y. Lai and Y. Xiao, Existence and asymptotic behavior of global solutions to chemorepulsion systems with nonlinear sensitivity, Electron. J. Differential Equations, (2017), pp. Paper No. 254, 9.
  * [18] I. R. Lapidus and R. Schiller, Model for the chemotactic response of a bacterial population, Biophysical journal, 16 (1976), pp. 779–789.
  * [19] A. Mayorcas, Stochastic PDEs and weakly interacting particle systems, PhD thesis, University of Oxford, Oxford, UK, Dec. 2020.
  * [20] A. Moinat and H. Weber, Local bounds for stochastic reaction diffusion equations, Electron. J. Probab., 25 (2020), p. 26 pp.
  * [21] J.-C. Mourrat and H. Weber, The dynamic $\Phi^{4}_{3}$ model comes down from infinity, Comm. Math. Phys., 356 (2017), pp. 673–753.
  * [22] , Global well-posedness of the dynamic $\phi^{4}$ model in the plane, The Annals of Probability, 45 (2017), pp. 2398–2476.
  * [23] J.-C. Mourrat, H. Weber, and W. Xu, Construction of $\Phi^{4}_{3}$ diagrams for pedestrians.  arxiv.org/abs/1610.08897, 2017.
  * [24] D. Nualart, Malliavin Calculus and Related Topics, Springer, 2010\.
  * [25] K. Osaki and A. Yagi, Finite dimensional attractor for one-dimensional Keller–Segel equations, Funkcial. Ekvac., 44 (2001), pp. 441–469.
  * [26] H. G. Othmer and A. Stevens, Aggregation, blowup, and collapse: the ABCs of taxis in reinforced random walks, SIAM J. Appl. Math., 57 (1997), pp. 1044–1081.
  * [27] F. Otto, H. Weber, and M. G. Westdickenberg, Invariant measure of the stochastic Allen–Cahn equation: the regime of small noise and large system size, Electron. J. Probab., 19 (2014), pp. no. 23, 76.
  * [28] C. S. Patlak, Random walk with persistence and external bias, Bull. Math. Biophys., 15 (1953), pp. 311–338.
  * [29] B. Perthame, PDE models for chemotactic movements: parabolic, hyperbolic and kinetic, Appl. Math., 49 (2004), pp. 539–564.
  * [30] D. Revuz and M. Yor, Continuous Martingales and Brownian Motion, Springer, 2008.
  * [31] Y. Tao, Global dynamics in a higher-dimensional repulsion chemotaxis model with nonlinear sensitivity, Discrete Contin. Dyn. Syst. Ser. B, 18 (2013), pp. 2705–2722.
  * [32] Y. Tao and M. Winkler, Boundedness in a quasilinear parabolic-parabolic Keller–Segel system with subcritical sensitivity, Journal of Differential Equations, 252 (2012), pp. 692–715.
  * [33] C. M.-R. Tomasz Cieślak, Philippe Laurençot, Global existence and convergence to steady states in a chemorepulsion system, Banach Center Publications, 81 (2008), pp. 105–117.
  * [34] P. Tsatsoulis and H. Weber, Spectral gap for the stochastic quantization equation on the 2-dimensional torus, Ann. Inst. H. Poincaré Probab. Statist., 54 (2018), pp. 1204–1249.


# Tiny-Inception-ResNet-v2: Using Deep Learning for Eliminating Bonded Labors
of Brick Kilns in South Asia

Usman Nazir, Numan Khurshid, Muhammad Ahmed Bhimra, Murtaza Taj
Department of Computer Science, Syed Babar Ali School of Science and
Engineering
Lahore University of Management Sciences (LUMS), Lahore, Pakistan
{17030059, 15060051, 17030015<EMAIL_ADDRESS>
###### Abstract

This paper proposes to employ a Inception-ResNet inspired deep learning
architecture called _Tiny-Inception-ResNet-v2_ to eliminate bonded labor by
identifying brick kilns within “Brick-Kiln-Belt” of South Asia. The framework
is developed by training a network on the satellite imagery consisting of 11
different classes of South Asian region. The dataset developed during the
process includes the geo-referenced images of brick kilns, houses, roads,
tennis courts, farms, sparse trees, dense trees, orchards, parking lots, parks
and barren lands. The dataset111https://cvlab.lums.edu.pk/?p=1779 is made
publicly available for further research. Our proposed network architecture
with very fewer learning parameters outperforms all state-of-the-art
architectures employed for recognition of brick kilns. Our proposed solution
would enable regional monitoring and evaluation mechanisms for the Sustainable
Development Goals.

##  1 Introduction

According to the Global Slavery Index of 2018 $24.9$ million people are
trapped in forced labor globally [1]. An estimated $12.7$ million are within
“Brick-Kiln-Belt” of South Asia ($1,551,997~{}km^{2}$ between Afghanistan,
Pakistan, India, Bangladesh and Nepal), because indebted labor is a major
contributing factor to the prevalence of slavery in the Asia-Pacific
constituting 55% of all forced labor-mapping brick kilns and estimating their
labor force is an essential step in eliminating slavery from the region. The
UN’s Sustainable Development Goal (SDGs) 8 specifically refers to forced labor
and governs the world over aims for ending modern slavery by 2025 [2]. They
are, however, faced with a lack of access to reliable, up-to-date and
actionable data on slavery activity. Such data needs to be spatially explicit
and scalable to allow governments to monitor activities and implement
strategies to emancipate individuals trapped within institutions of slavery.
The absence and insufficiency of data compromises evidence-based action and
policy formulation, as has been highlighted by various sustainable development
data-gap analysis reports. Thus, to meet this challenge, new and innovative
approaches are needed.

|  |
---|---|---
|  |
|  |
Figure 1: Example satellite imagery of brick kilns from different spatial
locations showing variation in quality, structure and color profile Figure 2:
Inception-ResNet-V2 contains 11 A-blocks, 20 B-blocks, and 10 C-blocks and our
proposed tiny version only contains 10 A-blocks, 3 B-blocks, and 3 C-blocks

One of the key challenges faced by all development efforts is selecting target
populations, or the beneficiaries of programs. In order to ensure that
interventions reach those who need them, and have an impact on economic and
human development, one must first identify those who are vulnerable and most
in need of help. This problem only increases in complexity when we seek to
identify individuals trapped in slavery or forced labor. In the past,
household surveys have generally been the norm in identifying those in need of
development assistance. However, surveys are costly and cumbersome, and rarely
have large enough samples to be representative at the village or town level.
Data is thus usually aggregated at the national, regional, provincial, or
district level. Statistics estimated at the national or district level,
however, may mask the fact that there are pockets of space which are severely
in need of assistance. In addition to this, slave owners may hinder data
collection in the areas that they control, or enslaved individuals may face
numerous pressures in disclosing themselves as enslaved to enumerators or
authorities.

Remote sensing is one such method that has been applied to various social and
developing issues to help better inform policy. Remote sensing refers to the
use of satellite or other sensor technologies to obtain information about the
physical characteristics of an area from afar without conducting expensive and
strenuous surveys. Moreover, new and improved techniques in computer vision
and deep learning use large datasets with somewhat similar proxy variables as
are collected in surveys to gauge vulnerability on a more minute scale and
provides solution with minimal human intervention.

To develop an efficient automated solution for brick kiln identification, we
propose Tiny-Inception-ResNet-v2 network applied to satellite imagery. We also
developed a new satellite image dataset comprising the images of the South
Asian region for 11 different categories. The dataset will be made publicly
available and will help inform the policy of governments within the region,
and allow them to gauge the success of past interventions, as well as serve as
a baseline by which to assess future progress. To the best of our knowledge,
this will be the first of its kind effort that goes beyond the geographic
boundaries of a specific country and aims to provide a comprehensive survey of
bonded labor at kilns along with their accurate geolocation.

###  1.1 Related Work

Availability of high resolution satellite imagery along with recent
advancements in machine learning particularly deep convolutional neural
networks (CNN) have paved a way for large scale analysis of wide variety of
parameters across the globe. High resolution satellite imagery has been used
by [3] to estimate per capita expenditure in Africa. Recently, remote sensing
images have also been used to analyze the extent of modern slavery. For
instance, [1] have mapped fish farms suspected of using child slave labor in
the Sundarbans Reserve Forest in Bangladesh. Similarly, [1],[4], through their
“Slavery from Space” project proposed a crowd-sourced procedure to manually
detect brick kilns from satellite imagery. However, they were only able to
manually annotate 320 geographic cells (i.e. only 2% of the entire Brick-Kiln-
Belt). This study also does not account for workforce size of each kiln.

The potential inherent within remote sensing and machine learning for social
science research and development purposes has been noted by many [1, 5, 6, 7].
Residual Neural Network [8] have been shown to produce favorable results on
satellite imagery for the problem of classifying imagery into known classes
such as road, houses, vegetation etc. Similarly, [9] proposed algorithm to
classify each pixel of satellite imagery into known classes thus obtaining an
exact boundary of the object of interest in the given satellite image.

One area that has especially benefited from such techniques is poverty
studies. Machine learning has been used to identify poverty stricken
populations and their geographic location through the analysis of high
resolution satellite images. One such technique utilizes the premise that
nighttime luminosity is strongly correlated (if noisily) with economic
activity and development [10]. Using satellite images, such efforts compare
daytime and nighttime images of a specified location (a country, a continent)
to identify which factors are associated with greater night lighting e.g.
paved road networks. [10] find in their study that their method of estimating
consumption expenditure and asset wealth through machine learning is able to
explain 75% of the variation in local level economic outcomes. Though such
methodologies cannot detect poverty at the level of the individual or
household, they can inform policy makers about which sub-populations or
communities are economically marginalized or vulnerable, and help them target
their interventions at particular geographical locations where they are most
needed. Recent contribution of [2] analyzed various machine learning
techniques to automate the process of brick kiln identification in the given
tile of images. However, their approach outputs very high false positive rate.
To cope with the scenario they proposed to train a two-staged R-CNN classifier
to achieve acceptable performance.

Our approach contributes by proposing Tiny-Inception-ResNet-v2, an improved
architecture inspired from Inception-ResNet-v2 [11] to identify kilns in the
satellite images. We showed that with fewer learnable parameters of Tiny-
Inception-ResNet-v2 we could easily classify brick kilns with about $95\%$
accuracy.

|  |
---|---|---
Barren Land | Sparse Trees | Brick Kiln
|  |
Dense Trees | Farms | Grass
|  |
Ground | Houses | Orchard
|  |
Parking | Roads |
Figure 3: Example satellite image from every class in dataset (Satellite image
courtesy Google Maps)

##  2 Challenges

Identifying brick kilns from satellite imagery involves the following
challenges:

###  2.1 Structural Variations

Brick kilns may significantly vary in structure, shape, and size as shown in
Fig. 1. They come in circular, oval, rectangular, and elongated structures.
These variations are usually enforced by geographical locations, environmental
conditions, manufacturing technologies, building material, and local building
regulatory authorities. This is considered as one of the obvious challenges,
for any machine learning approach, to learn a generic set of visual features
for each type of kilns.

###  2.2 Environmental Variations

Satellite images are prone to atmospheric variations including but not limited
to cloud cover, pollution, variation in luminosity, and seasonal changes in
the environment for which the imagery has been acquired. These variations in
images contributes in confusing the classifier. An obvious effect of these
changes could be observed between the top-left and top-right images of Fig. 1.

###  2.3 Sensor Variations

Satellite images are acquired through various imaging devices and sensors
installed on satellites. Some of the famous ones include WorldView, Pleiades,
GeoEye-1, and QuickBird. Sensor variations are usually caused by difference in
parameters of imaging devices capturing same imagery differently over time.
When spatial analysis is spread across multiple cities, these changes in
sensors can be observed as huge variations in the quality, resolution and
color profile of the imagery due to the different satellite data (see Fig. 1).

Figure 4: Inception Block

##  3 Proposed Method

###  3.1 Dataset Generation and Pre-processing

We developed a dataset of satellite imagery from South Asian region by
acquiring images at zoom level 20 with their geo-locations from OpenStreetMaps
(OSM). These images were then manually labelled through crowd-sourcing into 11
different classes each having around $600$ images. All these $6827$ images
were used to train the proposed classifier for categories of kilns, houses,
road, tennis courts, farms, sparse trees, dense trees, orchards, parking,
parks and barren lands (see Fig. 3).

To identify brick kilns in the ”Brick-Kiln-Belt” we used the same network,
testing it on an entirely different images from digital globe imagery
($787,100$ images) for the entire Afghanistan-to-Nepal region
$1,551,997~{}km^{2}$ for the year 2018. These images acquired at zoom-level
$17$ ($1.193~{}\frac{meters}{pixel}$ on Equator) were first converted to zoom-
level $20$ (about $50$ million images) using following pre-processing steps:

  1. 1.

Finding the midpoint _latitude-longitude_ (lat-lon) of each image tile of
zoom-level 17 using following equations:

$lat^{z_{17}}=\Bigg{[}\Big{(}(Y_{tile}^{z_{17}}-A)\times B\Big{)}+C\Bigg{]},$
(1)

where

     * •

$A=42295$ is a reference tile #. It can be chosen randomly; however, we
selected $42295$ throughout our experiments.

     * •

$B=0.0054931640625$ is difference between adjacent tiles of zoom-level 17.

     * •

$C=52.33612060546875$ is a relative lat and lon with respect to $A$.

$lon^{z_{17}}=\Bigg{[}\Big{(}(X_{tile}^{z_{17}}-A)\times B\Big{)}+C\Bigg{]},$
(2)

where $A$, $B$ and $C$ are given above.

  2. 2.

Converting zoom-level $17$ lat-lon to zoom-level 20 lat-lon using following
equations:

$cornerLat^{z_{20}}=lat^{z_{17}}+4*D,$ (3)
$cornerLon^{z_{20}}=lon^{z_{17}}-4*D,$ (4)

where $D=6.866455078125*10^{-4}$ is the difference between lat/lon of adjacent
zoom-level $20$ tiles.

In our work each image is associated with its geo-location which makes it very
easy for images classified as brick kilns class to be identified on the maps.
The geo-location of the identified brick kiln could then be linked to the
census data provided by Punjab-Brick-Kiln-Census-dataset
222http://202.166.167.115/brick_kiln_dashboard/ which is publicly available
containing all the relevant information regarding number of workers, number of
children, age of each child and schooling status associated with each kiln.

We proposed a classifier that is inspired from the well-known Inception-
ResNet-v2 netowrk. The performance gain of the network as compare to other
deep learning architectures encompassing the inception block incorporated in
the network.

###  3.2 Inception Block

The inception module consists of convolutions of different sizes that allow
the network to process features at different spatial scales. They are then
lumped and fed to the next layer for further processing as shown in Fig. 4.
For dimensionality reduction, 1x1 convolutions are used before the more
expensive 3x3 and 5x5 convolutions. In many remote-sensing problems, we need
the deeper network to process features at different spatial scales. To cope
with our challenges, such flexibility can be incorporated in convolutional
neural networks by introducing inception blocks.

###  3.3 Tiny-Inception-ResNet-v2

An improved form of Inception-ResNet-v2 called as ”Tiny-Inception-ResNet-v2”
with much less inception blocks have been proposed to classify satellite
images. We suppose that even fewer number of inception blocks used in
Inception-ResNet network are able to classify the satellite images in classes.
For this purpose optimal number of inception blocks are to be identified. Fig.
2 shows the detailed architecture of our proposed network with a block level
architecture at the top of the figure while the expended version at the
bottom. Three main blocks namely A, B, and C containing different number of
stacked inception blocks are shown. To identify the optimal number of
inception modules to be used in each block, we tried various combinations
given in Table 1 and found that network with $10$ inception block in A, $3$
blocks in B, and $3$ blocks in C outperforms the other combinations. It has
also been observed that this network has much fewer learnable parameters as
compare to other combinations. More inception blocks have been used in A than
B and C, to learn the features of the image responsible for visual appearance.
The later ones focus on discriminating the learned features.

Table 1: Comparison between different versions of Inception-ResNet-v2

A-block | B-block | C-block | Val Loss | Precison | Recall | F1 Score | Parameters
---|---|---|---|---|---|---|---
11 | 20 | 10 | 0.00037 | 0.9950 | 0.8550 | 0.9200 | 54,353,643
20 | 5 | 5 | 0.00016 | 0.9940 | 0.2990 | 0.4597 | 27,243,579
10 | 3 | 3 | 0.00022 | 0.9854 | 0.9052 | 0.9435 | 19,682,011
10 | 1 | 1 | 0.00041 | 0.993 | 0.3130 | 0.4760 | 13,354,523

Figure 5: Qualitative Analysis on Kot Radha Kishan, in Punjab Province,
Pakistan. It can be seen within the green region-of-interest; proposed network
correctly classified all brick kilns with threshold 0.5. (Satellite image
courtesy Google Maps) Figure 6: Qualitative Analysis on Province Punjab,
India. It can be seen within the green region-of-interest; proposed network
correctly classified all brick kilns with threshold 0.9. (Satellite image
courtesy Google Maps)

##  4 Results and Evaluation

We evaluated our trained network of Tiny-Inception-ResNet-v2 on unseen dataset
consisting of 700 images of zoom 20 for quantitative analysis and 787,100
images of zoom 17 acquired from Digital Globe Imagery for qualitative
analysis. The classifier with output kiln is considered as 1 whereas the rest
of the class outputs are considered as non-kiln (0). Following benchmark
performance metrics have been calculated to validate the performance of our
approach.

###  4.1 Performance Metrics

Precision, Recall, and F1-Score have been calculated for each of the networks
as:

$\small{Precision=\frac{True~{}Positives}{True~{}Positives+False~{}Positives}}$
(5)
$\small{Recall=\frac{True~{}Positives}{True~{}Positives+False~{}Negatives}}$
(6) $\small{F1~{}Score=\frac{2}{\frac{1}{Precision}+\frac{1}{Recall}}}$ (7)

Table 2: Table showing quantitative evaluation of the proposed network with state-of-the-art architectures Network Architectures | Precision | Recall | F1 Score | # of Parameters
---|---|---|---|---
Two Staged R-CNN [2] | 0.9494 | 0.9494 | 0.9494 | -
ResNet-152 | 0.9906 | 0.8166 | 0.8952 | 41,407,238
ResNet-50 | 0.9909 | 0.8416 | 0.9102 | 21,312,267
ResNet-34 | 0.9892 | 0.8841 | 0.9337 | 21,312,267
Inception-v3 | 0.9846 | 0.7413 | 0.8458 | 21,815,078
Inception-ResNet-v2 | 0.9955 | 0.8552 | 0.9200 | 54,345,958
Tiny-Inception-ResNet-v2 (Proposed) | 0.9854 | 0.9052 | 0.9435 | 19,682,011

###  4.2 Quantitative Results

The proposed method explained in the previous section to identify brick kilns
has been tried for various standard deep learning residual architectures.
Performance metrics with their number of learnable parameters have been shown
in Table. 2. It is evidently clear from the table that even though deep
ResNets and Inception networks may consists of many layers and much more
learnable parameters, still they fail to perform for our problem. Our proposed
Tiny-Inception-ResNet-v2 may lag behind in terms of precision value; however,
it leads the list in recall, F1-score, and number of learnable parameters. One
of the major problems with the two staged R-CNN approach [2] was, huge number
of false positives identified by the binary classifier. The issue was resolved
by adapting the two-staged classification approach. Our proposed multi-class
classifier has very less false positive results with a recall of about
$90.52\%$ and F1-Score of $94.35\%$ with 19m parameters only.

###  4.3 Qualitative Results

For the qualitative analysis of our approach we showed the results of Kot
Radha Kishan, Punjab, Pakistan in Figures 5 and Punjab, India in 6
respectively. The images clearly show that brick kilns with various structural
differences have been identified by our proposed classifier. Moreover, our
classifier efficiently recognizes the brick kiln in the image despite the
variation in size, spatial location and orientation of the kiln in the image.
The heatmap generated pin points each and every brick kiln identified in block
at zoom 17 adapting the proposed method.

##  5 Conclusion

We proposed a novel network architecture: “Tiny-Inception-ResNet-v2” to
identify the brick kilns in satellite images. The network was trained on zoom-
level $20$ images ($782$ images) of brick kilns of Lahore, Pakistan while
tested on the zoom-level $17$ images of South Asia (about $700$). Despite the
structural and color variation in the images of kilns, our proposed network,
with very fewer learning parameters, outperforms all the state-of-the-art
network architectures achieving F1-Score of around $94\%$. We also provided
detailed geo-referenced dataset and annotations for 11 classes, which will
serve as a valuable resource for further such analysis.

## References

  * [1] D. S. Boyd, B. Jackson, J. Wardlaw, G. M. Foody, S. Marsh, and K. Bales, “Slavery from space: Demonstrating the role for satellite remote sensing to inform evidence-based action related to un sdg number 8,” _ISPRS journal of photogrammetry and remote sensing_ , vol. 142, pp. 380–388, 2018.
  * [2] G. M. Foody, F. Ling, D. S. Boyd, X. Li, and J. Wardlaw, “Earth observation and machine learning to meet sustainable development goal 8.7: Mapping sites associated with slavery from space,” _Remote Sensing_ , vol. 11, no. 3, p. 266, 2019.
  * [3] M. Xie, N. Jean, M. Burke, D. Lobell, and S. Ermon, “Transfer learning from deep features for remote sensing and poverty mapping,” in _Thirtieth AAAI Conference on Artificial Intelligence_ , 2016.
  * [4] B. Jackson, K. Bales, S. Owen, J. Wardlaw, and D. S. Boyd, “Analysing slavery through satellite technology: How remote sensing could revolutionise data collection to help end modern slavery,” _J. Mod. Slavery_ , vol. 4, pp. 169–199, 2018.
  * [5] D. Blumberg and D. Jacobson, “New frontiers: remote sensing in social science research,” _The American Sociologist_ , vol. 28, no. 3, pp. 62–68, 1997\.
  * [6] Y. Huo, Z. Xu, S. Bao, A. Assad, R. G. Abramson, and B. A. Landman, “Adversarial synthesis learning enables segmentation without target modality ground truth,” in _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_. IEEE, 2018, pp. 1217–1220.
  * [7] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, p. 436, 2015.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [9] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask r-cnn,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 2961–2969.
  * [10] N. Jean, M. Burke, M. Xie, W. M. Davis, D. B. Lobell, and S. Ermon, “Combining satellite imagery and machine learning to predict poverty,” _Science_ , vol. 353, no. 6301, pp. 790–794, 2016.
  * [11] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. Alemi, “Inception-v4, Inception-ResNet and the impact of residual connections on learning,” in _Association for the Advancement of Artificial Intelligence_ , vol. 4, 2017, p. 12.