#  Asymptotic stability of viscous contact wave and rarefaction waves for the
system of heat-conductive ideal gas without viscosity

Lili Fan  Lili Fan: School of Mathematics and Computer Science, Wuhan
Polytechnic University, Wuhan 430023, P. R. China<EMAIL_ADDRESS>,  Guiqiong
Gong  Guiqiong Gong: School of Mathematics and Statistics, Wuhan University,
Wuhan 430072, P. R. China<EMAIL_ADDRESS>and  Shaojun Tang  Shaojun
Tang: School of Mathematics and Statistics, Wuhan University, Wuhan 430072, P.
R. China<EMAIL_ADDRESS>
*Corresponding author. Email address<EMAIL_ADDRESS>(S. J. Tang)

This paper is concerned with the Cauchy problem of heat-conductive ideal gas
without viscosity. We show that, for the non-viscous case, if the strengths of
the wave patterns and the initial perturbation are suitably small, the unique
global-in-time solution exists and asymptotically tends toward the
corresponding the viscous contact wave or the composition of a viscous contact
wave with rarefaction waves determined by the initial condition, which
extended the results by Huang-Li-Matsumura [13], where they treated the
viscous and heat-conductive ideal gas.

Keywords: Non-viscous; asymptotic behavior; viscous contact wave; rarefaction
waves.

Mathematics Subject Classification 2010: 35Q35, 35B35, 35L65

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries and Main Results
    1. 2.1 Viscous Contact Wave
    2. 2.2 Composition Waves
  3. 3 Stability Analysis
    1. 3.1 Reformed the System
    2. 3.2 A Priori Estimates

##  1\. Introduction

We consider the Cauchy problem for the equation of heat-conductive ideal gas
without viscosity

$\displaystyle\left\\{\begin{array}[]{ll}v_{t}-u_{x}=0,\\\\[5.69054pt]
u_{t}+p_{x}=0,\\\\[5.69054pt]
(e+\frac{u^{2}}{2})_{t}+(pu)_{x}=\kappa\big{(}\frac{\theta_{x}}{v}\big{)}_{x}\,,\quad\quad
x\in{\mathbb{R}},\ t>0\end{array}\right.$ (1.4)

with the following prescribed initial data and the far field state

$\displaystyle\left\\{\begin{array}[]{l}(v,u,\theta)(x,0)=(v_{0},u_{0},\theta_{0})(x),\quad
x\in{\mathbb{R}}\,,\\\\[5.69054pt]
(v,u,\theta)(\pm\infty,t)=(v_{\pm},u_{\pm},\theta_{\pm}),\quad
t>0\,.\end{array}\right.$ (1.7)

Here $v(x,t)>0$, $u(x,t)$, $\theta(x,t)>0$, $e(x,t)>0$ and $p(x,t)$ are the
specific volume, fluid velocity, internal energy, absolute temperature and
pressure of the gas, respectively, while $\kappa>0$ denotes the heat
conduction coefficient. $v_{\pm}(>0)$, $\theta_{\pm}(>0)$ and
$u_{\pm}(\in{\mathbb{R}})$ are given constants and the initial data
$(v_{0}(x),u_{0}(x),\theta_{0}(x))$ are assumed to satisfy
$\inf\limits_{x\in{\mathbb{R}}}v_{0}(x)>0$,
$\inf\limits_{x\in{\mathbb{R}}}\theta_{0}(x)>0$ and the compatibility
conditions
$(v_{0},u_{0},\theta_{0})(\pm\infty)=(v_{\pm},u_{\pm},\theta_{\pm})$.

Throughout this paper, we are concerned with the ideal and polytropic fluids
and in such a case, $p$ and $e$ are given by the following state equations

$p=\frac{R\theta}{v}=Av^{-\gamma}e^{\frac{\gamma-1}{R}s}\,,\quad
e=\frac{R}{\gamma-1}\theta+const.\,,$ (1.8)

where $s$ is the entropy, $\gamma>1$ is the adiabatic exponent and both $A$
and $R$ are positive constants.

If $\kappa=0$, we can rewrite the system (1.4) as

$\displaystyle\left\\{\begin{array}[]{ll}v_{t}-u_{x}=0,\\\\[5.69054pt]
u_{t}+p_{x}=0,\\\\[5.69054pt]
(e+\frac{u^{2}}{2})_{t}+(pu)_{x}=0.\end{array}\right.$ (1.12)

It is well-known that the above system has the three eigenvalues:
$\lambda_{1}=-\sqrt{\gamma p/v}<0,\ \lambda_{2}=0,\
\lambda_{3}=-\lambda_{1}>0$, where the second characteristic field is linear
degenerate and the others are genuinely nonlinear. The solutions of this
inviscous equations with the Riemann initial data

$\displaystyle(v,u,\theta)(x,0)=\left\\{\begin{array}[]{l}(v_{-},u_{-},\theta_{-}),x<0\,,\\\\[5.69054pt]
(v_{+},u_{+},\theta_{+}),x>0\,.\end{array}\right.$ (1.15)

has the basic Riemann solutions which are dilation invariant: shock wave,
rarefaction wave, contact discontinuity, and certain linear combinations of
these basic wave patterns [39]. The inviscid system (1.12) is a typical
example of hyperbolic conservation laws and is of great importance to study
the large-time asymptotic behavior of solutions of the corresponding viscous
system (1.4). In fact, if the unique global entropy to its Riemann problem
(1.12), (1.15) consists of shock wave, rarefaction wave, contact
discontinuity, and/or their linear superpositions, then the large time
behavior of the Cauchy problem (1.4), (1.7) of the corresponding viscous
conservation laws is expected to be precisely described by the corresponding
viscous shock wave, rarefaction wave, viscous contact discontinuity and/or
their linear superpositions.

The rigorous mathematical justifications of the above expectation is one of
the hottest topics in the field of nonlinear partial differential equations,
especially in the field of nonlinear dissipative hyperbolic conservation laws.
In fact, new phenomena have been discovered and new techniques, such as the
weighted characteristic energy method and the approximate Green function
method have been developed based on the intrinsic properties of the underlying
nonlinear waves, cf. [13, 22, 23, 25, 27, 40, 41] and the references cited
therein.

By combining these methods together with the fundamental energy method, some
interesting results have been obtained for the compressible Navier-Stokes
equations:

  * •

For the case when both viscosity coefficient $\mu$ and the heat conductivity
coefficient $\kappa$ are positive constants, many excellent results have been
obtained for the nonlinear stability of some basic wave patterns consisting of
diffusion waves, viscous shock waves, rarefaction waves, viscous contact
discontinuities and/or their certain linear superpositions with small
perturbation, cf. [25] for diffusion waves, [23, 22, 20, 31] for viscous shock
profiles, [24, 31] for rarefaction waves, [15, 17, 18] for viscous contact
discontinuities, [14] for the composition of viscous shock profiles of
different family, [13] for the composition a viscous contact wave and
rarefaction waves. For the corresponding results with large initial
perturbation, see [5, 7, 16, 19, 32, 33, 34, 36, 38] and the reference cited
therein.

  * •

For the case when $\mu>0$, $\kappa=0$, Liu and Zeng [26] studied the large-
time behavior of solutions around a constant state for this case. For the
three-dimensional case, the global existence and the temporal convergence rate
of solutions was obtained by Duan and Ma [6].

  * •

When $\mu=0$, $\kappa>0$, if the corresponding Riemann solution consists of
two shock waves form different families, the nonlinear stability of the
composition of viscous shock waves from different families was investigated by
Fan and Matsumura in [8], while the nonlinear stability of viscous contact
waves together with the temporal decay rates was obtained by Ma and Wang in
[29]. In all these results, both the strengths of the underlying wave patterns
and the initial perturbation are assumed to be small.

From the above results, a natural question is: _Whether can we get the
nonlinear stability of combination of viscous contact wave with rarefaction
waves for the case $\mu=0$ and $\kappa>0?$_ This is the motivation of our
work. As a continuation of [8], this manuscript showed that, if the strengths
of the viscous waves and the initial perturbation are suitably small, the
unique global-in-time solution exists and asymptotically tends toward the
corresponding viscous contact wave or the composition of a viscous contact
wave with rarefaction waves from different families. Our result generalizes
the corresponding results obtained by Huang, Li and Matsumura in [13] for the
case of $\mu>0,\kappa>0$ to the case of $\mu>0,\kappa=0$ and extends the
result of Ma and Wang [29] for the nonlinear stability of the viscous contact
waves to the nonlinear stability of the composition of a viscous contact wave
with rarefaction waves from different families.

Now we outline the main difficulties of the problem and our strategy to
overcome the difficulties. The first difficulty is due to the fact that the
system (1.4) is of less dissipation caused by the fact that $\mu=0$, thus we
need more subtle estimates to recover the regularity and dissipation for the
components of the hyperbolic part. We shall overcome this difficulty by
manipulating several new energy estimates and also looking for the perturbed
solution for the integrated system of (1.4) in $C([0,\infty),H^{2})$ to
control the nonlinearity of the hyperbolic part, instead of the usual
$C([0,\infty),H^{1})$ in [13]. Secondly, our stability results included the
superposition of the rarefaction waves with the viscous contact discontinuity,
which will lead to the degenerate characteristics. Thanks to the estimates on
the heat kernel function obtained by Huang, Li and Matsumura established in
[13], we can still close the energy type estimates for our non-viscous case
and get the desired results.

Before concluding this section, it is worth pointing out that there are many
results on the nonlinear stability of basic wave patterns to some hyperbolic
conservation laws with dissipation, cf. [9, 27, 40] for hyperbolic
conservation laws with artificial viscosity, [10, 21, 11, 12, 35, 37] for
compressible Navier-Stokes equations with density and/or temperature dependent
transportation coefficients, and [1, 2, 3, 4] for compressible fluid models of
Korteweg type, and so on.

The rest of the paper is arranged in the following way: in the next section,
we will give some elementary properties of the viscous contact wave and
rarefaction wave and state the main results. Main theorem will be proved in
the section 3.

##  2\. Preliminaries and Main Results

To show our main results, in this section, we will construct the two desired
viscous contact wave and viscous rarefaction waves for (1.4) and state the
main results. For each $z_{-}:=(v_{-},u_{-},\theta_{-})$, we can see our
situation takes place provided $z_{+}:=(v_{+},u_{+},\theta_{+})$ is located on
a quarter of a curved surface in a small neighborhood of $z_{-}$. In what
follows, as [8], the neighborhood of $z_{-}$ denoted by $\Omega_{-}$ is given
by

$\Omega_{-}=\\{(v,u,\theta)|\
|(v-v_{-},u-u_{-},\theta-\theta_{-})|\leq\bar{\delta}\\},$

where $\bar{\delta}$ is a positive constant depending only on $z_{-}$.

###  2.1. Viscous Contact Wave

we firstly recall the viscous contact wave
$(\widetilde{v},\widetilde{u},\widetilde{\theta})$ for the compressible system
(1.4) defined in [18]. For the Riemann problem (1.12),(1.15), it is known that
the contact discontinuity solution
$\widetilde{Z}(x,t):=(\widetilde{V},\widetilde{U},\widetilde{\Theta})(x,t)$
takes the form

$\displaystyle(\widetilde{V},\widetilde{U},\widetilde{\Theta})(x,t)=\left\\{\begin{array}[]{l}(v_{-},u_{-},\theta_{-}),\
x<0,\ t>0\,.\\\\[5.69054pt] (v_{+},u_{+},\theta_{+}),\ x>0,\
t>0\,.\end{array}\right.$ (2.3)

provided that

$\displaystyle u_{-}=u_{+},\quad
p_{-}=\frac{R\theta_{-}}{v_{-}}=\frac{R\theta_{+}}{v_{+}}=p_{+}\,.$ (2.4)

In the setting of the compressible Navier-Stokes system (1.4), the smooth
approximate wave $(\widetilde{v},\widetilde{u},\widetilde{\theta})$ to the
contact wave behaves as a diffusion wave due to the dissipation effect and we
call this wave ”viscous contact wave”. it can be constructed as follows. Since
the pressure is constant under the condition (2.11), we set

$\displaystyle p_{+}=\frac{R\widetilde{\theta}}{\widetilde{v}}\,,$ (2.5)

which indicates the leading part of the energy equation $(\ref{1.1})_{3}$ is

$\displaystyle\tfrac{R}{\gamma-1}\widetilde{\theta}_{t}+p_{+}\widetilde{u}_{x}=\kappa\Big{(}\tfrac{\widetilde{\theta}_{x}}{\widetilde{v}}\Big{)}_{x}\,.$
(2.6)

Meanwhile the equation $\widetilde{v}_{t}=\widetilde{u}_{x}$ leads to a
nonlinear diffusion equation

$\displaystyle\left\\{\begin{array}[]{l}\widetilde{\theta}_{t}=a\Big{(}\frac{\widetilde{\theta}_{x}}{\widetilde{\theta}}\Big{)}_{x},\quad
a=\tfrac{\kappa p_{+}(\gamma-1)}{\gamma R^{2}},\\\\[5.69054pt]
\Theta(\pm,t)=\theta_{\pm}\,.\end{array}\right.$ (2.9)

which has a unique self similarity solution
$\widetilde{\theta}(x,t)=\widetilde{\theta}(\xi),\ \xi=\frac{x}{\sqrt{1+t}}$.
Furthermore, on the one hand, $\widetilde{\theta}(\xi)$ is a monotone
function, increasing if $\theta_{+}>\theta_{-}$ and decreasing if
$\theta_{+}<\theta_{-}$; on the other hand, there exists some positive
constant $\delta_{0}$, such that for
$\delta=|\theta_{+}-\theta_{-}|\leq\delta_{0}(\leq\overline{\delta})$,
$\widetilde{\theta}$ satisfies

$\displaystyle(1+t)|\widetilde{\theta}_{xx}|+(1+t)^{\frac{1}{2}}|\widetilde{\theta}_{x}|+|\widetilde{\theta}-\theta_{\pm}|\lesssim\delta
e^{-\frac{x^{2}}{1+t}},\quad as\quad|x|\rightarrow\infty\,.$ (2.10)

Once $\widetilde{\theta}$ is defined, the viscous contact profile

$Z^{c}(x,t):=(V^{c},U^{c},\Theta^{c})(x,t)$

is determined as follows:

$\displaystyle\Theta^{c}=\widetilde{\theta},\quad
V^{c}=\frac{R}{p_{+}}\widetilde{\theta}\,,\quad
U^{c}=u_{-}+\frac{\kappa(\gamma-1)}{\gamma
R}\frac{\widetilde{\theta}_{x}}{\widetilde{\theta}}\,.$ (2.11)

It is straightforward to check that $Z^{c}(x,t)$ satisfies

$\displaystyle\|(Z^{c}-\widetilde{Z})(t)\|_{L^{p}}=O(\kappa^{\frac{1}{2p}})(1+t)^{\frac{1}{2p}},\
p\geq 1\,,$

which means the nonlinear diffusion wave $Z^{c}(x,t)$ approximates the contact
discontinuity $\widetilde{Z}(x,t)$ to the Euler system (1.12) in $L^{p}(p\geq
1)$ norm. Moreover, the viscous contact wave $Z^{c}(x,t)$ solves the
compressible Navier-Stokes system without viscosity (1.4) as

$\displaystyle\left\\{\begin{array}[]{ll}V^{c}_{t}-U^{c}_{x}=0\,,\\\\[5.69054pt]
U^{c}_{t}+\Big{(}\frac{R\Theta^{c}}{V^{c}}\Big{)}_{x}=U^{c}_{t}\,,\\\\[5.69054pt]
\tfrac{R}{\gamma-1}\Theta^{c}_{t}+p_{+}U^{c}_{x}=\kappa\Big{(}\frac{\Theta^{c}_{x}}{V^{c}}\Big{)}_{x}\,.\end{array}\right.$
(2.15)

Our first main result is as follow:

######  Theorem 2.1.

For any given $z_{-}$, assume that $z_{+}\in\Omega_{-}$ satisfies (2.4), let
$Z^{c}(x,t)$ is the viscous contact wave defined in (2.11) with strength
$\delta=|\theta_{+}-\theta_{-}|$. There exist positive constants
$\epsilon_{1}$ and $\delta_{1}(\leq\delta_{0})$, such that if
$\delta<\delta_{1}$ and the initial data satisfies

$\displaystyle\|(v_{0}(\cdot)-V(\cdot,0),u_{0}(\cdot)-U(\cdot,0),\theta_{0}(\cdot)-\Theta(\cdot,0))\|_{2}\leq\epsilon_{1},$
(2.16)

then the Cauchy problem (1.4), (1.8) admits a unique global solution
$(v,u,\theta)(t,x)$ satisfies

$\displaystyle\begin{aligned} (v-V^{c},u-U^{c},\theta-\Theta^{c})(t,x)\in
X[0,+\infty)\end{aligned}$

and

$\displaystyle\lim\limits_{t\rightarrow\infty}\sup\limits_{x\in{\mathbb{R}}}|(v-V^{c},u-U^{c},\theta-\Theta^{c})(x,t)|=0\,,$
(2.17)

where the solution space $X(0,t)$ is defined in (3.11).

###  2.2. Composition Waves

When the relation (2.4) is fails, the basic theory of conservation laws [39]
implies that for any given constant state $z_{-}$ and $z_{+}\in\Omega_{-}$ (
$\overline{\delta}$ is suitably small ), the Riemann problem (1.12), (1.15)
has a unique solution. Based on this, our second purpose is concerned with the
stability of superposition of a viscous contact wave with rarefaction waves.
In this situation, we assume that

$\displaystyle z_{+}\in R_{1}CR_{3}(z_{-})\subseteq\Omega_{-}\,,$ (2.18)

where

$\displaystyle
R_{1}CR_{3}(z_{-}):=\bigg{\\{}(v,u,\theta)\in\Omega_{-}\bigg{|}s\neq s_{-}$
(2.19) $\displaystyle u\geq
u_{-}-\int^{e^{\frac{\gamma-1}{R\gamma}(s-s_{-})v}}_{v_{-}}\lambda_{1}(\eta,s_{-})d\eta,\quad
u\geq
u_{-}-\int^{v_{-}}_{e^{\frac{\gamma-1}{R\gamma}(s-s_{-})v}}\lambda_{3}(\eta,s)d\eta\bigg{\\}}$

with the entropy $s$ in (1.8) is defined as follows:

$\displaystyle s=\frac{R}{\gamma-1}\ln\frac{R\theta}{A}+R\ln v\,,\quad
s_{\pm}=\frac{R}{\gamma-1}\ln\frac{R\theta_{\pm}}{A}+R\ln v_{\pm}\,.$ (2.20)

It is known that there exists some suitably small $\delta_{1}>0$ such that for

$|\theta_{-}-\theta_{+}|\leq\delta_{1},$

there exists a unique pair of points
$z^{m}_{-}:=(v^{m}_{-},u^{m},\theta^{m}_{-})$ and
$z^{m}_{+}:=(v^{m}_{+},u^{m},\theta^{m}_{+})$ in $\Omega_{-}$ satisfying

$\displaystyle\frac{R\theta^{m}_{-}}{v^{m}_{-}}=\frac{R\theta^{m}_{+}}{v^{m}_{+}}:=p_{m}$
(2.21)

and

$\displaystyle|v^{m}_{\pm}-v_{\pm}|+|u^{m}-u_{\pm}|+|\theta^{m}_{\pm}-\theta_{\pm}|\lesssim|\theta_{+}-\theta_{-}|\,.$
(2.22)

Moreover, the point $z^{m}_{-}$ belongs to the 1-rarefaction wave curve
$z^{r}_{1}:=(v^{r}_{1},u^{r}_{1},\theta^{r}_{1})(\frac{x}{t})$ which connected
with $z_{-}$, while $z^{m}_{+}$ belongs to the 3-rarefaction wave curve
$z^{r}_{3}:=(v^{r}_{3},u^{r}_{3},\theta^{r}_{3})(\frac{x}{t})$ which connected
with $z_{+}$. that is, the 1-rarefaction wave is the weak solution of Riemann
problem of the Euler system (1.12) with the following Riemann data

$\displaystyle
z^{r}_{1}(x,0)=\left\\{\begin{array}[]{l}(v_{-},u_{-},\theta_{-}),\
x<0,\\\\[5.69054pt] (v^{m}_{-},u^{m},\theta^{m}_{-}),\ x>0\end{array}\right.$

and the 3-rarefaction wave with Riemann data as

$\displaystyle
z^{r}_{3}(x,0)=\left\\{\begin{array}[]{l}(v^{m}_{+},u^{m},\theta^{m}_{+}),\
x<0\,,\\\\[5.69054pt] (v_{+},u_{+},\theta_{+}),\ x>0\,.\end{array}\right.$

Since rarefaction waves $z^{r}_{i}(\frac{x}{t})(i=1,3)$ are not smooth enough,
it is convenient to construct smooth approximations of them. Motivated by
[20], we start with the problem of the Burgers equation:

$\displaystyle\left\\{\begin{array}[]{l}w^{r}_{t}+w^{r}w^{r}_{x}=0,\quad
x\in{\mathbb{R}},\ t>0,\\\\[5.69054pt]
w^{r}(0,x)=w^{r}_{0}(x):=\tfrac{1}{2}\big{(}w_{r}+w_{l}\big{)}+\tfrac{1}{2}\big{(}w_{r}-w_{l}\big{)}\tanh(x),\end{array}\right.$
(2.27)

where $w_{l}=\lambda_{1}(v_{-},s_{-})$, $w_{r}=\lambda_{1}(v_{+},s_{-})$. Let
$w(x,t)$ be the unique global solution of (2.27), then
$Z^{r}_{1}(x,t):=(V^{r}_{1},U^{r}_{1},\Theta^{r}_{1})(x,t)$ are defined by

$\displaystyle\left\\{\begin{array}[]{ll}\lambda_{1}(V^{r}_{1},s_{-})=w(x,t),\\\\[5.69054pt]
U^{r}_{1}=u_{-}-\int^{V^{r}_{1}}_{v_{-}}\lambda_{1}(\eta,s_{-})d\eta,\\\\[5.69054pt]
\Theta^{r}_{1}=\theta_{-}(v_{-})^{\gamma-1}(V^{r}_{1})^{1-\gamma}\end{array}\right.$
(2.31)

with

$w_{l}=\lambda_{1}(v_{-},s_{-}),\quad w_{r}=\lambda_{1}(v^{m}_{-},s_{-})$

and the smooth approximation of
$(v^{r}_{3},u^{r}_{3},\theta^{r}_{3})(\frac{x}{t})$ is given by
$Z^{r}_{3}(x,t):=(V^{r}_{3},U^{r}_{3},\Theta^{r}_{3})(x,t)$ constructed by the
same manner as (2.31)

$\displaystyle\left\\{\begin{array}[]{ll}\lambda_{3}(V^{r}_{3},s_{+})=w(x,t),\\\\[5.69054pt]
U^{r}_{3}=u_{+}-\int^{V^{r}_{3}}_{v_{+}}\lambda_{3}(\eta,s_{+})d\eta,\\\\[5.69054pt]
\Theta^{r}_{3}=\theta_{+}(v_{+})^{\gamma-1}(V^{r}_{3})^{1-\gamma},\end{array}\right.$
(2.35)

where $w(x,t)$ is defined in (2.27) with

$w_{l}=\lambda_{3}(v^{m}_{+},s_{+}),\quad w_{r}=\lambda_{3}(v_{+},s_{+}).$

And since the condition (2.21), (2.22), $z^{m}_{\pm}$ is connected by the
viscous contact wave $Z^{c}(x,t)$ constructed in (2.9), (2.11). Finally, we
investigate some properties of $Z^{r}_{i}(x,t)(i=1,3)$ and $Z^{c}(x,t)$ as
[13], divided ${\mathbb{R}}\times[0,t]$ into three parts
${\mathbb{R}}\times[0,t]=\Omega_{1}\cup\Omega_{c}\cup\Omega_{3}$ with

$\displaystyle\begin{aligned}
&\Omega_{1}=\big{\\{}(x,t)|2x<\lambda_{1}(v_{-}^{m},s_{-})t\big{\\}},\\\\[5.69054pt]
&\Omega_{3}=\big{\\{}(x,t)|2x>\lambda_{3}(v_{+}^{m},s_{+})t\big{\\}},\\\\[5.69054pt]
&\Omega_{c}=\big{\\{}(x,t)|\lambda_{1}(v_{-}^{m},s_{-})t\leq
2x\leq\lambda_{3}(v_{+}^{m},s_{+})t\big{\\}},\end{aligned}$ (2.36)

it holds that

######  Lemma 2.1.

([13]) For any given $z_{-}$, assume that $z_{+}\in
R_{1}CR_{3}(z_{-})\subset\Omega_{-}$, then the smooth rarefaction wave
$Z^{r}_{i}(x,t)(i=1,3)$ and $Z^{c}$ satisfy:
$(i)$ $V^{r}_{it}=U^{r}_{ix}>0(i=1,3)$ for all $x\in{\mathbb{R}},\ t>0$.
$(ii)$ For $1\leq p\leq\infty$ and $i=1,3$, it holds

$\displaystyle\|(V^{r}_{ix},U^{r}_{ix},\Theta^{r}_{ix})(t)\|_{L^{p}}\lesssim\min\\{\delta,\delta^{\frac{1}{p}}(1+t)^{-1+\frac{1}{p}}\\},$
$\displaystyle\|(V^{r}_{ixx},U^{r}_{ixx},\Theta^{r}_{ixx})(t)\|_{L^{p}}\lesssim\min\\{\delta,(1+t)^{-1}\\}.$

$(iii)$ In $\Omega_{c}$, we have for $i=1,3$

$\displaystyle|Z^{r}_{ix}|+|Z^{r}_{1}-z^{m}_{-}|+|Z^{r}_{3}-z^{m}_{+}|\lesssim\delta
e^{-c(|x|+t)}$

and in $\Omega_{i}(i=1,3)$

$\displaystyle|Z^{c}_{x}|+|Z^{c}-z^{m}_{\pm}|\lesssim\delta e^{-c(|x|+t)},$
$\displaystyle|Z^{r}_{ix}|+|Z^{r}_{1}-z^{m}_{-}|+|Z^{r}_{3}-z^{m}_{+}|\lesssim\delta
e^{-c(|x|+t)}.$

$(iv)$ For the rarefaction wave $z^{r}_{i}\big{(}\frac{x}{t}\big{)}(i=1,3)$,
it holds

$\displaystyle\lim\limits_{t\rightarrow\infty}\sup\limits_{x\in{\bf{R}}}|Z^{r}_{i}(x,t)-z^{r}_{i}(\tfrac{x}{t})|=0.$

Setting $Z(x,t):=(V,U,\Theta)(x,t)$ is

$\displaystyle\left\\{\begin{array}[]{ll}V(x,t)=V^{r}_{1}(x,t)+V^{c}(x,t)+V^{r}_{3}(x,t)-v^{m}_{-}-v^{m}_{+},\\\\[5.69054pt]
U(x,t)=U^{r}_{1}(x,t)+U^{c}(x,t)+U^{r}_{3}(x,t)-2u^{m},\\\\[5.69054pt]
\Theta(x,t)=\Theta^{r}_{1}(x,t)+\Theta^{c}(x,t)+\Theta^{r}_{3}(x,t)-\theta^{m}_{-}-\theta^{m}_{+}.\end{array}\right.$
(2.40)

Then our second main result is as follow:

######  Theorem 2.2.

For any given $z_{-}$, assume that $z_{+}\in
R_{1}CR_{3}(z_{-})\subseteq\Omega_{-}$ with
$|\theta_{+}-\theta_{-}|\leq\delta_{1}$. There exist positive constants
$\epsilon_{2}$ and $\delta_{2}(\leq\min\\{\overline{\delta},\delta_{1}\\})$,
such that if $\delta<\delta_{2}$ and the initial data satisfies

$\displaystyle\|(v_{0}(\cdot)-V(\cdot,0),u_{0}(\cdot)-U(\cdot,0),\theta_{0}(\cdot)-\Theta(\cdot,0))\|_{2}\leq\epsilon_{2},$
(2.41)

then the Cauchy problem (1.4), (1.8) admits a unique global solution
$(v,u,\theta)(x,t)$ satisfying

$(v-V,u-U,\theta-\Theta)\in X([0,+\infty))$

and

$\displaystyle\lim\limits_{t\rightarrow\infty}\sup\limits_{x\in{\mathbb{R}}}\left\\{\begin{array}[]{ll}|(v-v^{r}_{1}-V^{c}-v^{r}_{3}+v^{m}_{-}+v^{m}_{+})(x,t)|\\\\[5.69054pt]
|(u-u^{r}_{1}-U^{c}-u^{r}_{3}+2u^{m})(x,t)|\\\\[5.69054pt]
|(\theta-\theta^{r}_{1}-\Theta^{c}-\theta^{r}_{3}+\theta^{m}_{-}+\theta^{m}_{+})(x,t)|\end{array}\right\\}=0$
(2.45)

where the solution space $X(0,t)$ is defined in (3.11).

######  Remark 2.1.

if $z^{m}_{-}=z^{m}_{+}=(v_{m},u_{m},\theta_{m})$, we can also get the
stability of the composition of two rarefaction waves.

##  3\. Stability Analysis

In this section, we will prove the asymptotic behavior of the solution for
nonvisous system (1.4), (1.7). If $z^{m}_{\pm}=z_{\pm}$, Theorem 2.2 will
coincide with the result of Theorem 2.1, here we will show the stability of
the composition wave only.

###  3.1. Reformed the System

$Z(x,t):=(V,U,\Theta)(x,t)$ defined in (2.40) satisfies

$\displaystyle\left\\{\begin{array}[]{ll}V_{t}-U_{x}=0,\\\\[5.69054pt]
U_{t}+P_{x}=-R_{1},\\\\[5.69054pt]
\tfrac{R}{\gamma-1}\Theta_{t}+PU_{x}=\kappa\big{(}\frac{\Theta_{x}}{V}\big{)}_{x}-R_{2},\end{array}\right.$
(3.4)

where

$\displaystyle\begin{aligned} P:=&\frac{R\Theta}{V},\quad
P_{i}=\frac{R\Theta^{r}_{i}}{V^{r}_{i}}\ (i=1,3),\\\\[5.69054pt]
R_{1}:=&-(P-P_{1}-P_{3}-p^{c})_{x}+U^{c}_{t}:=R^{1}_{1}+U^{c}_{t},\\\\[5.69054pt]
R_{2}:=&\big{\\{}(p_{m}-P)U^{c}_{x}+(P_{1}-P)U^{r}_{1x}+(P_{3}-P)U^{r}_{3x}\big{\\}}\\\\[5.69054pt]
&+\kappa\bigg{\\{}\Big{(}\frac{\Theta_{x}}{V}\Big{)}_{x}-\Big{(}\frac{\Theta^{c}_{x}}{V^{c}}\Big{)}_{x}\bigg{\\}}:=R^{1}_{2}+R^{2}_{2}.\end{aligned}$
(3.5)

Let the perturbation function is

$\displaystyle(\phi,\psi,\xi):=(v,u,\theta)-(V,U,\Theta),$

then the reformed equation is

$\displaystyle\left\\{\begin{array}[]{ll}\phi_{t}-\psi_{x}=0,\\\\[5.69054pt]
\psi_{t}+\Big{(}\frac{R\xi}{v}-\frac{P\phi}{v}\Big{)}_{x}=R_{1},\\\\[5.69054pt]
\tfrac{R}{\gamma-1}\xi_{t}+p\psi_{x}+(p-P)U_{x}=\kappa\Big{(}\frac{\xi_{x}}{v}-\frac{\Theta_{x}\phi}{vV}\Big{)}_{x}+R_{2}\end{array}\right.$
(3.9)

with the initial data

$\displaystyle(\phi,\psi,\xi)(x,0)$ $\displaystyle=$
$\displaystyle(\phi_{0},\psi_{0},\xi_{0})(x)$ $\displaystyle=$
$\displaystyle(v_{0}(x)-V(x,0),u_{0}(x)-U(x,0),\theta_{0}(x)-\Theta(x,0))\,.$

And the solution space is defined as

$\displaystyle\begin{aligned}
X(0,t):=&\big{\\{}(\phi,\psi,\xi)|(\phi,\psi,\xi)\in
C([0,t),H^{2}({\mathbb{R}})),\\\\[5.69054pt] &(\phi,\psi)_{x}\in
L^{2}([0,t),H^{1}({\mathbb{R}})),\xi_{x}\in
L^{2}([0,t),H^{2}({\mathbb{R}}))\big{\\}}\,.\end{aligned}$ (3.11)

The local existence is known in [34]

######  Proposition 3.1.

(Local existence) Under the assumptions stated in Theorem 2.1, the Cauchy
problem (3.9), (3.1) admits a unique smooth solution $(\phi,\psi,\xi)(x,t)\in
X(0,t_{1})$ for some sufficient small $t_{1}>0$, and $(\phi,\psi,\xi)(x,t)$
satisfies

$\displaystyle\sup\limits_{0\leq t\leq
t_{1}}\|(\phi,\psi,\xi)(t)\|^{2}_{2}\leq
2\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}\,.$ (3.12)

Suppose that $(\phi,\psi,\xi)(x,t)$ has been extended to the time $t>t_{1}$,
we want to get the following a priori estimates to obtain a global solution.

######  Proposition 3.2.

(A prior estimates) Under the assumptions stated in Theorem 2.1, there exist
positive constants $\epsilon_{2}\leq 1$,
$\delta_{2}\leq\min\\{\delta_{1},\overline{\delta},1\\}$ and C, such that
$(\phi,\psi,\xi)\in X([0,t])$ for some $t>0$ satisfying

$\displaystyle\begin{aligned} &N(t)=\sup\limits_{0\leq\tau\leq
t}\|(\phi,\psi,\xi)(\tau)\|_{2}\leq\epsilon_{2}\,,\\\\[5.69054pt]
&\delta=|\theta_{-}-\theta_{+}|<\delta_{2}\,,\end{aligned}$ (3.13)

it follows the estimate

$\displaystyle\begin{aligned} &\sup\limits_{0\leq\tau\leq
t}\|(\phi,\psi,\xi)(\tau)\|^{2}_{2}+\int^{t}_{0}\big{(}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}+\|\xi_{x}(\tau)\|^{2}_{2}\big{)}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+\delta_{2}\,.\end{aligned}$
(3.14)

Once Proposition 3.2 is proved, we can extend the unique local solution
$(\phi,\psi,\xi)(x,t)$ obtained in Proposition 3.1 to $t=\infty$, moreover,
estimate (3.14) implies that

$\displaystyle\int^{\infty}_{0}\big{(}\|(\phi_{x},\psi_{x},\xi_{x})(t)\|^{2}+\big{|}\frac{d}{dt}\|(\phi_{x},\psi_{x},\xi_{x})(t)\|^{2}\big{|}\big{)}d\tau\lesssim+\infty\,,$
(3.15)

which together with Sobolev inequality easily leads to the asymptotic behavior
(2.45), this concludes the proof of Theorem 2.1. In the rest of this section,
our main task is to show this a prior estimates.

###  3.2. A Priori Estimates

At first, we show the basic estimates.

######  Lemma 3.1.

Under the assumptions in proposition 3.2, we have

$\displaystyle\begin{aligned}
&\|(\phi,\psi,\xi)(t)\|^{2}+\int^{t}_{0}\int_{\mathbb{R}}\big{(}(|U^{r}_{1x}|+|U^{r}_{3x}|)(\phi^{2}+\xi^{2})+\xi^{2}_{x}\big{)}dxd\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}+\delta+\delta^{\frac{1}{2}}\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}d\tau\\\\[5.69054pt]
&+\delta\int^{t}_{0}\frac{1}{1+\tau}\int_{\mathbb{R}}(\phi^{2}+\xi^{2})e^{\frac{-cx^{2}}{1+\tau}}dxd\tau.\end{aligned}$
(3.16)

###### Proof.

The first step, by using lemma 2.1, we need get some estimates of $R_{1}$ and
$R_{2}$. Since $R\Theta^{c}=p_{m}V^{c}$, direct calculation yields that

$\displaystyle R^{1}_{1}$ $\displaystyle=$ $\displaystyle
R\Big{(}\frac{\Theta^{r}_{1}}{V^{r}_{1}}+\frac{\Theta^{r}_{3}}{V^{r}_{3}}+\frac{\Theta^{cd}}{V^{cd}}-\frac{\Theta}{V}\Big{)}_{x}$
$\displaystyle=$ $\displaystyle
R\Theta^{r}_{1x}((V^{r}_{1})^{-1}-V^{-1})+R\Theta^{r}_{3x}((V^{r}_{3})^{-1}-V^{-1})$
$\displaystyle+R\Theta^{cd}_{x}((V^{cd})^{-1}-V^{-1})+RV^{r}_{1x}\Big{(}\frac{\Theta}{V^{2}}-\frac{\Theta^{r}_{1}}{V^{r2}_{1}}\Big{)}$
$\displaystyle+RV^{r}_{3x}\Big{(}\frac{\Theta}{V^{2}}-\frac{\Theta^{r}_{3}}{V^{r2}_{3}}\Big{)}+RV^{cd}_{x}\Big{(}\frac{\Theta}{V^{2}}-\frac{\Theta^{c}}{V^{c2}}\Big{)}.$

It is easy to compute

$\displaystyle|\Theta^{r}_{1x}((V^{r}_{1})^{-1}-V^{-1})|\lesssim|\Theta^{r}_{1x}|(|V^{r}_{3}-v^{m}_{+}|+|V^{c}-v^{m}_{-}|)\lesssim\delta
e^{-c(|x|+t)}$

and we can treat the other terms on the righthand side of (3.2) in the same
way to obtain

$\displaystyle|R^{1}_{1}|\lesssim\delta e^{-c(|x|+t)}$ (3.18)

and

$\displaystyle|R_{1}|\lesssim|R^{1}_{1}|+|U^{c}_{t}|\lesssim\delta
e^{-c(|x|+t)}+\tfrac{\delta}{(1+t)^{\frac{3}{2}}}e^{\frac{-cx^{2}}{1+t}}.$
(3.19)

Similarly we have

$\displaystyle|R^{1}_{2}|\lesssim\delta e^{-c(|x|+t)}.$

Since

$\displaystyle
R^{2}_{2}=\kappa\Big{(}\frac{\Theta^{r}_{1x}}{V}+\frac{\Theta^{r}_{3x}}{V}\Big{)}_{x}+\kappa\Big{(}\frac{\Theta^{cd}_{x}}{V}-\frac{\Theta^{cd}_{x}}{V^{cd}}\Big{)}_{x}:=R^{2}_{21}+R^{2}_{22}$

and

$\displaystyle R^{2}_{21}$ $\displaystyle\lesssim$
$\displaystyle\Big{(}|\Theta^{r}_{1xx}|+|\Theta^{r}_{3xx}|+|\Theta^{r}_{1x}||V^{r}_{1x}|+|\Theta^{r}_{3x}||V^{r}_{3x}|\Big{)}$
$\displaystyle+|\Theta^{r}_{1x}|(|V^{r}_{3x}|+|V^{cd}_{x}|)+|\Theta^{r}_{3x}|(|V^{r}_{1x}|+|V^{cd}_{x}|),$

it follows from (2.10) and Lemma 2.1 that

$\displaystyle\begin{aligned}
|R^{2}_{21}|\lesssim&\delta^{\frac{1}{8}}(1+t)^{-\frac{7}{8}},\\\\[5.69054pt]
|R^{2}_{22}|\lesssim&\Big{(}|\Theta^{c}_{xx}|+|\Theta^{c}_{x}||V^{c}_{x}|\Big{)}\Big{(}|V^{r}_{3}-v^{m}_{3}|+|V^{r}_{1}-v^{m}_{-}|\Big{)}\\\\[5.69054pt]
&+|\Theta^{c}_{x}|(|V^{r}_{1x}|+|V^{r}_{3x}|)\lesssim\delta
e^{-c(|x|+t)}.\end{aligned}$ (3.21)

Thus we get

$\displaystyle|R_{2}|\lesssim\delta^{\frac{1}{8}}(1+t)^{-\frac{7}{8}}\,.$
(3.22)

And by the direct calculation, we can also obtain

$\displaystyle\begin{aligned} &|(R_{1x},R_{1xx})(x,t)|\lesssim\delta
e^{-c(|x|+t)}+\frac{\delta}{(1+t)^{\frac{3}{2}}}e^{\frac{-cx^{2}}{1+t}},\\\\[5.69054pt]
&|(R_{2x},R_{2xx})(x,t)|\lesssim\delta^{\frac{1}{8}}(1+t)^{-\frac{7}{8}}.\end{aligned}$
(3.23)

Then multiplying $(\ref{3.3})_{1}$ by
$-R\Theta\big{(}\frac{1}{v}-\frac{1}{V}\big{)}$, $(\ref{3.3})_{2}$ by $\psi$
and $(\ref{3.3})_{3}$ by $\frac{\xi}{\theta}$, respectively, and adding the
results together, we get

$\displaystyle\begin{aligned}
&\bigg{\\{}R\Theta\Phi\Big{(}\frac{v}{V}\Big{)}+\tfrac{1}{2}\psi^{2}+\tfrac{R\Theta}{\gamma-1}\Phi\Big{(}\frac{\theta}{\Theta}\Big{)}\bigg{\\}}_{t}+\frac{\kappa}{v\theta}\xi^{2}_{x}+H_{1x}+Q_{1}+Q_{2}\\\\[5.69054pt]
=&R_{1}\psi+R_{2}\frac{\xi}{\theta},\end{aligned}$ (3.24)

where

$\displaystyle\begin{aligned}
H_{1}&:=(p-P)\psi-\frac{\kappa\xi}{\theta}\Big{(}\frac{\xi_{x}}{v}-\frac{\Theta_{x}\phi}{vV}\Big{)},\\\\[5.69054pt]
Q_{1}&:=-R\Theta_{t}\Phi\Big{(}\frac{v}{V}\Big{)}+\frac{PU_{x}}{vV}\phi^{2}+\frac{R\Theta_{t}}{\gamma-1}\Phi\Big{(}\frac{\Theta}{\theta}\Big{)}+\frac{\xi}{\theta}(p-P)U_{x},\\\\[5.69054pt]
Q_{2}&:=-\kappa\frac{\theta_{x}}{\theta^{2}v}\xi\xi_{x}-\kappa\frac{\xi_{x}\phi}{\theta
vV}\Theta_{x}+\kappa\frac{\theta_{x}\xi\phi}{\theta^{2}vV}\Theta_{x},\end{aligned}$
(3.25)

here

$\Phi(s)=s-1-\ln s.$

In the process of the calculation, we have used the following equalities

$\displaystyle\begin{aligned}
-R\Theta\Big{(}\frac{1}{v}-\frac{1}{V}\Big{)}\phi_{t}&=\bigg{\\{}R\Theta\Phi\Big{(}\frac{v}{V}\Big{)}\bigg{\\}}_{t}+\frac{PU_{x}}{vV}\phi^{2}-R\Theta_{t}\Phi\Big{(}\frac{v}{V}\Big{)},\\\\[5.69054pt]
\frac{R\xi_{t}}{\gamma-1}\frac{\xi}{\theta}&=\bigg{\\{}\frac{R\Theta}{\gamma-1}\Phi\Big{(}\frac{\theta}{\Theta}\Big{)}\bigg{\\}}_{t}+\frac{R\Theta_{t}}{\gamma-1}\Phi\Big{(}\frac{\Theta}{\theta}\Big{)}\end{aligned}$
(3.26)

and the equations in (3.9).

Since the composite waves affected each other, we will try to divide it the
viscous contact wave $Z^{c}$ and viscous rarefaction wave $Z^{r}_{i}(i=1,3)$,

$\displaystyle\begin{aligned}
-R\Theta_{t}=&(\gamma-1)(P^{r}_{1}U^{r}_{1x}+P^{r}_{3}U^{r}_{3x})+(\gamma-1)\bigg{(}p^{m}U^{c}_{x}-\kappa\Big{(}\frac{\Theta^{c}_{x}}{V^{c}}\Big{)}_{x}\bigg{)}\\\\[5.69054pt]
=&(\gamma-1)P(U^{r}_{1x}+U^{r}_{3x})+(\gamma-1)(P^{r}_{1}-P)U^{r}_{1x}\\\\[5.69054pt]
&+(\gamma-1)(P^{r}_{3}-P)U^{r}_{3x}+(\gamma-1)\bigg{(}p^{m}U^{c}_{x}-\kappa\Big{(}\frac{\Theta^{c}_{x}}{V^{c}}\Big{)}_{x}\bigg{)},\end{aligned}$
(3.27)

thus

$\displaystyle\begin{aligned}
Q_{1}=&-R\Theta_{t}\Phi\Big{(}\frac{v}{V}\Big{)}+\frac{R\Theta_{t}}{\gamma-1}\Phi\Big{(}\frac{\Theta}{\theta}\Big{)}+\frac{P(U_{1x}+U_{3x}+U^{c}_{x})}{vV}\phi^{2}\\\\[5.69054pt]
&+\frac{\xi}{\theta}(p-P)(U_{1x}+U_{3x}+U^{c}_{x})\\\\[5.69054pt]
=&:(|U^{r}_{1x}|+|U^{r}_{3x}|)Q_{11}+Q_{12},\end{aligned}$ (3.28)

where

$\displaystyle\begin{aligned}
Q_{11}:=&(\gamma-1)P\Phi\Big{(}\frac{v}{V}\Big{)}-P\Phi\Big{(}\frac{\Theta}{\theta}\Big{)}+\frac{P}{vV}\phi^{2}+\frac{\xi}{\theta}(p-P)\\\\[5.69054pt]
=&P\bigg{(}\frac{\theta V}{\Theta
v}-1+\gamma\Big{(}\frac{\theta}{\Theta}\Big{)}-\Big{(}\log\frac{\theta}{\Theta}+(\gamma-1)\log\frac{v}{V}\Big{)}\bigg{)}\\\\[5.69054pt]
=&P\bigg{(}\Phi\Big{(}\frac{\theta V}{\Theta
v}\Big{)}+\gamma\Phi\Big{(}\frac{v}{V}\Big{)}\bigg{)}\geq
C(\phi^{2}+\xi^{2})\end{aligned}$ (3.29)

and

$\displaystyle\begin{aligned}
Q_{12}:=&U^{c}_{x}\bigg{(}\frac{P\phi^{2}}{vV}+(\gamma-1)p_{m}\Phi\Big{(}\frac{v}{V}\Big{)}+p_{m}\Phi\Big{(}\frac{\Theta}{\theta}\Big{)}+\frac{\xi}{\theta}(p-P)\bigg{)}\\\\[5.69054pt]
&+(P^{r}_{1}-P)U^{r}_{1x}\Big{[}(\gamma-1)\Phi\Big{(}\frac{v}{V}\Big{)}-\Phi\Big{(}\frac{v}{V}\Big{)}\Big{]}\\\\[5.69054pt]
&-\kappa\Big{(}\frac{\Theta^{c}_{x}}{V^{c}}\Big{)}_{x}\Big{[}(\gamma-1)\Phi\Big{(}\frac{v}{V}\Big{)}-\Phi\Big{(}\frac{v}{V}\Big{)}\Big{]}\\\\[5.69054pt]
&+(P^{r}_{3}-P)U^{r}_{3x}[(\gamma-1)\Phi\Big{(}\frac{v}{V}\Big{)}-\Phi\Big{(}\frac{v}{V}\Big{)}\Big{]}\\\\[5.69054pt]
\lesssim&|(U^{c}_{x},\Theta^{c}_{xx},\Theta^{c}_{x}V^{c}_{x})|(\phi^{2}+\xi^{2})+\delta(U^{r}_{1x}+U^{r}_{3x})(\phi^{2}+\xi^{2})\,.\end{aligned}$
(3.30)

Meanwhile,

$\displaystyle|Q_{2}|\lesssim\big{(}N(t)+\delta+\tfrac{1}{8}\big{)}\xi^{2}_{x}+\Theta^{2}_{x}\big{(}\phi^{2}+\xi^{2}\big{)},$
(3.31)

where

$\displaystyle\begin{aligned}
\Theta^{2}_{x}\lesssim&({\Theta^{r}_{1x}}^{2}+{\Theta^{r}_{3x}}^{2})+({\Theta^{c}_{x}}^{2})\\\\[5.69054pt]
\lesssim&{\Theta^{r}_{1x}}^{2}+{\Theta^{r}_{3x}}^{2}+\frac{\delta}{1+t}e^{\frac{-cx^{2}}{1+t}}\,.\end{aligned}$
(3.32)

Therefore, after integrating (3.25) on $[0,t]\times\mathbb{R}$ and using the
above estimates, we obtain

$\displaystyle\begin{aligned}
&\|(\phi,\psi,\xi)(t)\|^{2}+\int^{t}_{0}\int_{\mathbb{R}}\big{(}(|U^{r}_{1x}|+|U^{r}_{3x}|)(\phi^{2}+\xi^{2})+\xi^{2}_{x}\big{)}dxd\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}+\int^{t}_{0}\frac{\delta}{1+\tau}\int_{\mathbb{R}}(\phi^{2}+\xi^{2})e^{\frac{-cx^{2}}{1+\tau}}dxd\tau\\\\[5.69054pt]
&+\int^{t}_{0}\int_{\mathbb{R}}({\Theta^{r}_{1x}}^{2}+{\Theta^{r}_{3x}}^{2})(\phi^{2}+\xi^{2})dxd\tau+\int^{t}_{0}\int_{\mathbb{R}}(|\psi||R_{1}|+|\xi||R_{2}|)dxd\tau.\end{aligned}$
(3.33)

Noticing that
$\|(\Theta^{r}_{1x},\Theta^{r}_{3x})\|_{\infty}\lesssim\delta^{\frac{1}{8}}(1+t)^{-\frac{7}{8}}$,
one can easily get

$\displaystyle\int^{t}_{0}\int_{\mathbb{R}}({\Theta^{r}_{1x}}^{2}+{\Theta^{r}_{3x}}^{2})(\phi^{2}+\xi^{2})dxd\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{4}}\int^{t}_{0}\|(\phi,\xi)\|^{2}_{\infty}(1+\tau)^{-\frac{7}{4}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{4}}\int^{t}_{0}\|(\phi,\xi)\|\|(\phi_{x},\xi_{x})\|(1+\tau)^{-\frac{7}{4}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{4}}\int^{t}_{0}\|(\phi_{x},\xi_{x})\|^{2}d\tau+\delta^{\frac{1}{4}}\int^{t}_{0}\|(\phi,\xi)\|^{2}(1+\tau)^{-\frac{7}{2}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{4}}\int^{t}_{0}\|(\phi_{x},\xi_{x})\|^{2}d\tau+\delta^{\frac{1}{4}}.$

For the last term in (3.33), we have the following estimate

$\displaystyle\int^{t}_{0}\int_{\mathbb{R}}(|\psi||R_{1}|+|\xi||R_{2}|)dxd\tau$
$\displaystyle\lesssim$ $\displaystyle\delta\int^{t}_{0}\|\psi\
|_{\infty}(\int_{\mathbb{R}}e^{-c(|x|+\tau)}dx+\int_{\mathbb{R}}\frac{\delta}{(1+\tau)^{\frac{3}{2}}}e^{\frac{-cx^{2}}{1+\tau}}dx)d\tau$
$\displaystyle+\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi\|_{\infty}(1+\tau)^{-\frac{7}{8}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta\int^{t}_{0}\|\psi\|^{\frac{1}{2}}\|\psi_{x}\|^{\frac{1}{2}}(1+\tau)^{-1}d\tau+\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi\|^{\frac{1}{2}}\|\xi_{x}\|^{\frac{1}{2}}(1+\tau)^{-\frac{7}{8}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta\int^{t}_{0}\|(\psi_{x},\xi_{x})\|^{2}d\tau+\delta^{\frac{1}{8}}\int^{t}_{0}\|(\psi,\xi)\|^{\frac{2}{3}}(1+\tau)^{-\frac{7}{6}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta\int^{t}_{0}\|(\psi_{x},\xi_{x})\|^{2}d\tau+\delta.$

Inserting (3.2) and (3.2) into (3.33), we can get (3.16) and this completes
the proof of Lemma 3.1. ∎

Then we claim that there exists some positive constant $C$ such that

######  Lemma 3.2.

([13]) Under the assumptions in proposition 3.2, we have

$\displaystyle\begin{aligned}
&\int^{t}_{0}\frac{1}{1+\tau}\int_{\mathbb{R}}(\phi^{2}+\psi^{2}+\xi^{2})e^{-\frac{cx^{2}}{1+\tau}}dxd\tau\\\\[5.69054pt]
\lesssim&1+\int^{t}_{0}\|(\phi_{x},\psi_{x},\xi_{x})\|^{2}d\tau+\int^{t}_{0}\int_{\mathbb{R}}(|U^{r}_{1x}|+|U^{r}_{3x}|)(\phi^{2}+\xi^{2})dxd\tau.\end{aligned}$
(3.36)

For $\alpha>0$, we define the following heat kernel $\omega(x,t)$ which will
play an essential role in the later estimates. We define

$\displaystyle\omega(x,t)=(1+t)^{-\frac{1}{2}}e^{-\frac{\alpha
x^{2}}{1+t}}\,,\quad g(x,t)=\int_{-\infty}^{x}\omega(y,t)dy\,.$

It is easy to check that

$\displaystyle 4\alpha
g_{t}=g_{xx}\,,\quad\|g(\cdot,t)\|_{L^{\infty}}=\sqrt{\pi}\alpha^{-\frac{1}{2}}\,.$

The proof (3.36) is divided into two parts:

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}\omega^{2}((R\xi-P\phi)^{2}+\psi^{2})dxd\tau\\\\[5.69054pt]
\lesssim&1+\int^{t}_{0}\|(\phi_{x},\psi_{x},\xi_{x})(\tau)\|^{2}d\tau+\delta\int^{t}_{0}\int_{\mathbb{R}}\omega^{2}(\phi^{2}+\psi^{2}+\xi^{2})dxd\tau\\\\[5.69054pt]
&+\delta\int^{t}_{0}\int_{\mathbb{R}}(|U^{r}_{1x}|+|U^{r}_{3x}|)(\phi^{2}+\psi^{2}+\xi^{2})dxd\tau\end{aligned}$
(3.37)

and for any $\eta>0$,

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}\omega^{2}(R\xi+(\gamma-1)P\phi)^{2}dxd\tau\\\\[5.69054pt]
\lesssim&1+\int^{t}_{0}\|(\phi_{x},\psi_{x},\xi_{x})(\tau)\|^{2}d\tau+(\delta+\eta)\int^{t}_{0}\int_{\mathbb{R}}\omega^{2}(\phi^{2}+\psi^{2}+\xi^{2})dxd\tau,\\\\[5.69054pt]
&+\delta\int^{t}_{0}\int_{\mathbb{R}}(|U^{r}_{1x}|+|U^{r}_{3x}|)(\phi^{2}+\psi^{2}+\xi^{2})dxd\tau\,.\end{aligned}$
(3.38)

Here we used the same method as in [13], the only difference lied in the loss
of diffusion term of fluid velocity, but that’s not the main term. We omit the
detail for simplify.

We now turn to obtain the higher order estimates, further calculation yields
the following Lemma

######  Lemma 3.3.

Under the assumptions in proposition 3.2, we have

$\displaystyle\begin{aligned}
&\|(\phi_{x},\psi_{x},\xi_{x})(t)\|^{2}+\int^{t}_{0}\|\xi_{xx}(\tau)\|^{2}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0x},\psi_{0x},\xi_{0x})\|^{2}+\delta+(\delta+N(t))\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\,.\end{aligned}$
(3.39)

###### Proof.

Multiplying $(\ref{3.3})_{1x}$ by $\frac{P}{v}\phi_{x}$, $(\ref{3.3})_{2x}$ by
$\psi_{x}$ and $(\ref{3.3})_{3x}$ by $\frac{\xi_{x}}{\theta}$, respectively,
and adding the resulting equations together, then we have

$\displaystyle\begin{aligned}
&\bigg{\\{}\frac{P}{2v}\phi_{x}^{2}+\frac{\psi^{2}_{x}}{2}+\frac{R\xi^{2}_{x}}{2(\gamma-1)\theta}\bigg{\\}}_{t}+\frac{\kappa}{v\theta}\xi^{2}_{xx}+H_{2x}+J_{2}\\\\[5.69054pt]
=&R_{1x}\psi_{x}+R_{2x}\frac{\xi_{x}}{\theta}\,,\end{aligned}$ (3.40)

where

$\displaystyle\begin{aligned}
H_{2}=&(p-P)_{x}\psi_{x}+\frac{\xi_{x}}{\theta}\bigg{(}(p-P)U_{x}+\Big{(}\frac{\kappa\xi_{x}}{v}-\frac{\kappa\Theta_{x}\phi}{vV}\Big{)}_{x}\bigg{)}\,,\\\\[5.69054pt]
J_{2}=&\Big{(}\frac{\xi_{x}}{\theta}\Big{)}_{x}\bigg{(}\Big{(}\frac{\kappa\xi_{x}}{v}-\frac{\kappa\Theta_{x}\phi}{vV}\Big{)}_{x}-(p-P)U_{x}\bigg{)}-\frac{\kappa}{v\theta}\xi^{2}_{xx}-\Big{(}\frac{P}{2v}\Big{)}_{t}\phi^{2}_{x}\\\\[5.69054pt]
&-\frac{R}{2(\gamma-1)}\Big{(}\frac{1}{\theta}\Big{)}_{t}\xi^{2}_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi\psi_{xx}+\Big{(}\frac{P}{v}\Big{)}_{x}\phi\psi_{xx}+p_{x}\psi_{x}\tfrac{\xi_{x}}{\theta}\\\\[5.69054pt]
=&O(1)\big{(}N(t)+\delta+\eta\big{)}|(\phi_{x},\xi_{x},\psi_{xx},\xi_{xx})|^{2}+|(V_{x},U_{x},\Theta_{x},\Theta_{xx})|^{2}(\phi^{2}+\xi^{2})\,.\end{aligned}$
(3.41)

After integrating (3.40) on $[0,t]\times\mathbb{R}$, we get

$\displaystyle\begin{aligned}
&\|(\phi_{x},\psi_{x},\xi_{x})(t)\|^{2}+\int^{t}_{0}\|\xi_{xx}(\tau)\|^{2}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0x},\psi_{0x},\xi_{0x})\|^{2}+(\delta+N(t)+\eta)\int^{t}_{0}\|(\phi_{x},\psi_{x},\xi_{x})(\tau)\|^{2}_{1}dxd\tau\\\\[5.69054pt]
&+\int^{t}_{0}\int_{\mathbb{R}}(|\Theta_{xx}|+|\Theta_{x}|)^{2}(\phi^{2}+\xi^{2})dxd\tau+\int^{t}_{0}\int_{\mathbb{R}}(|R_{1x}\psi_{x}|+|R_{2x}\xi_{x}|)dxd\tau\,,\end{aligned}$
(3.42)

where $\eta>0$ is a constant suitably small, the last two terms on the right
hand side of (3.42) can be treated similarly as (3.2) and (3.2), respectively.
Then, with the help of the results of Lemma 3.1 and Lemma 3.2 we can deduce
the estimate (3.39), and this completes the proof of Lemma 3.3. ∎

######  Lemma 3.4.

Under the assumption in proposition 3.2, we have

$\displaystyle\begin{aligned}
&\|(\phi_{xx},\psi_{xx},\xi_{xx})(t)\|^{2}+\int^{t}_{0}\|\xi_{xxx}(\tau)\|^{2}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0xx},\psi_{0xx},\xi_{0xx})\|^{2}+\delta+(\delta+N(t))\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\,.\end{aligned}$
(3.43)

###### Proof.

Multiplying $(\ref{3.3})_{1xx}$ by $\frac{P}{v}\phi_{xx}$, $(\ref{3.3})_{2xx}$
by $\psi_{xx}$ and $(\ref{3.3})_{3xx}$ by $\frac{\xi_{xx}}{\theta}$,
respectively, and adding the results together, it is easily to get

$\displaystyle\begin{aligned}
&\bigg{\\{}\frac{P}{2v}\phi_{xx}^{2}+\frac{\psi^{2}_{xx}}{2}+\frac{R\xi^{2}_{xx}}{2(\gamma-1)\theta}\bigg{\\}}_{t}+\frac{\kappa}{v\theta}\xi^{2}_{xxx}+H_{3x}+J_{3}\\\\[5.69054pt]
=&R_{1xx}\psi_{xx}+R_{1xx}\frac{\xi_{xx}}{v}\,,\end{aligned}$ (3.44)

where

$\displaystyle\begin{aligned}
H_{3}=&(p-P)_{xx}\psi_{xx}+\frac{\xi_{xx}}{\theta}\bigg{(}(p-P)U_{x}-\Big{(}\frac{\kappa\xi_{x}}{v}-\frac{\kappa\Theta_{x}\phi}{vV}\Big{)}_{x}\bigg{)}_{x},\\\\[5.69054pt]
J_{3}=&\Big{(}\frac{\xi_{xx}}{\theta}\Big{)}_{x}\bigg{(}\Big{(}\frac{\kappa\xi_{x}}{v}-\frac{\kappa\Theta_{x}\phi}{vV}\Big{)}_{x}-(p-P)U_{x}\bigg{)}_{x}-\frac{\kappa}{v\theta}\xi^{2}_{xxx}\\\\[5.69054pt]
&-\Big{(}\frac{P}{2v}\Big{)}_{t}\phi_{xx}^{2}-\Big{(}\frac{R}{2(\gamma-1)\theta}\Big{)}_{t}\xi^{2}_{xx}+(2p_{x}\psi_{xx}+p_{xx}\psi_{x})\frac{\xi_{xx}}{\theta}\\\\[5.69054pt]
&+3\bigg{\\{}\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{\\}}_{x}\psi_{xx}-3\bigg{\\{}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}\bigg{\\}}_{x}\psi_{xx}+2J^{1}_{3}+J^{2}_{3}\,.\end{aligned}$
(3.45)

Here $J^{1}_{3}$, $J^{2}_{3}$ are the following equalities

$\displaystyle\begin{aligned}
J^{1}_{3}:&=\psi_{xxx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{)}\,,\\\\[5.69054pt]
J^{2}_{3}:&=\psi_{xxx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{xx}\phi-\Big{(}\frac{R}{v}\Big{)}_{xx}\xi\bigg{)}\,.\end{aligned}$

Meanwhile, we can get

$\displaystyle\begin{aligned}
J^{1}_{3}=&\bigg{\\{}\psi_{xx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{)}\bigg{\\}}_{x}-\psi_{xx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{)}_{x}\\\\[5.69054pt]
=&\bigg{\\{}\psi_{xx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{)}\bigg{\\}}_{x}+O(1)(N(t)+\delta)|(\phi_{x},\xi_{x},\phi_{xx},\psi_{xx},\xi_{xx})|^{2}\end{aligned}$
(3.46)

and

$\displaystyle\begin{aligned}
&\psi_{xxx}\Big{(}\frac{P}{v}\Big{)}_{xx}\phi\\\\[5.69054pt]
=&\psi_{xxx}\phi\bigg{(}\frac{P_{xx}}{v}+2P_{x}\Big{(}\frac{1}{v}\Big{)}_{x}+P\Big{(}\frac{2v^{2}_{x}}{v^{3}}-\frac{V_{xx}}{v^{2}}\Big{)}\bigg{)}-\frac{P\phi}{v^{2}}\phi_{xx}\psi_{xxx}\\\\[5.69054pt]
=&-\frac{P\phi}{v^{2}}\phi_{xx}\phi_{txx}+\bigg{\\{}\psi_{xx}\phi\bigg{(}\frac{P_{xx}}{v}+2P_{x}\Big{(}\frac{1}{v}\Big{)}_{x}+P\Big{(}\frac{2v^{2}_{x}}{v^{3}}-\frac{V_{xx}}{v^{2}}\Big{)}\bigg{)}\bigg{\\}}_{x}\\\\[5.69054pt]
&-\psi_{xx}\bigg{\\{}\phi\bigg{(}\frac{P_{xx}}{v}+2P_{x}\Big{(}\frac{1}{v}\Big{)}_{x}+P\Big{(}\frac{2v^{2}_{x}}{v^{3}}-\frac{V_{xx}}{v^{2}}\Big{)}\bigg{)}\bigg{\\}}_{x}\\\\[5.69054pt]
=&-\bigg{\\{}\frac{P\phi}{v^{2}}\frac{\phi^{2}_{xx}}{2}\bigg{\\}}_{t}+\bigg{\\{}\psi_{xx}\phi\bigg{(}\frac{P_{xx}}{v}+2P_{x}\Big{(}\frac{1}{v}\Big{)}_{x}+P\Big{(}\frac{2v^{2}_{x}}{v^{3}}-\frac{V_{xx}}{v^{2}}\Big{)}\bigg{)}\bigg{\\}}_{x}\\\\[5.69054pt]
&+O(1)(N(t)+\delta)|(\phi_{x},\xi_{x},\phi_{xx},\psi_{xx},\xi_{xx})|^{2}+|(\Theta_{x},\Theta_{xx})|^{2}|(\phi,\xi)|^{2}\,.\end{aligned}$
(3.47)

Similar to the estimate of $\psi_{xxx}\Big{(}\frac{P}{v}\Big{)}_{xx}\phi$, we
have

$\displaystyle\begin{aligned}
\psi_{xxx}\Big{(}\frac{R}{v}\Big{)}_{xx}\xi=&\bigg{\\{}\frac{R\xi}{v^{2}}\frac{\phi^{2}_{xx}}{2}\bigg{\\}}_{t}+\bigg{\\{}R\psi_{xx}\xi\Big{(}\frac{V_{xx}}{v^{2}}-\frac{2v^{2}_{x}}{v^{3}}\Big{)}\bigg{\\}}_{x}+|(\Theta_{x},\Theta_{xx})|^{2}|(\phi,\xi)|^{2}\\\\[5.69054pt]
&+O(1)(N(t)+\delta)|(\phi_{x},\xi_{x},\phi_{xx},\psi_{xx},\xi_{xx})|^{2}\,.\end{aligned}$

Therefore, it holds

$\displaystyle\begin{aligned}
J_{3}=&\Big{\\{}\frac{R\xi}{v^{2}}\frac{\phi^{2}_{xx}}{2}-\frac{P\phi}{v^{2}}\frac{\phi^{2}_{xx}}{2}\Big{\\}}_{t}+\bigg{\\{}\psi_{xx}\phi\bigg{(}\frac{P_{xx}}{v}+2P_{x}\Big{(}\frac{1}{v}\Big{)}_{x}+P\Big{(}\frac{2v^{2}_{x}}{v^{3}}-\frac{V_{xx}}{v^{2}}\Big{)}\bigg{)}\bigg{\\}}_{x}\\\\[5.69054pt]
&+\bigg{\\{}\psi_{xx}R\xi\Big{(}\frac{V_{xx}}{v^{2}}-\frac{2v^{2}_{x}}{v^{3}}\Big{)}\bigg{\\}}_{x}+\bigg{\\{}\psi_{xx}\bigg{(}\Big{(}\frac{P}{v}\Big{)}_{x}\phi_{x}-\Big{(}\frac{R}{v}\Big{)}_{x}\xi_{x}\bigg{)}\bigg{\\}}_{x}\\\\[5.69054pt]
&+O(1)(N(t)+\delta)|(\phi_{x},\xi_{x},\phi_{xx},\psi_{xx},\xi_{xx})|^{2}+|(\Theta_{x},\Theta_{xx})|^{2}|(\phi,\xi)|^{2}\,.\end{aligned}$
(3.48)

After integrating (3.44) on $[0,t]\times\mathbb{R}$, we get

$\displaystyle\begin{aligned}
&\|(\phi_{xx},\psi_{xx},\xi_{xx})(t)\|^{2}+\int^{t}_{0}\|\xi_{xxx}(\tau)\|^{2}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+(\delta+N(t))\int^{t}_{0}\|(\phi_{x},\psi_{x},\xi_{x})(\tau)\|^{2}_{1}d\tau\\\\[5.69054pt]
&+\int^{t}_{0}\int_{\mathbb{R}}(|\Theta_{xx}|+|\Theta_{x}|)^{2}|(\phi,\xi)|^{2}dxd\tau\\\\[5.69054pt]
&+\int^{t}_{0}\int_{\mathbb{R}}|R_{1xx}\psi_{xx}+R_{2xx}\xi_{xx}|dxd\tau\,.\end{aligned}$
(3.49)

For the estimate of the last term in (3.49), we have

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}|R_{1xx}\psi_{xx}|dxd\tau\\\\[5.69054pt]
\lesssim&\delta\int^{t}_{0}\big{\|}\psi_{xx}(\tau)\big{\|}\bigg{\\{}\Big{(}\int_{\mathbb{R}}e^{-2c|x|}e^{-2c\tau}dx\Big{)}^{\frac{1}{2}}+\Big{(}\frac{1}{(1+\tau)^{3}}e^{\frac{-2cx^{2}}{1+\tau}}dx\Big{)}^{\frac{1}{2}}\bigg{\\}}d\tau\\\\[5.69054pt]
\lesssim&\delta\int^{t}_{0}(1+t)^{-\frac{5}{4}}\|\psi_{xx}(\tau)\|d\tau\lesssim\delta\int^{t}_{0}\|\psi_{xx}(\tau)\|^{2}d\tau+\delta\end{aligned}$
(3.50)

and

$\displaystyle\int^{t}_{0}\int_{\bf R}|R_{2xx}\xi_{xx}|dxd\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi_{xx}\|_{\infty}(1+\tau)^{-\frac{7}{8}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi_{xx}\|^{\frac{1}{2}}\|\xi_{xxx}\|^{\frac{1}{2}}(1+\tau)^{-\frac{7}{8}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi_{xx}\|\|\xi_{xxx}\|d\tau+\delta^{\frac{1}{8}}\int^{t}_{0}(1+\tau)^{-\frac{7}{3}}d\tau$
$\displaystyle\lesssim$
$\displaystyle\delta^{\frac{1}{8}}\int^{t}_{0}\|\xi_{xx}\|^{2}_{1}d\tau+\delta^{\frac{1}{8}}\,,$

then using the results of Lemma 3.1-Lemma 3.3, we can get (3.43) and this
completes the proof of Lemma 3.4. ∎

Combining the results of Lemma 3.1-Lemma 3.4, we kown that

$\displaystyle\begin{aligned}
&\|(\phi,\psi,\xi)(t)\|^{2}_{2}+\int^{t}_{0}\|\xi_{x}(\tau)\|^{2}_{2}d\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+\delta+(\delta+N(t)+\eta)\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\,.\end{aligned}$
(3.52)

Based on these analysis, we now deal with the last term on the right hand in
(3.52).

######  Lemma 3.5.

Under the assumption in proposition 3.2, we have

$\displaystyle\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\lesssim\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+\delta\,.$
(3.53)

###### Proof.

Multiplying $(\ref{3.3})_{2}$ by $-\frac{P}{2}\phi_{x}$, and $(\ref{3.3})_{3}$
by $\psi_{x}$, respectively, and adding all the resultant equations, we have

$\displaystyle\begin{aligned}
&\Big{\\{}\tfrac{R}{\gamma-1}\xi\psi_{x}-\frac{P}{2}\phi_{x}\psi\Big{\\}}_{t}+\Big{\\{}\frac{P}{2}\phi_{t}\psi-\tfrac{R}{\gamma-1}\xi\psi_{t}\Big{\\}}_{x}+\frac{P^{2}}{2v}\phi^{2}_{x}+\frac{P}{2}\psi^{2}_{x}\\\\[5.69054pt]
=&\frac{P_{x}}{2}\psi\psi_{x}-\frac{P_{t}}{2}\phi_{x}\psi+\frac{P}{2}\phi_{x}\bigg{(}\Big{(}\frac{R\xi}{v}\Big{)}_{x}-\Big{(}\frac{P}{v}\Big{)}_{x}\phi+R_{1}\bigg{)}-\tfrac{R}{\gamma-1}\xi_{x}\psi_{t}\\\\[5.69054pt]
&+\kappa\Big{(}\frac{\xi_{x}}{v}-\frac{\Theta_{x}\phi}{vV}\Big{)}_{x}\psi_{x}-(p-P)(U_{x}+\psi_{x})\psi_{x}+R_{2}\psi_{x}\,.\end{aligned}$
(3.54)

Integrating (3.54) on $[0,t]\times\mathbb{R}$ and using the inequality (3.52),
it holds that

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}(\phi^{2}_{x}+\psi^{2}_{x})dxd\tau\\\\[5.69054pt]
\lesssim&\|(\phi,\psi,\xi)\|^{2}_{1}+\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{1}+\int^{t}_{0}\|\xi_{x}(\tau)\|^{2}_{1}d\tau\\\\[5.69054pt]
&+\big{(}\tfrac{1}{4}+\delta+N(t)\big{)}\int^{t}_{0}\int_{\mathbb{R}}(\phi^{2}_{x}+\psi^{2}_{x})dxd\tau+\int^{t}_{0}\int_{\mathbb{R}}(|\Theta_{xx}|+|\Theta_{x}|)^{2}(\phi^{2}+\psi^{2})dxd\tau\\\\[5.69054pt]
&+\int^{t}_{0}\int_{\mathbb{R}}|R_{1}\phi_{x}+R_{2}\psi_{x}|dxd\tau\,.\end{aligned}$
(3.55)

Similar as the estimates of (3.2) and (3.2), we can easily get

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}{(\phi^{2}_{x}+\psi^{2}_{x})}dxd\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+\delta+(\delta+N(t))\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\,.\end{aligned}$
(3.56)

Similarly, multiplying $(\ref{3.3})_{2x}$ by $-\frac{P}{2}\phi_{xx}$,
$(\ref{3.3})_{3x}$ by $\psi_{xx}$, respectively, and integrating the result on
$[0,t]\times\mathbb{R}$, then by using(3.52), we can also get

$\displaystyle\begin{aligned}
&\int^{t}_{0}\int_{\mathbb{R}}(\phi^{2}_{xx}+\psi^{2}_{xx})dxd\tau\\\\[5.69054pt]
\lesssim&\|(\phi_{0},\psi_{0},\xi_{0})\|^{2}_{2}+\delta+(\delta+N(t))\int^{t}_{0}\|(\phi_{x},\psi_{x})(\tau)\|^{2}_{1}d\tau\,.\end{aligned}$
(3.57)

Adding the results of (3.56) with (3.57), we can get (3.53) and this completes
the proof of Lemma 3.5. ∎

Inserting (3.53) into (3.52), then it yields (3.14), that is the result of
Proposition 3.2.

Acknowledgments: The authors are grateful to Professor A. Matsumura his
support and advice. This work work was supported by the Fundamental Research
Funds for the Central Universities and three grants from the National Natural
Science Foundation of China under contracts 11301405, 11671309 and 11731008,
respectively.

## References

  * [1] Z. Z. Chen, L. He and H. J. Zhao, Nonlinear stability of traveling wave solutions for the compressible fluid models of Korteweg type, J. Math. Anal. Appl. 422 (2015) 1213–1234.
  * [2] Z. Z. Chen and Q. H. Xiao, Nonlinear stability of viscous contact wave for the one-dimensional compressible fluid models of Korteweg type, Math. Methods Appl. Sci. 36 (2013) 2265–2279.
  * [3] Z. Z. Chen, L. J. Xiong and Y. J. Meng, Convergence to the superposition of rarefaction waves and contact discontinuity for the 1-D compressible Navier-Stokes-Korteweg system, J. Math. Anal. Appl. 412 (2014) 646–663.
  * [4] Z. Z. Chen and H. J. Zhao, Existence and nonlinear stability of stationary solutions to the full compressible Navier-Stokes-Korteweg system, J. Math. Pures Appl. 101 (2014) 330–371.
  * [5] R. J. Duan, H. X. Liu and H. J. Zhao, Nonlinear stability of rarefaction waves for the compressible Navier-Stokes equations with large initial perturbation, Trans. Amer. Math. Soc. 361(1) (2009) 453–493.
  * [6] R. J. Duan and H. F. Ma, Global existence and convergence rates for the 3-D compressible Navier-Stokes equations without heat conductivity, Indiana Univ. Math. J. 57 (2008) 2299- C2319.
  * [7] L. L. Fan, H. X. Liu, T. Wang and H. J. Zhao, Inflow problem for the one-dimensional compressible Navier-Stokes equations under large initial perturbation, J. Differential Equations 257 (2014) 3521–3553.
  * [8] L. L. Fan and A. Matsumura, Asymptotic stability of a composite wave of two viscous shock waves for a one-dimensional system of non-viscous and heat-conductive ideal gas, J. Differential Equations 258 (2015) 1129–1157.
  * [9] J. Goodman, Nonlinear asymptotic stability of viscous shock profiles for conservation laws, Arch.Rational Mech. Anal. 95 (1986) 325–344.
  * [10] L. He, S. J. Tang and T. Wang, Stability of viscous shock waves for the one-dimensional compressible Navier-Stokes equations with density-dependent viscosity, Acta Math. Sci. Ser. B Engl. Ed. 36 (2016) 34–48.
  * [11] B. K. Huang and Y. K. Liao, Global stability of viscous contact wave with rarefaction waves for compressible Navier-Stokes equations with temperature-dependent viscosity, Math. Models Methods Appl. Sci. 27 (2017) 2321–2379.
  * [12] B. K. Huang, L. S. Wang and Q. H. Xiao, Global nonlinear stability of rarefaction waves for compressible Navier-Stokes equations with temperature and density dependent transport coefficients, Kinet. Relat. Models 9 (2016) 469–514.
  * [13] F. M. Huang, J. Li and A. Matsumura, Asymptotic stability of combination of viscous contact wave with rarefaction waves for one-dimensional compressible Navier-Stokes system, Arch. Rational Mech. Anal. 197 (2010) 89–116.
  * [14] F. M. Huang and A. Matsumura, Stability of a composite wave of two viscous shock waves for the full compresible Navier-Stokes equation, Comm. Math. Phys. 289 (2009) 841–861.
  * [15] F. M. Huang, A. Matsumura and X. D. Shi, On the stability of contact discontinuty for compressible Navier-Stokes equations for free boundary, Osaka J. Math. 41 (2004) 193–210.
  * [16] F. M. Huang and T. Wang, Stability of superposition of viscous contact wave and rarefaction waves for compressible Navier-Stokes system, Indiana Univ. Math. J. 65(6) (2016) 1833–1875.
  * [17] F. M. Huang, A. Matsumura and Z. P. Xin, Stability of contact discontinuties for the 1-D compressible Navier-Stokes equations, Arch. Ration. Mech. Anal. 179 (2005) 55–77.
  * [18] F. M. Huang, Z. P. Xin and T. Yang, Contact discontinuity with general perturbation for gas motions, Adv. Math. 219 (2008) 1246–1297.
  * [19] F. M. Huang and H. J. Zhao, On the global stability of contact discontinuity for compressible Navier-Stokes equations, Rend. Semin. Mat. Univ. Padova 109(2003) 283–305.
  * [20] S. Kawashima and A. Matsumura, Asymptotic stability of travelling wave solutions of systems for one-dimensional gas motion, Comm. Math. Phys. 101 (1985) 97–127.
  * [21] H. X. Liu, T. Yang, H. J. Zhao and Q. Y. Zou, One-dimensional compressible Navier-Stokes equations with temperature dependent transport coefficients and large data, SIAM J. Math. Anal. 46 (2014) 2185–2228.
  * [22] T. P. Liu, Nonlinear stability of shock waves for viscous conservation laws, Mem. Am. Math. Soc. 56(1985) 1–108.
  * [23] T. P. Liu, Shock wave for compresible Navier-Stokes equations are stable, Comm. Pure Appl. Math. 39(1986) 565–594.
  * [24] T. P. Liu and Z. P. Xin, Nonlinear stability of rarefaction waves for compressible Navier-Stokes equations, Comm. Math. Phys. 118(1988) 451–465.
  * [25] T. P. Liu and Y. N. Zeng, Large time behavior of solutions for general quasilinear hyperbolic-parabolic systems of conservation laws, Mem. Amer. Math. Soc. 125(599) (1997) pp. viii+120
  * [26] T. P. Liu and Y. N. Zeng, Compressible Navier-Stokes equations with zero heat conductivity, J. Differential Equations 153(2) (1999) 225–291.
  * [27] T. P. Liu and Y. N. Zeng, Time-asymptotic behavior of wave propagation around a viscous shock profile, Comm. Math. Phys. 290(1) (2009) 23–82.
  * [28] T. P. Liu and Y. N. Zeng, Shock waves in conservation laws with physical viscosity, Mem. Amer. Math. Soc. 234(1105) (2015) pp. vi+168
  * [29] S. X. Ma and J. Wang, Decay rates to viscous contact waves for the compressible Navier-Stokes equations, J. Math. Phys. 57 (2016) 1–14.
  * [30] A. Matsumura and K. Nishihara, Asymptotics toward the rarefaction waves of a one-dimensional model system for compressible viscous gas, Japan J. Appl. Math. 3 (1986) 1–13.
  * [31] A. Matsumura and K. Nishihara, On the stability of travelling wave solutions of a one-dimensional model system for compressible viscous gas, Japan J. Appl. Math. 2(1) (1985) 17–25.
  * [32] A. Matsumura and K. Nishihara, Global stability of the rarefaction wave of a one-dimensional model system for compressible viscous gas, Comm. Math. Phys. 144 (1992) 325–335.
  * [33] A. Matsumura and K. Nishihara, Global asymptotics toward the rarefaction waves for solutions of viscous p-system with boundary effect, Quart. Appl. Math. 58 (2000) 69–83.
  * [34] K. Nishihara, T. Yang and H. J. Zhao, Nonlinear stability of strong rarefaction waves for compressible Navier-Stokes equations, SIAM J. Math. Anal. 35 (2004) 1561–1597.
  * [35] L. Wan and T. Wang, Symmetric flows for compressible heat-conducting fluids with temperature dependent viscosity coefficients, J. Differential Equations 262 (2017) 5939–5977.
  * [36] L. Wan, T. Wang and H. J. Zhao, Asymptotic stability of wave patterns to compressible viscous and heat-conducting gases in the half-space, J. Differential Equations 261 (2016) 5949–5991.
  * [37] T. Wang and H. J. Zhao, One-dimensional compressible heat-conducting gas with temperature-dependent viscosity, Math. Models Methods Appl. Sci. 26 (2016) 2237–2275.
  * [38] T. Wang, H. J. Zhao and Q. Y. Zou, One-dimensional compressible Navier-Stokes equations with large density oscillation, Kinet. Relat. Models 6 (2013) 649–670.
  * [39] J. Smoller, Shock Wave and Reaction-Diffusion Equations, 2nd Edn. ( Springer-Verlag, New York, 1994).
  * [40] A. Szepessy and Z. P. Xin, Nonlinear stability of viscous shock waves, Arch. Ration. Mech. Anal. 122(1993) 53–103.
  * [41] Z. P. Xin, On nonlinear stability of contact discontinuities. In: Hyperbolic problems: theory, numerics, applications (World Sci. Publishing, River Edge, NJ, 1996), pp. 249–257.


Simultaneously sorting overlapping quantum states of light

Suraj Goel
These authors contributed equally

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

Max Tyler
These authors contributed equally

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

Feng Zhu

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

Saroch Leedumrongwatthanakun

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

Mehul Malik

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

Jonathan Leach

School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, UK

The efficient manipulation, sorting, and measurement of optical modes and single-photon states is fundamental to classical and quantum science.  Here, we realise simultaneous and efficient sorting of non-orthogonal, overlapping states of light, encoded in the transverse spatial degree of freedom.  We use a specifically designed multi-plane light converter (MPLC) to sort states encoded in dimensions ranging from $d = 3$ to $d = 7$. Through the use of an auxiliary output mode, the MPLC simultaneously performs the unitary operation required for unambiguous discrimination and the basis change for the outcomes to be spatially separated. Our results lay the groundwork for optimal image identification and classification via optical networks, with potential applications ranging from self-driving cars to quantum communication systems.

Introduction:– The task of discriminating between a set of quantum states is a fundamental requirement in quantum information science, and in particular, quantum communication <cit.>. However, in general, two different quantum states can have a finite, non-zero overlap with respect to each other, making them non-orthogonal and therefore difficult to separate. It is theoretically impossible to perform a measurement that allows us to perfectly distinguish between such states 100% of the time. An important question now follows: given a set of quantum states with a non-zero overlap, what is the best measurement strategy to distinguish between them?

The answer to this question lies in quantum measurement theory, where strategies for the optimal measurement of non-orthogonal quantum states are known <cit.>.  In the extreme, we are left with a choice between a measurement strategy that is either efficient OR accurate (only orthogonal states can be sorted efficiently AND accurately), see Fig. <ref>. The efficient option is minimum error state discrimination (MESD) <cit.>. Here, one seeks to perform a set of measurements that categorise every input state.  The drawback to MESD is that errors are inevitable, and we have to accept that we will be incorrect with some probability relating to the overlap of the input states.  The accurate option is unambiguous state discrimination (USD) <cit.>. Here, one seeks to perform measurements that never incorrectly identify the input state.  The downside here is that state identification occurs with a reduced probability, i.e. a measurement does not always provide a result, but when it does, it is always correct.    In addition to the extremes, there are intermediate strategies that have been developed that interpolate between MESD and USD.  Several theoretical studies include discrimination with a known error margin <cit.>, discrimination with a fixed rate of inconclusive outcomes <cit.>, and partial state separation with a MESD procedure <cit.>.

Measuring non-orthogonal states using minimum error or unambiguous state discrimination. The MESD protocol is efficient in that every input state is always categorised, even if this leads to errors.  The USD protocol is accurate in that every input state is correctly identified, even if this does not happen 100% of the time.

The problem that we address in this work is practical high-dimensional unambiguous state discrimination, i.e. a positive operator-valued measure (POVM) for non-orthogonal, high-dimensional states with simultaneous outcomes.  High-dimensional quantum states (or qudits) allow for quantum information to be encoded in a $d$-dimensional space, enabling quantum communication protocols with increased information capacity and robustness to noise <cit.>.  While the theoretical foundation for the necessary measurement strategies has already been developed, their experimental realisation has proved to be a significant challenge.

(a) Schematic of the mode sorter for non-orthogonal states with colour used to represent different input modes.  The $d$ input modes pass through a multi-plane light converter that sorts them into $d+1$ orthogonal outputs of spatially separated Gaussian spots.  The one additional output mode $\ket{?}$ corresponds to the ambiguous outcome. (b) Example holograms used for the MPLC. (c) Amplitudes of non-orthogonal modes as they propagate through the MPLC. SLM: Spatial Light Modulator.

Unambiguous state discrimination was simulated for single-photon states of light encoded in high dimensions using single-outcome projective measurements <cit.>.  Such projective measurements provide a limited functionality in quantum communication systems, where multi-outcome measurements are necessary for maximising key rates and loophole-free tests of Bell non-locality <cit.>.  Recent work on state discrimination includes sorting high-dimensional states using optimal measurement strategies <cit.>, multi-state quantum discrimination through optical networks <cit.>, and quantum state elimination <cit.>.  Additionally, much of this recent work is related to the complementary field of classical deep optical networks used for information processing and image classification using diffractive optics <cit.>.  Methods for sorting and manipulating states using bulk optics <cit.>, two phase screens <cit.>, complex media <cit.>, and multi-plane light converters (MPLCs) <cit.> have been the topic of significant recent research.

The problem we set out to solve in this work is the simultaneous sorting of $d$ equally overlapping $d$-dimensional quantum states of light $\{\ket{\psi_1},...,\ket{\psi_d}\}$, i.e. with all pairwise fidelities of these states equal to each other $F = \left|\bra{\psi_i} \psi_{j \ne i}\rangle\right|^2$ <cit.>.  There are several requirements in order to perform this task. Firstly, every input mode in the set $\{\ket{\psi_1},...,\ket{\psi_d}\}$ is mapped to an individual measurement mode $\{ \ket{1},...,\ket{d} \}$ that uniquely identifies it, i.e. when a photon in state $\ket{\psi_1}$ passes through the system, only the measurement mode $\ket{1}$ can “click.”  The challenge is to perform this sorting process when every state in the set has a non-zero fidelity with respect to every other state in the set, i.e. $0< \left|\braket{\psi_i | \psi_{j \ne i}}\right|^2 < 1$.  Secondly, the system should not make any errors and incorrectly identify any input state, i.e. for the $\ket{\psi_1}$ input state, the probability of all output modes other than $\ket{1}$ clicking should be equal to zero.   Finally, the sorting process occurs simultaneously, and with the highest possible probability, which for USD is given by $\eta=1-|\bra{\psi_i} \psi_{j \ne i} \rangle|^2$ <cit.>.

MPLC for USD:– We realise unambiguous state discrimination (USD) within the framework of a multi-plane light converter (MPLC) as shown in Fig. <ref>.  This allows us to simultaneously transform a set of non-orthogonal states in any basis into a new basis of spatially separated measurement modes that can be simultaneously detected by a position-resolving detector. Here, the non-orthogonal states of light we sort are super-positions of Hermite-Gaussian (HG) modes, and we design the output modes of the MPLC to be spatially separated Gaussian spots that can be directly read by a camera or a single-photon detector array.  The MPLC is programmed to perform the required unitary (USD operation) and the necessary mode conversion (HG $\rightarrow$ Gaussian spots) at the same time.

The holograms used in MPLC devices are typically constructed using an inverse-design technique known as wavefront matching  <cit.>. This is an iterative algorithm where the optical fields for the input and output modes are forward and backward propagated respectively and overlapped at each plane of the MPLC. The phase of the MPLC at each plane is calculated in such a way that the entire set of input modes are phase-matched to the respective output modes. This process is repeated for each plane sequentially until the algorithm converges and the difference between the forward and backwards propagating light is minimised.  Each reflection from a mask performs a phase-only transformation $\hat{\Phi}_i$ of the input states of light, which is followed by free-space propagation ${\hat{H}}$ to the next plane. The total operator of the device after $n$ reflections is given by $\hat{U} = \hat{H} \prod_{i=n}^1\left(\hat{\Phi}_i\hat{H} \right)$. The free-space propagation operator ${\hat{H}}$ is calculated by simulating light propagation in free-space as discussed in detail in Supplementary Information.

Here we introduce an additional output state labelled $\ket{?}$ into the wavefront matching protocol, which means that for $d$ input modes, there are now $d+1$ output modes. The MPLC device supports a large number of modes, which is ultimately limited by the number of pixels and spatial resolution of the phase masks. The number of modes of the MPLC greatly exceeds the number of modes we sort, which allows us to include the additional auxiliary mode in a straightforward manner.  We are free to choose $d$ input modes (superpositions of HG modes) and $d+1$ output modes (Gaussian spots placed symmetrically on the circumference of a circle). The wavefront matching technique then ensures the correct mapping between the two sets of modes.

The purpose of the $\ket{?}$ output is that if a photon is detected in this mode, it provides no information about the input state. However, this also implies that we have not made any incorrect identification, as required by USD.  As we can successfully discriminate between any two states with a probability of $1-F$, the probability that $\ket{?}$ clicks is equal to $F$. The key to the success of this protocol is that this operator is designed to perform both a unitary operation that maps a set of input states onto the required unambiguous measurement states and simultaneously changes the basis for the measurement outcomes to be spatially separated. Both of these transformations are performed concurrently within the MPLC.

Measurement data for simultaneous USD in $d = 7$ using an MPLC. (a) Measurement matrices for a range of fidelities ($F$) between the input states.  The insets show the numerically modelled results. The results are normalised by the total power detected in each output state. (b) Probability of measuring an input state $\ket{\psi_1}$ in a given output state $\ket{x}$ as a function of inter-state fidelity. This is a cross section of the measurement matrices for the $\ket{\psi_1}$ input state.  The points are the measured values, which should be compared to the theoretical predictions, indicated by the dashed lines.

Evaluation of the MPLC for USD compared to the theoretical limit of MESD. The grey shaded area indicates the error-fidelity region accessible via MESD. Points outside this region represent USD with an error rate that is below than what is possible with MESD.

Results:– We performed unambiguous state discrimination for sets of symmetric non-orthogonal states constructed from modes in the Hermite-Gaussian (HG) basis.   In our experiment, we generated these modes with a HeNe laser and a spatial light modulator, and detected them with a CMOS camera.

In each realisation of USD for high-dimensional states of light, $d$ non-orthogonal modes were transformed into $d+1$ output modes and measured simultaneously on a CMOS camera.  The intensities recorded by the camera pixels located in the output modes were integrated and converted to a detection probability.  We achieved USD for sets of states in dimensions ranging from $d = 3$ to $7$ for fidelities in the range $F(\psi_i, \psi_{j \ne i}) \in (0, 1)$. Each state was generated as a complex superposition of Hermite-Gauss modes while varying the inter-state fidelity from 0 to 1 (see Supplementary Information for further details).   The results for USD in 7 dimensions are displayed in the form of correlation matrices in Fig. <ref>. Each row provides its detection probability in all possible outputs.   Fig. <ref>b compares the experimentally measured probabilities of successful USD to the theoretical predictions. This data clearly demonstrates that unambiguous sorting is achieved for 7-dimensional overlapping states of light.

To quantitatively assess our system, we analyse the performance of the MPLC compared to the theoretical limit of minimum error state discrimination (MESD) <cit.>. No auxiliary state is used in MESD in $d$ dimensions, and the error is instead distributed between the $d$ output states. This leads to a probability of error ($p_{\mathrm{err}}$) in the output, which is defined as the probability of measuring any output state $\ket{j\neq i}$ when given an input state $\ket{\psi_i}$. MESD for uniform-fidelity states has a minimum possible error probability given by $p_{\mathrm{err}}\ge\frac{1}{2}\left(1-\sqrt{1-F(\psi_i, \psi_{j \ne i})}\right)$ <cit.>. In Fig. <ref>, we plot the measured error probabilities of our USD protocol against this MESD threshold (grey area). We see that our system outperforms MESD and has a lower error rate than any strategy using MESD over a wide range of inter-state fidelities (overlaps). The error probabilities are higher for lower $F$ because the MPLC struggles more when sorting states that are farther apart. Additionally, as $d$ increases, so too does the total error. The performance reduction of the MPLC  for increasing $d$ can be explained as the accuracy of such transformations has been shown to reduce as the dimensionality increases for a fixed number of phase planes <cit.>.  Here, we use four phase plane while increasing the dimensionality of the set of modes that we are sorting.   We note that the error rate that we observe is approximately two times higher that that of Agnew et al. <cit.>, where measurements were performed one outcome at a time, which necessarily includes $(d-1)/d$ amount of loss.  However, in this work, all outcomes were measured simultaneously, which is a significant advance over the prior work in terms of practical applications.

Theoretically, the error should be equal to zero for USD, yet we see the experimental implementation using the MPLC performs better for states with a higher initial fidelity (overlap).  The reason for this is that as $F$ increases, a higher fraction of the energy is put into the single ambiguous outcome.  This mapping, where all input states are sorted to a single output state, is a simpler task to achieve for the MPLC than the case where all input states are sorted to individual outcomes.    Additionally, we are comparing our experimental USD implementation to the theoretical limit of MESD.  We see that, consequently, there are is a region where the error probability does not fall below what MESD could achieve.  However, any practical implementation of high-dimensional MESD would be subject to similar error as our experiment and would not perform at this theoretical limit.

Figure <ref> shows the extension of our method to the sorting of overlapping images.   We sort three images depicting a smiley face, sad face, and neutral face that have a large and symmetric overlap ($F=0.34$) with respect to each other—the eyes in the images are the same, while the mouth expressions are slightly different, connoting completely different emotions. In simple optical image classification, the three faces would be transformed directly to three spatially separated spots, and the wavefront matching algorithm would attempt to direct all input light into all output modes, leading to imperfect classification.  When using the extra mode $\ket{?}$, we have a place to direct any overlapping light, potentially leading to no errors in our measurement state outputs.

The correlation matrix for the images can be reformatted to a confusion matrix by removing the $\ket{?}$ mode and renormalising the rows. The success probability for image classification is then calculated as the average ratio of the light intensity in the outcomes of interest compared to the total light intensity of all other outcomes, excluding the ambiguous outcome.  Despite the large overlap between the input images, the USD protocol enables the sorting and classification with an average accuracy of 97.6%. We see that the success probability is not constant across all three input states.  As numerical simulations suggest that an equal success probability can be obtained, we believe that this asymmetry is due to the combined experimental error of the generation of the input modes and the misalignment between phase masks  in the MPLC.

Sorting overlapping images with the MPLC. The input images have a large fidelity ($F=0.34$) with each other, yet we can sort them with an accuracy of 97.6%.  (a) Shows the measured intensities of the input and output modes. (b) Shows the confusion matrix of the sorting and the associated success probabilities.

Conclusions:– In this work, we solve the measurement problem of simultaneous sorting of high-dimensional non-orthogonal states of light.  While previous works simulate a POVM with consecutive measurements, this work simultaneously realises the outcomes of a POVM  for high-dimensional non-orthogonal states of light, extending the use of MPLCs to include non-orthogonal input modes. Such simultaneous measurements can also be performed by engineering bases of light that are spatially separable on the detector <cit.>, however in this work we harness the abilities of MPLC to perform unitary operations in arbitrary spatial bases to perform this task experimentally for the first time.

The key to the success of this protocol is the additional output mode that provides extra flexibility to the system and enables perfect mapping from every $d$ non-orthogonal input state to a unique output. The method we adopt is the optimal strategy for minimising errors and correctly identifying the input states.  The consequence of USD is that at the single-photon level, the sorting does not provide an outcome 100% of the time, and for intense modes of light, there is a reduction in the total power that is transmitted into the known output modes.  Future work will focus on extending our method to sort quantum states and modes that have a non-uniform fidelity with respect to each other, as well as applying such generalised measurement strategies on high-dimensional entangled states of light.

Although this experiment was performed with coherent states of light, the theoretical formalism still applies at the single-photon level since, in the case of linear optics, there is an equivalence between the probabilities associated with single-photon detection and classical laser fields <cit.>.  The extension of our work, however, to include multi-photon states is of great interest, and in this case, multi-photon interference can lead to photon bunching or anti-bunching.  For multi-photon inputs, it would be necessary to measure the outcomes of our sorter in coincidence to fully reveal the photon statistics <cit.>.

Acknowledgements:–We thank Will McCutcheon for fruitful discussions regarding this work. This work was supported by EPSRC grants EP/T00097X/1 and EP/P024114/1, by QuantERA ERA-NET Co-fund (FWF Project I3773-N36), and the European Research Council (ERC) Starting grant PIQUaNT (950402).

Unambiguous State Discrimination:– In the process of unambiguous state discrimination (USD), our task is to accurately identify a random state from a set of states $S=
\left\{
\psik{i}
\right\}_i$. In general, this set of states $S$ are non-orthogonal with each other. If two states have an overlap along a vector, when a state measurement collapses into this vector, we get no information of which state we started with. Since this measurement needs to be accurate, we denote this outcome as an ambiguous outcome. Alternatively if the result of the collapse is a vector which only has overlap with one of our initial states, then we can be sure that was the state we started with. The vectors we measure with are denoted $\ket{i}$ for the state which only has overlap with $\psik{i}$ and $\ket{?_j}$ for the state with multiple non-zero overlaps. This set of states forms a projective operator-valued measure (POVM). In general, there could be different overlaps between different states in $S$ corresponding to different ambiguous outcomes $\ket{?_j}$, however in this work we looked at symmetric states which have equal overlap with each other corresponding to a single ambiguous outcome denoted by $\ket{?}$.

Generation of symmetric states:– A set of $d$ $d$-dimensional symmetric states $\left\{\psik{i}\right\}_i$ is defined such that
\braket{\psi_i | \psi_j}
\right|^2 = \delta_{ij} + (1-\delta_{ij})\left|\beta\right|^2$, where $\beta$ is the overlap between any two different states in the set. To be able to perform our USD measurements, we need a set of orthogonal measurement states ($\ket{i}$ and $\ket{?}$). To find these, we first construct the set of states $
\left\{
\ket{i}
\right\}_i$ where $\ket{i}$ is orthogonal to every state but $\psik{i}$ ($\braket{i|\psi_j}=\delta_{ij}\alpha$). These states form an orthogonal basis for the initial set of states in $d$ dimensions. We then increase the dimension of the space to $d' = d+1$ by adding in an extra basis element $\ket{?}$ which is orthogonal to every $\ket{i}$. We now only need to find the transformation $\hat{V}$ which takes $\psik{i}$ ($= \sum_j\psi_{ij}\ket{j}=\alpha\ket{i}+\sum_{j\neq i}\psi_{ij}\ket{j}$) to $\ket{\psi_i'}=\hat{V}\psik{i}=\alpha\ket{i} + \beta \ket{?}$.

For this experiment, we constructed a set of states with a parameter $\theta$ which allowed us to control the fidelity. We start by building a set of $d$ states with overlap:

\begin{equation}
\braket{\psi'_i|\psi'_j} = -\frac{1}{d-1}
\end{equation}

These are generated with an iterative technique similar to Gram-Schmidt orthonormalisation. This theoretical technique is a generalization of the problem in three dimensions <cit.> to $d$-dimensions.
To generate the set of $d$ vectors $\left\{\ket{\psi_i}\right\}_i$ in $d-1$ dimensions:
Let $\psi_i^j$ be the $j$'th component of the $i$'th vector.

Find the components as:

\begin{equation}
\psi_i^j =
\begin{cases}
1, & \text{when }~i=j=1 \\
\frac{-\frac{1}{d-1} - \sum\limits_{a < j}\psi_i^a\psi_j^{*a}}{\psi_j^j}, & \text{when } j < i \text{ or } j+1 = i=d \\
\sqrt{1-\sum\limits_{a < j} \psi_i^a\psi_i^{*a}}, & \text{when } j = i \\
0, & \text{otherwise}
\end{cases}
\end{equation}

where $\ket{\psi'_i} = \sum_{j = 1}^{d} \ket{j}$.
Once these states are made, we mix them into the additional dimension by:

\begin{equation}
\ket{\psi_i} = \sin(\theta)\ket{\psi'_i} + \cos(\theta)\ket{d+1}.
\end{equation}

Experimentally, the states $\{ \ket{1},\ket{2}...\ket{d+1}\}$ are chosen from the Hermite-Gauss basis. Thus each input state $\ket{\psi_i}$  is prepared as a superposition of $d+1$ Hermite-Gauss modes. For a given dimension of the extended Hilbert space $d+1$, we choose a certain set of  the Hermite-Gauss modes such that they follow $n+m+1=d+1$, where $m$ and $n$ are number of modes in either direction.

Sorting non-orthogonal states:– In order to perform USD measurements on these set of $d$ $d$-dimensional symmetric states $S = \left\{\psik{i}\right\}_i$, we need to define a set of measurement states. To do so, first calculate  $\ket{\psi_i^\bot}$ , which are orthogonal to every state but $\ket{\psi_i}$. This can be done by taking the SVD of all the $\ket{\psi_{j\neq i}}$ stacked as a matrix, and removing the components of $\ket{\psi}$ from the subspace this matrix spans. We then define the set of measurement states $\left\{D_i\right\}$ by expanding the Hilbert-space into a higher dimension using $\psik{d+1}$ which is orthogonal to every state in $S$ as

\begin{equation}
D_i = \ket{\psi_i^\bot} + \sqrt{-\braket{\psi_1^\bot|\psi_2^\bot}}\ket{d+1}
\end{equation}

and then normalise it. The overlap $\braket{\psi_i^\bot|\psi_{j\neq i}^\bot}=\braket{\psi_1^\bot|\psi_2^\bot}$ is the same for any two states $i≠j$. The unknown state is included by Gram-Schmidt orthogonalisation with this set of states. We also need to extend our original states $|ψ_i⟩$ into this space, and we do this just by adding the extra basis element $|ψ_d+1⟩$ with $0$ as its coefficient.

The final sorting is achieved using multi-plane light conversion that converts the set of states ${D_i}$  into spatially separated outcomes on a camera.

\emph{Multi Plane Light Conversion}:--  The multi-plane light converter (shown in Fig.~\ref{fig:mplc_diagram}) is a device which is trained to do a specific transformation on input light. It is formed of multiple planes in series which each impart a phase change on the light-field, and this phase change combined with the propagation of light between planes can approximate unitary transformations on the input light.

\begin{figure}
\centering
\begin{tikzpicture}[scale=2]

\tikzmath{\lighty = 0.7; \xwidth = 1;}
\foreach \xoffset in {0, 1, 2, 3}
\tikzmath{\planename = \xoffset -1;}
\begin{scope}[cm={0.4,0.3,0,1,(\xoffset * \xwidth,0)}]
\node[transform shape] at (0.2,0.3) {\includegraphics[width=.04\textwidth]{  fig_mask_\xoffset}};
\end{scope}

\begin{scope}[shift={(-1 * \xwidth, 0)}]
% \draw[] (0, 0) -- (0, 0.6) -- (0.1, 0.8) -- (0.1, 0.2) -- (0, 0);
\draw[->] (0.6, 0.4) -- (0.8, 0.4) node [midway, above] (TextNode) {$\hat{H}$};

\end{scope}

\foreach \xoffset in {1, 2}
\tikzmath{\planename = \xoffset -1;}
\begin{scope}[shift={(\planename * \xwidth, 0)}]
% \draw[] (0, 0) -- (0, 0.6) -- (0.1, 0.8) -- (0.1, 0.2) -- (0, 0);
\draw[->] (0.4, 0.4) -- (0.8, 0.4) node [midway, above] (TextNode) {$\hat{H}$};
\node at (0.1, 1) {$\hat{\Phi}_{\xoffset}$};

\end{scope}
\begin{scope}[shift={(2, 0)}]
% \draw[] (0, 0) -- (0, 0.6) -- (0.1, 0.8) -- (0.1, 0.2) -- (0, 0);
\draw[->] (0.4, 0.4) -- (0.8, 0.4) node [midway, above] (TextNode) {$\hat{H}$};
\node at (0.1, 1) {$\hat{\Phi}_{i}$};

\end{scope}

\begin{scope}[shift={(3 * \xwidth, 0)}]
% \draw[] (0, 0) -- (0, 0.6) -- (0.1, 0.8) -- (0.1, 0.2) -- (0, 0);
\draw[->] (0.4, 0.4) -- (0.6, 0.4) node [midway, above] (TextNode) {$\hat{H}$};

\end{scope}

\draw[decorate,thick,decoration={calligraphic brace,amplitude=3pt}] (1.9, -0.1) -- (1.9, 1.2);
\draw[decorate,thick,decoration={calligraphic brace,amplitude=3pt}] (2.8, 1.2) -- (2.8, -0.1);
node [end, right] (TextNode) {$i$};

\begin{scope}[shift={(3, 0)}]
% \draw[thin] (0, 0) -- (0, 0.6) -- (0.1, 0.8) -- (0.1, 0.2) -- (0, 0);
\node at (0.1, 1) {$\hat{\Phi}_{n}$};
\end{scope}

\node at (0, 1.3) (nodeA){};
\node at (1, 1.3) (nodeB){};
\draw[<->,dashed] (nodeA) -- (nodeB) node [midway, above] (TextNode) {plane spacing};
\draw [decorate,very thick,decoration={calligraphic brace,mirror,amplitude=10pt}] (-0.3, -0.3) -- (3.3, -0.3);

\node at (1.65, -0.7) {$\hat{U} =
\hat{H}
\prod_{i=n}^1\left(\hat{\Phi}_i\hat{H}
\right)$};
\end{tikzpicture}
\caption{ Illustrative diagram of MPLC: The unitary operation given as a product of multiple free-space propagation operators and phase planes.  }
\label{fig:mplc_diagram}
\end{figure}

Plane $i$ performs a phase-only transformation $Φ̂_i$ to the light, which is followed by propagation to the next plane $Ĥ$, giving the matrix of the device: $Û =

The calculation of the masks is done with the wavefront matching method.

\emph{Simulation of Light Propagation}:-- To simulate light propagation, we use the approximation of a 3-dimensional scalar field $ψ(𝐱)$ \cite{siegman1986lasers,khare2015fourier}, and define an initial field over a plane $P$ orthogonal to the light's main propagation direction\ $𝐧̂$: $
P = {
𝐱 |(𝐱 - 𝐚) ·𝐧̂ = 0
}$. The field is then propagated a distance $α$ along the vector $𝐧̂$ by the operator  $Ĥ(α)$
%$e^{-i\hat{H}\alpha}$ to give its value for all points $\mathbf{x}$.

To approximate this calculation on the computer, we first discretise the field $ψ(𝐱)$ over the plane P into the vector $ψ_ij$, where the indices $i$ and $j$ index two orthogonal directions. We take the Fourier transform of this vector in both of these indices to give $ψ̃=ℱ(ψ)$. The phase change over a distance $α$ along the $𝐧̂$ direction for light of spatial frequency $
| = 2π/λ$ is given by:

\begin{equation}
\label{eq:phase_change}
\exp
\left(
\right)
= \exp
\left(
\left(
\frac{2\pi}{\lambda}
\right)^2-
\left|
\right|^2}\alpha
\right)
\end{equation}
Where the value $|
|^2$ is sum of the squares of the coordinates $ψ̃$ is defined over. The inverse Fourier transform is then taken of this product to get the field at the new plane $α𝐧̂$ from the initial:
\begin{equation}
\label{eq:light_simulation}
\psi(\alpha) = \mathcal{F}^{-1}
\left\{
\exp
\left(
\right)\mathcal{F}
\left(
\psi_0
\right)
\right\}
\end{equation}

The propagation operator $Ĥ(α)$ is then given by:

\begin{equation}
\label{eq:propagation_operator}
\hat{H}(\alpha) = \hat{F_2}^\dagger \hat{D}(\alpha)\hat{F_2}
\end{equation}

where $F̂_̂2̂$ is the operator that takes 2D-Fourier transform of the field and  $D̂(̂α̂)̂$ is the operator that adds the distance dependent spatial phase to the field.

\emph{Wavefront Matching}:--  In order to encode a transformation within an MPLC, we need to calculate the right patterns to be displayed on each plane. This is done by a process called ``wavefront matching''.  It consists of propagating the input modes forward through the planes, and the output modes backward through the planes. At the plane where the modes meet, we create a new mask as the average overlap over all modes.

If our input modes at plane $p$ are $a_i, p(x)$, and our output modes are $b_i, p(x)$ then the expression for the next plane mask is given by:
\begin{widetext}

\begin{equation}
\Phi_p'(x) = \sum_i a_{i,p}(x)b_{i,p}^*(x)\exp
\left(
\left[
\int dx \Phi_p^*(x) a_{i,p}(x)b_{i,p}^*(x)
\right]
\right)
\end{equation}

\end{widetext}
The forward, and backward propagating fields at plane $p$ are given by $a_i,p(x)$ and $b_i,p(x)$ respectively.
The goal is to have the fields exactly matching at those points, so that $⟨a_i,p|Φ_p^†|b_i,p|=⟩1$ over all input modes $i$.
The change in the phase mask can be written in terms of its phase:
\begin{align}
\label{eq:new_phase_mask}
\Phi'(x) &= e^{i(\theta(x) + \delta \theta(x))}  \\
&\approx e^{i\theta(x)}
\left(
\right) \\
&= \Phi(x)
\left(
1 + i\delta\theta(x)
\right)
\end{align}

The new overlap is given by:
\begin{equation}
\label{eq:new_overlap}
\braket{a_{i,p}|\Phi'^{\dagger}|b_{i,p}}
\end{equation}

We want the new $Φ'$ to increase the overlap:
\begin{widetext}

\begin{align}
\label{eq:overlap}
\eta' &=
\left|
\braket{a_{i,p}|\Phi'^{\dagger}|b_{i,p}}
\right|^{2} \\
&\approx   \left|\braket{a_{i,p}|\Phi^{\dagger}(1-i\delta\theta^{\dagger})|b_{i,p}}
\right|^{2} \\
&=   \left|
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}} -i \braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}}
\right|^{2} \\
&= \left(
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}} -i \braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}}
\right)
\overline{\left(
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}} -i \braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}}
\right)} \\
&= \eta + 2\Re\left\{\
\overline{\braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}} }
\right\} + \mathcal{O}(\delta\theta^{2}) \\
&\approx \eta -2\Im
\left\{
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}}
\overline{\braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}} }
\right\}
\end{align}

\end{widetext}
And so to satisfy $η'>η$, we just need:
\begin{equation}
\label{eq:eta_greater}
\Im
\left\{
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}}
\overline{\braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}} }
\right\} < 0
\end{equation}
We let $α=⟨a_i,p|Φ^†|b_i,p|$⟩, and separate $\overline{\braket{a_{i,p}|\Phi^{\dagger}\delta\theta^{\dagger}|b_{i,p}} }$ into its components, requiring the inequality to be satisfied for each inidividual component:

\begin{equation}
\label{eq:component_wise_inequality}
\Im
\left\{
\alpha
\right\} < 0
\end{equation}

And letting $\gamma(x)=\alpha a_{i,p}(x)e^{i\theta(x)}b_{i,p}^{\dagger}(x)$:

\begin{equation}
\label{eq:gamma_sub}
\delta\theta(x)\Im
\left\{
\gamma(x)
\right\}<0
\end{equation}

By letting $\delta\theta(x)=-\Im
\left\{
\gamma(x)
\right\}$, we can satisfy this equation.
So, to update the phase mask:

\begin{equation}
\label{eq:phase_mask_update}
\Phi'(x) = \sum_i e^{i\theta(x)} e^{
\left(
\left\{
\braket{a_{i,p}|\Phi^{\dagger}|b_{i,p}}
\right\}
\right)}
\end{equation}

Addition of auxiliary mode:– Here we provide some explanation as to how the auxiliary mode is added to the system. There are an infinite number of free-space optical modes, but modes that are accessible in this experiment are limited to the Hilbert space of the MPLC.  This space is limited by the area and spatial frequencies of the masks on the SLM, which is set by the number of pixels on the SLM and the numerical aperture of the collection optics.  For example, if we want to sort a large number of input modes, the displayed masks on the SLM need to have a large range of spatial frequencies, which, in turn, requires many pixels.

It is the large number of modes provided by the MPLC enables the addition of the auxiliary mode.  We are free to choose where the outcomes of the MPLC are directed to, and we choose a spatially separated circular arrangement of spots.

In order to sort a greater number of mode modes, we would require larger masks and also larger spatial frequencies along with a greater number of phase planes.  However, there are physical constraints to this, for example, the propagation distance from one mask to the next, and the physical size the SLM that is used to display all the masks.   In our case, each mask is of dimension 250$\times$250 pixels, and we use four masks.  This enables us to sort up to seven non-orthogonal modes. For a comprehensive analysis of the trade-offs and scalability of multi-plane optical circuits in this context, see ref <cit.>.

Experimental Setup:–
As illustrated in  Fig. <ref>, the experimental setup consists of three sections: generation, conversion and measurement. The desired non-orthogonal states are generated at a first SLM are then unambiguously discriminated at the second SLM, which acts as a multi-plane light converter (MPLC). The intensity distribution of the output states are then measured using a camera.

For the generation of the input states, we first use a 5 mW He-Ne Laser to produce a linearly polarized beam at 633 nm. This is then passed though a single-mode fiber and is collimated with a pair of lenses with focal lengths 45 mm and 100 mm. This beam then passes through a half-wave plate (HWP) such that it is polarized in the optical axis of the generation SLM (Holoeye Pluto 2.1). Our input modes are generated using computer generated holograms (CGH), which are calculated using the type 3 method as detailed in <cit.>. This beam is imaged to the input plane of the MPLC using two lenses of focal length 750 mm and 500 mm.  We use a Fourier filter at the focus of the 750 mm lens to select generated beams.

The MPLC is build in the reflection mode with 4 phase planes comprising of a phase-only SLM (Meadowlark E19x12) and a mirror, separated by 17 mm. Inside the MPLC, light undergo a transformation due to multiple reflection from the SLM followed by free-space propagation. Each reflection within the MPLC changes the phase profile of the mode. A grating of 10-pixel period is displayed across the whole SLM in order to efficiently separate modulated light from the residual. The target input and output modes are assigned at $17~$mm propagation distance before the first and after the last reflection from the SLM.

Finally, the intensity distributions of the output modes are measured via the CMOS Camera (Thorlabs USB3.0) using a 4-$f$ system with a focal length of 250 mm. Since our beam is reflected though the same grating 4 times, at the focal plane of this imaging system, we filter the $4^th$ order diffraction that carries the fully transformed modes.

Schematic of the experiment: Input modes are generated using CGH holograms at SLM 1 (a).  The light from SLM 1 is filtered and imaged to the MPLC where 4 phase masks with a constant grating (b-e) are displayed. The output modes of the MPLC are filtered and imaged to a CMOS camera where light is sorted into Gaussian spots (f). Lenses L1-6 have focal length 45mm, 100mm, 750mm, 500mm, 250mm and 250mm, respectively. HWP 1 and HWP 2 are half-wave plates.

Data Processing:– We measure the experimental coupling matrices between the set of input modes and the output modes by measuring the intensity on each Gaussian spot. There are discrepancies between the ideal and measured intensities. This is due to a variety of factors like mode-dependent loss, imperfect wavefront matching, misalignments and subnormalised input state generation using computer-generated holograms (CGH).

We correct for these imperfections by multiplying the measured matrices with a “correction vector”. Experimentally, this can be interpreted as applying a non-uniform attenuation to different ouput modes. If our ideal USD matrix is $M$, and our experimentally measured USD matrix is $E$, we can construct a correction vector $\mathbf{v}$ with elements:

\begin{equation}
\label{eq:correction_vector_components}
\mathbf{v} =
\begin{pmatrix}
\frac{E_{11}}{M_{11}} & \frac{E_{22}}{M_{22}} & \dots & \frac{E_{dd}}{M_{dd}} & \frac{\sum_{i=1}^d{E_{i(d+1)}}}{\sum_{i=1}^d{M_{i(d+1)}}}
\end{pmatrix}
\end{equation}

Which can correct our measurements by multiplying column-wise:

\begin{equation}
\label{eq:measurement_correction}
M'_{ij} = M_{ij}\mathbf{v}_{j}
\end{equation}

Correlation matrices for all dimensions: Correlation matrices of USD measurements in dimensions $d=2$ to $d=8$ for various input state fidelities ranging from $F=0.2$ to $F=0.9$

Confusion matrices for all dimensions: Confusion matrices of USD measurements in dimensions $d=2$ to $d=8$ for various input state fidelities ranging from $F=0.2$ to $F=0.9$