Birth, Annexation, and Squeezing of Cities in a Prefecture: Can the Ranking of Competitive Areas of Municipalities Obey the Authentic Power Law?

As the first step for revealing potential rules inherent in cities that are closely squeezed in a sectioned domain, municipalities in the entire prefecture in Japan are considered and their distributions of the areas are analyzed in details according to a rank-size procedure. Computed results suggest that among the population, area, and population density, the last becomes the most important factor in finding the rank-size rule. Indeed, of the 47 Japanese prefectures the Metropolis of Tokyo and Fukuoka Prefecture exhibit the most typical rank-size rules, where the former possesses the exceptionally high population density as well as urbanized rate. The underlying mechanism of the rule can be supported by a toy model with a tournament game using a sequence of random numbers, where teams (municipalities) are highly competitive in gaining the final wins (broadest territory). A stability analysis implying perturbations due to global warming allows one to confirm unexpected robustness of the rank-size relation. Finally, the authenticity of the log-log relation in the rank-area data of Tokyo Metropolis is tested statistically.


INTRODUCTION
Cities, towns, and villages have long been regarded as an attractive object of scientific study [1][2][3]. In the context of applied statistical physics, their dynamics, configuration, geography, and population have been studied from a variety of viewpoints [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In an effort to reveal implicit rules that govern the organized whole, the possibility of bearing rank-size rules such as the Zipf's law and its variants has been constantly discussed for realistic cities that are scattered all over the world [7,8,13,15,18,19]. Besides the conventional approach to cities, in recent years, arguments from a new perspective have been assembled. The most typical examples can be seen in quantitative studies on the basis of social interaction, information entropy, environmental change, and sustainability [21][22][23][24]. It should be noted here that in most literature the key term "size" has been employed implicitly as the meaning of population. One should notice, however, that besides the population there are two meanings in the term, namely, area and population density. In this paper we focus our principal attention on the former. Specifically, areas of municipalities (not only cities but towns and villages being included) that are squeezed in a prefecture of a country are considered. For an example we consider Japan because this country has been divided into 47 prefectures; in each prefecture, typically tens of municipalities are closely packed with a unique configuration that has been arranged according to a self-organized process over a long period. The distributions of the areas are analyzed in detail according to a rank-size procedure. Through numerical results, one can identify, among the populations, areas, and population densities of prefectures, a key factor in finding the rank-size rule. Indeed, of the 47 Japanese prefectures, only four (specifically, the Metropolis of Tokyo; Hiroshima, Shimane, and Fukuoka Prefectures) are shown to exhibit the typical rank-size rules, where the first in the parentheses possesses the exceptionally high population density as well as urbanized rate. The underlying mechanism of the rule is explained by a toy model with a tournament game using a sequence of random numbers, where teams (municipalities) are highly competitive for gaining the final wins (broadest territory). A stability analysis implying perturbations due to global warming allows one to discuss robustness of the rank-size relation. Finally, the authenticity of the log-log relation in the rank-area data of Tokyo Metropolis is tested with the Durbin-Watson d statistics.

MOTIVATION AND METHODOLOGY
Suppose that you are an inhabitant in a square-shaped prefecture, who is asked by the king to divide the land into unequal squares as few as possible in number. It is assumed further that if the requirement is met, the king is pleased to give you the entire land and at the same time to assign you to the governor of the prefecture. Struggling with the problem for hours or days, finally you will find it too difficult to solve. Actually, this is termed the simplest perfect squared square problem, which, along with the Kepler's closest packing and the four colors problem, had long been one of the most difficult problems in geometry. In 1948, 23 years after the first explicit mention of division into unequal squares by Morón, a compound perfect square with 24 squares was realized by Willcocks [25] (see Figure 1A). Furthermore, after 30 years from his finding, in the long run the ultimately simple perfect squared square with only 21 elements was found by Duijvestijn [26,27] (see Figure 1B), Dutch geometrician, who had devoted himself to solving the extremely hard problem. He applied Kirchhoff's laws to polar nets with which the dissection was obtainable through calculating a current flow [26]. Here we note that for both squares there are pronounced inequalities in the areas of the elements, which may remind one of a realistic division inherent in lands that are constrained with a number of boundaries. Indeed, in 2003, 25 years after the great achievement by the Dutch scholar the areas y in the squared squares were found to obey the statistical rule [5] that corresponds to a q-extension of the logarithmic distribution where log abbreviates the common logarithm; x represents the rank in descending order; a, b, and q are positive real constants to be determined with the least squares fit. The analyzed results for the squares in Figure 1 are listed in Table 1 [5]. Here n is the number of elements; the accuracy of the regression model can be examined by the degree of fit, |r|, namely with the Pearson's coefficient (0 < |r| < 1), and with the Durbin-Watson ratio, for i 1 to n -1, where for i 1 to n; Y y q ; and the prime on Y indicates the point on the regression line. In Table 1  negative counterpart, i.e., 2 < d ≤ 4, d must be replaced with 4d. In contrast to non-ranking as well as rank-rank data (such as, e.g., Spearman's and Kendall's approaches), for the rank-size analysis, in most cases the positive correlations are included between the neighboring data. Note that along with detecting autocorrelations in the time-series data the Durbin-Watson d statistics was used also in other data; typical examples are seen in the Rietveld analysis of powder diffraction data [29][30][31][32] and in the Lorentzian fits in the 151 Eu Mössbauer spectroscopy of oxide glasses [33]. More recently has the statistics been applied to testing assumptions of independent residuals in the physics education research [34].
Besides the science of cities, to date, sustained efforts have been made to find nontrivial rules in the ranking of a variety of complex systems, not only in conventional physics but in  Frontiers in Physics | www.frontiersin.org February 2022 | Volume 9 | Article 789571 4 geometry, geography, demography, linguistics, and sciences on social phenomena [35][36][37][38][39][40][41][42][43]. More recently has ranking been regarded as a tool useful for condensing large-scale data that have been accumulating in contemporary sciences such as, e.g., computational metallurgy [44] and gravitational wave astronomy [45], though the results are not yet ready for finding a rule.

GEOMETRY OF PROBLEM
In this paper we consider the entire prefecture in Japan, which is divided into 47 prefectures that include Hokkaido and the Metropolis of Tokyo. In Figure 2, a complete map of A) Japan is given along with the ones of B) the Metropolis of Tokyo (n 53, excluding remote islands), and C) Fukuoka Prefecture (n 60). Here n indicates the number of municipalities squeezing in a prefecture, which consist of cities, towns, and villages. The two prefectures are to be cited in the rank-size analysis, and are marked with the blue and red ink, respectively, in Figure 2A. The colored dashed lines in Figures 2B,C indicate the boundaries among municipalities in the prefecture.
Since the opening of Japan to the West, this country has experienced three major municipal consolidations [46]

RESULTS
Below numerical results for the areas in the 47 Japanese prefectures [46] are given. Note that in Figures 3, 4, truncation is made for the exceptional datum of Ehime Prefecture (q 3.32), but it is included for the calculation of regression lines. For Fukui Prefecture, because of an ill-posed nature, there has been no optimal solution available: # The three frames in Figure 3 show a cross-sectional view of the three-dimensional scattergram of optimal solutions for the entire prefecture in Japan. The purple line in each individual cross section indicates the regression line. It can be seen that for all the cross sections the aggregation of dots is highly clustered around a "center of gravity": specifically, A) 0.96 ≲ |r| ≲ 1.00 and 0.5 ≲ d ≲ 2.0 in the d-|r| plane, B) 0.96 ≲ |r| ≲ 1.00 and 0.2 ≲ q ≲ 1.5 in the q-|r| plane, and C) 0.5 ≲ d ≲ 2.0 and 0.2 ≲ q ≲ 1.5 in the q-d plane. Correlation analysis has shown that irrespective of the projective directions the dots are randomly scattered in each cluster. It can be found in Figure 3A that, in contrast to the dots' aggregation in 0.99 < |r|, the d values exhibit concentrations in the vicinity of unity, indicating that, for reasons identical to quantifying serial correlations between adjacent least-squares residuals in Rietveld refinements of step-scan power diffraction data [29][30][31][32], the use of the Durbin-Watson d statistics is necessary for the rank-size analysis.
# Figure 4 shows the dependence of scaling exponent, q, in the left-hand side of Eq. 1 on the geographic data of the entire prefecture in Japan. First, as specific data we consider A) the population and B) the area (km 2 ). Through correlation analysis it has been found that for the former, there exists a negative correlation (R −0.3833 with d 2.113), in contrast to a positive correlation (R 0.4177 with d 2.197) for the latter. Here R represents the correlation coefficient between the optimal parameters and the geographic data. The results allow one to expect that there will exist a stronger correlation between the scaling exponent and the population density (km −2 ) that is defined with the population divided by the area (km 2 ). The plot on Figure 4C indicates that, as has been expected, there is a negative correlation with the maximum magnitude in the correlation |R| among the three (R −0.5522 with d 2.174). Frontiers in Physics | www.frontiersin.org February 2022 | Volume 9 | Article 789571 # Above we have confirmed the noticeable dependence of the scaling exponent, q, on the population density. Therefore, it appears interesting to examine how the other two parameters, |r| and d, depend on the population density (km −2 ). Figure 5A shows the dependence of the degree of fit, |r|. It can be seen that, although in comparison with Figure 4C the magnitude of the correlation between |r| and the density reduces substantially, the parameter is positively correlated with the density (R 0.2215 with d 2.528). For the Durbin-Watson ratio, however, as is apparent in Figure 5B, correlation vanishes (R -0.0095 with d 2.206).
# As has been seen in Figure 3A, there are four prefectures that meet the requirement: 0.99 < |r| and d U < d. The geographic data of the four prefectures are listed in Table 2. Figure 6A plots the rank dependence of the areas (km 2 ) of municipalities on the main land in the Metropolis of Tokyo. The purple declining line indicates the optimized fit (|r| 0.9961 with d 1.861 > d U for q 0.21 and n 53). As is found in Figures 3B,C, the prefecture shows minimum in the scaling exponent. The entire Metropolis includes nine outlying municipalities on islets in the Pacific Ocean [46]. Analysis of the extended region (n 62) has yielded the optimal   , which is consistent with our expectations that expanding territories give rise to the reduction of competitions. # Figure 6B plots the rank dependence of the areas (km 2 ) of municipalities in Fukuoka Prefecture. The purple declining line indicates the optimized fit (|r| 0.9942 with d 1.493 > d U for q 0.53 and n 60). In both plots (Figures 6A,B) the dots are found to be hierarchically clustered along the regression line, the feature of which will be reviewed in performing a modeling. It has been found empirically that without the hierarchy the positive correlation arises in the alignment of dots, which results in d < d U . Note again that in the rank-size analysis the serial correlation analysis is indispensable. In summary, along with the other two prefectures (Hiroshima and Shimane Prefectures) the numerical results are given in Table 3.

Constructing a Toy Model
Along with data on Tables 2, 3, the results of Figure 6 suggest that municipalities in a prefecture are in intense competition for the extension of their areas, which might remind one inevitably of the tournament game in sports. Indeed, in recent years sports games have become one of the fascinating topics in applied statistical physics [40,[47][48][49][50][51]. As a specific sport, here we consider a virtual rugby with 64 teams competing on the tournament diagram as depicted in Figure 7A. Because the sport is chosen solely as an easy example for finding an analogy to the territory expansion, one may arbitrarily replace this with another such as basketball, bowling, cricket, soccer, or cycle race. For the entire team, according to a Monte Carlo procedure, wins and losses of each game are determined stochastically by the use of two-digit random numbers. After the 62 games being finished the cumulative scores of the 64 teams are recorded. Subsequently, to enhance the success-breeds-success effect the scores are modified by calculating the (1/q)-th power of them. For instance, for a team (assuming team #i) that has survived till the final, the modified cumulative score, S i , can be written as while for a team (assuming team #k) that has lost in the first round, the counterpart, S k , can be given as Here s ij represents the score for team i (i 1-64) at round j (j 1-6). Consequently a non-square 64 × 6 matrix is generated. In comparison of the two teams in the scores it is apparent that for 0 < q < 1 the team (team #i) that has survived the final gains far greater advantage than the one (team #k) that has lost in the first round. Note that in the context of economics the effect due to a positive feedback mechanism is equivalent to the so-called richget-richer effect, and in the context of sociology it corresponds to the Mathew effect [51][52][53][54][55]. Besides, our model could be consistent with the one that was once portrayed by Calderón (1847-94), who modeled society as a thermodynamic machine [56]. There the temperature gradient was represented by the gap between rich and poor. Figure 7B shows a typical example of the rank dependence of the q-th power of cumulative scores of the entire team for a tournament using a sequence of two-digit random numbers. The purple declining line indicates the optimized fit (|r| 0.9963 with d 1.477 > d U with n 57; n < 64 because of 7 ties being included). In the plots the relation of Eq. 1 is found to be reproduced.
In recent years, efforts have been made to establish an analytical methodology for deriving stationary distributions of complex systems with unidirectional random growth with resetting [57][58][59]. To date the method using master equations has succeeded in finding extensive applications such as degree distributions in networks [60,61], distributions in citations [62], and income distributions [63,64]. It seems, however, that the possibility of application of the analytical approach to our case remains pending, because the highly complex morphologies ( Figures 2B,C) in the present system have been organized through a long-term blend of spontaneousness and intermittent stimulations by the administration. Moreover, even if the unidirectional requirement can be met, the interference based on an administrative law might be capricious, sporadic, and far from periodic.

Analyzing History of the Metropolis of Tokyo
To examine the formation of the complicated geometrical configuration we consider the history of the Metropolis of Tokyo. Computed results from 1955 to 2005 are given in Table 4. Here d L denotes the lower critical value in the Durbin-Watson d statistics [28]; note that with d < d L the null hypothesis is unconditionally rejected (α 0.01). It can be seen that the two parameters, |r| and d, increase constantly over the entire period. It should be noted here that the prefecture (the Metropolis of Tokyo) is a very rare exception in the entire Japanese prefecture, because of the absence of large-scale municipal consolidations [37].

Examining Effects Due to Global Warming: A Stability Analysis
Not to mention islets in oceans, a number of cities on the globe are now confronting a sign of the inundation crisis due to global warming [65][66][67][68]. As is well known, the Japanese Islands (see Figure 2A) are surrounded by the four seas: the Sea of Okhotsk, Pacific Ocean, East China Sea, and Japan Sea, indicating that along with Dhaka, Jakarta, and Miami, coastal megacities in Japan are far from exceptional in the susceptibility to the forthcoming aqua-disaster caused by the climate crisis [69]. Figure 8 shows the dependence of optimal parameters on the inundated rate on the supposition that effects due to global warming will get more and more inevitable in the future. In the Metropolis of Tokyo and Fukuoka Prefecture, respectively, there are 6 and 14 municipalities that are located along the coast. Here it is assumed that these are uniformly inundated. It is surprising to note that both the scaling exponent ( Figure 8A) and the degree of fit ( Figure 8B) preserve stabilities against the large-scale perturbation. The results of the Durbin-Watson ratio ( Figure 8C), however, exhibit a sharp contrast. Namely, for the Metropolis of Tokyo the value of d experiences a substantial reduction with increasing the inundated rate, whereas for Fukuoka Prefecture a phase transition occurs just in front of the 50% inundation, and the new phase is maintained to the complete inundation. Note that for the Metropolis of Tokyo the phase transition that occurs just in front of the complete inundation is due to the simultaneous vanishment of the 6 coastal cities (i.e., n 53 → n 47).

Testing Stretched Exponentiality
For the value of q much smaller than unity, Eq. 1 becomes compatible with its dual counterpart Here, for 0 < q < 1 the exponential decay is stretched along the x axis. The qualitative difference from Eq. 1 is seen in the appearance of a rank that characterizes the decay along the axis, where e is the Napier's constant. Figure 9 shows the stretched exponentiality for the Metropolis of Tokyo (n 53). Here we note that there are two solutions: A) Mode #1 (q 0.09; |r| 0.9948 with d 2.002), and B) Mode #2 (q 0.23; |r| 0.9953 with d 2.046). The reason why the bimodal solutions are possible can be explained as follows: because the dependence of d against q exhibits a convex curve, there exists a single maximum, assuming here (q 0 , d 0 ). If d 0 exceeds 2, the curve crosses the q axis twice (assuming q L and q U ) across the center q q 0 (i.e., q L < q 0 < q U ). Namely, at q q 0 zd zq 0, z 2 d zq 2 < 0, and d > 2.
Hence one obtains two solutions: (q L , 2) and (q U , 2), for which all the blue dots in Figure 9 are randomly distributed on the regression line (log y versus x q ). Along with the characteristic rank ρ, comparison of the two stretched exponential modes is given in Table 5.

Testing Power Law
To comprehend the link with the power law (i.e., log-log) relation [70], with the use of the Box-Cox transformation [71], Eq. 1 will be rewritten as log e y q − 1 q a′ − b′ log x, a′ log e (a − 1) q, b′ log e b q.
In Figure 10 the log-log plot is given for rank dependence of the areas of municipalities in the Metropolis of Tokyo (n 53). The purple declining line indicates the optimized fit (|r| 0.9912 with d 0.792 < d L 1.34). Because the Durbin-Watson ratio is much smaller than the lower critical value d L , it is judged that the resemblance to Figure 6A is solely apparent and therefore the pseudo-power law is far from authentic.

CONCLUSION
As the first step for revealing potential rules inherent in cities, towns, and villages that are closely squeezed in a sectioned domain, municipalities in the entire prefecture in Japan have been considered and their highly squeezed distributions of the areas have been analyzed according to a rank-size procedure. Through computed results it has been confirmed that among the population, area, and population density, the last becomes the most important factor in finding the rank-size rule. Indeed, of the 47 Japanese prefectures the Metropolis of Tokyo exhibits the exceptionally high density. The remark can be supported by the model using a tournament game, where teams (municipalities) are highly competitive for gaining the final wins (broadest territory). Finally, the study of this paper remains open. To seek all over the world for prefectures (including their analogs) bearing the authentic power law, international cooperations are needed. The most likely candidates are the ones with population densities comparable to or higher than the one of Tokyo Metropolis. If on the globe there is no prefecture that passes the statistical testing, the proposition on the interrogative subtitle of this paper will be truly turned down.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.