\section{Introduction}
The resilience of communities to crises such as flooding, wildfires, pandemics and other extreme events hinges on a deep understanding and effective simulation of the impacts of the crisis on human and physical environment \cite{Zhu2019}. Urban mobility plays a vital role in community resilience to crises by enabling populations to access critical facilities, such as healthcare, pharmacies, and grocery stores, for instance, \cite{Liu2019}; hence, the ability to simulate and examine the impacts of crises on urban mobility is essential to effectively improve the resilience of cities \cite{Sadri2020}. Studies, including Lei et al. \cite{Lei2020}, have examined the vulnerability and resilience of transportation systems during natural disasters and other crises. The majority of existing studies focus primarily on examining the vulnerability of physical road networks \cite{Fan2020e,Dong2020} and quantifying the effects of crises on traffic patterns \cite{Oh2020,Wang2020b}. A critical missing component is the capability to simulate crisis perturbation scenarios for predictive evaluation of the impacts on movement trajectories and traffic patterns at urban scale to inform emergency-response and resilience-planning decisions. To address this gap, the goal of this study is to create and test a deep learning model that can capture individuals’ movement patterns during normal situations and impose crisis-induced perturbations to simulate impacts on movement trajectories and traffic patterns.

With the wide use of smartphones and apps with location services, a person’s digital footprints can be generated based on daily trips, which allow researchers to gather fine-grained data related to individuals’ movement trajectories during normal situations for model training and prediction \cite{Olmos2018}. Prior studies have attempted to address mobility simulation through statistical models and machine-learning techniques at different scales \cite{Yan2017}. One stream of the existing studies focuses on quantitative assessments of the dynamical and statistical properties of human travel \cite{Kitamura2000}. With more prevalence of location data related to mobile phone users, González et al. \cite{Gonzalez2008} demonstrated that human movement trajectories have a high degree of temporal and spatial regularity, as they follow simple reproducible patterns. This finding has significant implications for the predictability of individual mobility patterns in normal conditions. Furthermore, measuring the entropy of an individual’s trajectory, Song et al. \cite{Song2010} showed a 93\% potential predictability of human mobility under normal situations. The statistical metrics provide essential knowledge about the travel distance, radius, and frequently visited locations of humans \cite{Song2010a}. Despite these advances in the statistical characterization of aggregated human mobility patterns, limited models exist that could capture the temporal sequence of activities and spatial distribution of activities of each individual. The majority of the standard mobility models have limited capability for predicting individual mobility at a local scale, such as a neighborhood, in a specific timestamp.

With advances in machine-learning techniques and enhanced computational power, it is now possible to learn and simulate finer-grained mobility patterns, such as the sequential patterns of individual location histories \cite{Tang2019}. In particular, deep neural networks (DNN) are prevalent in learning and simulating human dynamics and density from large data sets \cite{Li2017b}. For example, Wang et al. \cite{Wang2020a} designed a DNN architecture with an alternative-specific utility using behavioral knowledge to analyze human choices of their travel mode in normal circumstances. Hosseini et al. \cite{KhajehHosseini2019} introduced a deep learning-based methodology to predict traffic state using convolutional neural networks (CNN) by taking the individual vehicle-level data as inputs. The CNN model is further improved to capture citywide crowd activities with parameter efficiency and stability \cite{liang2020revisiting}. These deep learning techniques offer significant improvements in predicting individual mobility patterns with high resolution in terms of locations and timestamps \cite{Wang2020}.

Human mobility, including origin-destination matrixes and the trajectories for each  origin and destination pair, however, is strongly influenced by the condition of road networks, specific gathering events, friendship effects \cite{Cho2011}, and residents’ travel habits and lifestyles \cite{Zong2019}. In particular, crises such as natural disasters cause perturbations in physical infrastructure and roads that influence human movement trajectories. For example, inundated road segments and buildings would necessitate evacuation of affected residents, thus reducing travel demand in certain areas \cite{Fan2020e}. Human trajectory data in the context of crises for model training and mobility simulation are scarce, as this data is rarely recorded in historical record \cite{Deville2014}, necessitating adaptive models for simulating human mobility under crisis scenarios (e.g., flooding, wildfires, and blizzards.) Information related to movement trajectories and traffic patterns in such extreme and complex situations is rarely recorded in historical data. The lack of data makes the simulation of  human mobility during crisis conditions challenging \cite{Lu2012}.

To address this information gap, this study tested an adaptive reinforcement learning model that can (1) learn individual mobility patterns from data collected from regular daily activities; and (2) simulate mobility under extreme events. The model incorporates both time and location factors as the input features and allows users (e.g., practitioners and decision makers) to adjust the environmental parameters for various application objectives in depending upon scenarios. To illustrate the application as well as the performance of the model in both prediction and application, we trained the model using data related to human mobility activities in Houston during March and April 2017, and used the model to simulate movement trajectory densities and traffic patterns during the flooding caused by the Hurricane Harvey in August 2017. The potential of the model for applications in emergency response and pandemic prediction is also discussed.

\section{Related work}
\subsection{Resilience of urban mobility to crises}
The role of urban mobility in a crisis such as access by emergency responders has been emphasized in the existing literature \cite{Wang2015a}. Assessing strategies to enhance the resilience of urban mobility is essential to mitigate the negative effects of crises and to improve the effectiveness of response efforts. To this end, existing studies \cite{Dong2019b,Wang2016} have put forth methodologies to evaluate the vulnerability of physical road networks and have proposed appropriate strategies for enhancement of urban mobility resilience. In another stream of work, Fan et al. developed a mathematical contagion model to predict the spread of floodwaters over road networks that informs about the disruptions of the urban mobility networks in flooding events \cite{Fan2020e}. Dong et al. proposed a machine-learning model that coupled road networks and water channels to examine the impacts of flooding on road networks through the lens of network dependencies \cite{Dong2020}. Despite the success of prior models, the majority of their efforts are focused on physical road networks, such as the accessibility to critical facilities. To expand the understanding of collective mobility patterns in crises, Wang et al. \cite{Wang2016} used empirical data to examine the impact of natural disasters on population travel distances and radius of gyrations. The result of their study indicates inherent resilience of urban mobility in crises. The majority of the existing studies (such as Lu et al., 2012 \cite{Lu2012}) have used empirical approaches or analytical methods (such as network science-based models) to examine mobility patterns when a population is hit by a crisis. A critical missing piece, however, is predictive data-driven models that can realistically simulate the spatial-temporal patterns of urban mobility during a crisis. Such predictive models and their resulting information should inform about the spatial distribution of the vehicles on road networks and emerging traffic jams caused by crisis perturbations, such as road inundations. This information is crucial for emergency response and mobility resilience enhancement; however, using empirical data from past crises, one could only evaluate one scenario of crisis impacts. The key to addressing this gap is a predictive data-driven model that can learn the patterns of movements by individual vehicles during normal situations and which is capable of imposing crisis-related perturbations to simulated changes in mobility patterns. Such models for urban mobility simulation in a crisis situation could offer a unique and effective approach that enables city planners, emergency managers, and decision makers to capture the spatial and temporal patterns of population movements during extreme events, and facilitate contingency and hazard mitigation planning.

\subsection{Mobility destination and trajectory prediction}
Simulating urban mobility is an important problem with a wide range of applications, such as traffic control, accident warning, pandemic prediction, and urban planning \cite{Chao2020}. Therefore, there is a need for general and robust methods to predict the destinations and trajectories for drivers, especially in populated urban areas \cite{Zhao2018a}. A number of prior studies have explored multiple methods to quantify, model and simulate human mobility, ranging from destination (next place) prediction \cite{Noulas2012} to route selection. This section discusses existing research work to show the achievements in addressing the challenge and to highlight the need for adaptive models to simulate urban mobility under crises situations.

First, a number of existing studies have proposed and tested methods for next-place prediction \cite{Ma2013}. Next-place prediction, also called destination prediction in local urban scale, predicts the movement patterns of individuals and accordingly estimates the flow of population in both spatial and temporal manners. In next-place prediction, the majority of the methods deal with the pairs of origins and destinations based on individuals’ historical trip data. As mentioned earlier, statistical evidence shows the presence of regularity in human daily mobility patterns \cite{Song2010}. Hence, to detect and learn this regularity, one commonly adopted approach is to develop clustering algorithms such as k-means and DBSCAN (Density-based Spatial Clustering of Applications with Noise) algorithm based on historical travel data and to identify association rules between the origin and destination \cite{Chu2010}. Multiple studies \cite{Du2019} indicated that the clustering-based approaches have achieved very high performance in next-place prediction under normal circumstances. Location embedding is another method to capture location semantics with comprehensive numerical representations \cite{Cheng2013,WangandLi2017}. For example, Shimizu et al. proposed a place-embedding method that can learn fine-grained representation with spatial hierarchical information and which achieved higher accuracy in predicting the next place of human trips in urban areas \cite{shimizu2020learning}.

Second, in addition to single destination prediction, existing studies have also explored models for sequential patterns of individual travel to places \cite{Zhao2018}. In this task, Markov chain-based models, including Mobility Markov Chain \cite{gambs2012}, Mixed Markov Chain \cite{Asahara2011}, and Hidden Markov Model \cite{Mathew2012}, are commonly used to predict the location sequence of a given person/vehicle. Associated with location sequences but more related to route selection, other research is concerned with predicting the trajectory selected by an individual from an origin to a destination has garnered significant attention over the past few years \cite{Qiao2015}. The predicted trajectories allow researchers to not only estimate travel time \cite{Tang2018}, but also to analyze travel demand over urban road networks \cite{Dabiri2020}. Due to the complexity of trajectory prediction and required data scale, only a limited number of studies have examined this problem. The most commonly used approach is the shortest-path algorithm that constructs routable graphs from historical trajectory data and then computes the route based on travel frequencies \cite{wei2012}. This approach relies on traffic flow on specific road segments, but does not consider origin and destination pairs, that significantly influence route selection \cite{Yang2017}. More recent studies direct attention to the behavioral mobility patterns of individuals by inferring home and workplace, duration of activities. and other relevant temporal and spatial features for route prediction \cite{Alexander2015}. For example, Jiang et al. proposed a mechanistic modeling framework, called TimeGeo, that can generate mobility behaviors with a resolution of 10 minutes and hundreds of meters \cite{Jiang2016}.

In summary, recent works have significantly advanced methods of simulating urban mobility and trajectory prediction under normal conditions based on the evidence of regularity in human movement and activity patterns. Still unsolved, however, is the challenge of adaptability of a model for simulating urban mobility when people are exposed to disruptive conditions, such as life-threatening situations during crises such as flooding, wildfires, and even pandemics. Modeling approaches with limited adaptability would not be useful in assessing movement trajectories and traffic under crisis situations. Advances in adaptive reinforcement learning models provide opportunities for simulating mobility under crises situations. In the following section, we will elaborate the proposed adaptive reinforcement learning model and then show its application in the context of flood impact analysis in Houston during Hurricane Harvey in 2017.

\section{Methodology}
The proposed adaptive reinforcement learning model comprises three modules: destination prediction, trajectory prediction, and crisis scenario application (Figure \ref{fig:1}). The proposed model uses the human movement trajectory data as input, learns mobility patterns in regular situations, then simulates the trajectories and traffic conditions in crisis situations through adjusting the reward table. The mathematics and steps are formulated and elaborated in the following sub-sections.

\begin{figure*}
    \centering
    \includegraphics[scale=0.43]{figures/Figure1.pdf}
    \caption{A schema of the proposed adaptive reinforcement learning model for simulating urban mobility.}
    \label{fig:1}
\end{figure*}

\subsection{Preliminary definitions}
In this section, we first define the notations and terminologies that are used in this study.

\begin{itemize}[leftmargin=.3in]
  \item \textbf{Definition 1}. A \textbf{trajectory $T_i$}, also named as \textbf{waypoints}, \textbf{route} or \textbf{path}, is a spatial trace of a vehicle (or an individual) generated by their mobile apps in geographical space. $T_i$ contains a sequence of locations with specific latitudes and longitudes: $[(x_0^i,y_0^i,t_0^i),...,(x_{k^i}^i,y_{k^i}^i,t_{k^i}^i),...,(x_{n^i}^i,y_{n^i}^i,t_{n^i}^i)]$, where $x_{k^i}^i \in R$ and $y_{k^i}^i \in R$ are the latitude and longitude of a trace record $i$ at the $k^i_th$ location at time $t_{k^i}^i \in R$, and $k^i \in [0,1,...,n^i]$.
  \item \textbf{Definition 2}. A \textbf{trip} $t_i$, also called \textbf{O-D pair}, is a pair of origin and destination for a given trace record $i$. Trips can be extracted from trajectory data.
  \item \textbf{Definition 3}. A \textbf{provider} is a service provider which is allowed to anonymously collect locations from mobile devices. The providers are sorted into four categories: consumer vehicles, taxi/shuttle/town car services, field service/local delivery fleets, and for-hire/private trucking fleets. Each provider’s services fall exclusively into one category.
  \item \textbf{Definition 4}. A \textbf{grid world} is a grid based on a geographical map. The cells in the grid world are in the same shape and with equal area. Each cell represents a specific state (or location) showing is the location of a vehicle. This term is used to address the Markov decision process (MDP).
  \item \textbf{Definition 5}. A \textbf{policy}, $\pi$, is a set of choices of actions for an agent (people or vehicle in this study) at each state.
\end{itemize}

Based on these definitions, the destination prediction problem can be described thus: given an origin with spatial coordinate $(x_0^i,y_0^i,t_0^i)$ and a specific scenario, predict the spatial coordinate of the destination $(x_{n^i}^i,y_{n^i}^i)$. The trajectory prediction problem can further be described as: given a pair of origin $(x_0^i,y_0^i,t_0^i)$ and destination $(x_{n^i}^i,y_{n^i}^i)$ with time information, predict a trajectory, $[(x_0^i,y_0^i,t_0^i  ),...,(x_{k^i}^i,y_{k^i}^i),...,(x_{n^i}^i,y_{n^i}^i)]$, that a vehicle (or an individual) will go through from the origin to the destination. The specific time in a location on the path to the destination is uncertain due to the dynamic traffic states and data collection frequency. Hence, in this study, the trajectory prediction task will predict only the route a vehicle (or an individual) will select from an origin to a destination, while the specific time for the location of the vehicle (or the individual) on the route will not be estimated.

\subsection{Feature generation}
The first module of the model is a tool that can predict the destination of a vehicle, given an origin at a specific time. As suggested by existing studies, human movements are spatial- and time-dependent, meaning that each pair of origin and destination should include both temporal and spatial features to ensure the model is both accurate and general. To achieve this, this section discusses the process of generating relevant features.

For the temporal dimension, each trip sets the start time at the origin, showing the time of the day and day of the week. Empirical results show that urban mobility presents a significant weekly recurrent pattern \cite{Cho2011}. To incorporate this feature and to mitigate the effect of mobility variations across different days in a week, we converted the time of the day and day of the week into an integrated feature, called $t_n$. Specifically, we considered that any time within a week can be represented by a real number from 0 to 1 \cite{Garcia2017}. 0 represents the start of Monday (00:00 a.m., Monday), and 1 represents the end of Sunday (23:59 p.m., Sunday). Each day can be divided into 24 bins based on hours. Then, a specific time can be translated by:

\begin{equation}
   \centering
   t_n=\frac{hours+mins/60+24 \cdot weekday}{7 \cdot 24}
\end{equation}

where $hours$ represents the hours of the time in 24-hour system, $mins$ represents the minutes of the time, and $weekday$ represents the number of the days from Monday (for example, $weekday$=2 for Wednesday). In this way, we can specify the time feature of a trip in a continuous manner across different times in a week.

Since people tend to have recurrent movement patterns during the weekday and special travels during the weekend (Saturday and Sunday), we also considered this difference in the model by using the category of weekday and weekend as a separate feature, $w_n$. The value of the trips in this feature can be represented using the following rule:

\begin{equation}
    \centering
    w_n=
    \left\{ \begin{array}{ll}
            0 & weekday \\
            1 & weekend
        \end{array} \right.
\end{equation}

The location of the origin and destination is usually recognized by the latitude and longitude coordinates. The resolution of the coordinates tends to be in meters. This raises a significant challenge for the predictive models to accurately learn the locations from the training data and predict location for a given origin. Also considering extremely high-resolution coordinates will also increase computational cost for the model. To simplify this process and maintain fairly high resolution, we round the coordinates into numbers with three decimals places, (This threshold is flexible and can be selected based on the requirement of the resolution.) Similar to the idea of grid, the coordinates are associated to a cell and represented by a centroid of the corresponding cell. Furthermore, since the earth is an ellipsoid which should be described as a three-dimensional space, we projected the latitude and longitude to a Euclidean 3D space using the following formula \cite{Garcia2017}:

\begin{equation}
    x_{st} = \cos(x_0^i) \cos(y_0^i), \\
    y_{st} = \cos(x_0^i) \sin(y_0^i), \\
    z_{st} = \sin(x_0^i)
\end{equation}

Through this approach, we can obtain three location features to represent a location of an origin $(x_{st},y_{st},z_{st})$. We also did the same transformation for the coordinates of the destination in the training data and obtain $(x_{ed},y_{ed},z_{ed})$.

\subsection{K-nearest neighbor regression for destination prediction and evaluation}
Once the five features for an origin are prepared, we can train a supervised learning model to learn the movement patterns by extracting the association between origins and destinations from the historical data in regular conditions. Among the commonly used learning methods \cite{Jia2017}, k-nearest neighbor (k-NN) regression demonstrates high performance in predicting the location feature values of the destinations. The k-NN regression method assigns weights based on the contributions of the neighbors for the predicted outcomes (location features in this study). The nearer the neighbors, the more contribution they can offer on the outcome values.

Specifically, the k-NN model will take the list of trips $t_i=(O^i,D^i)$ as the training set, including all time and location features of origin and destination. $O^i=(t_n^i,w_n^i,x_{st}^i,y_{st}^i,z_{st}^i)$, and $D^i=(x_{ed}^i,y_{ed}^i,z_{ed}^i)$. Then, we defined the number ($k$) of nearest neighbors that would contribute to predicting the destination. The value of $k$ should be learned from the training data by minimizing the root mean square error (RMSE) which is defined below. The distance between every two neighbors (two origins: $O^i$ and $O^j$) is calculated using Euclidean norm ($L^2$ norm), $\parallel \cdot \parallel$, as shown below:

\begin{equation}
    \parallel O^i-O^j \parallel = \sqrt{(t_n^i-t_n^j)^2+(w_n^i-w_n^j)^2+(x_{st}^i-x_{st}^j)^2+(y_{st}^i-y_{st}^j)^2+(z_{st}^i-z_{st}^j)^2}
\end{equation}

The location feature values of the predicted destination are the average values of the location feature for the corresponding destinations of the k-nearest origins. This process generates three location feature values for the predicted destination. We then convert the location features back to the latitude and longitude with three decimal places of accuracy, analogous to the centroid of the destination cell. This coordinate is the predicted coordinates for the destination. By adjusting the value of $k$ and calculating corresponding RMSE based on all pairs of predicted coordinates and actual coordinates, we can find an optimized value of $k$ which can minimize the RMSE with high predictive performance.

This process is analogous to addressing a regression problem. Hence, we use the RMSE \cite{Chai2014} to measure the error between the actual and the predicted values:

\begin{equation}
    RMSE=\sqrt{\frac{1}{N} \sum_{i}^{N} \frac{1}{2} [(x_{n^i}^i-{x'}_{n^i}^i)^2 + (y_{n^i}^i-{y'}_{n^i}^i )^2])}
\end{equation}

where ${x'}_{n^i}^i$ and ${y'}_{n^i}^i$ are the latitude and longitude of the predicted destination for a trace record $i$, and $N$ is the number of traces in the test data. Since we calculate the RMSE for each trace, the value of $N$, here, should be 1. The RMSE takes the Euclidean distances between the predicted outcomes and the observations which are then normalized across distances among all test data.

\subsection{Markov decision process and reinforcement learning}
Once the destinations are predicted for given origins, the next step is to estimate the trajectories between the origins and destinations. An effective and commonly used approach is to structure the space of learned policies, which more generally called Markov decision process. The Markov decision process is a stochastic control process with a tuple of $(S,A,P,R)$, where $S$ is a set of states $s_i \in S$, here meaning the cell in a representing vehicle location; $A$ is a set of actions $a \in A$ which includes the four directions a vehicle can move from one cell to another in this context; $P({s'}_i | s_i,a)$ is the transition probability that action $a$ in state $s_i$ for agent $i$ leads to state ${s'}_i$ in the next timestamp. Following a transition matrix $T(s_i,a,{s'}_i)$; and $R(s_i,a,{s'}_i)$ is the reward function that specifies the reward a vehicle will immediately receive after transitioning from state $s_i$ to state ${s'}_i$ due to action $a$.

To simplify the state space and also to enable the matrix computation, we create a grid world that represents all possible states of an agent in the environment. With a start state (i.e., the origin of a vehicle) and the terminate state (i.e., the destination), the objective of this MDP is to learn an optimal policy $\pi:S \rightarrow A$, that can maximize the total rewards for an agent to move in the grid world \cite{ziebart2008maximum}. To learn a policy, we first need to specify the parameters in the model, the probability distribution over a trajectory $P(\zeta | T,\kappa)$ ($\zeta$ is a trajectory; and $\kappa$ is the reward weights), transition matrix $T$, and the reward function $R$. The calculation starts with the transition matrix $T$. In this study, we assume that the agents have the same probability to go in each direction a among the four directions. Hence, the transition probability is set to be equal (the value is 1 if the action can be taken) among the four directions. Since the agent is usually not allowed to move out of the grid and the destination, we set the transition probabilities of some actions at the border and destination cell to be 0.

Then, the probability $P_T$ that an agent will be in the state ${s'}_i$ at next step given a state $s_i$, an action $a$, and a transition matrix $T$, can be computed as $P_T({s'}_i | a,s_i)$. The probability distribution over a trajectory $P(\zeta | T,\kappa)$ would be proportional to the product of all probabilities over the path:

\begin{equation}
    P(\zeta | T,\kappa) \propto \prod_{s_i,a,{s'}_i \in \zeta}  P_T({s'}_i | a,s_i)
\end{equation}

To learn the reward function from historical trajectory data, we maximize the likelihood of the historical trajectory data under the maximum entropy distribution:

\begin{equation}
    \kappa^* = \arg\max_\kappa \sum_{data} \log{P(\widetilde{\zeta} | T, \kappa)}
\end{equation}

where $\widetilde{\zeta}$ specifically refers to the historical trajectories. Solving this function by gradient-based optimization method, we can obtain the reward weights for the grid world. The reward value of a trajectory is simply the sum of the state rewards, which is calculated by the product of the path feature and the reward weights.

\begin{equation}
    R(\mathbf{f}_{\widetilde{\zeta}})= \sum_{s_i \in \zeta} \kappa^\top \mathbf{f}_{s_i}
\end{equation}

where $\mathbf{f}_{\widetilde{\zeta}}$ is the path feature, which is the sum of the grid feature of each state, $\mathbf{f}_{s_i}$. The formula is shown as follows:

\begin{equation}
    \mathbf{f}_{\widetilde{\zeta}} = \sum_{s_i \in \widetilde{\zeta}} \mathbf{f}_{s_i}
\end{equation}

In implementing the model, the first challenge we encounter is the computational cost for training the agent on the grid world. That is because, in each step, when the agent wants to determine the action and next state, the model has to do matrix multiplications. Since we include the whole Houston metropolitan area to create the grid world, the matrix would be very large and, thus, calculating the reward would and informing the agent to make a decision would require an onerously long calculation time. This cost is usually not allowed for the model that wants to be applied to various contexts, especially in crisis settings. In fact, the cells distant from the origin and destination are not to make any impact on the trajectory selection of an agent. To reduce the computational cost, hence, in the implementation process, we draw only a small bounding box that can incorporate both origin and destination. Then, all cells in the bounding box are considered to be in the reward matrix to do reward calculations.

The second challenge would be the effect of the historical trajectories on the immediate reward values. As discussed earlier, the reward values in the reward table are learned from the historical trajectories that pass through the cell. For an origin and destination pair, we aim to find a trajectory to connect these two cells. The historical trajectories that connect other O-D pairs would, however, give weights to other cells in the reward table. Sometimes, the weights from the trajectories for other O-D pairs are quite high, inducing the agent to go other directions, moving around some high reward cells, and even leaving the actual destinations. To mitigate the effect of noisy historical trajectories, we include only the historical trajectories that pass through the given origin and destination to learn the reward table. This solution would effectively overcome the noise in the reward tables, reduce the training time, and finally improve the efficiency and performance of the model in learning reward tables.

\subsection{Transformation of the reward function}
The reward table for each pair of origin and destination is the key to predicting the trajectory that an agent would select under a specific circumstance. Through examining the reward function learned from the Markov decision process, some common problems hinder an agent to find the optimal policy.

First, due to the sparsity of the data, it is often the case that the immediate rewards among all cells are quite close to each other, meaning that an agent would be not able to effectively distinguish the benefits of two states to make an accurate choice of the next state. As a result, the model is usually unstable, generating different trajectories in different rounds of implementation. In addition, the reward values learned by the reinforcement learning model are always positive for all cells in the grid world. The agent tends to get stuck in some cells with relatively high rewards since the next state does not offer substantial reward to enable a moving action. To address these pitfalls, we conducted a normalization to project the values of the immediate rewards to a negative space with a range of $[-1,0]$. The normalized reward matrix is denoted as $R'$. The process can be formulated using the following equation:

\begin{equation}
    {r'}_{ij} = \frac{r_{ij}-r_{max}}{r_{max} - r_{min}}
\end{equation}

where $r_{ij}$ is the immediate reward in cell at $i^{th}$ row and $j^{th}$ column in the grid in according to reward function $R$; $r_{max}$ and $r_{min}$ are the maximum and minimum immediate reward respectively among all cells in the same grid. By doing so, the range of the reward values will be extended and distributed in a greater space. The negative values of the rewards can enable the agent to move out of some cells halfway on the trajectory and get a high reward only at the destination cell.

Second, in the reward table, the reward of an action given to the agent does not always compel the agent to move towards the destination due to a high density of the road networks and concentration of historical trajectories. That means an agent tends to move around the cells with high reward values, sometimes even moving back to the origin. The agent, however, should intentionally move towards the destination. The closer the next state to the destination, the higher reward the agent can obtain. To this end, we created a new matrix, $L$, that can be added up to the reward table to compensate for the effect of proximity to the destination on trajectory selection. Empirically, the reward distribution accounts for the closeness to the destination centers at the destination cell. The reward has to be reduced if the agent does not get to the destination or if the agent moves further from it. The idea of the reward distribution is similar to the Gaussian distribution. Hence, we employed the Gaussian distribution to assign rewards to the cells in the grid world. The reward regarding the closeness to a destination is denoted as $l_{ij}$ for the cell at $i^{th}$ row and $j^{th}$ column. The value of $l_{ij}$ can be obtained using the following formula:

\begin{equation}
    l_{ij}=f(i)+f(j)-\delta
\end{equation}

\begin{equation}
    f(i)=  \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2}(\frac{i-\mu}{\sigma})^2}
\end{equation}

\begin{equation}
    f(j)=  \frac{1}{{\sigma'} \sqrt{2\pi}} e^{-\frac{1}{2}(\frac{j-{\mu'}}{\sigma'})^2}
\end{equation}

where $f(i)$ and $f(j)$ are the reward values calculated based on the row and column; $\mu$ and ${\mu'}$ are the row and columns of the destination, respectively; $\sigma$ and ${\sigma'}$ are calculated as the variance of rows and columns. $\delta$ is a parameter that controls the range of the reward values. Regarding the contribution of the closeness to the destination, the value of $\delta$ can be tuned to improve the prediction performance of the model.

Third, the model usually does not converge at the destination cell because the reward at the destination is not significantly higher than other states. The agent would go back and forth around the destination cell. To address this problem, we simply assigned a large reward to the destination cell so that and the process can be converged when the agent arrives at and stops at the destination. The matrix can be created based on the formula shown as follows:

\begin{equation}
    d_{ij}=
    \left\{ \begin{array}{ll}
            500 & if \; ij \; is \; the \; destination \\
            0 & otherwise
        \end{array} \right.
\end{equation}

where $d_{ij}$ is the extra reward given to the destination cell, while other cell would receive no reward.

After all transformation matrices are prepared, finally, the reward table for a pair of origin and destination can be obtained by combining all of these matrixes:

\begin{equation}
    R_t = {R'} + L + D
\end{equation}

Using this transformed reward table, we can train agents for each origin and destination pair to effectively learn the optimal policy.

\subsection{Optimal trajectory prediction and evaluation}
Once the reward table for each pair of origin and destination is prepared, we define that the expected utility starting in $s_i$ and acting optimally is $V^* (s_i)$. The expected utility starting out having taken the action a from state $s_i$ and thereafter acting optimally is defined as $Q^*(s_i,a)$. Here $Q^*(s_i,a)$ is computed by:

\begin{equation}
    Q^*(s_i,a)= \sum_{{s'}_i} T(s_i,a, {s'}_i)[R_t (s_i,a,{s'}_i) + \gamma V^*({s'}_i)]
\end{equation}

where $\gamma$ is a hyperparameter that can be adjusted to balance the contributions of the reward function and the expected utility. Based on the historical data, we can maximize the value of a state by:

\begin{equation}
    V^*(s_i) = \max_a Q^*(s_i,a)
\end{equation}

Finally, we can learn the optimal policy:

\begin{equation}
    V^*(s_i)= \max_a \sum_{{s'}_i} T(s_i,a,{s'}_i) [R_t (s_i,a,{s'}_i) + \gamma V^* ({s'}_i)]
\end{equation}

\begin{equation}
    \pi^*(s_i) = \arg\max_a \sum_{{s'}_i} T(s_i,a,{s'}_i)[R_t(s_i,a,{s'}_i) + \gamma V^* ({s'}_i)]
\end{equation}

A deep Q-network (DQN) with prioritized experience replay is adopted to solve the behavioral policy and predict the optimal trajectory that the agent would select for a given O-D pair. Consider a historical trajectory where the model has learned and estimated the $Q^*$ value for an action. The empirical trajectories should be encouraged to be sampled in the prediction. To prioritize the empirical trajectories, we first measure the difference between the predicted $Q^*$ value and the experienced $Q^*$ value in the same state for the same action, which is represented by $\delta_i$:

\begin{equation}
    \delta_i = R_t(s_i,a,{s'}_i) + \gamma V_{\theta^-}^* ({s'}_i)-Q_\theta^*(s_i,a)
\end{equation}

where $\theta^-$ represents the target deep neural network, and $\theta$ represents the current deep neural network. The equation can further be reformulated as:

\begin{equation}
    \delta_i = R_t(s_i,a,{s'}_i) + \gamma Q_{\theta^-}^* ({s'}_i,\arg\max_a Q_\theta^*({s'}_i,a)) - Q_\theta^*(s_i,a)
\end{equation}

The left-hand side of the equation, $R_t(s_i,a,{s'}_i) + \gamma Q_{\theta^-}^*({s'}_i,\arg\max_a Q_\theta^*({s'}_i,a)))$, is the target value, and the right-hand side of the equation, $Q_\theta^*(s_i,a)$, is the predicted value. This equation provides a quantitative measure of how much the deep neural network can learn from the given experience sample $i$. This is notated as the Double-Q temporal difference (TD) error.

To find an appropriate $\theta$ and an optimal trajectory, the TD error is minimized by using expectation of the samples from the replay buffer $D$:

\begin{equation}
    \min_\theta E_{(s_i,a,R,{s'}_i)} \sim D[(R_t(s_i,a,{s'}_i) + \gamma Q_{\theta^-}^* ({s'}_i, \arg\max_a Q_\theta^*({s'}_i,a)) - Q_\theta^*(s_i,a))^2]
\end{equation}

By doing so, we can predict the optimal trajectory for each pair of origin and destination.

The performance of trajectory prediction is measured based on the overlaps between the actual and the predicted trajectory. Since we created the grid world for the Markov decision process and optimal policy learning, the trajectories are characterized by the cells. The intersection between the set of cells in the predicted trajectory and the set of cells in the actual trajectory indicates the precision and recall of the model for O-D pairs \cite{Buckland1994}. Hence, we adopted the following formulas for measuring the precision and recall for each O-D pair to assess the performance of the model:

\begin{equation}
    Precision = \frac{True \; positive}{True \; positive + False \; positive}
\end{equation}

\begin{equation}
    Recall= \frac{True \; positive}{True \; positive + False \; negative}
\end{equation}

where $True \; positive$ is the number passing cells that the model correctly predicts and that presents in the actual trajectory; $False \; positive$ is the number of passing cells that the model incorrectly predicts and that are not in the actual trajectory; $False \; negative$ is the number of cells that are incorrectly predicted to be out of the trajectory. The calculation of the precision and recall allows us to quantitatively assess the performance of the model for each O-D pair and to identify the specific correctly predicted cells.

\subsection{Crisis scenario application with contextual factors}
The adaptability of the model is an important feature that enables the model to be widely used in various contexts, especially in crises situations, such as flooding, hurricanes, wildfires, and pandemics. The major change in the model due to the variation of the contexts is the destination and immediate reward when agents make their decisions on the actions. Therefore, in this section, we discuss the steps and strategies to update the destination selection and the reward table. The detailed steps are shown below:

\begin{itemize}[leftmargin=.3in]
    \item \textit{Step 1}: sampling the origins based on the historical data (or the density of population);
	\item \textit{Step 2}: predicting the destinations corresponding to each sampled origin using the well-trained first module of the model. In particular contexts, such as crises, some destinations might have been damaged or closed. Based on the observed situations, these destinations must be removed from results and only reachable destinations retained.
	\item \textit{Step 3}: generating the situation matrix $F$ to represent the situation between the origin and destination and updating the transformed reward table with the situation matrix $F$. Take the flooding event as an example. The value of each cell $f_{ij}$ in the situation matrix $F$ can represent the extent of flooding in that cell, such as the number of flooded road segments or the area of the flooded cell. Then, multiplied by a parameter $1/\beta$ to normalize the values, we can add the matrix to the transformed matrix or subtract the transformed matrix with the situation matrix, following the formula shown below:
	\begin{equation}
	    R_t = R' + L + D - F/\beta
	\end{equation}
	where, the value of $\beta$ can be selected based on the contribution of the situation to the trajectory selection of the agents.
    \item \textit{Step 4}: predicting the trajectory of the given origin and destination using the updated reward table and the priority deep Q-network.
    \item \textit{Step 5}: augmenting the prediction results. It is important to note that the model can simulate only the trajectory for each origin and destination pair, instead of directly predicting the number of vehicles on the road. The model, however, has the capability to simulate traffic conditions on the road based on the simulated trajectories. In this step, the algorithm can first simulate the trajectory for an origin and destination pair 50 times. Then the algorithm selects the cells with high pass-through probability. We use the number of times, $n_{ij}$, that the cell is included in the predicted trajectory as the weight of the cell, $v_{ij}$, and other cells are considered as zero, as the equation shown below:
    \begin{equation}
        v_{ij}^{(k)}=
        \left\{ \begin{array}{ll}
            n_{ij}^{(k)} & if \; grids \; are \; selected \\
            0 & otherwise
        \end{array} \right.
    \end{equation}
    Each pair of origin and destination would have one vehicle matrix, $V^{(k)}$, that encompasses the entire grid world. Since a cell would likely to be on multiple trajectories for different pairs of origins and destinations, for all pairs of origins and destinations, we use $c_{ij}=\sum_{k=1}^{N} v_{ij}^{(k)}$  to represent the number of vehicles passing through this cell, where $N$ is the number of simulated pairs of origins and destinations.
\end{itemize}

Through the implementation of all these steps, we can arrive at a final traffic matrix $C$ to represent the number of vehicles simulated in the cells on the grid world. This matrix would not only show the common trajectories that people would choose to reach their destinations, but also can inform the emergency response and resource allocation based on perturbations in mobility and traffic patterns in the context of crises. The capabilities of the model in simulation of urban mobility during crises is illustrated in the Section 4.4.

\section{Results}
To demonstrate the performance and capabilities of the proposed model, we first introduced the data sets and then elaborated a specific use case of the model with an application in flooding impact analysis in the Houston metropolitan area in the context of Hurricane Harvey in 2017.

\subsection{Data collection and preprocessing}
The data set, which is used to train, tune, and evaluate the adaptive reinforcement learning model, comes from INRIX, a private location intelligence company providing location-based data and analytics. The reason why we employ these data are twofold. First, the INRIX data provides very detailed coordinates along with the trajectory of the vehicles. That is, INRIX collected the coordinates of the vehicles every few seconds. Hence, the temporal and spatial resolution of the data is high. Second, the INRIX data are collected over the entire Houston metropolitan area over a time period of two months (March through April 2017) at all times of the day. This yielded a dataset of more than 26 million trip records collected in a continuous timeframe. Third, the INRIX data contains the trips for hundreds of providers (defined in Section 3.1) with four vehicle types. The vast volume of the trips and the diversity of vehicle types allow us to capture the various activity patterns and to enhance the robustness of the proposed approach in learning and simulating large-scale urban mobility.

To prepare this data for training, we selected providers with large record sizes and discarded roughly 20\% of the trips that were too short (distance between origin and destination less than 5 miles) that might induce noise into the training process. Then, we randomly selected 80\% of the trips as the training set and 20\% of the remaining data as the testing set.

\subsection{Destination prediction}
The first task is to predict the destinations, given the origins. We mainly train the model to automatically find an optimal value for the number of neighbors, $k$, which indicates the number of neighbors that should be included for estimating the coordinates of the destinations. Since we have a large number of providers with different types, we trained the model on the dataset from each provider separately.

\begin{figure}
    \centering
    \includegraphics[scale=0.5]{figures/Figure2.pdf}
    \caption{The performance of the model for destination prediction in four example providers. The average RMSEs for providers 1 through 4 are 0.0288, 0.0701, 0.0780, and 0.0824, respectively.}
    \label{fig:2}
\end{figure}

Figure \ref{fig:2} shows the prediction performance of the models in four example providers. As shown in the figure, the RMSE goes from a very small value, about 0.005, and then grows gradually to 0.1 when 80\% of the O-D pairs are predicted. Only 20\% of the O-D pairs have an RMSE greater than 0.1. The result indicates that the model can accurately predict the destinations for given origins based on the parameters learned from training data. Since we filtered out short-distance trips, the lengths of the remaining trips in the training and testing data are relatively long. By measuring the distances between the actual destinations and the predicted destinations, we find that a few predicted destinations are less than 0.3-mile from the actual destination, about 10\% of the predicted destinations are less than 0.6-mile from the actual destination, and more than 50\% of the predicted destinations are within the 3-mile distance from the actual destination. Considering the length of the trips themselves, the differences between the actual and predicted destinations are much shorter. Hence, the accuracy of the results is acceptable. Figure \ref{fig:3} further shows the spatial accuracy of some example O-D pairs. Figure \ref{fig:3} also shows the high performance of the model. Comparing the results among all providers, we also find that the model is quite stable with respect to the prediction performance and is not affected by the types of providers. This result demonstrates that the model is robust and generalizable to be applied to different data providers for destination prediction.

\begin{figure}
    \centering
    \includegraphics[scale=0.47]{figures/Figure3.pdf}
    \caption{Example origin-destination pairs showing the performance of the model. (The blue polygon is the borderline of Harris County, Texas, USA.)}
    \label{fig:3}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.5]{figures/Figure4.pdf}
    \caption{Example reward tables for origin and destination pairs, learned from training data using reinforcement learning.}
    \label{fig:4}
\end{figure}

\subsection{Trajectory prediction}
Once we acquire the destinations for given origins, we further use them to predict the optimal trajectories between the origins and destinations. The first step is to compute the reward table for each O-D pair. As discussed in the methodology section, we include only the cells within the bounding box. Figure \ref{fig:4} shows examples of reward tables for some of O-D pairs. The bounding box is a bit larger than the actual box tightly covering the origin and destination. In the majority of the example tables, such as R1, R2, and R5, we can find that the some significantly highlighted cells have high reward values, and the agent can effectively determine the optimal trajectories from the origin to the destination. Although historical data do not have a significant footprint for some O-D pairs, such as R3 and R4, the model can still impulse the agent to find an optimal trajectory.

A quantitative assessment of the model performance is shown in Figure \ref{fig:5}. We tested the model for multiple providers and O-D pairs. As the results suggest, the majority of the predicted trajectories are similar to the actual trajectories selected by the vehicles. Both average precision and recall are fairly high (although they are negatively affected by some extreme values). Figure \ref{fig:6} shows some examples of actual trajectories, predicted trajectories, and their overlaps. Despite the complexity of the historical trajectories transited by the vehicles, the model can still identify optimal trajectories for each pair of origin and destination. The performance of the model is also stable when it is trained and tested on different types of providers and datasets. These results and findings indicate that the proposed model is robust and could be used for trajectory-finding tasks.

\begin{figure}
    \centering
    \includegraphics[scale=0.42]{figures/Figure5.pdf}
    \caption{Performance of the trajectory prediction using module 2 of the proposed model (average precision: 0.765; average recall: 0.766).}
    \label{fig:5}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.41]{figures/Figure6.pdf}
    \caption{Example trajectories showing the performance of the model. (The blue polygon is the borderline of Harris County.)}
    \label{fig:6}
\end{figure}

\subsection{Simulation for flooding impact analysis}
To illustrate the application and capabilities of the proposed adaptive reinforcement model, we implemented the model for the situation during Hurricane Harvey in late August 2017 in the Houston metropolitan area. Hurricane Harvey, a Category 4 tropical storm, landed in Houston on August 26, 2017, and dissipated inland August 30, 2017. August 27 is the date with extreme and sustained rainfall that subsequently caused large-scale flooding over the urban areas, especially on urban road networks. As reported in news articles and government reports, more than 115,000 buildings were damaged, and 290 roads and highways were flooded \cite{Ibrahim2017,Sebastian2017}. To illustrate the performance of the model, we selected a specific short time interval, 10:00 a.m. through 12:00 p.m., August 29, 2017. We simulated vehicles that would move during this time interval and estimated their trajectory to understand the traffic conditions under the perturbation of urban flooding. Here we measured the traffic conditions in different locations by using the number of vehicles that move across a cell.

To be consistent, the data set used to validate the capability of the proposed model for scenario application was also obtained from INRIX. For each road segment, the dataset includes the average speed in 5-minute intervals and the speed in free-flow conditions (estimated as the speed limit for the roads). The flood situation on the road segments can also be reflected by the INRIX data. In the INRIX average speed dataset, the flooded road segments are identified by a designation of NULL for the traffic speed since no traffic records and no vehicles are driving through flooded road segments. To test the accuracy and validity of the data for identifying flooded roads, we also did a comparison of the NULL records between the datasets collected before and during the Hurricane Harvey, and also checked the records with flooded roads from government reports \cite{Dong2020}. We found that the NULL records presented only during the flooding period and were consistent with the flooded locations.

Using this dataset, we identified the flooded road segments and quantified the extent of flooding in a cell. In the grid map for Houston, Harris County, each cell contains multiple road segments. To measure the extent of flooding in a specific cell, we considered the number of flooded road segments within a cell as a metric. As shown in Figure \ref{fig:7}(a), the flooded cells are concentrated on the central area of Houston and along the major roads such as highways.

The traffic conditions in a cell can be captured by the actual traffic speed on the road segment. Since our model is mainly learned from the vehicles driving on major roads, here, we take into account only the road segments on which the speed limit is equal or greater than 50 km/hour so that all simulated results can be on the same road basis. This step also contributes to eliminating the effect of speed limits (types) of roads on the results. Hence, we computed the average of the actual speed on these selected road segments as the metric for traffic conditions for each cell. Figure \ref{fig:7}(b) shows areas where traffic was quite heavy due to flooding.

To understand the impact of flooding on the population movement patterns, especially in the trajectory selection, we implemented the proposed adaptive reinforcement learning model by adopting the flood conditions to adjust the parameters in the model. The cells with more flooded road segments would have higher negative values added to the reward table. Accordingly, the agent would be less likely to transit the severely flooded cell to reach their destination. The simulated results are shown in Figure \ref{fig:7}(c) and (d). From these results, we found that the majority of the road segments hit hardest by the flood had fewer vehicles compared to the road segments where flooding was less severe (i.e., passable roads). This result indicates that, the model effectively simulated the trajectories selected by people who deliberately avoided flooded road segments. Hence, the road segments that were less flooded tend to have a great number of vehicles with low speed, or fewer vehicles with high speed. The results are quite realistic, which reveals the redistribution of vehicles and the utility of roads in urban road networks when flooding disrupted this area.

\begin{figure}
    \centering
    \includegraphics[scale=0.5]{figures/Figure7.pdf}
    \caption{The simulation results. (a) Number of flooded segments in Houston, Harris County (the blue borderline); (b) Actual average speed for major road segments during the investigated time period; (c) Predicted density of vehicles during the investigated time period; and (d) Relationships among three variables for simulation validation.}
    \label{fig:7}
\end{figure}

\section{Discussion and Concluding Remarks}
In this paper, we presented an adaptive reinforcement learning model that can simulate human mobility under crisis conditions based on the mobility patterns learned from historical data collected under normal circumstances. The proposed model can overcome the scarcity of mobility pattern data in crisis contexts, and also has the capability to simulate different crisis situations for impact analysis. The study showed the application of the proposed model using the digital trace data collected in Houston during regular situations and demonstrated the adaptability of the model in the case of urban flooding during Hurricane Harvey in 2017.

The presented model has multiple contributions. First, the proposed simulation model and its adaptability to different crisis contexts can contribute to hazard mitigation and resilience improvement of urban areas. Predictive emergency warning and situational awareness are key components of hazard mitigation, which would benefit from effective simulation of crisis situations. Simulating urban mobility during crises is an essential step to capture the crisis situation and further inform early warning and mitigation planning. For example, using the proposed simulation model, city planners, emergency managers, and decision-makers can simulate scenarios of flooding (such as 100-year and 500-year flood scenarios) and road inundations to enable examination of the effects of movement trajectories and traffic congestion. This output information could inform about areas that would experience significant mobility disruptions and traffic congestion under varying flood scenarios. With sufficient examination of such contingencies, first responders and decision-makers can effectively assess the extent to which reducing the flood vulnerability of a major road or highway would decrease traffic congestion and improve movement trajectories. Also, the proposed model could be used to simulate mobility disruptions between critical origins (such as areas where vulnerable populations, such as the elderly and low-income households, are located) and critical destinations (such as healthcare facilities). In addition, the outcomes of the proposed model can help with near-real-time prediction of mobility disruptions and traffic congestions based on the real-time information, such as the inundated roads, to improve the accessibility of critical facilities and efficiency of responder dispatch.

Second, the presented study also offers a new perspective for using deep learning techniques and digital trace data for urban mobility analysis. Although extremely big data related to human digital traces has been generated by daily mobility activity of individuals, existing deep learning techniques are limited to the prediction of regular mobility patterns. The proposed model enables consideration of contextual factors, such as road conditions, road network connectivity and population distribution, in the simulation process, which would significantly extend the capabilities of existing deep learning models gleaned from normal situations and applied to unfamiliar conditions, such as emergencies and traffic anomalies, to support decision-making and response strategies. For example, one can learn the mobility demand of the road networks from daily mobility data and predict the utility of each road segment by changing the layout of the networks. This can help design efficient transportation networks and traffic control systems to improve the performance of urban road networks.

Finally, although the proposed model is general and can be adaptive to various phenomena and applications, future research could address some limitations. First, although this study achieves good performance for predicting the destinations and trajectories at a relatively fine-grained scale, the prediction of local movements, such as visits to nearby neighborhoods and points of interest, still needs to be improved. Local movements provide essential information about the lifestyle of residents. This study predicts the pure coordinates of destination and trajectories regardless of the specific locations visited and the roads traversed. Inclusion of this information in creating the movement profile of populations and predicting local movements would be of interest and importance in applications related to city and emergency planning. Second, the model has a great potential to be used as an application for rapid prediction. Implementing the proposed algorithm is still time-consuming, which would take a great amount of computational cost for rapid prediction. Hence, there is a need of improving the efficiency of the algorithm so that it can be used for more real-time prediction and across more cities and crises contexts.

\section*{Acknowledgement}
This material is based in part upon work supported by the National Science Foundation under Grant Number CMMI-1846069 (CAREER), CMMI-1832662 (CRISP 2.0 Type 2), the National Academies’ Gulf Research Program Early-Career Research Fellowship, the Amazon Web Services (AWS) Machine Learning Award. The authors also would like to thank INRIX for providing the data for this research. Jan Gerston provided editorial services. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation and Amazon Web Services.

\section*{Data availability}
The data that support the findings of this study are available from INRIX and TTI (Texas A\&M Transportation Institute), but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available.

\section*{Competing interests}
The authors declare that they have no competing interests.

\section*{Author contribution}
C.F., X.J., and A.M. designed the study. X.J. and C.F. implemented the method and empirical case study. A.M. support with data acquisition. C.F. and A.M. wrote the main manuscript. All authors reviewed the manuscript.

\bibliographystyle{unsrt}


\section{ABSTRACT}
\begin{abstract}
We report emergence of photoluminescence at room temperature in trigonal $Cs_3Bi_2Br_9$ at high pressures. Enhancement in intensity with pressure is found to be driven by increase in distortion of $BiBr_6$ octahedra and iso-structural transitions. Electronic band structure calculations show the sample in the high pressure phase to be an indirect band gap semiconductor. The luminescence peak profile show signatures related to the recombination of free and self trapped excitons, respectively. Blue shift of the both peaks till about 4.4 GPa are due to the exciton recombination before relaxation due to the decrease in exciton lifetime with scattering from phonons
\end{abstract}
\maketitle

In recent years, halide perovskites of the formula $A_3B_2X_9$ ($A$:Cs, Rb; $B$:Sb or Bi; and $X$:halogen) have been in research forefront due to their potential applications in optoelectronic devices and solar cells\cite{Anna,Fan,Bin}. Since toxicity of lead has been a consistent concern for the environment and health safety, lead free halide perovskites have a great importance. Among several important materials under study, $Cs_3Bi_2Br_9$ has emerged as a model system for  detailed investigation as it has been found to be stable under atmospheric conditions. At ambient conditions $Cs_3Bi_2Br_9$  has trigonal crystal structure with space group $P\bar{3}m1$ and has been found to be ferroelectric\cite{Lazarini}. $Cs_3Bi_2Br_9$ is found to undergo a structural transition to monoclinic phase at 95 K \cite{Bator,IP}. Early studies on $Cs_3Bi_2Br_9$ single crystals by Timmermans and Blasse \cite{cwm} showed blue luminescence at 4K, which got quenched by 20K. A recent careful measurements by Bass et al.\cite{bass} on high purity nanocrystalline $Cs_3Bi_2Br_9$ show a weak but structured photoluminescence (PL) at room temperature and their feature is explained by excitons coupling to vibrons in the lattice.  All the above studies indicate large exciton binding energy leading to radiative recombination of excitons.

Pressure, being a unique thermodynamic variable, may dramatically alter the crystal structure, and hence optical properties of materials. High pressure (HP) studies on $CsPbBr_3$ nanocrystals have led to new structures with enhanced luminescence properties \cite{Yasutaka}. Ma et al.\cite{Zhiwei} have shown that by pressure tuning, initially non-luminescence $Cs_4PbBr_6$ nanocrystals exhibit a blue luminescence above 3 GPa  due to a structural transition.
The unit cell of $Cs_3Bi_2Br_9$ has a lot of empty space as one third of octahedral positions are vacant. Therefore, pressure is expected have a significant effect on  its structural behaviour. Recently Li et al.\cite{Qian} have reported pressure induced structural transition of 2D hybrid halide perovskite $(CH_3NH_3)_3Bi_2Br_9 $. In this work we have carried out a detailed pressure dependent study on structural, and optical  properties of $Cs_3Bi_2Br_9$. We observe an anomalous increase in PL intensity at about 2.9 GPa mediated by distortion and elongation of $BiBr_6$ octahedra due to iso-structural transition induced by high compressibility of the sample.

The sample is synthesized in the laboratory using sol-gel technique (details are given in Supplementary Information (SI)). Characterization of synthesized sample using x-ray diffraction (XRD) revealed trigonal structure ($P\bar{3}m1$) with lattice parameters: a=7.9689(2)$\AA$, c=9.8577(6)$\AA$, and are in good agreement with previous report \cite{Lazarini}.
We have recorded HP PL spectra of $Cs_3Bi_2Br_9$ using excitation wavelength of 488nm and shown in Fig.1. $Cs_3Bi_2Br_9$ does not show any observable PL at ambient condition.
Interestingly with increase in pressure, a luminescence appears at about 1.4 GPa as shown in the bottom of Fig.1(a). The intensity is too low compared to the ambient background. The PL intensity shows a maximum at about 2.9 GPa and then almost disappears above 11 GPa. The PL peak profiles show a long tail at lower energies with a kink at higher energy side(SI Fig.S2).
In Fig.1(b), we have shown the pressure evolution of peak positions of the two peaks. Both the peaks show a blue-shift till about 4.4 GPa followed by a redshift. Fig.1(c) shows the relative change in integrated intensities of both the peaks. The higher energy Peak-A shows much larger change, though its intensity is much lower compared to Peak-B. We shall discuss the results in later sections in detail after we discuss the structural evolution with pressure.

In Fig.2(a) we have shown the pressure evolution of Raman spectra (excitation $\lambda$ = 532 nm) of $Cs_3Bi_2Br_9$ up to about 16 GPa.
At ambient condition we observe six Raman modes (SI Fig.S3). The modes at $167cm^{-1}$ and $192cm^{-1}$ are indexed to $E_g$ and  $A_{1g}$, respectively and are associated with stretching vibrations of $Bi-Br$ bonds \cite{Bator}. Other four are labeled as $M_1$($63.9cm^{-1}$), $M_2$($67.3cm^{-1}$), $M_3$($78.8cm^{-1}$), $M_4$($92.2cm^{-1}$) and originate from deformation of $BiBr_6$ octahedra due to vibrations of Br atoms only\cite{ Bator, MYa}.
With pressure all modes are blue shifted. The first four low frequency modes broaden with pressure with reordering of their respective intensities.
At 4.7 GPa $E_g$, $A_{1g}$ modes split into two modes each (labeled as $M_5$, $M_6$, $M_7$, $M_8$) indicating a structural transition with lowering in symmetry. Gradual extensive broadening of all Raman modes with pressure indicate to amorphization of the sample.
Pressure evolution of full width half maximum (FWHM) of $A_{1g}$, $E_g$ are shown in Fig.2(b). Increase in the life time of both the phonon modes is evident from sharp decrease of FWHM till 2.4 GPa. Similarly in Fig.2(c) one can see an increase in the intensity of symmetric stretching mode $A_{1g}$ with respect to that of asymmetric stretching mode $E_g$.  Enhancement in a Raman mode intensity represents a strong coupling between the excited electrons and particular vibrational mode of the lattice as well as enhancement in the dipole-dipole polarizability of the lattice. The above anomalous behaviour of the Raman modes and enhancement of PL do indicate a strong exciton-phonon coupling.

We have carried out HP XRD measurements up to about 25 GPa on powder $Cs_3Bi_2Br_9$ at room temperature.
Pressure evolution of XRD patterns show appearance and disappearance of certain Bragg peaks (Fig.3(a)). We assign the ambient phase as Trigonal(I). At about 1.3 GPa a new peak appears at $2\theta$ = 8.72$^{\circ}$. Indexing results in a trigonal structure with lattice parameters: $a$=7.7350(2), $c$=28.934(1) $\AA$; space group $P\bar{3}m1$ (Trigonal(II)). In Trigonal(II) phase the $c$-axis is almost three times of that  of  Trigonal(I) phase.  With further increase in pressure to 3.1 GPa, another new Bragg peak appears at $2\theta$ = 9.93$^{\circ}$. The XRD pattern is indexed again to trigonal structure with lattice parameters: $a$=15.0956(4), $c$=9.4438(4)$\AA$; space group $P\bar{3}m1$ (Trigonal(III)) with $a$ and $b$-axes being double with respect to that of Trigonal(I) phase. Several new peaks emerge at 6 GPa indicating lowering in symmetry and the XRD pattern is indexed to a triclinic phase with structural parameters: $a$=9.455(1), $b$=7.4426(9), c=14.168(2) $\AA$, $\alpha$=59.160(7)$^{\circ}$, $\beta$=105.635(7)$^{\circ}$, $\gamma$=98.549(7)$^{\circ}$; space group $P\bar{1}$. Therefore the XRD data show two iso-structural transitions followed by a change in symmetry. The broadening of Bragg peaks and their disappearance above 11 GPa corroborate the pressure induced amorphization of the sample. Below 10 GPa, the sharp Bragg peaks of the Ag-pressure marker show the absence of non-hydrostatic stress on the sample. $3^{rd}$ order Birch-Murnaghan EOS is fitted to the volume {\it vs.} pressure data by normalizing the volume of all the trigonal phases with respect to the ambient phase (SI Fig.S6). The EOS fit results in low bulk moduli ($B_0$) values: $12\pm 1$ and $61.5\pm0.6$ GPa for trigonal and  triclinic phases, respectively.
Rietveld refinement of all trigonal phases are carried out starting with the model of the Trigonal(I) phase \cite{Lazarini}. The refined atom positions and XRD patterns are documented in the SI (Table-SI, SII, SIII and Fig.S7).
In Fig.3(b) we have shown the octahedral arrangement in the unit cell for all trigonal phases. Trigonal(II) and Trigonal(III) phase consists of  three and two types $BiBr_6$ octahedra, respectively. The large compressibility results in a increase of distortion index (DI) by about 67\% and quadratic elongation (QE) by about 1\% of $BiBr_6$ octahedra in the pressure range 0-0.9 GPa itself( SI Table-SIV).
It can be seen that Bi-atoms have moved significantly away from the central positions of the octahedra resulting in different $Bi-Br$ bond lengths and bond angles giving rise to significant rotation and distortion of the octahedra. The increase in distortion in $BiBr_6$ octahedra leads to the anomalous behaviour of Raman modes.

Using first principles DFT calculations we have computed the PBE+SOC band structure of $Cs_3Bi_2Br_9$ at  0.1 and 3.1GPa. Electronic band structure and partial density of states (PDOS) predicted by our calculations for the 0.1 GPa structure agree very well with those calculated using hybrid functionals (HSE06+SOC) for the same system at ambient  pressure\cite{bass} (SI Fig.S8).
It is to be noted that PBE is known to under-estimate the band-gap and a more accurate calculation involving computationally far more expensive hybrid functional is expected to give a better estimate of the gap.  Nonetheless, our PBE+SOC calculations are expected to provide a good qualitative picture of the nature of the band-structure of the material close to the Fermi level.
Fig.4 shows the PBE+SOC band structure and PDOS of this material at 3.1 GPa, calculated using our experimentally determined crystal structure. The band structure shows that the system is an indirect band gap semiconductor, with none of the low lying excitations occur at the same $k$-point. We find that the top of the valence band is dominated by $Br-p$ and $Bi-s$, whereas the bottom of the conduction band is dominated by $Bi-p$ indicating the importance of Bi-centers in the $BiBr_6$ octahedra for photoinduced activity. Also the electronic  band gap(3.1 GPa) is found to get reduced with increase in pressure.

The anomalous changes observed in Raman spectroscopy and XRD studies do indicate that application of pressure leads to several structural changes induced by the deformation of $BiBr_6$ polyhedra  affecting the electronic properties. The emergence of PL can be related to the phase transition  with its appearance at Trigonal(I) - Trigonal(II) transition. The same can be corroborated by  the fact that  the PL intensity reaches its maximum at about 2.9 GPa (Trigonal-III phase) and shows a $2^{nd}$ maximum at about 4.4 GPa(Triclinic phase) (Fig.1(b)). The PL peak feature shows a structure with two peaks suggesting coexistence of free exciton(FE) and self-trapped exciton (STE) recombination\cite{Shunran, Jiangjian, Ivan}. We have carried out an analysis of the obtained PL intensity ($I$) by changing the excitation laser power ($L$) using the power law behaviour $I\sim L^K$ law\cite{Schmidt} ({\bf SI Fig.S9}).The value of $K$ from the fit is found to be 1.17 for Peak-A and 1.04 for Peak-B. Schimdt and Lischka \cite{Schmidt} showed that for free and bound excitons the value of $K$ should be: $1<K<2$. This confirms the excitonic nature of pressure induced luminescence in $Cs_3Bi_2Br_9$. Peak-A with its narrow feature can be the free exciton. Iso-structural phase transition driven by structural distortion of $BiBr_6$ octahedra can give rise to photoinduced pseudo Jahn-Teller distortion of the individual $BiBr_6$ octahedra and  lead to the formation and radiative recombination of STE. The same is reflected in the asymmetric nature of Peak-B profile with a long tail in low-energy side. The initial blue shift of the PL peaks till about 4.4 GPa also can be related to the extensive structural distortion before transition to lowest symmetry triclinic phase. Due to the indirect band gap of the sample, the life time of the excitons shortens because of scattering with one or more phonon and recombination of exciton happens much before exciton relaxation. The red-shift of the PL peaks above 4.2 GPa is probably due to decrease of band gap \cite{Takuyao}. This shows that the structural distortion is a important factor to be considered for lead-free halide perovskites for exhibiting luminescence. We believe that this work will lead to further studies on Pb-free stable inorganic halide perovskites.

\section{Conclusions}
We have carried out a detailed high pressure study on $Cs_3Bi_2Br_9$ perovskite by PL, XRD and Raman spectroscopy measurements. HP XRD and Raman measurements show two isostructural transitions in ambient trigonal phase at about 0.9 GPa and 2.4 GPa, respectively. These isostructural transitions take place with extensive structural distortion due to large compressibility of the system. The sample shows appearance of PL at about 1.4 GPa driven by the isostructural transition. PL feature shows a kink at about 2.1 eV due to free exciton recombination and another broad peak at about 1.8 eV due to self trapped exciton (STE) recombination. Electronic density of states and band structure calculation show  that the valence band and conduction band consist of hybridization of $Br-p$, $Bi-s$ orbitals and $Bi-p$ orbitals. The enhancement of PL intensity and blue shift in PL peaks are related to the large distortion and quadratic elongation of $BiBr_6$ octahedra.

\section{Acknowledgments}
We acknowledge the financial support from the Ministry of Earth Sciences, Government of India, grant number MoES/16/25/10-RDEAS. The financial support from Department of Science and Technology, Government of India under Indo-Italy Executive Programme of Scientific and Technological Cooperation is gratefully acknowledged. DS, PS and BG also acknowledge the fellowship grant supported by the INSPIRE program, Department of Science and Technology, Government of India.

\section{Sample preparation}

Crystalline powders of $Cs_3Bi_2Br_9$ are prepared by dissolving 63.84mg of $CsBr$ and 89.74mg of $BiBr_3(\geq99\%)$ with $3:2$ molar ratio in 2 mL of $HBr(99.999\%)$ (48 wt. $\%$ in water $\geq$ 99.99 $\%$) in a capped vial and then the mixture are heated to $110^{o}$C for 2 hours. After that, the mixture are kept for an additional 2 hours at the same temperature without stirring. The obtained solution is allowed to cool to room temperature overnight, followed by washing with ethanol and drying.

\section{Experimental methods}
High pressure (HP) experiments at room temperature are performed using piston cylinder type diamond anvil cell (DAC) of culet size 300$ \mu$m. A small amount of powder sample along with pressure marker are loaded inside the central hole of diameter 100$\mu$m of steel gasket preindented to a thickness 60$\mu$m. For high pressure XRD and Raman measurements 4:1 methanol ethanol mixture are used as a pressure transmitting medium (PTM), whereas for PL measurements silicon oil is used. PL and Raman spectroscopy measurements are carried out using a confocal micro-Raman spectrometer(Monovista from SI GmBH) and Horiba Jobin-Yvon LabRAM HR-800 spectrometer with back scattering geometry respectively.
A micron sized ruby is loaded along with the sample in the DAC for determination of pressure using ruby fluorescence technique\cite{Mao}.
488nm and 532nm lasers are used to excite the sample for PL and Raman measurements, respectively.
Infinity corrected 20X objective having long working distance is used to focus the laser and collect the radiation.
Pressure dependent XRD measurements are carried out at  XPRESS beamline in the Elettra synchrotron radiation source using a monochromatic x-ray of wavelength 0.4957$\AA$ collimated to 20 $mu$m at the sample.  We have used $Ag$ as pressure marker and pressure is calculated using  equation of state (EOS) of $Ag$ \cite{Dewaele}.
FIT2D software\cite{Hammersley} is employed to convert 2D diffraction image to intensity versus 2$\theta$ plot.
All the XRD data are analyzed by CRYSFIRE\cite{Shirley} and GSAS \cite{Brian} program.
We have used VESTA\cite{Koichi} software for visualization of crystal structure. UV-absorption spectra are recorded using Jasco V-670 spectrophotometer.

\section{THEORETICAL METHODS}

First principles calculations are performed within the framework of density functional theory(DFT)\cite{Hohenberg, Sham} using projected augmented wave(PAW) \cite{PE, Kresse} method as implemented in the plane-wave based VASP code \cite{Hafner, Furth, JFurth}. The exchange-correlation functional is chosen to be the Perdew-Burke-Ernzerhof(PBE)\cite{John} implementation of the generalized gradient approximation (GGA). Spin-orbit coupling(SOC), which is crucial for understanding the electronic structure of the Bi-containing material has been included self-consistently in our calculation. An energy cut-off of 450 eV is used for plane wave expansions. A k-mesh of $9\times9\times9$ is used for the 0.1 GPa structure that contains 14 atoms. Whereas, a k-mesh of $3\times3\times5$ is used for the 3.1GPa crystal structure containing 56 atoms. The lattice parameters have been fixed at experimental values and the ionic positions have been relaxed using conjugate gradient algorithm, until the Helmann-Feynmann forces became less than 0.005 eV/$\AA$. The energy convergence with respect to all computational parameters have been carefully examined.

The upper panel of the Supplementary Fig.S8 shows the computed PBE+SOC band structure at 0.1GPa.  $Cs_3Bi_2Br_9$ has both direct(2.52eV) and indirect (2.4eV) band gap. Direct energy gap is at G whereas indirect energy gap takes place from G to A. It is to be noted that PBE is known to under-estimate the band-gap and a more accurate calculation involving computationally far more expensive hybrid functional is expected to give a better estimate of the gap. Nonetheless, our PBE+SOC calculations are expected to provide a good qualitative picture of the nature of the band-structure of the material close to the Fermi level. This is ascertained by the fact that the nature of the transitions as predicted by our calculations for the 0.1 GPa structure agree very well with those calculated using hybrid functionals (HSE06+SOC) for the same system at ambient  pressure\cite{bass}.  The calculated partial density of states (lower panel of Supplementary Fig.S8) shows that the top of the valence band is dominated by $Br-p$ and $Bi-s$ hybridized orbitals, whereas the bottom of the conduction band is dominated by $Bi-p$.


\section{Introduction}
Recently, the  LHCb measurements has confirmed the $X(4140)$ state with high statistic data~\cite{Aaij:2016iza,Aaij:2016nsc},  with a mass $4146.5\pm 4.5^{+4.6}_{-2.8}$~MeV and a width $ 83\pm 21^{+21}_{-14}$~MeV, much larger than the previous experimental measurements~\cite{PDG2018}, and the quantum numbers were determined to be $J^{PC}=1^{++}$.
There are many different suggestions about the structure of the $X(4140)$~\cite{Chen:2016iua,Wang:2018qpe}, especially because of the large discrepancy of the width.

Indeed, it is natural and necessary to exhaust the possible $q\bar{q}$ description of the observed states before restoring to the more exotic assignments.
 In this work, taking the meson wave functions by solving the relativistic/non-relativistic Schr\"{o}dinger equation, we investigate the decay properties of the $X(4140)$ as the assignment of charmonium state in the $^3P_0$ model, and provide more information about the decay modes for searching for the $X(4140)$, which could be useful to extract the more precision width.

\section{The $^3P_0$ decay model}
\label{sec:model}
The $^3P_0$ model, also known as the quark-pair creation model, has been widely applied to study strong decays of hadrons with considerable success~\cite{Micu:1968mk,Close:2005se,Li:2009rka,Lu:2016bbk}.
In the model, the transition operator $T$ of the decay $A\rightarrow BC$ in the $^3P_0$ model can be written,
\begin{eqnarray}
  T &=&-3\gamma\sum_m\langle 1m1-m|00\rangle\int
    d^3\boldsymbol{p}_3d^3\boldsymbol{p}_4\delta^3(\boldsymbol{p}_3+\boldsymbol{p}_4)\nonumber\\
&& \times   {\cal{Y}}^m_1\left(\frac{\boldsymbol{p}_3-\boldsymbol{p}_4}{2}\right
    )\chi^{34}_{1,-m}\phi^{34}_0\omega^{34}_0b^\dagger_3(\boldsymbol{p}_3)d^\dagger_4(\boldsymbol{p}_4),
\end{eqnarray}
where $\gamma$ is a dimensionless parameter corresponding the strength of quark-antiquark $q_3\bar{q}_4$ pair producing from the vacuum.
The decay width
$\Gamma(A\rightarrow BC)$ can be expressed in terms of the partial wave amplitude,
\begin{eqnarray}
\Gamma(A\rightarrow BC)= \frac{\pi
|\boldsymbol{P}|}{4M^2_A}\sum_{LS}|{\cal{M}}^{LS}(\boldsymbol{P})|^2, \label{width1}
\end{eqnarray}
where the explicit expressions for ${\cal{M}}^{LS}(\boldsymbol{P})$ can be found in Refs.~\citenum{Close:2005se,Li:2009rka,Lu:2016bbk}.

In this work, we  take into account two choices of the wave functions for mesons: Case A: the modified non-relativistic quark model (MNRQM), Case B: the modified  Godfrey-Isgur model (MGI). The details of model and the wave functions can be found in Ref.~\citenum{Hao:2019fjg}.

\section{Results and discussions}
\label{sec:results}

The mass spectra of the charmonium states predicted by two kinds of quark models are shown in Ref.~\citenum{Hao:2019fjg}. According to the PDG~\cite{PDG2018}, the mass of the $X(4140)$ [$I^G(J^{PC})=0^+(1^{++})$] is $4146.8\pm 2.4$~MeV, which is consistent with the predicted mass of $\chi_{c1}(3P)$ in both models within the uncertainties of the models. Next, we will calculate the strong decay widths of the $X(4140)$ state as the $\chi_{c1}(3P)$  assignment.

In the $^3P_0$ model, the strength of quark-antiquark pair creating from the vacuum, is taken to be $\gamma=4.52$ in Case A,  and $\gamma=5.90$ for Case B, by fitting to the total widths of the well established charmonium states,  $\psi(3770)$  ($1^3D_1$), $\psi(4040)$ ($3^3S_1$), $\psi(4160)$ ($2^3D_1$), and $\chi_{c2}(2P)$.
\begin{table
\begin{center}\label{tab:width}
\tbl{Decay widths of the ${\chi}_{c0}(3P)$, ${\chi}_{c1}(3P)$ and ${\chi}_{c2}(3P)$ states (in MeV). The mass of the ${\chi}_{c1}(3P)$ is taken to be the one of the $X(4140)$, and the masses of the ${\chi}_{c0}(3P)$,  and ${\chi}_{c2}(3P)$ are taken from Ref.~\citenum{Hao:2019fjg}, respectively for Case A and Case B.}
{\begin{tabular}{c|c|c|c|c}
\hline\hline
 State                &Channel                        &Mode                 &$\Gamma$ (Case A)  &$\Gamma$ (Case B)\\
\hline
${\chi}_{c0}(3P)$     &$0^+\rightarrow 0^-0^-$           &$D\bar{D}$           &10.58          &0.22    \\
    $$                &$$                             &$D^+_sD^-_s$           &0.37            &1.87    \\
    $$                &$0^+\rightarrow 1^-1^-$            &$D^*\bar{D}^*$       &16.28           &35.95  \\
Total Width           &$$                             &$$                    &27.23          &38.03  \\
\hline
${\chi}_{c1}(3P)$     &$1^+\rightarrow 0^-1^-$          &$D\bar{D}^*$         &4.54           &14.48     \\
    $$                &$$                             &$D_s\bar{D}_s^*$     &1.23            &0.70    \\
    $$                &$1^+\rightarrow 1^-1^-$           &$D^*\bar{D}^*$       &6.86           &16.17   \\
Total Width           &$$                             &$$                    &12.63          &31.34    \\
\hline
${\chi}_{c2}(3P)$     &$2^+\rightarrow 0^-0^-$           &$D\bar{D}$           &7.71            &8.79      \\
    $$                &$$                             &$D^+_sD^-_s$           &0.63            &0.10 \\
    $$                &$2^+\rightarrow 0^-1^-$           &$D\bar{D}^*$         &20.04            &11.34      \\
    $$                &$$                             &$D_s\bar{D}_s^*$     &0.17            &0.13     \\
    $$                &$2^+\rightarrow 1^-1^-$           &$D^*\bar{D}^*$       &11.33            &26.87    \\
Total Width          &$$                             &$$                    &39.89         &47.23    \\
\hline\hline
\end{tabular}}
\end{center}
\end{table}

\begin{figure}[htpb]
  \centering

  \includegraphics[scale=0.33]{c0MNR.eps}
  \includegraphics[scale=0.33]{c1MNR.eps}
  \includegraphics[scale=0.33]{c2MNR.eps}\\
  \includegraphics[scale=0.33]{c0MGI.eps}
\includegraphics[scale=0.33]{c1MGI.eps}
\includegraphics[scale=0.33]{c2MGI.eps}
  \caption{The dependence of the width of ${\chi}_{c0}(3P)$, ${\chi}_{c1}(3P)$ and ${\chi}_{c2}(3P)$ on the initial state mass with the wave functions of Case A (up panels) and Case B (low panels).}
  \label{fig:widthNRQM}
\end{figure}

With the above parameters, we have calculated the partial decay widths and total decay width, as shown in Table~\ref{tab:width} for both Case A and Case B. The total widths are 12.63~MeV for Case A, and $31.34$~MeV for Case B, which are consistent with the average value $\Gamma=22^{+8}_{-7}$~MeV of PDG, considering the uncertainties of the $^3P_0$ model.  It should be pointed out that the decay modes $D\bar{D}^*$ and $D^*\bar{D}^*$ have large decay widths. We suggest to search for this state in those two channels, and to measure the width precisely, which can be shed light on its nature.

We also show the dependence of the $\chi_{c1}(3P)$ decay width on the masses in Fig.~\ref{fig:widthNRQM}. Indeed, if the small width of the $X(4140)$ is confirmed in future high-precision measurements, the $X(4140)$ could be explained as the charmonium state $\chi_{c1}(3P)$. In addition, the processes of the $B^+\to  J/\psi \phi K^+$  and $e^+e^- \to \gamma J/\psi \phi K^+$ decay were investigated in Refs.~\citenum{Wang:2017mrt,Wang:2018djr} ,
where the existences of the $X(4140)$, with the small width, and the molecular state $X(4160)$ is supported by the low $J/\psi \phi$  invariant mass distributions ~\cite{Aaij:2016iza,Aaij:2016nsc,Ablikim:2017cbv}.

Finally, the strong decay width of the $\chi_{c0}(3P)$ and $\chi_{c2}(3P)$ states in Table~\ref{tab:width}, and the mass dependence of the total widths are also shown in Fig.~\ref{fig:widthNRQM}.

\section{SUMMARY}
\label{sec:summary}
We have investigated the strong decay properties of the $X(4140)$ with the assignment of the $\chi_{c1}(3P)$ states in the $^3P_0$ model, taking two cases for the wave functions.
The total decay width of the $\chi_{c1}(3P)$ is predicted to be $12.63$~MeV for Case A, and $31.34$~MeV for Case B, both of which are in agreement with the PDG average width of $X(4140)$. Thus, we conclude that, the $X(4140)$, with a small width, could be explained as the charmonium state $\chi_{c1}(3P)$, and the high-precision measurement about the $X(4140)$ could shed light on the nature of the $X(4140)$.

\section{Acknowledgements}

This work is partly supported by the National Natural Science Foundation of China under Grant Nos. 11505158, 11605158, the Key Research Projects of Henan Higher Education Institutions (No. 20A140027), and the Academic Improvement Project of Zhengzhou University.

\bibliographystyle{ws-procs9x6}