# So2Sat POP - A Curated Benchmark Data Set for Population Estimation from
Space on a Continental Scale

Sugandha Doda
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
Yuanyuan Wang
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
Matthias Kahl
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
Eike Jens Hoffmann
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
Kim Ouan
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
Hannes Taubenböck
German Remote Sensing Data Center,
German Aerospace Center,
Münchener Straße 20, Weßling, 82234, Germany
Institute for Geography and Geology,
Julius-Maximilians-Universitẗ Würzburg,
Würzburg 97074, Germany Xiao Xiang Zhu
Chair of Data Science in Earth Observation,
Department of Aerospace and Geodesy,
Technical University of Munich,
Arcisstraße 21, Munich,80333, Germany
corresponding author: Xiao Xiang Zhu<EMAIL_ADDRESS>
###### Abstract

Obtaining a dynamic population distribution is key to many decision-making
processes such as urban planning, disaster management and most importantly
helping the government to better allocate socio-technical supply. For the
aspiration of these objectives, good population data is essential. The
traditional method of collecting population data through the census is
expensive and tedious. In recent years, statistical and machine learning
methods have been developed to estimate population distribution. Most of the
methods use data sets that are either developed on a small scale or not
publicly available yet. Thus, the development and evaluation of new methods
become challenging. We fill this gap by providing a comprehensive data set for
population estimation in 98 European cities. The data set comprises a digital
elevation model, local climate zone, land use proportions, nighttime lights in
combination with multi-spectral Sentinel-2 imagery, and data from the Open
Street Map initiative. We anticipate that it would be a valuable addition to
the research community for the development of sophisticated approaches in the
field of population estimation.

##  1 Background & Summary

Rapid urbanization in cities is leading to environmental concerns such as
climatic changes, food and water scarcity, poor air quality, deforestation,
etc [1, 2, 3, 4, 5]. To understand the key trends in urbanization, population
estimation plays a crucial role. Traditionally, population estimation is done
through a census. In this procedure, the population data is systematically
collected and compiled over an administrative unit or a census unit. The
accuracy of the data depends upon the number and size of the administrative
unit, the collection method, the completeness of the survey, and varies
significantly within the rural and urban regions [6].

In recent years, statistical and machine learning methods have been applied
directly to remote sensing data to estimate the population distribution [7,
6]. In general, these methods have been applied either to a smaller region or
based on some country-specific data such as building footprints or other
detailed geospatial data sets which are not easily available in other
countries. Stevens et al. [8] used a random forest approach to estimate the
population at $\sim$ 100 m resolution for Vietnam, Cambodia, and Kenya. They
have incorporated a wide range of remotely-sensed and geospatial datasets,
such as distance to roads, health facility, elevation, land cover, vegetation,
settlements, and nighttime lights, and used the country-specific census data
collected from the National Institute of Statistics for Cambodia, the National
Statistics Office in Vietnam, and the National Bureaus of Statistics for
Kenya. Doupe et al. [9] proposed a new method that used a Convolutional Neural
Network (CNN) to estimate the population by combining Landsat-7 satellite
imagery with (DMSP/OLS) nighttime lights. They trained their model on data
from Tanzania at a 250 m satellite pixel resolution and estimated the
population for Kenya at an 8 km resolution. They have published the code to
reconstruct their data set for Tanzania and Kenya. Another similar CNN
approach has been proposed by Robinson et al. [10] They prepared their data
from US census summary grids in combination with Landsat imagery to estimate
the population in the US counties at a 1 km resolution. Hu et al. [11] also
suggested a deep learning approach by combining satellite imagery from
Landsat-8 and Sentinel-1 and used the Socio-Economic Caste Census survey to
derive the population density for India. In most of the methods above, either
the data is not available for download or could be reconstructed only for a
few cities. Other gridded population products include Global Human Settlement
Population Grid (GHS-POP) [12], WorldPop [13], Oak Ridge National Laboratory’s
LandScan [14], and High-Resolution Settlement Layer (HRSL) [15], and so on.
The difference in the estimation method and ancillary data used in the data
sets lead to different results [6]. Few studies have been conducted to assess
and compare the accuracy of gridded population products by comparing their
estimations with the actual population counts [16, 17]. However, these studies
require collecting and processing the census data. Therefore, it becomes
difficult and time-consuming to reproduce the results or compare the methods.

With our data set, we aim to fill these gaps by providing a systematic
regression and classification scheme for population estimation in 98 European
cities. The cities cover 28 European Union (EU) member states and the four
EFTA countries. It represents a wide range of topography, demography, and
architectural designs across the countries. It would save the cost of
collecting and processing of a new data set to develop and validate the
methods. The data set comprises digital elevation models (DEM), local climate
zone (LCZ), land use (LU), and nighttime lights (VIIRS) in combination with
multi-spectral Sentinel-2 imagery (SEN2), and data from the Open Street Map
initiative (OSM). This multi-data source combination has not been explored
before in the domain of population estimation. We expect that it will be a
valuable addition to the research community for developing sophisticated
approaches in the field of population estimation.

In this paper, we contribute to the current literature by providing a
benchmark data set created from publicly available data sets. We investigated
the fusion of multi-source data over a large number of cities. To demonstrate
the potential capability of our data set, we trained the Random Forest model
using the extracted features from the input data to estimate the population on
our test data set. The initial results indicate that there is a conceivable
potential for the development of powerful machine learning methods with the
So2Sat POP data set.

##  2 Methods

Our region of interest (ROI) is spread over Europe (Figure 1). Initially, we
select all the cities across Europe with 300,000 inhabitants or more in 2014
as per the UN World Urbanization Prospects - The 2014 Revision [18]. Out of
these cities, we select 106 cities based on the availability of population
data. Typically, a city can be described as a permanent large human settlement
defined by the administrative boundaries. However, defining an administrative
boundary could be very tricky because it changes as the census tracts merge or
split over time. The outward expansion of the cities is far beyond their
formal administrative boundaries [19, 20]. So, we employed an algorithm to
determine the extent of the city. The city center coordinates listed in the UN
World Urbanization Prospects - The 2014 Revision [18] have been used as a
starting point in combination with Global Urban Footprint (GUF) [21] which
provides a binary mask of urban vs. non-urban regions. A rectangle centered at
the coordinate of each city is adaptively grown outwards until half of the
area of the rectangle is not built up anymore according to the GUF. To take
the rapid urbanization into account, each side of the rectangles is expanded
by a factor of two (i.e. a factor of four in the area). Since the resulting
rectangles of two neighboring cities might intersect with another, a set of
rules is employed to allocate the intersecting area to one of the two cities
and ensure that each city’s extent covers a unique area. This set of rules is
summarized in Algorithm 1. The algorithm runs recursively and in descending
order, starting from the cities with the highest relative overlap. According
to the set of rules, $city_{a}$ is either merged into $city_{b}$ or the
overlapping are is allocated to $city_{a}$ and removed from $city_{b}$,
depending on the relative size of the overlapping area. The number of cities
is merged to 98 in order to handle these overlaps. As per the defined extent
of each city, we processed and prepared the following input data sources.

Algorithm 1 Allocation of intersecting areas - Pseudocode.

1:$city_{b}\geq city_{a}$

2:for $city_{a}\cap city_{b}$ in data set do

3: while Highest overlap > $0.5*city_{a}$ do

4: if $city_{a}\subset city_{b}$ then

5: Remove $city_{a}$

6: else

7: Merge $city_{a}$ into $city_{b}$

8: end if

9: end while

10: Remove overlapping area from $city_{b}$

11:end for

###  2.1 Input data sources

####  2.1.1 Population data

Generally, the country-specific Census Bureau provides the population data by
the hierarchy of administrative units. The heterogeneity in the size of the
administrative units makes it difficult to use the data directly for the
analysis. In such a scenario, gridded population products with uniformly sized
grid cells offer detailed and consistent population maps. The most common way
to produce the gridded population products is to redistribute the census count
from the administrative level to the fine grid cells, conditioned by ancillary
data such as land use and land cover maps, this method is also known as
dasymetric mapping [22].

The European Statistical System (ESSnet) project, in cooperation with the
European Forum for Geography and Statistics (EFGS), aimed to produce the high
resolution (1 km) population grids from the population census in Europe. The
methodology constitutes aggregation, disaggregation, and hybrid approaches
based on the availability of the data. Aggregation (bottom-up approach) is
assumed to be the best way to produce the population grids [23]. In this
project, approximately 18 countries are using the aggregation or at least a
hybrid method to produce the grid statistics [24]. Due to a lack of geocoded
microdata, the disaggregation (top-down approach) method has been employed in
some regions. In disaggregation, the difference between the modeled population
and the actual population of the region depends upon the size of the
administrative unit. Often higher misplacement of the people is observed in
the larger administrative unit. The overall quality of the product varies
depending upon the data availability ranging from positional accuracy of 0.1m
for each address and building in Austria to up to 100m in Estonia [25] (for
details see the document named ‘GEOSTAT_grid_POP_1K_2011_V2_0_QA.pdf’ in the
folder titled ‘GEOSTAT-grid-POP-1K-2011-V2-0-1’). The European Environment
Agency published the population grids for EU28+ EFTA countries that comprise
approximately 4.3 million km2 with 480 million inhabitants [23]. It is freely
available via Eurostat, for non-commercial purposes. Further details regarding
the product standards, methodology, and quality assessments can be found at
the GEOSTAT 1B project website [26].

####  2.1.2 Sentinel-2

The Sentinel-2 mission [27] was launched in June 2015 by the European Space
Agency (ESA) and consists of two identical satellites, 2A and 2B, phased at
180 degrees toward each other. Sentinel-2 satellites provide multi-spectral
optical images with a span of 13 spectral bands at a spatial resolution of 10
m, 20 m, 60 m. Thus, freely available Sentinel-2 imagery offers great
potential for the fine-scale mapping of human settlements. The analysis-ready
cloud-free mosaics of Sentinel-2 data have been achieved using the cloud
detection approach based on a pixel-wise analysis [28]. We have used all four
seasonal sets of Sentinel-2 images to capture the seasonal variation in the
data. The spring, summer, and autumn seasons are from 2017 and winter is from
2016.

####  2.1.3 TanDEM-X Digital Elevation Model

The goal of the TanDEM-X mission is to create a high-quality 3D topographic
map of the Earth that is homogeneous in quality and unparalleled in accuracy.
The data acquisition period for the global DEM product is between December
2010 and January 2015 and the global DEM has been produced in September 2016.
With the coverage of 150 Million km2 of the complete Earth’s landmasses and 10
m absolute height accuracy (90% linear error) [29], it is suitable for various
applications in Environmental research such as land cover and land use
analysis, urban planning, climate change, etc. Currently, the TanDEM-X digital
elevation model is the most accurate digital elevation model (DEM) available
globally [30]. We have used the freely available TanDEM-X 90 m (3 arcsec) DEM
global product [31] that contains the final, global Digital Elevation Model
(DEM) of the landmasses of the Earth [32].

####  2.1.4 Local climate zones

Local climate zones (LCZ) are a freely available systematic housing density
classification data, formally developed to standardize urban heat island
studies [33]. Based on the land surface and properties, it consists of 17
structural types of which 10 are describing built zones from compact high-rise
to open low-rise and 7 are natural zones ranging from dense vegetation to bare
lands. Thus, each zone is characterized by the built-up and land cover
properties. We have used the urban local climate zone classifications,
So2SatLCZ v1.0 at the resolution of 100 m, produced by fusing the freely
available satellite data from Sentinel-1 and Sentinel-2 satellites using deep
learning [34]. The patches in this data set are hand-labeled as per the local
climate zones classification scheme by 15 domain experts, followed by a visual
and quantitative evaluation process over six months. This benchmark data set,
therefore, is of potential use for urbanologists, demographers,
climatologists, and many other researchers.

####  2.1.5 Nighttime lights

Observational Nighttime lights (NTL) show a strong correlation with the
spatial distribution of the human population [35]. We studied the two widely
used Nighttime lights data, the DMPS-OLS and the NPP-VIIRS satellite images.
NPP-VIIRS with a better spatial resolution (15 arc-second, about 500 m) has a
higher potential in modeling the socio-economic indicators [36]. The annual
composites of VIIRS nighttime lights derived from monthly mean data are freely
available from 2012 to 2019 (https://eogdata.mines.edu/products/vnl/). We used
a masked average radiance version of VNL V2, which is a preprocessed version
free of outliers from fleeting events [37].

####  2.1.6 OSM

OpenStreetMap (OSM) (http://www.openstreetmap.org) is an open and crowd-
sourced platform for maps, available under the Open Data Commons Open Database
License (ODbL) (http://www.opendatacommons.org/licenses/odbl/1.0/). Geo-
referenced locations on a very detailed level can be entered as nodes, ways,
or relations and specified with informative tags. These locations include any
type of buildings, streets, boundaries, water bodies, etc.[38]. Our data set
contains low-level and high-level features derived from OSM data. At a low
level, a locale with a high number of certain nodes (e.g. supermarkets, gas
stations, residential buildings, schools, etc.) correlates strongly with the
population in the vicinity. A simple counter statistic of such node-types is a
strong indicator and can serve as a feature vector for population density
estimation. As high-level features we extracted the building functions from
OSM building tags to represent urban land use. This information is closely
linked to employment, social support, and population [39] and indicates the
interaction between human activities and the environment.

###  2.2 Data Preprocessing

We employed two-step preprocessing for all input data sources. In the first
step, input data for each city has been created. The second step includes the
creation of the 1x1 km patches for each city. Figure 2 outlines the step-by-
step preprocessing of all the data that has been employed to create the input
data for each city. All the input data has been cropped using the city
boundaries defined by our algorithm. DEM data is standardized by subtracting
the DEM mean and then scaled to unit variance. To match the spatial resolution
of the Sentinel-2 RGB bands, other data sources such as the digital elevation
model (DEM), local climate zone (LCZ), land use (LU), and VIIRS nighttime
lights (VIIRS) have been up-sampled to 10 m. Since the input data sets are
from various sources, they are in different Coordinate Reference Systems
(CRS). The VIIRS is in WGS84 (EPSG:4326), the LCZ, DEM, and SEN2 data are in
the Universal Transverse Mercator (UTM) zone and the population grid is in the
EPSG:3035 - ETRS89-extended / LAEA Europe. To align the input data with the
population grid, all the input data have been reprojected from their
corresponding coordinate system to the EPSG:3035 coordinate reference system.

####  2.2.1 Low-Level Features

The OSM planet dump of 2017-07-03 is being downloaded directly from the
OpenStreetMap archive (https://planet.osm.org/planet/2017/). To reduce the
computation time by touching the huge planet file only once, we extract the
bounding box of each corresponding city and subsequently 1x1 km patches from
these extracted city-dumps with the Osmosis command line
tool(https://github.com/openstreetmap/osmosis). The OSM node statistics are
extracted for each 1x1 km patch of the city using the OSMnx python Library
[40]. Table 1 shows the considered OSM tags over which the statistical
analysis has been done.

aerialway | boundary | geological | landuse | natural | public_transport | sport | waterway
---|---|---|---|---|---|---|---
aeroway | building | healthcare | leisure | office | railway | telecom | addr:housenumber
amenity | craft | highway | man_made | place | route | tourism | restrictions
barrier | emergency | historic | military | power | shop | water | other: True
Table 1: Nodes with these OSM tags are considered for the statistical
analysis/counting of the corresponding 1x1 km patch.

####  2.2.2 High-Level Features

To create the land use data we analyzed three different tags of buildings in
OSM: building, amenity, and shop. For each of them, OSM provides a guideline
on possible values. In total, there are 341 possible values for these three
tags, which are mapped to a homogenized and simplified land use classification
scheme: commercial, industrial, residential, and other. As the three tags can
occur jointly, we make sure that they are not contradicting each other and
omit buildings that have inconsistent values. Moreover, the tags are captured
as free-form text fields and hence, OSM contributors are not restricted to use
them but can enter any text. After homogenizing the semantic information, we
convert the vector data of building polygons into raster data. The rasterized
value represents the area covered by building polygons inside a raster pixel.
For scaling, we divided the area by the area of a pixel to yield a relative
number of how much of a pixel is covered by building polygons. Applying this
procedure for each land use class results in a four-band raster with
corresponding land use proportions. The output of the first step of data
preprocessing for Munich city can be seen in Figure 3.

In the second step, for each city, we created the patches using all the input
data processed in the first step. The population grid of a city is used as a
reference grid to crop all the other input data. The size of a grid cell in
the population grid is 1x1 km and each cell represents the population count
living per square km of the cell. The grid cells along the border of the
population grid might belong to two adjacent cities. To avoid this duplicity
such that one grid cell should belong to only one city, we applied an area
thresholding to all the cells. Only the cells with an area greater than
900,000 m2 or size $\sim$ 0.95 x 0.95 km have been considered as a part of the
city. This eliminates the grid cells from the edges which are not fully
included in the city boundary. Also, the reference Geostat population grid has
some missing cells. The missing cells mostly contain the uninhabited regions.
Such regions cover the green fields and water bodies. We included these
missing cells in our data set to allow for zero population predictions. The
patches cropped using the population grid have been assigned the population
count as of the corresponding grid cell. The rest of the missing/uninhabited
patches have been assigned as zero population count. For some applications
such as environmental impact assessments, land use analysis, climate change,
etc. it is sufficient to know the range of the people living in a region. So,
we further preprocessed the population grids by binning the population count
to a population class. We assigned a grid cell, Class 0, if the population
count of the cell is zero, Ccell=0 if Pcell=0 and subsequently Ccell=1 if 20
$\leq$ Pcell ${<}$ 21, Ccell=2 if 21 $\leq$ Pcell ${<}$ 22, Ccell=3 if 22
$\leq$ Pcell ${<}$ 23…..Ccell=${k+1}$ if 2k $\leq$ Pcell ${<}$ 2k+1 where k
$\in$ $\mathbb{N}$. For our data, the highest value of ${k}$ is 16. This
process of discretization is inspired by Robinson et al. [10]. Thus, each grid
cell has been assigned a population class in addition to the absolute
population count. It would give more flexibility to the end-users to develop
either regression or a classification model for the task considering the
requirements of the application.

Figure 4 illustrates the patch creation process. For every grid cell, a total
of 9 patches, one from each data source have been created. We called the 9
patches corresponding to one grid cell, a patch-set. Each patch-set represents
a population count as of the corresponding population grid cell and a
population class depending on which bin the population count of the grid cell
falls at the resolution of 1 km. Figure 5 depicts the odd-numbered class
samples from our data set with their corresponding patch-set, population class
and population count. The lower classes represent sparsely populated regions.
Patches belonging to lower classes mostly contain green fields, water bodies,
and bare land. As the class number goes higher, patches contain few low-rise
to dense high-rise built-up regions. In other words, patches from lower to
higher classes represent rural to urban regions.

##  3 Data Records

The final data set consists of two parts, So2Sat POP Part1 [41] and So2Sat POP
Part2 [42]. All the data patches except OSM data are available as GeoTiff
images. Along with the raw osm patches, we also provide the features extracted
from the OSM data as separate Comma Separated Value (CSV) files. So2Sat POP
Part1 consists of the patches from local climate zones, land use, nighttime
lights, Open Street Map features, and from all seasons(autumn, summer, spring,
winter) of Sentinel-2 imagery, a total of 1,104,688 patches. So2Sat POP Part2
consists of patches from the digital elevation model and Open Street Map only,
a total of 276,172 patches. So2Sat POP Part1 has the storage requirement of
$\sim$ 98 GB and So2Sat POP Part2 requires $\sim$ 5.20 GB.

###  3.1 Data set structure and naming convention

Both parts of the data set consist of a predefined train and test split. Out
of 98 cities, 80 cities ($\sim$ 80% of the data) have been randomly selected
as a training set and the rest of the 18 cities ($\sim$ 20% of the data)
constitucorrecte the test set. The city folder has been named as
xxxx_xxxxx_city_name, where xxxx_xxxxx constitutes randomly generated
identification number and the postal code of the city.

The digital elevation model (DEM), local climate zone (LCZ), land use (LU),
nighttime lights (VIIRS), Sentinel-2 autumn data (sen2autumn), Sentinel-2
spring data (sen2spring), Sentinel-2 summer data (sen2summer), Sentinel-2
winter data (sen2winter), Open Street Map (OSM) and its corresponding
extracted features (osm_features) are the utilized input sources. All city
folders in So2Sat POP Part1 contain 7 sub-folders, one for each input data
source except the OSM and DEM data, a separate folder for processed osm
features, and a comma-separated value (*.csv) file that contains the absolute
population count and population class for each patch. In So2Sat POP Part2, the
city folders contain the Open Street Map and digital elevation model data sub-
folder. All data folders have their class sub-folders. The class folders have
been named as Class_x where x denotes the class number. The number of class
folders in a city depends on its population distribution. For example, Malaga
has the highest class folder as 16 because the highest population count in the
1x1 km area of the city is 39535 while in Riga the highest population count is
15839 so the highest class folder in its city folder is 14. A patch has been
assigned a unique identification code using the naming convention of its
corresponding population grid cell. The naming of the grid cell (based on the
LAEA grid) begins with the size of the cell (1 km) followed by the coordinates
(in km) of the lower left-hand corner, starting with the letter N followed by
the latitude and E followed by the longitude, e.g. 1kmN4101E4453 [24]. The
patches that do not correspond to the grid cells in the population grid have
been simply given a numeric identification number.

###  3.2 Access and availability

The data set is available for download at the persistent links provided by the
official media library of the Technical University of Munich (TUM). So2Sat POP
Part 1 (https://mediatum.ub.tum.de/1633792) is distributed under the Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/) and So2Sat POP Part 2
(https://mediatum.ub.tum.de/1633795) is distributed under the Creative Commons
Attribution Share-Alike International License
(http://creativecommons.org/licenses/by-sa/4.0/). Kindly cite this paper when
the data set is used.

##  4 Technical Validation

To demonstrate the suitability of the data set for population estimation, we
implemented the popular Random Forest (RF) algorithm [43], because of its
flexibility, efficiency in handling the noisy input data, and relative
resistance to overfitting [44]. Another advantage of the random forest
algorithm is that it is easy to measure the relative importance of each
feature on the prediction. We implemented the supervised random forest
algorithm in Python with the scikit-learn library for both the regression and
classification tasks. We used grid search to automatically fine-tune the
number of trees to grow and the maximum number of features to consider
splitting a node and assess the performance through a 10-fold cross-
validation. To train the model, different features are constructed from all
the input data patches of 80 train-set cities. The constructed features
include min, max, mean, median, and standard deviation from only the RGB bands
of Sentinel-2 imagery, mean and max for digital elevation model and
nightlights, the total area covered by each class in land use, majority class
for local climate zone, and osm-based features such as street density,
presence of highways, railways, etc., extracted from osm patches. Using this
process, we calculated 125 features for each patch. For regression, we used
the absolute population count as the response variable while in
classification, class labels are used as ground truth.

The trained model has been evaluated on the 18 unseen test cities. Figure 6
shows only the twelve most relevant features that have been selected by the
Random Forest algorithm and used to estimate the population count and a
population class for the set of test data. The features extracted from the OSM
data, nightlights, and LCZ classes have been ranked as the most important
features for both tasks. For regression, we calculated the root-mean-square
error (RMSE) and the mean absolute error (MAE). Table 2 indicates the
performance of the regression model. Since the data set is imbalanced due to
the higher percentage of non-urban regions than the urban regions, we used
balanced accuracy to evaluate the classification performance. Also, we used
macro-averaged Precision, Recall, and F1-score metrics to treat each class
equally, regardless of any imbalance. Table 3 summarized the results from the
classification. To further describe the performance of the classification
model, we plotted a normalized confusion matrix on the set of test data.
Figure 7 (b) illustrates that the model is confident in predicting the higher
classes (urban regions) while it does not perform well on the lower classes,
which represent the sparsely populated regions. The initial three classes
(Class 1, 2, and 3) represent the regions where the population count is in the
range 1-8, and, likely, the features among these three classes are not
distinguishable enough. So, it becomes very difficult to differentiate the
patches from Class 1 to Class 3. One simple solution to it is the merging of
these initial classes, which represent the population count range of less than
8. It could significantly improve the results without the loss of any critical
information. To give more flexibility to the user, we provided the data set
without any such post-processing, and it could easily be rearranged as per the
specific needs. For regression, to visually assess the model fit, we plotted
the actual population count of the grid cells versus the population count
predicted for each grid cell of the test data. Figure 7 (a) indicates that the
model underestimates the actual values for the patches where the population
count is high while for the patches with a population count of less than
15.000, it is relatively a good fit. We believe that with more sophisticated
features and machine learning methods, a powerful model could be developed to
estimate the population using our data set.

##  5 Usage Notes

In this paper, we proposed a unique data set that combines multi-data sources
which have not been explored at the cross-country level in this domain before.
Our data set covers many cities across Europe thus offering diverse topography
and demography. Also, the fusion of distinct data sources gives much more
information about the landscape and other socio-demographic attributes of a
region.

The availability of precise and detailed population data varies from country
to country. The population grid used as a reference grid in our data set is
available throughout Europe at a consistent resolution of 1 km, however, its
precision may vary depending upon the quality of the data available in a
country. Also, the reference population grid is based on the 2011 population
and housing census data while other input data sources belong to a different
time frame, for example, Sentinel-2 data belongs to the year 2016 as the
mission itself has been launched in 2015. The time difference between the
collection of the population data and other corresponding input data might
introduce some noise to the evaluations. Usually, the population census has
been conducted once in a decade, so it becomes very difficult to collect the
other corresponding data from the same year, especially when collecting from
multiple data sources. We still believe that this data set will be helpful in
the development of the machine-learning based approaches in this domain. Also,
the evaluation of different methodologies has always been a challenge due to
the unavailability of the freely available benchmark data set. So, we hope
this data set will be useful for the comparison and evaluations of the
different approaches in the state of the art.

The data preparation is performed in Python using the common libraries for the
geospatial data such as Geopandas, Fiona, Rasterio, and Geospatial Data
Abstraction Library (GDAL). Our computing system consists of an HPC Server
with 2x AMD EPYC ROME 7402 CPUs (48 core) and needed 5-6 days to complete the
OSM patch files and their statistics. Along with the patches, we also provided
the Comma Separated Value (CSV) files for each city that contains the actual
population count and its corresponding population class for each grid cell.
Thus, the end user could rearrange the data set using its corresponding
population count according to the requirements of the use case.

We hope that this basic data set will make it possible to develop new
statistical and machine learning approaches that derive the population in
higher spatial resolution, but above in more consistency. It is precisely this
consistency that is often lacking in today’s population data sets across
countries. This data set is intended to lay the foundation for improved
comparative studies in various application domains. This improved database for
population distribution may be central information both, in the academic and
non-academic fields. Be it as necessary information for spatial or urban
planning, such as the provision of living space or socio-technical
infrastructure, be it in the locale of risk analysis or coordination in an
event, be it in the validation, support, or updating of censuses, or be it in
comparative studies on topics such as migration in general and, more
specifically, on urbanization or sub-urbanization trends.

##  6 Code availability

Python is used for all the analyses and implementations. The code to create
the features for each city and to run the baseline experiments is available on
our GitHub project (https://github.com/zhu-xlab/So2Sat-POP).

## Acknowledgements

This research was funded by the European Research Council (ERC) under the
European Unions Horizon 2020 research and innovation program with the grant
number ERC-2016-StG-714087 (Acronym: So2Sat, project website: www.so2sat.eu),
Helmholtz Association under the framework of the Helmholtz AI (grant number:
ZT-I-PF-5-01) –Local Unit “Munich Unit @Aeronautics, Space and Transport
(MASTr),” and Helmholtz Excellent Professorship “Data Science in Earth
Observation – Big Data Fusion for Urban Research (grant number: W2-W3-100) and
by the German Federal Ministry of Education and Research (BMBF) in the
framework of the international future AI lab "AI4EO – Artificial Intelligence
for Earth Observation: Reasoning, Uncertainties, Ethics and Beyond" (Grant
number: 01DD20001). Additionally, Sugandha Doda is supported by the Helmholtz
Association under the joint research school “Munich School for Data Science -
MUDS”.

## Author contributions statement

S.D. drafted the manuscripts, S.D., E.H., M.K. and K.O. undertook data
preparation, S.D. implemented and performed experiments, S.D. performed
technical validation of the results. S.D., X.Z., Y.W., H.T., and M.K. edited
the manuscript. All authors read and approved the final version of the
manuscript.

## Competing interests

The authors declare no competing of interest.

## Figures & Tables

Figure 1: The orange dots on the figure above indicate the location of
selected EU cities in our study. Figure 2: Step-by-step preprocessing of all
the input data sources to prepare the corresponding input data for each city.
Figure 3: All the input data for the Munich city, created using the first step
of data preprocessing. Figure 4: Patch creation process, second step of data
preprocessing. All input data sources have been cropped for each cell in the
population grid. The size of each patch is 1 x 1 km. Figure 5: Sample patches
from the odd numbered classes of our data set. Lower classes depicts sparsely
populated regions while higher classes depicts densely populated regions.

(a)

(b)

Figure 6: Random Forest feature importances based on mean decrease in impurity
(MDI). The higher the value the more important the feature. Plot shows only
the twelve most relevant features for both regression (a) and classification
(b)

(a)

(b)

Figure 7: (a) Predicted vs. Actual Values for regression, the model fits good
except for the high population counts where the points appeared dispersed from
regressed diagonal line (b) Confusion matrix for classification, normalized by
class support size (number of patches in each class). Confusion among the non-
urban classes is higher than among the urban classes. Regression
---
Method |  Internal OOB Score |  RMSE |  MAE
RF |  0.864 |  1276.26 |  463.35
Table 2: Evaluation of Random Forest model to estimate the population count on
test data set. Classification
---
Method |  Accuracy |  Balanced Accuracy |  F1 score |  Precision |  Recall
RF |  0.5913 |  0.3795 |  0.3833 |  0.4533 |  0.3795
Table 3: Evaluation of Random Forest model to predict the population class on
test data set.

## References

  * [1] Robert I McDonald, Pamela Green, Deborah Balk, Balazs M Fekete, Carmen Revenga, Megan Todd, and Mark Montgomery.  Urban growth, climate change, and freshwater availability.  Proceedings of the National Academy of Sciences, 108(15):6312–6317, 2011.
  * [2] Andrew J Tatem.  Mapping the denominator: spatial demography in the measurement of progress.  International health, 6(3):153–155, 2014.
  * [3] Gordon McGranahan, Deborah Balk, and Bridget Anderson.  The rising tide: assessing the risks of climate change and human settlements in low elevation coastal zones.  Environment and urbanization, 19(1):17–37, 2007.
  * [4] Xiaoyan Zhang, Lijian Han, Haiyan Wei, Xiaorui Tan, Weiqi Zhou, Weifen Li, and Yuguo Qian.  Linking urbanization and air quality together: A review and a perspective on the future sustainable urban development.  Journal of Cleaner Production, page 130988, 2022.
  * [5] Sylvia Szabo.  Urbanisation and food insecurity risks: Assessing the role of human development.  Oxford Development Studies, 44(1):28–48, 2016.
  * [6] Stefan Leyk, Andrea E Gaughan, Susana B Adamo, Alex de Sherbinin, Deborah Balk, Sergio Freire, Amy Rose, Forrest R Stevens, Brian Blankespoor, Charlie Frye, et al.  The spatial allocation of population: a review of large-scale gridded population data products and their fitness for use.  Earth System Science Data, 11(3):1385–1409, 2019.
  * [7] Shuo-sheng Wu, Xiaomin Qiu, and Le Wang.  Population estimation methods in gis and remote sensing: A review.  GIScience & Remote Sensing, 42(1):80–96, 2005.
  * [8] Forrest R Stevens, Andrea E Gaughan, Catherine Linard, and Andrew J Tatem.  Disaggregating census data for population mapping using random forests with remotely-sensed and ancillary data.  PloS one, 10(2):e0107042, 2015.
  * [9] Patrick Doupe, Emilie Bruzelius, James Faghmous, and Samuel G Ruchman.  Equitable development through deep learning: The case of sub-national population density estimation.  In Proceedings of the 7th Annual Symposium on Computing for Development, pages 1–10, 2016.
  * [10] Caleb Robinson, Fred Hohman, and Bistra Dilkina.  A deep learning approach for population estimation from satellite imagery.  In Proceedings of the 1st ACM SIGSPATIAL Workshop on Geospatial Humanities, pages 47–54, 2017.
  * [11] Wenjie Hu, Jay Harshadbhai Patel, Zoe-Alanah Robert, Paul Novosad, Samuel Asher, Zhongyi Tang, Marshall Burke, David Lobell, and Stefano Ermon.  Mapping missing population in rural india: A deep learning approach with satellite imagery.  In Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society, pages 353–359, 2019.
  * [12] Sergio Freire, Kytt MacManus, Martino Pesaresi, Erin Doxsey-Whitfield, and Jane Mills.  Development of new open and free multi-temporal global population grids at 250 m resolution, 2016.
  * [13] WorldPop.  School of geography and environmental science, university of southampton; department of geography and geosciences, university of louisville; departement de geographie, universite de namur) and center for international earth science information network (ciesin), columbia university.  Global High Resolution Population Denominators Project-Funded by The Bill and Melinda Gates Foundation (OPP1134076), 2018.
  * [14] Budhendra Bhaduri, Edward Bright, Phillip Coleman, and Jerome Dobson.  Landscan.  Geoinformatics, 5(2):34–37, 2002.
  * [15] High Resolution Settlement Layer.  Facebook connectivity lab and center for international earth science information network-ciesin-columbia university. source imagery for hrsl© 2016 digitalglobe, 2016.  Accessed on: 2022-06-21.
  * [16] Ruxia Chen, Huimin Yan, Fang Liu, Wenpeng Du, and Yanzhao Yang.  Multiple global population datasets: Differences and spatial distribution characteristics.  ISPRS International Journal of Geo-Information, 9(11):637, 2020\.
  * [17] Richard Sliuzas, Monika Kuffer, and Thomas Kemper.  Assessing the quality of global human settlement layer products for kampala, uganda.  In 2017 Joint Urban Remote Sensing Event (JURSE), pages 1–4. IEEE, 2017.
  * [18] United Nations United Nations.  World Urbanization Prospects: 2014 Revision.  United Nation, 2014.
  * [19] UN Habitat.  State of the world’s cities 2012/2013: Prosperity of cities.  Routledge, 2013.
  * [20] Hannes Taubenböck, Matthias Weigand, Thomas Esch, Jeroen Staab, Michael Wurm, Johannes Mast, and Stefan Dech.  A new ranking of the world’s largest cities—do administrative units obscure morphological realities?  Remote Sensing of Environment, 232:111353, 2019.
  * [21] Thomas Esch, Wieke Heldens, Andreas Hirner, Manfred Keil, Mattia Marconcini, Achim Roth, Julian Zeidler, Stefan Dech, and Emanuele Strano.  Breaking new ground in mapping human settlements from space–the global urban footprint.  ISPRS Journal of Photogrammetry and Remote Sensing, 134:30–42, 2017\.
  * [22] Mitchel Langford.  Rapid facilitation of dasymetric-based population interpolation by means of raster pixel maps.  Computers, Environment and Urban Systems, 31(1):19–32, 2007.
  * [23] Francisco Javier Gallego.  A population density grid of the european union.  Population and Environment, 31(6):460–473, 2010.
  * [24] Efgs - essnet project geostat 1b - final report.  https://www.efgs.info/wp-content/uploads/geostat/1b/GEOSTAT1B-final-technical-report.pdf.  Accessed on: 2022-10-05.
  * [25] Eurostat gisco geostat $1\,km^{2}$ population grid.  https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/population-distribution-demography/geostat, 2011\.  Accessed on: 2022-10-05.
  * [26] Efgs - essnet project geostat 1b - geostat 2011 quality assessment.  http://www.efgs.info/wp-content/uploads/geostat/1b/GEOSTAT1B-Appendix17-GEOSTAT-grid-POP-1K-ALL-2011-QA.pdf.  Accessed on: 2022-10-05.
  * [27] Matthias Drusch, Umberto Del Bello, Sébastien Carlier, Olivier Colin, Veronica Fernandez, Ferran Gascon, Bianca Hoersch, Claudia Isola, Paolo Laberinti, Philippe Martimort, et al.  Sentinel-2: Esa’s optical high-resolution mission for gmes operational services.  Remote sensing of Environment, 120:25–36, 2012.
  * [28] M Schmitt, Lloyd H Hughes, C Qiu, and X X Zhu.  Aggregating cloud-free Sentinel-2 images with Google Earth Engine.  In ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, volume IV-2/W7, pages 145–152, 2019.
  * [29] Birgit Wessel, Martin Huber, Christian Wohlfart, Ursula Marschalk, Detlev Kosmann, and Achim Roth.  Accuracy assessment of the global tandem-x digital elevation model with gps data.  ISPRS Journal of Photogrammetry and Remote Sensing, 139:171–182, 2018.
  * [30] Thomas Esch, Julian Zeidler, Daniela Palacios-Lopez, Mattia Marconcini, Achim Roth, Milena Mönks, Benjamin Leutner, Elisabeth Brzoska, Annekatrin Metz-Marconcini, Felix Bachofer, et al.  Towards a large-scale 3d modeling of the built environment—joint analysis of tandem-x, sentinel-2 and open street map data.  Remote Sensing, 12(15):2391, 2020.
  * [31] German aerospace center (dlr) (2018): Tandem-x - digital elevation model (dem) \- global, 90m.
  * [32] Birgit Wessel.  Tandem-x ground segment–dem products specification document, 2018.  EOC, DLR, Oberpfaffenhofen, Germany, Public Document TD-GS-PS-0021, Issue 3.2. [Online]. Available online: https://tandemx-science.dlr.de/.
  * [33] Ian D Stewart and Tim R Oke.  Local climate zones: Origins, development, and application to urban heat island studies.  In Proceedings of the Annual Meeting of the American Association of Geographers, Seattle, WA, USA, pages 12–16, 2011.  Accessed on: 2022-06-21.
  * [34] Xiao Xiang Zhu, Jingliang Hu, Chunping Qiu, Yilei Shi, Jian Kang, Lichao Mou, Hossein Bagheri, Matthias Haberle, Yuansheng Hua, Rong Huang, Lloyd Hughes, Hao Li, Yao Sun, Guichen Zhang, Shiyao Han, Michael Schmitt, and Yuanyuan Wang.  So2sat lcz42: A benchmark data set for the classification of global local climate zones [software and data sets].  IEEE Geoscience and Remote Sensing Magazine, 8(3):76–89, 2020.
  * [35] Qing Liu, Paul C Sutton, and Christopher D Elvidge.  Relationships between nighttime imagery and population density for hong kong.  Proc. Asia-Pac. Adv. Netw, 31:79, 2011.
  * [36] Kaifang Shi, Bailang Yu, Yixiu Huang, Yingjie Hu, Bing Yin, Zuoqi Chen, Liujia Chen, and Jianping Wu.  Evaluating the ability of npp-viirs nighttime light data to estimate the gross domestic product and the electric power consumption of china at multiple scales: A comparison with dmsp-ols data.  Remote Sensing, 6(2):1705–1724, 2014.
  * [37] Christopher D Elvidge, Mikhail Zhizhin, Tilottama Ghosh, Feng-Chi Hsu, and Jay Taneja.  Annual time series of global viirs nighttime lights derived from monthly averages: 2012 to 2019.  Remote Sensing, 13(5):922, 2021.
  * [38] Map features documentation wiki.  https://wiki.openstreetmap.org/wiki/Map_features.  Accessed on: 2021-08-21.
  * [39] Xiangmei Li, Ying Wang, Jiangfeng Li, and Bin Lei.  Physical and socioeconomic driving forces of land-use and land-cover changes: A case study of wuhan city, china.  Discrete Dynamics in Nature and Society, 2016, 2016.
  * [40] Geoff Boeing.  Osmnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks.  Computers, Environment and Urban Systems, 65:126–139, 2017.
  * [41] Sugandha Doda, Yuanyuan Wang, Matthias Kahl, Eike Jens Hoffmann, Kim Ouan, Hannes Taubenböck, and Xiao Xiang Zhu.  So2sat pop part 1, mediatum, 2022.
  * [42] Sugandha Doda, Yuanyuan Wang, Matthias Kahl, Eike Jens Hoffmann, Kim Ouan, Hannes Taubenböck, and Xiao Xiang Zhu.  So2sat pop part 2, mediatum, 2022.
  * [43] Leo Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [44] Taïs Grippa, Catherine Linard, Moritz Lennert, Stefanos Georganos, Nicholus Mboga, Sabine Vanhuysse, Assane Gadiaga, and Eléonore Wolff.  Improving urban population distribution models with very-high resolution satellite information.  Data, 4(1):13, 2019.


# Non-axisymmetric modes in ultrathin annular liquid films coating a
cylindrical fibre

D. Moreno-Boza<EMAIL_ADDRESS>Departamento de Ingeniería Térmica y de
Fluidos, Universidad Carlos III de Madrid. Avda. de la Universidad 30, 28911,
Leganés, Madrid, Spain.  A. Sevilla  Departamento de Ingeniería Térmica y de
Fluidos, Universidad Carlos III de Madrid. Avda. de la Universidad 30, 28911,
Leganés, Madrid, Spain.

(August 30, 2024)

###### Abstract

This paper presents a detailed analysis of the three-dimensional stability
properties of an annular liquid film coating a cylindrical fibre in the
presence of van der Waals (vdW) interactions, whose influence depends on the
wettability of the solid by the liquid. Under wetting conditions, vdW
interactions can stabilise a uniform annular film when its thickness is
smaller than a critical value that depends only on the fibre radius, the
Hamaker constant, and the surface tension coefficient [26]. In contrast, under
non-wetting conditions, both surface tension and vdW forces contribute to
destabilise the interface, and non-axisymmetric modes may become dominant
depending on the thickness of the film and the relative strength of the
surface tension and vdW forces. We perform temporal stability analyses of both
the Stokes and lubrication equations of motion, allowing us to reveal the
dominant azimuthal mode, as well as the optimal axial wavenumber and the
corresponding temporal growth rate, as a function of the relevant governing
parameters.

††preprint: Accepted in Phys. Fluids

##  I Introduction

A static annular liquid layer of uniform radius coating the external surface
of a cylindrical solid fibre is generically unstable due to the Plateau-
Rayleigh mechanism, evolving into a sequence of axisymmetric collars along the
fibre which are connected by thin annular lobes. The temporal linear stability
of the uniform-radius state was solved by Goren [9] using the full Navier-
Stokes equations, leading to an exact dispersion relation for the temporal
modes. In subsequent papers [10, 27] nonlinear equilibrium solutions under
quasi-static conditions were found, finding good agreement with the
experimentally recorded interface shapes at different times during the
unstable evolution.

Since the linear cut-off wavelength is equal to the circumference of the
coating, the typical axial length scales are much larger than the radial ones,
thus rendering lubrication theory as the natural framework in developing
models able to tackle the linear and non-linear dynamics in a simple and
accurate way. Indeed, Dumbleton and Hermans [6] showed that the linear
stability properties of the film are correctly captured by a leading-order
lubrication model, while its nonlinear long-time evolution was studied in a
series of more recent papers [12, 15, 34, 20], which provide a complete
picture of the flow for arbitrary times in the absence of van der Waals (vdW)
forces. In addition, many studies have been devoted to the annular coating
process leading to the type of liquid films considered in the present work
[see 25, and references therein], and have also considered other effects like
the gravitational fall of the film along a vertical fibre, which saturates the
Plateau-Rayleigh instability [8], and leads to different non-linear regimes
depending on the convective or absolute nature of the instability [16, 7, 28,
17]. Thermocapillary effects [21], the presence of surfactants [2] and the
influence of external electric fields [33] have also been addressed. A closely
related configuration is that of a thin annular film coating the inner side of
a tube, which has also been widely studied due to its relevance in physiology
[see, for instance, 12, 11, 5]. Many more examples can be found in the
excellent and thorough reviews due to Oron et al. [24] and Craster and Matar
[3], where annular films are discussed in the more general context of thin
liquid films.

When the thickness of the film is smaller than about $100\leavevmode\nobreak\
\textrm{nm}$, long-range van der Waals (vdW) forces become comparable to
surface tension forces, and their influence depends crucially on the
wettability of the solid by the liquid. In cases where the liquid wets the
solid, it was shown by Quéré et al. [26] that vdW interactions are able to
stabilise the annular film to a uniform thickness that depends only on the
fibre radius, the Hamaker constant associated with the disjoining pressure and
the surface tension coefficient of the interface separating the liquid layer
from the ambient air. In contrast, if the liquid does not wet the solid the
coating is always unstable, since both surface tension and vdW forces
contribute to destabilise the interface. The resulting unstable flow then
combines the effects of both forces, and it is expected that the film dynamics
be modified by vdW interactions when the film thickness is sufficiently small.
In particular, the linear stability properties of these annular ultrathin non-
wetting films, which to the best of our knowledge have not been reported
before in full detail, are of interest to describe the resulting dewetting
patterns. In contrast with the case of wetting films, where only axisymmetric
modes are unstable, it was concluded by [19] using an approximate model that
non-axisymmetric modes may also become unstable under non-wetting conditions,
and that they may even become dominant depending on the thickness of the film
and on the relative strength of surface tension and vdW forces.

Thus, our main objective is to provide a detailed account of the three-
dimensional stability properties of annular liquid films in the presence of
vdW interactions, including both the wetting and the non-wetting scenarios.
The mathematical models are formulated in §II, and the corresponding temporal
stability analyses are presented in §III distinguishing between wetting and
non-wetting cases. Finally, some conclusions are drawn in §IV.

##  II Formulation of the problem

An annular film of a Newtonian liquid of density $\rho$ and viscosity $\mu$
coating a cylindrical solid substrate of radius $R$ is considered. The local
film thickness will be referred to as $h$ such that the interface has the
radial position $r=R+h(x,\theta,t)$ in the cylindrical coordinate system
$(x,r,\theta)$. The film, of initial thickness $h_{o}$, is surrounded by
passive atmospheric air at constant pressure $p_{a}$. The coefficient of
surface tension $\sigma$ between the air and the liquid will be taken as a
constant, since the presence of temperature gradients and adsorbed surfactants
will be disregarded in the present work. The long-range dispersive van der
Waals forces will be modelled by means of a disjoining pressure
$\phi^{*}=A/(6\pi h^{3})$, where $A$ is the effective Hamaker constant, which
will be considered to be either positive or negative to contemplate,
respectively, both the non-wetting and the wetting cases in the analysis. The
reader is referred to the review of Craster and Matar [3] (Section V.B) or the
books by Israelachvili [14] and Blossey [1] for more details on the derivation
of the disjoining pressure. Also, short-range Born repulsion effects will be
neglected throughout the analysis. A detailed schematic of the configuration
is shown in Figure 1.

Figure 1: Schematic of the flow configuration illustrating the main governing
parameters as well as the cylindrical coordinate system employed in the
description. The perturbed air-liquid interface (dashed lines) is shown for an
axisymmetric mode of instability in the main panel, whereas different helical
modes $m=0,1,2$ are shown in the right-hand panel. Arbitrary values of $k$ and
$\eta$ are used for visualisation.

###  II.1 The Stokes and lubrication equations

The formulation to be presented below will disregard the relative effects of
liquid inertia and gravity. On the one hand, the relative importance of
inertia may be quantified through the Laplace number $\mbox{{La}}=\rho\sigma
h_{o}/\mu^{2}$, which is a Reynolds number based on the visco-capillary
velocity $v_{\sigma}=\sigma/\mu$, and is readily seen to be much smaller than
unity in typical configurations of very viscous liquids coating a fibre. On
the other hand, gravity effects may also be neglected upon comparing the
velocity $v_{g}=\rho gh_{o}^{2}/\mu$, stemming from a balance between viscous
and gravitational forces, with the velocity $v_{\mathrm{vdW}}=|A|/(6\pi\mu
h_{o}^{2})$ arising from a molecular-viscous balance, which yields the
rescaled Bond number $v_{g}/v_{\mathrm{vdW}}=6\pi\rho gh_{o}^{4}/|A|$.
Therefore, the governing equations of motion under the double limit $La\ll 1$
and $v_{g}/v_{\mathrm{vdW}}\ll 1$111 For a fibre coated by a
$1000\leavevmode\nobreak\ \mathrm{cSt}$ silicone oil film of initial thickness
$h_{o}=100\leavevmode\nobreak\ \mathrm{nm}$, $\rho=970\leavevmode\nobreak\
\mathrm{kg}/\mathrm{m}^{3}$, $\mu=0.97\leavevmode\nobreak\
\mathrm{Pa}\cdot\mathrm{s}$, $\sigma=27\times 10^{-3}\mathrm{N}/m$ we get
$v_{g}/v_{\mathrm{vdW}}=1.7\times 10^{-4}$ and $La=2.8\times 10^{-6}$. are
rendered dimensionless upon taking $h_{o}$, $\sigma/\mu$, $\mu h_{o}/\sigma$,
and $\sigma/h_{o}$ as the relevant scales for length, velocity, time, and
pressure, respectively, and read

$\bm{\nabla}\bm{\cdot}\bm{u}=0,\quad\bm{\nabla}\bm{\cdot}\mathsf{T}=\bm{0},$
(1)

where $\bm{u}=u\bm{e}_{x}+v\bm{e}_{r}+w\bm{e}_{\theta}$ is the liquid velocity
field,
$\mathsf{T}=-p\mathsf{I}+\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{\rm{T}}$ is
the liquid stress tensor, $p$ is the gauge liquid pressure and $\mathsf{I}$ is
the identity tensor. At the liquid-air interface, $r=\eta+h$, where
$\eta=R/h_{o}$ is the ratio of the cylinder radius to the unperturbed film
thickness and $h$ now stands for the dimensionless film thickness, both the
kinematic boundary condition,

$\frac{\partial h}{\partial t}=\bm{u}\cdot\bm{\nabla}(r-h),$ (2)

and the stress boundary condition,

$-\mathsf{T}\bm{\cdot}\bm{n}=\bm{n}(\bm{\nabla}\bm{\cdot}\bm{n})+\phi(h)\bm{n},$
(3)

are prescribed, where $\bm{n}=\bm{\nabla}(r-\eta-h)/|\bm{\nabla}(r-\eta-h)|$
is the outward-pointing normal unit vector, $\phi(h)=\mathcal{A}/h^{3}$ is the
dimensionless disjoining pressure, and $\mathcal{A}=A/(6\pi\sigma h_{o}^{2})$
is a dimensionless Hamaker constant that compares the relative importance of
vdW and surface tension forces. Note also that
$\mathcal{A}=\text{sgn}(A)(h_{o}/a)^{-2}$, where $a=\sqrt{|A|/(6\pi\sigma)}$
is the molecular length.[29, 31, 4]. Finally, no-slip conditions are imposed
at the surface of the solid cylinder $r=\eta$, where $\bm{u}=0$.

The lubrication model, derived in Appendix A, yields the following partial
differential equation for the interface position $h(x,\theta,t)$,

$\displaystyle(\eta+h)\frac{\partial h}{\partial t}+\frac{\partial}{\partial
x}\left(G\frac{\partial p}{\partial
x}\right)+\frac{\partial}{\partial\theta}\left(H\frac{\partial
p}{\partial\theta}\right)=0,$ (4)

where the mobility functions $G=G(\eta,h)$ and $H=H(\eta,h)$, derived in
Appendix A, are given by

$\displaystyle
G(h,\eta)=\frac{h}{16}\left(h+2\eta\right)\left(3h^{2}+6h\eta+2\eta^{2}\right)+\frac{1}{4}(h+\eta)^{4}\log\frac{\eta}{h+\eta},$
(5) $\displaystyle
H(h,\eta)=-\frac{h(h+2\eta)(h^{2}+2h\eta+2\eta^{2})}{8\eta^{2}}+\frac{1}{2}(h+\eta)^{2}\log\frac{h+\eta}{\eta},$
(6)

while the pressure is given by the normal-stress balance

$p={\mathcal{C}}+\mathcal{A}{h^{-3}},$ (7)

being $\mathcal{C}$ twice the mean curvature of the film (see Appendix A for
details). The recent study by Zhao et al. [35] on annular thin films coating a
fibre derived a similar model but in the absence of molecular forces,
$\mathcal{A}=0$, and including the effect of wall slip. The validity of the
lubrication model presented above relies in the approximation of the flow
being slender, $\epsilon=h_{o}/\lambda\ll 1$, where $\lambda$ is the
characteristic wavelength of the disturbances. We may anticipate that
$\lambda\simeq 2\pi h_{o}/\sqrt{3\mathcal{A}}$ is the dimensional cut-off
wavelength [29, 32]. Since the aspect ratio $\eta=R/h_{o}$ satisfies
$\epsilon=h_{o}/\lambda=h_{o}R/(R\lambda)=\eta^{-1}R/\lambda\sim\mathcal{A}^{1/2}/\eta^{2}$,
the lubrication approximation is expected to break down for combinations of
$\eta$ and $\mathcal{A}$ such that $\epsilon$ is no longer small, i.e. for
$\eta\simeq O(1)$ and $\mathcal{A}\lesssim O(1)$, as will be evidenced below.

Two dimensionless parameters govern the flow at hand, namely the dimensionless
Hamaker constant $\mathcal{A}$, and the ratio of the fibre radius to the
initial film thickness, $\eta$.

##  III Linear stability analysis

###  III.1 The unperturbed state

The base state consists of a static annular film of uniform thickness $h_{o}$
that depends on the coating process leading to the formation of the film, and
may in principle take any value in the context of the stability analysis
performed herein. The corresponding dimensional base pressure is given by
$p_{o}=\sigma/(R+h_{o})+A/(6\pi h_{o}^{3})$. In the non-wetting case, $A>0$,
no stable equilibrium is possible, since both surface tension and vdW forces
contribute to destabilise the interface. However, as shown by Quéré et al.
[26], a stable equilibrium is possible in the wetting case, $A<0$, given by a
balance between the surface tension and vdW forces, whereby films with
thicknesses smaller than $\sqrt{aR}$ are stable.

###  III.2 Critical conditions and cut-off wavenumber

Following the Plateau-like argument given e.g. by Scheludko [29] [see also 31,
32], in which periodic undulations of the free surface with negative extra
Gibbs energy tend to grow in time, let us consider the dimensionless excess
Gibbs energy $\Delta G$ for our particular configuration in which the initial
perturbation to the fibre surface may be written as

$r=r_{o}\equiv 1+\eta+\delta\cos\left(kx+m\theta\right),$ (8)

where $\delta\ll 1$ is the amplitude of the undulations. To that end, let
$\Delta G=\Delta G_{s}+\Delta G_{\textrm{molec}}$, where $\Delta G_{s}$ is the
Gibbs energy associated to the extra surface area due to the perturbation and
$\Delta G_{\textrm{molec}}$ is the extra interaction energy due to the
molecular interactions, which take the form

$\displaystyle\Delta
G_{s}=\int_{L}\int^{2\pi}_{0}\left|\bm{X}_{\theta}\wedge\bm{X}_{\varphi}\right|\,\mathrm{d}\theta\mathrm{d}x-S_{o},$
(9) $\displaystyle\Delta
G_{\textrm{molec}}=\int_{L}\int^{2\pi}_{0}\left[\phi(1)\Delta
h+\phi^{\prime}(1)(\Delta
h)^{2}/2+\ldots\right](1+\eta)\,\mathrm{d}\theta\mathrm{d}x,$ (10)

where $\bm{X}=\left(r_{o}\cos\theta,r_{o}\sin\theta,x\right)$ is a
parameterization of the corrugated surface, $S_{o}=2\pi(1+\eta)L$ is the
unperturbed interface area, $\Delta h=\delta\cos\left(kx+m\theta\right)$, and
the axial length $L$ is to be defined below. Note that
$\int_{h}\phi(h)\,\mathrm{d}h$ is the energy of interaction per unit area
among the molecules in the film, assumed to be essentially determined by vdW
forces. The incompressibility condition

$\Delta\mathcal{V}=\int_{L}\int^{2\pi}_{0}\int^{r_{o}}_{\eta}r\,\mathrm{d}r-V_{o}=0,$
(11)

where $V_{o}=\pi(1+2\eta)L$ is the unperturbed volume of the liquid fibre,
yields the following restriction for the first- and second-order terms

$(1+\eta)\delta\int_{L}\int^{2\pi}_{0}\cos(kx+m\theta)\,\mathrm{d}\theta\mathrm{d}x=-\frac{1}{2}\delta^{2}\int_{L}\int^{2\pi}_{0}\cos^{2}(kx+m\theta)\,\mathrm{d}\theta\mathrm{d}x.$
(12)

Introducing this restriction onto (9) and letting $L=2\pi/k$ be a wavelength
produces

$\Delta G=\Delta G_{s}+\Delta
G_{\textrm{molec}}=\frac{\delta^{2}\pi^{2}}{k(1+\eta)}\left[{k^{2}(1+\eta)^{2}+m^{2}-1}+(1+\eta)^{2}\phi^{\prime}(1)\right]+O(\delta^{4})$
(13)

Following Scheludko [29], Vrij [31], Vrij and Overbeek [32], any perturbation
is allowed to grow exponentially if the total excess energy $\Delta G$ becomes
negative. The cutoff is thus found when $\Delta G=0$ or, equivalently, to
second order,

$k^{2}(1+\eta)^{2}+m^{2}-1+(1+\eta)^{2}\phi^{\prime}(1)=0,$ (14)

for which the cutoff wavenumber $k_{c}$ is obtained as the positive solution
to (14), i.e.,

$k_{c}=\sqrt{-\phi^{\prime}(1)+\frac{1-m^{2}}{\left(1+\eta\right)^{2}}},$ (15)

where $\phi^{\prime}(1)=-3\mathcal{A}$ for the vdW potential employed
throughout this work. In the wetting case, $\mathcal{A}<0$, only the
axisymmetric modes $m=0$ are unstable, and $k_{c}>0$ for
$3|\mathcal{A}|<(1+\eta)^{-2}$. Thus, the film is stable if
$3|\mathcal{A}|>(1+\eta)^{-2}$, which corresponds to initial thicknesses
smaller than the Quéré thickness $h_{c}$, which is given by the algebraic
equation $h_{c}^{2}-\sqrt{3}\,a\,(R+h_{c})=0$ whose solution is
$h_{c}=3^{1/4}\sqrt{aR}+O(a)$ [26]. It is noteworthy that the same stability
condition as in Quéré et al. [26] is obtained, given that their argument is
based on the condition that the pressure is an increasing function of the film
thickness, but disregards the possibility of small-amplitude harmonic
disturbances.

In the non-wetting case, $\mathcal{A}>0$, both axisymmetric and non-
axisymmetric modes are potentially unstable. In particular, the axisymmetric
and first helical modes, $m=0,\pm 1$ are unstable for any values of
$\mathcal{A}$ and $\eta$, while the higher-order azimuthal modes $|m|\geq 2$
are unstable for $\mathcal{A}>(m^{2}-1)(1+\eta^{2})^{-2}/3$. Since different
modes compete in a generic low-noise scenario, it is crucial to elucidate the
dominant one depending on the values of $\mathcal{A}$ and $\eta$, which
requires obtaining the maximum growth rate from the dispersion relation.

###  III.3 The dispersion relation

The dispersion relation is obtained by introducing the usual normal-mode
decomposition

$(\bm{u},p,h)=\left(0,\frac{p_{o}}{\sigma/h_{o}},1\right)+\delta\left[{\hat{\bm{u}\mkern
3.0mu}\mkern-3.0mu}{}(r),{\hat{p\mkern 3.0mu}\mkern-3.0mu}{}(r),{\hat{h\mkern
3.0mu}\mkern-3.0mu}{}\right]\mathrm{e}^{\mathrm{i}(kx+m\theta)+\omega t},$
(16)

into the corresponding equations of motion and linearising about the base
flow, where $\delta$ is an arbitrarily small quantity, hats are used for
eigenfunctions, $\omega$ is the dimensionless growth rate, $|m|=0,1,\ldots$
and $k\in\mathcal{R}$ is the dimensionless wavenumber. In the case of the
Stokes equations (1), a closed-form solution for the growth rate
$\omega=\omega(k,\mathcal{A},\eta,m)$ can be found analytically in terms of
modified Bessel functions. The general process to arrive at such expression is
provided in Appendix B. In the case of the lubrication model (4), the solution
reads

$\omega=\frac{\varphi_{1}(k,\eta,m)+\varphi_{2}(k,\eta,m)}{16\,(1+\eta)}\underbrace{\left[k^{2}+\frac{m^{2}-1}{(1+\eta)^{2}}-3\mathcal{A}\right]}_{\textrm{cutoff
function}},$ (17)

where the auxiliary functions

$\displaystyle\varphi_{1}={(2\eta+1)}\left\\{[2\eta\,(\eta+3)+3]\,k^{2}-\frac{2m^{2}}{\eta^{2}}\,[2\eta\,(\eta+1)+1]\right\\},$
(18)
$\displaystyle\varphi_{2}=8(1+\eta)^{2}\left[m^{2}-\frac{k^{2}}{2}(1+\eta)^{2}\right]\log\frac{1+\eta}{\eta},$
(19)

have been defined. A large-$\eta$ expansion for the growth rate (17) is
available, namely,

$\displaystyle\omega=\,$ $\displaystyle
k^{2}\left(\mathcal{A}-\frac{k^{2}}{3}\right)+\frac{k^{2}\left(9\mathcal{A}-40\,m^{2}+20\right)+60\mathcal{A}\,m^{2}-3k^{4}}{60\,\eta^{2}}+O(\eta^{-3}),$
(20)

which, at leading order, is independent of the azimuthal number $m$ and
coincides with the classical result for a planar film [31]. Three-dimensional
effects appear in equation (20) in the next-to-leading, $O(\eta^{-2})$ term.
Also, in an effort to elucidate the large-$\eta$ behaviour of the dominant
mode, given by the pair $(k_{M},\omega_{M})$, one can find the solution to the
condition $\partial_{k}\omega=0$ from Eq. (17), wherefrom it is found

$\omega_{M}=\frac{3\mathcal{A}^{2}}{4}+\frac{9\mathcal{A}^{2}+40\mathcal{A}}{80\eta^{2}}+\frac{120\mathcal{A}m^{2}-21\mathcal{A}^{2}-160\mathcal{A}}{160\eta^{3}}+O(\eta^{-4}).$
(21)

###  III.4 The dominant azimuthal mode of non-wetting films, $\mathcal{A}>0$

The analysis of the Stokes and lubrication dispersion relations reveals that,
in both cases, there is only one relevant temporal branch for each combination
of $(\mathcal{A},\eta,m)$.

Figure 2: Temporal stability branches for $\mathcal{A}=0.01$ and $\eta=1$
(left panel), $\eta=10$ (central panel), and $\eta=20$ (right panel), for the
azimuthal modes $m=(0,\ldots,3)$. Solid and dashed lines correspond to the
Stokes and lubrication results, respectively.

Figure 2 shows the temporal branches for $\mathcal{A}=0.01$, $\eta=(1,10,20)$
in the right, middle and left panels, respectively, and different values of
$m=(0,1,2,3)$. The results obtained using the Stokes dispersion relation are
represented with solid lines, while those extracted from the lubrication model
(17) are plotted with dashed lines. Note first that only for relatively thick
films, illustrated by the case $\eta=1$, is there a significant quantitative
mismatch between the Stokes equations and lubrication theory, the latter
overestimating the growth rate considerably. In contrast, both results are
nearly identical for the cases $\eta=(10,20)$, which represent films much
thinner than the fibre radius. It is thus deduced that lubrication theory
provides a good quantitative description of the linearised dynamics of the
flow except for cases where $\eta\lesssim 1$. Secondly, as explained in
§III.2, both the axisymmetric, $m=0$, and first helical, $|m|=1$, modes are
unconditionally unstable, while azimuthal modes of increasingly higher order
become unstable for increasing values of $\eta$. Note that, according to the
Stokes equations, the axisymmetric mode is dominant for $\eta=1$, in that its
associated maximum growth rate is larger than that associated to the first
helical mode. In contrast, lubrication theory predicts the dominance of the
latter mode over the axisymmetric one, as evidenced in the left panel of
figure 2, thus indicating that lubrication theory must be used with great care
when dealing with relatively thick films. In addition, this result suggests
that lubrication theory may also fail when applied to the nonlinear flow
regimes that take place for large times. However, this issue shall not be
clarified in the present contribution, which is focused on the linearised
dynamics only. For a relatively thinner film with $\eta=10$, the $|m|=1$ mode
is the dominant one according to both the Stokes and lubrication dispersion
relations, and the $|m|=2$ mode is also unstable, although with a
comparatively smaller maximum growth rate. This trend continues as the film
becomes thinner, as can be seen in the left panel, which for $\eta=20$ reveals
that the dominant mode is $|m|=2$, and that the $|m|=3$ mode has a growth rate
comparable to that of the axisymmetric and $|m|=1$ ones.

Figure 3: Maximum growth rate $\omega_{M}$ (left panel) and corresponding
wavenumber $k_{M}$ (right panel) for $\mathcal{A}=0.01$ computed according to
the Stokes (blue dots) and lubrication (red dots) dispersion relations. The
associated azimuthal mode is indicated explicitly.

Figure 4: Contour of maximum growth rate (left panel) and corresponding wave
number (right panel) over the parametric plane ($\eta$,$\mathcal{A}$)
predicted by the linearised Stokes equations. The boundaries between different
dominant azimuthal modes are indicated with solid lines, with the associated
values of $m$ indicated explicitly.

Figure 3, which represents the maximum growth rate $\omega_{M}(\eta)$ (left
panel) and the associated axial wavenumber $k_{M}(\eta)$ (right panel) for
$\mathcal{A}=0.01$, clearly shows that the order of the dominant azimuthal
mode (coloured vertical regions) increases monotonically with $\eta$, starting
with $m=0$ for $0<\eta<1.8$, $|m|=1$ for $1.8<\eta<12.3$, $|m|=2$ for
$12.3<\eta<20.9$ and $|m|=4$ for $29.3<\eta<37.4$, the trend continuing for
larger values of $\eta$. It is noteworthy that lubrication theory (red dots)
fails to capture the transition between the axisymmetric and first azimuthal
mode, since it predicts that the dominant mode is always $|m|=1$ for
$\eta\lesssim 12$. The parametric behaviour of the optimal wavenumber $k_{M}$,
represented in the right panel of figure 3, reveals an interesting non-trivial
trend. For values of $\eta<1.8$ the dominant axisymmetric mode has a nonzero
associated axial wavenumber with values in the range $0.27<k_{M}<0.35$.
However, when the first azimuthal mode crosses over and becomes dominant, the
value of $k_{M}$ switches to zero in the range $1.8<\eta<5.6$, and becomes
strictly positive and increasing as a function of $\eta$ in the remaining
range, $5.6<\eta<12.3$. This behaviour of $k_{M}$ is seen to repeat each time
that a transition between the dominant azimuthal modes takes place.

The results of figure 3 correspond to a fixed value of the dimensionless
Hamaker constant, $\mathcal{A}=0.01$. To assess the role played by
$\mathcal{A}$, figure 4 shows contours of $\log_{10}(\omega_{M})$ (left panel)
and $k_{M}$ (right panel) in the $(\eta,\mathcal{A})$ parameter plane. The
examination of lines of constant $\mathcal{A}$ reveals that the trend
associated with the case $\mathcal{A}=0.01$ of figure 3, and discussed in the
previous paragraph, applies to any value of $\mathcal{A}$, except that the
transition to higher-order azimuthal modes takes place for smaller values of
$\eta$, i.e. relatively thicker films, when the value of $\mathcal{A}$
increases, i.e. for stronger van der Waals forces relative to surface tension
forces. Similarly, the transitions take place for larger values of $\eta$ when
the value of $\mathcal{A}$ is decreased. In particular, the boundaries between
the different azimuthal modes, represented with solid black lines, have
complicated stair-like shapes due to the branch switching behaviour explained
in the context of figure 3. In addition, the dominant wavenumber is zero in
the region surrounding each transition curve, a result that can be understood
in the light of the transitions shown in figure 3 for $\mathcal{A}=0.01$.
Indeed, except for the axisymmetric mode $m=0$ which always has a value of
$k_{M}>0$, all the non-axisymmetric modes have an associated value of
$k_{M}=0$ when they become dominant at each transition. It should be noted
that the numerical campaign to produce figure 4 only reaches $m=5$. Thus, the
parametric region corresponding to these higher-order modes is not explored
but its trend can be easily conceptualized.

##  IV Concluding remarks

We have presented a detailed analysis of the linear stability of static thin
annular films of liquid deposited on cylindrical fibres and subjected to van
der Waals interactions, using both the full Stokes equations of motion and a
leading-order lubrication description. In the wetting case only axisymmetric
modes are unstable, but only when the film thickness is larger than thickness
deduced by Quéré [25]. In the non-wetting case, our analysis has revealed that
both the axisymmetric and the first helical modes are unconditionally
unstable, while azimuthal modes of increasing order with $|m|\geq 2$ become
unstable when the film thickness is sufficiently smaller than the fibre
radius. The transition between dominant modes has been mapped in the parameter
plane spanned by the dimensionless strength of the van der Waals forces and
the film-to-fibre radius ratio, which are the only dimensionless parameters
governing the flow, revealing that the boundaries separating the dominant
modes present a surprisingly complicated structure in experimentally relevant
parametric ranges.

It is important to note that the maximum growth rates associated to the
different azimuthal modes are comparable, especially for large values of
$\eta$. For instance, it is deduced from the right panel of figure 2 that,
although the dominant mode is $m=2$ for $\eta=20$ and $\mathcal{A}=0.01$, the
values of $\omega_{M}$ are quite similar for all the unstable modes
$m\in\\{0,1,2,3\\}$. This fact may question the validity of the linear
stability analysis in predicting the resulting nonlinear dynamics and
dewetting patterns. However, it should be noticed that a small difference in
the growth rate of two different normal modes becomes large prior to the
nonlinear stages of the flow, provided that the initial amplitude of the noise
is small enough. Indeed, let $\omega_{M}^{(1)}>\omega_{M}^{(2)}$ be the
maximum growth rates of two different modes that share the same initial
amplitude $\epsilon$. After a time $t_{L}$ the relative amplitude is then
$\exp[(\omega_{M}^{(1)}-\omega_{M}^{(2)})t_{L}]$, which is of order unity
provided that $\omega_{M}^{(1)}-\omega_{M}^{(2)}\sim t_{L}^{-1}$. The key
factor for the observation of the dominant mode during the nonlinear regime is
then the time $t_{L}$ during which the linear regime prevails, which can be
estimated in terms of the initial amplitude $\epsilon\ll 1$ as
$\epsilon\exp(\omega_{M}^{(1)}t_{L})\sim 1\Rightarrow
t_{L}\sim\log(\epsilon^{-1})/\omega_{M}^{(1)}$. Consequently, the relative
amplitude at the end of the linear stage is of order unity provided that
$(\omega_{M}^{(1)}-\omega_{M}^{(2)})/\omega_{M}^{(1)}\sim(\log(\epsilon^{-1}))^{-1}\ll
1$. Thus, the smaller the amplitude of the initial noise, the smaller the
relative difference between the growth rates of the dominant and subdominant
modes that is needed to reach order-unity relative amplitudes at the end of
the linear stage.

As a final remark, our findings indicate the need for further theoretical and
experimental explorations. Based on our conclusions, non-axisymmetric modes
should be observed in annular dewetting experiments, and to the best of our
knowledge, this phenomenon has not been documented in prior studies. In
addition, a study of the nonlinear and large-time dynamics of the film is
mandatory, since nontrivial three-dimensional dewetting patterns are expected
to arise as a consequence of the non-axisymmetric stability modes reported
herein.

###### Acknowledgements.

This work has been supported by the Spanish MINECO under project
PID2020-115655GB-C22, partly financed through FEDER European funds. The
authors thank Dr. A. Martínez-Calvo for fruitful discussions.

Declaration of Interests
The authors report no conflict of interest.

##  Appendix A Derivation of the leading-order lubrication model

We present here the derivation of the Reynolds equation for the evolution of
the annular film, where the main validity assumption relies on the parameter
$\epsilon=h_{o}/\lambda=\mathcal{A}^{1/2}/\eta^{2}$ being small. We first
recall that the kinematic condition at $r=\eta+h$ reads

$\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial
x}+\frac{w}{r}\frac{\partial h}{\partial\theta}=v.$ (22)

Integration of the continuity equation across the liquid layer produces

$\int^{\eta+h}_{\eta}\left(r\frac{\partial u}{\partial
x}+\frac{\partial}{\partial r}\left(rv\right)+\frac{\partial
w}{\partial\theta}\right)\,\mathrm{d}r=0,$ (23)

which upon application of the Leibniz integral rule and use made of (22)
yields

$(\eta+h)\frac{\partial h}{\partial t}+\frac{\partial}{\partial
x}\int^{\eta+h}_{\eta}r\frac{\partial u}{\partial
x}\,\mathrm{d}r+\frac{\partial}{\partial\theta}\int^{\eta+h}_{\eta}\frac{\partial
w}{\partial\theta}\,\mathrm{d}r=0.$ (24)

The velocities $u$ and $w$ are found from integration of the momentum equation
under the assumption of negligible streamwise and azimuthal viscous diffusion,
respectively, and negligible radial pressure gradients compared with the axial
and azimuthal ones, as a result of the slender approximation. This yields222 A
fully detailed asymptotic expansion of the momentum equation and the
corresponding boundary conditions is available in the work by Zhao et al. [35]

$\frac{\partial p}{\partial x}=\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial u}{\partial r}\right),\quad\frac{1}{r}\frac{\partial
p}{\partial\theta}=\frac{\partial}{\partial
r}\left(\frac{1}{r}\frac{\partial}{\partial r}(rw)\right),$ (25)

together with the no-slip boundary condition $u=w=0$ at the cylinder wall
$r=\eta$ and vanishing leading-order tangential stresses
$\tau^{\prime}_{xr}=\partial_{r}u=0$ and $\tau^{\prime}_{\theta
r}=\partial_{r}w-w/r=0$ at the free surface $r=\eta+h$, to yield

$u=\frac{1}{4}\frac{\partial p}{\partial
x}\left[r^{2}-\eta^{2}-2(h+\eta)^{2}\log\frac{r}{\eta}\right],\,\,\,w=\frac{1}{4}\frac{\partial
p}{\partial\theta}\left[2r\log\frac{r}{\eta}+\frac{(h+\eta)^{2}(\eta^{2}-r^{2})}{r\eta^{2}}\right].$
(26)

Inserting these profiles into (24) produces the following evolution equation
for $h$

$\displaystyle(\eta+h)\frac{\partial h}{\partial t}+\frac{\partial}{\partial
x}\left(G\frac{\partial p}{\partial
x}\right)+\frac{\partial}{\partial\theta}\left(H\frac{\partial
p}{\partial\theta}\right)=0,$ (27)

where the mobility functions $G$ and $H$ are given by

$\displaystyle
G(h,\eta)=\frac{h}{16}\left(h+2\eta\right)\left(3h^{2}+6h\eta+2\eta^{2}\right)+\frac{1}{4}(h+\eta)^{4}\log\frac{\eta}{h+\eta},$
(28) $\displaystyle
H(h,\eta)=-\frac{h(h+2\eta)(h^{2}+2h\eta+2\eta^{2})}{8\eta^{2}}+\frac{1}{2}(h+\eta)^{2}\log\frac{h+\eta}{\eta}.$
(29)

The specification of the pressure through the leading-order normal stress
balance completes the description, i.e.,

$p={\mathcal{C}}+\mathcal{A}{h^{-3}},$ (30)

where $\mathcal{C}$ is (twice) the signed curvature of the liquid-air
interface,

$\mathcal{C}=\bm{\nabla}\cdot\left(\frac{\bm{\nabla}F}{|\bm{\nabla}F|}\right),\quad\textrm{with}\quad
F=r-\eta-h(x,\theta,t)=0,$ (31)

which, expressed in terms of $h(x,\theta,t)$, yields

$\mathcal{C}=\,\frac{1}{n(\eta+h)}-\frac{h_{xx}+h_{\theta\theta}/(\eta+h)^{2}}{n}+\frac{1}{n^{3}}\left[h_{x}^{2}h_{xx}+\frac{2h_{x}h_{\theta}h_{x\theta}}{(\eta+h)^{2}}+\frac{h_{\theta}^{2}\left(\eta+h+h_{\theta\theta}\right)}{(\eta+h)^{4}}\right]$
(32)

with $n=\sqrt{1+h_{x}^{2}+h_{\theta}^{2}/(1+\eta)^{2}}$ and subscripts now
indicating partial derivative.  We note that the mobility function $G$ admits
the large-$\eta$ expansion $G/\eta\simeq-h^{3}/3+O(\eta^{-1})$ and that the
leading-order curvature reads
$\mathcal{C}\simeq{1}/{\eta}-{h}/{\eta^{2}}-h_{xx},$ for $h\ll\eta$ and
$h_{x}\ll 1$. Thus, the model (27) reduces to

$\frac{\partial h}{\partial t}+\frac{\partial}{\partial
x}\left[\frac{h^{3}}{3}\left(h_{xxx}-\frac{h_{x}}{\eta^{2}}-3\mathcal{A}h_{x}h^{-4}\right)\right]=0$
(33)

which coincides with the axisymmetric models presented in Hammond [12] or
Lister et al. [20] for $\mathcal{A}=0$ up to a rescaling of the
$x$-coordinate.

##  Appendix B The Stokes dispersion relation

Introduction of the normal-mode decomposition into the Stokes equations (1)
furnishes the linear problem

$\displaystyle 0$
$\displaystyle=\mathrm{i}k\hat{u}+\frac{1}{r}\frac{\partial}{\partial
r}(r\hat{v})+\frac{\mathrm{i}m}{r}\hat{w},$ (34)
$\displaystyle\mathrm{i}k{\hat{p\mkern 3.0mu}\mkern-3.0mu}{}$
$\displaystyle=\bm{\nabla}^{2}_{k,m}\hat{u},$ (35)
$\displaystyle\frac{\partial{\hat{p\mkern 3.0mu}\mkern-3.0mu}{}}{\partial r}$
$\displaystyle=\bm{\nabla}^{2}_{k,m}\hat{v}-\frac{\hat{v}}{r^{2}}-\frac{2\mathrm{i}m}{r^{2}}\hat{w},$
(36) $\displaystyle\frac{\mathrm{i}m}{r}{\hat{p\mkern 3.0mu}\mkern-3.0mu}{}$
$\displaystyle=\bm{\nabla}^{2}_{k,m}\hat{w}-\frac{\hat{w}}{r^{2}}+\frac{2\mathrm{i}m}{r^{2}}\hat{v},$
(37)

where $\bm{\nabla}^{2}_{k,m}=\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial r}\right)-(k^{2}+m^{2}/r^{2})$ is the
normal-mode Laplacian. Since the perturbed velocity is solenoidal, it follows
that ${\hat{p\mkern 3.0mu}\mkern-3.0mu}{}$ satisfies the normal-mode Laplace
equation

$\bm{\nabla}^{2}_{k,m}{\hat{p\mkern 3.0mu}\mkern-3.0mu}{}=0,$ (38)

so it can be resolved into modified Bessel functions of the form

${\hat{p\mkern 3.0mu}\mkern-3.0mu}{}=A_{1}I_{m}(kr)+A_{2}K_{m}{(kr)},$ (39)

where the $A_{i}$’s are constants of integration, and $I_{m}$ and $K_{m}$ are
the modified Bessel functions of first and second kind of order $m$,
respectively. Consequently, the most general solution to (35) is

$\hat{u}=A_{3}I_{m}(kr)+A_{4}K_{m}(kr)+\frac{\mathrm{i}r}{2}\left[A_{1}I_{m+1}(kr)-A_{2}K_{m-1}(kr)\right].$
(40)

Upon introduction of (39) into linear combinations of (36) and (37) and
leveraging the solenoidality condition (34), it follows that

$\displaystyle\hat{v}=$
$\displaystyle-\frac{A_{1}}{2}\left[\frac{I_{m+1}(kr)}{k}-\frac{r}{2}\left\\{I_{m+2}(kr)+I_{m}(kr)\right\\}\right]-$
$\displaystyle\frac{A_{2}}{2}\left[\frac{2mK_{m+1}(kr)-K_{m+1}(kr)}{k}-\frac{r}{2}\left\\{K_{m+2}(kr)+K_{m}(kr)\right\\}\right]-$
$\displaystyle\mathrm{i}A_{3}I_{m+1}(kr)+\mathrm{i}A_{4}K_{m+1}(kr)+\frac{\mathrm{i}A_{5}}{2}\left[I_{m+1}(kr)-I_{m-1}(kr)\right]+$
$\displaystyle\frac{\mathrm{i}A_{6}}{2}\left[K_{m+1}(kr)-K_{m-1}(kr)\right],$
(41)

and

$\displaystyle\hat{w}=$
$\displaystyle\frac{\mathrm{i}A_{1}}{2}\left[\frac{I_{m+1}(kr)}{k}-\frac{r}{2}\left\\{I_{m+2}(kr)-I_{m}(kr)\right\\}\right]+$
$\displaystyle\frac{\mathrm{i}A_{2}}{2}\left[\frac{2mK_{m+1}(kr)-K_{m+1}(kr)}{k}-\frac{r}{2}\left\\{K_{m+2}(kr)-K_{m}(kr)\right\\}\right]-$
$\displaystyle
A_{3}I_{m+1}(kr)+A_{4}K_{m+1}(kr)+\frac{A_{5}}{2}\left[I_{m-1}(kr)+I_{m+1}(kr)\right]+$
$\displaystyle\frac{A_{6}}{2}\left[K_{m-1}(kr)+K_{m+1}(kr)\right].$ (42)

The six constants of integration $A_{1},\ldots,A_{6}$ need to be found upon
prescription of the linearised boundary conditions, which in normal-mode form
read

$\hat{u}=\hat{v}=\hat{w}=0,\quad\textrm{at}\quad r=\eta,$ (43)

and

$\displaystyle-$ $\displaystyle{\hat{p\mkern
3.0mu}\mkern-3.0mu}{}+2\hat{v}^{\prime}+\left(k^{2}-\frac{1-m^{2}}{(1+\eta)^{2}}-3\mathcal{A}\right){\hat{h\mkern
3.0mu}\mkern-3.0mu}{}=0,$ (44)
$\displaystyle\hat{u}^{\prime}+\mathrm{i}k\hat{v}=0,$ (45)
$\displaystyle\hat{w}^{\prime}-\hat{w}/(1+\eta)+\mathrm{i}m\hat{v}/(1+\eta)=0,$
(46) $\displaystyle\omega{\hat{h\mkern 3.0mu}\mkern-3.0mu}{}=\hat{v},$ (47)

at $r=1+\eta$, where primes indicate radial derivatives. Thus, the problem may
be recast into a matrix system of the form $\mathsf{M}\bm{\cdot}\bm{A}=0$,
where $\bm{A}=(A_{1},\ldots,A_{6})^{\mathrm{T}}$ and
$\mathsf{M}(\omega,k,\mathcal{A},\eta,m)$ is a $6\times 6$ matrix. The entries
$M_{ij}$ of the matrix $\mathsf{M}$ are listed here for completeness:

$\displaystyle M_{11}=\frac{\mathrm{i}}{2}\eta I_{m+1}(k\eta),\
M_{12}=-\frac{\mathrm{i}}{2}\eta K_{1-m}(k\eta),\ M_{13}=I_{m}(\eta k),\
M_{14}=K_{m}(\eta k),\ M_{15}=M_{16}=0,$ $\displaystyle
M_{21}=\frac{\eta}{4}I_{m}(k\eta)+\frac{\eta}{4}I_{m+2}(k\eta)-\frac{I_{m+1}(k\eta)}{2k},\
M_{22}=\frac{\eta}{4}K_{m}(k\eta)+\frac{\eta}{4}K_{m+2}(k\eta)-\frac{(2m-1)K_{m+1}(k\eta)}{2k},$
$\displaystyle M_{23}=-iI_{m+1}(k\eta),\ M_{24}=iK_{m+1}(k\eta),$
$\displaystyle M_{25}=\frac{1}{2}iI_{m+1}(k\eta)-\frac{1}{2}iI_{m-1}(k\eta),\
M_{26}=\frac{1}{2}iK_{m+1}(k\eta)-\frac{1}{2}iK_{m-1}(k\eta),$ $\displaystyle
M_{31}=\frac{\mathrm{i}}{4}\eta
I_{m}(k\eta)+\frac{iI_{m+1}(k\eta)}{2k}-\frac{\mathrm{i}}{4}\eta
I_{m+2}(k\eta),$ $\displaystyle M_{32}=\frac{\mathrm{i}}{4}\eta
K_{m}(k\eta)+\frac{i(2m-1)K_{m+1}(k\eta)}{2k}-\frac{\mathrm{i}}{4}\eta
K_{m+2}(k\eta),$ $\displaystyle M_{33}=-I_{m+1}(k\eta),\ M_{34}=K_{m+1}(\eta
k),$ $\displaystyle M_{35}=\frac{1}{2}(I_{m-1}(k\eta)+I_{m+1}(k\eta)),\
M_{36}=\frac{1}{2}K_{m-1}(k\eta)+\frac{1}{2}K_{m+1}(k\eta),$ $\displaystyle
M_{41}=\frac{1}{4}\left[\phi_{1}+\frac{\phi_{2}I_{m}(\tilde{k})+\phi_{3}I_{m+1}(\tilde{k})}{k(1+\eta)^{2}\omega}\right],\
M_{42}=\frac{1}{4}\left[-\phi_{1}+\frac{\phi_{4}K_{m}(\tilde{k})+\phi_{5}K_{m+1}(\tilde{k})}{k(1+\eta)^{2}\omega}\right],$
$\displaystyle
M_{43}=-\mathrm{i}\left[kI_{m}(\tilde{k})+kI_{m+2}(\tilde{k})+\frac{\left\\{m^{2}-1-(3\mathcal{A}-k^{2})(\eta+1)^{2}\right\\}I_{m+1}(\tilde{k})}{(\eta+1)^{2}\omega}\right],$
$\displaystyle
M_{44}=-\mathrm{i}\left[kK_{m}(\tilde{k})+kK_{m+2}(\tilde{k})-\frac{\left\\{m^{2}-1-(3A-k^{2})(\eta+1)^{2}\right\\}K_{m+1}(\tilde{k})}{(\eta+1)^{2}\omega}\right],$
$\displaystyle
M_{45}=-\frac{\mathrm{i}m}{k(1+\eta)^{3}\omega}\left[2(\eta+1)^{2}k\omega
I_{m-1}(\tilde{k})+\phi_{6}I_{m}(\tilde{k})\right],$ $\displaystyle
M_{46}=\frac{\mathrm{i}}{2}\left[kK_{m-2}(\tilde{k})-kK_{m+2}(\tilde{k})+\frac{2m\left\\{m^{2}-1-(3\mathcal{A}-k^{2})(\eta+1)^{2}\right)\\}K_{m}(\tilde{k})}{k(1+\eta)^{3}\omega}\right],$
$\displaystyle
M_{51}=-i\tilde{k}(m+1)I_{m-1}(\tilde{k})+\mathrm{i}\left[\tilde{k}^{2}+2m(m+1)\right]I_{m}(\tilde{k}),$
$\displaystyle
M_{52}=-i\tilde{k}(m-1)K_{m-1}(\tilde{k})+\mathrm{i}\left[\tilde{k}^{2}-(m-2)m\right]K_{m}(\tilde{k}),$
$\displaystyle M_{53}=2\tilde{k}kI_{m-1}(\tilde{k})-3mkI_{m}(\tilde{k}),\
M_{54}=-2\tilde{k}kK_{m-1}(\tilde{k})-3mkK_{m}(\tilde{k}),$ $\displaystyle
M_{55}=kmI_{m}(\tilde{k}),\ M_{56}=-kmK_{m}(\tilde{k})$ $\displaystyle
M_{61}=\mathrm{i}(m+1)\left[\tilde{k}\left(\tilde{k}^{2}+2m(m+2)\right)I_{m-2}(\tilde{k})-m\left(3\tilde{k}^{2}+4\left(m^{2}+m-2\right)\right)I_{m-1}(\tilde{k})\right],$
$\displaystyle
M_{62}=\mathrm{i}\left[\tilde{k}\left\\{\tilde{k}^{2}-2m\left(m^{2}-m-2\right)\right\\}K_{m-2}(\tilde{k})-m\left[\tilde{k}^{2}(m-3)+4(m-2)(m-1)(m+1)\right]K_{m-1}(\tilde{k})\right],$
$\displaystyle
M_{63}=4mk\left[\tilde{k}^{2}+2\left(m^{2}-1\right)\right]I_{m-1}(\tilde{k})-k\tilde{k}\left[\tilde{k}^{2}+4m(m+1)\right]I_{m-2}(\tilde{k}),$
$\displaystyle
M_{63}=-k\tilde{k}\left[\tilde{k}^{2}+4m(m+1)\right]K_{m-2}(\tilde{k})-4mk\left[\tilde{k}^{2}+2\left(m^{2}-1\right)\right]K_{m-1}(\tilde{k}),$
$\displaystyle
M_{63}=k\tilde{k}\left[\tilde{k}^{2}+2m(m+1)\right]I_{m-2}(\tilde{k})-2mk\left[\tilde{k}^{2}+2\left(m^{2}-1\right)\right]I_{m-1}(\tilde{k}),$
$\displaystyle
M_{63}=-k\tilde{k}\left[\tilde{k}^{2}+2m(m+1)\right]K_{m-2}(\tilde{k})-2mk\left[\tilde{k}^{2}+2\left(m^{2}-1\right)\right]K_{m-1}(\tilde{k}),$

where $\tilde{k}=k(1+\eta)$ and the definition of the following functions has
been made for conciseness:

$\displaystyle\phi_{1}=\tilde{k}I_{m-1}(\tilde{k}),\ \
\phi_{2}=2\tilde{k}\left(-3\mathcal{A}(\eta+1)^{2}+\tilde{k}^{2}+m^{2}-(\eta+1)(m+4)\omega-1\right),$
$\displaystyle\phi_{3}=6\mathcal{A}(\eta+1)^{2}(m+2)+\tilde{k}^{2}[3(\eta+1)\omega-2(m+2)]-2(m+1)(m+2)\left[m-2(\eta+1)\omega-1\right],$
$\displaystyle\phi_{4}=2\tilde{k}\left[\tilde{k}^{2}-3\mathcal{A}(\eta+1)^{2}+m^{2}+(\eta+1)(3m-4)\omega-1\right],$
$\displaystyle\phi_{5}=-6\mathcal{A}(\eta+1)^{2}(m-2)+\tilde{k}^{2}[3(\eta+1)\omega+2m-4]+2(m-2)(m+1)(-2(\eta+1)\omega+m-1),$
$\displaystyle\phi_{6}=m^{2}-1-(3\mathcal{A}-k^{2})(\eta+1)^{2}-2(\eta+1)(m+1)\omega,$

The dispersion relation is finally obtained upon requiring that
$\operatorname{det}(\mathsf{M})=0$. The necessary algebra needed to arrive at
the final expression for the growth rate
$\omega=\omega(k,\mathcal{A},\eta,m)$, which consists in finding the
determinant of the matrix $\mathsf{M}$ and solving for $\omega$, is carried
out using Mathematica [13] and, due to its length, it is not amenable for
transcription. By means of validation, the eigenproblem (34)–(37), (43)–(47)
is solved separately using a spectral collocation technique. To that end, the
physical domain $\eta<r<\eta+1$ is mapped onto $N$ Chebyshev nodes defined on
the domain $-1<\xi<1$ by means of the transformation $2r=1+2\eta-\xi$. The
differentiation matrices are obtained following Trefethen [30] upon
discretization of the Chebyshev domain, which eventually are used to recast
the eigenproblem into the generalized form
$\left(\mathsf{L}-\omega\mathsf{B}\right)\bm{q}=0$, where $\mathsf{L}$ and
$\mathsf{B}$ are the discretized version of their differential counterparts
and $\bm{q}=({\hat{\bm{u}\mkern 3.0mu}\mkern-3.0mu}{},{\hat{p\mkern
3.0mu}\mkern-3.0mu}{},{\hat{h\mkern 3.0mu}\mkern-3.0mu}{})^{\mathrm{T}}$ is
the vector of eigenfunctions. The problem is solved using the Matlab routine
eigs, which is based on the ARPACK [18] library. Below is a comparison for the
results obtained for the sample case $\mathcal{A}=0.01$, $\eta=10$, $m=1$,
$k=0.1$, where the reference computation by Mathematica is
$\omega=0.0001260252$ and the results provided by the spectral collocation
technique are

$N$=5, | 1.261314041e-04
---|---
$N$=10, | 1.260252407e-04
$N$=20, | 1.260252407e-04
$N$=50, | 1.260252379e-04
$N$=100, | 1.260252164e-04

showing that even for $N=5$ Chebyshev nodes, the relative error falls below
$10^{-3}$.

## References

  * Blossey [2012] R. Blossey.  _Thin liquid films: dewetting and polymer flow_.  Springer Science & Business Media, 2012.
  * Carroll and Lucassen [1974] B.J. Carroll and J. Lucassen.  Effect of surface dynamics on the process of droplet formation from supported and free liquid cylinders.  _J. Chem. Soc. Faraday Trans._ , 70:1228–1239, 1974.
  * Craster and Matar [2009] R.V. Craster and O.K. Matar.  Dynamics and stability of thin liquid films.  _Rev. Mod. Phys._ , 81(3):1131, 2009.
  * De Gennes [1985] P.-G. De Gennes.  Wetting: statics and dynamics.  _Rev. Mod. Phys._ , 57(3):827, 1985.
  * Dietze and Ruyer-Quil [2015] G.F. Dietze and C. Ruyer-Quil.  Films in narrow tubes.  _J. Fluid Mech._ , 762:68–109, 2015.
  * Dumbleton and Hermans [1970] J.H. Dumbleton and J.J. Hermans.  Capillary instability of a thin annular layer of liquid around a solid cylinder.  _Industrial & Engineering Chemistry Fundamentals_, 9(3):466–469, 1970.
  * Duprat et al. [2007] C. Duprat, C. Ruyer-Quil, S. Kalliadasis, and F. Giorgiutti-Dauphiné.  Absolute and convective instabilities of a viscous film flowing down a vertical fiber.  _Phys. Rev. Lett._ , 98(24):244502, 2007.
  * Frenkel et al. [1987] A.L. Frenkel, A.J. Babchin, B.G. Levich, T. Shlang, and G.I. Sivashinsky.  Annular flows can keep unstable films from breakup: nonlinear saturation of capillary instability.  _J. Colloid Interface Sci._ , 115(1):225–233, 1987.
  * Goren [1962] S.L. Goren.  The instability of an annular thread of fluid.  _J. Fluid Mech._ , 12(2):309–319, 1962.
  * Goren [1964] S.L. Goren.  The shape of a thread of liquid undergoing break-up.  _J. Colloid Sci._ , 19(1):81–86, 1964.
  * Grotberg [1994] J. B. Grotberg.  Pulmonary flow and transport phenomena.  _Annu. Rev. Fluid Mech._ , 26(1):529–571, 1994\.
  * Hammond [1983] P.S. Hammond.  Nonlinear adjustment of a thin annular film of viscous fluid surrounding a thread of another within a circular cylindrical pipe.  _J. Fluid Mech._ , 137:363–384, 1983.
  * Inc. [2021] Wolfram Research, Inc.  Mathematica, Version 12.3.1, 2021.  URL https://www.wolfram.com/mathematica.  Champaign, IL.
  * Israelachvili [1985] J.N. Israelachvili.  _Intermolecular and surface forces_.  Academic Press, 1985.
  * Johnson et al. [1991] M. Johnson, R.D. Kamm, L.W. Ho, A. Shapiro, and T.J. Pedley.  The nonlinear growth of surface-tension-driven instabilities of a thin annular film.  _J. Fluid Mech._ , 233:141–156, 1991.
  * Kalliadasis and Chang [1994] S. Kalliadasis and H.-C. Chang.  Drop formation during coating of vertical fibres.  _J. Fluid Mech._ , 261:135–168, 1994.
  * Kalliadasis et al. [2011] S. Kalliadasis, C. Ruyer-Quil, B. Scheid, and M.G. Velarde.  _Falling liquid films_ , volume 176.  Springer Science & Business Media, 2011.
  * Lehoucq et al. [1998] R.B. Lehoucq, D.C. Sorensen, and C. Yang.  _ARPACK users’ guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods_.  SIAM, 1998.
  * Lin et al. [2002] C.K. Lin, C.C. Hwang, and T.C. Ke.  Three-dimensional nonlinear rupture theory of thin liquid films on a cylinder.  _J. Colloid Interface Sci._ , 256(2):480–482, 2002.
  * Lister et al. [2006] J.R. Lister, J.M. Rallison, A.A. King, L.J. Cummings, and O.E. Jensen.  Capillary drainage of an annular film: the dynamics of collars and lobes.  _J. Fluid Mech._ , 552:311–343, 2006.
  * Moctezuma-Sánchez and Dávalos-Orozco [2015] M. Moctezuma-Sánchez and L.A. Dávalos-Orozco.  Azimuthal instability modes in a viscoelastic liquid layer flowing down a heated cylinder.  _Int. J. Heat Mass Transfer_ , 90:15–25, 2015.
  * Note [1] Note1.  For a fibre coated by a $1000\leavevmode\nobreak\ \mathrm{cSt}$ silicone oil film of initial thickness $h_{o}=100\leavevmode\nobreak\ \mathrm{nm}$, $\rho=970\leavevmode\nobreak\ \mathrm{kg}/\mathrm{m}^{3}$, $\mu=0.97\leavevmode\nobreak\ \mathrm{Pa}\cdot\mathrm{s}$, $\sigma=27\times 10^{-3}\mathrm{N}/m$ we get $v_{g}/v_{\mathrm{vdW}}=1.7\times 10^{-4}$ and $La=2.8\times 10^{-6}$.
  * Note [2] Note2.  A fully detailed asymptotic expansion of the momentum equation and the corresponding boundary conditions is available in the work by Zhao et al. [35] .
  * Oron et al. [1997] A. Oron, S. H. Davis, and S. G. Bankoff.  Long-scale evolution of thin liquid films.  _Rev. Mod. Phys._ , 69(3):931, 1997.
  * Quéré [1999] D. Quéré.  Fluid coating on a fiber.  _Annu. Rev. Fluid Mech._ , 31(1):347–384, 1999\.
  * Quéré et al. [1990] D. Quéré, J.-M. Di Meglio, and F. Brochard-Wyart.  Spreading of liquids on highly curved surfaces.  _Science_ , 249(4974):1256–1260, 1990.
  * Roe [1975] R.-J. Roe.  Wetting of fine wires and fibers by a liquid film.  _J. Colloid Interface Sci._ , 50(1):70–79, 1975\.
  * Ruyer-Quil et al. [2008] C. Ruyer-Quil, P. Treveleyan, F. Giorgiutti-Dauphiné, C. Duprat, and S. Kalliadasis.  Modelling film flows down a fibre.  _J. Fluid Mech._ , 603:431–462, 2008.
  * Scheludko [1962] A. Scheludko.  Sur certaines particularités des lames mousseuses. ii. stabilité cinétique, épaisseur critique et épaisseur d’équilibre.  _Proc Koninkl Nederl Akad Wetenschap_ , 65:87, 1962.
  * Trefethen [2000] L.N. Trefethen.  _Spectral methods in MATLAB_.  SIAM, 2000.
  * Vrij [1966] A. Vrij.  Possible mechanism for the spontaneous rupture of thin, free liquid films.  _Discussions of the Faraday Society_ , 42:23–33, 1966.
  * Vrij and Overbeek [1968] A. Vrij and J.T.G. Overbeek.  Rupture of thin liquid films due to spontaneous fluctuations in thickness.  _J. Am. Chem. Soc._ , 90(12):3074–3078, 1968\.
  * Wray et al. [2013] A.W. Wray, D.T. Papageorgiou, and O.K. Matar.  Electrostatically controlled large-amplitude, non-axisymmetric waves in thin film flows down a cylinder.  _J. Fluid Mech._ , 736:R2, 2013.
  * Yarin et al. [1993] A.L. Yarin, A. Oron, and P. Rosenau.  Capillary instability of thin liquid film on a cylinder.  _Phys. Fluids A_ , 5(1):91–98, 1993.
  * Zhao et al. [2023] C. Zhao, Y. Zhang, and T. Si.  Slip-enhanced Rayleigh–Plateau instability of a liquid film on a fibre.  _J. Fluid Mech._ , 954:A46, 2023.