Let $H:(f_{e})\mapsto(-f^{\prime\prime}_{e}+Uf_{e})$ be the self-adjoint
operator acting on
$(f_{e})\in\bigoplus_{e\in\mathcal{E}}{\mathcal{H}}^{2}[0,1]$ and satisfying
the conditions (5.10), for all $v\in\mathcal{V}$. Making some adaptations to
the argument set out above, it can be shown that the Theorem 5.1.1 is valid in
the same way if we replace the operator $\Delta$ with its magnetic version
$\Delta_{\beta}$ given by (5.7). In particular, this construction can be
applied to a two-dimensional lattice with a uniform magnetic field.

Furthermore, in [29] it was observed that the $H$ operator can be studied at
an abstract level using the boundary triple language and self-adjoint
extensions. More precisely, using boundary triples, Pankrashkin related
spectral properties of the continuous Laplacian operator $H$ in
$\mathcal{H}=\bigoplus\mathcal{H}^{2}[0,1]$ with the discrete operator
$\Delta_{\beta}$ in $l^{2}(G)$. Or, the author demonstrates that

$\sigma(H)\backslash\sigma_{D}=\Big{\\{}z\in\mathbb{R}\backslash\sigma_{D}:\mu_{\alpha}(z)\in\sigma(\Delta_{\beta})\Big{\\}}=\mu^{-1}_{\alpha}\big{(}\sigma(\Delta_{\beta})\big{)}\backslash\sigma_{D},$
(5.11)

where $\mu_{\alpha}$ o discriminante de Hill

$\mu_{\alpha}(z)=c(1,z)+\frac{\alpha}{2}s(1,z)\,.$ (5.12)

In other words, for $z\in\mathbb{R}\backslash\sigma_{D}$, Pankrashkin then
obtains the relationship between the spectra of $H$ and $\Delta_{\beta}$:

$z\in\sigma(H)\Longleftrightarrow\mu_{\alpha}(z)\in\sigma(\Delta_{\beta}),$
(5.13)

######  Observation 5.1.1.

As shown in [29], the analysis of the $H$ operator above can be done using
frontier triples; more precisely through the functions $c(\cdot,z),s(\cdot,z$)
associated with the operator $H$ and the constant $\alpha$ of the boundary
condition. The relationship determined in (5.13) is obtained under the
hypothesis that the symmetry condition remains and the coupling constants
comply with $\alpha(v)=\frac{\alpha}{2}\operatorname{deg}(v)$, with
$\alpha\in\mathbb{R}$. However, if one of these conditions fails (for example,
if the orientation of a single edge changes), one cannot relate the spectra of
the continuous and discrete problems as seen above.

In [29], Pankrashkin shows that the spectrum of $H$ has a band structure if
and only if the spectrum of $\Delta_{\beta}$ has a band structure, that is,
the spectrum and the union of a locally finite family of segments, while in
[30] he shows that the continuous Laplacian $H$ and the discrete Schrödinger
operator $\Delta_{\beta}$ are unitarily equivalent, that is, they have the
same spectra.

Based on these results, we will then apply them to the study of monolayers of
hexagonal structures, such as graphene and boron nitride, with rational
magnetic flux $\phi=2\pi\frac{p}{q}$ in a hexagonal lattice hexagon. We
reinforce that the quantum graph model for graphene with parameter
$\alpha(v)=0$ for every vertex $v$, in certain magnetic fields, has already
been studied in detail in [2, 3], and our contribution here it is the case of
graphene with $\alpha(v)$ constant, but not null, and also of hBN, always with
these fields.

###  5.2 Graphene and boron nitride with magnetic field

Let us consider the magnetic flux $\phi$ as a rational multiple of $2\pi$,
that is, $\phi=2\pi\frac{p}{q}$. Therefore, the operator $\Delta_{\beta}$ is
periodic [2]. Furthermore, given
$\theta=(\theta_{1},\theta_{2})\in[-\pi,\pi)\times[-\pi,\pi)$, we consider
$\Delta_{\beta}$ acting on $l^{2}(\hexagon)$ subject to pseudo-periodic
condition

$u(\gamma+qE_{k},v_{{}_{j-1}})=e^{i\theta_{k}}u(\gamma,v_{{}_{j-1}})\,,\,\,j,k=1,2\,,$
(5.14)

where $E_{1},E_{2}$ are base vectors of $\hexagon$ and $v_{0},v_{1}$ are
vertices of the fundamental domain $W_{\\!\hexagon}$.

We also know, from [2], that given the discrete operator $\Delta_{\beta}$ with
flow $\phi$ defined in (5.7) and acting on $l^{2}(\mathbb{Z}^{2},\ C^{2}),$
the Floquet matrix of $\Delta_{\beta}$ is given by

$M_{F}(\theta)=\frac{1}{3}\left(\begin{array}[]{cc}0&I_{q}+e^{i\theta_{1}}J_{p,q}+e^{i\theta_{2}}K_{q}\\\
I_{q}+e^{-i\theta_{1}}J^{*}_{p,q}+e^{-i\theta_{2}}K^{*}_{q}&0\end{array}\right)\,,$
(5.15)

where $J_{p,q}$ and $K_{q}$ are $q\times q$ matrices given by

$J_{p,q}=\operatorname{diag}\Big{(}\\{e^{i(j-1)\phi}\\}^{q}_{j=1}\Big{)}$

and

$(K_{q})_{jk}=\left\\{\begin{array}[]{ll}1\,,\,k=(j+1)\operatorname{mod}(q)\\\
0\,,\,\mbox{otherwise}\end{array}\right.\,.$ (5.16)

In particular, for $p=1,\,q=2$, we have

$M_{F}\big{(}\theta_{1},\theta_{2}\big{)}=\frac{1}{3}\left(\begin{array}[]{cccc}0&0&1+e^{i\theta_{1}}&e^{i\theta_{2}}\\\
0&0&e^{i\theta_{2}}&e^{i(\theta_{1}+\phi)}\\\
1+e^{-i\theta_{1}}&e^{-i\theta_{2}}&0&0\\\
e^{-i\theta_{2}}&e^{i(\theta_{1}+\phi)}&0&0\end{array}\right)\,.$ (5.17)

If the quasi-moment
$\theta=(\theta_{1},\theta_{2})\in[-\pi,\pi)\times[-\pi,\pi)$ is restricted to
$[0,\frac{\pi}{q})\times[-\frac{\pi}{q},\frac{\pi}{q})$, the authors
demonstrate that the scattering surface of $\Delta_{\beta}$ has Dirac at
energy level $0$, for any rational $p$ and $q$. To do so, they use the elegant
Chambers formula (see [20]) to help calculate the determinant of the matrix
$M_{F}(\theta)$ associated with the operator, since this matrix is of order
$q\times q$.

Note that, in the boundary conditions in (5.10), the same value $\alpha(v)$
was admitted at all vertices $v\in\mathcal{V}$ of the network; furthermore, in
the relation (5.13) such an assumption was also made, that is,
$\alpha(v)=\frac{\alpha}{2}\operatorname{deg}(v)$, with $\alpha\in\mathbb{R}$
constant, since, in a hexagonal network ⎔, $\ \operatorname{deg}(v)=3$ for any
vertex of the network. However, we will show that if we admit a certain
variation of the flow at the vertices, that is, $\alpha(v)$ is not constant,
the operator Dirac cones can disappear, as we present below.

  * •

If $\alpha(v)=0$ for every vertex $v$, we know that

$\sigma(H)\backslash\sigma_{D}=\\{z\notin\sigma_{D}:\mu_{0}(z)=\cos(\sqrt{z})\in\sigma(\Delta_{\beta})\\}\,,$

with $\sigma_{D}=\\{(\pi n)^{2}:n\in\mathbb{N}\\}$. With this, we see that
given a cone at a point in the spectrum of $H$, this cone is preserved by
applying $\mu_{0}(\cdot)=\cos(\cdot)$ in the spectrum of $H_{\beta}$ and, as
the Proposition 5.1.1, we conclude that, for any integer values $p,q$ such
that the magnetic flux has the form $\phi=2\pi\frac{p}{q}$, the operator
$\Delta_{\beta}$ has a Dirac cone, since its spectra are unitarily equivalent.

  * •

Similarly, for $\alpha(v)\neq 0$, and constant for all vertices, by Theorem
5.1.1, the presence of cones still remains for arbitrary values of $p,q$
(still with $\phi=2\pi\frac{p}{q}$). In fact, initially let us note that, for
$z\in\mathbb{R}\backslash\sigma_{D}$, the functions $c(\cdot,z),s(\cdot,z)$ in
(5.12) are of class $C^{1}$; therefore, the function $\mu_{\alpha}$ is also
$C^{1}$. Thus, taking a Dirac cone of $H$, simply show that its image by the
function

$\mu_{\alpha}(z)=c(\cdot,z)+\frac{\alpha}{2}\cdot s(\cdot,z)$

it is still a Dirac cone of $\Delta_{\beta}$. In fact, more generally, suppose
that $f(\theta)$ has a Dirac cone at $\theta_{0}$, that is, there is
$\gamma>0$, namely, $f^{\prime}(\theta_{0})\neq 0$, so that

$f(\theta)-f(\theta_{0})=\pm\gamma\cdot|\theta-\theta_{0}|+\mathcal{O}\big{(}|\theta-\theta_{0}|^{2}\big{)}\,.$
(5.18)

Now, if we perform the composition with another function $g$ of class $C^{1}$,
the cone remains if $g^{\prime}(f(\theta_{0}))\neq 0$. In fact, according to
the chain rule,

$g(y)-g(y_{0})=\pm
g^{\prime}(y_{0})|y-y_{0}|+\mathcal{O}\big{(}|y-y_{0}|^{2}\big{)}$

and choosing $y=f(\theta)$, $y_{0}=f(\theta_{0})$, we get

$g(f(\theta))-g(f(\theta_{0}))=\pm
g^{\prime}\big{(}f(\theta_{0})\big{)}\cdot|f(\theta)-f(\theta_{0})|+\mathcal{O}\big{(}|f(\theta)-f(\theta_{0})|^{2}\big{)}$

and using (5.18),

$g(f(\theta))-g(f(\theta_{0}))\\!=\\!\pm
g^{\prime}(f(\theta_{0}))\cdot\gamma\cdot|\theta-\theta_{0}|+\mathcal{O}(|\theta-\theta_{0}|^{2})\\!=\\!\gamma_{g}\cdot|\theta-\theta_{0}|+\mathcal{O}(|\theta-\theta_{0}|^{2})$

wint $\gamma_{g}:=g^{\prime}(f(\theta_{0}))\cdot\gamma\neq 0$, concluding that
the cone remains after composition with $g$, as long as
$g^{\prime}(f(\theta_{0}))\neq 0$; in particular, for $g=\mu_{\alpha}$, we
have $\gamma_{\mu_{\alpha}}\neq 0$ and the result follows.

  * •

However, if $\alpha(v)$ is not constant, corresponding to our modeling of
boron nitride, according to Observation 5.1.1, Pankrashkin does not guarantee
the relation (5.13), so that it is not possible to relate the spectra of
$\Delta_{\beta}$ in $l^{2}(\hexagon)$ and $H$ in
$\bigoplus_{e\in\mathcal{E}}{\mathcal{H}}^{2}[0,1]$, as this relationship was
obtained by supposing $\alpha(v)=\frac{\alpha}{2}\operatorname{deg}(v)$, with
$\alpha\in\mathbb{R}$ taking the same value on all vertices. We will then show
that for certain values of $\alpha$ the conical touches in the dispersion
relation of the operator $\Delta_{\beta}$ cease to exist so that we guarantee
the existence of gaps. To do so, we will use the Theorem 5.1.1 and the
relation (5.13) above, but we will base ourselves on the technique presented
in Chapter 1 (mainly due to Kuchment-Post [26] ), which we will adapt to the
case with a magnetic field in a specific configuration; the need to work with
a particular (and “small”) period comes from the lack of an analogue to the
Chambers Formula. Such adaptation and the example with rational flow in
relation to $2\pi$ without Dirac cones constitute our main contributions to
this model (and this chapter).

Let us then, from now on, fix $p=1$ and consider the pseudo-periodic condition

$u(\gamma+qE_{k},v_{j-1})=e^{i\theta_{k}}u(\gamma,v_{j-1})\,,\,\,j,k=1,2\,,$
(5.19)

for, initially, the period $q=1$ and, subsequently, $q=2$. Given an edge $e$
of the network $\hexagon$, let us take $\alpha(v)$ so that

$\alpha(v)=\left\\{\begin{array}[]{ll}\delta_{{}_{N}}\,,\,\mbox{if}\,v=i(e)=0\\\
\delta_{{}_{B}}\,,\,\mbox{if}\,v=t(e)=1\end{array}\right.\,.$ (5.20)

In the case $q=1$, the period of the hexagonal lattice ⎔ coincides with that
of the operator $H$ and, therefore, the modeling coincides with that already
covered in Chapter 1, for a sheet of nitride boron, so that the fundamental
domain is given by $W_{\\!\hexagon}=\\{f,g,h,v_{0},v_{1}\\}$, as shown in
Figure 5.1.

Figure 5.1: The fundamental domain $W_{\hexagon}$ and the generating vectors
$E_{1}$ and $E_{2}$ of the ⎔ network. For $q=1$, the period of the ⎔ network
coincides with that of the $H$ operator.

The continuity and flow equations at the vertices are given below.

$\mbox{Continuity at the
vertices:}\left\\{\begin{array}[]{cll}u_{f}(0)=u_{g}(0)=u_{h}(0)&=:&A_{1}\\\
e^{i\beta_{f}}u_{f}(1)=e^{i\beta_{g_{10}}}u_{g_{10}}(1)=e^{i\beta_{h_{01}}}u_{h_{01}}(1)&=:&B_{1}\end{array}\right.\,,$
(5.21)

so that, for the edges that are outside the fundamental domain, using the
cyclic Floquet conditions 111As the period of the operator and the lattice
coincide, the periodic Floquet relations and the pseudo-periodic ones in
(5.19) coincide., we obtain the relations
$u_{g_{10}}(1)=e^{i\theta_{1}}u_{g}(1)$ and
$u_{h_{01}}(1)=e^{i\theta_{2}}u_{h}(1)$;

$\mbox{Flux
Condition:}\left\\{\begin{array}[]{cll}u^{\prime}_{f}(0)+u^{\prime}_{g}(0)+u^{\prime}_{h}(0)&=&\delta_{{}_{N}}A_{1}\\\
-e^{i\beta_{f}}u^{\prime}_{f}(1)-e^{i\beta_{g_{10}}}u^{\prime}_{g_{10}}(1)-e^{i\beta_{h_{01}}}u^{\prime}_{h_{01}}(1)&=&\delta_{{}_{B}}B_{1}\end{array}\right.\,.$
(5.22)

Similar to what we saw in the initial chapters, for
$\lambda\notin\sigma(H^{D})$, there are functions $\varphi_{0},\varphi_{1}$
such that

$\left\\{\begin{array}[]{lll}u_{f}&=&A_{1}\varphi_{0}+B_{1}\varphi_{1}e^{-i\beta_{f}}\\\
u_{g}&=&A_{1}\varphi_{0}+B_{1}\varphi_{1}e^{-i(\theta_{1}+\beta_{g_{{}_{10}}})}\\\
u_{h}&=&A_{1}\varphi_{0}+B_{1}\varphi_{1}e^{-i(\theta_{2}+\beta_{h_{{}_{01}}})}\end{array}.\,\,\right.$
(5.23)

Combining the equations (5.22) and (5.23), with
$\alpha_{k}=\dfrac{\delta_{k}}{\varphi^{\prime}_{1}(0)}$, $k\in{B,N}$ , we
obtain the matrix representation

$M_{{}_{\hexagon}}(\theta)=\left(\begin{array}[]{cc}-3\mu-\alpha_{{}_{N}}&e^{-i\beta_{f}}+e^{-i(\theta_{2}+\beta_{h_{01}})}+e^{-i(\theta_{1}+\beta_{g_{10}})}\\\
e^{i\beta_{f}}+e^{i(\theta_{2}+\beta_{h_{01}})}+e^{i(\theta_{1}+\beta_{g_{10}})}&-3\mu-\alpha_{{}_{B}}\end{array}\right).$

Replacing222For simplicity, we will just use $u_{e}$ instead of
$u_{e_{{}_{00}}}$, to $e=f,g,h$; ditto for $\beta_{e}$.
$\beta_{f_{\gamma_{1},\gamma_{2}}}=\beta_{g_{\gamma_{1},\gamma_{2}}}=0$ and
$\beta_{h_{\gamma_{1},\gamma_{2}}}=-\gamma_{1}\phi$, with
$\phi=2\pi\frac{1}{1}=2\pi$, we get the polynomial
$P(\eta):=\det(M_{{}_{\hexagon}}(\theta))$ and the dispersion relation
$r^{\pm}(\theta)$ given by

$P(\eta)=9\eta^{2}+3(\alpha_{{}_{N}}+\alpha_{{}_{B}})\eta+\alpha_{{}_{B}}\alpha_{{}_{N}}-3-2\cos(\theta_{1})-2\cos(\theta_{2})-2\cos(\theta_{1}-\theta_{2})\,,$
$r^{\pm}(\theta)=\frac{-\alpha_{{}_{B}}-\alpha_{{}_{N}}\pm\sqrt{(\alpha_{{}_{B}}-\alpha_{{}_{N}})^{2}+4(3+2\cos(\theta_{1})+2\cos(\theta_{2})+2\cos(\theta_{1}-\theta_{2}))}}{6}\,,$

and we have already seen in the Chapter 1, Corollary 1.2.1, that this
dispersion relation has Dirac cones if $\alpha_{{}_{N}}=\alpha_{{}_{B}}$; see
Figure 5.2.

Figure 5.2: For $q=1$, the cone is maintained at
$(\theta_{1},\theta_{2})\in[0,\pi)\times[-\pi,\pi)$ for
$\alpha_{{}_{N}}=\alpha_{{}_{B}}$. Here $\alpha_{{}_{B}}=\alpha_{{}_{N}}=-1.$

In the case $q=2$, the period of the hexagonal lattice is 1, while the period
of the operator is $2$ and, therefore, we will have another modeling different
from that covered in Chapter 1, for a boron nitride sheet, as shown in Figure
5.3.

Figure 5.3: Period of the operator $\Delta_{\beta}$ for $q=2$.

In this case, the fundamental domain is given by

$W_{\hexagon}=\big{\\{}f,g,h,f_{10},g_{10},h_{10},v_{0},v_{1},v^{1,0}_{0},v^{1,0}_{1}\big{\\}}.$

The continuity and flow equations at the vertices are now given by

$\\!\\!\\!\mbox{\\!\\!\\!Continuity at the
vertices:}\left\\{\begin{array}[]{llll}\\!\\!\\!u_{f}(0)=u_{g}(0)=u_{h}(0)\\!\\!\\!&=:&\\!\\!\\!A_{1}\\\
\\!\\!\\!u_{f_{10}}(0)=u_{g_{10}}(0)=u_{h_{10}}(0)\\!\\!\\!&=:&\\!\\!\\!A_{2}\\\
\\!\\!\\!e^{i\beta_{f}}u_{f}(1)=e^{i\beta_{g_{10}}}u_{g_{10}}(1)=e^{i\beta_{h_{01}}}e^{i\theta_{2}}u_{h}(1)\\!\\!\\!&=:&\\!\\!\\!B_{1}\\\
\\!\\!\\!e^{i\beta_{f_{10}}}u_{f_{10}}(1)=e^{i\beta_{g_{20}}}e^{i\theta_{1}}u_{g}(1)=e^{i\beta_{h_{11}}}e^{i\theta_{2}}u_{h_{10}}(1)\\!\\!\\!&=:&\\!\\!\\!B_{2}\end{array}\right.$
(5.24)

and

$\mbox{Flux
Condition:}\left\\{\begin{array}[]{llll}u^{\prime}_{f_{00}}(0)+u^{\prime}_{g_{00}}(0)+u^{\prime}_{h_{00}}(0)&=&\delta_{{}_{N}}A_{1}\\\
-u^{\prime}_{f_{00}}(1)-u^{\prime}_{g_{10}}(1)-u^{\prime}_{h_{01}}(1)&=&\delta_{{}_{B}}B_{1}\\\
u^{\prime}_{f_{10}}(0)+u^{\prime}_{g_{10}}(0)+u^{\prime}_{h_{10}}(0)&=&\delta_{{}_{N}}A_{2}\\\
-u^{\prime}_{f_{10}}(1)-u^{\prime}_{g_{20}}(1)+u^{\prime}_{h_{11}}(1)&=&\delta_{{}_{B}}B_{2}\end{array}.\right.$
(5.25)

For $\lambda\notin\sigma(H^{D})$, we can write

$\left\\{\begin{array}[]{llllll}u_{f_{{}_{00}}}&=&A_{1}\varphi_{0}+B_{1}e^{-i\beta_{f_{{}_{00}}}}\varphi_{1}\\\
u_{g_{{}_{00}}}&=&A_{1}\varphi_{0}+B_{2}e^{-i(\theta_{1}+\beta_{g_{{}_{20}}})}\varphi_{1}\\\
u_{h_{{}_{00}}}&=&A_{1}\varphi_{0}+B_{1}e^{-i(\theta_{2}+\beta_{h_{{}_{01}}})}\varphi_{1}\\\
u_{f_{{}_{10}}}&=&A_{2}\varphi_{0}+B_{2}e^{-i\beta_{f_{{}_{10}}}}\varphi_{1}\\\
u_{g_{{}_{10}}}&=&A_{2}\varphi_{0}+B_{1}e^{-i\beta_{g_{{}_{10}}}}\varphi_{1}\\\
u_{h_{{}_{10}}}&=&A_{2}\varphi_{0}+B_{2}e^{-i(\theta_{2}+\beta_{h_{{}_{11}}})}\varphi_{1}\end{array}.\,\,\right.$
(5.26)

Note that the edges $g_{20},h_{01},h_{11}$ are outside the fundamental domain
$W_{\\!{}_{\hexagon}}$; in this case, we have the relations
$u_{g_{20}}(1)=e^{i\theta_{1}}u_{g}(1)$,
$u_{h_{01}}(1)=e^{i\theta_{2}}u_{h}(1)$ and
$u_{h_{11}}(1)=e^{i\theta_{2}}u_{h_{10}}(1)$. Using this fact, combining the
equations (5.25) and (5.26), with the same notations used in the previous
chapters for $\eta$ and $\alpha_{{}_{B}},\alpha_{{}_{N}}$, we obtain the
matrix representation

$M_{\hexagon}(\theta)\\!=\\!\left(\\!\\!\\!\begin{array}[]{cccc}-3\eta-\alpha_{{}_{N}}&e^{-i\beta_{f}}+e^{-i(\theta_{2}+\beta_{h_{01}})}&0&e^{-i(\theta_{1}+\beta_{g_{20}})}\\\
e^{i\beta_{f}}+e^{i(\theta_{2}+\beta_{h_{01}})}&-3\eta-\alpha_{{}_{B}}&e^{i\beta_{g_{10}}}&0\\\
0&e^{-i\beta_{g_{10}}}&-3\eta-\alpha_{{}_{N}}&e^{-i\beta_{f_{10}}}+e^{-i(\theta_{2}+\beta_{h_{11}})}\\\
e^{i(\theta_{1}+\beta_{g_{20})}}&0&e^{i\beta_{f_{10}}}+e^{i(\theta_{2}+\beta_{h_{11}})}&-3\eta-\alpha_{{}_{B}}\end{array}\\!\\!\\!\\!\right)$

and considering
$\beta_{f_{\gamma_{1},\gamma_{2}}}=\beta_{g_{\gamma_{1},\gamma_{2}}}=0$ and
$\beta_{h_{\gamma_{1},\gamma_{2}}}=-\gamma_{1}\phi$, with $\phi=\pi$, we
rewrite the matrix

$M_{\hexagon}\big{(}\theta_{1},\theta_{2}\big{)}=\left(\begin{array}[]{cccc}-3\eta-\alpha_{{}_{N}}&1+e^{-i\theta_{2}}&0&e^{-i\theta_{1}}\\\
1+e^{i\theta_{2}}&-3\eta-\alpha_{{}_{B}}&1&0\\\
0&1&-3\eta-\alpha_{{}_{N}}&1-e^{-i\theta_{2}}\\\
e^{i\theta_{1}}&0&1-e^{i\theta_{2}}&-3\eta-\alpha_{{}_{B}}\end{array}\right).$

Putting
$P(\eta):=\det\Big{(}M_{\hexagon}\big{(}\theta_{1},\theta_{2}\big{)}\Big{)}$,
we get

$\displaystyle P(\eta)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!81\eta^{4}+54(\alpha_{{}_{N}}+\alpha_{{}_{B}})\eta^{3}+9\big{[}\alpha_{{}_{N}}^{2}+\alpha_{{}_{B}}^{2}+4\alpha_{{}_{N}}\alpha_{{}_{B}}-6\big{]}\eta^{2}-6(\alpha_{{}_{N}}+\alpha_{{}_{B}})(3-\alpha_{{}_{N}}\alpha_{{}_{B}})\eta$
$\displaystyle+$
$\displaystyle\alpha_{{}_{B}}^{2}\alpha_{{}_{N}}^{2}-6\alpha_{{}_{N}}\alpha_{{}_{B}}+3+2\cos(\theta_{1}-2\theta_{2})-2\cos(\theta_{1})-2\cos(2\theta_{2})\,.$

Case $\alpha_{{}_{N}}=\alpha_{{}_{B}}$: Here, we get

$\displaystyle P_{B}(\eta)\\!\\!\\!$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\\!\\!\\!81\eta^{4}\\!+\\!108\alpha_{{}_{B}}\eta^{3}\\!+\\!54(\alpha_{{}_{B}}^{2}\\!-\\!1)\eta^{2}\\!+\\!12(\alpha_{{}_{B}}^{2}\\!-\\!3)\alpha_{{}_{B}}\eta\\!+\\!\alpha_{{}_{B}}^{4}\\!-\\!6\alpha_{{}_{B}}^{2}$
$\displaystyle\\!+\\!$
$\displaystyle\\!3\\!+\\!2\cos(\theta_{1}\\!-\\!2\theta_{2})\\!-\\!2\cos(\theta_{1})\\!-\\!2\cos(2\theta_{2}),$

whose roots are

$r^{\pm}_{\mp}(\theta_{1},\theta_{2})=\frac{\alpha_{{}_{B}}}{3}\pm\frac{\sqrt{3\mp\sqrt{2}\sqrt{3+\cos(\theta_{1})+\cos(2\theta_{2})-\cos(\theta_{1}-2\theta_{2})}\,\,}\,}{3}\,,$

which, according to its graphical representation in Figure 5.4 (with
$\alpha_{{}_{B}}=\alpha_{{}_{N}}=-1$), we see that the conic touch is
maintained for
$(\theta_{1},\theta_{2})\in[0,\frac{\pi}{2})\times[-\frac{\pi}{2},\frac{\pi}{2})$
for $\alpha_{{}_{B}}$ and $\alpha_{{}_{N}}$ equal; a fact that was already
expected, since $\alpha(v)$ would be a constant in this case.

Figure 5.4: For $q=2$, the cone is maintained at
$(\theta_{1},\theta_{2})\in[0,\frac{\pi}{2})\times[-\frac{\pi}{2},\frac{\pi}{2})$
for $\alpha_{{}_{N}}=\alpha_{{}_{B}}$, as expected, since in this case,
$\alpha(v)$ is a constant.

Case $\alpha_{{}_{N}}=-\alpha_{{}_{B}}$. In this case, we obtain a biquadratic
polynomial in $\eta$

$\displaystyle P_{B}(\eta)$ $\displaystyle\\!\\!=\\!\\!$ $\displaystyle
81\eta^{4}-6(3\alpha_{{}_{B}}^{2}+1)\eta^{2}+\alpha_{{}_{B}}^{4}+6\alpha_{{}_{B}}^{2}+3+2\cos(\theta_{1}-2\theta_{2})-2\cos(\theta_{1})-2\cos(2\theta_{2})$

whose roots are

$\eta^{\pm}_{\mp}(\theta_{1},\theta_{2})=\pm\frac{\sqrt{\alpha_{{}_{B}}^{2}+3\mp\sqrt{2}\sqrt{3+\cos(\theta_{1})+\cos(2\theta_{2})-\cos(\theta_{1}-2\theta_{2})}\,\,}\,}{3}\,,$

and from Figure 5.5 (with $\alpha_{{}_{B}}=1$ and $\alpha_{{}_{N}}=-1$), we
see that gaps (open) are formed at the points where conical touches occurred,
or that is, for non-constant $\alpha(v)$, the operator dispersion relation
does not have Dirac cones (even with rational magnetic flux with respect to
$2\pi$).

Figure 5.5: Graphical representation of the roots $\eta^{\pm}_{\mp}$ with
$\alpha_{{}_{B}}=1,\alpha_{{}_{N}}=-1$. For $q=2$, in
$(\theta_{1},\theta_{2})\in[0,\frac{\pi}{2})\times[-\frac{\pi}{2},\frac{\pi}{2})$
the cone is not maintained for $\alpha_{{}_{N}}=-\alpha_{{}_{B}}$, that is,
for non-constant $\alpha(v)$, the operator dispersion relation has an (open)
gap.

In general, for arbitrary values of real $\alpha_{{}_{B}}$ and
$\alpha_{{}_{N}}$, on each edge $e\in\hexagon$, we summarize these facts in
the following theorem.

######  Theorem 5.2.1 (Magnetic Graphene and hBN).

Let $\alpha_{{}_{B}},\alpha_{{}_{N}}$ be any real numbers,
$(\theta_{1},\theta_{2})\in[0,\frac{\pi}{2})\times[-\frac{\pi}{2},\frac{\pi}{2})$,
$a_{e}\in C^{1}[0,1]$ is a magnetic potential and $U\in L^{2}[0,1]$ is a
symmetric potential. Consider the operator

$H^{\phi}_{\hexagon}:(f_{e})\mapsto(-f^{\prime\prime}_{e}+Uf_{e}),\quad(f_{e})\in\bigoplus_{e\in\mathcal{E}}{\mathcal{H}}^{2}[0,1],$

with constant magnetic flux $\phi$ in the hexagons of the hexagonal lattice,
and boundary conditions

$\left\\{\begin{array}[]{ll}f_{a}(0)=e^{i\beta_{e}}f_{b}(1)=:f(v)\,,\forall\,\,a,b\in\mathcal{E}_{v}\,;\\\
\displaystyle\sum_{e;i(e)=v}f^{\prime}_{e}(0)-\sum_{e;t(e)=v}e^{\beta_{e}}f^{\prime}_{e}(1)=\alpha(v)f(v)\,,\,\,v\in\mathcal{V}\,,\end{array}\,\right.\,$

where $\beta_{e}:=\int^{1}_{0}\\!\\!a_{e}(s)ds$. Suppose $\phi=\pi$; then the
dispersion relation of the operator $H^{\phi}_{\hexagon}$ is given by

$\eta^{\pm}_{\mp}(\theta_{1},\theta_{2})=-\frac{\alpha_{{}_{N}}+\alpha_{{}_{B}}}{6}\pm\frac{\sqrt{(\alpha_{{}_{N}}-\alpha_{{}_{B}})^{2}+12\mp
4\sqrt{2}\sqrt{3+\cos(\theta_{1})+\cos(2\theta_{2})-\cos(\theta_{1}-2\theta_{2})\,\,}\,}}{6}\,.$

From this, it can be concluded that

  * (a)

If $\alpha(v)$ is constant at the vertices $v\in\hexagon$,
$H^{\phi}_{\hexagon}$ has Dirac cones;

  * (b)

If $\alpha(v)$ is not constant at the vertices $v\in\hexagon$, then
$H^{\phi}_{\hexagon}$ has gaps.

In particular, for $\phi=\pi$ or $\phi=2\pi$, magnetic graphene has a Dirac
cone, while magnetic hexagonal boron nitride has a gap between the bands.

Note that taking
$(\theta_{1},\theta_{2})\in[0,\frac{\pi}{q})\times[-\frac{\pi}{q},\frac{\pi}{q})$,
we saw that for $q=1$, according to the Chapter 1 and also in this section,
the Dirac Cones, if they exist, occur when $\theta_{1}=-\theta_{2}$, while for
$q=2$ conic touches, when they exist, occur for $\theta_{1}=-2\theta_{2}$. In
fact, putting

$F(\theta_{1},\theta_{2}):=3+\cos(\theta_{1})+\cos(2\theta_{2})-\cos(\theta_{1}-2\theta_{2}),$

we see that $F(\theta_{1},\theta_{2})=0$ if, and only if,
$\theta_{1}=\pi+2m\pi$ and $\theta_{2}=-\frac{\pi}{2}+n\pi$ ,
$m,n\in\mathbb{Z}$, that is, $F$ cancels only under the condition
$\theta_{1}=-2\theta_{2}$.

###  5.3 Analysis and conclusions

According to the development explained above, of these models with magnetic
field, if $\alpha(v)$ is not constant, it is not possible to guarantee the
relation (5.13), that is, it is not guaranteed that it is possible to relate
the spectra from $\Delta_{\beta}$ into $l^{2}(\hexagon)$ and $H$ into
$\bigoplus_{e\in\mathcal{E}}{\mathcal{H}}^{2}[0,1]$ , since such a relation
was obtained by supposing $\alpha(v)=\frac{\alpha}{2}\operatorname{deg}(v)$,
with $\alpha\in\mathbb{R}$ taking the same value in all the vertices
(according to our proposal, they would all be models with a single type of
atom at all vertices, like graphene).

We then show that there are integer values $p$ and $q$, and for certain values
of $\alpha$ taking two values on the hexagonal lattice $\hexagon$ (thus a
model for boron nitride), in which the conical touches on the dispersion
relation of the operator $H^{\phi}_{\hexagon}$ cease to exist and we also
guarantee the existence of gaps. To do so, we use the Theorem 5.1.1 and the
relation (5.13), however we are mainly based on the technique presented in the
Chapter 1, which we adapt to the case with a magnetic field in a specific
configuration in which we used a particular period of the operator, which
differed from the period of the hexagonal lattice. It was convenient to use
this technique due to the fact that we do not have an analogue to the Chambers
Formula used in [2], which also restricts us to small periods.

We then conclude that the hexagonal quantum graph model, with the continuous
Laplacian operator on the edges and with a magnetic field, which models the
magnetic hexagonal boron nitride, does not always have a Dirac cone if the
modeling considers the flow $\phi=2\pi\frac{p}{q}$ (rational multiple of
$2\pi$) constant in this network. In the case of boron nitride, we take for
the parameters in the boundary conditions the values
$\alpha(v)=\delta_{{}_{N}}$ if $v$ is the initial vertex of the edges, and
$\alpha(v)=\delta_{{}_{B}}$ if $v$ for final vertex. On the other hand, this
freedom of choosing different values for $\alpha(v)$ at the vertices does not
occur in the case of graphene, since this material has the same constant
$\alpha(v)=\delta_{{}_{C}}$ at all vertices of the hexagonal lattice and,
consequently, the operator proposed to model graphene will always present
Dirac cones, for any values of $p$ and $q$ such that the magnetic flux $\phi$
in each hexagon is a rational multiple of $2\pi$. These adaptations to
$\delta_{{}_{C}}\neq 0$ and the example with rational flow in relation to
$2\pi$ without Dirac cones constitute our main contributions to this model and
also to this last chapter of the thesis.

##  Chapter
6
Final considerations

In Chapter 1, we presented the modeling of a sheet of hexagonal material with
two types of alternating atoms at the vertices and concluded that its
dispersion relationship will have a Dirac cone if the atoms are equal and will
have gaps, otherwise. In particular, graphene monolayer always has a conical
touch, while hexagonal boron nitride always has gaps.

Then, in Chapter 2, we extend the study in Chapter 1 to stacked bilayers of
the $AA$ type; we conclude, according to our modeling, that in materials of
this nature there will be a Dirac cone only if the atoms are equal, and there
will be parabolic touches if the parameters comply with the relation
$\alpha_{b}=\pm\,2t_{0}^{2}\,+\,\alpha_{a}$; in particular, two sheets of hBN,
stacked like $AA$, do not have Dirac cones, but may have parabolic touches.

On the other hand, if in the configuration of Chapter 2 the stacking in type
$AA^{\prime}$ is considered, we conclude in Chapter 3 that this material will
not have conical or parabolic touches if the atoms are distinct, but they have
parabolic touch if the atoms are equal. In particular, in the $AA^{\prime}$
stacking, hBN bilayers always have holes, while graphene bilayers always have
Dirac cones.

Heterostructures were considered in Chapter 4; we confirmed the presence of
Dirac cones in the system composed of a layer of graphene and another of hBN,
a fact already covered in the physics literature. We also identified conical
touches or gaps in the graphene between boron nitride and boron nitride
between graphene sandwiches, respectively. In the former, graphene caused a
reduction in the openings of existing gaps in the hBN bilayers, while in this
case, the hBN induced a gap in some Dirac points of the graphene and a conical
touch is always maintained in this configuration.

Finally, in Chapter 5, we use quantum graphs in a proposal for modeling
hexagonal boron nitride under the effect of a magnetic field, through the
Schrödinger operator, and demonstrate that if the magnetic flux $\phi$ is
constant in the lattice hexagons and of the form $\phi=2\pi\dfrac{p}{q}$,
there are integer values $p$ and $q$ such that, for certain boundary
conditions $\alpha(v)$ at the vertices $v$ of the network, gaps are formed at
the conic points of the dispersion relation of the operator defined in this
network.

Although we present a modest and simplified mathematical modeling proposal for
hBN and graphene, we consider the results obtained analytically to be
consistent and satisfactory: we confirm some facts already known in the
physics literature obtained, in general, numerically or through laboratory
experiments, and in some In some cases (such as $AA$ stacking) we have
behaviors obtained from this type of model that would be interesting to
confirm.

## References

  * [1] Amovilli, C., Leys, F. E., March, N. H. Electronic energy spectrum of two-dimensional solids and a chain of C atoms from a quantum network model. J. Math. Chem., 36, 93–112 (2004).
  * [2] Becker, S., Han, R. Density of states and delocalization for discrete magnetic random Schrödinger operators. International Mathematics Research Notices, 5, 1–58 (2021).
  * [3] Becker, S.; Han, R.; Jitomirskaya, S. Cantor spectrum of graphene in magnetic fields. Invent. Math., 218, 979–1041 (2019).
  * [4] Berkolaiko, G., Comech, A. Symmetry and Dirac points in graphene spectrum. J. Spect. Theor., 8, 1099–1147 (2018).
  * [5] Berkolaiko, G., Kuchment, P. Introduction to quantum graphs. American Mathematical Society (2012).
  * [6] Brown, M. B., Eastham, M. S. P., Schmidt, K. M. Periodic differential operators. Birkhäuser Basel (2013).
  * [7] Castro Neto, A. H., Guinea, F., Peres, N. M. R., Novoselov, K. S., Geim, A. K. The electronic properties of graphene, Rev. Mod. Phys., 81, 109–162 (2014).
  * [8] de Oliveira, C. R., Rocha, V. L. Dirac cones for bi-and trilayer Bernal-stacked graphene in a quantum graph model. J. Phys. A: Math. Theor., 53, 505201 (2020).
  * [9] de Oliveira, C. R., Rocha, V. L. From multilayer $AA$-stacked graphene sheets to graphite: graph models and Dirac cones. Zeitschrift für Naturforschung A, 76, 371–384 (2021).
  * [10] de Oliveira, C. R., Rocha, V. L. Bidimensional honeycomb materials: a graph model through Dirac operator. Rep. Math. Phys., 89, 231–252 (2022).
  * [11] Do, N. T., Kuchment, P. Quantum graph spectra of a graphyne structure. Nanoscale systems MMTA, 2, 107–123 (2013).
  * [12] Eastham, M. S. P. The spectral theory of periodic differential equations. Scottish Acad. Press Ltd., Edinburgh-London (1973).
  * [13] Fang, W. Synthesis of bilayer graphene and hexagonal boron nitride by chemical vapor deposition method. Massachusetts Institute of Technology (2016).
  * [14] Fefferman, C. L., Weinstein, M. I. Honeycomb lattice potentials and Dirac points. J. Amer. Math. Soc., 25, 1169–1220 (2012).
  * [15] Fefferman, C. L., Weinstein, M. I. Wave packets in honeycomb structures and two-dimensional Dirac equations. Commun. Math. Phys., 326, 251–286 (2014).
  * [16] Feijó, T. O. Crescimento de heteroestruturas de van der Waals visando a aplicações em nanoeletrônica. Tese de Doutorado, Programa de Pós-Graduação em Microeletrônica, Universidade Federal do Rio Grande do Sul (2021).
  * [17] Fisher, L., Li, W., Shipman, S. P. Reducible Fermi Surface for Multi-layer Quantum Graphs Including Stacked Graphene. Commun. Math. Phys. 385, 1499–1534 (2021).
  * [18] Geim, A. K. Van der Waals heterostructures. Nature, 499, 419–425 (2013).
  * [19] Giovannetti, G., Khomyakov, P. A., Brocks, G., van den Brink, J., Kelly, P. J. Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations. Phys Review B, 76, 073103 (2007).
  * [20] Helffer, B., Kerdelhué, P., and Royo-Letelier, J. Chambers’s formula for the graphene and the Hou model with Kagome periodicity and applications. Annales Henri Poincaré, 17, 795–818 (2016).
  * [21] J. Brüning, V. Geyler, K. Pankrashkin. Cantor and band spectra for periodic quantum graphs with magnetic fields. Comm. Math. Phys., 269, 87–105 (2007).
  * [22] Katsnelson, M. I., The Physics of Graphene, 2nd Ed., Cambridge Univ. Press, Cambridge, 2012.
  * [23] Kopeliovich, D. Thermal conductivity of Boron nitride, SubsTech substances and technologies articles. USA (2011).
  * [24] Kuchment, P. Floquet theory for partial differential equations. Basel, Birkhäuser (1993).
  * [25] Kuchment, P. Quantum graphs: I. Some basic structures. Waves Random Media, 14, S107-S128 (2004).
  * [26] Kuchment, P., Post, O. On the spectra of carbon nano-structures. Commun. Math. Phys., 275, 805–826 (2006).
  * [27] Magnus, W., Winkler, S. Hill’s Equation. Wiley, New York (1966).
  * [28] Matos, M. J de Souza. Estudo teórico das propriedades de nitreto de boro hexagonal e grafeno: resposta dielétrica, empilhamento e estrutura eletrônica. Tese de doutorado. Departamento de Física, UFMG (2014).
  * [29] Pankrashkin, K. Spectra of Schrödinger operators on equilateral quantum graphs. Letters in Mathematical Physics, 77, 139–154 (2006).
  * [30] Pankrashkin, K. Unitary dimension reduction for a class of self-adjoint extensions with applications to graph-like structures. Elsevier, Paris (2012).
  * [31] Reed, M., Simon, B. Methods of Modern Mathematical Physics IV: Analysis of Operators. New York, Academic Press (1978).
  * [32] Rocha, Vinícius L. Multilayer graphene through quantum periodic graphs: Dirac cones. Tese de Doutorado, Departamento de Matemática, Universidade Federal de São Carlos (2020).
  * [33] Shipman, S. P. Reducible Fermi surfaces for non-symmetric bilayer quantum-graph operators. J. Spectr. Theory, 10, 33–72 (2020).
  * [34] Souza, E. dos Santos. Estudo teórico de bicamada de grafeno em superfície CU (111). Dissertação de mestrado, UFES (2014).
  * [35] Tarnopolsky, G., Kruchkov, A. J., Vishwanath, A. Origin of magic angles in twisted bilayer graphene. Phys. Rev. Lett., 122, 106405 (2019).
  * [36] Wallace, P. R., The band theory of graphite. Phys. Rev., 71, 622–634 (1947).
  * [37] Wang, E., Lu, X., Ding, S., Yao, W., Yan, M., Wan, G., Deng, K., Wang, S., Chen, G., Ma, L., Jung, J., Fedorov, A. V., Zhang, Y., Zhang, G., Zhou, S. Gaps induced by inversion symmetry breaking and second-generation Dirac cones in graphene/hexagonal boron nitride. Nature Physics, 12, 1111-1116 (2016).


# Similarity solutions for two-phase fluids models

Andronikos Paliathanasis
Institute of Systems Science, Durban University of Technology
PO Box 1334, Durban 4000, Republic of South Africa
Instituto de Ciencias Físicas y Matemáticas,
Universidad Austral de Chile, Valdivia, Chile Email<EMAIL_ADDRESS>
###### Abstract

The algebraic properties of drift-flux two-phase fluids models without
gravitational and wall friction forces are studied. More precisely, for the
two fluids we consider equation of states of polytropic gases. We perform a
classification scheme of the unknown parameters of the model such that to
determine all the possible admitted Lie symmetries. We find that in the most
general case the dynamical system of hyperbolic equations is invariant under
the action of a four dimensional Lie algebra, while the larger number of
admitted Lie symmetries is six. For each admitted Lie algebra the one-
dimensional optimal system is derived which is applied for the determination
of all the unique similarity transformations which lead to similarity
solutions. Our results are compared with that of previous studies from where
we see that most of the solutions presented in this study have not found
before in the literature.

Keywords: Lie symmetries; invariants; two-phase flows; similarity solutions;
hyperbolic equations

##  1 Introduction

A systematic approach for the determination of exact and analytical solutions
of nonlinear differential equations is Lie’s theory [1, 2, 3]. The novelty of
Lie’s consideration is that someone investigate the invariance properties of a
given differential equation under the action of infinitesimal transformations
in order to write the differential equation into an equivalent form through
algebraic representation [4, 5, 6, 7]. The generator of the infinitesimal
transformation which keeps a differential equation invariant is called Lie
symmetry of the right type. A main property for the admitted Lie symmetries
for a differential equation is that they form an algebra known as Lie algebra.

The existence of a large-dimensional Lie algebra of Lie symmetries for a
differential equation enables one to solve the differential equation by means
of repeated reduction of order with the use of similarity transformations, or
by means of the determination of a sufficient number of first integrals [8,
9]. In particular if an $n$th-order differential equation admit at least $n$
Lie point symmetries with a solvable algebra, one knows that the system is
reducible to quadratures and the system is integrable [10]. Thus, the absence
of the latter property does not immediately obviate the possibility of
integrability. For instance, the application of a similarity transformation
reduce the differential equation into a new differential equation which may
has different algebraic properties and additional (new) symmetries to follow,
[11]. On the other hand, there is a zoology on the nature of symmetries which
means that because an equation does not admit Lie point symmetries it does not
means that it is not invariant under the action of other kind transformations
which follows from generalized forms of the definition of symmetries for
instance see [12, 13, 14] and references therein. Another important
application of the symmetries is that they are applied for the classification
of equations and establish classes and families of equations which has as
common feature the common admitted Lie algebra [15, 16, 17, 18, 19, 20].

There are various applications in the literature on the application of Lie
symmetries for the study of nonlinear differential equations. Some results on
ordinary differential equations are presented in [21, 22], while there are
various studies and for elliptic partial differential equations, for instance
see [23, 24, 25] and references therein. As far as the application of Lie
symmetries on hyperbolic partial differential equations is concerned, there
are various of studies in all areas of applied mathematics [26, 27, 28, 29].
In theory of fluid dynamics Lie symmetries have played an important role on
the determination of exact and analytic solutions. The complete symmetry
classification of the shallow-water equations with or without a gravitational
field performed in [30, 31], while the case with a Coriolis force in [32, 33].
The case with various class of bottoms investigated in [34, 35, 36] while for
other studies on shallow-water equations we refer the reader in [37, 38, 39]
and references therein. A detailed discussion of the symmetry approach on
mechanical theories of continuous media is presented in [40]. The authors by
studying some systems of physical interests demonstrate the application of Lie
invariants for the derivation of similarity solutions.

Another important system of hyperbolic equations where Lie point symmetries
have been used for the determination of new solutions is the two-phase flow
model [41, 42]. The two-phase fluids model describes the evolution of two
fluids with different phases, such as liquid and gas, in a tube. The two-phase
models have many physical applications from oil extraction, underground water,
nuclear reactor and many others [43]. For some important results on two-phase
flow models we refer the reader in [44, 45, 46, 47, 48, 49, 50] and references
therein. The exact solution of the Riemann problem for the drift-flux
equations for a two-phase flow system was studied in [51] while numerical
simulations of wave propagation in compressible two-phase flow were presented
in [52].

In [41, 42] Lie’s theory applied for the simplest two-phase flow model where
there is not any mass transfer from the one fluid to the other while the
pressure and the energy density of the two fluids is given by an polytropic
equation of state as it is given by Lane-Emden equation. In this work we
revise the results of [41, 42] for a non-flip drift flux model of multi-phase
flow defined [53]. More specifically, we find new solutions which have not
been presented before, for linear and nonlinear equation of state parameters
for the two fluids. More precisely, we derive the admitted Lie point
symmetries of the two-phase fluids model and for the admitted Lie symmetries
we calculate the commutators and the Adjoint representation. By using these
results we are able to determine the one-dimensional optimal system such that
to perform all the independent similarity transformations. We find that the
similarity solutions are expressed either by closed-form expressions or by
quadratures. There are not many known exact solutions in two-phase models, see
for instance [54, 55, 56], hence the analysis presented here is important for
the determination of exact solutions. The plan of the paper is as follows.

The main mathematical properties and definitions on the symmetries of
differential equations are presented in Section 2. The main results of this
study are presented in Section 3. The drift-flux two-phase fluids model is
given where we determine the admitted Lie point symmetries. The unknown
parameters of the problem are the two polytropic exponents
$\gamma_{1},~{}\gamma_{2}$ of the two fluids. For arbitrary values of the
polytropic exponents such that $\gamma_{i}\geq 1,$ the set of three hyperbolic
equations admit four Lie point symmetries, while two additional Lie symmetries
exist when $\gamma_{1}=\gamma_{2}=\gamma$ which means that in the latter case
the dynamical system is invariant under a sixth dimensional Lie algebra. For
each different admitted Lie algebra we calculate the commutators and the
Adjoint representation. By using these results we are able to derive the one-
dimensional optimal system, necessary for the derivation of all the
independent similarity solutions. On the derivation of the similarity
solutions emphasis is given in the two cases in which
$\gamma_{1}=\gamma_{2}=2$ and $\gamma_{1}=\gamma_{2}=1$. Finally, in 4 we
discuss our results and we draw our conclusions.

##  2 Lie symmetries of differential equations

In this Section we present the basic properties and definitions of the Lie
point symmetries. Consider the system of differential equations
$H^{A}\left(x^{i},\Phi^{A},\Phi_{,i}^{A},...\right)=0$ where $x^{i}$ denotes
the indepedent variables and $\Phi^{A}$ are the independent variables while a
$\Phi_{,i}^{A}$ indicates derivative with respect to the variable $x^{i}$,
that is, $\Phi_{,i}^{A}=\frac{\partial}{\partial x^{i}}\Phi^{A}$.

We proceed by assuming the infinitesimal one-parameter point transformation of
the form

$\bar{x}^{i}=x^{i}\left(x^{j},\Phi^{B};\varepsilon\right),~{}~{}\bar{\Phi}^{A}=\Phi^{A}\left(x^{j},\Phi^{B};\varepsilon\right),$
(1)

where $\varepsilon~{}$is the infinitesimal parameter. Point transformation (1)
connects two different points $P\left(x^{j},\Phi^{B}\right)\rightarrow
Q\left(\bar{x}^{j},\bar{\Phi}^{B},\varepsilon\right)$, while at these two
points the system of differential equations is defined $H^{A}\left(P\right)$
and $\bar{H}^{A}\left(Q\right)$.

We shall say that the system $H^{A}$ remains invariant under the action of the
one-parameter transformation if and only if $\bar{H}^{A}=H^{A}$, that is [4,
5, 6, 7]

$\lim_{\varepsilon\rightarrow
0}\frac{\bar{H}^{A}\left(\bar{y}^{i},\bar{u}^{A},...;\varepsilon\right)-H^{A}\left(y^{i},u^{A},...\right)}{\varepsilon}=0.$
(2)

Which means that the solutions $\Phi^{A}~{}$of the system $H^{A}$ at the two
different points, i.e. $\Phi^{A}\left(P\right)$ and $\Phi^{A}\left(Q\right),$
are related through the point transformation (1).

By definition, expression (2) is the Lie derivative of $H^{A}$ along the
vector field $X$ of the one-parameter point transformation (1), in which $X$
is defined as

$X=\frac{\partial\bar{x}^{i}}{\partial\varepsilon}\partial_{i}+\frac{\partial\bar{\Phi}}{\partial\varepsilon}\partial_{A}.$

Hence, an equivalent form of the symmetry condition is

$\mathcal{L}_{X}\left(H^{A}\right)=0,$ (3)

where $\mathcal{L}$ denotes the Lie derivative with respect to the vector
field $X^{\left[n\right]}.$ Vector field $X^{\left[n\right]}$ the $n$th-
extension of generator $X~{}$of the transformation (1) in the jet space
$\left\\{x^{i},\Phi^{A},\Phi_{,i}^{A},\Phi_{,ij}^{A},...\right\\}$ where
$H^{A}$ is defined. The vector field $X^{\left[n\right]}$ is calculated by the
generator of the point transformation $X$ as follows

$X^{\left[n\right]}=X+\eta^{\left[1\right]}\partial_{\Phi_{i}^{A}}+...+\eta^{\left[n\right]}\partial_{\Phi_{i_{i}i_{j}...i_{n}}^{A}},$
(4)

where the new terms are

$\eta^{\left[n\right]}=D_{i}\eta^{\left[n-1\right]}-u_{i_{1}i_{2}...i_{n-1}}D_{i}\left(\frac{\partial\bar{x}^{j}}{\partial\varepsilon}\right)~{},~{}i\succeq
1~{},~{}\eta^{\left[0\right]}=\left(\frac{\partial\bar{\Phi}^{A}}{\partial\varepsilon}\right).$
(5)

The Lie symmetries for a given differential equation form a Lie algebra. Lie
symmetries can be used by different ways [5] in order to study a differential
equation. However, their direct application is on the determination of the so-
called similarity solutions. The steps which we follow to determine a
similarity solution is based on the determination and application of the Lie
invariant functions.

In order to determine all the possible independent solutions of a given
dynamical system we should derive the one-dimensional optimal system. Let the
$n$-dimensional Lie algebra $G_{n}$ with elements
$X_{1},~{}X_{2},~{}...~{}X_{n}$ admitted by the system $H^{A}$. The two vector
fields [7]

$Z=\sum\limits_{i=1}^{n}a_{i}X_{i}~{},~{}W=\sum\limits_{i=1}^{n}b_{i}X_{i}~{},~{}\text{\
}a_{i},~{}b_{i}\text{ are constants.}$ (6)

are equivalent if and only if

$\mathbf{W}=Ad\left(\exp\left(\varepsilon_{i}X_{i}\right)\right)\mathbf{Z}$
(7)

or

$W=cZ~{},~{}c=const.$ (8)

where the operator $Ad\left(\exp\left(\varepsilon
X_{i}\right)\right)X_{j}=X_{j}-\varepsilon\left[X_{i},X_{j}\right]+\frac{1}{2}\varepsilon^{2}\left[X_{i},\left[X_{i},X_{j}\right]\right]+...~{}$is
called the adjoint representation.

##  3 Two-phase flow model

Let us assume a mixture fluid stratified flow in a pipe, where the mixture is
consisting by two phases of the same fluid, i.e. liquid
$\left(\rho_{l}\right)$ and gas $\left(\rho_{g}\right)$. Furthermore, each
phase moves on the pipe at a local section with average velocity $u_{l}$ and
$u_{g}$ . Hence, the continuous equations for the two fluids in the one-
dimensional read [57]

$\displaystyle\left(\rho_{l}a_{l}\right)_{,t}+\left(\rho_{l}a_{l}u_{l}\right)_{,x}$
$\displaystyle=$ $\displaystyle\Gamma_{l}a_{l},$ (9)
$\displaystyle\left(\rho_{g}a_{g}\right)_{,t}+\left(\rho_{g}a_{g}u_{g}\right)_{,x}$
$\displaystyle=$ $\displaystyle\Gamma_{g}a_{g},$ (10)

where $\Gamma_{l},~{}\Gamma_{g}$ are the particle creation terms, where for a
closed system hold $\Gamma_{l}=-\Gamma_{g}$. Furthermore,
$a_{l},~{}a_{g}~{}$are the volume fractions for the two fluids with
$a_{l}+a_{g}=1$. We continue our analysis by assume the simplest model without
particle creation term, or any interaction between the two fluids, while we
omit any gravitational effect.

Hence, the momentum equations for the fluids reads [57]

$\left(\rho_{l}a_{l}u_{l}\right)_{,t}+\left(\rho_{l}a_{l}u_{l}^{2}+a_{l}p_{l}\right)_{,x}-p_{l}^{i}\left(a_{l}\right)_{,x}=0,$
(11)
$\left(\rho_{g}a_{g}u_{g}\right)_{,t}+\left(\rho_{g}a_{g}u_{g}^{2}+a_{g}p_{g}\right)_{,x}-p_{g}^{i}\left(a_{g}\right)_{,x}=0.$
(12)

in which $p_{l},~{}p_{g}$ are the pressure terms for the two fluids and
$p_{l}^{i},~{}p_{g}^{i}$ describe the interfacial pressures on the gas–liquid
interface on the side of the liquid and of the gas respectively.

Hence, from the momentum equations (11), (12) it follows

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(\rho_{l}a_{l}u_{l}+\rho_{g}a_{g}u_{g}\right)_{,t}+\left(\rho_{l}a_{l}u_{l}^{2}+\rho_{g}a_{g}u_{g}^{2}\right)_{,x}+$
(13)
$\displaystyle+a_{l}p_{l,x}+a_{g}p_{g,x}+\left(p_{l}-p_{l}^{i}\right)\left(a_{l}\right)_{,x}+\left(p_{g}-p_{g}^{i}\right)\left(a_{g}\right)_{,x}$

Therefore, with the use of the new variables $u_{r}=u_{g}-u_{l}$ and
$u=\frac{\rho_{g}a_{g}u_{g}+\rho_{l}a_{l}u_{l}}{\rho_{g}a_{g}+\rho_{l}a_{l}}+u_{r},$
and by considering that for the drift model that $u_{r}=0$, we end with the
finally system [42]

$\displaystyle\rho_{1,t}+\left(\rho_{1}u\right)_{,x}$ $\displaystyle=$
$\displaystyle 0,$ (14) $\displaystyle\rho_{2,t}+\left(\rho_{2}u\right)_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (15)
$\displaystyle\left(\left(\rho_{1}+\rho_{2}\right)u\right)_{,t}+\left(\left(\rho_{1}+\rho_{2}\right)u^{2}+p_{1}+p_{2}\right)_{,x}$
$\displaystyle=$ $\displaystyle 0.$ (16)

where
$\rho_{\left(1,2\right)}\left(t,x\right)=a_{\left(l,g\right)}\rho_{\left(l,g\right)}\left(t,x\right)$
are the energy density of the fluid of phase and
$p_{\left(1,2\right)}\left(t,x\right)=a_{\left(l,g\right)}p_{\left(l,g\right)}\left(t,x\right)$
denote the pressure of the fluids. Recall that for the drift flux model, by
definition, the pressure in the gas and the liquid surfaces are equal for the
same cross-sectional area [57].

Finally, for the fluids we assume polytropic equation of state parameters of
the form
$p_{\left(1,2\right)}=\kappa_{\left(1,2\right)}\rho^{\gamma_{\left(1,2\right)}}$
with $\gamma_{\left(1,2\right)}\geq 1$.

###  3.1 Arbitrary polytropic exponents

For the system (14)-(16) we apply Lie’s theory and we determine the Lie point
symmetries which are [42]

$X_{1}=\partial_{t}~{},~{}X_{2}=\partial_{x}~{},~{}X_{3}=t\partial_{t}+x\partial_{x}\text{
and }X_{4}=t\partial_{x}+\partial_{u}$

for arbitrary values of the polytropic exponents $\gamma_{1},~{}\gamma_{2}$.
The commutators of the admitted Lie symmetries by the system (14)-(16) and the
Adjoint representation of the Lie algebra consisted by the elements
$\left\\{X_{1},X_{2},X_{3},X_{4}\right\\}$are presented in Tables 1 and 2. The
admitted Lie algebra is a solvable four-dimensional Lie algebra classified as
$A_{4,6}$ in the Patera and Winternitz classification scheme [18].

From Tables 1 and 2 we infer that the one-dimensional optimal system is
consisted by the elements

$\displaystyle\left\\{X_{1}\right\\}~{},~{}\left\\{X_{2}\right\\}~{},~{}\left\\{X_{3}\right\\}~{},~{}\left\\{X_{4}\right\\}~{},~{}$
$\displaystyle\left\\{X_{1}+\beta X_{2}\right\\}~{},~{}\left\\{X_{3}+\beta
X_{4}\right\\}~{},\left\\{X_{1}+\beta X_{4}\right\\}.$

Table 1: Commutator table for the Lie point symmetries of the two-phase flow system (14)-(16) $\left[X_{I},X_{J}\right]$ | $\mathbf{X}_{1}$ | $\mathbf{X}_{2}$ | $\mathbf{X}_{3}$ | $\mathbf{X}_{4}$
---|---|---|---|---
$\mathbf{X}_{1}$ | $0$ | $0$ | $X_{1}$ | $X_{2}$
$\mathbf{X}_{2}$ | $0$ | $0$ | $X_{2}$ | $0$
$\mathbf{X}_{3}$ | $-X_{1}$ | $-X_{2}$ | $0$ |  $0$
$\mathbf{X}_{4}$ | $-X_{2}$ | $0$ | $0$ | $0$
Table 2: Adjoint representation for the Lie point symmetries of the two-phase flow system (14)-(16) $Ad\left(e^{\left(\varepsilon\mathbf{X}_{i}\right)}\right)\mathbf{X}_{j}$ | $\mathbf{X}_{1}$ | $\mathbf{X}_{2}$ | $\mathbf{X}_{3}$ | $\mathbf{X}_{4}$
---|---|---|---|---
$\mathbf{X}_{1}$ | $X_{1}$ | $X_{2}$ | $-\varepsilon X_{1}+X_{3}$ | $-\varepsilon X_{2}+X_{4}$
$\mathbf{X}_{2}$ | $X_{1}$ | $X_{2}$ | $-\varepsilon X_{2}+X_{3}$ | $X_{4}$
$\mathbf{X}_{3}$ | $e^{\varepsilon}X_{1}$ | $e^{\varepsilon}X_{2}$ | $X_{3}$ | $X_{4}$
$\mathbf{X}_{4}$ | $X_{1}+\varepsilon X_{2}$ | $X_{2}$ | $X_{3}$ | $X_{4}$

We proceed with the application of the Lie symmetries such that to find the
invariant transformations where the solution of the system (14)-(16) is
expressed by the solution of an ordinary differential equation.

####  3.1.1 Reduction with $X_{1}$

Application of the Lie point symmetry $X_{1}$ gives stationary solutions, that
is, $\rho_{i}=\rho_{i}\left(x\right),~{}u=u\left(x\right)$ where the reduced
system is

$\left(\rho_{1}u\right)_{,x}=0~{},~{}\left(\rho_{2}u\right)_{,x}=0~{},\left(\left(\rho_{1}+\rho_{2}\right)u^{2}+\kappa_{1}\rho_{1}^{\gamma_{1}}+\kappa_{2}\rho_{2}^{\gamma_{2}}\right)_{,x}=0,$
(17)

which means $\rho_{1}u=\rho_{10},~{}\rho_{2}u=\rho_{20}$ and

$\left(\rho_{10}+\rho_{20}\right)u+\kappa_{1}\left(\rho_{10}\right)^{\gamma_{1}}u^{-\gamma_{1}}+\kappa_{2}\left(\rho_{20}\right)^{\gamma_{2}}u^{-\gamma_{2}}=u_{10}.$
(18)

Thus, it follows, $u\left(x\right)=const$.

####  3.1.2 Reduction with $X_{2}$

From the Lie point symmetry $X_{2}$ the stationary similarity transformation
are $\rho_{i}=\rho_{i}\left(t\right),~{}u=u\left(t\right)$, where

$\rho_{1,t}=0~{},~{}\rho_{2,t}=0~{},~{}\left(\left(\rho_{1}+\rho_{2}\right)u\right)_{,t}=0,$
(19)

that is, $\rho_{i}=\rho_{i0}$ and $u=u_{0}$.

####  3.1.3 Reduction with $X_{3}$

From the scaling symmetry $X_{3}$ it follows the similarity transformation
$\rho_{i}=\rho_{i}\left(w\right),~{}u=u\left(w\right)$ where the new
independent variable is defined as $w=\frac{x}{t}$. The reduced system is

$\displaystyle w\rho_{1,w}-\left(\rho_{1}u\right)_{,w}$ $\displaystyle=$
$\displaystyle 0,$ (20) $\displaystyle
w\rho_{2,w}-\left(\rho_{2}u\right)_{,w}$ $\displaystyle=$ $\displaystyle 0,$
(21) $\displaystyle
w\left(\left(\rho_{1}+\rho_{2}\right)u\right)_{,w}+\left(\left(\rho_{1}+\rho_{2}\right)u^{2}+p_{1}+p_{2}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0.$ (22)

From (20) and (21) it follows

$\rho_{i}\left(w\right)=\rho_{i0}\exp\left(-\int\frac{u_{,w}}{u-w}dw\right),$
(23)

where by replacing in (22) we end with the equation

$U_{,ww}\left(\left(U_{,w}\right)^{-1}\left(\rho_{01}+\rho_{02}\right)\left(U_{0}-U\right)^{2}-\kappa_{1}\left(\rho_{01}\right)^{\gamma_{1}}\left(U_{,w}\right)^{\gamma_{1}}-\kappa_{1}\left(\rho_{02}\right)^{\gamma_{1}}\left(U_{,w}\right)^{\gamma_{2}}\right)=0,$
(24)

in which we have defined
$u\left(w\right)=w-\left(U_{0}-U\left(w\right)\right)\left(U_{,w}\right)^{-1}$.

Hence, function $U\left(w\right)$ is given by the non-static solution of the
first-order differential equation

$\left(U_{,w}\right)^{-1}\left(\rho_{01}+\rho_{02}\right)\left(U_{0}-U\right)^{2}-\kappa_{1}\left(\rho_{01}\right)^{\gamma_{1}}\left(U_{,w}\right)^{\gamma_{1}}-\kappa_{1}\left(\rho_{02}\right)^{\gamma_{1}}\left(U_{,w}\right)^{\gamma_{2}}=0.$
(25)

The latter equation is solved by quadratures, however there are some closed-
form solutions for specific values of the polytropic exponents $\gamma_{i}$.
Indeed, for $\gamma_{i}=1$ the closed-form solution is

$U_{\pm}\left(w\right)=U_{0}-U_{1}e^{\pm\lambda w}$ (26)

where now
$\lambda=\sqrt{\frac{\rho_{01}+\rho_{02}}{\kappa_{1}\rho_{0}+\kappa_{2}\rho_{02}}}$,
thus the similarity solution is

$u\left(w\right)=w-\frac{1}{\lambda}~{},~{}\rho_{i}\left(w\right)=\rho_{i0}e^{\lambda
w}.$ (27)

This is a new close-form solution which has not found before.

####  3.1.4 Reduction with $X_{4}$

The similarity transformation given by the symmetry vector $X_{4}$ is
$\rho_{i}=\rho_{i}\left(t\right)$, $u=\frac{x}{t}+v\left(t\right)$ where

$\displaystyle t\rho_{1,t}+\rho_{1}$ $\displaystyle=$ $\displaystyle 0,$ (28)
$\displaystyle t\rho_{2,t}+\rho_{2}$ $\displaystyle=$ $\displaystyle 0,$ (29)
$\displaystyle\left(t\left(\rho_{1,t}+\rho_{2,t}\right)+\left(\rho_{1}+\rho_{2}\right)\right)x+\left(\rho_{1}v+\rho_{2}v\right)_{,t}+2t\left(\rho_{1}+\rho_{2}\right)v$
$\displaystyle=$ $\displaystyle 0,$ (30)

from where it follows the exact solution

$\rho_{i}\left(t\right)=\rho_{i0}t^{-1}\text{ and
}v\left(t,x\right)=\frac{\left(x+u_{0}\right)}{t}\text{. }$ (31)

####  3.1.5 Travellin waves $X_{1}+\beta X_{2}$

From the Lie symmetry $X_{1}+\beta X_{2}$ it follows
$\rho_{i}=\rho_{i}\left(w\right),~{}u=u\left(w\right)$ where $w=x-\beta t$
while the reduced system is

$\displaystyle\beta\rho_{1,w}-\left(\rho_{1}u\right)_{,w}$ $\displaystyle=$
$\displaystyle 0,$ (32)
$\displaystyle\beta\rho_{2,w}-\left(\rho_{2}u\right)_{,w}$ $\displaystyle=$
$\displaystyle 0,$ (33)
$\displaystyle\left(\left(\rho_{1}+\rho_{2}\right)u\right)_{,w}-\beta\left(\left(\rho_{1}+\rho_{2}\right)u^{2}+p_{1}+p_{2}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (34)

which provides the constant solution $\rho_{i}=\rho_{i0}$ and $u=u_{0}$, or
for $u=\beta$ it follows

$\rho_{1}\left(w\right)=\left(\rho_{10}-\frac{\kappa_{2}}{\kappa_{1}}\rho_{2}\left(w\right)\right)^{\frac{1}{\gamma_{1}}}$
(35)

where $\rho_{2}\left(w\right)$ is an arbitrary function.

Recall that in order the latter solution to be physically accepted
$\rho_{2}\left(w\right)$ should be defined such that $\rho_{1},\rho_{2}$ to be
real positive functions. That is a more general solution that the one found
before in [42]. The only solution which was found in [42] provides that one of
the $\rho_{i}\left(w\right)$ will be negative, or one of the pressure
$p_{i}\left(\rho_{i}\left(w\right)\right)$ will be negative. That is not
physically accepted, since there are not known physically fluids with negative
energy density or negative pressure. Although that kind of fluids are used in
theoretical astrophysics as toy models, they have not observed yet.

####  3.1.6 Reduction with $X_{3}+\beta X_{4}$

Application of the symmetry vector $X_{3}+\beta X_{4}$ gives the similarity
transformation $\rho_{i}=\rho_{i}\left(w\right),~{}u=\beta\ln
t+v\left(w\right)$ with $w=x-\beta\ln t$. The reduced system provides

$\rho_{i}\left(w\right)=\rho_{i0}V_{,w}~{},~{}v\left(w\right)=\left(w+\beta\right)+\frac{V_{0}-V\left(w\right)}{V_{,w}}$
(36)

where $V\left(w\right)$ satisfies the second-order ordinary differential
equation

$0=\left(\rho_{10}+\rho_{20}\right)\left(\beta\left(V_{,w}\right)^{3}-V_{,ww}\left(V-V_{0}\right)^{2}\right)+V_{,w}V_{,ww}\left(\kappa_{1}\left(\rho_{10}V_{,w}\right)^{\gamma_{1}}+\kappa_{2}\left(\rho_{20}V_{,w}\right)^{\gamma_{2}}\right),$
(37)

Equation (37) is an autonomous equation, hence we can define the new dependent
variable $\phi=V_{,w}$ and the new independent variable $z=V$.

Therefore, it follows

$0=\left(\rho_{10}+\rho_{20}\right)\left(\beta\left(\phi\left(z\right)\right)^{3}-\phi\left(z\right)\phi_{,z}\left(z-V_{0}\right)^{2}\right)+\phi_{,z}\left(\phi\left(z\right)\right)^{2}\left(\kappa_{1}\left(\rho_{10}\phi\left(z\right)\right)^{\gamma_{1}}+\kappa_{2}\left(\rho_{20}\phi\left(z\right)\right)^{\gamma_{2}}\right),$

which can be solved by quadratures.

####  3.1.7 Reduction with $X_{1}+\beta X_{4}$

From the symmetry vector $X_{1}+\beta X_{4}$ it follows
$\rho_{i}=\rho_{i}\left(w\right)$, $u=\beta t+V\left(w\right)$ where now
$w=x-\frac{\beta}{2}t^{2}$. The reduced system is

$\left(\rho_{1}v\right)_{,w}=0~{},~{}\left(\rho_{2}v\right)_{,w}=0~{},$ (38)
$\beta\left(\rho_{1}+\rho_{2}\right)+\left(\left(\rho_{1}+\rho_{2}\right)u^{2}+p_{1}+p_{2}\right)_{,w}=0$
(39)

that is
$\rho_{1}\left(w\right)=\rho_{10}v^{-1},~{}~{}\rho_{2}\left(w\right)=\rho_{20}v^{-1}$.
Hence $v\left(w\right)$ is given by the equation

$\left(\rho_{10}+\rho_{20}\right)w+v_{0}+\frac{1}{2}\left(\rho_{10}+\rho_{20}\right)v^{2}-\frac{\kappa_{1}\left(\rho_{10}\right)^{\gamma_{1}}}{1-\gamma_{1}}v^{1-\gamma_{1}}-\frac{\kappa_{1}\left(\rho_{20}\right)^{\gamma_{2}}}{1-\gamma_{2}}v^{1-\gamma_{2}}=0$
(40)

for $\gamma_{i}\neq 1$, or

$\left(\rho_{10}+\rho_{20}\right)w+v_{0}+\frac{1}{2}\left(\rho_{10}+\rho_{20}\right)v^{2}-\kappa_{1}\rho_{10}\ln
v-\frac{\kappa_{1}\left(\rho_{20}\right)^{\gamma_{2}}}{1-\gamma_{2}}v^{1-\gamma_{2}}=0.$
(41)

for $\gamma_{1}=1,~{}\gamma_{2}\neq 1$ or

$\left(\rho_{10}+\rho_{20}\right)w+v_{0}+\frac{1}{2}\left(\rho_{10}+\rho_{20}\right)v^{2}-\kappa_{1}\rho_{10}\ln
v^{\kappa_{1}\rho_{10}+\kappa_{2}\rho_{29}}=0.$ (42)

for $\gamma_{i}=1$. This is a new solutions which has not found before in the
literature.

Nevertheless for specific values of the polytropic exponents the dynamical
system (14)-(16) admits additional symmetries.

###  3.2 Polytropic exponents $\gamma_{1}=\gamma_{2}=\gamma~{},~{}\gamma\neq
1$

Consider now the case where the polytropic exponents $\gamma_{i}$ are equal,
that is, $\gamma_{1}=\gamma_{2}=\gamma$. In such scenario the Lie point
symmetries of the two-phase model (14)-(16) are the vector fields
$X_{1},~{}X_{2},~{}X_{3},~{}X_{4},$ plus the additional symmetries

$~{}X_{5}=\frac{2}{\gamma-1}\left(\rho_{1}\partial_{\rho_{1}}+\rho_{2}\partial_{\rho_{2}}+u\partial
u-t\partial_{t}\right)~{},~{}X_{6}=f_{1}\left(\rho_{1},\rho_{2};\kappa_{1},\kappa_{2},\gamma\right)\partial_{\rho_{1}}+f_{2}\left(\rho_{1},\rho_{2};\kappa_{1},\kappa_{2},\gamma\right)\partial_{\rho_{2}}.$

Vector field $X_{5}$ is a scaling symmetry, while $X_{6}$ is a rotation
symmetry which indicates the invariance of the two-phase model if
$\left(\rho_{1},\rho_{2}\right)\rightarrow\left(\bar{\rho}_{1},\bar{\rho}_{2}\right)$.
However, functions
$f_{i}\left(\rho_{1},\rho_{2};\kappa_{1},\kappa_{2},\gamma\right)$ are not
expressed always into closed-form expressions. Thus for $\gamma=2$ it follows

$\displaystyle
f_{1}\left(\rho_{1},\rho_{2};\kappa_{1},\kappa_{2},\gamma\right)$
$\displaystyle=$
$\displaystyle\frac{\kappa_{1}\rho_{1}^{2}-2\kappa_{2}\rho_{1}\rho_{2}-\kappa_{2}\rho_{2}^{2}}{\kappa_{1}\rho_{1}-\kappa_{2}\rho_{2}},$
(43) $\displaystyle
f_{2}\left(\rho_{1},\rho_{2};\kappa_{1},\kappa_{2},\gamma\right)$
$\displaystyle=$
$\displaystyle\frac{\kappa_{1}\rho_{1}^{2}+2\kappa_{1}\rho_{1}\rho_{2}-\kappa_{2}\rho_{2}^{2}}{\kappa_{1}\rho_{1}-\kappa_{2}\rho_{2}}.$
(44)

For arbitrary polytropic index $\gamma$ we found that
$f_{i}=f_{i}\left(\frac{\rho_{2}}{f_{1}}\right)\rho_{i}$, which indicates that
$\left[X_{5},X_{6}\right]=0$.

Symmetry vector $X_{6}$ has not been derived before in the literature and it
is a new symmetry. In Table 3 the commutators of the admitted Lie symmetries
is presented while in Table 4 the corresponding adjoint representation is
derived.

From Tables 3 and 4 the one-dimensional optimal system is calculated consisted
by the vector fields

$\displaystyle\left\\{X_{1}\right\\}~{},~{}\left\\{X_{2}\right\\}~{},~{}\left\\{X_{3}\right\\}~{},~{}\left\\{X_{4}\right\\}~{},~{}\left\\{X_{5}\right\\}~{},~{}\left\\{X_{6}\right\\}~{},$
$\displaystyle\left\\{X_{1}+\beta X_{2}\right\\}~{},~{}\left\\{X_{3}+\beta
X_{4}\right\\}~{},\left\\{X_{1}+\beta X_{4}\right\\}~{},$
$\displaystyle\left\\{X_{2}+\beta X_{5}\right\\}~{},~{}\left\\{X_{3}+\beta
X_{5}\right\\}~{},~{}\left\\{X_{1}+\beta X_{6}\right\\}~{},$
$\displaystyle\left\\{X_{2}+\beta X_{6}\right\\}~{},~{}\left\\{X_{3}+\beta
X_{6}\right\\}~{},\left\\{X_{4}+\beta X_{6}\right\\}~{},$
$\displaystyle\left\\{X_{5}+\beta X_{6}\right\\}~{}~{},\left\\{X_{2}+\beta
X_{5}+\delta X_{6}\right\\}~{},~{}\left\\{X_{3}+\beta X_{5}+\delta
X_{6}\right\\}.$

We continue our analysis with the new reductions, while for simplicity on our
calculations and on the presentation we select the polytropic index
$\gamma=2$.

Table 3: Commutator table for the Lie point symmetries of the two-phase flow system (14)-(16) for equal indices $\gamma_{i}$ $\left[X_{I},X_{J}\right]$ | $\mathbf{X}_{1}$ | $\mathbf{X}_{2}$ | $\mathbf{X}_{3}$ | $\mathbf{X}_{4}$ | $\mathbf{X}_{5}$ | $\mathbf{X}_{6}$
---|---|---|---|---|---|---
$\mathbf{X}_{1}$ | $0$ | $0$ | $X_{1}$ | $X_{2}$ | $-X_{1}$ | $0$
$\mathbf{X}_{2}$ | $0$ | $0$ | $X_{2}$ | $0$ | $0$ | $0$
$\mathbf{X}_{3}$ | $-X_{1}$ | $-X_{2}$ | $0$ |  $0$ |  $0$ | $0$
$\mathbf{X}_{4}$ | $-X_{2}$ | $0$ | $0$ | $0$ | $X_{4}$ | $0$
$\mathbf{X}_{5}$ | $X_{1}$ | $0$ | $0$ | $-X_{4}$ | $0$ | $0$
$\mathbf{X}_{6}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
Table 4: Adjoint representation for the Lie point symmetries of the two-phase flow system (14)-(16) for equal indices $\gamma_{i}$ $Ad\left(e^{\left(\varepsilon\mathbf{X}_{i}\right)}\right)\mathbf{X}_{j}$ | $\mathbf{X}_{1}$ | $\mathbf{X}_{2}$ | $\mathbf{X}_{3}$ | $\mathbf{X}_{4}$ | $\mathbf{X}_{5}$ | $\mathbf{X}_{6}$
---|---|---|---|---|---|---
$\mathbf{X}_{1}$ | $X_{1}$ | $X_{2}$ | $-\varepsilon X_{1}+X_{3}$ | $-\varepsilon X_{2}+X_{4}$ | $\varepsilon X_{1}+X_{5}$ | $\mathbf{X}_{6}$
$\mathbf{X}_{2}$ | $X_{1}$ | $X_{2}$ | $-\varepsilon X_{2}+X_{3}$ | $X_{4}$ | $X_{5}$ | $\mathbf{X}_{6}$
$\mathbf{X}_{3}$ | $e^{\varepsilon}X_{1}$ | $e^{\varepsilon}X_{2}$ | $X_{3}$ | $X_{4}$ | $X_{5}$ | $\mathbf{X}_{6}$
$\mathbf{X}_{4}$ | $X_{1}+\varepsilon X_{2}$ | $X_{2}$ | $X_{3}$ | $X_{4}$ | $-\varepsilon X_{4}+X_{5}$ | $\mathbf{X}_{6}$
$\mathbf{X}_{5}$ | $e^{\varepsilon}X_{1}$ | $X_{2}$ | $X_{3}$ | $e^{\varepsilon}X_{4}$ | $X_{5}$ | $\mathbf{X}_{6}$
$\mathbf{X}_{6}$ | $X_{1}$ | $X_{2}$ | $X_{3}$ | $X_{4}$ | $X_{5}$ | $\mathbf{X}_{6}$

####  3.2.1 Reduction with $X_{5}$

Application of the symmetry vector $X_{5}$ provides the similarity
transformation $\rho_{i}=\rho_{i}\left(x\right)t^{-2}$ and
$u=u\left(x\right)t^{-2}$, where

$\rho_{i}\left(x\right)=\rho_{i0}\exp\left(\int\frac{2-u_{,x}}{u}dx\right)$
(45)

where $u=\left(2v\left(x\right)+v_{0}\right)\left(v_{,x}\right)^{-1}$ and
$v\left(x\right)$ is given by the differential equation

$v_{,xx}\left(\left(\kappa_{1}\rho_{10}^{2}+\kappa_{2}\rho_{20}^{2}\right)\left(v_{,x}\right)^{3}-\left(2v+v_{0}\right)^{2}\left(\rho_{10}+\rho_{20}\right)\right)+\left(v_{,x}\right)^{2}\left(\rho_{10}+\rho_{20}\right)\left(2v+v_{0}\right)=0,$
(46)

the later equation can be integrated by quadratures.

####  3.2.2 Reduction with $X_{2}+\beta X_{5}$

From the Lie symmetry vector $X_{2}+\beta X_{5}$ we find the similarity
transformation
$\rho_{i}=\rho_{i}\left(w\right)t^{-2}~{},~{}u=\bar{u}\left(w\right)t^{-1}$
where $w=x+\frac{1}{\beta}\ln t$. The reduced system provides

$\rho_{i}=\rho_{i0}\exp\left(\beta\int\frac{2-\bar{u}_{,w}}{1+\beta\bar{u}_{,w}}dw\right)$
(47)

where

$\bar{u}=-\frac{1}{\beta}+\left(2v\left(w\right)+v_{0}\right)\left(v_{,w}\right)^{-1}$
(48)

where now $v\left(x\right)$ is given by the second-order differential equation

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\beta
v_{,ww}\left(2\left(\kappa_{1}\rho_{10}^{2}+\kappa_{2}\rho_{20}^{2}\right)\left(v_{,w}\right)^{3}-\left(2v+v_{0}\right)^{2}\left(\rho_{10}+\rho_{20}\right)\right)$
(49)
$\displaystyle+\left(v_{,w}\right)^{2}\left(\rho_{10}+\rho_{20}\right)\left(\beta\left(2v+v_{0}\right)+v_{,w}\right).$

which again can be solved by quadratures.

####  3.2.3 Reduction with $X_{3}+\beta X_{5}$

The Lie symmetry vector $\left\\{X_{3}+\beta X_{5}\right\\}$ provides the
similarity transformation
$\rho_{i}=\rho_{i}\left(w\right)t^{-\frac{2\beta}{\beta-1}},~{}u=\bar{u}\left(w\right)t^{-\frac{\beta}{\beta-1}}$
with new independent variable $w=xt^{\frac{1}{\beta-1}}$, while the reduced
system provides

$\rho_{i}\left(w\right)=\rho_{i0}v_{x}~{},~{}\bar{u}=-\frac{1}{\beta}+\left(2v\left(w\right)+v_{0}\right)\left(v_{,w}\right)^{-1}$

where $v\left(w\right)$ satisfies the second-order differential equation

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
v_{,ww}\left(\left(\rho_{10}+\rho_{20}\right)\left(v\left(2\beta+1\right)-v_{0}\left(1-\beta\right)\right)-2\left(\kappa_{1}\rho_{10}^{2}+\kappa_{2}\rho_{20}^{2}\right)\left(v_{,w}\right)^{3}\left(1-\beta\right)^{2}\right)+$
(50)
$\displaystyle-w\left(v_{,w}\right)^{3}\beta\left(\rho_{10}+\rho_{20}\right)-\left(\rho_{10}+\rho_{20}\right)\beta\left(v_{,w}\right)^{2}\left(v\left(2\beta+1\right)-v_{0}\left(1-\beta\right)\right)$

which can be solved by quadratures.

####  3.2.4 Reduction with $X_{1}+\beta X_{6}$

In order to proceed with the reduction with the use of the rest of the
symmetry vectors we define the new variables

$\rho\left(t,x\right)=\rho_{1}\left(t,x\right)+\rho_{2}\left(t,x\right),$ (51)
$p\left(t,x\right)=\kappa_{1}\rho_{1}^{2}\left(t,x\right)+\kappa_{2}\rho_{2}^{2}\left(t,x\right),$
(52)

where now system (14)-(16) becomes

$\displaystyle\rho_{,t}+\left(\rho u\right)_{,x}$ $\displaystyle=$
$\displaystyle 0,$ (53) $\displaystyle p_{,t}+\left(pu\right)_{,x}+pu_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (54) $\displaystyle\left(\rho
u\right)_{,t}+\left(\rho u^{2}+p\right)_{,x}$ $\displaystyle=$ $\displaystyle
0.$ (55)

In the new coordinates, symmetries $X_{5,}~{}X_{6}$ are written
$X_{5}=-2\rho\partial_{\rho}+u\partial_{u}-t\partial_{t}$ and
$X_{6}=\rho\partial_{\rho}+p\partial_{p}\,$.

Thus, application of the Lie symmetry vector $X_{1}+\beta X_{6}$ gives

$\rho=e^{t}\rho\left(x\right)~{},~{}p=e^{t}p\left(x\right)~{},~{}u=u\left(x\right)$
(56)

while the reduced system is

$\displaystyle\rho+\left(\rho u\right)_{,x}$ $\displaystyle=$ $\displaystyle
0,$ (57) $\displaystyle p+\left(pu\right)_{,x}+pu_{,x}$ $\displaystyle=$
$\displaystyle 0,$ (58) $\displaystyle\rho u+\left(\rho u^{2}+p\right)_{,x}$
$\displaystyle=$ $\displaystyle 0$ (59)

which gives the solution

$\rho\left(x\right)=\frac{2pu_{,x}+\beta
p}{u^{2}u_{,x}}~{},~{}p=p_{0}\exp\left(-\int\frac{\beta+2u_{,x}}{u}dx\right)$
(60)

where $u\left(x\right)$ is given by the polynomial

$u_{0}\beta u^{4}-2\beta u-x-\beta x_{0}=0.$ (61)

This is also a new solution for the two-phase fluids model. As also all the
following similarity solutions are new and they have not been calculated
before.

####  3.2.5 Reduction with $X_{2}+\beta X_{6}$

The similarity transformation which follows from $X_{2}+\beta X_{6}$ is
$\rho=e^{\beta x}\bar{\rho}\left(t\right)~{},~{}p=e^{\beta
x}\bar{p}\left(t\right)$ and $u=u\left(t\right)$ while system (53)-(55)
becomes

$\displaystyle\bar{\rho}_{t}+\beta\bar{\rho}u$ $\displaystyle=$ $\displaystyle
0,$ (62) $\displaystyle\bar{p}_{t}+\beta\bar{p}u$ $\displaystyle=$
$\displaystyle 0,$ (63)
$\displaystyle\left(\bar{\rho}u\right)_{,t}+\beta\bar{\rho}u^{2}+\beta\bar{p}$
$\displaystyle=$ $\displaystyle 0,$ (64)

from where we infer the physically accepted solution

$\displaystyle u\left(t\right)$ $\displaystyle=$ $\displaystyle
u_{0}t+u_{1},~{}$ (65) $\displaystyle\bar{p}\left(t\right)$ $\displaystyle=$
$\displaystyle
p_{0}\exp\left(-\frac{\beta}{2}\lambda_{0}t^{2}-\beta\lambda_{1}t\right)~{},~{}$
(66) $\displaystyle\bar{\rho}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{\beta
p_{0}}{\lambda_{0}}\exp\left(-\frac{\beta}{2}\lambda_{0}t^{2}-\beta\lambda_{1}t\right)~{}.$
(67)

####  3.2.6 Reduction with $X_{3}+\beta X_{6}$

From the symmetry vector $X_{3}+\beta X_{6}$ we find
$\rho=t^{\beta}\bar{\rho}\left(w\right),~{}p=t^{\beta}\bar{p}\left(w\right)$
and $u=u\left(w\right)$ with $w=\frac{x}{t}$. Hence, from (53)-(55) we find
the solution

$\bar{\rho}\left(w\right)=\rho_{0}\exp\left(\frac{\beta+u_{,w}}{w-u}\right)~{},~{}\bar{p}=p_{1}\exp\left(\int\frac{\beta-2u_{w}}{w-u}\right),$
(68)

where $u\left(w\right)=w+v\left(w\right)$, while $v\left(w\right)$ is given by
the integral
$\int^{v\left(w\right)}\left(f\left(s\right)\right)^{-1}ds-\left(w-w_{0}\right)=0$,
where $f\left(s\right)$ is a solution of the polynomial equation

$-s^{-2-\frac{3}{2}\beta}\left(f\left(s\right)+1\right)^{-1-\frac{3}{4}\beta}\left(3f\left(s\right)+1\right)^{\frac{\beta}{4}}\left(2f\left(s\right)+\beta+2\right)^{\frac{\beta}{2}+1}+p_{0}=0.$
(69)

####  3.2.7 Reduction with $X_{4}+\beta X_{6}$

Application of the symmetry vector $X_{4}+\beta X_{6}$ gives the
transformation
$\rho=e^{\beta\frac{x}{t}}\bar{\rho}\left(t\right)~{},~{}p=e^{\beta\frac{x}{t}}\bar{p}\left(t\right)$
and $u=\frac{x}{t}+v\left(t\right)$, where now the reduced system is derived

$\displaystyle t\bar{\rho}_{,t}+\rho\left(1+\beta v\right)$ $\displaystyle=$
$\displaystyle 0,$ (70) $\displaystyle t\bar{p}_{,t}+p\left(2+\beta v\right)$
$\displaystyle=$ $\displaystyle 0,$ (71) $\displaystyle
t\left(\bar{\rho}v\right)_{,t}+\rho\left(\beta v^{2}+2v\right)+\beta\bar{p}$
$\displaystyle=$ $\displaystyle 0,$ (72)

that is

$\displaystyle v\left(t\right)$ $\displaystyle=$
$\displaystyle\frac{u_{0}}{t}+\frac{u_{1}}{t}\ln t~{},$ (73)
$\displaystyle\bar{p}\left(t\right)$ $\displaystyle=$ $\displaystyle
p_{0}t^{\frac{\beta}{t}-2}\exp\left(\frac{\beta}{t}\left(u_{0}+u_{1}\right)\right)~{},$
(74) $\displaystyle\bar{\rho}\left(t\right)$ $\displaystyle=$
$\displaystyle-\beta\frac{p_{0}}{u_{1}}t^{\frac{\beta}{t}-1}\exp\left(\frac{\beta}{t}\left(u_{0}+u_{1}\right)\right).$
(75)

Thus, in order the solution to be physically accepted, $p_{0}>0$ and $\beta
u_{1}<0$.

####  3.2.8 Reduction with $X_{5}+\beta X_{6}$

The Lie symmetry vector $X_{5}+\beta X_{6}$ provides
$\rho=t^{2-\beta}\bar{\rho}\left(x\right),~{}p=t^{-\beta}\bar{p}\left(x\right)$
and $u=t^{-1}\bar{u}\left(x\right)$, while the reduced system is

$\displaystyle\left(2-\beta\right)\bar{\rho}+\left(\bar{\rho}\bar{u}\right)_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (76)
$\displaystyle-\beta\bar{p}+\left(\bar{p}u\right)_{,x}+\bar{p}\bar{u}_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (77)
$\displaystyle\left(1-\beta\right)\bar{\rho}\bar{u}+\left(\bar{\rho}\bar{u}^{2}+\bar{p}\right)_{,x}$
$\displaystyle=$ $\displaystyle 0.$ (78)

that is,

$\displaystyle\bar{\rho}\left(x\right)$ $\displaystyle=$
$\displaystyle\rho_{0}\exp\left(-\int\frac{2-\beta+\bar{u}_{,x}}{\bar{u}}dx\right)~{},$
(79) $\displaystyle\bar{p}\left(x\right)$ $\displaystyle=$ $\displaystyle
p_{0}\exp\left(-\int\frac{2\bar{u}_{,x}-\beta}{\bar{u}}dx\right),$ (80)

where $\bar{u}\left(w\right)$ is given by the expression

$\int^{\bar{u}\left(x\right)}\left(g\left(s\right)\right)^{-1}ds=\left(x-x_{0}\right),$
(81)

where

$-2\left(sg\left(s\right)-1\right)^{\frac{4-3\beta}{\beta}}\left(3g\left(s\right)-2\right)^{-\frac{2}{\beta}\left(2-\beta\right)}\left(2g\left(s\right)-\beta\right)+g_{0}=0\text{.}$

####  3.2.9 Reduction with $X_{2}+\beta X_{5}+\delta X_{6}$

From $X_{2}+\beta X_{5}+\delta X_{6}$ it follows
$\rho=t^{2-\frac{\delta}{\beta}}\bar{\rho}\left(w\right),~{}p=t^{-\frac{\delta}{\beta}}\bar{p}\left(w\right),~{}u=t^{-1}\bar{u}\left(w\right)$,
where $w=x+\frac{1}{\beta}\ln t$. Thus, system (53)-(55) is reduced

$\displaystyle\bar{\rho}_{,w}+\left(2\beta-\delta\right)\bar{\rho}+\beta\left(\bar{\rho}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (82)
$\displaystyle\bar{p}_{,w}-\delta\bar{p}+\beta\left(\left(\bar{p}\bar{u}\right)_{,w}+\bar{p}\bar{u}_{,w}\right)$
$\displaystyle=$ $\displaystyle 0,$ (83)
$\displaystyle\left(\bar{\rho}\bar{u}\right)_{,w}+\left(\beta-\delta\right)\bar{\rho}\bar{u}+\beta\left(\bar{\rho}\bar{u}^{2}+\bar{p}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0.$ (84)

The latter system can be integrated further and its solution is expressed in
terms of quadratures.

####  3.2.10 Reduction with $X_{3}+\beta X_{5}+\delta X_{6}$

Similarly, from the symmetry vector $X_{3}+\beta X_{5}+\delta X_{6}$ it
follows
$\rho=t^{\frac{2\beta-\delta}{\beta-1}}\bar{\rho}\left(w\right),~{}p=t^{-\frac{\delta}{\beta-1}}\bar{p}\left(w\right),~{}u=t^{-\frac{\beta}{\beta-1}}\bar{u}\left(w\right)$
with $w=xt^{\frac{1}{\beta-1}}$. Therefore, the reduced system is

$\displaystyle\bar{\rho}_{,w}+\left(2\beta-\delta\right)\bar{\rho}+\left(\beta-1\right)\left(\bar{\rho}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (85)
$\displaystyle\bar{p}_{,w}-\delta\bar{p}+\left(\beta-1\right)\left(\left(\bar{p}\bar{u}\right)_{,w}+p\bar{u}_{,w}\right)$
$\displaystyle=$ $\displaystyle 0,$ (86)
$\displaystyle\left(\bar{\rho}\bar{u}\right)_{,w}+\left(\beta-\delta\right)\bar{\rho}\bar{u}+\left(\beta-1\right)\left(\bar{\rho}\bar{u}^{2}+\bar{p}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0.$ (87)

which again can be integrated by quadratures.

We continue our analysis with the special case of the polytropic exponents
$\gamma_{1}=\gamma_{2}=1$.

###  3.3 Polytropic exponents $\gamma_{1}=\gamma_{2}=1$

In the special case where the two polytropic exponents are equal with one,
that is, $\gamma_{1}=\gamma_{2}=1$. Similarly with the previous case we define
the new dependent variables

$\rho\left(t,x\right)=\rho_{1}\left(t,x\right)+\rho_{2}\left(t,x\right),$ (88)
$p\left(t,x\right)=\kappa_{1}\rho_{1}\left(t,x\right)+\kappa_{2}\rho_{2}\left(t,x\right),$
(89)

where the dynamical system (14)-(16) takes the following form

$\displaystyle\rho_{,t}+\left(\rho u\right)_{,x}$ $\displaystyle=$
$\displaystyle 0,$ (90) $\displaystyle p_{,t}+\left(pu\right)_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (91) $\displaystyle\left(\rho
u\right)_{,t}+\left(\rho u^{2}+p\right)_{,x}$ $\displaystyle=$ $\displaystyle
0.$ (92)

The set of equations (90)-(92) admits a six dimensional Lie algebra consisted
by the symmetry vectors
$\left\\{X_{1},X_{2},X_{3},X_{4},X_{5},X_{6}\right\\}$. Thus, there is not any
different between the symmetry classification of the two cases $\gamma=1$ and
$\gamma\neq 1$, as it was found before in [42]. We proceed with the
application of the symmetry vectors for the derivation of similarity
solutions.

####  3.3.1 Reduction with $X_{5}$

Application of the symmetry vector $X_{5}$ gives
$\rho=t^{-2}\bar{\rho}\left(x\right),~{}p=p\left(x\right)$ and
$u=t^{-1}\bar{u}\left(x\right)$ where the reduced system is

$\displaystyle\left(\bar{\rho}\bar{u}\right)_{,x}+2\bar{\rho}$
$\displaystyle=$ $\displaystyle 0,$ (93)
$\displaystyle\left(p\bar{u}\right)_{,x}$ $\displaystyle=$ $\displaystyle 0,$
(94)
$\displaystyle\bar{\rho}\bar{u}+\left(\bar{\rho}\bar{u}^{2}+p\right)_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (95)

from where we find $\bar{\rho}=v_{,x}$, $p=p_{0}u^{-1}$ and
$\bar{u}=\frac{2\left(v_{0}-v\right)}{v_{,x}}$, where $v\left(x\right)$
satisfies the differential equation

$6\rho_{0}\left(v-v_{0}\right)+\frac{p_{0}}{2\left(v-v_{0}\right)^{2}}\left(\left(v_{0}-v\right)v_{,xx}+v_{,x}^{2}\right)-\frac{4\rho_{0}}{v_{,x}^{2}}\left(v-v_{0}\right)^{2}v_{,xx}=0,$
(96)

which can be solved by quadratures.

####  3.3.2 Reduction with $X_{2}+\beta X_{5}$

The vector field $X_{2}+\beta X_{5}$ provides
$\rho=t^{2}\rho\left(w\right),~{}p=p\left(w\right)$ and
$u=t^{-1}u\left(w\right)$ where now the independent variable is defined as
$w=x+\frac{1}{\beta}\ln t$. Hence, by replacing in (90)-(92) we find

$\rho\left(w\right)=\rho_{0}v_{,w}~{},~{}p\left(w\right)=p_{0}\left(1+\beta
u\right)^{-1}~{},~{}u\left(w\right)=-\frac{1}{\beta}+\frac{2\left(v_{0}-v\right)}{v_{,x}}$
(97)

where $v\left(x\right)$ satisfies the second-order ode

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
v_{xx}\left(v_{0}-v\right)\left(8\beta\rho_{0}\left(v-v_{0}\right)^{3}+p_{0}\left(v_{,x}\right)^{2}\right)+p_{0}\left(v_{,x}\right)^{4}+$
(98)
$\displaystyle+2\rho_{0}\left(v-v_{0}\right)^{2}\left(v_{x}\right)^{3}+12\beta\rho_{0}\left(v_{,x}\right)^{2}\left(v-v_{0}\right)^{3}$

which can be solved by quadratures.

####  3.3.3 Reduction with $X_{3}+\beta X_{5}$

From the vector field $X_{3}+\beta X_{5}$ it follows
$\rho=t^{\frac{2\beta}{\beta-1}}\bar{\rho}\left(w\right),~{}p=p\left(w\right)$
and $u=t^{-\frac{\beta}{\beta-1}}\bar{u}\left(w\right)$ where
$w=xt^{\frac{1}{\beta-1}}$. The reduced system is

$\displaystyle\left(x+\left(\beta-1\right)\bar{u}\right)\bar{\rho}_{,x}+2\beta\bar{\rho}+\left(\beta-1\right)\bar{\rho}\bar{u}_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (99)
$\displaystyle\left(x+\left(\beta-1\right)\bar{u}\right)p_{,x}+\left(\beta-1\right)p\bar{u}_{,x}$
$\displaystyle=$ $\displaystyle 0,$ (100) $\displaystyle
x\bar{u}\bar{\rho}_{x}+x\bar{\rho}\bar{u}_{,x}+\beta\bar{u}\bar{\rho}+\left(\beta-1\right)\left(\bar{\rho}\bar{u}^{2}+p\right)_{,x}$
$\displaystyle=$ $\displaystyle 0.$ (101)

For $\beta=1$, the closed-form solution of the later system is
$p\left(w\right)=p_{0},~{}\bar{\rho}=\rho_{0}w^{-2}$ and $\bar{u}=u_{0}w$.

####  3.3.4 Reduction with $X_{1}+\beta X_{6}$

The similarity transformation which follows from the symmetry vector
$X_{1}+\beta X_{6}$ is $\rho=e^{\beta t}\bar{\rho}\left(x\right),~{}p=e^{\beta
t}\bar{p}\left(x\right),~{}u=\bar{u}\left(x\right)$, where by replacing in
(90)-(92) we find

$\bar{\rho}=\rho_{0}v_{,x}~{},~{}\bar{p}=p_{0}v_{,x},~{}\bar{u}=\frac{v_{0}-\beta
v}{v_{,x}}$ (102)

where $v\left(x\right)$ is a solution of the second-order differential
equation

$v_{,xx}\left(p_{0}\left(v_{,x}\right)^{2}-\rho_{0}\left(\beta
v-v_{0}\right)^{2}\right)+\beta\rho_{0}v_{,x}^{2}\left(\beta
v-v_{0}\right)=0\text{.}$ (103)

which can be solved by quadratures.

####  3.3.5 Reduction with $X_{2}+\beta X_{6}$

From the symmetry vector $X_{2}+\beta X_{6}$ it follows $\rho=e^{\beta
x}\bar{\rho}\left(t\right),~{}p=e^{\beta x}\bar{p}\left(t\right)$ and
$u=u\left(t\right)$, where by replacing in the system (90)-(92) it follows

$\displaystyle u\left(x\right)$ $\displaystyle=$ $\displaystyle u_{1}x+u_{0},$
(104) $\displaystyle\bar{p}\left(x\right)$ $\displaystyle=$ $\displaystyle
p_{0}\exp\left(-\beta\left(\frac{1}{2}u_{1}x^{2}+u_{0}x\right)\right),$ (105)
$\displaystyle\bar{\rho}$ $\displaystyle=$ $\displaystyle-\frac{\beta
p_{0}}{u_{1}}\exp\left(-\beta\left(\frac{1}{2}u_{1}x^{2}+u_{0}x\right)\right).$
(106)

####  3.3.6 Reduction with $X_{3}+\beta X_{6}$

The Lie symmetry vector $X_{3}+\beta X_{6}$ provides
$\rho=t^{\beta}\bar{\rho}\left(w\right),~{}p=t^{\beta}\bar{p}\left(w\right),~{}u=u\left(w\right)$
where $w=\frac{x}{t}$. Hence, from the system (90)-(92) we find

$\bar{\rho}=\rho_{0}\exp\left(-\int\frac{\beta+u_{,w}}{u-w}dw\right)~{},~{}\bar{p}=p_{0}\exp\left(-\int\frac{\beta+u_{,w}}{u-w}dw\right),$
(107)

and $u\left(w\right)$ is given by the algebraic equation

$\arctan
h\left(\frac{e^{\frac{p_{0}}{2}\beta}\left(u\left(w\right)-w\right)}{\sqrt{\beta+1}}\right)\beta+\sqrt{\left(\beta+1\right)e}\left(u_{0}-u\left(w\right)\right)=0.$
(108)

In the special case where $\beta=-1\,,$ function $u\left(w\right)$ is
expressed as follows

$u_{\pm}\left(w\right)=\frac{1}{2}\left(\left(w-w_{0}\right)\pm\sqrt{\left(w-w_{0}\right)^{2}-4p_{0}}\right).$
(109)

####  3.3.7 Reduction with $X_{4}+\beta X_{6}$

The similarity transformation which correspond to the symmetry vector
$X_{4}+\beta X_{6}$ is
$\rho=e^{\beta\frac{x}{t}}\bar{\rho}\left(t\right),~{}p=e^{\beta\frac{x}{t}}\bar{p}\left(t\right),~{}u=\frac{x}{t}+v\left(t\right)$.
Therefore by replacing in the original system (90)-(92) the closed-form
solution it follows

$\displaystyle v\left(t\right)$ $\displaystyle=$
$\displaystyle\frac{u_{1}}{t}+u_{0}~{},~{}$ (110)
$\displaystyle\bar{p}\left(t\right)$ $\displaystyle=$ $\displaystyle
t^{-1-u_{0}\beta}\exp\left(\frac{u_{1}\beta}{t}\right)~{},~{}$ (111)
$\displaystyle\bar{\rho}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{\beta
p_{0}}{u_{0}}x^{-1-u_{0}\beta}\exp\left(\frac{u_{1}\beta}{t}\right).$ (112)

####  3.3.8 Reduction with $X_{5}+\beta X_{6}$

From the Lie symmetry $X_{5}+\beta X_{6}$ we calculate
$\rho=t^{2-\beta}\bar{\rho}\left(x\right),~{}p=t^{-\beta}\bar{p}\left(x\right)$
and $u=t^{-1}\bar{u}\left(x\right)$. Thus, from (90)-(92) we calculate

$\bar{p}\left(x\right)=p_{0}\exp\left(-\int\frac{\bar{u}_{,x}-\beta}{\bar{u}}dx\right)~{},~{}\bar{\rho}\left(x\right)=\frac{\bar{u}_{,x}\bar{p}-\beta\bar{p}}{\bar{u}_{,x}\bar{u}^{2}-\bar{u}^{2}}.$
(113)

while $\bar{u}\left(x\right)$ is expressed in quadratures, that is

$\int_{0}^{\bar{u}\left(x\right)}\left(e^{r\left(s\right)}+\beta\right)^{-1}ds=x-x_{0},$

where $r\left(s\right)$ is a solution of the algebraic equation

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
e^{r\left(s\right)}\beta\left(Z+p_{0}\right)+\left(\beta-1\right)\left(r\left(s\right)\beta+p_{0}\beta+1\right)$
(114)
$\displaystyle-\beta\left(e^{r\left(s\right)}+\beta-1\right)\ln\left(e^{r\left(s\right)}+\beta-1\right)-2\left(\beta-1\right)\left(e^{r\left(s\right)}+\beta-1\right)\ln
s.$

In the special case where $\beta=1$ the closed-form solution follow

$\bar{\rho}=\frac{\rho_{0}}{\left(x-x_{0}\right)^{2}}~{},~{}\bar{u}=x-x_{0}~{},~{}\bar{p}=\rho_{0}.$
(115)

####  3.3.9 Reduction with $X_{2}+\beta X_{5}+\delta X_{6}$

The similarity transformation which correspond to the vector field
$X_{2}+\beta X_{5}+\delta X_{6}$ is
$\rho=t^{2-\frac{\delta}{\beta}}\bar{\rho}\left(w\right),~{}p=t^{-\frac{\delta}{\beta}}\bar{p}\left(w\right),~{}u=t^{-1}\bar{u}\left(w\right)$,
where $w=x+\frac{1}{\beta}\ln t$. The reduced system is

$\displaystyle\bar{\rho}_{,w}+\left(2\beta-\delta\right)\bar{\rho}+\beta\left(\bar{\rho}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (116)
$\displaystyle\bar{p}_{,w}-\delta\bar{p}+\beta\left(\bar{p}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (117)
$\displaystyle\left(\bar{\rho}\bar{u}\right)_{,w}+\left(\beta-\delta\right)\bar{\rho}\bar{u}+\beta\left(\bar{\rho}\bar{u}^{2}+\bar{p}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (118)

which can be solved by quadratures.

In the limit where $\beta=\frac{\delta}{2}$, the solution is

$\bar{\rho}=\frac{\rho_{0}}{2+\delta
u}~{},~{}\bar{p}\left(x\right)=\rho_{0}\frac{\left(\delta\bar{u}+2\right)\bar{u}_{,w}-\delta\bar{u}}{\delta^{2}\left(\bar{u}-2\right)}$
(119)

where $\bar{u}\left(w\right)$ is given by the first-order differential
equation

$\frac{df\left(\bar{u}\right)}{d\bar{u}}f-\frac{2\left(f\left(\bar{u}\right)-2\right)\left(\left(\bar{u}\delta+2\right)f\left(\bar{u}\right)-\left(2f\left(\bar{u}\right)-1\right)\delta\bar{u}-3f\left(\bar{u}\right)\right)}{\left(\bar{u}\delta+4\right)\left(\bar{u}\delta+2\right)},~{}f\left(\bar{u}\right)=\bar{u}_{,w}.$
(120)

####  3.3.10 Reduction with $X_{3}+\beta X_{5}+\delta X_{6}$

The symmetry vector $X_{3}+\beta X_{5}+\delta X_{6}$ provides the similarity
transformation
$\rho=t^{\frac{2\beta-\delta}{\beta-1}}\bar{\rho}\left(w\right),~{}p=t^{-\frac{\delta}{\beta-1}}\bar{p}\left(w\right),~{}u=t^{-\frac{\beta}{\beta-1}}\bar{u}\left(w\right)$
with $w=xt^{\frac{1}{\beta-1}}$. Therefore, the reduced system is

$\displaystyle\bar{\rho}_{,w}+\left(2\beta-\delta\right)\bar{\rho}+\left(\beta-1\right)\left(\bar{\rho}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (121)
$\displaystyle\bar{p}_{,w}-\delta\bar{p}+\left(\beta-1\right)\left(\bar{p}\bar{u}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0,$ (122)
$\displaystyle\left(\bar{\rho}\bar{u}\right)_{,w}+\left(\beta-\delta\right)\bar{\rho}\bar{u}+\left(\beta-1\right)\left(\bar{\rho}\bar{u}^{2}+\bar{p}\right)_{,w}$
$\displaystyle=$ $\displaystyle 0.$ (123)

which can be integrated by quadratures.

When $\beta=1$ and $\delta=-1$ the similarity transformation is
$\rho=x^{-3}\bar{\rho}\left(t\right),~{}p=\bar{p}\left(t\right)$ and
$u=x\bar{u}\left(t\right)$, where now the exact solution is

$\bar{\rho}\left(x\right)=-\frac{\delta\bar{u}\left(t\right)}{\bar{p}^{2}+\bar{p}_{,t}}~{},~{}\bar{u}=u_{0}e^{-\left(\delta+1\right)\int\bar{p}dt},$

where

$\int^{p\left(t\right)}\left(e^{W\left(-\frac{s}{e^{p_{0}}}+p_{0}\right)}-s^{2}\right)ds=t-t_{0}$
(124)

in which $W\left(t\right)$ is the Lambert function.

##  4 Conclusions

In this work we considered a drift-flux two-phase fluids model without
gravitational and wall friction forces where the equation of state parameter
for the fluids is that of a polytropic gas. The system of three hyperbolic
differential equations studied with the use of Lie’s theory. In particular we
investigate the algebraic properties of the two-phase fluids model by
calculate the one-parameter point transformations in which the dynamical
system is invariant.

We found that in the general scenario where the polytropic exponents
$\gamma_{1},~{}\gamma_{2}$ of the two fluids are arbitrary the admitted Lie
symmetries from a algebra of dimension fourth, while when in the special case
where the polytropic exponents are equal, that is, $\gamma_{1}=\gamma_{2}$ the
admitted Lie symmetries form a sixth dimensional Lie algebra. That result is
different from that previously found in the literature where it was found that
a sixth dimensional Lie algebra it is admitted only when
$\gamma_{1},~{}\gamma_{2}$ are equal with one, that is,
$\gamma_{1}=\gamma_{2}=1$ [41].

Thus, the missed symmetry vector by the previous study where
$\gamma_{1},~{}\gamma_{2}$ are equal it is important for the determination of
new similarity transformations. Indeed, because the symmetry vector is defined
only in the two-dimensional space of the dependent variables and it commutes
with the rest symmetries it provides a large number of independent similarity
transformations which lead to similarity solutions which can not connect
through an Adjoint transformation.

Furthermore, the one-dimensional optimal systems were determined for all the
cases which followed by the classification of the Lie symmetries, the
knowledge of the one-dimensional optimal system is essential because all the
unique similarity solutions can be classified. The results of this work
includes and new similarity solutions for the two-phase fluids model of our
study.

In this work we have not studied the initial value problem for the two-phase
fluids model. Indeed, physical problems are defined with a set of boundary and
initial conditions. Consequently, not all similarity solutions found in this
work will satisfy the initial value problems for all the physical states.
However, we were able to find all the possible similarity solutions, where
these can be constrained according to the initial conditions. Such analysis
extends the scopus of this work and will be investigated in a future study.

This study contributes on the subject of the study of the algebraic properties
of hyperbolic equations in fluid dynamics. From the result of this work it is
clear that Lie symmetries play an important role on the determination of exact
solutions in the two-phase fluids models. Although in this work we studied the
simplest drift-flux two-phase fluids model, in the future we plan to extend
our analysis in more general models.

Conflict of interest: This work does not have any conflicts of interest

Funding information: There are no funders to report for this submission

## References

  * [1] S. Lie, Theorie der Transformationsgruppen: Vol I, Chelsea, New York (1970)
  * [2] S. Lie, Theorie der Transformationsgruppen: Vol II, Chelsea, New York (1970)
  * [3] S. Lie, Theorie der Transformationsgruppen: Vol III, Chelsea, New York (1970)
  * [4] N.H. Ibragimov, CRC Handbook of Lie Group Analysis of Differential Equations, Volume I: Symmetries, Exact Solutions, and Conservation Laws, CRS Press LLC, Florida (2000)
  * [5] G.W. Bluman and S. Kumei, Symmetries of Differential Equations, Springer-Verlag, New York, (1989)
  * [6] H. Stephani, Differential Equations: Their Solutions Using Symmetry, Cambridge University Press, New York, (1989)
  * [7] P.J. Olver, Applications of Lie Groups to Differential Equations, Springer-Verlag, New York, (1993)
  * [8] W. Sarlet and M. Crampin, J. Phys. A: Math Gen. 18, L563 (1985)
  * [9] P.G.L. Leach, Applicable Analysis and Discrete Mathematics 6, 238 (2012)
  * [10] A. Paliathanasis and P.G.L. Leach, Int. J. Geom. Meth. Mod. Phys. 13, 1630009 (2016)
  * [11] B. Abraham-Shrauner, J. Math. Phys. 34, 4809 (1993)
  * [12] S. Anco and G. Bluman, J. Math. Phys. 37, 2361 (1996)
  * [13] G. Bluman, Math. Comput. Modelling, 8, 1 (1993)
  * [14] P.J. Olver, J.A. Sanders and J.P. Wang, J. Nonl. Math. Phys. 9, 164 (2002)
  * [15] G.M. Mubarakzyanov, On solvable Lie algebras Izvestia Vysshikh Uchebn Zavendeniĭ Matematika 32, 114 (1963)
  * [16] G.M. Mubarakzyanov,Classification of real structures of five-dimensional Lie algebras Izvestia Vysshikh Uchebn Zavendeniĭ Matematika 34, 99 (1963)
  * [17] G.M. Mubarakzyanov, Classification of solvable six-dimensional Lie algebras with one nilpotent base element Izvestia Vysshikh Uchebn Zavendeniĭ Matematika 35, 104 (1963)
  * [18] J. Patera, P. Winternitz and H. Zassenhaus J. Math. Phys. 16, 1597 (1975)
  * [19] J. Patera, R.T. Sharp and P. Winternitz, J. Math. Phys, 17, 986 (1976)
  * [20] J. Patera and P. Winternitz,, J. Math. Phys. 18, 1449 (1977)
  * [21] P.G.L. Leach, J. Aust. Math. Soc. B 20, 97 (1977)
  * [22] F.M. Mahomed, Math. Meth. Appl. Sci. 30, 1995 (2007)
  * [23] R. Cherniha, V. Davydovych and J.R. King, Symmetry 10, 171 (2018)
  * [24] M. Craddock, J. Diffential Equations 166, 107 (2000)
  * [25] A. Paliathanasis, Math. Meth. Appl. Sci. 43, 716 (2020)
  * [26] G.M. Webb, J. Phys. A: Math. Gen. 23, 3885 (1990)
  * [27] S. Kontogiorgis, R.O. Popovych and C. Sophocleous, Acta Applicandae Mathematicae 163, 91 (2019)
  * [28] S. Moyo and P.G.L. Leach, Proceedigns of Institute of Mathematics of NAS of Ukraine 50, 204 (2004)
  * [29] A. Paliathanasis, K. Krishnakumar, K.M. Tamizhmani and P.G.L. Leach, Mathematics 4, 28 (2016)
  * [30] A.A. Chesnokov, J. Appl. Mech. Techn. Phys. 49, 737 (2008)
  * [31] A.A. Chesnokov, Eur. J. Appl. Math. 20, 461 (2009)
  * [32] A. Paliathanasis, Symmety 11, 1115 (2019)
  * [33] A. Paliathanasis, International Journal of Nonlinear Sciences and Numerical Simulation, 20, 739 (2020)
  * [34] S.V. Meleshko and N.F. Samatova, AIP Conference Proceedings 2164, 050003 (2019)
  * [35] A. Bihlo, N. Poltavets and R.O. Popovych, Chaos 30, 073132 (2020)
  * [36] A.V. Aksenov and K.P. Druzhkov, J. Phys.: Conf. Ser. 722, 012001 (2016)
  * [37] V.A. Dorodnitsym and E.I. Kaptsov, Commun. Nonl. Sci. Num. Sim. 89, 105343 (2020)
  * [38] J.-G. Liu, Z.-F. Zeng, Y. He and G.-P. Ai, Int. J. Nonl. Sci. Num. Sim. 16, 114 (2013)
  * [39] A. Hematulin, S.V. Meleshko and S.G. Gavrilyuk, Math. Methods Appl. Sci. 30 2101 (2007)
  * [40] E. Pucci, G. Saccomandiand R. Vitolo, Int. J. Eng. Sci. 99, 13 (2016)
  * [41] B. Bira and T.R. Sekhar, Appl. Math. Mech. -Engl. Ed. 36, 1105 (2015)
  * [42] B. Bira, T.R. Sekhar and D. Zeidan, Comput. Math. Appl. 71, 46 (2016)
  * [43] G.F. Hewitt, Heat Transfer Enginnering 4, 67 (1983)
  * [44] D. Zeidan, P. Bähr, P. Farber, J. Gräbel and P. Ueberholz, Computer and Fluids 181, 90 (2019)
  * [45] E. Goncalves and D. Zeidan, Commun. Comput. Phys. 24, 167 (2018)
  * [46] J. Zhai, W. Liu and L. Yuan, Computers and Fluids 134, 90 (2016)
  * [47] D. Zeidan, L.T. Zhang and E. Goncalves, Int. J. Appl. Mech. 12, 2050049 (2020)
  * [48] D. Zeidan and B. Bira, Math. Meth. Appl. Sci. 42, 4679 (2019)
  * [49] A. Kozakevicius, D. Zeidan, A.A. Schmidt and S. Jakobsson, International Journal of Numerical Methods for Heat and Fluid Flow 28, 2052 (2018)
  * [50] E. Goncalves and D. Zeidan, International Journal of Numerical Methods for Heat and Fluid Flow, 27, 1487 (2017)
  * [51] S. Kuila, T.R. Sekhar and D. Zeidan, International Journal of Computational Methods 13, 1650009 (2016)
  * [52] D. Zeidan, E. Romenski, A. Slaouti and E.F. Toro, International Journal for Numerical Methods in Fluids, 54, 393 (2007)
  * [53] M.K. Banda, M. Herty, J.-M. T. Ngnotghouye, SIAM J. Sci. Comput. 31, 4633 (2010)
  * [54] L. Pan, S.W. Webb and C.M. Oldenburg, Adv. Water Resources 34, 1656 (2011)
  * [55] G.C. Sander, J.-Y. Paralange, I.G. Lisle and S.W. Weeks, Adv. Water Resources 28, 1112 (2005)
  * [56] M.D. Thanh, Nonlinear Analysis: Real World Appications 13, 987 (2012)
  * [57] Y. Taitel and D. Barena, Encyclopedia of Two-Phase Heat Transfer and Flow I, Fundamentals and Methods Volume 1 :Modeling of Gas Liquid Flow in Pipes, edited by J.R. Thome, World Scientiffic Publishing Co. Pte. Ltd., Singapure, (2015)