# Intermunicipal Travel Networks of Mexico (2020-2021)

Oscar Fontanelli1 Plinio Guzmán2 Amilcar Meneses3 Alfredo Hernández4 Marisol
Flores-Garrido5 Maribel Hernández-Rosales1∗ Guillermo de Anda-Jáuregui6,7,8,∗

(1 CINVESTAV, Irapuato, Mexico
6 Computational Genomics Division, National Institute of Genomic Medicine,
Mexico City, Mexico
7 Cátedras Conacyt Para Jóvenes Investigadores, National Council on Science
and Technology, Mexico City, Mexico
8 Centro de Ciencias de la Complejidad (C3), Universidad Nacional Autónoma de
México , Mexico City, Mexico
2 Astronomer LTD
3 CINVESTAV CDMX, Mexico
4 Centro de Ciencias Genómicas, UNAM
5 ENES Morelia, UNAM
∗<EMAIL_ADDRESS>
∗<EMAIL_ADDRESS>
)

###### Abstract

We present a collection of networks that describe the travel patterns between
municipalities in Mexico between 2020 and 2021. Using anonymized mobile device
geo-location data we constructed directed, weighted networks representing the
(normalized) volume of travels between municipalities. We analysed changes in
global (graph total weight sum), local (centrality measures), and mesoscale
(community structure) network features. We observe that changes in these
features are associated with factors such as Covid-19 restrictions and
population size. In general, events in early 2020 (when initial Covid-19
restrictions were implemented) induced more intense changes in network
features, whereas later events had a less notable impact in network features.
We believe these networks will be useful for researchers and decision makers
in the areas of transportation, infrastructure planning, epidemic control and
network science at large.

## Introduction

Intermunicipal mobility is a kind of medium and large scale mobility within a
country where millions of individuals daily travel from one county or
municipality to another, either going from home to work, shopping, accessing
public services, cargo loading, vacation, etc. These movements and travels
generate complex structures and dynamics of socio-economic interactions
between different areas both at regional and national levels.

Given the nature of these mobility systems, complex networks have been widely
adopted to model this commuting phenomena [6][10][23][25][35]. Characterizing
and understanding the properties of these mobility networks is crucial for
decision-making, urban planning, traffic engineering and, as has become clear
with the Covid-19 pandemic, designing, implementing and evaluating mobility
restrictions and lockdowns in order to contain or control the epidemic spread
[3][17][26][34][38][40][45]. Therefore the need for public mobility network
datasets, aiming to researchers and decision-makers in the areas of mobility,
urban planning, epidemic control, etc.

Traditionally these networks had been obtained from mobility surveys. However,
the emergence in recent years of cell phone and GPS data have facilitated the
acquisitions of accurate and large sets of human mobility data, thus allowing
the construction and characterization of large and detailed human mobility
networks.

Here we introduce a new public dataset of daily intermunicipal origin-
destination networks in Mexico for 2020 and 2021, which were directly
constructed from large datasets of geolocation data. With this dataset we hope
to contribute to research and decision-making communities from diverse
interests, from pure network theory to those studying human mobility, urban
planning, national scale social and economic relations, epidemic control, etc.
As far as we know, this is the first public dataset of its kind for a country
in Latin America.

There are other works that have built origin-destination networks aiming to
describe mobility patterns, measure the volume of public transportation and to
plan public transportation, among other things. These mobility patters include
large-scale and long-range commuting patterns [20][36]; spatio-temporal
patterns for different socioeconomic strata [25]; patterns in bike-sharing
systems [24], etc. Origin-destination matrices have also been used to measure
the volume of use of public transport [2][44][28] and public transport
planning [29][1].

There are several approaches for modeling and generating origin-destination
matrices. These approaches include gravity models [42][39][16],Bayesian model
[31][32][43], linear assignment matrix approximation [41], Principal
Components Analysis [12] and gradient approximation method [18], among others.

Origin-destination networks have been elaborated from different types of data.
Most of these works use data from the different public transport systems or
from road side interview [42]. In works such as [2] and [28] authors use smart
card fare data to estimate origin-destination networks. In [25] the data
considered for these matrices are those of the bike-shared system. In [44]
authors conduct experiments on two datasets generated by ride-hailing
applications. In [41], authors estimate the matrices using data from traffic
counts on the network links.

Recently, other works have attempted to utilize social network data from
Twitter or Facebook [30][33][5][4][15] [19] or another data set from mobile
phone data such ACAPS dataset or SafeGraph [37][7][11]. In particular, Edsberg
et al. [15] show that origin-destination networks based on data from cell
phones and social networks present high-quality results.

This article is organized as follows: Results presents and describes the
collection of mobility networks and it shows analysis of changes in global
(sum of weights), local (centrality measures), and mesoscale (community
structure) network features. In Discussion we argue how events in early 2020
(when initial Covid-19 restrictions were implemented) induced more intense
changes in network features, whereas later events had a less notable impact in
network features. In Methods we describe the methodology we utilized to
collect data and show the algorithm for network construction.

##  1 Results

### A travel network dataset

We release a public dataset of 731 intermunicipal origin-destination networks
in Mexico. These networks were constructed from a large and anonymized mobile
location dataset (see [9] for more information about this dataset). Each
network is the intermunicipal origin-destination network in Mexico for each
day during the 2020-2021 period. Nodes represent municipalities (third level
administrative division) or official metropolitan zones, see Methods section
below. These are weighted and directed networks, where the weight of edge
$a_{ij}$ is equal to the total number of observed travels from node $i$ to
node $j$ normalized by the different number of mobile devices we recorded on
that day. The data set with these 731 networks is freely available in a OSF
repository http://dx.doi.org/10.17605/OSF.IO/42XQZ.

For analysis and visualization purposes we chose nine representative dates
capturing different important events during the evolution of the pandemic in
Mexico; these are shown in Table 1. As a first visualization, we show in Fig.
1 mobility networks over the Mexico map for this set of dates, marking only
1$\%$ of edges with highest weight.

Date | Event
---|---
2020-02-24 | First reported Covid-19 case.
2020-03-23 | Beginning of official lockdowns (National Program of Social Distance).
2020-06-01 | Beginning of Epidemiological Stoplight Program.
2020-07-30 | First national peak of daily contagions.
2020-09-21 | Local minimum of daily contagions (between first and second wave).
2021-01-19 | Second national peak of daily contagions.
2021-05-24 | Local minimum of daily contagions (between second and third wave).
2021-08-16 | Third national peak of daily contagions
2021-12-27 | Local minimum of daily contagions (just before fourth wave)
Table 1: Set of dates for the analysis. We chose these dates to be either
Monday or the Monday closest to the referred event. All of these make
reference to events in Mexico. Figure 1: Network visualizations for selected
days.

### Changes in graph total weight sums were observed throughout the first two
years of the pandemic

As a measure to quantify the total observed movement in each network of the
collection, we consider the total sum of weights in the network $S_{G}=\sum
w_{i}$, where $i$ runs over all nodes of the network. In our context, a higher
value of $S_{G}$ can be understood as higher mobility between municipalities
in the country, which, in turn, is associated with people’s decisions to move
outside their locality.

Fig. 2 shows the time series related to the $S_{G}$ parameter, also indicating
the set of dates described in Table 1 and official school vacation periods
(summer breaks, winter breaks and Easter holidays, shaded in gray). It can be
observed that the decay in mobility that happens in January 2020, probably due
to being the post-holiday season, is prolonged after the start of the
lockdown, reaching a local minimum point shortly before the beginning of the
summer holiday season. Following summer break mobility continues to decay,
reaching its lowest point again shortly before the beginning of winter break,
when there was a pronounced mobility rise. For the first half of 2021 we
observe a sustained rise in mobility until July 2020, when it reaches a
relatively high plateau.

Figure 2: Time series of smoothed edge weight sum (weighted density). Vertical
dotted lines correspond to dates shown on Table 1. Shaded areas are official
school holiday periods.

Since the networks exhibit a different mobility pattern each day, centrality
measures associated with each node (locality) also change. In this section we
explore the variability of centrality measures over the period we studied. We
choose three different centrality measures: (a) degree centrality, (b) node
strength, (c) betweenness centrality. Variability is measured through the
coefficient of variation (cv).

Figures 3-5 show the variation of centrality for nine different nodes
(municipalities or metropolitan zones). These nodes were chosen to illustrate
different behaviours across the studied time period.

The regions of Valle de México, Guadalajara, Monterrey and Puebla-Tlaxcala
correspond to densely populated regions. In these areas, the coefficients of
variation for the three different centrality measures are: (30.0, 53.1, 9.7)
in Valle de México, (24.3, 50.3, 17.8) in Guadalajara, (31.3, 49.2, 25.9) in
Monterrey, and (29.0, 44.0, 34.7) in Puebla-Tlaxcala. Morelia, a city in west-
central Mexico, is included because it maintains a highly stable mobility
pattern with coefficients of variation 24.0, 40.3 and 16.0 for the considered
centrality measures. The municipality of San Pedro Topiltepec is an example of
a locality with high degree-centrality variability. The coefficients of
variation for this location are 320.4, 410.1 and 285.8. Tijuana is an
important border city and changes in centrality are expected whenever there
are changes in the regulations on the border with the USA. The coefficients of
variation are 34.9, 41.5 and 89.6. Finally, Acapulco and Cancun are cities
with a beach that represents a popular tourist destination. The coefficients
of variation for these cities are (41.4, 41.1, 60.2) in Acapulco, and (33.3,
41.8, 62.1) in Cancún.

Figure 3: Time series of total degree for nine representative nodes. Figure 4:
Time series of total strength for nine representative nodes. Figure 5: Time
series of betweenness centrality for nine representative nodes.

In order to explore the variability observed in node centrality measures with
the population of the area they represent, the relationship between the three
centrality measures considered and the population reported in the 2020
national census was analyzed
https://www.inegi.org.mx/programas/ccpv/2020/default.html#Datos_abiertos.
Table 2 shows estimated parameters when fitting a linear regression model,
while Fig. 6 shows the identified relationship on a logarithmic scale.

| Node metric | Statistic | Model | R squared
---|---|---|---|---
1 | strength_in | mean | lineal | 0.94
2 | strength_all | mean | lineal | 0.93
3 | strength_out | mean | lineal | 0.93
4 | degree_out | mean | loglog | 0.84
5 | degree_all | mean | loglog | 0.84
6 | degree_in | mean | loglog | 0.84
7 | betweenness_undirected | mean | lineal | 0.74
8 | betweenness_directed | mean | lineal | 0.69
9 | betweenness_undirected | cv | loglog | 0.75
10 | betweenness_directed | cv | loglog | 0.74
11 | strength_out | cv | loglog | 0.50
12 | degree_out | cv | loglog | 0.50
13 | strength_in | cv | loglog | 0.50
14 | strength_all | cv | loglog | 0.48
15 | degree_in | cv | loglog | 0.48
16 | degree_all | cv | loglog | 0.44
Table 2: Results for linear regressions for node metrics (mean and coefficient
of variation) against population. Some of them where in lineal scale and some
other in log scale. Figure 6: left: log-log relationship between mean of daily
centrality measures and population (degree, strength, and betweenness
centrality); right: log-log relationship between coefficient of variation of
daily centrality measures and population (degree, strength, and betweenness
centrality)

### Formation and evolution of communities in the network.

Recalling that a community in a network is a set of nodes with a larger
density of connections between them than external to the set, communities in
these mobility networks correspond to groups of municipalities or with a high
internal (within the group) mobility and a relative low external mobility
(from the region to the outside or the other way around). Community detection
algorithms in mobility or commuting networks have been widely utilized to
detect or delimitate geographical functional regions [13][14][21][22][27].

We utilized the label-propagation algorithm to detect communities in our
networks. Alluvial plot in Fig. 7 shows in schematic way the time evolution of
community structure on the networks. In this diagram each line represents a
municipality in Mexico and they are grouped according to the network community
to which they belong. We show this structure for the chosen set dates. Color
of each line represents the state to which they belong.

Figure 7: Alluvial plot for network communities for the set of nine relevant
dates. Each line represents a node (municipality or metropolitan zone), they
are coloured according to state and for each day they are grouped according to
their community. For each day communities are shown in the vertical bars.
Figure 8: Alluvial plots for network community inside four states in Mexico.
Each line represents a municipality or metropolitan zone, they are coloured
according to their initial community and they are grouped according to the
community to which they belong.

## Discussion

Sum of total (in and out) weights over all edges of these mobility networks is
a gross indicator of the total intermunicipality mobility that we observe on
each day. We can see in Fig 2 how there is a pre-pandemic peak at late
February - early March 2020 (first confirmed case in Mexico was on February
24, 2020). This peak is followed by a persistent decay, whose rate does not
seem to be affected by the beginning of lockdowns. In fact, this decay
continues until June-July of 2020, when we observe a peak that somewhat
coincides with official summer break. After summer break we still see a
downward trend until winter break and, agian, we see a peak that coincides
with official vacation period. In late February 2021 we observe the lowest
point of mobility and from here there is sustained upward train, apparently
unaffected by Easter holidays until it reaches a new maximum that is even
higher than the pre-pandemic peak. At this point mobility reaches a high
plateau and, even though mobility still fluctuates, these fluctuations does
not seem to be correlated neither with the pandemic third wave nor with
vacation periods. From July 2021 total intermunicipal mobility stays at high
levels, comparable with pre-pandemic levels.

Centrality measures evolve over time. How they change seem to be weakly
associated to the evolution of the pandemic. We also observe that not all
nodes exhibit the same behavior. While some nodes, such as Guadalajara,
Puebla-Tlaxcala and Valle de Mexico show larger fluctuations in their degree
centrality, other nodes such as San Pedro Topiltepec or Morelia exhibit much
flatter time series. Changes in betweenness centrality for Valle de Mexico are
very different to those for the other nodes. Evolution for degree and strength
in Valle de Mexico seem to be correlated, but they differ from evolution of
betweenness.

A relatively simple intuition is that municipalities with larger populations
will have higher values of centrality measures regardless of time. We do see
that all of the analysed centrality measures are correlated with population
size. However, we observe that this correlation scales in the same way; Some
of these metrics show a better lineal correlation with population on lineal
scale while some other show a better lineal correlation on logarithmic scale
(see Appendix for model fits).

A related question is whether variations in these centrality measure are also
correlated with population size. We observe a negative correlation between
coefficient of variation and population, showing that centrality variation is
higher in smaller municipalities. Again, the best fit for this correlation is
not necessarily found in the linear scale (see Appendix). In any case, it
should be noted that while population size is clearly a relevant factor, the
dispersion of the point cloud in the range of midsized municipalities indicate
other factors that may be in play to fully explain centrality measures in
these networks.

Regarding formation of communities in the network, we observe that, in a very
broad way, for each day in this analysis there is a “giant” community which
includes about half or more of the nodes in the network and that this
structure seems to preserve for all dates (see Fig. 7. Therefore, if there are
changes in the community structure of inter-municipal mobility, these have to
occur at more local levels. In this diagram lines (nodes in the network) are
colored according to their state. For example, all light-blue lines starting
at the bottom-left corner of the diagram correspond to Oaxaca municipalities,
while green lines just above are Puebla municipalities. Notice here how all
nodes in Puebla tend to stay inside the same community (the largest one), but
municipalities in Oaxaca move to different communities, indicating a change in
the structure of inter-municipal level mobility within this state.

Alluvial plot on Fig. 7 shows community structure and its time evolution for
the hole country. We show in Fig. 8 local versions of these plots for four
different states. Campeche and Nuevo León are among top three states in the
country for Gross Domestic Product per capita (the other being Mexico City,
which is only one node in this analysis, therefore it lacks true intra-state
community structure, since all municipalities are collapsed in the same node);
on the other hand, Oaxaca and-Chiapas are the two states with the lowest GDP
per capita.

## Methods

### Mobile device location data

We used mobile device location data for the time period between 2020-01-01 and
2021-12-31 within Mexican territory provided by Veraset, a company that
aggregates anonymized mobile device location data. This source dataset is
provided as a table in which each record (called a ping) contains the position
(latitude and longitude) of a given (anonymized) device for a given timestamp
(with temporal resolution up to seconds). The set of all unique device ids for
a given day is called the device panel.

### Intermunicipal travel network construction

We define an Intermunicipal Travel Network (IMTN) as a directed, weighted
graph G(V,E), for a given day, such that:

  * •

nodes represent localities (either municipalities or metropolitan zones, as
defined by the national geographic agency; see Note 1)

  * •

links represent mobility from the source node to the target node, defined by
observing at least one device that moved from node i to node j.

  * •

Link weights represent the total fraction of observed devices that moved from
node i to node j, out of the total number of observed devices; this acts as a
normalized measure of flow between nodes (see Note 2).

#### Note 1: on the aggregation of metropolitan area nodes

The political division of Mexico has municipalities as the smallest unit.
Generally, a population center is contained within a municipality; however,
there are large urban areas in which a single population center extends
through many different municipalities, such that the movement between
municipal boundaries is capturing the urban mobility and not travel between
different locations. The National Geography and Statistics Institute (INEGI)
defines 74 metropolitan areas in Mexico, based on measurements from 2015
https://www.inegi.org.mx/contenido/productos/prod_serv/contenidos/espanol/bvinegi/productos/nueva_estruc/702825006792.pdf.
Since intracity mobility is beyond the scope of this manuscript (and would
greatly skew the mobility metrics, as the volume of intra-city mobility is way
larger than that of true travels between different locations), we decided to
aggregate the municipalities that form these metropolitan areas into single
nodes in the network.

#### Note 2: on the interpretation of edge weights in the network

We defined edge weights in the network as follows:

$W_{ij}=\frac{|D_{ij}|}{|D|}$ (1)

where $|D_{ij}|$ is the number of devices that were observed to move from i to
j (that is, were observed in i and their next immediate ping was in j) and
$|D|$ is the total number of devices in the day’s device panel. In this way,
the weight represents a normalized measure of flow between regions; we may
observe that in limit cases, the weight will be zero when there is no movement
observed from one municipality to the other, and the weight would be 1 if all
observed devices within the country travelled from region i to region j (which
would be a virtually impossible scenario). An advantage of using this approach
is that it controls variability in the number of observed devices each day,
allowing for comparisons between days.

#### Network construction algorithm

The following pseudocode shows the the algorithm used for network
construction, starting with the daily mobile device location data table:

Algorithm 1 Intermunicipal Travel Network Construction

1:inputs

2:pings = data.frame(device, timestamp, lat, long)

3:panel = unique device set

4:mun = municipality and metropolitan zone polygon

5:

6:intermediate objects

7:transitions = data.frame(from, to)

8:

9:for ping in pings do

10: pings.mun = locate latitude, longitude in mun

11:end for

12:

13:for k in panel do

14: k.pings = pings[k]

15: sort k.pings by timestamp

16: if k.pings[t].mun != k.pings[t+1].mun then

17: transitions.append(from = k.pings[t].mun, to = k.pings[t+1].mun)

18: end if

19:end for

20:

21:imtn =

22:for i,j in transitions[from,to] do

23: imtn[i,j].weight = $\frac{count(imtn[i,j])}{size(panel)}/$

24:end for

25:

26:RETURN imtn

### Network analysis

Networks were analyzed using the igraph library, version 1.2.7 [8], for the R
programming language, version 4.1.0.

## Data availability

The collection of 731 intermunicipal networks is publicly available on a OSF
repository http://dx.doi.org/10.17605/OSF.IO/42XQZ.

## References

  * [1] A. Ait-Ali and J. Eliasson, The value of additional data for public transport origin–destination matrix estimation, Public Transport, (2021), pp. 1–21.
  * [2] A. A. Alsger, M. Mesbah, L. Ferreira, and H. Safi, Use of smart card fare data to estimate public transport origin–destination matrix, Transportation Research Record, 2535 (2015), pp. 88–96.
  * [3] S. Bai, J. Jiao, and Y. Chen, Mapping the intercounty transmission risk of covid-19 in new york state, Available at SSRN 3582774, (2020).
  * [4] G. Bonaccorsi, F. Pierri, M. Cinelli, A. Flori, A. Galeazzi, F. Porcelli, A. L. Schmidt, C. M. Valensise, A. Scala, W. Quattrociocchi, et al., Economic and social consequences of human mobility restrictions under covid-19, Proceedings of the National Academy of Sciences, 117 (2020), pp. 15530–15535.
  * [5] P. Bonnel, M. Fekih, and Z. Smoreda, Origin-destination estimation using mobile network probe data, Transportation Research Procedia, 32 (2018), pp. 69–81.
  * [6] S. Caschili and A. De Montis, Accessibility and complex network analysis of the us commuting system, Cities, 30 (2013), pp. 4–17.
  * [7] S. Chang, E. Pierson, P. W. Koh, J. Gerardin, B. Redbird, D. Grusky, and J. Leskovec, Mobility network models of covid-19 explain inequities and inform reopening, Nature, 589 (2021), pp. 82–87.
  * [8] G. Csardi and T. Nepusz, The igraph software package for complex network research, InterJournal, Complex Systems (2006), p. 1695.
  * [9] G. de Anda-Jáuregui, P. Guzmán, O. Fontanelli, A. Meneses, A. Hernández, J. de Anda-Gil, M. Flores Garrido, and M. Hernández-Rosales, The contact and mobility networks of mexico city, arXiv e-prints, (2020), pp. arXiv–2007.
  * [10] A. De Montis, A. Chessa, M. Campagna, S. Caschili, and G. Deplano, Modeling commuting systems through a complex network analysis: A study of the italian islands of sardinia and sicily, Journal of Transport and Land Use, 2 (2010), pp. 39–55.
  * [11] H. Deng, J. Du, J. Gao, and Q. Wang, Network percolation reveals adaptive bridges of the mobility network response to covid-19, PloS one, 16 (2021), p. e0258868.
  * [12] T. Djukic, G. Flötteröd, H. Van Lint, and S. Hoogendoorn, Efficient real time od matrix estimation based on principal component analysis, in 2012 15th International IEEE Conference on Intelligent Transportation Systems, IEEE, 2012, pp. 115–121.
  * [13] S. Drobne, A. Garre, E. Hontoria, and M. Konjar, Comparison of two network-theory-based methods for detecting functional regions, Business Systems Research: International journal of the Society for Advancing Innovation and Research in Economy, 11 (2020), pp. 21–35.
  * [14] G. Duranton, Delineating metropolitan areas: Measuring spatial labour market networks through commuting patterns, in The economics of interfirm networks, Springer, 2015, pp. 107–133.
  * [15] P. Edsberg Møllgaard, S. Lehmann, and L. Alessandretti, Understanding components of mobility during the covid-19 pandemic, Philosophical Transactions of the Royal Society A, 380 (2022), p. 20210118.
  * [16] I. Ekowicaksono, F. Bukhari, and A. Aman, Estimating origin-destination matrix of bogor city using gravity model, in IOP Conference Series: Earth and Environmental Science, vol. 31, IOP Publishing, 2016, p. 012021.
  * [17] P. D. Fajgelbaum, A. Khandelwal, W. Kim, C. Mantovani, and E. Schaal, Optimal lockdown in a commuting network, American Economic Review: Insights, 3 (2021), pp. 503–22.
  * [18] R. Frederix, F. Viti, R. Corthout, and C. M. Tampère, New gradient approximation method for dynamic origin–destination matrix estimation on congested networks, Transportation Research Record, 2263 (2011), pp. 19–25.
  * [19] A. Galeazzi, M. Cinelli, G. Bonaccorsi, F. Pierri, A. L. Schmidt, A. Scala, F. Pammolli, and W. Quattrociocchi, Human mobility in response to covid-19 in france, italy and uk, Scientific reports, 11 (2021), pp. 1–10.
  * [20] A. Hadachi, M. Pourmoradnasseri, and K. Khoshkhah, Unveiling large-scale commuting patterns based on mobile phone cellular network data, Journal of Transport Geography, 89 (2020), p. 102871.
  * [21] M. He, J. Glasser, N. Pritchard, S. Bhamidi, and N. Kaza, Demarcating geographic regions using community detection in commuting networks with significant self-loops, PloS one, 15 (2020), p. e0230941.
  * [22] Y. Hong and Y. Yao, Hierarchical community detection and functional area identification with osm roads and complex graph theory, International Journal of Geographical Information Science, 33 (2019), pp. 1569–1587.
  * [23] M. Lenormand, S. Huet, F. Gargiulo, and G. Deffuant, A universal model of commuting networks, (2012).
  * [24] D. Loaiza-Monsalve and A. Riascos, Human mobility in bike-sharing systems: Structure of local and non-local dynamics, PLoS One, 14 (2019), p. e0213106.
  * [25] L. Lotero, R. G. Hurtado, L. M. Floría, and J. Gómez-Gardeñes, Rich do not rise early: spatio-temporal patterns in the mobility networks of different socio-economic classes, Royal Society open science, 3 (2016), p. 150654.
  * [26] H. P. Melo, J. Henriques, R. Carvalho, T. Verma, J. P. Da Cruz, and N. Araújo, Heterogeneous impact of a lockdown on inter-municipality mobility, Physical Review Research, 3 (2021), p. 013032.
  * [27] X. Mu and A. G.-O. Yeh, Regional delineation of china based on commuting flows, Environment and Planning A: Economy and Space, 52 (2020), pp. 478–482.
  * [28] M. A. Munizaga and C. Palma, Estimation of a disaggregate multimodal public transport origin–destination matrix from passive smartcard data from santiago, chile, Transportation Research Part C: Emerging Technologies, 24 (2012), pp. 9–18.
  * [29] E. Nasiboglu, U. Kuvvetli, M. Ozkilcik, and U. Eliiyi, Origin-destination matrix generation using smart card data: Case study for izmir, in 2012 IV International Conference” Problems of Cybernetics and Informatics”(PCI), IEEE, 2012, pp. 1–4.
  * [30] J. Osorio-Arjona and J. C. García-Palomares, Social media and urban mobility: Using twitter to calculate home-work travel matrices, Cities, 89 (2019), pp. 268–280.
  * [31] K. Perrakis, D. Karlis, M. Cools, D. Janssens, K. Vanhoof, and G. Wets, A bayesian approach for modeling origin–destination matrices, Transportation Research Part A: Policy and Practice, 46 (2012), pp. 200–212.
  * [32] A. R. Pitombeira-Neto, C. F. G. Loureiro, and L. E. Carvalho, A dynamic hierarchical bayesian model for the estimation of day-to-day origin-destination flows in transportation networks, Networks and Spatial Economics, 20 (2020), pp. 499–527.
  * [33] N. Pourebrahim, S. Sultana, A. Niakanlahiji, and J.-C. Thill, Trip distribution modeling with twitter data, Computers, Environment and Urban Systems, 77 (2019), p. 101354.
  * [34] G. Pullano, E. Valdano, N. Scarpa, S. Rubrichi, and V. Colizza, Population mobility reductions during covid-19 epidemic in france under lockdown, MedRxiv, 29 (2020), p. 2020.
  * [35] J. J. Ramasco, V. Colizza, and P. Panzarasa, Using the weighted rich-club coefficient to explore traffic organization in mobility networks, in International Conference on Complex Sciences, Springer, 2009, pp. 680–692.
  * [36] A. Riascos and J. L. Mateos, Networks and long-range mobility in cities: A study of more than one billion taxi trips in new york city, Scientific Reports, 10 (2020), pp. 1–14.
  * [37] F. Schlosser, B. F. Maier, O. Jack, D. Hinrichs, A. Zachariae, and D. Brockmann, Covid-19 lockdown induces disease-mitigating structural changes in mobility networks, Proceedings of the National Academy of Sciences, 117 (2020), pp. 32883–32890.
  * [38] C. Seto, A. Khademi, C. Graif, and V. G. Honavar, Commuting network spillovers and covid-19 deaths across us counties, arXiv preprint arXiv:2010.01101, (2020).
  * [39] F. Simini, M. C. González, A. Maritan, and A.-L. Barabási, A universal model for mobility and migration patterns, Nature, 484 (2012), pp. 96–100.
  * [40] Q. Sun, Y. Pan, W. Zhou, C. Xiong, and L. Zhang, Quantifying the influence of inter-county mobility patterns on the covid-19 outbreak in the united states, arXiv preprint arXiv:2006.13860, (2020).
  * [41] T. Toledo and T. Kolechkina, Estimation of dynamic origin–destination matrices using linear assignment matrix approximations, IEEE Transactions on Intelligent Transportation Systems, 14 (2012), pp. 618–626.
  * [42] R. Tolouei, S. Psarras, and R. Prince, Origin-destination trip matrix development: Conventional methods versus mobile phone data, Transportation Research Procedia, 26 (2017), pp. 39–52.  Emerging technologies and models for transport and mobility.
  * [43] N. Wang and P. Mirchandani, Sensor location model to optimize origin–destination estimation with a bayesian statistical procedure, Transportation research record, 2334 (2013), pp. 29–39.
  * [44] Y. Wang, H. Yin, H. Chen, T. Wo, J. Xu, and K. Zheng, Origin-destination matrix prediction via graph convolution: a new perspective of passenger demand modeling, in Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2019, pp. 1227–1235.
  * [45] H. Yilmazkuday, Covid-19 spread and inter-county travel: Daily evidence from the us, Transportation Research Interdisciplinary Perspectives, 8 (2020), p. 100244.


# Liouville equations on complete surfaces with nonnegative Gauss curvature

Xiaohan Cai  School of Mathematical Sciences, Shanghai Jiao Tong University
<EMAIL_ADDRESS>and  Mijia Lai  School of Mathematical Sciences,
Shanghai Jiao Tong University<EMAIL_ADDRESS>
###### Abstract.

We study finite total curvature solutions of the Liouville equation $\Delta
u+e^{2u}=0$ on a complete surface $(M,g)$ with nonnegative Gauss curvature. It
turns out that the asymptotic behavior of the solution separates two extremal
cases: on the one end, if the solution decays not too fast, then $(M,g)$ must
be isometric to the standard Euclidean plane; on the other end, if $(M,g)$ is
isometric to the flat cylinder $S^{1}\times\mathbb{R}$, then solutions must
decay linearly and are completely classified.

Lai’s research is supported by NSFC No. 12031012, No. 12171313 and the
Institute of Modern Analysis-A Frontier Research Center of Shanghai. Cai’s
research is supported by NSFC No.12171313.

##  1\. Introduction

In the seminal work [CL1] of Chen and Li, they obtained the radial symmetry of
the solution of

(1.1) $\displaystyle\Delta u+e^{2u}=0$

on $\mathbb{R}^{2}$, provided that $\int_{\mathbb{R}^{2}}e^{2u}dx<\infty$. Put
the center of symmetry at origin and up to a rescaling, then

$u(x)=\ln\left(\frac{2}{1+|x|^{2}}\right).$

The geometric meaning of above equation is that the conformal metric
$g=e^{2u}g_{0}$ has constant Gauss curvature $1$. It is tempting to think that
$g$ is isometric to the standard round sphere. It is indeed true as the
solution is the pull back of the round metric via stereogrpahic projection.
Nevertheless this line of reasoning is valid only if one establishes the
precise asymptotic behavior of $u$ at $\infty$, so that the metric extends to
a smooth metric on the sphere from $\mathbb{R}^{2}$. The reader is refereed to
[LiTa] for this line of reasoning, see also [GL1] for more general equations.

The assumption $\int_{\mathbb{R}^{2}}e^{2u}dx<\infty$ is natural since there
are infinitely many solutions of (1.1) with
$\int_{\mathbb{R}^{2}}e^{2u}dx=\infty$. One way to obtain such solution is to
pull back the spherical metric via a univalent holomorphic map from
$\mathbb{C}$ to $\overline{\mathbb{C}}$. Recently, there appeared some
interesting studies on (1.1) subject to
$\int_{\mathbb{R}^{2}}e^{2u}dx=\infty$. Eremenko-Gui-Li-Xu [EGLX] give a
complete classification of solutions of (1.1) which are bounded from above. We
also refer to [GL2], [BEL], [L] for some studies on (1.1) from geometric point
of view.

The story in higher dimension was accomplished even earlier. For $n\geq 3$,
let $u$ be a positive solution of

(1.2) $\displaystyle\Delta u+u^{\frac{n+2}{n-2}}=0.$

We refer it as the scalar curvature equation as the conformal metric
$g=u^{\frac{4}{n-2}}g_{0}$ has positive constant scalar curvature. Gidas-Ni-
Nirenberg [GNN] first proved the radial symmetry of the solution under the
assumption $u(x)\sim O(|x|^{2-n})$ as $|x|\to\infty$. This can be viewed as an
analytical proof of a famous result of Obata on classification of constant
scalar curvature metrics which are conformal to an Einstein metric. In a
remarkable paper [CGS], Caffarelli-Gidas-Spruck established the radial
symmetry of the solution without any assumption on the asymptotic behavior of
$u$.

The scalar curvature equation for conformal metrics has critical Sobolev
power. In the subcritical case,

(1.3) $\displaystyle\Delta u+u^{p}=0,\quad 1<p<\frac{n+2}{n-2}.$

Gidas-Spruck [GS] showed any nonnegative solution must be trivial. In a recent
paper [CM], Catino-Monticelli carried out a systematic study of above
mentioned equations (1.1,1.2, 1.3) on complete manifolds with nonnegative
Ricci curvature. Among many results, one particular case is a full extension
of Caffarelli-Gidas-Spruck’s result in dimension three to complete manifolds
with nonnegative Ricci curvature.

Inspired by Catino-Monticelli’s work, we aim to study the Liouville equation
(1.1) on complete surfaces with nonnegative Gauss curvature, we are able to
connect the asymptotic behavior of the solution with the underlying manifold.

To be more precise, let $(M,g)$ be a complete surface with nonnegative Gauss
curvature. We study the Liouville equation

(1.4) $\displaystyle\Delta_{g}u+e^{2u}=0$

on $M$. A solution is called to have finite total curvature if
$\int_{M}e^{2u}dg<\infty$.

In view of the Cohn-Vossen splitting Theorem, a complete surface $(M,g)$ with
nonnegative Gauss curvature is

  * •

either isometric to the flat cylinder $S^{1}\times\mathbb{R}$ or the flat
Möbius band (unorientable),

  * •

or diffeomorphic to $(\mathbb{R}^{2},g_{0})$.

In the latter case, by Huber’s theorem [Hu], $(M,g)$ is conformal to
$(\mathbb{R}^{2},g_{0})$.

Without loss of generality, we assume from now on that $M$ is orientable. In
the former case, we have the following classification of solutions to (1.4).

######  Theorem 1.

Let $u$ be a solution of (1.4) with finite total curvature on the flat
cylinder $(S^{1}\times\mathbb{R},g_{prod})$. Then there exists
$\mu\in[0,\infty)$ and $\beta\in(-1,\infty)$, such that either $\beta$ is an
integer or $\mu=0$, and up to a rescaling, we have

$\displaystyle e^{2u(z)}=\frac{(2\beta+2)^{2}|z|^{2\beta+2}}{(|1+\mu
z^{\beta+1}|^{2}+|z|^{2\beta+2})^{2}}\quad\text{on
}(\mathbb{C}-\\{0\\},\frac{1}{|z|^{2}}g_{0}).$

The classification result is in fact not new. Since the Gauss curvature for
the flat cylinder is identically zero, (1.4) has a geometric meaning that the
conformal metric $e^{2u}g_{prod}$ has Gauss curvature $1$. Note the flat
cylinder is conformal to $(\mathbb{R}^{2}\setminus\\{0\\},g_{0})$, thus (1.4)
can be translated to the Liouville equation on
$\mathbb{R}^{2}\setminus\\{0\\}$. Then the theorem follows from a combination
of results of Chou-Wan [CW, Theorem 5] , Chen-Li [CL2, Theorem 3.1] and
Troyanov [Tr, Theorem II].

Our main theorem is the following rigidity result.

######  Theorem 2.

Let $u$ be a solution of (1.4) with finite total curvature on a complete
surface $(M,g)$ with nonnegative Gauss curvature. Let $r(x)$ be the distance
function on $M$ with respect to a fixed point. If $u(x)\geq-2\ln r(x)+o(\ln
r(x))$, for $r(x)$ large, then $(M,g)$ must be isometric to
$(\mathbb{R}^{2},g_{0})$. Moreover, $-2$ is optimal in the sense that there
exists non flat $(M,g)$ which admits solutions verifying $u(x)\sim\gamma\ln
r(x)$ for any $\gamma<-2.$

A similar result has been proved in [CM, Theorem 1.10]. Our contribution here
has two-fold. On the one hand, our assumption on $u$ is weaker than that in
[CM] and our treatment emphasizes the analysis of asymptotic behavior of the
solution which helps to identify the threshold where the rigidity occurs. On
the other hand, by setting the stage on the complete surfaces with nonnegative
Gauss curvature, we unite two works of Chen-Li [CL1] and [CL2].

As mentioned above, we focus on the asymptotic behavior of the solution. If
$(M,g)$ is conformal to $(\mathbb{R}^{2},g_{0})$, we assume $g=e^{2f}g_{0}$,
then (1.4) becomes

(1.5) $\displaystyle\Delta u+e^{2f}e^{2u}=0\quad\text{on $\mathbb{R}^{2}$.}$

This is the so-called the prescribing Gauss curvature equation on
$\mathbb{R}^{2}$, which has been investigated intensively over past few
decades. Under suitable decay assumption of $e^{2f}$ near infinity, Cheng-Lin
[ChL, Theorem 1.1] showed that the solution $u$ of (1.5) has the following
asymptotic behavior

$\lim_{x\to\infty}\frac{u(x)}{\ln|x|}=-\frac{1}{2\pi}(\int_{\mathbb{R}^{2}}e^{2f}e^{2u}dx)$

if and only if $\int_{\mathbb{R}^{2}}e^{2f}e^{2u}dx<\infty$. However, a
priori, there is no any decay control for $e^{2f}$. In fact, $f$ satisfies an
equation of similar type

$\Delta f+K_{g}e^{2f}=0,$

where $K_{g}$ is the Gauss curvature of $g$. The only information here is that
$K_{g}\geq 0$. Nevertheless, using Arsove-Huber’s result [AH], there exists an
$m\in[0,1]$ and an exceptional set $E$ which is thin at infinity such that

(1.6) $\displaystyle\lim_{x\to\infty,x\notin
E}\frac{u(x)}{\ln|x|}=\liminf_{x\to\infty}\frac{u(x)}{\ln|x|}=-m.$

Here the thinness of a set at infinity is a concept concerning the logarithmic
capacity. For a conformal metric $e^{2f}g_{0}$ on $\mathbb{R}^{n}(n\geq 3)$
with nonnegative Ricci curvature, Ma-Qing [MQ] obtained a similar asymptotic
behavior for the conformal factor $f$.

While Cheng-Lin and Arsove-Huber’s works are main analytical inspirations for
our work, we also benefit from two interesting geometric ingredients: the
first is Li-Tam’s work [LT] on comparison between the intrinsic distance and
the Euclidean distance on $(\mathbb{R}^{2},e^{2f}g_{0})$ (see Lemma 2.2) and
the second is an isoperimetric inequality on surfaces with nonnegative Gauss
curvature established recently by Brendle [B] (see Lemma 2.3).

We present proofs in the next section. The natural analog on the study of
equation (1.2) on higher dimensional complete locally conformally flat
manifolds with nonnegative Ricci curvature will appear in a future work.

Acknowledgement: Both authors wish to thank Prof. Shiguang Ma for helpful
discussions.

##  2\. Proof of Main theorem

###### Proof of Theorem 1.

Note that the flat cylinder $\mathbb{S}^{1}\times\mathbb{R}$ is conformal to
$(\mathbb{R}^{2}\setminus\\{0\\},g_{0})$ since

$\displaystyle
dt^{2}+d\theta^{2}=\frac{1}{r^{2}}dr^{2}+d\theta^{2}=\frac{1}{r^{2}}g_{0},$

by setting $t=\ln r$. Let $e^{2w(x)}=\frac{1}{|x|^{2}}e^{2u(x)}$, then
equation $\Delta_{g}u+e^{2u}=0$ is equivalent to

(2.1) $\displaystyle\begin{cases}\Delta w+e^{2w}=0\quad\mathrm{on}\
\mathbb{R}^{2}\setminus\\{0\\},\\\
\int_{\mathbb{R}^{2}}e^{2w(x)}dx<\infty.\end{cases}$

Chou-Wan’s complex analysis argument [CW, Theorem 5] shows that

$\displaystyle w(x)=\beta_{1}\ln|x|+O(1)\quad\mathrm{as}\ x\to 0,\ \text{for
some}\ \beta_{1}>-1.$

Let $\tilde{w}(x)=w(\frac{x}{|x|^{2}})-2\ln|x|$, it is easy to see that
$\tilde{w}$ satisfies

$\displaystyle\begin{cases}\Delta\tilde{w}+e^{2\tilde{w}}=0\quad\mathrm{on}\
\mathbb{R}^{2}\setminus\\{0\\},\\\
\int_{\mathbb{R}^{2}}e^{2\tilde{w}(x)}dx<\infty.\end{cases}$

Apply Chou-Wan’s asymptotic result [CW, Theorem 5] to $\tilde{w}$ and trace
back to $w$, we get

$\displaystyle w(x)=\beta_{2}\ln|x|+O(1)\quad\mathrm{as}\ x\to\infty,\
\text{for some}\ \ \beta_{2}<-1.$

Therefore, $w(x)$ is a solution of (2.1) with conical singularities at $x=0$
and $x=\infty$. Hence the classification result of Troyanov [Tr, Theorem II]
yields, there exists $\mu\in[0,\infty)$ and $\beta\in(-1,\infty)$ such that
either $\beta$ is an integer or $\mu=0$, and up to a rescaling, we have

$\displaystyle e^{2w(z)}=\frac{(2\beta+2)^{2}|z|^{2\beta}}{(|1+\mu
z^{\beta+1}|^{2}+|z|^{2\beta+2})^{2}}\quad\text{on }\,\mathbb{C}-\\{0\\}.$

Then the desired result follows since $e^{2u(z)}=|z|^{2}e^{2w(z)}$. Note if
both cone angles are less than $2\pi$ ($\beta\in(-1,0)$), Chen-Li [CL2,
Theorem 3.1] also obtained such classification. ∎

In this following, we give the complete proof of Theorem 2.

First we exclude the case of flat cylinder in Theorem 2. Suppose $u$ is a
finite total curvature solution of (1.4) on the flat cylinder, then Theorem 1
implies

$\displaystyle u(x)\sim-(\beta+1)r(x),\quad\text{for $r(x)$ large,}$

where $\beta>-1$ is a constant. This contradicts with the assumption that
$u(x)\geq-2\ln r(x)+o(\ln r(x))$ for $r(x)$ large. In conclusion, $(M,g)$ can
not be the flat cylinder and thus is conformal to $(\mathbb{R}^{2},g_{0})$ by
Cohn-Vossen splitting theorem and Huber’s theorem.

Now we write $g=e^{2f}g_{0}$, then the finite total curvature solution $u$ of
(1.4) becomes

(2.2) $\displaystyle\begin{cases}\Delta u+e^{2f}e^{2u}=0\quad\text{on
$\mathbb{R}^{2}$,}\\\ \int_{\mathbb{R}^{2}}e^{2f+2u}dx<\infty.\end{cases}$

To fix the notation, we consider the quantity

(2.3) $\displaystyle\alpha:=-\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{2f+2u}dx.$

The strategy of our proof is as follows: using the asymptotic lower bound
assumption of the solution $u$, we establish a lower bound of $\alpha$ by
analyzing carefully the asymptotic upper bound of the solution to (2.2). On
the other hand, with the help of Brendle’s isoperimetric inequality, we prove
that the reversed inequality still holds. Hence the equality is obtained and
the rigidity part of the isoperimetric inequality brings the rigidity of the
underlying manifold.

First, we aim at getting the lower bound of $\alpha$. It is tempting to obtain
a pointwise upper bound of the solution $u$ to (2.2) im terms of $\alpha$ so
that the lower bound assumption on $u$ could imply immediately the lower bound
of $\alpha$. However, due to the lack of a uniform asymptotic behavior of the
conformal factor $f$, it’s impossible to derive such a pointwise bound for
$u$. Instead, we shall give an upper bound of the integral average of $u$ on
small balls. The argument is based on that of [ChL].

######  Lemma 2.1.

Let $(M,g)=(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete surface with
nonnegative Gauss curvature. Assume $u\in C^{2}(\mathbb{R}^{2})$ is a solution
to (2.2). Then for any $\epsilon>0$, $\sigma>0$, there exists $R>0$ such that
for $|x|\geq R$ and $\rho=|x|^{-\sigma}$, there holds

$\displaystyle\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}u(y)dy\leq(\alpha+\epsilon)\ln|x|+C,$

where $\alpha$ is given by (2.3) and $C$ is a constant depends on
$\epsilon,\sigma,R$.

###### Proof.

Construct an auxiliary function

$\displaystyle
v(x)=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\psi(y)\ln(\frac{|x-y|}{|y|})dy,$

where $\psi(y)=e^{2f(y)+2u(y)}$.

The proof is mainly composed of three claims:

  1. (1)

$v(x)\leq-\alpha\ln|x|+C$ for $|x|\geq 2$.

  2. (2)

$u+v$ is a constant.

  3. (3)

For any $\epsilon>0,\sigma>0$, there exists $R>0$ such that for $|x|\geq R$
and $\rho=|x|^{-\sigma}$, there holds

(2.4) $\displaystyle
u(x)\leq(\alpha+\epsilon)\ln|x|+\frac{1}{2\pi}\int_{B_{\rho}(x)}\psi(y)\ln(\frac{|y|}{|x-y|})dy+C,$

Proof of claim (1): For fixed $x$ with $|x|\geq 2$,

$\displaystyle 2\pi v(x)$
$\displaystyle=\int_{T_{1}}\psi(y)\ln(\frac{|x-y|}{|y|})dy+\int_{T_{2}}\psi(y)\ln(\frac{|x-y|}{|y|})dy+\int_{T_{3}}\psi(y)\ln(\frac{|x-y|}{|y|})dy$
$\displaystyle\overset{\text{def}}{=}I_{1}+I_{2}+I_{3},$

where $T_{1}=\\{y\\!:|y|\leq
2\\},T_{2}=\\{y\\!:|y-x|\leq\frac{|x|}{2},|y|\\!\geq
2\\},T_{3}=\\{y\\!:|y-x|\geq\frac{|x|}{2},|y|\geq 2\\}$.

Note that for $|x|\geq 2$ and $y\in T_{1}$, we have
$\ln|x-y|\leq\ln(|x|+2)\leq\ln|x|+\ln 2$. Thus

$\displaystyle I_{1}$
$\displaystyle=\int_{T_{1}}\psi(y)\ln|x-y|dy-\int_{T_{1}}\psi(y)\ln|y|dy$
$\displaystyle\leq\int_{T_{1}}\psi(y)(\ln|x|+\ln
2)dy-\int_{T_{1}}\psi(y)\ln|y|dy$
$\displaystyle=(\ln|x|)\int_{T_{1}}\psi(y)dy+C.$

Now for $y\in T_{2}$, we have $|x-y|\leq\frac{|x|}{2}\leq|y|$. Hence

$\displaystyle I_{2}\leq 0.$

For $y\in T_{3}$ and $|x|\geq 2$, there holds $|x-y|\leq|x|+|y|\leq|x||y|$.
Therefore

$\displaystyle I_{3}\leq(\ln|x|)\int_{T_{3}}\psi(y)dy.$

We conclude that

$\displaystyle 2\pi v(x)=I_{1}+I_{2}+I_{3}\leq-2\pi\alpha\ln|x|+C.$

The proof of claim (1) is finished.

Proof of claim (2): It is easy to see that $\Delta v=e^{2f+2u}$ and $u+v$ is a
harmonic function on $\mathbb{R}^{2}$. Hence there exists an entire function
$f(z)$ such that Re$f=2(u+v)$. Let $F(z)=e^{f(z)}$. Clearly, by claim $(1)$ we
get

$\displaystyle|F(z)|=e^{2u+2v}\leq C|z|^{-2\alpha}e^{2u},$

for $|z|\geq 2$. Using the lower bound (1.6) for the conformal factor $f$
($e^{2f}\geq|z|^{-2m}$), we get that for some $R_{0}$ large enough,

$\displaystyle\int_{|z|\geq R_{0}}|F(z)||z|^{2\alpha}|z|^{-2m}dx\leq
C\int_{|z|\geq R_{0}}e^{2u}e^{2f}dx<\infty.$

Let $M(\rho)=\max_{|z|=\rho}|F(z)|$, we shall show that $M(\rho)\leq
C\rho^{2m-2\alpha}$ for $\rho\geq R_{0}+1$. In fact, assume $|z_{0}|=\rho$ and
$M(\rho)=|F(z_{0})|$. The mean value property implies

$\displaystyle|F(z_{0})|\leq\frac{1}{\pi}\int_{B_{1}(z_{0})}|F(z)|dx\leq\frac{1}{\pi}\int_{\rho-1\leq|z|\leq\rho+1}|F(z)|dx.$

Hence we get

$\displaystyle M(\rho)\rho^{2\alpha-2m}$
$\displaystyle\leq\frac{1}{\pi}\int_{\rho-1\leq|z|\leq\rho+1}|F(z)|\rho^{2\alpha-2m}dx$
$\displaystyle\leq\frac{1}{\pi}\int_{\rho-1\leq|z|\leq\rho+1}|F(z)||z|^{2\alpha-2m}(\frac{\rho}{\rho+1})^{2\alpha-2m}dx$
$\displaystyle\leq\frac{2^{2m-2\alpha}}{\pi}\int_{|z|\geq\rho-1}|F(z)||z|^{2\alpha-2m}dx<\infty.$

Therefore, the order of the entire function $F(z)$ is

$\displaystyle\lambda:=\limsup_{\rho\to\infty}\frac{\ln\ln
M(\rho)}{\ln\rho}=0.$

By a theorem of Hadamard (see Theorem 8 of p. 209 in [A]), we conclude that
the genus of $F(z)$ is zero and $F(z)$ is a constant since $F$ has no zeros.
The proof of claim (2) is completed.

Proof of claim (3): For any $\epsilon>0,\sigma>0$, choose $R>0$ large enough
such that

$\displaystyle(\sigma+1)\int_{|y|\geq R}\psi(y)dy\leq\pi\epsilon,$

where $\psi(y)=e^{2f(y)+2u(y)}$. By claim (2), we have

$\displaystyle 2\pi u(x)$
$\displaystyle=C+\int_{\mathbb{R}^{2}}\psi(y)\ln(\frac{|y|}{|x-y|})dy$
$\displaystyle=C+\int_{\tilde{T}_{1}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+\int_{\tilde{T}_{2}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+\int_{\tilde{T}_{3}}\psi(y)\ln(\frac{|y|}{|x-y|})dy$
$\displaystyle\overset{\text{def}}{=}\tilde{I}_{1}+\tilde{I}_{2}+\tilde{I}_{3},$

where $\tilde{T}_{1}=\\{y\\!:|y|\\!\leq
R\\},\tilde{T}_{2}=\\!\\{y\\!:|y-x|\leq\frac{|x|}{2},|y|\\!\geq
R\\},\tilde{T}_{3}=\\!\\{y\\!:|y-x|\geq\frac{|x|}{2},|y|\\!\geq R\\}$.

Now for $|x|\geq\frac{R^{2}}{R-1}$ and $y\in\tilde{T}_{1}$, we have
$\ln|x-y|\geq\ln(|x|-R)\geq\ln|x|-\ln R$. Thus

$\displaystyle\tilde{I}_{1}$
$\displaystyle=\int_{\tilde{T}_{1}}\psi(y)\ln|y|dy-\int_{\tilde{T}_{1}}\psi(y)\ln|x-y|dy$
$\displaystyle\leq C-(\ln|x|)\int_{\tilde{T}_{1}}\psi(y)dy+(\ln
R)\int_{\tilde{T}_{1}}\psi(y)dy$
$\displaystyle\leq-(\ln|x|)\int_{\tilde{T}_{1}}\psi(y)dy+C.$

To estimate $\tilde{I}_{2}$, let
$\tilde{T}^{\sigma}=\\{y:|y-x|\leq|x|^{-\sigma},|y|\geq R\\}$. Then we have

$\displaystyle\tilde{I}_{2}$
$\displaystyle=\int_{\tilde{T}^{\sigma}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+\int_{|x|^{-\sigma}\leq|y-x|\leq\frac{|x|}{2},|y|\geq
R}\psi(y)\ln(\frac{|y|}{|x-y|})dy$
$\displaystyle\leq\int_{\tilde{T}^{\sigma}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+\int_{|y|\geq
R}\psi(y)\ln(\frac{\frac{3}{2}|x|}{|x|^{-\sigma}})dy$
$\displaystyle\leq\int_{|y-x|\leq|x|^{-\sigma}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+(\sigma+1)\int_{|y|\geq
R}\psi(y)dy+C.$

Now for $y\in\tilde{T}_{3}$, one easily gets $|y|\leq 4|x-y|$. Therefore

$\displaystyle\tilde{I}_{3}=\int_{\tilde{T}_{3}}\psi(y)\ln(\frac{|y|}{|x-y|})dy\leq(\ln
4)\int_{\tilde{T}_{3}}\psi(y)dy\leq C.$

In conclusion, there holds

$\displaystyle 2\pi u(x)$
$\displaystyle=\tilde{I}_{1}+\tilde{I}_{2}+\tilde{I}_{3}$ $\displaystyle\leq
C-(\ln|x|)\int_{|y|\leq
R}\psi(y)dy+\int_{|y-x|\leq|x|^{-\sigma}}\psi(y)\ln(\frac{|y|}{|x-y|})dy+(\sigma+1)\int_{|y|\geq
R}\psi(y)dy$ $\displaystyle\leq
C+2\pi(\alpha+\epsilon)\ln|x|+\int_{|y-x|\leq|x|^{-\sigma}}\psi(y)\ln(\frac{|y|}{|x-y|})dy,$

for $|x|\geq\frac{R^{2}}{R-1}$. The proof of claim (3) is completed.

Finally, we give the upper bound of the integral average of $u$. By Green’s
formula, we get

$\displaystyle
u(x)=\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}u(y)dy+\frac{1}{2\pi}\int_{B_{\rho}(x)}\psi(y)\ln(\frac{\rho}{|x-y|})dy,$

for every $x\in\mathbb{R}^{2}$ and $\rho>0$. Combined with (2.4), we have for
any $\epsilon>0,\sigma>0$, there exists $R>0$ such that for $|x|\geq R$ and
$\rho=|x|^{-\sigma}$,

(2.5)
$\displaystyle\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}u(y)dy\leq(\alpha+\epsilon)\ln|x|+\frac{1}{2\pi}\int_{B_{\rho}(x)}\psi(y)\ln(\frac{|y|}{\rho})dy+C.$

Note that
$\frac{|y|}{\rho}\leq\frac{|x|+\rho}{\rho}=|x|^{\sigma+1}+1\leq|x|^{\sigma+2}$
for $|x|$ large enough, the second term in the right hand side of (2.5) could
be estimated as

$\displaystyle\frac{1}{2\pi}\int_{B_{\rho}(x)}\psi(y)\ln(\frac{|y|}{\rho})dy\leq\frac{\sigma+2}{2\pi}(\ln|x|)\int_{|y|\geq
R/2}\psi(y)dy\leq\epsilon\ln|x|,$

for $|x|\geq R$ provided $R$ is large enough. Inserting this into (2.5), the
proof of the lemma is completed. ∎

To derive the lower bound of $\alpha$, we need a lower bound of $u$ in terms
of the Euclidean distance $\ln|x|$ rather than the intrinsic distance $\ln
r(x)$ appeared in the hypotheses of Theorem 2. Fortunately, the comparison of
these two distances is established by Li-Tam [LT, Corollary 3.3]. Moreover,
Hartman [Ha, Theorem 7.1] revealed the connection between this limit and the
asymptotic volume ratio of the manifold. Their results are combined as
follows.

######  Lemma 2.2 (Hartman, Li-Tam).

Let $(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete manifold with nonnegative
Gauss curvature $K$. Then

$\displaystyle\lim_{x\to\infty}\frac{\ln
r(x)}{\ln|x|}=1-\frac{1}{2\pi}\int_{\mathbb{R}^{2}}Kdg=\beta,$

where

$\displaystyle\beta:=\lim_{t\to\infty}\frac{Area(B(p,t))}{\pi t^{2}}\in[0,1]$

is the asymptotic volume ratio of the manifold $(\mathbb{R}^{2},e^{2f}g_{0})$.

Given this asymptotic behavior of $r(x)$, the a prior assumption on $u$ could
be applied to obtain the lower bound of $\alpha$ in terms of the asymptotic
volume ratio.

######  Proposition 2.1.

Let $(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete surface with nonnegative
Gauss curvature. Let $u$ be a solution of (2.2). Assume

$\displaystyle u(x)\geq-2\ln r(x)+o(\ln r(x)),$

for $r(x)$ large. Then

$\displaystyle\alpha\geq-2\beta,$

where $\alpha$ is given by (2.3) and $\beta$ is the asymptotic volume ratio of
$(\mathbb{R}^{2},e^{2f}g_{0})$.

###### Proof.

By our assumption on $u$ and Lemma 2.2, we get for any $\epsilon>0$, there
exits $R>0$ such that for $r(x)\geq R$,

$\displaystyle u(x)\geq-2\ln r(x)+o(\ln
r(x))\geq(-2\beta-2\epsilon)\ln|x|+o(\ln|x|).$

While Lemma 2.1 yields that for any $\epsilon>0,\sigma>0$, there exists $R>0$
such that for $|x|\geq R$ and $\rho=|x|^{-\sigma}$,

$\displaystyle\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}u(y)dy\leq(\alpha+\epsilon)\ln|x|+C,$

where $C$ is a constant depends on $\epsilon,\sigma,R$.

We conclude that for any $\epsilon>0,\sigma>0$, there exists $R>0$ such that
for $|x|\geq R$ and $\rho=|x|^{-\sigma}$,

$\displaystyle(\alpha+\epsilon)\ln|x|+C$
$\displaystyle\geq\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}u(y)dy$
$\displaystyle\geq(-2\beta-2\epsilon)\frac{1}{\pi\rho^{2}}\int_{B_{\rho}(x)}\ln|y|dy+o(\ln|x|)$
$\displaystyle\geq(-2\beta-2\epsilon)(\ln|x|-\epsilon)+o(\ln|x|).$

Let $x\to\infty$, we get $\alpha+\epsilon\geq-2\beta-2\epsilon$. Since
$\epsilon$ could be arbitrarily small, we get

$\displaystyle\alpha\geq-2\beta.$

∎

We shall see that the reversed inequality also holds and thus the equality is
obtained. For this, we need the isoperimetric inequality on nonnegatively
curved surfaces established by Brendle [B, Corollary 1.3], and it also helps
to get the rigidity of the underlying manifold in our setting.

######  Lemma 2.3 (Brendle).

Let $(M^{2},g)$ be a complete noncompact manifold with nonnegative Gauss
curvature. Let $D$ be a compact domain in $M$ with boundary $\partial D$. Then

$\displaystyle L(\partial D)^{2}\geq 4\pi\beta\ A(D),$

where $L(\partial D)$ and $A(D)$ represent the length of $\partial D$ and the
area of $D$, respectively, and $\beta$ is the asymptotic volume ratio of
$(M,g)$. Moreover, the equality holds if and only if $(M,g)$ is isometric to
Euclidean space and $D$ is a ball.

Now with the help of Lemma 2.3, one could mimic the argument in [CL1] to give
the upper bound of $\alpha$.

######  Proposition 2.2.

Let $(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete surface with nonnegative
Gauss curvature. Let $u$ be a solution of (2.2). Then

$\displaystyle\alpha\leq-2\beta,$

where $\alpha$ is given by (2.3) and $\beta$ is the asymptotic volume ratio of
$(\mathbb{R}^{2},e^{2f}g_{0})$.

###### Proof.

Consider $F(t):=\int_{\Omega_{t}}e^{2u}dg$, where $\Omega_{t}=\\{x:u(x)>t\\}$
is the upper level set of $u$.

The finite total curvature assumption $\int_{M}e^{2u}dg<\infty$ and the
Minkowski inequality yield

$\displaystyle A(\Omega_{t})<\infty,$

where $A(\Omega_{t})$ represents the area of $\Omega_{t}$ in
$(\mathbb{R}^{2},g=e^{2f}g_{0})$.

It follows from the equation (1.4) and the divergence theorem that

$\displaystyle F(t)=\int_{\Omega_{t}}e^{2u}dg=-\int_{\Omega_{t}}\Delta
udg=-\int_{\partial\Omega_{t}}<\nabla
u,\eta>dS_{g}=\int_{\partial\Omega_{t}}|\nabla u|dS_{g}.$

By the coarea formula, there holds

$\displaystyle F^{\prime}(t)=-\int_{\partial\Omega_{t}}\frac{e^{2u}}{|\nabla
u|}dS_{g}=-e^{2t}\int_{\partial\Omega_{t}}\frac{1}{|\nabla u|}dS_{g}.$

Then the Hölder inequality and the isoperimetric inequality (Lemma 2.3) imply

(2.6)
$\displaystyle\begin{split}(F^{2}(t))^{\prime}&=-2e^{2t}\int_{\partial\Omega_{t}}|\nabla
u|dS_{g}\int_{\partial\Omega_{t}}\frac{1}{|\nabla u|}dS_{g}\\\
&\leq-2e^{2t}L(\partial\Omega_{t})^{2}\\\ &\leq-8\pi\beta\
e^{2t}A(\Omega_{t}).\end{split}$

Note that the isoperimetric inequality still holds for noncompact regions
whose area are finite, since the length of its boundary must be infinite by
the completeness of $(\mathbb{R}^{2},e^{2f}g_{0})$.

Finally integrating $(\ref{ineq.ode})$ from $-\infty$ to $\infty$ yields

$\displaystyle-(\int_{M}e^{2u}dg)^{2}$
$\displaystyle\leq-8\pi\beta\int_{-\infty}^{\infty}e^{2t}A(\\{x:e^{2u(x)}>e^{2t}\\})dt$
$\displaystyle=-4\pi\beta\int_{0}^{\infty}A(\\{x:e^{2u(x)}>\lambda\\})d\lambda$
$\displaystyle=-4\pi\beta\int_{M}e^{2u}dg.$

Thus the desired inequality holds. ∎

###### Proof of Theorem 2.

By Proposition 2.1 and Proposition 2.2, we get

$\displaystyle\alpha=-2\beta.$

Inspecting the proof of Proposition 2.2 shows that
$L(\partial\Omega_{t})^{2}=4\pi\beta\ A(\Omega_{t})$ for every
$t\in\mathbb{R}$. Hence Lemma 2.3 tells $(\mathbb{R}^{2},e^{2f}g_{0})$ must be
isometric to the Euclidean space $(\mathbb{R}^{2},g_{0})$.

To see the sharpness of the coefficient $-2$ in the assumption $u(x)\geq-2\ln
r(x)+o(\ln r(x))$, consider the following conformal flat manifolds:

Choose the conformal factor
$e^{2f(x)}=\frac{\gamma}{(1+|x|^{2})^{2-2\gamma}}.$ Then for
$\gamma\in[\frac{1}{2},1)$, $(\mathbb{R}^{2},g=e^{2f}g_{0})$ is a complete
surface with nonnegative Gauss curvature
$K_{g}=\frac{4(1-\gamma)}{\gamma(1+|x|^{2})^{2\gamma}}$.

Take $e^{2u(x)}=\frac{4}{(1+|x|^{2})^{2\gamma}}$, it is easy to see that
$\Delta u+e^{2f}e^{2u}=0$. In other words,

$\displaystyle\Delta_{g}u+e^{2u}=0.$

Moreover,
$\int_{\mathbb{R}^{2}}e^{2u}dg=\int_{\mathbb{R}^{2}}\frac{4\gamma}{(1+|x|^{2})^{2}}dx=4\pi\gamma<\infty.$

Direct computation shows

$\displaystyle\lim_{x\to\infty}\frac{\ln
r(x)}{\ln|x|}=2\gamma-1,\quad\mathrm{for}\ \gamma\in(\frac{1}{2},1).$
$\displaystyle\lim_{x\to\infty}\frac{r(x)}{\ln|x|}=1,\quad\mathrm{for}\
\gamma=\frac{1}{2}$

Thus for $\gamma\in(\frac{1}{2},1)$, we have

$\displaystyle u(x)\sim-2\gamma\ln|x|\sim-\frac{2\gamma}{2\gamma-1}\ln r(x),$

where $-\frac{2\gamma}{2\gamma-1}\in(-\infty,-2)$.

In conclusion, for any $k<-2$, there exists a complete surface
$(\mathbb{R}^{2},e^{2f}g_{0})$ with nonnegative Gauss curvature which admits a
finite total curvature solution $u$ of (1.4) with $u(x)\sim k\ln r(x)$.

∎

We conclude this paper with following remark. The main theorem states the
rigidity of the underlying manifold under the assumption $u(x)\geq-2\ln
r(x)+o(\ln r(x))$. However, on the other end, we cannot expect such rigidity
as illustrated by examples above. More precisely, when $\gamma=\frac{1}{2}$,
it readily follows that the solution $u$ decays linearly with respect to the
distance induced by the metric. Hence one cannot distinguish the flat cylinder
by linear decay condition on the solution.

Nevertheless, when the solution decays sufficiently fast, we can get the
volume growth control of the underlying manifold. We record here as a result
of independent interest.

######  Proposition 2.3.

Let $(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete surface with nonnegative
Gauss curvature. Let $u$ be a solution of $(\ref{Efinite})$ satisfying

$\displaystyle\liminf_{x\to\infty}\frac{u(x)}{\ln r(x)}=-\infty.$

Then the asymptotic volume ratio of $(\mathbb{R}^{2},e^{2f}g_{0})$ is zero.

###### Proof.

Suppose the asymptotic volume ratio $\beta$ is positive. According to the
claims (1) and (2) in Lemma 2.1, there holds

$\displaystyle u(x)\geq\alpha\ln|x|+C,\quad\mathrm{for}\ |x|\geq 2,$

where $\alpha=-\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{2f}e^{2u}dx$. Combined
with Lemma 2.2 one get

$\displaystyle\liminf_{x\to\infty}\frac{u(x)}{\ln
r(x)}=\liminf_{x\to\infty}\frac{u(x)}{\ln|x|}\frac{\ln|x|}{\ln
r(x)}\geq\frac{\alpha}{\beta}>-\infty.$

This contradicts to the hypothesis. Hence we get $\beta=0$. ∎

We also have a partial converse to Proposition 2.3.

######  Proposition 2.4.

Let $(\mathbb{R}^{2},e^{2f}g_{0})$ be a complete surface with nonnegative and
bounded Gauss curvature. Suppose the asymptotic volume ratio $\beta=0$, then
there exists a solution of (2.2) satisfying

$\displaystyle\lim_{x\to\infty}\frac{u(x)}{\ln r(x)}=-\infty.$

###### Proof.

Since $f$ satisfies

$\Delta f+Ke^{2f}=0,$

where $0\leq K\leq C$ by assumption. Based on a work of Taliaferro [Ta],
Bonini-Ma-Qing [BMQ, Lemma 4.2] showed that

$e^{2f}\sim|x|^{-2(1-\beta)}=|x|^{-2}\quad\text{as $|x|\to\infty$}.$

In view of the existence theorem of McOwen [Mc, Theorem 1], for any
$\alpha\in(-2,0)$, there exists a solution $u$ of (2.2) verifying

$u(x)\sim\alpha\ln|x|+O(1)\quad\text{at $\infty$}.$

Since Lemma 2.2 still holds for $\beta=0$, the conclusion readily follows.

∎

## References

  * [A] L. Ahlfors, Complex analysis—an introduction to the theory of analytic functions of one complex variable, 3rd edn. McGraw-Hill Book Co., New York, 1978. xi+331 pp.
  * [AH] M. Arsove and A. Huber, Local behavior of subharmonic functions, Indiana Univ. Math. J.22 (1972/73/1973), 1191-1199.
  * [BEL] W. Bergweiler, A. Eremenko and J. Langley, On conformal metrics of constant positive curvature in the plane, J. Math. Phys. Anal. Geom. 19 (2023), no. 1, 59-73
  * [BMQ] V. Bonini, S. Ma and J. Qing, On nonnegatively curved hypersurfaces in $\mathbb{H}^{n+1}$. Math. Ann. 372(2018), no.3-4, 1103-1120.
  * [B] S. Brendle, Sobolev inequalities in manifolds with nonnegative curvature, Comm. Pure Appl. Math. 76 (2023), no. 9, 2192–2218.
  * [CGS] L. Caffarelli, B. Gidas and J. Spruck, Asymptotic symmetry and local behavior of semilinear elliptic equations with critical Sobolev growth, Comm. Pure Appl. Math. 42 (1989), no. 3, 271–297.
  * [CM] G. Catino and D. Monticelli, Semilinear elliptic equations on manifolds with nonnegative Ricci curvature, arXiv: 2203.03345.
  * [CL1] W. Chen and C. Li, Classification of solutions of some nonlinear elliptic equations, Duke Math. J. 63 (1991), no. 3, 615–622.
  * [CL2] W. Chen and C. Li, What kinds of singular surfaces can admit constant curvature?, Duke Math. J. 78 (1995), no. 2, 437–451.
  * [ChL] K.S. Cheng and C.S. Lin, On the asymptotic behavior of solutions of the conformal Gaussian curvature equations in $\mathbb{R}^{2}$, Math. Ann. 308 (1997), no. 1, 119–139.
  * [CW] K. Chou and T. Wan, Asymptotic radial symmetry for solutions of $\Delta u+e^{u}=0$ in a punctured disc, Pacific J. Math. 163 (1994), 269-276.
  * [EGLX] A. Eremenko, C. Gui, Q. Li, and L. Xu, Rigidity results on Liouville equation, arXiv: 2207.05587.
  * [GNN] B. Gidas, W. M. Ni and L. Nirenberg, Symmetry of positive solutions of nonlinear elliptic equations in $\mathbb{R}^{n}$, Mathematical analysis and applications, Part A, pp. 369-402, Adv. in Math. Suppl. Stud., 7a, Academic Press, New York-London, 1981.
  * [GS] B. Gidas, J. Spruck, Global and local behavior of positive solutions of nonlinear elliptic equations, Comm. Pure Appl. Math. 34 (1981), no. 4, 525-598.
  * [GL1] C. Gui and Q. Li, Completion of $\mathbb{R}^{2}$ with a conformal metric as a closed surface, Anal. Theory Appl. 37 (2021), no. 1, 59–73.
  * [GL2] C. Gui and Q. Li, Some geometric inequalities related to Liouville equation, arXiv: 2208.03612.
  * [Ha] P. Hartman, Geodesic parallel coordinates in the large. Amer. J. Math.86(1964), 705–727.
  * [Hu] A. Huber, On subharmonic functions and differential geometry in the large, Comment. Math. Helv. 32, 13-72 (1957)
  * [LT] P. Li and L.F. Tam, Complete surfaces with finite total curvature, J. Differential Geom. 33 (1991), no. 1, 139–168.
  * [LiTa] Y. Li and H. Tang, Remark on Gauss curvature equations on punctured disk, Front. Math. China 15(2020), no.4, 701-707.
  * [L] A.Lytchak, On conformal planes of finite area, arXiv:2211.01347
  * [MQ] S. Ma and J. Qing, On $n$-superharmonic functions and some geometric applications Calc. Var. Partial Differential Equations 60 (2021), no. 6, Paper No. 234, 42 pp.
  * [Mc] R. McOwen, Conformal metrics in $\mathbb{R}^{2}$ with prescribed Gaussian curvature and positive total curvature, Indiana Univ. Math. J.34(1985), no.1, 97–104.
  * [Ta] S. Taliaferro, On the growth of superharmonic functions near an isolated singularity I, J. Differ. Equ. 158, 28-47 (1999)
  * [Tr] M. Troyanov, Metrics of constant curvature on a sphere with two conical singularities, in Differential Geometry, Lecture Notes in Math. 1410, Springer-Verlag, Berlin, 1989, 296-306.


April 17, 2019

# Characterization of SiPM Optical Crosstalk and its Dependence on the
Protection-Window Thickness

Yuki Nakamura Akira Okumura Hiroyasu Tajima Nobuhito Yamane and Anatolii Zenin
Institute for Space–Earth Environmental ResearchInstitute for Space–Earth
Environmental Research Nagoya University Nagoya University Furo-cho Furo-cho
Chikusa-ku Chikusa-ku Nagoya 464-8601 Nagoya 464-8601 Japan Japan
<EMAIL_ADDRESS>(Y. N<EMAIL_ADDRESS>(A. O.)

###### Abstract

Owing to their high photon detection efficiency, compactness, and low
operating voltage, silicon photomultipliers (SiPMs) have found widespread
application in many fields, including medical imaging, particle physics, and
high-energy astrophysics. However, the so-called optical crosstalk (OCT)
phenomenon of SiPMs is a major drawback to their adoption. Secondary infrared
photons are emitted inside the silicon substrate spontaneously after the
avalanche process caused by the primary incident photons, and they can be
detected by the surrounding photodiodes. As a result large output pulses that
are equivalent to multiple photoelectrons are observed with a certain
probability (OCT rate), even for single-photon events, making the charge
resolution worse and increasing the rate of accidental triggers by single-
photon events in applications such as atmospheric Cherenkov telescopes. In our
previous study, we found that the OCT rates of single-channel SiPMs was
dependent on the thickness of their protection resin window, which may be
explained by photon propagation inside the resin. In the present study, we
measured the OCT rate of a multichannel SiPM and those of neighboring channels
caused by photon propagation. Both OCT rates were found to be dependent on the
protection-window thickness. We report our OCT measurements of a multichannel
SiPM and comparisons with a ray-tracing simulation.

Silicon photomultipliers (SiPMs), ray-tracing simulation, Cherenkov telescopes

##  1 Introduction

Silicon photomultipliers (SiPMs) are a gathering of multiple avalanche
photodiode cells that operate in a Geiger mode (G-APD). While a single G-APD
cell produces a constant charge signal that does not depend on the number of
photons hitting the cell, SiPMs are capable of counting the number of incident
photons by arranging multiple G-APD cells in parallel. Nowadays SiPMs are used
in a variety of applications, such as in medical imaging, particle physics,
and high-energy astrophysics owing to their robustness, low cost per channel,
low operating voltage, and high tolerance to bright illumination. However,
SiPMs have a detrimental phenomenon, the so-called optical crosstalk (OCT),
that can negatively affect the system performance in some applications.

OCT is a phenomenon whereby the avalanche process in a single G-APD cell emits
secondary infrared photons that are then detected by the surrounding cells
with a certain probability (OCT rate). As a result, the output chage from the
SiPM channel is not constant, even for a single-photoelectron (1-p.e.) event,
and it is proportional to one plus the number of detected secondary photons.
For very-high-energy ground-based gamma-ray telescopes (also referred to as
Cherenkov telescopes), which detect UV–visible photons from Cherenkov
radiation in atmospheric showers, OCT can degrade the energy resolution of the
primary particles, because the number of Cherenkov photons is approximately
proportional to the energy. Furthermore, OCT can increase the rate of
accidental triggers, which are caused by photons randomly coming from the
night sky. This is because a single photon may generate a large output pulse
equivalent to multiple photons, which cannot be distinguished from temporally
correlated Cherenkov photons. Thus the level-1 trigger threshold needs to be
set higher than preferred, resulting in a higher energy threshold and lower
detection sensitivity to celestial gamma-ray sources.

The use of SiPMs as the focal-plane camera pixels of Cherenkov telescopes was
first realized by the FACT telescope[2]. As of 2019, several SiPM camera
prototypes are also being developed for the Cherenkov Telescope Array (CTA)[3]
with the aim to build compact and less expensive cameras with thousands of
pixels that can have a high tolerance to bright moon nights.

In our previous study, in which an extensive characterization of various SiPM
products was conducted for the CTA, we found that the OCT rate of single-
channel SiPMs had a coating-thickness dependence, as shown in Fig. 1[5]. This
dependence may be explained by photon propagation in the protection window, as
discussed in the paper. However, a firm conclusion could not be derived
because the OCT rate comparison was made for different products and resin
materials, and because escaping photons (illustrated as path E in Fig. 2)
could not be confirmed with single-channel SiPMs.

Figure 1: The OCT rates of five different SiPM products as functions of over
voltage (defined to be the supplied bias voltage minus the breakdown voltage).
The thicker the protection window, the lower OCT rate; but, those with no
protection coating exhibit the lowest OCT rate. The figure is taken from [5]
(reproduced under CC BY 4.0, see reference for details). Figure 2: A schematic
diagram of possible OCT photon paths inside a multichannel SiPM. After the
avalanche process is triggered by the first incident photon (solid-line
arrow), secondary photons may propagate and be detected by surrounding cells
via several paths. Path A is due to Fresnel reflection at the boundary between
the silicon substrate and the air. Path B is a direct path to an adjacent
G-APD cell but mostly blocked by a trench. Path C is by Fresnel reflection at
the boundary between the protection window and the air. Path D is similar to
Path C, but the photon is reflected back to the first G-APD cell, which is
already saturated. Path E is also similar to Path C, but it may be detected by
neighboring channels in the case of multichannel SiPMs.

Fig. 2 illustrates the possible paths of secondary photons propagating in the
protection window and in the silicon substrate. It supports our photon
propagation scenario and explains Fig. 1 as follows.

Secondary photons that directly enter the neighboring cells (path B) can be
blocked by the so-called trench structure developed in recent SiPM products
[1]. When the protection window is very thin or does not exist, the number of
photons propagating in path D increases (i.e., path C or E does not occur),
and the OCT rate decreases because the first G-APD cell is already saturated.
By increasing the protection window thickness, we can make the number of
photons propagating in path E larger than that in path C, and thus the OCT
rate of a single-channel SiPM is reduced. However, if this scenario is
correct, path E may generate OCT pulses in neighboring channels for
multichannel SiPMs that are planned to be used in CTA cameras.

To characterize and understand the OCT rate of multichannel SiPMs and to
confirm the possible photon paths in Fig. 2, we measured the OCT rate of
multichannel SiPM S13361-6805 ($8\times 8$ channels, $3\times 3\
\mathrm{mm}^{2}$ each, 50-$\mu$m cells, no protection coating, Hamamatsu
Photonics). To eliminate the possible systematic uncertainty due to different
resin materials, we used thin glass plates and optical grease to emulate the
different protection-window thickness.

In Section 2, we show the OCT rate for different (glass) window thickness to
confirm the results of the previous study. In Section 3, we present OCT events
occurring in neighboring channels and their rates for the first time. In
Section 4, we also compare our measurements with a simple ray-tracing
simulation and provide a conclusion in Section 5.

##  2 Optical Crosstalk Rate of a Single SiPM Channel

Figure 3: Example SiPM waveforms recorded with TARGET 7 chips: (a) 1-p.e.
pulse in the reference channel is seen at around $155$ ns. (b) Same as (a) but
for Adjacent Channel 1. The counts in the reference channel and this channel
are not coincident. (c) Same as (a) but for Adjacent Channel 2. This may be
temporally correlated with the pulse in the reference channel.

The OCT rate of a particular single channel (hereafter referred to as the
reference channel, see Fig. 7) can be measured by taking the charge
distribution of the channel output in a certain time window. Deriving the
ratio between the numbers of 1-p.e. and 2-p.e. events under an assumption that
the 2-p.e. events do not follow a pure Poisson distribution, the excess number
of 2-p.e. events is considered to be due to OCT. Output waveforms were
recorded with a 1 GSa/s sampling circuit named TARGET 7, a member of TARGET
family developed for CTA [6, 7]. A single TARGET 7 chip is able to record 16
input channels. To record all the 64 output channels simultaneously, a TARGET
7 camera module, equipped with four TARGET 7 chips, was used for this study.
Example waveforms can be seen in Fig 3.

Fig. 4 shows a charge distribution of the reference channel, to which
multiGaussian are fitted to derive the OCT rate,

$R_{\mathrm{OCT}}=1-\frac{1-\frac{N(\geq 2~{}\mathrm{p.e.})}{N(\geq
1~{}\mathrm{p.e.})}}{\exp{(-f_{\mathrm{DCR}}\Delta t_{\mathrm{ps}})}},$ (1)

where $N(\geq 2~{}\mathrm{p.e.})$ and $N(\geq 1~{}\mathrm{p.e.})$ are the
numbers of events with a charge equivalent to 2 p.e. or more and 1 p.e. or
more, respectively, $f_{\mathrm{DCR}}$ is the dark count rate and $\Delta
t_{\mathrm{ps}}$ is the maximum time window in which two accidental dark
counts cannot be temporally resolved in our deconvolution analysis[8].

Figure 4: The charge distribution of the reference channel (solid black line),
the best-fit multiGaussian (dashed green lines), and their sum (solid red
line). Five obvious peaks can be seen from left to right: pedestal, 1 p.e., 2
p.e., and 3 p.e.. Systematic uncertainty in calculating the OCT rate due to
the non-Gaussianity in this distribution is shown with vertical brackets in
Fig. 5.

The OCT rate of the reference channel with various dummy protection-window
thickness is shown in Fig 5 (solid red squares). It is at the minimum and at
the maximum when the thickness is $0$ $\mu$m (i.e., no glass plate is put) and
$\sim 100$ $\mu$m, respectively, and gradually decreases as the window becomes
thicker. This confirms the similar trend found in the previous study (see Fig.
1).

Figure 5: (a) OCT rate vs. window thickness. The OCT rate in the reference
channel (OCT self), that in neighboring channels at “$+$” and “$\times$”
directions (OCT+ and OCT×), and the total. Ray-tracing simulations are shown
with solid curves. Open symbols show the OCT rates of the same SiPM product
with a 100-$\mu$m thick silicone resin to verify our window emulation method
using thin glass plates. (b) Definition of “$+$” and “$\times$” directions.

##  3 Optical Crosstalk to Neighboring Channels

We measured the OCT rates in neighboring channels to test the OCT photon
propagation pathways depicted in Fig. 2. If OCT propagation to neighboring
channels according to path E exists, coincident dark counts may be observed in
the adjacent channels, as shown in Figs. 3(a) and 3(c). Assuming that the
first dark count happened in the reference channel, the 1-p.e. pulse in
Adjacent Channel 2 is considered to be another dark count in the latter
channel or an OCT photon detected right after the dark count in the reference
channel. For individual pulses, we cannot determine which of the two scenarios
is correct. On the contrary, for the pulse in Fig. 3(b), it is not temporally
correlated with the other two pulses and can be considered to be a dark count
occurring in Adjacent Channel 1.

The two possibilities can be later separated by making distributions of the
time difference of photoelectron detection timings ($\Delta t$ ($\equiv
t_{0}-t_{1}$) distribution, $t_{0}$ is the detection timing in the reference
channel, $t_{1}$ is in a neighboring channel). Fig. 6 shows two $\Delta t$
distributions: The red distribution was made from only waveforms in which a
photoelectron was detected in a certain time window in the reference channel.
The black distribution was made from waveforms without any pulse in the
reference channel and $t_{0}$ was randomly chosen in $\pm 100$ ns, resulting
in an exponential-like distribution.

The excess around $\Delta t=0$ in the red histogram was caused by OCT photons
propagating to the neighboring channels. Integrating this excess within $\pm
32$ ns (wider than the timing resolution in our peak search), the OCT rate in
the neighboring channel, which originates from dark counts in the reference
channel, can be calculated.

Figure 6: $\Delta t$ distributions of the reference channel and a neighboring
channel. The red histogram shows $\Delta t$ for the two closest pulses
appeared in the reference and neighboring channels. The black histogram was
made from waveforms without any pulse in the reference channel ($t_{0}$ was
randomly chosen in $\pm 100$ ns).

Figure 7: (a) $8\times 8$ channel SiPM array (Hamamatsu Photonics
S13361-6805). (b) An example 2D distribution of neighboring-channel OCT rates.
The blank 15 channels were excluded in the analysis to avoid common electrical
noise of the TARGET ASIC connected to the reference channel.

Fig. 7 shows the $8\times 8$ channel SiPM and a 2D distribution of
neighboring-channel OCT rates originating from dark counts in the reference
channel. Repeating measurements with different window thicknesses, we obtained
the window-thickness dependence of the OCT rates (Fig 5). Neighboring-channel
OCT occurs in adjacent channels with non-negligible rates ($>\sim 1$%) when
the window is $\sim 100$ $\mu$m or thicker and the total OCT rate summed from
the reference and adjacent channels is roughly constant in this thickness
region. These findings confirmed that OCT photons indeed propagate to
neighboring channels due to reflection at the media boundary between the
window and the ambient air.

##  4 Ray-tracing Simulation

The measurement results shown in Fig. 5 and OCT photon paths shown in Fig. 2
were also verified by ray-tracing simulations. We used the ROBAST library[9]
to build the multichannel SiPM and window geometries, inside which the OCT
photon propagation was simulated assuming a simple Fresnel reflection at the
media boundaries and isotropic and uniform OCT photon emission from all the
APD cells in the reference channel. As shown in Fig. 5, the simulation results
after rough normalization by eye are consistent with the measurements.

##  5 Conclusion

We measured the neighboring-channel OCT rates of a multichannel SiPM and its
window-thickness dependence for the first time. The “self OCT” rate was the
highest at around $100$ $\mu$m and it can be decreased by removing the window
or by increasing the thickness of the window. The neighboring-channel OCT
rates increased as the window thickness increased and neighboring-channel OCT
could only be reduced by removing the window entirely. Thus, removing the
protection window is the preferred way to significantly suppress the total OCT
rate in multichannel SiPM applications like CTA unless the SiPM surface
handling is problematic.

## References

  * [1] D. McNally, V. Golovin, Nucl. Instrum. Methods Phys. Res., Sect. A 610, 150–153 (2009).
  * [2] H. Anderhub, et al., J. Instrum. 8 P06008 (2013).
  * [3] B. S. Acharya, et al., Astropart. Phys. 40, 3–18, (2013).
  * [4] Hamamatsu Photonics. in Opto-Semiconductor handbook, Chap. 3.
https://www.hamamatsu.com/resources/pdf/ssd/03_handbook.pdf

  * [5] A. Asano, et al., Nucl. Instrum. Methods Phys. Res., Sect. A 912, 177–181 (2018).
  * [6] K. Bechtol, et al., Astropart. Phys. 36, 156–165 (2012).
  * [7] A. Albert, et al., Astropart. Phys. 92 49–61 (2017).
  * [8] H. Otono, et al., Nucl. Instrum. Methods Phys. Res., Sect. A, 610, 397–399 (2009).
  * [9] A. Okumura, K. Noda, and C. Rulten, Astropart. Phys. 76 38–47 (2016).