# Causal Analysis and Prediction of Human Mobility in the U.S. during the
COVID-19 Pandemic

Subhrajit Sinha<EMAIL_ADDRESS>Meghna Chakraborty<EMAIL_ADDRESS>
###### Abstract

Since the increasing outspread of COVID-19 in the U.S., with the highest
number of confirmed cases and deaths in the world as of September 2020, most
states in the country have enforced travel restrictions resulting in sharp
reductions in mobility. However, the overall impact and long-term implications
of this crisis to travel and mobility remain uncertain. To this end, this
study develops an analytical framework that determines and analyzes the most
dominant factors impacting human mobility and travel in the U.S. during this
pandemic. In particular, the study uses Granger causality to determine the
important predictors influencing daily vehicle miles traveled and utilize
linear regularization algorithms, including Ridge and LASSO techniques, to
model and predict mobility. State-level time-series data were obtained from
various open-access sources for the period starting from March 1, 2020 through
June 13, 2020 and the entire data set was divided into two parts for training
and testing purposes. The variables selected by Granger causality were used to
train the three different reduced order models by ordinary least square
regression, Ridge regression, and LASSO regression algorithms. Finally, the
prediction accuracy of the developed models was examined on the test data. The
results indicate that the factors including the number of new COVID cases,
social distancing index, population staying at home, percent of out of county
trips, trips to different destinations, socioeconomic status, percent of
people working from home, and statewide closure, among others, were the most
important factors influencing daily VMT. Also, among all the modeling
techniques, Ridge regression provides the most superior performance with the
least error, while LASSO regression also performed better than the ordinary
least square model.

###### keywords:

COVID-19 Pandemic , Time-series Data , Time-series Analysis , Causality ,
Mobility , Vehicle Miles Traveled , Granger Causality , Linear Regression ,
Regularization Algorithms , LASSO , Ridge Regression

††journal: Journal of LaTeX Templatesfn1fn1footnotetext: Pacific Northwest
National Laboratory. 902 Battelle Blvd, Richland, Washington,
99354fn2fn2footnotetext: Michigan State University. 428 S Shaw Ln. East
Lansing, Michigan, 48910

##  1 Introduction and Background

The novel coronavirus (COVID-19) pandemic is delineating our times’ global
health crisis and has had a significant impact on the way we understand our
everyday lives. Since its emergence in Asia in late 2019, all continents in
the world except Antarctica have been fighting the virus in earnest. The World
Health Organization (WHO) [1] has declared COVID-19 a global pandemic on March
11, 2020 and the United States declared a national emergency on March 13, 2020
[2]. As of September 15, 2020, almost 30 million cases of COVID-19 were
confirmed in 215 countries around the globe [3] and among all countries, the
U.S. has the highest number of confirmed cases and fatalities in the world due
to COVID-19 [3]. Several countries have closed their borders, exercised
lockdowns, curfews, stay-at-home orders, and social distancing protocols,
resulting in sharp reductions in mobility and travel demand at local,
regional, national, and even international levels. By March 24, 2020, more
than 20 percent of the world’s population has been ordered to remain at home
as governments, health, and administrative organizations take extreme measures
to protect their communities from the spread of the virus [4].

In the U.S., more than 40 states have already enforced stay-at-home order, the
earliest being in California effective from March 19, 2020 [5]. Only in the
U.S., the mobility restrictions during the COVID-19 pandemic, in the form of
travel bans, stay-at-home mandates, and lockdown policies, have impacted
millions of people. Overall, human mobility has been severely impacted due to
the travel restrictions and individual concerns to avoid public gatherings,
resulting in tremendous economic impacts in transportation sectors. However,
the overall influence and the long-term implications of this pandemic to
mobility and transportation systems still remain unknown at this point in
time. Against the background of this unprecedented global crisis, questions
remain as to how the different factors during the pandemic affect human
mobility and travel.

With the increasing availability of high-quality data related to COVID-19,
analyzing transportation and mobility during and after this crisis is
imperative. Although there are still many unknowns, statistical models and
analytical tools would help produce evidence-based research and policy
interventions after COVID-19 outbreak. At the national level, Zhang et al.
(2020) developed a COVID-19 impact analysis platform [6] that can inform users
about the spread of COVID-19 in the U.S., and the effects of the virus spread
and government orders on mobility and social distancing in the country, using
privacy-protected smartphone device location data, coupled with the census
information. The platform gets updated daily and goes back to January 1, 2020,
for benchmarking, and the results are scaled and aggregated to the entire
population for both state and county levels [7]. Gao et al. (2020) reported on
the interactive web-based mapping platform [8] developed by the GeoDS Lab at
the University of Wisconsin, Madison, with the support of the National Science
Foundation RAPID program [9]. This platform provides information on how people
in different counties and states in the U.S. responded to the social
distancing directives and guidelines. The platform integrates geographic
information systems (GIS) and the daily updated human mobility patterns
obtained from large-scale anonymous and aggregated smartphone data at the
county-level [10].

In recent times, research in several domains is heavily focusing on large-
scale data analysis utilizing sophisticated computing capabilities and machine
learning [11, 12, 13, 14, 15, 16, 17]. However, causal analysis of data for
prediction purposes has received limited attention in the literature.
Causality has been one of the oldest and under-examined questions in all of
science. However, conceptually, its importance has been widely acknowledged in
prior studies and used in various scientific disciplines including social
media research [18], neuroscience [19], biological networks [20, 21], and
economics and finance [22, 23, 24], among others. Causal analysis and
influence characterization were initially geared towards time-series analysis
and among the different measures, Granger causality [22, 24], directed
information [25, 26], and Schreiber’s transfer entropy [27] have gained the
most attention. In the realm of dynamical systems, a new definition and
measure of causal characterization, called information transfer, has been
proposed recently [28, 29, 30], where the authors show that the existing
definitions of causality, namely, Granger causality, directed information, and
transfer entropy fail to capture the correct causal structure in a dynamical
system. Additionally, some recent studies [31, 32] provided a data-driven
approach to infer the causal structure of a dynamical system from time-series
data.

Although there are different measures of causality, for time-series data
analysis, Granger causality is one of the simplest and most widely used
methods. Broadly, in econometric studies, Granger causality test has been
popular in analyzing time-series data for identifying influential factors and
prediction purposes. Konstantakis et al. (2017) employed Granger causality,
along with other quantitative techniques, to investigate the factors that
affect automobile sales in Greece [33]. More recently, Homolka et al. (2020)
determine macro and socio-economic indicators that may significantly predict
car registrations across European countries with the help of Granger causality
test and Vector Autoregressive models [34]. Also, in the domain of
transportation research, particularly in predictive analysis with big data,
Granger causality has been employed widely. Beyzatlar et al. (2014) analyzed
panel data from fifteen European countries (EU-15) using Granger causality to
investigate the relationship between income and transportation [35]. In order
to accurately build the traffic flow prediction models, Li et al. (2015)
utilized Granger causality to determine the potential dependence among the
pool of predictor variables in the time-series big data collected by different
sensors [36]. In a relevant study, McMullen and Eckstein (2012) tested Granger
causality between vehicle miles traveled (VMT) and various measures of the
national economic activity over time [37]. The authors rightly argued that it
is imperative to properly understand the relationship between VMT and other
pertinent transportation and economic factors, as VMT trend is one of the key
components in transportation policies, in general.

Moreover, regression analysis is commonly used in data analysis and machine
learning [38], and is applied extensively in various domains, including
transportation [14, 39, 40, 41, 42]. Among the various regression algorithms,
linear regression is computationally one of the most efficient techniques and
is often used as a starting point for many different problems, although in
many cases, the linear regression model may be suboptimal. In the analysis of
time-series data coming from a dynamical system, linear regression arises
naturally in data-driven analysis of dynamical systems using the transfer
operators, namely Perron-Frobenius and Koopman operators [43, 44, 45, 46].
More recently, studies [47, 48] used robust optimization techniques to compute
the aforementioned operators from noisy data sets, and the resulting
optimization problem was a variation of the ordinary least squares (OLS),
namely, least squares with regularization. Regularization is a standard
technique used in data analysis to overcome some of the limitations of OLS,
including the overfitting and susceptibility to noise in the data [38].

The literature on transportation analysis utilizing regularization is scant.
Recently Polson and Sokolov (2017) developed a deep learning model utilizing a
linear model that is fitted with $\ell_{1}$-regularization to predict traffic
flows. The study showed that deep learning architecture was capable of
capturing non-linear spatio-temporal effects in traffic and providing short-
term predictions of traffic flow [49]. Tan et al. (2011) proposed a semi-
supervised Elastic Net regression method for pedestrian counting by utilizing
sequential information between unlabelled samples and their temporally
neighboring samples as a regularization term. The developed model was able to
attain superior prediction performance and select representative features from
the original set of features without losing their interpretability [50]. Hasan
et al. (2017) proposed statistical techniques to identify spatial
relationships among road links in an urban road network to select predictors
for a short-term traffic prediction model for a given road link. The study
used a time-lagged multiple linear regression method and utilized two
analytical methods, including the Granger causality test and Elastic Net
regularization, using one year of traffic flow and speed data from the
selected road network in Brisbane, Australia. For a given target link, the
relevant predictors obtained by the Granger causality and Elastic Net were
used separately to build the respective traffic prediction models. The results
showed that Granger causality-based traffic prediction model provided superior
prediction accuracy than that using the Elastic Net regression [51]. More
recently, Battifarano and Qian (2019) explored the spatio-temporal
correlations between the urban environment, traffic flow characteristics, and
surge multipliers and proposed a general framework for predicting the short-
term evolution of surge multipliers in real-time using a log-linear model with
$\ell_{1}$-regularization, integrated with pattern clustering. The modeling
algorithm was validated by using Uber and Lyft data from Pittsburgh [52].

While there is a plethora of information available related to mobility and
travel during this pandemic, it is critical to develop a robust methodological
framework to accurately identify and analyze the key factors influencing human
mobility and subsequently predict mobility using the selected factors in the
time of such health crisis. To this end, this study develops an analytical
framework that helps determine the most significant factors affecting mobility
by utilizing Granger causality followed by predicting mobility using linear
regularization algorithms including the Ridge, and Least Absolute Shrinkage
and Selection Operator (LASSO) modeling techniques.

The paper is organised as follows. Section 2 describes the theoretical
concepts of Granger causality test, while Section 3 explains the analytical
methodology of ordinary least square method, and Ridge and LASSO linear
regularization algorithms. Additionally, Section 4 describes the time-series
data analyzed in this study along with its descriptive statistics. Moreover,
Sections 5 and 6 discuss the results of the causal analysis and the prediction
of the regression models developed. Finally, a summary and conclusion of this
study along with its limitations and future research scope is included in
Section 7.

##  2 Granger Causality

In this section, we briefly discuss the concept of Granger causality with its
two different methodological approaches. Granger causality [22, 24] is a
quantitative measure for inferring causal relationships between variables of a
time-series data. It is based on the following two principles.

###### Assumption 1

Cause happens before the effect and the cause has unique information about the
future of the effect.

The intuition behind the definition of Granger causality is as follows.
Suppose the goal is to predict the future of a variable $Y$. If it happens
that the prediction of $Y$ improves by considering the past of variables $Y$
and $X$ as opposed to considering the past of only $Y$, then it is said that
$X$ “Granger causes” $Y$.

###### Definition 2

Under assumption 1, the hypothesis for testing Granger causality of $X$ on $Y$
is

$\displaystyle P[Y(t+1)\in\Omega|I(t)]\neq P[Y(t+1)\in\Omega|I_{\not{X}}(t)]$
(1)

where $P[Y(t+1)\in\Omega|I(t)]$ is the probability of $Y(t+1)$ belonging to
the set $\Omega$ when the entire information till time $t$ is considered
$(I(t))$ and $P[Y(t+1)\in\Omega|I_{\not{X}}]$ is the probability of $Y(t+1)$
belonging to the set $\Omega$ when $X$ is removed from the information set
(denoted by $I_{\not{X}}(t)$). When the above hypothesis (1) is satisfied,
then we say $X$ Granger causes $Y$.

###  2.1 Bivariate Granger Causality

Let $X_{t}$ and $Y_{t}$ be two time-series data, individually each of which
can be represented by the following regressive models:

$\displaystyle\begin{aligned}
&X_{t}=\sum_{j=1}^{\infty}\alpha_{xj}X_{t-j}+\epsilon_{xt},\quad\textnormal{var}(\epsilon_{xt})=\sigma_{x}\\\
&Y_{t}=\sum_{j=1}^{\infty}\beta_{yj}Y_{t-j}+\epsilon_{yt},\quad\textnormal{var}(\epsilon_{yt})=\sigma_{y}.\end{aligned}$
(2)

When considered jointly, the bivariate autoregressive models are:

$\displaystyle\begin{aligned}
&X_{t}=\sum_{j=1}^{\infty}\alpha_{xj}^{\prime}X_{t-j}+\sum_{j=1}^{\infty}\alpha_{xyj}^{\prime}Y_{t-j}+\epsilon_{xt}^{\prime}\\\
&Y_{t}=\sum_{j=1}^{\infty}\beta_{yxj}^{\prime}X_{t-j}+\sum_{j=1}^{\infty}\beta_{yj}^{\prime}Y_{t-j}+\epsilon_{yt}^{\prime},\end{aligned}$
(3)

where $\epsilon_{xt}^{\prime}$ and $\epsilon_{yt}^{\prime}$ are uncorrelated
over time, such that their covariance matrix $\Sigma_{xy}$ is

$\displaystyle\Sigma_{xy}^{\prime}=\begin{pmatrix}\sigma_{xx}^{\prime}&\sigma_{xy}^{\prime}\\\
\sigma_{yx}^{\prime}&\sigma_{yy}^{\prime}\end{pmatrix}.$ (4)

With this, the Granger causality from each of the variable to the other is
given by

$\displaystyle G_{X\to Y}=\ln\frac{\sigma_{y}}{\sigma_{yy}^{\prime}};\quad
G_{Y\to X}=\ln\frac{\sigma_{x}}{\sigma_{xx}^{\prime}}.$ (5)

Furthermore, the interdependence between the variables $X_{t}$ and $Y_{t}$ is
given by

$\displaystyle
G_{X,Y}=\ln\frac{\sigma_{x}\sigma_{y}}{|\Sigma_{xy}^{\prime}|},$ (6)

where $|\cdot|$ denotes the determinant of a matrix. Note that when $X_{t}$
and $Y_{t}$ are independent,
$|\Sigma_{xy}^{\prime}|=\sigma_{xx}^{\prime}\sigma_{yy}^{\prime}$ and hence
$G_{X,Y}$ is zero. When the interdependence is non-zero, in [53], it was shown
that the total interdependence can be decomposed as

$\displaystyle G_{X,Y}=G_{X\to Y}+G_{Y\to X}+G_{X\cdot Y},$ (7)

where $G_{X\cdot
Y}=\ln\frac{\sigma_{xx}^{\prime}\sigma_{yy}^{\prime}}{|\Sigma_{xy}^{\prime}|}$
is the instantaneous causality between $X_{t}$ and $Y_{t}$.

###  2.2 Conditional Granger Causality

Equations (5) give the Granger causality values for a bivariate time-series
data. However, in most real-life applications, the obtained time-series data
has more than two variables. In such cases, computing pairwise dependence
using bivariate Granger causality may lead to ambiguous results [54, 55]
because there may be direct or indirect causal links and in such
circumstances, conditional Granger causality [54] is more appropriate to infer
causality. For simplicity, the case with three variables is discussed and the
general case with more variables follows directly.

Consider three time-series data $X_{t}$, $Y_{t}$ and $Z_{t}$ and suppose
$Y_{t}$ has a pairwise causal influence on $X_{t}$. We now consider the causal
influence of $Y_{t}$ on $X_{t}$ that is mediated through $Z_{t}$. Let the
joint autoregressive model of $X_{t}$ and $Z_{t}$ be

$\displaystyle\begin{aligned}
&X_{t}=\sum_{j=1}^{\infty}\alpha_{xj}X_{t-j}+\sum_{j=1}^{\infty}\alpha_{xzj}Z_{t-j}+\epsilon_{xt}\\\
&Z_{t}=\sum_{j=1}^{\infty}\gamma_{zxj}X_{t-j}+\sum_{j=1}^{\infty}\gamma_{zj}Z_{t-j}+\epsilon_{zt},\end{aligned}$
(8)

with residual covariance matrix
$\Sigma_{xz}=\begin{pmatrix}\sigma_{xx}&\sigma_{xz}\\\
\sigma_{zx}&\sigma_{zz}\end{pmatrix}$.

Again, the joint autoregressive model for all the variables $X_{t}$, $Y_{t}$
and $Z_{t}$ is

$\displaystyle\begin{aligned}
&X_{t}=\sum_{j=1}^{\infty}\alpha_{xj}^{\prime}X_{t-j}+\sum_{j=1}^{\infty}\alpha_{xyj}^{\prime}Y_{t-j}+\sum_{j=1}^{\infty}\alpha_{xzj}^{\prime}Z_{t-j}+\epsilon_{xt}^{\prime}\\\
&Y_{t}=\sum_{j=1}^{\infty}\beta_{yxj}^{\prime}X_{t-j}+\sum_{j=1}^{\infty}\beta_{yj}^{\prime}Y_{t-j}+\sum_{j=1}^{\infty}\beta_{yzj}^{\prime}Z_{t-j}+\epsilon_{yt}^{\prime}\\\
&Z_{t}=\sum_{j=1}^{\infty}\gamma_{zxj}^{\prime}X_{t-j}+\sum_{j=1}^{\infty}\gamma_{zyj}^{\prime}Y_{t-j}+\sum_{j=1}^{\infty}\gamma_{zj}Z_{t-j}+\epsilon_{zt}^{\prime},\end{aligned}$
(9)

with the residual covariance being

$\displaystyle\Sigma_{xyz}^{\prime}=\begin{pmatrix}\sigma_{xx}^{\prime}&\sigma_{xy}^{\prime}&\sigma_{xz}^{\prime}\\\
\sigma_{yx}^{\prime}&\sigma_{yy}^{\prime}&\sigma_{yz}^{\prime}\\\
\sigma_{zx}^{\prime}&\sigma_{zy}^{\prime}&\sigma_{zz}^{\prime}\\\
\end{pmatrix}.$ (10)

With this, the Granger causality of $Y_{t}$ on $X_{t}$ conditioned on $Z_{t}$
is defined as

$\displaystyle G_{Y\to X|Z}=\ln\frac{\sigma_{xx}}{\sigma_{xx}^{\prime}}.$ (11)

With this, if $G_{Y\to X|Z}>0$ and bivariate Granger analysis shows that
$G_{Y\to X}\neq 0$, then the inclusion of $Y$ results in an improved
prediction of $X$ and one can conclude that $Y$ influences $X$ directly. If
$G_{Y\to X|Z}=0$ and bivariate Granger causality analysis results in $G_{Y\to
X}\neq 0$, then the influence of $Y_{t}$ on $X_{t}$ is entirely through
$Z_{t}$.

In the cases where there are more than three variables, conditional Granger
causality can be defined similarly with the scalar $\sigma_{ij}$ and
$\sigma_{ij}^{\prime}$ replaced by corresponding elements of the residual
covariance matrices. However, it is to note that although Granger causality
identifies the most influential variables, it does not provide the direction
of the association between the explanatory and the dependent variables.

##  3 Linear Models

In this section, we discuss the simple linear regression (ordinary least
square or OLS) model and two of its variants, namely Ridge and LASSO
regression models.

###  3.1 Linear Regression

Consider the data set $\\{y_{i},x_{i1},x_{i2},\cdots,x_{iN}\\}_{i=1}^{n}$
where $y_{i}$ is the $i^{th}$ observation of the dependent variable $y$ and
$x_{j}$, $j=1,2,\cdots,N$ are the $N$ independent variables. In case of linear
regression, the model tries to fit a straight line by minimizing the
residuals. In particular, it assumes that the dependent variable can be
expressed as a linear combination of the independent variables, as given by
the following,

$\displaystyle
y_{i}=\alpha_{0}+\alpha_{1}x_{i1}+\alpha_{2}x_{i2}+\cdots+\alpha_{N}x_{iN}+\epsilon_{i},\quad
i=1,2,\cdots,n,$ (12)

where $\epsilon_{i}$ is the residual.

In matrix form, the equation (12) can be written as

$\displaystyle{\bf Y}={\bf X}{\bf\alpha}+{\bf\epsilon},$ (13)

where

$\displaystyle{\bf Y}=[y_{1},y_{2},\cdots,y_{n}]^{\top},{\bf
X}=\begin{pmatrix}1&x_{11}&\cdots&x_{1N}\\\ 1&x_{21}&\cdots&x_{2N}\\\
\vdots&\vdots&\ddots&\vdots\\\
1&x_{n1}&\cdots&x_{nN}\end{pmatrix},{\bf\alpha}=[\alpha_{0},\alpha_{1},\cdots,\alpha_{N}]^{\top},{\bf\epsilon}=[\epsilon_{1},\epsilon_{2},\cdots,\epsilon_{n}]^{\top}.$

The linear regression selects the parameters $\alpha_{j}$’s $(j=0,\cdots,N)$
such that the norm of the residual for every $y_{i}$ $(i=1,\cdots,n)$ is
minimized. Hence the optimal ${\bf\alpha}$ is obtained as a solution of the
following optimization problem,

$\displaystyle\min_{\bf\alpha}\parallel{\bf Y}-{\bf
X}{\bf\alpha}\parallel_{2}.$ (14)

where $\parallel\cdot\parallel_{2}$ is the 2-norm of a vector. The
optimization problem (14) is convex and can be solved efficiently either using
convex optimization techniques or analytically, such that the optimal
${\bf\alpha}^{\star}$ is given by

$\displaystyle{\bf\alpha}^{\star}={\bf Y}{\bf X}^{\dagger},$ (15)

where ${\bf X}^{\dagger}$ is the Moore-Penrose inverse of $\bf X$.

###  3.2 Ridge Regularization

A major drawback of linear regression is that this algorithm has low bias and
high variance [38]. This means that the linear regression may perform well on
the train data, but it may not generalize well to the test data set, thereby
making the model performance unsatisfactory. In machine learning literature,
this phenomenon is known as Bias-Variance trade-off [38]. The intuition of
bias-variance trade-off is explained in Figure 1(a). Usually, with a highly
complex model, it is possible to fit the training data as closely as possible.
In this case, the training error is zero and the model is said to have low
bias. However, the highly complex model may not generalize well to the test
data, thus making the test error large. This is due to the overfitting of the
training data. The complex model, which overfits the training data and
produces high test error, is said to have high variance. This situation is
often reversed if the model considered is fairly simple. Ridge regression,
which puts a 2-norm ($\ell_{2}$-norm) constraint on the set of coefficients,
is able to overcome this challenge efficiently.

Figure 1: (a) Bias-variance trade-off. (b) Effect of outlier on linear
regression.

Another drawback of OLS is that the obtained model is highly influenced by the
outliers in the training data set. For example, as in Figure 1(b), the outlier
data point $y_{k}$ results in the linear model fit as represented by the red
line. However, it is obvious by looking at the overall data that the model fit
depicted by the black line is the more appropriate linear fit to the data.

Additionally, on many occasions, the real-life data is noisy or uncertain.
When the OLS attempts to fit that noise in the data, it eventually results in
overfitting and consequently degrading its performance for model prediction.
As stated earlier, the data for this study were obtained mostly from
smartphone devices, and the chances of acquiring this data may also be subject
to individual user’s discretion, so it is reasonable to assume that the data
utilized in this study may contain some noise or uncertainty in it. To account
for the noise in the data, it is assumed that there is some uncertainty,
$\Delta\bf Y$ and $\Delta\bf X$, in the dependent and independent variables,
respectively. It is assumed that the uncertainties in both $\bf Y$ and $\bf X$
are bounded, i.e. there exists some positive real number $\rho>0$ such that
$\parallel\Delta\bf X\parallel_{2}\leq\rho$ and $\parallel\Delta\bf
Y\parallel_{2}\leq\rho$. With this, the optimization problem (14) is modified
to a min-max optimization problem [47, 48] given by,

$\displaystyle\min_{\bf\alpha}\max_{\begin{subarray}{c}\parallel\Delta\bf
X\parallel_{2}\leq\rho,\\\ \parallel\Delta\bf
Y\parallel_{2}\leq\rho\end{subarray}}\parallel(\bf Y+\Delta\bf Y)-(\bf
X+\Delta X)\bf\alpha\parallel_{2}.$ (16)

Min-max optimization problems are generally hard to solve, but in this case,
the optimization problem (16) can be equivalently expressed as a convex
optimization problem as follows,

###### Theorem 3

The optimization problem

$\displaystyle\min_{\bf\alpha}\max_{\begin{subarray}{c}\parallel\Delta\bf
X\parallel_{2}\leq\rho,\\\ \parallel\Delta\bf
Y\parallel_{2}\leq\rho\end{subarray}}\parallel(\bf Y+\Delta\bf Y)-(\bf
X+\Delta X)\bf\alpha\parallel_{2}.$ (17)

is equivalent to the following,

$\displaystyle\min_{\bf\alpha}\parallel\bf Y-\bf
X\bf\alpha\parallel_{2}+\lambda_{2}\parallel\bf\alpha\parallel_{2},$ (18)

where $\lambda_{2}$ is a positive real number, depending on the uncertainty
bound $\rho$.

Proof. For proof, see [47, 48]. $\hfill\square$

The optimization problem (18), known as Ridge regression, is a convex problem
and can be solved efficiently using any of the available convex optimization
problem solvers. The parameter $\lambda_{2}$ is called the regularization
parameter and it acts as a trade-off between the OLS cost and the cost on the
coefficients $\bf\alpha$.

###  3.3 LASSO Regularization

Apart from the 2-norm regularization, another popular regularization is the
$\ell_{1}$-norm regularization, which is also known as Least Absolute
Shrinkage and Selection Operator (LASSO) regression. In particular, instead of
putting a 2-norm cost on the coefficients of the linear model, LASSO employs a
1-norm cost on the coefficients. Hence, the LASSO regression model is obtained
by solving the following optimization problem:

$\displaystyle\begin{aligned} \min_{\bf\alpha}\qquad&\parallel{\bf Y}-{\bf
X}{\bf\alpha}\parallel_{2}+\lambda_{1}\parallel\bf\alpha\parallel_{1}\\\
\textnormal{subject to}\quad&\parallel\bf\alpha\parallel_{1}\leq
t,\end{aligned}$ (19)

where $\parallel\cdot\parallel_{1}$ is the $1$-norm of a vector and the bound
$t$ is the tuning parameter. If $t$ is large, it has no effect on the
regression coefficients $\alpha_{i}$s and in this case, the solution to the
optimization problem (19) approach the solution of normal linear regression
optimization problem (14) in the limit of large $t$. However, when the bound
$t$ is small, the parameters $\alpha_{i}$s are constrained and hence are
shrunk and are smaller versions of the original least squares estimates. The
$1$-norm minimization puts constraints on parameters that shrink coefficients
towards zero and thus leads to a sparse solution for the linear model.

##  4 Data Description

Data for this study were collected and combined from multiple web-based open-
access sources. The majority of data used in this analysis were requested and
obtained from the COVID-19 Impact Analysis Platform developed at the Maryland
Transportation Institute of the University of Maryland (UMD) [6]. This
platform provides both state and county-based information for 50 states in the
U.S. and the District of Columbia. To match with the data available from other
sources, for the purpose of this study, state-wise daily time-series data were
requested from this platform. The relevant statewide data obtained and
analyzed from this source include the daily number of new COVID-19 confirmed
cases per 1,000 people, social distancing index, percent of out of county and
state trips, transit mode share, population, percent of people older than 60
years, percent of African American or Hispanic Americans, median income,
percent of male population, number of hot spots per 1,000 people, unemployment
rate, percent of people working from home, among others. Social distancing
index in the data set indicates the increasing space between individuals and
decreasing frequency of contact and is represented as an integer from 0 to
100, where 0 indicates no social distancing in the state and 100 indicates all
residents are staying at home.

Figure 2: The average vehicle miles traveled during the analysis period.

Additional information was collected and appended with the data obtained from
the UMD platform. The data for daily vehicle miles traveled (VMT) starting
from March 1, 2020 through June 19, 2020 for 48 states (excluding Alaska and
Hawaii), and the District of Columbia were requested and obtained from the
Streetlight Data [56]. The movement trends over time by geography, across
different categories of places such as retail and recreation, groceries and
pharmacies, parks, transit stations, workplaces, and residential, reported by
Google, were further joined with the study data set [57]. Additionally, the
percent of the population at different education levels and gross domestic
product (GDP) information were collected from the U.S. Census Bureau [58].
Information on the percent of population staying at home during this time was
obtained from the Bureau of Transportation Statistics [59]. Furthermore,
several states have exercised travel restrictions in the form of stay at home
order, limitations on gatherings, domestic travel limitations, or school
closures. This information was obtained from the COVID-19 State and Territory
Action Tracker [60] provided by ESRI (Environmental Systems Research
Institute) [61].

Table 1: Descriptive Statistics of the Data Parameter | Minimum | Maximum | Mean | Std. Dev
---|---|---|---|---
Daily vehicle miles traveled ($10^{6}$ miles) | 0.36 | 2,519.35 | 152.48 | 200.31
Number of new COVID-19 cases (per 1,000 people) | 0 | 0.72 | 0.06 | 0.07
Social distancing index | 10 | 83 | 38.05 | 16.07
Percent of out of county trips | 6.9 | 52.1 | 28.01 | 7.17
Percent of out of state trips | 0.4 | 52.1 | 6.8 | 6.51
Total population | 577,737 | 39,557,045 | 6,632,847 | 7,336,262
Transit mode share | 0.29 | 34.83 | 3.66 | 6.34
Percent of population older than 60 years | 15 | 27 | 21.94 | 2.24
Median income (US dollars) | 44,445 | 84,342 | 61,137 | 10,219.58
Percent of Hispanic or African American population | 1.5 | 48.5 | 11.85 | 10.32
Percent of male | 47.46 | 51.34 | 49.32 | 0.73
Number of hot spots (per 1,000 people) | 47.46 | 51.34 | 49.32 | 0.73
Unemployment rate | 2.2 | 46.6 | 16.67 | 9.29
Percent of population working from home | 2.3 | 55.7 | 22.88 | 9.8
Number of tests conducted (per 1,000 people) | 0 | 181.14 | 22.07 | 25.63
Change in the number of trips to transit stations from baseline | -81 | 44 | -27.32 | 22.55
Change in the number of trips to retail and recreation places from baseline | -77 | 24 | -23.9 | 19.15
Change in the number of trips to grocery and pharmacy from baseline | -62 | 48 | -2.84 | 13.01
Change in the number of trips to parks from baseline | -77 | 388 | 38.184 | 60.2
Change in the number of trips to workplaces from baseline | -78 | 18 | -32.08 | 17.55
Change in the number of trips to residential places from baseline | -5 | 33 | 11.87 | 7.31
Population staying at home | 108,677 | 14,180,383 | 1,646,749 | 1,932,225
Socioeconomic status | 0 | 1 | 0.51 | 0.29
State employee travel restriction | 0 | 1 | 0.71 | 0.45
Stay-at-home order or guidance | 0 | 1 | 0.76 | 0.43
School closures | 0 | 1 | 0.96 | 0.20
Closure of some or all facilities | 0 | 1 | 0.86 | 0.35
Mandatory or recommended domestic travel limitations | 0 | 1 | 0.43 | 0.49
Mandatory statewide mask policy | 0 | 1 | 0.69 | 0.46
Median age | 31 | 45.1 | 38.66 | 2.38
GDP in 2019 ($10^{9}$ U.S. dollars) | 30.5 | 2,792.03 | 382.32 | 488.51
Percent of population aged 25 years or over with high school or higher | 83.8 | 93.9 | 89.74 | 2.58

Following the joining of the data from various aforementioned sources, a
thorough screening and quality check of the data was performed for any missing
values. The final data set includes daily data starting from March 1, 2020
through June 13, 2020 from 48 U.S. states excluding Alaska and Hawaii and
including the District of Columbia, and consists of a total of 5,145
observations for further analysis. In the final data set, the two states,
Alaska and Hawaii, were excluded as the VMT data was not available for the
same.

The preliminary analysis examined several predictor variables of interest that
were initially included in the analysis. The distribution of the dependent
variable i.e., the daily VMT during the analysis period (from March 1, 2020
through June 13, 2020) is shown in Figure 2. As can be clearly seen from
Figure 2, the average VMT dropped substantially from the beginning of March,
2020, around the time when the U.S. started experiencing rapid community
outbreaks of the virus and the country declared a national emergency, and it
continued to remain low until around end of April, 2020. However, it is
interesting to see that the daily VMT gradually increased from around early
May, 2020, although the number of daily COVID-19 cases continued to grow
considerably over time.

As stated earlier, the compiled data set includes a large number of variables
of interest that may impact the daily VMT. Table 1 presents the descriptive
statistics, including the minimums, maximums, means, and standard deviations
of thirty-one variables that were initially considered as the set of potential
determinants. For the purpose of scaling the data, VMT, population, median
income, and GDP were included in their natural log forms.

##  5 Causal Analysis of Data

In any multivariate time-series data, there is causal influence between the
variables involved. Correlations between the variables quantify the extent of
their interdependence, but they do not specifically identify the cause and the
effect. In particular, the correlations are symmetric in the variables and
lack a directional sense. However, to better understand the relationships
among the variables in the data, identification of the causal structure and
influential predictor variables is critical. For example, for a large time-
series data with a massive number of predictor variables, it is advisable to
consider a reduced order model, which can be done by identifying the most
influential variables. In this study, the goal of determining the most
important predictor variables from the large set of time-series data was
accomplished by employing the Granger causality test.

Figure 3: Bivariate Granger causality test values of the predictor variables
to dependent variable.

As discussed earlier in section 2, Granger causality between two variables can
be analyzed by either doing a pairwise study by developing bivariate VAR model
or by conditional Granger causality by considering multivariate VAR model. For
multivariate data, bivariate Granger causality analysis often provides
ambiguous results [54, 55]. For validation purposes, the study initially
carries out a bivarite Granger causality test whose results are shown in
Figure 3. In this figure, the Granger causality test values of all predictor
variables on the dependent variable (daily VMT) using bivariate models are
displayed. The lags in this VAR model were chosen by Akaike Information
Criterion (AIC) and the optimal lag was found to be 4 for the bivariate VAR
models. From Figure 3, it can be seen that the model was not efficient in
determining causality for most of the variables and several apparently
important predictors were not correctly captured in explaining the causal
relationships between the explanatory and dependent variable. For example, the
most influential predictor variable was identified as _the change in the
percent of trips to grocery and pharmacy from the average baseline_ , whereas
factors like _social distancing index_ , or _population staying at home_ or
_percent of people working from home_ etc. had very little to no causal
influence on the dependent variable. This observation is counter-intuitive and
supports the previous theoretical argument that bivariate Granger causality
tests provide incoherent results for multivariate time-series data.

Figure 4: Multivariate Granger causality test values of the predictor
variables to dependent variable.

For conditional Granger causality test, a multivariate VAR model was developed
and similar to the bivariate model, the optimal lag for the VAR model was
determined by the Akaike Information Criterion and this optimal lag was
obtained as 10. While determining the conditional Granger causality (11) of
the $i^{th}$ predictor variable on the dependent variable, the conditional set
was chosen as all the predictor variables except the $i^{th}$ predictor
variable. Figure 4 presents the results of the conditional Granger causality
test by showing the causality values for each of the predictor variables. As
can be seen from Figure 4, the conditional Granger causality efficiently
identified the important predictors from a pool of a large number of
explanatory variables. For example, unlike the bivariate causality test,
_social distancing index_ or _population staying at home_ were identified as
the important variables in explaining the daily vehicle miles traveled.

##  6 Regression Models and Prediction

In the previous section, the influence of the predictor variables on the VMT
was computed and was rank-ordered according to the influence each predictor
variable has on the VMT. In this section, the most influential predictor
variables identified by Granger causality test were used to develop three
different reduced order linear models, namely, ordinary least squares, LASSO,
and Ridge models. In particular, the largest value of the causal influence was
around 0.03 and the cut-off was selected to be 0.003. This cut-off resulted in
the selection of seventeen predictor variables and with these selected
variables, reduced order models for predicting VMT were developed.

###  6.1 Model Training

The final data set used in this study included information from March 1, 2020
through June 13, 2020. After splitting the data set into two parts for
training and testing, the train and test data sets ranged from March 1, 2020
through May 18, 2020 and from May 19, 2020 through June 13, 2020,
respectively. The first (train) part was utilized for training the regression
models, while the latter part (test) was used to test the efficiency of the
prediction of the obtained regression models. As mentioned earlier, the
seventeen most influential variables were considered for the reduced order
models. It is important to note that, except for the OLS model, the other
computation of Ridge and LASSO models (optimization problems (18) and (19))
involve one regularization parameter and in all the optimization problems the
$\lambda_{i}$s were chosen such that, $0<\lambda_{i}\leq 1$.

####  6.1.1 Interpretation and Discussion of the Model Results

Table 2 compares the optimal coefficients of the final set of explanatory
variables between all three regression modeling techniques. As can be seen
from Table 2, the coefficients of the predictors were fairly comparable across
all modeling techniques. As expected, when the number of new COVID cases per
1,000 people increased, the daily vehicle miles traveled decreased.
Expectedly, with the increase in the social distancing index, percent of
people working from home, statewide closure, and tests done per 100 people,
the daily VMT decreased. Additionally, the vehicle miles traveled per day
increased with the increase in population, unemployment rate, person of out of
county trips, socio-economic status, and the increase in percent of trips to
transit stations, retail and recreational places, grocery and pharmacy,
workplaces, and residences. The predictor for population staying at home was
rightly captured by regularization methods, namely LASSO and Ridge
regressions, and showed a negative association with daily VMT. Lastly,
although the negative association between the increase in the percent of trips
to parks and daily VMT seemed counter-intuitive, this could be partially due
to the reason that people might also rely on non-motorized transport to go to
parks and such trips would not be captured by the daily VMT.

Furthermore, the log-lambda plots shown in Figures 5 and 6 depicts how the
predictor variables that enter in the model, varied across the LASSO and Ridge
regression techniques as the regularization parameter change. When the
regularization parameter $\lambda$ is small, the contribution of the
regularization part to the cost functions in the optimization problems (18)
and (19) is small and as such all these optimization problems are reduced to
the simple linear regression or OLS. However, as $\lambda$ is increased, the
weight on the regularization component in the optimization problems increases
and as such, the coefficients of the independent variables become smaller and
approach zero. It is clearly seen from Figures 5 and 6 that all the
coefficients did not approach zero at the same time. In particular, the
important variables remained non-zero for larger values of $\lambda$ as
compared to the relatively non-important variables. Based on the order of
importance, the coefficients of the explanatory variables differed slightly
between the LASSO and Ridge regularization techniques. This difference could
be due to the characteristic of the Ridge regularization that reduces the norm
of the coefficients more uniformly, while the LASSO model attempts to set as
many coefficients to zero as possible.

Table 2: Comparison of Model Coefficients between the Linear, LASSO, and Ridge Regressions Parameter | Linear | LASSO | Ridge
---|---|---|---
Intercept | 3.119 | -1.319 | -1.049
New COVID-19 cases per 1,000 people | -0.010 | -0.148 | -0.128
Social distancing index | -0.844 | -0.266 | -0.169
Percent of out of county trips | 2.126 | 1.064 | 1.139
Percent of out of state trips | -1.854 | -0.177 | -0.263
Ln(Population) | 0.607 | 2.504 | 2.351
Unemployment rate | 0.005 | 0.005 | 0.005
Percent of population working from home | -0.221 | 0.356 | 0.307
Tests done per 100 people | -0.015 | -0.021 | -0.019
Trips to transit stations-percent change from baseline | 1.081 | 1.920 | 1.899
Trips to retail and recreation places-percent change from baseline | 0.707 | 0.225 | 0.145
Trips to grocery and pharmacy-percent change from baseline | 0.717 | 1.265 | 1.199
Trips to parks-percent change from baseline | -0.073 | -0.066 | -0.069
Trips to workplaces-percent change from baseline | 0.387 | 0.046 | 0.033
Trips to residence-percent change from baseline | 0.327 | 0.847 | 0.788
Population staying at home | 0.432 | -1.368 | -1.216
Socioeconomic status | 0.731 | 0.486 | 0.501
Statewide closure of some or all facilities | -0.175 | -0.136 | -0.139
Figure 5: Relative importance of the predictors in terms of log-lambda for
LASSO regression. Figure 6: Relative importance of the predictors in terms of
log-lambda for Ridge regression.

###  6.2 Prediction Performance

Following the development of the models using the training data, the
prediction of the dependent variable (daily VMT) was tested and compared
between the three regression techniques. Essentially, the performance of the
three different models on the test data set was evaluated, and the predicted
values were compared with the observed values. Ultimately, the root mean
square errors (RMSEs) from all models provided the measure of performance and
efficiency of the models.

Table 3: Root mean square error (RMSE) results of training and test data Modeling Technique | RMSE in Training | RMSE in Test
---|---|---
Linear | 0.3733 | 0.4330
LASSO | 0.3346 | 0.3998
Ridge | 0.3348 | 0.3952

The RMSEs of the three different models utilizing both the train and test data
set are presented in Table 3. The comparison of the RMSEs between the three
modeling techniques clearly shows that the Ridge regression performed the best
for both train and test data by having the least RMSE among all models. This
is expected, as the Ridge regression provides superior prediction by
overcoming the issue of overfitting with low variance and better
generalization to the test data compared to other regularization methods
(refer to Figure 1). Additionally, based on the RMSEs, also the LASSO model
provided better prediction performance compared to the OLS. However, the
higher error in the test data in all models could partially be due a much
smaller sample size and the consistent gradual increase in VMT in the test
period as opposed to the change in VMT from the average baseline in both
directions during the training data period.

For the purpose of showing how the different modeling techniques perform at an
individual state level, a comparative graphical representation of RMSEs in the
prediction of daily VMT for the three regression techniques are presented
separately for the states of Florida, New Jersey, New York, and Texas as
examples in Figure 7. These figures clearly show that even at the individual
state level, normal linear model provided the poorest performance in terms of
large errors, while Ridge and LASSO models showed comparable prediction
performance. Among the three models, Ridge showed the most superior
performance, followed by LASSO regression with a slightly higher prediction
error.

Figure 7: Error in prediction of daily VMT in the states of (a) Florida, (b)
New Jersey, (c) New York and (d) Texas .

Figure 8: (a) Root mean square error (RMSE) in prediction of the different
models (Eq. (21)) over the analysis period. (b) Average date-wise root mean
square error (RMSE).

Additionally, the date-wise RMSEs in the prediction, utilizing the three
developed models over the analysis period are presented in Figure 8. Let
$y_{d_{i}}$ be the observed value of daily VMT on the $d^{th}$ day for the
state $i$. For the test data set, $d=1,2,\cdots,26$ and $i=1,2,\cdots,49$. Let
$\hat{y}_{d_{i}}^{p}$ be the predicted daily VMT on the $d^{th}$ day for
$i^{th}$ state when using the $p^{th}$ modeling technique. Here $p$ is one of
linear, LASSO or Ridge regression techniques. Therefore, the error in the
prediction using the $p^{th}$ model on the $d^{th}$ day for $i^{th}$ state is,

$\displaystyle\epsilon_{d_{i}}^{p}=y_{d_{i}}-\hat{y}_{d_{i}}^{p}.$ (20)

The RMSE in the prediction of the $p^{th}$ model for the $d^{th}$ day over all
the states is,

$\displaystyle
r_{d}^{p}=\sqrt{\frac{\sum_{i=1}^{n_{s}}(\epsilon_{d_{i}}^{p})^{2}}{n_{s}}}=\sqrt{\frac{\sum_{i=1}^{n_{s}}(y_{d_{i}}-\hat{y}_{d_{i}}^{p})^{2}}{n_{s}}},$
(21)

where $n_{s}=49$ is the number of the states and D.C.

In Figure 8(a), $r_{d}^{p}$ is plotted for all modeling techniques over the
period ranging between May 19, 2020 through June 13, 2020 (test period). From
Figure 8(a), it can once again be seen that the Ridge regression performed the
best for most of the test period, as its RMSE was usually lower compared to
the other methods. The average of $r_{d}^{p}$ across the test period, which
can be expressed as,

$\displaystyle
r^{p}=\frac{\sum_{d=1}^{n_{d}}r_{d}^{p}}{n_{d}}=\left(\sum_{d=1}^{n_{d}}\sqrt{\frac{\sum_{i=1}^{n_{s}}(y_{d_{i}}-\hat{y}_{d_{i}}^{p})^{2}}{n_{s}}}\right)\Bigg{/}n_{d},$
(22)

where $n_{d}$ is the number of days (equals to 26 in this study), is plotted
in Figure 8(b). This plot also confirms that the Ridge regression model has
the least error in prediction.

##  7 Summary and Conclusions

Since the emergence and rapid growth of novel coronavirus (COVID-19) in
December 2019, countries worldwide are taking extreme measures to prevent the
spread of the virus. The U.S. is greatly hit by the pandemic and currently
(September 2020) has the highest number of confirmed COVID cases and deaths in
the world. Since the national emergency declared by the White House on March
13, 2020, most states in the U.S. implemented travel restrictions and social
distancing protocols to combat the crisis, causing drastic reductions in
mobility and travel demand at all levels. However, the overall impact and the
long-term implications of this crisis to mobility and travel still remain
unknown at this point in time. In order to understand these implications
better, statistical models and analytical tools utilizing the increasingly
available open-access data is the need of the hour. To that end, this study
developed an analytical framework that helped determine and analyze the most
important factors impacting and predicting human mobility and travel in the
U.S. during the pandemic by employing Granger causality and linear
regularization algorithms, including the Ridge, and LASSO modeling techniques.

Data for this study were obtained for 48 states (excluding Alaska and Hawaii)
and the District of Columbia from various databases created and maintained to
analyze the impacts of this pandemic. The data obtained and analyzed in this
study included information on mobility and movement trends, travel
restrictions and social distancing, and health and demographics of the
population. The compiled data set included daily time-series data starting
from March 1, 2020 through June 13, 2020.

Evaluating a large-scale data set requires advanced analytical techniques to
identify the most important factors in explaining the response variable.
Commensurate with analyzing such rich data, this study employed Granger
causality and linear regularization techniques, including the Ridge, and LASSO
models, along with ordinary least square regression. The entire data set was
split into two parts, where approximately 75 percent of the data was used for
training the models, and the remaining 25 percent was used to test the
prediction performance. Determining the set of most important predictors
impacting the daily VMT from the pool of several potential determinants was
accomplished using the Granger causality. The seventeen selected variables
were further used to develop reduced order models, employing the linear,
Ridge, and Lasso regression techniques on the training data. Finally, the
performance of the prediction was tested by feeding the test data into the
models for all regression techniques.

The results of this study revealed that the coefficients of the predictors
were comparable across all modeling techniques. When factors including the
number of new COVID cases, the social distancing index, percent of people
working from home, statewide closure, and tests done per 100 people increased,
the daily VMT decreased. Conversely, the vehicle miles traveled per day
increased with the increase in population, unemployment rate, person of out of
county trips, socio-economic status, and the increase in percent of trips to
transit stations, retail and recreational places, grocery and pharmacy,
workplaces, and residences. The population staying at home was rightly
captured by regularization methods, namely LASSO and Ridge regressions, and
shows a negative association with the daily VMT.

Furthermore, the developed models were used to predict the daily VMT for all
the states for a period of 26 days (from May 19, 2020 through June 13, 2020).
Although all the developed models compare favorably, the Ridge regression
model performed the best by having the least root mean square error (RMSE) in
prediction among all models. This result makes sense because the Ridge
regression is robust in overcoming the issue of overfitting and thus
generalizes better to the test data set, resulting in lesser prediction error.
Also, LASSO regularization technique performed superior to the ordinary least
square regression.

The study is only the starting point to help understand the associations
between different factors and human mobility during the COVID-19 pandemic. The
authors of this study intend to expand this study to utilize county-based data
to understand these associations from a more granular level. Moreover, it
would be insightful to include additional variables into the analysis as
potential predictors. From the modeling perspective, it is reasonable to argue
that the available data is subjected to some uncertainties and future research
should be carried out to explicitly take the uncertainties into account to
derive at more precise models. Furthermore, as the crisis is moving on to the
greater peaks in terms of the number of confirmed cases and deaths over time
in the U.S., subsequent analysis is warranted with data from the following
months (post June 13, 2020).

##  8 Acknowledgement

None. This research did not receive any specific grant from funding agencies
in the public, commercial, or not-for-profit sectors.

##  9 CRediT Author Statement

Subhrajit Sinha: Conceptualization, Methodology, Validation, Formal Analysis,
Investigation, Data Curation, Original Draft Preparation, Review and Editing.
Meghna Chakraborty: Methodology, Validation, Formal Analysis, Investigation,
Data Curation, Original Draft Preparation, Review and Editing.

## References

  * [1] World Health Organization (2020 (accessed July 19, 2020)).
URL https://www.who.int/emergencies/diseases/novel-coronavirus-2019

  * [2] Proclamation on Declaring a National Emergency Concerning the Novel Coronavirus Disease (COVID-19) Outbreak (2020 (accessed July 19, 2020)).
URL https://www.whitehouse.gov/presidential-actions/proclamation-declaring-
national-emergency-concerning-novel-coronavirus-disease-covid-19-outbreak/

  * [3] COVID-19 Coronavirus Pandemic (2020 (accessed July 19, 2020)).
URL https://www.worldometers.info/coronavirus/

  * [4] Helen Davidson, Around $20\%$ of global population under coronavirus lockdown (2020 (accessed July 19, 2020)).
URL https://www.theguardian.com/world/2020/mar/24/nearly-20-of-global-
population-under-coronavirus-lockdown

  * [5] COVID19.CA.GOV, Stay at home $Q\&A$ (2020 (accessed July 19, 2020)).
URL https://covid19.ca.gov/stay-home-except-for-essential-needs/

  * [6] Maryland Transportation Institute, University of Maryland COVID-19 Impact Analysis Platform, University of Maryland, College Park, USA (2020 (accessed July 19, 2020)).
URL https://data.covid.umd.edu

  * [7] L. Zhang, S. Ghader, M. L. Pack, C. Xiong, A. Darzi, M. Yang, Q. Sun, A. Kabiri, S. Hu, An interactive covid-19 mobility impact and social distancing analysis platform, medRxiv (2020).
  * [8] Gao, Song and Rao, Jinmeng and Kang, Yuhao and Liang, Yunlei and Kruse, Jake, Mapping Mobility Changes in Response to COVID-19 (2020 (accessed July 19, 2020)).
URL https://geods.geography.wisc.edu/covid19/physical-distancing/

  * [9] National Science Foundation, National Science Foundation awards rapid response grants to support coronavirus (COVID-19) research (2020 (accessed July 19, 2020)).
URL https://www.nsf.gov/news/special_reports/announcements/041720.jsp

  * [10] S. Gao, J. Rao, Y. Kang, Y. Liang, J. Kruse, Mapping county-level mobility pattern changes in the united states in response to covid-19, SIGSPATIAL Special 12 (1) (2020) 16–26.
  * [11] J. Golbeck, C. Robles, K. Turner, Predicting personality with social media, in: CHI’11 extended abstracts on human factors in computing systems, 2011, pp. 253–262.
  * [12] X. Zheng, W. Chen, P. Wang, D. Shen, S. Chen, X. Wang, Q. Zhang, L. Yang, Big data for social transportation, IEEE Transactions on Intelligent Transportation Systems 17 (3) (2015) 620–630.
  * [13] H. Nguyen, L.-M. Kieu, T. Wen, C. Cai, Deep learning methods in transportation domain: a review, IET Intelligent Transport Systems 12 (9) (2018) 998–1004.
  * [14] M. Chakraborty, M. S. Mahmud, T. Gates, S. Sinha, Linear regularization-based analysis and prediction of human mobility in the us during the covid-19 pandemic (2020).
  * [15] M. Chakraborty, S. Stapleton, M. Ghamami, T. Gates, Safety effectiveness of all-electronic toll collection systems., Advances in Transportation Studies (2020).
  * [16] M. Chakraborty, H. Singh, P. T. Savolainen, T. J. Gates, Examining correlation and trends in seatbelt use among occupants of the same vehicle using a bivariate probit model, Transportation Research Record (2021) 0361198121995487\.
  * [17] B. J. Erickson, P. Korfiatis, Z. Akkus, T. L. Kline, Machine learning for medical imaging, Radiographics 37 (2) (2017) 505–515.
  * [18] G. Ver Steeg, A. Galstyan, Information transfer in social media, in: Proceedings of the 21st international conference on World Wide Web, 2012, pp. 509–518.
  * [19] O. Sporns, The networks of the brain, MIT Press, 2010.
  * [20] D. J. S. A. Rao, A. O. Hero, J. D. Engel, Motif discovery in tissue-specific regulatory sequences using directed information., in: EURASIP J. on Bioinformatics and Systems Biology, 2007.
  * [21] D. J. S. A. Rao, A. O. Hero, J. D. Engel, Inference of biologically relevant gene influence networks using the directed information criterion., in: In proc. ICASSP, Toulouse, France, 2006.
  * [22] C. W. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica: journal of the Econometric Society (1969) 424–438.
  * [23] C. A. Sims, Money, income, and causality, The American economic review 62 (4) (1972) 540–552.
  * [24] C. W. Granger, Causality, cointegration, and control, Journal of Economic Dynamics and Control 12 (2-3) (1988) 551–559.
  * [25] J. L. Massey, Causality, feedback and directed information., in: Proc. Intl. Symp. on Info. th. and its Applications, Waikiki, Hawai, USA, 1990.
  * [26] G. Kramer, Directed information for channels with feedback, in: PhD Thesis, Swiss Federal Institute of Technology Zurich, 1998.
  * [27] T. Schreiber, Measuring information transfer, Physical Review Letters 85, no. 2 (July, 2000) 461–464.
  * [28] S. Sinha, U. Vaidya, Formalism for information transfer in dynamical network, in: 2015 54th IEEE Conference on Decision and Control (CDC), IEEE, 2015, pp. 5731–5736.
  * [29] S. Sinha, U. Vaidya, Causality preserving information transfer measure for control dynamical system, in: 2016 IEEE 55th Conference on Decision and Control (CDC), IEEE, 2016, pp. 7329–7334.
  * [30] S. Sinha, U. Vaidya, On information transfer in discrete dynamical systems, in: 2017 Indian Control Conference (ICC), IEEE, 2017, pp. 303–308.
  * [31] S. Sinha, U. Vaidya, Data-driven approach for inferencing causality and network topology, in: 2018 Annual American Control Conference (ACC), IEEE, 2018, pp. 436–441.
  * [32] S. Sinha, U. Vaidya, On data-driven computation of information transfer for causal inference in discrete-time dynamical systems, Journal of Nonlinear Science (2020) 1–26.
  * [33] K. N. Konstantakis, C. Milioti, P. G. Michaelides, Modeling the dynamic response of automobile sales in troubled times: A real-time vector autoregressive analysis with causality testing for greece, Transport Policy 59 (2017) 75–81.
  * [34] L. Homolka, V. M. Ngo, D. Pavelková, B. T. Le, B. Dehning, Short-and medium-term car registration forecasting based on selected macro and socio-economic indicators in european countries, Research in Transportation Economics 80 (2020) 100752.
  * [35] M. A. Beyzatlar, M. Karacal, H. Yetkiner, Granger-causality between transportation and gdp: A panel data approach, Transportation Research Part A: Policy and Practice 63 (2014) 43–55.
  * [36] L. Li, X. Su, Y. Wang, Y. Lin, Z. Li, Y. Li, Robust causal dependence mining in big data network and its application to traffic flow predictions, Transportation Research Part C: Emerging Technologies 58 (2015) 292–307.
  * [37] B. S. McMullen, N. Eckstein, Relationship between vehicle miles traveled and economic activity, Transportation research record 2297 (1) (2012) 21–28.
  * [38] T. M. Mitchell, et al., Machine learning. 1997, Burr Ridge, IL: McGraw Hill 45 (37) (1997) 870–877.
  * [39] M. Chakraborty, T. J. Gates, Association between driveway land use and safety performance on rural highways, Transportation research record 2675 (1) (2021) 114–124.
  * [40] M. Chakraborty, M. S. Mahmud, T. Gates, Analysis of trends and correlation in child restraint use and seating position of child passengers in motor vehicles: Application of a bivariate probit model (2021).
  * [41] M. Chakraborty, T. Gates, Assessing safety performance on urban and suburban roadways of lower functional classification: A comparison of minor arterial and collector roadway segments (2021).
  * [42] M. Chakraborty, T. J. Gates, Relationship between horizontal curve density and safety performance on rural two-lane road segments by road jurisdiction and surface type, Tech. rep. (2021).
  * [43] C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter, D. HENNINGSON, et al., Spectral analysis of nonlinear flows, Journal of fluid mechanics 641 (1) (2009) 115–127.
  * [44] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data–driven approximation of the koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (6) (2015) 1307–1346.
  * [45] S. Sinha, U. Vaidya, E. Yeung, On computation of koopman operator from sparse data, in: 2019 American Control Conference (ACC), IEEE, 2019, pp. 5519–5524.
  * [46] S. Sinha, S. P. Nandanoori, E. Yeung, Koopman operator methods for global phase space exploration of equivariant dynamical systems, IFAC-PapersOnLine 53 (2) (2020) 1150–1155.
  * [47] S. Sinha, B. Huang, U. Vaidya, Robust approximation of koopman operator and prediction in random dynamical systems, in: 2018 Annual American Control Conference (ACC), IEEE, 2018, pp. 5491–5496.
  * [48] S. Sinha, B. Huang, U. Vaidya, On robust computation of koopman operator and prediction in random dynamical systems, Journal of Nonlinear Science (2019) 1–34.
  * [49] N. G. Polson, V. O. Sokolov, Deep learning for short-term traffic flow prediction, Transportation Research Part C: Emerging Technologies 79 (2017) 1–17.
  * [50] B. Tan, J. Zhang, L. Wang, Semi-supervised elastic net for pedestrian counting, Pattern Recognition 44 (10-11) (2011) 2297–2304.
  * [51] M. M. Hasan, J. Kim, C. Prato, et al., Spatial variable selection methods for network-wide short-term traffic prediction, in: 39 th Australasian Transport Research Forum (ATRF) Proceedings, 2017.
  * [52] M. Battifarano, Z. S. Qian, Predicting real-time surge pricing of ride-sourcing companies, Transportation Research Part C: Emerging Technologies 107 (2019) 444–462.
  * [53] J. Geweke, Measurement of linear dependence and feedback between multiple time series, Journal of the American statistical association 77 (378) (1982) 304–313.
  * [54] J. F. Geweke, Measures of conditional linear dependence and feedback between time series, Journal of the American Statistical Association 79 (388) (1984) 907–915.
  * [55] Y. Chen, S. L. Bressler, M. Ding, Frequency decomposition of conditional granger causality and application to multivariate neural field potential data, Journal of neuroscience methods 150 (2) (2006) 228–237.
  * [56] Streetlight Data, Big Data for Mobility (2020 (accessed July 22, 2020)).
URL https://www.streetlightdata.com/covid-transportation-metrics/

  * [57] COVID-19 Community Mobility Reports, Google (2020 (accessed July 22, 2020)).
URL https://www.google.com/covid19/mobility/

  * [58] U.S. Census Bureau Data (2020 (accessed July 19, 2020)).
URL https://www.census.gov/data.html

  * [59] Bureau of Transportation Statistics, U.S. Department of Transportation (2020 (accessed July 19, 2020)).
URL https://www.bts.gov/

  * [60] National Governors Association, COVID-19 State and Territory Actions Tracker (2020 (accessed July 19, 2020)).
URL https://experience.arcgis.com/experience/fd73fbfbe3e84e53bcbc45f7ca18db2e/

  * [61] Environmental Systems Research Institute (2020 (accessed July 19, 2020)).
URL https://www.esri.com/en-us/home


11institutetext: Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France

# Impact of clustering inside compact tetraquarks

Hagop Sazdjian 11<EMAIL_ADDRESS>
###### Abstract

Due to the reducibility of tetraquark operators into mesonic clusters, a
strong interplay exists in tetraquarks between compact and molecular
structures. This issue is studied within an effective field theory approach,
where the compact tetraquark is treated as an elementary particle. Under the
influence of the coupling to the mesonic clusters, an initially formed compact
tetraquark bound state is deformed towards a new structure of the molecular
type, having the attributes of a shallow bound state.

##  1 Introduction

There has been growing experimental evidence, during the last two decades,
about the existence of exotic hadrons, also called multiquark states,
containing more valence quarks than ordinary mesons ($\overline{q}q$) and
baryons ($qqq$). Prototypes are tetraquarks with valence quark structure
$\bar{q}\bar{q}qq$, pentaquarks, with structure $\overline{q}qqqq$,
hexaquarks, with structure $qqqqqq$ or $\bar{q}\bar{q}\bar{q}qqq$ Choi:2003ue
; Aubert:2003fg ; Besson:2003cp ; Aubert:2005rm ; Ablikim:2013mio ;
Liu:2013dau ; Ablikim:2013wzq ; Aaij:2014jqa ; Aaij:2015tga ; Aaij:2020fnh ;
Aaij:2020ypa ; Wu:2020hmk ; SpadaroNorella:2021yje ; LHCb:2021auc ;
Gershon:2022xnn .

However, contrary to ordinary hadrons, multiquark states are not color-
irreducible, in the sense that they can be decomposed along a finite number of
combinations of ordinary mesonic or baryonic clusters Jaffe:2008zz ;
Lucha:2019cdc . Schematically, ignoring here flavor and spin indices, a
tetraquark interpolating operator, $(\bar{q}\bar{q}qq)$, which is globally
color invariant, can be decomposed by means of Fierz rearrangements into a
form where clusters of mesonic operators have emerged:

$(\bar{q}\bar{q}qq)=\sum(\bar{q}q)(\bar{q}q),$ (1)

where the parentheses signify color invariance of the included operator.
Similar decompositions can be done with the pentaquark and hexaquark
interpolating operators:

$\displaystyle(\overline{q}qqqq)=\sum(\overline{q}q)(qqq),$ (2)
$\displaystyle(qqqqqq)=\sum(qqq)(qqq),\ \ \ \ \
(\bar{q}\bar{q}\bar{q}qqq)=\sum(\bar{q}\bar{q}\bar{q})(qqq)+\sum(\bar{q}q)(\bar{q}q)(\bar{q}q).$
(3)

Hadronic clusters, being color-singlets, mutually interact by means of short-
range forces, like meson-exchanges or contacts. They would then form loosely
bound states, in similarity with atomic molecules. These are called hadronic
molecules or molecular states Voloshin:1976ap ; DeRujula:1976zlg ;
Tornqvist:1993ng .

In contrast, multiquark states, formed directly from confining interactions
acting on all quarks, would form compact bound states, called compact
multiquark states Jaffe:2003sg ; Shuryak:2003zi ; Maiani:2004vq ;
Maiani:2005pe .

Multiquark states can thus be formed by two different mechanisms, each leading
to a different structure. The issue is to find, by theoretical justification,
also guided by experimental data, the most faithful description.

In the following, we concentrate on tetraquark states, ignoring details coming
from quark flavors and spins, which do not play a fundamental role.

##  2 Energy balance

A first hint is provided by the study of the energy balance of the system
Lucha:2021mwx . This is most easily done in the static limit of the theory,
with very heavy quarks, fixed at spatial positions. The system would choose
configurations with minimal energy.

The problem is analytically solved in the strong coupling limit of lattice
theory Dosch:1982ep  and confirmed by direct lattice numerical calculations
Alexandrou:2004ak ; Okiharu:2004wy ; Okiharu:2004ve ; Suganuma:2011ci ;
Cardoso:2011fq ; Bicudo:2017usw . In the strong coupling limit, the potential
energy is concentrated on the Wilson lines, which are the path-ordered gluon
field phase factors, with constant linear energy density.

Figure 1 represents schematic geometric configurations of the two types of
structure. In Fig. 1a, one has a compact tetraquark, formed by confining
interactions, through diquark–antidiquark global interaction. In Fig. 1b, one
has two meson clusters.

Figure 1: (a): Compact tetraquark, formed by confining interactions, through
diquark–antidiquark global interaction. (b): Two meson clusters.

The compact tetraquark formation is energetically favored if the diquark
interdistance $d$ is much smaller than the quark-antiquark interdistance
$\ell$:

$d\ll\ell.$ (4)

This shows that both structures have nonzero probabilities to be formed.
However, quarks have finite masses and move in space. Even if a compact
tetraquark has been formed, there is a sizable probability that the quarks,
during their motion, reach the two-meson-cluster configuration, in which case
the system dislocates or decays.

An interesting and nontrivial case is when the compact tetraquark mass lies
below the two-meson threshold, which prevents dislocation. Nevertheless, there
will be, through fluctuations, constant transitions to the two-meson virtual
states. This will have the tendancy to expand the compact system into a more
loosely bound system, close to a molecular-type state. This is a dynamical
mechanism, whose precise outcome necessitates the resolution of the four-body
bound state problem, in the presence of the confining forces. For the time
being, this problem has not yet been satisfactorily solved. One is obliged, to
analyze the problem with sufficient accuracy, to have recourse to
approximation schemes. We shall consider below the framework provided by the
effective field theory approach. A detailed account of this work can be found
in Sazdjian:2022kaf .

##  3 Effective field theory approach

According to the energy balance analysis, there is always a probability that a
compact tetraquark be formed from the action of the confining forces (actually
an infinite tower of such states). We assume that in first approximation,
because of the compactness of the state, the latter can be assimilated to a
pointlike object and treated as an elementary particle.

We are mainly interested in the case where the mass of the latter particle
lies below the lowest possible two-meson threshold. This might be the case
when the compact tetraquark is the ground state of the corresponding spectrum.

The two mesons are designated by $M_{1}$ and $M_{2}$, with masses $m_{1}$ and
$m_{2}$, respectively. The bare mass of the tetraquark is $m_{t0}$.

Quark flavors and spins will be ignored, as not playing a fundamental role
here.

Because of the initial general structure of the tetraquark, containing two
meson clusters, the pointlike tetraquark has necessarily a coupling to the
mesons $M_{1}$ and $M_{2}$. The coupling is assumed scalar and is designated
by $g^{\prime}$; the latter is dimensionless and is factored by the mass term
$(m_{1}+m_{2})$.

The coupling $g^{\prime}$ generates, through meson loops, radiative
corrections inside the tetraquark propagator, modifying the parameters of the
bare propagator. In particular, the bare mass $m_{t0}$ gets changed into a
physical mass $m_{t}$. This is graphically represented in Fig. 2.

Figure 2: The full tetraquark propagator, including the radiative corrections,
coming from the coupling to the two meson clusters.

Mutual interactions of mesons are neglected (in first approximation); they
play a less important role than the direct tetraquark-two-meson coupling.

The full tetraquark propagator becomes

$D_{t}(s)=\frac{i}{s-m_{t0}^{2}+i(m_{1}+m_{2})^{2}g^{\prime 2}J(s)},$ (5)

where $s$ stands for $p^{2}$ and $J$ is the standard loop function. The
divergence of $J$ is absorbed by the bare mass term, yielding a renormalized
mass $m_{t1}$:

$m_{t1}^{2}=m_{t0}^{2}-i(m_{1}+m_{2})^{2}g^{\prime 2}J^{\mathrm{div}}.$ (6)

The renormalized tetraquark propagator is now

$D_{t}(s)=\frac{i}{s-m_{t1}^{2}+i(m_{1}+m_{2})^{2}g^{\prime
2}J^{\mathrm{r}}(s)},$ (7)

where $J^{\mathrm{r}}$ is the renormalized finite part of $J$. Notice that
$g^{\prime}$ does not undergo any renormalization. The mass term $m_{t1}$ does
not yet represent the physical mass of the tetraquark. The latter is
determined from the pole position of the propagator.

We stick here to the case of heavy quarks and heavy mesons, treating the
problem in its nonrelativistic limit, referred to the two-meson threshold. The
nonrelativistic energy $E$ is introduced through the standard definition

$\sqrt{s}=(m_{1}+m_{2})+E.$ (8)

The nonrelativistic energy corresponding to the renormalized mass $m_{t1}$ of
the tetraquark is defined similarly:

$E_{t1}=m_{t1}-(m_{1}+m_{2}),\ \ \ \ \ E_{t1}<0.$ (9)

Since we are considering the case of stable tetraquarks (under the strong
interactions), the mass of the tetraquark is expected to lie below the two-
meson threshold. The physical nonrelativistic energy of the tetraquark will be
designated by $E_{t}$ and is also expected to be negative.

To simplify notations, we shall use henceforth the reduced dimensionless
energy variables $e$ through the definitions

$e\equiv\frac{E}{2m_{r}},\ \ \ \ e_{t1}\equiv\frac{E_{t1}}{2m_{r}},\ \ \ \
e_{t}\equiv\frac{E_{t}}{2m_{r}},\ \ \ \
m_{r}=\frac{m_{1}m_{2}}{(m_{1}+m_{2})}.$ (10)

The quantity $-e_{t}$ represents the nonrelativistic binding energy of the
tetraquark by reference to the two-meson threshold, although it is different
from the binding energy defined from the confining forces by reference to the
quark masses.

One finds for the nonrelativistic energy $e_{t}$ of the tetraquark the
equation

$-e_{t}+e_{t1}+\frac{g^{\prime 2}}{16\pi}\sqrt{-e_{t}}=0,$ (11)

whose solution is

$\sqrt{-e_{t}}=\frac{1}{2}\Big{[}-\frac{g^{\prime
2}}{16\pi}+\sqrt{\Big{(}\frac{g^{\prime
2}}{16\pi}\Big{)}^{2}-4e_{t1}}\Big{]}.$ (12)

Figure 3: Variation of the square-root of the binding energy as a function of
$g^{\prime 2}/(16\pi)$.

The binding energy $-e_{t}$ is a decreasing function of $g^{\prime 2}/(16\pi)$
and comes out smaller than $-e_{t1}$, reaching the value 0 when
$g^{\prime}\rightarrow\infty$. This is represented in Fig. 3.

We find a very rapid decrease of the binding energy. With ordinary values of
$g^{\prime 2}/(16\pi)$, of the order of 1, the binding energy decreases by a
factor of 1/100. The state takes the appearance of a shallow bound state.

##  4 Compositeness

The comparison of the molecular and compact schemes is reminiscent of a
general problem, already raised in the past in the case of the deuteron state,
denoted under the term of compositeness Weinberg:1965zz .

Weinberg has shown that this question can receive, in the nonrelativistic
limit, a precise and model-independent answer, by relating the probability of
a state as being elementary (or compact) to observable quantities, represented
by the scattering length and the effective range of the two constituents of
the molecular scheme in the $S$-wave state of their scattering amplitude.
Designating by $Z$ this probability, one has the following relations for the
scattering length $a$ and the effective range $r_{e}$, adapted to the
tetraquark problem:

$a=\frac{2(1-Z)}{(2-Z)}R,\ \ \ \ \ r_{e}=-\frac{Z}{(1-Z)}R,\ \ \ \ \
R=(-2m_{r}E_{t})^{-1/2},$ (13)

where $R$ is the radius of the bound state, $m_{r}$ the reduced mass of the
two-meson system and $E_{t}$ the tetraquark nonrelativistic energy.

Investigations about the compositeness criterion and its applicability to
various tetraquark candidates, as well as to ordinary hadrons, can be found in
Baru:2003qq ; Cleven:2011gp ; Hanhart:2011jz ; Hyodo:2011qc ; Aceti:2012dd ;
Sekihara:2014kya ; Guo:2015daa ; Kang:2016ezb ; Meissner:2015mza ;
Oller:2017alp ; Guo:2020vmu ; Esposito:2021vhu ; Li:2021cue ; Baru:2021ldu ;
Kinugawa:2021ykv ; Song:2022yvz .

The contribution of the tetraquark state, in the $s$-channel, to the two-meson
elastic scattering amplitude is obtained by inserting the tetraquark
propagator between two tetraquark-two-meson couplings (cf. Fig. 4).

Figure 4: The tetraquark contribution, in the $s$-channel, to the two-meson
elastic scattering amplitude.

From the scattering length and the effective range one obtains $Z$:

$Z=\frac{\sqrt{-e_{t}}}{\sqrt{-e_{t}}+\frac{1}{2}\frac{g^{\prime 2}}{16\pi}}.$
(14)

When $Z=1$, we have the case of a pure compact tetraquark, while for $Z=0$, we
have a pure molecular state. The variation of $Z$ with respect to $g^{\prime
2}/(16\pi)$ is represented in Fig. 5.

Figure 5: Variation of $Z$ as a function of $g^{\prime 2}/(16\pi)$.

Like the binding energy, $Z$ is a rapidly decreasing function of $g^{\prime
2}/(16\pi)$. For values of the latter of the order of 1, $Z$ is very close to
zero. This confirms the interpretation that the coupling of the compact
tetraquark to its internal mesonic clusters qualitatively deforms its initial
structure, bringing it into a form closer to a molecular-type state.

##  5 Conclusion

A compact tetraquark, formed from the confining forces acting between quarks
and gluons, rapidly evolves, under the influence of the clustering phenomenon,
towards a molecular-type state. The origin of the state is, however, of
compact nature: nowhere, in the present model, did we consider direct
interactions between mesons.

The experimental test for this phenomenon is provided by the measure of the
elementariness coefficient $Z$. Pure molecular states are characterized by the
value $Z=0$. A value $Z\neq 0$, beyond uncertainties, reflects the existence
of an original compact tetraquark. Many tetraquark candidates fall in this
category.

Because of the shrinking of the tetraquark binding energy to values close to
zero, shallowness of many bound states may receive a natural explanation from
the above mechanism.

The present study can also be extended to the case of resonances and to the
case when direct meson-meson interactions are incorporated.

Funding. This research has received financial support from the EU research and
innovation programme Horizon 2020, under Grant agreement No. 824093, and from
the joint CNRS/RFBR Grant No. PRC Russia/19-52-15022.

## References

  * (1) S. Choi et al. (Belle), Phys. Rev. Lett. 91, 262001 (2003), hep-ex/0309032
  * (2) B. Aubert et al. (BaBar), Phys. Rev. Lett. 90, 242001 (2003), hep-ex/0304021
  * (3) D. Besson et al. (CLEO), Phys. Rev. D 68, 032002 (2003), [Erratum: Phys. Rev. D 75 (2007) 119908], hep-ex/0305100
  * (4) B. Aubert et al. (BaBar), Phys. Rev. Lett. 95, 142001 (2005), hep-ex/0506081
  * (5) M. Ablikim et al. (BESIII), Phys. Rev. Lett. 110, 252001 (2013), 1303.5949
  * (6) Z. Liu et al. (Belle), Phys. Rev. Lett. 110, 252002 (2013), [Erratum: Phys. Rev. Lett. 111 (2013) 019901], 1304.0121
  * (7) M. Ablikim et al. (BESIII), Phys. Rev. Lett. 111, 242001 (2013), 1309.1896
  * (8) R. Aaij et al. (LHCb), Phys. Rev. Lett. 112, 222002 (2014), 1404.1903
  * (9) R. Aaij et al. (LHCb), Phys. Rev. Lett. 115, 072001 (2015), 1507.03414
  * (10) R. Aaij et al. (LHCb), Sci. Bull. 65, 1983 (2020), 2006.16957
  * (11) R. Aaij et al. (LHCb), Phys. Rev. D 102, 112003 (2020), 2009.00026
  * (12) L. Wu, Nucl. Part. Phys. Proc. 312-317, 15287 (2021), 2012.15473
  * (13) E. Spadaro Norella (LHCb), PoS LHCP2021, 051 (2021)
  * (14) R. Aaij et al. (LHCb), Nature Commun. 13, 3351 (2022), 2109.01056
  * (15) T. Gershon (LHCb) (2022), 2206.15233
  * (16) R.L. Jaffe, Nucl. Phys. A 804, 25 (2008)
  * (17) W. Lucha, D. Melikhov, H. Sazdjian, Phys. Rev. D 100, 094017 (2019), 1908.10164
  * (18) M.B. Voloshin, L.B. Okun, JETP Lett. 23, 333 (1976), [Pisma Zh. Eksp. Teor. Fiz. 23, 369 (1976)]
  * (19) A. De Rújula, H. Georgi, S.L. Glashow, Phys. Rev. Lett. 38, 317 (1977)
  * (20) N.A. Törnqvist, Z. Phys. C 61, 525 (1994), hep-ph/9310247
  * (21) R. Jaffe, F. Wilczek, Phys. Rev. Lett. 91, 232003 (2003), hep-ph/0307341
  * (22) E. Shuryak, I. Zahed, Phys. Lett. B 589, 21 (2004), hep-ph/0310270
  * (23) L. Maiani, F. Piccinini, A.D. Polosa, V. Riquer, Phys. Rev. D 71, 014028 (2005), hep-ph/0412098
  * (24) L. Maiani, V. Riquer, F. Piccinini, A.D. Polosa, Phys. Rev. D 72, 031502 (2005), hep-ph/0507062
  * (25) W. Lucha, D. Melikhov, H. Sazdjian, Prog. Part. Nucl. Phys. 120, 103867 (2021), 2102.02542
  * (26) H.G. Dosch, Phys. Rev. D 28, 412 (1983)
  * (27) C. Alexandrou, G. Koutsou, Phys. Rev. D 71, 014504 (2005), hep-lat/0407005
  * (28) F. Okiharu, H. Suganuma, T.T. Takahashi, Phys. Rev. Lett. 94, 192001 (2005), hep-lat/0407001
  * (29) F. Okiharu, H. Suganuma, T.T. Takahashi, Phys. Rev. D 72, 014505 (2005), hep-lat/0412012
  * (30) H. Suganuma, T. Iritani, F. Okiharu, T.T. Takahashi, A. Yamamoto, AIP Conf. Proc. 1388, 195 (2011), 1103.4015
  * (31) N. Cardoso, M. Cardoso, P. Bicudo, Phys. Rev. D 84, 054508 (2011), 1107.1355
  * (32) P. Bicudo, M. Cardoso, O. Oliveira, P.J. Silva, Phys. Rev. D 96, 074508 (2017), 1702.07789
  * (33) H. Sazdjian, Symmetry 14, 515 (2022), 2202.01081
  * (34) S. Weinberg, Phys. Rev. 137, B672 (1965)
  * (35) V. Baru, J. Haidenbauer, C. Hanhart, Y. Kalashnikova, A.E. Kudryavtsev, Phys. Lett. B 586, 53 (2004), hep-ph/0308129
  * (36) M. Cleven, F.K. Guo, C. Hanhart, U.G. Meissner, Eur. Phys. J. A 47, 120 (2011), 1107.0254
  * (37) C. Hanhart, Y.S. Kalashnikova, A.V. Nefediev, Eur. Phys. J. A 47, 101 (2011), 1106.1185
  * (38) T. Hyodo, D. Jido, A. Hosaka, Phys. Rev. C 85, 015201 (2012), 1108.5524
  * (39) F. Aceti, E. Oset, Phys. Rev. D 86, 014012 (2012), 1202.4607
  * (40) T. Sekihara, T. Hyodo, D. Jido, PTEP 2015, 063D04 (2015), 1411.2308
  * (41) Z.H. Guo, J.A. Oller, Phys. Rev. D 93, 096001 (2016), 1508.06400
  * (42) X.W. Kang, Z.H. Guo, J.A. Oller, Phys. Rev. D 94, 014012 (2016), 1603.05546
  * (43) U.G. Meissner, J.A. Oller, Phys. Lett. B 751, 59 (2015), 1507.07478
  * (44) J.A. Oller, Annals Phys. 396, 429 (2018), 1710.00991
  * (45) Z.H. Guo, J.A. Oller, Phys. Rev. D 103, 054021 (2021), 2012.11904
  * (46) A. Esposito, L. Maiani, A. Pilloni, A.D. Polosa, V. Riquer, Phys. Rev. D 105, L031503 (2022), 2108.11413
  * (47) Y. Li, F.K. Guo, J.Y. Pang, J.J. Wu, Phys. Rev. D 105, L071502 (2022), 2110.02766
  * (48) V. Baru, X.K. Dong, M.L. Du, A. Filin, F.K. Guo, C. Hanhart, A. Nefediev, J. Nieves, Q. Wang, Phys. Lett. B 833, 137290 (2022), 2110.07484
  * (49) T. Kinugawa, T. Hyodo, EPJ Web Conf. 262, 01019 (2022), 2112.00249
  * (50) J. Song, L.R. Dai, E. Oset, Eur. Phys. J. A 58, 133 (2022), 2201.04414