# Discrepancies in Epidemiological Modeling of Aggregated Heterogeneous Data

Anna L. Trella1 Peniel N. Argaw1,111Equal contribution Michelle M.
Li2,11footnotemark: 1 James A. Hay3 1Harvard John A. Paulson School Of
Engineering and Applied Sciences, 2Harvard Medical School, 3Harvard T.H. Chan
School of Public Health {annatrella, peniel<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Within epidemiological modeling, the majority of analyses assume a single
epidemic process for generating ground-truth data. However, this assumed data
generation process can be unrealistic, since data sources for epidemics are
often aggregated across geographic regions and communities. As a result,
state-of-the-art models for estimating epidemiological parameters, e.g.
transmission rates, can be inappropriate when faced with complex systems. Our
work empirically demonstrates some limitations of applying epidemiological
models to aggregated datasets. We generate three complex outbreak scenarios by
combining incidence curves from multiple epidemics that are independently
simulated via SEIR models with different sets of parameters. Using these
scenarios, we assess the robustness of a state-of-the-art Bayesian inference
method that estimates the epidemic trajectory from viral load surveillance
data. We evaluate two data-generating models within this Bayesian inference
framework: a simple exponential growth model and a highly flexible Gaussian
process prior model. Our results show that both models generate accurate
transmission rate estimates for the combined incidence curve at the cost of
generating biased estimates for each underlying epidemic, reflecting highly
heterogeneous underlying population dynamics. The exponential growth model,
while interpretable, is unable to capture the complexity of the underlying
epidemics. With sufficient surveillance data, the Gaussian process prior model
captures the shape of complex trajectories, but is imprecise for periods of
low data coverage. Thus, our results highlight the potential pitfalls of
neglecting complexity and heterogeneity in the data generation process, which
can mask underlying location- and population-specific epidemic dynamics.

††Our Python implementation of the Bayesian inference method using viral load
surveillance data is open source and available at
https://github.com/pargaw/pyCt

##  1 Introduction

Reliable testing and reporting of new cases play key roles in early detection
and management of an epidemic. For example, testing data are often used to
generate predictions regarding the dynamics of an epidemic Tolles and Luong
(2020) to guide interventions enforced by governments Brauner et al. (2021).
Data sources for epidemics are usually aggregated across geographic regions
and spatial scales Dong et al. (2020). However, aggregating datasets can mask
spatial heterogeneity in underlying incidence trends Chang et al. (2021). In
fact, it has been shown that individual-level data can better inform risk
models’ predictions Xu et al. (2020), especially for disadvantaged communities
Chang et al. (2021).

Still, due to limited access to data (e.g. during early detection stages),
epidemiologists must simulate epidemics to generate scenarios or predictions
for possible outcomes of interventions Tolles and Luong (2020). Because of the
complexities in disease and population dynamics, routinely used
epidemiological models (e.g. compartmental models) make assumptions that can
be oversimplifying Tolles and Luong (2020).

For the scope of this work, we adopt the recent Ct model to evaluate our
simulated complex outbreak scenarios. The Ct model estimates epidemiological
dynamics from viral load distributions, namely cycle threshold (Ct) values
from real-time PCR assays Hay et al. (2021). Ct values represent the viral
load in a person at a given point in time, and the distribution of detectable
viral loads under random cross-sectional surveillance has been shown to be
informative of the epidemic growth rate Hay et al. (2021). Specifically, low
Ct values ($\text{Ct}\leq 29$) correspond to high viral loads observed,
indicating recent epidemic growth, while high Ct values ($38\leq\text{Ct}\leq
40$; Ct values above 40 cannot be measured due to technical limitations of
real-time PCR assays) correspond to low viral loads observed, indicating
recent epidemic decline Hay et al. (2021). Still, the Ct model, like other
state-of-the-art (SoTA) models for estimating epidemiological dynamics,
assumes a single epidemic process that generates the data. We use the Ct model
as a demonstration—any epidemiological model can be applied instead—to
evaluate an epidemiological model’s robustness in differentiating complex
population structures and epidemiological dynamics.

Our evaluation process comprises of two main components: simulating three
scenarios in which two distinct populations—with different susceptibility,
infectious, incubation, and recovery rates—interact (Sections 2 and 3); and
measuring the fit of the Ct model on the combined populations’ viral load
distributions (Sections 2 and 4). We also empirically demonstrate the model’s
sensitivity to the choice of a posterior sampling method (Section 4).

Ultimately, our results highlight some of the challenges in implementing and
deploying epidemiological models. While it is clearly naive to assume a single
data-generating process for a dataset constructed by two underlying epidemic
processes, we hypothesize that estimates from an epidemiological model can
only reflect one of the epidemics and/or populations, depending on the stage
of the two underlying epidemics and how well the data represent the two
populations. For example, suppose we have an abundance of surveillance data
from population #1 that indicate an epidemic decline, and sparse surveillance
data from population #2, which is entering an epidemic growth phase. After
combining the datasets, we must be careful not to incorrectly conclude that
population #2 is also in the decline phase of the epidemic. Blindly
aggregating data can lead to flawed policy decisions. Thus, as shown through
this work, our recommendation is to evaluate the level of aggregation in the
data before drawing conclusions (e.g. Given only state-aggregated data, how
well-equipped are we to make county-level decisions?).

##  2 Methods

We begin with an overview of the models used to simulate and evaluate the
outbreak scenarios described in Section 3. We utilize these simulated data to
assess the accuracy of the inference approach in estimating the true epidemic
trajectory. We note that our simulation settings are set to broadly reflect a
SARS-CoV-2 outbreak, but our findings can be considered pathogen-nonspecific.

###  2.1 Simulations via SEIR model

We simulate our complex scenarios with a compartmental model, specifically the
Susceptible-Exposed-Infectious-Recovered (SEIR) model. We use a basic SEIR
model Ma (2020), where the population size $N=S+E+I+R$, and the transition
rates between the compartments $S$, $E$, $I$, and $R$ are defined by

$\displaystyle\frac{dS}{dt}$ $\displaystyle=-\beta SI$ (1)
$\displaystyle\frac{dE}{dt}$ $\displaystyle=\beta SI-\sigma E$ (2)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=\sigma E-\gamma I$ (3)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=\gamma I$ (4)

where $\beta$ is the transmission rate (or rate in which susceptible
individuals get exposed), $\sigma$ is the incubation rate (or rate at which
exposed individuals become infectious), and $\gamma$ is the recovery rate (or
rate at which infectious individuals recover).

Table 1: Default parameters for the simulated scenarios. We set the following default values for the parameters in our SEIR model. Parameter | Value
---|---
$N$ | 1,000,000
$R_{0}$ | 2.5
$\gamma$ | 1/8
$\sigma$ | 1/4
$S_{initial}$ | $0.9999N$
$E_{initial}$ | 0
$I_{initial}$ | $0.0001N$
$R_{initial}$ | 0
Table 2: Scaling factors for the simulated outbreak scenarios. For our simulations, we define scaling factors for the SEIR model parameters $N$, $R_{0}$, and $\gamma$, where $R_{0}=\beta/\gamma$. We only report the scaling factors for population #1 because the parameters are scaled relative to population #2, whose parameters are scaled by 1. Scenario | $N$ | $R_{0}$ | $\gamma$
---|---|---|---
Out-of-state travel (Section 3.2) | 1 | 2 | 1.5
Seasonal travel (Section 3.3) | 0.5 | 0.8 | 0.8
NPIs (Section 3.4) | 1 | 0.8 | 1

###  2.2 Simulations via Ct model

After generating epidemic incidence curves for each of the scenarios using the
SEIR model, we simulate observable Ct values under random cross-sectional
surveillance, assuming that infected individuals are tested at an unknown
point in their infection. The frequency and number of tests simulated are
specified in Table S1. To capture the fact that measured Ct values depend on
the age of an infection, $a$, we use a previously described model for modal
viral load (or Ct) kinetics, $C_{mode}(a)$, assuming that observed Ct values
follow a Gumbel distribution Hay et al. (2021)

$C(a)\sim Gumbel(C_{mode}(a),\sigma(a))$ (5)

where $\sigma(a)$ captures the observation that the Ct values exhibit less
variability in older infections. We also assume that the parameters underlying
the Ct model are known precisely based on previous estimates Hay et al.
(2021).

###  2.3 Single Cross-Section Model

We infer epidemic growth rates by fitting an exponential growth model given a
single cross-section of observed Ct values. Briefly, the aim is to infer the
recent growth rate of the epidemic using the distribution of Ct values
obtained at a single point in time. An exponential growth model is simple yet
interpretable because it has only one free parameter, which denotes the
direction and magnitude of the epidemic trajectory from a single cross-section
of data. We assume that the daily incidence over the $A_{max}$ days before the
testing day $t$ grows or declines exponentially with rate $\beta$ such that

$\pi_{t-a}=\pi_{0}e^{\beta(t-a)}$ (6)

where $\beta$ is the logarithm of the daily growth rate in incidence over the
days $t-A_{max}$ and $t-1$, $\pi_{0}$ is a nuisance parameter that can be
ignored, and we assume a prior $\beta\sim\mathcal{N}(0,0.25)$.

###  2.4 Multiple Cross-Section Model

We also infer epidemic growth rates by fitting a Gaussian process model given
multiple cross-sections. Gaussian Processes (GPs) Rasmussen (2003) are a
powerful non-parametric framework for performing Bayesian inference. The
framework pairs a GP prior on daily incidence rate Xu et al. (2016) with a
custom Ct likelihood function. The flexibility of the GP prior allows for
practitioners to analyze multiple test days, as the priors do not assume
specific parameters. Our choices of a covariance function and hyperparameters
are manually specified by domain experts. Hay et al. (2021) uses a modified
version of the squared exponential (SE) kernel:

$k(x,x^{\prime})=\eta^{2}\exp\left[-\rho^{2}(x-x^{\prime})^{2}\right]$ (7)

with hyperparameters $\eta=1.5,\rho=0.03.$ Therefore, daily incidence is
defined as

$\pi_{t-a}=(1+\exp(-f))^{-1}$ (8)

where $f\sim\mathcal{N}(\mathbf{0},\mathbf{K})$ and $\mathbf{K}$ is a
covariance matrix generated by $k$. For sampling, we use Hamilton Monte Carlo
Duane et al. (1987) with a No-U-Turn sampler Hoffman and Gelman (2014) because
it is a more efficient MCMC method than Metropolis-Hastings Robert and Casella
(1999), the sampler used in the implementation by Hay et al. (2021).

##  3 Simulated Scenarios

(

(a)

(

(b)

(

(c)

(

(d)  Figure 1: Incidence rates across individual and aggregated populations
for each outbreak scenario. (a) The baseline scenario consists of only one
population (Section 3.1), whose incidence rates are shown in blue. Each of the
(b) out-of-state travel (Section 3.2), (c) seasonal travel (Section 3.3), and
(d) NPIs (Section 3.4) scenarios has two panels: on the left are two incidence
rate curves from populations #1 (blue) and #2 (green); and on the right is the
incidence rate curve of the aggregated dataset. The red dots indicate times
when testing is performed.

a) b) c) d)

We next introduce our scenarios for evaluating the robustness of the Ct model.
We create the scenarios by independently simulating outbreaks in two
populations, and then combining their simulated Ct values observed through
surveillance testing (Figure S1a). In each scenario, we arbitrarily define
scaling factors for SEIR model parameters (Table 2) to generate outbreaks with
different transmission dynamics (Figure 1a). The choices of parameters and
scaling factors are designed to qualitatively demonstrate the desired
epidemiological dynamics within the following scenarios, and we note that
other parameter choices may give different dynamics but qualitatively similar
results.

###  3.1 Baseline

As a baseline, we simulate a simple scenario with only one population using
the default parameters in Tables 1 and S1.

###  3.2 Out-of-State Travel

We simulate a scenario in which population #1 consists of younger individuals
(e.g. college students) who travel into a community (population #2) from out-
of-state (e.g. for spring break). We assume that the rate of recovery
($\gamma$) in the incoming younger individuals is higher on average than that
of the community, broadly reflecting the epidemiological situation for SARS-
CoV-2, where younger individuals are less likely to experience symptomatic
infections Davies et al. (2020). We also assume that the potential for
transmission ($R_{0}$) is higher in population #1 to reflect higher numbers
and intensity of contacts amongst younger individuals. Further, the Ct model
is built to incorporate variation across a population, including in age
structures Hay et al. (2021). Here, we assume that individuals in the two
populations follow the same viral kinetics model (Table S1), and only vary the
SEIR model parameters in Table 2. We combine the populations by adding their
Ct values at corresponding times $t\in[0,250]$ (Figure 1a).

###  3.3 Seasonal Travel

We simulate the dynamics of an epidemic in a community (e.g. college town) as
a result of seasonal travel (e.g. students in town for the semester or leaving
during the breaks). Population #1 consists of individuals found in the town
during the breaks; population #2 consists of those present during the
semester. In this scenario, the effect of population density can be modeled
using a density-dependent compartment model (Table 2). The transmission rate
is defined as $\beta=pC$, where $p$ is the probability of infection per
contact, and $C$ is the number of contacts per unit time. We assume that
population #1 has a lower transmission potential since the number of contacts
per unit time decreases ($R_{0}$ is scaled by $0.8$). As the younger
population leaves the town, population size and recovery rate also decrease.
We combine the two populations by specifying the outbreaks in populations #1
and #2 to occur at $t\in[0,250]$ and $t\in[251,500]$, respectively (Figure
1a).

###  3.4 Non-Pharmaceutical Interventions

We simulate a scenario in which non-pharmaceutical interventions (NPIs) (e.g.
curfews, mask mandates) are enforced in one population (population #1) but not
in another (population #2). We assume that the enforcement of NPIs results in
a lower transmission rate ($R_{0}$) (Table 2). We combine the populations as
described in Section 3.2 (Figure 1a).

##  4 Results

The exponential growth and GP prior models output average growth rates for a
set of 8 defined time points (or 16 in the seasonal travel scenario by the
exponential growth model). At each time point $t$, we show the predicted
growth rate (namely, $\beta$ for the exponential model and $f$ for the GP
prior model) and the true growth rate derived from the true incidence values.
To evaluate the robustness of the two models, we compare the predicted growth
rates against the true growth rates within a 95% credible interval.
Predictions are generated using only observed Ct values as input data (i.e. Ct
$<40$), so only data from infected individuals are included in the inference.
Thus, estimates for the epidemic trajectory give relative, rather than
absolute, incidence rates since we do not include information on the
proportion of individuals who are infected.

###  4.1 Baseline

We first compare both methods’ performance on a single standard incidence rate
trajectory generated from a SEIR model (Section 3.1). Both models are able to
estimate the correct direction of the epidemic, and the GP prior model can
even capture the full shape of the trajectory as we increase the number of
cross-sections (Figures 2aa and 3aa). However, despite the simple trajectory,
the exponential growth model has difficulty predicting growth rates near the
inflection point of the curve (Figure S2), and the GP prior model needs—in our
case—at least two cross-sections to capture the shape (Figure 4).

###  4.2 Out-of-State Travel

In the out-of-state travel scenario, our exponential model is able to fit the
combined incidence curve within all testing times except for when there is
high variability across the two populations (namely, $t\in\\{55,69\\}$)
(Figure S3). This is due to the fact that the aggregated dataset’s incidence
rate is the mean of the two populations’ incidence rates (Figure 1aa). Since
the rate of decrease in population #1 is greater than the rate of increase in
population #2 at $t\in\\{55,69\\}$, the predicted growth rate shows a sharper
decline than the true growth rate (Figures 2ab and S3). Interestingly, the
mean predicted growth rate for $t=83$ is approximately two times greater than
the true growth rate, although it is still within the 95% credible interval
(Figure 2ab).

The GP prior model is able to fully capture the shape of the trajectory of
population #2, but not that of the trajectory of population #1 (Figure 3ab).
The peak of population #1’s epidemic occurs before the first testing time
point (Figure 1ab), so the GP prior model is only able to infer the epidemic
decline in population #1 (Figure 3ab). As a result of this limitation in
surveillance data, even the GP prior model cannot accurately fit the initial
incidence curve of population #1 despite its overconfident uncertainty
intervals (Figure 3ab).

###  4.3 Seasonal Travel

As there is no overlap between the two populations in the seasonal travel
scenario (Figure 1ab), there is very little variability in the aggregated
dataset. Thus, the exponential model seems to be robust in most of the sampled
times (Figures 2ac and S4). The GP prior model is able to capture the full
trajectory, including both peaks, as well because their testing time points
are densely scattered (Figure 3ac). We also note that the GP prior model is
able to handle the aleatoric uncertainty well, as shown by the large
uncertainty intervals between times $t\in[175,250]$, when no surveillance data
is available.

###  4.4 Non-Pharmaceutical Interventions

The exponential model performs slightly better in the NPIs scenario than in
the out-of-state travel scenario (Figures 2ad and S5). Although the
populations have opposing trajectories at $t\in\\{97,111\\}$, the fit of the
exponential model falls within our 95% credible boundary across all time
points.

The GP prior model is also able to capture the full trajectory when given
enough cross-sections (Figure 3ad). Similarly with our baseline experiment,
the GP prior model struggles when there are too few testing time points (i.e.
$t<50$). However, it is important to note that the GP prior model captures a
fair amount of uncertainty for the first cross-section.

(

(a)

(

(b)

(

(c)

(

(d)  Figure 2: Evaluation of the exponential growth model for each outbreak
scenario. For the (a) baseline, (b) out-of-state travel, (c) seasonal travel,
and (d) NPIs scenarios, we show the predicted growth rates (average taken
between $t-35$ and $t$) and the 95% credible interval for the 8 time points
$t$ (dot and vertical line), and their corresponding true growth rates (gray
curve). Note that the true growth rates exhibit high variability during time
periods with very few new infections (i.e. at the start or end of an
outbreak).

a) b) c) d)

(

(a)

(

(b)

(

(c)

(

(d)  Figure 3: Evaluation of the Gaussian process prior model for each
outbreak scenario. For the (a) baseline, (b) out-of-state travel, (c) seasonal
travel, and (d) NPIs scenarios, we show the posterior draws of $t$ (purple)
and the true incidence rate in the overall combined population (black).

a) b) c) d)

Figure 4: Limitations of the Gaussian process prior model. We use the baseline
scenario to demonstrate the limitations of the GP prior model. As the number
of cross-sections increase (the leftmost panel shows one cross-section; the
rightmost panel, four), the GP prior model’s ability to capture the trajectory
of the epidemic improves.

##  5 Conclusion

In this work, we demonstrate the limitations of SoTA epidemiological models in
capturing the heterogeneity in underlying aggregated datasets’ population
structures and epidemiological dynamics. While both the exponential growth and
GP prior models are able to infer the direction of the combined epidemic
across most time points, the models’ estimated trajectories tend to reflect
the population and/or epidemic with the highest abundance in surveillance data
(e.g. in the out-of-state travel scenario). Although the exponential growth
model seems to accurately estimate simple epidemic trajectories with limited
data at time points far from inflection points (e.g. in the baseline and
seasonal travel scenarios), it fails to capture the complexities of two
conflicting trajectories (e.g. in the out-of-state travel scenario). On the
other hand, the highly flexible GP prior model can only accurately fit
aggregated trajectories when provided sufficient cross-sections that are
densely scattered (e.g. unlike in the seasonal travel scenario).

Future work regarding the epidemiological models includes exploring the
effects of model selection (e.g. kernels $k$, hyperparameters $\eta$ and
$\rho$) for the GP prior model in order to address its limitations in
estimating aggregated datasets. Our Python implementation of the Ct,
exponential growth, and Gaussian process prior models are open source and
available. Given the potentially high-impact repercussions of biased and
inaccurate models (e.g. when being used to advise government interventions),
we must explore strategies for data aggregations such that models are both
robust and fair towards heterogeneous population structures and epidemic
dynamics. Moreover, epidemiological models’ performance could possibly provide
insights on the allocation of testing resources to ensure that adequate data
is obtained for fitting the models.

##  6 Broader Impact Statement

Our work is primarily motivated by the need to account for spatial
heterogeneity when estimating epidemiological parameters. This is particularly
important due to the use of epidemiological models to inform public health
policies (e.g. travel restrictions, mask-wearing) Brauner et al. (2021). Some
potential implications of an inaccurate model include (1) inadvertently
forcing the population to stay longer than necessary in quarantine or
lockdown, which could negatively affect the economy, overall psychological
health, and interpersonal relationships; or (2) releasing interventions when
incidence is still at the early stage of growth. Nevertheless, our models are
unable to make direct recommendations; we are simply highlighting the
discrepancies of epidemiological models.

As our results are drawn directly from simulated data, we must consider the
advantages and disadvantages of generating our own data rather than using
real-world data. One major advantage is that it enables us to preserve the
privacy of individuals El Emam et al. (2021). Moreover, because it is
relatively straightforward to manipulate the parameters of compartmental
models (e.g. SEIR model), we can simulate a wide range of desired scenarios
for investigation. While many studies have validated synthetic data as
adequate proxies for real-world data El Emam et al. (2021), there may still be
artificial biases introduced by the simulations.

Finally, our criticism of epidemiological models should be considered in the
context of current state-of-the-art models’ capabilities. Although standard
models often make oversimplifying assumptions, it is important to recognize
that modeling real world dynamics is highly complex and challenging. The
ability to capture even some of the key contributing factors of an epidemic
using a simple model is valuable, as it can enable further extrapolation to
more robust models.

## Acknowledgements

We are grateful for the guidance and support from Dr. Michael Mina and Dr.
Milind Tambe throughout our project. M.M.L. is supported by T32HG002295 from
the National Human Genome Research Institute and a National Science Foundation
Graduate Research Fellowship.

## References

  * Brauner et al. [2021] Jan M Brauner, Sören Mindermann, Mrinank Sharma, David Johnston, John Salvatier, Tomáš Gavenčiak, Anna B Stephenson, Gavin Leech, George Altman, Vladimir Mikulik, et al.  Inferring the effectiveness of government interventions against covid-19.  Science, 2021.
  * Chang et al. [2021] Serina Chang, Emma Pierson, Pang Wei Koh, Jaline Gerardin, Beth Redbird, David Grusky, and Jure Leskovec.  Mobility network models of covid-19 explain inequities and inform reopening.  Nature, 2021.
  * Davies et al. [2020] Nicholas G Davies, Petra Klepac, Yang Liu, Kiesha Prem, Mark Jit, and Rosalind M Eggo.  Age-dependent effects in the transmission and control of covid-19 epidemics.  Nature Medicine, 2020.
  * Dong et al. [2020] Ensheng Dong, Hongru Du, and Lauren Gardner.  An interactive web-based dashboard to track covid-19 in real time.  The Lancet Infectious Diseases, 2020.
  * Duane et al. [1987] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth.  Hybrid monte carlo.  Physics Letters B, 1987.
  * El Emam et al. [2021] Khaled El Emam, Lucy Mosquera, Elizabeth Jonker, and Harpreet Sood.  Evaluating the utility of synthetic covid-19 case data.  JAMIA Open, 2021.
  * Hay et al. [2021] James A. Hay, Lee Kennedy-Shaffer, Sanjat Kanjilal, Niall J. Lennon, Stacey B. Gabriel, Marc Lipsitch, and Michael J. Mina.  Estimating epidemiologic dynamics from cross-sectional viral load distributions.  Science, 2021.
  * Hoffman and Gelman [2014] Matthew D Hoffman and Andrew Gelman.  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  JMLR, 2014.
  * Ma [2020] Junling Ma.  Estimating epidemic exponential growth rate and basic reproduction number.  Infectious Disease Modelling, 2020.
  * Rasmussen [2003] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer School on Machine Learning, 2003.
  * Robert and Casella [1999] Christian P Robert and George Casella.  The metropolis—hastings algorithm.  In Monte Carlo Statistical Methods. 1999.
  * Tolles and Luong [2020] Juliana Tolles and ThaiBinh Luong.  Modeling epidemics with compartmental models.  JAMA, 2020.
  * Xu et al. [2016] Xiaoguang Xu, Theodore Kypraios, and Philip D O’Neill.  Bayesian non-parametric inference for stochastic epidemic models using gaussian processes.  Biostatistics, 2016.
  * Xu et al. [2020] Bo Xu, Bernardo Gutierrez, Sumiko Mekaru, Kara Sewalk, Lauren Goodwin, Alyssa Loskill, Emily L Cohn, Yulin Hswen, Sarah C Hill, Maria M Cobo, et al.  Epidemiological data from the covid-19 outbreak, real-time case information.  Scientific Data, 2020.

##  Appendix A Inference Models

The inference models from Hay et al. [2021] can be decomposed into two
components: (1) a custom Ct likelihood and (2) a prior on the growth rate
parameter(s). The custom likelihood for a single-cross section at time $t$ is
as follows:

$\mathcal{L}(X_{1},...,X_{n}|\pi_{t-A_{max}},...,\pi_{t-a})=\prod_{i=1}^{n}\bigg{[}\bigg{(}\sum_{a=1}^{A_{max}}p_{a}(X_{i})\phi_{a}\pi_{t-a}\bigg{)}^{\mathbf{I}(X_{i}\leq
C_{LOD})}\bigg{(}1-\sum_{a=1}^{A_{max}}\phi_{a}\pi_{t-a}\bigg{)}^{\mathbf{I}(X_{i}>C_{LOD})}\bigg{]}$

where, for testing day $t$, $\pi_{t-a}:=$ probability that a randomly-selected
individual in the population was infected on day $t-a$, $p_{a}(x):=$
probability that Ct value is $x$ for a test conducted $a$ days after
infection, and $\phi_{a}:=$ probability of a Ct value being detectable $a$
days after infection. Since we only consider detectable Ct values for a cross-
section, the likelihood simplifies into:

$\mathcal{L}(X_{1},...,X_{n}|\pi_{t-A_{max}},...,\pi_{t-a})=\frac{\prod_{i=1}^{n}(\sum_{a=1}^{A_{max}}p_{a}(X_{i})\phi_{a}\pi_{t-a}])}{(\sum_{a=1}^{A_{max}}\phi_{a}\pi_{t-a}])^{n}}$

where $n$ is the number of detectable Ct values in the cross-section. Since
the likelihood can lead to overflow values, we use log-likelihood. Using the
log-likelihood does not disrupt the optimization problem because log is
monotonic.

$\displaystyle\implies
l(X_{1},...,X_{n}|\pi_{t-A_{max}},...,\pi_{t-a})=\log\mathcal{L}$
$\displaystyle=\sum_{i=1}^{n}\log(\sum_{a=1}^{A_{max}}p_{a}(X_{i})\phi_{a}\pi_{t-a})-n\log(\sum_{a=1}^{A_{max}}\phi_{a}\pi_{t-a})$

###  A.1 Single Cross-Section Model

The single-cross section model assigns a normal prior on a single growth rate
parameter $\beta\sim\mathcal{N}(0,0.25)$. We assume a fixed growth rate for
$A_{max}=35$ days. Daily incidence within that interval is
$\pi_{t-a}=\exp(\beta(t-a))$. With both components, we obtain a posterior
distribution:

$p(\beta|X_{1},...,X_{n})\propto l(X_{1},...,X_{n}|\beta)p(\beta)$

We sample from the posterior and obtain draws of $\beta$. We then use these
$\beta$’s to construct a posterior predictive, which we compare against the
true incidence rates. Note that when we compare against incidence rates, we
need to modify the likelihood function such that
$\pi=[\pi_{t-a}/\sum_{a=1}^{A_{max}}\pi_{t-a}]$ reflects an incidence rate.

###  A.2 Multiple Cross-Section Model

For the multiple cross section model, we leverage the function space view of
Gaussian process priors Rasmussen [2003]. Instead of sampling a parameter, we
directly sample functions that represent the incidence rate curve. A GP prior
is defined entirely by its kernel function, which are usually specified
directly by domain experts. Accounting for multiple cross sections is a
natural extension from the single cross section model and is given by:

$\mathcal{L}(X_{1}^{t_{1}},...,X_{n_{1}}^{t_{1}},...,X_{1}^{t_{T}},...,X_{n_{T}}^{t_{T}}|\pi_{t-A_{max}^{(T)}},...,\pi_{t_{T}-1})=\prod_{j=1}^{T}\bigg{[}\frac{\prod_{i=1}^{n_{j}}(\sum_{a=1}^{A_{max}^{(j)}}p_{a}(X_{i}^{t_{j}})\phi_{a}\pi_{t_{j}-a})}{(\sum_{a=1}^{A_{max}^{(j)}}\phi_{a}\pi_{t_{j}-a})^{n_{j}}}\bigg{]}$

The log likelihood is derived below:

$\log\mathcal{L}(X_{1}^{t_{1}},...,X_{n_{1}}^{t_{1}},...,X_{1}^{t_{T}},...,X_{n_{T}}^{t_{T}}|\pi_{t-A_{max}^{(T)}},...,\pi_{t_{T}-1})=\sum_{j=1}^{T}\bigg{[}\sum_{i=1}^{n_{j}}\bigg{[}\log[(\sum_{a=1}^{A_{max}^{(j)}}p_{a}(X_{i}^{t_{j}})\phi_{a}\pi_{t_{j}-a})]-{n_{j}}\log(\sum_{a=1}^{A_{max}^{(j)}}\phi_{a}\pi_{t_{j}-a})\bigg{]}\bigg{]}$

Again, we need to modify the likelihood function to replace
$\pi=[\pi_{t-a}/\sum_{a=1}^{A_{max}}\pi_{t-a}]$ such that it generates
incidence rates.

##  Appendix B Supplementary Tables

Table S1: Testing times and counts from Ct simulations. We present the number of Ct values collected at each testing time for all scenarios. Scenario | $t=55$ | $t=69$ | $t=83$ | $t=97$ | $t=111$ | $t=125$ | $t=139$ | $t=153$
---|---|---|---|---|---|---|---|---
Baseline | 49 | 172 | 466 | 777 | 612 | 342 | 140 | 50
Out-of-state travel | 1148 | 466 | 523 | 710 | 616 | 334 | 140 | 61
Seasonal travel | 6 | 9 | 27 | 71 | 120 | 240 | 402 | 463
NPIs | 52 | 221 | 553 | 883 | 959 | 787 | 573 | 390
| $t=305$ | $t=319$ | $t=333$ | $t=347$ | $t=361$ | $t=375$ | $t=389$ | $t=403$
Seasonal travel | 36 | 149 | 412 | 737 | 619 | 328 | 142 | 50

##  Appendix C Supplementary Figures

(

(a)

(

(b)

(

(c)

(

(d)  Figure S1: Ct values of simulated data. For the (a) baseline, (b) out-of-
state travel, (c) seasonal travel, and (d) NPIs scenarios, we show the
distribution of Ct values across the combined populations at each observed
time point.

a) b) c) d)

Figure S2: Evaluation of the exponential growth model on the baseline
scenario. For each of the 8 time points, we show the predicted growth rates
$\beta$ from the exponential growth model (blue) and the true incidence rate
(black). The value for $\beta$ shown in each figure legend is the posterior
mean estimate. The model is able to estimate the correct direction of the
epidemic; however, despite the simple trajectory, it has difficulty predicting
growth rates near the inflection point of the curve. Figure S3: Evaluation of
the exponential growth model on the out-of-state travel scenario. For each of
the 8 time points, we show the predicted growth rates $\beta$ from the
exponential growth model (blue) and the true incidence rate (black). The value
for $\beta$ shown in each figure legend is the posterior mean estimate. The
model is able to correctly fit the combined incidence curve within all testing
times except at $t\in\\{55,69,111\\}$, of which $t\in\\{55,69\\}$ are
locations of high variability across the populations, and $t=111$ is an
inflection point. Figure S4: Evaluation of the exponential growth model on the
seasonal travel scenario. For each of the 16 time points, we show the
predicted growth rates $\beta$ from the exponential growth model (blue) and
the true incidence rate (black). The value for $\beta$ shown in each figure
legend is the posterior mean estimate. As there is no overlap between the two
populations in this scenario, we observe very little variability in the
aggregated dataset. The exponential growth model is robust in all sampled
times except at $t=361$, the inflection point of the curve. Figure S5:
Evaluation of the exponential growth model on the NPIs scenario. For each of
the 8 time points, we show the predicted growth rates $\beta$ from the
exponential growth model (blue) and the true incidence rate (black). The value
for $\beta$ shown in each figure legend is the posterior mean estimate. The
model performs relatively well in this scenario. Although the populations have
opposing trajectories at $t\in\\{97,111\\}$, the fit of the exponential model
falls within the 95% credible interval across all time points.


# Strong decays of newly observed $D_{sJ}$ states in a constituent quark model
with effective Lagrangians

Xian-Hui Zhong1,3111E-mail<EMAIL_ADDRESS>and Qiang Zhao2,3222E-mail:
<EMAIL_ADDRESS>1) Department of Physics, Hunan Normal University, and Key
Laboratory of Low-Dimensional Quantum Structures $\&$ Quantum Control of
Ministry of Education, Changsha 410081, P.R. China  2) Institute of High
Energy Physics, Chinese Academy of Sciences, Beijing 100049, P.R. China  3)
Theoretical Physics Center for Science Facilities, CAS, Beijing 100049, P.R.
China

###### Abstract

The strong decay properties of the newly observed states $D_{sJ}(3040)$,
$D_{sJ}(2860)$ and $D_{sJ}(2710)$ are studied in a constituent quark model
with quark-meson effective Lagrangians. We find that the $D_{sJ}(3040)$ could
be identified as the low mass physical state $|2{P_{1}}\rangle_{L}$
($J^{P}=1^{+}$) from the $D_{s}(2^{1}P_{1})$-$D_{s}(2^{3}P_{1})$ mixing. The
$D_{sJ}(2710)$ is likely to be the low-mass mixed state $|(SD)\rangle_{L}$ via
the $1^{3}D_{1}$-$2^{3}S_{1}$ mixing. In our model, the $D_{sJ}(2860)$ cannot
be assigned to any single state with a narrow width and compatible partial
widths to $DK$ and $D^{*}K$. Thus, we investigate a two-state scenario as
proposed in the literature. In our model, one resonance is likely to be the
$1^{3}D_{3}$ ($J^{P}=3^{-}$), which mainly decays into $DK$. The other
resonance seems to be the $|1{D_{2}}^{\prime}\rangle_{H}$, i.e. the high-mass
state in the $1^{1}D_{2}$-$1^{3}D_{2}$ mixing with $J^{P}=2^{-}$, of which the
$D^{*}K$ channel is its key decay mode. We also discuss implications arising
from these assignments and give predictions for their partner states such as
$|(SD)^{\prime}\rangle_{H}$, $|2{P^{\prime}_{1}}\rangle_{H}$, $2^{3}P_{0}$ and
$2^{3}P_{2}$, which could be helpful for the search for these new states in
future experiment.

###### pacs:

12.39.Fe, 12.39.Jh, 13.25.Ft, 13.25.Hw

##  I Introduction

Experimental progress on the study of $D$ and $D_{s}$ states in the past few
years provides a great opportunity for theory development. Recently a new
broad resonance $D_{sJ}(3040)$ with a mass of $(3044\pm
8_{\mathrm{stat}}(^{+30}_{-5})_{\mathrm{syst}})$ MeV and a width of
$\Gamma=(239\pm 35_{\mathrm{stat}}(^{+46}_{-42})_{\mathrm{syst}})$ MeV is
reported in the $D^{*}K$ channel Aubert:2009 . Apart from the $D_{sJ}(3040)$
another two states $D_{sJ}(2710)$ and $D_{sJ}(2860)$, which were observed by
BABAR and Belle two years ago Aubert:2006mh ; jb:2007aa , are also examined.
Their branching ratio fractions between $D^{*}K$ and $DK$ are measured
Aubert:2009 ,

$\displaystyle\frac{D_{sJ}(2710)^{+}\rightarrow
D^{*}K}{D_{sJ}(2710)^{+}\rightarrow DK}=0.91\pm 0.13_{\mathrm{stat}}\pm
0.12_{\mathrm{syst}},$ (1) $\displaystyle\frac{D_{sJ}(2860)^{+}\rightarrow
D^{*}K}{D_{sJ}(2860)^{+}\rightarrow DK}=1.10\pm 0.15_{\mathrm{stat}}\pm
0.19_{\mathrm{syst}}.$ (2)

These new observations stimulate great interest in the understanding of their
nature and strong coupling properties in theory. Different theoretical
approaches for the study of the strong coupling properties of heavy-light
mesons can be found in the literature, such as the heavy quark effective field
theory approach (HQEFT) Nowak:1992um ; Deandrea:1999pa ; Deandrea:1998uz ;
Dai:1997df ; Eichten:1993ub ; Dai:1998ve ; Hiorth:2002pp ; Fajfer:2006hi ;
Fajfer:2008zz ; Kamenik:2008zza ; Bardeen:1993ae ; Bardeen:2003kt ;
Colangelo:2004vu ; Colangelo:2007ds ; Fazio , QCD sum rules Colangelo:1995ph ;
Dai:1997df ; Wang:2008tm , ${}^{3}P_{0}$ model Close:2006gr ; Close:2005se ;
Zhang:2006yj ; Wei:2006wa ; Liuxiang , and chiral quark model Pierro .

In this work, we present an analysis of these $D_{s}$ states in a constituent
quark model with effective Lagrangians for the quark-meson couplings, and try
to clarify the following issues: (i) To gain information about the structure
of the newly observed state $D_{sJ}(3040)$ according to its strong decay
properties. (ii) With the new data for the $D_{sJ}(2710)$ and $D_{sJ}(2860)$,
we reanalyze their strong decays and examine their structures again. The
quantum numbers of these two states remain controversial. The $D_{sJ}(2710)$
is identified as a state of $J^{P}=1^{-}$ in $B$ decays jb:2007aa , while it
is explained by various models as the $2^{3}S_{1}$, $1^{3}D_{1}$, the
admixtures of $2^{3}S_{1}$-$1^{3}D_{1}$, molecular structure, or tetraquark
state Zhang:2006yj ; Wei:2006wa ; Close:2006gr ; Colangelo:2007ds ;
Wang:2007nfa . There are also a lot of solutions proposed for the
$D_{sJ}(2860)$. The assignments of $1^{3}D_{3}$ or $2^{3}P_{0}$ have been
discussed in Refs. Zhang:2006yj ; Close:2006gr ; vanBeveren:2006st ;
Wei:2006wa ; Colangelo:2006rq ; Koponen:2007nr ; Zhong:2008kd . A recent
comment by Ref. eef:2009  suggests a two-state structure for the
$D_{sJ}(2860)$ in order to understand the controversial aspects arising from
its strong decays. (iii) The quark-model assignment of these states will
result in implications of their partner multiplets. We discuss some of those
relevant states, for which the experimental observations would be able to
clarify some of those theoretical and experimental issues.

By treating the light mesons (pseudoscalar and vector mesons) as effective
fields, we introduce constituent-quark-meson couplings to describe the charmed
meson strong decays into a charmed meson plus a light pseudoscalar or vector
meson in the final state. The quark-pseudoscalar-meson coupling is given by
the chiral quark model at the leading order as proposed by Manohar and Georgi
Manohar:1983md . Its application to pseudoscalar meson photoproduction in the
quark model turns out to be promising and many low-energy phenomena can be
highlighted in such a framework Li:1994cy ; Li:1997gda ; qk3 ; Zhong:2007fx ;
Zhong:2009 . In particular, the axial current conservation allows one to
extract the axial coupling in terms of the meson decay constant and a form
factor arising from the microscopic quark model wavefunctions Riska:2000gd .
With an effective quark-vector-meson coupling, one can also extract the vector
couplings in a similar way zhao:1998fn ; zhao-k1 ; zhao-kstar ; Riska:2000gd .

A natural extension of this picture is to apply this effective Lagrangian
approach to heavy-light meson strong decays involving light pseudoscalar or
vector mesons, which would be a good place to examine the validity of the
light axial and vector fields in such a transition. On the one hand, the
quark-meson coupling is the same as that defined in meson photoproduction
which is proportional to the meson decay constant. On the other hand, the
heavy-light meson in the initial and final state would provide information
about the coupling form factor and can be calculated in the quark model
framework. Thus, one can study the heavy-light meson strong decays by
combining dynamical information from meson photoproduction off nucleons.

The paper is organized as follows. In Sec. II, a brief review of the quark-
meson effective Lagrangian approach is given. The numerical results are
presented and discussed in Sec. III. Finally, a summary is given in Sec. IV.

##  II framework

Figure 1:  Diagrams (A) and (B) stand for the quark-meson couplings in meson-
baryon interactions and light-meson production in heavy-light meson strong
decays, respectively.

In Fig. 1, we illustrate the similarity of the quark-meson coupling in meson-
baryon and light-meson production in heavy-light meson strong decays. It
should be pointed out that since the flavor symmetries beyond the SU(3) are
badly broken, the contributions from transitions of treating the final-state
heavy-light meson as an effective field are strongly suppressed. We thus can
neglect those contributions safely in our approach. An early study of the
charmed meson strong decays can be found in Ref. Zhong:2008kd .

In the chiral quark model Manohar:1983md , the low energy quark-pseudoscalar-
meson interactions in the SU(3) flavor basis are described by the effective
Lagrangian Li:1997gda ; qk3 ; qk4 ; Zhong:2007fx ; Zhong:2009

${\cal
L}_{Pqq}=\sum_{j}\frac{1}{f_{m}}\bar{\psi}_{j}\gamma^{j}_{\mu}\gamma^{j}_{5}\psi_{j}\partial^{\mu}\phi_{m}.$
(3)

where $\psi_{j}$ represents the $j$-th quark field in the hadron, and
$\phi_{m}$ is the pseudoscalar meson field.

The effective Lagrangian for quark-vector-meson interactions in the SU(3)
flavor basis is zhao:1998fn ; zhao-k1 ; zhao-kstar

${\cal
L}_{Vqq}=\sum_{j}\bar{\psi}_{j}(a\gamma^{j}_{\mu}+\frac{ib}{2m_{j}}\sigma_{\mu\nu}q^{\nu})V^{\mu}\psi_{j}\
,$ (4)

where $V^{\mu}$ represents the vector meson field with four-vector moment $q$.
Parameters $a$ and $b$ denote the vector and tensor coupling strength,
respectively.

As follows, we provide the quark-pseudoscalar and quark-vector-meson coupling
operators in a non-relativistic form Li:1997gda ; qk3 ; qk4 ; Zhong:2007fx ;
Zhong:2009 ; zhao:1998fn ; zhao-k1 ; zhao-kstar . Considering light meson
emission in a heavy-light meson strong decays, the effective quark-
pseudoscalar-meson coupling operator in the center-of-mass (c.m.) system of
the initial meson is

$\displaystyle
H_{m}=\sum_{j}\left[-\left(1+\frac{\omega_{m}}{E_{f}+M_{f}}\right)\mbox{\boldmath$\sigma$\unboldmath}_{j}\cdot\textbf{q}+\frac{\omega_{m}}{2\mu_{q}}\mbox{\boldmath$\sigma$\unboldmath}_{j}\cdot\textbf{p}_{j}\right]I_{j}\varphi_{m}\
.$ (5)

In a case that a light vector meson is emitted, the transition operators for
producing a transversely or longitudinally polarized vector meson are as
follows

$\displaystyle
H_{m}^{T}=\sum_{j}\left\\{i\frac{b^{\prime}}{2m_{q}}\mbox{\boldmath$\sigma$\unboldmath}_{j}\cdot(\mathbf{q}\times\mathbf{\epsilon})+\frac{a}{2\mu_{q}}\mathbf{p}_{j}\cdot\mathbf{\epsilon}\right\\}I_{j}\varphi_{m},$
(6)

and

$\displaystyle H_{m}^{L}=\sum_{j}\frac{aM_{v}}{|\mathbf{q}|}I_{j}\varphi_{m}\
.$ (7)

In the above three equations, q and $\omega_{m}$ are the three-vector momentum
and energy of the final-state light meson, respectively. $\textbf{p}_{j}$ is
the internal momentum operator of the $j$-th quark in the heavy-light meson
rest frame. $\mbox{\boldmath$\sigma$\unboldmath}_{j}$ is the spin operator on
the $j$-th quark of the heavy-light system and $\mu_{q}$ is a reduced mass
given by $1/\mu_{q}=1/m_{j}+1/m^{\prime}_{j}$ with $m_{j}$ and
$m^{\prime}_{j}$ for the masses of the $j$-th quark in the initial and final
mesons, respectively. Here, the $j$-th quark is referred to the active quark
involved at the quark-meson coupling vertex. $M_{v}$ is the mass of the
emitted vector meson. The plane wave part of the emitted light meson is
$\varphi_{m}=e^{-i\textbf{q}\cdot\textbf{r}_{j}}$, and $I_{j}$ is the flavor
operator defined for the transitions in the SU(3) flavor space Zhong:2007gp ;
Li:1997gda ; qk3 ; qk4 ; Zhong:2009 ; Zhong:2007fx ; Zhong:2008kd ;
zhao:1998fn . Parameters $a$ and $b$ are the vector and tensor coupling
strengths of the quark-vector-meson couplings, respectively. Studies of vector
meson photoproduction zhao-k1 ; zhao-kstar ; JLab-Kstar ; elsa-kstar  suggest
that $a=g_{\omega qq}=g_{\rho qq}\simeq-3$ and $b^{\prime}\equiv b-a\simeq 5$.
Because of vector current conservation, one has $a=g_{\rho NN}=g_{\omega
NN}/3$ Riska:2000gd ; zhao:1998fn ; zhao-kstar .

The heavy-light meson wavefunctions have been given in Ref. Zhong:2008kd , and
some of the decay amplitudes have also been deduced there. In the charmed
meson decays, the SU(4) flavor symmetry is broken. Thus, the charm quark is
treated as a spectator and the transition amplitude is proportional to the
final-state light meson decay constant associated with a form factor arising
from the convolution of the initial and final-state charmed meson
wavefunctions.

In the calculation, the standard quark model parameters are adopted. Namely,
we set $m_{u}=m_{d}=330$ MeV, $m_{s}=450$ MeV, and $m_{c}=1700$ MeV for the
constituent quark masses. The harmonic oscillator parameter $\beta$ is usually
adopted in the range of (0.4–0.5) GeV, in this work we take it as $\beta=0.45$
GeV. The decay constants for $K$ and $\eta$ mesons are $f_{K}=f_{\eta}=160$
MeV. As shown in Refs. Zhong:2007gp ; Zhong:2008kd , the flavor symmetry
breaking will lead to corrections to the quark-pseudoscalar-meson coupling
vertex, for which an additional global parameter $\delta$ is introduced. Here,
we fix its value the same as that in Refs. Zhong:2008kd ; Zhong:2007gp , i.e.
$\delta=0.557$. For the quark-vector-meson coupling strength which still
suffers relatively large uncertainties, we adopt the values extracted from
vector meson photoproduction as mentioned earlier, i.e. $a\simeq-3$ and
$b^{\prime}\simeq 5$. The masses of the mesons used in the calculations are
adopted from the PDG PDG .

Justification of the non-relativistic formulation is not obvious for the light
quark sector in the heavy-light meson transitions. This is similar to the case
of a non-relativistic quark model for baryons, where the results would rely on
the experimental data to tell how far they deviate from reality. Treating the
light meson as a chiral field somehow assumes that the light meson is produced
at short distance, and the spectators (i.e. the two spectator quarks inside a
baryon or the heavy quark in the heavy-light meson transitions) do not respond
to the internal structure of the light meson. Instead, the propagation of the
light quark pair would feel the hadronic environment from the convolution of
initial and final-state heavy-light mesons. Such an implicated assumption
means that only the processes with relatively small momentum transfers between
the light quarks inside the light meson would dominantly contribute to the
transition matrix element. This empirically supports the validity of the non-
relativistic formulation as a leading order approximation.

##  III results and discussions

###  III.1 $D_{sJ}(3040)$

The $D_{sJ}(3040)$ is observed in the $D^{*}K$ mode, while there is no sign of
$D_{sJ}(3040)\rightarrow DK$ in experiment Aubert:2009 . This allows its
quantum number to be $J^{P}=0^{-}$, $1^{+}$, $2^{-}$ etc. The $J^{P}=0^{-}$
state $2^{1}S_{0}$ seems not a good candidate since its predicted mass, $\sim
2.7$ GeV Pierro ; Godfrey:1985xj ; Close:2006gr ; Ebert:2009ua , is much less
than $3.04$ GeV. The predicted masses of $J^{P}=1^{+}$ and $2^{-}$ are close
to $3.04$ GeV. We hence discuss these two possibilities for the $D_{sJ}(3040)$
in this work.

First, we considered it as the $J^{P}=1^{+}$ states $2^{1}P_{1}$ and
$2^{3}P_{1}$. These two states also can decay into $D^{*}K$, $DK^{*}$,
$D^{*}_{s}\eta$, $D_{s}\phi$, $D_{0}(2400)K$, $D_{1}(2430)K$, $D_{1}(2420)K$,
$D_{2}(2460)K$, $D_{s}(2317)\eta$, $D_{s}(2460)\eta$. With a mass of 3.04 GeV,
we calculate their decay widths, which are listed in Tab. 1. From the table,
it is found that the decay width of $2^{1}P_{1}$ and $2^{3}P_{1}$ are $\sim
115$ MeV and $\sim 93$ MeV, respectively, which are too small to compare with
the data, although the decay mode, dominated by the $D^{*}K$, is consistent
with the observation Aubert:2009 . Thus, the $D_{sJ}(3040)$ may not be
considered as pure $2^{1}P_{1}$ or $2^{3}P_{1}$ state.

Figure 2:  (Color online) The partial decay widths and total width of
$|2P_{1}\rangle_{L}$ with a mass of 3040 MeV as functions of mixing angle
$\phi$. The data are from BABAR Aubert:2009 .

Since the heavy-light mesons are not charge conjugation eigenstates, state
mixing between spin $\mathbf{S}=0$ and $\mathbf{S}=1$ states with the same
$J^{P}$ can occur via the spin-orbit interactions Godfrey:1986wj ;
Close:2005se ; Swanson . The physical states with $J^{P}=1^{-}$ can then be
described as

$\displaystyle|2P_{1}\rangle_{L}=+\cos(\phi)|2^{1}P_{1}\rangle+\sin(\phi)|2^{3}P_{1}\rangle,$
(8)
$\displaystyle|2{P_{1}}^{\prime}\rangle_{H}=-\sin(\phi)|2^{1}P_{1}\rangle+\cos(\phi)|2^{3}P_{1}\rangle\
,$ (9)

where the subscripts $L$ and $H$ stand for the low mass and high mass of the
physical states after the mixing.

Usually, the low mass state has a broad width while the high mass state has a
narrow width. We set the mass of $|2P_{1}\rangle_{L}$ with $3.04$ GeV, and
plot its decay width as a function of the mixing angle $\phi$, which is shown
in Fig. 2. It shows that when the mixing angle is in the range
$\phi\simeq-(40\pm 12)^{\circ}$, the total decay width, $\Gamma=(162\sim 170)$
MeV, is in the range of the experimental data (close to the lower limit of the
data) Aubert:2009 . The mixing angle predicted here is consistent with the
result $\phi\simeq-55^{\circ}$ in the heavy quark limit Godfrey:1986wj ;
Close:2005se ; Swanson . The $D^{*}K$ governs the decays of
$|2P_{1}\rangle_{L}$, while the $DK$ channel is forbidden. This is also in
agreement with the observations. These results suggest that the $D_{sJ}(3040)$
favors the $|2P_{1}\rangle_{L}$ classification.

Apart from the $D^{*}K$ mode, the $D_{1}(2430)K$, $D_{2}(2460)K$,
$D_{0}(2400)K$, $DK^{*}$, and $D_{s}^{*}\eta$ are also important in the decays
of $|2P_{1}\rangle_{L}$ as shown by Fig. 2. In particular, the partial widths
of $D_{1}(2420)K$, $D_{s}(2317)\eta$ and $D^{*}K^{*}$ turn out to be sizable.
A search for those channels would be useful for clarifying the property of the
$D_{sJ}(3040)$. With the mixing angle $\phi\simeq-55^{\circ}$, the relative
decay ratios among those decay channels are
$D^{*}K:D_{1}(2430)K:D_{2}(2460)K:D_{0}(2400)K:D_{1}(2420)K:DK^{*}:D_{s}(2317)\eta:D_{s}^{*}\eta:D^{*}K^{*}=78:17:19:13:4:8:4:11:2$.

Figure 3:  (Color online) The partial decay widths and total width of
$|2P_{1}\rangle_{L}$ as functions of mass. The data are from BABAR Aubert:2009
.

Since the mass of $D_{s}(3040)$ still has a large uncertainty, it may bring
uncertainties to the theoretical predictions on the decay widths. To
investigate this effect, we plot the decay widths as a function of the mass in
Fig. 3 with the mixing angle fixed at $\phi=-50^{\circ}$. It shows that the
mass uncertainty gives rise to an uncertainty of about $\sim 70$ MeV in the
total decay width. The predicted widths are much closer to the central value
of the data with the increasing mass. The sensitivity of different decay modes
to the mass can also be seen clearly in the plot.

Table 1: The decay widths (MeV) for the $D_{sJ}(3040)$ as $1^{1}D_{2}$, $1^{3}D_{2}$, $2^{1}P_{1}$ and $2^{3}P_{1}$ candidates. | $D^{*}K$ | $DK^{*}$ | $D^{*}K^{*}$ | $D^{*}_{s}\eta$ | D(2430)K | D(2420)K | $D_{s}(2460)\eta$ | $D_{s}\phi$ | D(2400)K | D(2460)K | $D_{s}(2317)\eta$ | total
---|---|---|---|---|---|---|---|---|---|---|---|---
$1^{1}D_{2}$ | 197 | 27 | 2 | 25 | 3 | 2 | 0.4 | 4 | 3 | 345 | 4 | 608
$1^{3}D_{2}$ | 256 | 21 | 33 | 34 | 1 | 18 | 0.01 | 0.05 | 3.4 | 512 | 1.6 | 879
$2^{1}P_{1}$ | 44 | 9 | 0.3 | 5.5 | 0.02 | 0.01 | $7.5\times 10^{-5}$ | 0.1 | 33 | 12 | 11 | 115
$2^{3}P_{1}$ | 41 | 2 | 2.5 | 7.5 | 24 | 7 | 0.5 | 0.002 | 0.09 | 9 | 0.06 | 93

Finally, we discuss the possibilities of $D_{sJ}(3040)$ as a $J^{P}=2^{-}$
candidate. There are two states, $1^{1}D_{2}$ and $1^{3}D_{2}$, with
$J^{P}=2^{-}$. If $1^{1}D_{2}$ and $1^{3}D_{2}$ have a mass of 3.04 GeV, they
can decay into the following channels, $D^{*}K$, $DK^{*}$, $D^{*}K^{*}$,
$D^{*}_{s}\eta$, $D_{s}\phi$, $D_{0}(2400)K$, $D_{1}(2430)K$, $D_{1}(2420)K$,
$D_{2}(2460)K$, $D_{s}(2317)\eta$, and $D_{s}(2460)\eta$. We calculate these
partial decay widths and list the results in Tab. 1. It shows that $D^{*}K$
and $D_{2}(2460)K$ are the two main decay channels. The total widths for both
$1^{1}D_{2}$ and $1^{3}D_{2}$ are very broad, i.e. $\Gamma\sim 608$ MeV and
$\sim 879$ MeV, respectively. They are too large to compare with the data
$\Gamma=(239\pm 35)$ MeV Aubert:2009 . Nevertheless, it shows that the
admixtures between $1^{1}D_{2}$ and $1^{3}D_{2}$ are unable to give a
reasonable explanation of the decay properties of $D_{sJ}(3040)$ as well.
Thus, the $D_{sJ}(3040)$ as a $J^{P}=2^{-}$ candidate is not favored.

In brief, the $D_{sJ}(3040)$ seems to favor a $|2P_{1}\rangle_{L}$ state with
$J^{P}=1^{+}$, which is an admixture of $2^{1}P_{1}$ and $2^{3}P_{1}$ with a
mixing angle $\phi\simeq-(40\pm 12)^{\circ}$. Our conclusion is in agreement
with that of a ${}^{3}P_{0}$ model analysis Liuxiang . The semiclassical flux
tube model Ailin  and relativistic quark model Ebert:2009ua  mass calculations
also support this picture.

###  III.2 $D_{sJ}(2710)$

The $D_{sJ}(2710)$ was first reported by BABAR Aubert:2006mh , and its quantum
number $J^{P}=1^{-}$ was determined by Belle jb:2007aa . Recently, the decay
ratios of the $D_{sJ}(2710)$ have also been reported Aubert:2009 , which is
very useful for understanding its nature. According to the classification of
the quark model, only two states $2^{3}S_{1}$ and $1^{3}D_{1}$ with the
quantum number $J^{P}=1^{-}$ are located around the mass range $(2.7\sim 2.8)$
GeV. This state is studied by various models, e.g. as a $2^{3}S_{1}$ state
Ebert:2009ua ; Fazio , $1^{3}D_{1}$ state Wei:2006wa , or admixture of
$2^{3}S_{1}$-$1^{3}D_{1}$ Close:2006gr . It should mention that in our
previous work Zhong:2008kd  an error occurred in the partial decay amplitude
of $1^{3}D_{1}\rightarrow DK$, which led to a rather small width for the
assignment of the admixture of $2^{3}S_{1}$-$1^{3}D_{1}$. Here we correct the
formulation and reanalyze the mixing scenario for the $D_{sJ}(2710)$.

Table 2: The decay widths (MeV) for the $D_{sJ}(2710)$ as $1^{3}D_{1}$ and $2^{3}S_{1}$ candidates. | $D^{0}K^{+}$ | $D^{+}K^{0}$ | $D^{*+}K^{0}$ | $D^{*0}K^{+}$ | $D_{s}\eta$ | $D^{*}_{s}\eta$ | total | $\Gamma(D^{*}K)/\Gamma(DK)$
---|---|---|---|---|---|---|---|---
$1^{3}D_{1}$ | 75 | 73.6 | 17.8 | 18.5 | 14 | 0.9 | 200 | 0.24
$2^{3}S_{1}$ | 5.4 | 5.6 | 9.0 | 9.1 | 1.7 | 0.7 | 31 | 1.65

We first assign the $D_{sJ}(2710)$ as the $2^{3}S_{1}$ and $1^{3}D_{1}$ states
and calculate its decay widths. The results are listed in Tab. 2,
respectively. For the assignment of the $2^{3}S_{1}$ state, the total decay
width and the decay branching ratio fraction between $D^{*}K$ and $DK$
channels are

$\displaystyle\Gamma\simeq 31~{}\mathrm{MeV},\ \
\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 1.65.$ (10)

It shows that the predicted width $\Gamma\simeq 31$ MeV is too narrow to
compare with the data, and the predicted decay ratio $D^{*}K/DK\simeq 1.67$ is
much larger than the measurement $D^{*}K/DK\simeq 0.91\pm 0.13\pm 0.12$
Aubert:2009 . The calculations of Ref. Zhang:2006yj  also tend to give a small
width $\Gamma\simeq 32$ MeV for the $2^{3}S_{1}$ configuration. The predicted
branching ratio fraction is also inconsistent with the observations
Aubert:2009 . In Ref. Close:2006gr , it is also found that the large branching
ratio fraction $D^{*}K/DK\simeq 3.55$ does not support the $D_{sJ}(2710)$ as a
pure $2^{3}S_{1}$ state.

On the other hand, if the $D_{sJ}(2710)$ is considered as a $1^{3}D_{1}$
state, the decay width and branching ratio fraction will be

$\displaystyle\Gamma\simeq 200~{}\mathrm{MeV},\ \
\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 0.24.$ (11)

In this case, the branching ratio fraction $\Gamma(D^{*}K)/\Gamma(DK)\simeq
0.24$ is too small though the decay width $\Gamma\simeq 200$ MeV is roughly
consistent with the upper limit of the data jb:2007aa ; Aubert:2009 . These
results suggest that either $1^{3}D_{1}$ or $2^{3}S_{1}$ is not a good
assignment for the $D_{sJ}(2710)$.

Thus, we consider the possibilities of the $D_{sJ}(2710)$ as a mixed state of
$2^{3}S_{1}$-$1^{3}D_{1}$, for which the physical states can be expressed as
Close:2006gr

$\displaystyle|(SD)_{1}\rangle_{L}=+\cos(\phi)|2^{3}S_{1}\rangle+\sin(\phi)|1^{3}D_{1}\rangle,$
(12)
$\displaystyle|(SD)^{\prime}_{1}\rangle_{H}=-\sin(\phi)|2^{3}S_{1}\rangle+\cos(\phi)|1^{3}D_{1}\rangle\
,$ (13)

where the physical partner in the mixing is included. Assuming that the low
mass state $|(SD)_{1}\rangle_{L}$ corresponds to the $D_{sJ}(2710)$
Close:2006gr , we plot the decay properties of $|(SD)_{1}\rangle_{L}$ as
functions of the mixing angle $\phi$ in Fig. 4. It shows that with the mixing
angle $\phi\simeq(-54\pm 7)^{\circ}$, the decay width and branching ratio
fraction are

$\displaystyle\Gamma\simeq(133\pm 22)~{}\mathrm{MeV},\ \
\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 0.91\mp 0.25\ ,$ (14)

which are in a good agreement with the data jb:2007aa ; Aubert:2009 .

Following this scheme, one can examine the high-mass partner
$|(SD)^{\prime}_{1}\rangle_{H}$, of which the expected mass is $\sim 2.81$ GeV
Close:2006gr . Taking into account the mass uncertainties of a region
$M\simeq(2.71\sim 2.88)$ GeV, we plot the mass-dependence of the partial and
total widths in Fig. 5. It shows that the $|(SD)^{\prime}_{1}\rangle_{H}$ also
has a broad width $\sim(120\pm 10)$ MeV, and the $DK$ channel is dominant over
others. In contrast, the partial width of $D_{s}\eta$ is also sizable, while
the $D^{*}_{s}\eta$ width is negligible. Around $M=2.81$ GeV, the predicted
branching ratio fractions are

$\displaystyle\frac{\Gamma(D_{s}\eta)}{\Gamma(DK)}\simeq
0.15,~{}\frac{\Gamma(D^{*}K)}{\Gamma(D_{s}\eta)}\simeq 0.06.$ (15)

The above mixing scheme is consistent with Ref. Close:2006gr  for the low-mass
state while the predicted suppression of the $D^{*}K$ decay mode is different
from that of Ref. Close:2006gr . In Ref. Close:2006gr  a very broad high-mass
state is predicted and would dominantly decay into both $DK$ and $D^{*}K$. In
our scheme, the predicted decay width for $|(SD)^{\prime}_{1}\rangle_{H}$ is
$\sim(120\pm 10)$ MeV. As a consequence, one would expect that it should
appear in the $DK$ spectrum similar to the $D_{s}(2710)$ signal. Taking into
account the still undetermined mass for $|(SD)^{\prime}_{1}\rangle_{H}$, one
possible explanation would be that the $|(SD)^{\prime}_{1}\rangle_{H}$ mass
may be larger than $M\simeq 2.88$ GeV. If so, its total width would be larger
than we estimated above and become much broader, thus, cannot be easily
identified in the present $DK$ spectrum. Interestingly, a recent study of the
$D_{s}$ spectrum suggests a larger mass for the $1^{3}D_{1}$ state
Ebert:2009ua .

It should be noted that different methods seem to lead to different
conclusions on the $D_{sJ}(2710)$ state. In Refs. Colangelo:2007ds ; Fazio ,
both the decay width and branching ratio fraction of the $D_{sJ}(2710)$ as the
$2^{3}S_{1}$ assignment can be well explained. In Ref. Ebert:2009ua , the mass
calculation also suggests that the $D_{sJ}(2710)$ is $2^{3}S_{1}$. However,
the recent study of a ${}^{3}P_{0}$ model tends to conclude that the
$D_{sJ}(2710)$ is a mixture of $2^{3}S_{1}$ and $1^{3}D_{1}$ Li:2009qu .
Therefore, additional information for $D_{sJ}(2710)\to D_{s}\eta$ and
$D_{s}^{*}\eta$, as well as a search for the $|(SD)^{\prime}_{1}\rangle_{H}$
partner in experiment would be useful for understanding the property of the
$D_{sJ}(2710)$.

Figure 4:  (Color online) The partial decay widths, total width, and the decay
branching ratio fraction $\Gamma(D^{*}K)/\Gamma(D^{*}K)$ of
$|(SD)_{1}\rangle_{L}$ as functions of mass, respectively. The data are from
BABAR Aubert:2009 . Figure 5:  (Color online) The partial decay widths and
total width of $|(SD)^{\prime}_{1}\rangle_{H}$ as functions of mass.

###  III.3 $D_{sJ}(2860)$

The situation about the $D_{sJ}(2860)$ is still controversial and different
solutions have been proposed in the literature. In Ref. vanBeveren:2006st ,
the $D_{sJ}(2860)$ is assigned as a $J^{P}=0^{+}$ state. However, the recent
observation of $D_{sJ}(2860)\rightarrow D^{*}K$ does not support this picture.
It is also proposed to be a $J^{P}=3^{-}$ state Colangelo:2006rq ;
Zhang:2006yj ; Zhong:2008kd . However, although the decay width and decay mode
are consistent with the observation, the predicted ratio $D^{*}K/DK\simeq 0.4$
is too small to compare with the data $D^{*}K/DK\simeq 1.1$ Aubert:2009 .

Table 3: The decay widths (MeV) for the $D_{sJ}(2860)$ as $1^{3}D_{3}$, $2^{3}P_{2}$, and $1^{3}F_{2}$ candidates. | $D^{0}K^{+}$ | $D^{+}K^{0}$ | $D^{*+}K^{0}$ | $D^{*0}K^{+}$ | $D_{s}\eta$ | $D^{*}_{s}\eta$ | $DK^{*}$ | total | $\Gamma(D^{*}K)/\Gamma(DK)$
---|---|---|---|---|---|---|---|---|---
$1^{3}D_{3}$ | 12.3 | 11.8 | 5 | 4.7 | 1.7 | 0.3 | 0.2 | 36 | 0.40
$2^{3}P_{2}$ | 1.3 | 1.3 | 2.1 | 1.9 | 0.01 | 1.7 | 0.02 | 8 | 1.53
$1^{3}F_{2}$ | 21.9 | 21.4 | 0.1 | 0.1 | 5.5 | 0.02 | 0.005 | 49 | 0.005

Since the $D_{sJ}(2860)$ is observed in both $D^{*}K$ and $DK$ channels, the
allowed quantum numbers would be $1^{3}D_{3}$, $2^{3}P_{2}$ and $1^{3}F_{2}$.
We calculate the total and partial widths for these configurations and list
the results in Tab. 3.

More specifically, as the $1^{3}D_{3}$ state, the predicted width and
branching ratio fraction between the $D^{*}K$ and $DK$ channel are

$\displaystyle\Gamma\simeq
36~{}\mathrm{MeV},~{}\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 0.4.$ (16)

The predicted ratio $\Gamma(D^{*}K)/\Gamma(DK)$ differs from the measurement
$D^{*}K/DK\simeq 1.1$ Aubert:2009  at the level of three standard deviations,
although the decay width is in agreement with the data. Our predictions are
consistent with those of Refs. Colangelo:2006rq ; Fazio . It should be
mentioned that the QCD-motivated relativistic quark model can not well explain
the mass of $D_{sJ}(2860)$ if it is considered as the $1^{3}D_{3}$ state
Ebert:2009ua . This could be a signal indicating the chiral symmetry in
association with the heavy quark symmetry in the heavy-light meson
transitions.

As a candidate of the $2^{3}P_{2}$ state, the decay width and branching ratio
fraction of $D_{sJ}(2860)$ are

$\displaystyle\Gamma\simeq
8~{}\mathrm{MeV},~{}\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 1.53,$ (17)

where both the predicted width and ratio are inconsistent with the data. It is
interesting to mention that our predicted ratio agrees with the estimation of
Ref. eef:2009 .

If the $D_{sJ}(2860)$ is a $1^{3}F_{2}$ state, the predicted width and
branching ratio fraction are

$\displaystyle\Gamma\simeq
49~{}\mathrm{MeV},~{}\frac{\Gamma(D^{*}K)}{\Gamma(DK)}\simeq 0.005\ ,$ (18)

where the decay mode of $D^{*}K$ turns out to be negligible in comparison with
the $DK$ mode, and disagrees with the experimental observation.

Figure 6:  (Color online) The partial decay widths and total width of
$2^{3}P_{0}$ as functions of mass. Figure 7:  (Color online) The partial decay
widths and total width of $|1D_{2}\rangle_{H}$ with a mass of 2860 MeV as
functions of mixing angle $\phi$. The data are from BABAR Aubert:2009 ;
Aubert:2006mh .

It can be seen from the above analysis that a simple assignment of the
$D_{sJ}(2860)$ to be a pure $2^{3}P_{0}$, $1^{3}D_{3}$, $2^{3}P_{2}$ or
$1^{3}F_{2}$ cannot well explain the data. We also point out that the
$2^{3}P_{2}$ and $1^{3}F_{2}$ mixing is unable to overcome the problem either
because of the narrow width of the $2^{3}P_{2}$ state or small branching ratio
fraction $\Gamma(D^{*}K)/\Gamma(DK)\simeq 0.005$ of $1^{3}F_{2}$.

In Ref. eef:2009 , van Beveren and Rupp recently proposed an alternative
solution that there might exist two largely overlapping resonances at about
2.86 GeV, i.e. a radially excited tensor ($2^{+}$) and a scalar ($0^{+}$)
$c\bar{s}$ state. Following this two-state assumption, one would expect that
one state $D_{sJ_{1}}(2860)$ dominantly decays into $DK$, while the other one
$D_{sJ_{2}}(2860)$ dominantly decays into $D^{*}K$. Both states have a mass
around 2.86 GeV, and comparable width $\Gamma\sim 50$ MeV. This idea may shed
some light on the controversial issues. As follows, we shall investigate such
a possibility in our approach.

It shows that the decays of $2^{3}P_{0}$, $1^{3}D_{3}$ and $1^{3}F_{2}$ is
dominated by the $DK$ channel, while the decay of $1^{3}D_{2}$, $1^{1}D_{2}$,
$2^{3}P_{2}$ is dominated by the $D^{*}K$ channel. We shall identify which
states are more appropriate candidates in the two-state scenario.

First, we analyze the states dominated by $DK$ decays, i.e. $2^{3}P_{0}$,
$1^{3}D_{3}$ and $1^{3}F_{2}$. In Fig. 6 the total and partial decay widths
for the $2^{3}P_{0}$ state are revealed. It shows that the $2^{3}P_{0}$
possesses a broad decay width $\Gamma\simeq 115$ MeV at about 2.86 GeV, which
is inconsistent with the data. The $1^{3}F_{2}$ is not considered as a good
candidate of $D_{sJ_{1}}(2860)$ as well since its mass is expected to be
larger than 3.1 GeV Pierro ; Ebert:2009ua . Furthermore, our earlier analysis
suggests that the $D_{sJ}(3040)$ may favor a configuration of
$|2P_{1}\rangle_{L}$ such that the mass of the $1^{3}F_{2}$ should be larger
than the $P$ wave state $|2P_{1}\rangle_{L}$ as a consequence. In contrast, we
find that the $1^{3}D_{3}$ could be a good candidate for $D_{sJ_{1}}(2860)$
since it is dominated by the $DK$ decay mode and has a narrow width
$\Gamma\simeq 36$ MeV. The calculation results for the total and partial decay
widths have been listed in Tab. 3.

Candidates for the $D_{sJ_{2}}(2860)$ could be $1^{3}D_{2}$, $1^{1}D_{2}$, or
$2^{3}P_{2}$ which dominantly decay into $D^{*}K$. As discussed earlier in
this section and shown in Tab. 3, the $2^{3}P_{2}$ is not a good candidate
since its total width is too small to compare with the data. Nevertheless, its
expected mass should be larger than 2.86 GeV Pierro ; Ebert:2009ua .

If the $D_{sJ_{2}}(2860)$ is considered as pure $1^{3}D_{2}$ or $1^{1}D_{2}$
state, their decay widths would be $\Gamma\simeq 170$ MeV and $\Gamma\simeq
130$ MeV, respectively, which are inconsistent with the data as well. In fact,
the physical states should be the admixtures between $1^{3}D_{2}$ and
$1^{1}D_{2}$ due to the presence of the spin-orbit interactions Godfrey:1986wj
; Close:2005se ; Swanson . Thus, the mixed states can be expressed as

$\displaystyle|1D_{2}\rangle_{L}=+\cos(\phi)|1^{1}D_{2}\rangle+\sin(\phi)|1^{3}D_{2}\rangle,$
(19)
$\displaystyle|1{D_{2}}^{\prime}\rangle_{H}=-\sin(\phi)|1^{1}D_{2}\rangle+\cos(\phi)|1^{3}D_{2}\rangle\
,$ (20)

where the subscripts $L$ and $H$ denote the low-mass and high-mass state due
to the mixing. Usually, the $|1{D_{2}}^{\prime}\rangle_{H}$ have a narrow
width Godfrey:1986wj ; Close:2005se ; Swanson . We thus consider the
$|1{D_{2}}^{\prime}\rangle_{H}$ as the $D_{sJ_{2}}(2860)$ in the calculation.
In Fig. 7 the decay properties as a function of the mixing angle are plotted.
We see that around $\phi=-65^{\circ}$ or $\phi=-35^{\circ}$ the decay width is
$\Gamma\simeq 40$ MeV, which is compatible with the observation, and the decay
mode is dominated by the $D^{*}K$. With $\phi=-35^{\circ}$, the corresponding
decay branching ratio fractions are

$\displaystyle\frac{\Gamma(D^{*}K)}{\Gamma(D^{*}_{s}\eta)}\simeq
1.2,~{}\frac{\Gamma(D^{*}K)}{\Gamma(DK^{*})}\simeq 13\ ,$ (21)

which fit in the experimental data quite well. This result turns out to
support the $|1{D_{2}}^{\prime}\rangle_{H}$ to be a candidates of
$D_{sJ_{2}}(2860)$ in the two-state scenario. In the range of
$\phi=-65^{\circ}\sim-35^{\circ}$ the partial widths do not change drastically
with the mixing angle. In contrast, the suggested value is consistent with
that ($\phi=-50.7^{\circ}$) obtained in the heavy quark effective theory
Godfrey:1986wj ; Close:2005se ; Swanson ; Ebert:2009ua .

In brief, it seems likely that the abnormal property with the $D_{sJ}(2860)$
arises from two overlapping resonances with the same mass but different decay
modes. One is $1^{3}D_{3}$ and the other is $|1{D_{2}}^{\prime}\rangle_{H}$
from the $1^{3}D_{2}$ and $1^{1}D_{2}$ mixing. The $1^{3}D_{3}$ state mainly
decays into $DK$ and the $|1{D_{2}}^{\prime}\rangle_{H}$ into $D^{*}K$. With
these two largely overlapping resonances at about 2.86 GeV, we can understand
both the observed decay widths and branching ratio fractions of the
$D_{sJ}(2860)$. It shows that the $1^{3}D_{3}$ has a sizable partial width in
the $D_{s}\eta$ channel, while the $|1{D_{2}}^{\prime}\rangle_{H}\to
D^{*}_{s}\eta$ also turns out to be measurable. Further measurements of
$\Gamma(D^{*}_{s}\eta)/\Gamma(D^{*}K)$ and $\Gamma(DK)/\Gamma(D_{s}\eta)$ may
be able to distinguish the $1^{3}D_{3}$ and $|1{D_{2}}^{\prime}\rangle_{H}$
and test the two-state scenario in experiment.

Figure 8:  (Color online) The partial decay widths and total width of the
$|2P^{\prime}_{1}\rangle_{H}$ state as functions of mass.

###  III.4 $D_{sJ}(2|P^{\prime}_{1}\rangle_{H})$, $D_{sJ}(2^{3}P_{0})$ and
$D_{sJ}(2^{3}P_{2})$

In this subsection we discuss the implications of other states following the
consequence of the assignments for the $D_{sJ}(3040)$, $D_{sJ}(2710)$ and
$D_{sJ}(2860)$.

Since the $D_{sJ}(3040)$ seems to favor a $P$ wave with $J^{P}=1^{+}$
($|2P_{1}\rangle_{L}$), experimental evidences for the other $P$ waves,
$D_{sJ}(2|P^{\prime}_{1}\rangle_{H})$, $D_{sJ}(2^{3}P_{0})$ and
$D_{sJ}(2^{3}P_{2})$, would be important to establish the spectrum. In
particular, its high-mass partner $|2P^{\prime}_{1}\rangle_{H}$ should be
searched in experiments. Supposing that the $|2P^{\prime}_{1}\rangle_{H}$ has
a mass in the range of $(3.04\sim 3.2)$ GeV, we plot in Fig. 8 the decay
widths as functions of the mass with the mixing angle $\phi=-50^{\circ}$ fixed
by $D_{sJ}(3040)$. It shows that the $|2P^{\prime}_{1}\rangle_{H}$ width is
indeed relative narrower around $M=3.04$ GeV, although we should note that the
decay width increases fast with the increasing mass. The decay channels,
$D_{0}(2400)K$, $D_{2}(2460)K$ and $D_{s}(2317)\eta$, are predicted to be the
dominant ones, which can be investigated in experiments. In contrast, the
$D^{*}K$ channel plays a less important role in the decays.

We further study the $D_{sJ}(2^{3}P_{0})$ in detail here. The decay widths as
a function of the possible mass range $M=(2.8\sim 2.9)$ GeV are plotted in
Fig. 6. In this range the total decay width is $\Gamma\simeq(90\sim 140)$ MeV,
and increases with the increasing mass. It shows that the $DK$ channel
dominates its decays. Taking the mass of the $D_{sJ}(2^{3}P_{0})$ as
$(2.82\sim 2.84)$ GeV Close:2006gr ; Matsuki , the total width and branching
ratio fractions between $D_{s}\eta$ and $DK$ are

$\displaystyle\Gamma\simeq(101\pm
5)~{}\mathrm{MeV},~{}\frac{\Gamma(D_{s}\eta)}{\Gamma(DK)}\simeq 0.08.$ (22)

It should be pointed out that the decay properties of $2^{3}P_{0}$ are similar
to those of $|(SD)^{\prime}_{1}\rangle_{H}$ in the mass range $M<2.9$ GeV (see
Fig. 5 and Fig. 6). Both of them have comparable decay widths $\Gamma\sim 100$
MeV, and mainly decay into $DK$. To distinguish them from each other, the
measurements of their decay ratio $\Gamma(D_{s}\eta)/\Gamma(DK)$ are
important. We also note that a recent calculation suggests a larger mass of
$M\simeq 3.054$ GeV for $2^{3}P_{0}$ Ebert:2009ua . As a consequence of this
scenario, its total decay width would become much broader than we estimated
above. Thus, it may not be easily isolated in experiment.

Figure 9:  (Color online) The partial decay widths and total width of
$2^{3}P_{2}$ as functions of mass.

As discussed earlier the $D_{sJ}(2860)$ does not favor the assignment of
$2^{3}P_{2}$. Thus, we investigate its decay properties and implications of
experimental measurement. We also plot its total and partial decay widths as
functions of the mass in the possible range $M=(3.04\sim 3.2)$ GeV in Fig. 9.
If $2^{3}P_{2}$ has a mass larger than 3.04 GeV, decay channels, $DK$,
$D^{*}K$, $DK^{*}$, $D^{*}K^{*}$, $D^{*}_{s}\eta$, $D_{s}\phi$, $D_{s}\eta$,
$D_{1}(2430)K$, $D_{1}(2420)K$, $D_{2}(2460)K$, $D_{s}(2460)\eta$, will open
in which $D_{1}(2430)K$ and $DK$ channels are dominant. In Fig. 9, we do not
show the results for the $D^{*}_{s}\eta$, $D_{s}\phi$ and $D_{s}(2460)\eta$
channels since they are negligibly small ($<1$ MeV). If we adopt the mass
$\sim 3.15$ GeV as predicted by Refs. Matsuki ; Ebert:2009ua ; Pierro , the
predicted width is $\Gamma\simeq 140$ MeV, and the relative decay strengths
are
$DK:D^{*}K:D_{1}(2430)K:D_{1}(2420)K:D_{2}(2460)K:DK^{*}:D^{*}K^{*}:D_{s}\eta\simeq
41:9:50:13:11:6:7:4$. It suggests that the $DK$, $D_{1}(2430)K$,
$D_{1}(2420)K$ channels may be the optimal ones for searching for the
$D_{sJ}(2^{3}P_{2})$ state in experiment.

###  III.5 Sensitivity to the harmonic oscillator parameter

It should be mentioned that model-dependent feature of our model arises from
the simple treatment of harmonic oscillator potential for the heavy-light
quark system. Therefore, uncertainties with the theoretical results are
present in the choice of the quark model parameter values. The most important
parameter in our model should be the harmonic oscillator strength $\beta$,
which controls the size effect or coupling form factor from the convolution of
the heavy-light meson wavefunctions. The commonly adopted range of this
quantity is $\beta=(0.4\sim 0.5)$ GeV, and we apply $\beta=0.45$ GeV in the
above calculations.

In order to examine the sensitivity of the calculation results to $\beta$, we
plot the decay widths and ratios of $2^{3}S_{1}$, $1^{3}D_{3}$, mixed state
$|(SD)\rangle_{L}$ of $2^{3}S_{1}$-$1^{3}D_{1}$, and mixed state
$|2{P_{1}}\rangle_{L}$ of $2^{1}P_{1}$-$2^{3}P_{1}$ as a function of $\beta$
in Fig. 10. It shows that the decay widths of these excited $D_{s}$ states
exhibit some sensitivities to the parameter $\beta$. Within the range of
$\beta=(0.45\pm 0.05)$ GeV, about $30\%$ uncertainties of the decay widths
would be expected. This is a typical order of accuracy for the constituent
quark model, and can be regarded as reasonable.

The ratio $\Gamma(D^{*}K)/\Gamma(DK)$ appears to behave differently. For the
$2^{3}S_{1}$, the sensitivity of the ratio to $\beta$ is apparent. In
contrast, the ratios of $|(SD)\rangle_{L}$ and $1^{3}D_{3}$ are quite
insensitive to $\beta$. The ratio $\Gamma(D^{*}K)/\Gamma(DK)$ is not shown for
$|2{P_{1}}\rangle_{L}$ since its decay into $DK$ is forbidden.

In brief, although the harmonic oscillator parameter $\beta$ can bring some
uncertainties to the final results, within the range of $\beta=(0.4\sim 0.5)$
GeV, our major conclusions will still hold.

Figure 10:  (Color online) The decay widths (upper panel) and ratios (lower
panel) of different configuration assignments as a function of $\beta$: the
solid lines are for the mixed state $|(SD)\rangle_{L}$ of
$2^{3}S_{1}$-$1^{3}D_{1}$ with a mass of 2710 MeV and mixing angle
$-55^{\circ}$; the dotted lines are for $2^{3}S_{1}$ with a mass of 2710 MeV;
the dot-dashed lines for $1^{3}D_{3}$ with a mass of 2860 MeV; and the dashed
line in the upper panel is for the mixed state $|2{P_{1}}\rangle_{L}$ of
$2^{1}P_{1}$-$2^{3}P_{1}$ with a mass of 3040 MeV and mixing angle
$-50^{\circ}$.

##  IV summary

In this work we investigate the strong decays of several newly observed
charmed mesons in a constituent quark model with effective Lagrangians for the
quark-meson interactions. The decay amplitudes are extracted for light
pseudoscalar meson or vector meson productions via axial or vector current
conservation between the quark-level and hadronic level couplings. The quark-
meson couplings can then be determined by independent measurements such as
meson photoproduction and meson-baryon scatterings.

We find that the new state $D_{sJ}(3040)$ can be identified as the low mass
physical state $|2{P_{1}}\rangle_{L}$ from the
$D_{s}(2^{1}P_{1})$-$D_{s}(2^{3}P_{1})$ mixing with a mixing angle
$\phi\simeq-(40\pm 12)^{\circ}$. Further experimental search for decay modes
of $D_{1}(2430)K$, $D_{2}(2460)K$, $D_{0}(2400)K$, $DK^{*}$, and
$D_{s}^{*}\eta$ should be able to disentangle its property and test our model
predictions.

The $D_{sJ}(2710)$ seems to favor a low mass physical state $|(SD)\rangle_{L}$
from the $2^{3}S_{1}$-$1^{3}D_{1}$ mixing with a mixing angle
$\phi\simeq(-54\pm 7)^{\circ}$. Both the ratio and width are in a good
agreement with the data. The decay properties of its heavy partner
$|(SD)^{\prime}_{1}\rangle_{H}$ are also discussed. It has a broad width
$\Gamma\simeq(110\sim 140)$ MeV at the $2.8$ GeV mass region, and dominated by
the $DK$ mode. We also point out that the $|(SD)^{\prime}_{1}\rangle_{H}$
state may be searched in the $DK$ spectrum as the $D_{sJ}(2710)$ if its mass
is $\sim 2.8$ GeV. Whether the present data have contained its signal could be
a crucial criteria for various model predictions.

The $D_{sJ}(2860)$ cannot be easily explained by a single configuration of
$2^{3}P_{0}$, $2^{3}P_{2}$, $1^{3}F_{2}$ or $1^{3}D_{3}$. To overcome this
problem we follow the proposal of a two-state picture by Ref. eef:2009  and
assume that two narrow resonances may have been observed around 2.86 GeV with
a width $\Gamma\simeq(40\sim 50)$ MeV. It shows that one resonance seems to be
the $1^{3}D_{3}$, which mainly decays into $DK$. The other resonance could be
the $|1{D_{2}}^{\prime}\rangle_{H}$, which is the high-mass state from the
$1^{1}D_{2}$-$1^{3}D_{2}$ mixing, and dominantly decays into $D^{*}K$. Further
theoretical and experimental efforts are needed to disentangle the mysterious
properties about this state.

We also study the implications arising from the assignments for those observed
resonances, e.g. their partner states in the mixing. In particular, if the
$D_{sJ}(3040)$ is indeed a $P$-wave state $|2{P_{1}}\rangle_{L}$, the other
three $P$-wave states $D_{sJ}(|2P^{\prime}_{1}\rangle_{H})$,
$D_{sJ}(2^{3}P_{0})$ and $D_{sJ}(2^{3}P_{2})$ may also have measurable effects
in experiment. Their strong decay properties are predicted, which could be
useful for future experimental studies.

##  Acknowledgements

This work is supported, in part, by the National Natural Science Foundation of
China (Grants 10675131 and 10775145), Chinese Academy of Sciences
(KJCX3-SYW-N2), and Ministry of Science and Technology of China
(2009CB825200).

## References

  * (1) B. Aubert et al. [BABAR Collaboration], Phys. Rev. D 80, 092003 (2009) [arXiv:0908.0806 [hep-ex]].
  * (2) B. Aubert [BABAR Collaboration], Phys. Rev. Lett. 97, 222001 (2006).
  * (3) J. Brodzicka et al. [Belle Collaboration], Phys. Rev. Lett. 100, 092001 (2008)
  * (4) E. J. Eichten, C. T. Hill and C. Quigg, Phys. Rev. Lett. 71, 4116 (1993) [arXiv:hep-ph/9308337].
  * (5) Y. B. Dai and S. L. Zhu, Phys. Rev. D 58, 074009 (1998) [arXiv:hep-ph/9802224].
  * (6) A. Hiorth and J. O. Eeg, Phys. Rev. D 66, 074001 (2002) [arXiv:hep-ph/0206158].
  * (7) S. Fajfer and J. F. Kamenik, Phys. Rev. D 74, 074023 (2006) [arXiv:hep-ph/0606278].
  * (8) S. Fajfer and J. F. Kamenik, Int. J. Mod. Phys. A 23, 3196 (2008).
  * (9) J. Kamenik and S. Fajfer, J. Phys. Conf. Ser. 110, 122015 (2008).
  * (10) W. A. Bardeen and C. T. Hill, Phys. Rev. D 49, 409 (1994) [arXiv:hep-ph/9304265].
  * (11) W. A. Bardeen, E. J. Eichten and C. T. Hill, Phys. Rev. D 68, 054024 (2003) [arXiv:hep-ph/0305049].
  * (12) P. Colangelo, F. De Fazio and R. Ferrandes, Mod. Phys. Lett. A 19, 2083 (2004) [arXiv:hep-ph/0407137].
  * (13) M. A. Nowak, M. Rho and I. Zahed, Phys. Rev. D 48, 4370 (1993) [arXiv:hep-ph/9209272].
  * (14) A. Deandrea, R. Gatto, G. Nardulli and A. D. Polosa, JHEP 9902, 021 (1999) [arXiv:hep-ph/9901266].
  * (15) A. Deandrea, N. Di Bartolomeo, R. Gatto, G. Nardulli and A. D. Polosa, Phys. Rev. D 58, 034004 (1998) [arXiv:hep-ph/9802308].
  * (16) P. Colangelo, F. De Fazio, S. Nicotri and M. Rizzi, Phys. Rev. D 77, 014012 (2008) [arXiv:0710.3068 [hep-ph]].
  * (17) F. De Fazio, arXiv:0910.0412 [hep-ph] (2009).
  * (18) Y. B. Dai, C. S. Huang, M. Q. Huang, H. Y. Jin and C. Liu, Phys. Rev. D 58, 094032 (1998) [Erratum-ibid. D 59, 059901 (1999)] [arXiv:hep-ph/9705223].
  * (19) P. Colangelo, F. De Fazio, G. Nardulli, N. Di Bartolomeo and R. Gatto, Phys. Rev. D 52, 6422 (1995) [arXiv:hep-ph/9506207].
  * (20) Z. G. Wang, Phys. Rev. D 77, 054024 (2008) [arXiv:0801.0267 [hep-ph]].
  * (21) F. E. Close, C. E. Thomas, O. Lakhina and E. S. Swanson, Phys. Lett. B 647, 159 (2007) [arXiv:hep-ph/0608139].
  * (22) F. E. Close and E. S. Swanson, Phys. Rev. D 72, 094004 (2005) [arXiv:hep-ph/0505206].
  * (23) Zhi-Feng Sun and Xiang Liu, Phys. Rev. D 80, 074037 (2009) [arXiv:0909.1658].
  * (24) B. Zhang, X. Liu, W. Z. Deng and S. L. Zhu, Eur. Phys. J. C 50, 617 (2007) [arXiv:hep-ph/0609013].
  * (25) W. Wei, X. Liu and S. L. Zhu, Phys. Rev. D 75, 014013 (2007) [arXiv:hep-ph/0612066].
  * (26) M. Di Pierro and E. Eichten, Phys. Rev. D 64, 114004 (2001).
  * (27) Z. G. Wang, Chin. Phys. C 32, 797 (2008) [arXiv:0708.0155 [hep-ph]].
  * (28) E. Van Beveren and G. Rupp, Phys. Rev. Lett. 97, 202001 (2006) [arXiv:hep-ph/0606110].
  * (29) P. Colangelo, F. De Fazio and S. Nicotri, Phys. Lett. B 642, 48 (2006) [arXiv:hep-ph/0607245].
  * (30) J. Koponen, Phys. Rev. D 78, 074509 (2008) [arXiv:0708.2807 [hep-lat]].
  * (31) X. H. Zhong and Q. Zhao, Phys. Rev. D 78, 014029 (2008) [arXiv:0803.2102 [hep-ph]].
  * (32) E. Van Beveren and G. Rupp, arXiv:0908.1142 [hep-ph].
  * (33) A. Manohar and H. Georgi, Nucl. Phys. B 234, 189 (1984).
  * (34) Z. P. Li, Phys. Rev. D 50, 5639 (1994) [arXiv:hep-ph/9404269].
  * (35) Z. P. Li, H. X. Ye and M. H. Lu, Phys. Rev. C 56, 1099 (1997).
  * (36) Q. Zhao, J. S. Al-Khalili, Z. P. Li and R. L. Workman, Phys. Rev. C 65, 065204 (2002).
  * (37) X. H. Zhong, Q. Zhao, J. He and B. Saghai, Phys. Rev. C 76, 065205 (2007).
  * (38) X. H. Zhong, Q. Zhao, Phys. Rev. C 79, 045202 (2009).
  * (39) D. O. Riska and G. E. Brown, Nucl. Phys. A 679, 577 (2001) [arXiv:nucl-th/0005049].
  * (40) Q. Zhao, Z. P. Li and C. Bennhold, Phys. Rev. C 58, 2393 (1998); Phys. Lett. B 436, 42 (1998).
  * (41) Q. Zhao, Phys. Rev. C 63, 025203 (2001).
  * (42) Q. Zhao, J. S. Al-Khalili, and C. Bennhold, Phys. Rev. C 64, 052201(R) (2001).
  * (43) Q. Zhao, B. Saghai and Z. P. Li, J. Phys. G 28, 1293 (2002).
  * (44) X. H. Zhong and Q. Zhao, Phys. Rev. D 77, 074008 (2008) [arXiv:0711.4645 [hep-ph]].
  * (45) I. Hleiqawi et al. [CLAS Collaboration], Phys. Rev. C 75, 042201 (2007) [Erratum-ibid. C 76, 039905 (2007)].
  * (46) M. Nanova et al. [CBELSA/TAPS Collaboration], Eur. Phys. J. A 35, 333 (2008) [arXiv:0803.2146 [nucl-ex]].
  * (47) C. Amsler et al. [Partical Data Group], Phys. Lett. B 667, 1 (2008).
  * (48) S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985).
  * (49) D. Ebert, R. N. Faustov and V. O. Galkin, arXiv:0910.5612 [hep-ph].
  * (50) S. Godfrey and R. Kokoski, Phys. Rev. D 43, 1679 (1991).
  * (51) E. S. Swanson, Phys. Rep. 429, 243 (2006).
  * (52) B. Chen, D. X. Wang and A. Zhang, Phys. Rev. D 80, 071502(R) (2009) [arXiv:0908.3261 ].
  * (53) D. M. Li and B. Ma, arXiv:0911.2906 [hep-ph].
  * (54) T. Matsuki, T. Morii and K. Sudoh, Eur. Phys. J.A 31,701 (2007).