# Evolution of waves in liquid films on moving substrates

Tsvetelina Ivanova  von Karman Institute for Fluid Dynamics,
Waterloosesteenweg 72, Sint-Genesius-Rode, Belgium  Fabio Pino  von Karman
Institute for Fluid Dynamics, Waterloosesteenweg 72, Sint-Genesius-Rode,
Belgium  Benoit Scheid  TIPs lab, Université Libre de Bruxelles, Belgium
Miguel A. Mendez<EMAIL_ADDRESS>von Karman Institute for Fluid Dynamics,
Waterloosesteenweg 72, Sint-Genesius-Rode, Belgium

###### Abstract

Accurate and computationally accessible models of liquid film flows allow for
optimizing coating processes such as hot-dip galvanization and vertical slot-
die coating. This paper extends the classic three-dimensional integral
boundary layer (IBL) model for falling liquid films (FF) to account for a
moving substrate (MS). We analyze the stability of the liquid films on
vertically moving substrates in a linear and in a nonlinear setting. In the
linear analysis, we derive the dispersion relation and the temporal growth
rates of an infinitesimal disturbance using normal modes and linearized
governing equations. In the nonlinear analysis, we consider disturbances of
finite size and numerically compute their evolution using the set of nonlinear
equations in which surface tension has been removed. We present the region of
(linear) stability of both FF and MS configurations, and we place the
operating conditions of an industrial galvanizing line in these maps. A wide
range of flow conditions was analyzed and shown to be stable according to
linear and nonlinear stability analyses. Moreover, the nonlinear analysis,
carried out in the absence of surface tension, reveals a nonlinear stabilizing
mechanism for the interface dynamics of a liquid film dragged by an upward-
moving substrate.

††preprint: AIP/123-QED

##  I Introduction

The evolution of waves in liquid films plays a fundamental role in many
coating processes. Their occurrence delimits the range of operating conditions
and influences the quality of final products in the coating industry.

In falling liquid films, waves naturally develop and evolve through various
patterns due to a fascinating interaction between inertial, viscous,
gravitational, and capillary forces[1]. One approach for investigating these
waves is by integral boundary layer (IBL) models (see Shkadov [2]). These
models of lower-dimensionality proved to be useful in simulating industrial
processes such as hot dip galvanization[3], which are not yet accessible by
direct high fidelity simulations because of the prohibitive computational cost
(see Aniszewski  _et al._ [4] and Barreiro-Villaverde, Gosset, and Mendez
[5]). In these models, the dynamics of the liquid film flow is described in
terms of film thickness and streamwise and spanwise flow rates, as opposed to
the Navier-Stokes equations where the film thickness, velocity and pressure
fields must be computed.

The literature on the integral modeling of falling liquid films is vast, with
pioneering two-dimensional formulation proposed by Kapitza [6], Shkadov [2]
and later extended to three-dimensional models by Demekhin and Shkadov [7]
(see also Demekhin and Kalaidin [8], Demekhin  _et al._ [9]). Improvements
over the classic self-similar formalism have been proposed by Ruyer-Quil and
Manneville [10, 11], who solved the inconsistency in the prediction of the
stability threshold by using the method of weighted residuals (see also
Scheid, Ruyer-Quil, and Manneville [12] for a three-dimensional extension). An
extensive review of the modeling of falling liquid film is proposed by Ruyer-
Quil  _et al._ [13] and by Kalliadasis  _et al._ [1].

Integral models enable analytical insight into the flow’s stability ranges and
enable computationally inexpensive simulations of their nonlinear dynamics.
The numerical advantages of using integral models for falling films are also
illustrated in Dietze  _et al._ [14] and Rohlfs, Rietz, and Scheid [15].
Recently, integral models have been extended by Mendez  _et al._ [3] to the
problem of liquid films evolving on a moving substrate in the presence of
pressure gradient and shear stress at the interface. This configuration is
encountered in the jet wiping process in hot-dip galvanization (see Buchlin,
J.M. [16],Gosset, Mendez, and Buchlin [17] Mendez, Gosset, and Buchlin [18]).
Integral models allowed for analyzing the liquid film response to various
disturbances in the process (see also Hocking  _et al._ [19] and Barreiro-
Villaverde, Gosset, and Mendez [5]). Although the configuration was limited to
2-D models, the simulations suggest that thin films are more stable on an
upward-moving substrate than on a fixed one.

This work aims to investigate the reasons for this difference and analyze how
the kinematic condition at the wall influences the stability of the liquid
interface. Moreover, we extend the 2-D models in Mendez  _et al._ [3] to a 3-D
configuration similar to the model by Demekhin and Kalaidin [8], here
generalized to account for the substrate motion. The formulation of the model
includes the presence of shear stress and pressure gradient exerted at the
interface by an external flow, even though these are disregarded in our
analysis. We first use the derived model to perform a classic linear stability
analysis via normal modes and compare the results with the well-known case of
falling liquid films. We then investigate the flow stability numerically in a
nonlinear framework by studying how disturbances of finite size evolve in the
film according to the (nonlinear) set of equations in the absence of surface
tension. Surprisingly, the results show that nonlinearities have a stabilizing
effect since conditions that are linearly unstable are nonlinearly stable.

The rest of the article is structured as follows. Section II reviews the
reference quantities used to scale the dynamics of a liquid film in both the
case of fixed or moving substrate. Section III introduces the integral models
while section IV presents the linear stability analysis for both cases.
Section V describes the numerical methods implemented in an in-house finite
volume solver for the nonlinear partial differential equations (PDEs) derived
in section III. The numerical implementations for a fixed and a moving
substrate are validated in V.1 and V.2, respectively. Section V.3 introduces
the test cases analyzed in this work. The results of the stability analysis
are collected in Section VI for both the linear (VI.1) and the nonlinear
(VI.2) analysis. Conclusions and perspectives are in Section VII.

##  II Scaling Laws

The configuration of interest is three-dimensional and it is illustrated in
Fig. 1. The liquid is assumed to be incompressible with kinematic viscosity
$\nu$, density $\rho$, dynamic viscosity $\mu=\rho\nu$, and surface tension
$\sigma$. As illustrated in the figure, we consider gravity directed towards
$x>0$, we set $y$ orthogonal to the substrate and $z$ in the spanwise
direction. The pressure in the liquid is denoted by $p$, and the three
velocity components are denoted as $\vec{u}=(u,v,w)$, oriented as shown in the
figure. The flow is bounded by the substrate at $y=0$, and the dynamic liquid
interface is at $y=h(x,z,t)$ ($h$ is also referred to as film thickness).

Figure 1: Sketch of the flow domain for a liquid film on a vertically moving
substrate. Gravity is in the opposite direction of the substrate motion.

In the classic falling film problem, herein denoted as ‘FF’, the substrate is
fixed. In the moving substrate problem, here denoted as ‘MS’, it is moving at
a velocity $U_{p}$ against gravity. The FF and MS problems are governed by
different scaling laws.

In the FF problem, one imposes the thickness (or the flow rate) of the film
and sets the scales for the velocity from the viscous-gravity balance in
steady conditions. Denoting as $h_{N}$ the (given) thickness in steady state
conditions, and using square brackets to denote reference quantities such that
$\hat{a}=a/[a]$ is the scaling of $a$ with respect to the reference $[a]$, for
a falling film one has (see Kalliadasis  _et al._ [1]):

$[h]=h_{N},\,\,\,\,\,\,[u]=\frac{g[h]^{2}}{\nu},\,\,\,\,\,\,[q]=[u][h]=\frac{gh^{3}_{N}}{\nu}\,.$
(1)

Accordingly, the Reynolds number in the FF problem is defined as
$Re=[q]/\nu=gh^{3}_{N}/\nu^{2}$.

In the MS problem, one imposes the velocity of the substrate and it is thus
natural to set $[u]=U_{p}$. The scale for the film thickness can also be
computed from the viscous-gravity balance in steady state conditions, leading
to:

$[u]=U_{p},\,\,\,\,\,\,[h]=\sqrt{\frac{\nu
U_{p}}{g}}\,\,\,\,\,\,[q]=[u][h]=\sqrt{\frac{\nu U_{p}^{3}}{g}}\,.$ (2)

Accordingly, the Reynolds number in the MS problem is defined as
$\text{Re}=[q]/\nu=\sqrt{U^{3}_{p}/g\nu}$.

In both problems, it is convenient to scale the streamwise direction $x$ such
that $\varepsilon=[h]/[x]\ll 1$ and the capillary forces
($\sim\sigma[h]/[x]^{3}$) balance the gravitational ones ($\sim\rho g$). This
is known as Shkadov’s scaling, (Shkadov [20], Kalliadasis  _et al._ [1]) and
leads to

$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\Biggl{(}\frac{[h]^{2}\rho
g}{\sigma}\Biggr{)}^{1/3}=\text{We}^{-1/3}\quad\mbox{ (in FF)},$ (3a)
$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\Biggl{(}\frac{[h]^{2}\rho
g}{\sigma}\Biggr{)}^{1/3}=\text{Ca}^{1/3}\qquad\mbox{(in MS)}\,,$ (3b)

having introduced the Weber number $\text{\text{We}}=\sigma/(\rho g[h]^{2})$
for the FF problem, as in Kalliadasis  _et al._ [1], and the capillary number
$\text{Ca}=\mu U_{p}/\sigma$ for the MS problem, as in Mendez  _et al._ [3].
Moreover, it is convenient to introduce a reduced Reynolds number, defined as
$\delta=\varepsilon\text{Re}$, and a dimensionless number that combines Re and
We (or Ca) and depends only on the liquid properties. This is the Kaptiza
number $\text{Ka}=\sigma/(\rho\nu^{4/3}g^{1/3})$, which weighs the importance
of surface tension over viscosity (see also Mendez, Scheid, and Buchlin [21]).
The remaining set of reference quantities is finally $[v]=\varepsilon[u]$ and
$[t]=[x]/[u]$, taking the appropriate references for the FF and the MS
problems.

##  III Integral Boundary Layer Models

In both the FF and the MS problems, the integral boundary layer (IBL) models
can be derived from the Navier-Stokes equations, scaled according to the
reference quantities in Section II and retaining only terms up to
$\mathcal{O}(\varepsilon)$ (see Appendix A for details). This results in the
boundary layer equations:

$\displaystyle\partial_{\hat{x}}\hat{u}+\partial_{\hat{y}}\hat{v}+\partial_{\hat{z}}\hat{w}$
$\displaystyle=$ $\displaystyle 0,$ (4a)
$\displaystyle\delta\big{(}\partial_{\hat{t}}\hat{u}+\hat{u}\partial_{\hat{x}}\hat{u}+\hat{v}\partial_{\hat{y}}\hat{u}+\hat{w}\partial_{\hat{z}}\hat{u}\big{)}$
$\displaystyle=$ $\displaystyle-\partial_{\hat{x}}\hat{p_{x}}$ (4b)
$\displaystyle+\partial^{2}_{\hat{y}\hat{y}}\hat{u}+1,$
$\displaystyle\partial_{\hat{y}}\hat{p_{y}}$ $\displaystyle=$ $\displaystyle
0,$ (4c)
$\displaystyle\delta\big{(}\partial_{\hat{t}}\hat{w}+\hat{u}\partial_{\hat{x}}\hat{w}+\hat{v}\partial_{\hat{y}}\hat{w}+\hat{w}\partial_{\hat{z}}\hat{w}\big{)}$
$\displaystyle=$
$\displaystyle-\partial_{\hat{z}}\hat{p_{z}}+\partial^{2}_{\hat{y}\hat{y}}\hat{w},$
(4d)

where (4a) is the continuity equation, and (4b), (4c), (4d) are the momentum
equations along $\hat{x}$, $\hat{y}$, and $\hat{z}$ respectively. The hat
indicates dimensionless quantities. The dimensionless kinematic boundary
conditions at the wall and at the interface are

$\displaystyle\vec{\hat{v}}\big{|}_{\hat{y}=0}$ $\displaystyle=$
$\displaystyle(\hat{u},\hat{v},\hat{w})\big{|}_{\hat{y}=0}=(\alpha,0,0),$ (5a)
$\displaystyle\hat{v}\big{|}_{\hat{y}=\hat{h}}$ $\displaystyle=$
$\displaystyle\partial_{\hat{t}}\hat{h}+\hat{u}\big{|}_{\hat{y}=\hat{h}}\partial_{\hat{x}}\hat{h}+\hat{w}\big{|}_{\hat{y}=\hat{h}}\partial_{\hat{z}}\hat{h},$
(5b)

where $\alpha=0$ for the FF problem and $\alpha=-1$ for the MS problem. This
parameter is introduced to link the derivation of the two models, but it is
worth stressing that these problems have different scaling laws, as described
in Section II. At $\mathcal{O}(\varepsilon)$, the dynamic boundary conditions
formulating the force balance at the free surface is:

$\displaystyle\hat{p}\big{|}_{\hat{y}=\hat{h}}$ $\displaystyle=$
$\displaystyle\hat{p}_{g}-(\partial_{\hat{x}\hat{x}}\hat{h}+\partial_{\hat{z}\hat{z}}\hat{h}),$
(6a) $\displaystyle\partial_{\hat{y}}\hat{u}\big{|}_{\hat{y}=\hat{h}}$
$\displaystyle=$ $\displaystyle\hat{\tau}_{g,x},$ (6b)
$\displaystyle\partial_{\hat{y}}\hat{w}\big{|}_{\hat{y}=\hat{h}}$
$\displaystyle=$ $\displaystyle\hat{\tau}_{g,z},$ (6c)

where $\hat{p}_{g}$, $\hat{\tau}_{g,x}$ and $\hat{\tau}_{g,z}$ are the gas
pressure and the shear stress components along $x$ and $z$ respectively,
imposed by an external air flow.

To derive the integral model, we integrate (4) along $y$ assuming a self-
similar parabolic velocity profile for both the streamwise $\hat{u}$ and the
spanwise $\hat{w}$ velocity components, as in Demekhin and Shkadov [7]. Using
the local flow rate definitions, the substrate motion and the interface shear
stress, the profiles for the MS case read:

$\displaystyle\hat{u}(\hat{h},\hat{q}_{x},\hat{q}_{z})$ $\displaystyle=$
$\displaystyle\frac{3}{4\hat{h}^{3}}\big{(}\hat{\tau}_{g,x}\hat{h}^{2}-2\hat{h}-2\hat{q}_{x}\big{)}\hat{y}^{2}$
(7a)
$\displaystyle+\frac{6\hat{h}+6\hat{q}_{x}-\hat{\tau}_{g,x}\hat{h}^{2}}{2\hat{h}^{2}}\hat{y}-1,$
$\displaystyle\hat{w}(\hat{h},\hat{q}_{x},\hat{q}_{z})$ $\displaystyle=$
$\displaystyle\frac{3}{4\hat{h}^{3}}\big{(}\hat{\tau}_{g,z}\hat{h}^{2}-2\hat{q}_{z}\big{)}\hat{y}^{2}$
(7b)
$\displaystyle+\frac{6\hat{q}_{z}-\hat{\tau}_{g,z}\hat{h}^{2}}{2\hat{h}^{2}}\hat{y}.$

The main hypothesis under this assumption is that the balance of viscosity and
gravity is not significantly altered by inertia and surface tension.

The integration results in a system of nonlinear partial differential
equations for the liquid film height $\hat{h}$, the streamwise $\hat{q}_{x}$
and spanwise $\hat{q}_{z}$ flow rates. In conservative form, this reads:

$\displaystyle\partial_{\hat{t}}\vec{U}+\nabla\cdot\mathbf{F}=\vec{S},$ (8)

with the state vector $\vec{U}$ consisting of the liquid film height and the
volumetric flow rates, $\vec{U}=(\hat{h},\hat{q}_{x},\hat{q}_{z})^{\text{T}}$.
The source vector is denoted by $\vec{S}=(S_{1},S_{2},S_{3})^{\text{T}}$, and
reads

$\displaystyle\vec{S}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}0\\\
\frac{1}{\delta}\Big{[}\hat{h}\Big{(}-\partial_{\hat{x}}\hat{p}_{x}+\partial_{\hat{x}\hat{x}\hat{x}}\hat{h}+\partial_{\hat{x}\hat{z}\hat{z}}\hat{h}+1\Big{)}+\Delta\hat{\tau}_{x}\Big{]}\\\
\frac{1}{\delta}\Big{[}\hat{h}\Big{(}-\partial_{\hat{z}}\hat{p}_{z}+\partial_{\hat{z}\hat{z}\hat{z}}\hat{h}+\partial_{\hat{z}\hat{x}\hat{x}}\hat{h}\Big{)}+\Delta\hat{\tau}_{z}\Big{]}\end{pmatrix}$
(9)

where the terms with third derivatives of $\hat{h}$ correspond to the
capillary pressure gradients, and the terms
$\Delta\hat{\tau}_{x}=\hat{\tau}_{g,x}+\hat{\tau}_{w,x}$ and
$\Delta\hat{\tau}_{z}=\hat{\tau}_{g,z}+\hat{\tau}_{w,z}$ result from the
integration of the viscous terms in (4b) and (4d). These represent the
difference in shear stress between the interface (terms $\hat{\tau}_{g,x}$ and
$\hat{\tau}_{g,z}$) and the wall (terms $\hat{\tau}_{w,x}$ and
$\hat{\tau}_{w,z}$). The shear stress at the wall, using the self-similar
assumption for the velocity profiles, reads:

$\displaystyle\hat{\tau}_{w,x}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{\tau}_{g,x}-\frac{3\hat{q}_{x}}{\hat{h}^{2}}+\alpha\frac{3}{\hat{h}},$
(10a) $\displaystyle\hat{\tau}_{w,z}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{\tau}_{g,z}-\frac{3\hat{q}_{z}}{\hat{h}^{2}}.$
(10b)

The flux matrix $\mathbf{F}$ in (8) is

$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}F_{11}&F_{12}&F_{13}\\\
F_{21}&F_{22}&F_{23}\end{pmatrix}=\begin{pmatrix}\hat{q}_{x}&\int_{0}^{\hat{h}}\hat{u}^{2}d\hat{y}&\int_{0}^{\hat{h}}\hat{u}\hat{w}d\hat{y}\\\
\hat{q}_{z}&\int_{0}^{\hat{h}}\hat{u}\hat{w}d\hat{y}&\int_{0}^{\hat{h}}\hat{w}^{2}d\hat{y}\end{pmatrix}$
(11)

and has the following components:

$\displaystyle F_{11}$ $\displaystyle=$
$\displaystyle\int_{0}^{\hat{h}}\hat{u}d\hat{y}\eqqcolon\hat{q}_{x},$ (12a)
$\displaystyle F_{21}$ $\displaystyle=$
$\displaystyle\int_{0}^{\hat{h}}\hat{w}d\hat{y}\eqqcolon\hat{q}_{z},$ (12b)
$\displaystyle F_{12}$ $\displaystyle=$
$\displaystyle\frac{1}{120\hat{h}}\Big{(}144\hat{q}_{x}^{2}+6\hat{\tau}_{g,x}\hat{h}^{2}\hat{q}_{x}+\hat{\tau}_{g,x}\hat{h}^{4}$
(12c)
$\displaystyle-\alpha\big{(}48\hat{h}\hat{q}_{x}+6\hat{\tau}_{g,x}\hat{h}^{3}+24\hat{h}^{2}\big{)}\Big{)},$
$\displaystyle F_{22}$ $\displaystyle=$
$\displaystyle\frac{1}{120\hat{h}}\Big{(}144\hat{q}_{x}\hat{q}_{z}+3\hat{\tau}_{g,x}\hat{h}^{2}\hat{q}_{z}+3\hat{\tau}_{g,z}\hat{h}^{2}\hat{q}_{x}$
(12d)
$\displaystyle+\hat{\tau}_{g,x}\hat{\tau}_{g,z}\hat{h}^{4}-\alpha\big{(}24\hat{h}\hat{q}_{z}+3\hat{\tau}_{g,z}\hat{h}^{3}\big{)}\Big{)},$
$\displaystyle F_{13}$ $\displaystyle=$ $\displaystyle F_{22},$ (12e)
$\displaystyle F_{23}$ $\displaystyle=$
$\displaystyle\frac{144\hat{q}_{z}^{2}+6\hat{\tau}_{g,z}\hat{h}^{2}\hat{q}_{z}+\hat{\tau}_{g,z}^{2}\hat{h}^{4}}{120\hat{h}}.$
(12f)

This model recovers the 3-D model for falling liquid films by Demekhin and
Kalaidin [8] if $\alpha=0$,
$\partial_{\hat{x}}\hat{p}_{x}=\partial_{\hat{x}}\hat{p}_{z}=0$,
$\hat{\tau}_{g,x}=\hat{\tau}_{g,z}=0$, and its two-dimensional version by
Shkadov [22] if also $\hat{q}_{z}=0$ and $\partial_{z}\rightarrow 0$.
Moreover, the model recovers the 2-D liquid film model in jet wiping by Mendez
_et al._ [3] if $\alpha=-1$, $\hat{q}_{z}=0$, $\hat{\tau}_{g,z}=0$ and
$\partial_{z}\rightarrow 0$. If $\alpha=-1$ is introduced in system (8), we
obtain the first 3-D formulation of an integral boundary layer model for the
jet wiping process.

##  IV Linear Stability Analysis in 2-D

We consider the 2-D linear stability analysis of the IBL model for FF and MS
conditions, i.e. setting $\hat{q}_{z}=0$ in the system (8) and assuming
$\partial_{\hat{z}}\rightarrow 0$, no shear stress at the interface
($\hat{\tau}_{g,x}=\hat{\tau}_{g,z}=0$) and no pressure gradient
($\partial_{\hat{x}}\hat{p}_{g}=0$). We then introduce

$\hat{h}=\hat{h}_{0}+\tilde{h},\quad\hat{q}_{x}=\hat{q}_{0}+\tilde{q},$ (13)

in the governing Eq. (8), with $\hat{h}_{0},\hat{q}_{0}$ denoting the
thickness and flow rates at an equilibrium solution and
$\tilde{h}<<\hat{h}_{0}$, $\tilde{q}<<\hat{q}_{0}$ some small perturbations.
Linearizing around $\hat{h}_{0},\hat{q}_{0}$ yields the perturbation
equations:

$\displaystyle\partial_{\hat{t}}\tilde{h}+\partial_{x}\tilde{q}$
$\displaystyle=$ $\displaystyle 0,$ (14a)
$\displaystyle\delta\Big{(}\partial_{\hat{t}}\tilde{q}+\partial_{\hat{x}}\tilde{F}\Big{)}$
$\displaystyle=$
$\displaystyle(\hat{h}_{0}+\tilde{h})+(\hat{h}_{0}+\tilde{h})\partial_{\hat{x}\hat{x}\hat{x}}\tilde{h}+\Delta\tilde{\tau},$
(14b)

with

$\displaystyle\tilde{F}$ $\displaystyle=$
$\displaystyle\frac{6(\hat{q}_{0}^{2}+2\hat{q}_{0}\tilde{q})}{5(\hat{h}_{0}+\tilde{h})}-\alpha\Biggl{[}\frac{2}{5}(\hat{q}_{0}+\tilde{q})+\frac{1}{5}(\hat{h}_{0}+\tilde{h})\Biggr{]}$
(15a) $\displaystyle\Delta\tilde{\tau}$ $\displaystyle=$
$\displaystyle-\frac{3(\hat{q}_{0}+\tilde{q})}{\hat{h}_{0}^{2}+2\hat{h}_{0}\tilde{h}}+\frac{3\alpha}{\hat{h}_{0}+\tilde{h}}\,.$
(15b)

We now consider a perturbation in the form of a normal mode, hence:

$\tilde{h}=h_{\varepsilon}\exp\bigl{[}i(\hat{k}\hat{x}-\hat{\omega}\hat{t})\bigr{]},\;\tilde{q}=q_{\varepsilon}\exp\bigl{[}i(\hat{k}\hat{x}-\hat{\omega}\hat{t})\bigr{]},$
(16)

where $\hat{k}$ is the dimensionless wave number and
$\omega=\hat{\omega}_{r}+i\hat{\omega}_{i}$ is the complex dimensionless
angular frequency. Substituting Eq. (16) into Eq. (14b), noticing that the
base state leads to $\hat{q}_{0}=\hat{h}^{3}_{0}/3+\alpha\hat{h}_{0}$ and
separating real and imaginary parts gives an algebraic system of nonlinear
equations:

$\displaystyle\delta\Big{[}2\hat{\omega}_{r}\hat{\omega}_{i}\hat{h}_{0}^{2}-\frac{12}{5}\hat{q}_{0}\hat{h}_{0}\hat{k}\hat{\omega}_{i}+\alpha\frac{2}{5}\hat{h}_{0}^{2}\hat{k}\hat{\omega}_{i}\Big{]}$
$\displaystyle\qquad-3k\hat{h}_{0}^{2}+3\hat{\omega}_{r}-3\alpha\hat{k}=0,$
(17a)
$\displaystyle\delta\Big{[}(\hat{\omega}_{r}^{2}-\hat{\omega}_{i}^{2})\hat{h}_{0}^{2}-\frac{12}{5}\hat{q}_{0}\hat{h}_{0}\hat{k}\hat{\omega}_{r}+\frac{6}{5}k^{2}\hat{q}_{0}^{2}\Big{]}$
$\displaystyle\qquad+\delta\alpha\Big{(}\frac{2}{5}\hat{h}_{0}^{2}\hat{k}\hat{\omega}_{r}+\frac{\hat{k}^{2}\hat{h}_{0}^{2}}{5}\Big{)}$
$\displaystyle\qquad-\hat{k}^{4}\hat{h}_{0}^{3}-3\hat{\omega}_{i}=0.$ (17b)

Since $\hat{q}_{0}$ and $\hat{h}_{0}$ are linked by
$\hat{q}_{0}=\hat{h}^{3}_{0}/3+\alpha\hat{h}_{0}$, for a given $\hat{h}_{0}$
and a pair of $(\hat{k},\delta)$, these equations can be solved for
$\hat{\omega}_{i},\hat{\omega}_{r},\hat{k}$. Moreover, by setting
$\hat{\omega}_{i}=0$, we find the dispersion relation and the neutral curves
(i.e. the loci of conditions in which disturbances neither grow nor decay).

In the FF problem, one has $\alpha=0$ and $\hat{h}_{0}=1,\hat{q}_{0}=1/3$. The
neutral curve is

$\displaystyle\hat{\omega}_{r}=\hat{k}$ $\displaystyle\quad(\text{Dispersion
Relation})$ (18a) $\displaystyle\hat{k}=\sqrt{\frac{\delta}{3}}$
$\displaystyle\quad(\text{Neutral Curve})\,.$ (18b)

In the MS problem, one has $\alpha=-1$ and any $\hat{h}_{0}\in[0,\sqrt{3}]$ is
a possible steady state solution (see Mendez  _et al._ [3]). The neutral curve
is

$\displaystyle\hat{\omega}_{r}=\Big{(}\hat{h}_{0}^{2}-1\Big{)}\hat{k}$
$\displaystyle\quad(\text{Dispersion Relation})$ (19a)
$\displaystyle\hat{k}=\sqrt{\frac{h_{0}^{3}\delta}{3}}$
$\displaystyle\quad(\text{Neutral Curve})\,.$ (19b)

Interestingly, at $\hat{h}_{0}=1$ one has $\hat{\omega}_{r}=0$ for all wave
numbers.

##  V Numerical Methods

We developed an in-house finite volume solver in Python. This is a 3-D
extension of the 2-D solver in Mendez  _et al._ [3]. More specifically, we
blend the two-steps Lax-Wendroff and the two-steps Lax-Friedrichs schemes by
Shampine [23]. The blending is carried out using flux limiters. These allow
for switching between second-order (Lax-Wendroff) and first-order (Lax-
Friedrichs) accuracy depending on the steepness of the solution. More details
on the solver and the discretization schemes can be found in Appendix B.

###  V.1 Validation of the solver for the FF problem

We begin by considering a test case of the FF problem. The test case is from
Doro and Aidun [24], who presented a numerical investigation of falling liquid
films using the Volume of Fluid (VOF) solver in OpenFOAM. In this test case,
the liquid is an aqueous solution of dimethylsulfoxide (DMSO) at Reynolds
number $\text{Re}=15$ and Kapitza number $\text{Ka}=509$. The liquid
properties are $\rho=1098.3$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$,
$\nu=2.85\times 10^{-6}$ ${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ and
$\sigma=0.0484$ $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$. The domain length in
our simulation covers $0.153$ $\mathrm{m}$. In dimensionless units, this
yields $L_{x}=140$, with cell size set to $dx=dz=0.1$. The time step is set to
$dt=0.01$. At the domain’s inlet, perturbations to the flow rate are
introduced at a frequency $f=16$ $\mathrm{Hz}$. Thus, the inlet conditions
are:

$\displaystyle\hat{q}_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\hat{q}_{A}\sin(2\pi\hat{f}\hat{t})+\frac{1}{3},$
(20) $\displaystyle\hat{h}$ $\displaystyle=$
$\displaystyle(3\hat{q}_{x})^{1/3},$

where $\hat{q}_{A}=0.05$ is the perturbation amplitude, $\hat{f}=12[t]=0.048$
is the dimensionless frequency, and $\hat{t}=ndt$ is the dimensionless time
stepping of the simulation.

The disturbances grow over the domain until they produce the classic wave
train observed in forced flows. Sufficiently far from the inlet, the waves are
developed and their shape and phase is nearly invariant to the streamwise
location, as discussed in Doro and Aidun [24]. For this region, the comparison
between the IBL simulations and the VOF simulations (from Fig. 5 in Doro and
Aidun [24]) is shown in Fig. 2.

We observe a good agreement between the two simulations, although it is over-
simplifying to assume self-similarity and $\mathcal{O}(\varepsilon)$ accuracy
in the modeling of a problem that has $\varepsilon=0.29$. Nevertheless,
considering that the computational cost of a (1D) IBL simulation is several
orders of magnitude lower than the computational cost of a (2-D) VOF
simulation, the result is particularly encouraging.

Figure 2: Simulation with the IBL solver of a DMSO falling film to validate it
with results from Doro and Aidun [24] for the FF problem.

###  V.2 Validation of the solver for the MS problem

For the MS problem, we consider the same validation test case as in Mendez
_et al._ [3] to validate our solver. This test case consists of a 2-D wave
train over a moving substrate and was also simulated using high-fidelity VOF
simulations in OpenFOAM. The liquid is water with a Reynolds number $Re=319$
and the substrate moves at $U_{p}=1$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.

The computational domain is rectangular, with a dimensionless length
$L_{x}=8400h_{0}$ in the streamwise direction, and $L_{z}=7.8h_{0}$ in the
spanwise direction, small enough to keep 2-D waves, i.e. stable with respect
to spanwise perturbations. A perturbation with a dimensionless frequency
$\hat{f}=0.05$ is introduced at the inlet’s streamwise flow rate:

$\displaystyle\hat{q}_{x}=\big{[}\frac{1}{3}\hat{h}_{0}^{3}-\hat{h}_{0}\big{]}\big{[}1+\hat{q}_{A}\sin(2\pi\hat{f}\hat{t})\big{]},$
(21)

where $\hat{q}_{A}$ is the perturbation amplitude.

This test case was also used to perform a mesh sensitivity analysis of our
solver. This analysis was based on solutions we obtained with the Lax-
Friedrichs scheme for three different cell sizes $d\hat{x}$, namely
$d\hat{x}=0.0138$, $d\hat{x}=0.0275$ and $d\hat{x}=0.0550$. The results for
the thickness evolution in these three cases are shown in Fig. 3a and compared
to the results obtained by the VOF simulations in OpenFOAM (in which
$d\hat{x}=0.0275$). Fig. 3b further compares the obtained solution for
$d\hat{x}=0.0275$ with the OpenFOAM validation case from Mendez  _et al._ [3].
While numerical dissipation is visible and it is larger for coarser meshes as
expected, its impact can be considered minor within the investigated domain.

(a) Comparison of the solutions obtained with the Lax-Friedrichs scheme for
different cell sizes.

(b) Validation of the current solver (with its solution in a dash-dotted line)
against the OpenFOAM case, which is also used in Mendez  _et al._ [3].

Figure 3: Validation of the solver for the MS case. The inlet is at
$\hat{x}=80$ and the substrate moves against the direction of gravity, which
is $\hat{x}\to\infty$. MAPE is the mean absolute percentage error between the
solution with Lax-Friedrichs (dot-dashed line) and the OpenFOAM data.

###  V.3 Investigated test cases

We analyzed the propagation of nonlinear waves in the MS problem keeping the
same configuration used for the solver validation in Section V.2, hence
introducing perturbations at the inlet flow rate as in Eq. (21). We consider
three dimensionless thicknesses $\hat{h}_{0}=0.1$, $0.2$, $0.3$, two reduced
Reynolds numbers $\delta_{1}=76$ and $\delta_{2}=153$, and fourteen
dimensionless frequencies in the range $\hat{f}=[0.005,0.2]$. This leads to 84
simulations. The simulations were carried out without surface tension to focus
on the impact of nonlinearities in the interface instability. Without surface
tension linear stability theory predicts that all configurations are unstable,
as shown in Section VI.1. However, this was found not to be the case if
nonlinearities are considered.

The dimensionless conditions are representative for hot dip galvanizing lines
as well as the laboratory model ESSOR at the von Karman institute (see
Buchlin, J.M. [16], Gosset, Mendez, and Buchlin [17], Mendez  _et al._ [3]).
The similarity between water and molten zinc in the Skhadov-like scaling used
in this work was discussed in Mendez  _et al._ [3]. For a plate moving at
$U_{p}=1$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, taking $\rho=1000$
$\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$, $\nu=1$
${\mathrm{mm}}^{2}\text{\,}{\mathrm{s}}^{-1}$, $\sigma=0.074$
$\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ for water, and $\rho=6500$
$\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$, $\nu=0.46$
${\mathrm{mm}}^{2}\text{\,}{\mathrm{s}}^{-1}$, $\sigma=0.78$
$\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ for molten zinc leads to $\delta\approx
76$ for both fluids.

In all test cases, waves propagate in the direction of the strip motion
($\hat{x}\rightarrow-\infty$, cf. Fig. 1). Therefore, the domain was set as
$\hat{x}\in[0,40]$, with initial disturbance placed at $\hat{x}=40$. Both the
size of the domain and the mesh size vary from test case to test case.
Specifically, considering that waves propagate at about
$\hat{u}_{\text{w}}\approx 1$ and their wavelength is of the order of
$\hat{\lambda}\approx 1/\hat{f}$, the domain length in the streamwise
direction is taken as $L_{x}=8\lambda$ while the width is taken as
$L_{z}\approx L_{x}/10$. The grid spacing is taken as $dx=\lambda/363$ and
$dz=1/100$, since $363$ points per wavelength proved to give a good compromise
between accuracy and computational cost. The time step is taken such that the
Courant–Friedrichs–Lewy (CFL) number is 0.3, using $|u_{w}|=1$ as an estimate
of the wave velocity. This yields a numerical viscosity
$\mu_{n}=d\hat{x}^{2}/d\hat{t}\propto d\hat{x}$.

Finally, we analyzed the spanwise propagation of three-dimensional
disturbances by considering an inlet flow rate $\hat{q}_{x}$ consisting of an
harmonic term modulated by a Gaussian function $G$ along $\hat{z}$

$\hat{q}_{x}=\big{[}\frac{1}{3}\hat{h}_{0}^{3}-\hat{h}_{0}\big{]}\Big{[}1+\hat{q}_{A}\sin(2\pi\hat{f}\hat{t})\sin\Big{(}\frac{2\pi}{\lambda_{z}}\hat{z}\Big{)}\Big{]}G(\hat{z}),$
(22)

where $\hat{\lambda}_{z}=1$ and the Gaussian modulation with a standard
deviation $\sigma=0.4$ was taken as

$G(\hat{z})=\frac{1}{\sigma\sqrt{2\pi}}\exp\Big{(}\frac{-(z-z_{\text{mean}})^{2}}{2\sigma^{2}}\Big{)}.$
(23)

For the inlet flow rate along $z$, we consider $\hat{q}_{z}=0$. This yields
$\hat{q}_{z}=0$ everywhere and at all times.

##  VI Results and Discussions

###  VI.1 Linear stability

We first consider the result from a linear stability analysis following
equations (17) for both FF and MS. Given a computational grid of 1000
dimensionless wave numbers $\hat{k}$ and 1000 reduced Reynolds numbers
$\delta$, Eq. (17) were solved for $\hat{\omega}_{r}$ and $\hat{\omega}_{i}$
for each pair $(\hat{k},\delta)$. Because multiple solutions exist, we only
focus on the largest $\hat{\omega}_{i}$ (regardless of the sign) and plot a
contour of $\log(\hat{\omega}_{i})$ on the grid $(\hat{k},\delta)$.

Figure 4: Color map of $\log(\hat{\omega}_{i})$ with the natural stability
curve (black thick line) for the FF case.

For the FF problem one has $\alpha=0$, $\hat{h}_{0}=1$ and $\hat{q}_{0}=1/3$.
The associated amplification factors are shown in Fig. 4. The region in white
in the $(k,\delta)$-plane is the region where $\hat{\omega}_{i}<0$, hence
where the film is stable according to the linear stability analysis. The line
separating the stable and unstable regions is the neutral stability curve. For
$k<k_{c}$, the amplification factors show a non-monotonic dependence on
$\delta$: while the region of maximum amplification occurs in $\delta\approx
16$, increasing the Reynolds number leads to reduction of $\omega_{i}$ for all
wavelengths. This trend is due to the scaling of the problem, since
$[\omega]=1/[t]\propto U_{p}^{5/6}\propto\delta^{5/11}$. For later reference,
at $\delta=76$ and $\delta=153$ the critical wave-numbers are $k_{c}\approx 5$
and $k_{c}\approx 7.1$ respectively.

(a) Case for $\hat{h}_{0}=0.1$.

(b) Case for $\hat{h}_{0}=0.2$.

(c) Case for $\hat{h}_{0}=0.3$.

Figure 5: Color map of $\log(\hat{\omega}_{i})$ with the natural stability
curve (black thick lines) for the MS case for $\hat{h}_{0}=0.1,0.2,0.3$. The
conditions with $\delta_{1}=76$ and $\delta_{2}=153$ analyzed in Sec. VI.2 are
located in the white zone further above these stability curves.  Figure 6:
Plot of the imaginary part in the dispersion relation
$\hat{\omega}_{i}(\hat{k})$ with and without surface tension for the case with
$\delta=76$ and $\hat{h}_{0}=0.2$.

The same plots are then produced for the MS problem, and shown in Fig. 5 for
$\hat{h}_{0}=0.1,0.2,0.3$. We recall that in the MS problem one has
$\alpha=-1$, hence $\hat{q}_{0}=\hat{h}_{0}^{3}/3-\hat{h}_{0}$. The location
of the conditions analyzed in the nonlinear setting in section VI.2 are not
show since these are far away from the neutral curve. We have $k$ in the range
$15-17$ for the case $\delta=76$ and in the range $21-23$ for the case
$\delta=153$. The critical wave-number at these Reynolds number (see eq. 19b)
is one order of magnitude lower than the one in the investigated conditions.
Therefore, according to the linear stability analysis, all the investigated
points should be in stable conditions.

The region of maximum amplification is located at much larger $\delta$
compared to the FF and the region of largest amplification moves towards lower
$\delta$.

###  VI.2 Nonlinear analysis of 2-D waves

We here move to the nonlinear analysis of the test cases introduced in Section
V.3. The perturbations are not infinitesimally small and the governing
equations are not linearized. We recall that in this numerical investigation
we do not include the contribution of the surface tension. Therefore, the
dynamics of the liquid film is lacking the stabilizing effect: while all the
investigated test cases are linearly stable if surface tension is included
(cf. Fig. 5), these are (linearly) unstable in absence of surface tension. We
further illustrate this in Fig. 6, which shows the imaginary part of the
dispersion relation $\hat{\omega}_{i}(\hat{k})$ with and without surface
tension for the case with $\delta=76$ and $\hat{h}_{0}=0.2$. We recall, from
Eq. (19a), that the dispersion relation is linear regardless of the surface
tension and hence waves are non-dispersive. For the illustrated case,
$\partial_{\hat{k}}\hat{\omega}_{r}=-0.96$, i.e. waves move approximately at
the substrate speed.

(a) Evolution of nonlinear waves at $\hat{h}_{0}=0.2$ and $\delta=76$.

(b) Evolution of nonlinear waves at $\hat{h}_{0}=0.2$ and $\delta=153$.

Figure 7: Two-dimensional waves evolving without surface tension on a liquid
film on a moving substrate with $\hat{h}_{0}=0.2$. The perturbation frequency
is $\hat{f}=0.2$ ($\hat{\omega}_{r}=1.26$), with amplitude $\hat{q}_{A}=0.2$.
The substrate motion is from left to right; gravity is from right to left. The
inlet is at $\hat{x}=40$.

Moving to the nonlinear analysis in absence of surface tension, Fig. 7 shows
an instantaneous of the film thickness profile, together with the maxima and
minima envelopes computed in the time for each position and considering
$\hat{h}_{0}=0.2$ and $\delta_{1}=76$ (top) and $\delta_{2}=153$ (bottom).
Despite the large disturbance ($\hat{q}_{A}=0.2$ in (21)) at the inlet
($\hat{x}=40$) and despite the initial growth of the waves, the flow remains
convectively stable: the wave amplitude decreases in the stream-wise
direction. These waves are nonlinear, with a steep front and a long tail
similar to the waves observed in the FF problem (Fig. 2) in the first portion
of the domain. However, the interplay of inertia, gravity and viscosity is
different in the MS and FF cases. Contrary to the FF problem, the relative
velocity between the flow and the wall _decreases_ at larger thicknesses if
$\hat{h}_{0}\ll 1$. This implies that, contrary to the FF problem, the crest
of the wave is _slower_ than the substrate film on which they travel.
Therefore, waves tend to level out as they flow and this tendency is more
pronounced for larger waves.

The initial growth produced in $30<\hat{x}<40$ is most probably due to the
mechanism through which the perturbations are injected (see Section V.2). The
boundary condition for injecting the perturbation simulate a manifold and the
flow needs a certain distance to adjust back to the equations governing the
film thickness and flow rate (see Section III). Once this occurs, a clear
decay of the wave amplitude is observed in all investigated configurations. We
can thus fit an exponential decay $h_{M}e^{-\beta\hat{x}}$ to the maximum
thickness temporal envelope and extract a spatial decay rate $\beta$ for each
of the $84$ investigated simulations. The results are shown in Fig. 8 for
$\delta_{1}=76$ (top) and $\delta_{2}=153$ (bottom) over the range of
perturbation frequencies ($\hat{f}\in[0.005,0.2]$) and for three thicknesses
($\hat{h}_{0}=0.1,0.2,0.3$).

(a) Decay amplitudes for a reduced Reynolds number $\delta_{1}=76$.

(b) Decay amplitudes for a reduced Reynolds number $\delta_{2}=153$.

Figure 8: Dependency of the water wave amplitude decay rates on the
perturbation frequencies for three different initial heights $h_{0}$ and two
different reduced Reynolds numbers $\delta_{1}=76$ and $\delta_{2}=153$.

While one would expect the decay rate of the waves to become stronger at
larger frequencies, somewhat less expected is the fact that thicker films
(within the investigated cases) lead to more substantial damping. In line with
the damping mechanism previously described (shown in Fig. 7), thicker films
result in larger waves at a short distance from the inlet, and larger waves
are slower and thus characterized by stronger damping. This mechanism is
promoted by the nonlinearities of the problem and is not captured in the
linear stability framework (which predicts unstable waves in the absence of
surface tension).

(a) spatiotemporal maps of the evolution of 2-D waves in a liquid film over a
moving substrate.

(b) Neither wave merging, nor frequency cross-talk is observed.

Figure 9: The operating conditions in these plots are the same as for the
water wave in Fig. 7 (top) with $\hat{h}_{0}=0.2,\delta=76,\hat{f}=0.2$.

Finally, we highlight that during the downstream propagation, the phase
velocity is approximately constant and equal to $\approx-1$, i.e., substrate
velocity. This can be seen from Fig. 9a, which shows the spatiotemporal
evolution of the waves for a case with $\hat{h}_{0}=0.2$, $\delta=76$ and
$\hat{f}=0.3$. The characteristic lines are straight, showing no appreciable
acceleration or deceleration during their evolution, and their velocity is
line with the linear stability analysis. Fig. 9b shows the evolution of the
frequency content of the thickness evolution at various locations. The
nonlinearities do not produce wave-merging mechanisms within the investigated
test cases: the same frequency is propagated downstream, and higher harmonics
(linked to the initial stiff front of the waves) gently vanishes as the wave
amplitude decreases and their profile becomes more harmonic.

###  VI.3 Analysis of 3-D waves

We conclude this investigation with the analysis of a three-dimensional test
case, with inlet flow rate $\hat{q}_{x}$ prescribed as in (22), and inlet flow
rate $\hat{q}_{z}=0$. The main interest was to analyze if and how three-
dimensional perturbations grow in the spanwise direction $\hat{z}$ (see also
Fig. 1). A snapshot of the liquid film surface for this test case is shown in
Fig. 10a, while Fig. 10b shows the film thickness contour plot.

The results show that the motion of the substrate dominates the direction of
the propagation and the region of influence of the disturbance is particularly
narrow in the $\hat{z}$ direction. This is in contrast to what happens, for
example, in the waves of shallow or deep waters where a perturbation in
relative motion with respect to the substrate (e.g. a ship) produces a
V-shaped wake envelope in case of shallow (non-dispersive) waves and the well-
known Kelvin wedge in the case of deep (dispersive) waves Fitzpatrick [25],
Lighthill [26]. Although we leave the analysis of more general kinds of 3-D
disturbances to future works, these results highlight the distinctive role of
nonlinearities in the MS problem compared to other cases of nonlinear waves in
fluid dynamics.

(a) Dimensionless liquid film height of a 3-D wave generated by perturbations.

(b) Contour plot of the height of the 3-D wave.

Figure 10: Dimensionless liquid film height (top) and its contour plot
(bottom) for 3-D waves generated by the flow rate perturbation of
$\hat{q}_{x}$ in Eq. (22). The substrate moves toward
$\hat{x}\rightarrow-\infty$ while gravity is directed toward
$\hat{x}\rightarrow\infty$ (see Fig. 1). The disturbances are introduced on
the left, at $\hat{x}=160$. The dimensionless perturbation frequency is
$\hat{f}=0.05$, with initial film thickness $\hat{h}_{0}=0.2$ and $\delta=76$.

##  VII Conclusions

We presented integral boundary layer models to describe the evolution of
interface waves in liquid films dragged along upward-moving walls (MS problem)
and compared the problem with the well-known case of liquid films falling
along fixed walls (FF problem). We introduced a dimensionless 3-D integral
boundary layer model that extends classic models for falling liquid films to
account for substrate motion, external pressure gradient and interface shear
stress. These extensions allow for modeling the jet wiping process in hot-dip
galvanization.

The interest in integral models is twofold: they allow for performing
computationally inexpensive numerical simulations of liquid film flows and
enable analytical insights on their stability. We analyzed the stability of
the MS problem in a linear and nonlinear setting. For the linear setting,
focusing on 2-D disturbances, we derived dispersion relations and neutral
curves, and compared the results with those of the FF problem. For the
nonlinear setting, we analyzed the flow response to 2-D disturbances of
various frequencies, Reynolds number and baseline film thickness numerically.
These conditions are relevant to industrial applications.

In the linear stability analysis, the dispersion relation shows the usual
stabilizing effect of surface tension. It was shown that the critical wave
number is smaller than in the FF problem if $\hat{h}_{0}<1$. In the nonlinear
stability analysis, a nonlinear stabilizing mechanism was identified even in
the absence of surface tension. The wave amplitude decays because the
interface velocity is lower than the substrate, and this leveling effect is
larger for thicker films. Finally, we presented a numerical test case with 3-D
disturbances. This showed no disturbance growth in the span-wise direction.

Future work will extend the current analysis to the evolution of other kind of
3-D disturbances, other integral models (e.g., the Weighted IBL formulation)
and the Orr-Sommerfeld problem from the full Navier-Stokes equations.

###### Acknowledgements.

This work is supported by ArcelorMittal in the framework of the ’Ondule VII’
project. Fabio Pino is supported by a FRIA grant from FNRS. Benoit Scheid
thanks the F.R.S.-FNRS for financial support. The solver and post-processing
scripts have been developed using free and open-source software: we thank the
community for building powerful tools from which everyone can benefit.
Finally, we thank David Barreiro from University of A Coruña for providing
valuable feedback and correcting a mistake in the derivation of Eqs. (17a) -
(17b).

##  Appendix A The full (dimensional) problem

The goal is to obtain an integral model for the 3-D liquid film on a moving
substrate, since this approach reduces the number of independent variables and
the dimension of the domain. This leads to a significantly lower computational
cost when performing numerical analysis as shown in Mendez  _et al._ [3]. The
developed finite volume solver for the purposes of this research is a 3-D
extension of the one in Mendez  _et al._ [3].

The starting point for the formulation of the integral model is the Navier-
Stokes equations for a divergence-free Newtonian liquid (Graebel [27]):

$\displaystyle\frac{\partial\vec{v}}{\partial
t}+\vec{v}\cdot\nabla\vec{v}=-\frac{1}{\rho}\nabla
p+\nu\nabla^{2}\vec{v}+\vec{f_{v}},$ (24)

where $\rho$ is the density of the liquid, $\nu$ is its kinematic viscosity,
$\vec{v}=(u,v,w)$ is the velocity field, and $p$ is the pressure field.

The continuity equation for a divergence-free fluid in a Cartesian coordinate
system states

$\displaystyle\nabla\cdot\vec{v}=0.$ (25)

The kinematic boundary conditions at the wall consist of the non-slip and non-
permeability conditions for this problem:

$\displaystyle\vec{v}\big{|}_{y=0}=(-U_{p},0,0),$ (26)

where $U_{p}$ is the speed of the substrate. The continuity of the interface
$h(x,z,t)$ is ensured by the kinematic boundary condition

$\displaystyle\frac{\partial}{\partial
t}(y-h(x,z,t))+\vec{v}\cdot\nabla(y-h(x,z,t))=0$

at $y=h(x,z,t)$ rewritten as

$\displaystyle v=\partial_{t}h+u\partial_{x}h+w\partial_{z}h.$ (27)

The dynamic boundary conditions at the interface represents the force balance
in a local Cartesian reference frame along the free surface normal $\vec{n}$,
and the streamwise and spanwise tangential directions, denoted by
$\vec{t}_{x}$ and $\vec{t}_{z}$, respectively. These vectors have the
following components:

$\displaystyle\vec{n}$ $\displaystyle=$
$\displaystyle\frac{1}{|\vec{n}|}\mbox{$\nabla$}(y-h(x,z,t))$ (28a)
$\displaystyle=$
$\displaystyle\frac{1}{|\vec{n}|}(-\partial_{x}h(x,z,t),1,-\partial_{z}h(x,z,t)),$
(28b) $\displaystyle\vec{t}_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{|\vec{t}_{x}|}(1,\partial_{x}h(x,z,t),0),$ (28c)
$\displaystyle\vec{t}_{z}$ $\displaystyle=$
$\displaystyle\frac{1}{|\vec{t}_{z}|}(0,\partial_{z}h(x,z,t),1).$ (28d)

To denote quantities referring to the gas, the subscript $g$ is introduced to
the variables notations. At the interface $y=h(x,z,t)$, the force balances
along these directions are:

$\displaystyle\vec{n}\cdot\hat{T}\cdot\vec{n}-\vec{n}\cdot\hat{T}_{g}\cdot\vec{n}$
$\displaystyle=$ $\displaystyle\sigma\mbox{$\nabla\cdot\,$}\vec{n},$ (29a)
$\displaystyle\vec{n}\cdot\hat{T}\cdot\vec{t}_{x}-\vec{n}\cdot\hat{T}_{g}\cdot\vec{t}_{x}$
$\displaystyle=$ $\displaystyle 0,$ (29b)
$\displaystyle\vec{n}\cdot\hat{T}\cdot\vec{t}_{z}-\vec{n}\cdot\hat{T}_{g}\cdot\vec{t}_{z}$
$\displaystyle=$ $\displaystyle 0.$ (29c)

The stress tensor for the two fluids is:

$\displaystyle\hat{T}=-p\hat{I}+2\mu\hat{E},$ (30)

where $\hat{I}$ is the identity tensor, and $\hat{E}$ is the strain of rate
tensor that is

$\displaystyle\hat{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\mbox{$\nabla$}\vec{v}+\mbox{$\nabla$}\vec{v}^{T}).$
(31)

It is worth remarking that $\vec{n}\cdot(-p\hat{I})\cdot\vec{n}=-p$ and
$\vec{n}\cdot(-p\hat{I})\cdot\vec{\tau}_{x,z}=0$ both for the liquid and the
gas. Substituting the stress tensor definition (30) into equations (29)
yields:

$\displaystyle-p+\vec{n}\cdot(2\mu\hat{E})\cdot\vec{n}+p_{g}-\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{n}$
$\displaystyle=$ $\displaystyle\sigma\mbox{$\nabla\cdot\,$}\vec{n},$
$\displaystyle\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{x}-\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{t}_{x}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{z}-\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{t}_{z}$
$\displaystyle=$ $\displaystyle 0.$

By introducing the following notations

$\displaystyle p_{g}-\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{n}$
$\displaystyle=$ $\displaystyle p_{g}(x,z,t),$
$\displaystyle\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{t}_{x}$ $\displaystyle=$
$\displaystyle\tau_{g,x}(x,z,t),$
$\displaystyle\vec{n}\cdot(2\mu\hat{E}_{g})\cdot\vec{t}_{z}$ $\displaystyle=$
$\displaystyle\tau_{g,z}(x,z,t),$

we reach the three scalar equations representing the projections of the force
balance along the normal and tangential vectors:

$\displaystyle-p+\vec{n}\cdot(2\mu\hat{E})\cdot\vec{n}+p_{g}(x,z,t)$
$\displaystyle=$ $\displaystyle\sigma\mbox{$\nabla\cdot\,$}\vec{n},$ (32a)
$\displaystyle\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{x}-\tau_{g,x}(x,z,t)$
$\displaystyle=$ $\displaystyle 0,$ (32b)
$\displaystyle\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{z}-\tau_{g,z}(x,z,t)$
$\displaystyle=$ $\displaystyle 0.$ (32c)

The computation of the expressions
$\vec{n}\cdot(2\mu\hat{E})\cdot\vec{n},\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{x}$
and $\vec{n}\cdot(2\mu\hat{E})\cdot\vec{t}_{z}$ in a 3-D Cartesian coordinate
system is a lengthy procedure and its details are therefore omitted in this
paper. It results in the dynamic boundary conditions that represent the force
balance at the interface with projections along $\vec{n},\vec{t}_{x}$, and
$\vec{t}_{z}$, which is shown in Section III in a first-order boundary layer
approximation followed by the integration along the wall-normal $y$-axis.

The 3-D liquid film on a moving substrate is therefore represented by the
dimensional equations (24) and (25) for an incompressible Newtonian fluid with
boundary conditions (26), (27), and (32). The scaling procedure of these
equations and boundary conditions with appropriate reference quantities is
described in Section II.

##  Appendix B Solver and numerical schemes

The numerical solution of the hyperbolic problem (8), a blending between the
Lax-Wendroff and the Lax-Friedrichs schemes is applied. This blending is
achieved with flux limiter functions. They are used to help preventing
oscillations near sharp changes in the solution.

Let the discretized state vector be denoted by $\vec{U}$. In order to evaluate
the state at the next time step, estimations of the quantities at mid-points
in space and time are necessary:

$\displaystyle\vec{U}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}=\frac{1}{2}\big{(}\vec{U}_{i,j}^{n}+\vec{U}_{i+1,j}^{n}\big{)}-\frac{\Delta
t}{2\Delta
x}\big{(}\vec{\mathcal{F}_{x,}}_{i+1,j}^{n}-\vec{\mathcal{F}_{x,}}_{i,j}^{n}\big{)},$
(33a)
$\displaystyle\vec{U}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}=\frac{1}{2}\big{(}\vec{U}_{i,j}^{n}+\vec{U}_{i,j+1}^{n}\big{)}-\frac{\Delta
t}{2\Delta
z}\big{(}\vec{\mathcal{F}_{z,}}_{i,j+1}^{n}-\vec{\mathcal{F}_{z,}}_{i,j}^{n}\big{)}.$
(33b)

The fluxes can be represented by low and high-resolution schemes and a flux
limiter can switch between these schemes depending on the gradients of the
solutions. Values of the fluxes at the half steps are evaluated in the
following way. The high-resolution fluxes are:

$\displaystyle\vec{\mathcal{F}_{x,}}_{i+\frac{1}{2},j}^{\text{high}}$
$\displaystyle=$
$\displaystyle\vec{\mathcal{F}_{x}}(\vec{U}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}),$
(34a) $\displaystyle\vec{\mathcal{F}_{z,}}_{i,j+\frac{1}{2}}^{\text{high}}$
$\displaystyle=$
$\displaystyle\vec{\mathcal{F}_{z}}(\vec{U}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}).$
(34b)

The low-resolution fluxes are:

$\displaystyle\vec{\mathcal{F}_{x,}}_{i+\frac{1}{2},j}^{\text{low}}=\vec{\mathcal{F}_{x,}}(\vec{U}_{i+1,j}^{n})+\frac{\Delta
t}{2\Delta
x}\Big{(}\vec{U}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}-\vec{U}_{i+1,j}^{n}\Big{)},$
(35a)
$\displaystyle\vec{\mathcal{F}_{z,}}_{i,j+\frac{1}{2}}^{\text{low}}=\vec{\mathcal{F}_{z,}}(\vec{U}_{i,j+1}^{n})+\frac{\Delta
t}{2\Delta
z}\Big{(}\vec{U}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-\vec{U}_{i,j+1}^{n}\Big{)}.$
(35b)

The blended fluxes are:

$\displaystyle\vec{F_{x,}}_{i+\frac{1}{2},j}=\vec{\mathcal{F}_{x,}}_{i+\frac{1}{2},j}^{\text{low}}-\phi_{x}(r_{x,i})\Big{(}\vec{\mathcal{F}_{x,}}_{i+\frac{1}{2},j}^{\text{low}}-\vec{\mathcal{F}_{x,}}_{i+\frac{1}{2},j}^{\text{high}}\Big{)}$
(36a)
$\displaystyle\vec{F_{x,}}_{i-\frac{1}{2},j}=\vec{\mathcal{F}_{x,}}_{i-\frac{1}{2},j}^{\text{low}}-\phi_{x}(r_{x,i-1})\Big{(}\vec{\mathcal{F}_{x,}}_{i-\frac{1}{2},j}^{\text{low}}-\vec{\mathcal{F}_{x,}}_{i-\frac{1}{2},j}^{\text{high}}\Big{)}$
(36b)
$\displaystyle\vec{F_{z,}}_{i,j+\frac{1}{2}}=\vec{\mathcal{F}_{z,}}_{i,j+\frac{1}{2}}^{\text{low}}-\phi_{z}(r_{z,i})\Big{(}\vec{\mathcal{F}_{z,}}_{i,j+\frac{1}{2}}^{\text{low}}-\vec{\mathcal{F}_{z,}}_{i,j+\frac{1}{2}}^{\text{high}}\Big{)}$
(36c)
$\displaystyle\vec{F_{z,}}_{i,j-\frac{1}{2}}=\vec{\mathcal{F}_{z,}}_{i,j-\frac{1}{2}}^{\text{low}}-\phi_{z}(r_{z,i-1})\Big{(}\vec{\mathcal{F}_{x,}}_{i,j-\frac{1}{2}}^{\text{low}}-\vec{\mathcal{F}_{z,}}_{i,j-\frac{1}{2}}^{\text{high}}\Big{)}$
(36d)

where $\phi_{x},\phi_{z}$ are the flux limiter functions in $x$\- and
$z$-directions, and $r$ represents the ratio of successive gradients on the
mesh:

$\displaystyle r_{i}$ $\displaystyle=$
$\displaystyle\frac{\vec{U}_{i}-\vec{U}_{i-1}}{\vec{U}_{i+1}-\vec{U}_{i}}.$
(37)

Finally, using an explicit scheme, the time stepping for the solution vector
is:

$\displaystyle\vec{U}_{i,j}^{n+1}=\quad\vec{U}_{i,j}^{n}$ $\displaystyle-$
$\displaystyle\frac{\Delta t}{\Delta
x}\Big{(}\vec{F_{x,}}_{i+\frac{1}{2},j}-\vec{F_{x,}}_{i-\frac{1}{2},j}\Big{)}$
(38) $\displaystyle-$ $\displaystyle\frac{\Delta t}{\Delta
z}\Big{(}\vec{F_{z,}}_{i,j+\frac{1}{2}}-\vec{F_{z,}}_{i,j-\frac{1}{2}}\Big{)}$
$\displaystyle+$ $\displaystyle\Delta t\vec{S}_{i,j}^{n}.$

Depending on the flux limiter functions, different schemes can be obtained
from the blended scheme (38). More specifically, if all values of the flux
limiters are set to $1$, then the solution is smooth and the fluxes are
represented by a high-resolution scheme. Substituting all values of the flux
limiter functions with $1$ retrieves the Lax-Wendroff scheme. Another example
is the Lax-Friedrichs scheme which can be obtained if all flux limiter values
are $0$, which means that a low-resolution approximation of the fluxes is
needed.

The preformed simulations in this work use also the minmod limiter functions
which are of the following kind:

$\displaystyle\phi_{x}=\max[0,\min(1,r_{x})],\quad\lim_{r_{x}\to\infty}\phi_{x}(r_{x})=1,$
(39)
$\displaystyle\phi_{z}=\max[0,\min(1,r_{z})],\quad\lim_{r_{z}\to\infty}\phi_{z}(r_{z})=1.$
(40)

The numerical stencil for these two-dimensional schemes is presented in Fig.
11. The horizontal plane defined by $i$\- and $j$-directions (the $xz$-domain)
physically represents the liquid film height. The third dimension of the
stencil $n$ represents the time.

Figure 11: Stencil for the implemented numerical schemes. The spatial
dimensions are represented by $i$ and $j$, and the time steps are indicated by
$n$.

## References

  * Kalliadasis  _et al._ [2012] S. Kalliadasis, C. Ruyer-Quil, B. Scheid, and M. G. Velarde,  _Falling Liquid Films_ (Springer London, 2012).
  * Shkadov [1968] V. Shkadov, “Wave flow theory for a thin viscous liquid layer,” Fluid Dynamics - FLUID DYN 3, 12–15 (1968).
  * Mendez  _et al._ [2021] M. Mendez, A. Gosset, B. Scheid, M. Balabane, and J.-M. Buchlin, “Dynamics of the jet wiping process via integral models,” Journal of Fluid Mechanics 911 (2021), 10.1017/jfm.2020.1075.
  * Aniszewski  _et al._ [2020] W. Aniszewski, Y. Saade, S. Zaleski, and S. Popinet, “Planar jet stripping of liquid coatings: Numerical studies,” International Journal of Multiphase Flow 132, 103399 (2020).
  * Barreiro-Villaverde, Gosset, and Mendez [2021] D. Barreiro-Villaverde, A. Gosset, and M. A. Mendez, “On the dynamics of jet wiping: Numerical simulations and modal analysis,” Physics of Fluids 33, 062114 (2021).
  * Kapitza [1949] P. Kapitza, “Wave flow of thin layers of a viscous fluid, i. the free flow,” Zh. Eksperim. Theor. Fiz. 18, 541–546 (1949).
  * Demekhin and Shkadov [1984] E. Demekhin and V. Shkadov, “Three-dimensional waves in a liquid flowing down a wall,” Fluid Dynamics 19, 689–695 (1984).
  * Demekhin and Kalaidin [2007] E. A. Demekhin and E. N. Kalaidin, “Three-dimensional localized coherent structures of surface turbulence. i. scenarios of two-dimensional–three-dimensional transition,” Physics of Fluids 19 (2007), https://doi.org/10.1063/1.2793148.
  * Demekhin  _et al._ [2010] E. A. Demekhin, E. N. Kalaidin, S. Kalliadasis, and S. Y. Vlaskin, “Three-dimensional localized coherent structures of surface turbulence: Model validation with experiments and further computations,” Physical Review E 82 (2010), https://doi.org/10.1103/PhysRevE.82.036322.
  * Ruyer-Quil and Manneville [1998] C. Ruyer-Quil and P. Manneville, “Modeling film flows down inclined planes,” The European Physical Journal B 6, 277–292 (1998).
  * Ruyer-Quil and Manneville [2000] C. Ruyer-Quil and P. Manneville, “Improved modeling of flows down inclined planes,” The European Physical Journal B 15, 357–369 (2000).
  * Scheid, Ruyer-Quil, and Manneville [2006] B. Scheid, C. Ruyer-Quil, and P. Manneville, “Wave patterns in film flows: modelling and three-dimensional waves,” Journal of Fluid Mechanics 562, 183 (2006).
  * Ruyer-Quil  _et al._ [2014] C. Ruyer-Quil, N. Kofman, D. Chasseur, and S. Mergui, “Dynamics of falling liquid films,” The European physical journal. E, Soft matter 37, 9982 (2014).
  * Dietze  _et al._ [2014] G. Dietze, W. Rohlfs, K. Nährich, R. Kneer, and B. Scheid, “Three-dimensional flow structures in laminar falling liquid films,” Journal of Fluid Mechanics 743, 75–123 (2014).
  * Rohlfs, Rietz, and Scheid [2018] W. Rohlfs, M. Rietz, and B. Scheid, “Wavemaker : The three-dimensional wave simulation tool for falling liquid films,” SoftwareX 7, 211–216 (2018).
  * Buchlin, J.M. [1997] Buchlin, J.M., “Modelling of Gas Jet Wiping,” in  _Thin Liquid Films and Coating Processes,VKI Lecture Series_ (von Karman Institute for Fluid Dynamics, 1640 Rhode-Saint-Genese, Belgium, 1997).
  * Gosset, Mendez, and Buchlin [2019] A. Gosset, M. A. Mendez, and J.-M. Buchlin, “An experimental analysis of the stability of the jet wiping process: Part i – characterization of the coating uniformity,” Experimental Thermal and Fluid Science 103, 51–65 (2019).
  * Mendez, Gosset, and Buchlin [2019] M. Mendez, A. Gosset, and J.-M. Buchlin, “Experimental analysis of the stability of the jet wiping process, part II: Multiscale modal analysis of the gas jet-liquid film interaction,” Experimental Thermal and Fluid Science 106, 48–67 (2019).
  * Hocking  _et al._ [2010] G. C. Hocking, W. L. Sweatman, A. D. Fitt, and C. Breward, “Deformations during jet-stripping in the galvanizing process,” Journal of Engineering Mathematics 70, 297–306 (2010).
  * Shkadov [1977] V. Y. Shkadov, “Solitary waves in a layer of viscous liquid,” Akademiia Nauk SSSR Izvestiia Mekhanika Zhidkosti i Gaza 12, 63–66 (1977).
  * Mendez, Scheid, and Buchlin [2017] M. Mendez, B. Scheid, and J.-M. Buchlin, “Low kapitza falling liquid films,” Chemical Engineering Science 170, 122–138 (2017).
  * Shkadov [1970] V. Y. Shkadov, “Wave flow regimes of a thin layer of viscous fluid subject to gravity,” Fluid Dynamics 2, 29–34 (1970).
  * Shampine [2005] L. F. Shampine, “Solving hyperbolic PDEs in MATLAB,” Applied Numerical Analysis & Computational Mathematics 2, 346–358 (2005).
  * Doro and Aidun [2013] E. O. Doro and C. K. Aidun, “Interfacial waves and the dynamics of backflow in falling liquid films,” Journal of Fluid Mechanics 726, 261–284 (2013).
  * Fitzpatrick [2019] R. Fitzpatrick,  _Oscillations and waves : an introduction_ (CRC Press/Taylor and Francis Group, Boca Raton, FL, 2019).
  * Lighthill [2007] M. J. Lighthill,  _Waves in Fluids_ (Cambridge University Press, 2007).
  * Graebel [2007] W. P. Graebel,  _Advanced fluid mechanics_ (Academic Press, Burlington, MA, 2007).


# Time Evolution of a Decaying Quantum State:
Evaluation and Onset of Non-Exponential Decay

Markus Nöth<EMAIL_ADDRESS>Mathematisches Institut, Ludwig-Maximilians-
Universität München, Theresienstr. 39, D-80333 München, Germany

###### Abstract

The method developed by van Dijk, Nogami and Toyama [1, 2, 3, 4] for obtaining
the time-evolved wave function of a decaying quantum system is generalized to
potentials and initial wave functions of non-compact support. The long time
asymptotic behavior is extracted and employed to predict the timescale for the
onset of non-exponential decay. The method is illustrated with a Gaussian
initial wave function leaking through Eckart’s potential barrier on the
halfline.

††preprint: APS/123-QED

##  I Introduction

Quantum tunneling, the phenomenon of a particle traversing a potential barrier
higher than its kinetic energy, has received a great deal of attention since
the early works of Gamow [5]. Despite the time and effort spent on this issue
(see [6, 7, 8, 9] for important contributions to the discussion and reviews
and [10] for a recent experiment), there still remain open questions about
tunneling. One reason for this is that while tunneling is fundamentally a
process in time, much of the analysis has been restricted to static aspects of
the relevant setups. While important properties of tunneling systems such as
exponential decay rates describing a decaying system for intermediate times
can be found in such an analysis, a more thorough understanding seems
desirable.

For this reason studying the time dependent wave function describing the
tunneling particle is essential to complete the picture. For the setup of a
decaying quantum state a great deal of progress in this direction has been
achieved by van Dijk, Nogami and Toyama [1, 2, 3, 4]. For a precursor to the
technique see [11]. They found a representation of the time dependent wave
function well suited for studying early, intermediate and late times. This
makes it well suited for analytic as well as numerical endeavors. The method
was only applied to compactly supported potentials and initial wave functions.
The compact support of the initial wave function played an important role, as
will be explained in the next section. In this article I extend this method
and apply it to a Gaussian initial wave packet leaking out of Eckart’s
potential barrier; neither of which has a compact support. The necessary
modification of the technique is introduced first.

Following a reminder of the tricks employed and the generalization of the
present work I analyze the late-time asymptotics. While the power-law for the
decay of the wave function

$|\psi(r,t)|\overset{t\rightarrow\infty}{\rightarrow}\frac{|\psi_{\infty}(r)|}{t^{3/2}},$
(1)

in the generic case is well-known, I derive a simple and fairly general
expression for the constant involved cf. [12]. Furthermore, an estimate of the
time $t_{\text{alg}}$ after which the transition from exponential to algebraic
decay takes place is given in (29). This result should be compared to the one
for the survival amplitude by [13]. For a complicated barrier it may not be
feasible to apply the full technique while the expressions for
$|\psi_{\infty}(r)|$ and $t_{\text{alg}}$ can be evaluated more easily and may
therefore be of greater importance when comparing with experiments.

##  II Recap and Extension

###  II.1 Main Argument

We consider the Schrödinger equation (setting $\hbar=2m=1$)

$i\partial_{t}\psi=-\partial_{r}^{2}\psi+V\psi,$ (2)

for the scalar wave function $\psi$ subject to the potential $V$ and boundary
condition $\psi(0,t)=0$ for all $t$, chosen such that there are no bound
states 111See [3] for how to modify the technique in the presence of bound
states. This theoretical setup may describe the radial part of a system
describing alpha-decay with zero angular momentum [15] or it may describe the
motion of an ion along the axis of a quantum optical trap inducing a back wall
[16, 17] and an extra trapping potential $V$. We will make use of the Jost
solutions $f(k,r)$ of (2) which are the generalized eigenfunctions of the
Hamiltonian of (2) of energy $k^{2}$ that satisfy

$\lim_{r\rightarrow\infty}e^{-ikr}f(k,r)=1.$ (3)

Jost solutions exist for potentials $V$ such that
$\int_{0}^{\infty}dr\,r|V(r)|<\infty$. For fixed real $r$ the Jost solution
$f(k,r)$ is analytic on $\\{\operatorname{Im}(k)\geq 0\wedge k\neq
0\\}\subset\mathbb{C}$ and fulfills $\overline{f(k,r)}=f(-\overline{k},r)$
[18, Thm XI.57]. Furthermore, I assume that 0 is not an eigenvalue or
resonance of the Hamiltonian involved, which is typically true.

The generalized eigenfunction $u(k,r)$ of the same energy satisfying the
proper boundary condition at the origin and the additional normalization

$u(k,0)=0,\quad u^{\prime}(k,0)=1,$ (4)

where the prime denotes a derivative with respect to the second argument, is
given by

$u(k,r)=\big{(}f(k,0)f(-k,r)-f(-k,0)f(k,r)\big{)}\frac{1}{2ik}.$ (5)

Furthermore, $u(k,r)$ is real for real $r$ and $k$ and even in $k$ and for $r$
fixed $u(k,r)$ is an entire function of $k$ [18, Thm XI.56]. The solutions $u$
form a complete set and are normalized to [1]

$\int_{0}^{\infty}\overline{u(k,r)}u(k^{\prime},r)dr=\frac{\pi}{2k^{2}}|f(k,0)|^{2}\delta(k-k^{\prime}),$
(6)

for $k,k^{\prime}>0$. So we may expand any solution $\psi(x,t)$ of equation
(2):

$\displaystyle\psi(r,t)=\frac{2}{\pi}\int_{0}^{\infty}\frac{k^{2}}{f(k,0)f(-k,0)}C(k)e^{-ik^{2}t}u(k,r)dk,$
(7)

with

$C(k)=\int_{0}^{\infty}u(k,r)\psi(r,0)dr,$ (8)

where I eliminated complex conjugates by the above mentioned properties of $u$
and $f$. Because of the symmetry properties of the integrand I split $u$ up as
in equation (5) to obtain

$\psi(r,t)=\int_{-\infty}^{\infty}e^{-ik^{2}t+ikr}\frac{C(k)}{i\pi}\frac{k\,e^{-ikr}f(-k,r)}{f(-k,0)}dk,$
(9)

where I added an additional plane wave so that the last factor approaches a
finite limit for $r\rightarrow\infty$. Next I introduce an auxiliary factor
$h(k)$

$\psi(r,t)=\frac{1}{i\pi}\int_{-\infty}^{\infty}e^{-ik^{2}t+ikr}h(k)\overbrace{\frac{k\,C(k)e^{-ikr}f(-k,r)}{h(k)f(-k,0)}}^{=:g(k,r)}dk,$
(10)

to achieve the following properties:

  * •

The fraction $g(k,r)$ is still meromorphic in $k$.

  * •

The expansion of $g(k,r)$ by the Mittag-Leffler [19] theorem is of the form

$g(k,r)=\sum_{n}\frac{a_{n}(r)}{k-k_{n}},$ (11)

where $a_{n}$ are the residues of $g$ at it’s poles $k_{n}$. A sufficient
criterion for (11) is that $g(k,r)$ falls off at infinity in the complex
$k-$plane (excluding small circles around the poles) and the sum in (11)
converges uniformly on compact sets.

  * •

The integral that remains after exchange of integral and sum which is of the
form

$\frac{1}{i\pi}\int_{-\infty}^{\infty}\frac{e^{-ik^{2}t+ikr}h(k)}{k-k_{n}},$
(12)

may be evaluated in closed form.

In [1, 2, 3, 4], the initial wave functions $\psi(r,0)$ under consideration
vanishes for any $r>R$ for some $R>0$, this results in growth of the fractions
of (9) like $e^{-\mathrm{Im}(k)R}$ in negative imaginary direction. The choice
$h=e^{ikR}$ exactly cancelles this growth so that Mittag-Leffler can be
applied.

The central innovation of the current paper is a different choice of $h$ so as
to allow for more generality:

$h_{\alpha}(k)=e^{i\alpha k^{2}}+e^{-i\alpha k^{2}}+e^{-\alpha k^{2}},$ (13)

where $\alpha>0$ is a free parameter.

To give an idea of which initial wave functions can be treated with this
method a quick estimate of the growth of $C(k)$ in imaginary direction may be
helpful. From general scattering theory we know that the eigenfunctions $u$
generically grow exponentially in the complex plane see [20, equation (12.8)],
[21, equations (6.4.13) and (6.4.17)]

$|u(k,r)|\leq\frac{c_{1}}{|k|}e^{|\operatorname{Im}(k)|r}.$ (14)

If we assume that $\psi(r,0)=\mathcal{O}\left(e^{-c_{2}r^{c_{3}}}\right)$ for
some $c_{2}>0,c_{3}>1$ for $r\rightarrow\infty$ big enough a stationary phase
argument yields

$|C(k)|\leq
c_{4}e^{(c_{2}c_{3})^{\frac{-1}{c_{3}-1}}(\left|\mathrm{Im}(k)\right|)^{\frac{c_{3}}{c_{3}-1}}(1-1/c_{3})}.$
(15)

If this growth can be controlled by a Gaussian the choice (13) will work for
some $\alpha>0$. Hence for any initial wave function with $c_{3}\geq 2$ the
method will succeed.

The choice (13) results in the following expression for the wave function

$\displaystyle\psi(r,t)=\sum_{n}a_{n}(r)\sum_{\gamma\in\\{\alpha,\pm
i\alpha\\}}M(k_{n},r,\gamma$ $\displaystyle+it),$ (16)

where $M$ is given by

$\displaystyle
M(k,r,\beta)=\int_{-\infty}^{\infty}\frac{dp}{i\pi}\frac{e^{-\beta
p^{2}+irp}}{p-k}$ (17) $\displaystyle=e^{ikr-\beta
k^{2}}\bigg{[}\text{sgn}(\text{Im}(k))+\operatorname{erf}\big{(}(r+2ik\beta)/(2\sqrt{\beta})\big{)}\bigg{]}$

and is closely related to the Moshinsky function [22]. In order to apply the
method one has to include all the poles of $g$, including the zeros of
$h_{\alpha}$, which are located close to the three sets

$\displaystyle\left\\{\pm e^{-\frac{3\pi
i}{8}}\sqrt{\pi(2n-1)/(\sqrt{2}\alpha)}\,\bigg{|}n\in\mathbb{N}\right\\},$
(18) $\displaystyle\left\\{\pm e^{\frac{3\pi
i}{8}}\sqrt{\pi(2n-1)/(\sqrt{2}\alpha)}\,\bigg{|}n\in\mathbb{N}\right\\},$
(19)
$\displaystyle\left\\{\pm\sqrt{\pi(2n-1)/(2\alpha)}\,\bigg{|}n\in\mathbb{N}\right\\}.$
(20)

###  II.2 Late Time Asymptotics

Besides the apparent use of the representation (16) to calculate $\psi(r,t)$
numerically with small error for early as well as for late times it is also
useful for theoretical considerations which is illustrated by this chapter. I
will discuss the late time asymptotics of a decaying quantum system recovering
the well known result

$\psi(r,t)\xrightarrow{t\rightarrow\infty}\frac{\psi_{\infty}(r)}{t^{3/2}},$
(21)

valid vor generic wave functions. Please note, that for special initial wave
functions a different power law may apply [23]. The asymptotic behavior (21)
was e.g. also obtained by [3, 8, 12] for rigorous bounds see [21, 24, 25]. The
amplitude $\psi_{\infty}(r)$ will be obtained for any system to which the
method of the last section is applicable and to arbitrary accuracy in the
expansion (16). This expression is simpler than the one of [12]. For $t$ large
relative to $r$ and $\alpha$, I employ the asymptotic expression for
$\operatorname{erf}$ [26, page 109-112]222formulated in terms of lower
incomplete gamma function

$\displaystyle\operatorname{erf}(\pm z)$ $\displaystyle=\pm
1-\frac{e^{-z^{2}}}{\sqrt{\pi}}\sum_{k=0}^{M}\frac{(-1)^{k}\left(\frac{1}{2}\right)_{k}}{(\pm
z)^{2k+1}}+\mathcal{O}(e^{-z^{2}}z^{-2M-3}),$ (22)

valid for $|\operatorname{ph}(\pm z)|<\frac{3}{4}\pi$ where the phase of
$k_{n}$ decides which of the two expansions is to be used and $(c)_{k}$ is the
Pochhammer-symbol. Taking $\gamma\in\\{\alpha,\pm i\alpha\\}$ and using (22) I
get

$\displaystyle\\!M(k_{n},r,\gamma+it)\approx
e^{ik_{n}r-(it+\gamma)k_{n}^{2}}(\operatorname{sgn}(\operatorname{Im}(k_{n}))+\nu)$
(23)
$\displaystyle-\frac{\sqrt{it+\gamma}e^{\frac{-r^{2}}{4(it+\gamma)}}}{\sqrt{\pi}(\frac{r}{2}+ik_{n}(it+\gamma))}\left(1-\frac{(it+\gamma)/2}{\big{(}\frac{r}{2}+ik_{n}(it+\gamma)\big{)}^{2}}\right),$

with $\nu=-1$ for $k_{n}=|k_{n}|e^{i\varphi}$ and
$-\frac{1}{4}\pi\leq\varphi\leq\frac{3}{4}\pi$ and $\nu=1$ otherwise. I
omitted errors of order $\mathcal{O}(t^{-5/2})$. There are several remarks in
order:

  * •

Except for the case $k_{n}\in\mathbb{R}$ the first summand in (23) is either
identically zero or decays exponentially in time for their respective region
of $\varphi$ and may therefore be absorbed into the error.

  * •

For a $k_{n}\in\mathbb{R}$, the first term is oscillatory. However, in this
case $h_{\alpha}(k_{n})=0$ as the Jost solution $f(-k,0)$ has no zeros and
$C(k)$ no poles on the real line. Hence the first term of (23) vanishes in the
sum over $\gamma\in\\{\alpha,\pm i\alpha\\}$.

  * •

The complicated fractions involving $k_{n},r$ and $t$ may be represented by a
series in inverse powers of $t$, which will also be a series in inverse powers
of $k_{n}$. We may then exploit the absence of a constant term in (11) to
obtain

$\displaystyle 0$ $\displaystyle=\sum_{n}\frac{a_{n}(r)}{k_{n}},$ (24)
$\displaystyle\frac{C(0)}{3}\frac{f(0,r)}{f(0,0)}$
$\displaystyle=-\sum_{n}\frac{a_{n}(r)}{k_{n}^{2}}$ (25)
$\displaystyle\frac{C(0)}{3}\partial_{k}\left.\frac{e^{-ikr}f(-k,r)}{f(-k,0)}\right|_{k=0}$
$\displaystyle=-\sum_{n}\frac{a_{n}(r)}{k_{n}^{3}}$ (26)

by inserting $k=0$ into (11) and its derivatives. This yields a simpler
expression for the asymptotically late wave function.

As a next step I employ the binomial series to express the fractions in (23)
and sum over the three possibilities of $\gamma\in\\{\alpha,\pm i\alpha\\}$
and over all poles $k_{n}$, this results in

$\displaystyle\psi_{\infty}(r)=-\frac{C(0)}{2\sqrt{i\pi}}\partial_{k}\left.\frac{f(-k,r)}{f(-k,0)}\right|_{k=0}.$
(27)

Corrections to (21) are of order $\mathcal{O}\left(t^{-5/2}\right)$. If
$C(0)=0$ the decay of $\psi(r,t)$ for large $t$ will be faster, as observed in
[23]. Higher orders can be worked out analogously to extract an asymptotic
series for $\psi(r,t)$ in $1/\sqrt{t}$, which implies a series of the same
kind for

$P(t)=\int_{0}^{\rho}|\psi(r,t)|^{2}dr,$ (28)

the probability of finding the particle inside the trapping region of length
$\rho$. The use of the binomial series is only justified for late times, where
what “late” means depends on $r$. This onset of the algebraic decay can be
estimated by equating the exponential term with slowest decay of (23) with the
term (27), yielding:

$t_{\text{alg}}=\frac{-3/2}{|\operatorname{Im}k_{0}^{2}|}W_{-1}\\!\\!\left(\frac{-|\operatorname{Im}k_{0}^{2}|}{3/2^{1/3}}\left|\frac{\psi_{\infty}(r)\partial_{k_{0}}f(-k_{0},0)}{k_{0}C(k_{0})f(-k_{0},r)}\right|^{2/3}\right),$
(29)

where $W_{-1}$ is the $-1$-branch of the Lambert W function and $k_{0}$ is the
wave vector belonging to the slowest exponential decay, i.e. the zero of
$f(-k,0)$ in the fourth quadrant of the complex plane closest to the origin.

Besides (28) the survival probability

$S(t)=\left|\int_{0}^{\infty}\psi^{*}(r,0)\psi(r,t)dr\right|^{2},$ (30)

is often used to study decaying quantum states. The long time limit of $S(t)$
can be found pluggin in (7) and (6):

$S(t)=\left|\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{k^{2}C^{\prime}(k)C(k)}{f(k,0)f(-k,0)h_{\alpha}(k)}e^{-ik^{2}t}h_{\alpha}(k)dk\right|^{2},$
(31)

with

$C^{\prime}(k)=\int_{0}^{\infty}u(k)\psi^{*}(r)dr.$ (32)

So that Mittag-Leffler can be applied again to obtain

$S(t)=\left|\sum_{n}b_{n}\sum_{\gamma\in\\{\alpha,\pm
i\alpha\\}}M(k_{n},0,\gamma+it)\right|^{2}.$ (33)

Using the same expansion and manipulations as for the the long time limit of
$\psi$ then yields:

$S(t)\xrightarrow{t\rightarrow\infty}\frac{1}{t^{3}}~{}\frac{|C(0)|^{4}}{4\pi(f(0,0))^{4}},$
(34)

showing that $S$ has the same power law as $P$ and also decays faster if
$C(0)=0$. Evaluating expansions (16) and (33) without closed form expressions
for $f$ is difficult, this is not the case for equations (27) and (34). They
may be evaluated numerically for any potential $V$ such that the growth
conditions on the Jost solution $f$ and the integrals over the initial wave
function $C$ are fulfilled. Hence these considerations of late times may be of
some use to test non-relativistic models of alpha-decay when compared to
experiment.

##  III Application to Eckart’s potential

Eckart’s potential

$V(r)=\frac{Ae^{r-\rho}}{\left(1+e^{r-\rho}\right)^{2}},$ (35)

where $\rho\in\mathbb{R}$ and $A\in\mathbb{R}$ are free parameters together
with the initial wave function

$\psi_{0}(r)=\frac{2^{5/2}}{\rho^{3/2}\sqrt[4]{\pi}}re^{-2\left(\frac{r}{\rho}\right)^{2}}$
(36)

provides a nice example for which an exponential choice of $h$ would not
suffice. This potential has been used to study tunneling [28] from the
perspective of theoretical physics as well as from the perspective of
chemistry [29]. The associated Jost solution is given by [8, 30]

$f(k,r)=e^{-ikr}{}_{2}F_{1}\Big{(}\frac{1}{2}-i\delta,\frac{1}{2}+i\delta;1+2ik;\frac{1}{1+e^{r-\rho}}\Big{)},$
(37)

with $\delta=\sqrt{A-1/4}$, where ${}_{2}F_{1}$ is the Gauss hypergeometric
function.

Figure 1:  Numerical logarithmic comparison of $|\psi(r,t)|$,
$|\psi_{\infty}(r)|/t^{1.5}$ and the slowest exponential decay corresponding
to the zero $k_{0}$ of $f(-k,0)$ closest to the origin in the fourth quadrant.
The thin vertical line marks $t_{\text{alg}}$. The parameters are
$A=49.25,\rho=1,r=0.5$ and $\alpha=1.25$. To compute $\psi(r,t)$ all summands
of (16) corresponding to poles $|k_{n}|<100$, which results in approximately
18000 summands. The code that produced this figure can be found here:
%****␣script.tex␣Line␣450␣****https://gitlab.com/therunomask/tunneling, it
uses mpmath [31].

In this setup the Mittag-Leffler theorem implies (11)

$g(k,r)=\frac{kC(k)e^{-ikr}f(-k,r)}{h_{\alpha}(k)f(-k,0)}=\sum_{n}\frac{a_{n}(r)}{k-k_{n}}.$
(38)

The main reason why my choice of $h_{\alpha}$ works while an exponential one
does not is as follows:

The generalized Fourier transform of of the initial wave function, $C(k)$,
satisfies the bound

$|C(k)|\leq\frac{\sqrt{2}\sqrt{\rho}c_{1}}{|k|\sqrt[4]{\pi}}+4c_{1}\sqrt[4]{\pi}\rho^{3/2}e^{\rho^{2}|\operatorname{Im}(k)|^{2}},$
(39)

which can be derived using (14). The $\leq$ becomes $\approx$ for $|k|$ large
enough and $k\not\in\mathbb{R}$. This can be compensated for by my choice of
$h_{\alpha}$ for $\alpha$ big enough, while clearly an exponential $h$ does
not suffice.

A detailed estimate on the real line shows

$C(k)=\mathcal{O}(k^{-6}),\quad k\rightarrow\infty$ (40)

implying that the sum (16) converges. In order to evaluate (16), the pole
structure of $g$ should be studied. Here, in addition to the artificial poles
close to (18), (19) and (20) there are poles due to the Jost solution $f$. It
has poles at $k=in/2$ for $n\in\mathbb{N}$, however the poles cancel in the
expression for $g$, as they appear in the numerator as well as the
denominator. The function $f(-k,0)$ has zeros for complex $k$, symmetric about
the imaginary axis. I conjecture the right arm of which approaches

$k_{n}=\frac{n-\frac{1}{4}+\frac{i}{2\pi}\ln\Big{(}e^{\pi\delta}+e^{-\pi\delta}\Big{)}}{i+\frac{\rho}{\pi}}$
(41)

for large $n$. This conjecture is supported by numerical evidence, but I do
not have a proof yet.

In figure 1 all of the above is used to check numerically how accurately
$t_{\text{alg}}$ reflects the transition from exponential to algebraic decay
for the arbitrary choice $r=1/2$. The behavior seen in figure 1 is consistent
with the analysis in [3].

Figure 2:  Comparison of evaluation of (16) with $A=49.25$, $\rho=1$ and
$\alpha=1.25$ with the Crank-Nicolson (CN) method at $t=3$ and $t=100$. All
poles with $|k_{n}|\leq 120$ are taken into account for (16). The
discretization parameters for CN are $\Delta r=0.01953125$ and $\Delta
t=\Delta r^{2}/4$, the length of the simulation box is $L=650$ for $t=3$ and
$L=4480$ for $t=100$. The CN method is implemented in Julia [32], the code
that produces this figure can also be found here
https://gitlab.com/therunomask/tunneling.

In figure 2 the evaluation of (16) is compared with the well established
Crank-Nicolson method at an early ($t=3$) and moderately late ($t=100$) time.
One can see clearly that the two methods agree reasonably well at early times
and the absolute value of $|\psi|$ still agrees at $t=100$ while substantial
phase difference has accumulated leading a large difference.

##  IV Conclusion

In the present paper I have generalized the method of van Dijk, Nogami, and
Toyama [1, 2, 3, 4] to a large class of non-compactly supported potentials and
initial wave functions which includes Gaussian initial states. This method was
used to derive the asymptotic amplitude $|\psi_{\infty}(r)|$ as well as
estimate the onset of algebraic decay $t_{\text{alg}}$ characteristic of large
times. Finally, the method was applied to a Gaussian wave packet leaking
through Eckart’s potential against a hard-backwall potential barrier.

##  V Disclaimer

This article reflects the author’s personal work and does not represent
scientific standpoints of the author’s employer.

##  VI Acknowledgments

I wholeheartedly thank Christian Bild for meticulously checking my
calculations as well as Siddhant Das for helpful discussions.

## References

  * van Dijk and Nogami [1999] W. van Dijk and Y. Nogami, Novel expression for the wave function of a decaying quantum system, Physical review letters 83, 2867 (1999).
  * Nogami  _et al._ [2000] Y. Nogami, F. Toyama, and W. van Dijk, Bohmian description of a decaying quantum system, Physics Letters A 270, 279 (2000).
  * van Dijk and Nogami [2002] W. van Dijk and Y. Nogami, Analytical approach to the wave function of a decaying quantum system, Physical Review C 65, 024608 (2002).
  * van Dijk and Toyama [2019] W. van Dijk and F. M. Toyama, Decay of a quasistable quantum system and quantum backflow, Physical Review A 100, 052101 (2019).
  * Gamow [1928] G. Gamow, Zur quantentheorie des atomkernes, Zeitschrift für Physik 51, 204 (1928).
  * Hartman [1962] T. E. Hartman, Tunneling of a wave packet, Journal of Applied Physics 33, 3427 (1962).
  * Muga  _et al._ [2007] G. Muga, R. S. Mayato, and I. Egusquiza,  _Time in quantum mechanics_ , Vol. 734 (Springer Science & Business Media, 2007).
  * Razavy [2003] M. Razavy,  _Quantum theory of tunneling_ (World Scientific, 2003).
  * Hauge and Støvneng [1989] E. Hauge and J. Støvneng, Tunneling times: a critical review, Reviews of Modern Physics 61, 917 (1989).
  * Ramos  _et al._ [2020] R. Ramos, D. Spierings, I. Racicot, and A. M. Steinberg, Measurement of the time spent by a tunnelling atom within the barrier region, Nature 583, 529 (2020).
  * Garcıa-Calderón  _et al._ [1996] G. Garcıa-Calderón, J. Mateos, and M. Moshinsky, Survival and nonescape probabilities for resonant and nonresonant decay, annals of physics 249, 430 (1996).
  * Jensen and Kato [1979] A. Jensen and T. Kato, Spectral properties of schrödinger operators and time-decay of the wave functions, Duke mathematical journal 46, 583 (1979).
  * Peres [1980] A. Peres, Nonexponential decay law, Annals of Physics 129, 33 (1980).
  * Note [1] See [3] for how to modify the technique in the presence of bound states.
  * Lovas  _et al._ [1998] R. G. Lovas, R. Liotta, A. Insolia, K. Varga, and D. Delion, Microscopic theory of cluster radioactivity, Physics reports 294, 265 (1998).
  * Das and Dürr [2019] S. Das and D. Dürr, Arrival time distributions of spin-1/2 particles, Scientific reports 9, 1 (2019).
  * Das  _et al._ [2019] S. Das, M. Nöth, and D. Dürr, Exotic bohmian arrival times of spin-1/2 particles: an analytical treatment, Physical Review A 99, 052124 (2019).
  * Reed and Simon [1979] M. Reed and B. Simon, Scattering theory. methods of modern mathematical physics, vol. _3., Scattering theory. Methods of modern mathematical physics, Vol. _3., by Reed, M.; Simon, B.. New York, NY (USA): American Press, 15+ 463 p. (1979).
  * Jeffreys  _et al._ [1999] H. Jeffreys, B. Jeffreys, and B. Swirles,  _Methods of mathematical physics_ (Cambridge university press, 1999).
  * Newton [1966] R. G. Newton,  _Scattering theory of waves and particles_ (Springer Science & Business Media, 1966).
  * Vona [2014] N. Vona, On time in quantum mechanics, arXiv preprint arXiv:1403.2496 (2014).
  * Moshinsky [1952] M. Moshinsky, Diffraction in time, Physical Review 88, 625 (1952).
  * Miyamoto [2004] M. Miyamoto, Initial wave packets and the various power-law decreases of scattered wave packets at long times, Physical Review A 69, 042704 (2004).
  * Journé  _et al._ [1991] J.-L. Journé, A. Soffer, and C. D. Sogge, Decay estimates for schrödinger operators, Communications on Pure and Applied mathematics 44, 573 (1991).
  * Schlag [2007] W. Schlag, Dispersive estimates for schrödinger operators: a survey, Mathematical aspects of nonlinear dispersive equations 163, 255 (2007).
  * Olver [1997] F. Olver,  _Asymptotics and special functions_ (CRC Press, 1997).
  * Note [2] Formulated in terms of lower incomplete gamma function.
  * de la Cal  _et al._ [2022] X. G. de la Cal, M. Pons, and D. Sokolovski, Speed-up and slow-down of a quantum particle, Scientific reports 12, 1 (2022).
  * Migliore  _et al._ [2014] A. Migliore, N. F. Polizzi, M. J. Therien, and D. N. Beratan, Biochemistry and theory of proton-coupled electron transfer, Chemical reviews 114, 3381 (2014).
  * Eckart [1930] C. Eckart, The penetration of a potential barrier by electrons, Physical Review 35, 1303 (1930).
  * Johansson  _et al._ [2013] F. Johansson _et al._ ,  _mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18)_ (2013), http://mpmath.org/.
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59, 65 (2017), https://doi.org/10.1137/141000671 .