# Heterogeneous Coefficients, Control Variables, and Identification of
Treatment Effects

Whitney K. Newey† and Sami Stouli§

###### Abstract.

Multidimensional heterogeneity and endogeneity are important features of
models with multiple treatments. We consider a heterogeneous coefficients
model where the outcome is a linear combination of dummy treatment variables,
with each variable representing a different kind of treatment. We use control
variables to give necessary and sufficient conditions for identification of
average treatment effects. With mutually exclusive treatments we find that,
provided the generalized propensity scores (Imbens, 2000) are bounded away
from zero with probability one, a simple identification condition is that
their sum be bounded away from one with probability one. These results
generalize the classical identification result of Rosenbaum and Rubin (1983)
for binary treatments.

† Department of Economics, MIT<EMAIL_ADDRESS>
$\S$ Department of Economics, University of Bristol<EMAIL_ADDRESS>
Keywords: Treatment effect; Multiple treatments; Heterogeneous coefficients;
Control variable; Identification; Conditional nonsingularity; Propensity
score.

##  1\. Introduction

Models that allow for multiple treatments are important for program evaluation
and the estimation of treatment effects (Cattaneo 2010; Heckman, Ichimura,
Smith, and Todd 1998; Imai and van Dyk 2004; Imbens 2000; Graham and Pinto
2018; Lechner 2001; Wooldridge 2004). A general class is heterogeneous
coefficient models where the outcome is a linear combination of dummy
treatment variables and unobserved heterogeneity. These models allow for
multiple treatment regimes, with each variable representing a different kind
of treatment. These models also feature multidimensional heterogeneity, with
the dimension of unobserved heterogeneity being determined by the number of
treatment regimes.

Endogeneity is often a problem in these models because we are interested in
the effect of treatment variables on an outcome, and the treatment variables
are correlated with heterogeneity. Control variables provide an important
means of controlling for endogeneity with multidimensional heterogeneity. For
treatment effects, a control variable is an observed variable that makes
heterogeneity and treatment variables independent when it is conditioned on
(Rosenbaum and Rubin, 1983).

We use control variables to give necessary and sufficient conditions for
identification of average treatment effects based on conditional
nonsingularity of the second moment matrix of the vector of dummy treatment
variables given the controls. This result is familiar in the binary treatment
case, but its generalization to multiple treatments appears to be new. With
mutually exclusive treatments we find that, provided the generalized
propensity scores (Imbens, 2000) are bounded away from zero with probability
one, a simple condition for identification is that their sum be bounded away
from one with probability one. These results provide an important
generalization of Rosenbaum and Rubin (1983)’s classical identification result
for binary treatments.

##  2\. Modeling of Treatment Effects

Let $Y$ denote an outcome variable of interest, and $X$ a vector of dummy
variables $X(t)$, $t\in\mathcal{T}\equiv\\{1,\dots,T\\}$, taking value one if
treatment $t$ occurs and zero otherwise, and $\varepsilon$ a structural
disturbance vector of finite dimension. We consider a heterogeneous
coefficients model of the form

(2.1) $Y=p(X)^{\prime}\varepsilon,\quad\quad
p(X)=(1,X(1),\dots,X(T))^{\prime}.$

This model is linear in the treatment dummy variables, with coefficients
$\varepsilon$ that need not be independent of $X.$ We assume that the vector
$\varepsilon$ is mean independent of the endogenous treatments $X$,
conditional on an observable control variable denoted $V$.

######  Assumption 1.

For the model in (2.1), there exists a control variable $V$ such that
$E[\varepsilon|X,V]=E[\varepsilon|V]$.

The Rosenbaum and Rubin (1983) treatment effects model is included as a
special case where $X\in\\{0,1\\}$ is a treatment dummy variable that is equal
to one if treatment occurs and equals zero without treatment, and

$p(X)=(1,X)^{\prime}.$

In this case $\varepsilon=(\varepsilon_{1},\varepsilon_{2})^{\prime}$ is two
dimensional with $\varepsilon_{1}$ giving the outcome without treatment, and
$\varepsilon_{2}$ being the treatment effect. Here the control variables in
$V$ would be observable variables such that Assumption 1 holds, i.e., the
coefficients $(\varepsilon_{1},\varepsilon_{2})$ are mean independent of
treatment conditional on controls; this is the unconfoundedness assumption of
Rosenbaum and Rubin (1983).

A central object of interest in model (2.1) is the average structural function
given by $\mu(X)\equiv p(X)^{\prime}E[\varepsilon]$; see Blundell and Powell
(2003) and Wooldridge (2005). This function is also referred to as the dose-
response function in the statistics literature (e.g., Imbens, 2000). When
$X\in\\{0,1\\}$ is a dummy variable for treatment, $\mu(0)$ gives the average
outcome if every unit remained untreated and $\mu(1)$ the average outcome if
every unit were treated, with $\mu(1)-\mu(0)$ being the average treatment
effect. In general, the average effect of some treatment $t\in\mathcal{T}$ is

$\mu(e_{t})-\mu(0_{T}),$

with $e_{t}=(0,\ldots,0,1,0,\ldots,0)^{\prime}$ defined as a $T$-vector with
all components equal to zero, except the $t$th, which is one, and $0_{T}$ a
$T$-vector of zeros.

The conditional mean independence assumption and the form of the structural
function $p(X)^{\prime}\varepsilon$ in (2.1) together imply that the control
regression function of $Y$ given $(X,V)$, $E[Y|X,V]$, is a linear combination
of the treatment variables:

(2.2) $E[Y\mid X,V]=p(X)^{\prime}E[\varepsilon\mid
X,V]=p(X)^{\prime}E[\varepsilon\mid V]=p(X)^{\prime}q_{0}(V),\;q_{0}(V)\equiv
E[\varepsilon\mid V].$

The average structural function can thus be expressed as a known linear
combination of $E[q_{0}(V)]$ from equation (2.2). By iterated expectations,

$p(X)^{\prime}E[q_{0}(V)]=p(X)^{\prime}E[E[\varepsilon\mid V]]=\mu(X).$

We use the varying coefficient structure of the control regression function
(2.2) and the implied linear form of $\mu(X)$ to give conditions that are
necessary as well as sufficient for identification.

##  3\. Identification Analysis

A sufficient condition for identification of the average structural function
is nonsingularity of the second moment matrix of the treatment dummies given
the controls,

$E[p(X)p(X)^{\prime}\mid V],$

with probability one. Under the maintained assumption that
$E[p(X)p(X)^{\prime}]$ is nonsingular, this condition is also necessary.

Theorem 1 states our first main result. The proofs of all formal results are
given in the Appendix.

######  Theorem 1.

Suppose that $E[\left\|\varepsilon\right\|^{2}]<\infty$,
$E[p(X)p(X)^{\prime}]$ is nonsingular, and Assumption 1 holds. Then:
$E[p(X)p(X)^{\prime}|V]$ is nonsingular with probability one if, and only if,
$\mu(X)$ is identified.

When $X\in\\{0,1\\}$ and $p(X)=(1,X)^{\prime}$, the identification condition
becomes the standard condition for the treatment effect model

$Y=\varepsilon_{1}+\varepsilon_{2}X,\quad E[\varepsilon\mid
X,V]=E[\varepsilon\mid
V],\quad\varepsilon\equiv(\varepsilon_{1},\varepsilon_{2})^{\prime}.$

The identification condition is that the conditional second moment matrix of
$(1,X)^{\prime}$ given $V$ is nonsingular with probability one, which is the
same as

(3.1) $\text{var}(X\mid V)=P(V)[1-P(V)]>0,\quad P(V)\equiv\Pr[X=1\mid V],$

with probability one, where $P(V)$ is the propensity score. Here we can see
that the identification condition is the same as $0<P(V)<1$ with probability
one, which is the standard identification condition.

With multiple treatments, because $p(X)$ includes an intercept, the
identification condition is the same as nonsingularity of the variance matrix
$\text{var}(X|V)$ with probability one. This result generalizes (3.1).

######  Theorem 2.

$E[p(X)p(X)^{\prime}|V]$ is nonsingular with probability one if, and only if,
the variance matrix $\text{var}(X|V)$ is nonsingular with probability one.

Considerable simplification occurs with mutually exclusive treatments, which
allows for the formulation of an equivalent condition for nonsingularity of
$E[p(X)p(X)^{\prime}|V]$ solely in terms of the generalized propensity scores
(Imbens, 2000). This result generalizes the standard identification condition
for binary $X$.

######  Theorem 3.

Suppose that $\Pr[X(t)=1|V]>0$ for each $t\in\mathcal{T}$ with probability
one. With mutually exclusive treatments, $E[p(X)p(X)^{\prime}|V]$ is
nonsingular with probability one if, and only if,

(3.2) $\Sigma_{s=1}^{T}\Pr[X(s)=1\mid V]<1,$

with probability one.

##  4\. Discussion

The heterogeneous coefficients formulation we propose for multiple treatment
effects reveals the central role of the conditional nonsingularity condition
for identification. Because this condition is in principle testable,
establishing that it is also necessary demonstrates testability of
identification (e.g., Breusch, 1986). With mutually exclusive treatments, the
formulation of the equivalent condition (3.2) thus relates testability of
identification to the generalized propensity scores. This is a generalization
of the relationship between testability of identification and the propensity
score in the binary treatment case.

Conditions that are both necessary and sufficient are also important for the
determination of minimal conditions for identification. In unpublished work
Wooldridge (2004) considers a restricted version of our model with
$E[X|\varepsilon,V]=E[X|V]$ and
$E[p(X)p(X)^{\prime}|\varepsilon,V]=E[p(X)p(X)^{\prime}|V]$, and shows that
$q_{0}(V)$ is identified if $E[p(X)p(X)^{\prime}|V]$ is invertible. The
additional conditional second moments assumption implies that his
identification condition differs from ours. Thus his result and proof do not
apply in our setting which only assumes conditional mean independence
$E[\varepsilon|X,V]=E[\varepsilon|V]$, and our results show that conditional
second moments independence is not necessary for identification in multiple
treatment effect models. Graham and Pinto (2018) consider a related approach
in work independent of the first version of this paper (Newey and Stouli,
2018) where the identification result (Lemma 1 in the Appendix) was derived.
The conditional nonsingularity condition we propose is weaker than their
identification condition, and we study necessity as well as sufficiency for
identification of average treatment effects.

The identification results we obtain here are of general interest for the vast
treatment effects literature (e.g., Imbens, 2004, Imbens and Wooldridge, 2009,
and Athey and Imbens, 2017, for reviews) and complement existing results on
identification of treatment effects.

##  Appendix A Proofs

###  A.1. Preliminary result

######  Lemma 1.

Suppose that $E[\left\|\varepsilon\right\|^{2}]<\infty$ and Assumption 1
holds. If $E[p(X)p(X)^{\prime}|V]$ is nonsingular with probability one then
$q_{0}\left(V\right)$ is identified.

###### Proof.

Let $\lambda_{\min}(V)$ denote the smallest eigenvalue of
$E[p(X)p(X)^{\prime}|V]$. Suppose that $\bar{q}\left(V\right)\neq
q_{0}\left(V\right)$ with positive probability on a set
$\widetilde{\mathcal{V}}$, and note that $\lambda_{\min}(V)>0$ on
$\mathcal{V}$ by assumption. Then

$\displaystyle
E\left[\left\\{p\left(X\right)^{\prime}\left\\{\bar{q}\left(V\right)-q_{0}\left(V\right)\right\\}\right\\}^{2}\right]$
$\displaystyle=$ $\displaystyle
E\left[\left\\{\bar{q}\left(V\right)-q_{0}\left(V\right)\right\\}^{\prime}E[p\left(X\right)p\left(X\right)^{\prime}\mid
V]\left\\{\bar{q}\left(V\right)-q_{0}\left(V\right)\right\\}\right]$
$\displaystyle\geq$ $\displaystyle
E\left[\left\|\bar{q}\left(V\right)-q_{0}\left(V\right)\right\|^{2}\lambda_{\min}\left(V\right)\right]$
$\displaystyle\geq$ $\displaystyle
E\left[1(V\in\mathcal{V}\cap\widetilde{\mathcal{V}})\left\|\bar{q}\left(V\right)-q_{0}\left(V\right)\right\|^{2}\lambda_{\min}\left(V\right)\right]$

By definition $\Pr[\widetilde{\mathcal{V}}]>0$ and
$\widetilde{\mathcal{V}}\subseteq\mathcal{V}$ so that
$\widetilde{\mathcal{V}}\cap\mathcal{V}=\widetilde{\mathcal{V}}$. Thus the
fact that
$\left\|\bar{q}\left(V\right)-q_{0}\left(V\right)\right\|^{2}\lambda_{\min}\left(V\right)$
is positive on $\widetilde{\mathcal{V}}\cap\mathcal{V}$ implies

$E\left[1(V\in\mathcal{V}\cap\widetilde{\mathcal{V}})\left\|\bar{q}\left(V\right)-q_{0}\left(V\right)\right\|^{2}\lambda_{\min}\left(V\right)\right]>0.$

We have shown that, for $\bar{q}\left(V\right)\neq q_{0}\left(V\right)$ with
positive probability on a set $\widetilde{\mathcal{V}}$,

$E\left[\left\\{p\left(X\right)^{\prime}\left\\{\bar{q}\left(V\right)-q_{0}\left(V\right)\right\\}\right\\}^{2}\right]>0,$

which implies $p\left(X\right)^{\prime}\bar{q}\left(V\right)\neq
p\left(X\right)^{\prime}q_{0}\left(V\right)$. Therefore, $q_{0}\left(V\right)$
is identified from $E[Y|X,V]$. ∎

###  A.2. Proof of Theorem 1

We first show that nonsingularity of $E[p(X)p(X)^{\prime}|V]$ with probability
one implies identification of $\mu(X)$. By Lemma 1, if
$E[p(X)p(X)^{\prime}|V]$ is nonsingular with probability one then $q_{0}(V)$
is identified, and hence $E[q_{0}(V)]$ also is. By $p(X)$ being a known
function, $p(X)^{\prime}E[q_{0}(V)]=\mu(X)$ is identified.

We now establish that nonsingularity of $E[p(X)p(X)^{\prime}|V]$ with
probability one is necessary for identification of $\mu(X)$. It suffices to
show that singularity of $E[p(X)p(X)^{\prime}|V]$ with positive probability
implies that $\mu(X)$ is not identified, i.e., there exists an observationally
equivalent $\overline{q}(V)\neq q_{0}(V)$ with positive probability such that
$p(X)^{\prime}E[\overline{q}(V)]\neq p(X)^{\prime}E[q_{0}(V)]$ with positive
probability. By nonsingularity of $E[p(X)p(X)^{\prime}]$ and linearity of
$\mu(X)$, the conclusion holds if, and only if, there exists an
observationally equivalent $\overline{q}(V)\neq q_{0}(V)$ with positive
probability such that $E[\overline{q}(V)]\neq E[q_{0}(V)]$.

Suppose that $E[p(X)p(X)^{\prime}|V]$ is singular with positive probability
and let $\Delta(V)$ be such that $E[p(X)p(X)^{\prime}|V]\Delta(V)=0$. We have
that $\Delta(V)\neq 0$ on a set $\widetilde{\mathcal{V}}$ with
$\Pr[\widetilde{\mathcal{V}}]>0$. For $J=T+1$, define
$\widetilde{\mathcal{V}}_{j}=\\{v\in\widetilde{\mathcal{V}}:\Delta_{j}(v)\neq
0\\}$, $j\in\\{1,\dots,J\\}$. Then
$\cup_{j=1}^{J}\widetilde{\mathcal{V}}_{j}=\\{v\in\widetilde{\mathcal{V}}:\Delta(v)\neq
0\\}=\widetilde{\mathcal{V}}$. Hence

$0<\Pr[\widetilde{\mathcal{V}}]=\Pr[\cup_{j=1}^{J}\widetilde{\mathcal{V}}_{j}]\leq\sum_{j=1}^{J}\Pr[\widetilde{\mathcal{V}}_{j}],$

which implies that $\Pr[\widetilde{\mathcal{V}}_{j^{*}}]>0$ for some
$j^{*}\in\\{1,\dots,J\\}$.

Set $\widetilde{\Delta}(v)=\Delta(v)$ for
$v\in\widetilde{\mathcal{V}}_{j^{*}}$, and $\widetilde{\Delta}(v)=0$
otherwise. By construction $\widetilde{\Delta}_{j^{*}}(V)\neq 0$, and letting

$\widetilde{\widetilde{\Delta}}(V)=\textrm{sign}\\{\widetilde{\Delta}_{j^{*}}(V)\\}\frac{\widetilde{\Delta}(V)}{||\widetilde{\Delta}(V)||},$

we have that $\widetilde{\widetilde{\Delta}}_{j^{*}}(V)>0$ on
$\widetilde{\mathcal{V}}_{j^{*}}$ and
$||\widetilde{\widetilde{\Delta}}(V)||=1$, and hence
$E[|\widetilde{\widetilde{\Delta}}(V)|]<\infty$ and
$E[\widetilde{\widetilde{\Delta}}_{j^{*}}(V)]\neq 0$. Therefore
$E[\widetilde{\widetilde{\Delta}}(V)]\neq 0$, which implies that
$E[q_{0}(V)+\widetilde{\widetilde{\Delta}}(V)]\neq E[q_{0}(V)]$. The result
follows. ∎

###  A.3. Proof of Theorem 2

The matrix $E[p(X)p(X)^{\prime}|V]$ is of the form

(A.1) $E[p(X)p(X)^{\prime}\mid V]=\begin{bmatrix}1&E[X^{\prime}\mid V]\\\
E[X\mid V]&E[XX^{\prime}\mid V]\end{bmatrix},$

and is positive definite if, and only if, the Schur complement of $1$ in (A.1)
is positive definite (Boyd and Vandenberghe, 2004, Appendix A.5.5.), i.e., if,
and only if,

$E[XX^{\prime}\mid V]-E[X\mid V]E[X^{\prime}\mid V]=\text{var}(X\mid V),$

is positive definite with probability one, as claimed.∎

###  A.4. Proof of Theorem 3

For a vector $w\in\mathbb{R}^{T}$, let $\text{diag}(w)$ denote the $T\times T$
diagonal matrix with diagonal elements $w_{1},\dots,w_{T}$. For mutually
exclusive treatments, the matrix $E[p(X)p(X)^{\prime}|V]$ is of the form

(A.2) $E[p(X)p(X)^{\prime}\mid V]=\begin{bmatrix}1&E[X^{\prime}\mid V]\\\
E[X\mid V]&\text{diag}(E[X\mid V])\end{bmatrix}.$

The matrix $\text{diag}(E[X|V])$ has diagonal elements
$E[X(t)|V]=\Pr[X(t)=1|V]>0$, for each $t\in\mathcal{T}$, and hence is positive
definite. Therefore, $E[p(X)p(X)^{\prime}|V]$ is positive definite if, and
only if, the Schur complement of $\text{diag}(E[X|V])$ in (A.2) is positive
definite (Boyd and Vandenberghe, 2004, Appendix A.5.5.), i.e., if, and only
if,

$\displaystyle 0<1-E[X^{\prime}\mid V]\text{diag}(E[X\mid V])^{-1}E[X\mid V]$
$\displaystyle=1-\Sigma_{s=1}^{T}E[X(s)\mid V]$
$\displaystyle=1-\Sigma_{s=1}^{T}\Pr[X(s)=1\mid V],$

with probability one, as claimed.∎

## References

  * Athey and Imbens (2017) Athey, S. and Imbens, G. (2017).  The state of applied econometrics: causality and policy evaluation.  Journal of Economic Perspectives. 31, 3–32.
  * Blundell and Powell (2003) Blundell, R., and Powell, J. L. (2003).  Endogeneity in nonparametric and semiparametric regression models.  Econometric society monographs. 36, 312–357.
  * Boyd and Vandenberghe (2004) Boyd, S. P. and Vandenberghe, L. (2004).  Convex Optimization.  Cambridge University Press.
  * Breusch (1986) Breusch, T. S. (1986).  Hypothesis testing in unidentified models.  Review of Economic Studies. 53, 635–651.
  * Cattaneo (2010) Cattaneo, M. (2010).  Efficient semiparametric estimation of multi-valued treatment effects under ignorability.  Journal of Econometrics. 155, 138–154.
  * Graham and Pinto (2018) Graham, B. S. and Pinto, C. C. D. X. (2018)  Semiparametrically efficient estimation of the average linear regression function.  eprint arXiv:1810.12511.
  * Heckman, Ichimura, Smith, and Todd (1998) Heckman, J. J., Ichimura, H., Smith, J. and Todd, P. (1998).  Characterizing selection bias using experimental data.  Econometrica. 66, 1017–1098.
  * Imai and van Dyk (2004) Imai, K. and van Dyk, D. A. (2004).  Causal inference with general treatment regimes: Generalizing the propensity score.  Journal of the American Statistical Association. 99, 854–866.
  * Imbens (2000) Imbens, G. W. (2000).  The role of the propensity score in estimating dose-response functions.  Biometrika. 87, 706–710.
  * Imbens (2004) Imbens, G. (2004).  Nonparametric estimation of average treatment effects under exogeneity.  Review of Economics and Statistics. 86, 4–29.
  * Imbens and Wooldridge (2009) Imbens, G. and Wooldridge, J. M. (2009).  Recent developments in the econometrics of program evaluation.  Journal of Economic Literature. 47, 5–86.
  * Lechner (2001) Lechner, M. (2001).  Identification and estimation of causal effects of multiple treatments under the conditional independence assumption.  In Econometric evaluation of labour market policies in Europe, 43–58. Physica, Heidelberg.
  * Newey and Stouli (2018) Newey, W. K. and Stouli, S. (2018).  Heterogenous coefficients, discrete instruments, and identification of treatment effects.  eprint arXiv:1811.09837.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika. 70, 41–55.
  * Wooldridge (2004) Wooldridge, J. M. (2004).  Estimating average partial effects under conditional moment independence assumptions.  Cemmap working paper CWP03/04.
  * Wooldridge (2005) Wooldridge, J. M. (2005).  Unobserved heterogeneity and the estimation of average partial effects. In  Identification and inference for econometric models, pp.27-55.


capsideboxfigure[][] [capbesidefigure]style=plain, heightadjust=object,
floatrowsep=none, capbesideposition=left,top, capbesidewidth=5mm,
capbesidesep=none

# Explore missing flow dynamics by physics-informed deep learning: the
parameterised governing systems

Hui Xu  School of Aeronautics and Astronautics, Shanghai Jiao Tong University,
Shanghai, 200240, China  Wei Zhang<EMAIL_ADDRESS>Science and Technology
on Water Jet Propulsion Laboratory, Marine and Research Institute of China,
Shanghai, China  Yong Wang<EMAIL_ADDRESS>Max Planck Institute for
Dynamics and Self-Organization, 37077 Göttingen, Germany

###### Abstract

Gaining and understanding the flow dynamics have much importance in a wide
range of disciplines, e.g. astrophysics, geophysics, biology, mechanical
engineering and biomedical engineering. As a reliable way in practice,
especially for turbulent flows, regional flow information such as velocity and
its statistics, can be measured experimentally. Due to the poor fidelity or
experimental limitations, some information may not be resolved in a region of
interest. On the other hand, detailed flow features are described by the
governing equations, e.g. the Navier-Stokes equations for viscous fluid, and
can be resolved numerically, which is heavily dependent on the capability of
either computing resources or modelling. Alternatively, we address this
problem by employing the physics-informed deep learning, and treat the
governing equations as a parameterised constraint to recover the missing flow
dynamics. We demonstrate that with limited data, no matter from experiment or
others, the flow dynamics in the region where the required data is missing or
not measured, can be reconstructed with the parameterised governing equations.
Meanwhile, a richer dataset, with spatial distribution of the control
parameter (e.g. eddy viscosity of turbulence modellings), can be obtained. The
method provided in this paper may shed light on data-driven scale-adaptive
turbulent structure recovering and understanding of complex fluid physics, and
can be extended to other parameterised governing systems beyond fluid
mechanics.

##  I Introduction

In fluid mechanics, acquiring the flow data is one of the key tasks in science
and engineering Duraisamy et al. (2018); Brunton et al. (2019). Analysing or
controlling the flow dynamics generally needs a rich data set or a continuous
representation of the flow field. Conventionally, the flow dynamics is often
conducted by either numerics or experiments. Numerically, the Navier-Stokes
equations (NSEs), which govern the flow dynamics in a form of partial
differential equations (PDEs), can be solved under typical assumptions or
specific conditions within a typical range of precision. For high Reynolds
number flows, or turbulent flows, directly solving the NSEs is extremely
challenging due to limited computer resources and large amounts of degrees of
freedom Moin and Mahesh (1998); Launder and Spalding (1974); Duraisamy et al.
(2018). To annihilate the restriction, various turbulence models, which are
coupled with the NSEs, are proposed to save computational cost Durbin (2018).
Meanwhile, suspicion around correctness and applicability of modelling arises,
which leads to unremitting pursuit of researches in turbulence modelling Ling
et al. (2016b); Xiao et al. (2016); Ling et al. (2016a); Wang et al. (2017),
although machine learning approaches have been used to enhance the predictive
capability Taghizadeh et al. (2020).

In contrast, nowadays, gain of data from experiments is reliable, which
generally has no restriction from the Reynolds number. However, due to the
precision and uncertainty of experimental equipment Xiao et al. (2016),
measured data can be obtained in a finite set of locations, and are precise
only in a limited-order statistical sense, which are not as sufficiently
accurate as the solutions described by the NSEs. In a statistical sense, the
experimental data can be generally regarded as an approximation of the
solutions of the NSEs. It turns out that leveraging the limited experimental
data to the underlying governing equations appears as a critical question in
modelling and recovering the flow dynamics. In order to address this problem,
we here consider the reformed governing equations where the control parameter
is regarded as an unknown, and employ the physics-informed learning technique
to preserve fundamental physical principles Taghizadeh et al. (2020); Shukla
et al. (2020); Raissi et al. (2020), by which the parameterised governing
equations and the missing flow dynamics can be physically determined at the
same time. To achieve the above, the parameterised governing equations are
treated as a physical constraint, which is used to form a loss function in the
training of neural networks Raissi et al. (2019a).

A data-driven discovery of PDEs, named PDE-FIND, was proposed by Rudy and his
collaborators. Especially, for the NSEs, the PDE-FIND was used to identify the
Reynolds number in the vorticity transport equation for fluid mechanicsRudy et
al. (2017). The research group of Karniadakis has made a great contribution in
physics-informed neural networks Raissi et al. (2019a); Meng and Karniadakis
(2020); Meng et al. (2020) and their applications to fluid flows. Generally,
physics-informed neural networks can be used for both forward and inverse
problems. For forward problem, Raissi et al. proposed the hidden fluid
mechanics Raissi et al. (2020), which provides a new way to learn velocity and
pressure fields from flow visualizations. Most recently Jin et al. proposed
the NSFNets to solve two different, velocity-pressure and vorticity-velocity,
formulations of the incompressible NSEs, in the sense of direct numerical
simulation for Kovasznay flow and turbulent channel flow Jin et al. (2021).
Other applications are for high speed inviscid flow Mao et al. (2020) and
vortex-induced vibration Raissi et al. (2019b). For inverse problems, physics-
informed neural networks were used to identify the Reynolds number from flow
data for the low Reynolds number flow past a cylinder Raissi et al. (2019a),
and to discover universal variable-order fractional model for turbulent
Couette flow Mehta et al. (2019).

Built on the above evidence, the generalisation of physics-informed neural
networks forms the basis of the parameterised governing system in exploring
missing flow physics and exploiting the essential of modelling parameters,
specially for turbulent flows. To overcome numerical stiffness, a recently
developed learning rate annealing algorithm Wang et al. (2020) is adopted in
this work. The proposed idea can be extended to a broad range of applications
where PDEs appear as physical constraints and the control parameter of the
PDEs plays a role in dominating the physics.

##  II Methodology

Mathematically, we consider the following system as a general description of
physical phenomena in a region $\Omega$:

$\mathcal{L}^{t}(\vec{p};\vec{p}_{0},\vec{g}_{\Gamma},\nu)=\vec{f},$ (1)

where the operator $\mathcal{L}^{t}$ represents the evolution of $\vec{p}$,
and $\vec{p}_{0}$ is the initial condition. $\vec{g}_{\Gamma}$ is the boundary
conditions where $\Gamma$ denotes the boundary $\partial\Omega$ of $\Omega$.
$\nu$ denotes the control parameter. For a steady system, the corresponding
representation can be written as
$\mathcal{L}^{\circ}(\vec{p};\vec{g}_{\Gamma},\nu)=\vec{f}$. In real-world
applications, in the sense of temporal (or ensemble) averaging or spatial
filtering, the control parameter can be represented as an effective one,
$\nu_{\rm eff}$, which is the summation of the physical $\nu$ and the modelled
$\nu_{t}$, while the original expression of the governing system is kept the
same. For instance, in turbulence modelling, at any level of turbulence,
closure of the averaged or filtered NSEs requires the parameter, eddy(subgrid)
viscosity, which is introduced in large-eddy simulations (LES), scale-
resolving simulations (SRS) and Reynolds-averaged Navier-Stokes (RANS) methods
Spalart (2009); Durbin (2018). Once the effective control parameter $\nu_{\rm
eff}$ is introduced as an unknown parameter,
$\mathcal{L}^{t}(\vec{p};\vec{p}_{0},\vec{g}_{\Gamma},\nu_{\rm eff})=\vec{f}$
or $\mathcal{L}^{\circ}(\vec{p};\vec{g}_{\Gamma},\nu_{\rm eff})=\vec{f}$ can
be regarded as the parameterised physical constraint which is used to explore
missing flow dynamics under the framework of the physics-informed neural
networks.

Let us consider the Leonardo’s self-portrait, shown in Figure 1 (A), as a
physical field. It is assumed that a part of Leonardo’s hair is missing. In
order to fix the portrait, we need to find a pattern which is consistent with
the missing one. Without referring the potential governing system of the
Leonardo’s portrait, by learning the hair pattern beyond the information-
missing region, it is possible to create a pattern which matches with the
pattern around the region. It is clear that in principle, there should exist
infinite possibilities. So, by virtue of limited information, how to
physically make it unique as possible as we can? This is termed as “the
parameterised physics-informed deep learning strategy”.

(A) Leonardo’s hair (B) Domain $\Omega$ and sub-domain $\Omega^{\prime}$

Figure 1: Illustration of missing information: (A) Leonardo’s self-portrait
with a part of curly hair missing; [Generated from
https://en.wikipedia.org/wiki/Leonardo_da_Vinci] (B) Mathematical
representation of the domain $\Omega$ and the sub-domain $\Omega^{\prime}$ on
which the information is missing.

As illustrated in Fig. 1(B), $\vec{p}$ is known on
$\Omega\setminus\Omega^{\prime}$ and unknown on $\Omega^{\prime}$.
Theoretically, in fluid mechanics, if the effective eddy viscosity $\nu_{\rm
eff}$ can be estimated or physically formulated, the missing $\vec{p}$ on
$\Omega^{\prime}$ can be solved from Eq. (1). In the framework of the
parameterised physics-informed learning, solving Eq. (1) is avoided and Eq.
(1) only leads us to directly find a physical unique solution on
$\Omega^{\prime}$ in terms of $\vec{p}$ on $\Omega\setminus\Omega^{\prime}$.
Let us turn back to the physical world of the fluid dynamics which is governed
by the PDEs, e.g. the incompressible NSEs as follows

$\mathcal{L}^{t}(\\{\mathbf{u},p\\};\mathbf{g},\nu)\equiv\left\\{\begin{array}[]{rl}\displaystyle\frac{\partial\mathbf{u}}{\partial
t}-\nu\nabla^{2}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}+\nabla
p&=\mathbf{f},\mbox{on }\Omega\\\\[10.00002pt]
\nabla\cdot\mathbf{u}&=0,\\\\[10.00002pt]
\mathbf{u}|_{\Gamma=\partial\Omega}&=\mathbf{g}.\end{array}\right.$ (2)

Theoretically, Eq. (2) can be used for all kinds of incompressible flows.
However, with extremely small viscosity $\nu$, and thus very high Reynolds
number, solving Eq. (2) is extra challenging in practical applications. In
another aspect, although some data can be obtained experimentally, it is also
challenging, if possible, to establish the direct connection between Eq. (2)
and the gained experimental data. In engineering applications, the RANS
equations, which are derived by the Reynolds decomposition of flow fields,
have been well validated for engineering applications by experimental data and
direct simulations of the NSEs Spalart (2009); Durbin (2018). The key question
in RANS is to model the Reynolds stress $\tau_{ij}$, which is represented as
follows

$\tau_{ij}=-\rho\overline{u_{i}^{\prime}u_{j}^{\prime}},$ (3)

where $\overline{(\cdot)}$ denotes the Reynolds average or the spatial
filtering, and $u^{\prime}_{i}=u_{i}-\overline{u}_{i}$. Similar to Fick’s law,
$\rho\overline{u_{i}^{\prime}u_{j}^{\prime}}$ can be modelled as

$\rho\overline{u_{i}^{\prime}u_{j}^{\prime}}=-\rho\nu_{t}(\frac{\partial\bar{u}_{i}}{\partial
x_{j}}+\frac{\partial\bar{u}_{j}}{\partial x_{i}}),$ (4)

where $\nu_{t}$ denotes the eddy viscosity. Although it is named as viscosity,
$\nu_{t}$ is an artificial parameter different from the molecular viscosity in
the original equations. Similar idea has been adopted to the large eddy
simulation (LES), where sub-grid eddy viscosity $\nu_{\rm sgs}$ is introduced
to model the subgrid scale turbulent flow. For hybrid RANS/LES such as detach
eddy simulation (DES) and scale adaptive simulation (SAS), $\nu_{t}$ in the
control system has different meaning in different region of the computational
domain, depending on the resolution of the turbulence scales. Normally,
$\nu_{t}$ is flow dependent. In the turbulence community, various modeling of
the parameter $\nu_{t}$ has been studied for several decades. With the big
data and machine learning, model constants of the RANS equations are tuned by
fitting Ling et al. (2016b); Poroseva et al. (2016); Edeling et al. (2017);
Duraisamy et al. (2018). However, compared with small data sets, observed
discrepancies are not explained any better with big ones Edeling et al.
(2017), and modelling $\nu_{t}$ so as the $\nu_{\rm sgs}$ in a universal sense
still remains as a persistent challenging because of the case-by-case
dependence. Luckily, engineering practices demonstrate that predictions of the
RANS equations for complex flows are mostly quite satisfactory. It indicates
that utilization of $\nu_{t}$ as an extra parameter can be a reliable way to
establish the fitting between averaged/filtered experiment data and solutions
of the corresponding control equations.

As indicated in the pioneering work Durbin (2018), representation of turbulent
transport by an eddy viscosity underlies most of the related theoretical work.
With this idea, the classical RANS equations can be rewritten as (without
considering the auxiliary equations)

$\mathcal{L}^{t}(\\{\bar{\mathbf{u}},p\\};g,\nu+\nu_{\eta})\equiv\left\\{\begin{array}[]{rl}\displaystyle\frac{\partial\bar{\mathbf{u}}}{\partial
t}-(\nu+\nu_{\eta})\nabla^{2}\bar{\mathbf{u}}+\bar{\mathbf{u}}\cdot\nabla\bar{\mathbf{u}}+\nabla\bar{p}&=\mathbf{f},\mbox{on
}\Omega\\\\[10.00002pt] \nabla\cdot\bar{\mathbf{u}}&=0,\\\\[10.00002pt]
\bar{\mathbf{u}}|_{\Gamma=\partial\Omega}&=\mathbf{g},\end{array}\right.$ (5)

where $\nu_{\rm\eta}$ is an artificial viscosity term introduced here for
parameterisation similar to that described above. From the control parameter
point of view, the data-based determination of $\nu_{\rm\eta}$ indicates
“scale adaption” of recovering turbulent structures, same as $\nu_{\rm t}$ in
DES or SAS. The above equation is the system of the operator equation, Eq.
(1), in a concrete way. For steady cases, the unsteady term
$\partial\bar{\mathbf{u}}/{\partial t}$ is neglected. It should be mentioned
that reformulating the NSEs does not imply employment of the RANS approach
since parameterisation of $\nu_{\rm eff}=\nu+\nu_{\rm\eta}$ appears only as an
undetermined parameter. The parameterisation is significantly distinguished
from the conventional tuning of model constants. The question remains as “how
to determine the parameter $\nu_{\rm eff}$ (or $\nu_{\rm\eta}$) regarded as a
spatial distributed variable?”.

With the paramaterised NSEs, assuming only limited data are available on a set
$\mathcal{X}\in\Omega\setminus\Omega^{\prime}$ (cardinality
$\mathcal{N}=|\mathcal{X}|$) of finite scattered measurement locations, we are
interested in the missing flow dynamics on $\Omega^{\prime}$ where the
specified point set $\mathcal{X}^{\prime}$ with cardinality
$\mathcal{N}^{\prime}=|\mathcal{X}^{\prime}|$ is chosen to represent the
dynamics. On $\mathcal{X}$, the set of measured averaged/filtered $\mathbf{u}$
is denoted by $\bar{\mathcal{S}}_{\mathcal{X}}$ and on $\Omega^{\prime}$, the
set of missing averaged/filtered $\mathbf{u}$ is denoted by
$\bar{\mathcal{S}}_{\mathcal{X}^{\prime}}$. Further, the undetermined
$\nu_{\rm eff}$ appears as a spatial variable defined on $\Omega$, and will be
inferred later. Here, the target is to approximate the functional relation
$(t,x,y,z)\mapsto(u,v,w,p,\nu_{\rm eff})$ by physics-informed neural network
$(t,x,y,z)\mapsto(e_{u},e_{v},e_{w},e_{\rm div})$. The parameterised NSEs are
encoded in $e_{u}$, $e_{v}$, $e_{w}$, and $e_{\rm div}$, which represent the
residuals of the momentum equations and the continuity equation, respectively.
To overcome numerical stiffness that lead to unbalanced back-propagated
gradients during model training, a recently developed learning rate annealing
algorithm is adopted Wang et al. (2020). The physics-informed residuals
$e_{u}$, $e_{v}$, $e_{w}$, and $e_{\rm div}$ are minimized with the mean
squared error loss function on the time set $\mathcal{T}$

${\rm Loss}=\mathbb{E}_{\mathcal{T}}[{\rm MSE}_{\rm Data}+{\rm MSE}_{\rm
PDEs}]$ (6)

where

${\rm MSE}_{\rm
Data}=\mathbb{E}_{\mathcal{X}}\left[\|\mathbf{u}(t;\mathcal{X})-\mathbf{u}_{(t;\mathcal{X})}\|^{2}\right],$
(7)

and

${\rm MSE}_{\rm
PDEs}=\mathbb{E}_{\hat{\mathcal{X}}}\left[e_{u}^{2}+e_{v}^{2}+e_{w}^{2}+e_{\rm
div}^{2}\right].$ (8)

Here $\hat{\mathcal{X}}\in\Omega$. The operator $\mathbb{E}_{\\#}[\cdot]$
means the average on the set ‘#’ in a distribution sense.
$\mathbb{E}_{\mathcal{T}}[\cdot]$ is an identity operator for stationary
modellings. Generally, the location cluster can be generated in a probability
sense. Once the above minimisation problem is solved, the extra variables $p$
and $\nu_{\rm eff}$ can be discovered.

##  III Numerical configuration

The diagram of the physics-informed neural network is illustrated in Fig. 2.
We here treat the parameter $\nu_{\rm\eta}$ as an output of the deep neutral
network. During the training process, the required weights of neural network
are determined by the loss function formed by the residual of the governing
equations. Once the training convergence is achieved, the variables $u$, $v$,
$w$, $p$ and $\nu_{\rm\eta}$ can be inferred. The training datasets were
generated by computational fluid dynamics (CFD) and the open-source software
OpenFOAM Weller et al. (1998) was used to implement required computations. In
two-dimension, only values of the velocity components $u$ and $v$ at mesh
nodes are considered as training data.

To check the ability of the method, we considered three different cases,
related to the parameterised governing equations of RANS, URANS and SAS for
two different turbulent flows, respectively. SAS denotes “scale-adaptive
simulation” Spalart (2009), and typically allows simulation of unsteady flows
with both RANS and LES content. It should be noted that only regional velocity
is used as training data, while velocity, pressure and artificial viscosity in
the entire domain are inferred.

Figure 2: Illustration of the physics-informed neural networks. The parameter
$\nu_{\rm\eta}$ is regarded as an output variable of the neural network.

Two flows, including a steady turbulent flow over a back-ward facing step
(RANS) and an unsteady turbulent flow past a cylinder (URANS and SAS), were
involved. For the steady one, there are $1.2\times 10^{4}$ cell nodes in
total, while the random sampling nodes are $3\times 10^{3}$ for training data
(on $\Omega\setminus\Omega^{\prime}$) and $4.3\times 10^{3}$ constraints of
the parameterised governing system (on $\Omega$), respectively. A neural
network with 6 hidden layers and 50 neurons per layer is considered. $2\times
10^{5}$ epochs are used. For the unsteady flow past a cylinder, there are 351
snapshots and $5.1\times 10^{5}$ cell nodes on every snapshot, thus
$1.7901\times 10^{8}$ in total. The random sampling nodes are $1.2\times
10^{7}$ and $1.4\times 10^{7}$ respectively, on all time snapshots.
Corresponding neural network has 10 hidden layers and 96 neurons per layer.
The learning rate is gradually decreased from $1.0\times 10^{-3}$ to
$1.0\times 10^{-6}$ with 400 epochs at each stage. The Adam optimizer is
adopted for all cases.

##  IV Results

We first consider a steady turbulent flow over a backward facing step,
experimentally studied by Pitz and Daily date back to 1981 Pitz and Daily
(1981). The Reynolds number is $5\times 10^{5}$ based on the step height. The
SST $k-\omega$ model was adopted to generate the data. Here SST denotes
“shear-stress transport” Spalart (2009). The contours of velocity components
along the streamwise and vertical directions are given in Figs. 3(a) and 4(a),
respectively. It is assumed that in the dash-lined rectangular box, denoted by
$\Omega^{\prime}$, the flow information therein is missing and thus not used
as training data.

After the neural network is trained with the data outside $\Omega^{\prime}$,
together with constrains of the governing equations, the flow variables in the
whole domain $\Omega$ are inferred. In Figs. 3(b) and 4(b), the inferred
velocity components by the physics-informed deep learning are given. It is
clear that in the box region, we get the similar velocity fields as
illustrated in Figs. 3(a) and 4(a), demonstrating that the averaged flow
dynamics is mimicked correctly with the proposed parameterised NSEs. In Fig.
5, the inferred pressure field is in good agreement with the reference one
from the SST $k-\omega$ model. It is interesting, on the other hand, to
compare the inferred artificial viscosity $\nu_{\rm\eta}$ with eddy viscosity
$\nu_{t}$ empirically determined by the SST $k-\omega$ model. Although they
have different meaning thus are different in both magnitude and spatial
distribution, as shown in Fig. 6, their distributions are related to the flow
structure and both are relatively large in the region with separation bubble.

Figure 3: Contours of the velocity component $u$ in the streamwise direction
for the steady turbulent flow over a backward facing step. (Top) from the SST
$k-\omega$ model; (Bottom) inferred by the physics-informed neural network.
The dash-lined box indicates the region where the flow information is missing.

Figure 4: Contours of the velocity component $v$ in the vertical direction for
the steady turbulent flow over a backward facing step. (Top) from the SST
$k-\omega$ model; (Bottom) inferred by the physics-informed neural network.
The dash-lined box indicates the region where the flow information is missing.

Figure 5: The pressure fields for the steady turbulent flow over a backward
facing step. (Top) from the SST $k-\omega$ model; (Bottom) inferred by the
physics-informed neural network.

Figure 6: Comparison between the eddy viscosity and the artificial viscosity
for the steady turbulent flow over a backward facing step. (Top) eddy
viscosity $\nu_{t}$ from the SST $k-\omega$ model; (Bottom) artificial
viscosity $\nu_{\rm\eta}$ inferred by the physics-informed neural network.

For the same steady turbulent flow, we specially consider a case where the
parameter $\nu_{\rm\eta}$ is removed and thus the original NSEs are used as a
physics constraint during the learning. The energy contours of the inferred
velocity field $\mathbf{u}_{\rm NSEs}$ and the ‘perturbation’ field
($\mathbf{u}_{\rm NSEs}-\bar{\mathbf{u}}$) are given in Fig. 7. Special
structures are observed around the tip of the separation bubble. Compared with
the solution based on the parameterised NSEs, this significant inconsistency
provides some potential information about the region where the most unstable
modes or the energetic turbulent structures emerge around the flow separation
region.

Figure 7: The total energy and ‘perturbation’ energy of the inferred velocity
for the steady turbulent flow over a backward facing step. (Top) the inferred
energy with the original NSEs; (Bottom) the ‘perturbation’ energy. The white
dashed-line indicates the separation region of the averaged field, while the
red dash-lined box is the region where the flow information is missing.

As a further demonstration, we consider the time-dependent turbulent flow past
a cylinder. The Reynolds number is, $Re=u_{in}D/\nu=10^{5}$, where $u_{in}$
and $D$ are the inlet velocity and diameter of the cylinder, respectively. The
reference flows are generated respectively by the SST $k-\omega$ and SST-SAS
turbulence models. It should be emphasized that, unlike other traditional RANS
models (e.g. the SST $k-\omega$ model), the SST-SAS model has a pure RANS
nature but achieves LES behavior. With this time-dependent turbulent flow, we
demonstrate that the proposed methodology has a great adaptability to
turbulence models for unsteady flow dynamics. In practice, the reference flow
data obtained by the SST-SAS and SST $k-\omega$ turbulence models can be
regarded as the “filtered" experimental data.

It is assumed, again, that in the system there exists a local region where the
flow dynamics is missing, as indicated by a dashed rectangular box in Fig. 8.
In order to employ the established methodology, we consider a series of flow
field snapshots obtained by the SST-SAS model, and only a part of the velocity
field beyond the dashed-rectangular box on every snapshot is employed for
training. In Fig. 8(Mid), the reconstruction of the missing flow dynamics at
$t=11.96$ in the dashed-rectangular box is presented. Surprisingly, compared
with Fig. 8(Top), the reconstructed missing flow coincides with the reference
flow very well, even with complex vortex structure in the box. Qualitatively,
we show the differences between the reference fields and the inferred ones in
Fig. 8(Bottom). Relative low deviation is observed. Meanwhile, the profiles of
the lift and drag coefficients are given in Fig. 9. The inferred $C_{L}$
profile is consistent with the reference one. However, there exist some
differences between the inferred $C_{D}$ and the reference one. As observed
from the comparison, the unsteadiness (or “irregularity” ) of the $C_{D}$
profile from the inferred data is reconstructed fairly well. Further, we show
the distribution of the viscosity in Fig. 10. It is observed that the patterns
of the viscosity from both the SST-SAS and neural network, are consistent with
the vortex shedding pattern, but with different magnitudes.

Figure 8: Unsteady turbulent flow past a cylinder: (Top) the instantaneous
reference velocity and pressure fields at $t=11.96$, generated by the SST-SAS
turbulence model; (Mid) the inferred fields at the same time instant with the
missing dynamics; (Bottom) the difference between the reference fields and the
inferred fields. The red dashed-line indicates the region where the flow
dynamics is missing.

Figure 9: Lift and drag coefficients over the cylinder: (Top) profiles of the
lift coefficient; (Bottom) profiles of the drag coefficient. The solid curves
are generated from the reference data by the SST-SAS turbulence model, while
the dashed curves generated from the inferred fields.

Figure 10: Comparison between the eddy viscosity and the artificial viscosity
for the unsteady turbulent flow past a cylinder: (Left) eddy viscosity
$\nu_{t}$ from SST-SAS turbulence model; (Right) the parameter $\nu_{\rm\eta}$
inferred by the physics-informed neural network.

Accordingly, for the flow past a cylinder, instantaneous results at $t=7$ with
the SST $k-\omega$ model are given in Figs. 11-13. We see that using the
methodology with the physics-informed deep learning, the missing flow dynamics
whose location is assigned arbitrarily, is recovered properly with good
accuracy.

Figure 11: Unsteady turbulent flow past a cylinder: (Top) the instantaneous
reference velocity and pressure fields at $t=7$, generated by the SST
$k-\omega$ turbulence model; (Mid) the inferred fields at the same time
instant with the missing dynamics; (Bottom) the difference between the
reference fields and the inferred fields. The red dashed-line indicates the
region where the flow dynamics is missing.

Figure 12: Lift and drag coefficients over the cylinder: (Top) profiles of the
lift coefficient; (Bottom) profiles of the drag coefficient. The solid curves
are generated from the reference case by the the SST $k-\omega$ turbulence
model, while the dashed curves generated from the inferred fields.

Figure 13: Comparison between the eddy viscosity and the artificial viscosity
for the unsteady turbulent flow past a cylinder: (Left) eddy viscosity
$\nu_{t}$ from SST $k-\omega$ turbulence model; (Right) the parameter
$\nu_{\rm\eta}$ inferred by the physics-informed neural network.

With the cases investigated, we demonstrate that when the parameter
$\nu_{\rm\eta}$ is introduced as an output variable of the neural network, the
parameterised governing system of flow dynamics has a strong adaptivity to
approximate the data from the turbulence models. Significantly, the missing
data in a randomly specified region can be inferred by the physics-informed
deep learning. Given the reality that reasonable accurate velocity and
pressure are predicted by the methodology, the inferred $\nu_{\rm\eta}$
indicates that the eddy viscosity in the conventional turbulence models
probably has an over-dissipation effect on the flow structures and a different
spatial distribution from the inferred one. This discovery implies that the
physics-preserved data driven techniques may have a transformative influence
on the turbulence closure modelling.

##  V Conclusion and Discussion

The methodology proposed in this paper provides a way to explore the missing
flow dynamics, which is well validated by the data from the turbulence models.
We show that from the single parameterisation of the Navier-Stokes equations,
the feasibility and potential of the method are significant for the “averaged”
or “filtered” data, no matter obtained experimentally or numerically.
Surprisingly, the current methodology allows to use the parameterised
governing equations without extra modellings, which is shaping a perspective
in exploring the dynamics beyond the description of the original governing
system, although the deviation of exploring missing dynamics exists. From this
point of view, the parameterisation, in fact, provides a solid foundation of
accurately describing practical complex flow dynamics, which may not be
accurately described by the traditional turbulence modellings where the
artificial viscosity $\nu_{\rm t}$ needs to be further tuned.

The proposed methodology provides a new way for further calibrations of the
turbulence closures, which is probably difficult for traditional case-by-case
machine learning with simply tuning the scalar modelling parameters. As
indicated from the theory of turbulence modelling, due to the significance of
the parameter $\nu_{\rm\eta}$, the proposed method is naturally termed as
“data-driven scale-adaptive turbulent structure recovering”.

Although only 2D incompressible flows were considered as examples in this
work, the proposed methodology can be extended to 3D incompressible or
compressible flows. While learning from the parameterised Navier-Stokes
equations is successful in coping with the traditional fluid dynamics
problems, this, in turn, implies generalisation of the employed methodology
could enable exploration and exploitation in a broad range of complex systems,
e.g. weather prediction, ocean recirculation and biofluid mechanics.

## Acknowledgments

The authors would like to acknowledge support from National Numerical Wind
Tunnel Project (Grant NO. NNW2019ZT4-B09) and the National Natural Science
Foundation of China (Grant NOs. 91852106).

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## references

## References

  * Brunton et al. [2019] S. L. Brunton, B. R. Noack, and P. Koumoutsakos.  Machine learning for fluid mechanics.  _Annual Review of Fluid Mechanics_ , 52:477 – 598, 2019\.
  * Duraisamy et al. [2018] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence modeling in the age of data.  _Annual Review of Fluid Mechanics_ , 51:357 – 77, 2018\.
  * Durbin [2018] P. A. Durbin.  Some recent developments in turbulence closure modeling.  _Annual Review of Fluid Mechanics_ , 50:77 – 103, 2018\.
  * Edeling et al. [2017] W. N. Edeling, G. Iaccarino, and P. Cinnella.  Data-free and data-driven rans predictions with quantified uncertainty.  _Flow, Turbulence and Combustion_ , 100:593 – 616, 2017\.
  * Jin et al. [2021] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  _Journal of Computational Physics_ , 426:109951, 2021.
  * Launder and Spalding [1974] B.E. Launder and D.B. Spalding.  The numerical computation of turbulent flows.  _Computer Methods in Applied Mechanics and Engineering_ , 3(2):269 – 289, 1974.
  * Ling et al. [2016a] J. Ling, A. Kurzawski, and J. Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807:155 – 166, 2016a.
  * Ling et al. [2016b] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807:155–166, 2016b.
  * Mao et al. [2020] Zhiping Mao, Ameya D. Jagtap, and George E Karniadakis.  Physics-informed neural networks for high-speed flows.  _Computer Methods in Applied Mechanics and Engineering_ , 360:112789, 2020.
  * Mehta et al. [2019] Pavan Pranjivan Mehta, Guofei Pang, Fangying Song, and George Em Karniadakis.  Discovering a universal variable-order fractional model for turbulent couette flow using a physics-informed neural network.  _Fractional Calculus and Applied Analysis_ , 22(6):1675–1688, 2019.
  * Meng and Karniadakis [2020] Xuhui Meng and George E Karniadakis.  A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems.  _Journal of Computational Physics_ , 401:109020, 2020.
  * Meng et al. [2020] Xuhui Meng, Zhen Li, Dongkun Zhang, and George E Karniadakis.  PPINN : Parareal physics-informed neural network for time-dependent PDEs.  _Computer Methods in Applied Mechanics and Engineering_ , 370:113250, 2020.
  * Moin and Mahesh [1998] P. Moin and K. Mahesh.  Direct numerical simulation: A tool in turbulence research.  _Annual Review of Fluid Mechanics_ , 39:538 – 578, 1998\.
  * Pitz and Daily [1981] R. W. Pitz and J. W. Daily.  An experimental study of combustion: The turbulent structure of a reacting shear layer formed at a rearward-facing step.  _Tech. rep., NASA, Lewis research center (contract number 165427)_ , 1981.
  * Poroseva et al. [2016] S. Poroseva, F. J. D. Colmenares, and S. Murman.  On the accuracy of rans simulations with dns data.  _Physics of Fluids_ , 28:115102, 2016.
  * Raissi et al. [2019a] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686 – 707, 2019a.
  * Raissi et al. [2019b] Maziar Raissi, Zhicheng Wang, Michael S. Triantafyllou, and George Em Karniadakis.  Deep learning of vortex-induced vibrations.  _Journal of Fluid Mechanics_ , 861:119–137, 2019b.
  * Raissi et al. [2020] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  _Science_ , 367(6481):1026–1030, 2020.
  * Rudy et al. [2017] Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Data-driven discovery of partial differential equations.  _Science Advances_ , 3(4):e1602614, 2017.
  * Shukla et al. [2020] Khemraj Shukla, Patricio Clark Di Leoni, James Blackshire, Daniel Sparkman, and George Em Karniadakis.  Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks.  _Journal of Nondestructive Evaluation_ , 39:61, 2020.
  * Spalart [2009] P. R. Spalart.  Detached-eddy simulation.  _Annual Review of Fluid Mechanics_ , 41:181 – 202, 2009\.
  * Taghizadeh et al. [2020] S. Taghizadeh, F. D. Witherden, and S. S. Girimaji.  Turbulence closure modeling with data–driven algorithms: physical compatibility and consistency considerations.  _New Journal of Physics_ , 2020.
  * Wang et al. [2017] JianXun Wang, JinLong Wu, and Heng Xiao.  Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data.  _Physical Review Fluids_ , 2(3):1–22, Mar 2017\.
  * Wang et al. [2020] Sifan Wang, Yujun Teng, and Paris Perdikaris.  Understanding and mitigating gradient pathologies in physics-informed neural networks.  _arXiv preprint arXiv:2001.04536_ , 2020.
  * Weller et al. [1998] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  _Computers in Physics_ , 12(6):620 – 631, 1998\.
  * Xiao et al. [2016] H. Xiao, J.L. Wu, J.X. Wang, R. Sun, and C.J. Roy.  Quantifying and reducing model-form uncertainties in reynolds-averaged navier–stokes simulations: A data-driven, physics-informed bayesian approach.  _Journal of Computational Physics_ , 324:115 – 136, 2016\.


# Some identities involving degenerate Stirling numbers arising from normal
ordering

Taekyun Kim  Department of Mathematics, Kwangwoon University, Seoul 139-701,
Republic of Korea<EMAIL_ADDRESS>and  DAE SAN KIM  Department of Mathematics,
Sogang University, Seoul 121-742, Republic of Korea<EMAIL_ADDRESS>
###### Abstract.

In this paper, we derive some identities and recurrence relations for the
degenerate Stirling numbers of the first kind and of the second kind which are
degenerate versions of the ordinary Stirling numbers of the first kind and of
the second kind. They are deduced from the normal oderings of degenerate
integral powers of the number operator and their inversions, certain relations
of boson operators and from the recurrence relations of the Stirling numbers
themselves. Here we note that, while the normal ordering of an integral power
of the number operator is expressed with the help of the Stirling numbers of
the second kind, that of a degenerate integral power of the number operator is
represented by means of the degenerate Stirling numbers of the second kind.

###### Key words and phrases:

normal odering; degenerate Stirling numbers of the first kind; degenerate
Stirling numbers of the second kind

###### 2010 Mathematics Subject Classification:

11B73; 11B83

##  1\. Introduction

The Stirling number of the second $S_{2}(n,k)$ is the number of ways to
partition a set of $n$ objects into $k$ nonempty subsets (see (6)). The
(signed) Stirling number of the first kind $S_{1}(n,k)$ is defined in such a
way that the number of permutations of $n$ elements having exactly $k$ cycles
is the nonnegative integer $(-1)^{n-k}S_{1}(n,k)=|S_{1}(n,k)|$ (see (5)). The
degenerate Stirling numbers of the second kind $S_{2,\lambda}(n,k)$ (see (9))
and of the first kind $S_{1,\lambda}(n,k)$ (see (7)) appear most naturally
when we replace the powers of $x$ by the generalized falling factorial
polynomials $(x)_{k,\lambda}$ in the defining equations (see (5), (6), (7),
(9)).

Carlitz initiated a study of degenerate versions of some special numbers and
polynomials in [2], where the degenerate Bernoulli and Euler numbers were
investigated. It is remarkable that in recent years intensive studies have
been done for degenerate versions of quite a few special polynomials and
numbers and have yielded many interesting results (see [6,7] and the
references therein). They have been explored by various methods, including
mathematical physics, combinatorial methods, generating functions, umbral
calculus techniques, $p$-adic analysis, differential equations, special
functions, probability theory and analytic number theory. It turns out that
the degenerate Stirling numbers play an important role in this exploration for
degenerate versions of many special numbers and polynomials. The normal
ordering of an integral power of the number operator $a^{+}a$ in terms of
boson operators $a$ and $a^{+}$ can be written in the form

$(a^{+}a)^{k}=\sum_{l=0}^{k}S_{2}(k,l)(a^{+})^{l}a^{l}.$

The normal ordering of the degenerate $k$th power of the number operator
$a^{+}a$, namely $(a^{+}a)_{k,\lambda}$, in terms of boson operators $a,a^{+}$
can be written in the form

(1) $(a^{+}a)_{k,\lambda}=\sum_{l=0}^{k}S_{2,\lambda}(k,l)(a^{+})^{l}a^{l},$

where the generalized falling factorials $(x)_{n,\lambda}$ are given by (1).

By inversion, from (1) we obtain

(2) $(a^{+})^{k}a^{k}=\sum_{l=0}^{k}S_{1,\lambda}(k,l)(a^{+}a)_{l,\lambda}.$

The aim of this paper is to derive some identities and recurrence relations
for the degenerate Stirling numbers of the first kind and of the second kind
by using the normal ordering in (1) and its inversion in (2), certain
relations of boson operators and the recurrence relations of the Stirling
numbers themselves. In more detail, our main results are as follows. Firstly,
we derive some recurrence relation and identity involving the degenerate
Stirling numbers of the second kind in Theorems 3 and 5 by using (1) and
certain relations for boson operators. Secondly, we find some recurrence
relation and identity involving the degenerate Stirling numbers of the first
kind in Theorems 4 and 6 by exploiting (2) and certain relations for boson
operators. Thirdly, we investigate some recurrence relations for the
degenerate Stirling numbers of the second kind in Theorems 7 and 9 by using
the recurrence relation in (15), and those for the degenerate Stirling numbers
of the first kind in Theorems 8 and 11 by using the recurrence relation in
(8). For the rest of this section, we recall the facts that are needed
throughout this paper.

For any $\lambda\in\mathbb{R}$, the degenerate exponential functions are
defined by

(3)
$e_{\lambda}^{x}(t)=\sum_{k=0}^{\infty}(x)_{n,\lambda}\frac{t^{k}}{k!},\quad(\mathrm{see}\
[6,7]),$

where

(4)
$(x)_{0,\lambda}=1,\quad(x)_{k,\lambda}=x(x-\lambda)\cdots\big{(}x-(k-1)\lambda\big{)},\quad(k\geq
1).$

When $x=1$, we use the notation $e_{\lambda}(t)=e_{\lambda}^{1}(t)$. It is
well known that the Stirling numbers of the first kind are defined by

(5) $(x)_{n}=\sum_{k=0}^{n}S_{1}(n,k)x^{k},\quad(n\geq 0),\quad(\mathrm{see}\
[1,2,3,9,10]),$

where $(x)_{0}=1,\ (x)_{n}=x(x-1)\cdots(x-n+1),\quad(n\geq 1)$.

As the inversion formula of (5), the Stirling numbers of the second kind are
defined by

(6) $x^{n}=\sum_{k=0}^{n}S_{2}(n,k)(x)_{k},\quad(n\geq 0).$

In [6], the degenerate Stirling numbers of the first kind are defined by

(7) $(x)_{n}=\sum_{k=0}^{n}S_{1,\lambda}(n,k)(x)_{k,\lambda},\quad(n\geq 0).$

From (7), we note that

(8)
$S_{1,\lambda}(k+1,l)=S_{1,\lambda}(k,l-1)-(k-l\lambda)S_{1,\lambda}(k,l),$

where $k,l\in\mathbb{N}$ with $k\geq l$.

As the inverse formula of (7), the degenerate Stirling numbers of the second
kind are defined by

(9) $(x)_{n,\lambda}=\sum_{k=0}^{n}S_{2,\lambda}(n,k)(x)_{k},\quad(n\geq
0),\quad(\mathrm{see}\ [6]).$

Note that $\displaystyle\lim_{\lambda\rightarrow
0}S_{2,\lambda}(n,k)=S_{2}(n,k)\displaystyle$ and
$\displaystyle\lim_{\lambda\rightarrow
0}S_{1,\lambda}(n,k)=S_{1}(n,k)\displaystyle$.

By (7) and (9), we easily obtain the next proposition.

######  Proposition 1.

The following orthogonality and inverse relations hold true.

$\displaystyle\sum_{k=l}^{n}S_{1,\lambda}(n,k)S_{2,\lambda}(k,l)=\delta_{n,l},\,\,\sum_{k=l}^{n}S_{2,\lambda}(n,k)S_{1,\lambda}(k,l)=\delta_{n,l},\,\,(0\leq
l\leq n),$ $\displaystyle
a_{n}=\sum_{k=0}^{n}S_{2,\lambda}(n,k)b_{k}\Longleftrightarrow
b_{n}=\sum_{k=0}^{n}S_{1,\lambda}(n,k)a_{k},$ $\displaystyle
a_{n}=\sum_{k=n}^{m}S_{2,\lambda}(k,n)b_{k}\Longleftrightarrow
b_{n}=\sum_{k=n}^{m}S_{1,\lambda}(k,n)a_{k},$

where $\delta_{n,l}$ is the Kronecker’s symbol.

Let $a$ and $a^{+}$ be the boson annihilation and creation operators
satisfying the commutation relation

(10) $[a,a^{+}]=aa^{+}-a^{+}a=1,\quad(\mathrm{see}\ [4,5,8]).$

It has been known for some time that the normal ordering of $(a^{+}a)^{n}$ has
a close relation to the Stirling numbers of the second kind ([4,5,8]). Indeed,
the normal ordering of an integral power of the number operator $a^{+}a$ in
terms of boson operators $a$ and $a^{+}$, that satisfy the commutation
$[a,a^{+}]=aa^{+}-a^{+}a=1$, can be written in the form

(11)
$(a^{+}a)^{k}=\sum_{l=0}^{k}S_{2}(k,l)(a^{+})^{l}a^{l},\quad(\mathrm{see}\
[4,5,8]).$

The number states $|m\rangle,\ m=1,2,\dots$, are defined as

(12) $a|m\rangle=\sqrt{m}|m-1\rangle,\quad
a^{+}|m\rangle=\sqrt{m+1}|m+1\rangle.$

By (12), we get $a^{+}a|m\rangle=m|m\rangle$, (see [4,5,8]). The coherent
states $|z\rangle$, where $z$ is a complex number, satisfy
$a|z\rangle=z|z\rangle,$ $z|z\rangle=1$. To show a connection to coherent
states, we recall that the harmonic oscillator has Hamiltonian
$\hat{n}=a^{+}a$ (neglecting the zero point energy) and the usual eigenstates
$|n\rangle$ (for $n\in\mathbb{N}$) satisfying $\hat{n}|n\rangle=n|n\rangle$
and $\langle m|n\rangle=\delta_{m,n}$, where $\delta_{m,n}$ is the Kronecker’s
symbol. In this paper, we derive some identities involving the degenerate
Stirling numbers arising from the normal ordering of degenerate integral
powers of the number operator $a^{+}a$ in terms of boson operators $a$ and
$a^{+}$, which are given by

$(a^{+}a)_{k,\lambda}=\sum_{l=0}^{k}S_{2,\lambda}(k,l)(a^{+})^{l}a^{l}.$

In this paper, we study the degenerate Stirling numbers associated with the
number operator $a^{+},a$ that satisfy $[a,a^{+}]=aa^{+}-a^{+}a=1$. Indeed, we
give some new formulae of the degenerate Stirling numbers which are derived
from an integral power of number operator $a^{+}a$ in terms of boson operators
$a$ and $a^{+}$ with

$(a^{+}a)_{k,\lambda}=\sum_{l=0}^{k}S_{2,\lambda}(k,l)(a^{+})^{l}a^{l}.$

##  2\. Some identities involving degenerate Stirling numbers

We recall that the bosonic commutation relation $[a,a^{+}]=aa^{+}-a^{+}a=1$
can be realized formally in a suitable space of functions by letting
$a=\frac{d}{dx}$ and $a^{+}=x$ (the operator of multiplication by $x$). It is
known that

(13)
$\bigg{(}x\frac{d}{dx}\bigg{)}_{n,\lambda}f(x)=\sum_{k=0}^{n}S_{2,\lambda}(n,k)x^{k}\bigg{(}\frac{d}{dx}\bigg{)}^{k}f(x),\quad(\mathrm{see}\
[7]),$

where $n$ is a nonnegative integer.

The equation (13) can be written as

(14) $(a^{+}a)_{n,\lambda}=\sum_{k=0}^{n}S_{2,\lambda}(n,k)(a^{+})^{k}a^{k}.$

From (9), we note that

(15)
$S_{2,\lambda}(k+1,l)=S_{2,\lambda}(k,l-1)+(l-k\lambda)S_{2,\lambda}(k,l),$

where $k,l\in\mathbb{N}$ with $k\geq l$. It is easy to show that

$[a,\hat{n}]=a,\quad[\hat{n},a^{+}]=a^{+}.$

Inverting (13) by using Proposition 1, we have

(16)
$x^{k}\bigg{(}\frac{d}{dx}\bigg{)}^{k}f(x)=\sum_{m=0}^{k}S_{1,\lambda}(k,m)\bigg{(}x\frac{d}{dx}\bigg{)}_{m,\lambda}f(x).$

In view of (16), the normal ordering of a degenerate integral power of the
number operator $a^{+}a$ in terms of boson operators $a,a^{+}$ can be
rewritten in the form

(17) $\displaystyle(a^{+})^{k}a^{k}$
$\displaystyle=\sum_{m=0}^{k}S_{1,\lambda}(k,m)(a^{+}a)_{m,\lambda}$
$\displaystyle=\sum_{m=0}^{k}S_{1,\lambda}(k,m)(\hat{n})_{m,\lambda}=(\hat{n})_{k},$

where $k$ is a positive integer.

######  Proposition 2.

For $k\in\mathbb{N}$, we have

$\displaystyle(\hat{n})_{k}=(a^{+})^{k}a^{k}$
$\displaystyle=\sum_{m=0}^{k}S_{1,\lambda}(k,m)(a^{+}a)_{m,\lambda}=\sum_{m=0}^{k}S_{1,\lambda}(k,m)(\hat{n})_{m,\lambda}.$

We note that

(18) $\displaystyle a^{+}(\hat{n}+1-\lambda)_{k,\lambda}a$
$\displaystyle=a^{+}\sum_{l=0}^{k}\binom{k}{l}(\hat{n})_{l,\lambda}(1-\lambda)_{k-l,\lambda}a$
$\displaystyle=\sum_{l=0}^{k}\binom{k}{l}(1-\lambda)_{k-l,\lambda}\sum_{m=0}^{l}S_{2,\lambda}(l,m)(a^{+})^{m+1}a^{m+1}$
$\displaystyle=\sum_{m=0}^{k}\bigg{(}\sum_{l=m}^{k}\binom{k}{l}(1-\lambda)_{k-l,\lambda}S_{2,\lambda}(l,m)\bigg{)}(a^{+})^{m+1}a^{m+1}.$

By (10), we also have

(19) $\displaystyle(\hat{n})_{k+1,\lambda}$
$\displaystyle=(\hat{n}-\lambda)_{k,\lambda}\hat{n}=a^{+}\big{(}(\hat{n}+1-\lambda)\cdots(\hat{n}+1-k\lambda)\big{)}a$
$\displaystyle=a^{+}(\hat{n}+1-\lambda)_{k,\lambda}a.$

By (19) and (14), we get

$\displaystyle a^{+}(\hat{n}+1-\lambda)_{k,\lambda}a$
$\displaystyle=(\hat{n})_{k+1,\lambda}=(a^{+}a)_{k+1,\lambda}$ (20)
$\displaystyle=\sum_{m=0}^{k+1}S_{2,\lambda}(k+1,m)(a^{+})^{m}a^{m}$
$\displaystyle=\sum_{m=1}^{k+1}S_{2,\lambda}(k+1,m)(a^{+})^{m}a^{m}$
$\displaystyle=\sum_{m=0}^{k}S_{2,\lambda}(k+1,m+1)(a^{+})^{m+1}a^{m+1}.$

Therefore, by (18) and (20), we obtain the following theorem.

######  Theorem 3.

For $m,k\in\mathbb{Z}$ with $k\geq m\geq 0$, we have

$\sum_{l=m}^{k}\binom{k}{l}(1-\lambda)_{k-l,\lambda}S_{2,\lambda}(l,m)=S_{2,\lambda}(k+1,m+1).$

From Proposition 2, we note that

(21) $\displaystyle(a^{+})^{k+1}a^{k+1}$
$\displaystyle=\sum_{m=0}^{k+1}S_{1,\lambda}(k+1,m)(a^{+}a)_{m,\lambda}$
$\displaystyle=\sum_{m=1}^{k+1}S_{1,\lambda}(k+1,m)(\hat{n})_{m,\lambda}$
$\displaystyle=\sum_{m=0}^{k}S_{1,\lambda}(k+1,m+1)(\hat{n})_{m+1,\lambda},\quad(k\geq
0).$

The degenerate rising factorial sequence is defined by

(22) $\langle x\rangle_{0,\lambda}=1,\quad\langle
x\rangle_{n,\lambda}=x(x+\lambda)\cdots(x+(n-1)\lambda),\quad(n\geq 1).$

Now, from (10) we note that

(23)
$a^{+}(\hat{n})_{k,\lambda}a=(\hat{n}-1)_{k,\lambda}\hat{n}=\hat{n}(\hat{n}-1)_{k,\lambda},\quad(k\in\mathbb{N}).$

On the other hand, by (17) and (23), we get

(24)
$\displaystyle(a^{+})^{k+1}a^{k+1}=\sum_{l=0}^{k}S_{1,\lambda}(k,l)a^{+}(\hat{n})_{l,\lambda}a=\sum_{l=0}^{k}S_{1,\lambda}(k,l)(\hat{n}-1)_{l,\lambda}\hat{n}$
$\displaystyle=\sum_{l=0}^{k}S_{1,\lambda}(k,l)\sum_{m=0}^{l}\binom{l}{m}(-1)^{l-m}\langle
1\rangle_{l-m,\lambda}(\hat{n})_{m,\lambda}\big{(}\hat{n}-m\lambda+m\lambda)$
$\displaystyle=\sum_{l=0}^{k}S_{1,\lambda}(k,l)\sum_{m=0}^{l}\binom{l}{m}(-1)^{l-m}\langle
1\rangle_{l-m,\lambda}(\hat{n})_{m+1,\lambda}$
$\displaystyle\quad\quad+\lambda\sum_{l=0}^{k}S_{1,\lambda}(k,l)\sum_{m=1}^{l+1}\binom{l}{m}(-1)^{l-m}\langle
1\rangle_{l-m,\lambda}(\hat{n})_{m,\lambda}m$
$\displaystyle=\sum_{m=0}^{k}\sum_{l=m}^{k}S_{1,\lambda}(k,l)\binom{l}{m}(-1)^{l-m}\langle
1\rangle_{l-m,\lambda}(\hat{n})_{m+1,\lambda}$
$\displaystyle\quad\quad+\sum_{m=0}^{k}\sum_{l=m}^{k}S_{1,\lambda}(k,l)(m+1)\lambda\binom{l}{m+1}(-1)^{l-m-1}\langle
1\rangle_{l-m-1,\lambda}(\hat{n})_{m+1,\lambda}.$

Therefore, by (21) and (24), we obtain the following theorem.

######  Theorem 4.

For $m,k\in\mathbb{Z}$ with $k\geq m\geq 0$, we have

$\displaystyle S_{1,\lambda}(k+1,m+1)$
$\displaystyle\quad=\sum_{l=m}^{k}S_{1,\lambda}(k,l)\bigg{\\{}\binom{l}{m}\langle
1\rangle_{l-m,\lambda}(-1)^{l-m}+l\lambda\binom{l-1}{m}(-1)^{l-m-1}\langle
1\rangle_{l-m-1,\lambda}\bigg{\\}}.$

From (23), we get

(25) $\displaystyle
a^{+}(\hat{n})_{k,\lambda}a=\hat{n}(\hat{n}-1)_{k,\lambda}=\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}\hat{n}(\hat{n})_{m,\lambda}$
$\displaystyle=\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}(\hat{n}-m\lambda+m\lambda)(\hat{n})_{m,\lambda}$
$\displaystyle=\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}\big{(}(\hat{n})_{m+1,\lambda}+m\lambda(\hat{n})_{m,\lambda}\big{)}$
$\displaystyle=\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}(\hat{n})_{m+1,\lambda}$
$\displaystyle\quad+\sum_{m=0}^{k}\binom{k}{m+1}(-1)^{k-m-1}\langle
1\rangle_{k-m-1,\lambda}(m+1)\lambda(\hat{n})_{m+1,\lambda}$
$\displaystyle=\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}\sum_{p=1}^{m+1}S_{2,\lambda}(m+1,p)(a^{+})^{p}a^{p}$
$\displaystyle\quad+\lambda\sum_{m=0}^{k}\binom{k}{m+1}(-1)^{k-m-1}\langle
1\rangle_{k-m-1,\lambda}(m+1)\sum_{p=1}^{m+1}S_{2,\lambda}(m+1,p)(a^{+})^{p}a^{p}.$
$\displaystyle=\sum_{p=1}^{k+1}\bigg{(}\sum_{m=p-1}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}S_{2,\lambda}(m+1,p)\bigg{)}(a^{+})^{p}a^{p}$
$\displaystyle\quad+\lambda\sum_{p=1}^{k+1}\sum_{m=p-1}^{k}\binom{k}{m+1}(-1)^{k-m-1}\langle
1\rangle_{k-m-1,\lambda}(m+1)S_{2,\lambda}(m+1,p)(a^{+})^{p}a^{p}$
$\displaystyle=\sum_{p=0}^{k}\bigg{\\{}\sum_{m=p}^{k}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}S_{2,\lambda}(m+1,p+1)\bigg{\\}}(a^{+})^{p+1}a^{p+1}$
$\displaystyle\quad+\lambda\sum_{p=0}^{k}\sum_{m=p}^{k}k\binom{k-1}{m}(-1)^{k-m-1}\langle
1\rangle_{k-m-1,\lambda}S_{2,\lambda}(m+1,p+1)(a^{+})^{p+1}a^{p+1}.$

We observe that

(26) $\displaystyle a^{+}(\hat{n})_{k,\lambda}a$
$\displaystyle=a^{+}\bigg{(}\sum_{p=0}^{k}S_{2,\lambda}(k,p)(a^{+})^{p}a^{p}\bigg{)}a$
$\displaystyle=\sum_{p=0}^{k}S_{2,\lambda}(k,p)(a^{+})^{p+1}a^{p+1}.$

Therefore, by (25) and (26), we obtain the following theorem.

######  Theorem 5.

For $p,k\in\mathbb{Z}$ with $0\leq p\leq k$, we have

$\displaystyle S_{2,\lambda}(k,p)=$
$\displaystyle\sum_{m=p}^{k}\bigg{\\{}\binom{k}{m}(-1)^{k-m}\langle
1\rangle_{k-m,\lambda}S_{2,\lambda}(m+1,p+1)$ $\displaystyle\quad+\lambda
k\binom{k-1}{m}(-1)^{k-m-1}\langle
1\rangle_{k-m-1,\lambda}S_{2,\lambda}(m+1,p+1)\bigg{\\}}.$

By (17), we obtain

(27)
$\displaystyle(a^{+})^{k+1}a^{k+1}=\sum_{p=0}^{k+1}S_{1,\lambda}(k+1,p)(a^{+}a)_{p,\lambda}=\sum_{p=1}^{k+1}S_{1,\lambda}(k+1,p)(\hat{n})_{p,\lambda}$
$\displaystyle=\sum_{p=1}^{k+1}S_{1,\lambda}(k+1,p)\hat{n}(\hat{n}-1+1-\lambda)_{p-1,\lambda}$
$\displaystyle=\sum_{p=1}^{k+1}S_{1,\lambda}(k+1,p)\hat{n}\sum_{l=0}^{p-1}\binom{p-1}{l}(1-\lambda)_{p-1-l,\lambda}(\hat{n}-1)_{l,\lambda}$
$\displaystyle=\sum_{l=0}^{k}\bigg{(}\sum_{p=l+1}^{k+1}S_{1,\lambda}(k+1,p)\binom{p-1}{l}(1-\lambda)_{p-1-l,\lambda}\bigg{)}\hat{n}(\hat{n}-1)_{l,\lambda}$
$\displaystyle=\sum_{l=0}^{k}\bigg{(}\sum_{p=l}^{k}S_{1,\lambda}(k+1,p+1)\binom{p}{l}(1-\lambda)_{p-l,\lambda}\bigg{)}\hat{n}(\hat{n}-1)_{l,\lambda}.$

Therefore, from the first line of (24) and (27), we obtain the following
theorem.

######  Theorem 6.

For $k,l\in\mathbb{Z}$ with $k\geq l\geq 0$, we have

$S_{1,\lambda}(k,l)=\sum_{p=l}^{k}S_{1,\lambda}(k+1,p+1)\binom{p}{l}(1-\lambda)_{p-l,\lambda}.$

From (9), we note that

(28)
$\frac{1}{k}\big{(}e_{\lambda}(t)-1\big{)}^{k}=\sum_{n=k}^{\infty}S_{2,\lambda}(n,k)\frac{t^{n}}{n!},\quad(k\geq
0).$

Thus, by (28), we get

(29) $\displaystyle\sum_{n=k}^{\infty}S_{2,\lambda}(n,k)\frac{t^{n}}{n!}$
$\displaystyle=\frac{1}{k!}\big{(}e_{\lambda}(t)-1\big{)}^{k}=\frac{1}{k!}\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}e_{\lambda}^{m}(t)$
$\displaystyle=\sum_{n=0}^{\infty}\bigg{(}\frac{1}{k!}\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}(m)_{n,\lambda}\bigg{)}\frac{t^{n}}{n!}.$

Comparing the coefficients on both sides of (29), we have

$\sum_{m=0}^{k}\binom{k}{m}(-1)^{k-m}(m)_{n,\lambda}=\left\\{\begin{array}[]{ccc}k!S_{2,\lambda}(n,k),&\textrm{if
$n\geq k,$}\\\ 0,&\textrm{if $0\leq n<k$}.\end{array}\right.$

By (15), we get

(30) $\displaystyle
S_{2,\lambda}(k+1,l+1)=S_{2,\lambda}(k,l)+(l+1-k\lambda)S_{2,\lambda}(k,l+1)$
$\displaystyle=S_{2,\lambda}(k,l)+(l+1-k\lambda)$
$\displaystyle\quad\quad\times\big{\\{}S_{2,\lambda}(k-1,l)+(l+1-(k-1)\lambda)S_{2,\lambda}(k-1,l+1)\big{\\}}$
$\displaystyle=S_{2,\lambda}(k,l)+(l+1-k\lambda)S_{2,\lambda}(k-1,l)$
$\displaystyle\quad\quad+\langle
l+1-k\lambda\rangle_{2,\lambda}S_{2,\lambda}(k-1,l+1)$
$\displaystyle=S_{2,\lambda}(k,l)+(l+1-k\lambda)S_{2,\lambda}(k-1,l)+\langle
l+1-k\lambda\rangle_{2,\lambda}S_{2,\lambda}(k-2,l)$
$\displaystyle\quad+\cdots+\langle
l+1-k\lambda\rangle_{k-l,\lambda}S_{2,\lambda}(k,l)$
$\displaystyle=\sum_{m=l}^{k}\langle
l+1-k\lambda\rangle_{k-m,\lambda}S_{2,\lambda}(l,l).$

Therefore, by (30), we obtain the following theorem.

######  Theorem 7.

For $l,k\in\mathbb{Z}$ with $0\leq l\leq k$, we have

$S_{2,\lambda}(k+1,l+1)=\sum_{m=l}^{k}\langle
l+1-k\lambda\rangle_{k-m,\lambda}S_{2,\lambda}(m,l).$

From (8), we note that

(31) $\displaystyle
S_{1,\lambda}(k+1,m+1)=S_{1,\lambda}(k,m)-\big{(}k-(m+1)\lambda\big{)}S_{1,\lambda}(k,m+1)$
$\displaystyle=S_{1,\lambda}(k,m)-\big{(}k-(m+1)\lambda\big{)}$
$\displaystyle\quad\quad\times\Big{(}S_{1,\lambda}(k-1,m)-((k-1)-(m+1)\lambda)S_{1,\lambda}(k-1,m+1)\Big{)}$
$\displaystyle=S_{1,\lambda}(k,m)-(k-(m+1)\lambda)S_{1,\lambda}(k-1,m)$
$\displaystyle\quad\quad+(-1)^{2}(k-(m+1)\lambda)_{2}S_{1,\lambda}(k-1,m+1)$
$\displaystyle=S_{1,\lambda}(k,m)-(k-(m+1)\lambda)S_{1,\lambda}(k-1,m)$
$\displaystyle\quad\quad+(-1)^{2}(k-(m+1)\lambda)_{2}S_{1,\lambda}(k-1,m+1)$
$\displaystyle\quad\quad+\cdots+(-1)^{k-m}(k-(m+1)\lambda)_{k-m}S_{1,\lambda}(m,m)$
$\displaystyle=\sum_{l=m}^{k}(-1)^{k-l}\big{(}k-(m+1)\lambda)_{k-l}S_{1,\lambda}(l,m).$

Therefore, by (31), we obtain the following theorem.

######  Theorem 8.

For $m,k\in\mathbb{Z}$ with $k\geq m\geq 0$, we have

$S_{1,\lambda}(k+1,m+1)=\sum_{l=m}^{k}(-1)^{k-l}\big{(}k-(m+1)\lambda\big{)}_{k-l}S_{1,\lambda}(l,m).$

By (15), we get

(32) $\displaystyle
S_{2,\lambda}(m+k+1,m)=S_{2,\lambda}(m+k,m-1)+\big{(}m-(m+k)\lambda\big{)}S_{2,\lambda}(m+k,m)$
$\displaystyle=S_{2,\lambda}(m+k-1,m-2)+\big{(}(m-1)-(m+k-1)\lambda\big{)}S_{2,\lambda}(m+k-1,m-1)$
$\displaystyle\quad\quad+\big{(}m-(m+k)\lambda\big{)}S_{2,\lambda}(m+k,m)$
$\displaystyle=\big{(}m-(m+k)\lambda\big{)}S_{2,\lambda}(m+k,m)+\big{(}(m-1)-(m+k-1)\lambda\big{)}$
$\displaystyle\quad\quad\times
S_{2,\lambda}(m+k-1,m-1)+\cdots+(0-k\lambda)S_{2,\lambda}(k,0)$
$\displaystyle=\sum_{l=0}^{m}\big{(}l-(k+l)\lambda\big{)}S_{2,\lambda}(k+l,l).$

Therefore, by (32), we obtain the following theorem.

######  Theorem 9.

For $k,m\in\mathbb{Z}$ with $k,m\geq 0$, we have

$S_{2,\lambda}(m+k+1,m)=\sum_{l=0}^{m}\big{(}l-(k+l)\lambda\big{)}S_{2,\lambda}(k+l,l).$

From Theorem 3, we note that

$S_{2,\lambda}(k+1,m+1)=\sum_{l=m}^{k}\binom{k}{l}(1-\lambda)_{k-l,\lambda}S_{2,\lambda}(l,m).$

Now, by using Propositon 1 we have

(33) $\displaystyle\sum_{l=m}^{k}S_{2,\lambda}(k+1,l+1)S_{1,\lambda}(l,m)$
$\displaystyle=\sum_{l=m}^{k}\sum_{p=l}^{k}\binom{k}{p}(1-\lambda)_{k-p,\lambda}S_{2,\lambda}(p,l)S_{1,\lambda}(l,m).$
$\displaystyle=\sum_{p=m}^{k}\binom{k}{p}(1-\lambda)_{k-p,\lambda}\sum_{l=m}^{p}S_{2,\lambda}(p,l)S_{1,\lambda}(l,m)$
$\displaystyle=\binom{k}{m}(1-\lambda)_{k-m,\lambda}.$

Therefore, by (33), we obtain the following theorem.

######  Theorem 10.

For $m,k\in\mathbb{Z}$ with $k\geq m\geq 0$, we have

$\binom{k}{m}=\frac{1}{(1-\lambda)_{k-m,\lambda}}\sum_{l=m}^{k}S_{2,\lambda}(k+1,l+1)S_{1,\lambda}(l,m).$

By (8), we get

(34) $\displaystyle
S_{1,\lambda}(m+k+1,m)=S_{1,\lambda}(m+k,m-1)-(m+k-m\lambda)S_{1,\lambda}(m+k,m)$
$\displaystyle=S_{1,\lambda}(m+k-1,m-2)-(m+k-1-(m-1)\lambda))S_{1,\lambda}(m+k-1,m-1)$
$\displaystyle\quad\quad-(m+k-m\lambda)S_{1,\lambda}(m+k,m)$
$\displaystyle=-(m+k-m\lambda)S_{1,\lambda}(m+k,m)-(m+k-1-(m-1)\lambda)$
$\displaystyle\quad\quad\times S_{1,\lambda}(m+k-1,m-1)-\cdots-
kS_{1,\lambda}(k,0)$
$\displaystyle=-\sum_{l=0}^{m}(k+l-l\lambda)S_{1,\lambda}(k+l,l).$

Therefore, by (34), we obtain the following theorem.

######  Theorem 11.

For $m,k\in\mathbb{Z}$ with $m,k\geq 0$, we have

$S_{1,\lambda}(m+k+1,m)=-\sum_{l=0}^{m}(k+l-l\lambda)S_{1,\lambda}(k+l,l).$

##  3\. Conclusion

In recent years, studying degenerate versions of some special numbers and
polynomials have drawn the attention of many mathematicians with their
regained interests not only in combinatorial and arithmetical properties but
also in applications to differential equations, identities of symmetry and
probability theory. These degenerate versions include the degenerate Stirling
numbers of the first and second kinds, degenerate Bernoulli numbers of the
second kind and degenerate Bell numbers and polynomials.

As a degenerate version of the well known normal ordering of an integral power
of the number operator, we considered the normal ordering of a degenerate
integral power of the number operator in terms of boson operators and its
inversion as well. We derived some identities and recurrence relations for the
degenerate Stirling numbers of the first kind and of the second kind by using
the normal ordering in (1) and its inversion in (2), certain relations of
boson operators and the recurrence relations of the Stirling numbers
themselves.

It is one of our future projects to continue to explore various degenerate
versions of many special polynomials and numbers by using various methods
mentioned in the Introduction.

Conflict of Interest

We have no conflicts of interest to disclose.

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * [1] Araci, S. _A new class of Bernoulli polynomials attached to polyexponential functions and related identities._ Adv. Stud. Contemp. Math. (Kyungshang) 31 (2021), no. 2, 195-204.
  * [2] Carlitz, L. _Degenerate Stirling, Bernoulli and Eulerian numbers._ Utilitas Math. 15 (1979), 51-88.
  * [3] Comtet, L. _Advanced combinatorics._ The art of finite and infinite expansions. Revised and enlarged edition. D. Reidel Publishing Co., Dordrecht, 1974. xi+343 pp.
  * [4] Katriel, J. _Bell numbers and coherent states._ Phys. Lett. A 273 (2000), no. 3, 159-161.
  * [5] Katriel, J. _Coherent states and combinatorics._ PJ. Opt. B.: Quantum Semiclass. Opt. 4 (2002), S200-S203.
  * [6] Kim, D. S.; Kim, T. _A note on a new type of degenerate Bernoulli numbers._ Russ. J. Math. Phys. 27 (2020), no. 2, 227-235.
  * [7] Kim, T.; Kim, D. S. _On some degenerate differential and degenerate difference operator._ Russ. J. Math. Phys. 29 (2022), no. 1, 37-46.
  * [8] Perelomov, A. _Generalized coherent states and their applications._ Texts and Monographs in Physics. SpringerVerlag, Berlin, 1986. xii+320 pp.
  * [9] Roman, S. _The umbral calculus._ Pure and Applied Mathematics, 111. Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York, 1984. x+193 pp.
  * [10] Simsek, Y. _Construction of generalized Leibnitz type numbers and their properties._ Adv. Stud. Contemp. Math. (Kyungshang) 31 (2021), no. 3, 311-323.