Current address: ]Departamento de Física, Universidade Estadual de Maringá,
87020-900 Maringá, PR, Brazil

# Magnetic field-induced weak-to-strong-link crossover in patterned
superconducting films

Davi A. D. Chaves<EMAIL_ADDRESS>Departamento de Física, Universidade
Federal de São Carlos, 13565-905 São Carlos, SP, Brazil  M. I. Valerio-Cuadros
[  Departamento de Física, Universidade Federal de São Carlos, 13565-905 São
Carlos, SP, Brazil  L. Jiang  School of Aeronautics, Northwestern
Polytechnical University, Xi’an 710072, China  Experimental Physics of
Nanostructured Materials, Q-MAT, CESAM, Université de Liège, B-4000 Sart
Tilman, Belgium  E. A. Abbey  Departamento de Física, Universidade Federal de
São Carlos, 13565-905 São Carlos, SP, Brazil  F. Colauto  Departamento de
Física, Universidade Federal de São Carlos, 13565-905 São Carlos, SP, Brazil
A. A. M. Oliveira  Instituto Federal de Educação, Ciência e Tecnologia de São
Paulo, Campus São Carlos, 13565-905, São Carlos, SP, Brazil  A. M. H. Andrade
Instituto de Física, Universidade Federal do Rio Grande do Sul, 91501-970
Porto Alegre, RS, Brazil.  L. B. L. G. Pinheiro  Departamento de Física,
Universidade Federal de São Carlos, 13565-905 São Carlos, SP, Brazil
Instituto Federal de Educação, Ciência e Tecnologia de São Paulo, Campus São
Carlos, 13565-905, São Carlos, SP, Brazil  T. H. Johansen  Department of
Physics, University of Oslo, P.O. Box 1048 Blindern, 0316 Oslo, Norway  C. Xue
School of Mechanics, Civil Engineering and Architecture, Northwestern
Polytechnical University, Xi’an 710072, China  Y.-H. Zhou  Department of
Mechanics and Engineering Sciences, Lanzhou University, Key Laboratory of
Mechanics on Disaster and Environment in Western China, Ministry of Education
of China, Lanzhou 730000, China  School of Aeronautics, Northwestern
Polytechnical University, Xi’an 710072, China  A. V. Silhanek  Experimental
Physics of Nanostructured Materials, Q-MAT, CESAM, Université de Liège, B-4000
Sart Tilman, Belgium  W. A. Ortiz  Departamento de Física, Universidade
Federal de São Carlos, 13565-905 São Carlos, SP, Brazil  M. Motta
<EMAIL_ADDRESS>Departamento de Física, Universidade Federal de São
Carlos, 13565-905 São Carlos, SP, Brazil

###### Abstract

Ubiquitous in most superconducting materials, weak-links are well-known for
their limiting effect on the critical current. Still, they are at the root of
important features of superconductivity. Using quantitative magneto-optical
imaging and the thermomagnetic model, we correlate local maxima in
magnetization loops of patterned Nb films to a magnetic field-induced weak-to-
strong-link crossover. This tunable phenomenon can be exploited in
superconducting electronics, e.g., flux valves and injectors, or as a simple
strategy to adjust the critical current in specialized devices.

Nanoscale patterning of superconducting thin films enables the optimization
and control of distinct material properties [1, 2, 3, 4, 5, 6, 7, 8, 9].
Moreover, patterning allows exploring different phenomena rising from the rich
physics of superconducting weak-links (WLs) [10, 11, 12, 13]. Today, patterned
superconducting structures find applications in memories [14, 15, 16, 17],
diodes [18, 19], transistors [20, 21, 22], quantum batteries [23], and phase
shifters [24, 25], serving as the backbone of superconducting electronics.

Meanwhile, focused-ion beam (FIB) milling is a prominent nanofabrication
technique with maximum spatial resolution in the nanometer range [26, 27, 28].
However, FIB unavoidably introduces defects along the patterned regions,
arising from ionic implantation [29, 30, 31], and locally increasing pinning
in superconducting films [3, 8, 32]. Recently, we investigated the effects of
a single FIB-milled WL on the properties of prototypical low-critical
temperature Nb films, revealing the existence of unexpected local maxima in
their magnetization loops [33].

Weak-links are also an important challenge in the optimization of large-scale
applications of high-temperature superconductors (HTSC) [34, 35].
Particularly, the overall critical current density ($J_{c}$) of HTSCs is
affected by an angular-dependent deterioration of $J_{c}$ along grain
boundaries (GB) [36, 37]. This obstacle led to the development of specialized
fabrication techniques to minimize the GB impact on $J_{c}$, such as IBAD and
RABiTS [38, 39, 40, 41]. Nevertheless, the improvement of transport carrying
capacities, both in HTSCs and low-critical temperature materials, remains an
active research topic [42, 43, 44, 45, 46, 47].

In this study, we experimentally and numerically investigate the magnetic flux
penetration and shielding currents in FIB-patterned Nb films containing a
single macroscopic WL. Our results shine light on the origin of unexpected
local maxima in the magnetization loops, revealing a field-induced crossover
from a weak-link to a strong-link behavior by the enhancement of transport
capacities across the WL.

Three 180-nm-thick Nb films were patterned with a single shallow groove across
their center using FIB doses of 0.1, 0.2, and 0.3 nC/µm2. The grooved films
are labeled GF01, GF02, and GF03, accordingly. GF01 has a surface area of 3
$\times$ 1 mm2, the others 2.5 $\times$ 0.8 mm2. This process creates WLs
comprised of thinner Nb regions pervaded by defects arising from Ga+
implantation, as shown in the Supplemental Material (SM)
[[SeeSupplementalMaterialat][formoreinformationonthemorphology, structure,
andsuperconductingaspectsofGF01.]supp]. The groove depths range from 4.2 nm to
14.6 nm. A pristine 180-nm Nb film with critical temperature ($T_{c}$) of 9.0
K is also studied. Details of the sample preparation are given in Ref. [33].
The global DC magnetic response of the samples is obtained in a SQUID-based
Quantum Design MPMS-5S system in the perpendicular geometry. Quantitative
magneto-optical imaging (MOI) is used to reveal the local out-of-plane
magnetic flux density ($B$) distribution in the films [49]. The local
shielding current density ($\bm{J}$) is obtained as described in Ref. [50].

Figure 1: Representation and MOI at 7 K of fully flux-penetrated (a) pristine
and (b) grooved GF01 films under a perpendicular $H$ after ZFC. Red arrows
represent the direction of $\bm{J}$, $\alpha$ is the angle between the d-line
and the edge. (c) The positive field-decreasing branch of $m(H)$ for GF01 and
GF03. Inset shows $m(H)$ for pristine and GF02 samples. Red arrows indicate
the local maxima in $m(H)$. (d) MOI of GF01 at 7 K for values of $H$ indicated
by vertical arrows in the main panel in (c). Colored rectangles represent
areas used in the quantitative study in Fig. 2(b-c).

The left panel of Fig. 1(a) schematically represents a superconducting film
fully penetrated by magnetic flux under a perpendicular applied magnetic field
$H$ after zero-field cooling (ZFC) [51, 52]. We observe five dark
discontinuity lines, or d-lines, where $\bm{J}$ changes its direction,
shielding flux more efficiently. They define two trapezoidal and two
triangular current domains in the sample. The magneto-optical (MO) image at
11.0 Oe demonstrates that flux penetrates uniformly from the edges toward the
center of the pristine film as ${H}$ is increased. This is manifested by
bright, flux-filled regions surrounding the dark shielded central portion of
the sample. At 62.0 Oe, we observe the four expected domains for the fully
penetrated film.

In contrast, a groove across the shortest symmetry line of the rectangular
film creates a WL and defines two apparently disconnected pristine regions
(PR). Since the WL has a lower $J_{c}$ than the PR, $\bm{J}$ needs to bend
away from the groove, resulting in the central diamond-shaped domain
represented in Fig. 1(b) [53, 54, 55]. In the limiting case when the depicted
angle $\alpha$ = 45∘, no current is able to flow through the WL, i.e.,
$J_{c}^{\text{WL}}=0$. MOI for GF01 at 13.9 Oe reveals that, although the flux
front has reached only a small portion of the PR, the WL is fully penetrated—a
consequence of its weaker magnetic shielding capacity. At 145.4 Oe, the
triangular and diamond-shaped domains can be clearly distinguished. Contrary
to the schematic representation, MOI for GF01 reveals $\alpha>45^{\circ}$ and
a faded dark vertical d-line inside the diamond-shaped domain. These effects
happen because $J_{c}^{\text{WL}}>0$, meaning that a fraction of $\bm{J}$ is
able to flow through the WL, as we will discuss ahead. The dark scratches
above the diamond-shaped and bottom triangular domains are defects on the MO
indicator and do not interfere with the flux penetration into the sample or
the presented data analysis.

The inset of Fig. 1(c) shows complete $m(H)$ hysteresis loops at 7 K for the
pristine film and GF02, where $m$ is the magnetic moment. For the pristine
film, the observed behavior matches that expected for a type-II superconductor
presenting a flux-dependent critical current density, $J_{c}(B)$ [56, 57]. In
contrast, the nanostructured sample GF02 exhibits unexpected local maxima in
the positive and negative decreasing-field branches of the $m(H)$ loop. As
shown in the main panel of Fig. 1(c), the same feature is observed around 50
Oe for GF01 and GF03. Notably, $J_{c}$ is proportional to the height of the
$m(H)$ curve [58].

We elect sample GF01 to illustrate our investigation on what leads to this
phenomenon. GF01 contains a 4.2-nm-deep groove with a 800 nm width at half
depth—see SM. Figure 1(d) shows a series of MO images of GF01 at 7 K as ${H}$
is decreased from a state of full penetration. For 140.8 Oe, we observe flux
penetration similar to that of Fig 1(b). The WL now appears in dark contrast,
because $H$ was reduced from 315.6 Oe, and the WL shows a lower, but positive
flux density due to its weaker shielding capability—see SM. However, as ${H}$
approaches 50 Oe, the images reveal that the d-lines forming the central
diamond-shaped domain move toward the groove—as if they were closing. The
image at 50.2 Oe shows great resemblance to the case of a pristine film, as
the diamond-shape practically vanishes. If $H$ is reduced beyond the local
maximum of $m(H)$, the d-lines now move away from the groove, as if they were
reopening, reestablishing the diamond-shaped domain, as seen at 25.6 Oe. This
analysis demonstrates that MOI is an ideal tool for our study, as it reveals
details of the local flux penetration and current distribution that are hidden
in global magnetic measurements, such as DC magnetization.

The observed d-line movement is associated to an enhancement of the transport
properties in the WL. This is demonstrated by the MO results in Fig. 2(a). As
the d-lines close, $\alpha$ increases and $\bm{J}$ changes
accordingly—streamlines that initially did not cross the WL become straighter
and have a higher density across the WL. This trend is maximized at 48.0 Oe,
when the diamond-shaped domain almost fades out. This result shows that, after
the sample is fully penetrated and $H$ is reduced to a specific value, the WL
behaves as a strong-link, as its transport properties are greatly enhanced and
current may flow through it largely unaffected. At 46.0 Oe and as $H$ is
further reduced, the d-lines reappear, reestablishing the WL behavior. This
notion explains the local peak in $m(H)$ as an increase in the overall $J_{c}$
due to an increase in $J_{c}^{\text{WL}}$ as the sample approaches the
pristine behavior.

Figure 2: (a) MOI of GF01 at 7 K as $H$ is reduced from full penetration. The
$\bm{J}$ distribution is shown as red streamlines. White dashed lines define
the angle $\alpha$ shown in the panels. (b) $m(H)$ for GF01 evaluated by MOI.
(c) $m_{l}(H)$ averaged over different regions of the sample. (d) Upper panel
shows $\Phi_{\text{out}}$ besides the WL and the PR. Lower panel highlights
the difference in $\Phi_{\text{out}}$ between the two regions. The presented
quantities are obtained for the color-matched regions shown in Fig. 1(d).

In Fig. 2(b), we evaluate the magnetic moment of GF01 via MOI as
$m=\sum_{y}m_{l}(y)l_{px}=l_{px}\sum_{y}2\int_{0}^{w/2}xj_{y}(x)dx$, where
$m_{l}(y)$ is $m$ per unit length at a position $y$, $l_{px}$ is the image
pixel size, $w$ is the film width, $j_{y}=J_{y}d$ is the sheet current flowing
perpendicular to the WL, and $d$ is the film thickness [59, 60]. The results
bare great resemblance to those in Fig. 1(c). Figure 2 also presents
quantitative results obtained along the four color-matched regions in Fig.
1(d). In Fig. 2(c), we average $m_{l}$ over the WL (blue) and on a selected
portion of the PR (red), confirming that the peak in $m(H)$ manifests an
effect tied to the WL—as a peak in $m_{l}$ appears for the WL, but not for the
PR. Moreover, we are able to gauge the flux behavior just outside the sample,
$\Phi_{\text{out}}$, by the numerical integration of $B$. Figure 2(d) shows
that, as $H$ is reduced, $\Phi_{\text{out}}$ is everywhere higher just besides
the WL (pink) than besides the PR (yellow) for equal areas. However, as the
d-lines are closing, $\Phi_{\text{out}}$ decreases close to the WL while it is
unaffected besides the PR, lessening the difference between
$\Phi_{\text{out}}$ along the two regions (green points). This indicates that
less flux is now able to escape the WL. When the d-lines reopen below 50 Oe,
$\Phi_{\text{out}}$ is significantly enhanced besides the groove, indicating a
large amount of flux is expelled through the edges of the WL.

Figure 3 locally resolves the behavior of $B$ and $J$ for GF01. Figure 3(a)
depicts the differential flux density, $B_{\text{dif}}$. It captures the $B$
variation due to the variation of $H$ by subtracting the target MO image $n$
by the one taken at the previous field step, i.e.
$B_{\text{dif}}~{}=~{}B_{n}-B_{n-1}$. The field step is kept around $-2$ Oe in
the $H$ range depicted in Fig. 3. In the first three $B_{\text{dif}}$ images,
we visualize the closing of the d-lines by the shrinking of a brighter
diamond-shaped inner region. This indicates that the $B$ decrease in that
region is lower than that outside the d-lines. Therefore, flux is pushed
toward the WL as the d-lines close. This agrees with the fact that flux is
unable to cross d-lines, as these are regions of high shielding capacity [61].
At 46.0 Oe, a dark-colored region anchored in the central part of the groove
appears in $B_{\text{dif}}$, revealing flux is intensely expelled from the WL,
as also seen in Fig. 2(d). At 43.8 Oe, flux is still being strongly pushed
away from the inner part of the sample and is directed toward the edges and
the PR—as confirmed by the bright halo centered in the groove, signaling a
lower variation in $B_{\text{dif}}$ although $H$ decreases uniformly [62].
This demonstrates that the reappearance of the d-lines begins from the central
part of the film, not from the edges. For 39.9 Oe and further, the flux pushed
out of the WL is systematically reduced.

Figure 3: Magnetic flux and shielding currents for GF01 at 7 K at different
$H$. (a) $B_{\text{dif}}$ distribution around the region of d-line movement.
Thin vertical lines represent the edges of the film. (b) $B$ and (c) $J_{y}$
profiles along the WL and PR, evaluated respectively at the black and red
dashed lines in panel (a).

Figure 3(b) shows the evolution of $B$ along the WL (in black) and the PR
(red). From 92.6 to 48.0 Oe, the black profile decreases less and becomes
increasingly similar to the red one, indicating that the flux pushed inward by
the d-line movement is partially retained by the WL. This is related to an
increased pinning in the WL due to its reduced thickness and FIB-induced
defects, which act as pinning centers that also locally suppress $T_{c}$
[33]—see SM. The relative increase in $B$ corroborates the increase in $M$
across the WL seen in Fig. 2(b). Starting at 46.0 Oe, trapped flux is
vigorously pushed away from the WL and the $B$ profiles become distinct again.
Interestingly, at 46.0 Oe $B$ is positive along the WL. These observations
suggest that the maximum in $m(H)$ is not related to a net neutral flux in the
WL. Here, we argue that the increased flux pinning in the WL reduces energy
dissipation coming from flux line movement, increasing $J_{c}^{\text{WL}}$
and, therefore, acting to equalize $J_{c}(B)$ for the PR and WL for a specific
$H$ value. A similar behavior exists in HTSCs [63, 64, 60, 65, 66, 67, 68].
The fact that the associated d-line movement was not observed in previously
MOI-studied YBCO bicrystalline thin films [53, 69, 70, 68] may be related to
the absence of the FIB-added pinning centers and the relatively higher $H$ for
which the peak in $m(H)$ occurs in HTSCs, rendering $J_{c}(B)$-dependent
effects less prominent [57].

Figure 3(c) presents an analysis of the current density flowing
perpendicularly to the WL, $J_{y}$. Initially, there is a gradual increase in
the current able to cross the WL, highlighted by an enhancement of the black
curve magnitude. At 48.0 Oe, $J_{y}$ across the WL bares a remarkable
resemblance to that flowing in the PR—as the sample behaves almost as a
pristine film due to a weak-to-strong-link behavior crossover. The trend is
reversed when the d-lines reappear and a sharp decrease in $J_{y}$ is observed
beginning in the center of the WL. This occurs due to pinning centers created
by the FIB. The $B$ profile along the WL at 48.0 Oe displays the highest flux
density at the center of the film. Therefore, as the density of vortices
increases, the electromagnetic repulsion between them eventually overcomes the
pinning force acting to keep the flux contained in the WL, thus expelling the
flux, inducing a strong-link to weak-link crossover. The SM depicts the
behavior of $B$ and $J_{y}$ for a large number of images at different applied
fields, corroborating their relationship to the d-line movement.

With the experimental input, we now turn to numerical simulations based on the
thermomagnetic model (TM), which allow for the comprehension of subtle
distinctions in specimen behavior arising from a $B$-dependent $J_{c}$ [71,
72, 73]—see SM for further information. The simulated Nb sample shares the
geometry of GF01, except for the groove width, set to 45 µm to improve
computational performance. Notably, the higher density of pinning sites
introduced by Ga implantation in the WL is not included in the simulations. A
reduction of $\sim$15% in the zero-field $J_{c}^{\text{WL}}$ is induced by
reducing $T_{c}$ in the WL, as experimentally observed—see SM and Ref. [33].

Figure 4: (a) $B$ distribution at 7 K captured by TM simulations as $H$ is
reduced from 250 Oe. Simulation parameters follow Ref. [73]. (b)
$\Delta\alpha(H)$ for the simulated sample. Red arrows indicate the direction
of field variation. Inset: experimental $\Delta\alpha(H)$ for GF01. Error lies
within symbol size.

Figure 4 describes the TM results obtained when considering
$j_{c}=j_{c}(T,B)=j_{c0}(1-T/T_{c})(B_{0}/(B+B_{0}))$, with
$B_{0}=\mu_{0}j_{c0}/\pi$ [74]. Disorder is introduced by lowering $j_{c0}$
for randomly selected simulation grid points [75]. First, the magnetic field
is increased up to 250 Oe after ZFC to 7 K. The resulting flux distribution in
Fig. 4(a) qualitatively matches that of GF01 in the full penetration state.
Then, $H$ is progressively reduced. At 50 Oe, we visually notice the
disappearance of the diamond-shaped domain, which is maximum at 20 Oe. This is
quantitatively corroborated by Fig. 4(b), which depicts
$\Delta\alpha(H)=\alpha(H=250~{}\text{Oe})-\alpha(H)$. The inset of Fig. 4(b)
confirms that the TM simulations reproduce the d-line movement observed
experimentally for GF01. The lesser closure of the d-lines is likely due to
the larger width of the simulated groove in comparison to GF01. Finally, the
d-lines reappear, as demonstrated in Fig. 4(a) at $-50$ Oe. Additional
simulations, conducted for $J_{c}$ independent of $B$, i.e.,
$J_{c}=J_{c0}(1-T/T_{c})$, did not reproduce the d-line movement—reinforcing
the $J_{c}(B)$-dependency role in the observed phenomenon.

In summary, combining quantitative magneto-optical imaging and thermomagnetic
model simulations, we demonstrate that the presence of local maxima in the
magnetization loop of FIB-patterned Nb films is associated with a field-
induced weak-link to strong-link to weak-link crossover. Such an effect may be
explored in the development of novel sensing applications, as well as for
different implementations in superconducting electronics. For instance, the
tunable flux escape observed just outside the groove may be employed in a flux
valve or injector. Mastering this phenomenon may also lead to the improvement
of operating conditions for different superconducting technologies, resulting
in effectively higher critical currents for films in which the presence of
weak-links is essential or unavoidable. We additionally envision that, if the
proper $J_{c}(B)$-dependency, pinning center distribution, and field
conditions are met, the reported mechanism could play an important role in
HTSCs—solidifying the understanding of a post-fabrication route for $J_{c}$
optimization in these materials.

###### Acknowledgements.

This work was partially supported by Coordenação de Aperfeiçoamento de Pessoal
de Nível Superior – Brasil (CAPES) – Finance Code 001, the São Paulo Research
Foundation (FAPESP, Grant No. 2021/08781-8), and the National Council for
Scientific and Technological Development (CNPq, Grants No. 431974/2018-7, No.
316602/2021-3, and No. 309928/2018-4). C. X. acknowledges the support by the
National Natural Science Foundation of China (Grants No. 11972298). L. J. was
supported by the China Scholarship Council. The authors thank the Laboratory
of Structural Characterization (LCE/DEMa/UFSCar) for the EDS and AFM
measurements reported in the SM.

## References

  * Harada  _et al._ [1996] K. Harada, O. Kamimura, H. Kasai, T. Matsuda, A. Tonomura, and V. V. Moshchalkov, Direct observation of vortex dynamics in superconducting films with regular arrays of defects, Science 274, 1167 (1996).
  * Silhanek  _et al._ [2003] A. V. Silhanek, L. Van Look, S. Raedts, R. Jonckheere, and V. V. Moshchalkov, Guided vortex motion in superconductors with a square antidot array, Phys. Rev. B 68, 214504 (2003).
  * Dobrovolskiy  _et al._ [2012] O. V. Dobrovolskiy, E. Begun, M. Huth, and V. A. Shklovskij, Electrical transport and pinning properties of Nb thin films patterned with focused ion beam-milled washboard nanostructures, New J. Phys. 14, 113027 (2012).
  * Shaw  _et al._ [2012] G. Shaw, P. Mandal, B. Bag, S. Banerjee, T. Tamegai, and H. Suderow, Properties of nanopatterned pins generated in a superconductor with FIB, Appl. Surf. Sci. 258, 4199 (2012).
  * Wang  _et al._ [2013] Y. L. Wang, M. L. Latimer, Z. L. Xiao, R. Divan, L. E. Ocola, G. W. Crabtree, and W. K. Kwok, Enhancing the critical current of a superconducting film in a wide range of magnetic fields with a conformal array of nanoscale holes, Phys. Rev. B 87, 220501 (2013).
  * Dobrovolskiy [2017] O. V. Dobrovolskiy, Abrikosov fluxonics in washboard nanolandscapes, Physica C 533, 80 (2017).
  * Kalcheim  _et al._ [2017] Y. Kalcheim, E. Katzir, F. Zeides, N. Katz, Y. Paltiel, and O. Millo, Dynamic control of the vortex pinning potential in a superconductor using current injection through nanoscale patterns, Nano Lett. 17, 2934 (2017).
  * Dobrovolskiy  _et al._ [2017] O. V. Dobrovolskiy, V. A. Shklovskij, M. Hanefeld, M. Zörb, L. Köhs, and M. Huth, Pinning effects on flux flow instability in epitaxial Nb thin films, Supercond. Sci. Technol. 30, 085002 (2017).
  * Colauto  _et al._ [2020] F. Colauto, M. Motta, and W. A. Ortiz, Controlling magnetic flux penetration in low-$T_{c}$ superconducting films and hybrids, Supercond. Sci. Technol. 34, 013002 (2020).
  * Likharev [1979] K. K. Likharev, Superconducting weak links, Rev. Mod. Phys. 51, 101 (1979).
  * Moshchalkov  _et al._ [2010] V. Moshchalkov, R. Woerdenweber, and W. Lang, eds.,  _Nanoscience and engineering in superconductivity_ (Springer Berlin, Heidelberg, 2010).
  * Fornieri  _et al._ [2019] A. Fornieri, A. M. Whiticar, F. Setiawan, E. Portolés, A. C. Drachmann, A. Keselman, S. Gronin, C. Thomas, T. Wang, R. Kallaher, _et al._ , Evidence of topological superconductivity in planar Josephson junctions, Nature 569, 89 (2019).
  * Holzman and Ivry [2019] I. Holzman and Y. Ivry, Superconducting nanowires for single-photon detection: Progress, challenges, and opportunities, Adv. Quantum Technol. 2, 1800058 (2019).
  * Murphy  _et al._ [2017] A. Murphy, D. V. Averin, and A. Bezryadin, Nanoscale superconducting memory based on the kinetic inductance of asymmetric nanowire loops, New J. Phys. 19, 063015 (2017).
  * Chen  _et al._ [2020] L. Chen _et al._ , Miniaturization of the superconducting memory cell via a three-dimensional Nb nano-superconducting quantum interference device, ACS Nano 14, 11002 (2020).
  * Ilin  _et al._ [2021] E. Ilin, X. Song, I. Burkova, A. Silge, Z. Guo, K. Ilin, and A. Bezryadin, Supercurrent-controlled kinetic inductance superconducting memory element, Appl. Phys. Lett. 118, 112603 (2021).
  * Chaves  _et al._ [2023] D. A. D. Chaves, L. Nulens, H. Dausy, B. Raes, D. Yue, W. A. Ortiz, M. Motta, M. J. Van Bael, and J. Van de Vondel, Nanobridge SQUIDs as multilevel memory elements, Phys. Rev. Appl. 19, 034091 (2023).
  * Lyu  _et al._ [2021] Y.-Y. Lyu _et al._ , Superconducting diode effect via conformal-mapped nanoholes, Nat. Commun. 12, 1 (2021).
  * Golod and Krasnov [2022] T. Golod and V. M. Krasnov, Demonstration of a superconducting diode-with-memory, operational at zero magnetic field with switchable nonreciprocity, Nat. Commun. 13, 3658 (2022).
  * Giazotto  _et al._ [2010] F. Giazotto, J. T. Peltonen, M. Meschke, and J. P. Pekola, Superconducting quantum interference proximity transistor, Nat. Phys. 6, 254 (2010).
  * Paolucci  _et al._ [2018] F. Paolucci, G. De Simoni, E. Strambini, P. Solinas, and F. Giazotto, Ultra-efficient superconducting dayem bridge field-effect transistor, Nano Lett. 18, 4195 (2018).
  * Ligato  _et al._ [2022] N. Ligato, F. Paolucci, E. Strambini, and F. Giazotto, Thermal superconducting quantum interference proximity transistor, Nat. Phys. , 1 (2022).
  * Strambini  _et al._ [2020] E. Strambini _et al._ , A Josephson phase battery, Nat. Nanotechnol. 15, 656 (2020).
  * Yamashita  _et al._ [2020] T. Yamashita, S. Kim, H. Kato, W. Qiu, K. Semba, A. Fujimaki, and H. Terai, $\pi$ phase shifter based on NbN-based ferromagnetic Josephson junction on a silicon substrate, Sci. Rep. 10, 1 (2020).
  * Golod  _et al._ [2021] T. Golod, R. A. Hovhannisyan, O. M. Kapran, V. V. Dremov, V. S. Stolyarov, and V. M. Krasnov, Reconfigurable Josephson phase shifter, Nano Lett. 21, 5240 (2021).
  * Volkert and Minor [2007] C. A. Volkert and A. M. Minor, Focused ion beam microscopy and micromachining, MRS Bull. 32, 389 (2007).
  * Wyss  _et al._ [2022] M. Wyss, K. Bagani, D. Jetter, E. Marchiori, A. Vervelaki, B. Gross, J. Ridderbos, S. Gliga, C. Schönenberger, and M. Poggio, Magnetic, thermal, and topographic imaging with a nanometer-scale SQUID-on-lever scanning probe, Phys. Rev. Appl. 17, 034002 (2022).
  * Sigloch  _et al._ [2022] F. Sigloch, S. Sangiao, P. Orús, and J. M. de Teresa, Direct-write of tungsten-carbide nanoSQUIDs based on focused ion beam induced deposition, Nanoscale Adv. 4, 4628 (2022).
  * Datesman  _et al._ [2005] A. M. Datesman, J. C. Schultz, T. W. Cecil, C. M. Lyons, and A. W. Lichtenberger, Gallium ion implantation into niobium thin films using a focused-ion beam, IEEE Trans. Appl. Supercond. 15, 3524 (2005).
  * De Leo  _et al._ [2016] N. De Leo, M. Fretto, V. Lacquaniti, C. Cassiago, L. D’Ortenzi, L. Boarino, and S. Maggi, Thickness modulated niobium nanoconstrictions by focused ion beam and anodization, IEEE Trans. Appl. Supercond. 26, 1 (2016).
  * Singh  _et al._ [2018] M. Singh, R. Chaujar, S. Husale, S. Grover, A. P. Shah, M. M. Deshmukh, A. Gupta, V. N. Singh, V. N. Ojha, D. K. Aswal, and R. K. Rakshit, Influence of fabrication processes on transport properties of superconducting niobium nitride nanowires, Current Science 114, 1443 (2018).
  * Aichner  _et al._ [2019] B. Aichner, B. Müller, M. Karrer, V. R. Misko, F. Limberger, K. L. Mletschnig, M. Dosmailov, J. D. Pedarnig, F. Nori, R. Kleiner, D. Koelle, and W. Lang, Ultradense tailored vortex pinning arrays in superconducting YBa2Cu3O7-δ thin films created by focused He ion beam irradiation for fluxonics applications, ACS Appl. Nano Mater. 2, 5108 (2019).
  * Valerio-Cuadros  _et al._ [2021] M. I. Valerio-Cuadros, D. A. D. Chaves, F. Colauto, A. A. M. Oliveira, A. M. H. Andrade, T. H. Johansen, W. A. Ortiz, and M. Motta, Superconducting properties and electron scattering mechanisms in a Nb film with a single weak-link excavated by focused ion beam, Materials 14, 7274 (2021).
  * Larbalestier  _et al._ [2001] D. Larbalestier, A. Gurevich, D. M. Feldmann, and A. Polyanskii, High-$T_{c}$ superconducting materials for electric power applications, Nature 414, 368 (2001).
  * Foltyn  _et al._ [2007] S. R. Foltyn, L. Civale, J. L. MacManus-Driscoll, Q. X. Jia, B. Maiorov, H. Wang, and M. Maley, Materials science challenges for high-temperature superconducting wire, Nat. Mater. 6, 631 (2007).
  * Hilgenkamp and Mannhart [2002] H. Hilgenkamp and J. Mannhart, Grain boundaries in high-${T}_{c}$ superconductors, Rev. Mod. Phys. 74, 485 (2002).
  * Graser  _et al._ [2010] S. Graser, P. J. Hirschfeld, T. Kopp, R. Gutser, B. M. Andersen, and J. Mannhart, How grain boundaries limit supercurrents in high-temperature superconductors, Nat. Phys. 6, 609 (2010).
  * Reade  _et al._ [1992] R. P. Reade, P. Berdahl, R. E. Russo, and S. M. Garrison, Laser deposition of biaxially textured yttria‐stabilized zirconia buffer layers on polycrystalline metallic alloys for high critical current Y‐Ba‐Cu‐O thin films, Appl. Phys. Lett. 61, 2231 (1992).
  * Iijima  _et al._ [1993] Y. Iijima, K. Onabe, N. Futaki, N. Tanabe, N. Sadakata, O. Kohno, and Y. Ikeno, Structural and transport properties of biaxially aligned YBa2Cu3O7-x films on polycrystalline Ni‐based alloy with ion‐beam‐modified buffer layers, J. Appl. Phys. 74, 1905 (1993).
  * Wu  _et al._ [1995] X. D. Wu, S. R. Foltyn, P. N. Arendt, W. R. Blumenthal, I. H. Campbell, J. D. Cotton, J. Y. Coulter, W. L. Hults, M. P. Maley, H. F. Safar, and J. L. Smith, Properties of YBa2Cu3O7-δ thick films on flexible buffered metallic substrates, Appl. Phys. Lett. 67, 2397 (1995).
  * Goyal  _et al._ [1996] A. Goyal, D. P. Norton, J. D. Budai, M. Paranthaman, E. D. Specht, D. M. Kroeger, D. K. Christen, Q. He, B. Saffian, F. A. List, D. F. Lee, P. M. Martin, C. E. Klabunde, E. Hartfield, and V. K. Sikka, High critical current density superconducting tapes by epitaxial deposition of YBa2Cu3Ox thick films on biaxially textured metals, Appl. Phys. Lett. 69, 1795 (1996).
  * Díez-Sierra  _et al._ [2020] J. Díez-Sierra _et al._ , High critical current density and enhanced pinning in superconducting films of YBa2Cu3O7-δ nanocomposites with embedded BaZrO3, BaHfO3, BaTiO3, and SrZrO3 nanocrystals, ACS Appl. Nano Mater. 3, 5542 (2020).
  * Ahmad  _et al._ [2021] I. Ahmad, S. Sarangi, and P. Sarun, Enhanced critical current density and flux pinning of anthracene doped magnesium diboride superconductor, J. Alloys Compd. 884, 160999 (2021).
  * Rocci  _et al._ [2020] M. Rocci, D. Suri, A. Kamra, G. Vilela, Y. Takamura, N. M. Nemes, J. L. Martinez, M. G. Hernandez, and J. S. Moodera, Large enhancement of critical current in superconducting devices by gate voltage, Nano Lett. 21, 216 (2020).
  * Caffer  _et al._ [2021] A. M. Caffer, D. A. D. Chaves, A. L. Pessoa, C. L. Carvalho, W. A. Ortiz, R. Zadorosny, and M. Motta, Optimum heat treatment to enhance the weak-link response of Y123 nanowires prepared by solution blow spinning, Supercond. Sci. Technol. 34, 025009 (2021).
  * Algarni  _et al._ [2021] R. Algarni, M. A. Almessiere, Y. Slimani, E. Hannachi, and F. B. Azzouz, Enhanced critical current density and flux pinning traits with Dy2O3 nanoparticles added to YBa2Cu3O7-d superconductor, J. Alloys Compd. 852, 157019 (2021).
  * Chaves  _et al._ [2021] D. A. D. Chaves, I. M. de Araújo, D. Carmo, F. Colauto, A. A. M. Oliveira, A. M. H. Andrade, T. H. Johansen, A. V. Silhanek, W. A. Ortiz, and M. Motta, Enhancing the effective critical current density in a Nb superconducting thin film by cooling in an inhomogeneous magnetic field, Appl. Phys. Lett. 119, 022602 (2021).
  * [48] URL_will_be_inserted_by_publisher.
  * Shaw  _et al._ [2018] G. Shaw, J. Brisbois, L. B. G. L. Pinheiro, J. Müller, S. Blanco Alvarez, T. Devillers, N. M. Dempsey, J. E. Scheerder, J. Van de Vondel, S. Melinte, P. Vanderbemden, M. Motta, W. A. Ortiz, K. Hasselbach, R. B. G. Kramer, and A. V. Silhanek, Quantitative magneto-optical investigation of superconductor/ferromagnet hybrid structures, Rev. Sci. Instrum. 89, 023705 (2018).
  * Meltzer  _et al._ [2017] A. Y. Meltzer, E. Levin, and E. Zeldov, Direct reconstruction of two-dimensional currents in thin films from magnetic-field measurements, Phys. Rev. Applied 8, 064030 (2017).
  * Schuster  _et al._ [1994a] T. Schuster, M. V. Indenbom, M. R. Koblischka, H. Kuhn, and H. Kronmüller, Observation of current-discontinuity lines in type-II superconductors, Phys. Rev. B 49, 3443 (1994a).
  * Schuster  _et al._ [1995] T. Schuster, H. Kuhn, and M. V. Indenbom, Discontinuity lines in rectangular superconductors with intrinsic and extrinsic anisotropies, Phys. Rev. B 52, 15621 (1995).
  * Polyanskii  _et al._ [1996] A. A. Polyanskii, A. Gurevich, A. E. Pashitski, N. F. Heinig, R. D. Redwing, J. E. Nordman, and D. C. Larbalestier, Magneto-optical study of flux penetration and critical current densities in [001] tilt YBa2Cu3O7-δ thin-film bicrystals, Phys. Rev. B 53, 8687 (1996).
  * Johansen  _et al._ [2019] T. H. Johansen, F. Colauto, A. M. H. de Andrade, A. A. M. Oliveira, and W. A. Ortiz, Transparency of planar interfaces in superconductors: a critical-state analysis, IEEE Trans. Appl. Supercond. 29, 1 (2019).
  * Colauto  _et al._ [2021] F. Colauto, D. Carmo, A. M. H. Andrade, A. A. M. Oliveira, W. A. Ortiz, Y. M. Galperin, and T. H. Johansen, Measurement of critical current flow and connectivity in systems of joined square superconducting plates, Physica C 589, 1353931 (2021).
  * Kim  _et al._ [1963] Y. B. Kim, C. F. Hempstead, and A. R. Strnad, Magnetization and critical supercurrents, Phys. Rev. 129, 528 (1963).
  * Shantsev  _et al._ [2000] D. V. Shantsev, Y. M. Galperin, and T. H. Johansen, Thin superconducting disk with field-dependent critical current: Magnetization and ac susceptibilities, Phys. Rev. B 61, 9699 (2000).
  * Johansen and Bratsberg [1995] T. H. Johansen and H. Bratsberg, Critical‐state magnetization of type‐II superconductors in rectangular slab and cylinder geometries, J. Appl. Phys. 77, 3945 (1995).
  * Brandt and Indenbom [1993] E. H. Brandt and M. Indenbom, Type-II-superconductor strip with current in a perpendicular magnetic field, Phys. Rev. B 48, 12893 (1993).
  * Shantsev  _et al._ [1999] D. V. Shantsev, M. R. Koblischka, Y. M. Galperin, T. H. Johansen, L. Půst, and M. Jirsa, Central peak position in magnetization loops of high- ${T}_{c}$ superconductors, Phys. Rev. Lett. 82, 2947 (1999).
  * Schuster  _et al._ [1994b] T. Schuster, M. V. Indenbom, M. R. Koblischka, H. Kuhn, and H. Kronmüller, Observation of current-discontinuity lines in type-II superconductors, Phys. Rev. B 49, 3443 (1994b).
  * Pinheiro  _et al._ [2020] L. B. L. G. Pinheiro, M. Caputo, C. Cirillo, C. Attanasio, T. H. Johansen, W. A. Ortiz, A. V. Silhanek, and M. Motta, Magnetic flux avalanches in Nb/NbN thin films, Low Temp. Phys. 46, 365 (2020).
  * Evetts and Glowacki [1988] J. E. Evetts and B. A. Glowacki, Relation of critical current irreversibility to trapped flux and microstructure in polycrystalline YBa2Cu3O7, Cryogenics 28, 641 (1988).
  * Dimos  _et al._ [1990] D. Dimos, P. Chaudhari, and J. Mannhart, Superconducting transport properties of grain boundaries in YBa2Cu3O7 bicrystals, Phys. Rev. B 41, 4038 (1990).
  * Hogg  _et al._ [2001] M. J. Hogg, F. Kahlmann, Z. H. Barber, and J. E. Evetts, Angular hysteresis in the critical current of YBa2Cu3O7 low-angle grain boundaries, Supercond. Sci. Technol. 14, 647 (2001).
  * Thompson  _et al._ [2004] J. R. Thompson, H. J. Kim, C. Cantoni, D. K. Christen, R. Feenstra, and D. T. Verebelyi, Self-organized current transport through low-angle grain boundaries in YBa2Cu3O7-δ thin films studied magnetometrically, Phys. Rev. B 69, 104509 (2004).
  * Palau  _et al._ [2004] A. Palau, T. Puig, X. Obradors, E. Pardo, C. Navau, A. Sanchez, A. Usoskin, H. C. Freyhardt, L. Fernández, B. Holzapfel, and R. Feenstra, Simultaneous inductive determination of grain and intergrain critical current densities of YBa2Cu3O7-x coated conductors, Appl. Phys. Lett. 84, 230 (2004).
  * Palau  _et al._ [2007] A. Palau, T. Puig, X. Obradors, and C. Jooss, Simultaneous determination of grain and grain-boundary critical currents in YBa2Cu3O7-coated conductors by magnetic measurements, Phys. Rev. B 75, 054517 (2007).
  * Guth  _et al._ [2004] K. Guth, V. Born, and C. Jooss, Inhomogeneous current distribution in wide high-temperature superconducting small-angle grain boundaries, Eur. Phys. J. B 42, 239 (2004).
  * Born  _et al._ [2004] V. Born, K. Guth, H. C. Freyhardt, and C. Jooss, Self-enhanced flux penetration into small angle grain boundaries in YBCO thin films, Supercond. Sci. Technol. 17, 380 (2004).
  * Jing  _et al._ [2016] Z. Jing, H. Yong, and Y. Zhou, Influences of non-uniformities and anisotropies on the flux avalanche behaviors of type-II superconducting films, Supercond. Sci. Technol. 29, 105001 (2016).
  * Jing  _et al._ [2017] Z. Jing, H. Yong, and Y. Zhou, Numerical simulation on the flux avalanche behaviors of microstructured superconducting thin films, J. Appl. Phys. 121, 023902 (2017).
  * Jiang  _et al._ [2020] L. Jiang, C. Xue, L. Burger, B. Vanderheyden, A. V. Silhanek, and Y.-H. Zhou, Selective triggering of magnetic flux avalanches by an edge indentation, Phys. Rev. B 101, 224505 (2020).
  * Kim  _et al._ [1962] Y. B. Kim, C. F. Hempstead, and A. R. Strnad, Critical persistent currents in hard superconductors, Phys. Rev. Lett. 9, 306 (1962).
  * Vestgården  _et al._ [2011] J. I. Vestgården, D. V. Shantsev, Y. M. Galperin, and T. H. Johansen, Dynamics and morphology of dendritic flux avalanches in superconducting films, Phys. Rev. B 84, 054537 (2011).


]Complexity Science Hub Vienna, Vienna, 1080, Austria

# Maximum Entropy Approach for the Prediction of Urban Mobility Patterns

Simone Daniotti  [<EMAIL_ADDRESS>Bernardo Monechi  Enrico Ubaldi

###### Abstract

The science of cities is a relatively new and interdisciplinary topic. It
borrows techniques from agent-based modeling, stochastic processes, and
partial differential equations. However, how the cities rise and fall, how
they evolve, and the mechanisms responsible for these phenomena are still open
questions. Scientists have only recently started to develop forecasting tools,
despite their importance in urban planning, transportation planning, and
epidemic spreading modeling. Here, we build a fully interpretable statistical
model that, incorporating only the minimum number of constraints, can predict
different phenomena arising in the city. Using data on the movements of car-
sharing vehicles in different Italian cities, we infer a model using the
Maximum Entropy (MaxEnt) principle. With it, we describe the activity in
different city zones and apply it to activity forecasting and anomaly
detection (e.g., strikes, and bad weather conditions). We compare our method
with different models explicitly made for forecasting: SARIMA models and Deep
Learning Models. We find that MaxEnt models are highly predictive,
outperforming SARIMAs and having similar results as a Neural Network. These
results show how relevant statistical inference can be in building a robust
and general model describing urban systems phenomena. This article identifies
the significant observables for processes happening in the city, with the
perspective of a deeper understanding of the fundamental forces driving its
dynamics.

††preprint: APS/123-QED

## Introduction

In recent years, pressing societal and environmental problems such as
population growth, migrations, and climate change, boosted the research on the
Science of Cities and the related study of mobility. The availability of large
and detailed datasets covering the mobility of individuals at different
granularity contributed largely to enhance the interest of researchers in this
field[1, 2]. Being multi-disciplinary, Science of Cities studies embrace
diverse areas of research. For example, statistical methods have been applied
to city growth [3], multi-layer networks to urban resilience [4] and spatial
networks to describe and characterize the structure and the evolution of
phenomena arising on it [5]. On the modeling side, co-evolution models [6] and
agent-based simulations [7] have been used to model stylized facts and take
into account policy making. Moreover, other frameworks such as ranking
dynamics [8] have been used to study urban environments and universal laws
arising in them.

The mobility patterns of individuals diffusing in an urban environment are
determined by the interplay of different mechanisms, such as the daily
routines of the individuals and environmental constraints, both again
regulated by even more fundamental and interrelated phenomena, like wealth,
trends, economic relations, socio-political disparities, and cultural
movements [9, 10, 11]. Urban environments are complex systems [12], and
describing their growth and relation with the surrounding cities is not
unanimously understood by the scientific community [13]. Being that the
commuting phenomenon inside and between cities may be driven by different
events and may be related to different causes, different tools and
understandings must be developed [14, 15].

In this work, we study urban mobility patterns, building a model based on only
a few constraints driven by data observations. The model can predict different
events in the urban environment with high precision and can be generalized to
other dynamical processes unfolding in urban spaces. Here, we propose a
Statistical Inference (Maximum Entropy) approach to study and predict urban
mobility patterns. Phenomena and relative observables that happen inside the
city are complex, they entangle with various indicators and are difficult to
predict. To build a powerful, general, and robust statistical model, in
principle we need to understand what are the important variables that play a
central role in the phenomenon we want to model. To solve it, we need to
analyze the dataset we want to study, identify the important dynamical
properties and then model it solving the problem of optimizing the resulting
entropy.

The data represents the $30-minute-binned$ multi-variate time-series of the
activity of different zones inside the city.

Being that urban systems are notoriously complex and that the fundamental
causes of the observed mobility patterns are various and interrelated, our
methodology is novel in this field in the sense that builds the most general
model constrained to reproduce the correlations observed.

In section Results we present the results for the Metropolitan City of Milan:
multi-variate forecasting of the zones activity and forecasting outliers
events, discussed in Discussion. We also settle the basis for the formal
derivation of the model. We give an introduction to Maximum Entropy models,
the formalism used in the following sessions and the optimization algorithm
used. Afterward,we introduce the formalism to compute an approximated Log-
Likelihood. In Additional Informations can be found similar results for the
cities of Rome, Florence, and Turin.

In section Methods we describe the data and the variables in use. To find out
which are the important variables to represent, we do a historical analysis
and observe a dynamics of contraction and dilation of the time-shifted
correlations between different zones. With this, we understand the important
effects and and finally we define and derive the formulae for the ME model
describing our data.

We use MaxEnt inference to obtain the parameters (using gradient ascent
algorithm) and reproduce lag-correlations of definite positive time series. We
derive a model that is highly predictive at least as more sophisticated models
that take into account nonlinear correlations and have more parameters. We
also use the obtained statistical model to find extrema events (outliers, such
as strikes and bad weather days). As a result, we find the dynamics of cars’
presence in urban areas to be extremely rich and complex. We infer the
couplings parameters between the activity profiles of different areas and use
them to project the cars’ locations in time efficiently.

## Results

In this section, given a brief introduction to the method and the data
studied, we present the results. For a more in depth description, we refer to
sec. Methods.

### Data

The data we use is a normalised multivariate time-series representing the
parking activity of the different municipalities of the city. The procedure to
obtain the final form of the time-series in described in sec. Methods. In Fig.
1 we show an example of activity for two zones of our dataset for the
Metropolitan City of Milan.

Figure 1: Time series activity data for the city center (a) and for one of the
suburbs (b).

### Maximum Entropy Principle

The principle of maximum entropy states that the probability distribution
which best represents the current state of knowledge is the one with the
largest entropy, in the context of precisely stated prior data. According to
this principle, the distribution with maximal information entropy is the best
choice. The principle was first shown by E. T. Jaynes in two papers in the
late fifties where he emphasized a natural correspondence between statistical
mechanics and information theory [16, 17].

In particular, Jaynes offered a new and very general rationale on why the
Gibbsian method of statistical mechanics works. He showed that statistical
mechanics, and in particular Ensemble Theory, can be seen just as a particular
application of information theory, hence there is a strict correspondence
between the entropy of statistical mechanics and the Shannon information
entropy.

Maximum Entropy models unveiled interesting results throughout the years for a
large variety of systems, like flocking birds [18], proteins [19], brain [20]
and social systems [21].

We will then implement this approach to define the model of our real-world
system in the following sections. A more general introduction to the maximum
entropy formalism is out of scope here and we refer to the existing literature
for details [22, 23, 24, 25].

The probability distribution with Maximum Entropy $P_{ME}$ results from the
extreme condition of the so-called Generalized Entropy:

$\mathcal{S}\big{[}P\big{]}=S\big{[}P\big{]}+\sum_{k=1}^{K}\theta_{k}(\langle
O_{k}\rangle_{P(\underline{X})}-\langle O_{k}\rangle_{obs}),$ (1)

where

$S\big{[}P\big{]}=-\sum_{\underline{X}}P(\underline{X})\log(P(\underline{X}))$
(2)

is the Shannon Entropy of the probability distribution $P(\underline{X})$. The
maximum of the Generalized Entropy is the maximum of the Entropy of the model
when it is subject to constraints. Computing the functional-derivative (1)
with respect to $P(\underline{X})$ and equating to zero results in:

$P_{me}(\underline{X})=\frac{1}{Z(\underline{\theta})}\exp\Big{[}-\sum_{k=1}^{K}\theta_{k}O_{k}(\underline{X})\Big{]},$
(3)

where

$Z(\underline{\theta})=\int\limits_{\Omega}d\underline{X}\exp\Big{[}-\sum_{k=1}^{K}\theta_{k}O_{k}(\underline{X})\Big{]}$
(4)

is the normalization (making a parallel with statistical physics, can be
called Partition Function). $Z(\underline{\theta})$ is written as a sum if
$\Omega$ is discrete. Hence, the maximum entropy probability distribution is a
Boltzmann distribution in the canonical ensemble at temperature $T=1K$, with
effective Hamiltonian
$\mathcal{H}(\underline{X})=-\sum_{k=1}^{K}\theta_{k}O_{k}(\underline{X})$.

It must be noticed that the minimization of the generalized entropy is
equivalent to the maximization of the experimental average of the likelihood:

$\mathcal{S}\big{[}P\big{]}=\log
Z(\underline{\theta})-\sum_{k=1}^{K}\theta_{k}\langle
O_{k}\rangle_{e}=-\langle\log P_{me}\rangle_{e}=\frac{1}{M}\sum_{m=1}^{M}\log
P(\underline{X}^{(m)}).$ (5)

In other words, the $\theta_{k}$ are chosen by imposing the experimental
constraints on Entropy or, equivalently, by maximizing the global,
experimental likelihood according to a model with the constraints cited above.
Focusing on this last sentence, we can say that the optimal parameters of
$\theta$ (called effective couplings) can be obtained through Maximum
Likelihood, but only once one has assumed (by the principle of Maximum
Entropy) that the most probable distribution has the form of $P_{ME}$.

Given the generative model probability distribution of configurations
$P(\underline{X}\mid\underline{\theta})$ and its corresponding partition
function by $\log Z(\underline{\theta})$, the estimator of $\theta$ can be
found by maximizing the log-likelihood:

$\mathcal{L}(\underline{\theta})=\langle\log(P(\underline{X}\mid\underline{\theta}))\rangle_{data}=-\langle\mathcal{H}(\underline{X};\underline{\theta})\rangle_{data}-\log
Z(\underline{\theta}).$ (6)

Having chosen the log-likelihood as our cost function, we still need to
specify a procedure to maximize it with respect to the parameters.

One common choice that is widely employed when training energy-based models is
Gradient Descent [25] or its variations: optimization is taken w.r.t. the
gradient direction. Once one has chosen the appropriate cost function
$\mathcal{L}$, the algorithm calculates the gradient of the cost function
concerning the model parameters. The update equation is:

$\theta_{ij}\leftarrow\theta_{ij}-\eta_{ij}\frac{\partial\mathcal{L}}{\partial\theta_{ij}}.$
(7)

Typically, the difficulty in these kind of problems is to evaluate $\log
Z(\underline{\theta})$ and its derivatives. The reason is that the partition
function is rarely an exact integral, and it can be calculated exactly only in
a few cases. However, it is still possible to find ways to approximate it and
compute approximated gradients.

### Pseudo-Log-Likelihood (PLL) Maximization

Pseudo-likelihood is an alternative method compared to the likelihood function
and leads to the exact inference of model parameters in the limit of an
infinite number of samples [26, 27]. Let’s consider the log-likelihood
function $\mathcal{L}(\underline{\theta})=\langle\log
P(\underline{X}\mid\underline{\theta})\rangle_{data}$. In some cases we are
not able to compute the partition function $Z(\underline{\theta})$, but it is
possible to derive exactly the conditional probability of one component of
$\underline{X}$ with respect to the others, i.e.
$P(X_{j}|\underline{X_{-j}},\underline{\theta})$ where $\underline{X_{-j}}$
indicates the vector $\underline{X}$ without the $j$-th component.

In this case, we can write an approximated likelihood called Pseudo-log-
likelihood, which takes the form:

$\mathcal{L}(\underline{\theta})_{\textit{pseudo}}=\sum_{j}\langle\log
P(X_{j}|\underline{X_{-j}},\underline{\theta})\rangle_{data}.$ (8)

The model we will introduce in this work does not have an explicit form for
the partition function, but Eq. (8) and its derivatives can be calculated
exactly. Thus, the Pseudo-log-likelihood is a convenient cost function for our
problem.

### Definition of MaxEnt model

We have seen that the Maximum Entropy inference scheme requires the definition
of some observables that are supposed to be relevant to the system under
study. Being that one aims to predict the evolution of the different
$x_{i}(t)$, the most simple choice is to consider their correlations.

As a preliminary analysis, we study the correlation between the activity
$x_{i}(t)$ of most central zone within the city (highlighted in red) and all
the other zones with a shift in time (i.e. we correlate $x_{i}(t)$ with
$x_{j}(t-\delta)$ for al $j\neq i$ and some $\delta>0$). The measure of
correlation between two vectors $u$ and $v$ represented is $\frac{u\cdot
v}{{||u||}_{2}{||v||}_{2}}$, where ${||.||}_{2}$ measures the euclidean norm
defined in $\mathbb{R}^{n}$ as
$\left\|\bm{x}\right\|_{2}:=\sqrt{x_{1}^{2}+\cdots+x_{n}^{2}}$. We see that
the areas with large correlations vary with $\delta$ so that they clustered
around the central area for small values, and they become peripherals when
$\delta\sim 31$. When $\delta\sim 48$ they start to cluster around the central
area again. An opposite behavior is shown when repeating the same procedure
focusing on a peripheral zone. We perform a historical analysis of time-
shifted correlations, observing contraction and dilation w.r.t. different zone
in a one-day periodicity. Due to the correlation dependency on the time shift,
we have to include it in the observables used to describe the system. Hence,
we chose as observables all the shifted correlations between all the couples
of zones defined as,

$\langle
x_{i}(t)x_{j}(t-\delta)\rangle_{data}=\frac{1}{T-\Delta}\sum_{t=\Delta}^{T}x_{i}(t)x_{j}(t-\delta),$
(9)

for $i,j=1,....,N$ (with $N$ as the number of zones we took in consideration)
and $\delta=0,....,\Delta$. Another common choice is to also fix the average
value of all the variables of the system. We took $\langle
x_{i}\rangle_{data}=\frac{1}{T-\Delta}\sum_{t=\Delta}^{T}x_{i}(t)$ and
$\langle
x^{2}_{i}\rangle_{data}=\frac{1}{T-\Delta}\sum_{t=\Delta}^{T}x^{2}_{i}(t)$.

From these, we obtain the equation for the probability:

$P(x(t),\dots,x(t-\Delta))=\frac{1}{Z}\exp\left[-\sum_{t=\Delta}^{T}\sum_{i}a_{i}x_{i}^{2}(t)+\sum_{t=\Delta}^{T}\sum_{i}h_{i}x_{i}(t)+\sum_{t=\Delta}^{T}\sum_{\delta=1}^{\Delta}J_{ij}^{\delta}x_{i}(t)x_{j}(t-\delta)\right],$
(10)

where $a_{i}$ is the $i$-th component’s standard deviation, $h_{i}$ is its
mean and $J_{ij}^{\delta}$ are the time shifted interactions.

Writing $v_{i}(t)=h_{i}+\sum_{\delta}J_{ij}^{\delta}x_{j}(t-\delta)$ one can
obtain:

$P(x(t),\dots,x(t-\Delta))=\frac{1}{Z}\prod_{t=\Delta}^{T}\exp\left[-\sum_{i}a_{i}x_{i}^{2}(t)+\sum_{i}v_{i}(t)x_{i}(t)\right].$
(11)

### Time series Forecasting

We trained the model for $100000$ steps using several hyperparamenters
($\Delta=\\{24,36,48,72\\}$ and $\lambda=\\{0.001,0.004,0.005,0.006,0.01\\}$,
see tab.1).

To quantify the predictive capabilities of our models, we used the Mean
Absolute Error (MAE), the Mean Squared Error (MSE) and the Coefficient of
Determination, also known as $R^{2}$ [28].

Table 1: $R^{2}$ values w.r.t. $\Delta$ and $\lambda$ values, resulting from $100000$ steps of training.  $\Delta$ |  $\lambda$ |  $R^{2}$
---|---|---
24 |  0.001 |  0.855316
24 |  0.004 |  0.854723
24 |  0.005 |  0.860069
24 |  0.006 |  0.853689
24 |  0.01 |  0.832542
36 |  0.001 |  0.854458
36 |  0.004 |  0.85442
36 |  0.005 |  0.853946
36 |  0.006 |  0.853447
36 |  0.01 |  0.832672
48 |  0.001 |  0.860072
48 |  0.004 |  0.861741
48 |  0.005 |  0.869466
48 |  0.006 |  0.861626
48 |  0.01 |  0.833099
72 |  0.001 |  0.85655
72 |  0.004 |  0.859633
72 |  0.005 |  0.867582
72 |  0.006 |  0.859632
72 |  0.01 |  0.858497

We typically find good predictive capabilities on the test set ( as in fig. 4
and the exact values for each predictor are found in fig 1). To check for the
best value, we used the average value at each time step during the day in the
training set as a baseline model. For $\Delta=48$ (exactly the day-
periodicity) and $\lambda=0.005$, we find our model’s best performance, so we
will use these values.

To assess the predictive power of our approach against more complex models, we
compared it with standard statistical inference and Machine Learning methods.
In particular:

  1. 1.

a SARIMA model [29]: we performed a grid search over the parameters using AIC
[30] as metrics, the best predictive model is, referring to standard
litherature on the topic, $\\{p=2,d=0,q=1,P=2,D=1,Q=2,m=48\\}$. In tab.2 the
table of the grid search.

  2. 2.

NN and LSTM: each of them has 3 hidden layers with a
$N_{nodes}=\\{64,128,256,512\\}$. Between the layers, we performed dropout
techniques and standard non-linear activation functions (ReLU). Performing a
grid-search over the other Machine Learning hyperparameters, the layers with
$128$ and $256$ nodes perform better.

Fig. 2 shows that our model outperforms SARIMA and performs as well as the
Machine Learning ones. Being that the LSTM one always performs slightly better
than the FeedForward, we will use the first as a reference. Fig. 3 shows this
in more detail: SARIMA fails to reproduce variations and it is mostly focused
on seasonality. Our model works as well as NNs but with two orders of
magnitude of fewer parameters ($109$ w.r.t., for the $128$-nodes
configuration, $39040$ for the FeedForward).

In fig. 7 we show data and results for the cities of Rome, Florence, and
Turin.

### Extreme events prediction

As already pointed out, prediction is a possible application of our procedure,
and we compared the results with state-of-the-art forecasting techniques.
Another possible application is the detection of outliers in the data. We see
from Fig. 1 that our car sharing data exhibits quite regular patterns on
different days. These patterns can be perturbed by a wide variety of external
events, influencing the demand for car-sharing vehicles as well as the traffic
patterns in the city. If we compute the log-likelihood from equation (17)
restricting to one of the perturbed days, we would expect it to be
particularly low, implying that the model has lower predictive power than
usual. Hence, we can use the log-likelihood $\mathcal{L}$ from equation (17)
as a score in a binary classification task. After training the model with a
specific parameter choice, we consider a set of days in the test dataset,
$d_{1}\dots d_{D}$. Each one is a set of consecutive time steps:
$d_{i}=[t^{d_{i}}_{1},t^{d_{i}}_{2}]$ with $t^{d_{i}}_{1}<t^{d_{i}}_{2}$. The
log-likelihood of a single day is just

$\mathcal{L}_{pseudo}(d_{i})=\frac{1}{t^{d_{i}}_{2}-t^{d_{i}}_{1}-\Delta}\sum_{t=t^{d_{i}}_{1}-\Delta}^{t^{d_{i}}_{2}}\log
P(x(t)|x(t-1),\dots,x(t-\Delta)).$ (12)

We can now assign a label $L(d_{i})=0$ if $d_{i}$ is a standard day and
$L(d_{i})=1$ if it is a day where a certain external event occurred. We
considered two types of external events:

  1. 1.

Weather conditions: on that day there was fog, a storm or it was raining.

  2. 2.

Strikes: in that day there was a strike in the city capable of affecting car
traffic. We considered Taxi strikes, Public Transportation strikes, and
general strikes involving all working categories.

We had $50\%$ of days labeled as 1 among the days in the test dataset. In fig.
5 we show the ROC curve for the classification of such days, using the log-
likelihood trained with the parameters $\\{\Delta=48,\lambda=0.006\\}$. The
Area under the ROC curve (AuROC) is 0.81, indicating that $\mathcal{L}$ is a
good indicator of whether an external event occurred on a specific day.

Figure 2:  Model comparison using different metrics. We see similar results
between the MaxEnt model and the Deep Learning ones.  Figure 3: Model
comparison: taking one zone’s time series, we see that SARIMA models are good
predictors of periodicity, but the other models perform better when the change
between one period and another is high.  Figure 4: R2 MaxEnt prediction over
the test set. To have r2=1, the point should align on the bisector line,
having a perfect predictor.  Figure 5: ROC curve for the detection of outliers
(bad weather conditions and strikes) for $\lambda=0.006$ and $\Delta=48$.

## Discussion

In this work, we addressed the problem of building statistical and machine
learning models to predict mobility patterns within the Milan Metropolitan
Area. We focused on the prediction of how many car-sharing vehicles are parked
within a specific area at different times. To do this, we analyzed a
longitudinal dataset of parking events coming from the APIs of the Italian
car-sharing company “Enjoy”. The processed data consisted of a collection of
time series representing the number of parked vehicles in a given area at a
given time.

To predict the evolution of these time series, we used a Maximum Entropy (ME)
approach that requires the identification of relevant observables from the
data and use these observables to define probability density functions
depending on some parameters.

Maximum Entropy modeling has proven to be very effective at studying and
determining the activity patterns inside the city (Fig. 3).

We compared our model with other models specifically built for time series
forecasting. In particular, we used a SARIMA model, a Feed-Forward Neural
Network, and a Long Short-Term Memory Neural Network.

Maximum Entropy models outperformed SARIMA models w.r.t. all the metrics used
in the evaluation. Our model is as predictive as LSTMs Neural Networks, but
using two orders of magnitude fewer parameters , interpretability of the
parameters, and the possibility of performing different studies using the same
model.

Finally, we used the statistical model also to identify extreme events
affecting urban traffic, such as bad weather conditions or strikes.

In conclusion, this work is a first attempt at modeling the inherently complex
interactions and causalities that shape the fluxes of people and vehicles
within urban spaces. We showed that the interplay of the activity of different
areas at a different time is far from trivial. Indeed, we found the models
derived within the ME framework to accurately reproduce the empirical time
series using a mixture of correlations between different areas at different
times.

Moreover, our finding that linear correlations are enough to predict the
mobility patterns (wrt neural networks) is an interesting result that requires
additional investigations.

In Additional Information we show the prediction ability on other three
cities: Rome, Florence, and Turin.

Given our results, several research directions can be taken to improve and
extend the results:

  * •

a more extensive study on the effects of seasonality could help in building
better models. Specific season-dependent models could be built by taking into
account larger datasets.

  * •

Evaluate how the prediction ability of the Maximum Entropy models is dependent
upon the structure of the city or the resolution of the time series.

  * •

Entangling mobility patterns with other city indicators, such as socio-
political disparities and economic indicators, can lead to a better model that
depicts human distribution.

  * •

The inclusion of non-linear interaction in the ME models could be hard if made
by defining standard observables. Instead, hidden variables approaches could
be take into account, e.g. Restricted Boltzmann Machines [31].

  * •

The ME models could be adapted to perform other anomaly detection tasks, i.e.
identifying parts of the time series which are not standard fluctuations of
the system [32, 33].

## Methods

### Dataset

The dataset used in this work is a log of the cars’ locations, belonging to
Enjoy, a popular Italian car-sharing service. Although the service is active
in six major Italian cities, we focus our analysis on the city of Milan. The
data has been collected via the Enjoy APIs, scraping the information through a
client account. We obtained information about the area of service (i.e., the
limit where people can start and end their ride), the location of vehicles,
and the points of interest (POIs) related to the service (such as fuel
stations and reserved parkings) from the API endpoints. The endpoints used
have base URL https://enjoy.eni.com/ajax/ and the following endpoints: missing
for the area, missing for cars, and missing for points of interests.

We collected the data recording all the movements of Enjoy cars within the
city of Milan in $2017$. Through the scraping procedure, we collected a series
of events, that are divided into two categories: parking events and travel
events of each Enjoy car we could observe. Each event comes with a set of
information such as the car plate, the time at which the event occurred (with
the precision of seconds), and the latitude and longitude of the parking spot
in the case of parking events. Starting and arrival points in the case of
travel events are also recorded. In the rest of this article, we will focus
mainly on parking events, but our methods and analyses could be also applied
to study the volume of travel events. The final aim of the work was to predict
the number of vehicles parked at a given time in each and every district of
the city. To this end, we divided the area of service into the different
municipalities (we will call them zones). Different tessellation strategies
(e.g., using hexagons with a $900m$ side) have been considered, but the
municipality division provides more statistichally stable training data, and
thus more precise models.

For each type of event (parked car) that we collected, we defined the activity
of a zone as the number of events (i.e., cars parked) occurring there in a
certain fixed amount of time $\delta t$. Indicating each zone with the index
$i$, we have obtained for each and every one of them a discrete time series
$x_{i}(t)$ representing the activity of the zone $i$ in the time frame
$[t,t+\delta t]$. We will indicate with $t=T$ the last time frame observed in
the data. In the following, we will use a $\delta t=1800\,$s, corresponding to
$30$ minutes. This time bin width has been chosen so to have a stable average
activity $\langle x_{i}(t)\rangle$ and a stable variance $\sigma(x_{i}(t))$
over all the zones. Indeed, narrower or larger time bins feature unstable mean
and/or variance in time: the $30-minutes$ binning is significantly stable
throughout all the observation period. This characteristic helps the model to
generalize and predict in a better way as the distributions of the modeled
quantities are not changing too much in time.

Some of the models we would like to implement require real variables
distributed over $\mathbb{R}$. However, the $x_{i}(t)$ activity we have
defined so far belongs to $\mathbb{N}$ by definition. For this reason, we
defined another kind of variable dividing the activity by the standard
deviation:

$z=\frac{x}{\sigma},$ (13)

where $x$ is the original activity data, and $\sigma$ is the standard
deviation of the activity in time. In this case, for the population we take
into account the typical day: the std is taken w.r.t. the same time bin for
every day. So, Eq. (13) gets into:

$z_{i}(t)=\frac{x_{i}(t)}{\sigma_{i}(t_{\delta t})},$ (14)

where $t_{\delta t}=t\mod(\delta t)$ is the integer division of t by the
$\delta t$ time bin width, $\mu_{i}(t_{\delta t})$ is the average over all the
activities of the zone $i$ at times with the same value of $t\mod\delta t$,
and $\sigma_{i}(t_{\delta t})$ is their standard deviation. To keep the
notation formalism, we will indicate $z_{i}(t)$ with $x_{i}(t)$, keeping in
mind that now $x$ refers to a normalized activity.

From now on, we will work on the normalized $x_{i}(t)$: these indicate how
much a given area is more ”active” - i.e., has more cars parked - concerning
the average volume of that time bin $t$ (typical day activity) by weighting
this fluctuation with the standard deviation observed for that zone-hour
volume. In this way, we can compare the signal of areas with high fluctuations
and high activities with less frequented areas around the city.

The final step when processing data for inference and Machine Learning is to
divide the data into a train set and a test set. The models will be then
trained using the first dataset, and their precision will be tested on the
second one so to check their ability to generalize to unseen data. Different
kinds of splittings have been done, like random splitting or taking the first
part of the time series as training and the last one as a test. Similar
results are obtained.

As a final remark, in the following we will indicate with $t$ the time bin of
activity, i.e. $x_{i}(t)$ will indicate the activity of the zone $i$ in a time
range $[t\delta t,(t+1)\delta t]$.

### Pseudo-Log-Likelihood Maximization

The Boltzmann probability is defined over the whole time series of all the
zones, i.e.

$P\Big{(}x(t),x(t-1),...,x(t-\Delta)\Big{)}\propto\exp\Big{(}\mathbf{H}(x(t),x(t-1),...,x(t-\Delta))\Big{)}.$
(15)

From this, it is straightforward to define the conditional probability of one-
time step $x(t)$ concerning all the previous ones:

$P(x(t)\mid
x(t-1),...,x(t-\Delta))=\frac{P\Big{(}x(t),x(t-1),...,x(t-\Delta)\Big{)}}{P\Big{(}x(t-1),...,x(t-\Delta)\Big{)}}.$
(16)

Using equation (8), we can define the Pseudo-Log-Likelihood as:

$\mathcal{L}_{pseudo}=\frac{1}{T-\Delta}\sum_{t=\Delta}^{T}\log
P(x(t)|x(t-1),\dots,x(t-\Delta)).$ (17)

Here, using Eq. (16) and substituting the functional form of the two total
probabilities, we obtain

$P(x(t)|x(t-1),\dots,x(t-\Delta))=\prod_{i}\frac{1}{Z_{i}(t)}\exp(-a_{i}x_{i}^{2}(t)+x_{i}(t)v_{i}(t)),$
(18)

with

$\displaystyle
Z_{i}(t)=\frac{1}{\sqrt{a_{i}}}\exp\left(\frac{v_{i}(t)^{2}}{4a_{i}}\right)\int_{-\infty}^{v_{i}(t)/2\sqrt{a_{i}}}dze^{-z^{2}}=\frac{1}{\sqrt{a_{i}}}\exp\left(\frac{v_{i}(t)^{2}}{4a_{i}}\right)I\left(\frac{v_{i}(t)}{2\sqrt{a_{i}}}\right),$
(19) $\displaystyle
v_{i}(t)=\sum_{n=1}^{d}\sum_{\delta}(J_{n}^{\delta}x^{n}(t-\delta))_{i}+h_{i}.$

Substituting in eq. (17), we get:

$\mathcal{L}_{pseudo}=\frac{1}{T-\Delta}\sum_{t=\Delta}^{T}\sum_{i}-a_{i}x_{i}^{2}(t)+x_{i}(t)v_{i}(t)-\log
Z_{i}(t)$ (20)

and we can calculate the gradients of the $\mathcal{L}_{pseudo}$ w.r.t. the
parameters:

$\displaystyle\frac{\partial\mathcal{L}_{pseudo}}{\partial
a_{i}}=\frac{1}{T-\Delta}\sum_{t}-x_{i}^{2}(t)-\frac{\partial\log
Z_{i}(t)}{\partial a_{i}},$ (21)
$\displaystyle\frac{\partial\mathcal{L}_{pseudo}}{\partial
J^{\delta}_{ij}}=\frac{1}{T-\Delta}\sum_{t}x_{i}(t)x_{j}(t-\delta)-\langle
x_{i}(t)\rangle x_{j}(t-\delta),$
$\displaystyle\frac{\partial\mathcal{L}_{pseudo}}{\partial
h_{i}}=\frac{1}{T-\Delta}\sum_{t}x_{i}(t)-\langle x_{i}(t)\rangle,$

where

$\langle
x_{i}(t)\rangle=\frac{v_{i}(t)}{2a_{i}}+\frac{1}{2\sqrt{a_{i}}I(\frac{v_{i}}{2\sqrt{a_{i}}})}\exp\left(-\frac{v_{i}(t)^{2}}{4a_{i}}\right)$
(22)

and

$\frac{\partial\log Z_{i}(t)}{\partial
a_{i}}=-\frac{1}{2a_{i}}+\frac{-v_{i}^{2}(t)}{4a_{i}^{2}}+\frac{1}{I(\frac{v_{i}}{2\sqrt{a_{i}}})}\exp\left(-\frac{v_{i}(t)^{2}}{4a_{i}}\right)\left(\frac{-v_{i}(t)}{4a_{i}^{3/2}}\right).$
(23)

The fact that the gradients and the cost function can be computed exactly
makes the inference of the parameters relatively easy with respect to other
cases where they need to be approximated.

Once the parameters of the model have been inferred with some method, it is
possible to use it to predict the temporal evolution of the normalized
activities of the system. Given some certain state of the system unit time
$t-1$, i.e. $(x(t-1),\dots x(t-\Delta))$ (past time steps further than
$\Delta$ from $t$ are not relevant), we can use equation (18) to predict the
next step $x(t)$. Since the probability in (18) is a normal distribution whose
average is completely defined by $(x(t-1),\dots x(t-\Delta))$, the best
prediction of $x_{i}(t)$ is the average of the distribution $\langle
x_{i}(t)\rangle=\frac{1}{2a_{i}}\left(v_{i}(t)+\frac{1}{Z_{i}(t)}\right)$.

In other words, we are using the generative model to make a discrimination
task. This makes possible the comparison of this model with standard machine
learning ones, by comparing their precision in the prediction of the time
series. To avoid over-fitting, we used L1-regularization [34, 35]. An in-depth
description of the technique can be found in the references. In practice, the
cost function that has to be optimized is:

$C(\theta)=\log P(\theta|X)=\log
P(X|\theta)-\lambda\sum_{i}|\theta_{i}|+cost.$ (24)

The first term of this sum is the Log-likelihood; instead, the second term is
the regularization term.

If the gradient is performed, we obtain:

$\frac{\partial C(\theta)}{\partial\theta_{i}}=\frac{\partial\log
P(X|\theta)}{\partial\theta_{i}}-\lambda sign(\theta_{i}),$ (25)

where $sign$ is the sign function.

The training curves show no sign of overfitting,as the log-likelihood
asymptotically stabilizes for the validation and train sets (see Fig. 6 in
Appendix).

##  Appendix A Appendixes

Table 2: Grid search over the SARIMA parameters p |  d |  q |  P |  D |  Q |  m |  AIC
---|---|---|---|---|---|---|---
1 |  0 |  1 |  0 |  1 |  1 |  48 |  inf
0 |  0 |  0 |  0 |  1 |  0 |  48 |  1243.245
1 |  0 |  0 |  1 |  1 |  0 |  48 |  359.487
0 |  0 |  1 |  0 |  1 |  1 |  48 |  inf
0 |  0 |  0 |  0 |  1 |  0 |  48 |  1241.606
1 |  0 |  0 |  0 |  1 |  0 |  48 |  615.697
1 |  0 |  0 |  2 |  1 |  0 |  48 |  285.295
1 |  0 |  0 |  2 |  1 |  1 |  48 |  270.777
1 |  0 |  0 |  1 |  1 |  1 |  48 |  279.520
1 |  0 |  0 |  2 |  1 |  2 |  48 |  237.403
1 |  0 |  0 |  1 |  1 |  2 |  48 |  inf
0 |  0 |  0 |  2 |  1 |  2 |  48 |  inf
2 |  0 |  0 |  2 |  1 |  2 |  48 |  236.616
2 |  0 |  0 |  1 |  1 |  2 |  48 |  inf
2 |  0 |  0 |  2 |  1 |  1 |  48 |  270.899
2 |  0 |  0 |  1 |  1 |  1 |  48 |  280.519
2 |  0 |  1 |  1 |  1 |  2 |  48 |  inf
2 |  0 |  1 |  2 |  1 |  1 |  48 |  273.366
2 |  0 |  1 |  1 |  1 |  1 |  48 |  276.647
1 |  0 |  1 |  2 |  1 |  2 |  48 |  237.017
2 |  0 |  2 |  2 |  1 |  2 |  48 |  237.784
1 |  0 |  2 |  2 |  1 |  2 |  48 |  236.860
2 |  0 |  1 |  2 |  1 |  2 |  48 |  233.925

Figure 6: Log-Likelihood training for different values of $\Delta$. Figure 7:
Time series activity data for the cities of Rome (a) having $R^{2}=0.816849$ ,
Florence (b) having $R^{2}=0.748687$ and Turin (c) having $0.727507$ .

## References

  * Bettencourt  _et al._ [2007] L. M. Bettencourt, J. Lobo, D. Helbing, C. Kühnert, and G. B. West, Growth, innovation, scaling, and the pace of life in cities, Proceedings of the national academy of sciences 104, 7301 (2007).
  * Bettencourt  _et al._ [2010] L. M. Bettencourt, J. Lobo, D. Strumsky, and G. B. West, Urban scaling and its deviations: Revealing the structure of wealth, innovation and crime across cities, PloS one 5, e13541 (2010).
  * Verbavatz and Barthelemy [2020a] V. Verbavatz and M. Barthelemy, The growth equation of cities, Nature 587, 397 (2020a).
  * Rogov and Rozenblat [2018] M. Rogov and C. Rozenblat, Urban resilience discourse analysis: towards a multi-level approach to cities, Sustainability 10, 4431 (2018).
  * Barthélemy [2011] M. Barthélemy, Spatial networks, Physics reports 499, 1 (2011).
  * Raimbault [2020] J. Raimbault, Hierarchy and co-evolution processes in urban systems, arXiv preprint arXiv:2001.11989 (2020).
  * Wise  _et al._ [2016] S. Wise, A. Crooks, and M. Batty, Transportation in agent-based urban modelling, in  _International Workshop on Agent Based Modelling of Urban Systems_ (Springer, 2016) pp. 129–148.
  * Noulas  _et al._ [2012a] A. Noulas, S. Scellato, R. Lambiotte, M. Pontil, and C. Mascolo, A tale of many cities: universal patterns in human urban mobility, PloS one 7, e37027 (2012a).
  * Davies and Maré [2021] B. Davies and D. C. Maré, Relatedness, complexity and local growth, Regional Studies 55, 479 (2021).
  * Tainter [1988] J. Tainter,  _The collapse of complex societies_ (Cambridge university press, 1988).
  * Schläpfer  _et al._ [2021] M. Schläpfer, L. Dong, K. O’Keeffe, P. Santi, M. Szell, H. Salat, S. Anklesaria, M. Vazifeh, C. Ratti, and G. B. West, The universal visitation law of human mobility, Nature 593, 522 (2021), number: 7860 Publisher: Nature Publishing Group.
  * Bettencourt [2021] L. M. Bettencourt,  _Introduction to urban science: evidence and theory of cities as complex systems_ (MIT Press, 2021).
  * Verbavatz and Barthelemy [2020b] V. Verbavatz and M. Barthelemy, The growth equation of cities, Nature 587, 397 (2020b).
  * Dong  _et al._ [2022] L. Dong, P. Santi, Y. Liu, S. Zheng, and C. Ratti, The universality in urban commuting across and within cities (2022), arXiv:2204.12865 [physics, q-fin].
  * Noulas  _et al._ [2012b] A. Noulas, S. Scellato, R. Lambiotte, M. Pontil, and C. Mascolo, A Tale of Many Cities: Universal Patterns in Human Urban Mobility, PLoS ONE 7, e37027 (2012b).
  * Jaynes [1957a] E. T. Jaynes, Information theory and statistical mechanics, Physical review 106, 620 (1957a).
  * Jaynes [1957b] E. T. Jaynes, Information theory and statistical mechanics. ii, Physical review 108, 171 (1957b).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Statistical mechanics for natural flocks of birds, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Ekeberg  _et al._ [2013] M. Ekeberg, C. Lövkvist, Y. Lan, M. Weigt, and E. Aurell, Improved contact prediction in proteins: using pseudolikelihoods to infer potts models, Physical Review E 87, 012707 (2013).
  * Tang  _et al._ [2008] A. Tang, D. Jackson, J. Hobbs, W. Chen, J. L. Smith, H. Patel, A. Prieto, D. Petrusca, M. I. Grivich, A. Sher, _et al._ , A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro, Journal of Neuroscience 28, 505 (2008).
  * Gresele and Marsili [2017] L. Gresele and M. Marsili, On maximum entropy and inference, Entropy 19, 642 (2017).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to machine learning for physicists, Physics reports (2019).
  * Giffin [2009] A. Giffin, Maximum entropy: the universal method for inference, arXiv preprint arXiv:0901.2987 (2009).
  * Ibáñez-Berganza  _et al._ [2019] M. Ibáñez-Berganza, G. L. Lancia, A. Amico, B. Monechi, and V. Loreto, Unsupervised inference approach to facial attractiveness, arXiv preprint arXiv:1910.14072 (2019).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep learning_ (MIT press, 2016).
  * Arnold and Strauss [1991] B. C. Arnold and D. Strauss, Pseudolikelihood estimation: some examples, Sankhyā: The Indian Journal of Statistics, Series B , 233 (1991).
  * Nguyen  _et al._ [2017] H. C. Nguyen, R. Zecchina, and J. Berg, Inverse statistical problems: from the inverse ising problem to data science, Advances in Physics 66, 197 (2017).
  * Devore [2011] J. L. Devore,  _Probability and Statistics for Engineering and the Sciences_ (Cengage learning, 2011).
  * Tseng and Tzeng [2002] F.-M. Tseng and G.-H. Tzeng, A fuzzy seasonal arima model for forecasting, Fuzzy Sets and Systems 126, 367 (2002).
  * Burnham and Anderson [2004] K. P. Burnham and D. R. Anderson, Multimodel inference: understanding aic and bic in model selection, Sociological methods & research 33, 261 (2004).
  * Sutskever  _et al._ [2009] I. Sutskever, G. E. Hinton, and G. W. Taylor, The recurrent temporal restricted boltzmann machine, in  _Advances in neural information processing systems_ (2009) pp. 1601–1608.
  * Fiore  _et al._ [2013] U. Fiore, F. Palmieri, A. Castiglione, and A. De Santis, Network anomaly detection with the restricted boltzmann machine, Neurocomputing 122, 13 (2013).
  * Sun  _et al._ [2004] B. Sun, F. Yu, K. Wu, and V. C. Leung, Mobility-based anomaly detection in cellular mobile networks, in  _Proceedings of the 3rd ACM workshop on Wireless security_ (2004) pp. 61–69.
  * Bühlmann and Van De Geer [2011] P. Bühlmann and S. Van De Geer,  _Statistics for high-dimensional data: methods, theory and applications_ (Springer Science & Business Media, 2011).
  * Bishop [2006] C. M. Bishop,  _Pattern recognition and machine learning_ (springer, 2006).