# Equator to Pole Solar Chromospheric Differential Rotation using Ca-K
Features Derived from Kodaikanal Data

Hema Kharayat Department of Physics
M. L. K. P. G. College
Balrampur, 271201, Uttar Pradesh, India  Indian Institute of Astrophysics
Bangalore, 560034, India  Jagdev Singh  Indian Institute of Astrophysics
Bangalore, 560034, India  Muthu Priyal Indian Institute of Astrophysics
Bangalore, 560034, India  B. Ravindra Indian Institute of Astrophysics
Bangalore, 560034, India

###### Abstract

Differential rotation is one of the basic characteristics of the Sun, and it
plays an important role in generating the magnetic fields and its activities.
We investigated rotation rate using chromospheric features such as plages,
enhanced network (EN), active network (AN), and quiet network (QN) separately
(for the first time). The digitized Ca-K images from Kodaikanal Observatory
for 1907–1996 are used to study rotation over 0–80 degrees latitudes at an
interval of 10∘. We find that plages and all types of networks exhibit the
differential rotation of the chromosphere. Furthermore, the rotation rate
shows a decreasing pattern as one move from the equator to the higher polar
latitudes for all the features used in the study. At the equator the rotation
rate (rotation period) is obtained to be $\sim$13.98 deg day-1 (25.74 days),
$\sim$13.91 deg day-1 (25.88 days), $\sim$13.99 deg day-1 (25.74 days), and
$\sim$14.11 deg day-1 (25.51 days) for plage, EN, AN, and QN areas,
respectively. By analyzing how the area of chromospheric features varies over
time, we can effectively map the Sun’s rotation rate at all latitudes,
including the polar regions. Interestingly, both plages and small-scale
networks exhibit similar differential rotation rate. This suggests these
features likely rooted at the same layer below the visible surface of the Sun.
Therefore, the long-term Ca-K data is very useful to study the solar rotation
rate at all latitudes including the polar regions.

Solar chromosphere (1479) — plages(1240) — Solar rotation(1524) — Solar
differential rotation(1996)

††journal: APJL

##  1 Introduction

Solar differential rotation is a topic of investigation among researchers due
to its close connection with solar dynamo theory (Babcock, 1961; Dikpati &
Gilman, 2006; Chu et al., 2010; Javaraiah, 2020; Hotta & Kusano, 2021). The
interaction of the Sun’s differential rotation with the magnetic field plays a
fundamental role in generating all solar activities (Javaraiah et al., 2005a).
Therefore, studying solar differential rotation is essential for better
understanding the physical mechanisms behind the various solar activities and
cycles.

Solar differential rotation has been studied for a long time using various
forms of solar activity indexes. Sunspots are one of the oldest recorded solar
parameters and have been widely used for the measurement of the rotation rate
of the solar atmosphere (Lustig, 1983; Gupta et al., 1999; Javaraiah et al.,
2005b, and references therein). Ruždjak et al. (2017) analysed sunspot groups
to find a relationship between the solar rotation and activity. They found
that the Sun rotates more differentially at the minimum than at the maximum of
activity. By analysing the sunspot data from the Royal Greenwich Observatory
and the USAF/NOAA, Zhang et al. (2013) found evidence of an anticorrelation of
the rotation in the two hemispheres and oscillatory behavior of the asymmetry
at a period of 80–90 years. They also found that the north-south asymmetry of
solar rotation have an inverse relationship with the area of large sunspots.
Jha et al. (2021) measured the solar rotation profile using the white-light
sunspot data from Kodaikanal Observatory. They found no variations in rotation
rate between activity extremes, i.e., solar maxima and minima. To study the
coronal rotation, Insley et al. (1995) considered coronal holes as tracers of
the differential rotation. They noticed that the mid-latitude corona rotates
more rigidly than the photosphere, but still exhibits significant differential
rotation. Vats et al. (2001) derived the solar coronal rotation from the disk-
integrated simultaneous daily measurements of solar flux at eleven radio
frequencies in the range of 275–2800 MHz. They found that the sidereal
rotation period decreases with the altitude. Wöhl et al. (2010) determined the
solar rotation velocity by tracing small bright coronal structures in SOHO-EIT
images. They found that it has a rotation velocity similar to those obtained
by small photospheric magnetic features. The differential rotation rate of the
solar corona at low-latitudes is found to be consistent with the rotational
profile of the near-surface convective zone of the Sun (Mancuso et al., 2020).
By analysing the observations from the Atmospheric Imaging Assembly (AIA)
telescope on board the Solar Dynamics Observatory (SDO) at different
wavelengths, Sharma et al. (2020) found that the sidereal rotation periods of
different coronal layers decrease with increasing temperature (or height).

The Ca II K line images serve as an excellent diagnostic tool to investigate
the rotation of the solar chromosphere. A wide range of previous studies are
available for the chromospheric differential rotation (Livingston, 1969;
Belvedere et al., 1977; Antonucci et al., 1979a, b; Singh & Prabhu, 1985;
Ternullo, 1986, 1987; Li et al., 2020; Xu et al., 2020; Wan & Gao, 2022; Li &
Xu, 2024, and references therein). Using the total plage area integrated over
the entire visible hemisphere of the Sun, Singh & Prabhu (1985) found a
dominant periodicity of 7 years in rotation rate. Bertello et al. (2020)
analyzed the full-disk Ca II K images obtained from Mount Wilson Observatory.
They found that temporal variations in full-disk chromospheric activity show
the signature of the 11 year solar cycle. Also, they noticed that there is no
indication of any detectable periodicity in the temporal behavior of the
orthogonalized rotation rate coefﬁcients. Schroeter & Woehl (1975) and
Belvedere et al. (1977) investigated differential rotation of Ca II Networks
and Plages, respectively. They found that the chromosphere rotates 5$\%$
faster than the photosphere. Li et al. (2020) and Xu et al. (2020) utilized
data in He I and Mg II lines and reported a faster rotation of the
chromosphere. In an attempt to measure the chromospheric rotation rate, Wan &
Li (2022) utilize Ca II K filaments and suggested faster rotation of
chromosphere. Recently, Mishra et al. (2024) utilized Ca II K data obtained
from Kodaikanal Solar Observatory and reported that the chromospheric plages
exhibit an equatorial rotation rate 1.59$\%$ faster than the photosphere. At
the chromospheric layer, the short-lived features of Ca II K rotate at the
same rate as the chromospheric plasma (Antonucci et al., 1979a). Recently, Wan
et al. (2023a) analyzed synoptic maps of Ca II K-normalized intensity to
investigate the long-term variation of the quiet chromospheric differential
rotation. They found that the rotation rate is smaller for the quiet
chromosphere than for the chromosphere on the whole.

Hathaway & Rightmire (2011) measure differential rotation of the small
magnetic elements using the data from Michelson Doppler Interferometer (MDI)
instrument on board Solar and Heliospheric Observatory (SOHO) spacecraft
(Scherrer et al., 1995). They found that the differential rotation vary
systematically over the solar cycle and the differential rotation is weaker at
maximum than at minimum. By tracking the motions of individual magnetic
features from a single month of observations from a Helioseismic and Magnetic
Imager (HMI) aboard SDO, Lamb (2017) has arrived at high-precision
measurements of solar rotation. Imada & Fujiyama (2018) studied the dependence
of the surface flow velocity on the magnetic field strength. They found that
magnetic elements with a strong (active region remnants) and weak (solar
magnetic networks) magnetic field show a faster and slower rotation speed,
respectively. Hathaway et al. (2022) measured differential rotation and
meridional flow in the Sun’s surface shear layer by tracking the motions of
the magnetic network using magnetic pattern tracking on magnetograms from the
MDI on board SOHO spacecraft and from HMI on board SDO during 1996-2022. They
found that both, differential rotation and meridional flow, vary in strength
with depth. The rotation rate increases inward while the meridional flow
weakens inward.

In this way, many authors have studied the solar rotation rate by using the
sunspot numbers and others using the whole disk solar chromospheric and
coronal features but not the small-scale networks. In a study of differential
rotation by tracking the motions of the magnetic network, Hathaway et al.
(2022) found slower flows on the poleward sides of the active latitudes and
faster flows equatorward. In our study, we also found a decreasing trend in
the differential rotation, of both plages and small-scale networks, towards
the poles. Further, most of the investigations made in the past have
considered only up to mid-latitudes (up to $\pm$ 50∘) for the study of
differential rotation. Recently, Singh et al. (2021) have generated a uniform
time series of Ca-K images using the Equal Contrast technique (ECT), which
permitted the extraction of the reliable values of areas of small-scale
features in daily available Ca-K data. The uniqueness of the present work lies
in the study of solar differential rotation by using various small-scale solar
chromospheric features, viz plage, Enhanced Network (EN), Active Network (AN),
and Quiet Network (QN) areas. Furthermore, with this data, we are fortunate
enough to extend the study of solar rotation profile up to higher polar
latitudes (up to $\pm$ 80∘).

In this paper, we present the long-term measurements (1907–1996) of average
solar rotation rate as a function of heliographic latitude with the help of
the fast Fourier transform method. We utilize the digitized data of plage area
and small-scale networks obtained from the Kodaikanal Observatory. The paper
is structured as follows. Section 2 provides the information about the data
and method used for the present study. In Section 3, we discuss the results of
the investigation. The paper ends with conclusions about the results obtained.

##  2 data and Method

We utilize digitized data from the Kodaikanal Observatory (KO) for the present
work. The daily Ca-K spectroheliograms obtained at KO with some data gaps for
1907–2007, were digitized with a pixel resolution of 0.86 arcsec and 16 bit
readout using a CCD detector of 4K$\times$4K format. Priyal et al. (2019) have
generated a time series of Ca-K images after correcting the limb darkening and
instrumental intensity vignetting. The contrast of the images in the time
series was made uniform by applying the ECT methodology (Singh et al., 2021).
The Ca-K features, such as plages, EN, AN, and QN regions, were segregated
based on intensity and area threshold values in the different latitude belts
with an interval of 10∘ up to 90∘ (Priyal et al., 2023). In panel (a) of
Figure 1, we present an example of the analyzed Ca-K image taken on 1936 April
11. The identified plages, EN, AN, and QN, in the binary format, are presented
in panels (b), (c), (d), and (e), respectively. These digitized Ca-K data are
used from 1907 to 1996 to study the differential rotation. After 1996, there
were significant data gaps. Therefore, we excluded the later part (1997–2007)
of the data for the present study. Distribution of daily fractional area of
plage, EN, AN, and QN is shown in the four panels of Figure 2 for the period
1907–1996.

In the absence of Ca-K data, the plages, EN, AN, and QN areas are interpolated
using the IDL subroutine ‘INTERPOL.PRO’ to fill the data gaps and obtain the
continuous time series. We have analyzed the data using the power spectral
analysis technique. Raw power spectra are obtained through a fast Fourier
transform (FFT) for data string of 512 days (Singh & Prabhu, 1985). The linear
trend in the data is removed before performing the FFT. To form successive
data string, the data for 256 days are taken common to successive data
intervals. The total data length of 32759 days (1907 January 01 to 1996
October 01) thus gave us total 127 epochs. This is done to assure the
continuous and at short interval information about the chromospheric rotation
period. The power spectra for one epoch encompassing 512 days from 1951 March
10 to 1952 August 03 in the 10–20∘ latitude belt of the southern hemisphere is
shown in Figure 3. The power spectrum shows a prominent peak at $\sim$27 days,
corresponding to the solar rotation rate. However, power spectra show double
peaks at $\sim$27 days, suggesting the temporal variation of the rotation
period over a time scale shorter than the total data length. However, we have
selected the most dominant peak only. A similar power spectral analysis is
also performed for other epochs, a total of 127 epochs, and the rotation
period corresponding to the maximum power is noted.

##  3 Results and Discussions

We identify the features such as plages, EN, AN, and QN based on intensity
threshold and computed areas of each feature, separately, were grouped into
different latitude belts considering the 0∘–10∘, 10∘–20∘, 20∘–30∘, 30∘–40∘,
40∘–50∘, 50∘–60∘, 60∘–70∘, and 70∘–80∘ belts in both the north and south
hemispheres. Hereafter, we refer to them by the mean latitude as 5∘, 15∘, 25∘,
35∘, 45∘, 55∘, 65∘, and 75∘ latitude belts. The derived rotation periods are
noted as a function of time for various latitudes. In this way, we have
obtained 127 rotation period values for each latitude belt. We have computed
the mean and standard deviation of the rotation period for each latitude belt
from the pool of data points. In some cases, the power spectrum shows two
peaks of similar strength separated by a few days. This can lead to data
scatter, particularly at higher latitudes. To address this, we filtered the
data by removing values that fell outside one and half times standard
deviation. This approach helps concentrate the data points around a central
value and reduce the influence of outliers. Finally, we calculated the average
and standard deviation of these refined rotation periods to determine the
final values for our analysis. Afterward, the rotation period is converted
into rotation rate by using the relationship as, Rotation Rate (deg day-1)=
360∘/Period (days).

The average synodic rotation rate measured using plage areas is plotted as a
function of heliographic latitude in Figure 4. Plages are the magnetically
active chromospheric structures prominently visible in the Ca II K line
(3933.67 $\AA$) and mainly occur up to mid-latitudes. Therefore, its rotation
rate is measured only up to $\pm$55∘ latitudes (positive for northern
hemisphere and negative for southern hemisphere). The figure shows average
rotation rates (represented by dots) for different latitude bands. The solid
line represents the best-fit model capturing the overall trend. Error bars
indicate the standard error associated with each data point. This model
fitting was achieved using a fourth-order polynomial equation. The best-fit to
the rotation rate ($\omega$) values is a fourth order polynomial given by the
equation,

$\omega=A+B\phi+C\phi^{2}+D\phi^{3}+E\phi^{4}$ (1)

Where, $\phi$ is the heliographic latitude (in degrees). Coefficient A
represents the rotation rate at the equator. Here, the value of coefficients
A, B, C, D, and E are obtained to be 13.983$\pm$0.044, -0.0013$\pm$0.0017,
-0.0003$\pm$8.24$\times$10-5, 1.54$\times 10^{-6}\pm$7.56$\times$10-7, and
-1.38$\times 10^{-7}\pm$2.65$\times$10-8, respectively. Using this equation,
we have calculated the rotation rate at different latitudes as shown in Table
1. At the equator, the rotation rate is obtained to be 13.983$\pm$0.044 deg
day-1 for the plage area. As we move towards the higher polar latitudes, the
rotation rate values decrease accordingly in both the north and south
hemispheres. This trend of rotation rate compares well with previous studies
(Lamb, 2017; Bertello et al., 2020), but they have used different data sets.
Lamb (2017) used the photospheric magnetic features from high-cadence HMI
aboard the SDO, while Bertello et al. (2020) used full-disk Ca II K images
from Mount Wilson Observatory. A comparison of the results of rotation rate is
presented in Beck (2000), Bertello et al. (2020), and Jha et al. (2021).

We extended our study to the small-scale network areas. Figures 5 and 6
present the rotation rate values as a function of heliographic latitude
measured by using EN and AN, respectively. We note that EN and AN show a
similar variation in rotation rate as for the plage area in Figure 4. Here the
best-fit model for the rotation rate is a fourth order polynomial given by
equation 1. Dots indicate the rotation rate data point at each latitude belt.
For EN, the value of coefficients is obtained to be 13.912$\pm$0.054,
-0.0043$\pm$0.0015, -0.0008$\pm$5.50$\times$10-5, 5.65$\times
10^{-7}\pm$3.82$\times$10-7, and 3.42$\times 10^{-8}\pm$9.80$\times$10-9. For
AN, the coefficients A, B, C, D, and E are measured as 13.985$\pm$0.053,
-0.0029$\pm$0.0015, -0.0009$\pm$5.44$\times$10-5, 6.13$\times
10^{-7}\pm$3.78$\times$10-7, and 5.19$\times 10^{-8}\pm$9.70$\times$10-9,
respectively. The equatorial rotation rate is obtained to be 13.91$\pm$0.05
deg day-1 for EN and 13.99$\pm$0.05 deg day-1 for AN. These values are
comparable to the value obtained for the plage area.

The same analysis is performed on the QN areas. Figure 7 shows a plot between
the rotation rate and the heliographic latitude for the quiet network area. We
find that plage, EN, AN, and QN show a similar variation of rotation rate with
latitudes. Coefficients A, B, C, D, and E of the fourth-order polynomial fit
(given by equation 1) are measured to be 14.111$\pm$0.083, -0.0031$\pm$0.0024,
-0.0009$\pm$8.41$\times$10-5, 5.36$\times 10^{-7}\pm$5.85$\times$10-7, and
5.54$\times 10^{-8}\pm$1.49$\times$10-8, respectively. The rotation rate for
QN at the equator is 14.11$\pm$0.08 deg day-1 which again compares well with
the rotation rate value at the equator obtained using other chromospheric
features.

In Table 2, we compare our results with the previous studies carried out using
Doppler shift, magnetic features, sunspot area, and Ca K plage area at the
equator and at 40∘ latitude. The rotation rate value at 40∘ latitude is
calculated with the help of best-fit equation obtained in the respective
studies. We found that our results are comparable with the previous studies.
Komm et al. (1993) measured sidereal rotation rate of small magnetic features
(without active regions). They found the equatorial rotation rate as 2.913
$\mu$rad s-1 or 14.420 deg day-1 and at 40∘ latitude as 13.236 deg day-1.

Gupta et al. (1999) measured rotation rate for three groups of sunspot,
sunspot area less than 5 millionth of hemisphere, sunspot area larger than 5,
but less than 15 millionth of hemisphere, and sunspot area larger than 15
millionth of hemisphere and found that the smaller sunspots (with small area)
rotate faster than the bigger ones (with large area) (see Table 2 for rotation
rate values). A similar result was obtained by Jha et al. (2021). They have
classified sunspots area into two groups: smaller sunspots with sunspot area
less than 200 millionth of hemisphere and the bigger sunspots with sunspot
area larger than 400 millionth of hemisphere and then measured the equatorial
rotation rate as 14.399 deg day-1 for small area group and 14.351 deg day-1
for large area group. Therefore, they have concluded that the bigger sunspots
group rotates slower than the smaller ones. One plausible explanation for
these results is that sunspots of different sizes are rooted at different
depths (Gilman & Foukal, 1979; Beck, 2000). It is found that the large size
sunspots suppress the differential rotation in the photosphere (Brajša et al.,
2006; Li et al., 2020). The smaller rotation rates of larger sunspot thus may
indicate the deeper anchoring depth of the parent ﬂux tubes (Livingston et
al., 2006). A detailed and comprehensive analysis of the physical mechanisms
underlying the differential rotation rates of different sized sunspots is
presented by Wan et al. (2023b).

A plot for the comparison of rotation rates of plage area and small-scale
networks is presented in Figure 8. Rotation rate of plages are slightly higher
than that of small-scale networks at mid-latitudes, slightly lower at higher
latitudes and comparable to small-scale network near the equator (Table 1).
However, the difference in the rotation rate of two populations is small.
Therefore, we can infer that the small scale chromospheric features observed
in Ca-K line show a similar rotation rate as large structures such as plages.
This indicates that large and small-scale features have the foot points at the
same layer of the Sun. Hale (1908) obtained that the rotation rate of plages
are consistent with that of sunspots or surface magnetic field. The rotation
rates obtained by the quiet background network emission pattern in the Ca II K
line are in close agreement with those of small sunspots (Schroeter et al.,
1978). Belvedere et al. (1977) investigated that the large evolved
chromospheric plages rotate more rigidly than the smaller ones. However,
Antonucci et al. (1979b) found that the chromospheric features with sizes
24–300 $\times$ 103 km, such as network elements and active regions, exhibit
the same differential rotation. Recent studies by Li et al. (2020) and Li et
al. (2023) reported that the quiet chromosphere rotates at the same rate as
the small-scale magnetic elements do. In the chromosphere, in contrast to
photosphere, the rotation rate is strengthened by large-scale magnetic field
(Li et al., 2013; Wan et al., 2023a). Further, it has been suggested that
magnetic elements in quiet regions rotate marginally faster than sunspots in
active regions (Xiang et al., 2014; Xu & Gao, 2016). Therefore, the driving
force for the chromosphere to rotate in this way is expected from the solar
interior, not from the photosphere (Li et al., 2020).

##  4 Conclusions

In this study, the digitized data of four chromospheric features viz plage,
AN, EN, and QN area obtained from Kodaikanal Observatory for the period
1907–1996 are used to investigate the differential rotation at different
latitude belt from 0 to 80 degree with a step size of 10∘ in both hemispheres.
We find that plage areas as well as chromospheric networks are showing the
pattern of solar differential rotation. Even QNs (most fainter structures) are
found to show the differential rotation profile. From this, we concluded that
chromospheric features like plage, EN, AN, and even QN areas are reliable
parameters for the study of differential rotation of the solar chromosphere.
Moreover, this time, we have used the data of higher polar latitudes up to
$\pm$80∘ for the study of differential rotation. We observed that rotation
rate presents a decreasing trend from equator to poles regardless of solar
structure type used to measure the rotation rate. The equatorial rotation rate
for plage, EN, AN, and QN area are estimated to be 13.983$\pm$0.044 deg day-1,
13.912$\pm$0.054 deg day-1, 13.985$\pm$0.053 deg day-1, and 14.111$\pm$0.083
deg day-1, respectively. The results of this study will help to understand the
differential rotation and so the dynamo theory of the Sun. In this way,
further analyses can be done in future by using the data of other
observatories and by other methods; and a comparision of such studies will
definitely enhance our understanding toward the important phenomena of solar
differential rotation.

We thank the anonymous referee for providing valuable comments and suggestions
on the paper. We would like to thank all the observers at Kodaikanal
Observatory who made the observations, maintained the data at the observatory,
and digitized all the images. The data digitization was supervised by Jagdev
Singh and supported by F. Gabriel and P U. Kamath.

## References

  * Antonucci et al. (1979a) Antonucci, E., Azzarelli, L., Casalini, P., Cerri, S., & Denoth, F. 1979a, Sol. Phys., 61, 9, doi: 10.1007/BF00155441
  * Antonucci et al. (1979b) —. 1979b, Sol. Phys., 63, 17, doi: 10.1007/BF00155691
  * Babcock (1961) Babcock, H. W. 1961, ApJ, 133, 572, doi: 10.1086/147060
  * Beck (2000) Beck, J. G. 2000, Sol. Phys., 191, 47, doi: 10.1023/A:1005226402796
  * Belvedere et al. (1977) Belvedere, G., Godoli, G., Motta, S., Paterno, L., & Zappala, R. A. 1977, ApJ, 214, L91, doi: 10.1086/182449
  * Bertello et al. (2020) Bertello, L., Pevtsov, A. A., & Ulrich, R. K. 2020, ApJ, 897, 181, doi: 10.3847/1538-4357/ab9746
  * Brajša et al. (2006) Brajša, R., Ruždjak, D., & Wöhl, H. 2006, Sol. Phys., 237, 365, doi: 10.1007/s11207-006-0076-7
  * Chu et al. (2010) Chu, Z., Zhang, J., Nie, Q. X., & Li, T. 2010, Sol. Phys., 264, 1, doi: 10.1007/s11207-010-9561-0
  * Dikpati & Gilman (2006) Dikpati, M., & Gilman, P. A. 2006, ApJ, 649, 498, doi: 10.1086/506314
  * Gilman & Foukal (1979) Gilman, P. A., & Foukal, P. V. 1979, ApJ, 229, 1179, doi: 10.1086/157052
  * Gupta et al. (1999) Gupta, S. S., Sivaraman, K. R., & Howard, R. F. 1999, Sol. Phys., 188, 225, doi: 10.1023/A:1005229124554
  * Hale (1908) Hale, G. E. 1908, ApJ, 28, 315, doi: 10.1086/141602
  * Hathaway & Rightmire (2011) Hathaway, D. H., & Rightmire, L. 2011, ApJ, 729, 80, doi: 10.1088/0004-637X/729/2/80
  * Hathaway et al. (2022) Hathaway, D. H., Upton, L. A., & Mahajan, S. S. 2022, Frontiers in Astronomy and Space Sciences, 9, 419, doi: 10.3389/fspas.2022.1007290
  * Hotta & Kusano (2021) Hotta, H., & Kusano, K. 2021, Nature Astronomy, 5, 1100, doi: 10.1038/s41550-021-01459-0
  * Howard & Harvey (1970) Howard, R., & Harvey, J. 1970, Sol. Phys., 12, 23, doi: 10.1007/BF02276562
  * Imada & Fujiyama (2018) Imada, S., & Fujiyama, M. 2018, ApJ, 864, L5, doi: 10.3847/2041-8213/aad904
  * Insley et al. (1995) Insley, J. E., Moore, V., & Harrison, R. A. 1995, Sol. Phys., 160, 1, doi: 10.1007/BF00679089
  * Javaraiah (2020) Javaraiah, J. 2020, Sol. Phys., 295, 170, doi: 10.1007/s11207-020-01740-x
  * Javaraiah et al. (2005a) Javaraiah, J., Bertello, L., & Ulrich, R. K. 2005a, ApJ, 626, 579, doi: 10.1086/429898
  * Javaraiah et al. (2005b) —. 2005b, Sol. Phys., 232, 25, doi: 10.1007/s11207-005-8776-y
  * Jha et al. (2021) Jha, B. K., Priyadarshi, A., Mandal, S., Chatterjee, S., & Banerjee, D. 2021, Sol. Phys., 296, 25, doi: 10.1007/s11207-021-01767-8
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys., 143, 19, doi: 10.1007/BF00619094
  * Lamb (2017) Lamb, D. A. 2017, ApJ, 836, 10, doi: 10.3847/1538-4357/836/1/10
  * Li et al. (2023) Li, K. J., Wan, M., & Feng, W. 2023, MNRAS, 520, 5928, doi: 10.1093/mnras/stad478
  * Li et al. (2013) Li, K. J., Xie, J. L., & Shi, X. J. 2013, ApJS, 206, 15, doi: 10.1088/0067-0049/206/2/15
  * Li & Xu (2024) Li, K. J., & Xu, J. C. 2024, MNRAS, 528, 1438, doi: 10.1093/mnras/stae044
  * Li et al. (2020) Li, K. J., Xu, J. C., Xie, J. L., & Feng, W. 2020, ApJ, 905, L11, doi: 10.3847/2041-8213/abcb84
  * Livingston et al. (2006) Livingston, W., Harvey, J. W., Malanushenko, O. V., & Webster, L. 2006, Sol. Phys., 239, 41, doi: 10.1007/s11207-006-0265-4
  * Livingston (1969) Livingston, W. C. 1969, Sol. Phys., 9, 448, doi: 10.1007/BF02391670
  * Lustig (1983) Lustig, G. 1983, A&A, 125, 355
  * Mancuso et al. (2020) Mancuso, S., Giordano, S., Barghini, D., & Telloni, D. 2020, A&A, 644, A18, doi: 10.1051/0004-6361/202039094
  * Mishra et al. (2024) Mishra, D. K., Routh, S., Jha, B. K., et al. 2024, ApJ, 961, 40, doi: 10.3847/1538-4357/ad1188
  * Priyal et al. (2019) Priyal, M., Singh, J., Ravindra, B., & Shekar B, C. 2019, Sol. Phys., 294, 131, doi: 10.1007/s11207-019-1522-7
  * Priyal et al. (2023) Priyal, M., Singh, J., Ravindra, B., & Sindhuja, G. 2023, ApJ, 944, 218, doi: 10.3847/1538-4357/acaf60
  * Ruždjak et al. (2017) Ruždjak, D., Brajša, R., Sudar, D., Skokić, I., & Poljančić Beljan, I. 2017, Sol. Phys., 292, 179, doi: 10.1007/s11207-017-1199-8
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129, doi: 10.1007/BF00733429
  * Schroeter & Woehl (1975) Schroeter, E. H., & Woehl, H. 1975, Sol. Phys., 42, 3, doi: 10.1007/BF00153278
  * Schroeter et al. (1978) Schroeter, E. H., Woehl, H., Soltau, D., & Vazquez, M. 1978, Sol. Phys., 60, 181, doi: 10.1007/BF00152342
  * Sharma et al. (2020) Sharma, J., Kumar, B., Malik, A. K., & Vats, H. O. 2020, MNRAS, 492, 5391, doi: 10.1093/mnras/staa188
  * Singh & Prabhu (1985) Singh, J., & Prabhu, T. P. 1985, Sol. Phys., 97, 203, doi: 10.1007/BF00152989
  * Singh et al. (2021) Singh, J., Priyal, M., & Ravindra, B. 2021, ApJ, 908, 210, doi: 10.3847/1538-4357/abd021
  * Ternullo (1986) Ternullo, M. 1986, Sol. Phys., 105, 197, doi: 10.1007/BF00156387
  * Ternullo (1987) —. 1987, Sol. Phys., 112, 143, doi: 10.1007/BF00148494
  * Vats et al. (2001) Vats, H. O., Cecatto, J. R., Mehta, M., Sawant, H. S., & Neri, J. A. C. F. 2001, ApJ, 548, L87, doi: 10.1086/318924
  * Wan & Gao (2022) Wan, M., & Gao, P.-x. 2022, ApJ, 939, 111, doi: 10.3847/1538-4357/ac930d
  * Wan et al. (2023a) Wan, M., Gao, P. X., Xu, J. C., et al. 2023a, MNRAS, 520, 988, doi: 10.1093/mnras/stad192
  * Wan et al. (2023b) Wan, M., Gao, P.-X., Zeng, S.-G., & Deng, L.-h. 2023b, ApJ, 956, 53, doi: 10.3847/1538-4357/acf75a
  * Wan & Li (2022) Wan, M., & Li, K. 2022, Sol. Phys., 297, 126, doi: 10.1007/s11207-022-02059-5
  * Wöhl et al. (2010) Wöhl, H., Brajša, R., Hanslmeier, A., & Gissot, S. F. 2010, A&A, 520, A29, doi: 10.1051/0004-6361/200913081
  * Xiang et al. (2014) Xiang, N. B., Qu, Z. N., & Zhai, Q. 2014, AJ, 148, 12, doi: 10.1088/0004-6256/148/1/12
  * Xu & Gao (2016) Xu, J. C., & Gao, P. X. 2016, ApJ, 833, 144, doi: 10.3847/1538-4357/833/2/144
  * Xu et al. (2020) Xu, J.-C., Gao, P.-X., & Shi, X.-J. 2020, ApJ, 902, 64, doi: 10.3847/1538-4357/abb5b7
  * Zhang et al. (2013) Zhang, L., Mursula, K., & Usoskin, I. 2013, A&A, 552, A84, doi: 10.1051/0004-6361/201220693

Table 1: Differential rotation rate measured using plage, EN, AN, and QN area. Latitude |  | Rotation rate (Rotation period) |  |
---|---|---|---|---
(Degree) |  | deg day-1 (days) |  |
| Plage | EN | AN | QN
-80 | - | 10.146 (35.481) | 10.012 (35.953) | 10.211 (35.252)
-70 | - | 10.843 (33.201) | 10.618 (33.904) | 10.772 (33.418)
-60 | 10.534 (34.173) | 11.554 (31.155) | 11.315 (31.813) | 11.445 (31.453)
-50 | 12.018 (29.954) | 12.231 (29.432) | 12.028 (29.928) | 12.147 (29.636)
-40 | 12.960 (27.777) | 12.831 (28.056) | 12.691 (28.365) | 12.808 (28.106)
-30 | 13.518 (26.629) | 13.320 (27.026) | 13.252 (27.164) | 13.371 (26.922)
-20 | 13.819 (26.049) | 13.673 (26.327) | 13.671 (26.332) | 13.794 (26.096)
-10 | 13.955 (25.797) | 13.873 (25.947) | 13.920 (25.861) | 14.046 (25.628)
0 | 13.983 (25.744) | 13.912 (25.876) | 13.985 (25.741) | 14.111 (25.511)
10 | 13.931 (25.841) | 13.787 (26.110) | 13.862 (25.968) | 13.984 (25.742)
20 | 13.790 (26.105) | 13.508 (26.650) | 13.563 (26.542) | 13.676 (26.321)
30 | 13.520 (26.625) | 13.089 (27.503) | 13.109 (27.462) | 13.211 (27.249)
40 | 13.049 (27.588) | 12.554 (28.675) | 12.534 (28.720) | 12.624 (28.516)
50 | 12.268 (29.344) | 11.936 (30.160) | 11.887 (30.283) | 11.965 (30.086)
60 | 11.038 (32.613) | 11.275 (31.928) | 11.227 (32.063) | 11.298 (31.863)
70 | - | 10.619 (33.900) | 10.627 (33.875) | 10.698 (33.649)
80 | - | 10.026 (35.906) | 10.170 (35.397) | 10.255 (35.101)
Table 2: Comparision of rotation rate measured for different features. Author | Features | Rotation rate (deg day-1) |
---|---|---|---
|  | at equator | at 40∘ latitude
Howard & Harvey (1970) | Doppler Shift | 13.761 | 12.675
Komm et al. (1993) | small magnetic features | 14.420 | 13.236
Gupta et al. (1999) | Sunspot area$<$5$\mu$Hem | 14.491 | 13.313
Gupta et al. (1999) | 5$\mu$Hem$<$Sunspot area$<$15$\mu$Hem | 14.380 | 13.206
Gupta et al. (1999) | Sunspot area$>$15$\mu$Hem | 14.279 | 13.109
Lamb (2017) | magnetic feature | 14.296 | 13.086
Bertello et al. (2020) | Ca II K plage | 14.286 | 13.024
Jha et al. (2021) | Sunspot area | 14.381 | 13.257
Wan & Gao (2022) | Ca II K plage | 13.496 | 12.960
Present work | Plage | 13.983 | 12.960
Present work | EN | 13.912 | 12.831
Present work | AN | 13.985 | 12.691
Present work | QN | 14.111 | 12.808
Figure 1: Panel (a) of the figure shows the analyzed Ca-K image taken on 1936
April 11, 07:46:00 UT. Panels (b), (c), (d), and (e) show the identified
plages, EN, AN, and QN in the binary format, respectively. Figure 2: (Top to
bottom) Distribution of daily fractional area of plages, EN, AN, and QN for
the period 1907-1996 for 10-20∘ latitudinal belt of northern hemisphere.
Figure 3: A sample FFT power spectra of 512 days of data (from 1951 March 10
to 1952 August 03) is shown for a plage area of 10-20∘ latitude belt. Figure
4: Average synodic solar rotation rate measured using the plage area during
1907–1996 as a function of latitude. Error bars are standard error values.
Solid line represents the fourth order polynomial fit as
13.983-0.0013$\phi$-0.0003($\phi$)2+1.54$\times$10-6($\phi$)3-1.38$\times$10-7($\phi$)4,
where $\phi$ is the heliographic latitude. Figure 5: Same as Figure 4 but for
EN area. Solid line represents the fourth order polynomial fit as
13.912-0.0043$\phi$-0.0008($\phi$)2+5.65$\times$10-7($\phi$)3+3.42$\times$10-8($\phi$)4,
where $\phi$ is the heliographic latitude. Figure 6: Same as Figure 4 but for
AN area. Solid line represents the fourth order polynomial equation as
13.985-0.0029$\phi$-0.0009($\phi$)2+6.13$\times$10-7($\phi$)3+5.19$\times$10-8($\phi$)4,
where $\phi$ is the heliographic latitude. Figure 7: Same as Figure 4 but for
QN area. Solid line represents the fourth order polynomial equation as
14.111-0.0031$\phi$-0.0009($\phi$)2+5.36$\times$10-7($\phi$)3+5.54$\times$10-8($\phi$)4,
where $\phi$ is the heliographic latitude. Figure 8: Comparision of
differential rotation rate profiles obtained for plage area (blue dashed
line), EN (black dotted line), AN (red dash-dotted line), and QN (green solid
line).


††thanks: Third version

# Challenges in identifying simple pattern-forming mechanisms
in the development of settlements using demographic data

Bartosz Prokop  Corresponding author<EMAIL_ADDRESS>Laboratory
of Dynamics in Biological Systems,
Department of Cellular and Molecular Medicine,
KU Leuven, Leuven, Belgium  Peter F. Pelz  Chair of Fluid Systems
TU Darmstadt, Darmstadt, Germany  Lendert Gelens  Laboratory of Dynamics in
Biological Systems,
Department of Cellular and Molecular Medicine,
KU Leuven, Leuven, Belgium  John Friesen  Chair of Fluid Systems
TU Darmstadt, Darmstadt, Germany

###### Abstract

Abstract

The rapid increase of population and settlement structures in the Global South
during recent decades motivates the development of suitable models to describe
their formation and evolution. Such settlement formation has been previously
suggested to be dynamically driven by simple pattern-forming mechanisms. Here,
we explore the use of a data-driven white-box approach, called SINDy, to
discover differential equation models directly from available spatiotemporal
demographic data for three representative regions of the Global South. We show
that the current resolution and observation time of the available data is
insufficient to uncover relevant pattern-forming mechanisms in settlement
development. Using synthetic data generated with a generic pattern-forming
model, the Allen-Cahn equation, we characterize what the requirements are on
spatial and temporal resolution, as well as observation time, to successfully
identify possible model system equations. Overall, the study provides a
theoretical framework for the analysis of large-scale geographical/ecological
systems, and it motivates further improvements in optimization approaches and
data collection.

††preprint: APS/123-QED

##  I Introduction

The fraction of the population living in urban or settlement structures has
grown exponentially over the last several decades, especially in the Global
South [1, 2]. This trend poses one of the main challenges in our world [3] as
the rising population in such structures is in need of vital infrastructure
[4] while simultaneously affecting (mostly negative) climatic developments [5,
6]. Consequently, there is an urgent need to understand underlying processes
of urbanization and anticipate the emergence of these structures.

Urbanization and development of settlement structures depends on several
mechanisms based on repulsion and attraction [7, 8]. Such interactions can
lead to three major settlement distributions (see Fig. 1). Their existence has
been confirmed in recent settlement pattern studies of different regions in
the Global South, in which regular distributions are dominating [9, 10, 11,
12, 13].

Figure 1: Possible settlement arrangements following the phenomenological
interpretation as point processes, resulting in either clustered (i), random
(ii) or regular distributions (iii).

The existence of regularly-patterned distributions in settlements, and in
other spatial systems, can be an indicator for the existence of instability-
driven pattern-forming mechanisms [14]. A similar concept of linking spatial
distributions with specific driving mechanisms has been successfully applied
in the field of plant and animal ecology for decades [15, 16, 17]. To
understand the emergence of these patterns from underlying interactions,
different modeling approaches have been developed. For example, urban
development can be modeled by agent- or cellular-automata-based approaches
[18, 19] which include detailed interactions at the level of individuals.
Despite their accuracy in specific cases, such models have several drawbacks.
They lack generalization, require large and detailed datasets to be fit on and
can become computationally expensive. Another approach is to use reaction-
diffusion models. Such systems consisting of differential equations are
simpler and directly interpretable, yet they can also lead to highly complex
patterns and dynamics [20, 21]. In the context of urban structures, Pelz et
al. [22] have developed a theoretical framework describing the formation of
informal settlements (so-called slums) in the Global South. Furthermore, this
framework was extended to describe the morphogenesis of urban systems in the
United States as a reaction-diffusion system in [23].

Deriving such models is typically done by suggesting sets of possible models
from experimental measurements of specific interactions paired with scientific
intuition and determining which one provides the best fit to measured data.
However, as data structures of urban system are complex, suggesting certain
model structures can be difficult. One way to address this is by using
recently developed data-driven model discovery approaches, such as the system
identification approach called ”Sparse Identification of Nonlinear Dynamics”
(SINDy) [24]. SINDy has recently gained attention in many fields, such as
engineering [25], physics [26], chemistry [27] and biology [28]. The method
has shown success in identifying interpretable models in the form of ordinary
or partial (through the extension PDE-FIND [29]) differential equations from
synthetic data of e.g. pattern-forming mechanisms in the form of reaction-
diffusion equations [29, 30]. However, the literature on model discovery from
real temporal or spatiotemporal data with SINDy is scarce (e.g. on generic
benchmark problems in [31] or [32]), as the method struggles with commonly
encountered high-noise and low-data situations.

In this work, we will theoretically and practically investigate if pattern-
forming processes can be responsible for settlement development in the Global
South and attempt to identify such mechanisms using SINDy directly from
existing spatiotemporal data of population patterns. In section II, we
motivate the potential relevance of pattern-forming equations from data
analysis of satellite images. In section III, we apply the SINDy algorithm to
satellite images and try to identify mathematical equations in the form of
one-component pattern-forming mechanisms. As model identification from the
real data sets turns out to be challenging, we identify and study two main
challenges, namely data availability and quality and determine their influence
on model discovery from spatiotemporal data in section IV. Lastly, we discuss
our findings and elaborate on these challenges and what is required to
overcome them in section V.

##  II Settlement development as a pattern-forming process

We start our study by selecting three representative regions of emerging
countries – the Punjab region in India, the Nile delta in Egypt and the Kano
State region in Nigeria (see Table 1). The regions are chosen as they lie in
countries which can be considered representative of the Global South: all
countries have had steady population as well as steady economic growth over
the last 20 years [1]. Despite and because of the great cultural differences,
all three societies are in transition from agricultural countries to
industrialized nations.

Table 1:  Investigated regions with their geopolitical location and attributes. Region | Attributes
---|---
Punjab, India | northwest India, border region with
| Pakistan, agricultural region called
| Granary of India
Nile delta, Egypt | north Egypt, densely populated, fast
| growing agricultural region
Kano State, Nigeria | north Nigeria, border region to Niger,
| agricultural region in one of the fastest
| growing economies

Furthermore, the spatial distribution of settlements has been studied in these
regions and their regularity characterized on different spatial scales [10,
12, 13]. This allows us to select subregions (here called regions of interest
(ROI)) which show a regular distribution. This is illustrated for the Punjab
region in India in Fig. 2i). For a detailed explanation of how we used
satellite data, in combination with the Global Artificial Impervious Area
(GAIA, temporal resolution $\Delta t=1\text{ year}^{-1}$, spatial resolution
$\Delta x=\Delta y\approx 30$ m) [33, 34] data set, to select the ROIs, we
refer to Appendix D. For each ROI, we then track the settlement evolution over
a period of about 15 years using the data set WorldPop [35] depicting spatial
population distributions, see Fig. 2ii. The WorldPop data set is used later
for the model identification process as described in Fig. 2iii 111The
computation of the ANN, the resulting selection of the ROIs, and the later
computation of the feature lengths was done with GAIA because it provides
discrete spatial data (1 - built-up, 0 - no built-up area). As such, and in
contrast to WorldPop, which distributes population quantities across
administrative regions, leading to non-existent values of 0, the
identification of settlement locations and sizes is straightforward. However,
for model identification WorldPop is used, since most pattern forming
mechanisms describe concentration distributions of different (chemical or
other) species and do not create discrete patterns [37, 39], similar to
WorldPop..

Figure 2: Overview of the workflow of this study i) Selection of ROIs,
starting from the identification of suitable regions of the Global South, here
Punjab, India, selection of an excerpt of the region in the data set GAIA and
determination of ROIs from regular settlement patterns with the ANN (see
Appendix D). ii) Data of population density patterns (on logarithmic scale)
from the data set WorldPop from ROI 2 of Punjab, India, for four time points
with a spatial resolution of $\Delta x,\Delta y\approx 100$m (enlarged in
figure) at the equator and temporal resolution of $\Delta t=1\text{
year}^{-1}$ iii) Work flow to investigate possible pattern-forming mechanisms
in the development of settlement structures. Starting from the time series of
population density patterns, we extract spatiotemporal features and apply
SINDy to discover models in the form of one-component partial differential
equations. This can help us understand the role of such mechanisms in
settlement structures and estimate future urban scenarios.

Previous work has shown that the observed emergence of rural settlement
structures could be caused by simple reaction-diffusion pattern-forming
mechanisms [22]. Following the example of Pelz et al., we divide the system
into a population living in the rural settlements, and a supply potential of
agricultural land which is complementary to the population. In other words,
areas that are in agricultural use are uninhabited, and vice-versa. We also
assume that there is a limit to agricultural exploitability and urban
densification. The population density at a spatial point $\bm{x}=(x,y)$ and
time $t$ is given by $u(\bm{x},t)$ and the corresponding supply potential by
$v(\bm{x},t)$.

The change of the concentrations of $u$ and $v$ is defined by three
complementary global contributions: (i) birth or death of the population or
the cultivation or sealing of agricultural space within a domain of size $A$,
(ii) migration to and from other cities outside of $A$ leading to a de- or
increase in supply potential and (iii) migration to areas with higher supply
potential over the boundary $C$ of $A$ where a settlement exists or is
created. These mechanisms lead to the formulation of a typical reaction-
diffusion model (for more detail, see Appendix A and [22]):

$\displaystyle\begin{split}&u_{t}=\nabla^{2}u+Rf(u,v)\\\
&v_{t}=D\nabla^{2}v+Rg(u,v),\\\
&\text{with}\quad(\bm{n}\cdot\nabla)\begin{pmatrix}u\\\ v\end{pmatrix}=0\text{
on }C,\end{split}$ (1)

with $D,R$ being the diffusion and reaction coefficients alongside the
respective reaction terms $f,g$, which are being evaluated under no-flux
boundary conditions.

The linear stability analysis of the homogeneous system (no diffusion)
provides the Jacobian matrix $\bm{J}$. It is well known that specific sign-
combinations of components of the Jacobian allow for a initial homogeneous
distribution of population and supply potential to be stable – either through
substrate inhibition or as an activator-inhibitor system [37]. In the case of
settlement structures, the only physically sensible choice is substrate
inhibition, as the resulting concentration patterns are out of phase [22, 37]:

Figure 3: (i) Example of multiple effects found in spatiotemporal data sets of
population densities indicating the existence of suggested reaction-diffusion
equations (population density data from data set WorldPop [35], ROI 2, Punjab,
India). The observable behaviors are internal migration into existing
settlements (A), invasion of not occupied agricultural space (B) or local
migration triggered by competition over available agricultural space (C). (ii)
Calculated characteristic length of settlement patterns resulting from the
mean and the median heavy sided feature size distribution (shown in inlet, ROI
2, Punjab, India in 2018). The characteristic lengths obey the power law when
growing, with $\ell(t)\sim t^{1/7}$ and $\ell(t)\sim t^{1/3}$ respectively
(with R2 scores). iii) Calculated characteristic length of the pattern
(example of a feature size shown in inlet) resulting from the Allen-Cahn
equation. Here, both the mean and median characteristic length obey the power
law with $\ell(t)\sim t^{1/2}$ when growing (with R2 scores).

$\bm{J}=\begin{pmatrix}f_{u}\big{|}_{0}&f_{v}\big{|}_{0}\\\
g_{u}\big{|}_{0}&g_{v}\big{|}_{0}\end{pmatrix}=\begin{pmatrix}+&+\\\
-&-\end{pmatrix}.$ (2)

$f_{u}\big{|}_{0}>0,$ people attract other people

The population $u$ increases due to self-reproduction of $u$ in an environment
of sufficient sustenance.

$f_{v}\big{|}_{0}>0,$ supply attracts people

The amount of available agricultural area $v$ attracts people, increasing the
population $u$ of the settlement.

$g_{u}\big{|}_{0}<0,$ people inhibit supply

The higher the population $u$, the less agricultural area $v$ is available,
especially due to the limitation by the agricultural needs of surrounding
settlements. This is the case until the maximum amount of agricultural area is
used and does not suffice to supply the population of a settlement, leading to
a decrease of population $u$.

$g_{v}\big{|}_{0}<0,$ supply inhibits additional supply

Due to limited resources, the supply production decreases, when agricultural
efficiency plateaus.

In the presence of diffusion, the linear stability shows that an equally
dispersed population and supply potential densities can destabilize to form
spatial patterns when the following condition is met (see Appendix B):

$f_{u}\big{|}_{0}>-\frac{g_{v}\big{|}_{0}}{D}.$ (3)

This is the case when the attraction of people to $A$ dominates the inhibition
of supply potential due to the emergence of settlements. If the domain size
$A$ and diffusion coefficient $D$ of the system are sufficiently large,
different settlements can emerge that constantly compete against each other
over the supply potential, eventually leading to a regular distribution of
settlements.

Assuming that the total population density $u$ and supply potential $v$ are
conserved, $u+v=c_{max}$, the system reduces to a simpler one-component
equation [38, 39]:

$\displaystyle\begin{split}u_{t}(x,y,t)&=\check{D}\nabla^{2}u(x,y,t)+R\check{f}(u(x,y,t)),\\\
\text{with}\quad\check{D}&=\frac{D+1}{2},\\\ \text{and
}\quad\check{f}&=\frac{1}{2}\left[f(u(x,y,t))+g(u(x,y,t))\right].\end{split}$
(4)

Based on this proposed reaction-diffusion description consistent with pattern
formation [22], we first set out to see whether the satellite data (Fig. 2ii)
contained any clear signatures of such developing rural settlement patterns.
Visual inspection of the WorldPop data set indicates the possible presence of
three key processes: local migration through attraction of bigger settlements
(e.g. between 2003 and 2008, Fig. 3i/A), invasion or occupation of available
agricultural space (e.g. between 2008 and 2013, Fig. 3i/B), and local
migration induced by competition of settlements over available space leading
to changes in settlement patterns (e.g. between 2013 and 2018, Fig. 3i/C).

We then analyzed the characteristic lengths $\ell(t)$ of settlements (see Fig.
3ii) using the GAIA data set [33, 34]. We chose the characteristic length to
be the size of the respective features, evaluating the distribution feature
sizes and determining the predominant size following [40] and [41] (see inset
in Fig. 3ii, for more details see Appendix E). This analysis shows that when
the characteristic size of settlements is growing, this growth is well
described by a power law $ct^{\beta}$, with an exponent of $\beta=1/7$ for the
mean and $\beta=1/3$ median characteristic lengths.

Interestingly, from literature we know that the time evolution of
characteristic lengths of patterns driven by coarsening mechanisms also follow
such a power law. For example, for the Allen-Cahn equation, it was
theoretically shown that the change of $\ell(t)$ is described by the power law
with $\beta=1/2$ [40, 42, 41], while for the Cahn-Hilliard equation [38],
$\beta=1/3$ [40, 41]. Indeed, we calculated the characteristic length for a
simulated pattern for the Allen-Cahn equation of the form,

$\displaystyle\begin{split}u_{t}(x,y,t)=\>&\alpha\nabla^{2}u+\beta u+\gamma
u^{2}-u^{3},\end{split}$ (5)

using the same method as for the settlement patterns, which revealed that the
mean and median characteristic length $\ell(t)$ follow the power law with
$\beta=1/2$ (see Fig. 3iii). Even though the exponent $\beta$ of the power law
for settlement patterns is not the same as for the Allen-Cahn or Cahn-Hilliard
equation, it indicates that settlement patterns could be a product of a
simple, coarsening pattern-forming mechanism in the proposed form of Eq. (4).

##  III Model identification from settlement data

###  III.1 Model identification using SINDy

The SINDy method to identify differential equation models from spatiotemporal
datasets has been increasingly applied in many fields, i.e. in fluid mechanics
[24]. However, it has, to our knowledge, never been used for studying large
scale, geo-sociological questions. Therefore, we asked ourselves whether we
could use this method to discover a partial differential equation (PDE) that
provides a good description of the measured time evolution of rural
settlements in the Global South. If successful, such a PDE would also provide
new insights in potential pattern-forming mechanisms in settlement
development.

The main idea behind the SINDy method is the assumption that dynamic systems
can be described through either ordinary, or in our case partial differential
equations (using PDE-FIND) [29], with sparse structure in the following form:

$u_{t}=N(u(\bm{x},t),u_{x},u_{y},...,\bm{x},\bm{\xi})$ (6)

The temporal change of $u$, $u_{t}$, is a function of the variable $u$ itself,
its spatial derivatives and a set of coefficients $\bm{\xi}$. Differential
equations of this form can be linearly combined:

$u_{t}=\xi_{1}+\xi_{2}u+\xi_{3}u^{2}+\xi_{4}u_{x}+\xi_{5}u_{xx}+...$ (7)

This equation can be rewritten as a row vector containing all combinations and
derivatives of the quantity, called the term library and a coefficient vector
$\bm{\xi}$ containing all coefficients:

$u_{t}=\begin{pmatrix}1&u&u^{2}&u_{x}&u_{xx}&...\end{pmatrix}\cdot\bm{\xi}.$
(8)

The values of each term in the library can be calculated from a single shot at
a given point in time. If this system is extended to all available time
points, a linear system of equations with the unknown parameter vector
$\bm{\xi}$ and the term library matrix $\bm{\Theta}$ is formed:

$\displaystyle\begin{split}\begin{pmatrix}\\\ \\\ \bm{u}_{t}\\\ \\\ \\\
\end{pmatrix}&=\begin{pmatrix}&&&&&\\\ &&&&&\\\ 1&u&u^{2}&u_{x}&u_{xx}&...\\\
&&&&&\\\ &&&&&\end{pmatrix}\cdot\begin{pmatrix}\\\ \bm{\xi}\\\ \\\
\end{pmatrix}=\bm{\Theta}\cdot\bm{\xi}\end{split}$ (9)

We assume that the settlement evolution could be captured by a PDE similar to
Eqs. (4) and (5). Therefore, we use the term library given in Tab. 2, which
includes derivatives up to the fourth order.

Table 2: Terms included in the library for the one-component, two-dimensional equation sorted by combinations of $u$ and its derivatives. |  Terms
---|---
Combinations |  $1$, $u$, $u^{2}$,$u^{3}$
Derivatives |  $u_{x}$,$u_{y}$, $u_{xx}$, $u_{yy}$, $u_{xy}$, $u_{xxx}$, $u_{yyy}$, $u_{xxy}$, $u_{yyx}$,
|  $u_{xxxx}$, $u_{yyyy}$, $u_{xxyy}$, $u_{xxxy}$, $u_{yyyx}$

This system poses an over-determined optimization problem for values of
$\bm{\xi}$ and can be solved using regression algorithms (for more detailed
information on regression algorithms for SINDy, see [43]). In contrast to the
original work in which the method PDE-FIND was introduced [29], we apply a
sparsity-promoting algorithm with the SR3 algorithm developed by [44]. This
method includes the additional variable $\bm{w}$, which is forced to be close
to the coefficient parameter and therefore relaxes the optimization problem.

$\displaystyle\begin{split}\min_{\bm{\xi},\bm{w}}&\frac{1}{2}\|\bm{u}_{t}-\bm{\Theta}\bm{\xi}\|^{2}_{2}+\lambda\|\bm{w}\|_{1}+\frac{\alpha}{2}\|\bm{w}-\bm{\xi}\|^{2}_{2}\\\
&\text{with}\quad\lambda=\frac{l^{2}}{2\alpha}\end{split}$ (10)

Here, two hyper-parameters of the optimization have to be set: the threshold
$l$ and the parameter of the optimization $\alpha$, which provides the
penalizing parameter $\lambda$ of the regularization.

After model identification, we analyze the discovered models with the Akaike
Information Criterion (AIC). The AIC is a measure of parsimony [45]. It
compares the goodness of fit of a given model to other proposed models and
weighs it with the model’s complexity aiming on maximizing the information
provided by the simplest-as-possible model. For our analysis, we apply the
corrected formulation for finite sample sizes of the AIC ($\text{AIC}_{c}$)
proposed by [46], resulting from [47]:

$\text{AIC}_{c}=\text{AIC}+\frac{2(k+1)(k+2)}{m-k-2}.$ (11)

The AIC is described by the likelihood function of average error over time and
space $\epsilon$ as follows:

$\displaystyle\begin{split}\text{AIC}&=m\ln{(\epsilon/m)}+2k\\\ \text{with
}\quad\epsilon&=\frac{\left|\sum_{i=1}^{m_{t}}y_{i}-N(x_{i},\bm{\xi})\right|}{m}\\\
\text{and}\quad m&=m_{s}n_{\text{ROI}}\end{split}$ (12)

Figure 4: Parameter sweep of all ROIs and regions for sets of threshold $l$
and optimization parameter $\alpha$. (i) The analysis of the ROIs shows
different combinations of parameters where the AIC falls beneath 2 for the
respective ROI. As we search for a regionally valid solution, we select only
unique solutions at combinations of $l,\alpha$ where the $\text{AIC}<2$ for
all ROIs overlap. The selected unique equations and their respective
optimization parameter sets where they were found first are shown. Red
bordered markers depict the best identified models of each ROI. We see that
the AIC provides the best equations with the lowest error and complexity
(except for ROI 2, from which we show later through analysis of contributions
that only the production terms are significant) (ii) The respective error and
complexity of the found equations are depicted. The red bordered markers show
the best model for each ROI in each region following the AIC. (iii) The
identified coefficients $\bm{\xi}$ (as in Eq. (9)) are shown for each
equation. All identified models contain most of the derivatives (diffusive
terms) with large coefficient values. If found the coefficient values of the
reaction terms are multiple magnitudes smaller than of the diffusive terms.

The AIC depends on the number of observations $m$ (size of region $m_{s}=XY$
and amount of included ROIs $n_{\text{ROI}}=3$, where we interpret an
observation as the time series at every spatial point) and the number of terms
($k$) describing the complexity of an identified model.

Using this approach, we determine the most parsimonious model among all
potential models and study its properties. In order to compare the identified
models, we further normalize the AIC by the minimal value of the respective
analysis $\text{AIC}_{min}$. Here, following [46] and [47], a model that has
an $\text{AIC}_{c}-\text{AIC}_{min}<2$ has strong support for being the
correct underlying system, while the ones with
$\text{AIC}_{c}-\text{AIC}_{min}<8$ have weak support. Hereafter, we always
refer to the corrected $\text{AIC}_{c}$ when the AIC is mentioned.

###  III.2 Application to settlement data

Using the outlined approach, we look for potential models to describe the
available settlement data from the WorldPop data set in the different regions
(India, Egypt, Nigeria). We scan sets of thresholds from $l=10^{-6}$ to
$l=10^{2}$ and optimization hyper-parameter $\alpha=10^{-3}$ to
$\alpha=10^{3}$. The SR3 algorithm was applied with a tolerance of $10^{-2}$
and using 200 iterations (for details see [44]).

From the parameter scan, we determine for which combinations of $l$ and
$\alpha$ the identified models provide an $\text{AIC}<2$. In order to pre-
select possible solutions, we only evaluate equations with this AIC at
overlapping parameter pairs of $l$ and $\alpha$. The reason for this is, that
if the dynamics in the whole region follow the same rules (or dynamical
behavior), all best identified models should have the same mechanistic form.
Therefore, in Fig. 4i, we show for which sets of ($l,\alpha$) the low-AIC
regions of each ROI overlap (in yellow).

Next, we selected unique model equations at their lowest parameter values for
$l$ and $\alpha$, respectively. The parameter combinations for unique
solutions of each ROIs are shown in Fig. 4i. These unique solutions are
further compared with the use of the AIC, where we also depict the error and
the complexity of the found equations in Fig. 4ii. The optimal selected
equations are highlighted by markers with red borders and the values of
coefficients are shown in Fig. 4iii.

Figure 5: (i) The simulation of found models shows that the concentration in
the whole region increases with time, but no patterns form. Despite being
trained on the data set, the models are not able to reproduce the training
data (Here ROI 2, India). (ii) Calculated contributions to the change of $u$.
Even though the magnitude of coefficients is in the order of $10^{6}$ for
derivatives, their contribution is negligible and lies between order of
$10^{-17}$ and $10^{-21}$.

For the identified models, we analyze if they indeed reproduce the same
spatiotemporal dynamics and patterns as present in the original data they were
trained on. Fig. 5i shows a representative time evolution of the original and
simulated data for ROI 2 in India. One can see that while the spatial patterns
become more pronounced in the original data, this is not the case in the
simulated model. This is also seen by explicitly plotting the difference
between the original data and simulation. These observations are general for
all analyzed ROIs (full analysis in the shared repository).

To better understand why the identified model equations do not accurately
capture the original data, we look more closely at the coefficients of each
model term, as well as the overall contribution of each term. We find that the
assigned coefficients vary over multiple orders of magnitude, ranging from
$10^{-3}$ for some production terms up to $10^{6}$ for some higher order
derivatives of the diffusive terms. We calculated the actual contribution
$c_{j}$ of each term as the value of each term $\theta_{ij}$ at a set time
point $t$ (we arbitrarily chose $t=10$) and averaged over all spatial points,
and then multiplied with the respective coefficient $\xi_{j}$,

$c_{j}=\left|\xi_{j}\cdot\bar{\theta_{j}}\right|\quad\text{with}\quad\bar{\theta_{j}}=\frac{\sum_{x=0}^{X}\sum_{y=0}^{Y}\theta_{i=10\text{
}j}(x,y)}{XY}$ (13)

Fig. 5ii shows that even though the coefficients of derivative terms were
multiple orders of magnitude larger than those of the production terms, their
contribution is in fact negligible compared to production terms. The dominance
of the terms proportional to $(1,u,u^{2},u^{3})$ over all terms with spatial
derivatives leads to the observed overall increase in concentration while
preserving the initial pattern. In conclusion, the identified models do not
capture any pattern-forming mechanism as the interaction of reaction and
diffusion processes is crucial (see Section II).

##  IV Low-data limits in model identification of pattern-forming processes

One possible reason is that the SINDy method was unable to recover a reaction-
diffusion model that correctly describes the settlement data due to a lack of
spatiotemporal resolution and/or insufficient observation time. Indeed, it is
known that model discovery with SINDy is dependent on the amount of temporal
points and the size of the time step $\Delta t$ [48, 44]. Therefore, we
decided to study the limits of the SINDy method in recovering a reaction-
diffusion model for low spatial and/or temporal resolution, as well as short
observation times. As we found that the observed settlements followed
coarsening dynamics (see section II), we decided to study the recovery of the
Allen-Cahn (AC) equation (5, [49, 39]) using SINDy in the low-data limit.

$\displaystyle\begin{split}u_{t}(x,y,t)=\>&\alpha\nabla^{2}u+\beta u+\gamma
u^{2}-u^{3}\\\ \text{with}\quad&\alpha=0.1,\beta=0.5,\gamma=-0.01.\end{split}$
(14)

With this set of coefficients, the initial condition
$u_{init}(x,y,t=0)\sim\mathcal{N}(0,0.01)$ (which was created once and used
for all simulations), and zero-flux boundary conditions, Eq. (14) creates
coarsening labyrinth patterns as shown in Fig. 6i (Original).

Figure 6: Recovery of the AC equation with SINDy: sensitivity to
spatiotemporal resolution. (i) Time simulation of the AC equation, sub-sampled
with different spatial and temporal resolutions, as indicated in (ii). (ii)
Diagram showing for which resolution SINDy is able to recover the AC equation.
Recovery is only successful for sufficiently high spatial and temporal
resolution ($214\times 214$, $N_{t}>80$ frames) of the simulated data set of
$t=100$.

We then first sub-sampled this dataset generated by simulating the AC equation
by imposing a spatial resolution of $\Delta x=\Delta y=0.39$
($N_{x}=N_{y}=N=256$ lattice points with $\Delta x=L/(N-1)$ and $L=100$) and a
temporal resolution of $\Delta t=0.08$ ($N_{t}=1250$ frames with
$N_{t}=t/\Delta t$ and $t=100$). This sub-sampled dataset was then used as
input to the SINDy algorithm, which was able to recover the original AC
equation (14) in the correct form with a maximum error of 3% in the
coefficients.

Figure 7: Recovery of the AC equation with SINDy: sensitivity to observation
time. (i) Time series of the simulated AC equation with temporal resolution of
$\Delta t=0.16$ and observation time $t=100$. (ii)-(iii) Diagrams showing for
which resolution SINDy is able to recover the AC equation with changing
observation time and spatial resolution. The observed time window used for
model recovery is chosen at the end (ii) or start (iii) of the time series in
(i). (iv) For a spatial resolution of $N=256$ lattice points varying the
temporal resolution shows that identification is less sensitive to resolution
then to the amount of available data (for case (iii), the transition points
are plotted with offset for visualisation).

Next, we further sub-sampled the data temporally and spatially in order to
identify the limits of the SINDy approach. We decreased the temporal
resolution from $\Delta t=0.08$ ($N_{t}=1250$ frames) in 26 steps to $\Delta
t=10$ ($N_{t}=10$ frames). Similarly, the spatial resolution was reduced from
$\Delta x=\Delta y=0.39$ ($N=256$ lattice points) in 18 steps to $\Delta
x=\Delta y=7$ ($N=14$ lattice points). We then compared the identified model
to the original AC model, characterizing if the form of the model (correct
terms) was correct, and whether it had identified the correct coefficients,
see Fig. 6. This analysis shows that the SINDy algorithm is sensitive to
spatial and temporal resolution. Only for sufficiently high spatial ($\Delta
x=\Delta y<0.46$, $N>214$ lattice points) and temporal resolution ($\Delta
t<1.25$, $N_{t}>80$ frames) we correctly recover the AC equation with the
proper coefficients (yellow region in Fig. 6ii). For lower resolution, we are
still able to recover the mechanistic form of the equation (but not the
correct coefficients) in the green region in Fig. 6ii. In Fig. 6i we show how
the detected mechanistic forms in this region are still able to recover the
overall dynamics of the system (for low-spatial/high-temporal and high-
spatial/low-temporal resolutions). However, the reduced resolution leads to
coefficients smaller than those in the original AC equation, which slows down
the dynamics. When decreasing the resolution even further, the optimization
algorithm includes and overestimates higher order derivatives. Such models are
no longer able to capture the dynamics of a system accurately (blue region in
Fig. 6ii; Fig. 6i low-spatial/low-temporal resolution).

We then wondered how the observation time affects the possibility to recover
the correct model. This question is especially relevant as we only observe
relatively slow changes of settlement structures over 20 years using the
WorldPop data set. We repeated the study of the recovery of the AC equation
(Fig. 6) for a fixed time resolution ($\Delta t=0.16$, $N_{t}=625$ frames) for
which we could successfully identify the AC equation (provided the spatial
resolution was sufficiently high). We then redid this analysis for varying
spatial resolution and observation time. The observation time was
progressively reduced by only considering the first or last fraction of the
original time series as input to the SINDy optimization (see sketch in Fig.
7i). Note that in this way we are also reducing the total amount of input
data. The results of this analysis are shown in Fig. 7ii-iii.

As expected, the AC equation can no longer be recovered when reducing the
observation time and/or spatial resolution below a certain threshold.
Interestingly, we find that the system is very sensitive to observation time
when using the later stages of the dynamical evolution (Fig. 7ii), while this
is much less the case when using the initial part of the dynamical evolution
of the AC equation (Fig. 7iii). This illustrates that the observation time
required to correctly recover the underlying model equation strongly depends
on when one measures the system dynamics. In particular, our analysis shows
that it is best to capture as much as possible of the dynamical changes at the
relevant time scales. In this case, much of the initial patterns form quickly
at the start, while later the patterns coarsen only slowly.

Finally, we then also investigated how sensitive these findings related to
observation time were with respect to time resolution. We fixed the spatial
resolution at $\Delta x=\Delta y=0.39$ or $N=256$ lattice points and
determined the critical thresholds in terms of observation time for correct
model identification for varying time resolution $\Delta t$. Fig. 7iv shows
that successful model recovery is less sensitive to time resolution (can be
varied over 2 orders of magnitude) than to observation time (both duration and
exact timing).

This also shows that there is a limit to the additional information that can
be provided by higher spatial and temporal resolution if the duration and/or
timing of the observed time window is not well chosen for successful model
recovery. In the case of the settlement data set under study here, this
analysis suggests that it is plausible that the observed changes in population
density in the provided data set are inadequate for proper model
identification due to the too short observation time compared to the relevant
time scales over which settlements develop.

##  V Discussion and Conclusion

The goal of this work was to not only theoretically describe the possible role
of simple pattern-forming mechanisms in the development of urban structures
(in our case settlements) of the Global South, as has been done before by [22,
23], but to provide an unbiased approach to identify such models directly from
data.

In order to do this, we selected three representative regions of the Global
South from India, Egypt and Nigeria and analyzed the occurrence of regularity
in settlement structures in these. Using this data, we selected smaller
regions of interest (ROIs) in the spatiotemporal data set WorldPop [35] of
population density distributions.

Following this, we extended the ideas of [22], motivating and providing a new
theoretical point of view on pattern-forming mechanisms in rural,
agriculturally dominated settlement structures. We argued that together with
features of regularity (as suggested by [12, 17, 14]), pattern-forming
mechanisms could be responsible for the emergence of settlement structures.
This extension can be a starting point in critically evaluating urban modeling
approaches that strive for more complexity over generalization. Here, we
substantiated this claim as we observed the suggested behavior in actual
population density patterns, while also showing that the characteristic length
of patterns resulting from settlements follows a power law, similarly to
coarsening patterns in e.g. Allen-Cahn or Cahn-Hilliard models.

We then introduced the SINDy [24, 29] method together with the AIC [46, 45],
allowing us to derive and investigate spatiotemporal models for the dynamics
of population density patterns. However, using the SINDy method, we were not
able to identify simple pattern-forming mechanisms directly from selected ROIs
of regions in the Global South. The found equations were neither sparse nor
represented known pattern-forming mechanisms from literature. The assigned
coefficients differed in multiple orders of magnitude between production terms
($10^{-3}$) and diffusion terms ($10^{6}$).

As a result, it seems necessary to change the target of optimization from
solely evaluating the coefficients to targeting the actual contribution of
terms, which e.g. has been recently suggested by [50]. Additionally, the
configuration of our SINDy approach does not include any time or space
dependency of parameters, as suggested by [51], which can prevent us from
capturing important dynamical behavior of settlement systems in the Global
South. The found models show that the models are not able to recreate training
data. An analysis of term contributions has revealed that the model dynamics
are dominated by the production terms and cannot be understood as pattern-
forming mechanisms.

Following this unsuccessful application, we identified and studied challenges
of model identification in spatiotemporal data sets considering the quality
and availability of data. We suggested that the used data set had too low
spatial and temporal resolution or an insufficient observation time
Subsequently we studied this question with a sensitivity study of SINDy
towards low-data limits. We show that SINDy (here PDE-FIND) is sensitive to
spatial and temporal resolution while identifying that observation time and
the observed dynamics have a significant influence on the recovery of
underlying dynamics. We see that when fast dynamics of pattern formation are
captured lower spatial and temporal resolution and a shorter observation time
are required to correctly identify the model. When only observing slow
dynamics (as seen in our settlement data) model identification is more
sensitive towards limited spatial resolution and requires longer observations
times in order to correctly identify the AC equation. Hence, we need to
closely follow the rapid improvements in data acquisition with satellite
imagery, which would provide us with sufficiently good data in resolution,
while simultaneously and more importantly provide us with longer observation
times. If such challenges are overcome this work provides a ready-to-use
framework to discover pattern-forming mechanisms in settlement development.

Moreover, the structure of the data should be adjusted. Currently the WorldPop
data set does not allow for uninhabited areas with a population density of 0.
Here, WorldPop itself is developing an improved data set, where population
densities and built-up areas are mapped. At the moment of publication this
data set only contains a single time point but, when extended, it will provide
new opportunities to study our question. Furthermore, the available data sets
limit us to only a single observed variable introducing a strong assumption
when considering models. Here, the application of methods coming from Koopman
theory, e.g. delay embedding [52, 43], could pose an interesting line of work
that can provide ’hidden’ variables, adding information for the optimization
and resulting in model of higher dimensions (as has been recently attempted
for synthetic data from a shallow-water model [53] or spatiotemporal Lotka-
Volterra model [54]).

In conclusion, we have provided an initial framework for the evaluation and
identification of the role of simple pattern-forming mechanisms in the
development of settlement structures. So far, the efforts were unfruitful to
provide model equations describing such mechanisms. However, we developed a
possible theoretical motivation and were able to identify the major challenges
of model identification in low-data limits in pattern-formation.

## Data and Code Availability

All calculations, simulations and graphs are done in Python. For SINDy, we use
the package PySINDy [55, 56] and for simulations, we developed a simple
forward-Euler solver. All algorithms are available in our Gitlab repository
(https://gitlab.kuleuven.be/gelenslab/publications/settlements.git).
Furthermore, raw data and algorithms are archived via RDR by KU Leuven under
the link XXX.

###### Acknowledgements.

The work of the author J.F. is funded by the LOEWE Program of Hesse State
Ministry for Higher Education, Research and the Arts within the project
”Uniform detection and modeling of slums to determine infrastructure needs”.
We also want to thank Nikita Frolov for his input and constructive
discussions.

##  Appendix A Formulation of reaction-diffusion equations

With the definitions from Section II we formulate balance equations for the
respective agents $u^{\prime},v^{\prime}$:

$\displaystyle\begin{split}\dot{N}_{u}^{\prime}&=\frac{\partial}{\partial
t}\int_{A}u^{\prime}\,dA=\int_{A}\hat{U}R\>f(u^{\prime},v^{\prime})-\oint_{C}\bm{J}_{u}^{\prime}\cdot\bm{n}\,dC,\\\
\dot{N}_{v}^{\prime}&=\frac{\partial}{\partial
t}\int_{A}v^{\prime}\,dA=\int_{A}\hat{V}R\>g(u^{\prime},v^{\prime})-\oint_{C}\bm{J}_{v}^{\prime}\cdot\bm{n}\,dC.\end{split}$
(15)

Here $N$ describes the amount of population or agriculturally used area in the
finite area $A$ and accordingly $\dot{N}$ describes the change in the whole
area, whereas $u^{\prime}$ and $v^{\prime}$ describe local changes. The long-
distance effects are a product of the reaction terms $f(u,v)$ or $g(u,v)$ and
the reaction rates $\hat{U}R$, $\hat{V}R$. Here, $u:=u^{\prime}/\hat{U}$ and
$v:=v^{\prime}/\hat{V}$ are dimensionless by division with reference or
maximum densities $\hat{U},\hat{V}$.

Similarly to [22], the short-distance effects are also driven by a density
gradient which can be modeled with Fick’s first law. By applying Gauss’
theorem, we get the two reaction diffusion equations:

$\displaystyle\begin{split}u_{t}&=\hat{U}R\>f(u,v)+D_{u}\nabla^{2}u,\\\
v_{t}&=\hat{V}R\>g(u,v)+D_{v}\nabla^{2}v.\end{split}$ (16)

With the additional dimensionless transformations $t:=Rt^{\prime}$,
$\bm{x}=\bm{x^{\prime}}\sqrt{R/D_{u}}$,$D:=D_{v}/D_{u}$ we derive the
dimensionless standard form of reaction-diffusion equations:

$\displaystyle\begin{split}u_{t}&=\nabla^{2}u+Rf(u,v)\\\
v_{t}&=D\nabla^{2}v+Rg(u,v)\end{split}$ (17)

##  Appendix B Linear stability analysis

As done in [22] and following [39] we perform a linear stability analysis
around the linearized state of Eq. (LABEL:eq:RD_general) with $u=U+\delta
u\,,u=V+\delta v$ with the homogeneous solutions $U,V$, With the the
perturbation ansatz $\delta u=\mathcal{R}[\delta\hat{u}\,\text{exp}(\sigma
t+i\bm{k}\bm{x})]$ or vice-versa with $v$, we derive an eigenvalue problem
with the eigenvalue $\sigma$, the Kronecker delta $\bm{\delta}$,
$\bm{u}=(u,v)$, the Jacobi $\bm{J}(f,g)$ and $\bm{D}=0$:

$\displaystyle\left(\sigma\bm{\delta}-\bm{J}(f,g)\right)\delta\bm{\hat{u}}$
$\displaystyle=0,$ (18)
$\displaystyle\rightarrow\quad\sigma^{2}-\bm{J}(f,g)\bm{I}\,\sigma+\text{det}\left(\bm{J}(f,g)\right)$
$\displaystyle=0.$ (19)

Solving the eigenvalue problem results in two conditions for the Jacobi matrix
$\bm{J}(f,g)$,

$\displaystyle\bm{J}=\begin{pmatrix}f_{u}\big{|}_{0}&f_{v}\big{|}_{0}\\\
g_{u}\big{|}_{0}&g_{v}\big{|}_{0}\end{pmatrix}.$ (20)

that lead to instability,

$\displaystyle f_{u}\big{|}_{0}+g_{v}\big{|}_{0}$ $\displaystyle<0,$ (21)
$\displaystyle\text{det}\left(\bm{J}(f,g)\right)=g_{v}\big{|}_{0}\,f_{u}\big{|}_{0}-g_{u}\big{|}_{0}\,f_{v}\big{|}_{0}$
$\displaystyle>0.$ (22)

As described in [22] the only reasonable formulation of the Jacobi matrix is:

$\bm{J}=\begin{pmatrix}+&+\\\ -&-\end{pmatrix}.$ (23)

Other forms where the column-wise signs are the same results in concentrations
spatially in phase and the form with row-wise same signs only the shown can be
suitably used as shown in Section II.

As Turing patterns can arise due to diffusion, we as well study the short-
distance effects. We can reformulate diffusion as a product of specific energy
$k_{B}T$, with the Boltzmann constant $k_{B}$ and temperature $T$, and the
mobility $\mu$. At constant $T$, the ratio $D=\mu_{v}/\mu_{v}$ results in,

$\displaystyle\bm{B}:=\bm{J}(f,g)-\bm{D}k\quad\text{with}\quad\bm{D}=\begin{pmatrix}1&0\\\
0&D\end{pmatrix}$ (24)

and allows to rewrite the eigenvalue problem in Eq. (18) to:

$\displaystyle\left(\sigma\bm{\delta}-\bm{B}(f,g)\right)\delta\bm{\hat{u}}=0.$
(25)

This results again in two conditions for the Jacobi:

$\displaystyle(f_{u}\big{|}_{0}+g_{v}\big{|}_{0}-k^{2}(1+D)$
$\displaystyle<0,$ (26)
$\displaystyle\text{det}(\bm{B})=(f_{u}\big{|}_{0}-k^{2})(g_{v}\big{|}_{0}-k^{2})-g_{u}\big{|}_{0}\,f_{v}\big{|}_{0}$
$\displaystyle>0.$ (27)

Turing instability is achieved when the condition Eq. (27) is violated,
resulting in the necessary condition for the diffusion induced instability:

$D\,f_{u}\big{|}_{0}+g_{v}\big{|}_{0}>0\quad\rightarrow\quad
f_{u}\big{|}_{0}>-\frac{g_{v}\big{|}_{0}}{D}$ (28)

##  Appendix C Region Information

Here we attach the geographical data of the regions used to demonstrate the
workflow of our method, see Table 3. The coordinates are given in decimal
degrees in reference system WGS84.

Table 3: Coordinates of the regions of interest.  Region |  | West | South | East | North
---|---|---|---|---|---
India |  | 75.3855 | 28.8265 | 77.4804 | 30.6380
| ROI 1 | 76.0086 | 29.4155 | 75.8528 | 29.5319
| ROI 2 | 76.7096 | 29.7066 | 76.5538 | 29.8230
| ROI 3 | 77.2547 | 30.2305 | 77.0990 | 30.3469
Egypt |  | 29.8650 | 29.9671 | 32.1010 | 31.8803
| ROI 1 | 30.3377 | 30.7896 | 30.1804 | 30.9044
| ROI 2 | 31.1260 | 30.6174 | 30.9686 | 30.7322
| ROI 3 | 31.5200 | 31.2488 | 31.3627 | 31.3636
Nigeria |  | 7.3774 | 11.1881 | 9.3336 | 13.1056
| ROI 1 | 8.2020 | 11.9172 | 8.0646 | 12.0492
| ROI 2 | 8.2708 | 12.5114 | 8.1333 | 12.6434
| ROI 3 | 8.9580 | 12.1152 | 8.8205 | 12.2473

##  Appendix D Calculation of ANN with moving windows of different sizes

In order to select suitable excerpt sizes from our selected WorldPop data sets
from Tab. 3, we follow a similar approach as in [12]. With varying window
sizes, starting from square windows with side lengths $L$ of 5 km up to 50 km,
we scan over the respective data set calculating the ANN while moving the
windows in North-South or East-West direction by $L/2$ (see Fig. 8). We
calculate the ANN as in [57],

$\text{ANN}=\frac{\frac{\sum_{i=1}^{N}d_{i}}{N}}{\sqrt{\frac{S}{N}}}$ (29)

with $d_{i}$ the distance of a settlement to the next nearest settlement, $N$
the total amount of settlements in a window area $S=L^{2}$. The ANN evaluates
the regularity of a point pattern and assigns a value between 0 and 2.14
describing the distribution of settlements (0– clustered, 1– random, 2.14–
regular). With this we generated contour diagrams of ANN in the selected
regions over time allowing us to estimate a ’characteristic’ window size where
regularity is dominating. By visual inspection we first selected the suitable
window size to be 15 km and three respective ROIs of this size which have the
most regular distribution over the observation time from our contour diagrams
(see Fig. 8 for a part of the obtained results).

Figure 8: Moving window method for ANN calculation with changing $L$ the
window side length and a part of the results for the Punjab region in India
from year 2003. Similar figures for all regions, years and 10 window sizes
between $L=5$km and $L=50$km can be found in the repository.

##  Appendix E Calculation of the characteristic length

We calculated the characteristic length following [40] and [41] which define
it as follows,

$\ell(t)=\frac{2\pi}{\int qp(q,t)dq}.$ (30)

Here, $q$ describes the modes or waves lengths of a Fourier analysis of a
spatial system with their respective probability of occurrence $p(q,t)$,
evaluated in all directions or over $2\pi$ at every point. To simplify the
analysis, we only scan four directions: horizontal, vertical and two diagonal
directions for every point.

## References

  * UN [2019] UN, World population prospects - average annual rate of population change (percentage) (2019).
  * UN [2020] UN, Population facts: Policies on spatial distribution and urbanization have broad impacts on sustainable development (2020).
  * Retief  _et al._ [2016] F. Retief, A. Bond, J. Pope, A. Morrison-Saunders, and N. King, Environmental Impact Assessment Review 61, 52 (2016).
  * Adams  _et al._ [2020] E. A. Adams, J. Stoler, and Y. Adams, American journal of human biology : the official journal of the Human Biology Council 32, e23368 (2020).
  * Nagendra  _et al._ [2018] H. Nagendra, X. Bai, E. S. Brondizio, and S. Lwasa, Nature Sustainability 1, 341 (2018).
  * Thacker  _et al._ [2019] S. Thacker, D. Adshead, M. Fay, S. Hallegatte, M. Harvey, H. Meller, N. O’Regan, J. Rozenberg, G. Watkins, and J. W. Hall, Nature Sustainability 2, 324 (2019).
  * Christaller [1933] W. Christaller,  _Die zentralen Orte in Süddeutschland (German), Central places in Southern Germany_ , 3rd ed. (Wissenschaftliche Buchgesellschaft, 1933).
  * Hudson [1969] J. C. Hudson, Annals of the Association of American Geographers 59, 365 (1969).
  * Yang  _et al._ [2016] R. Yang, Q. Xu, and H. Long, Journal of Rural Studies 47, 413 (2016).
  * AbouKorin [2018] A. A. AbouKorin, Ain Shams Engineering Journal 9, 1819 (2018).
  * Friesen  _et al._ [2018] J. Friesen, H. Taubenböck, M. Wurm, and P. F. Pelz, Habitat International 73, 79 (2018).
  * Henn  _et al._ [2020] K. Henn, J. Friesen, J. Hartig, and P. F. Pelz, ISPRS International Journal of Geo-Information 9, 541 (2020).
  * Prokop and Friesen [2021] B. Prokop and J. Friesen, Preprints 10.20944/preprints202104.0752.v1 (2021).
  * Pringle and Tarnita [2017] R. M. Pringle and C. E. Tarnita, Annual review of entomology 62, 359 (2017).
  * Theraulaz  _et al._ [2002] G. Theraulaz, E. Bonabeau, S. C. Nicolis, R. V. Solé, V. Fourcassié, S. Blanco, R. Fournier, J.-L. Joly, P. Fernández, A. Grimal, P. Dalle, and J.-L. Deneubourg, Proceedings of the National Academy of Sciences of the United States of America 99, 9645 (2002).
  * Grohmann  _et al._ [2010] C. Grohmann, J. Oldeland, D. Stoyan, and K. E. Linsenmair, Insectes Sociaux 57, 477 (2010).
  * Tarnita  _et al._ [2017] C. E. Tarnita, J. A. Bonachela, E. Sheffer, J. A. Guyton, T. C. Coverdale, R. A. Long, and R. M. Pringle, Nature 541, 398 (2017).
  * Losiri  _et al._ [2016] C. Losiri, M. Nagai, S. Ninsawat, and R. Shrestha, Sustainability 8, 686 (2016).
  * Batty and Milton [2021] M. Batty and R. Milton, Urban Studies 58, 3071 (2021).
  * May [1976] R. M. May, Nature 261, 459 (1976).
  * Turing [1990] A. M. Turing, Bulletin of Mathematical Biology 52, 153 (1990).
  * Pelz  _et al._ [2019] P. F. Pelz, J. Friesen, and J. Hartig, Physical Review E 99, 022302 (2019).
  * Friesen  _et al._ [2019] J. Friesen, R. Tessmann, and P. F. Pelz, Proceedings of the 5th International Conference on Geographical Information Systems Theory, Applications and Management - GISTAM (2019).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Proceedings of the National Academy of Sciences 113, 3932 (2016).
  * Reinbold  _et al._ [2021] P. A. K. Reinbold, L. M. Kageorge, M. F. Schatz, and R. O. Grigoriev, Nature Communications 12, 3219 (2021).
  * Ermolaev  _et al._ [2022] A. V. Ermolaev, A. Sheveleva, G. Genty, C. Finot, and J. M. Dudley, Scientific Reports 12, 12711 (2022).
  * Hoffmann  _et al._ [2019] M. Hoffmann, C. Fröhner, and F. Noé, Journal of Chemical Physics 150, 10.1063/1.5066099 (2019).
  * Mangan  _et al._ [2016] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz, IEEE Transactions on Molecular, Biological and Multi-Scale Communications 2, 52 (2016), 1605.08368 .
  * Rudy  _et al._ [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Science advances 3, e1602614 (2017).
  * Schaeffer [2017] H. Schaeffer, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 20160446 (2017).
  * Fasel  _et al._ [2022] U. Fasel, J. N. Kutz, B. W. Brunton, and S. L. Brunton, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478, 10.1098/rspa.2021.0904 (2022).
  * Hirsh  _et al._ [2022] S. M. Hirsh, D. A. Barajas-Solano, and J. N. Kutz, Royal Society Open Science 9, 10.1098/rsos.211823 (2022).
  * Gong  _et al._ [2020] P. Gong, X. Li, J. Wang, Y. Bai, B. Chen, T. Hu, X. Liu, B. Xu, J. Yang, W. Zhang, and Y. Zhou, Remote Sensing of Environment 236, 111510 (2020).
  * Liu  _et al._ [2020] X. Liu, Y. Huang, X. Xu, X. Li, X. Li, P. Ciais, P. Lin, K. Gong, A. D. Ziegler, A. Chen, P. Gong, J. Chen, G. Hu, Y. Chen, S. Wang, Q. Wu, K. Huang, L. Estes, and Z. Zeng, Nature Sustainability 3, 564 (2020).
  * Worldpop [2019] Worldpop, Worldpop - population density (2019).
  * Note [1] The computation of the ANN, the resulting selection of the ROIs, and the later computation of the feature lengths was done with GAIA because it provides discrete spatial data (1 - built-up, 0 - no built-up area). As such, and in contrast to WorldPop, which distributes population quantities across administrative regions, leading to non-existent values of 0, the identification of settlement locations and sizes is straightforward. However, for model identification WorldPop is used, since most pattern forming mechanisms describe concentration distributions of different (chemical or other) species and do not create discrete patterns [37, 39], similar to WorldPop.
  * Murray [2003] J. D. Murray,  _Mathematical biology. II Spatial models and biomedical applications_ (Springer, 2003) p. 811.
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, The Journal of Chemical Physics 28, 258 (1958).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Reviews of Modern Physics 65, 854 (1993).
  * Puri [2004] S. Puri, Phase Transitions 77, 407 (2004).
  * König  _et al._ [2021] B. König, O. J. Ronsin, and J. Harting, Physical Chemistry Chemical Physics 23, 24823 (2021), 2107.07234 .
  * Christiansen  _et al._ [2020] H. Christiansen, S. Majumder, M. Henkel, and W. Janke, Physical Review Letters 125, 180601 (2020).
  * Champion  _et al._ [2020] K. Champion, P. Zheng, A. Y. Aravkin, S. L. Brunton, and J. N. Kutz, IEEE Access 8, 169259 (2020).
  * Zheng  _et al._ [2019] P. Zheng, T. Askham, S. L. Brunton, J. N. Kutz, and A. Y. Aravkin, IEEE Access 7, 1404 (2019).
  * Akaike [1973] H. Akaike, In proceedings of 2nd International Symposium on Information Theory , 267 (1973).
  * Mangan  _et al._ [2017] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 10.1098/rspa.2017.0009 (2017).
  * Burnham and Anderson [2004] K. P. Burnham and D. R. Anderson, Sociological Methods & Research 33, 261 (2004).
  * Thaler  _et al._ [2019] S. Thaler, L. Paehler, and N. A. Adams, Journal of Computational Physics 397, 108851 (2019).
  * Allen and Cahn [1979] S. M. Allen and J. W. Cahn, Acta Metallurgica 27, 1085 (1979).
  * Naozuka  _et al._ [2022] G. T. Naozuka, H. L. Rocha, R. S. Silva, and R. C. Almeida, Nonlinear Dynamics , 1 (2022).
  * Rudy  _et al._ [2019] S. Rudy, A. Alla, S. L. Brunton, and J. N. Kutz, SIAM Journal on Applied Dynamical Systems 18, 643 (2019).
  * Bakarji  _et al._ [2022] J. Bakarji, K. Champion, J. N. Kutz, and S. L. Brunton, ArXiv (2022), arXiv:2201.05136 .
  * Ouala  _et al._ [2023] S. Ouala, S. L. Brunton, B. Chapron, A. Pascual, F. Collard, L. Gaultier, and R. Fablet, Physica D: Nonlinear Phenomena 446, 133630 (2023).
  * Lu  _et al._ [2022] P. Y. Lu, J. Ariño Bernad, and M. Soljačić, Communications Physics 5, 206 (2022), 2107.10879 .
  * de Silva  _et al._ [2020] B. de Silva, K. Champion, M. Quade, J.-C. Loiseau, J. Kutz, and S. Brunton, Journal of Open Source Software 5, 2104 (2020).
  * Kaptanoglu  _et al._ [2022] A. Kaptanoglu, B. de Silva, U. Fasel, K. Kaheman, A. Goldschmidt, J. Callaham, C. Delahunt, Z. Nicolaou, K. Champion, J.-C. Loiseau, J. Kutz, and S. Brunton, Journal of Open Source Software 7, 3994 (2022).
  * Clark and Evans [1954] P. J. Clark and F. C. Evans, Ecology (1954).