††thanks: Corresponding author<EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding author<EMAIL_ADDRESS>
# Effects of the initial perturbations on the Rayleigh-Taylor-Kelvin-Helmholtz
instability system

Feng Chen  School of Aeronautics, Shan Dong Jiaotong University, Jinan 250357,
China.  Aiguo Xu  Laboratory of Computational Physics, Institute of Applied
Physics and Computational Mathematics, P. O. Box 8009-26, Beijing 100088,
China.  HEDPS,Center for Applied Physics and Technology, and College of
Engineering, Peking University, Beijing 100871, China.  State Key Laboratory
of Explosion Science and Technology, Beijing Institute of Technology, Beijing
100081, China.  Yudong Zhang  School of Mechanics and Safety Engineering,
Zhengzhou University, Zhengzhou 450001, China.  Yanbiao Gan  Hebei Key
Laboratory of Trans-Media Aerial Underwater Vehicle, North China Institute of
Aerospace Engineering, Langfang 065000, China.  Bingbing Liu  Naval
Architecture and Port Engineering College, Shan Dong Jiaotong University,
Weihai 264200, China.  Shuang Wang  School of science, Shandong Jianzhu
University, Jinan 250101, China.

###### Abstract

In the paper, the effects of initial perturbations on the Rayleigh-Taylor
instability (RTI), Kelvin-Helmholtz instability (KHI), and the coupled
Rayleigh-Taylor-Kelvin-Helmholtz instability (RTKHI) systems are investigated
using a multiple-relaxation-time discrete Boltzmann model. Six different
perturbation interfaces are designed to study the effects of the initial
perturbations on the instability systems. It is found that the initial
perturbation has a significant influence on the evolution of RTI. The sharper
the interface, the faster the growth of bubble or spike. While the influence
of initial interface shape on KHI evolution can be ignored. Based on the mean
heat flux strength $D_{3,1}$, the effects of initial interfaces on the coupled
RTKHI are examined in detail. The research is focused on two aspects: (i) the
main mechanism in the early stage of the RTKHI, (ii) the transition point from
KHI-like to RTI-like for the case where the KHI dominates at earlier time and
the RTI dominates at later time. It is found that the early main mechanism is
related to the shape of the initial interface, which is represented by both
the bilateral contact angle $\theta_{1}$ and the middle contact angle
$\theta_{2}$. The increase of $\theta_{1}$ and the decrease of $\theta_{2}$
have opposite effects on the critical velocity. When $\theta_{2}$ remains
roughly unchanged at $90$ degrees, if $\theta_{1}$ is greater than $90$
degrees (such as the parabolic interface), the critical shear velocity
increases with the increase of $\theta_{1}$, and the ellipse perturbation is
its limiting case; If $\theta_{1}$ is less than $90$ degrees (such as the
inverted parabolic and the inverted ellipse disturbances), the critical shear
velocities are basically the same, which is less than that of the sinusoidal
and sawtooth disturbances. The influence of inverted parabolic and inverted
ellipse perturbations on the transition point of the RTKHI system is greater
than that of other interfaces: (i) For the same amplitude, the smaller the
contact angle $\theta_{1}$, the later the transition point appears; (ii) For
the same interface morphology, the disturbance amplitude increases, resulting
in a shorter duration of the linear growth stage, so the transition point is
greatly advanced.

††preprint: AIP/123-QED

##  I Introduction

Hydrodynamic instability-driven mixing processes are ubiquitous in nature, and
of great significance both for fundamental research and for engineering
application in a number of research fields, including Inertial Confinement
Fusion (ICF), astrophysical phenomena, and supersonic combustion. In the last
several decades, much attention has been paid to the challenging problem
Zhou2017PR-I ; Zhou2017PR-II ; Zhou2019PP ; LiHF2021PRE ; WangLF2009POP ;
epl2010 ; Liang2016PRE ; Liang2019POF ; Liang2021POF ; Zou2017JMES ;
Zou2017PRE ; Zou2019JFM ; RavidNatureC ; chen2021jfm ; Lai2020 ; Lai2021 ;
Tang2021 ; LL2019 ; Tao2020 ; Tao2021 .

In view of its direct consequence in various applications, the effects of
initial perturbations are being explored extensively Miles2004 ; Dimonte2004 ;
Ramaprabhu2005 ; Olson2009 ; Gowardhan2011 ; Doron2011 ; Wei2012 ;
Kuchibhatla2013 ; Liu2013 ; McFarland2013 ; ZhaiJFM2014 ; ZhaiJFM2015 ;
Dell2015 ; Xiao2016 ; Xie2017 ; Kord2019 ; Ding2020 ; Si2021 . Miles et al.
Miles2004  investigated the effect of initial conditions on two-dimensional
(2D) Rayleigh-Taylor instability (RTI) and transition to turbulence in planar
blast-wave-driven systems, and found that the initial conditions have a strong
effect on the time to transition to the quasi-self-similar regime. Dimonte
Dimonte2004  studied the dependency of the self-similar Rayleigh-Taylor bubble
acceleration constant on the initial perturbation amplitude. Ramaprabhu et al.
Ramaprabhu2005  investigated the dependence of the RTI growth coefficient and
the self-similar parameter on the amplitude, the spectral shape, the longest
wavelength imposed, and mode-coupling effects. Gowardhan et al. Gowardhan2011
demonstrated that the initial material interface morphology controls the
evolution characteristics of Richtmyer-Meshkov instability (RMI). Wei et al.
Wei2012  investigated the role of the initial perturbation shape on the RTI
development. It is found that the time when the instability reaches a mean
quadratic growth depends on the initial perturbation shape, and the subsequent
vortical interactions are also very sensitive to details of the initial
perturbation shape. Liu et al. Liu2013  found that the temporal evolution of
the bubble tip velocity is sensitively dependent on the Atwood numbers, the
initial perturbation amplitude and the initial perturbation velocity in
classical RTI. McFarland et al. McFarland2013  investigated the effects of
inclination angle and incident shock Mach number on the inclined interface
RMI. Zhai et al. studied the RMI problems of three light gas interfaces
ZhaiJFM2014  and six heavy gas interfaces ZhaiJFM2015 , including the
different shock wave refraction on the interface and the influence of initial
interface shape on the interface features. Dell et al. Dell2015
systematically studied the effect of the initial perturbation on RMI, and
observed that the initial growth rate of RMI is a non-monotone function of the
initial perturbation amplitude. Xiao et al. Xiao2016  studied the effects of
an initial perturbation on RMI, and found that the evolution of the interface
with large initial amplitude in a low-density nonuniform area is fastest,
while that with a small initial amplitude in a high-density nonuniform area is
slower. Xie et al. Xie2017  showed that the growth of the RTI mixing zone is
substantially retarded by superimposing an optimized additional mode on its
random initial perturbations. Liang et al. Liang2019POF  analyzed the effects
of the initial conditions in terms of the perturbation wavelength and
amplitude. It is found that the instability undergoes a faster growth at the
intermediate stage for a larger wavelength, while the late-time bubble and
spike growth rates are insensitive to the changes of the initially perturbed
wavelength and amplitude. Kord et al. Kord2019  investigated the sensitivity
of objective functions to the initial perturbation amplitudes at various
stages of the RTI, such as the mole fraction, kinetic energy norm. Ding et al.
Ding2020  demonstrated that controlling the layer thickness is an effective
way to modulate the late-stage instability growth, which may be useful for the
target design. Si et al. Si2021  investigated the mode-composition effect on
the three-dimensional (3D) interface RMI development, and found that the mode-
coupling has an evident influence on the bubble evolution.

It is generally believed that in the ICF, in the late stage of RTI, the
materials on both sides penetrate each other to form shear flow, and Kelvin-
Helmholtz instability (KHI) occurs at the interface. However, in nature or
practical engineering applications, such as in the atmosphere and oceans,
air/fuel mixing in combustion chambers, the outer region of supernovae, and
the compression of the fuel capsule in ICF, there is rarely a single
instability. Due to the complexity of the interface, KHI often occurs
simultaneously when RTI develops. It would be interesting to study the effect
of the combined occurrence of two instability mechanisms. A complete
understanding of the relation between the RTI and the KHI is important in
understanding the hydrodynamic instability-driven mixing processes. In
addition to the growth rate and final state, some studies have been made on
the early linear and early nonlinear evolution of the coupled instability
Wang2010 ; Ye2011 ; Mandal2011 ; Olson2011 ; Akula2013 ; Vadivukkarasan2016 ;
Vadivukkarasan2017 ; Vadivukkarasan2020 ; Sarychev2019 ; Brizzolara2021 ;
Chen2020 . Both the instability modes and devolution processes that may be
manifested when two such mechanisms occur simultaneously are of current
interest. Ye et al. Wang2010 ; Ye2011  investigated the competitions between
RTI and KHI in 2D incompressible fluids within a linear growth regime. It is
found that the competition between the RTI and the KHI is dependent on the
Froude number, the density ratio of the two fluids, and the thicknesses of the
density transition layer and the velocity shear layer. Mandal et al.
Mandal2011  investigated the nonlinear evolution of bubble and spike due to
the combined action of RTI and KHI. Olson et al. Olson2011  studied the
coupled RTKHI in the early nonlinear regime, and indicated a complex and non-
monotonic behavior where small amounts of shear in fact decrease the growth
rate. Akula et al. Akula2013  showed that the superposition of shear on RTI at
small Atwood number increases the mixing width and growth rate at early times.
Vadivukkarasan et al. Vadivukkarasan2016 ; Vadivukkarasan2017 ;
Vadivukkarasan2020  described the 3D destabilization characteristics of
cylindrical and annular interfaces under the combined RTKHI, and the effects
of various parameters on the most unstable wavenumbers were studied. Sarychev
et al. Sarychev2019  found that an undulating topography on the interface
coating/base material is resulted from the coupled RTKHI. Brizzolara et al.
Brizzolara2021  observed that the cross-over time can correctly predict the
transition from shear- to buoyancy-driven turbulence, in terms of turbulent
kinetic energy production, energy spectra scaling and mixing layer thickness.
In 2020, we investigated the coupled RTKHI system with a Multiple-Relaxation
Time (MRT) Discrete Boltzmann Model (DBM) Chen2020 . To quantitatively analyse
the coupled RTKHI process, we resort to morphological and non-equilibrium
analysis techniques. We found that both the total boundary length $L$ of the
condensed temperature field and the mean heat flux strength $D_{3,1}$ can be
used to quantitatively judge the main mechanism in the early stage of the
RTKHI system. For the case where the KHI dominates at earlier time and the RTI
dominates at later time, the ending point of linear increasing $L$ or
$D_{3,1}$ can work as a geometric or physical criterion for discriminating the
two stages.

In this paper, we investigate the effects of initial perturbations on the RTI,
KHI and the coupled RTKHI systems with the MRT-DBM. The paper is organized as
follows: Section II briefly introduces the methodology of Discrete Boltzmann
Modeling (DBM) method DBM . Systematic numerical simulations on the effects of
initial perturbations are shown and analyzed in Section III. A brief
conclusion is given in Section IV.

##  II Discrete Boltzmann modeling method

The DBMxu2012 ; xu2015aps ; xu2016 ; Xu2018-Chapter2 ; Xu2021aas ; Xu2021cjcp
; Xu2021aaas  is developed from a hybrid of Lattice Boltzmann method (LBM)
ss2001 ; ss1992 ; ShanChen1 ; Qin2005PRE ; SofoneaSpringer ; SXW2011 ; li2012
; Qian2020 ; Shu2018 ; Tian2011 ; Tian2018 ; Liangpre2014 ; Zhong2020 ;
Cbx2020pre ; QiuRF2020PoF ; QiuRF2021PRE ; SunDK-AML ; SunDK-HMT ; Chai-AMC
and the phase space description method of statistical physics. It is one of
the concrete applications of statistical physics coarse-grained modeling
method in the field of fluid mechanics, and is a further development of the
description method of statistical physical phase space in the form of discrete
Boltzmann equation. Its idea originated from a research review published by Xu
et al. in 2012 xu2012 . In the process of development, it was inspired by the
morphological phase space description method Xu2021cjcp ; Xu2021aaas ; Xu2009
; Xu2010sc ; Xu2016sc . The methodology of DBM is as follows: It selects a
perspective to study a set of kinetic properties of the system according to
research requirements. Therefore, the kinetic moments describing this set of
properties are required to maintain their values in the model simplification.
Based on the independent components of the non-conserved kinetic moments of
$(f-f^{eq})$, construct the phase space, and the phase space and its subspaces
are used to describe the non-equilibrium behaviors of the system. The research
perspective and modeling accuracy should be adjusted as the research
progresses, where $f^{eq}$ is the corresponding equilibrium distribution
function of $f$.

The establishment of a DBM model needs to go through three steps. The first
step is the linearization of the collision term, and the second step is the
discretization of the particle velocity space. The principle for coarse-
grained modeling is that the physical quantities we choose to measure the
system must keep the same values after simplification. The third step is the
purpose and core of DBM modeling, and the specific scheme of non-equilibrium
state and behavior description should be given. DBM is mainly aimed at the
“mesoscale" and “dilemma" situations where continuum modeling fails or
physical functions are insufficient and molecular dynamics method is unable to
do due to the limited applicable scale. The physical information provided by
DBM is between macroscopic continuous description and microscopic molecular
dynamics. Compared with macroscopic description, DBM observes the system from
a wider perspective. The necessity and benefit of non-conserved moment
description in DBM increase with the increase of non-equilibrium degree.

Generally, different from the traditional LBM, whose main function is to
recover macroscopic fluid equations, the DBM does not aim to solve the macro
fluid equations, but mainly to describe and measure both hydrodynamic and
thermodynamic nonequilibrium effects beyond the macro fluid equations. DBM
only puts forward a physical requirement for the used discrete velocities and
does not include specific discrete format. The physical requirement is that
the kinetic moments concerned must keep the value unchanged after being
converted into sums for calculation. DBM does not use the “lattice gas" image
of “virtual particle propagation $+$ collision". This is also one reason why
it should no longer be called the “lattice" Boltzmann method. A more detailed
introduction can be seen in our recent review articles Xu2021aas ; Xu2021cjcp
; Xu2021aaas .

The DBM equation with external force term can be written as :

$\frac{\partial f_{i}}{\partial t}+v_{i\alpha}\frac{\partial f_{i}}{\partial
x_{\alpha}}=-M_{il}^{-1}[\hat{S}_{lk}(\hat{f}_{k}-\hat{f}_{k}^{eq})+\hat{A}_{l}]-g_{\alpha}\frac{(v_{i\alpha}-u_{\alpha})}{RT}f_{i}^{eq}\text{,}$
(1)

where the variable $t$ and $x_{\alpha}$ are the time and spatial coordinates,
$T$ is the temperature, $g_{\alpha}$ and $u_{\alpha}$ denote the acceleration
and velocity in the $x_{\alpha}$ direction, $v_{i\alpha}$ is the discrete
particle velocity, $i=1$, $\ldots$, $N$, and the subscript $\alpha$ indicates
the $x$, $y$, or $z$ component. $f_{i}$ and $\hat{f}_{i}$ ($f_{i}^{eq}$ and
$\hat{f}_{i}^{eq}$) are the particle (equilibrium) distribution functions in
velocity space and kinetic moment space, respectively; the mapping between
moment space and velocity space is defined by the linear transformation
$M_{ij}$, i.e., $\hat{f}_{i}=M_{ij}f_{j}$ and $f_{i}=M_{ij}^{-1}\hat{f}_{j}$.
The matrix $\hat{S}={\rm diag}(s_{1},s_{2},\cdots,s_{N})$ is the diagonal
relaxation matrix. $\hat{A}_{l}$ is the correction term used for recovering
reasonable hydrodynamic behaviors.

In the paper, the following two-dimensional discrete velocity model is used:

$\displaystyle\left(v_{ix,}v_{iy}\right)=\left\\{\begin{array}[]{cc}\mathbf{cyc}\\!:c\left(\pm
1,0\right),&\text{for }1\leq i\leq 4,\\\ c\left(\pm 1,\pm 1\right),&\text{for
}5\leq i\leq 8,\\\ \mathbf{cyc}\\!:2c\left(\pm 1,0\right),&\text{for }9\leq
i\leq 12,\\\ 2c\left(\pm 1,\pm 1\right),&\text{for }13\leq i\leq
16,\end{array}\right.$ (6)

and $\eta_{i}=\eta_{0}$ for $i=1$, …, $4$, and $\eta_{i}=0$ for $i=5$, …,
$16$, which is introduced to control the specific-heat-ratio $\gamma$. “cyc"
indicates the cyclic permutation, and $c$ and $\eta_{0}$ are two free
parameters, which are adjusted to optimize the properties of the model. The
specific forms of the correction term, transformation matrix and the
corresponding equilibrium distribution functions in kinetic moment space can
be seen in the appendix for details.

The kinetic moments of $f_{i}-f_{i}^{eq}$ can be used to check and measure the
Thermodynamic Non-Equilibrium (TNE) state and effects.
$\Delta_{i}^{\ast}=M_{ij}^{\ast}(f_{j}-f_{j}^{eq})$, where
$M_{ij}^{\ast}(f_{i})$ represents kinetic central moments of $f_{i}$, and its
specific expression can be obtained by replacing the variable $v_{i\alpha}$ by
$v_{i\alpha}-u_{\alpha}$ in $M_{ij}$. At present, the commonly used non-
equilibrium characteristic quantities are $\Delta_{2\alpha\beta}^{\ast}$,
$\Delta_{(3,1)\alpha}^{\ast}$, $\Delta_{3\alpha\beta\gamma}^{\ast}$ and
$\Delta_{(4,2)\alpha\beta}^{\ast}$. $\Delta_{2\alpha\beta}^{\ast}$ indicates
the viscous stress tensor and $\Delta_{(3,1)\alpha}^{\ast}$ indicates the heat
flux. $\Delta_{3\alpha\beta\gamma}^{\ast}$ and
$\Delta_{(4,2)\alpha\beta}^{\ast}$ indicate the flux of viscous stress and of
heat flux, respectively. Where, the subscripts “2" and “3" indicate tensors of
second order and third order, the subscript “3,1" means a first-order tensor
contracted from a third-order tensor, and the subscript “4,2" means a second-
order tensor contracted from a fourth-order tensor.

In the phase space opened by the non-conserved kinetic moments and their
subspaces, the origin of coordinates corresponds to a thermodynamic
equilibrium state, and any other point corresponds to a specific thermodynamic
non-equilibrium state. The corresponding non-equilibrium strength can be
defined by means of the distance from the origin. For example, the mean heat
flux strength can be defined as
$D_{(3,1)}=\overline{\sqrt{\Delta_{(3,1)\alpha}^{\ast 2}}}$. These more
abundant but previously poorly understood characteristics of non-equilibrium
behavior contain a large number of physical functions to be explored. At the
same time, it should be pointed out that in addition to the idea of DBM
modeling, a series of data processing techniques, characteristic scale
extraction technology and structure analysis techniques developed in DBM
application process can be directly used for data processing and feature
analysis of various complex configuration and dynamic physical field.

##  III Numerical Simulations

In recent decades, numerous investigations have shown that the initial
conditions have important effects on the evolution of interface instability.
In this section, we investigate the effects of initial perturbations on the
RTI, KHI and the coupled RTKHI systems with the multiple-relaxation-time DBM
model, which has been validated by some well-known benchmark testsChen2016 ;
Chen2018 ; Chen2020 .

Figure 1: The schematic diagram of the initial disturbance interface. (a)
sinusoidal perturbation, (b) sawtooth perturbation, (c) parabolic $y=x^{2}$,
(d) elliptic disturbance, (e) inverted parabolic $y=-x^{2}$, and (f) inverted
elliptic disturbance interface. Figure 2: Schematic diagram. (a) Definition of
the contact angle, $\theta_{1}$ is the bilateral contact angle, $\theta_{2}$
is the middle contact angle. (b) Definition of the bubble amplitude and the
spike amplitude.

Six different initial conditions, shown in Fig. 1, are designed to study the
effects of the initial perturbation shape on the instability system, including
(a) sinusoidal perturbation, (b) sawtooth perturbation, (c) parabolic
$y=x^{2}$, (d) elliptic disturbance, (e) inverted parabolic $y=-x^{2}$, and
(f) inverted elliptic disturbance interface. The angle between the tangent at
the end of the disturbance interface and the boundary is called the bilateral
contact angle $\theta_{1}$, and the angle between the centerline and the
tangent at the center point is called the middle contact angle $\theta_{2}$,
shown in Fig. 2(a). The two contact angles together describe the basic
morphology of the disturbance interface. Except for the sawtooth perturbation,
the middle contact angle $\theta_{2}$ of the other five interfaces is
approximately equal to $90^{o}$. These six interfaces basically include the
continuous curve shape from the boundary to the center line determined by
$\theta_{1}$ and $\theta_{2}$. Their related understanding and recognition can
be helpful to analyze the multimode disturbance in the future.

In the simulation, a two-dimensional computational domain with height $H=80$
and width $W=20$ is adopted, and divided into $NX\times NY=101\times 400$
mesh-cells. The computational domain can be divided into upper and lower parts
by the disturbed interface. In the two parts, two different homogeneous
temperatures are fixed, and the corresponding hydrostatic density profiles
$\rho(y)$ follow the regularity, $\partial_{y}p(y)=-g_{y}\rho(y)$, and
$p(y)=T\rho(y)$ in each part. With fixed $T$, the solution has an
exponentially decaying behavior in the two half volumes. The initial
hydrostatic unstable configuration is therefore given by

$\displaystyle T(y)=T_{u},\;\rho(y)=\rho_{u}\exp(-g_{y}(y-y_{s})/T_{u}),$ (7)
$\displaystyle u_{x}(y)={u_{0}},\;u_{y}(y)=0,$ $\displaystyle\text{for }y\geq
y_{s},$ $\displaystyle
T(y)=T_{b},\;\rho(y)=\rho_{b}\exp(-g_{y}(y-y_{s})/T_{b}),$
$\displaystyle\;u_{x}(y)=-u_{0},\;u_{y}(y)=0,$ $\displaystyle\text{for
}y<y_{s},$

where $\rho_{b}$ and $T_{b}$ ($\rho_{u}$ and $T_{u}$) are the density and
temperature of the light (heavy ) fluid, $u_{0}$ is the shear velocity,
$g_{y}$ is the gravitational acceleration, $y_{s}$ is the initial disturbance
interface. For these six interfaces, the specific expressions of initial
disturbances are as follows:

$y_{s}^{a}=H/2+a_{0}\cos(2\pi x/\lambda),$ (8a)
$y_{s}^{b}=H/2-a_{0}+2a_{0}\frac{|{x-\lambda/2|}}{\lambda/2},$ (8b)
$y_{s}^{c}=H/2-a_{0}+2a_{0}\frac{(x-\lambda/2)^{2}}{(\lambda/2)^{2}},$ (8c)
$y_{s}^{d}=H/2+a_{0}-2a_{0}(1-\frac{(x-\lambda/2)^{2}}{(\lambda/2)^{2}})^{0.5},$
(8d) $y_{s}^{e}=H/2+a_{0}-2a_{0}\frac{(x-\lambda/2)^{2}}{(\lambda/2)^{2}},$
(8e)
$y_{s}^{f}=H/2-a_{0}+2a_{0}(1-\frac{(x-\lambda/2)^{2}}{(\lambda/2)^{2}})^{0.5},$
(8f)

where $\lambda=W$ is the wavelength of the perturbation, $a_{0}$ is the
initial amplitude. To have a finite width of the initial interface, a smooth
interpolation between the two half-volumes is plused on the initial
configuration. When the relative velocity $u_{0}$ is set to zero, this is a
pure RTI system; when the acceleration $g_{y}$ is set to zero, it is a pure
KHI; with various tangential velocities and accelerations, the coupled RTKHI
systems are gotten. The fifth-order weighted essentially non-oscillatory
(WENO) scheme is adopted for space discretization, and the third-order Runge-
Kutta scheme is adopted for time discretization. The left and right boundaries
are periodic boundary conditions. In the bottom and the top boundaries, we
apply the solid wall boundary conditions, i.e.,
$(\rho,U_{x},T)|_{ix,-2}=(\rho,U_{x},T)|_{ix,-1}=(\rho,U_{x},T)|_{ix,0}=(\rho,U_{x},T)|_{ix,1}$,
$U_{y}|_{ix,-2}=U_{y}|_{ix,-1}=U_{y}|_{ix,0}=0$, where $(ix,-2)$, $(ix,-1)$,
and $(ix,0)$ are the indexes of ghost nodes out of the bottom boundary
$(ix,1)$ when the fifth-order WENO scheme is used. The distribution functions
$f_{i}$ on the ghost nodes are approximated by theirs equilibrium distribution
functions $f_{i}^{eq}$. On the top side, we can operate in a similar way.

We first investigate the influence of the initial interface shape on the
evolution of RTI. Figure 3 shows the temperature patterns of the above six
different initial disturbances evolving to $t=300$ under the action of RTI.
Figure 4 shows the time evolutions of bubbles, spikes, disturbance amplitudes
and morphological boundary lengths for different initial interface systems.
The definitions of bubble and spike amplitudes are shown in Fig. 2(b). The
disturbance amplitude of RTI is defined as half of the sum of the bubble
amplitude and the spike amplitude. The morphological boundary length $L$ is
the total length of the interface between high and low temperature (light and
heavy) fluids. The methodology of morphological analysis technique is as
follows: By defining a temperature threshold $T_{th}$, the value higher than
the threshold is defined as the high temperature area (white area), and the
value below the threshold is defined as the low temperature area (black area).
In this way, the original continuous temperature field is transformed into a
temperature pattern of black and white pixels, the total length of the
dividing lines between white and black regions is defined as the morphological
boundary length $L$. The specific calculation method of $L$ can refer to our
previous work Xu2021cjcp ; Xu2009 ; Xu2010sc ; Chen2020 . Both disturbance
amplitude and morphological boundary length are important parameters that
describe the degree of instability development and material mixing. They have
different perspectives, but they are related and cannot replace each other. In
the simulation, the dimensionless parameters are set as follows: $\rho_{b}=1$,
$T_{b}=1.4$, $\rho_{u}=2.33333$, $T_{u}=0.6$, $g_{x}=0$, $g_{y}=0.005$,
$u_{0}=0$, $c=1$, $\eta_{0}=3$, temperature threshold $T_{th}=1.0$, the
specific-heat-ratio $\gamma=1.4$, the initial amplitude $a_{0}=2$, and spatial
step $dx=dy=0.2$, temporal step $dt=10^{-3}$.

It can be seen from Figs. 3 and 4 that the initial disturbance interface shape
has a great influence on the evolution of RTI. In the evolution of RTI, the
heavy and light fluids gradually penetrate into each other as time progresses,
with the light fluid rising to form a bubble and the heavy fluid falling to
generate a spike. The sharper the interface (the interface has a higher
curvature), the faster the growth of bubble or spike. For example, in the
system with an inverted elliptic interface, the heavy fluid falls along both
sides of the computing domain to form spikes. Compared with other interfaces
where spikes are formed, the interface in this place is the sharpest, and the
spike falls at the fastest rate; The light fluid rises in the middle part of
the interface to form bubbles, compared with other interfaces to form bubbles,
the interface here is the bluntest and the bubble rises at the slowest speed.
The difference of amplitudes of RTI systems with different initial interfaces
is small (Fig. 4(c)), but the morphological boundary length differs greatly in
the later stage (Fig. 4(d)), indicating that the initial interface shape has
substantial influence on the mixing degree of different media in the RTI
system.

Figure 3: (Color online) The influence of the initial interface shape on the
evolution of RTI. (a) to (f) are the temperature patterns of the six systems
with different initial disturbances described in Fig. 1 at time $t=300$. The
horizontal and vertical coordinates are the mesh numbers in the computing
domain. For clarity, only half of the computational domains are shown. From
blue to red corresponds to an increase of temperature. Each figure follows the
same legend. Figure 4: (Color online) Effects of initial interface shape on
RTI evolution, (a) bubble amplitude, (b) spike amplitude, (c) total
disturbance amplitude (= half of the sum of bubble amplitude and spike
amplitude), (d) morphological boundary length. Figure 5: (Color online) The
influence of the initial interface shape on the evolution of KHI. (a) to (f)
are the temperature patterns of the six systems with different initial
disturbances described in Fig. 1 at time $t=300$. From blue to red corresponds
to an increase of temperature. Each figure follows the same legend. Figure 6:
(Color online) Effects of initial interface shape on KHI evolution, (a)
amplitude (1/2 width of the mixing layer), (b) morphological boundary length
$L$ of the temperature field.

Figure 5 shows the temperature patterns of the six different initial
disturbances evolving to $t=300$ under the action of KHI. Figure 6 shows the
time evolutions of amplitudes and morphological boundary lengths for different
KHI systems. The amplitude of KHI systems is defined as half of the mixing
width. The initial parameters are $g_{x}=0$, $g_{y}=0$, $u_{0}=0.1$, and the
other parameters are the same as Figure 3. It can be seen from Figs. 5 and 6
that the initial interface shape has little influence on the evolution of KHI.
There is no significant difference in the final degree of material mixing in
KHI systems with different initial interface shapes.

Figure 7: (Color online) Temperature patterns of the RTKHI. From (a) to (f),
the initial perturbations are sinusoidal, sawtooth, parabolic, ellipse,
inverted parabolic and inverted elliptic perturbation interfaces,
respectively. The images from left to right correspond to
$t=0,50,100,150,200$, respectively. Each figure follows the same legend.

Figure 7 shows the evolutions of the coupled RTKHI systems with different
initial disturbances at times $t=0$, $50$, $100$, $150$, and $200$. The
initial conditions are $g_{x}=0$, $g_{y}=0.005$, $u_{0}=0.15$, and the other
parameters remain the same. As shown in the figure, under the action of
initial perturbation and tangential velocity, the perturbations gradually grow
to rolled-up vortexes, which are the main characteristic of free shear flows,
so KHI plays a major role at the initial stage. For the inverted parabolic and
elliptic interfaces, the vortexes are rolled up on the right side of the
computational domain, which is different from the other four interfaces. As
time progresses, more fluid is entrained in the vortical structure, the
secondary RTI develops along the vortex arms, and plays a major role. It
should be pointed out that in the systems with different disturbance
interfaces, the secondary RTI develops at different times and with different
speeds. Among them, the first three systems develop relatively fast, while the
last three develop relatively slowly. In other words, the shape of the
perturbed interface has a certain degree of influence on the evolution of the
coupled RTKHI system. The focal points of the following research are: (i) the
influence of perturbed interface on the transition point of the coupled
system; (ii) the influence of perturbed interface on the main mechanism in the
early stage of the RTKHI system.

Figure 8: (Color online) Influence of the initial interface shapes on the
early mechanism judgment of the coupled system. Figure 9: (Color online)
Influence of the initial interface shapes on the transition points of the
coupled system. From (a) to (f), the initial perturbations are sinusoidal,
sawtooth, parabolic, ellipse, inverted parabolic, and inverted elliptic
perturbation interfaces, respectively.

In the previous work, we have shown that the main mechanism of the coupled
RTKHI system in the early stage depends on the comparison of buoyancy and
shear strength, i.e., gravity acceleration $g$ and shear velocity $u_{0}$, and
both the morphological boundary length $L$ and the mean heat flux strength
$D_{3,1}$ can be used to quantitatively judge the main mechanism in the early
stage. The judgment method is as follows: Comparing the morphological boundary
length $L$ or the non-equilibrium strength $D_{3,1}$ of pure RTI with
acceleration $g$ and KHI systems with shear velocity $u_{0}$, if
$L^{KHI}=L^{RTI}$ ($D_{3,1}^{KHI}=D_{3,1}^{RTI}$), the buoyancy and shear
effects are balanced in the early stages of the corresponding coupled system;
The shear velocity $u_{0}$ is defined as the critical velocity $u_{C}$ of the
coupled system. It stands to reason that, keeping $g$ constant, if the shear
velocity increases, the shear effect is stronger than the buoyancy effect;
when the shear velocity decreases, the buoyancy effect is stronger than the
shear effect.

Figure 8 shows the influence of different initial disturbance interface shapes
on the early mechanism judgment of the coupled system. In the previous work,
we have proved that, the two quantities, $L$ and $D_{3,1}$, always show a high
correlation, especially in the early stage. Therefore, the simulation results
of $D_{3,1}$ are only shown. In figures 8 (a) and (b), the black curve of the
pure RTI system (g = 0.005) and the green curve of the pure KHI system with
shear velocity $u_{0}=0.1$ intersect at a point, and remain approximately
coincidence until this point, i.e., these two systems have the highest similar
degree of medium mixing in the early stages. The critical shear velocity is
$u_{C}=0.1$. In figures 8 (c) and (d), as the interface shape changes, the
bilateral contact angle $\theta_{1}$ increases, the black curve of the pure
RTI system deviates from the green curve, which indicates that the critical
shear velocity $u_{C}$ increases. For the elliptic structure (Fig. 8 (d)), the
critical shear velocity increases to $u_{C}=0.12$. Figures 8 (e) and (f) show
similar trends. The black curve of the pure RTI system with $g=0.005$ is
approximately coincident with the red curve of the pure KHI system with shear
velocity $u_{0}=0.08$ in the early stage, and then gradually deviates from it,
and intersects with the green and blue lines (shear velocity
$u_{0}=0.1,0.12$). In general, for a given gravitational acceleration $g$, the
critical shear velocity $u_{C}$ is related to the shape of the initial
interface.

For the case where the KHI dominates at earlier time and the RTI dominates at
later time, the evolution process can be roughly divided into two stages. In
the previous work, we have pointed out that, before the transition point of
the two stages, both $L$ and $D_{3,1}$ initially increase exponentially, and
then increase linearly. Hence, the ending point of linear increasing $L$ or
$D_{3,1}$ can work as a geometric or physical criterion for discriminating the
two stages.

Figure 9 shows the effects of initial disturbance shapes on the transition
points of the coupled system. The green box denotes the mean heat flux
strength $D_{3,1}$, the red circle represents the $dD_{3,1}/dt$, and the blue
vertical line marks the position of the transition point. As can be seen from
the figure, for the initial disturbance interface whose bilateral contact
angle $\theta_{1}$ is greater than $90^{o}$ (Figs. 9 (b), (c), and (d)), the
delay of transition point is relatively slow and the magnitude is almost
negligible, as the bilateral contact angle increases; For the inverted
parabolic and elliptic interfaces (Figs. 9 (e) and (f)), the transition point
has a relatively rapid and substantial delay, as the bilateral contact angle
decreases; The inverted parabolic and elliptic interfaces (Figs. 9 (e) and
(f)) have a longer linear growth stage, and the transition point is later than
that of the parabolic and elliptic interfaces (Figs. 9 (c) and (d)).

Figure 10: (Color online) Effects of disturbance amplitude on the early main
mechanism of the coupled system. (a), (b), (c), and (d) correspond to
sinusoidal, parabolic, inverted parabolic and sawtooth disturbances,
respectively. Rows 1 to 3 correspond to the initial amplitudes of $2$, $2.8$,
$3.6$, respectively. Figure 11: (Color online) Effects of the initial
disturbance amplitude on the transition point of the coupled system.

Figure 10 shows the effect of disturbance amplitude on the early main
mechanism of the coupled system. From top to bottom, the corresponding initial
disturbance amplitudes are $a_{0}=2$, $2.8$, and $3.6$, respectively. Figures
10 (a), (b), (c), and (d) correspond to sinusoidal, parabolic, inverted
parabolic and sawtooth disturbances, respectively. From the curves we can see
that:

(i) For the sinusoidal disturbance, the contact angles ($\theta_{1}$ and
$\theta_{2}$) of the interface are basically unchanged with the increase of
amplitude, so the change of amplitude has no significant impact on the
judgment of the early mechanism of the RTKHI system (Figs. 10 (a1)-(a3)); For
ellipse and inverted ellipse perturbation interfaces, there are consistent
conclusions.

(ii) For the parabolic disturbance, with the increase of the disturbance
amplitude, the middle contact angle $\theta_{2}$ almost remains unchanged, but
the bilateral contact angle $\theta_{1}$ further increases, the black curve of
pure RTI system with $g=0.005$ gradually approaches the blue curve of pure KHI
system with $u_{0}=0.12$, that is, the critical velocity tends to increase
with the increase of bilateral contact angle $\theta_{1}$ (Figs. 10(b1)-(b3)).
This is consistent with the conclusion of Figs. 8(c) and 8(d).

(iii) For the inverted parabolic disturbance ($\theta_{2}=90^{o}$ and
$\theta_{1}<90^{o}$), the amplitude increases and the bilateral contact angle
$\theta_{1}$ decreases, but it does not affect the judgment of the early
mechanism of the RTKHI system (Figs. 10(c1)-(c3)), and the images are similar
to the system with inverted ellipse interface. This is also consistent with
the conclusion of Figs. 8 (e) and (f).

(iv) For the sawtooth disturbance, the amplitude increases, the bilateral
contact angle $\theta_{1}$ increases, but the middle contact angle
$\theta_{2}$ decreases, and the critical velocity is basically unchanged due
to the combined action of the two factors (Figs. 10 (d1)-(d3)).

Figure 11 shows the effect of the initial disturbance amplitude on the
transition point of the coupled system. The green and red curves correspond to
the cases with $a_{0}=2$ and $a_{0}=3.6$, respectively. The virtual vertical
lines indicate the positions of transition points. It can be seen from Fig. 11
that, (i) For the RTKHI systems with different initial interfaces, the effects
of amplitude on the transition point are not the same; (ii) For the sinusoidal
and sawtooth (Figs. 11 (a) and (b)), the increase of amplitude has no
significant effect on the transition point (marked by blue vertical lines).
(iii) For the parabolic and ellipse interfaces (Figs. 11 (c) and (d)), as the
amplitude increases, the transition points will be brought forward by a
smaller amplitude. The results of $a_{0}=2$ and $a_{0}=3.6$ are marked with
green and red vertical lines, respectively. (iv) For the inverted parabolic
and inverted ellipse perturbation interfaces Figs. 11 (e) and (f)), the
increase of amplitude has a significant effect on the evolution process of the
system and makes the transition point of the coupled system advance
significantly. The effect of amplitude on the inverted ellipse interface is
the greatest.

##  IV Conclusions

In this paper, we design six typical initial disturbance interfaces, and
compare their effects on the RTI, KHI, and the coupled RTKHI systems. It is
found that the initial perturbation has a great influence on the evolution of
RTI. The sharper the interface, the faster the growth of bubbles or spikes.
The influence of initial interface shape on KHI evolution can be ignored.
Based on the mean heat flux strength $D_{3,1}$, the effects of initial
interfaces on the coupled RTKHI system are investigated, and the researches
focus on two aspects: (i) how the perturbed interface affects the transition
point of the coupled system; (ii) how the perturbed interface affects the main
mechanism in the early stage of the RTKHI system.

It is found that, in the coupled RTKHI system, for a given gravitational
acceleration $g$, the magnitude of the critical shear velocity $u_{C}$ is
closely related to the shape of the initial disturbance interface. The
interface shape can be described by two contact angles, the bilateral contact
angle $\theta_{1}$ and the middle contact angle $\theta_{2}$. For the
interface that $\theta_{1}$ and $\theta_{2}$ are basically constant, the
critical velocity is basically constant. The increase of $\theta_{1}$ and the
decrease of $\theta_{2}$ have opposite effects on the critical velocity. For
the parabolic perturbation interface ($\theta_{2}=90^{o}$,
$\theta_{1}>90^{o}$), the critical shear velocity increases as the increase of
bilateral contact angle $\theta_{1}$. The ellipse perturbation interface is
its limiting case. For the inverted parabolic and the inverted ellipse
disturbance ($\theta_{2}=90^{o}$, $\theta_{1}<90^{o}$), the critical shear
velocity is basically the same, which is less than that of the sinusoidal and
sawtooth disturbances. The influence of inverted parabolic and inverted
ellipse perturbations on the transition point of the RTKHI system is greater
than that of other interfaces. When the amplitude is constant, the bilateral
contact angle $\theta_{1}$ of the interface decreases, the transition point
will be greatly delayed, and the transition point of the inverted ellipse
structure system appears at the latest. For the same interface morphology, the
disturbance amplitude increases, resulting in a shorter duration of the linear
growth stage, and the transition point is greatly advanced.

## Acknowledgments

This work was supported by the Natural Science Foundation of Shandong Province
(under Grant Nos. ZR2020MA061, ZR2019PA021), and Shandong Province Higher
Educational Youth Innovation Science and Technology Program (under Grant No.
2019KJJ009), the National Natural Science Foundation of China (under Grant
Nos. 11772064, 11875001, 12102397), CAEP Foundation (under Grant No.
CX2019033), the opening project of State Key Laboratory of Explosion Science
and Technology (Beijing Institute of Technology) (under Grant No. KFJJ21-16M),
the China Postdoctoral Science Foundation (under Grant No. 2019M662521),
Science Foundation of Hebei Province (under Grant No. A2021409001), “Three,
Three and Three Talent Project" of Hebei Province (under Grant No.
A202105005).

## Appendix

$\hat{A}_{l}$ is the $l$th element of
$\hat{\mathbf{A}}=(0,\cdots,0,\hat{A}_{8},\hat{A}_{9},0,\cdots,0)$ and is a
modification to the collision operator
$\hat{S}_{lk}(\hat{f}_{k}-\hat{f}_{k}^{eq})$, where

$\displaystyle\hat{A}_{8}=(s_{T}/s_{v}-1)\rho Tu_{x}\left(2\frac{\partial
u_{x}}{\partial x}-\frac{2}{b}\frac{\partial u_{x}}{\partial
x}-\frac{2}{b}\frac{\partial u_{y}}{\partial y}\right)$
$\displaystyle\qquad+(s_{T}/s_{v}-1)\rho Tu_{y}\left(\frac{\partial
u_{y}}{\partial x}+\frac{\partial u_{x}}{\partial y}\right)\text{,}$
$\displaystyle\hat{A}_{9}=(s_{T}/s_{v}-1)\rho Tu_{x}\left(\frac{\partial
u_{y}}{\partial x}+\frac{\partial u_{x}}{\partial y}\right)$
$\displaystyle\qquad+(s_{T}/s_{v}-1)\rho Tu_{y}\left(2\frac{\partial
u_{y}}{\partial y}-\frac{2}{b}\frac{\partial u_{x}}{\partial
x}-\frac{2}{b}\frac{\partial u_{y}}{\partial y}\right)\text{,}$

$s_{v}=s_{5}=s_{6}=s_{7}$, and $s_{T}=s_{8}=s_{9}$. This is to maintain the
isotropy constraint of viscous stress tensor and heat conductivity.

The transformation matrix and the corresponding equilibrium distribution
functions in kinetic moment space are constructed according to the moment
relations. Specifically, the transformation matrix is

$\mathbf{M}=(m_{1},m_{2},\cdots,m_{16})^{T}\text{,}\;$
$m_{1}=1\text{,}\;m_{2}=v_{ix}\text{,}\;m_{3}=v_{iy}\text{,}\;m_{4}=(v_{i\alpha}^{2}+\eta_{i}^{2})/2\text{,}\;$
$m_{5}=v_{ix}^{2}\text{,}\;m_{6}=v_{ix}v_{iy}\text{,}\;m_{7}=v_{iy}^{2}\text{,}\;$
$m_{8}=(v_{i\beta}^{2}+\eta_{i}^{2})v_{ix}/2\text{,}\;m_{9}=(v_{i\beta}^{2}+\eta_{i}^{2})v_{iy}/2\text{,}\;$
$m_{10}=v_{ix}^{3}\text{,}\;m_{11}=v_{ix}^{2}v_{iy}\text{,}\;m_{12}=v_{ix}v_{iy}^{2}\text{,}\;m_{13}=v_{iy}^{3}\text{,}\;$
$m_{14}=(v_{i\chi}^{2}+\eta_{i}^{2})v_{ix}^{2}/2\text{,}\;m_{15}=(v_{i\chi}^{2}+\eta_{i}^{2})v_{ix}v_{iy}/2\text{,}\;$
$m_{16}=(v_{i\chi}^{2}+\eta_{i}^{2})v_{iy}^{2}/2\text{.}\;$

The corresponding equilibrium distribution functions in KMS are

$\hat{f}_{1}^{eq}=\rho\text{,}\;\hat{f}_{2}^{eq}=\rho
u_{x}\text{,}\;\hat{f}_{3}^{eq}=\rho
u_{y}\text{,}\;\hat{f}_{4}^{eq}=e\text{,}\;$ $\hat{f}_{5}^{eq}=P+\rho
u_{x}^{2}\text{,}\;\hat{f}_{6}^{eq}=\rho
u_{x}u_{y}\text{,}\;\hat{f}_{7}^{eq}=P+\rho u_{y}^{2}\text{,}\;$
$\hat{f}_{8}^{eq}=(e+P)u_{x}\text{,}\;\hat{f}_{9}^{eq}=(e+P)u_{y}\text{,}\;\hat{f}_{10}^{eq}=\rho
u_{x}(3T+u_{x}^{2})\text{,}\;$ $\hat{f}_{11}^{eq}=\rho
u_{y}(T+u_{x}^{2})\text{,}\;\hat{f}_{12}^{eq}=\rho
u_{x}(T+u_{y}^{2})\text{,}\;\hat{f}_{13}^{eq}=\rho
u_{y}(3T+u_{y}^{2})\text{,}\;$
$\hat{f}_{14}^{eq}=(e+P)T+(e+2P)u_{x}^{2}\text{,}\;\hat{f}_{15}^{eq}=(e+2P)u_{x}u_{y}\text{,}$
$\hat{f}_{16}^{eq}=(e+P)T+(e+2P)u_{y}^{2}\text{,}$

where pressure $P=\rho RT$ and energy $e=b\rho RT/2+\rho u_{\alpha}^{2}/2$.
$R$ is the specific gas constant and $b$ is a constant related to the
specific-heat-ratio $\gamma$ by $\gamma=(b+2)/b$.

## References

## References

  * (1) Y. Zhou, “Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. I," Phys Rep, 720, 1-136 (2017).
  * (2) Y. Zhou, “Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. II," Phys Rep, 723, 1-160 (2017).
  * (3) Y. Zhou, T. T. Clark, D. S. Clark, G. S. Gail, S. M. Aaron, C. M. Huntington, O. A. Hurricane, A. M. Dimits, B. A. Remington, “Turbulent mixing and transition criteria of flows induced by hydrodynamic instabilities," Phys. Plasmas, 26(8), 080901 (2019).
  * (4) H. Li, B. Tian, Z. He, Y.Zhang, “Growth mechanism of interfacial fluid-mixing width induced by successive nonlinear wave interactions, " Phys. Rev. E, 103, 053109 (2021).
  * (5) L. F. Wang, C. Xue, W. H. Ye, Y. J. Li, “Destabilizing effect of density gradient on the Kelvin-Helmholtz instability," Phys. Plasmas, 16, 112104 (2009).
  * (6) F. Chen, A. G. Xu, G. C. Zhang, Y. J. Li, S. Succi, “Multiple-relaxation-time lattice Boltzmann approach to compressible flows with flexible specific-heat ratio and Prandtl number," Europhys. Lett., 90: 54003 (2010).
  * (7) H. Liang, Q. X. Li, B. C. Shi, Z. H. Chai, “Lattice Boltzmann simulation of three-dimensional Rayleigh-Taylor instability," Phys. Rev. E, 93, 033113 (2016).
  * (8) H. Liang, X. L. Hu, X. F. Huang, J. R. Xu, “Direct numerical simulations of multi-mode immiscible Rayleigh-Taylor instability with high Reynolds numbers," Phys. Fluids, 31(11), 112104 (2019).
  * (9) H. Liang, Z. H. Xia, H. W Huang, “Late-time description of immiscible Rayleigh-Taylor instability: A lattice Boltzmann study," Phys. Fluids, 33, 082103 (2021).
  * (10) Z. Zhai, L. Zou, Q. Wu and X. Luo, “Review of experimental Richtmyer-Meshkov instability in shock tube: From simple to complex," Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 232(16), 2830 (2018).
  * (11) L. Zou, J. Liu, S. Liao, and X. Zheng,“Richtmyer-Meshkov instability of a flat interface subjected to a rippled shock wave," Phys. Rev. E, 95, 013107 (2017).
  * (12) L. Zou, M. Al-Marouf, W. Cheng, R. Samtaney, J. Ding and X. Luo, “Richtmyer-Meshkov instability of an unperturbed interface subjected to a diffracted convergent shock," J. Fluid Mech., 879, 448 (2019).
  * (13) A. Ravid, R. I. Citron, R. Jeanloz, “Hydrodynamic instability at impact interfaces and planetary implications," Nature Communications, 12, 1 (2021).
  * (14) Y. W. Bin, M. J. Xiao, Y. P. Shi, Y. S. Zhang, S. Y. Chen, “A new idea to predict reshocked Richtmyer-Meshkov mixing: constrained large-eddy simulation," J. Fluid Mech., 918 (2021).
  * (15) H. Y. Ye, H. L. Lai, D. M. Li, Y. B. Gan, C. D. Lin, L. Chen, A. G. Xu, “Knudsen number effects on two-dimensional Rayleigh-Taylor instability in compressible fluid: based on a discrete Boltzmann method," Entropy, 22(5), 500 (2020).
  * (16) L. Chen, H. L. Lai, C. D. Lin, D. M. Li, “Specific heat ratio effects of compressible Rayleigh-Taylor instability studied by discrete Boltzmann method," Front. Phys. 16(5), 52500 (2021).
  * (17) J. G. Tang, F. Zhang, X. S. Luo, Z. G. Zhai, “Effect of Atwood number on convergent Richtmyer-Meshkov instability," Acta Mech. Sin. 37, 434-446 (2021).
  * (18) C. D. Lin, K. H. Luo, Y. B. Gan, Z. P. Liu, “Kinetic simulation of nonequilibrium Kelvin-Helmholtz instability," Commun. Theor. Phys., 71(1), 132-142 (2019).
  * (19) R. H. Zeng, J. J. Tao, Y. B. Sun, “Three-dimensional viscous Rayleigh-Taylor instability at the cylindrical interface," Phys. Rev. E, 102, 023112 (2020).
  * (20) Y. B. Sun, R. H. Zeng, J. J. Tao, “Effects of viscosity and elasticity on rayleigh-taylor instability in a cylindrical geometry," Phys. Plasmas, 28(6), 062701 (2021).
  * (21) G. Dimonte, “Dependence of turbulent Rayleigh-Taylor (RT) instability on initial perturbations," Phys. Rev. E, 69, 056305 (2004).
  * (22) A. R. Miles, M. J. Edwards, and J. A. Greenough, “Effect of initial conditions on two-dimensional Rayleigh-Taylor instability and transition to turbulence in planar blast-wave-driven systems," Phys. Plasmas, 11, 5278 (2004).
  * (23) P. Ramaprabhu, G. Dimonte, M. J. Andrews, “A numerical study of the influence of initial perturbations on the turbulent Rayleigh-Taylor instability," J. Fluid Mech., 536, 285-319 (2005).
  * (24) D. H. Olson, J. W. Jacobs, “Experimental Study of Rayleigh-Taylor Instability With a Complex Initial Perturbation," Phys. Fluids, 21(3), 034103 (2009).
  * (25) A. A. Gowardhan, J. R. Ristorcelli, F. F. Grinstein, “The bipolar behavior of the Richtmyer-Meshkov instability," Phys. Fluids, 23, 071701 (2011).
  * (26) Y. Doron, A. Duggleby, “Optical Density Measurements and Analysis for Single-Mode Initial-Condition Buoyancy-Driven Mixing," J. Fluids Eng., 133(10): 101204 (2011).
  * (27) T. Wei, D. Livescu, “Late-time quadratic growth in single-mode Rayleigh-Taylor instability," Phys. Rev. E, 86, 046405 (2012).
  * (28) S. Kuchibhatl, D. Ranjan, “Effect of initial conditions on Rayleigh-Taylor mixing: modal interaction," Phys. Scr. T, 155, 014057 (2013).
  * (29) W. H. Liu, L. F. Wang, W. H. Ye, X. T. He, “Temporal evolution of bubble tip velocity in classical Rayleigh-Taylor instability at arbitrary Atwood numbers," Phys. Plasmas 20, 062101 (2013).
  * (30) J. A. McFarland, J. A. Greenough, D. Ranjan, “Investigation of the initial perturbation amplitude for the inclined interface Richtmyer-Meshkov instability," Phys. Scr., 2013(T155), 014014 (2013).
  * (31) Z. G. Zhai, M. H. Wang, T. Si, X. S. Luo, “On the interaction of a planar shock with a light polygonal interface, J. Fluid Mech., 757, 800-816 (2014).
  * (32) X. S. Luo, M. H. Wang, T. Si, Z. G. Zhai, “On the interaction of a planar shock with an SF6 polygon, J. Fluid Mech., 773, 366-394 (2015).
  * (33) Z. Dell, R. F. Stellingwerf, S. I. Abarzhi, “Effect of initial perturbation amplitude on Richtmyer-Meshkov flows induced by strong shocks," Phys. Plasmas, 22, 092711 (2015).
  * (34) J. X. Xiao, J. S. Bai, T. Wang, “Numerical study of initial perturbation effects on Richtmyer-Meshkov instability in nonuniform flows," Phys. Rev. E, 94, 013112 (2016).
  * (35) C. Y. Xie, J. J. Tao, Z. L. Sun, J. Li, “Retarding viscous Rayleigh-Taylor mixing by an optimized additional mode," Phys. Rev. E, 95(2), 023109 (2017).
  * (36) A. Kord, J. Capecelatro, “Optimal perturbations for controlling the growth of a Rayleigh-Taylor instability," J. Fluid Mech., 876, 150-185 (2019).
  * (37) R. Sun, J. C. Ding, Z. G. Zhai, T. Si, X. S. Luo, “Convergent Richtmyer-Meshkov instability of heavy gas layer with perturbed inner surface," J. Fluid Mech., 902, A3 (2020).
  * (38) Y. Liang, L. L. Liu, Z. G. Zhai, T. Si, X. S. Luo, “Universal perturbation growth of Richtmyer-Meshkov instability for minimum-surface featured interface induced by weak shock waves," Phys. Fluids, 33, 032110 (2021).
  * (39) L. F. Wang, W. H. Ye, and Y. J. Li, “Combined effect of the density and velocity gradients in the combination of Kelvin-Helmholtz and Rayleigh-Taylor instabilities," Phys. Plasmas 17, 042103 (2010).
  * (40) W. H. Ye, L. F. Wang, C. Xue, Z. F. Fan and X. T. He, “Competitions between Rayleigh-Taylor instability and Kelvin-Helmholtz instability with continuous density and velocity profifiles," Phys. Plasmas 18, 022704 (2011).
  * (41) L. Mandal, S. Roy, R. Banerjee, M. Khan, M. R. Gupta, “Evolution of nonlinear interfacial structure induced by combined effect of Rayleigh-Taylor and Kelvin-Helmholtz instability," Nucl. Instr. Meth. in Phys. Res. A 653, 103-108 (2011).
  * (42) B. J. Olson, J. Larsson, S. K. Lele, and A. W. Cook, “Nonlinear effects in the combined Rayleigh-Taylor/Kelvin-Helmholtz instability," Phys. Fluids 23, 114107 (2011).
  * (43) B. Akula, M. J. Andrews, D. Ranjan, “Effect of shear on Rayleigh-Taylor mixing at small Atwood number," Phys. Rev. E, 87(3): 033013 (2013).
  * (44) M. Vadivukkarasan and M. V. Panchagnula, “Helical modes in combined Rayleigh-Taylor and Kelvin-Helmholtz instability of a cylindrical interface ," Int. J. Spray Combust. 8, 219-234 (2016).
  * (45) M. Vadivukkarasan and M. V. Panchagnula, “Combined Rayleigh-Taylor and Kelvin-Helmholtz instabilities on an annular liquid sheet," J. Fluid Mech. 812, 152-177 (2017).
  * (46) M. Vadivukkarasan, “Temporal instability characteristics of Rayleigh-Taylor and Kelvin-Helmholtz mechanisms of an inviscid cylindrical interface," Meccanica, 56, 117-124 (2021).
  * (47) V. D. Sarychev, S. A. Nevskii, A. Y. Granovskii, S. V. Konovalov, and V. E. Gromov, “Combined Rayleigh-Taylor-Kelvin-Helmholtz instability and its role in the formation of the surface relief of the coating/substrate," AIP Conf. Proc. 2167, 020307 (2019).
  * (48) S. Brizzolara, J. Mollicone, M. Van Reeuwijk, A. Mazzino, M. Holzner, “Transition from shear-dominated to Rayleigh-Taylor turbulence," J. Fluid Mech., 924, A10, (2021).
  * (49) F. Chen, A. G. Xu, Y. D. Zhang, Q. K. Zeng, “Morphological and non-equilibrium analysis of coupled Rayleigh-Taylor-Kelvin-Helmholtz instability," Phys. Fluids 32, 104111 (2020).
  * (50) Without causing misunderstanding, DBM is used as an abbreviation of discrete Boltzmann model/Modeling/Method.
  * (51) A. G. Xu, G. C. Zhang, Y. B. Gan, F. Chen, X. Yu, “Lattice Boltzmann modeling and simulation of compressible flows," Front. Phys. 7, 582-600 (2012).
  * (52) A. G. Xu, G. C. Zhang, Y. J. Ying, “Progess of discrete Boltzmann modeling and simulation of combustion system," Acta. Phys. Sin. 64, 184701 (2015).
  * (53) A. Xu, G. Zhang, Y. Gan, “Progress in studies on discrete Boltzmann modeling of phase separation process," Mech. Eng. 38, 361-374 (2016).
  * (54) A. G. Xu, G. C. Zhang, Y. D. Zhang, _Discrete Boltzmann Modeling of Compressible Flows, Chapter 2 in Kinetic Theory_ , edited by G. Z. Kyzas and A. C. Mitropoulos (Rijeka: InTech, 2018).
  * (55) A. G. Xu, J. Chen, J. H. Song, D. W. Chen, Z. H. Chen, “Progress of discrete Boltzmann study on multiphase complex flows," Acta Aerodynamica Sinica, 39(3): 138-169 (2021).
  * (56) A. G. Xu, J. H. Song, F. Chen, K. Xie, Y. J. YING, “Modeling and analysis methods for complex fields based on phase space," Chinese J. Comput. Phys., (2021-05-25). https://kns.cnki.net/kcms/detail/11.2011.O4.20210524.1535.002.html
  * (57) A. G. Xu, Y. M. Shan, F. Chen, Y. B. Gan, C. D. Lin,“Progress of m esoscale modeling and investigation of combustion multiphase flow," Acta Aeronautica et Astronautica Sinica, 2021, 42(12), 625842 (2021).
  * (58) S. Succi, _The Lattice Boltzmann Equation for Fluid Dynamics and Beyond_ (New York: Oxford University Press, 2001)
  * (59) R. Benzi, S. Succi, and M. Vergassola, “The lattice Boltzmann equation: Theory and applications," Phys. Rep. 222, 145-197 (1992).
  * (60) X. Shan, and H. Chen, “Lattice Boltzmann model for simulating flows with multiple phases and components," Phys. Rev. E 47(3), 1815 (1993).
  * (61) Y. Zhang, R. Qin, and D. Emerson, “Lattice Boltzmann simulation of rarefied gas flows in microchannels," Phys. Rev. E 71, 047702 (2005).
  * (62) V. E. Ambruş, V. Sofonea, _Quadrature-Based Lattice Boltzmann Models for Rarefied Gas Flow_ , edited by F. Toschi, M. Sega (Springer, 2019)
  * (63) Y. B. Li, X. W. Shan, “Lattice Boltzmann method for adiabatic acoustics," Phil. Trans. R. Soc. A, 369, 2371-2380 (2011).
  * (64) Q. Li, K. H. Luo, Y. J. Gao, and Y. L. He, “Additional interfacial force in lattice Boltzmann models for incompressible multiphase flows," Phys. Rev. E 85(2), 026704 (2012).
  * (65) Z. Wang, Y. Wei, Y. Qian, “A simple direct heating thermal immersed boundary-lattice Boltzmann method for its application in incompressible flow," Comput. Math. Appl. 80, 1633-1649 (2020).
  * (66) Z. Chen, C. Shu, D. Tan, “Highly accurate simplified lattice Boltzmann method," Phys. Fluids 30, 103605 (2018).
  * (67) F. B. Tian, H. Luo, L. Zhu, J. C. Liao, X. Y. Lu, “An efficient immersed boundary-lattice Boltzmann method for the hydrodynamic interaction of elastic filaments," J. Comput. Phys. 230, 7266-7283 (2011).
  * (68) F. B. Tian, Y. Wang, H. Liu, Y. Zhang, “The lattice Boltzmann method and its applications in complex flows and fluid-structure interactions," Inst. Mech. Eng. C J. Mech. Eng. Sci. 232, 403-404 (2018).
  * (69) H. Liang, B. C. Shi, Z. L. Guo, Z. H. Chai, “Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows," Phys. Rev. E, 89, 053320 (2014).
  * (70) Y. Wang, C. Zhong, C. Zhuo, S. Liu, “A simplified finite volume lattice Boltzmann method for simulations of fluid flows from laminar to turbulent regime, Part I: Numerical framework and its application to laminar flow simulation," Comput. Math. Appl. 79, 1590-1618 (2020).
  * (71) K. Pasieczynski, B. X. Chen, “Multipseudopotential interaction models for thermal lattice Boltzmann method simulations," Phys. Rev. E, 102, 013311 (2020).
  * (72) R. Qiu, Y. Bao, T. Zhou, H. Che, R. Chen, and Y. You, “Study of regular reflection shock waves using a mesoscopic kinetic approach: Curvature pattern and effects of viscosity," Phys. Fluids, 32, 106106 (2020).
  * (73) R. Qiu, T. Zhou, Y. Bao, K. Zhou, H. Che, and Y. You, “Mesoscopic kinetic approach for studying nonequilibrium hydrodynamic and thermodynamic effects of shock wave, contact discontinuity, and rarefaction wave in the unsteady shock tube," Phys. Rev. E, 103, 053113 (2021).
  * (74) D. K. Sun, “A discrete kinetic scheme to model anisotropic liquid-solid phase transitions," Appl. Math. Lett., 103, 106222 (2020).
  * (75) D. K. Sun, H. Xing, X. L. Dong, Y. S. Han, “An anisotropic lattice Boltzmann - Phase field scheme for numerical simulations of dendritic growth with melt convection," Int. J. Heat Mass Tran., 133, 1240-1250 (2019).
  * (76) C. J. Zhan, Z. H. Chai, B. C. Shi, “A lattice Boltzmann model for the coupled cross-diffusion-fluid system," Appl. Math. Comput., 400, 126105, (2021).
  * (77) A. Xu, G. Zhang, X. Pan, P. Zhang, J. Zhu, “Morphological characterization of shocked porous material," J. Phys. D: Appl. Phys. 42, 075409 (2009).
  * (78) A. G. Xu, G. C. Zhang, H. Li, Y. J. Ying, X. J. Yu, J. S. Zhu, “Temperature pattern dynamics in shocked porous materials," Sci. China-Phys. Mech. Astron., 53(8): 1466-1474 (2010).
  * (79) A. G. Xu, G. C. Zhang, Y. J. Ying, C. Wang, “Complex fields in heterogeneous materials under shock: modeling, simulation and analysis," Sci. China-Phys. Mech. Astron. 59, 650501 (2016).
  * (80) F. Chen, A. Xu, G. Zhang, “Viscosity, heat conductivity, and Prandtl number effects in the Rayleigh Taylor Instability," Front. Phys. 11(6), 114703 (2016).
  * (81) F. Chen, A. Xu, G. Zhang, “Collaboration and competition between Richtmyer-Meshkov instability and Rayleigh-Taylor instability," Phys. Fluids, 30(10), 102105 (2018).


# Semi-Supervised Translation with MMD Networks

Mark Hamilton

###### Abstract

This work aims to improve semi-supervised learning in a neural network
architecture by introducing a hybrid supervised and unsupervised cost
function. The unsupervised component is trained using a differentiable
estimator of the Maximum Mean Discrepancy (MMD) distance between the network
output and the target dataset. We introduce the notion of an $n$-channel
network and several methods to improve performance of these nets based on
supervised pre-initialization, and multi-scale kernels. This work investigates
the effectiveness of these methods on language translation where very few
quality translations are known a priori. We also present a thorough
investigation of the hyper-parameter space of this method on both synthetic
data.

Semi-supervised Learning, Maximum Mean Discrepency, MMD, Translation, SMT,
Regression

##  1 Introduction

Often in data analysis, one has a small set of quality labeled data, and a
large pool of unlabeled data. It is the task of semi-supervised learning to
make as much use of this unlabeled data as possible. In the low-data regime,
the aim is to create models that perform well after seeing only a handful of
labeled examples. This is often the case with machine translation and
dictionary completion, as it can be difficult to construct a large number of
labeled instances or a sufficiently large parallel corpora. However, this
domain offers a huge number of monolingual corpora to make high quality
language embeddings (Tiedemann, 2012; Al-Rfou et al., 2013). The methods
presented in this paper are designed to take into consideration both labeled
and unlabeled information when training a neural network. The supervised
component uses the standard alignment-based loss functions and the
unsupervised component attempts to match the distribution of the network’s
output to the target data’s distribution by minimizing the Maximum Mean
Discrepancy (MMD) “distance” between the two distributions. This has the
effect of placing a prior on translation methods that preserve the
distributional structure of the two datasets. This limits the model space and
increases the quality of the mapping, allowing one to use less labeled data.

Related methods such as Auto-Encoder pre-initialization (Erhan et al., 2010),
first learn the structure of the input, then learn a mapping. In this setup,
unsupervised knowledge enters through learning good features to describe the
dataset. The MMD method of unsupervised training directly learns a mapping
between the two spaces that aligns all of the moments of the mapped data and
the target data. This method can be used to improve any semi-supervised
mapping problem, such as mappings between languages (Dinu et al., 2014), image
labeling, FMRI analysis (Mitchell et al., 2008), and any other domains where
transformations need to be learned between data. This investigation aims to
study these methods in the low data regime, with the eventual goal of studying
of dying or lost languages, where very few supervised training examples exist.

##  2 Background

###  2.1 Maximum Mean Discrepancy

The Maximum Mean Discrepancy (MMD) put forth by (Gretton et al., 2012a) is a
measure of distance between two distributions $p,q$. More formally, letting
$x$, $y$ be variables defined on a topological space $\mathcal{X}$ with Borel
measures $p,q$, and $\mathcal{F}$ be a class of functions from
$\mathcal{X}\to\mathbb{R}$. The MMD semi-metric is defined as:

$MMD_{\mathcal{F}}\left(p,q\right)=sup_{f\in\mathcal{F}}\Bigl{(}\mathbb{E}_{x\sim
p}f\left(x\right)-\mathbb{E}_{y\sim q}f\left(y\right)\Bigr{)}$ (1)

Where $\mathbb{E}$ is the first raw moment defined as:

$\mathbb{E}_{x\sim p}f(x)=\int_{\mathcal{X}}f(x)dp$ (2)

Intuitively, the MMD is a measure of distance which uses a class of functions
as a collection of “trials” to put the two distributions through. The
distributions pass a trial if the function evaluated on both distributions has
the same expectation or mean. Two distributions fail a trial if they yield
different means, the size of the difference measures how much the
distributions fail that trial. Identical distributions should yield the same
images when put through each function in $\mathcal{F}$, so the means (first
moments) of the images should also be identical. Conversely, if the function
class is “large enough” this method can distinguish between any two
probability distributions that differ, making the MMD a semi-metric on the
space of probability distributions. A unit ball in a Reproducing Kernel
Hilbert Space (RKHS) is sufficient to discern any two distributions provided
the kernel, $k$, is universal. (Cortes et al., 2008) If $\mathcal{F}$ is equal
to a unit ball in kernel space, Gretton $\it{et.al.}$ showed that the
following is an unbiased estimator of the MMD: (Gretton et al., 2012a)

$MMD_{u}^{2}(X,Y)=\frac{1}{m(m-1)}\sum_{i=1}^{m}\sum_{j\neq
i}^{m}k(x_{i},x_{j})+\\\ \frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\neq
i}^{n}k(y_{i},y_{j})-\frac{2}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}k(x_{i},y_{j})$
(3)

If the kernel function is differentiable, this implies that the estimator of
the MMD is differentiable, allowing one to use it as a loss function that can
be optimized with gradient descent.

###  2.2 MMD Networks

The differentiability of the MMD estimator allows it to be used as a loss
function in a feed-forward network. Li $\it{et.al.}$ showed that by using the
MMD distance as a loss function in a neural net, $\mathcal{N}$, one can learn
a transformation that maps a distribution of points $X=(x_{i})_{1}^{n}$ in
$\mathbb{R}^{d}$ to another distribution $Y=(y_{i})_{1}^{m}$ in
$\mathbb{R}^{n}$ while approximately minimizing the MMD distance between the
image of $X$, $\mathcal{N}(X)$, and $Y$. (Li et al., 2015)

$l_{MMD}(X,Y,\mathcal{N})=MMD_{u}^{2}(\mathcal{N}(X),Y)$ (4)

This loss function allows the net to learn transformations of probability
distributions in a completely unsupervised manner. Furthermore, the MMD-net
can also be used to create generative models, or mappings from a simple
distribution to a target distribution.(Li et al., 2015) Where simple usually
means easy to sample from, or a maximum entropy distribution. Often, a
multivariate uniform or Gaussian source distribution is used in these
generative models. This loss function can be optimized via mini-batch
stochastic gradient descent, though the samples from X and Y need not be
paired in any way. To avoid over-fitting, the minibatches for X and Y should
be sampled independently, which this paper refers to as “unpaired”
minibatching.

##  3 Methods

###  3.1 $n$-Channel Networks

This work introduces a generalization of a feed forward net, called an
$n$-Channel net. This architecture allows an unsupervised loss term that
requires unpaired mini-batching and a paired mini-batching scheme of a
standard feed forward network to be mixed.

An $n$-channel net is a collection of $n$ networks with tied weights that
operate on $n$ separate datasets $(X_{i},Y_{i})_{1}^{n}$. More formally, an
$n$-channel net is a mapping:

$\mathcal{N}_{n}:\left(\mathbb{R}^{d}\right)^{n}\to\left(\mathbb{R}^{e}\right)^{n}$
(5)

defined as:

$\mathcal{N}_{n}\Bigl{(}\left(X_{i}\right)_{1}^{n}\Bigr{)}\equiv\Bigl{(}\mathcal{N}(X_{i})\Bigr{)}_{1}^{n}$
(6)

where where $\mathcal{N}:\mathbb{R}^{d}\to\mathbb{R}^{e}$ is a feed forward
network. Each channel of the network can have it’s own loss function and be
fed with a separate data source. Most importantly, these separate data sources
can be trained in a paired or unpaired manner.

###  3.2 A Semi-Supervised MMD-Net

In many applications where one is interested in estimating a transformation
between data spaces, one has a small labeled dataset $(X,Y)$, and large,
unlabeled datasets $(S,T)$. Throughout the literature, MMD networks have only
been applied to the case of unpaired data.(Li et al., 2015) We expand on this
work by augmenting the completely unsupervised MMD distance with a semi-
supervised alignment term. More formally, if one has a collection of k paired
vectors $(x_{i},y_{i})_{1}^{k}$ with $x_{i}\in X$ and $y_{i}\in Y$ that should
be aligned through the transformation $\mathcal{N}$, one can use the standard
loss function:

$l_{alignment}(X,Y,\mathcal{N})=\sum_{i=1}^{k}\lVert\mathcal{N}(x_{i})-y_{i}\rVert$
(7)

Where $\lVert\cdot\rVert$ is any differentiable norm in $\mathbb{R}^{d}$. This
work uses the standard $l_{2}$ vector norm. This is the standard norm used in
regression, where the goal of the network is to minimize the distance between
the network output$\mathcal{N}(x_{i})$, and the observed responses $y_{i}$.

Using a hyperparameter, we can blend the cost functions of the supervised
alignment loss and the unsupervised MMD loss. The full cost function for the
MMD network then becomes:

$l(X,Y,S,T,\mathcal{N})=\alpha_{pair}l_{alignment}(X,Y,\mathcal{N})+\\\
(1-\alpha_{pair})l_{MMD}(S,T,\mathcal{N})$ (8)

###  3.3 Supervised Pre-Initialization

The MMD term of the cost function scales as $\mathcal{O}\left(M^{2}\right)$
where $M$ is the size of the mini-batch. This significantly increases training
time for large batch sizes slowing convergence in wall-time. To mitigate this
effect, we first train the network until convergence with only the supervised
term of the cost function. Once converged, we then switch to the semi-
supervised cost function.

This also helps the network avoid local minima as it already starts close to
the optimal solution. Because the MMD cost function is inherently unpaired, it
is susceptible to getting stuck in local minima when there are multiple ways
to map the mass of one probability distribution into another distribution. We
say that a mapping from the supports, $f:\mathcal{X}\to\mathcal{Y}$, is a MMD-
mode from distributions $p$ to $q$ if $f(p)\sim q$. Here $f(p)$ is the
distribution formed by sampling from $p$ and then applying $f$. These modes
coincide with critical points of the $MMD_{u}^{2}$ cost function and are
therefore tough to escape with gradient descent methods. As the class of
functions represented by the network increases, the more distinct MMD-modes
arise. This increases the number of critical points, though these probably
tend to be saddle points rather than local minima as the dimensionality of the
function space increases. (Dauphin et al., 2014)

One can escape these local minima, by increasing $\alpha_{pair}$ to the point
where the signal from the supervised term overcomes that the signal from the
unsupervised cost function. However, if the network is within the pull of the
correct minima, it is often better to rely on the robust unsupervised signal
than the noisy supervised signal, which requires a small $\alpha_{pair}$. We
found that supervised pre-training helped guide the network parameters to
within the basin of attraction for the correct unsupervised minima. From here
the unsupervised signal was much more reliable and led to better results on
synthetic and language datasets. Furthermore on all datasets, the supervised
warm-start greatly reduced fitting time, as convergence of the expensive MMD
cost function needed fewer optimization steps. Future work could involve
annealing the supervised term to a small number, though this would eliminate
the aforementioned computational speedup.

To demonstrate the effect of pre-initialization, we show the unbiased MMD
estimator of a simple synthetic experiment. We generate two datasets of two
dimensional points. The first, shown in Figure 1 left is sampled from a
uniform distribution on the unit square support centered at $(0,0)$. To
generate a simple target shown in Figure 1 middle, we rotate the source cloud
of points by an angle $\theta^{*}=255^{\circ}$ and add a small Gaussian noise
term. Figure 1 right shows the that MMD loss as a function of angle of
rotation transformation has several modes caused by the symmetries of the
square. To simulate a very noisy MSE, we use the MSE of one randomly sampled
point and its respective pair. The noisy MSE loss function has two local
minima and the global minima $\hat{\theta}$ is within the correct basin of
attraction of the unsupervised cost function. This basin of attraction of the
unsupervised cost has a minima that is indistinguishable from the correct
value of theta and much more accurate than the supervised loss term.

Figure 1: Left: Initial dataset $X$ sampled uniformly from the unit square.
Colors indicate how points are mapped through the transform. Middle:
$Y=X_{255^{\circ}}+Gaussian(\mu=0,\sigma=.1)$ Where $X_{\theta}$ denotes a
rotation clockwise by $\theta$. Right: Unit scaled
$MMD_{u}^{2}(X_{\theta},Y)$, and unit scaled $MSE(X_{\theta,1},Y_{1})$ as a
function of $\theta$. Where $X_{1}$ denotes the first element of $X$.

###  3.4 Choice of Kernel

The $MMD_{\mathcal{F}}$ is able to differentiate between any two distributions
if the function class, $\mathcal{F}$, is a unit ball in the reproducing kernel
Hilbert space (RKHS) of a universal kernel.(Cortes et al., 2008) One of the
simplest and most commonly used universal kernels is the Gaussian or radial
basis function kernel, which excels at representing smooth functions.

$k_{\sigma}(x,y)=exp\left(-\frac{\lVert x-y\rVert^{2}}{2\sigma^{2}}\right)$
(9)

The parameter $\sigma$ controls the width of the Gaussian, and needs to be set
properly for good performance. If $\sigma$ is too low, each point’s local
neighborhood will be effectively empty, and the gradients will vanish. If it
is too high, every point will be in each point’s local neighborhood and the
kernel will not have enough resolution to see the details of the distribution.
In this scenario, the gradients vanish. We found that $\sigma$ was one of the
most important hyper-parameters for the success of the method. In both our
synthetic data and natural language examples, we found that the method
performed well in a small window of kernel scale settings.

To improve the robustness of this method, this investigation used the
following multi-scale Gaussian kernel:

$k(x,y)=\sum_{i=0}^{n}c_{i}k_{\sigma_{i}}(x,y)$

Where $c_{i}=1$, $\sigma_{i}=s10^{w(i/n)-w/2}$, $w=4$, $n=10$. The scalar $s$
is the average scale of the multi-scale and the width, $w$, controls the width
of the frequency range covered by the kernel. $n$ controls how many samples
are taken from this range. Choosing a larger $n$ improves performance as there
are more scales in the kernel, but increases computation time. By including
multiple scales in the kernel, the gradients from the larger kernels will move
the parameters to a region where the distributions are aligned at a large
scale, they will then begin to vanish and the smaller scale gradients will
become more relevant. Setting $w=4$ allows the kernel to be sensitive to
functions with scales that are within $2$ orders of magnitude of the average
scale $s$. We find that choosing this kernel significantly broadens the areas
of parameter space where the method succeeds, without hurting the performance.

Many have investigated the kernel scale problem and there are several
heuristics available for choosing the scale based on optimizing statistical
power or median distances to nearest neighbors. (Gretton et al., 2012b) For
clarity, we explicitly investigated and set the kernel scale based on a grid
search evaluating on a held out validation set. Figure 2 demonstrates that the
method was fairly robust to settings of average kernel scale on synthetic data
and language data.

Figure 2: Left: Performance comparison on word embeddings in the $0-5k$
frequency bin as a function of the average kernel scale $s$. Middle:
Performance comparison on synthetically generated data in $\mathbb{R}^{30}$ as
a function of $\alpha_{pair}$. Right: Performance comparison on synthetically
generated data in $\mathbb{R}^{300}$ as a function of $\alpha_{pair}$.

###  3.5 Globally Corrected (GC) Retrieval

In this analysis, performance of translation methods are compared on their
ability to infer the correct translation on a held out test set. More
specifically, we use the precision at $N$, which is the fraction of examples
where the correct word was in the top $N$ most likely translations of model.
This is a natural choice for translation, as it estimates the probability of
translating a word correctly when $N=1$.

To generate the list of $N$ most likely translations for a given word, one can
use nearest neighbor (NN) retrieval. In this method, one uses the $N$ closest
neighbors in the target space of the mapped word vector as the list of best
guesses. We find that it is always better to use cosine distance for nearest
neighbor calculations. Finding the first nearest neighbor to a point $\hat{y}$
can be more formally expressed as:

$NN_{1}(\hat{y})=argmin_{y\in T}Rank_{T}(\hat{y},y)$ (10)

Where $\hat{y}$ is our mapped word vector, $T$ is our target space, and
$Rank_{T}(\hat{y},y)$ is a function that returns the rank of $y$ in the sorted
list of distances between $\hat{y}$ and the points in $T$.

If the space of word embeddings is not uniformly distributed, there will be
areas where word embeddings bunch together in higher densities. The points
towards the center of these bunches act as hub points, and may be the nearest
neighbors of many other points. Dinu et. al. 2014 have shown that naive NN
retrieval results in over-weighting these hub points as they are more
frequently the neighbors of points. They called this the “Hubness Problem” and
introduced a corrected form of the nearest-neighbor retrieval called the
globally corrected neighbor retrieval method (GC). In this method, instead of
using distance to select translates as in $NN_{1}$, one uses:

$GC_{1}(\hat{y})=argmin_{y\in
T}\left(Rank_{P}\left(y,\hat{y}\right)-cos\left(\hat{y},y\right)\right)$ (11)

Where $P$ is a random sampling of points from $T$ and $cos(x,y)$ is the cosine
distance between $x$ and $y$. Instead of returning the nearest neighbor of
$\hat{y}$, GC returns the point in $T$ that has $\hat{y}$ ranked the highest.
The cosine distance term breaks ties. GC retrieval has been shown to
outperform the nearest neighbor retrieval in all frequency bins when the
transformation is a linear mapping.(Dinu et al., 2014) Figure 4 shows that it
also improves the performance of the semi-supervised translation task.

###  3.6 Neural Network Implementation

This work implemented the network in Theano, (Bergstra et al., 2011) an
automatic differentiation software written in python. The net was trained with
RMSProp (Tieleman & Hinton, 2012) on both the unpaired and paired batches with
a batch size of 200 for each set. The unregularized pre-initialization was
trained for $4000$ epochs and the regularized network was trained for $250$
epochs, which gave ample time for convergence. Hyperparameter optimization was
perfomed through parallel grid searches a TORQUE Cluster, where each job ran
for $\sim 20$ hours. A validation set consisting of a random sample of $10\%$
of the training set was used to choose the parameters for the final reported
results.

##  4 Data

###  4.1 Synthetic Data

Several synthetic datasets were used to demonstrate the method’s ability to
accurately learn linear transformations using a very small paired dataset.
Furthermore, we used this synthetic data to investigate the effects of the
network’s hyper-parameters.

Two datasets were created, one with the dimension of the source and target
equal to $30$ and the other $300$, the same dimensionality as the embeddings.
The datasets contained $100,000$ points and various sized paired subsets were
used to calculate the supervised alignment loss in the experiments.

Source data was generated as a multivariate Gaussian with zero mean and unit
variance. A ground truth mapping was generated by sampling the entries of a
$d\times d$ matrix of independent Gaussians with zero mean and unit variance.
The target data was generated by applying the ground truth transformation to
the source data and adding Gaussian noise with zero mean and a variance of
$0.1$.

###  4.2 Embedding Data

This analysis used 300 dimensional English (EN) and Italian (IT) monolingual
word embeddings from (Dinu et al., 2014). These embeddings were trained with
word2vec’s CBOW model on $2.8$ billion tokens as input (ukWaC + Wikipedia +
BNC) for English and the 1.6 billion itWaC tokens for Italian.(Dinu et al.,
2014) The embeddings contained the top 200,000 words in each language.
Supervised training and testing sets were constructed from a dictionary built
from Europarl, available at http://opus.lingfil.uu.se/. (Tiedemann, 2012)

Two training sets consisted of the $750$ and $5,000$ most frequent words from
the source language (English) which had translations in the gold dictionary.
Five disjoint test sets were created consisting of roughly 400 translation
pairs randomly sampled from the frequency ranked words in the intervals 0-5k,
5k-20k, 20k-50k, 100k-200k.

##  5 Results

###  5.1 Synthetic Data

Adding the MMD term to the loss function dramatically improved the ability to
learn the transformation on all synthetic datasets. The synthetic data also
provided a clean environment to see the effect of varying hyper-parameters.
The experiment used a “linear network” which is equivalent to learning a
linear transformation between the spaces. In general, if the hyper-parameters
are set correctly, the MMD assisted learner can approach the true
transformation with significantly less paired data.

Our first investigation aimed to understand the effect and robustness of the
kernel scale parameter. As one can see from Figure 2, the performance of the
method is robust to a setting of the average kernel scale within $+/-2$ orders
of magnitude of the optimal scale. This empirically confirms the intuition
behind the width parameter of the multi-scale kernel. As the width parameter
decreases, this valley of good performance becomes narrower by the expected
amount. A similar pattern arose in the $300$ dimensional dataset.

In order to simulate the environment of the embedding experiment that required
a validation set of $\sim 10\%$ of the data, we also removed $\sim 10\%$ of
our data. The plots in Figure 2 demonstrate that even with the data removed
for a validation set, the method still significantly beats linear regression
trained on the training and validation set, justifying the use of data for
parameter tuning. The models in $d=30$ and $d=300$ both reach error rates
comparable to the ground truth regressor learned on all $100,000$ data points.

Figure 3: Performance of methods on synthetically generated data in
$\mathbb{R}^{300}$ as a function of $\alpha_{pair}$, $s=10$.

Figure 3 investigates various settings of $\alpha_{pair}$ and shows that
decreasing $\alpha_{pair}$ drives the performance down to the ground truth
level. This trend appears in both the low and high dimensional data and
suggests that the supervised pre-initialization yields a configuration that is
within the basin of attraction of the true parameters in vector field $\nabla
l_{MMD}$. Thus, only the unsupervised term is needed as the supervised
initialization has already eliminated the ambiguity of the MMD loss function
modes.

###  5.2 Embedding Data

Figure 4 shows that the semi-supervised MMD-Net was able to significantly
outperform the standard linear regression on a paired dataset of $750$ and
$5000$ word-translation pairs in every frequency bin . Furthermore, this
dominance over linear regression follows a similar pattern in the precisions
@5 and @10. The method also outperformed several other linear and nonlinear
methods as shown in Table 1.

Figure 4: Model performance as a function of English word frequency bins using
the top 5000 (left) and 750 (right) EN-IT word pairs as training data.
Precision@1 refers to the fraction of words correctly translated by the method
on held out testing sets. Table 1: Comparison of Precision@1 across different
algorithms and dimensionality reduction schemes. PCA S and PCA T refers to
projecting the source and target respectively onto their first 270 principal
vectors. KR refers to Kernel Ridge Regression an RBF refers to the radial
basis function kernel with heuristically set scale
| 0-5k | 5k-20k | 20k-50k | 50k-100k | 100k-200k
---|---|---|---|---|---
Linear | 0.228 | 0.052 | 0.028 | 0.015 | 0.011
Linear + PCA S | 0.236 | 0.057 | 0.031 | 0.036 | 0.019
Linear + PCA T | 0.207 | 0.044 | 0.031 | 0.028 | 0.011
Linear + PCA S + T | 0.212 | 0.072 | 0.033 | 0.043 | 0.029
Random Forrest | 0.008 | 0.000 | 0.000 | 0.000 | 0.000
KR 2-deg Poly | 0.057 | 0.003 | 0.008 | 0.010 | 0.008
KR 3-deg Poly | 0.049 | 0.005 | 0.003 | 0.013 | 0.008
KR RBF | 0.057 | 0.003 | 0.010 | 0.010 | 0.008
Linear + MMD | 0.347 | 0.129 | 0.099 | 0.094 | 0.035

##  6 Discussion and Future Work

The addition of the MMD cost function term significantly improves the results
of regression in the low data regime. Furthermore, to the best knowledge of
the authors, this method achieves state of the art results on the embeddings
of (Dinu et al., 2014). The authors also experimented with deeper nets, but
did not observe significant performance improvements, an observation
consistent with the observations of (Mikolov et al., 2013).

###  6.1 Adversarial Distribution Matching

One promising future direction involves replacing the MMD unsupervised term
with a Generative Adversarial Network (GAN) (Goodfellow et al., 2014). Like
the MMD, the GAN also involves a maximization over a function class of a
measure of dissimilarity. Similarly, the GAN loss function can be used for
unsupervised learning of probability distributions. However, the GAN is
usually optimized directly by stochastic gradient descent, trading the
quadratic time dependence on minibatch size with a linear one. In practice
however, the maximization over the function class (the discriminator) is
usually done in $k$ gradient descent steps for every one step of training the
distribution matching net (the generator). Furthermore, the GAN cost function
does not have a dependence on kernel scale.

Analogous to the discriminator in the GAN, we can also adversarially learn the
MMD. In this setup, the function class takes the the form of a parametrized
network. Instead of estimating the supremum of the mean discrepancy over a
ball in RKHS, we would be finding the supremum through gradient ascent on the
network. This would also have the effect of eliminating the quadratic compute
and the dependence on kernel scale. This formulation of the MMD would allow
for a more direct comparison between the GAN and MMD loss functions, and
warrants future investigation. These two loss functions are in-equivalent, as
the only intersection between $f$-divergences, like the Jensen-Shannon
Divergence which is equivalent to the GAN, and integral measures like the MMD
is the total variation distance. (Mohamed & Lakshminarayanan, 2016) Thus, one
might be able to leverage more diverse information by combining the two.

###  6.2 Bi-Directional Networks

In the case of translation between two spaces of equal dimension, the inverse
of the translation transformation should also be a translation from the target
to the source space. We can capitalize on this observation to further
constrain our set of possible translations. This allows the transformation to
also draw information from the structure of the source space. More
specifically one can minimize:

$L=\alpha_{target}\|RT-S\|_{target}^{2}+\\\
(1-\alpha_{target})\|R-ST^{-1}\|_{source}^{2}$ (12)

where $T\in GL_{d}$, $\alpha_{target}\in[0,1]$ and $R,S\in\mathbb{R}^{d\times
n_{pair}}$. This would result in twice as much supervisory signal and maintain
the same number of parameters. Furthermore, this can also be applied in
conjunction with the GAN loss. It is also compatible with the pre-
initialization scheme. In the case of a more complex nonlinear network where
an inverse transformation cannot be easily calculated, the architecture could
include an encoder network which maps from the source to the target and a
decoding network which maps from the target to the source. These two mappings
could then be constrained to be close to mutual inverses through a
reconstruction loss penalty.

## References

  * Al-Rfou et al. (2013) Al-Rfou, Rami, Perozzi, Bryan, and Skiena, Steven.  Polyglot: Distributed word representations for multilingual nlp.  _arXiv preprint arXiv:1307.1662_ , 2013.
  * Bergstra et al. (2011) Bergstra, James, Bastien, Frédéric, Breuleux, Olivier, Lamblin, Pascal, Pascanu, Razvan, Delalleau, Olivier, Desjardins, Guillaume, Warde-Farley, David, Goodfellow, Ian, Bergeron, Arnaud, et al.  Theano: Deep learning on gpus with python.  In _NIPS 2011, BigLearning Workshop, Granada, Spain_ , volume 3. Citeseer, 2011.
  * Cortes et al. (2008) Cortes, Corinna, Mohri, Mehryar, Riley, Michael, and Rostamizadeh, Afshin.  Sample selection bias correction theory.  In _International Conference on Algorithmic Learning Theory_ , pp. 38–53. Springer, 2008.
  * Dauphin et al. (2014) Dauphin, Yann N, Pascanu, Razvan, Gulcehre, Caglar, Cho, Kyunghyun, Ganguli, Surya, and Bengio, Yoshua.  Identifying and attacking the saddle point problem in high-dimensional non-convex optimization.  In _Advances in neural information processing systems_ , pp. 2933–2941, 2014.
  * Dinu et al. (2014) Dinu, Georgiana, Lazaridou, Angeliki, and Baroni, Marco.  Improving zero-shot learning by mitigating the hubness problem.  _arXiv preprint arXiv:1412.6568_ , 2014.
  * Erhan et al. (2010) Erhan, Dumitru, Bengio, Yoshua, Courville, Aaron, Manzagol, Pierre-Antoine, Vincent, Pascal, and Bengio, Samy.  Why does unsupervised pre-training help deep learning?  _The Journal of Machine Learning Research_ , 11:625–660, 2010.
  * Goodfellow et al. (2014) Goodfellow, Ian, Pouget-Abadie, Jean, Mirza, Mehdi, Xu, Bing, Warde-Farley, David, Ozair, Sherjil, Courville, Aaron, and Bengio, Yoshua.  Generative adversarial nets.  In _Advances in neural information processing systems_ , pp. 2672–2680, 2014.
  * Gretton et al. (2012a) Gretton, Arthur, Borgwardt, Karsten M, Rasch, Malte J, Schölkopf, Bernhard, and Smola, Alexander.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012a.
  * Gretton et al. (2012b) Gretton, Arthur, Sejdinovic, Dino, Strathmann, Heiko, Balakrishnan, Sivaraman, Pontil, Massimiliano, Fukumizu, Kenji, and Sriperumbudur, Bharath K.  Optimal kernel choice for large-scale two-sample tests.  In _Advances in neural information processing systems_ , pp. 1205–1213, 2012b.
  * Li et al. (2015) Li, Yujia, Swersky, Kevin, and Zemel, Richard.  Generative moment matching networks.  _arXiv preprint arXiv:1502.02761_ , 2015.
  * Mikolov et al. (2013) Mikolov, Tomas, Le, Quoc V, and Sutskever, Ilya.  Exploiting similarities among languages for machine translation.  _arXiv preprint arXiv:1309.4168_ , 2013.
  * Mitchell et al. (2008) Mitchell, Tom M, Shinkareva, Svetlana V, Carlson, Andrew, Chang, Kai-Min, Malave, Vicente L, Mason, Robert A, and Just, Marcel Adam.  Predicting human brain activity associated with the meanings of nouns.  _science_ , 320(5880):1191–1195, 2008.
  * Mohamed & Lakshminarayanan (2016) Mohamed, Shakir and Lakshminarayanan, Balaji.  Learning in implicit generative models.  _arXiv preprint arXiv:1610.03483_ , 2016.
  * Tiedemann (2012) Tiedemann, Jörg.  Parallel data, tools and interfaces in opus.  2012\.
  * Tieleman & Hinton (2012) Tieleman, Tijmen and Hinton, Geoffrey.  Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude.  _COURSERA: Neural networks for machine learning_ , 4(2), 2012.