Analysis of the decay $Y(4500)\to D^{*}\bar{D}^{*}\pi$ with the light-cone QCD
sum rules

Zhi-Gang Wang 111E-mail<EMAIL_ADDRESS>
Department of Physics, North China Electric Power University, Baoding 071003,
P. R. China

###### Abstract

In this work, we tentatively assign the $Y(4500)$ as the
$[uc]_{\tilde{A}}[\overline{uc}]_{V}+[uc]_{V}[\overline{uc}]_{\tilde{A}}+[dc]_{\tilde{A}}[\overline{dc}]_{V}+[dc]_{V}[\overline{dc}]_{\tilde{A}}$
tetraquark state with the quantum numbers $J^{PC}=1^{--}$, and study the
three-body strong decay $Y(4500)\to D^{*-}D^{*0}\pi^{+}$ with the light-cone
QCD sum rules. It is the first time to use the light-cone QCD sum rules to
calculate the four-hadron coupling constants, the approach can be extended to
study other three-body strong decays directly and diagnose the $X$, $Y$ and
$Z$ states.

PACS number: 12.39.Mk, 12.38.Lg

Key words: Tetraquark state, QCD sum rules

##  1 Introduction

In last two decays, several vector charmonium-like states have been observed,
they cannot be accommodated comfortably in the traditional charmonia, we have
to introduce additional quark or gluon degrees of freedom in assignments [1].
For example, the $Y(4260)$ observed in the $J/\psi\pi^{+}\pi^{-}$ invariant
mass spectrum by the BaBar collaboration [2], the $Y(4220)$ and $Y(4390)$
($Y(4320)$) observed in the $h_{c}\pi^{+}\pi^{-}$ ($J/\psi\pi^{+}\pi^{-}$)
invariant mass spectrum by the BESIII collaboration [3, 4], and the $Y(4360)$
and $Y(4660)$ ($Y(4630)$) observed in the $\psi^{\prime}\pi^{+}\pi^{-}$
($\Lambda_{c}^{+}\bar{\Lambda}_{c}^{-}$) invariant mass spectrum by the Belle
collaboration [5, 6, 7] are excellent candidates for the vector tetraquark
states.

In 2022, the BESIII collaboration explored the $e^{+}e^{-}\to
K^{+}K^{-}J/\psi$ cross sections at center-of-mass energies from 4.127 to
4.600 GeV based on $15.6\,\rm{fb}^{-1}$ data, and observed two resonant
structures, one is consistent with the established $Y(4230)$; the other was
observed for the first time with a significance larger than $8\sigma$ and
denoted as $Y(4500)$, its Breit-Wigner mass and width are $4484.7\pm 13.3\pm
24.1\,\rm{MeV}$ and $111.1\pm 30.1\pm 15.2\,\rm{MeV}$, respectively [8].

Recently, the BESIII collaboration explored the Born cross sections of the
process $e^{+}e^{-}\to D^{*-}D^{*0}\pi^{+}$ at center-of-mass energies from
4.189 to 4.951 GeV using the data samples corresponding to an integrated
luminosity of $17.9\,\rm{fb}^{-1}$ and observed three enhancements, whose
masses are $4209.6\pm 4.7\pm 5.9\,\rm{MeV}$, $4469.1\pm 26.2\pm 3.6\,\rm{MeV}$
and $4675.3\pm 29.5\pm 3.5\,\rm{MeV}$, respectively, and widths are $81.6\pm
17.8\pm 9.0\,\rm{MeV}$, $246.3\pm 36.7\pm 9.4\,\rm{MeV}$ and $218.3\pm 72.9\pm
9.3\,\rm{MeV}$, respectively. The first and third resonances are consistent
with the $Y(4230)$ and $Y(4660)$ states, respectively, while the second
resonance is compatible with the $Y(4500)$ [9].

In fact, analogous decays were already observed in the process $e^{+}e^{-}\to
Y\to\pi^{+}D^{0}D^{*-}$ for the center-of-mass energies from 4.05 to 4.60 GeV
by the BESIII collaboration in 2018, and the two enhancements $Y$ lie around
4.23 and 4.40 GeV, respectively [10].

In the scenario of tetraquark states, the calculations based on the QCD sum
rules have given several reasonable assignments of the $Y$ states [11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22]. For example, in Ref.[22], we take the
scalar, pseudoscalar, axialvector, vector and tensor (anti)diquarks to
construct vector and tensor four-quark currents without introducing explicit
P-waves, and explore the mass spectrum of the vector hidden-charm tetraquark
states via the QCD sum rules in a comprehensive way. At the energy about
$4.5\,\rm{GeV}$, we obtain three hidden-charm tetraquark states with the
$J^{PC}=1^{--}$, the tetraquark states with the symbolic structures
$[uc]_{\tilde{V}}[\overline{dc}]_{A}-[uc]_{A}[\overline{dc}]_{\tilde{V}}$,
$[uc]_{\tilde{A}}[\overline{dc}]_{V}+[uc]_{V}[\overline{dc}]_{\tilde{A}}$ and
$[uc]_{S}[\overline{dc}]_{\tilde{V}}-[uc]_{\tilde{V}}[\overline{dc}]_{S}$ have
the masses $4.53\pm 0.07\,\rm{GeV}$, $4.48\pm 0.08\,\rm{GeV}$ and $4.50\pm
0.09\,\rm{GeV}$, respectively [22]. Thus we have three candidates for the
$Y(4500)$, the best assignment of the symbolic structure is
$[uc]_{\tilde{A}}[\overline{dc}]_{V}+[uc]_{V}[\overline{dc}]_{\tilde{A}}=Y(4500)$
comparing with the BESIII experimental data $M_{Y(4500)}=4469.1\pm 26.2\pm
3.6\,\rm{MeV}$ [9], where we have taken the isospin limit, the tetraquark
states with the valence quark structures,

$\displaystyle\bar{c}c\bar{d}u,\,\,\bar{c}c\bar{u}d,\,\,\bar{c}c\frac{\bar{u}u-\bar{d}d}{\sqrt{2}},\,\,\bar{c}c\frac{\bar{u}u+\bar{d}d}{\sqrt{2}}\,,$
(1)

have degenerated masses and pole residues. As there exist three four-quark
currents with the $J^{PC}=1^{--}$, which couple potentially to the vector
hidden-charm tetraquark states with almost degenerated masses [22], we can
also tentatively say that there only exists one vector tetraquark state with
three different Fock components
$[uc]_{\tilde{V}}[\overline{dc}]_{A}-[uc]_{A}[\overline{dc}]_{\tilde{V}}$,
$[uc]_{\tilde{A}}[\overline{dc}]_{V}+[uc]_{V}[\overline{dc}]_{\tilde{A}}$ and
$[uc]_{S}[\overline{dc}]_{\tilde{V}}-[uc]_{\tilde{V}}[\overline{dc}]_{S}$, we
can choose either Fock component (in other words, current) to explore the
hadronic properties. At the first step, we choose the optimal current
corresponding to the optimal Fock component
$[uc]_{\tilde{A}}[\overline{dc}]_{V}+[uc]_{V}[\overline{dc}]_{\tilde{A}}$,
i.e. the current $J_{\mu}^{Y}(0)$ in Eq.(2) in the isospin limit. However, we
cannot assign a hadron unambiguously with the mass alone, we have to explore
the decay width to make more robust assignment. If we want to investigate the
three-body strong decays $Y\to J/\psi\pi^{+}\pi^{-}$,
$\psi^{\prime}\pi^{+}\pi^{-}$, $J/\psi K^{+}K^{-}$, $h_{c}\pi^{+}\pi^{-}$,
$D^{0}D^{*-}\pi^{+}$ and $D^{*-}D^{*0}\pi^{+}$ with the QCD sum rules
directly, we have to introduce four-point correlation functions, the hadronic
spectral densities are complex enough to destroy the reliability of the
calculations.

In this work, we tentatively assign the $Y(4500)$ as the
$[uc]_{\tilde{A}}[\overline{uc}]_{V}+[uc]_{V}[\overline{uc}]_{\tilde{A}}+[dc]_{\tilde{A}}[\overline{dc}]_{V}+[dc]_{V}[\overline{dc}]_{\tilde{A}}$
tetraquark state with the $J^{PC}=1^{--}$, and extend our previous works to
study the three-body strong decay $Y(4500)\to D^{*-}D^{*0}\pi^{+}$ with the
light-cone QCD sum rules, where only three-point correlation function is
needed. It is the first time to investigate the three-body strong decays with
the light-cone QCD sum rules. In our previous works, we have obtained rigorous
quark-hadron duality for the three-point correlation functions, which work
very well. There are other procedures in dealing with the three-point QCD sum
rules exploring the hadronic coupling constants [30, 31, 32], for detailed
discussions about the differences, one can consult Refs.[23, 24].

The article is arranged as follows: we derive the light-cone QCD sum rules for
the $YD^{*}\bar{D}^{*}\pi$ coupling constants in section 2; in section 3, we
present numerical results and discussions; section 4 is reserved for our
conclusion.

##  2 Light-cone QCD sum rules for the $YD^{*}\bar{D}^{*}\pi$ coupling
constants

Firstly, we write down the three-point correlation function
$\Pi_{\mu\alpha\beta}(p)$ in the light-cone QCD sum rules,

$\displaystyle\Pi_{\mu\alpha\beta}(p,q)$ $\displaystyle=$ $\displaystyle
i^{2}\int d^{4}xd^{4}y\,e^{-ip\cdot x}e^{-iq\cdot y}\,\langle
0|T\left\\{J_{\mu}^{Y}(0)J_{\alpha}^{D^{*+}}(x)J^{\bar{D}^{*0}}_{\beta}(y)\right\\}|\pi(r)\rangle\,,$
(2)

where the currents

$\displaystyle J_{\mu}^{Y}(0)$ $\displaystyle=$
$\displaystyle\frac{\varepsilon^{ijk}\varepsilon^{imn}}{2}\Big{[}u^{T}_{j}(0)C\sigma_{\mu\nu}\gamma_{5}c_{k}(0)\bar{u}_{m}(0)\gamma_{5}\gamma^{\nu}C\bar{c}^{T}_{n}(0)+u^{T}_{j}(0)C\gamma^{\nu}\gamma_{5}c_{k}(0)\bar{u}_{m}(0)\gamma_{5}\sigma_{\mu\nu}C\bar{c}^{T}_{n}(0)$
$\displaystyle+d^{T}_{j}(x)C\sigma_{\mu\nu}\gamma_{5}c_{k}(0)\bar{d}_{m}(0)\gamma_{5}\gamma^{\nu}C\bar{c}^{T}_{n}(0)+d^{T}_{j}(0)C\gamma^{\nu}\gamma_{5}c_{k}(0)\bar{d}_{m}(0)\gamma_{5}\sigma_{\mu\nu}C\bar{c}^{T}_{n}(0)\Big{]}\,,$
$\displaystyle J_{\alpha}^{D^{*+}}(y)$ $\displaystyle=$
$\displaystyle\bar{d}(x)\gamma_{\alpha}c(x)\,,$ $\displaystyle
J_{\beta}^{\bar{D}^{*0}}(x)$ $\displaystyle=$
$\displaystyle\bar{c}(y)\gamma_{\beta}u(y)\,,$ (3)

interpolate the mesons $Y(4500)$, $\bar{D}^{*}$ and $D^{*}$, respectively
[22], the $|\pi(r)\rangle$ is the external $\pi$ state. The physical process
is shown explicitly in Fig.1.

In the present work, we take the isospin limit, the current $J_{\mu}^{Y}(x)$
in Eq.(2) and the current $J_{-,\mu}^{\widetilde{A}V}(x)$ chosen in Ref.[22]
couple potentially to the vector tetraquark states with the same masses and
pole residues, where

$\displaystyle J_{-,\mu}^{\widetilde{A}V}(x)$ $\displaystyle=$
$\displaystyle\frac{\varepsilon^{ijk}\varepsilon^{imn}}{\sqrt{2}}\Big{[}u^{Tj}(x)C\sigma_{\mu\nu}\gamma_{5}c^{k}(x)\bar{d}^{m}(x)\gamma_{5}\gamma^{\nu}C\bar{c}^{Tn}(x)$
(4)
$\displaystyle+u^{Tj}(x)C\gamma^{\nu}\gamma_{5}c^{k}(x)\bar{d}^{m}(x)\gamma_{5}\sigma_{\mu\nu}C\bar{c}^{Tn}(x)\Big{]}\,.$

Figure 1:  The decay $Y(4500)\to\bar{D}^{*}D^{*}\pi^{+}$.

At the hadron side, we insert a complete set of intermediate hadronic states
having non-vanishing couplings with the interpolating currents into the three-
point correlation function, and isolate the ground state contributions
clearly,

$\displaystyle\Pi_{\mu\alpha\beta}(p,q)$ $\displaystyle=$
$\displaystyle\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}\frac{-iG_{\pi}r_{\tau}+iG_{Y}p^{\prime}_{\tau}}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}\varepsilon^{\rho\sigma\lambda\tau}\left(-g_{\mu\rho}+\frac{p^{\prime}_{\mu}p^{\prime}_{\rho}}{p^{\prime
2}}\right)$ (5)
$\displaystyle\left(-g_{\alpha\sigma}+\frac{p_{\alpha}p_{\sigma}}{p^{2}}\right)\left(-g_{\lambda\beta}+\frac{q_{\lambda}q_{\beta}}{q^{2}}\right)+\cdots\,,$

where $p^{\prime}=p+q+r$, the decay constants $\lambda_{Y}$, $f_{D^{*}}$,
$f_{\bar{D}^{*}}$ and hadronic coupling constants $G_{\pi}$, $G_{Y}$ are
defined by,

$\displaystyle\langle 0|J_{\mu}^{Y}(0)|Y_{c}(p^{\prime})\rangle$
$\displaystyle=$ $\displaystyle\lambda_{Y}\varepsilon_{\mu}\,\,,$
$\displaystyle\langle 0|J_{\alpha}^{D^{*}{\dagger}}(0)|\bar{D}^{*}(p)\rangle$
$\displaystyle=$ $\displaystyle
f_{\bar{D}^{*}}M_{\bar{D}^{*}}\xi_{\alpha}\,\,,$ $\displaystyle\langle
0|J_{\beta}^{\bar{D}^{*{\dagger}}}(0)|D^{*}(q)\rangle$ $\displaystyle=$
$\displaystyle f_{D^{*}}M_{D^{*}}\zeta_{\beta}\,\,,$ (6) $\displaystyle\langle
Y_{c}(p^{\prime})|\bar{D}^{*}(p)D^{*}(q)\pi(r)\rangle$ $\displaystyle=$
$\displaystyle
G_{\pi}\varepsilon^{\rho\sigma\lambda\tau}\varepsilon^{*}_{\rho}\xi_{\sigma}\zeta_{\lambda}r_{\tau}-G_{Y}\varepsilon^{\rho\sigma\lambda\tau}\varepsilon^{*}_{\rho}\xi_{\sigma}\zeta_{\lambda}p^{\prime}_{\tau}\,\,,$
(7)

the $\varepsilon_{\mu}$, $\xi_{\alpha}$ and $\zeta_{\beta}$ are polarization
vectors of the $Y(4500)$, $\bar{D}^{*}$ and $D^{*}$, respectively. In the
isospin limit, $m_{u}=m_{d}$, $f_{D^{*}}=f_{\bar{D}^{*}}$ and
$M_{D^{*}}=M_{\bar{D}^{*}}$.

We multiply Eq.(5) with the tensor
$\varepsilon_{\theta\omega}{}^{\alpha\beta}$ and obtain

$\displaystyle\widetilde{\Pi}_{\mu\theta\omega}(p,q)$ $\displaystyle=$
$\displaystyle\varepsilon_{\theta\omega}{}^{\alpha\beta}\,\Pi_{\mu\alpha\beta}(p,q)$
(8) $\displaystyle=$
$\displaystyle\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}\frac{iG_{\pi}\left(g_{\mu\omega}r_{\theta}-g_{\mu\theta}r_{\omega}\right)-iG_{Y}\left(g_{\mu\omega}p^{\prime}_{\theta}-g_{\mu\theta}p^{\prime}_{\omega}\right)}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\cdots\,.$

Again, we take the isospin limit, then
$\widetilde{\Pi}_{\mu\theta\omega}(p,q)=\widetilde{\Pi}_{\mu\theta\omega}(q,p)$,
such a relation can simplify the calculations at the QCD side greatly, and we
write down the relevant components,

$\displaystyle\widetilde{\Pi}_{\mu\theta\omega}(p,q)$ $\displaystyle=$
$\displaystyle\left[i\Pi_{\pi}(p^{\prime 2},p^{2},q^{2})-i\Pi_{Y}(p^{\prime
2},p^{2},q^{2})\right]\left(g_{\mu\omega}r_{\theta}-g_{\mu\theta}r_{\omega}\right)$
(9) $\displaystyle+i\Pi_{Y}(p^{\prime
2},p^{2},q^{2})\left(g_{\mu\omega}q_{\theta}-g_{\mu\theta}q_{\omega}\right)+\cdots\,,$

where

$\displaystyle\Pi_{\pi}(p^{\prime 2},p^{2},q^{2})$ $\displaystyle=$
$\displaystyle\frac{\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}G_{\pi}}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\cdots\,,$
$\displaystyle\Pi_{Y}(p^{\prime 2},p^{2},q^{2})$ $\displaystyle=$
$\displaystyle\frac{\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}G_{Y}}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\cdots\,.$ (10)

Then we choose the tensor structures
$g_{\mu\omega}r_{\theta}-g_{\mu\theta}r_{\omega}$ and
$g_{\mu\omega}q_{\theta}-g_{\mu\theta}q_{\omega}$ to study the hadronic
coupling constants $G_{\pi}$ and $G_{Y}$, respectively. And we obtain the
hadronic spectral densities $\rho_{H}(s^{\prime},s,u)$ through triple
dispersion relation,

$\displaystyle\Pi_{H}(p^{\prime 2},p^{2},q^{2})$ $\displaystyle=$
$\displaystyle\int_{\Delta_{s}^{\prime
2}}^{\infty}ds^{\prime}\int_{\Delta_{s}^{2}}^{\infty}ds\int_{\Delta_{u}^{2}}^{\infty}du\frac{\rho_{H}(s^{\prime},s,u)}{(s^{\prime}-p^{\prime
2})(s-p^{2})(u-q^{2})}\,,$ (11)

where the $\Delta_{s}^{\prime 2}$, $\Delta_{s}^{2}$ and $\Delta_{u}^{2}$ are
the thresholds, and we add the subscript $H$ to represent the hadron side.

We carry out the operator product expansion up to the vacuum condensates of
dimension 5 and neglect the tiny gluon condensate contributions [23, 24],

$\displaystyle\Pi_{\pi}(p^{2},q^{\prime 2},q^{2})$ $\displaystyle=$
$\displaystyle
f_{\pi}m_{c}\int_{0}^{1}du\varphi_{\pi}(u)\left[\int_{0}^{1}dxx\bar{x}\frac{\Gamma(\epsilon-1)}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon-1}}-\frac{2m_{c}\langle\bar{q}q\rangle}{3(p^{2}-m_{c}^{2})}\right.$
(12) $\displaystyle\left.+\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{3(p^{2}-m_{c}^{2})^{3}}\right]\frac{1}{(q+ur)^{2}-m_{c}^{2}}$
$\displaystyle+\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\int_{0}^{1}du\varphi_{5}(u)\bar{u}\left[\int_{0}^{1}dxx\bar{x}\frac{\Gamma(\epsilon-1)}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon-1}}-\frac{2m_{c}\langle\bar{q}q\rangle}{3(p^{2}-m_{c}^{2})}\right.$
$\displaystyle\left.+\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{3(p^{2}-m_{c}^{2})^{3}}\right]\frac{1}{(q+ur)^{2}-m_{c}^{2}}$
$\displaystyle-\frac{f_{\pi}m_{c}^{2}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36}\int_{0}^{1}du\varphi_{\pi}(u)\frac{1}{(p^{2}-m_{c}^{2})((q+ur)^{2}-m_{c}^{2})^{2}}$
$\displaystyle+\frac{f_{\pi}m_{\pi}^{2}m_{c}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36(m_{u}+m_{d})}\int_{0}^{1}du\varphi_{5}(u)\bar{u}\frac{1}{(p^{2}-m_{c}^{2})((q+ur)^{2}-m_{c}^{2})^{2}}\,,$

$\displaystyle\Pi_{Y}(p^{2},q^{\prime 2},q^{2})$ $\displaystyle=$
$\displaystyle\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\int_{0}^{1}du\varphi_{5}(u)\left[\int_{0}^{1}dxx\bar{x}\frac{\Gamma(\epsilon-1)}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon-1}}-\frac{2m_{c}\langle\bar{q}q\rangle}{3(p^{2}-m_{c}^{2})}\right.$
(13) $\displaystyle\left.+\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{3(p^{2}-m_{c}^{2})^{3}}\right]\frac{1}{(q+ur)^{2}-m_{c}^{2}}$
$\displaystyle+\frac{f_{\pi}m_{\pi}^{2}m_{c}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36(m_{u}+m_{d})}\int_{0}^{1}du\varphi_{5}(u)\frac{1}{(p^{2}-m_{c}^{2})((q+ur)^{2}-m_{c}^{2})^{2}}$
$\displaystyle+f_{3\pi}m_{\pi}^{2}\int_{0}^{1}dx\bar{x}\left[\frac{3\Gamma(\epsilon)}{8\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon}}-\frac{p^{2}}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})}\right]\frac{1}{q^{2}-m_{c}^{2}}$
$\displaystyle-
f_{3\pi}m_{\pi}^{2}\int_{0}^{1}dxx\bar{x}\left[\frac{\Gamma(\epsilon-1)}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon-1}}+\frac{p^{2}\Gamma(\epsilon)}{2\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon}}\right]\frac{1}{(q^{2}-m_{c}^{2})^{2}}$
$\displaystyle-
f_{3\pi}m_{\pi}^{2}\int_{0}^{1}dxx\left[\frac{3\Gamma(\epsilon)}{8\pi^{2}(p^{2}-\tilde{m}_{c}^{2})^{\epsilon}}+\frac{p^{2}}{4\pi^{2}(p^{2}-\tilde{m}_{c}^{2})}\right]\frac{1}{q^{2}-m_{c}^{2}}\,,$

where $q^{\prime}=q+r$, $\bar{u}=1-u$, $\bar{x}=1-x$,
$\tilde{m}_{c}^{2}=\frac{m_{c}^{2}}{x}$,
$(q-ur)^{2}-m_{c}^{2}=(1-u)q^{2}+u(q+r)^{2}-u\bar{u}m_{\pi}^{2}-m_{c}^{2}$.
And we have used the definitions for the $\pi$ light-cone distribution
functions [33],

$\displaystyle\langle 0|\bar{d}(0)\gamma_{\mu}\gamma_{5}u(x)|\pi(r)\rangle$
$\displaystyle=$ $\displaystyle if_{\pi}r_{\mu}\int_{0}^{1}due^{-iur\cdot
x}\varphi_{\pi}(u)+\cdots\,,$ $\displaystyle\langle
0|\bar{d}(0)\sigma_{\mu\nu}\gamma_{5}u(x)|\pi(r)\rangle$ $\displaystyle=$
$\displaystyle\frac{i}{6}\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\left(r_{\mu}x_{\nu}-r_{\nu}x_{\mu}\right)\int_{0}^{1}due^{-iur\cdot
x}\varphi_{\sigma}(u)\,,$ $\displaystyle\langle
0|\bar{d}(0)i\gamma_{5}u(x)|\pi(r)\rangle$ $\displaystyle=$
$\displaystyle\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\int_{0}^{1}due^{-iur\cdot
x}\varphi_{5}(u)\,,$ (14)

and the approximation,

$\displaystyle\langle
0|\bar{d}(x_{1})\sigma_{\mu\nu}\gamma_{5}g_{s}G_{\alpha\beta}(x_{2})u(x_{3})|\pi(r)\rangle$
$\displaystyle=$ $\displaystyle
if_{3\pi}\left(r_{\mu}r_{\alpha}g_{\nu\beta}+r_{\nu}r_{\beta}g_{\mu\alpha}-r_{\nu}r_{\alpha}g_{\mu\beta}-r_{\mu}r_{\beta}g_{\nu\alpha}\right)\,,$
(15)

for the twist-3 quark-gluon light-cone distribution functions. Such terms
proportional to $m_{\pi}^{2}$ and their contributions are greatly suppressed
[34, 35], the approximation in Eq.(15) works well. However, the terms
proportional to $\frac{m_{\pi}^{2}}{m_{u}+m_{d}}$ are Chiral enhanced due to
the Gell-Mann-Oakes-Renner relation
$\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}=-\frac{2\langle\bar{q}q\rangle}{f_{\pi}}$,
and we take account of those contributions fully in Eqs.(12)-(13). In the
following, we list out the light-cone distribution functions explicitly,

$\displaystyle\varphi_{\pi}(u)$ $\displaystyle=$ $\displaystyle
6u\bar{u}\left[1+A_{2}\frac{3}{2}\left(5t^{2}-1\right)+A_{4}\frac{15}{8}\left(21t^{4}-14t^{2}+1\right)\right]\,,$
$\displaystyle\varphi_{5}(u)$ $\displaystyle=$ $\displaystyle
1+B_{2}\frac{1}{2}\left(3t^{2}-1\right)+B_{4}\frac{1}{8}\left(35t^{4}-30t^{2}+3\right)\,,$
$\displaystyle\varphi_{\sigma}(u)$ $\displaystyle=$ $\displaystyle
6u\bar{u}\left[1+C_{2}\frac{3}{2}\left(5t^{2}-1\right)\right]\,,$ (16)

where $t=2u-1$, and the coefficients $A_{2}=0.44$, $A_{4}=0.25$, $B_{2}=0.43$,
$B_{4}=0.10$, $C_{2}=0.09$, and the decay constant
$f_{3\pi}=0.0035\,\rm{GeV}^{2}$ at the energy scale $\mu=1\,\rm{GeV}$ [33,
36]. In the present work, we neglect the twist-4 light-cone distribution
functions due to their small contributions. In Fig.2, we draw the lowest order
Feynman diagrams as example to illustrate the operator product expansion.

In the soft limit $r_{\mu}\to 0$, $(q+r)^{2}=q^{2}$, we can set
$\Pi_{\pi/Y}(p^{2},q^{\prime 2},q^{2})=\Pi_{\pi/Y}(p^{2},q^{2})$, then we
obtain the QCD spectral densities $\rho_{QCD}(s,u)$ through double dispersion
relation,

$\displaystyle\Pi^{QCD}_{\pi/Y}(p^{2},q^{2})$ $\displaystyle=$
$\displaystyle\int_{\Delta_{s}^{2}}^{\infty}ds\int_{\Delta_{u}^{2}}^{\infty}du\frac{\rho_{QCD}(s,u)}{(s-p^{2})(u-q^{2})}\,,$
(17)

again the $\Delta_{s}^{2}$ and $\Delta_{u}^{2}$ are the thresholds, we add the
superscript or subscript $QCD$ to stand for the QCD side.

Figure 2:  The lowest order Feynman diagrams, where the dashed (solid) lines
denote the heavy (light) quark lines, the ovals denote the external $\pi^{+}$
meson.

We match the hadron side with the QCD side bellow the continuum thresholds
$s_{0}$ and $u_{0}$ to acquire rigorous quark-hadron duality [23, 24],

$\displaystyle\int_{\Delta_{s}^{2}}^{s_{0}}ds\int_{\Delta_{u}^{2}}^{u_{0}}du\frac{\rho_{QCD}(s,u)}{(s-p^{2})(u-q^{2})}$
$\displaystyle=$
$\displaystyle\int_{\Delta_{s}^{2}}^{s_{0}}ds\int_{\Delta_{u}^{2}}^{u_{0}}du\left[\int_{\Delta_{s}^{\prime
2}}^{\infty}ds^{\prime}\frac{\rho_{H}(s^{\prime},s,u)}{(s^{\prime}-p^{\prime
2})(s-p^{2})(u-q^{2})}\right]\,,$ (18)

and we carry out the integral over $ds^{\prime}$ firstly, then

$\displaystyle\Pi_{H}(p^{\prime 2},p^{2},q^{2})$ $\displaystyle=$
$\displaystyle\frac{\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}G_{\pi/Y}}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\int_{s^{\prime}_{0}}^{\infty}ds^{\prime}\frac{\tilde{\rho}_{H}(s^{\prime},M_{\bar{D}^{*}}^{2},M_{D^{*}}^{2})}{(s^{\prime}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}$ (19)
$\displaystyle+\cdots\,,$ $\displaystyle=$
$\displaystyle\frac{\lambda_{Y}f_{D^{*}}^{2}M_{D^{*}}^{2}G_{\pi/Y}}{(M_{Y}^{2}-p^{\prime
2})(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\frac{C_{\pi/Y}}{(M_{\bar{D}^{*}}^{2}-p^{2})(M_{D^{*}}^{2}-q^{2})}+\cdots\,,$

where
$\rho_{H}(s^{\prime},s,u)=\tilde{\rho}_{H}(s^{\prime},s,u)\delta(s-M_{\bar{D}^{*}}^{2})\delta(u-M_{D^{*}}^{2})$,
and we introduce the parameters $C_{\pi/Y}$ to parameterize the contributions
concerning the higher resonances and continuum states in the $s^{\prime}$
channel,

$\displaystyle C_{\pi/Y}$ $\displaystyle=$
$\displaystyle\int_{s^{\prime}_{0}}^{\infty}ds^{\prime}\frac{\tilde{\rho}_{H}(s^{\prime},M_{\bar{D}^{*}}^{2},M_{D^{*}}^{2})}{s^{\prime}-p^{\prime
2}}\,.$ (20)

As the strong interactions among the ground states $\pi$, $D^{*}$,
$\bar{D}^{*}$ and excited $Y^{\prime}$ states are complex, and we have no
knowledge about the corresponding four-hadron contact vertex. In practical
calculations, we can take the unknown functions $C_{\pi/Y}$ as free parameters
and adjust the values to acquire flat platforms for the hadronic coupling
constants $G_{\pi/Y}$ with variations of the Borel parameters. Such a method
works well in the case of three-hadron contact vertexes [23, 24, 25, 26, 27,
28, 29], and we expect it also works in the present work.

In Eq.(5) and Eq.(8), there exist three poles in the limit $p^{\prime 2}\to
M_{Y}^{2}$, $p^{2}\to M_{\bar{D}^{*}}^{2}$ and $q^{2}\to M_{D^{*}}^{2}$.
According to the relation $M_{Y}\approx M_{\bar{D}^{*}}+M_{D^{*}}$, we can set
$p^{\prime 2}=4q^{2}$ in the correlation functions $\Pi_{H}(p^{\prime
2},p^{2},q^{2})$, and perform double Borel transform in regard to the
variables $P^{2}=-p^{2}$ and $Q^{2}=-q^{2}$ respectively, then we set the
Borel parameters $T_{1}^{2}=T_{2}^{2}=T^{2}$ to acquire two QCD sum rules,

$\displaystyle\frac{\lambda_{YD^{*}D^{*}}G_{\pi}}{4\left(\widetilde{M}_{Y}^{2}-M_{D^{*}}^{2}\right)}\left[\exp\left(-\frac{M_{D^{*}}^{2}}{T^{2}}\right)-\exp\left(-\frac{\widetilde{M}_{Y}^{2}}{T^{2}}\right)\right]\exp\left(-\frac{M_{\bar{D}^{*}}^{2}}{T^{2}}\right)+C_{\pi}\exp\left(-\frac{M_{D^{*}}^{2}+M_{\bar{D}^{*}}^{2}}{T^{2}}\right)$
$\displaystyle=f_{\pi}m_{c}\int_{m_{c}^{2}}^{s_{0}}ds\int_{0}^{1}du\varphi_{\pi}(u)\left[\frac{1}{2\pi^{2}}\int_{x_{i}}^{1}dxx\bar{x}(s-\tilde{m}_{c}^{2})-\left(\frac{2m_{c}\langle\bar{q}q\rangle}{3}-\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{6T^{4}}\right)\delta(s-m_{c}^{2})\right]$
$\displaystyle\exp\left(-\frac{s+m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)$
$\displaystyle+\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\int_{m_{c}^{2}}^{s_{0}}ds\int_{0}^{1}du\varphi_{5}(u)\bar{u}\left[\frac{1}{2\pi^{2}}\int_{x_{i}}^{1}dxx\bar{x}(s-\tilde{m}_{c}^{2})-\left(\frac{2m_{c}\langle\bar{q}q\rangle}{3}-\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{6T^{4}}\right)\delta(s-m_{c}^{2})\right]$
$\displaystyle\exp\left(-\frac{s+m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)+\frac{f_{\pi}m_{c}^{2}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36T^{2}}\int_{0}^{1}du\varphi_{\pi}(u)\exp\left(-\frac{2m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)$
$\displaystyle-\frac{f_{\pi}m_{\pi}^{2}m_{c}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36(m_{u}+m_{d})T^{2}}\int_{0}^{1}du\varphi_{5}(u)\bar{u}\exp\left(-\frac{2m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)\,,$
(21)

$\displaystyle\frac{\lambda_{YD^{*}D^{*}}G_{Y}}{4\left(\widetilde{M}_{Y}^{2}-M_{D^{*}}^{2}\right)}\left[\exp\left(-\frac{M_{D^{*}}^{2}}{T^{2}}\right)-\exp\left(-\frac{\widetilde{M}_{Y}^{2}}{T^{2}}\right)\right]\exp\left(-\frac{M_{\bar{D}^{*}}^{2}}{T^{2}}\right)+C_{Y}\exp\left(-\frac{M_{D^{*}}^{2}+M_{\bar{D}^{*}}^{2}}{T^{2}}\right)$
$\displaystyle=\frac{f_{\pi}m_{\pi}^{2}}{m_{u}+m_{d}}\int_{m_{c}^{2}}^{s_{0}}ds\int_{0}^{1}du\varphi_{5}(u)\left[\frac{1}{2\pi^{2}}\int_{x_{i}}^{1}dxx\bar{x}(s-\tilde{m}_{c}^{2})-\left(\frac{2m_{c}\langle\bar{q}q\rangle}{3}-\frac{m_{c}^{3}\langle\bar{q}g_{s}\sigma
Gq\rangle}{6T^{4}}\right)\delta(s-m_{c}^{2})\right]$
$\displaystyle\exp\left(-\frac{s+m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)-\frac{f_{\pi}m_{\pi}^{2}m_{c}\langle\bar{q}g_{s}\sigma
Gq\rangle}{36(m_{u}+m_{d})T^{2}}\int_{0}^{1}du\varphi_{5}(u)\exp\left(-\frac{2m_{c}^{2}+u\bar{u}m_{\pi}^{2}}{T^{2}}\right)$
$\displaystyle-\frac{f_{3\pi}m_{\pi}^{2}}{2\pi^{2}}\int_{m_{c}^{2}}^{s_{0}}ds\int_{x_{i}}^{1}dx\bar{x}\left[\frac{3}{4}+s\,\delta(s-\tilde{m}_{c}^{2})\right]\exp\left(-\frac{s+m_{c}^{2}}{T^{2}}\right)$
$\displaystyle-\frac{f_{3\pi}m_{\pi}^{2}}{2\pi^{2}T^{2}}\int_{m_{c}^{2}}^{s_{0}}ds\int_{x_{i}}^{1}dxx\bar{x}\,\tilde{m}_{c}^{2}\exp\left(-\frac{s+m_{c}^{2}}{T^{2}}\right)$
$\displaystyle+\frac{f_{3\pi}m_{\pi}^{2}}{4\pi^{2}}\int_{m_{c}^{2}}^{s_{0}}ds\int_{x_{i}}^{1}dxx\left[\frac{3}{2}-s\,\delta(s-\tilde{m}_{c}^{2})\right]\exp\left(-\frac{s+m_{c}^{2}}{T^{2}}\right)\,,$
(22)

where $\lambda_{YD^{*}D^{*}}=\lambda_{Y}f_{D^{*}}^{2}M^{2}_{D^{*}}$,
$\widetilde{M}_{Y}^{2}=\frac{M_{Y}^{2}}{4}$ and $x_{i}=\frac{m_{c}^{2}}{s}$.
In numerical calculations, we take the $C_{\pi}$ and $C_{Y}$ as free
parameters, and search for the best values to acquire stable QCD sum rules.

##  3 Numerical results and discussions

We take the standard values of the vacuum condensates,
$\langle\bar{q}q\rangle=-(0.24\pm 0.01\,\rm{GeV})^{3}$,
$\langle\bar{q}g_{s}\sigma Gq\rangle=m_{0}^{2}\langle\bar{q}q\rangle$,
$m_{0}^{2}=(0.8\pm 0.1)\,\rm{GeV}^{2}$ at the energy scale $\mu=1\,\rm{GeV}$
[37, 38, 39], and take the $\overline{MS}$ mass $m_{c}(m_{c})=(1.275\pm
0.025)\,\rm{GeV}$ from the Particle Data Group [1]. We set $m_{u}=m_{d}=0$ and
take account of the energy-scale dependence of the input parameters,

$\displaystyle\langle\bar{q}q\rangle(\mu)$ $\displaystyle=$
$\displaystyle\langle\bar{q}q\rangle({\rm 1GeV})\left[\frac{\alpha_{s}({\rm
1GeV})}{\alpha_{s}(\mu)}\right]^{\frac{12}{33-2n_{f}}}\,,$
$\displaystyle\langle\bar{q}g_{s}\sigma Gq\rangle(\mu)$ $\displaystyle=$
$\displaystyle\langle\bar{q}g_{s}\sigma Gq\rangle({\rm
1GeV})\left[\frac{\alpha_{s}({\rm
1GeV})}{\alpha_{s}(\mu)}\right]^{\frac{2}{33-2n_{f}}}\,,$ $\displaystyle
m_{c}(\mu)$ $\displaystyle=$ $\displaystyle
m_{c}(m_{c})\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}(m_{c})}\right]^{\frac{12}{33-2n_{f}}}\,,$
$\displaystyle\alpha_{s}(\mu)$ $\displaystyle=$
$\displaystyle\frac{1}{b_{0}t}\left[1-\frac{b_{1}}{b_{0}^{2}}\frac{\log
t}{t}+\frac{b_{1}^{2}(\log^{2}{t}-\log{t}-1)+b_{0}b_{2}}{b_{0}^{4}t^{2}}\right]\,,$
(23)

where $t=\log\frac{\mu^{2}}{\Lambda_{QCD}^{2}}$,
$b_{0}=\frac{33-2n_{f}}{12\pi}$, $b_{1}=\frac{153-19n_{f}}{24\pi^{2}}$,
$b_{2}=\frac{2857-\frac{5033}{9}n_{f}+\frac{325}{27}n_{f}^{2}}{128\pi^{3}}$,
$\Lambda_{QCD}=210\,\rm{MeV}$, $292\,\rm{MeV}$ and $332\,\rm{MeV}$ for the
flavors $n_{f}=5$, $4$ and $3$, respectively [1, 40], and we choose $n_{f}=4$.

At the hadron side, we take the parameters as $m_{\pi}=0.13957\,\rm{GeV}$,
$f_{\pi}=0.130\,\rm{GeV}$ [39], $M_{D^{*}}=2.01\,\rm{GeV}$,
$f_{D^{*}}=263\,\rm{MeV}$, $s^{0}_{D^{*}}=6.4\,\rm{GeV}^{2}$ [41],
$M_{Y}=4.48\,\rm{GeV}$, $\lambda_{Y}=9.47\times 10^{-2}\,\rm{GeV}^{5}$ [22],
and $f_{\pi}m^{2}_{\pi}/(m_{u}+m_{d})=-2\langle\bar{q}q\rangle/f_{\pi}$ from
the Gell-Mann-Oakes-Renner relation.

In calculations, we fit the free parameters to be
$C_{\pi}=0.00101(T^{2}-3.6\,\rm{GeV}^{2})\,\rm{GeV}^{4}$ and
$C_{Y}=0.00089(T^{2}-3.2\,\rm{GeV}^{2})\,\rm{GeV}^{4}$ to acquire uniform flat
Borel platforms $T^{2}_{max}-T^{2}_{min}=1\,\rm{GeV}^{2}$ (just like in our
previous works [23, 24, 25, 26, 27, 28, 29]), where the max and min represent
the maximum and minimum values, respectively. The Borel windows are
$T^{2}_{\pi}=(4.6-5.6)\,\rm{GeV}^{2}$ and $T^{2}_{Y}=(4.4-5.4)\,\rm{GeV}^{2}$,
where the subscripts $\pi$ and $Y$ represent the corresponding channels, the
uncertainties $\delta G_{\pi/Y}$ come from the Borel parameters $T^{2}$ are
less than $0.01\,(\rm GeV^{-1})$. In Fig.3, we plot the hadronic coupling
constants $G_{\pi}$ and $G_{Y}$ with variations of the Borel parameters. In
the Borel windows, there appear very flat platforms indeed, it is reasonable
and reliable to extract the $G_{\pi}$ and $G_{Y}$.

If we take the symbol $\xi$ to stand for the input parameters, then the
uncertainties $\bar{\xi}\to\bar{\xi}+\delta\xi$ result in the uncertainties
$\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}\to\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}+\delta\,\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}$,
$\bar{C}_{\pi/Y}\to\bar{C}_{\pi/Y}+\delta C_{\pi/Y}$,

$\displaystyle\delta\,\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}$
$\displaystyle=$
$\displaystyle\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}\left(\frac{\delta
f_{D^{*}}}{\bar{f}_{D^{*}}}+\frac{\delta
f_{\bar{D}^{*}}}{\bar{f}_{\bar{D}^{*}}}+\frac{\delta\lambda_{Y}}{\bar{\lambda}_{Y}}+\frac{\delta
G_{\pi/Y}}{\bar{G}_{\pi/Y}}\right)\,,$ (24)

where the short overline $\bar{}$ on all the input parameters represents the
central values. In calculation, we observe that the uncertainties $\delta
C_{\pi/Y}$ are very small, and set $\delta C_{\pi/Y}=0$ and $\frac{\delta
f_{D^{*}}}{\bar{f}_{D^{*}}}=\frac{\delta
f_{\bar{D}^{*}}}{\bar{f}_{\bar{D}^{*}}}=\frac{\delta\lambda_{Y}}{\bar{\lambda}_{Y}}=\frac{\delta
G_{\pi/Y}}{\bar{G}_{\pi/Y}}$ approximately. Now we obtain the hadronic
coupling constants routinely,

$\displaystyle G_{\pi}$ $\displaystyle=$ $\displaystyle 15.9\pm
0.5\,\rm{GeV}^{-1}\,,$ $\displaystyle G_{Y}$ $\displaystyle=$ $\displaystyle
10.4\pm 0.6\,\rm{GeV}^{-1}\,,$ (25)

by setting

$\displaystyle\delta\,\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}$
$\displaystyle=$
$\displaystyle\bar{\lambda}_{Y}\bar{f}_{D^{*}}\bar{f}_{\bar{D}^{*}}\bar{G}_{\pi/Y}\frac{4\delta
G_{\pi/Y}}{\bar{G}_{\pi/Y}}\,.$ (26)

Figure 3:  The hadronic coupling constants with variations of the Borel
parameters $T^{2}$, where the (I) and (II) denote the $G_{\pi}$ and $G_{Y}$,
respectively, the regions between the two vertical lines are the Borel
windows.

Then it is direct to obtain the partial decay width by taking the hadron
masses $M_{D^{*-}}=2.01026\,\rm{GeV}$, $M_{D^{*0}}=2.00685\,\rm{GeV}$ and
$m_{\pi}=0.13957\,\rm{GeV}$ from the Particle Data Group [1] and
$M_{Y(4500)}=4.4691\,\rm{GeV}$ from the BESIII collaboration [9],

$\displaystyle\Gamma\left(Y(4500)\to D^{*}\bar{D}^{*}\pi^{+}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{24\pi M_{Y}}\int
dk^{2}(2\pi)^{4}\delta^{4}(p^{\prime}-k-p)\frac{d^{3}\vec{k}}{(2\pi)^{3}2k_{0}}\frac{d^{3}\vec{p}}{(2\pi)^{3}2p_{0}}$
(27)
$\displaystyle(2\pi)^{4}\delta^{4}(k-q-r)\frac{d^{3}\vec{q}}{(2\pi)^{3}2q_{0}}\frac{d^{3}\vec{r}}{(2\pi)^{3}2r_{0}}\Sigma|T|^{2}$
$\displaystyle=$ $\displaystyle 6.43^{+0.80}_{-0.76}\,\rm{MeV}\,,$

where $T=\langle Y_{c}(p^{\prime})|\bar{D}^{*}(p)D^{*}(q)\pi(r)\rangle$
defined in Eq.(7).

The partial decay width $\Gamma\left(Y(4500)\to
D^{*}\bar{D}^{*}\pi^{+}\right)=6.43^{+0.80}_{-0.76}\,\rm{MeV}$ is much smaller
than the total width $\Gamma=246.3\pm 36.7\pm 9.4\,\rm{MeV}$ from the BESIII
collaboration [9], which is consistent with our naive expectation that the
main decay channels of the vector tetraquark states are two-body strong decays
$Y\to D\bar{D}$, $D^{*}\bar{D}^{*}$, $D\bar{D}^{*}$, $D^{*}\bar{D}$,
$J/\psi\pi$, $\eta_{c}\rho$. The observations of the $Y(4500)$ in the channels
$D\bar{D}$, $D^{*}\bar{D}^{*}$, $D\bar{D}^{*}$, $D^{*}\bar{D}$, $J/\psi\pi$,
$\eta_{c}\rho$ would shed light on the nature of the $Y(4500)$, and we would
explore those two-body strong decays in our next work in a comprehensive way.

We choose the process $Y(4500)\to D^{*-}D^{*0}\pi^{+}$ to explore whether or
not the four-hadron coupling constants can be calculated directly using the
(light-cone) QCD sum rules, as this process is not expected to be the dominant
decay channel, which only servers as a powerful constraint to examine the
calculations, i.e. the partial decay width should be small enough to satisfy
the BESIII experimental data. We should admit that it would be better to find
a tetraquark candidate, whose dominant decay mode is the three-body strong
decay, to examine the present approach (or procedure), however, at the present
time, we cannot find such a tetraquark candidate. In short, the present work
supports assigning the $Y(4500)$ to be the
$[uc]_{\tilde{A}}[\overline{uc}]_{V}+[uc]_{V}[\overline{uc}]_{\tilde{A}}+[dc]_{\tilde{A}}[\overline{dc}]_{V}+[dc]_{V}[\overline{dc}]_{\tilde{A}}$
hidden-charm tetraquark state with the quantum numbers $J^{PC}=1^{--}$. It is
the first time to use the light-cone QCD sum rules to study the four-hadron
coupling constants, the approach can be used to explore the $Y\to
J/\psi\pi^{+}\pi^{-}$, $\psi^{\prime}\pi^{+}\pi^{-}$, $J/\psi K^{+}K^{-}$,
$h_{c}\pi^{+}\pi^{-}$, $D^{0}D^{*-}\pi^{+}$, and diagnose the nature of the
$X$, $Y$ and $Z$ states.

##  4 Conclusion

In this work, we tentatively assign the $Y(4500)$ as the
$[uc]_{\tilde{A}}[\overline{uc}]_{V}+[uc]_{V}[\overline{uc}]_{\tilde{A}}+[dc]_{\tilde{A}}[\overline{dc}]_{V}+[dc]_{V}[\overline{dc}]_{\tilde{A}}$
tetraquark state with the quantum numbers $J^{PC}=1^{--}$, and extend our
previous works to study the three-body strong decay $Y(4500)\to
D^{*-}D^{*0}\pi^{+}$ with the light-cone QCD sum rules, the partial width is
consistent with the experimental data from the BESIII collaboration. It is the
first time to use the light-cone QCD sum rules to study the four-hadron
coupling constants, we choose the process $Y(4500)\to D^{*-}D^{*0}\pi^{+}$ to
explore whether or not the (light-cone) QCD sum rules can be used to calculate
the four-hadron coupling constants directly, as the process is not the main
decay channel, which servers as a powerful constraint to testify the approach,
i.e. the partial decay width should be small enough to be match the
experimental data. The approach can be used to investigate the three-body
strong decays $X/Y\to J/\psi\pi^{+}\pi^{-}$, $\psi^{\prime}\pi^{+}\pi^{-}$,
$J/\psi K^{+}K^{-}$, $h_{c}\pi^{+}\pi^{-}$, $D^{0}D^{*-}\pi^{+}$ directly, and
shed light on the nature of the $X$, $Y$ and $Z$ states.

## Acknowledgements

This work is supported by National Natural Science Foundation, Grant Number
12175068.

## References

  * [1] R. L. Workman et al, Prog. Theor. Exp. Phys. 2022 (2022) 083C01.
  * [2] B. Aubert et al, Phys. Rev. Lett. 95 (2005) 142001.
  * [3] M. Ablikim et al, Phys. Rev. Lett. 118 (2017) 092002.
  * [4] M. Ablikim et al, Phys. Rev. Lett. 118 (2017) 092001.
  * [5] X. L. Wang et al, Phys. Rev. Lett. 99 (2007) 142002.
  * [6] X. L. Wang et al, Phys. Rev. D91 (2015) 112007.
  * [7] G. Pakhlova et al, Phys. Rev. Lett. 101 (2008) 172001.
  * [8] M. Ablikim et al, Chin. Phys. C46 (2022) 111002.
  * [9] M. Ablikim et al, Phys. Rev. Lett. 130 (2023) 121901.
  * [10] M. Ablikim et al, Phys. Rev. Lett. 122 (2019) 102002.
  * [11] R. M. Albuquerque and M. Nielsen, Nucl. Phys. A815 (2009) 532009; Erratum-ibid. A857 (2011) 48.
  * [12] W. Chen and S. L. Zhu, Phys. Rev. D83 (2011) 034010.
  * [13] Z. G. Wang, Eur. Phys. J. C78 (2018) 518.
  * [14] Z. G. Wang, Eur. Phys. J. C74 (2014) 2874.
  * [15] Z. G. Wang, Eur. Phys. J. C76 (2016) 387.
  * [16] J. R. Zhang and M. Q. Huang, Phys. Rev. D83 (2011) 036005.
  * [17] J. R. Zhang and M. Q. Huang, JHEP 1011 (2010) 057.
  * [18] Z. G. Wang, Eur. Phys. J. C78 (2018) 933.
  * [19] Z. G. Wang, Eur. Phys. J. C79 (2019) 29.
  * [20] Z. G. Wang, Commun. Theor. Phys. 71 (2019) 1319.
  * [21] H. Sundu, S. S. Agaev and K. Azizi, Phys. Rev. D98 (2018) 054021.
  * [22] Z. G. Wang, Nucl. Phys. B973 (2021) 115592.
  * [23] Z. G. Wang and J. X. Zhang, Eur. Phys. J. C78 (2018) 14.
  * [24] Z. G. Wang, Eur. Phys. J. C79 (2019) 184.
  * [25] Z. G. Wang and Z. Y. Di, Eur. Phys. J. C79 (2019) 72.
  * [26] Z. G. Wang, Acta Phys. Polon. B51 (2020) 435.
  * [27] Z. G. Wang, Int. J. Mod. Phys. A34 (2019) 1950110.
  * [28] Z. G. Wang, Chin. Phys. C46 (2022) 103106.
  * [29] Z. G. Wang, Chin. Phys. C46 (2022) 123106.
  * [30] J. M. Dias, F. S. Navarra, M. Nielsen and C. M. Zanetti, Phys. Rev. D88 (2013) 016004.
  * [31] W. Chen, T. G. Steele, H. X. Chen and S. L. Zhu, Eur. Phys. J. C75 (2015) 358.
  * [32] H. Sundu, S. S. Agaev and K. Azizi, Eur. Phys. J. C79 (2019) 215.
  * [33] P. Ball, JHEP 9901 (1999) 010.
  * [34] Z. G. Wang and S. L. Wan, Phys. Rev. D74 (2006) 014017.
  * [35] Z. G. Wang, J. Phys. G34 (2007) 753.
  * [36] V. M. Braun and I. E. Filyanov, Z. Phys. C48 (1990) 239.
  * [37] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B147 (1979) 385; Nucl. Phys. B147 (1979) 448.
  * [38] L. J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rept. 127 (1985) 1.
  * [39] P. Colangelo and A. Khodjamirian, hep-ph/0010175.
  * [40] S. Narison and R. Tarrach, Phys. Lett. 125 B (1983) 217.
  * [41] Z. G. Wang, Eur. Phys. J. C75 (2015) 427.


# Early Stage Growth of Amorphous Thin Film:
Average Kinetics, Nanoscale Dynamics and Pressure Dependence

Chenyu Wang<EMAIL_ADDRESS>Department of Physics, Boston University, Boston,
MA 02215  Karl F. Ludwig, Jr<EMAIL_ADDRESS>Department of Physics and Division
of Materials Science and Engineering, Boston University, Boston, MA 02215
Christa Wagenbach  Current Address: Thermo Fisher Scientific Inc.  Meliha G.
Rainville  Division of Materials Science and Engineering, Boston University,
Boston, MA 02215  Suresh Narayanan  Hua Zhou  Advanced Photon Source, Argonne
National Lab, Lemont, IL 60439  Jeffrey G. Ulbrandt  Randall L. Headrick
Department of Physics and Materials Science Program, University of Vermont,
Burlington, VT 05405

###### Abstract

We used the Coherent Grazing Incidence Small Angle X-Ray Scattering (Co-
GISAXS) technique to study the average kinetics and nanoscale dynamics during
early-stage a-WSi2 sputter deposition. The kinetic and dynamic properties are
examined as a function of pressure, which is known to be a critical factor in
determining final surface roughness. Surface growth kinetics and dynamics are
characterized by time parameters extracted from the height-height structure
factor and correlation functions. The roughness at a given length scale
reaches a maximum before relaxing down to a steady state. The lateral length
scale dependence and pressure dependence are then compared among measured
kinetics and dynamics time parameters. Surfaces grown at lower pressures are
smoother, leading to longer correlation times. The time to reach a dynamic
steady state and a kinetic steady state show contrasting pressure dependence.
A dynamic steady state is reached earlier than the kinetic steady state at
high pressure. A more random deposition direction and lower kinetic energy at
higher pressures can explain these phenomena, along with the hypothesis that
larger nanoclusters form in vapor before arriving at the surface. A continuum
model is applied to simulate the overall behavior with mixed success.

††preprint: APS/123-QED

##  I Introduction

The high level of materials science development in recent years enables the
“materials by design” concept [1]. In-depth understanding of the structure and
surface growth is a prerequisite to design materials with novel functions to
meet the increasing industrial and societal needs. Non-epitaxial thin film
growth exhibits qualitatively universal morphologies which have long been
encoded in structure zone models based upon empirical relationships between
growth conditions and post-facto observations of film structure [2, 3]. During
a thin film growth process, the final structure arises from the kinetics of
the average surface morphology evolution, which is in turn determined by local
dynamical processes. In this context, the average kinetics means the spatially
averaged evolution over the nanoscale surface structure, and the local
dynamics means the nanoscale fluctuation in temporal evolution [4]. At the
atomic scale these dynamical processes include adatom arrival, diffusion,
incorporation, and desorption [5]. During ongoing growth of the thin film, we
distinguish here between obtaining a kinetic steady state, when there is no
longer evolution of the average structure at a given length scale, and
obtaining the dynamic steady state, when the spectrum of local dynamic
processes reaches a steady state.

When used in a surface sensitive mode, X-ray Photon Correlation Spectroscopy
(XPCS) enables the analysis of dynamical process in the near surface
structure. Coherent X-ray scattering can produce a speckle pattern, which
reflects the temporal configuration in the coherent illuminated volume [6].
Under coherent illumination, scattering from surface height fluctuations forms
a speckle pattern on the detector, and one can infer the dynamical processes
leading to such evolution [7]. In this paper, we use XPCS to examine the
relationship between average surface structure and the underlying dynamics
during early stage sputter deposition of a-WSi2 as a function of growth
pressure. The technique of Grazing Incidence Small Angle X-Ray Scattering
(GISAXS) is applied to improve the surface sensitivity. It is known that
growth pressure plays a critical role during the growth of a-WSi2, with a
transition occurring between relatively smooth growth surfaces at low pressure
and rough growth surfaces at higher pressure [8].

A number of theories of amorphous thin film growth have been developed to
explain the morphology evolution [9, 10, 11, 12, 13, 14, 15, 16, 17, 18].
Nonlinear continuum models are the most concise method to simulate growth
surface evolution. The data obtained in these XPCS experiments provide
considerable information about evolving kinetic and dynamic properties of the
growing film in $k$-space. Therefore a model of surface growth was simulated
and the results were compared with the X-ray scattering. With the
implementation of further analysis, we can then discuss the ability of the
continuum model to reproduce the experimental results.

The rest of the paper is organized as follows. In Section II, we introduce the
main experimental technique we used, Co-GISAXS, and give a detailed
description of the sputter deposition experiments under different pressure
conditions. In Section III, the data taken under three pressure conditions are
analyzed in both kinetics and dynamics aspects. The height-height structure
factor, two-time correlation function and $g_{2}$ autocorrelation function are
calculated to extract different time parameters which characterize the early
stage behavior. In Section IV, we discuss a continuum model of amorphous thin
film growth and the simulation results. The experimental and simulation
results are compared in this section. In Section V comprehensive conclusions
are given.

##  II Experiments

###  II.1 XPCS and GISAXS

Photon Correlation Spectroscopy (PCS) provides information on the dynamic
behavior of condensed matter across a wide range of length- and time-scales.
Since the wavelength of hard X-rays lies in the subnanometer range, XPCS can
provide such information at the nanoscale. XPCS measures spontaneous density
fluctuations of a material. In the kinetic steady state, conventional X-ray
scattering fails to provide such information about a material growth since
incoherent X-ray beam measures only the ensemble-averaged density correlation
function [19].

Fig.1 is a schematic diagram of Co-GISAXS geometry. $\mathbf{k}_{i}$ and
$\mathbf{k}_{f}$ are the wave vectors of incident and scattered beams, making
angles $\alpha_{i}$ and $\alpha_{f}$ relative to the surface. The change of
wave vector during the scattering is

$\mathbf{q}=\mathbf{k}_{f}-\mathbf{k}_{i}$ (1)

and the magnitude of $\mathbf{q}$ is denoted as $q$. In order to obtain the
scattering pattern with surface information, the use of a grazing incidence
angle is required. During the diffuse scattering process, there are X-rays
propagating along the surface. The beam exits at an angle of total external
reflection, forming the band of enhanced scattering in Fig.1, known as the
Yoneda wing [20]. When $\alpha_{i}\leqslant\alpha_{c}$ and
$\alpha_{f}\leqslant\alpha_{c}$, where $\alpha_{c}$ is the critical angle of
total reflection, the scattering depth is at the order of a few nanometers.
Under such circumstances, the scattered beam is surface sensitive [21], so it
is not mixed with much bulk information during the growth. Co-GISAXS is a non-
destructive technique which can be performed as a real-time measurement.

###  II.2 Experimental details

Figure 1: The experimental setup for GISAXS studies during WSi2 sputter
deposition

Coherent X-ray experiments were performed in a GISAXS geometry at the 8-ID-I
beamline of the Advanced Photon Source (APS) (Fig.1) at Argonne National
Laboratory. The scattering in the horizontal direction on the detector is
primarily in the plane of the sample, and for simplicity the wavenumber of
scattering in this direction will be referred to as $q_{\parallel}$. The
energy of the incoming X-rays was $7.38~{}\mathrm{k}\mathrm{eV}$ and the
focused beam was $20~{}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ horizontal and
$4~{}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ vertical. The X-ray flux at the
sample was $\sim 7\times 10^{10}$ photons/ s . Scattering was captured by a
two-dimensional Princeton Instruments LCX-1300 CCD detector, located
$4067~{}\mathrm{m}\mathrm{m}$ from the sample. The detector has 1340 pixels
horizontal and 1300 pixels vertical with each pixel being $20\times
20~{}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$. In the experimental geometry,
it covered a range of $q_{\parallel}$ of $0.23\sim
0.47~{}\mathrm{n}\mathrm{m}^{-1}$. Data collection scans were 1024 frames in
length, with each frame having 1 second of data and followed by a 1 second
readout time of the detector. The X-ray incidence angle, $\alpha_{i}$, was
varied around $0.45\mathrm{\SIUnitSymbolDegree}$, always remaining close to
and above the critical angle $\alpha_{c}=0.4\mathrm{\SIUnitSymbolDegree}$ of
the sample. The relationship between $\alpha_{i}$ and $\alpha_{c}$ could
ensure scattering from the surface instead of the bulk of the film. For the
data analyzed here, the X-ray exit angle, $\alpha_{f}$, was centered around
the Yoneda peak.

The base pressure of the custom UHV chamber was $7\times
10^{-7}~{}\mathrm{Torr}$ and it used $2^{\prime\prime}$ MeiVac MAK DC Sputter
source at a working distance of $10~{}\mathrm{c}\mathrm{m}$. The sputter
deposition source is located at a position normal to the sample surface. The
target was pre-sputtered for at least $30~{}\mathrm{min}$ prior to deposition
with the sample shutter closed to remove contamination from the target.
Deposition was performed at room temperature onto Si substrates with a
$500~{}\mathrm{n}\mathrm{m}$ thermal SiO2 layer. Pressures of 16 mTorr (Sample
1), 10 mTorr (Sample 2) and 8 mTorr (Sample 3) were used with ultra-high
purity argon (99.999%) as the working gas. It is known that there is a
transition around 5 mTorr growth pressure between smooth films below this
value to films which exhibit increasing roughness with increasing pressure
[22]. However, due to the decreased X-ray scattering from smoother films,
growth at 8 mTorr was the lowest for which we could obtain sufficient signal
for XPCS experiments.

For a given deposition power, the deposition rate was approximately
independent of pressure. Sample 2 and 3, grown at $25~{}\mathrm{W}$ of
deposition power, had a growth rate of
$0.18~{}\mathrm{n}\mathrm{m}/\mathrm{s}$. Sample 1 was grown at
$50~{}\mathrm{W}$. Doubling the deposition power increases the deposition rate
by a factor of two and likewise decreases the growth time constants by a
factor of two. Real-time GISAXS measurements show that surface evolution stops
when deposition ends, suggesting that the dominant surface evolution at this
temperature is driven by the deposition itself and that thermal effects at
this temperature are minimal. We therefore do not expect the difference in
deposition rates to affect the fundamental evolution of the surface structure,
but only its rate. To aid the reader in comparability of the samples, the
timescales from Sample 1 will be scaled by a factor of two so that samples are
directly comparable using “effective time”. The “effective time” is used as
the default time scale except in Fig.3.

##  III Data Analysis

###  III.1 Height-height structure factor

In the experimental geometry (Fig.1), the projected direction of the incident
beam on the sample surface is taken to be the $x$-axis, and the $z$-axis
points normally out of the surface. The $y$-axis is taken to be perpendicular
to the $x$ direction in the sample surface. Therefore, the Yoneda wing extends
along the $y$ direction in Fig.1 and the detector is in the $yz$-plane. When
surface features are not too high, the intensity $I(\mathbf{q},t)$ on the
Yoneda wing is approximately proportional to the square modulus of the Fourier
transformation of the surface height [23], i.e., $I(\mathbf{q},t)\propto
S(\mathbf{q},t)$, where the height-height structure factor $S(\mathbf{q},t)$
is defined as

$S(\mathbf{q},t)=h(\mathbf{q},t)h^{*}(\mathbf{q},t)$ (2)

where $h(\mathbf{q},t)$ is the Fourier transformation of $h(x,y,t)$.

(a)

(b)

(c)

Figure 2: Intensity (arb. unit) on the detector of Sample 1 (a) before
deposition; the intensity is too low to be seen on this scale; (b) 78 s after
the start of deposition; and (c) 758 s after the start of deposition (kinetic
steady state). The $x$ and $y$ axes are the pixel numbers. Figure 3: Speckle
averaged intensity vs. time in the early stage deposition. “Effective time” is
not used here.

Fig.2 shows how the scattering pattern for Sample 1 changes with time. It is
seen that the scattering intensity increases to a maximum (Fig.2b) before
reaching a steady state (Fig.2c). In order to demonstrate the kinetics
behavior, the evolution of speckle averaged intensity along the Yoneda wings
as a function of time are plotted in Fig.3. Deposition starts at about 40 s.
Each dot in Fig.3 is calculated as the averaged intensity over a square of
pixels centered at the corresponding $\mathbf{q}$ at the Yoneda wing.
Therefore, each curve in Fig.3 represents one of the $q_{\parallel}$ positions
along the Yoneda wing. In the regime of the detector,

$q_{\parallel}\approx q_{y}=\frac{2\pi}{\lambda}\cos(\alpha_{f})\sin(\psi)$
(3)

because $q_{x}\ll q_{y}$. In each curve, the intensity increases rapidly in
time up to a peak and then descends to a steady-state value. Moreover, the
peak intensity is reached at a later time for lower $q_{\parallel}$ values
than for higher $q_{\parallel}$’s. When a curve with position $q_{\parallel}$
reaches the peak, the roughness on the length scale
$\frac{2\pi}{q_{\parallel}}$ is maximum. Focusing on the kinetic steady-state
behavior, we find that positions with lower $q_{\parallel}$ stabilize later
than the higher $q_{\parallel}$ positions. This behavior is consistent with
mound formation and coarsening, so that the kinetic properties reach a steady
state later at longer length scales. Comparing the results among the three
pressures, we see that the time to reach the kinetic steady state increases
with decreasing pressure, while the final scattering intensity decreases with
decreasing pressure. This is consistent with the known result that the grown
surface is smoother at lower pressures [8].

The time to reach the peak intensity is denoted as $T_{max}$ and the time to
reach the steady-state intensity (kinetic steady state) is called $T_{f}$.
These two time parameters, $T_{max}$ and $T_{f}$, characterize the kinetic
evolution at each wavenumber during the thin film growth. On account of the
early-stage asymmetric peak and the late-stage steady state, the $I(t)$ curves
can be fit with an empirical function consisting of the sum of a lognormal
distribution plus a hyperbolic tangent function

$I(t)=I_{1}(t)+I_{2}(t)$ (4a) where $\displaystyle I_{1}(t)$
$\displaystyle=I_{max}\frac{1}{\sqrt{2\pi}\sigma
t}e^{-\frac{(\ln(t)-\mu)^{2}}{2\sigma^{2}}}$ (4b) $\displaystyle I_{2}(t)$
$\displaystyle=I_{\infty}\tanh(\frac{x}{T_{I}})$ (4c)

Figure 4: An example of speckle averaged intensity $I$ as a function of film
thickness and curve fitting. Thickness is calculated as a product of
deposition rate and deposition time.

All fitting parameters in Eq.(4) are $q_{\parallel}$ dependent. $T_{max}$ is
read directly as the time when the maximum intensity at a given wavenumber is
reached, and $T_{f}$ can be calculated from the curve fitting (Fig.4)

$T_{f}=e^{\mu+2\sigma}$ (5)

On the other hand, the evolution of intensity can be replotted in a way that
the $I$ vs. $q_{\parallel}$ curves evolve with time (Fig.5). At all
$q_{\parallel}$ positions, the intensity increases rapidly with time at the
very beginning. At low $q_{\parallel}$, the intensities rise to their peaks
later than at high $q_{\parallel}$, so that crossovers between curves are
observed. Ultimately the intensities reach their steady state values. The
crossovers of different curves is consistent with coarsening of lateral
structure during growth [24].

Figure 5: Speckle averaged Yoneda peak intensities of Sample 1 as a function
of $q_{\parallel}$ at different time (thicknesses) Figure 6: Speckle averaged
intensity $I$ vs. $q_{\parallel}$ plots and fits in the kinetic steady state
of all three samples

Scans of a broader $q_{\parallel}$ range are taken after the early stage
evolution, and the results are shown in Fig.6. The curves are fit by a power
law plus a Gaussian tail

$I(q_{\parallel})=I_{p}q_{\parallel}^{-m}+I_{g}e^{-\frac{q_{\parallel}^{2}}{2\sigma^{2}}}$
(6)

The power law behavior in low $q_{\parallel}$ conforms to fractal growth [5],
and the Gaussian tail in high $q_{\parallel}$ ranges has been observed in
previous experiments [23]. As discussed in Ref.[23], it is believed that the
Gaussian component contains the bulk structure information from the near
surface region. The lower intensity of the low pressure sample indicates a
smoother surface at all length scales examined.

###  III.2 Two-time correlation function

The kinetic behavior only gives us the information on the average evolution.
In order to go beyond the average structure, the two-time correlation function
(2TCF) is applied to explore the dynamical correlation process. The 2TCF is
calculated as

$C(\mathbf{q},t_{1},t_{2})=\frac{\langle
I(\mathbf{q},t_{1})I(\mathbf{q},t_{2})\rangle_{\mathbf{q}}}{\langle
I(\mathbf{q},t_{1})\rangle_{\mathbf{q}}\langle
I(\mathbf{q},t_{2})\rangle_{\mathbf{q}}}$ (7)

The $\langle I(\mathbf{q},t)\rangle_{\mathbf{q}}$ means the averaged intensity
over a rectangle of pixels centered at $\mathbf{q}$. In practice, however, the
Savitzky-Golay filter is applied to better estimate the normalized scattering
[25]

$C(\mathbf{q},t_{1},t_{2})=\langle
I^{\prime}(\mathbf{q},t_{1})I^{\prime}(\mathbf{q},t_{2})\rangle_{\mathbf{\mathbf{q}}}$
(8a) where
$I^{\prime}(\mathbf{q},t)=\frac{I(\mathbf{q},t)}{\tilde{I}(\mathbf{q},t)}$
(8b)

$\tilde{I}(\mathbf{q},t)$ is composed of pixels smoothed by the 2D Savitzky-
Golay filter.

Figure 7: 2TCF of Sample 2 at
$q_{\parallel}=0.463~{}\mathrm{n}\mathrm{m}^{-1}$ at the Yoneda wing. $t_{1}$
and $t_{2}$ are deposition time. $t_{1,2}<0$ before deposition starts. The
whole figure is symmetric about the main diagonal line $t_{1}=t_{2}$ by
definition. The white arrows demonstrate the direction of slicing along the
off-diagonal lines keeping $T$ as a constant in Eq.(9). Figure 8: Off-diagonal
line slices of 2TCF along the three white arrows in Fig.7 and the KWW fits.
$T=127~{}\mathrm{s}$ lies in the early stage, and the other two curves lie in
the dynamic steady state.

The 2TCF of Sample 2 for $q_{\parallel}=0.463~{}\mathrm{n}\mathrm{m}^{-1}$ at
the Yoneda wing is given in Fig.7. The red band near the diagonal line shows
the high correlation of the surface when $t_{1}$ is close to $t_{2}$. The
‘width’ of the central band tells us how rapidly the surface structure changes
during the growth of layers. The dynamic steady state is attained when the
correlation width becomes stable, implying that the dynamical processes at the
scale of $\frac{2\pi}{q_{\parallel}}$ has reached a dynamical equilibrium. We
need to go further, performing off-diagonal slices through the $t_{1}=t_{2}$
line in order to characterize the dynamical properties. Keeping

$T=\frac{1}{2}(t_{1}+t_{2})$ (9)

as a constant, three typical off-diagonal slices are shown in Fig.8, where
$\Delta t=\left|t_{2}-t_{1}\right|$ is the $x$-axis. The curves are fit by the
Kohlrausch-Williams-Watts (KWW) form [26]

$C_{\mathbf{q},T}(\Delta t)=B+\beta e^{-2(\frac{\Delta
t}{\tau(\mathbf{q};T)})^{n}}$ (10)

where $\beta$ is the contrast, $\tau(\mathbf{q},T)$ is the correlation time,
$n$ determines the compressed ($n\geqslant 1$) or stretched ($n\leqslant 1$)
exponential form. Since we fix $q_{z}$ at the Yoneda wing,
$\tau(\mathbf{q},T)$ is actually $\tau(q_{\parallel},T)$ in this context. In
our experiments, $n$ is always greater than $1$, indicating a compressed
exponential relaxation. $B$ is the baseline, which is controlled close to 1
(varied from 0.9 to 1.1 in the fittings). The compressed exponential result
suggests nonlinear dynamics [27] which will be discussed further in Section
IV.

Figure 9: $\tau$-$T$ curve at
$q_{\parallel}=0.346~{}\mathrm{n}\mathrm{m}^{-1}$ of Sample 2. The
oscillations in the dynamic steady state are due to challenges from
normalizations of the constantly evolving speckle patterns.

Once we get all fitting parameters along the $T$ axis, the function of
correlation time $\tau$ vs. $T$ is obtained (Fig.9). The curves can be fit
empirically using an exponential decay function:

$\tau_{q_{\parallel}}(T)=\tau_{\infty}(1-e^{-\frac{T}{T_{\tau}}})+\tau_{0}e^{-\frac{T}{T_{\tau}}}$
(11)

The $q_{\parallel}$-dependent fitting parameters $T_{\tau}$ and
$\tau_{\infty}$ are used as dynamical time parameters. The time $2T_{\tau}$
approximately characterizes the time from the beginning of the deposition to
reaching the dynamic steady state, and $\tau_{\infty}$ characterizes the
correlation time in the dynamic steady state. Note that $T_{\tau}$ is the time
for the surface to fundamentally change its morphology on a given length
scale. Thus $2T_{\tau}$ is the time it takes for the local surface morphology
evolution to reach a steady state rate of change. $\tau_{\infty}$
characterizes the time scale of that steady state rate of change on a given
length scale.

###  III.3 $g_{2}$ function

The steady state correlation property can also be calculated according to the
autocorrelation function

$g_{2}(\mathbf{q};\Delta t)=\langle\frac{\langle
I(\mathbf{q},t)I(\mathbf{q},t+\Delta t)\rangle_{t}}{\langle
I(\mathbf{q},t)\rangle_{t}^{2}}\rangle_{\mathbf{q}}$ (12)

$g_{2}$ analysis makes up the basis of dynamic light scattering (DLS) [7].
There are two brackets of average. $\langle\cdots\rangle_{\mathbf{q}}$ has the
same meaning as in Eq.(7). $\langle\cdots\rangle_{t}$ takes the average over
time. Symmetrical normalization [28] is applied to avoid violent oscillations
at large lag times, and the multi-tau method is applied for the consideration
of calculation efficiency [29, 30]. The curves are fit with the KWW form
(Eq.(10)) with the baseline constrained to determine $\tau$.

Figure 10: $\tau$-$q_{\parallel}$ curves of three pressures given by $g_{2}$
analysis

The resulting $\tau$ values are shown in Fig.10. The high $q_{\parallel}$ data
of Sample 3 is missing because the intensity on the detector is too weak to
calculate $g_{2}$ functions. $\tau$ decreases with increasing $q_{\parallel}$
and pressure, demonstrating that under lower pressure, the smoother surfaces
remain correlated for a longer period of time during the growth process. The
$\tau$’s obtained in $g_{2}$ functions are used as the best estimated
$\tau_{\infty}$ since the $g_{2}$ functions are based on a large amount of
steady-state data, which reflect the steady-state dynamical behavior more
accurately.

###  III.4 Results and discussions

Time Params | K/D | Characterization | Physical Interpretation
---|---|---|---
$T_{max}(q_{\parallel})$ | Kinetics |  $I$-$t$ curves reach the max. | The mounds have coarsened through the length scale $\frac{2\pi}{q_{\parallel}}$.
$T_{f}(q_{\parallel})$ | Kinetics |  $I$-$t$ curves reach the steady state. | Local morphology reaches steady state at the scale.
$T_{\tau}(q_{\parallel})$ | Dynamics | 2TCFs reach the steady state. | Local dynamics reaches steady state at the scale.
|  |  | Surface structure at the corresponding length scale
$\tau_{\infty}(q_{\parallel})$ | Dynamics | The width of 2TCFs | remains correlated over this time (thickness).
Table 1: Illustration of four time parameters Figure 11: The $q_{\parallel}$
dependence of the four experimentally determined time parameters for the 3
samples

The overall goal of the X-ray data analysis has been to extract two kinetic
parameters, $T_{max}$ and $T_{f}$, and two dynamic parameters, $T_{\tau}$ and
$\tau_{\infty}$. A summary is given in Table 1 and the final results are given
in Fig.11. The $T_{\tau}$ for Sample 3 is absent because the surface is too
smooth to extract signals with meaningful contrast in the 2TCF. Effective
deposition time as mentioned in Section II.2 has been applied to balance the
different deposition rates.

In order to demonstrate that the kinetic energy and the velocity direction of
adatoms are strongly affected by pressure, let us estimate the mean free paths
through which silicon and tungsten molecules travel in argon atmosphere. The
mean free path of a gas molecule is calculated as:

$\mathrm{bar}{\lambda}=\frac{k_{B}T}{\sqrt{2}\sigma p}$ (13)

where $k_{B}$ is the Boltzmann constant, $T$ is the gas temperature, $\sigma$
is the scattering cross section, and $p$ is the gas pressure. The cross
section can be calculated as $\sigma=\pi d^{2}$, where $d$ is the distance at
which there starts to be interatomic forces between two atoms. We now estimate
the scattering cross section of a tungsten or silicon atom scattered by an
argon atom. The radius of one silicon atom is
$r_{\mathrm{Si}}=0.132~{}\mathrm{n}\mathrm{m}$, the radius of one tungsten
atom is $r_{\mathrm{W}}=0.141~{}\mathrm{n}\mathrm{m}$, and the radius of one
argon atom is $r_{\mathrm{Ar}}=0.071~{}\mathrm{n}\mathrm{m}$ [31]. With the
atomic radii, we have

$\displaystyle d_{\mathrm{W}}$
$\displaystyle=r_{\mathrm{W}}+r_{\mathrm{Ar}}=0.212~{}\mathrm{n}\mathrm{m}$
(14a) $\displaystyle d_{\mathrm{Si}}$
$\displaystyle=r_{\mathrm{Si}}+r_{\mathrm{Ar}}=0.203~{}\mathrm{n}\mathrm{m}$
(14b)

The experiment was done under room temperature so $T=298~{}\mathrm{K}$.
Substituting the experimental pressures in Eq.(13), we have the mean free
paths of Si and W atoms in Table 2.

$p$ | 16 mTorr | 10 mTorr | 8 mTorr
---|---|---|---
$\mathrm{bar}{\lambda}_{\mathrm{W}}$(cm) | 0.966 | 1.55 | 1.93
$\mathrm{bar}{\lambda}_{\mathrm{Si}}$(cm) | 1.05 | 1.69 | 2.11
Table 2: Mean free paths of W and Si atoms at three different pressures

Since the mean free path of a molecule is much less than the working distance
$10~{}\mathrm{c}\mathrm{m}$, the adatoms will experience multiple collisions
on their way to the surface. It has been hypothesized that atom collisions in
the vapor play an increasingly important role with increasing pressure in the
experimental pressure regime [8]. Furthermore, the mean free path varies from
about 1 cm at 16 mTorr to about 2 cm at 8 mTorr. Thus increased pressure can
potentially lead to more random deposition direction and to nanoclusters
forming in the vapor. Collisions also reduce the kinetic energy of the
adatoms, presumably leading to smaller diffusion lengths. Less diffusion and
greater cluster sizes could then result in a rougher surface with a morphology
which stabilizes more quickly. We note that simulations and calculations have
shown the importance of adatom kinetic energy to diffusion lengths on the
surface [32, 33]. However, unlike previous works which examined the effect of
this behavior on cluster formation, here we evoke kinetic energy-enhanced
diffusion as a mechanism to smoothen the amorphous growth surface. Another
effect of increased adatom scattering in the vapor above the surface is to
randomize their arrival directions. Atom steering due to surface interactions
is a potentially important mechanism of surface roughening [34, 35]. It is
more effective for lower energy adatoms and grazing incidence approach to the
surface, so that it would increase with increasing pressure.

Comparing the results among three pressure conditions, the following
observations can be made:

  1. 1.

At 16 mTorr, it is seen that $T_{f}<2T_{\tau}$. Physically, this means that
the kinetic steady state is reached earlier than the fluctuation dynamics
reaches a steady state. At the measured length scales, this phenomenon
physically means that dynamical processes driven by adatom absorption and
diffusion continue to evolve after the average morphology is no longer
changing during the thin film growth. Zhou et al. [8] have suggested that the
deposited particles are in clusters at high pressures. Also, according to the
mean free path calculation above, particles at higher pressures are deposited
with lower kinetic energy. Under such circumstances, we hypothesize that the
effects of adatom diffusion on the surface structure are weaker, causing the
local morphology to become saturated earlier than in the case at lower
pressures. While deposited nanoclusters are still contributing to evolution of
the dynamic processes, the surface has already reached a steady state kinetic
roughness morphology.

  2. 2.

$q_{\parallel}$ dependence. Apart from $T_{\tau}$, the other three parameters
show obvious $q_{\parallel}$ dependence, and $T_{f}$ shows the strongest
$q_{\parallel}$ dependence. Therefore, the time scale to reach the kinetic
steady state varies significantly with lateral length scale. From the strong
$q_{\parallel}$ dependence of $T_{f}$ and weak $q_{\parallel}$ dependence of
$T_{\tau}$, we could deduce that, if this behavior continues beyond the
observed $q_{\parallel}$ window, there will be a crossing point at
$q_{\parallel}<0.23~{}\mathrm{n}\mathrm{m}^{-1}$ in the $16$ mTorr experiment
and at $q_{\parallel}>0.47~{}\mathrm{n}\mathrm{m}^{-1}$ in the $10$ mTorr
experiment, at which length scale the kinetic and dynamic steady states are
reached simultaneously. The $q_{\parallel}$ independence of $T_{\tau}$
suggests that the fluctuation dynamics reaches the steady state at the same
time on a wide range of length scales.

  3. 3.

Pressure dependence. Pressure-dependent behavior is most strongly seen in
$\tau_{\infty}$ and $T_{f}$. The time scales decrease with increasing
pressure. As the surface becomes rougher and the kinetic energy of adatoms
becomes lower with increasing pressure, the steady-state local morphology will
form faster. Under low pressure conditions, a smooth surface with sufficient
dynamical relaxation preserves the existing morphology over a greater time of
continuing film deposition.

##  IV Theory and Simulation

In this section we introduce a continuum model for comparison with the kinetic
and dynamic properties discussed in Section III. $x$ and $y$ axes are
equivalent in the experimental sputter deposition geometry, so an isotropic
equation should be applied. The Kuramoto-Sivashinsky (KS) equation was
invented to describe spontaneous instabilities in nonequilibrium physical and
chemical systems, such as chemical turbulence and burning fronts [9, 10].
Besides isotropy of $\mathbf{x}$, the translational symmetry of $t$ and $h$,
and inverse symmetry of $t$ should also be taken into consideration. All
physically possible terms up to the fourth order are included in the following
differential equation [13]:

$\begin{split}&\frac{\partial H}{\partial
t}=a_{1}\nabla^{2}H+a_{2}\nabla^{4}H\\\ &+a_{3}\nabla^{2}(\nabla
H)^{2}+a_{4}(\nabla H)^{2}+a_{5}M+F+\eta\end{split}$ (15)

with

$M=\begin{vmatrix}\partial_{x}^{2}H&\partial_{y}\partial_{x}H\\\
\partial_{x}\partial_{y}H&\partial_{y}^{2}H\end{vmatrix}$ (16)

In Eq.(15), $H$ is the actual height of the surface, $F$ is the average growth
velocity and noise $\eta$ conforms to stochastic uniform distribution.

The term with $M$ can be neglected compared to the $\nabla^{2}H$ term [13]. By
defining $h=H-Ft$ as the height measured from the average surface level, the
equation we adopted for simulation is expressed as

$\frac{\partial h}{\partial
t}=a_{1}\nabla^{2}h+a_{2}\nabla^{4}h+a_{3}\nabla^{2}(\nabla
h)^{2}+a_{4}(\nabla h)^{2}+\eta$ (17)

where the random noise $\eta$ satisfies

$\displaystyle\langle\eta(\mathbf{x},t)\rangle$ $\displaystyle=0$ (18a)
$\displaystyle\langle\eta(\mathbf{x},t)\eta(\mathbf{x^{\prime}},t^{\prime})\rangle$
$\displaystyle=2D\delta^{d}(\mathbf{x}-\mathbf{x^{\prime}})\delta(t-t^{\prime})$
(18b)

We used a square grid of size $L\times L$ with $L=1024$, and set the time step
to be $\Delta t=0.01$. One length unit is set as $1$ nm and one time unit
represents $1$ s.

Physical backgrounds of the five terms are discussed in Ref.[5, 13, 12, 11,
16, 14]. The signs of coefficients in Eq.(17) are

$a_{1},a_{2},a_{3}<0\text{ and }a_{4}>0$

The coefficient $a_{1}<0$ because incident particles on their way to the
surface will be deflected toward mounds by atomic attractions [35]. Thus, more
particles arrive at the surface where $\nabla^{2}h<0$ (Fig.2 in Ref.[13]),
contrary to the effect of a positive Laplacian term in the Edwards-Wilkinson
(EW) equation. In the EW equation, a positive $\nabla^{2}$ term comes from
surface desorption, but this mechanism is highly suppressed in the present
case because the experiment is done under room temperature. The $\nabla^{4}h$
term results from the surface diffusion in Mullin’s equation [36]. The third
term $a_{3}\nabla^{2}(\nabla h)^{2}$ arises from the equilibration of the
inhomogeneous surface concentration, and $a_{4}(\nabla h)^{2}$ can be derived
from the variation of surface inclinations. $(\nabla h)^{2}$ is usually
referred to as the Kardar-Parisi-Zhang (KPZ) term [37] and $\nabla^{2}(\nabla
h)^{2}$ is referred to as the conserved KPZ term. The KPZ term is non-
conserved so it can lead to excess velocity.

Figure 12: Schematic diagram of scattering inside the material

We investigated a wide range of parameters for Eq.(17). According to Ref.[12],
$a_{1}=-Fb$, where $b$ is the difference between the typical range of
interatomic forces and the equilibrium distance of the deposited particles.
Even if taking the ten times of the atomic radius to be $b$, that is to say,
$b=10r_{\mathrm{W}}=1.41~{}\mathrm{n}\mathrm{m}$, $a_{1}$ will only reach a
value of $-0.25~{}\mathrm{n}\mathrm{m}^{2}\cdot\mathrm{s}^{-1}$, whose
absolute value is too small to generate a peak in the early stage $I$-$t$
curves even if the other parameters are also small in magnitude. In order to
resemble the $I$-$t$ behavior in Fig.3, $\left|a_{1}\right|$ cannot be too
small because a negative $a_{1}$ generates the initial instability. There will
not be a peak in early stage $I$-$t$ curves unless a sufficiently negative
$a_{1}$ is applied to compete with the terms $a_{2}\nabla^{4}h$ and
$a_{4}(\nabla h)^{2}$ that stabilize the morphology evolution and cause the
scattering intensity to saturate. Correspondingly, although some theories were
developed to predict the other parameters [13], the core features in the
experiment must be reproduced with a proper selection of the other
coefficients compatible with $a_{1}$. Ultimately, we found the best general
agreement with experiment using the following parameters:

$\begin{split}a_{1}=-0.64~{}\mathrm{n}\mathrm{m}^{2}\cdot\mathrm{s}^{-1}&\text{;~{}~{}~{}}a_{2}=-1~{}\mathrm{n}\mathrm{m}^{4}\cdot\mathrm{s}^{-1}\text{;}\\\
a_{3}=-2~{}\mathrm{n}\mathrm{m}^{4}\cdot\mathrm{s}^{-1}&\text{;~{}~{}~{}}a_{4}=0.1~{}\mathrm{n}\mathrm{m}\cdot\mathrm{s}^{-1}\text{;~{}and~{}}\\\
D=3.8\times
10^{-3}&~{}\mathrm{n}\mathrm{m}^{4}\cdot\mathrm{s}^{-1}\end{split}$ (19)

Real space result of $t=2000~{}\mathrm{s}$ is shown in Fig.13. Mounds are
found on top of longer wavelength roughness, and these mounds coalesce and
split during the growth.

After simulating the surface in real space, a transformation to coherent X-ray
scattering data is required. The scattering pattern is therefore obtained from
the simulated surface using the formula [38, 17]

$\begin{split}I(q_{x},q_{y})&\propto\frac{1}{q_{z}^{\prime
2}L^{2}}\left|\int\int\,\mathrm{d}x\,\mathrm{d}ye^{-iq_{z}^{\prime}h(x,y)}e^{-i(q_{x}x+q_{y}y)}\right|^{2}\\\
&=\frac{1}{q_{z}^{\prime
2}L^{2}}\left|\mathcal{F}\\{e^{-iq_{z}^{\prime}h(x,y)}\\}\right|^{2}\end{split}$
(20)

where $q_{z}^{\prime}$ is the $z$ component of the wave vector change
$\mathbf{q^{\prime}}$ inside the material (Fig.12):

$\mathbf{q^{\prime}}=\mathbf{k_{f}^{\prime}}-\mathbf{k_{i}^{\prime}}$ (21)

$\mathbf{k}_{i}$ and $\mathbf{k}_{f}$ are incident and scattered beams outside
the material. $\mathbf{k}_{i}^{\prime}$ is refracted beam entering the
material and $\mathbf{k}_{f}^{\prime}$ is the scattered beam inside the
material. For scattering at the Yoneda wing, $\mathbf{k}_{f}^{\prime}$
propagates along (parallel to) the surface. According to Snell’s law and
scattering geometry,

$\begin{cases}n=\cos\alpha_{c}\\\ n\cos\alpha_{i}^{\prime}=\cos\alpha_{i}\\\
q_{z}^{\prime}=k_{i}^{\prime}\sin\alpha_{i}^{\prime}\\\ \end{cases}$ (22)

With small-angle approximation, $\alpha_{i}^{\prime}$ can be expressed as

$\alpha_{i}^{\prime 2}=\alpha_{i}^{2}-\alpha_{c}^{2}$ (23)

Therefore,

$q_{z}^{\prime}\approx\frac{2\pi\alpha_{i}^{\prime}}{\lambda}=\frac{2\pi\sqrt{\alpha_{i}^{2}-\alpha_{c}^{2}}}{\lambda}$
(24)

Substituting $\alpha_{i}=0.45\mathrm{\SIUnitSymbolDegree}$,
$\alpha_{c}=0.40\mathrm{\SIUnitSymbolDegree}$ and $\lambda=\frac{hc}{E}$ where
$E=7.38~{}\mathrm{k}\mathrm{eV}$ in Eq.(24), we have
$q_{z}^{\prime}=0.134~{}\mathrm{n}\mathrm{m}^{-1}$.

Figure 13: Simulation result $h(x,y)$ in the real space at $t=2000\mathrm{s}$.
The unit of height is nm. Figure 14: Simulation result: speckle averaged
intensity $I$-$t$ in the early stage Figure 15: Simulation result: speckle
averaged intensity $I$-$q_{\parallel}$ in the early stage evolving with time
Figure 16: 2TCF from simulation at
$q_{\parallel}=0.264~{}\mathrm{n}\mathrm{m}^{-1}$ Figure 17: The
$q_{\parallel}$ dependence of the four time parameters from the simulation

Now we are able to characterize the kinetic properties from the calculated
intensity evolution. The early stage kinetic behavior is shown in Fig.14 and
Fig.15. To compare the dynamics between the simulation and the experiment,
speckle correlation analysis was performed for the simulation in the same way
as for the X-ray data. Fig.16 gives an example of what the 2TCF looks like,
and the final results of the four time parameters are plotted in Fig.17.

The simulation is successful to some extent because many features have
resemblances to what was observed in the experiment. In Fig.14, it is firstly
observed in the continuum model that a peak develops in the structure factor
at early times. The peak in $I(q_{\parallel})$ (Fig.15) shifts to smaller
$q_{\parallel}$ with time and there is overlapping of subsequent curves. The
$q_{\parallel}$ dependence of $T_{f}$ and $\tau_{\infty}$ resembles Fig.11.

However, the simulation results also display a few disagreements with the
experiment. For the kinetics, peaks are observed in $I(t)$ in only part of the
$q_{\parallel}$ range instead of the entire $q_{\parallel}$ range from $0.23$
to $0.47~{}\mathrm{n}\mathrm{m}^{-1}$, and the peaks in Fig.15 are seen at
greater $q_{\parallel}$’s compared with the experimental result. Focusing on
the dynamic behavior, $2T_{\tau}$ decreases so significantly with
$q_{\parallel}$ that the surface reaches the dynamic steady state very quickly
after the deposition starts. Therefore $T_{f}$ is always greater than
$2T_{\tau}$. Furthermore, the strong $q_{\parallel}$ dependence of $T_{max}$
and $T_{\tau}$ does not accord with the experiments.

##  V Conclusion

In conclusion, the Co-GISAXS technique is used to study the kinetic and
dynamic behavior during a-WSi2 thin film deposition as a function of growth
pressure. After analyzing height-height structure factors and correlation
functions, we characterize the kinetic and dynamic processes by extracting
four time parameters, $T_{f}$, $T_{max}$, $T_{\tau}$ and $\tau_{\infty}$ with
physical interpretations in Table 1. The peaks in the $I$-$t$ curves lead to
$T_{max}$, implying that at each length scale, the roughness reaches a maximum
before relaxing down to a steady state. Kinetic properties reach a steady
state later at longer length scales because of mound formation and coarsening.
Comparing the kinetic behavior at different pressures, we can conclude that
the surface is rougher at higher pressure, as is expected because of
decreasing kinetic energy and increasingly random directions of adatoms. This
phenomenon can also be explained by the hypothesis that clusters of greater
sizes aggregate in the vapor at higher pressures, so the deposited particles
are actually nanoclusters instead of single atoms [8]. With such an
assumption, the mean free path would be even smaller. The surface at 8 mTorr
is too smooth to calculate the autocorrelation function. Therefore we could
not extract $T_{\tau}$ value for 8 mTorr.

We then discuss the $q_{\parallel}$ and pressure dependence in Fig.11. Apart
from $T_{\tau}$, the other three parameters show obvious $q_{\parallel}$
dependence under each pressure condition, and $T_{f}$ shows the strongest
$q_{\parallel}$ dependence. $T_{f}$ varies significantly with $q_{\parallel}$
because the average structure at greater lateral sizes of course takes longer
time to saturate, but $T_{\tau}$ does not show such strong $q_{\parallel}$
dependence, presumably because dynamical processes driven by atomic arrival,
diffusion and incorporation do not rely on the length scale we detect.
$T_{max}$ increases with increasing length scale because it takes longer to
form larger mounds, and the surface structure at longer length scales stays
correlated longer than at shorter length scales.

As for pressure dependence, a marked finding is the inversion between $T_{f}$
and $2T_{\tau}$. At pressures as high as 16 mTorr, the kinetic steady state is
found to be reached even earlier than the fluctuation dynamics reaches a
steady state. At high pressures, the adatoms are deposited with lower kinetic
energy, so they diffuse over a short length scale before they stop, and it
will take a shorter time for the local morphology to reach a steady state.
However, if local dynamical processes depend on surface roughness, then they
will take longer to stabilize under higher pressures because of the higher
roughness that develops. At low pressures, the surface roughness is small, so
the dynamical processes do not change much from what they were on the original
flat surface. Coherent X-ray scattering essentially reflects the height
fluctuation, so $2T_{\tau}$ decreases a little bit from 16 to 10 mTorr, and
$2T_{\tau}>T_{f}$ at 16 mTorr. $T_{max}$ hardly changes with pressure, which
shows that early stage coarsening through a definite lateral length scale does
not depend on pressure. A smoother surface structure remains correlated for a
longer time, so $\tau_{\infty}$ increases with decreasing pressure.

In the late stages of the experiment (Fig.6), the power law behavior of
$I(q_{\parallel})$ indicates a self-affine structure on the measured length
scales, which may be described by a continuum model. Moreover, the non-
exponential intensity evolution and the compressed exponential decay in Fig.8
require nonlinear terms. Therefore, Eq.(17) was adopted to describe the early
stage kinetic and dynamic behavior. The set of parameter in Eq.(19) represents
the best overall behavior among all of our simulations. This simulation result
resembles many features in the experiment. In part of the $q_{\parallel}$
ranges, $I(t)$ reaches a peak before reaching a lower steady-state value. The
peak in $I(q_{\parallel})$ shifts to smaller $q_{\parallel}$ with time and the
subsequent curves overlap, although the peaks are seen at greater
$q_{\parallel}$’s compared with the experiment. The $q_{\parallel}$ dependence
of $T_{f}$ and $\tau_{\infty}$ is consistent with the experiment. However, the
strong $q_{\parallel}$ dependence of $T_{max}$ is not consistent with the
experiment. The least agreement with experiment comes from $T_{\tau}$, which
almost vanishes in high $q_{\parallel}$. In addition, attempts to simulate the
pressure dependence behavior by simply varying the simulation coefficients
were not successful.

Improvements to the continuum model could be possible. According to Ref.[12],
$a_{5}=Fb^{2}$. If we adopt the ‘nanoclusters forming in vapor’ hypothesis,
$b$ will be in the order of 1 nm or above, so $a_{5}M$ in Eq.(15) cannot be
neglected. The nonuniformity of cluster sizes at high pressures may require a
more sophisticated treatment of the noise term. Besides, higher order
derivatives may also contribute to kinetics and dynamics, and different
magnitudes of discretization could be applied to different pressure conditions
if we could confirm the hypothetical dependence of the mean cluster sizes of
deposition on pressures.

The nature of the continuum model has itself limited the range of simulated
dynamic behavior and surface evolution. Eq.(17) has contained all possible
terms to $\mathcal{O}(\nabla^{4},(\nabla h)^{2})$ for an isotropic growth
process with definite physical sense. The shortcomings noted above were
observed for all simulation parameter sets examined. Restricted by local
behavior, local continuum models cannot describe nonlocal mechanisms like
shadowing. Furthermore, continuum models fail to track short length-scale
processes of particles or clusters because of discretization, losing part of
the correlation behavior especially at higher $q_{\parallel}$’s. In other
words, local morphology below a nanometer is lost because of the
discretization. An alternative solution could be atomistic models, which might
overcome some of these limitations.

###### Acknowledgements.

We thank Peco Myint for use of his XPCS analysis program and for helpful
conversations. The component of this work at BU was partly supported by the
National Science Foundation (NSF) under Grant No. DMR-1709380. At UVM, J.G.U.
and R.L.H. were partly supported by the U.S. Department of Energy (DOE) Office
of Science under Grant No. DE-SC0017802. This research used resources of the
Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science
User Facility operated for the DOE Office of Science by Argonne National
Laboratory under Contract No. DE-AC02-06CH11357.

## References

  * Olson [2000] G. B. Olson, Designing a new material world, Science 288, 993 (2000).
  * Movchan and Demchishin [1969] B. Movchan and A. Demchishin, Fiz met metalloved, J 28, 653 (1969).
  * Thornton [1977] J. A. Thornton, High rate thick film growth, Annual Review of Materials Science 7, 239 (1977).
  * Myint  _et al._ [2021] P. Myint, K. F. Ludwig Jr, L. Wiegart, Y. Zhang, A. Fluerasu, X. Zhang, and R. L. Headrick, Nanoscale dynamics during self-organized ion beam patterning of $\mathrm{Si}$. $\mathrm{I}$. $\mathrm{Ar}^{+}$ bombardment, Physical Review B 103, 195423 (2021).
  * Barabási  _et al._ [1995] A.-L. Barabási, H. E. Stanley, _et al._ ,  _Fractal concepts in surface growth_ (Cambridge University Press, 1995).
  * Headrick  _et al._ [2019] R. L. Headrick, J. G. Ulbrandt, P. Myint, J. Wan, Y. Li, A. Fluerasu, Y. Zhang, L. Wiegart, and K. F. Ludwig, Coherent $\mathrm{X}$-ray measurement of step-flow propagation during growth on polycrystalline thin film surfaces, Nature Communications 10, 1 (2019).
  * Sutton [2008] M. Sutton, A review of $\mathrm{X}$-ray intensity fluctuation spectroscopy, Comptes Rendus Physique 9, 657 (2008).
  * Zhou  _et al._ [2010] L. Zhou, Y. Wang, H. Zhou, M. Li, R. L. Headrick, K. MacArthur, B. Shi, R. Conley, and A. T. Macrander, Pressure-dependent transition from atoms to nanoparticles in magnetron sputtering: Effect on $\mathrm{WSi}_{2}$ film roughness and stress, Physical Review B 82, 075408 (2010).
  * Kuramoto [2003] Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Sivashinsky [1983] G. I. Sivashinsky, Instabilities, pattern formation, and turbulence in flames, Annual Review of Fluid Mechanics 15, 179 (1983).
  * Raible  _et al._ [2000a] M. Raible, S. J. Linz, and P. Hänggi, Amorphous thin film growth: Minimal deposition equation, Physical Review E 62, 1691 (2000a).
  * Raible  _et al._ [2000b] M. Raible, S. Mayr, S. J. Linz, M. Moske, P. Hänggi, and K. Samwer, Amorphous thin-film growth: Theory compared with experiment, EPL (Europhysics Letters) 50, 61 (2000b).
  * Raible  _et al._ [2001] M. Raible, S. J. Linz, and P. Hänggi, Amorphous thin film growth: Effects of density inhomogeneities, Physical Review E 64, 031506 (2001).
  * Pelliccione and Lu [2008] M. Pelliccione and T.-M. Lu, Evolution of thin film morphology, in  _Springer Series in Materials Science_ , Vol. 108 (Springer, 2008).
  * Villain [1991] J. Villain, Continuum models of crystal growth from atomic beams with and without desorption, Journal de Physique I 1, 19 (1991).
  * Lai and Sarma [1991] Z.-W. Lai and S. D. Sarma, Kinetic growth with surface relaxation: Continuum versus atomistic models, Physical Review Letters 66, 2348 (1991).
  * Mokhtarzadeh and Ludwig Jr [2017] M. Mokhtarzadeh and K. F. Ludwig Jr, Simulations of $\mathrm{Co}$-$\mathrm{GISAXS}$ during kinetic roughening of growth surfaces, Journal of Synchrotron Radiation 24, 1187 (2017).
  * Sneppen  _et al._ [1992] K. Sneppen, J. Krug, M. Jensen, C. Jayaprakash, and T. Bohr, Dynamic scaling and crossover analysis for the $\mathrm{K}$uramoto-$\mathrm{S}$ivashinsky equation, Physical Review A 46, R7351 (1992).
  * Sandy  _et al._ [2018] A. R. Sandy, Q. Zhang, and L. B. Lurio, Hard $\mathrm{X}$-ray photon correlation spectroscopy methods for materials studies, Annual Review of Materials Research 48, 167 (2018).
  * Yoneda [1963] Y. Yoneda, Anomalous surface reflection of $\mathrm{X}$ rays, Physical Review 131, 2010 (1963).
  * Renaud  _et al._ [2009] G. Renaud, R. Lazzari, and F. Leroy, Probing surface and interface morphology with grazing incidence small angle $\mathrm{X}$-ray scattering, Surface Science Reports 64, 255 (2009).
  * Salditt  _et al._ [1996] T. Salditt, D. Lott, T. Metzger, J. Peisl, R. Fischer, J. Zweck, P. Høghøj, O. Schärpf, and G. Vignaud, Observation of the $\mathrm{H}$uygens-principle growth mechanism in sputtered $\mathrm{W/Si}$ multilayers, EPL (Europhysics Letters) 36, 565 (1996).
  * Rainville  _et al._ [2015] M. G. Rainville, C. Wagenbach, J. G. Ulbrandt, S. Narayanan, A. R. Sandy, H. Zhou, R. L. Headrick, and K. F. Ludwig Jr, $\mathrm{Co}$-$\mathrm{GISAXS}$ technique for investigating surface growth dynamics, Physical Review B 92, 214102 (2015).
  * Arslan  _et al._ [2012] L. Ç. Arslan, C. Sanborn, E. Anzenberg, and K. F. Ludwig Jr, Evidence for $\mathrm{F}$amily-$\mathrm{M}$eakin dynamical scaling in island growth and coalescence during vapor phase deposition, Physical Review Letters 109, 106102 (2012).
  * Savitzky and Golay [1964] A. Savitzky and M. J. Golay, Smoothing and differentiation of data by simplified least squares procedures, Analytical Chemistry 36, 1627 (1964).
  * Williams and Watts [1970] G. Williams and D. C. Watts, Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function, Transactions of the Faraday Society 66, 80 (1970).
  * Ulbrandt  _et al._ [2016] J. G. Ulbrandt, M. G. Rainville, C. Wagenbach, S. Narayanan, A. R. Sandy, H. Zhou, K. F. Ludwig Jr, and R. L. Headrick, Direct measurement of the propagation velocity of defects using coherent $\mathrm{X}$-rays, Nature Physics 12, 794 (2016).
  * Schätzel  _et al._ [1988] K. Schätzel, M. Drewel, and S. Stimac, Photon correlation measurements at large lag times: improving statistical accuracy, Journal of Modern Optics 35, 711 (1988).
  * Culbertson and Burden [2007] M. J. Culbertson and D. L. Burden, A distributed algorithm for multi-tau autocorrelation, Review of Scientific Instruments 78, 044102 (2007).
  * Magatti and Ferri [2001] D. Magatti and F. Ferri, Fast multi-tau real-time software correlator for dynamic light scattering, Applied Optics 40, 4011 (2001).
  * [31] www.princeton.edu/ maelabs/mae324/glos324/index.htm.
  * Morales-Cifuentes  _et al._ [2014] J. R. Morales-Cifuentes, T. Einstein, and A. Pimpinelli, How “hot precursors” modify island nucleation: A rate-equation model, Physical Review Letters 113, 246101 (2014).
  * Gao  _et al._ [2012] D. Z. Gao, M. B. Watkins, and A. L. Shluger, Transient mobility mechanisms of deposited metal atoms on insulating surfaces: $\mathrm{Pd}$ on $\mathrm{MgO}$ (100), The Journal of Physical Chemistry C 116, 14471 (2012).
  * Van Dijken  _et al._ [1999] S. Van Dijken, L. C. Jorritsma, and B. Poelsema, Steering-enhanced roughening during metal deposition at grazing incidence, Physical Review Letters 82, 4038 (1999).
  * Shevchik [1973] N. J. Shevchik, Growth instabilities in the deposition of amorphous films, Journal of Non-Crystalline Solids 12, 141 (1973).
  * Mullins [1957] W. W. Mullins, Theory of thermal grooving, Journal of Applied Physics 28, 333 (1957).
  * Kardar  _et al._ [1986] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Physical Review Letters 56, 889 (1986).
  * Sinha  _et al._ [1988] S. Sinha, E. Sirota, S. Garoff, and H. Stanley, X-ray and neutron scattering from rough surfaces, Physical Review B 38, 2297 (1988).