\section{Introduction}

Flash floods are among the most destructive surface processes around the world, especially in mountainous areas \citep{au1998,borga2011,gomez2005a,jonkman2005}. They are mainly initiated by rapid and intense rainfall, often discharged in few hours \citep[e.g.,][]{borga2007,Bout2018,he2018,loczy2012}, and by complex interactions of the climatic conditions with topography and hydrology \citep[e.g.,][]{hatheway2005}. Because of the very rapid raise in water levels caused by flash floods, it is challenging to take timely and effective actions to contain the associated damage. Flash flood disasters are essentially flash floods that have caused losses either in terms of human lives or economy, or both \citep{gaume2009,jonkman2005analysis,kelman2004}. In China, approximately 70\% of the total area is covered by mountains and hills, which exposes a substantial surface of the national territory to flash flood disasters' risk \citep{liu2018}. Additionally, the more frequent extreme precipitation associated with climate change has increased the number of flash flood disasters in recent decades \citep{sampson2015}.

Historical inventories of flash flood disasters are a precious source of information allowing to investigate their spatiotemporal pattern distribution and evolution. Furthermore, this information can be related with the geomorphological setting of the area and the climatic/meteorological conditions to detect triggering factor, highlight the more vulnerable areas, and to prevent and forecast their effects in the future.

The susceptibility to hydro-geomorphological processes is commonly assessed by considering only the spatial distribution of observed events \citep{cama2015predicting,cama2017improving,santangelo2012,zaharia2017}. However, this is purely a convenient assumption from the modeling perspective. Recently, a growing amount of evidence indicates that these events tend to aggregate in space conditioned by the temporal variability, attesting for an interaction between space and time on event frequency and distribution  \citep{GARIANO2016227,kouli_landslide_2010,nhess-15-2127-2015,merz2016, tonini2019}. In other words, when an event occurs at a specific location, a temporary increase in the probability that other events will cluster at nearby locations should be accounted for. This increase in probability can be captured through clustering analyses and various examples already exist in literature where this has been done at different spatial and temporal scales and via different analytical approaches. Notably, this type of application spans in many areas of natural hazards and have become mainstream in case of seismicity \citep[e.g.,][]{fischer2003,GEORGOULAS20134183, varga2012,woodward2018, yang_spatial-temporal_2019}, joint sets and their orientation in rock outcrops \citep[e.g.,][]{tokhmechi2011,zhan2017comprehensive}, groundwater monitoring \citep{chambers2015}, wildfires \citep[e.g.,][]{orozco2012, costafreda2016,fuentes2013,tonini2017}, and landslides \citep[e.g.,][]{Lombardo.etal:2018,lombardo2019geostatistical,tonini2019}. In the specific case of flooding, \citet{zhao2014} used the projection pursuit theory to cluster spatial data and to build a dynamic risk assessment model for flood disasters. Moreover, \citet{renard2017} detected flood vulnerability accounting for clustering effects in key areas with high flood risk. \citet{pappada2018} also investigated the flood risks in a given region and identified clusters where the floods show a similar behavior with respect to multivariate criteria. Another example can be found in \citet{merz2016} where the authors analyzed the inter-annual and intra-annual flood clustering in Germany. All these examples confirm a substantial scientific interest in recent years dedicated to investigate the clustering behaviors of flash floods and the associated risk; and, more generally, to concurrently analyze their spatial and temporal persistence. However, despite the scientific efforts, detecting flash flood patterns at long temporal scale is still scarce in literature, mainly because of technical limitations. In fact, limited information and records are available in digital form reporting locations and dates of flash flood (and flash flood disasters), especially over long periods. Nevertheless, very recent advances in data collection and sharing techniques are gradually filling this gap, and an increasing number of databases are being published and made available to the scientific community with the records of historical and hydro-geomorphological disasters at the global, continental, or regional scale over long periods \citep{archer2019,de2019,gourley2013,haigh2017,nowicki2020,vennari2016,wood2020}.

Typically, flash flood disasters (as many other hydro-geomorphological disasters) can be considered as a stochastic point processes \citep{stoyan2006} acting in both spatial and temporal dimensions \citep[e.g.,][]{lombardo2019space}. Point patterns can be analyzed in terms of their random distribution, dispersion and clustering behaviour \citep{merz2016,tonini2019}. Several methods can be implemented to deal with stochastic properties. Some classic models, such as Moran’s I \citep{moran1950}, Ripley’s K-function \citep{ripley1977}, fractal dimension \citep{lovejoy1986}, and Allan factor \citep{allan1966}, have been used to detect clustering behaviour in space and in time. Representative models for local clustering analysis (i.e. allowing to detect clusters and their specific location) include Geographical Analysis Machine \citep[GAM,][]{openshaw1987}, Turnbull’s Cluster Evaluation Permutation Procedure \citep[CEPP,][]{turnbull1990}, Scan Statistics \citep{kulldorff1997}, and DBSCAN \citep{ester1996}. For flash floods, which are triggered by storms, the temporal dependency among persistent events is mainly driven by climatic and meteorological conditions. However, global cluster indicators only take into consideration one dimension, disregarding the interaction between space and time. In this sense, spatiotemporal Scan Statistics is a good tool to detect clusters since it allows to identify statistically significant excess of observations thanks to a moving cylindrical window that scans all locations both in space and time \citep{kulldorff_evaluating_1998}. Therefore, it is especially useful to investigate hydro-geomorphological processes such as flash floods. For such phenomena, the detection of events aggregated over a given region and in a specific period, generally yields more informative results than the purely spatial or temporal analyses. Furthermore, understanding the magnitude of the persistence for flash floods is an important requirement to predict where, when and how their highest probability to occur distributes in the future.

In this study, we explored the spatiotemporal pattern distribution of flash flood disasters which have caused either or both life and economic losses in China over the period 1950-2015. Firstly, the deviation of flash flood disasters from a spatiotemporal random process is explored by applying the spatiotemporal Ripley's K-function. Then, the Scan Statistics was applied to detect statistically significant spatiotemporal clusters. Finally, the possible relationship between the detected clusters and local climatic proxy factors is discussed. To the best of our knowledge, it is the first time that such a long-term inventory is analysed to explore the spatiotemporal patterns of flash flood disasters, especially in China. This study provide useful insights on flood dynamics over a large spatiotemporal domain. Moreover, because of the long time-span, results can be useful to indicate how flash flood disasters have evolved in response to climate changes.

\section{Material and methods}
\subsection{Data description}
\subsubsection{Study area}
China lies between latitudes 18$^\circ$ and 54$^\circ$ N, and longitudes 73$^\circ$ and 135$^\circ$ E. With an area of about 9.6 million square kilometers, it is the world's third-largest country. The landscape varies significantly across this vast area, ranging from the Gobi and Taklamakan deserts in the north to the subtropical forests in the wetter south. The eastern plains and southern coasts are the location of most of China's agricultural land and settlements. The southern areas consist of hilly and mountainous terrain. The west and north of the country are dominated by sunken basins (such as the Gobi and the Taklamakan desert), towering massifs and rolling plateaus, including part of the highest tableland on earth, the Tibetan Plateau. Based on its topography, China can be divided into six homogeneous geomorphological macro-regions \citep{wang2020geomorphological}: eastern plain, southeastern hills, southwestern mountains, north-central plains, northwestern basins and Tibetan Plateau. Mountains (33\% of the territory), plateaus (26\%) and hills (10\%) account together for nearly 70\% of the entire surface.

Influenced by the East Asian summer monsoon and the geomorfologic settings, the climatic condition across the whole country varies considerably. In general, the wet season in China lasts from May to September \citep{song2011rain}. In the Eastern area, the annual rainfall decreases from south to north with an average annual precipitation that ranges from 250 to 750 $mm$ \citep{zhang2007climate}. In the west and central part of North China, due to its far distance away from ocean, the climate tends to be more arid and the landscape transitions to large deserts. The Tibetan plateau is characterized by wet and humid summers with cool and dry winters. More than 60–90\% of the annual total precipitation falls between June and September \citep{xu2008decadal}.

\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE1.jpg}
	\caption{Distribution of flash flood disasters and background setting of China.}
	\label{Figure1}
\end{figure}

\subsubsection{Flash flood disasters inventory}
The dataset used in this study has been collated and made accessible for the present research as part of a national effort carried out by the Chinese Institute of Water Resources and Hydropower Research \citep{liu2018}. It reports flash flood occurrences in China since 1950 until 2015 together with available information, namely longitude and latitude, date, fatalities and economic losses. Due to the lack of specific terminology and/or detailed descriptions of the disaster process in the database, the data does not differentiate the initial mechanism, be it water floods or debris floods/flows \citep[e.g.,][]{fernandez2010,gartner2014}. The only common information is that for each specific case, a large amount of overland flows, mixed with an unspecified solid fraction, rapidly flooded a given area with disastrous effects \citep[e.g.,][]{chang2011,pierson1987}.

To better understand the spatiotemporal dynamics of flash flood and associated disasters, as well as the relationship with the triggering factors, the date of occurrence is of vital importance. Therefore, for consistency reasons, we considered only the records whose metadata contained a full temporal description (year-month-day) resulting in a subset of 32,473 flash flood disasters (accounting for 68\% of the entire dataset) precisely located in space and time (Figure \ref{Figure1}). We further defined the impact of flash flood disasters as the combination of fatalities and economic losses (see Table \ref{Table1}), and we refer to this classification throughout the manuscript.

\begin{table}[H]
	\centering
	\caption{Impact of flash flood disasters (RMB = renminbi, the official currency of China).}
	\includegraphics[width=0.7\linewidth]{Table1.pdf}
	\label{Table1}
\end{table}

\subsection{Methodological overview}\label{methods}

\subsubsection{Spatiotemporal K-function}

The Ripley’s K-function ($K_{(s)}$) is largely applied in environmental studies to analyse the pattern distribution of spatial point processes and to detect deviation from spatial randomness. $K_{(s)}$ allows to determine if a set of mapped punctual events show a random, dispersed or cluster distribution pattern over increasing distance values \citep{ripley1977}. It is computed as the ratio between the expected number of  events falling at a distance $r$ from an arbitrary event and the average number of points per unit area, corresponding to the intensity of the spatial point process ($\lambda$). In the same way, it is possible to define the temporal K-function ($K_{(t)}$) allowing to asses for the randomness of events in time. The spatiotemporal K-function ($K_{(s,t)}$) is a generalization of the univariate Repley’s K-function which allows to test for the independence between two variables, space ($s$) and time ($t$). Therefore, the $K_{(s,t)}$ is a suitable tool to investigate the clustering behaviour of a set of events occurred in a given area at a given time.  For a point process $X$ with intensity $\lambda$, according to equation \ref{eq:EQ1}, it is defined as the number of expected further events ($E$) occurring within a distance $r$ and time $t$ from an arbitrary event $u$, where $a$ define the contouring circle.

\begin{equation}\label{eq:EQ1}
K_{(s,t)} = 1/\lambda \times E[n(X \cap a(u,r,t){u})| u \in X]
\end{equation}

To illustrate the interaction between space and time, it can be useful to evaluate the value $D_{(s,t)}$, defining the difference between the spatiotemporal K-function and the product of the purely spatial and the purely temporal K-function (see equation \ref{eq:EQ2}).
\begin{equation}\label{eq:EQ2}
D_{(s,t)} = K_{(s,t)} - K_{(s)} \times K_{(t)}
\end{equation}

If space and time are independent variables, this value equals to zero. Otherwise, positive values of $D_{(s,t)}$ indicates the interaction among events in space and in time. In other words, events closer in space are more likely to occur in a closer time. On the contrary, the negative values means a dispersed pattern.

In this study, spatiotemporal K-function analyses were performed with the package ``Spatial and Space-Time Point Pattern Analysis'' \citep[splancs,][]{rowlingson2017} in R \citep{team2019r}.

\subsubsection{Spatiotemporal scan statistics}

Scan statistic was originally developed by \citet{naus_clustering_1965,naus_distribution_1965} to detect cluster in a one-dimensional point process. Subsequently \citet{kulldorff1997} extended this approach to multi-dimensional point process, introducing the use of scanning windows. The procedure was implemented into a free software, SaTScan\textsuperscript{TM} (\href{satscan.org}{satscan.org}) which can handle a purely spatial, purely temporal or spatiotemporal datasets and includes different probability models depending on the nature of the data and the scope of the research (e.g. for prospective or retrospective cluster detection). In the purely spatial case, the aim of scan statistics is the early detection of clusters, allowing one to map them and to assess their statistical significance. Moving windows scan the region increasing their radius up to a fixed limit ($R_{max}$) and count the number of events falling inside and outside the area. The probability that a window contains more observations than expected is assessed via the likelihood ratio, by comparing with the background population. Then, the null hypothesis of randomness is tested by Monte Carlo experiments, based on repeated random sampling. The spatiotemporal scan statistic use cylinders instead of circular windows, where the height of the cylinder account for the temporal dimension.
In order to deal with flash foods, the retrospective spatiotemporal permutation scan statistics \citep[STPSS,][]{kulldorff_spacetime_2005} seems to be the most adequate model. Indeed, for environmental processes, the definition of the background population at risk needed for the statistical significance assessment of the detected clusters is quite problematic. STPSS assesses the expected number of cases using only the observed cases by permutation, supposing that each event has the same probability for all the times. Computationally, if $C$ is the total number of observer cases and $c_{zd}$ the number of cases observed in a zone $z$ and a day $d$, the expected number of cases per zone and day ($\mu_{zd}$) is:

\begin{equation}\label{eq:mu_zd}
\mu_{zd} = \frac1C\left(\sum_z c_{zd}\right) \left(\sum_d c_{zd}\right)
\end{equation}

It follows that, for a spatiotemporal cylinder $A$, the expected number of cases $\mu_A$ can be estimated as the sum of each $\mu_{zd}$ inside the cylinder $A$:

\begin{equation}\label{eq:muA}
\mu_A = \sum_{z\,,d\in A} \mu_{zd}
\end{equation}

If $C_A$ is the number of observed cases in $A$, considered as Poisson-distributed with mean $\mu_A$, the Poisson generalized likelihood ratio ($GLR$) can be computed as:

\begin{equation}\label{eq:GLR}
GLR = \left(\frac{c_{A}}{\mu_{A}}\right)^{c_{A}}\left(\frac{C - c_{A}}{C - \mu_{A}}\right)^{C - c_{A}}
\end{equation}

This ratio is calculated and maximized for every possible scanning cylinder. Then, the Monte Carlo simulations are performed and the statistical significance of the detected clusters is assigned by ranking the clusters according to their $GLR$-value.

\section{Results}

\subsection{Deviation from a random process}

In the present study the spatiotemporal K-function is used to assess the interaction between the two variables, space and time, in generating clusters at increasing  distances. Figure \ref{Figure2} (panel a) shows the 3D-plot of $D_{(s,t)}$ with a zoom up to 2000 $km$ (panel b). Positive values indicate that space and time interact in generating clusters: in other words, events closer in space are also closer in time. This is the case at any increasing distance, from hundred meters to thousands meters and from few years to decades.

In addition, we computed the spatiotemporal K-function separately for the eastern and western side of China (Figure \ref{Figure3}). We did this because the southeastern area, which is the rainiest part of the country, is highly affected by flash floods, while the northwestern area is predominantly desert and flash floods are less frequent. It results that, although events are clustered in both the areas, in the southeastern area (panel a) clusters arise at a shorter spatial distance and closer in time than in the northwestern area (panel b). More specifically, in the southeast China the spatiotemporal interaction generates clusters starting from 200 $km$ and a plateau is reached at about 1800 $km$. In Northwest China the global cluster behaviour is more evident from about 1000 $km$ to higher distances. As regards the temporal dimension, the two part of the country show a similar cluster behaviour, with a strong attraction among events during the first 20 years lasting in time with a more relaxed clustering behaviour.

To summarize, the spatiotemporal K-function reveals a deviation of flash flood disasters and associated spatiotemporal pattern distribution from a random process at specific scales, measured and quantified both in space, as distances-values, and in time, as yearly periods. These values can provide a useful indication to set up the parameters for further clustering algorithms, acting at local scale such as, for example, the spatiotemporal scan statistics.
\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE2.jpg}
	\caption{Three dimensional summary of flash flood disasters in China during 1950-2015.}
	\label{Figure2}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE3.jpg}
	\caption{Three dimensional summary of flash flood disasters in China, separated between two eastern and western sectors and with a maximum spatial bandwidth of 2000 $km$.}
	\label{Figure3}
\end{figure}

\subsection{Spatiotemporal clusters}

\subsubsection{Cluster detection and spatial distribution}\label{Sec3.2.1}

Scan statistics was performed to detect spatiotemporal clusters of flash flood disasters. The size and the duration of the detected clusters are influenced by the input parameters of the scanning windows, namely the maximum radius ($R_{max}$), the maximum temporal duration ($T_{max}$), and the time aggregation ($T_{agg}$). Indeed, values of $R_{max}$ exceeding the 50\% of the total area or, for $T_{max}$, the 50\% of the entire study period, can result in an exceptionally low rate outside the scanning window rather than detecting an exceptionally high rate inside. $T_{agg}$ is used to adjust the aggregation of the data over time and allows adjusting for cyclic temporal trends: for example, a time aggregation of one year automatically adjusts for the seasonal variability, while the contrary happen with monthly aggregations. Moreover, both spatial and temporal aggregations can highly reduce the computer processing time. Following the results obtained by the spatiotemporal K-function and discussed above, few radii for each area (southeast and northwest China) were tested. Performed analyses indicated that the effect onto the detected clusters were negligible and finally we considered the spatiotemporal distribution of flash flood disasters as a whole rather than splitting the Chinese territory in two areas. We opted for a set of possible combinations of $R_{max}$ and $T_{max}$, keeping $T_{agg}$ fixed to one year. More specifically, to compare the combination of these parameters, and to obtain reasonable clusters, we tested three $R_{max}$ values equal to 100, 200 and 300$km$, and three $T_{max}$ values equal to 1, 3 and 5 years. The choice for $R_{max}$ is corroborated by \citet{zhang2010} who report measurements constantly less than 500 $km$ for the radius of typical convective storms in the Chinese mainland, which can trigger flash floods. Results of STPSS for each of the nine combinations of these parameters are shown in Figure \ref{Figure4}.

\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE4.jpg}
	\caption{Significant (p$<$0.005) spatiotemporal clusters of flash flood disasters in China during 1950-2015. }
	\label{Figure4}
\end{figure}

\begin{table}[H]
	\centering
	\caption{Number of detected spatiotemporal clusters of flash flood disasters in China during 1950-2015 using different parameters.}
	\includegraphics[width=0.4\linewidth]{Table2.pdf}
	\label{Table2}
\end{table}

The largest variation in the number of detected clusters is mainly associated with $R_{max}$ -- as $R_{max}$ increases, the detected flash flood disaster clusters exhibit a clear decrease -- rather than with $T_{max}$. This result is summarized in Table \ref{Table2}. More specifically, large $R_{max}$ values affect the detection of clusters acting at a fine scale, which tend to be missed or merged into larger ones; conversely, very large clusters, acting at a coarse spatial scale, are still detected. This is geographically visible in the south-easternmost sector of China (Figure \ref{Figure4}). Changes on $T_{max}$ have almost no effect on the number of clusters since, even allowing for a maximum duration of 5 years, almost all the clusters do not exceed one year.
As complementary information, Table \ref{Table3} presents the temporal duration, expresses as start and end date, for the first ten clusters of flash flood disasters using $T_{max}$ equals to 3 years and for increasing values of $R_{max}$, equal to 100, 200 and 300 $km$. Results confirm that the cluster duration does not exceed one year. The most significant cluster was detected in 1975, while the rating for the following clusters can change in the three cases. Nevertheless, it is important to notice that the top-ten clusters are well distributed over the entire study period, with the oldest one detected between 1963 and 1969.

\begin{table}[!htbp]
	\centering
	\caption{Temporal duration of the first 10 clusters of flash flood disasters detected via three different models (left: $R_{max}=100 km$; center: $R_{max} = 200 km$; right: $R_{max} = 300 km$)}
	\includegraphics[width=1\linewidth]{Table3.pdf}
	\label{Table3}
\end{table}

\subsubsection{Clusters characterization}

Detected clusters where further analyzed by considering the impact of flash flood disasters. To this end, we examined only clusters detected by using $R_{max}=200 km$ and $T_{max}=3 years$.  The choice of a $T_{agg}=1 year$ was originally meant to focus our analyses on effects that may exhibit a yearly cycle. However, this would have smoothed nested effects acting at the seasonal scale. For this reason, we opted to carry out additional analyses using a $T_{agg}$ of three months (hereafter referred as \emph{monthly model}). Results are shown in Figure \ref{Figure5} where information on the spatial distribution of the detected clusters is combined with the impact related to the single flash flood events (see Table \ref{Table1}). Overall, the clusters chiefly appear along the main river systems in China, namely the Yangtze, the Yellow, the Pearl and the Yarlung Zangbo Rivers. In addition, some clusters stand out on high mountains such as the Qinling-Daba and the Changbai Mountains.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=\linewidth]{FIGURE5.jpg}
	\caption{Significant (p$<$0.005) spatiotemporal clusters of flash flood disasters in China during 1950-2015 ($R_{max}=200 km$, $T_{max}=3 years$, $T_{agg}=3 months$). Each event belonging to a single cluster is further resized as a function of its impact, in accordance to Table \ref{Table1}.}
	\label{Figure5}
\end{figure}

\begin{table}[!htbp]
	\centering
	\caption{Temporal duration of the first 10 clusters of flash flood disasters during 1950-2015 ($R_{max}=200km$, $T_{max}=1 year$, $T_{agg} = 3 months$).}
	\includegraphics[width=0.5\linewidth]{Table4.pdf}
	\label{Table4}
\end{table}

Forcing the model parameterization to aggregate the time over a fraction of the year (three months) allows us to investigate potential seasonal effects. Indeed, even if the maximum temporal duration is still of one year, looking at the ten most significant clusters detected under the \emph{monthly model} (Table \ref{Table4}), it results that all of them have a duration of three (six clusters) or six (four clusters) months. Notably, almost every cluster (nine clusters) encompass the period from July to September, with an earlier start date (in April) for the ones which have a longer duration.

\subsubsection{Temporal duration of detected clusters}

The temporal variation in the duration of the detected clusters could have been driven by the precipitation regime. In additional, space-time dependency may have been induced by the geomorphological setting of the area and by anthropogenic pressures, but these last factors should have a minor effect compared to the rainfall pattern, which acts as the primary triggering factor of flash floods. Therefore, in the present study we assume the precipitation as the main driver for flash floods detected clusters, and results are interpreted and discussed on the basis of this hypothesis.
Allowing for $T_{max}=3 years$ in the parameterization of the \emph{yearly models}, the temporal duration of the detected clusters ranges from one to three years (see Figure \ref{Figure6}). The cluster detection pattern appears quite clear and well defined. However, since 1980 some clusters partially overlap. This can lead to two separate interpretations. Firstly, the relative small number of clusters detected between 1950 and 1980 may imply that the data acquisition and report in the Chinese database of hydro-morphological disasters was not fully operational at the time. Conversely, from 1980 to present days the Chinese database has evolved into a mature and detailed geographic information system. Secondly, the same pattern can be justified as a result of climatic changes. In fact, overlapping clusters of one, two and three years duration essentially appear only after 1980. These concurrent clusters may reflect similar synchronous variations of the climate settings and rainfall regimes across China in the recent period.

\begin{figure}[H]
	\centering
	\includegraphics[width=\linewidth]{FIGURE6.jpg}
	\caption{Temporal duration of flash flood disasters clusters in China during 1950-2015 ($R_{max}=200km$, $T_{agg}=1 year$, $T_{max}=3 years$).}
	\label{Figure6}
\end{figure}

We summarized the same results for the \emph{monthly model} in Figure \ref{Figure7}. To better visualize the seasonality trend, we opted for a cyclic representation of the detected clusters, plotting their pattern in four temporal duration classes of 3, 6 and 9 months as well as one year. Most clusters show a 3-months duration, concentrated in the period between July and October, and an increasing density after 1980. Furthermore, clusters of 6-months temporal duration are most likely to occur from January to July or from April to October. As for clusters with 9-months temporal duration, these mostly cover the period of July-August-September, irrespective of the starting month. Ultimately, as noticed for the \emph{yearly model}, also in the \emph{monthly model} much more clusters were detected in the late period, mainly from 2000. Moreover, the vast majority of flash flood disasters clusters happened between July and October, a period coinciding with the wet season in China.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.6\linewidth]{FIGURE7.jpg}
	\caption{Seasonal effect of flash flood disasters clusters in China during 1950-2015 ($R_{max} = 200km$, $T_{max} = 1 year$, $T_{agg} = 3 months$).}
	\label{Figure7}
\end{figure}

\subsubsection{Recurrence of clusters at decades-scale}

The analyses run in the previous sections were all voted to search for clusters in a relatively small temporal window. However, environmental changes, and especially those related to climate change, usually act on a longer time span. To better investigate this effect, we considered a temporal subdivision of the dataset into six subsets, each one lasting ten years (starting from 1956). For each decade (1956-1965, 1966-1975, 1976-1985, 1986-1995, 1996-2005, 2006-2015) the following parameter for the scanning widows were imposed: $R_{max} = 200km$, $T_{max}= 2 years$ and $T_{agg} = 1 year$. As shown in Figure \ref{Figure8}, the number of detected clusters increases from the early to recent periods. These are compared with the rainfall distribution, derived from the daily rainfall data provided by the China Meteorological Administration (\href{http://data.cma.cn/}{http://data.cma.cn/}). In the present study, only the weather stations (a total of 699 rain gauges) with complete data for the period 1955-2015 were considered. The mean monthly and annual rainfall were computed for each station and this data were then regionalized on a 2$km$ $\times$ 2$km$ lattice, via Ordinary Kriging interpolator. It results that flash floods detected clusters are mainly located in the southeastern most humid regions in every period. However, in the last two decades, clusters appear also in the northwestern arid regions. Even if the rainfall distribution, averaged over each decades, does not allow to discover clear changes along the subsequent periods, these newly detected clusters can be due to the intensification of the extreme rainfall events occurring in the area in recent periods. This assumption is confirmed by the statistics on clusters duration (Figure \ref{Figure9}). From the boxplot summarizing the descriptive statistics it is evident that the median values of clusters duration tends to slightly decrease from 46 days (1956-1965), to 17 days (1986-1995), to stabilise at a value around 20 days in the two last decades. At the same time, the overall duration, measured as difference between the maximum and the minimum value, is higher in the late periods (140 days in 1956-1965, and 93 and 74 days respectively in the two following decades) than in the early periods (about 65 days for the last two decades). This is even more evident looking at the interquantile ranges, which decrease with time. To resume, from these analyses, the number of detected clusters globally increase in time, but their duration drastically decreases in the recent period, indicating a possible activation induced by short-duration extreme rainfall events.

\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE8.jpg}
	\caption{Significant (p $<$ 0.005) spatiotemporal clusters of flash flood disasters in China every ten years. The size of the circles indicates the impact of flash flood disasters according to the classification proposed in Table \ref{Table1}.}
	\label{Figure8}
\end{figure}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=0.6\linewidth]{FIGURE9.jpg}
	\caption{Boxplots summarizing the descriptive statistics of the duration of clusters reported on Figure \ref{Figure8}.}
	\label{Figure9}
\end{figure}

 We further explored how many times the clusters detected from the previous investigation overlap, considering the catchment level. Results provide an useful information on the recurrence of clusters of flash flood disasters every ten years. To perform this analysis, the centroid of each cluster (with reference to Figure \ref{Figure8}) was extracted and intersected with the catchment boundaries. In a second step, the number of repeated clusters per catchment was computed and  their distribution investigated through the time. Results are shown in Figure \ref{Figure10}, where panel (a) reports the number of repeated clusters and panel (b) reports the information on their relative occurrence across time (similarly to the concept of return time but in the context of spatiotemporal clustering).

 Figure \ref{Figure10}a shows that, as for the spatial trends of the detected clusters, the catchments with recurrent clusters are mainly located in the southeast sector and essentially in the coastal mountains. From Figure \ref{Figure10}b it emerges that, on average, most of the repeated cluster occur with an interval between 10-20 and 20-50 years.
\begin{figure}[t!]
	\centering
	\includegraphics[width=\linewidth]{FIGURE10.jpg}
	\caption{Catchments with clusters detected more than once (a) and return period for the clusters (b).}
	\label{Figure10}
\end{figure}

\section{Discussions}\label{Discuss}

The present study aims at exploring the spatiotemporal clustering characteristics of flash flood disasters in China. For this purpose, we analyzed the official historical inventory, which covers a very long period (from 1950 to 2015). Results are interpreted with a particular regard to the rainfall distribution, being these two processes highly related \citep{wei_rainfall_2018}.
The spatiotemporal K-function was fist computed to assess the deviation of flash flood pattern distribution from a random process. This revealed a clustering behavior at specific spatial distances and yearly periods. Scan Statistics, the spatiotemporal permutation model we adopted, was then performed to identify statistically significant clusters together with their duration (start and end date). This allowed us to detect areas and periods more susceptible to flash flood disasters. We opted for a set of possible combinations for the maximum spatial and temporal extension of the scanning windows, while the data were aggregated both at yearly and at seasonal scale. More specifically, we tested three $R_{max}$ values equal to 100, 200 and 300 $km$, and three $T_{max}$ values equal to 1, 3 and 5 years for the \emph{yearly model}, with an aggregation of three months for the \emph{monthly model}. The most significant cluster resulting from the yearly model was detected in 1975, while the rating for the following clusters can change by varying $R_{max}$; nevertheless, it is important to note that the top-ten clusters are well distributed over the entire study period, with the oldest one detected in 1963-1969. Results of the monthly model show that the top-ten detected clusters have a duration of three (six clusters) or six (four clusters) months. Notably, almost every cluster encompasses the period from July to September, a period coinciding with the wet season in China, with an earlier start date (in April) for the ones which have a longer duration. Globally, much more clusters were detected in the late period, mainly from 2000.
Overall, clusters are chiefly located along the main river systems in China (the Yangtze, the Yellow, the Pearl and the Yarlung Zangbo Rivers). In addition, some clusters stand out on high mountains such as the Qinling-Daba and the Changbai Mountains.

Finally, to detect changes acting at a larger temporal scale, dates were grouped each ten years over the last six decades (from 1956 to 2015). As for the previous analyses, detected clusters are mainly located in the southeastern most humid regions in every period. However, in the last two decades, clusters appear also in the northwestern  arid regions. These newly detected clusters can be due to the intensification the extreme rainfall events occurring in the area in recent periods, as a consequence of climate changes \citep{Song_2011}. This important fact is confirmed by checking the descriptive statistics of the duration of clusters: globally, the number of detected clusters increases in time, but the duration drastically decreases in recent periods, indicating a possible activation induced by short-duration extreme rainfall events.
Our analyses reveled that the catchments with recurrent clusters are mainly located in the southeast sector and essentially in the coastal mountains. We can assume these catchments to be exposed at the highest potential risk across the whole Chinese territory also in the short to long term future. Nevertheless, catchments with repeated clusters in a shorter timespan (5 to 10 years) may also pose a relevant threat, especially in the near future.

In the present study spatiotemporal clusters of flash floods were detected chiefly on the basis of two parameters ($R_{max}$ and $T_{max}$), without featuring terrain attributes, precipitation regimes and anthropogenic pressure. However, these factors may have played and still play a significant role to explain the distribution of flash flood disasters. For instance, the approach we adopted may over-rely on spatial distances to detect clusters. In fact, the natural landscape has mountain belts that can act as orographic barriers to the incoming cloudbursts, effectively limiting the rainfall distribution -- hence flash flood occurrences -- on one or the other side of a catchment divide (at various scales). As for the temporal scale, due to the large timespan, the detected temporal patterns may reflect more information due to long-term climatic variations rather than specific conditions. For this reason, we are planning to extend our spatiotemporal cluster analyses to more complex models, which can concurrently capture multivariate contributions featuring environmental effects, even at the latent level \citep{Lombardo.etal:2018,lombardo2019geostatistical}.

\section{Conclusion}

In this work, we explore the national archive of flash flood disasters in China from 1950 to 2015. The term disaster is meant to describe the destructiveness of the flash floods, since each record in this archive has produced economic, life losses, or both.

The clustering procedure highlighted distinct spatial and temporal patterns at different scales. For instance, flash flood disasters cluster in specific regions and closely follow the mean rainfall distribution. Additionally, we were also able to distinguish seasonal, yearly and even long-term flash flood persisting behaviors. The persistence of disasters is a crucial information because it indicates the risk that a community may undergo in response to a flash flood. Moreover, we studied the cycle of such disasters with particular emphasis on their repeated occurrence per catchment. This complementary information can be further used in relation to engineering and structural design. In fact, infrastructure is usually built to sustain the damage of an event of certain return time. In our analyses, we show that the very same area has been hit and incurred losses up to six times in the last 66 years. This may suggest locally-tailored structural improvements which may lengthen the life expectancy of specific infrastructure as well as reduce the number of victims.

We would like to stress that, as advanced as it may be, our clustering framework is essentially a descriptive tool. And yet, the amount of information one can draw from a descriptive tool can be extremely valuable. Nowadays, the hazard community's effort is mainly dedicated to predictive modeling of various natures and purposes, thus leaving under-explored or even unexplored some basic concepts and interpretative conclusions that data description and visualization can provide. Long time series of national hazard phenomena are one of these examples where studying variations over space and time can highlight very important environmental dynamics, even in the direction of climate change and its implications.

\section*{Acknowledgement}

 This work was supported by the China National Flash Flood Disasters Prevention and Control Project. The authors are grateful for financial support from the China Institute of Water Resources and Hydropower Research (IWHR), grant number No. SHZH-IWHR-57 and National Natural Science Foundation of China, grant number No. 41571388.

\newpage
\bibliographystyle{CUP}


\section{Introduction}

The epoch of reionization (EoR) leaves footprints in the observed cosmic microwave background (CMB) as the photons Thomson scatter off free electrons.
These include damping the primary temperature anisotropies,
inducing secondary anisotropies from the bulk motion of ionized gas
(i.e. the kinetic Sunyaev–Zel'dovich effect),
and prompting curl-less (i.e. E-mode) polarization from the CMB quadrupole (e.g. \citealt{Sunyaev1980MNRAS.190..413S,Vishniac1987ApJ...322..597V,Hu2000ApJ...529...12H,Hu2002ARA&A..40..171H,McQuinn05,Dvorkin2009PhRvD..79j7302D}).
Of these, the large-scale E-mode polarization anisotropies are a particularly powerful probe of the EoR as they are less plagued by degeneracies  and systematics (e.g. \citealt{Reichardt_2016}).
Reionization models can therefore constrain their largely uncertain parameters that describe the ionizing emissivity of the early Universe, through forward-modelling the EE autocorrelation power spectra (PS) and comparing against measurements from the {\it Planck} satellite (\citealt{Planck2016A&A...596A.108P,Planck2019arXiv190712875P}; e.g. \citealt{Hu2003PhRvD..68b3001H, Mortonson2008ApJ...672..737M, Miranda2017MNRAS.467.4050M, Hazra19}).

Nevertheless, most reionization models do not directly use the CMB PS to constrain their parameters.
Instead, they use a summary statistic which has become one of the de-facto standard cosmological parameters -- the direction-averaged Thomson scattering optical depth, $\tau_e$.
In going from the observed PS to $\tau_e$, one needs to adopt a basis set for the EoR history -- a parametrization of the redshift evolution of the comoving number density of free electrons ($n_{e}$).
Early works adopted a simple step function reionization at a given redshift $z_{\rm re}$ \citep{Page2007ApJS..170..335P,Spergel2007ApJS..170..377S}.  Currently, the most common choice has a hyperbolic tangent functional form ({\it Tanh}; \citealt{Lewis_2008}) parametrized by a reionization midpoint ($z_{\rm re}$)  and a redshift duration ($\Delta_{\rm re}$; see equation \ref{eq:tanh}).
For example, the latest constraints on $\tau_e$
published by {\it Planck}, $\tau_e = 0.0522\pm0.0080$ (the {\it TT+lowE} reconstruction in \citealt{Planck2018arXiv180706209P}), were generated by fixing a width of $\Delta_{\rm re}=0.5$,  sampling a flat prior over $z_{\rm re}$, and comparing against the observed large-scale E-mode PS (and also the temperature autocorrelation).

However, the redshift-symmetric evolution given by {\it Tanh}
differs in shape from both physical and empirical models of EoR history, based on the growth of cosmic structure and/or fit to observed galaxy luminosity functions (LFs; e.g. \citealt{Choudhury2006MNRAS.371L..55C,Raskutti2012,Koh2016,Greig2017,Qin2017a,Kulkarni2018,Gorce2018A&A...616A.113G,Roy2018JCAP...05..014R,Qin2020arXiv200304442Q}).
Therefore, {\it computing a likelihood using the $\tau_e$ summary statistic instead of directly forward-modelling the CMB PS can bias EoR model constraints} (e.g. \citealt{Douspis2015A&A...580L...4D,Planck2016A&A...596A.108P,Miranda2017MNRAS.467.4050M, Hazra19}).
For example, \citet{Miranda2017MNRAS.467.4050M} and \citet{Heinrich_2018} claimed (though see \citealt{Millea2018A&A...617A..96M} and \citealt{Planck2018arXiv180706209P}) that the {\it Planck} 2015 E-mode PS prefers reionization histories with an extended tail of partial reionization towards very high redshifts ($z{>}15$).  This would have significant implications for our understanding of the very first galaxies.  However, quantifying such claims is difficult without directly performing inference on galaxy model parameters.

In this work, we use a physically-motivated EoR model to quantify how inference using the $\tau_e$ summary statistic (instead of directly forward-modelling the E-mode PS)
impacts the resulting constraints on galaxy properties and EoR histories.  Using the last {\it Planck} data release, we show that the marginalized constraints on the EoR history are far more sensitive to the choice of the basis set ({\it Tanh} vs physical model) than to the CMB likelihood statistic ($\tau_e$ vs PS).
Specifically, we use the latest v3.0.0 release\footnote{\url{https://github.com/21cmfast/21cmFAST}}
(Murray et al. in prep.) of {{21\textsc{cmFAST}}} \citep{Mesinger2007ApJ...669..663M,Mesinger2011MNRAS.411..955M}, whose parametrization for galaxy properties is informed by high-redshift UV LFs \citep{Park2019MNRAS.484..933P}.
To calculate the CMB PS for a given reionization history and compute the corresponding likelihood,
we add the Cosmic Linear Anisotropy Solving System
(\textsc{class}\footnote{\url{https://github.com/lesgourg/class_public}};
\citealt{Lesgourgues2011arXiv1104.2932L}) Boltzmann solver and the {\it Planck}
likelihood codes (\textsc{plc}\footnote{\url{http://pla.esac.esa.int/pla}}; \citealt{Planck2016A&A...594A..11P,Planck2019arXiv190712875P}) to the upcoming
v1.0.0 release\footnote{\url{https://github.com/21cmfast/21CMMC}} of the public {21\textsc{CMMC}}
Markov chain Monte Carlo (MCMC) framework \citep{Greig2015MNRAS.449.4246G,Greig2017MNRAS.472.2651G}.
All codes developed here are publicly available.

This paper is organized as follows. We present our analysis of the {\it Planck} data and
discuss the difference between the 2015 and 2018 results in Section \ref{sec:model_cmb}.
In Section \ref{sec:model_eor}, we briefly introduce our EoR model (Sec. \ref{subsec:21cmfast})
and show the resulting constraints inferred from CMB and
other observations.  We quantify the bias from the choice of basis set for EoR histories in Sec. \ref{subsec:tanhvs21cmfast} and the choice of likelihood statistics in Sec. \ref{subsec:tauvsps}.  We summarize our results and conclusions in Section \ref{sec:conclusion}.

\section{Modelling CMB observables}\label{sec:model_cmb}

\textsc{class} \citep{Lesgourgues2011arXiv1104.2932L} computes
	CMB anisotropies including temperature and polarization, and calculates their autocorrelation and cross-correlation
	PS. In its default configuration,
	\textsc{class} computes the ionization history from $z{=}10^4$ and
	throughout recombination (e.g. via the code \textsc{recfast};
	\citealt{Seager1999ApJ...523L...1S,Chluba2010MNRAS.402.1195C}), and includes a parametrized function for the EoR history.
	The EoR is assumed to have the {\it Tanh} form as follows
	\begin{equation}\label{eq:tanh}
	n_e {=} \frac{n_{\rm H}{+}n_{\rm He}}{2} \left\{1{+}\tanh\left[\left(1{-}\left(\frac{1{+}z}{1{+}z_{\rm re}}\right)^{1.5}\right)\frac{1{+}z_{\rm re}}{1.5\Delta_{\rm re}}\right]\right\}
	\end{equation}
	with ${n}_{\rm H}$ and ${n}_{\rm He}$ representing the average comoving number density of hydrogen and helium, respectively.
	Instead of using the default {\it Tanh} parametrization,
	here we forward-model the CMB observables by passing
	any given ionization history directly to \textsc{class}.
	We then compare the theoretical PS against observations using the {\it Planck} likelihood codes (\citealt{Planck2016A&A...594A..11P,Planck2019arXiv190712875P}; see more in Table \ref{tab:planck}).
	On the other hand,	the integrated history of reionization
	can also be summarized using the Thomson scattering optical depth
	\begin{equation}\label{eq:tau_e}
	\tau_e = \int_{0}^{z_{\rm d}} c H^{-1}\left(1{+}z\right)^2 n_e\sigma_{\rm T},
	\end{equation}
	where $z_{\rm d}{\sim}1100$, $c$, $H(z)$ and $\sigma_{\rm T}{=}6.65\times10^{-25}{\rm cm^2}$ are the
	redshift of the last scattering surface at recombination,
	the speed of light,
	the Hubble constant at redshift $z$, and
	Thomson scattering cross-section,
	respectively. Below we also
	perform EoR inference with a Gaussian likelihood ($\mathscr{L}$) computed using $\tau_e$ from {\it Planck} \citep{Planck2016A&A...594A..16P,Planck2018arXiv180706209P}. This allows us to compare the resulting EoR parameter constraints to those obtained from using directly the EE PS for the likelihood.

\subsection{Preference for an early reionization?}
The mapping of CMB PS to the $\tau_e$ summary statistic depends on the chosen basis set and corresponding priors for the EoR history.  The sensitivity of the resulting $\tau_e$ constraints to this choice has been debated extensively in the literature (e.g. \citealt{Douspis2015A&A...580L...4D, Planck2016A&A...596A.108P, Heinrich2016,Reichardt_2016,Planck2018arXiv180706209P, Millea2018A&A...617A..96M, Hazra19}).
As mentioned in the introduction, using different EoR basis sets to infer $\tau_e$ from the CMB PS  and then from $\tau_e$ to astrophysical parameters could result in biases.

Constraining astrophysics directly using the CMB PS bypasses this issue.  For instance,
\citet{Miranda2017MNRAS.467.4050M} adopted a simple EoR galaxy model, parametrized by the ionizing efficiencies of Pop-II or Pop-III dominated galaxies, and constrained their model parameters directly from the {\it Planck} 2015 \citep{Planck2016A&A...596A.108P} PS observations (see also \citealt{Hazra19}).
The top panel of Fig. \ref{fig:miranda} shows their best-fit models, where
\begin{enumerate}
	\item \textit{Pop-II} only considers UV ionizing photons from Pop-II star-dominated galaxies;
	\item \textit{Pop-III} considers UV ionizing photons from both Pop-II and Pop-III star-dominated galaxies; and
	\item \textit{Pop-III, self-regulated} is similar to the previous model but assumes a significant
	contribution from Pop-III stars in the early universe ($z{\sim20}$) before their formation
	becomes completely quenched when $n_e/n_{\rm H}$ reaches 0.2.
\end{enumerate}
Also shown is a \textit{Tanh} model that assumes a fairly sharp transition at
$z_{\rm re}{\sim}10$ (see equation \ref{eq:tanh}). Comparing the corresponding large-scale E-mode polarization PS of these
models against the {\it Planck} 2015 measurement \citep{Planck2016A&A...596A.108P},
\citet{Miranda2017MNRAS.467.4050M} concluded {\it Planck} 2015 might favour a significant
UV ionizing photon contribution from Pop-III star-dominated galaxies in the early Universe and
the resulting optical depth is much higher than
what the default {\it Tanh} parametrization suggests.
This example illustrates how a reionization model informed only by the CMB optical depth could result in biased constraints, compared to using the PS directly in the inference.

\subsubsection{From Planck 2015 to 2018}\label{subsubsec:planck}

From their 2015 to 2018 data release, the {\it Planck} collaboration has made
tremendous efforts in improving the characterization and removal of systematic uncertainties
affecting the polarization data of the {\it Planck} High Frequency Instrument (HFI) on large angular
scales \citep{Planck2016A&A...596A.107P}. With such improvements, the {\it Planck} collaboration has
shown that the mean value and uncertainties of the optical depth deduced from the low-$\ell$
data assuming the {\it Tanh} EoR model have significantly decreased from $\tau =0.078 \pm 0.010$ (the {\it TT+lowP} reconstruction in \citealt{Planck2016A&A...594A..13P})
to $\tau=0.0522\pm0.0080$ ({\it TT+lowE} in \citealt{Planck2018arXiv180706209P}).

Before discussing the difference between EoR inference from the CMB optical depth
and PS, we revisit the 4 best-fit models from \citet{Miranda2017MNRAS.467.4050M} mentioned above,
using the updated measurement from {\it Planck} 2018 \citep{Planck2019arXiv190712875P} to
see if the latest data still supports an earlier reionization and larger optical depths.

\begin{figure}
	\begin{minipage}{\columnwidth}
		\begin{center}
			\includegraphics[width=0.97\textwidth]{figs/Miranda.pdf}\vspace*{-2mm}
		\end{center}
	\end{minipage}
	\caption{\label{fig:miranda} Reanalysing models from \citet{Miranda2017MNRAS.467.4050M}
		with cosmological parameters given by the {\it TT+lowP} reconstruction in
		\citet{Planck2018arXiv180706209P}. \textit{Top panel:} the best-fit EoR models from
		\citet{Miranda2017MNRAS.467.4050M}. \textit{Bottom panel:} corresponding PS of
		the E-mode polarization anisotropies with {\it Planck} 2015 and 2018 measurements
		shown as squares and circles, respectively.  Only the large-scale ($2{\le} \ell{\le}29$) PS is considered when
		evaluating the likelihoods as these scales are most relevant for the EoR.}
\end{figure}
\begin{table}
	\caption{Revisiting the models from \citet{Miranda2017MNRAS.467.4050M}, using both {\it Planck} 2015 and 2018 data. The models correspond to those shown in Fig. \ref{fig:miranda}; the parameters are fixed and are not allowed to vary so as to better fit the 2018 data.  For each model, we list the optical depth and the difference in the reduced $\chi^2$ with respect to the {\it Pop-II} model (computed using the {\it Planck} likelihood code).
}
	\vspace{-3mm}
	\begin{threeparttable}
		\label{tab:planck}
		\begin{tabular}{l|cc|cc}
			\hline \hline

			\multirow{2}{*}{Model} & \multicolumn{2}{c|}{{\it Planck} 2015\tnote{a}}& \multicolumn{2}{c}{{\it Planck} 2018\tnote{b}}\\\cline{2-5}
			& $\tau_e$&$\chi^2-\chi^2_{\rm PopII}$& $\tau_e$\tnote{c}&$\chi^2-\chi^2_{\rm PopII}$\\\hline
			Tanh&0.0792&1.13&0.0788&-1.14\\%\hline
			Pop-II\tnote{d}&0.0832&0&0.0827&0\\%\hline
			Pop-III&0.0926&-0.96&0.0921&5.47\\%\hline
			Pop-III,self-&\multirow{2}{*}{0.1049}&\multirow{2}{*}{-2.27}&\multirow{2}{*}{0.1043}&\multirow{2}{*}{16.10}\\
			regulated&&&&\\%\hline
			\hline
		\end{tabular}
		\begin{tablenotes}
			\item[a] This column assumes cosmological parameters (i.e. the density and Hubble constant) from the {\it TT+lowP} reconstruction in \citet{Planck2016A&A...594A..13P} and calculates $\chi^2$ using the {\it \scriptsize low\_l/bflike/lowl\_SMW\_70\_dx11d\_2014\_10\_03\_v5c\_Ap} likelihood in \textsc{plc} 2.0 \citep{Planck2016A&A...594A..11P}. This likelihood considers both EE, TT, TE and BB components.
			\item[b] Results in this column, similarly to the {\it Planck} 2015 column, assume cosmological parameters from the {\it TT+lowE} reconstruction in \citet{Planck2018arXiv180706209P}, and use the {\it \scriptsize low\_l/simall/simall\_100x143\_offlike5\_EE\_Aplanck\_B} likelihood in \textsc{plc} 3.0 \citep{Planck2019arXiv190712875P}. This likelihood allows users to consider only the EE PS, which is used in this work.
			\item[c] For each model, there are slight changes in the resulting optical depth when comparing the {\it Planck} 2015 and 2018 columns. This is due to the variation in the Hubble constant (see equation \ref{eq:tau_e}) when different cosmological parameters are adopted.
			\item[d] $\chi^2_{\rm PopII} {=} 10492.45$ (407.50) for {\it Planck} 2015 (2018).
		\end{tablenotes}
	\end{threeparttable}
\end{table}

We show the PS of the E-mode polarization anisotropies using the
\citet{Miranda2017MNRAS.467.4050M} EoR models in the bottom panel of
Fig. \ref{fig:miranda} assuming {\it Planck} 2018 {\it TT+lowE} cosmology
($\Omega_{\mathrm{m}}, \Omega_{\mathrm{b}}, \Omega_{\mathrm{\Lambda}}, h, \sigma_8$ =
0.321, 0.04952, 0.679, 0.6688, 0.8118; \citealt{Planck2018arXiv180706209P}). We present the
relative $\chi^2{\equiv}-2\ln \mathscr{L}$ and the corresponding optical depth\footnote{Note that the updated optical depths from the 4 models presented in
		Fig. \ref{fig:miranda} using {\it Planck} 2018 cosmology are all higher than the reported
		value ($0.0522\pm0.0080$; \citealt{Planck2018arXiv180706209P}). This is due to the fact that the model parameters were chosen to fit the 2015 data, not the 2018. Here we do not vary the parameters of the
		\citet{Miranda2017MNRAS.467.4050M} models but instead use a different parametrization introduced below.} in Table \ref{tab:planck}, together
with the results using {\it Planck} 2015 {\it TT+lowP} cosmology
($\Omega_{\mathrm{m}}, \Omega_{\mathrm{b}}, \Omega_{\mathrm{\Lambda}}, h, \sigma_8$ =
0.315, 0.04904, 0.685, 0.6731, 0.829; \citealt{Planck2016A&A...594A..13P}) for
comparison. We conclude that, unlike {\it Planck} 2015, the
2018 measurement no longer prefers an earlier reionization or a significant contribution from
Pop-III stars to early reionization -- the likelihood decreases when ionization
starts at earlier times (see also, e.g. \citealt{Millea2018A&A...617A..96M}).
This is mainly driven by the reduced amplitude at multipole
$\ell{\sim}10 - 20$ in the updated {\it Planck} E-mode PS (see the bottom panel of Fig. \ref{fig:miranda}).
Consequently, models with a better fit to the
updated PS measurement also return an optical depth closer to the reported value from
\citet{Planck2018arXiv180706209P}, indicating that EoR inference biases
from using the CMB optical depth might be insignificant in {\it Planck} 2018 (see also v2 of \citealt{Hazra19}).

For the remainder of the paper, we focus only on the {\it Planck} 2018 data.  Using a flexible, physical EoR model, we forward-model the EE PS, and quantify the bias from using the $\tau_e$ summary statistic instead of the PS.

\section{EoR inference from the CMB}\label{sec:model_eor}

To compute the impact of realistic reionization histories on the CMB PS,
we connect {{21\textsc{cmFAST}}} (\citealt{Mesinger2007ApJ...669..663M,Mesinger2011MNRAS.411..955M}; Murray et al. in prep)
to \textsc{class} \citep{Lesgourgues2011arXiv1104.2932L}.
Specifically, for a given set of cosmological and astrophysical parameters, {{21\textsc{cmFAST}}} performs a 3D EoR simulation\footnote{We do not explicitly model helium reionization. Instead we assume
helium to be singly ionized following the same rate as hydrogen
before becoming fully ionized at $z{=}3$
\citep{Hogan1997AJ....113.1495H,Worseck2011ApJ...733L..24W,Worseck2016ApJ...825..144W}.}.  The volume-averaged EoR history of that simulation is passed to \textsc{class}, which then computes the corresponding CMB PS.

 With this interface, we use the MCMC driver {21\textsc{CMMC}} \citep{Greig2015MNRAS.449.4246G, Greig2017MNRAS.472.2651G}
to sample two EoR models: (i) the astrophysical galaxy-based parametrization native to {{21\textsc{cmFAST}}} (see below); and (ii) the commonly used, two-parameter {\it Tanh} model (see equation \ref{eq:tanh}).  In addition to varying the EoR parameterization, we also check how the inference is affected by the choice of likelihood statistic: (i) $\tau_e$; or (ii) the EE PS.  Specifically, for (i) we use $\tau_e=0.0522\pm0.0080$ from the {\it TT+lowE} result in \citet{Planck2018arXiv180706209P}; for (ii) we use the large-scale ($2{\le} \ell {\le} 29$) E-mode polarization measurements and the likelihood from \citet{Planck2019arXiv190712875P}.

\subsection{Modelling the EoR}\label{subsec:21cmfast}
Our EoR model describes galaxy properties mostly using power-law scaling relations with respect
to their host halo masses \citep{Park2019MNRAS.484..933P} and calculates the 3D reionization evolution following an excursion-set approach based on the cumulative number density of ionizing photons and recombinations (e.g. \citealt{Furlanetto2004ApJ...613....1F, Sobacchi2014MNRAS.440.1662S}).
We
introduce some basic characteristics of our model below and refer interested readers
to the aforementioned references for more details.

We start from an initial Gaussian realization of the density and velocity fields
\citep{Mesinger2007ApJ...669..663M} in a large-volume (250Mpc), grid-based,
high-resolution (${\sim}$0.65Mpc; i.e 250Mpc/384) simulation box assuming periodic
boundary conditions. These fields are then evolved according to second-order Lagrangian
perturbation theory \citep{Scoccimarro1998MNRAS.299.1097S}, and re-gridded to a lower resolution
(${\sim}$1.95Mpc; i.e. 250Mpc/128) for the sake of computing efficiency.
Then, for each cell centred at a spatial position and redshift of $({\bf r}, z)$,
we compare the cumulative number per baryon of ionizing photons\footnote{In this work, we do not consider ionization by X-ray photons or earlier objects such as minihalo-hosted galaxies, and their corresponding parameters.  Although efficient at heating the IGM before reionization, X-rays and minihalos are expected to have a very minor contribution to the EoR for reasonable galaxy models (e.g. \citealt{RO04, McQuinn12, Mesinger2013MNRAS.431..621M, Ross17, Eide18,Qin2020arXiv200304442Q}).} ($\bar{n}_{\rm ion}$) to that of recombinations ($\bar{n}_{\rm rec}$) in spheres with decreasing radii \citep{Furlanetto2004ApJ...613....1F, Sobacchi2014MNRAS.440.1662S}. A cell is ionized if at any radius $R$,
\begin{equation}\label{eq:ionization}
\bar{n}_{\rm ion} \ge (1+\bar{n}_{\rm rec}).
\end{equation}
Unresolved $\mathrm{H}\textsc{ii}$ regions (smaller than the cell size) are accounted for according to \citet{Zahn2011MNRAS.414..727Z}.

The cumulative number of ionizing photons per baryon is obtained with
\begin{equation}\label{eq:n_ion}
\bar{n}_{\rm ion}\left({\bf r},z|R,\delta_{\rm R|_{{\bf r},z}}\right) {=}  \int {\rm d}M_{\rm vir} {\phi}f_{\rm duty} f_* \frac{\Omega_{\mathrm{b}}}{\Omega_{\mathrm{m}}}\frac{M_{\rm vir}}{\rho_{\rm b}}n_{\gamma} f_{\rm esc}
\end{equation}
where $\delta_{\rm R|_{{\bf r},z}}{\equiv}\rho_{\rm b}/\bar{\rho}_{\rm b} {-}1 $ is the
average overdensity within the spherical region, $\rho_{\rm b}$
and $\bar{\rho}_{\rm b}$ represent the baryonic density and its cosmic mean. In equation (\ref{eq:n_ion})
\begin{enumerate}
	\item $M_{\rm vir}$ and
$\phi\left(M_{\rm vir},z|R, \delta_{\rm R|_{{\bf r},z}}\right)$ are the halo mass and halo
mass function;
	\item $f_{\rm duty}\left(M_{\rm vir}\right) {=} \exp\left(-{M_{\rm turn}}/{M_{\rm vir}}\right)$, with a characteristic mass ($M_{\rm turn}$) as a free parameter, accounts for a decreasing occupation fraction of star forming galaxies inside smaller halos due to inefficient cooling, photo-heating feedback
\citep{Efstathiou1992MNRAS.256P..43E,Shapiro1994ApJ...427...25S,Thoul1996ApJ...465..608T,Hui1997MNRAS.292...27H,Sobacchi2013MNRAS.432L..51S,Sobacchi2014MNRAS.440.1662S} or supernovae feedback \citep{Hopkins2014MNRAS.445..581H, Hopkins2017, Wyithe2013MNRAS.428.2741W,Sun2015,Mutch2016};
 \item $f_*\left(M_{\rm vir}\right) {=} \min\left[1, f_{*,10} \left({M_{\rm vir}}/{10^{10}{{\rm M}_\odot}}\right)^{\alpha_*}\right]$ is the fraction of galactic gas in stars and is assumed to scale with the host halo mass
\citep{Moster2013MNRAS.428.3121M,Sun2015,Mutch2016,Tacchella2018ApJ...868...92T,Behroozi2019MNRAS.488.3143B,Yung2019MNRAS.490.2855Y}
according to the two free parameters, $f_{*,10}$ and $\alpha_*$, representing
the normalization and power-law index;
\item $n_\gamma$ is the number of ionizing
photons intrinsically emitted per stellar baryon; and
\item $f_{\rm esc}\left(M_{\rm vir}\right) = \min\left[1, f_{{\rm esc},10} \left(\frac{M_{\rm vir}}{10^{10}{{\rm M}_\odot}}\right)^{\alpha_{\rm esc}}\right]$
is the ionizing escape fraction defined as the number ratio of photons reaching
the IGM to those emitted in the galaxy, and is also assumed to scale with the
halo mass
\citep{Ferrara2013MNRAS.431.2826F,Kimm2014,Paardekooper2015MNRAS.451.2544P,Xu2016ApJ...833...84X}.
\end{enumerate}

Inside the $\mathrm{H}\textsc{ii}$ regions, we estimate the local, average photoionization rate with
\citep{Sobacchi2014MNRAS.440.1662S}:
\begin{equation}\label{eq:gamma}
{\bar{\Gamma}_{\rm ion}\left({\bf r},z\right)} = \left(1+z\right)^2 R\sigma_{\rm H}\frac{\alpha_{\rm UVB}}{\beta_{\rm H}+\alpha_{\rm UVB}}\frac{\bar{\rho}_{\rm b}}{m_{\rm p}}{\dot{\bar{n}}_{\rm ion}},
\end{equation}
where $\sigma_{\rm H}=6.3\times10^{-18}{\rm cm}^2$ and $\beta_{\rm H}{\sim}2.75$ are the photoionization cross-section at Lyman limit and its spectral dependence; $\alpha_{\rm UVB}{\sim}5$, $m_{\rm p}$ and $\dot{\bar{n}}_{\rm ion}$ are
the spectral indices of a stellar-driven UV ionizing
background \citep{Thoul1996ApJ...465..608T}, the mass of a proton, and the local production rate of ionizing photons. Assuming the
typical star formation time-scale is $t_*H^{-1}\left(z\right)$, with $t_*$ being a free parameter,
we calculate the average
star formation rate (SFR) of galaxies in halos of a given mass at a given redshift by
\begin{equation}\label{eq:sfr}
{\rm SFR}\left(M_{\rm vir},z\right) = \frac{M_*}{t_* H(z)^{-1}},
\end{equation}
convert it to the non-ionizing UV luminosity via $L_{1500}/{\rm SFR}=8.7\times10^{27}{\rm erg\ s^{-1}Hz^{-1}{{\rm M}_\odot}^{-1}yr}$ \citep{Madau2014ARA&A..52..415M},
and estimate $\dot{\bar{n}}_{\rm ion}$ using equation (\ref{eq:n_ion}) with $M_*$
being replaced by the SFR.

We follow \citet{Sobacchi2014MNRAS.440.1662S} and estimate the recombination
rate in each cell with a spatial position and redshift of ({\bf r, $z^\prime$}) as well
as an overdensity of $\Delta_{\rm cell}$ by
\begin{equation}
\dot{n}_{\rm rec}\left({\bf r},z^\prime \right){=}\int{\rm d}\Delta_{\rm sub}\phi_{\rm sub}\alpha_{\rm B}f_{\rm H} \frac{\bar{\rho}_{\rm b}}{m_{\rm p}} \frac{\Delta_{\rm sub}^2}{\Delta_{\rm cell}} \left(1{-}x_{\rm \mathrm{H}\textsc{i}, sub}\right)^2.
\end{equation}
Here $\Delta_{\rm sub}$,
$\phi_{\rm sub}\left(z^\prime, \Delta_{\rm sub}| \Delta_{\rm cell}\right)$, $\alpha_{\rm B}$,
$f_{\rm H}$ and $x_{\mathrm{H}\textsc{i}, {\rm sub}}\left(z^\prime, \Delta_{\rm sub}, T_{\rm g}, \bar{\Gamma}_{\rm ion}\right)$
are the sub-grid (unresolved) overdensity, its probability distribution function (PDF; e.g. \citealt{Miralda2000ApJ...530....1M}) within our
large-scale simulation cell (${\sim}2$Mpc), case-B recombination coefficient evaluated at $10^4{\rm K}$,
number fraction of hydrogen in the Universe, and the neutral hydrogen fraction of the sub-grid gas element, respectively.
We assume photoionization
equilibrium in the ionized IGM, accounting for the attenuation of the
local photoionization rate, $\bar{\Gamma}_{\rm ion}$, according to the radiative transfer
simulations from \citet{Rahmati2013MNRAS.430.2427R}.
We then compute the cumulative number of recombinations for each cell (see equation \ref{eq:ionization}) by integrating $\dot{n}_{\rm rec}$ from the time the cell was ionized to the redshift of interest.

In summary, our model consists of the following six astrophysical parameters that we sample within the MCMC:
\begin{itemize}
\item $f_{*,10}$, the normalisation of the stellar-to-galactic gas mass relation at $M_{\rm vir}=10^{10}{{\rm M}_\odot}$, sampled with a flat prior in log space between $10^{-3}$ and 1;
\item $\alpha_*$, the power-law index of the stellar-to-galactic gas mass relation, sampled with a flat prior between -0.5 and 1;
\item $f_{\rm esc,10}$, the normalisation of the ionizing escape fraction to halo mass relation at $M_{\rm vir}=10^{10}{{\rm M}_\odot}$, sampled with a flat prior in log space
between $10^{-3}$ and 1;
\item $\alpha_{\rm esc}$, the power-law index of the  ionising escape fraction to halo mass relation, sampled with a flat prior between -1 and 0.5;
\item $M_{\rm turn}$, the turnover halo mass below which the number density of halos hosting star-forming galaxies drops exponentially, sampled with a flat prior in log space between $10^8$ and $10^{10}{{\rm M}_\odot}$;
\item $t_*$, the star-formation timescale as a fraction of the Hubble time, sampled with a flat prior between 0 and 1.
\end{itemize}
These prior ranges are chosen based on the physical meaning of the parameters.  For example, fractions must range from 0 to 1, and we observe galaxies inside halos with masses around $10^{10}{{\rm M}_\odot}$ thus setting an upper limit on $M_{\rm turn}$.  More detailed discussion on the parameters and corresponding observational constraints can be found in \citet{Park2019MNRAS.484..933P}.

We stress that our EoR model is both {\it flexible} and {\it physical}; both properties are important for useful inference.  It is {\it flexible} in that it is capable of reproducing (see e.g. \citealt{Park2019MNRAS.484..933P}) the bulk properties and scalings of high-redshift galaxy observations (e.g. \citealt{Bouwens2015a,Bouwens2016,Oesch2018ApJ...855..105O}),
as well as results from more sophisticated semi-numerical models and hydrodynamic simulations (e.g. \citealt{Mutch2016,Xu2016ApJ...833...84X,Tacchella2018ApJ...868...92T,Behroozi2019MNRAS.488.3143B,Yung2019MNRAS.490.2855Y}). It is {\it physical} in that the equations and parameters have a straightforward interpretation in terms of galaxy evolution, allowing us to set physically-motivated priors for the free parameters.

For computational convenience, we fix the cosmology to the best-fit values of {\it Planck} 2018 {\it TT+lowE}
($\Omega_{\mathrm{m}}, \Omega_{\mathrm{b}}, \Omega_{\mathrm{\Lambda}}, h, \sigma_8$ =
0.321, 0.04952, 0.679, 0.6688, 0.8118). In practice, one should co-vary astrophysical and cosmological parameters when performing inference. Fixing cosmological parameters effectively assumes that they are mostly constrained by the
temperature and higher-$\ell$ polarization PS {\color{black}(which are not considered in the
likelihoods in this work)}, while the $\ell\le29$ E-mode PS constrains the EoR.
{\color{black} There are less than 5 percent differences ($<$ 1/3 $\sigma$) in the recovered optical
depth from
ignoring the
temperature component in the likelihood.}
However, the most important
degeneracy affecting the determination of $\tau_e$ is the known ``$A_s\exp(-2\tau_e)$''
degeneracy \citep{Planck2016A&A...596A.107P}.  Using the {\it Tanh} model, we show in
Appendix \ref{sec:appendix} that co-varying $\sigma_8$ following that degeneracy has no
impact on the reconstructed EoR history.

\subsection{Hyperbolic tangent vs. physically motivated EoR model}\label{subsec:tanhvs21cmfast}

We perform MCMC simulations for the two different EoR models ({\it Tanh} and {{21\textsc{cmFAST}}}), both constrained using the low-$\ell$ EE PS. As an additional constraint, we also include the model-independent upper limit on the neutral hydrogen fraction at $z\sim5.9$, $\bar{x}_{\mathrm{H}\textsc{i}}<0.06{+}0.05(1\sigma)$, measured from the dark fraction in QSO spectra \citep{McGreer2015MNRAS.447..499M} and modelling the associated likelihood as a one-sided Gaussian \citep{Greig2017}.

In Fig. \ref{fig:MCMC_tanh}, we show the resulting posterior for the {\it Tanh} model using orange shaded regions, while the posterior corresponding to the astrophysical model is denoted with purple lines ([14, 86] and [2.3, 97.7] percentiles).  The recovered EE PS, $\tau_e$, and EoR history are shown clockwise from the top panel.

In the top panel, we see that the differences in the recovered PS are negligible between the two models.  Neither model is able to recover the excess power at $20\lsim \ell \lsim 30$. However, the probability-to-exceed of these data as computed in \cite{Planck2019arXiv190712875P} points to a statistical fluctuation rather than a residual systematic error (or new physics not captured by our EoR modelling).

\begin{figure*}
	\begin{minipage}{\textwidth}
		\begin{center}
			\includegraphics[width=0.90\textwidth]{./figs/posterior_tanh.pdf
		\end{center}
	\end{minipage}
	\caption{\label{fig:MCMC_tanh} Posteriors of our astrophysical EoR model ({\it purple}) and the {\it Tanh} EoR model ({\it orange}), recovered from the E-mode polarization PS on large scales ($2\le \ell \le 29$; \citealt{Planck2019arXiv190712875P}; grey circles shown in panel 1) and the upper limit on the neutral hydrogen fraction at $z\sim5.9$ (\citealt{McGreer2015MNRAS.447..499M}; grey shaded regions shown in panel 2a).
		The 14 to 86 (2.3 to 97.7) percentiles of (1) the E-mode polarization PS and (2) the mean neutral hydrogen fraction evolution in the posteriors are presented using thick (thin) lines or dark (light) shaded regions. PDFs of (2a) the neutral hydrogen fraction at $z=5.9$ and (3) the Thompson optical depth are also shown for the two MCMC results.
		Although the recovered PS and $\tau_e$ are comparable between the two models, the astrophysical model (based on the hierarchical structural growth) results in asymmetric EoR histories with tails towards high redshifts.
	}
\end{figure*}

In the bottom right panel, we see that the distributions of $\tau_e$ peak on similar scales, though the shape of the PDFs are qualitatively different.  The recovered median and [14,86] percentiles are comparable (i.e. $\tau_e = 0.0533^{+0.0078}_{-0.0068}$ for {\it Tanh}; $\tau_e = 0.0556^{+0.0080}_{-0.0064}$ for the astrophysical model).  However, we note that the astrophysical model results in an asymmetric PDF of $\tau_e$, with a tail extending towards high values.

The reason for this is apparent looking at the recovered EoR histories in the bottom left panel.  The two distributions are comparable around the midpoint of reionization, where a large fraction of the EE power is imprinted.  However, the astrophysical model (based on the growth of structure) results in asymmetric EoR histories with a tail towards higher redshifts.  For example, at $z=10$, the astrophysical model recovers $\bar{x}_{\mathrm{H}\textsc{i}} \gsim 0.849 (1 \sigma)$, while the {\it Tanh} model recovers $\bar{x}_{\mathrm{H}\textsc{i}} \gsim 0.925 (1 \sigma)$.
Thus, although the preference for a higher median optical depth is reduced in {\it Planck} 2018 compared to the 2015 data (see section \ref{subsubsec:planck}), the {\it natural shape} of the EoR history implied by the growth of structure does result in a (modest) high-redshift tail.

\subsection{Optical depth vs. power spectra}\label{subsec:tauvsps}

In the previous section we discussed how inference from the {\it Planck} 2018 E-mode PS is affected by the choice of EoR models -- {\it Tanh} and astrophysical.  In this section, we only use the astrophysical EoR model, and instead investigate the impact of the choice of likelihood statistic -- using the EE PS directly vs using the $\tau_e$ summary statistic.  In addition to the choice of CMB statistics, we also account for the (model independent) upper limit on the neutral hydrogen fraction at $z\sim5.9$ from the dark fraction in QSO spectra \citep{McGreer2015MNRAS.447..499M} as well as the galaxy UV luminosity functions (LFs) at $z{=}6{-}10$ from \citet{Bouwens2015a,Bouwens2016,Oesch2018ApJ...855..105O}; both of these independent data sets are included as priors in the total likelihood, as described in \citet{Park2019MNRAS.484..933P}.

In Fig. \ref{fig:MCMC}, we show the marginalized posteriors for the astrophysical parameters (panels in {\it lower left corner}), together with the corresponding: (1) EE PS; (2) EoR history; and (3) $\tau_e$  ({\it upper right}).  We run three different MCMC simulations corresponding to different combinations of observational data sets used for the likelihood:
\begin{enumerate}
	\item \textit{DarkFraction\_LF}, including high-redshift LFs and the QSO dark fraction upper limit on $\bar{x}_{\mathrm{H}\textsc{i}}(z=5.9)$. This run does not consider any CMB observations;
     \item \textit{DarkFraction\_tau\_LF}, based on \textit{DarkFraction\_LF}, but including an additional constraint on $\tau_e = 0.0522\pm0.0080$, which is taken from the {\it TT+lowE} reconstruction in {\it Planck} 2018 \citep{Planck2018arXiv180706209P} and was generated using a {\it Tanh} basis set; and
     \item \textit{DarkFraction\_EE\_LF}, based on \textit{DarkFraction\_LF},  but also forward-modelling the low-$\ell$ E-mode PS using \textsc{class} \citep{Lesgourgues2011arXiv1104.2932L} and includes {\color{black}the low-$\ell$ EE likelihood from {\it Planck} 2018} \citep{Planck2019arXiv190712875P}. These are our ``flagship'' constraints.
\end{enumerate}

\begin{figure*}
	\begin{minipage}{\textwidth}
		\begin{center}
			\includegraphics[width=\textwidth]{./figs/posterior.pdf
		\end{center}
	\end{minipage}
	\caption{\label{fig:MCMC}Marginalized posterior distributions of the astrophysical model
		parameters with different observational constraints:
		(i) \textit{DarkFraction\_LF} (brown shaded regions) uses high-redshift LFs and the QSO dark fraction measurements;
		(ii) \textit{DarkFraction\_tau\_LF} (red dashed lines) uses high-redshift LFs, the QSO dark fraction measurements and $\tau_e$ derived using a {\it Tanh} EoR model; and
		(iii) \textit{DarkFraction\_EE\_LF} (blue solid lines) uses high-redshift LFs, the QSO dark fraction measurements, and the low-$\ell$ EE PS.
		The 2D distributions correspond to 68th (dark regions or thick lines) and 95th (light regions or thin lines) percentiles.
		The upper-right three sub-panels
		present the [14, 86] (and [2.3, 97.7]) percentiles of (1) the E-mode
		polarization PS and (2) evolution of the mean neutral hydrogen fraction
		($\bar{x}_{\mathrm{H}\textsc{i}}$), as well as the PDFs of (2a) $x_{\mathrm{H}\textsc{i}}$
		at $z{=}5.9$ and (3) the Thompson scattering optical depth.
		The median and [14, 86] percentiles of the inferred
		optical depth are presented in the lower right box (4).
		Observational constraints are indicated in
		grey. Overall, we see {\it insignificant difference} in the posteriors
		between
		\textit{DarkFraction\_tau\_LF} and \textit{DarkFraction\_EE\_LF}, indicating a negligible bias in inference when using the $\tau_e$ summary statistic generated with a different EoR basis set (instead of directly forward-modelling the EE PS).}
\end{figure*}

Looking at the posterior of \textit{DarkFraction\_LF}, we see that some of our astrophysical parameters are already constrained by galaxy and QSO observations.  As pointed out in \citet{Park2019MNRAS.484..933P} (see also e.g. \citealt{Tacchella2018ApJ...868...92T,Behroozi2019MNRAS.488.3143B,Yung2019MNRAS.490.2855Y}), the observed high-redshift UV LFs already constrain the SFR-to-halo mass relation to within factors of ${\sim}2$ (i.e. $f_{*, 10}/t_\ast$ and $\alpha_\ast$) and provide an upper limit on the characteristic halo mass below which the galaxy occupancy fraction starts decreasing ($M_{\rm turn}\lsim 10^{9.5} {{\rm M}_\odot}$; 2$\sigma$).
Additionally, the dark fraction measurement of QSO spectra sets a lower limit for the ionizing escape fraction normalization ($f_{\rm esc, 10} > 0.02$; 2$\sigma$), requiring the bulk of reionization to occur before $z\sim6$ (panel 2) and setting a lower limit on $\tau_e$ (panel 3). However, without the CMB, the early
stages of reionization are unconstrained (c.f. \citealt{Greig2017}).  This is evident from the broad range of EoR histories allowed beyond $z\gtrsim8$
(panel 2), as well as the broad distributions of $\tau_e$ (panel 3) and the PS (panel 1).  These early EoR models generally correspond to the high $f_{\rm esc}$ + low $\alpha_{\rm esc}$ corner of astrophysical parameter space.\footnote{We note that if our prior ranges were extended even further, reionization would be allowed at even higher redshifts.  However, an important benefit of using an astrophysical EoR model is that it allows us to place physically-motivated priors on the parameters.  For example, an ionizing escape fraction cannot be higher than unity, nor can star formation occur efficiently inside halos whose virial temperature is smaller than available gas cooling channels.  This is not the case for non-physical or so-called model-independent constraints, for which it can be difficult to choose reasonable priors on the model parameters.}

Including CMB observations rules out early reionizing models.  Both \textit{DarkFraction\_tau\_LF} and \textit{DarkFraction\_EE\_LF} posteriors disfavour the high $f_{\rm esc}$ and low $\alpha_{\rm esc}$ corner of parameter space. The EE PS, EoR history and $\tau_e$ distributions all shrink.

Comparing the  \textit{DarkFraction\_tau\_LF} and \textit{DarkFraction\_EE\_LF} posteriors in Fig. \ref{fig:MCMC} quantifies the bias of using the $\tau_e$ summary statistic, generated with a different EoR model, instead of directly forward-modelling the CMB PS.  There is a small difference in the EE PS with \textit{DarkFraction\_EE\_LF} allowing for a slightly earlier EoR. In general however the two posteriors are nearly identical.  This indicates that, with the {\it Planck} 2018 data, the bias in using $\tau_e$ for the likelihood instead of the EE PS directly is negligible.

{\color{black}
Finally, Fig. \ref{fig:gamma} presents the constrained photonionization rate ($\Gamma_{\rm ion}$; see equation \ref{eq:gamma}).  This serves to further illustrate how a physical model allows us to predict additional IGM properties.  We only show {\it DarkFraction\_tau\_LF}, as {\it DarkFraction\_EE\_LF} is nearly identical, while not including the CMB observations allows $\Gamma_{\rm ion}$ to become unrealistically large.
The rise in $\Gamma_{\rm ion}$ with redshift is determined by the formation of structure in our galaxy model, with the ``flattening'' seen at $z\lsim$6--7 being due to photo-heating suppression of gas accretion onto galaxies following reionization
(e.g. \citealt{Sobacchi2013MNRAS.432L..51S}).
We see that while our prediction is consistent with the measured UV ionizing background at $z\sim 5$ to 6 \citep{Bolton2007MNRAS.382..325B,Calverley2011MNRAS.412.2543C,Wyithe2011MNRAS.412.1926W}, observational data lies on the lower boundary.
Additionally including $z\sim5$--6 Ly$\alpha$ transmission statistics constrains the upper envelope of $\Gamma_{\rm ion}$ significantly, though the results are more model dependent than those presented here (Qin et al. in prep).
}

\begin{figure}
	\begin{minipage}{\columnwidth}
		\begin{center}
			\includegraphics[width=\textwidth]{./figs/gamma.pdf
		\end{center}
	\end{minipage}
	\caption{\label{fig:gamma}{\color{black}The [14, 86] (and [2.3, 97.7]) percentiles of the photonionization rate from models used in the {\it DarkFraction\_tau\_LF} posterior in Fig. \ref{fig:MCMC}. Observational constraints \citep{Bolton2007MNRAS.382..325B,Wyithe2011MNRAS.412.1926W,Calverley2011MNRAS.412.2543C} are indicated in grey.}
	}
\end{figure}

We should caution however that our findings here are valid in the context of a flat $\Lambda$CDM Universe. An important follow-up to this study will be to generalize these trends to alternative cosmologies (e.g. \citealt{Paoletti2020arXiv200512222P}).
It is possible that in some exotic cosmologies, the correlation between parameters describing exotic physics and reionization exist at the level of $\tau_e$ but are broken once the full power of the CMB PS is included, rendering the necessity of joint analysis of reionization data and CMB observations.

\section{Conclusions}\label{sec:conclusion}

In this work, we develop an interface between {{21\textsc{cmFAST}}} and \textsc{class}, allowing us to forward-model the large-scale E-mode polarization power spectra inside the {21\textsc{CMMC}} sampler.  With this setup, we study how the choice of (i) EoR model (astrophysical vs. a {\it Tanh}), and (ii) likelihood statistic ($\tau_e$ vs EE PS), impacts EoR parameter inference.

The marginalized posteriors of $\tau_e$ for the {\it Planck} 2018 data \citep{Planck2019arXiv190712875P} are fairly insensitive to the parametrization of the EoR history: {\it Tanh} vs a galaxy model.  This is contrary to claims based on the earlier 2015 release (e.g. \citealt{Miranda2017MNRAS.467.4050M,Planck2018arXiv180706209P,Millea2018A&A...617A..96M,Hazra19}). However, the galaxy model (based on standard hierarchical growth of structure) results in asymmetric EoR histories, with the early stages extending to higher redshifts (see also \citealt{Choudhury2006MNRAS.371L..55C,Wise2014,Price2016,Qin2017a,Gorce2018A&A...616A.113G,Kulkarni2018}).  As a result, the lower limit on the volume averaged neutral hydrogen fraction at $z=10$ changes from $\bar{x}_{\mathrm{H}\textsc{i}} \gsim 0.93 (1 \sigma)$ using the {\it Tanh} model to $\bar{x}_{\mathrm{H}\textsc{i}} \gsim 0.85 (1 \sigma)$ using an astrophysical model.

Using our galaxy EoR model, we quantify the bias in inference when the likelihood is computed from the $\tau_e$ summary statistic, compared with directly using the
the measured E-mode polarization PS. We perform MCMC simulations
taking into account the current observational constraints from high-redshift
galaxy UV LFs \citep{Bouwens2015a,Bouwens2016,Oesch2018ApJ...855..105O}
and the model-independent constraints from the dark fraction in QSO spectra \citep{McGreer2015MNRAS.447..499M}.
Additionally including either CMB statistic helps constrain the posterior, ruling out models that have a high escape fraction for faint galaxies and hence an early reionization.  However, the difference between using $\tau_e$ for the likelihood, compared with the EE PS, is negligible.
Our flagship constraints, based on the QSO dark fraction + UV LFs + EE PS, result in an optical depth of
$\tau_e = 0.0569_{-0.0066}^{+0.0081} (1\sigma)$, with asymmetric EoR histories.

\section*{Acknowledgements}
The authors thank V. Miranda for providing data from \citet{Miranda2017MNRAS.467.4050M}.
This work was supported by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (AIDA -- \#638809). The results presented here
reflect the authors’ views; the ERC is not responsible for their use. Parts of this research
were supported by the Australian Research Council Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO 3D), through project \#CE170100013.

\section*{Data availability}
The data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mn2e}