# DeepAL for regression using $\epsilon$-weighted Hybrid Query strategy

Harsh Vardhan<EMAIL_ADDRESS>
Janos Sztipanovits<EMAIL_ADDRESS>
Department of Computer Science
Vanderbilt University
Nashville, TN 98195-4322, USA

###### Abstract

Designing an inexpensive approximate surrogate model that captures the salient
features of an expensive high-fidelity behavior is a prevalent approach in
design optimization. In recent times, Deep Learning (DL) models are being used
as a promising surrogate computational model for engineering problems.
However, the main challenge in creating a DL-based surrogate is to
simulate/label a large number of design points, which is time-consuming for
computationally costly and/or high dimensional engineering problems. In the
present work, we propose a novel sampling technique by combining the active
learning (AL) method with DL. We call this method $\epsilon$-weighted hybrid
query strategy ($\epsilon$-HQS), which focuses on the evaluation of the
surrogate at each learning iteration and provides an estimate of the failure
probability of the surrogate in the Design Space. By reusing already collected
training and test data, the learned failure probability guides the next
iteration’s sampling process to the region of the high probability of failure.
We empirically evaluated this method in two different engineering design
domains, finite element based static stress analysis of submarine pressure
vessel (computationally costly process) and second submarine propeller design
( high dimensional problem). Empirical evaluation shows better performance of
our technique as compared to other methods and in one case was able to achieve
same level of prediction accuracy as random sample selection process with 58%
less number of samples.
https://github.com/vardhah/epsilon_weighted_Hybrid_Query_Strategy.

Keywords: Deep Learning, Active Learning, Sampling, Propeller design , Finite
element analysis.

##  1 Problem Formulation

In the context of this paper, Design Space Exploration (DSE) is defined as a
process to characterize the response surface of a simulator’s input–output
behavior, where the input describes the relevant system parameters, and the
output describes quantities of interest. These simulators can be parameterized
and the values of a parameter can be changed in a specified bounded range.
These all possible combinations of different parameter values construct the
design space. We will first introduce our notations for formalization. Let
$DS$ be a bounded subspace of $\mathbb{R}^{d}$, i.e.,
$DS\subset\mathbb{R}^{d}$ and it is represented by an unlabeled set of $n$
samples. $d$ is the number of input parameters that describes the relevant
system properties. The training data is $D_{train}=(X,Y^{*})$, where $X$ and
$Y^{*}$ are the possible collection of all $x_{i}$ (input samples to the
simulator) and $y_{i}^{*}$ (corresponding simulation’s output). The role of a
simulator is to map the input space to the solution space. i.e.
$Sim:X\rightarrow Y^{*}$. We view these simulators as a nonlinear function
which can be modeled with a $k$–parameter configuration nonlinear regression
model i.e. $Sim=Model(\theta_{1},\theta_{2},...,\theta_{k})$. Accordingly, the
surrogate modeling problem is the conditioning of a parametric model ($h$) on
$D_{train}$ to find its parameters ($\theta$) which results into:
$h(x,\theta|D_{train})=y,\;\forall\,x\in DS:y\approx y^{*}$ where $y^{*}$
represents the ground truth. Since our evaluation is costly and/or has high
dimensional design space, we only want to take strategic samples with the goal
to obtain a certain level of prediction accuracy with the least amount of
simulation queries. In this case, our problem can be formalized as:

$\displaystyle min\;|{D_{train}}|$ $\displaystyle\textit{subject
to}\;\;h(x,\theta|D)=y,\,\forall\,x\in DS:y\approx y^{*}$

This is an optimization problem, whose direct solution is NP-hard (Cunningham
and Schrijver, 1997). However, the heuristic-based iterative approach has been
successful and this field of research is known as Active Learning (Settles,
2009). The interpolant model $h$ is an instance of a computational
architecture $A$, which defines the functional dependence of the interpolant
on its parameter $\theta$. We choose the architecture $A$ to be a Deep Neural
Network(DNN), which has a proven track of powerful learning capability, cheap
parallel prediction, and automatic feature extraction. Deployment of active
learning strategies is tightly bounded with the class of chosen architecture.
Active learning is a semi-supervised machine learning approach for sequential
sampling with the goal of identifying strategic samples. The general flow of
the AL procedures works as follows:

  1. 1.

Select some initial randomised samples via methods like Uniform random, Latin
hypercube sampling and evaluate these to create $D^{train}\leftarrow
D_{train}^{init}$.

  2. 2.

Train a machine learning model ($h$) using $D^{train}$ for constraining the
sampling procedure in next iteration.

  3. 3.

Design a query strategy ($Q$) to compute a score ($s\in\mathbb{R^{+}}$) using
the trained model ($h$) on entire design space ($DS$).

  4. 4.

Select new samples in next iteration based on score ($s$), and evaluate these
samples to create new additional training data ($D_{train}^{new}$)

  5. 5.

Update $D^{train}:D^{train}\leftarrow D_{train}^{new}$

  6. 6.

Iterate 2-5 until budget of training sample ($B=|{D_{train}}|$)

By formulating the AL framework in a pool-based sampling scheme, where
initially the design space $DS$ has an unlabeled set of $n$ samples and the
current labeled training set $D^{train}$ is the labeled set of $m$ samples,
our goal is to design a query strategy $Q$ such that $Q:DS\mapsto D^{train}$.
The combined goals in a supervised environment can be expressed as follows:

$\displaystyle\underset{D^{train}\subseteq DS,(x,y)\in
D^{train}}{\mathrm{argmin}}\mathbb{E}_{(x,y)}[l(h(x,\theta|D^{train})]$
$\displaystyle min\;|{D_{train}}|$ $\displaystyle\textit{subject
to}\;\;h(x,\theta|D)=y,\,\forall\,x\in DS:y\approx y^{*}$

where $l(h(x,\theta|D^{train})$ is the loss function and $l(.)\in R^{+}$ and
the first goal is to reduce loss on $D^{train}$ and the second goal is to
reduce $m$ ($m=|{D_{train}}|$) as small as possible while ensuring a
predetermined level of accuracy on $DS$. Therefore, the query strategy
$Q:DS\mapsto D^{train}$ is crucial. Also at each acquisition step, we want to
select a batch of data points $\\{x_{1}^{sel},x_{2}^{sel},...,x_{n}^{sel}\\}$
rather than just one point $x^{sel}$ because retraining the learning model
($h$) on every single additive data-point would be computationally infeasible
and even may not give any performance gain. So, based on our current trained
model’s state $h(w|D_{train})$, we need to design the query strategy ($Q$)
which scores a candidate batch of unlabelled data points ($L_{candidate}$)
from $DS$. An Active learning framework has various query strategies like
(uncertainty-based, variance-based, and disagreement-based) however directly
applying these strategies to regression problems out of the Bayesian framework
is less understood/studied. Accordingly, the standard AL approach is to
estimate their uncertainty directly, using Bayesian methods (Rasmussen, 2003),
or estimate uncertainty indirectly, through the variance of the outputs of an
ensemble of models. Bayesian methods are computationally expensive when the
number of samples in DS is more than $10000$ and/ or input dimensions are more
than $10$. Ensemble-like models require independent training of several
models, which is also time-consuming due to training a family of models. In
this work, we tried to solve this problem out of the Bayesian framework and
ensemble models to make it scalable and practically useful for a large class
of problems. In the next section, we will explain our approach.

##  2 Approach

To integrate the DL framework with the AL framework, we propose to use two
neural networks in a student-teacher fashion, first network is called the
student network (whose goal is to learn the input-output behavior) and the
second network is called the teacher network (whose goal is to guide next
iteration’s sampling process). For training the teacher network, we use
already evaluated historical training and test data from prior iterations and
label it with a binary indicator ($C$), where $C:X\mapsto\\{0,1\\}$, where
$C=1$, indicates student network is able to predict on the sample within an
acceptable level of accuracy otherwise $C=0$. The objective of training the
teacher network is to make an inference on the failure probability of the
student network on the entire design space. This failure probability can be
used to identify the regions/ samples in $DS$ on which student network has
poor performance/ high likelihood to fail. We continue refining the student
and teacher networks on further evaluated samples at each iteration until we
exhaust our budget for iterations. (refer figure 1)

Figure 1: Student teacher architecture: the process

Ideally, the true failure probability function $f^{*}$ is unknown but for
sampling purposes, the approximation would work. The benefit of choosing a
binary failure indicator is, that we can have a very small DNN because
creating a coarse approximation of the decision boundary is a much simpler job
than creating a nonlinear regression model. By leveraging the GPU
implementation of DNN, we can make inferences on thousands/millions of samples
at a cost negligible to running the costly simulation. With the knowledge of
the failure probability of student network, it is possible to identify at each
iteration, what are the samples $x\in DS$ on which student network has a high
likelihood to fail.

###  2.1 Batching

Using the learned probability function, the active learning framework ideally
would select one sample per iteration which has maximum probability to fail.
However, this one additive sample per iteration does not go well with the DL
framework because of retraining time. So, at each iteration, a batch of data
points ($L_{candidate}$) should be selected. For a given batch size ($B$) ,
the goal is to design an query strategy $Q$ such that
$Q:L_{candidate}\rightarrow DS$. The most obvious approach for batching is to
select top $B$ samples sorted in descending order by failure probability, we
call this approach as top-B, where $B$ is the batch size at each iteration.
However, this method is highly likely to select a set of information-rich but
similar samples because the query strategy considers each sample independently
and ignores the correlation between samples. This may result in a wastage of
the evaluation budget. Our goal is not only to select information-rich samples
but also to maintain diversity in the selected samples. One way to solve it is
by combining probability information with unsupervised clustering. Zhdanov
(2019) proposed the batched active learning using the clustering approach and
it utilizes both the amount of information in the probability function as well
as the diversity of information. We call this approach Diverse Batched Active
Learning (DBAL).

###  2.2 Debiasing scheme

Empirical evaluation of both DBAL and top-B sampling approaches performed
poorer than randomized sampling. The underlying reason is, in the initial
phase of training, the teacher network which learns the probability function
($f:x\rightarrow[0,1]\;\forall\;x\in DS$) has only limited information about
the design space $DS$ (a maximum up to training and testing data).
Accordingly, the “trained model does not know what it does not know” and it
creates a biased estimate of the probability function. If $f^{*}$ is the true
failure probability function then $f$ can approximate it only up to $\beta-
approximation$. The value of $\beta$ is an open research question, but we can
increase the robustness of the teacher network and can reduce the sensitivity
of mismatch between the distribution of $f$ and $f^{*}$. For this purpose, we
introduce a hyper-parameter $\epsilon$ ($0<\epsilon\leq 1$) and adopt
filtering of evaluated samples using rejection threshold and then running
$\epsilon$-weighted Monte Carlo sampling on filtered samples. The $\epsilon$
also controls the exploration similar to used in Reinforcement Learning(Sutton
and Barto, 2018)(refer algorithm 1). The value of $\epsilon$ reflects our
belief in the accuracy of the failure predictor. In the next section, we will
empirically evaluate our approach with all above mentioned acquisition
functions: top-B, DBAL-10 , DBAL-50, $\epsilon$-weighted Hybrid Query
Strategy, Uniformly Random Monte Carlo.

Algorithm 1 $\epsilon$-weighted Sampling and de-biasing scheme

1:$Input:{trained\;teacher\;model(f),Design\;space(DS),Budget(B),\&\;\epsilon}$

2:$Output:{Candidate\;Samples(L_{candidate})}$

3:$Initialization:S\leftarrow\\{\\}$

4:for $size(S)\leq\epsilon*B$ do

5: $Sample\;proposal\;instance/full\;DS,X_{t}$

6: $Calculate\;f(X_{t})$ $\triangleright$ Estimate failure probability on
samples in DS

7: $X_{f}\subseteq X_{t}:X_{f}=X_{t}\;\forall\;X_{t},\;where\;f(X_{t})\geq
0.5$ $\triangleright$ filtered high failure probability samples $X_{f}$

8: $X_{sel}=Choose\\_Uniform\\_rand(X_{f},\epsilon*B)$ $\triangleright$ Choose
$\epsilon*B$ samples randomly from $X_{f}$

9: $S\leftarrow S+X_{sel}$

10:$X_{rest}=Choose\\_Uniform\\_rand(X_{t}-X_{sel},(1-\epsilon)*B)$

11:$S\leftarrow S+X_{rest}$

##  3 Related Work

There are many ongoing and historical efforts for creating learning-based
surrogate models in both propeller design and the FEA domain.(Vesting and
Bensow, 2014)(Vardhan et al., 2021)(Liang et al., 2018). However, without any
strategic sampling, all work is done for limited design space. Active
learning-based sampling by utilizing information-based objective function has
also a long history(Fedorov, 2013)(MacKay, 1992)(Settles, 2009) and has been
extensively used for regression problems under the Bayesian framework
especially by using Kriging and Gaussian Process-based learning
models(Forrester et al., 2008)(Rasmussen, 2003). With the advent of deep
learning(LeCun et al., 2015), it is found advantageous to integrate DL
models(which is a powerful learning model) with AL to get the best of both
worlds however the traditional AL literature has little to offer when the
parameters of learning models are large. There is a plethora of ongoing
research work in integrating DL with AL and it is called DeepAL but most work
is done in the context of classification problems where direct probabilistic
output makes creating information-based objective function much easier(Ren et
al., 2021). Our work is more inspired by (Uesato et al., 2018)(Vardhan and
Sztipanovits, 2021), in which an attempt made to learn the failure probability
in the context of rare event failures in RL agents.

##  4 Experiments and Evaluation

For empirical evaluation, we select two engineering domains–propeller design
domain and the Finite Element Analysis (FEA) domain. In both cases, we want to
replace a physics-based simulation process with a computationally cheap DNN-
based surrogate model. In the propeller domain, the goal is to create a
surrogate model to predict the efficiency of the propeller on the given
requirement and geometric design. In the case of the FEA domain, the goal is
to create a surrogate model of the FEA process to predict the maximum Von
mises stress (which determines the integrity of the vessel) in the capsule-
shaped pressure vessel used in an Underwater vehicle (refer figures 3 and 3).

Figure 2: Propeller design surrogate

Figure 3: Pressure vessel Surrogate

For data generation and evaluation, we used physics-based simulators -openprop
(Epps, Chalfant, Kimball, Techet, Flood, and Chryssostomidis, 2009) and
freeCAD (Riegel, Mayer, and van Havre, 2016) for the propeller and FEA domain
respectively. In propeller cases, we generated 275000 labeled data in $14$
dimensional input parameters space because the average evaluation time per
sample was small ($\approx\;0.5\;sec$). However, in the case of pressure
vessels, the evaluation time is high ($\approx\;202\;sec$) and we generated
11131 labeled samples in the $4$ dimensional input parameter space. For
experiment I in both domains, the total budget of training data is fixed to
2500 (50 iterations of the AL process, with a budget of 50 per iteration). The
leftover labeled data is used as the test data (to measure the accuracy of the
surrogate). For evaluating the consistency of results, we ran multiple
experiments and averaged the prediction accuracy.

We observed that in both domains, DBAL and top-b strategies perform poorer
than the random sampling and the underlying reason is the complete reliance on
a biased estimation of failure probability.

Figure 4: Propeller Domain : Experiment I

Figure 5: FEA domain :Experiment II

Figure 7 shows the accuracy of trained DL-based surrogates using different
sampling strategies (DL models and other learning parameters are kept the
same). Our accuracy metrics is fraction of test predictions which are out of
acceptable accuracy(+/- 5% error) i.e.,
$\frac{\sum_{i=1}^{T}\mathbb{I}(|y_{i}-y_{i}^{*}|\leq 0.05*|y_{i}^{*}|)}{T}$.
Here $y_{i}$ is a prediction on $i^{th}$ test sample having ground truth
$y_{i}^{*}$ and $T$ is the total number of samples in the test data. The
DBAL50 and DBAL10 are two different DBAL (Diverse Batched Active learning)
strategies with hyper-parameter $\beta=50$ and 10 respectively(refer Zhdanov
(2019)). In random strategy, at each iteration samples are selected uniformly
randomly for unlabeled pool data($DS$). Five $\epsilon-HQS$ approaches are
tested labeled as $ep\\_g,ep\\_025,ep\\_05,ep\\_075$ and $ep\\_1$, called as
logarithmic increasing $\epsilon$, $\epsilon$ with weights 0.25,0.50, 0.75 and
1.0 respectively. The epsilon-greedy(Sutton and Barto, 2018) approach
increases belief in the prediction of teacher network logarithmic, while other
fixed $\epsilon$ strategies assign linear belief to the teacher network. For
experiment II, we used only FEA domain and reduced the batch size per
iteration to 20 and the increased the number of iterations to 400 (i.e., a
total budget of training is increased from 2500 to 8000), and the rest of the
data is used for testing. We observed the improvement in accuracy of surrogate
using $epsilon$-greedy HQS at a higher rate than randomized sampling as the
$epsilon$-greedy HQS can find strategic samples, which random sampling cannot.
The accuracy increment was 4.14 %(86.47% to 90.61%) for $epsilon$-greedy HQS
in comparison to 1.99%(83.76% to 85.75%) in the case of random sampling when
training data increased from 2500 samples to 8000 samples. Another result that
we want to share is for reaching the same level of accuracy as the randomized
sampling strategy, the $epsilon$-greedy HQS saves $57.5\%$ of samples and
consequently $250$ hours of data generation time (includes the samples
simulation time minus training and prediction time from teacher network).

Figure 6: Prediction accuracies

Figure 7: FEA domain: Experiment II

##  5 Conclusion and Future Work

In this work, we developed a new method of data efficient surrogate modeling
by combining Deep Learning and Active Learning for regression problem. When
compared this methods like DBAL, VMC and top-B strategy, it outperformed all
of them.

##  6 Acknowledgments

This work is supported by DARPA’s Symbiotic Design for CPS project
(FA8750-20-C-0537).

## References

  * Cunningham and Schrijver (1997) William J Cook William H Cunningham and William R Pulleyblank Alexander Schrijver.  Combinatorial optimization, 1997.
  * Epps et al. (2009) Brenden Epps, Julie Chalfant, Richard Kimball, Alexandra Techet, Kevin Flood, and Chrysssostomos Chryssostomidis.  Openprop: An open-source parametric design and analysis tool for propellers.  In _Proceedings of the 2009 grand challenges in modeling & simulation conference_, pages 104–111, 2009.
  * Fedorov (2013) Valerii Vadimovich Fedorov.  _Theory of optimal experiments_.  Elsevier, 2013.
  * Forrester et al. (2008) Alexander Forrester, Andras Sobester, and Andy Keane.  _Engineering design via surrogate modelling: a practical guide_.  John Wiley & Sons, 2008.
  * LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Liang et al. (2018) Liang Liang, Minliang Liu, Caitlin Martin, and Wei Sun.  A deep learning approach to estimate stress distribution: a fast and accurate surrogate of finite-element analysis.  _Journal of The Royal Society Interface_ , 15(138):20170844, 2018.
  * MacKay (1992) David JC MacKay.  Information-based objective functions for active data selection.  _Neural computation_ , 4(4):590–604, 1992.
  * Rasmussen (2003) Carl Edward Rasmussen.  Gaussian processes in machine learning.  In _Summer school on machine learning_ , pages 63–71. Springer, 2003\.
  * Ren et al. (2021) Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Brij B Gupta, Xiaojiang Chen, and Xin Wang.  A survey of deep active learning.  _ACM Computing Surveys (CSUR)_ , 54(9):1–40, 2021\.
  * Riegel et al. (2016) Juergen Riegel, Werner Mayer, and Yorik van Havre.  Freecad, 2016.
  * Settles (2009) Burr Settles.  Active learning literature survey.  2009\.
  * Sutton and Barto (2018) Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Uesato et al. (2018) Jonathan Uesato, Ananya Kumar, Csaba Szepesvari, Tom Erez, Avraham Ruderman, Keith Anderson, Nicolas Heess, Pushmeet Kohli, et al.  Rigorous agent evaluation: An adversarial approach to uncover catastrophic failures.  _arXiv preprint arXiv:1812.01647_ , 2018.
  * Vardhan and Sztipanovits (2021) Harsh Vardhan and Janos Sztipanovits.  Rare event failure test case generation in learning-enabled-controllers.  In _2021 6th International Conference on Machine Learning Technologies_ , pages 34–40, 2021.
  * Vardhan et al. (2021) Harsh Vardhan, Peter Volgyesi, and Janos Sztipanovits.  Machine learning assisted propeller design.  In _Proceedings of the ACM/IEEE 12th International Conference on Cyber-Physical Systems_ , pages 227–228, 2021.
  * Vesting and Bensow (2014) Florian Vesting and Rickard E Bensow.  On surrogate methods in propeller optimisation.  _Ocean engineering_ , 88:214–227, 2014.
  * Zhdanov (2019) Fedor Zhdanov.  Diverse mini-batch active learning.  _arXiv preprint arXiv:1901.05954_ , 2019.


#  $D$ meson mixing as an inverse problem

Hsiang-nan Li  Institute of Physics, Academia Sinica, Taipei, Taiwan 115,
Republic of China  Hiroyuki Umeeda  Institute of Physics, Academia Sinica,
Taipei, Taiwan 115, Republic of China  Fanrong Xu  Department of Physics,
Jinan University, Guangzhou 510632, People’s Republic of China  Fu-Sheng Yu
School of Nuclear Science and Technology, Lanzhou University, Lanzhou 730000,
People’s Republic of China

###### Abstract

We calculate the parameters $x$ and $y$ for the $D$ meson mixing in the
Standard Model by considering a dispersion relation between them. The
dispersion relation for a fictitious charm quark of arbitrary mass squared $s$
is turned into an inverse problem, via which the mixing parameters at low $s$
are solved with the perturbative inputs $x(s)$ and $y(s)$ from large $s$. It
is shown that nontrivial solutions for $x$ and $y$ exist, whose values around
the physical charm scale agree with the data in both CP-conserving and CP-
violating cases. We then predict the observables $|q/p|-1\approx 2\times
10^{-4}$ and $Arg(q/p)\approx 6\times 10^{-3}$ degrees associated with the
coefficient ratio for the $D$ meson mixing, which can be confronted with more
precise future measurements. Our work represents the first successful
quantitative attempt to explain the $D$ meson mixing parameters in the
Standard Model.

How to understand the large $D$ meson mixing in the Standard Model has been a
long-standing challenge. Previous evaluations based on box diagrams Cheng ;
BSS ; Datta:1984jx , and on heavy quark effective field theory Georgi:1992as ;
Ohl:1992sr  led to the mixing parameters $x$ and/or $y$ far below the current
data. The updated inclusive analysis Golowich:2005pt , including next-to-
leading-order QCD corrections, still gave small $x\sim y\simeq 6\times
10^{-7}$. Some authors Bigi:2000wn ; Falk:2001hx ; Bobrowski:2010xg  claimed
that higher dimensional operators, for which the strong Glashow-Iliopoulos-
Maiani (GIM) suppression Glashow:1970gm  might be circumvented, yielded
dominant contributions. This claim has not been verified quantitatively, which
requires information on a large number of nonperturbative matrix elements.
Another uncertainty in the heavy quark expansion originates from violation of
the quark-hadron duality, which represents an error in the analytic
continuation from deep Euclidean to Minkowskian domains. A simple
phenomenological argument Jubb:2016mvq  indicated that 20% duality violation
could explain the width difference in the presence of the GIM cancellation.

On the other hand, the exclusive approach, where the $D$ meson mixing is
extracted from hadronic processes, led to an enhancement by relevant long-
distance effects Wolfenstein:1985ft ; Donoghue:1985hh ; Colangelo:1990hj ;
Buccella:1994nf ; Kaeding:1995zx ; Falk:2001hx ; Falk:2004wg ; Cheng:2010rv ;
Gronau:2012kq ; Jiang:2017zwr . Modern works along this direction, e.g.,
Cheng:2010rv ; Jiang:2017zwr  showed that a half value of $y$ was accounted
for roughly with contributions from two-body decays, albeit the difficulty in
taking account of other multi-body channels. Thus, the quantitative
understanding is still not attained in this data-driven approach, while the
order of magnitude of the mixing parameters was properly described.

The complexities are attributed to the notorious difficulty of charm physics:
the charm scale is too heavy to apply the chiral perturbation theory and
possibly too light to apply the heavy quark expansion. Moreover, the $D$ meson
mixing, strongly suppressed by the GIM mechanism, is sensitive to
nonperturbative SU(3) breaking effects Kingsley:1975fe  characterized by the
strange and down quark mass difference, and to CKM-suppressed diagrams with
bottom quarks in the loop. On the contrary, the heavy quark expansion
accommodates the data for the $B_{d,s}$ meson mixings satisfactorily
Artuso:2015swg ; Jubb:2016mvq .

In this letter we will analyze the $D$ meson mixing in a novel approach based
on a dispersion relation, which relates $x$ and $y$ for a fictitious $D$ meson
of an arbitrary mass. The dispersion relation is separated into a low mass
piece and a high mass piece, with the former being treated as an unknown, and
the latter being input from reliable perturbative results. We then turn the
study of the $D$ meson mixing into an inverse problem: the mixing parameters
at low mass are solved as ”source distributions”, which produce the
”potential” observed at high mass. It will be demonstrated that nontrivial
correlated solutions for $x$ and $y$ exist, whose values around the physical
charm quark mass $m_{c}\approx 1.3$ GeV match the data in both cases with and
without CP violation. Our observation implies that resonance properties can be
extracted from asymptotic QCD by solving an inverse problem.

Consider the analytical transition matrix element for a $D$ meson formed by a
fictitious charm quark of invariant mass squared $s$,

$\displaystyle M_{12}(s)-\frac{i}{2}\Gamma_{12}(s)=\langle D^{0}(s)|{\cal
H}_{w}^{\Delta C=2}|\bar{D}^{0}(s)\rangle,$ (1)

whose branch cut runs from the threshold $s=4m_{\pi}^{2}$ to infinity with the
pion mass $m_{\pi}$. The effective weak Hamiltonian ${\cal H}_{w}^{\Delta
C=2}$ contains two four-fermion operators $(V-A)(V-A)$ and $(S-P)(S-P)$, which
will be abbreviated to $V-A$ and $S-P$ below, respectively. The right hand
side of Eq. (1) starts with the evaluation of box diagrams, whose dispersive
and absorptive contributions give rise to $M_{12}$ and $\Gamma_{12}$,
respectively. The dispersive part $M_{12}$ and the absorptive part
$\Gamma_{12}$ then obey the dispersion relation Falk:2004wg

$\displaystyle
M_{12}(s)=\frac{P}{2\pi}\int_{0}^{\infty}ds^{\prime}\frac{\Gamma_{12}(s^{\prime})}{s-s^{\prime}},$
(2)

where $P$ denotes the principal value prescription, and the lower bound of the
integration variable $s^{\prime}$, being of $O(m_{\pi}^{2})$, has been
approximated by zero. Equation (1) governs the time evolution of the $D^{0}$
and $\bar{D}^{0}$ mesons, whsoe diagonalization yields the mass eigenstates
$D_{1,2}=pD^{0}\pm q\bar{D}^{0}$ as linear combinations of the weak
eigenstates $D^{0}$ and $\bar{D}^{0}$. The mass and width differences of
$D_{1,2}$ define the mixing parameters

$\displaystyle
x\equiv\frac{m_{1}-m_{2}}{\Gamma}=\frac{2M_{12}}{\Gamma},\;\;\;\;y\equiv\frac{\Gamma_{1}-\Gamma_{2}}{2\Gamma}=\frac{\Gamma_{12}}{\Gamma},$
(3)

in the CP-conserving case with the total decay width $\Gamma$.

Figure 1:  $\lambda_{s}^{2}$, $\lambda_{s}\lambda_{b}$ and $\lambda_{b}^{2}$
contributions to the real parts of $M_{12}$ and $\Gamma_{12}$ from the $S-P$
and $V-A$ operators.

The elements $M_{12}$ and $\Gamma_{12}$, extracted from the evaluation of box
diagrams Cheng ; BSS , can be applied to the mixing of a heavy meson with
arbitrary mass, and will be adopted directly below. The $b$ quark mass $m_{b}$
should remain constant in the evaluation of $\Gamma_{12}$, so that the
fictitious $D$ meson can decay into a $b$ quark, as its mass crosses the $b$
quark threshold. The right hand side of Eq. (2) then contains heavy quark
contributions to be consistent with the heavy quark dynamics involved in
$M_{12}$. The $V-A$ contribution $\Gamma_{12}^{V-A}$ is given, for
$s>4m_{b}^{2}$, by

$\displaystyle\Gamma_{12}^{V-A}\propto\lambda_{s}^{2}(B_{dd}^{(a)}-2B_{ds}^{(a)}+B_{ss}^{(a)})+2\lambda_{s}\lambda_{b}(B_{dd}^{(a)}-B_{ds}^{(a)}-B_{db}^{(a)}+B_{sb}^{(a)})+\lambda_{b}^{2}(B_{dd}^{(a)}-2B_{db}^{(a)}+B_{bb}^{(a)}),$
(4)

where $\lambda_{k}\equiv V_{ck}V^{*}_{uk}$, $k=s,b$, are the products of the
Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, and the functions
$B_{ij}^{(a)}$ BSS  with the internal quarks $i,j=d,s,b$ arise from the
absorptive contributions of the box diagrams for Eq. (1). The terms up to
$B_{ss}^{(a)}$ ($B_{db}^{(a)}$, $B_{sb}^{(a)}$) are kept in the range
$s<(m_{b}+m_{d})^{2}$ ($(m_{b}+m_{d})^{2}\leq s<(m_{b}+m_{s})^{2}$,
$(m_{b}+m_{s})^{2}\leq s<4m_{b}^{2}$). The expression of the $S-P$
contribution $\Gamma_{12}^{S-P}$ is similar but with $B_{ij}^{(a)}$ in Eq. (4)
being replaced by $C_{ij}^{(a)}$ BSS . Equation (4) shows clearly that the $D$
meson mixing results from the flavor symmetry breaking. We have confirmed that
$\Gamma_{12}$ deceases like $1/s^{2}$ at large $s$, so the integral on the
right hand side of Eq. (2) converges.

Figure 2:  Dispersive and absorptive contributions to the real part of
$\omega(s)$ from the $S-P$ and $V-A$ operators.

We rewrite the dispersion relation as

$\displaystyle\int_{0}^{\Lambda}ds^{\prime}\frac{y(s^{\prime})}{s-s^{\prime}}=\pi
x(s)-\int_{\Lambda}^{\infty}ds^{\prime}\frac{y(s^{\prime})}{s-s^{\prime}}\equiv\omega(s),$
(5)

where both sides have been divided by the measured total width $\Gamma_{\rm
exp}=1.61\times 10^{-12}$ GeV PDG  to get the variables $x$ and $y$. The
separation scale $\Lambda$ is arbitrary, but should be large enough to justify
the perturbative calculation of $y$ on the right hand side, and below the $b$
quark threshold to avoid the $b$ quark contribution to the left hand side. The
product $f_{D}^{2}m_{D}$ appearing in the expressions of $M_{12}$ and
$\Gamma_{12}$ BSS  on the right hand side of Eq. (5), with the $D$ meson decay
constant $f_{D}$ and its mass $m_{D}$, scales like a constant in the heavy
quark limit. Here we adopt the value for a $B_{s}$ meson PDG , ie.,
$f_{D}^{2}m_{D}\sim 0.3$ GeV3. The behaviors of $M_{12}(s)$ and
$\Gamma_{12}(s)$ from the $S-P$ and $V-A$ operators with the masses $m_{d}=5$
MeV, $m_{s}=109.9$ MeV, $m_{b}=4.8$ GeV and $m_{W}=80.379$ GeV, the separation
scale $\Lambda=m_{b}^{2}/2\approx 12$ GeV2, and the bag parameters equal to
unity are displayed in Fig. 1, which have been decomposed into three pieces
proportional to the real parts of $\lambda_{s}^{2}$, $\lambda_{s}\lambda_{b}$
and $\lambda_{b}^{2}$. The above choice of $\Lambda$ can be regarded as being
of $O(m_{b}^{2})$, so the perturbation theory is applicable to the mixing of
the fictitious $D$ meson for $s>\Lambda$. The choice of $m_{s}=109.9$ MeV is
within the range of the strange quark mass $m_{s}=108^{+13}_{-6}$ MeV given
for the renormalization scale $\mu=m_{c}$ in PDG . It is seen in Fig. 2 that
both terms on the right-hand side of Eq. (5) exhibit cusps as $s$ crosses the
$b$ quark and $b$ quark pair thresholds. Their sum $\omega(s)$ behaves
smoothly and, furthermore, turns out to be independent of $m_{b}$. This
feature, existent for the two four-fermion operators, indicates that $y$ in
the low mass region $s<\Lambda$ decouples from the $b$ quark dynamics as
expected.

In principle, we can have separate dispersion relations associated with the
three CKM products. However, it is reasonable to combine all the terms in Eq.
(4) into a single dispersion relation due to the dominance of the
$\lambda_{s}^{2}$ contribution to the real part of $\omega(s)$. The $S-P$ and
$V-A$ contributions are opposite in sign, and the corresponding bag parameters
are roughly equal. The significant cancellation between these two pieces
causes sensitivity to the bag parameters, which have not yet been computed
precisely enough in lattice QCD. To reduce the sensitivity to this potential
cancellation, we consider separate dispersion relations for these two
operators. Equation (5) will be treated as an inverse problem, in which
$\omega(s)$ for $s>\Lambda$ from Fig. 2 is an input, and $y(s)$ in the range
$s<\Lambda$ is solved with the boundary condition $y(0)=0$ and the continuity
of $y$ at $s=\Lambda$. That is, the ”source distribution” $y(s)$ will be
inferred from the ”potential” $\omega(s)$ observed outside the distribution.

For such an ill-posed inverse problem, the ordinary discretization method to
solve an integral (Fredholm) equation does not work. The discretized version
of Eq. (5) is in the form $\sum_{i}A_{ij}y_{j}=\omega_{i}$ with $A_{ij}\propto
1/(i-j)$. It is easy to find that any two adjacent rows of the matrix $A$
approach to each other as the grid becomes infinitely fine. Namely, $A$ tends
to be singular, and has no inverse. We stress that this singularity, implying
no unique solution, should be appreciated actually. If $A$ is not singular,
the solution to Eq. (5) will be unique, which must be the tiny perturbative
result obtained in the literature. It is the existence of multiple solutions
that allows possibility to explain the observed large $D$ meson mixing.

Figure 3:  $s$ dependence of $1/\omega$.

We notice that the smooth curves of $\omega(s)$ can be well described by
simple functions proportional to $1/(s-m^{2})$, as indicated by the almost
straight lines for $1/\omega(s)$ down to $s=15$ GeV2 in Fig. 3. These straight
lines, as extrapolating to the low $s$ region, cross the horizontal axis at
some small scale $s=m^{2}$. The power-law behavior is understandable, since
only the effect from the monopole component of the distribution dominates at
large $s$, which decreases like $1/s$. The meaning of the scale $m^{2}$ will
become clear later. If $\omega(s)$ followed the power law exactly, the
solution to Eq. (5) would be a $\delta$-function,
$y(s)\propto\delta(s-m^{2})$. The slight deviation from the power-law behavior
suggests mild broadening of $y(s)$ into a resonance-like distribution located
at $s\approx m^{2}$, if $m^{2}>0$.

Viewing the difficulty to solve an inverse problem with multiple solutions and
the qualitative resonance-like behavior of a solution, we propose the
parametrization

$\displaystyle
y(s)=\frac{Ns[b_{0}+b_{1}(s-m^{2})+b_{2}(s-m^{2})^{2}]}{[(s-m^{2})^{2}+d^{2}]^{2}},$
(6)

for $0\leq s\leq\Lambda$, and determine the free parameters $b_{0}$, $b_{1}$,
$b_{2}$, $m^{2}$ and $d$ from the best fit to the input $\omega(s)$. The
normalization $N$ respects $Ns/[(s-m^{2})^{2}+d^{2}]^{2}\to\delta(s-m^{2})$ in
the vanishing width limit $d\to 0$. Equation (6) with the completely free
parameters is general enough, which can also describe a nonresonant behavior
with $m^{2}<0$ and a flat behavior with large $d$. The convergence of the
expansion in the numerator will be verified, so keeping terms up to
$(s-m^{2})^{2}$ is sufficient. Equation (6) obeys the boundary condition
$y(0)=0$. The continuity of $y(s)$ at $s=\Lambda$, ie., the equality of
$y(\Lambda)$ to the perturbative input imposes a constraint among the five
parameters. We emphasize that a systematic expansion of $y(s)$ in terms of a
complete basis of orthogonal functions also works, but the numerical analysis
is more tedious, and will be performed in a forthcoming paper.

The separation scale $\Lambda$ introduces an end-point singularity to the
integral on the right hand side of Eq. (5), as $s\to\Lambda$. To reduce the
effect caused by this artificial singularity, we consider $\omega(s)$ from the
range 30 GeV${}^{2}<s<$ 250 GeV2, in which 200 points $s_{i}$ are selected. We
have checked the cases with 100, 200 and 300 points, and confirmed that the
results have little dependence on these numbers. For each point $(m^{2},d)$ on
the $m^{2}$-$d$ plane, we search for $b_{0}$ and $b_{1}$, that minimize the
deviation

$\displaystyle\sum_{i=1}^{200}\left|\int_{0}^{\Lambda}ds^{\prime}\frac{y(s^{\prime})}{s_{i}-s^{\prime}}-\omega(s_{i})\right|^{2}.$
(7)

The above definition, characterizing the relative quality of solutions, is
referred to as the goodness-of-fit (GOF) hereafter. The value of $b_{2}$ is
fixed by the continuity constraint at $s=\Lambda$. The scanning on the
$m^{2}$-$d$ plane generates the arc-shaped distribution of the GOF minima
associated with the $S-P$ operator in Fig. 4, which ranges roughly in -0.2
GeV${}^{2}<m^{2}<1.8$ GeV2. The minima along the arc, having similar GOF about
$10^{-21}$-$10^{-22}$ relative to $10^{-17}$ from outside the arc, hint the
existence of multiple solutions. If a resonance-like solution with $m^{2}\sim
m_{c}^{2}$ and small $d$ exists, ie., obeys the dispersion relation, it will
be revealed by the scanning, and indeed it is as shown in Fig. 4.

Figure 4:  Distributions of GOF minima in the common logarithmic scale on the
$m^{2}$-$d$ plane.

Evaluating $y(s)$ at low $s$ in perturbation theory with a finite running
coupling constant $\alpha_{s}$, we get different results at various orders.
These different results lead to almost identical $\omega(s)$ in the large $s$
limit, where $\alpha_{s}$ diminishes. The solutions from $m^{2}$ away from
$m_{c}^{2}$ might correspond to fixed-order results, since they generate tiny
$y$ at the physical scale $m_{c}^{2}$, while those near $m^{2}\approx
m_{c}^{2}$ correspond to nonperturbative results. We select a typical
nonresonant solution for $y(s)$ with $m^{2}=0$ and $d=0.38$ GeV2 from the arc
associated with the $S-P$ operator, and compare it to a resonance-like
solution with $m^{2}=1.713$ GeV2 and $d=3.876\times 10^{-2}$ GeV2 from the
same arc in Fig. 5. The dramatic distinction in the shape and in the order of
magnitude between these two cases supports that Eq. (6) is general enough to
exhibit very different behaviors. The observation that the above perturbative
and nonperturbative solutions give the same $\omega(s)$ at large $s$ realizes
the concept of the global quark-hadron duality postulated in QCD sum rules
Shifman:2000jv . The arc-shaped distribution from the $V-A$ operator is also
displayed in Fig. 4, where a solution with $m^{2}\approx m_{c}^{2}$ has a
large $d$, so its contribution to $y$ is negligible.

Figure 5:  Comparison of a typical nonresonant solution and a resonance-like
solution for $y$. The horizontal band represents the data with $1\sigma$
errors, and the vertical line corresponds to $s=m_{c}^{2}$.

Selecting a point $(m^{2},d)$ on the arc, we get a solution of $y(s)$.
Substituting the obtained $y(s)$, ie., $\Gamma_{12}(s)$ in the whole range of
$s$ into Eq. (2), we calculate the corresponding $x(s)$. The values
$x(m_{c}^{2})$ and $y(m_{c}^{2})$ are then compared with the data. It is seen
in the left plot of Fig. 6 that the data $x=(0.50^{+0.13}_{-0.14})\%$ and
$y=(0.62\pm 0.07)\%$ in the CP-conserving case Amhis:2019xyh  can be
accommodated simultaneously by the $S-P$ contribution with the parameters

$\displaystyle m^{2}=1.713\;{\rm GeV}^{2},\;\;d=3.876\times 10^{-2}\;{\rm
GeV}^{2},$ $\displaystyle b_{0}=-3.296\times 10^{-5}\;{\rm
GeV}^{2},\;\;b_{1}=-3.234\times 10^{-2},\;\;b_{2}=5.617\times 10^{-2}\;{\rm
GeV}^{-2}.$ (8)

Equation (8) justifies that the lower bound of the integral in Eq. (2), being
of $O(m_{\pi}^{2})$, can be set to zero safely, because $y(s)$ takes
substantial values only around $s\sim O(m_{c}^{2})$. We remind that the values
in Eq. (8) are representative, and their slight variations are allowed for
explaining the data of $x$ and $y$ within 1$\sigma$. For instance, the width
$d$ is allowed to vary by 20%. The uncertainties from the fitting procedure
and from the parametrization for $y(s)$ will be investigated rigorously in a
subsequent publication.

It has been concluded Jiang:2017zwr  that two-body modes in $D$ meson decays
are insufficient for understanding $y$, and multi-particle modes play a
crucial role for this purpose. When $s$ increases, single strange quark
channels with destructive contributions, like $KKK\pi$, are enhanced by phase
space, and double strange quark channels with constructive contributions, like
$KKKK$, are opened. This tendency fits the behavior of $y(s)$ in Fig. 6, which
first decreases from a positive value expected in the two-body analysis
Jiang:2017zwr  to a negative value, and then increases with $s$. It also
explains why the width $d$, within which the above oscillation occurs, is of
$O(m_{s}^{2})$. As a single resonance around $s\approx m_{c}^{2}$ accommodates
the data, it would hint that the multi-particle channels with the total rest
mass around $m_{c}$ give dominant contributions to $y$. Certainly, our
observation does not exclude other shorter peaks at lower $s$ but above the
threshold for two-body channels. That is, the curve in Fig. 6 has caught the
major features of $y(s)$, though its true behavior might be more complicated.
We have also examined that the $b_{1}$ term dominates, and the $b_{2}$ term
contributes only about 10% of $x(m_{c}^{2})$ and $y(m_{c}^{2})$. The
convergence of the parametrization in Eq. (6) is verified. To test whether
$x(m_{c}^{2})$ and $y(m_{c}^{2})$ exhibit a quadratic rise with $m_{s}$, as
expected from the SU(3) symmetry breaking Falk:2001hx , we fix
$\Lambda=12~{}\mathrm{GeV}^{2}$ and $m^{2}=1.713\;{\rm GeV}^{2}$ in Eq. (8),
and then derive $x(m_{c}^{2})$ and $y(m_{c}^{2})$ from the dispersion relation
for various $m_{s}$. The quadratic increase with a vanishing slope at small
$m_{s}$ is indeed observed.

Figure 6:  $x(s)$ and $y(s)$ in the cases without and with CP violation.

As CP violation is allowed, both $M_{12}$ and $\Gamma_{12}$ become complex due
to the weak phase in the CKM matrix elements, but Eq. (2) still holds. The
expressions of $x$ and $y$ in terms of the complex $M_{12}$ and $\Gamma_{12}$
are referred to PDG . In this case the same parameters $\Lambda=12$ GeV2 and
$m_{s}=109.9$ MeV are chosen, and the $\lambda_{b}^{2}$ contribution is found
to dominate the imaginary part of $\omega(s)$. An additional parametrization
similar to Eq. (6) but with primed parameters is proposed. The imaginary part
of $\omega(s)$ is fitted by the primed parametrization independently of the
fitting to its real part. The scanning on the $m^{2}$-$d$ planes yield the
arc-shaped distributions of the GOF minima similar to Fig. 4. Taking a common
value for $m^{2}$ and $m^{\prime 2}$, one finds that the real part of the
$S-P$ contribution still dominates $x$ and $y$. The parameters in Eq. (8) and

$\displaystyle m^{\prime 2}=m^{2}=1.713\;{\rm
GeV}^{2},\;\;d^{\prime}=4.970\;{\rm GeV}^{2},$ $\displaystyle
b^{\prime}_{0}=-8.238\times 10^{-7}\;{\rm
GeV}^{2},\;\;b^{\prime}_{1}=4.355\times
10^{-7},\;\;b^{\prime}_{2}=-7.192\times 10^{-8}\;{\rm GeV}^{-2},$ (9)

for the $S-P$ imaginary contribution, and those for the $V-A$ contribution,
which are not presented for simplicity, accommodate the data
$x=(0.39^{+0.11}_{-0.12})\%$ and $y=(0.651^{+0.063}_{-0.069})\%$ Amhis:2019xyh
simultaneously, as illustrated in the right plot of Fig. 6. Given the
parameters in Eqs. (8) and (9) and those of the $V-A$ contribution, we then
derive $|q/p|-1\approx 2.2\times 10^{-4}$ and $Arg(q/p)\approx(6.2\times
10^{-3})^{\circ}$ associated with the coefficient ratio as predictions, which
are comparable to the data
$q/p=(0.969^{+0.050}_{-0.045})e^{i(-3.9^{+4.5}_{-4.6})^{\circ}}$ Amhis:2019xyh
, and can be confronted with more precise measurements in future.

To examine the uncertainty from the theoretical input, we increase $m_{s}$ to,
say, 130.6 MeV, for which $\Lambda$ needs to increase to 14 GeV2 accordingly
to accommodate the observed $x$ and $y$. That is, a positive correlation
between $m_{s}$ and $\Lambda$ is observed. In this case, the representative
parameters for the real and imaginary parts of the $S-P$ contribution from the
fit are

$\displaystyle m^{2}=m^{\prime 2}=1.720~{}\mathrm{GeV},\;\;d=5.467\times
10^{-2}~{}\mathrm{GeV}^{2},\;\;d^{\prime}=5.267~{}\mathrm{GeV}^{2},$
$\displaystyle b_{0}=7.989\times 10^{-6}~{}{\rm
GeV}^{2},\;\;b_{1}=-3.025\times 10^{-2},\;\;b_{2}=3.011\times 10^{-2}~{}{\rm
GeV}^{-2},$ $\displaystyle b_{0}^{\prime}=-1.438\times 10^{-6}~{}{\rm
GeV}^{2},\;\;b_{1}^{\prime}=6.680\times
10^{-7},\;\;b_{2}^{\prime}=-8.870\times 10^{-8}~{}{\rm GeV}^{-2}.$ (10)

Compared to Eqs. (8) and (9), the result of $m^{2}$ varies slightly, the
dominant coefficient $b_{1}$ changes by 10% roughly, and $d$ exhibits about
30% uncertainty. The above parameters, together with those for the $V-A$
contribution, lead to $|q/p|-1\approx 3.2\times 10^{-4}$ and
$Arg(q/p)\approx(7.1\times 10^{-3})^{\circ}$. It is seen that our predictions
for $|q/p|-1$ and $Arg(q/p)$ are quite stable with respect to the variation of
$m_{s}$, which change by only $\sim$10%-30%. That is, $q/p$ can serve as an
ideal observable for constraining new physics effects.

This work represents the first successful quantitative attempt in the sense
that definite values have been presented for the $D$ meson mixing parameters
$x$ and $y$ in both the CP-conserving and CP-violating cases in the Standard
Model. The key is to transform the dispersion relation between $x$ and $y$
into an inverse problem, in which the nonperturbative observables at low mass
are solved with the perturbative inputs from high mass. It is nontrivial to
find a solution under the analyticity constraint from the perturbative inputs
that explains the data of $x$ and $y$. The accommodation of the data by a
single resonance around the charm mass hints that multi-particle channels of
$D$ meson decays give dominant contributions to $y$. If such a solution does
not exist, it would be a strong indication that the large mixing parameters
are attributed to new physics. The obtained solution has been employed to
predict the coefficient ratio $q/p$ in the CP-violating case. To improve the
precision of our results, high-power corrections to the inputs can be included
systematically. Theoretical uncertainties in this approach will be
investigated in detail in the future. Once the $D$ meson mixing is understood,
relevant data, especially those for the coefficient ratio $q/p$, can be used
to constrain new physics effects appearing in the box diagrams. Our approach
will be developed into a fundamental nonperturbative QCD formalism, with the
insight that resonance properties are extractable from asymptotic QCD as
demonstrated in the $D$ meson mixing case.

Acknowledgement

This work was supported in part by MOST of R.O.C. under Grant No.
MOST-107-2119-M-001-035-MY3, and by NSFC under Grant Nos. U1932104, 11605076,
U173210, and 11975112.

## References

  * (1) H. Y. Cheng, Phys. Rev. D 26, 143 (1982).
  * (2) A. J. Buras, W. Slominski and H. Steger, Nucl. Phys. B245, 369 (1984).
  * (3) A. Datta and D. Kumbhakar, Z. Phys. C 27, 515 (1985).
  * (4) H. Georgi, Phys. Lett. B 297, 353-357 (1992).
  * (5) T. Ohl, G. Ricciardi and E. H. Simmons, Nucl. Phys. B 403, 605-632 (1993).
  * (6) E. Golowich and A. A. Petrov, Phys. Lett. B 625, 53 (2005).
  * (7) I. I. Bigi and N. G. Uraltsev, Nucl. Phys. B 592, 92-106 (2001).
  * (8) A. F. Falk, Y. Grossman, Z. Ligeti and A. A. Petrov, Phys. Rev. D 65, 054034 (2002).
  * (9) M. Bobrowski, A. Lenz, J. Riedl and J. Rohrwild, [arXiv:0904.3971 [hep-ph]]; JHEP 1003, 009 (2010).
  * (10) S. L. Glashow, J. Iliopoulos and L. Maiani, Phys. Rev. D 2, 1285 (1970).
  * (11) T. Jubb, M. Kirk, A. Lenz and G. Tetlalmatzi-Xolocotzi, Nucl. Phys. B 915, 431-453 (2017).
  * (12) L. Wolfenstein, Phys. Lett. B 164, 170-172 (1985).
  * (13) J. F. Donoghue, E. Golowich, B. R. Holstein and J. Trampetic, Phys. Rev. D 33, 179 (1986).
  * (14) P. Colangelo, G. Nardulli and N. Paver, Phys. Lett. B 242, 71-76 (1990).
  * (15) F. Buccella, M. Lusignoli, G. Miele, A. Pugliese and P. Santorelli, Phys. Rev. D 51, 3478-3486 (1995).
  * (16) T. A. Kaeding, Phys. Lett. B 357, 151-155 (1995).
  * (17) A. F. Falk, Y. Grossman, Z. Ligeti, Y. Nir and A. A. Petrov, Phys. Rev. D 69, 114021 (2004).
  * (18) H. Y. Cheng and C. W. Chiang, Phys. Rev. D 81, 114020 (2010).
  * (19) M. Gronau and J. L. Rosner, Phys. Rev. D 86, 114029 (2012).
  * (20) H. Y. Jiang, F. S. Yu, Q. Qin, H. n. Li and C. D. Lü, Chin. Phys. C 42, 063101 (2018).
  * (21) R. Kingsley, S. Treiman, F. Wilczek and A. Zee, Phys. Rev. D 11, 1919 (1975).
  * (22) M. Artuso, G. Borissov and A. Lenz, Rev. Mod. Phys. 88, 045002 (2016).
  * (23) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * (24) M. A. Shifman, hep-ph/0009131.
  * (25) Y. S. Amhis et al. [HFLAV], [arXiv:1909.12524 [hep-ex]].


# The standard model of star formation applied to massive stars: accretion
disks and envelopes in molecular lines

Eric Keto1 and Qizhou Zhang1
1Harvard-Smithsonian Center for Astrophysics, 160 Garden St, Cambridge, MA
02420, USA E-mail<EMAIL_ADDRESS>(EK<EMAIL_ADDRESS>(QZ)

(May 8, 2009)

###### Abstract

We address the question of whether the formation of high-mass stars is similar
to or differs from that of solar-mass stars through new molecular line
observations and modeling of the accretion flow around the massive protostar
IRAS20126+4104. We combine new observations of NH3 (1,1) and (2,2) made at the
Very Large Array, new observations of CH3CN(13-12) made at the Submillimeter
Array, previous VLA observations of NH3(3,3), NH3(4,4), and previous Plateau
de Bure observations of C34S(2-1), C34S(5-4), and CH3CN(12-11) to obtain a
data set of molecular lines covering 15 to 419 K in excitation energy. We
compare these observations against simulated molecular line spectra predicted
from a model for high-mass star formation based on a scaled-up version of the
standard disk-envelope paradigm developed for accretion flows around low-mass
stars. We find that in accord with the standard paradigm, the observations
require both a warm, dense, rapidly-rotating disk and a cold, diffuse
infalling envelope. This study suggests that accretion processes around 10 M⊙
stars are similar to those of solar mass stars.

###### keywords:

Keywords.

##  1 Introduction

Does the formation of massive stars differ significantly from that of solar
mass stars? As far as we know, stars of all masses form in gravitationally
unstable regions of molecular clouds and gain their mass by accretion. A
standard model developed for accretion flows around low-mass stars consists of
two-components, a rotationally-supported disk inside a freely-falling envelope
(Shu, Adams, Lizano, 1987; Hartmann, 2001). This model has been particularly
successful in explaining infrared observations of low-mass star formation. The
disk and envelope produce an excess of long-wavelength infrared emission that
has been adopted as the identifying signature of accreting protostars in
Galactic (Whitney et al., 2003) and extragalactic (Whitney et al., 2008) star-
forming regions. Furthermore, because the disk and envelope have different
densities and temperatures, the evolutionary state of the protostars can be
identified by the shape of the infrared spectral energy distribution: class 0
(envelope dominated) and class I (disk dominated) (Lada, 1987; Andre et al.,
1993).

Does this standard two component accretion model developed for low-mass stars
also describe the accretion flows around massive stars? There are some doubts.
The more massive stars are luminous enough to generate radiation pressure and
hot enough to ionize their own accretion flows such that the outward radiative
and thermal pressures rival the inward pull of the stellar gravity (Larson &
Starrfield, 1971; Kahn, 1974; Keto, 2002; Keto & Wood, 2006). Do these outward
pressures result in accretion flows that are different around more massive
stars? In this paper we compare new and previous molecular line observations
of an accretion flow around one massive star against the standard disk-
envelope paradigm for accretion flows around low-mass stars.

The previous molecular line observations of the massive protostar
IRAS20126+4104 111 IRAS20126 is located in the Cygnus-X region at a distance
of 1.7 kpc (Wilking et al. 1989). , suggest an accretion disk and bipolar
outflow around a 7 to 15 M⊙ protostar embedded in a dense molecular envelope
(Cesaroni et al., 1997; Zhang, Hunter, & Sridharan, 1998; Cesaroni et al.,
1999; Kawamura et al., 1999; Zhang et al., 1999; Shepherd et al., 2000;
Cesaroni et al., 2005; Lebron et al., 2006; Su et al., 2007; Qiu et al.,
2008). Previous infrared observations of absorption and scattering also reveal
a disk and outflow cavity immediately around the star (Sridharan, Williams, &
Fuller, 2005; deBuizer, 2007).

In this study, we assemble a suite of observations of molecular lines of
different excitation temperature in order to compare with molecular spectra
predicted from the disk-envelope model. Lines of different excitation
temperature are useful because massive stars heat the surrounding molecular
gas to observationally significant temperatures ($>100K$) at observationally
significant distances from the star, ($T_{gas}\sim
T_{\ast}(R/R_{\ast})^{\alpha}$), and we can exploit the relationship between
temperature and radius to distinguish emission from gas at different radii
around the star. We expect to identify the emission from the higher excitation
temperature lines with gas in the flow closer to the star and also to separate
the emission from the disk and envelope components. In previous observations,
Keto, Ho & Haschick (1987) and Cesaroni et al. (1994) used this technique,
observing several lines of NH3 to study the accretion flows around very high
mass stars associated with HII regions.

We present new observations of the NH3 (1,1) and (2,2), inversion transitions
made with the National Radio Astronomy Observatory’s Very Large Array
(VLA)222The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated
Universities, Inc. and new observations of CH3CN(13-12) made with the
Submillimeter Array (SMA). The new observations of NH3(1,1) and NH3(2,2) have
a factor of 3 better sensitivity than the earlier observations of Zhang,
Hunter, & Sridharan (1998). Combined with previous observations of NH3(3,3)
and NH3(4,4) (Zhang et al., 1999) the 4 NH3 lines span a range of excitation
temperatures from 23 K to 200 K. We also have additional lines of lower and
higher excitation temperature from previous observations of the C34S(2-1) and
C34S(5-4) lines with energies of 7 K and 35 K respectively, and the ladder of
CH3CN (J=12-11) lines with energies ranging from 69 K for (K=0) to 419 K for
(K=7). The C34S and CH3CN observations were previously presented in Cesaroni
et al. (1999) and Cesaroni et al. (2005).

We specify the two-component accretion model in terms of 6 parameters: the
scale factors for the density and temperature of the envelope and of the disk,
the angular momentum of the envelope, and the stellar mass. We use our
molecular line emission code MOLLIE to predict molecular line spectra from the
parameterized model, and we use a least-squares fitting procedure to adjust
the parameters to fit the observations.

We find that the two-component, disk-envelope model can successfully describe
the observations, but a single component model cannot. A warm, dense
rotationally-supported disk is required to obtain sufficient brightness and
width in the high-excitation lines, and a cold, large-scale envelope is
required to match the emission from the lower excitation lines. We find no
evidence that the accretion flow around IRAS20126 is profoundly altered by the
outward force of radiation pressure or by ionization. At 10 M⊙, the star might
simply not be luminous enough or hot enough for its radiation pressure or
ionization to significantly affect the accretion process. Based on this
example, the accretion flows around 10 M⊙ stars are quite similar to flows
around lower-mass stars.

##  2 The Data

###  2.1 Ammonia inversion lines

We observed the IRAS 20126+4104 region, R.A. (J2000) = 20:14:26.06, decl.
(J2000) = 41:13:31.50, in the inversion transitions of NH3 (J,K) = (1,1),
(2,2), (3,3) and (4,4) on 1999 March 27 and 1999 May 29 with the VLA in its D
configuration. We used the 2IF correlator mode to sample both the right and
left polarizations, a spectral bandwidth of 3.13 MHz, and a channel width of
24 kHz or 0.3 kms-1. The primary beam of the VLA at the NH3 line frequencies
is about $2^{\prime}$. Quasars 3C48 and 3C286 were used for flux calibration,
3C273 and 3C84 for bandpass calibration, and quasar 2013+370 for phase
calibrations. The NH3(3,3) and NH3(4,4) observations were previously discussed
in Zhang et al. (1999).

The VLA data were processed using the NRAO Astronomical Image Processing
System (AIPS) package. Data from the two days were combined to achieve an rms
of 3 mJy per 0.3 km s-1 channel for the NH3 (1,1), (2,2) lines and 6 mJy and 4
mJy per 0.6 km s-1 channel for the NH3 (3,3) and (4,4) lines respectively.
Figure 1 shows two images of IRAS20126 in the integrated intensity of NH3(1,1)
and (2,2). The protostar is located within the bright core of the larger scale
cloud. Also shown in figure 1 is the image of the bipolar outflow in SiO(2-1)
from Cesaroni et al. (1999). The NH3 spectra used in the analysis are taken at
positions offset from the phase center by -0.23, 1.00 (center), -3.10, 3.25
(left), and 2.50,-0.94 (right) in arc seconds of RA and dec. The corresponding
absolute positions are R.A. (J2000) = 20:14:26.02, decl. (J2000) = 41:13:32.8
(center) , R.A. (J2000) = 20:14:25.79, decl. (J2000) = 41:13:30.8 (right) ,
and R.A. (J2000) = 20:14:26.28, decl. (J2000) = 41:13:35.0 (left) , The
locations of the spectra are shown as red stars on figure 1.

###  2.2 Methyl Cyanide and Carbon Sulfide lines

The CH3CN(13-12) observations were made with the SMA on July 10, 2006 with a
maximum baseline of $\sim 0.5$ km, a beam size of $0.40\times 0.36$”, and a
phase center of R.A. (J2000) = 20:14:26.02, decl. (J2000) = 41:13:32.7. The
channel width is 0.8 kms-1 at the observing frequency of 239 GHz. The rms
noise of the data is 0.05 Jy beam-1 channel-1 or 8 K channel-1 in brightness
temperature.

Data in the CH3CN J = 12-11, and C34S J = 5-4 and J = 2-1 transitions were
obtained from the IRAM interferometer between 1996 December through 1997
March, and 2002 January through March. The phase center of the observations is
R.A. (J2000) = 20:14:26.03, decl. (J2000) = 41:13:32.7. The primary beam of
the array at the frequency of the observed CH3CN transitions (230 GHz) and
C34S(5-4) (241 GHz) is approximately $20^{\prime\prime}$ with a spatial
resolution of about $0.8^{\prime\prime}$. At the lower frequency of the
C34S(2-1) transtion (96 GHz) the primary beam is about $50^{\prime\prime}$ and
the spatial resolution is and $2.4^{\prime\prime}$. The full details of these
observations can be found in Cesaroni et al. (1999, 2005).

##  3 The Standard Model for Star Forming Accretion Flows

###  3.1 The density of the accretion flow

The standard model for a star-forming accretion flow consists of an inflowing
envelope around a rotationally supported disk. Similar to the approach in
Whitney et al. (2003), we use the model of Ulrich (1976) to describe the gas
density in the infalling envelope and standard thin-disk theory (Pringle,
1981) to describe the disk,

The model for the envelope assumes that the gas flows toward a gravitational
point source along ballistic trajectories. The flow conserves angular momentum
and ignores the pressure and self-gravity of the gas. This is the same
envelope model that Whitney et al. (2003) adopted for low-mass stars and
similar to the model used to describe the accretion flow onto a cluster of
high-mass stars, G10.6-0.07 (Keto, 1990; Keto & Wood, 2006). Thus despite, or
perhaps because of, its simplicity, the model has found application to
accretion flows from low-mass stars to star clusters. The model is
particularly useful in analyzing observational data because only 3 parameters
are required to describe individual cases: 1) the density of the gas at a
single radius (the gas density elsewhere follows from mass conservation), 2)
the mass of the point source, and 3) the specific angular momentum.

The model for the envelope is fully described in Ulrich (1976). We use the
equations as presented in Keto (2007), and Mendoza, Canto, Raga (2004). The
envelope density is (equation 10 of Keto (2007)),

$\displaystyle\rho_{\rm
env}=\rho_{e0}(r/R_{d})^{-3/2}\bigg{(}1+{{\cos{\theta}}\over{\cos{\theta_{0}}}}\bigg{)}^{-1/2}$
$\quad\quad\quad\times\quad(1+(r/R_{d})^{-1}(3\cos^{2}{\theta_{0}}-1))^{-1}$
(1)

where $\rho_{e0}=\rho_{\rm env}(R_{d})$ is the density in the mid-plane,
$\theta=\pi/2$, at radius, $R_{d}$, $\theta_{0}$ is the initial polar angle of
the streamline, and $\theta$ and $r$ are the polar angle and radius at each
point along a streamline. The angle, $\theta$ is related to the polar radius,
$r$, and the initial angle $\theta_{0}$ (equation 7 of Keto (2007)),

$r={R_{d}\cos{\theta_{0}}\sin^{2}{\theta_{0}}\over{\cos{\theta_{0}}-\cos{\theta}}}$
(2)

The density is related to the mass accretion rate (equation 11 of Keto (2007))

$\dot{M}=\rho_{e0}(R_{d})4\pi R_{d}^{2}v_{k}$ (3)

where $v_{k}$ is the Keplerian velocity at $R_{d}$.

The model for the disk assumes that a rotationally-supported disk forms at the
radius where the centrifugal force in the rotating envelope equals the
gravitational force of the point mass,

${{\Gamma^{2}}\over{R^{3}_{d}}}={{GM_{\ast}}\over{R^{2}_{d}}}$ (4)

and $\Gamma=v_{\phi}R$. The disk is truncated at $R_{d}$.

The gas density in the disk is (equation 3.14 of Pringle (1981)),

$\rho_{\rm disk}(z,R)=\rho_{0}(R)\exp(-z^{2}/2H^{2})$ (5)

where $\rho_{0}(R)$ is the density in the mid-plane, $\theta=\pi/2$, at any
radius $R$. We use $R=\sqrt{x^{2}+y^{2}}$ to denote the cylindrical radius,
and $r=\sqrt{x^{2}+y^{2}+z^{2}}$ for the polar radius. In the thin disk
theory, the scale height, $H$, is

$H^{2}=c_{s}^{2}R^{3}/GM_{\ast}$ (6)

where $c_{s}$ is the sound speed. We follow Whitney et al. (2003) and use a
modified version of this equation,

$H=H_{0}(R/R_{\ast})^{1.25}$ (7)

with $H_{0}=0.01R_{\ast}$. The density in the mid-plane of the accretion disk,
$\rho_{0}(R)$, is

$\rho_{0}=\rho_{d0}({{R_{d}}/{R}})^{2.25}$ (8)

where $\rho_{d0}$ is the density in the mid-plane at $R_{d}$. Again following
Whitney et al. (2003), we use an exponent of 2.25 in equation 8 rather than 2
which is derived in the thin-disk theory. As explained in Whitney et al.
(2003) the modifications to the equations for the scale height and mid-plane
densities are based on fits to numerical models of disk structure.

The densities in the disk and envelope at radius $R_{d}$ are related by a
factor, $A_{\rho}$,

$\rho_{d0}=A_{\rho}\rho_{e0}$ (9)

For example, in the steady-state flow, the disk would have a higher density
if, as should be the case, the inward velocities in the rotationally-supported
disk were lower than in the freely-falling envelope. Because the inward
velocities in the disk are not described by the thin-disk theory, we leave
$A_{\rho}$ as an adjustable parameter. The total gas density at any point in
the model is the sum of the densities in the disk and envelope,

$\rho_{\rm gas}=\rho_{\rm disk}+\rho_{\rm env}$ (10)

###  3.2 The temperature

We assume that the envelope is heated by the star. The envelope temperature is
(equation 7.36 of Lamers & Cassinelli (1999))

$T_{\rm env}=T_{\ast}(R_{\ast}/2r)^{2/(4+p)}$ (11)

where $p$ in our model is an adjustable parameter related to the dust opacity
and the geometry of the flow. If the density structure were spherically
symmetric, then $p$ would be the exponent in the frequency dependence,
$\nu^{p}$, of the Planck mean opacity of the dust. In a flattened flow, $p$
can be negative if the geometrical dilution of the radiation in the accretion
flow is greater than $r^{-2}$. This can happen if the dust at each radius
absorbs and isotropically re-emits the outward flowing radiation. Radiation
that is emitted perpendicular to the disk escapes from the flow. Only the
radiation that is emitted in the direction along the flattened flow continues
to heat the dust at larger radii. The result is a decrease in the radiation in
the disk faster than $r^{-2}$.

In the thin-disk theory, the disk is heated by dissipation related to the
accretion rate. The disk temperature is (equation 3.23 of Pringle (1981))

$T_{\rm disk}=B_{T}\bigg{[}\bigg{(}{{3GM_{\ast}\dot{M}}\over{4\pi
R^{3}\sigma}}\bigg{)}\bigg{(}1-\sqrt{{{R_{\ast}}\over{R}}}\bigg{)}\bigg{]}^{1/4}$
(12)

We include an adjustable factor, $B_{T}$, to allow for additional, or possibly
less, disk heating. For example, the observations constrain the gas density
through the observed optical depth of NH3. This means that the density and
therefore the accretion rate, $\dot{M}$, are dependent on the assumed NH3
abundance which is not well known. The adjustable factor, $B_{T}$, decouples
the disk temperature from the assumed molecular abundance. Also the disk
temperature may be raised by stellar radiation (passive heating). The factor,
$B_{T}$, allows for these effects in an approximate way.

The gas temperature in the model is the density weighted average of the
envelope and disk temperatures,

$T_{\rm gas}={{T_{\rm disk}\rho_{\rm disk}+T_{\rm env}\rho_{\rm
env}}\over{\rho_{\rm disk}+\rho_{\rm env}}}$ (13)

###  3.3 The velocity

The 3 components of the gas velocity in the envelope are given in spherical
coordinates by equations 4,5,6 of Keto (2007), but there are errors in earlier
papers. In particular, equations 5 and 6 of Keto (2007) and equation 8 of
Ulrich (1976) are incorrect. Equation 8 of Ulrich (1976) (same as equation 5
of Keto (2007)) does not result in energy conservation, $mv^{2}/2=GM/r$, when
combined with the other velocity components. The error is small, 1 part in
$10^{4}$. The following equations, same as in Mendoza, Canto, Raga (2004), are
exact.

$v_{r}(r,\theta)=-\bigg{(}{{GM_{\ast}}\over{r}}\bigg{)}^{1/2}\bigg{(}1+{{\cos{\theta}}\over{\cos{\theta_{0}}}}\bigg{)}^{1/2}$
(14)
$v_{\theta}(r,\theta)=\bigg{(}{{GM_{\ast}}\over{r}}\bigg{)}^{1/2}\bigg{(}{{\cos{\theta_{0}}-\cos{\theta}}\over{\sin\theta}}\bigg{)}\bigg{(}1+{{\cos\theta}\over{\cos{\theta_{0}}}}\bigg{)}^{1/2}$
(15)
$v_{\phi}(r,\theta)=\bigg{(}{{GM_{\ast}}\over{r}}\bigg{)}^{1/2}{\sin{\theta_{0}}\over{\sin{\theta}}}\bigg{(}1-{{\cos{\theta}}\over{\cos{\theta_{0}}}}\bigg{)}^{1/2}$
(16)

The velocity in the disk is simply the Keplerian velocity,

$v_{\rm disk}(R)=\sqrt{GM_{\ast}/R}$ (17)

where the velocity in the disk is purely azimuthal. We assume that the radial
velocity in the rotationally-supported disk is comparatively small. The gas
velocity in the model is the density weighted average of the envelope and disk
velocities,

$\vec{v}_{\rm gas}={{\vec{v}_{\rm disk}\rho_{\rm disk}+\vec{v}_{\rm
env}\rho_{\rm env}}\over{\rho_{\rm disk}+\rho_{\rm env}}}$ (18)

where $\vec{v}_{\rm env}$ is given by equations 14 through 16 and
$\vec{v}_{\rm disk}$ is given by 17.

###  3.4 The density singularities in the accretion model

The density in the Ulrich model is singular in the mid-plane of the disk at
the centrifugal radius and also at the origin. These singularities are caused
by the convergence of the streamlines in the simple mathematical description
of the flow. This convergence is not expected on physical grounds because gas
pressure, neglected in the Ulrich model, would prevent it. We handle the
singularities in two ways. We define the computational grid to have an even
number of cells so that the centers of the middle cells are above and below
the midplane and around the origin. Second, we smooth the density in the
radial direction with a Gaussian with a width of $R_{d}/2$.

###  3.5 The model parameters

Based on the analyses of the previous observations cited in the introduction,
we assume that the disk-envelope is viewed edge-on. The model contains 6
adjustable parameters: 1) $\rho_{0}$ sets the density of the envelope
(equation 1) and the mass accretion rate (equation 3), 2) $p$ sets the
exponent of the power law decrease of the temperature in the envelope
(equation 11), 3) $\Gamma=v_{\phi}r$ is the specific angular momentum
(equation 4) of the envelope flow, 4) $M_{\ast}$ is the stellar mass 5)
$A_{\rho}$ sets the ratio (equation 9) of the disk density to the envelope
density at $R_{d}$, 6) $B_{T}$ is factor multiplying the disk temperature
(equation 12).

##  4 Fitting the model to the NH3 data

The disk enters into the model additively, and we can test for the presence of
a disk by fitting models to the data with and without the disk. In the first
case, with the disk, we adjust all 6 model parameters for both the envelope
and the disk, and in the second case, without the disk, we adjust only the
first 4 parameters for the envelope. The fitting is done independently in each
case, and the 4 envelope parameters are therefore different in the 2 cases.
The procedure for fitting the data is the same as described in Keto et al.
(2004). We use a fast simulated annealing algorithm to adjust the model
parameters to minimize the summed squared difference ($\chi^{2}$) between the
data and the model spectra.

The model spectra for each particular set of model parameters are generated by
our radiative transfer code MOLLIE (Keto, 1990; Keto et al., 2004). We assume
LTE conditions for the NH3 and CH3CN lines. The LTE approximation is
appropriate for NH3 because in the absence of strong infrared radiation the
level populations are expected to be mostly in the ”metastable” states which
are the lowest J-state of each K-ladder. Since radiative transitions between
K-ladders are forbidden, the coupling between the metastable states is purely
collisional and the population approximates a Boltzmann distribution as in LTE
(Ho & Townes, 1983). CH3CN is also a symmetric top and transitions across the
K-ladders are similarly forbidden; however, unlike NH3, the upper states are
easily populated in warm gas ($>100$ K). While the justification for the LTE
approximation is not as strong for CH3CN as for NH3, we find that CH3CN always
traces hot, very dense gas where collisional transitions should be important.
For C34S we use the accelerated lambda iteration algorithm (ALI) of Rybicki &
Hummer (1991) to solve for the non-LTE level populations.

In comparing the model to the data, we take into account the spatial averaging
of the brightness by the width of the observing beam. We compute spectra over
a grid of positions and smooth the result by convolution with a Gaussian with
the FWHM equivalent to the observing beam of each observation, §2.

We simultaneously fit 12 NH3 spectra, 4 transitions at 3 locations. The 3
locations are marked on figure 1. We ran twenty thousand trial models for each
of the two cases with and without the disk. The spectra of the 2 best-fit
models (with and without the disk) are shown in figures 2 and 3. Parameters
for both cases are listed in table 1. Figure 4 shows the temperature and
density in the mid-plane of our best fitting models. Figure 5 shows the
density and velocity of the disk-envelope model on planes parallel and
perpendicular to the rotation.

We also have CH3CN and C34S data from the observations of Cesaroni et al.
(1999). The way the radiative transfer simulation program operates, we cannot
use our automated search algorithm on more than one molecule at a time. The
program is recompiled for each molecule. We also have not implemented the
automated search for CH3CN. We could fit the model to the C34S data, but this
line does not have the hyperfine structure of NH3 that is so useful in
constraining the optical depth and temperature. Therefore we opted for a
different strategy. We do not use the CH3CN or the C34S data in the fitting.
We use the data on these other lines as a check on the model derived from the
NH3 data. We simulate the CH3CN and C34S emission using the same two models
(with and without the disk) previously derived from the NH3 observations for
comparison with the predicted spectra against the data. The models derived
from the NH3 fitting fix the temperature, densities, and velocities, but we
still need to assume abundances for CH3CN and C34S. For CH3CN we assume an
abundance of $6\times 10^{-8}$. The assumed abundance of C34S is $9.0\times
10^{-11}$, chosen to match the brightness of the (2-1) line.

###  4.1 Comparison of the model with the observed spectra

The comparison of figures 2 and 3 shows that the higher temperature and
density disk and the lower temperature and density envelope are both required
to fit both the high-excitation and low-excitation NH3 spectra. Without the
disk, there is not enough warm, dense gas to reproduce the (4,4) line
brightness. Without the large-scale envelope, there is not enough cold gas to
fit the observed ratios of the NH3(1,1) hyperfine lines.

Further evidence for the disk component comes from the CH3CN spectrum. The
results again show that the warm, high-density, disk component is required to
get enough brightness in lower frequency (higher velocity) K transitions.
Based on the observed line width, the CH3CN(12-11) emission comes from gas
very close to the star. The observed CH3CN(12-11) lines that are not blended
with other molecular lines, have line widths of 10.0 kms-1. The 2 lines that
appear much broader in the data contain contaminating emission from
transitions of CH13CN(12-11) and HNCO (Cesaroni et al., 1999). Rotational
velocities $>10$ kms-1 are only found at radii closer than 600 AU
($<GM/v^{2}$) around a 10 M⊙ star. Figure 4 shows that the gas in the disk has
a temperature of greater than 150 K within this radius. Thus the CH3CN
emission derives mostly from warm gas that is very close to the star and is
thus a good molecular tracer of the inner disk.

The CH3CN(13-12) spectra is too noisy to help discriminate between the two
cases. The peak signal-to-noise ratio is about 6 after smoothing by every
other channel. Nonetheless, the models are at least consistent with these
observations.

The brightness ratio of the C34S (5-4) and (2-1) lines also suggests the
presence of a warm, dense disk (figure 7). The model with the disk reproduces
the observed brightness ratio, 6:1 for the 2 lines, although line width is
less than observed. The model without the disk is not able to generate
sufficient brightness in the higher excitation line, and the shapes of the
line profiles do not match the data. In particular, the strong splitting that
is seen in the model profiles without the disk is not seen in the data. This
difference is also seen in the NH3 and CH3CN spectra of the 2 models, although
not as prominently. We usually associate split line profiles with spatially
unresolved observations of Keplerian disks (Beckwith & Sargent, 1993) whereas
here the disk produces a triangular profile. What happens here is that the
disk component is spatially resolved and has a very high density. Thus within
the beam through the center of the model there is a lot of high density gas
from the outer part of the disk with very low velocity projected along the
line-of-sight. This gas disk creates the peak in the spectrum around zero
velocity. In the envelope-only model without the high density disk, the model
is optically thin in C34S. In this case, we get the usual split line profiles
from the rotation and infall in the envelope.

###  4.2 Goodness of fit

How well do the data constrain the model parameters? Some appreciation can be
gained by plotting $\chi^{2}$ obtained from all the trial models versus each
model parameter. Each panel shows the fits obtained (ordinate) for each value
of a single parameter (abscissa) as all the other parameters are varied. The
lowest value of $\chi^{2}$ (ordinate) at each parameter value (abscissa) is
the best fit that can be obtained for that parameter value, for any
combination of all the other parameters. The formal error of a model parameter
is proportional to the second derivative of $\chi^{2}$ with respect to the
model parameter. Thus the curvature or width of the lower boundary of the
collection of points in each figure is a qualitative measure of the
sensitivity of the model to the parameter.

##  5 Discussion

###  5.1 The density

The gas density is constrained by the optical depths of the NH3 lines, which
can be determined from the brightness ratios of the hyperfine lines. The
density, $\rho_{0}$, at $R_{d}$ is $7.9\times 10^{4}$ (table 1) assuming an
NH3 abundance of $10^{-7}$. With this abundance, the mass of the disk is 2.5
M⊙, the mass of the envelope is 12.6 M⊙, and the accretion rate (equation 3),
asssumed to be the same in the envelope and the disk, is $\dot{M}=7.6\times
10^{-5}$ M⊙ yr-1. The total mass in the accretion flow is mass of the disk and
envelope together which is 15.1 M⊙. This mass estimate is the total mass
within the model boundary of 26000 AU. Cesaroni et al. (1999) estimates that
there is between 0.6 and 8.0 M⊙ within a radius of 5000 AU.

The mass estimate derived from molecular line observations is subject to a
large uncertainty. The radiative transfer modeling determines the column
density of NH3 rather than H2. Therefore, the gas density, the masses of the
disk and envelope, and the accretion rate depend inversely on the assumed
abundance of NH3 which is not well known. Estimates from models and
observations of similar clouds range from $10^{-9}$ to $10^{-6}$ (Herbst &
Klemperer, 2003; Keto, 1990; Estalella et al., 1993; Caproni et al., 2000;
Galvan-Madrid et al., 2009). Thus the abundance of NH3 may be uncertain by
more than an order of magnitude, yet a factor of 2 in the masses of the disk
and envelope is significant in an interpretation of the accretion dynamics.
Furthermore, the abundance of NH3 could be different in the disk and envelope.
Estimates of the NH3 abundance are generally lowest in colder clouds and
highest in warm gas around massive stars. Some NH3 may be frozen onto dust
grains in colder gas and sublimated into the gas phase as the temperature
rises. Therefore, it is possible that the NH3 abundance is higher in the warm
disk than the cold envelope. If so, the mass of the envelope could be higher
than 12.6 M⊙, assuming a lower abundance, without necessarily implying a
higher mass for the disk.

###  5.2 The disk density factor

In the model with the flared disk, the density in the mid-plane of the disk at
the radius of the disk boundary $R_{d}$ is 5 times more dense than the
smoothed density of the envelope at the same point. Considering that the gas
density changes by 6 orders of magnitude across the model, there is little
discontinuity between the envelope and disk densities (figure 4). In the model
without the disk the smoothed density is approximately constant from $R_{d}$
to the origin (figure 4).

###  5.3 The stellar mass

The stellar mass, 10.7 M⊙, is constrained primarily by the disk velocities
(equation 17 ) required to match the observed linewidths (about 10 kms-1) and
the gas temperature that determines the brightness ratios of the low and high
excitation lines. The gas velocities and therefore line widths due to
unresolved motions within the beam depend on the stellar mass because the
velocities are proportional to the square root of the mass. The gas
temperature depends on the mass of the star because the stellar temperature is
a strong function of the stellar mass and because the accretion rate depends
on the stellar mass. In the disk, heating by both the star and by accretion
(dissipation) are important (equations 11, 12 and 13 ). In the best-fit model,
the disk temperature would be about 1/3 lower without the passive heating from
the star, The temperature of the infalling envelope is determined entirely by
the stellar temperature (equation 11).

###  5.4 The angular momentum of the envelope

The angular momentum of the envelope is determined from the widths of the
spectral lines and from the VLSR of the NH3 spectra at the locations left and
right of the center. If there were too much rotation then both of the off-
center spectra would have an incorrect VLSR; too little, and all the NH3 line
widths would be too narrow. The disk radius, $R_{d}=6900$, is set by the
angular momentum in the envelope, $\Gamma$, and the stellar mass (equation 4).
The velocity at $R_{d}$, $v_{k}=1.2$ kms-1, increases inward as $R^{-1/2}$. We
assume a microturbulent broadening of 1 kms-1, added in quadrature to the
thermal broadening. Since the observed line widths are about 10 kms-1 most of
the width of the spectral lines comes from rotation and infall in the model.

The model radius of 6900 AU seems large, and we would regard this as an upper
limit. First, the Ulrich flow conserves angular momentum whereas real
accretion flows probably involve some braking of the spin-up. If the spin-up
were slower, then the disk would be found at a smaller radius. Second, the
high density of the disk is helpful in providing optical depth to strengthen
the NH3(1,1) outer hyperfines. The Ulrich flow has a density singularity in
the mid-plane which we have handled as desribed in §3.4 by a combination of
gridding and smoothing. Maybe our mid-plane density in the envelope is too
low, and to compensate, the thin disk is bigger.

###  5.5 The exponent of the temperature power law

The modelling suggests that the temperature of the envelope falls off very
quickly away from the star, as $r^{-1}$. This implies that the dilution of the
radiation is faster than spherical, $r^{-2}$. The rapid dilution is consistent
with the rotationally flattened geometry. At each radius, radiation is
absorbed and re-emitted by dust in the flow. In a flattened flow, much of this
reprocessed radiation escapes vertically out of the flow. In contrast, in a
spherical flow, all the reprocessed radiation continues to interact with the
flow at larger radii. If the gas temperature decreases outward fast enough
that most of the envelope is at the minimum gas temperature, then a larger
value of the exponent (more negative value of the parameter $p$ in equation
11) would have no further effect. We assume a minimum temperature of 10 K, a
typical temperature of molecular gas, We derive an upper limit $p<-1$.

###  5.6 The disk temperature multiplier

The disk temperature multiplier determines the relative importance of active
disk heating, owing to accretion and dissipation (equation 12), and passive
heating from the protostar (equation 11). With active disk heating, the
temperature in the disk is about 50% above that from passive heating alone.

##  6 Conclusions

This investigation shows that a standard model of a freely-falling,
rotationally-flattened accretion envelope around a rotationally supported disk
is able to reproduce the NH3, CH3CN, and C34S spectral line observations of
IRAS20126. Both the disk and envelope components are required to fit the
observed brightness of both the low and high excitation lines, and to fit the
observed line widths.

This disk-envelope model was developed for low-mass stars and is quite
successful in explaining many of their observable characteristics. The success
of this model in explaining the molecular line observations of the massive
star IRAS20126 suggests that at least up to 10 M⊙, the accretion processes of
massive stars are similar to those of solar mass stars.

There are some differences. Although the mass of the disk is uncertain owing
to its dependence on the NH3 abundance, the disk mass is a significant
fraction of the mass of the star. Furthermore, the extent of the disk is quite
large, $\sim$6900 AU. This suggests that self-gravity in the disk is
dynamically important, and therefore, the disk may be unstable to local
fragmentation and the formation of companion stars.

The physical model in combination with the molecular line radiative transfer
presented in this paper has a further application to other massive star
forming regions. The spatial resolution and sensitivity of the present
generation of interferometers cannot spatially resolve accretion disks in
massive protostellar objects at multi-kpc distances. However, spectral lines
sampling a wide range of densities and temperatures can still provide
constraints to the physical structure of the core and disk. As shown in
IRAS20126, the spectral lines formed in the higher density and temperature
regions confirm the presence of an accretion disk. While future observatories
such as ALMA will be able to spatially resolve the flow on the disk scale,
before the science commissioning of ALMA, this method is a promising technique
to probe the spatially unresolved regions of the flow and disk using data from
the current interferometers.

## Acknowledgments

The authors thank Riccardo Cesaroni and T.K. Sridharan for the use of their
observational data.

## References

  * Andre et al. (1993) Andre, P., Ward-Thompson, D., & Barsony, M., 1993, 406, 122
  * Caproni et al. (2000) Caproni, A., Abraham, Z., & Vilas-Boas, J.W.S., 2000, AA, 361, 685
  * Cesaroni et al. (1997) Cesaroni, R., Felli, M., Walmsley, C.M., & Olmi, L, 1997, AA, 325, 725
  * Cesaroni et al. (1999) Cesaroni, R., Felli, M., Jenness, T., Neri, R., Olmi, L., Robberto, M., Testi, L., Walmsley, C.M., 1999, AA, 345, 949
  * Cesaroni et al. (2005) Cesaroni, R., Neri, R., Olmi, L., Testi, L., Walmsley, C.M., & Hofner, P., 2005, AA, 434, 1039
  * Estalella et al. (1993) Estalella, R., Mauersberger, R., Torrelles, J.M., Anglada, G., Gomez, J.F., Lopez, R., Murders, D., 1993, ApJ, 419,698
  * deBuizer (2007) deBuizer J.M., 2007, ApJ, 654, L147
  * Galvan-Madrid et al. (2009) Galvan-Madrid, R., Keto, E., Zhang, Q., Kurtz, S., Rodriguez,L.F., Ho, P.T.P., 2009, ApJ, submitted
  * Herbst & Klemperer (2003) Herbst, E. & Klemperer, W., 1973, ApJ, 616, 301
  * Keto, Ho & Haschick (1987) Keto, E., Ho, P., & Haschick, A., 1987, ApJ, 318, 712
  * Cesaroni et al. (1994) Cesaroni, R., Churchwell, E., Hofner, P., Walmsley, C.M., & Kurtz, S., 1994, AA, 288, 903
  * Hartmann (2001) Hartmann, L., Accretion Processes in Star Formation Cambridge, UK: Cambridge University Press
  * Ho & Townes (1983) Ho, P.T.P., Townes, C.H., 1983, ARAA, 21, 239
  * Kahn (1974) Kahn, F.D., 1974, AA, 37, 149
  * Kawamura et al. (1999) Kawamura, J.H., Hunter, T.R. , Tong, C.-Y.E., Blundell, R., Zhang, Q., Katz, C., Papa, D.C. & T.K. Sridharan, 1999, PASP, 111, 1088
  * Keto (1990) Keto, E., 1990, ApJ, 355, 190
  * Keto (2002) Keto, E. 2002b, ApJ, 580, 980
  * Keto et al. (2004) Keto, E., Rybicki, G.B., Bergin, E.A., Plume, R., 2004, ApJ, 613, 355
  * Keto & Wood (2006) Keto, E. & Wood, K., 2006, ApJ, 637, 850
  * Keto (2007) Keto, E., 2007, ApJ, 666, 976
  * Kumar & Grave (2007) Kumar M.S.N. & Grave, J.M.C., 2007, AA, 472, 155
  * Lada (1987) Lada, C.J., in IAU Symp. 115, Star Forming Regions, ed. M Peimbert & J.Jugaku (Dordrecht:Kluwer), 1
  * Lamers & Cassinelli (1999) Lamers, H.J.G.L.M. & Cassinelli, J.P., 1999, Introduction to Stellar Winds Cambridge, UK: Cambridge University Press
  * Larson & Starrfield (1971) Larson, R.B. & Starrfield, S., 1971, AA, 13, 190
  * Lebron et al. (2006) Lebron, M., Beuther, H., Schilke, P. & Stanke, T. 2006, AA, 448, 1037
  * Mendoza, Canto, Raga (2004) Mendoza, S., Canto, J., Raga, A.C., 2004 Rev. Mex. AA, 40, 147
  * Pringle (1981) Pringle J., 1981, ARAA, 19, 137
  * Qiu et al. (2008) Qiu, Keping, Zhang, Qizhou, Megeath, S Thomas, Gutermuth, Robert A., Beuther, Henrik, Shepherd, Debra S., Testi, L., and De Pree, C. G. 2008, ApJ, 685, 1005
  * Rybicki & Hummer (1991) Rybicki, G.B. & Hummer, D.G., 1991, AA, 245, 171
  * Beckwith & Sargent (1993) Beckwith, S.V.W. & Sargent, A.I., 1993, ApJ, 402, 280
  * Shepherd et al. (2000) Shepherd, D. S., Yu, K. C., Bally, J. & Testi, L. 2000, ApJ, 535, 833
  * Shu, Adams, Lizano (1987) Shu, F.H., Adams, F.C., Lizano, S., 1987, ARAA, 25, 23
  * Sridharan, Williams, & Fuller (2005) Sridharan, T.K., Williams, S.J., & Fuller, G.A., 2005, ApJ, L73
  * Su et al. (2007) Su, Yu-Nung, Liu, Sheng-Yuan, Chen, Huei-Ru, Zhang, Qizhou, and Cesaroni, Riccardo, 2007, ApJ, 671, 571
  * Ulrich (1976) Ulrich J., 1976, ApJ, 210, 377
  * Whitney et al. (2003) Whitney, B.A., Wood, K., Bjorkman, J.E., Wolff, M.J., 2003, ApJ, 591, 1049
  * Whitney et al. (2008) Whitney, B.A., et al., 2008, AJ, 136, 18
  * Wilking, Lada, Young (2007) Wilking, Bruce A, Lada, Charles J., & Young, Eric T. 1989, ApJ, 340, 823
  * Wolfile & Cassinelli (1987) Wolfire, M.G., Cassinelli, J.P., 1987, ApJ, 319, 850
  * Zhang, Hunter, & Sridharan (1998) Zhang, Q., Hunter, T.R., Sridharan, T.K., 1998, ApJ, 505, L151
  * Zhang et al. (1999) Zhang, Q., Hunter, T. R., Sridharan, T. K., & Cesaroni, R., 1999, Astrophysical Journal, 527, L117

Table 1: Model Parameters. Parameter | Symbol | Disk | No Disk
---|---|---|---
Env. Density at $R_{d}$ (cm-3) | $\rho_{0}$ | $7.9\times 10^{4}$ | $7.9\times 10^{5}$
Temperature power law exp. | $p$ | $<-1$ | 0.4
Angular Momentum (AU kms-1) | $\Gamma$ | 8100 | 3500
Stellar Mass (M⊙) | M∗ | 10.7 | 7.3
Disk Density Ratio | $A_{\rho}$ | 5.1 | …
Disk Temperature factor | $B_{T}$ | 15.0 | …
Centrifugal radius (AU) | $R_{d}$ | 6900 | 1900
Velocity at $R_{d}$ (kms-1) | $v_{k}$ | 1.2 | 1.8
Total mass1 within 0.128 pc ($M_{\odot}$) |  | 12.6 | 10.1
Disk mass ($M_{\odot}$) |  | 2.5 | …
Envelope accretion rate ($M_{\odot}$ yr-1) |  | $7.6\times 10^{-5}$ | $1.0\times 10^{-4}$
1 The mass estimates include the envelope and the disk for the case |  |  |
with a disk, and the envelope only for no disk. |  |  |
Both cases assume an NH3 abundance of $1.0\times 10^{-7}$. |  |  |
Figure 1: Images of IRAS20126 in NH3(1,1), NH3(2,2), and SiO(2-1). The red
stars in the NH3 images show the locations of the NH3 spectra in figures 2 and
3. See text for absolute coordinates. The single star in the SiO image
corresponds to the middle star in the other images. The angular resolution
(FWHM) is indicated by the shaded ellipse, lower-left corner of each panel.
The contour interval and first contour is 35 mJy kms-1 for the NH3(1,1) and
NH3(2,2) emission, and 130 mJy kms-1 for the SiO(2-1). The SiO image is from
Cesaroni et al. (1999).  Figure 2: Data (red) and model (blue) spectra of
NH3(1,1) (top row) NH3(2,2) (2nd row) NH3(3,3) (3rd row) and NH3(4,4) (bottom
row) at 3 positions across the mid-plane of IRAS20126 for the accretion flow
including both an infalling envelope and a flared disk. The spectra in the
middle column are toward the center of the flow, and the the columns to the
left and right show spectra on either side of the center, 6200 AU to the left,
and 5700 au to the right. The 3 locations are shown in figure 1.

Figure 3: Spectra in the same format as figure 2 for the accretion flow with
an envelope only without a flared disk.

Figure 4: The temperature (red) and density (blue) for the model with both an
envelope and flared disk (solid line) and for the model without a disk (dashed
line). The temperatures are about the same in both models, but the disk adds
considerable density. The increase in density in the envelope-only model
(dashed) at 3300 AU is the smoothed density singularity at $R_{d}$ in the
Ulrich (1976) model.

Figure 5: Density (color) and velocities (arrows) for the model with both an
envelope and flared disk. Left: The XZ plane (perpendicular to the rotation)
and right: the XY plane (the mid-plane of the disk and rotationally-flattened
envelope). In the left panel, the color scale is logarithmic between
$9.2\times 10^{4}$ and $3.7\times 10^{6}$ cm-3. Contours (white) indicate the
density at 0.006, 0.01, and 0.03, of the peak density. The longest arrow
corresponds to 2.5 kms-1 In the right panel the color scale is logarithmic
between $1.1\times 10^{5}$ and $8.1\times 10^{7}$ cm-3. Contours (white) show
the density at 0.03, 0.10, and 0.3 of the peak density. The linear scale in
both panels is $R_{d}=6900$ AU.

Figure 6: Data (red) and model (blue) spectra of CH3CN(12-11). The left panel
shows the spectrum for the model with both an accretion envelope and flared
disk. The model is the same as shown in figures 2, 4, 5. The right panel is
for the model with only an accretion envelope and without a disk (figures 3,
4). In the observed spectrum, the line that is anomalously wide (K=6, 2nd from
right) is contaminated by HNCO emission. Additionally, contamination from weak
CH${}_{3}^{13}$CN emission shows up to the red (right) of 4 of the lines. The
velocity of the K=6 hyperfine line has been set to zero. The observed spectra
have been shifted in velocity to match.

Figure 7: C34S(2-1) (blue) and C34S(5-4)(red) spectra of the data (solid
lines) and the model (dashed lines). Left: Model with both the accretion
envelope and flared disk. Right: Model with only an accretion envelope and
without a flared disk. The brightness of the C34S(5-4) spectra, both observed
and modeled, have been divided by 6. The comparison shows that the model with
the disk produces the correct brightness ratio. Without the disk, the
C34S(5-4) line is not bright enough and both lines have the wrong profile
shape. The zero velocity refers to the model. The observed spectra have been
shifted to match.

Figure 8:  Spectra of CH3CN(13-12) in the same format as figure 6 except that
we have only the 4 lowest K transitions. The K=0 and K=1 transitions are
blended together in the bluest peak. The left panel shows the spectrum for the
model with both an accretion envelope and flared disk. The right panel is for
the model with only an accretion envelope and without a disk (figures 3, 4).
The signal-to-noise ratio of the observed spectra, about 6 at the peak, is not
high enough to discriminate between the two models. At least the models are
consistent with the data. The data are from T.K. Sridharan (private
communication).

Figure 9: The value of $\chi^{2}$ as a function of a single parameter. The
plots show $\chi^{2}$ for the better fitting disk-envelope models as a
function of a single parameter. The width of the collection of points in each
figure is a qualitative measure of the sensitivity of the models to each
parameter. Values of the parameters outside the range shown by the lower
$\chi^{2}$ values produce very poor fits. The units of the density are the log
of $cm^{-3}$, the units of stellar mass are $M_{\odot}$, and the units of
angular momentum are AU $\times$ kms-1 divided by 1000. The parameters are
listed in Table 1.