# On explicit solutions for coupled reaction-diffusion and Burgers-type
equations with variable coefficients through a Riccati System

José M. Escorcia  Escuela de Ciencias Aplicadas e Ingeniería, Universidad
EAFIT, Carrera 49 No. 7 Sur-50, Medellín 050022, Colombia.
<EMAIL_ADDRESS>and  E. Suazo  School of Mathematical and Statistical
Sciences, University of Texas Rio Grande Valley, 1201 W. University Drive,
Edinburg, Texas, 78539-2999<EMAIL_ADDRESS>
(Date: August 28, 2024)

###### Abstract.

This work is concerned with the study of explicit solutions for generalized
coupled reaction-diffusion and Burgers-type systems with variable
coefficients. Including nonlinear models with variable coefficients such as
diffusive Lotka-Volterra model, the Gray-Scott model, the Burgers equations.
The equations’ integrability (via the explicit formulation of the solutions)
is accomplished by using similarity transformations and requiring that the
coefficients fulfill a Riccati system. We present traveling wave type
solutions, as well as solutions with more complex dynamics and relevant
features such as bending. A Mathematica file has been prepared as
supplementary material, verifying the Riccati systems used in the construction
of the solutions.

Keywords: Coupled reaction-diffusion equations, Coupled Burgers equations,
Traveling wave solution, Similarity transformations, and Riccati system.

###### 1991 Mathematics Subject Classification:

Primary 81Q05, 35C05. Secondary 42A38

##  1\. Introduction

A large number of problems in engineering and applied sciences are modeled by
partial differential equations (PDEs) or systems of PDEs. In this regard, the
study of such PDEs or systems of PDEs is a constantly evolving area due to
their importance in several fields of knowledge. One type of differential
equations that has attracted the attention of researchers in recent years are
those of reaction-diffusion. In biology, the reaction-diffusion systems
explain the pattern of stripes of leopards, snakes, zebras, and jaguars [31,
53, 65] and creates a model for replicating the hepatitis B virus [20, 36,
67]. Reaction-diffusion equations can also be used to study tumor growth [8,
22, 69] and fission wave dynamics [42], among other topics.

From a theoretical standpoint, aspects of reaction-diffusion models, such as
well-posedness and asymptotic behavior of the solutions, have also been
thoroughly explored [18, 35, 39, 56, 57], allowing for significant progress in
our understanding of these systems. A hot topic in the study of PDEs and with
great impact on both theoretical and applied domains is the development of
techniques that allow exact solutions [47, 48, 50] . In the case of reaction-
diffusion equations, there is a wealth of literature describing various
strategies; a few instances include the Hirota method [45, 64], similarity
transformations [46, 60], and the tanh-method [33, 51] among others [3, 5, 47,
49], that have enabled the construction of explicit solutions.

The aforementioned methods have been successfully applied to classical models
such as the linear coupled reaction-diffusion equation [14, 55], the diffusive
Lotka-Volterra system [9, 10, 63], the Gray-Scott model [21, 26, 52], and the
coupled Burgers system [1, 37, 68] (Section 2 provides additional information
on these systems). Essentially, the type of solutions most reported in the
literature for these models are the so-called traveling waves, due to their
fascinating mathematical and physical properties. However, there is still a
lot of interest in finding solutions with more complicated structures
(compared to solitary waves) that allow us to describe phenomena that the
former cannot explain. The main goal of this research is to introduce some
generalizations of these models as well as a mechanism for finding explicit
solutions.

To be more specific, we primarily study the explicit solutions for the general
coupled equations of the type

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
H(\psi)+F_{1}(\psi,\varphi),$ (1.1) $\displaystyle\varphi_{t}$
$\displaystyle=$ $\displaystyle H(\varphi)+F_{2}(\psi,\varphi),$ (1.2)

where the variable operator $H$ has a quadratic form

$H(\psi):=a(t)\psi_{xx}-b(t)x^{2}\psi+c(t)x\psi_{x}+d(t)\psi+xf(t)\psi-g(t)\psi_{x}$

and $F_{i}(\psi,\varphi)$ ($i=1,2$) will be given. Here, the coefficients
$a(t),$ $b(t),$ $c(t),$ $d(t),$ $f(t),$ $g(t)$ are time-dependent real
functions. By selecting appropriate $F_{i}$, the system (1.1)-(1.2) can be
viewed as generalizations of the classical reaction-diffusion models and the
Burgers system discussed above. The integrability of equations (1.1)-(1.2) is
accomplished by requiring the coefficients to fulfill a Riccati system and
extending, some classical model solutions via similarity transformations.
Although the solutions discovered are primarily derived from traveling waves,
we emphasize that our findings are considerably more general and applicable to
any sort of solution. However, the employment of traveling waves did not
constitute a limitation since it was clear from these examples that the system
(1.1)-(1.2) allowed solutions of the same kind, as well as others with more
fascinating dynamics.

The Riccati system (3.10)-(3.17) has been applied to the analysis of several
PDEs with variable coefficients, including the nonlinear Schrödinger equation
[2, 4, 12, 15, 59, 61], diffusion-type equations [60], and reaction-diffusion
equations [46]. In our recent paper [16], we use the Riccati system to
investigate the existence (by explicit construction) of Rogue wave and dark-
bright soliton-like solutions for a coupled system of Schrödinger equations
(generalizing the Manakov system). And now in the present work, we extend the
ideas of [46] by investigating a coupled reaction-diffusion system and coupled
Burgers equations with the structure indicated in (1.1)-(1.2).

After the system (3.10)-(3.17) was introduced in [12] it has brought attention
of different researchers. In [7] it was considered a new example of a Lie
system that is not a Hamiltonian system. It is also one of the few Lie systems
connected to PDEs until now [7]. System (3.10)-(3.17) describes integral
curves of the t-dependent vector field of the form

$X_{t}=a(t)X_{1}+b(t)X_{2}+c(t)X_{3}-2d(t)X_{4}+f(t)X_{5}+g(t)X_{6}$

where $X_{1},...,X_{6}$ satisfy certain commutations rules showing that
(3.10)-(3.17) is a Lie system associated to a Vessiot Guldberg Lie algebra of
Hamiltonian vector fields with respect to a presymplectic form, see [7].
Systems such as (3.10)-(3.17) originated the introduction of a particular
class of Lie systems on Dirac manifolds, denominated Dirac–Lie systems, see
[7]. It is part of the objectives of this work show more applications of this
Lie system in the construction of transformations for nonlinear
reaction–diffusion equations.

The consideration of parameters in this work is inspired by the work of
Marhic, who, in 1978 [34] introduced (probably for the first time) a one-
parameter $\\{\alpha(0)\\}$ family of solutions for the linear Schrödinger
equation of the one-dimensional harmonic oscillator. The solutions presented
by Marhic constituted a generalization of the original Schrödinger wave packet
with oscillating width.

We have also introduced solutions for a modified Riccati system (3.34)-(3.40),
see the appendix. We have prepared a Mathematica file as supplementary
material verifying the Riccati systems used in this work.

This work makes the following contributions:

  * •

The introduction of _new coupled reaction-diffusion and Burgers-type equations
with variable coefficients_ and, by the use of similarity transformations and
Riccati system, we provide a mechanism to find explicit solutions. These new
models generalize some classical systems such as the linear coupled reaction-
diffusion equation, the diffusive Lotka-Volterra system, the Gray-Scott model,
and the coupled Burgers system.

  * •

The report of _new solutions to coupled reaction-diffusion and Burgers-type
equations_ , include traveling wave solutions as well as solutions with more
advanced kinetics determined by bending properties.

  * •

The introduction of a _modified Riccati system_ that allows us to generate
exponential-type solutions for the general linear reaction-diffusion system.
The solution for this Riccati system is provided. All solutions of the Riccati
systems utilized in the current research were verified using the computer
algebra system Mathematica.

The rest of the paper is organized as follows: In Section 2, we present a
brief overview of classical models, such the linear coupled reaction-diffusion
equation, the diffusive Lotka-Volterra system, the Gray-Scott model, and the
coupled Burgers system. Some of the solutions found in the literature are also
presented in this section. Section 3 is the core of the current research.
Generalizations of the systems discussed in the preceding section are
presented here, as well as sufficient conditions for the construction of exact
solutions in terms of the solutions of the constant coefficient models (see
Theorems 1-5 and Proposition 1). In particular, Theorem 2 employs the modified
Riccati system. As a result of the coefficient balance requirements (Riccati
system), traveling wave type solutions and solutions with noticeable bending
dynamics are shown. A Mathematica file has been prepared as supplementary
material verifying the Riccati systems used in the construction of the
solutions in this section. Conclusions and some final remarks are given in
Section 4. Section 5 corresponds to the appendix, in which the solutions of
the Riccati systems are presented.

##  2\. Classical Coupled Reaction-Diffusion and Burgers Systems

In this section, we briefly describe some well-known reaction-diffusion
equations, such as: the linear reaction-diffusion model, the diffusive Lotka-
Volterra, and the Gray-Scott model. Likewise, the coupled Burgers equations
are presented. We mainly exhibit the traveling wave solutions of these
equations. However, the findings of the paper are valid for all solutions of
these classical systems.

###  2.1. The linear reaction-diffusion system

Several processes during the embryonic state are associated with the migration
and proliferation of cells within growing tissues. A canonical model for these
processes is the linear reaction-diffusion system [14, 55]

$u_{\tau}=a_{1}u_{\xi\xi}-b_{1}u+v,$ (2.1) $v_{\tau}=a_{1}v_{\xi\xi}-b_{2}v,$
(2.2)

where the coefficients $a_{1},\ b_{i}$ $(i=1,2)$ are nonnegative. Another
application is the study of the motion and biodegradation of dissolved organic
contaminants in a saturated porous medium [54]. The following explicit
solutions are reported in the literature [11, 38]:

$u(\xi,\tau)=\left[e^{-(b_{1}+a_{1})\tau}+e^{-(b_{2}+a_{1})\tau}\right]\cos\xi,$
(2.3) $v(\xi,\tau)=(b_{1}-b_{2})e^{-(b_{2}+a_{1})\tau}\cos\xi.$ (2.4)

A family of exact solutions for (2.1)-(2.2) were constructed in the framework
of a growing domain (time-dependence domain) [55]. But, for the papers
purpose, we will only consider the solutions described above.

###  2.2. The diffusive Lotka-Volterra system

The diffusive Lotka-Volterra (DLV) system is a generalization of the classical
Lotka-Volterra system [32, 66], in which the space diffusion is taking into
account in order to have a more realistic population dynamics description.
This system is used as model in different processes in biology [6], ecology
[41], and medicine [28].

The two-component DLV system

$u_{\tau}=u_{\xi\xi}+u(a_{1}-b_{1}u-c_{1}v),$ (2.5)
$v_{\tau}=v_{\xi\xi}+v(a_{2}-b_{2}u-c_{2}v)$ (2.6)

with $a_{i},b_{i},c_{i}\geq 0$ $(i=1,2)$ describes competition between two
populations $u$ and $v$, in which diffusivity is the same for both
populations. This DLV system admits the traveling wave solutions [9, 10, 63]:

$u(\xi,\tau)=\frac{A}{4B}\left[1-\tanh\left(\frac{\sqrt{A}}{24}\xi-\frac{5A}{12}\tau\right)\right]^{2},$
(2.7)
$v(\xi,\tau)=\nu_{0}+\nu_{1}\frac{A}{4B}\left[1-\tanh\left(\frac{\sqrt{A}}{24}\xi-\frac{5A}{12}\tau\right)\right]^{2}.$
(2.8)

Here, the parameters $A,\ B,\ \nu_{0},\ \nu_{1}$, are defined by the
relations:

$A=\left\\{\begin{array}[]{lcc}a_{1}=a_{2},&if&\nu_{0}=0,\\\ \\\
a_{1}-a_{2}\frac{c_{1}}{c_{2}},&if&\nu_{0}=\frac{a_{2}}{c_{2}}\end{array}\right.,\quad
B=\left\\{\begin{array}[]{lcc}\frac{c_{1}b_{2}-b_{1}c_{2}}{c_{1}-c_{2}},&if&\nu_{0}=0,\\\
\\\ b_{1}+c_{1}\nu_{1},&if&\nu_{0}=\frac{a_{2}}{c_{2}}\end{array}\right.,$

$\nu_{1}=\left\\{\begin{array}[]{lcc}\frac{b_{1}-b_{2}}{c_{2}-c_{1}},&if&c_{1}\not=c_{2},\
b_{1}\not=b_{2},\\\ \\\ -\frac{a_{2}b_{1}}{a_{1}c_{1}},&if&c_{1}=c_{2},\
b_{1}=b_{2}.\end{array}\right.$

It is well known that these solutions possess the following properties [9, 10,
63]: If $\nu_{0}\not=0$, the solutions have the asymptotical behavior

$(u,v)\rightarrow\left(\frac{a_{1}}{b_{1}},0\right),\quad\quad\mbox{as}\quad\quad\tau\rightarrow\infty,$
(2.9)

provided the condition $\hat{A}>\max\\{\hat{B},\hat{C}\\}$ where
$\hat{A}=\frac{a_{1}}{a_{2}},$ $\hat{B}=\frac{b_{1}}{b_{2}},$ and
$\hat{C}=\frac{c_{1}}{c_{2}}$ is satisfied. In contrast, if $\nu_{0}=0$ (with
$a_{1}=a_{2}$), then the solutions converge to the stationary state

$(u,v)\rightarrow\left(\frac{a_{1}(\hat{C}-1)}{b_{2}(\hat{C}-\hat{B})},\frac{a_{1}(1-\hat{B})}{c_{2}(\hat{C}-\hat{B})}\right),\quad\quad\mbox{as}\quad\quad\tau\rightarrow\infty.$
(2.10)

The asymptotic behavior is valid if one of these relations is satisfied:
$\hat{B}>\hat{A}=1>\hat{C}$ or $\hat{C}>\hat{A}=1>\hat{B}$. In a biological
context, the first property can describe the competition of two species, in
which one of them dominates and the other eventually dies out. Behavior (2.10)
represents an arbitrary long coexistence of the species.

In order to describe a larger number of population interaction such as a
predator-predator-prey model and a predator-prey-competition model, the three-
component DLV system is commonly assumed in the literature. This DLV system
reads [9, 43, 44]:

$u_{\tau}=u_{\xi\xi}+u(a_{1}-b_{1}u-c_{1}v-e_{1}w),$ (2.11)
$v_{\tau}=v_{\xi\xi}+v(a_{2}-b_{2}u-c_{2}v-e_{2}w),$ (2.12)
$w_{\tau}=w_{\xi\xi}+w(a_{3}-b_{3}u-c_{3}v-e_{3}w),$ (2.13)

where $a_{i},b_{i},c_{i},e_{i},$ $(i=1,2,3)$ are constants. As in the previous
case, system (2.11)-(2.13) possesses solutions in a traveling wave form [9,
23]:

$u(\xi,\tau)=\left(2+\theta-\frac{a_{1}}{4}\right)\left[1-\tanh\left(\xi-\theta\tau\right)\right]^{2},$
(2.14)
$v(\xi,\tau)=\frac{a_{1}}{4}\left[1+\tanh\left(\xi-\theta\tau\right)\right]^{2},$
(2.15)
$w(\xi,\tau)=\left(a_{1}-\theta-2\right)\left[1-\tanh\left(\xi-\theta\tau\right)\right],$
(2.16)

as long as $a_{1}=a_{2}=a_{3},$ $b_{1}=1,\
b_{2}=-\frac{a_{1}-24}{8-a_{1}+4\theta},\
b_{3}=-\frac{a_{1}-4-2\theta}{8-a_{1}+4\theta},$ $c_{1}=-\frac{4\theta-
a_{1}-16}{a_{1}},\ c_{2}=1,\ c_{3}=-\frac{2\theta-a_{1}-4}{a_{1}},\
e_{1}=-\frac{a_{1}-4-2\theta}{2+\theta-a_{1}},\
e_{2}=-\frac{a_{1}-4+2\theta}{2+\theta-a_{1}},\ e_{3}=1,$ with
$\theta+2<a_{1}<4(\theta+2)$. The last inequality is used to guarantee the
positivity of the solutions.

###  2.3. The Gray-Scott model

The following is one the most important reaction-diffusion systems due to the
diversity of phenomena that can be described by the model. More specifically,
we consider the system of equations

$u_{\tau}=u_{\xi\xi}-uv^{2}+b_{1}(1-u),$ (2.17)
$v_{\tau}=v_{\xi\xi}+uv^{2}-b_{1}v,$ (2.18)

where $b_{1}$ is the diffusion constant [21, 26, 52]. Gray and Scott [21]
introduced this system in 1984 as a mathematical model for the cubic
autocatalytic reaction. One of the exact solutions reported in the literature
has the form [26, 52]

$u(\xi,\tau)=\frac{3-\sqrt{1-4b_{1}}}{4}-\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(\xi-\theta\tau)\right),$
(2.19)
$v(\xi,\tau)=\frac{1+\sqrt{1-4b_{1}}}{4}+\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(\xi-\theta\tau)\right),$
(2.20)

with $\theta=\frac{\sqrt{2}(1-3\sqrt{1-4b_{1}})}{4}.$ Travelling waves
solutions allow us to describe a large number of chemically reacting systems
[17, 30]. Recently in [25], numerical simulations of pattern formation for
more complex solutions revealed the potential applicability in 3D printing.

###  2.4. The coupled Burgers system

The coupled Burgers model is a system of partial differential equations as
follows:

$u_{\tau}=u_{\xi\xi}-b_{1}uu_{\xi}-c_{1}(uv)_{\xi},$ (2.21)
$v_{\tau}=v_{\xi\xi}-b_{2}vv_{\xi}-c_{2}(uv)_{\xi},$ (2.22)

with $b_{i},c_{i}$ $(i=1,2)$ as constants. The system (2.21)-(2.22) has been
studied extensively from theoretical and numerical point of view [1, 13, 29,
37, 68]. As in the previous models, the coupled Burgers equations admit a
travelling wave solution in the form [1, 29]:

$u(\xi,\tau)=B-2A(2c_{1}-1/4c_{1}c_{2}-1)\tanh\left[A\left(20\xi-10-2A\tau\right)\right],$
(2.23)
$v(\xi,\tau)=B(2c_{1}-1/4c_{1}c_{2}-1)\tanh\left[A\left(20\xi-10-2A\tau\right)\right],$
(2.24)

where $A=\left(4c_{1}c_{2}-1/4c_{1}-2\right).$ Multiple-kink solutions are
also reported for this system [24]. In the field of fluid dynamics, the
coupled Burgers equation plays a crucial role in the study of turbulence
because it can be considered as a simple case of the Navier–Stokes equations.

##  3\. The Generalized Coupled Reaction-Diffusion and Burgers-type Systems

This section demonstrates how to generate explicit solutions for general
variable coefficients systems. Each of these general models can be viewed as
extensions of the systems discussed in the preceding section. In our findings,
the integrability of such systems is obtained on the premise of an equilibrium
in the coefficients, which essentially assumes that they satisfy a Riccati
system; as a result, the solutions of models with variable coefficients, which
have more complicated dynamics, are obtained by transforming those solutions
of models corresponding to systems of constant coefficients.

###  3.1. Solutions for the generalized linear reaction-diffusion system

Let’s start with the first result of this section, which allows the
construction of explicit solutions of the general linear reaction-diffusion
model.

######  Theorem 1 (Generalized Linear Reaction-Diffusion System).

The variable coefficients reaction-diffusion system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{1}(t)-xf(t))\psi-(g(t)-c\left(t\right)x)\psi_{x}+h(t)\varphi,$
(3.1) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{2}(t)-xf(t))\varphi-(g(t)-c\left(t\right)x)\varphi_{x}$
(3.2)

can be transformed into the constant coefficients system

$u_{\tau}=u_{\xi\xi}-b_{1}u+v,$ (3.3) $v_{\tau}=v_{\xi\xi}-b_{2}v.$ (3.4)

###### Proof.

We will consider solutions for (3.1)-(3.2) of the form

$\psi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}u(\xi,\tau),\qquad\xi=\beta(t)x+\varepsilon(t),\qquad\tau=\gamma(t),$
(3.5)

and

$\varphi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}v(\xi,\tau).\qquad$
(3.6)

Let $S=\alpha(t)x^{2}+\delta(t)x+\kappa(t).$ Then, the derivatives of the
function $\psi$ (similarly, derivatives of $\phi$ are obtained by changing the
function $u$ by $v$) are

$\psi_{x}=\mu^{-1/2}e^{S}\left[u(2\alpha x+\delta)+\beta u_{\xi}\right],$
(3.7) $\psi_{xx}=\mu^{-1/2}e^{S}\left[u(2\alpha x+\delta)^{2}+2\beta
u_{\xi}(2\alpha x+\delta)+u(2\alpha x+\delta)+\beta^{2}u_{\xi\xi}\right],$
(3.8)
$\psi_{t}=\mu^{-1/2}e^{S}\left[-\frac{1}{2}\mu^{-1}\mu^{\prime}u+u(\alpha^{\prime}x^{2}+\delta^{\prime}x+\kappa^{\prime})+\gamma^{\prime}u_{\tau}+u_{\xi}(\beta^{\prime}x+\varepsilon^{\prime})\right].$
(3.9)

After substituting (3.7)-(3.9) into (3.1)-(3.2), we obtain the Riccati system
[12, 15, 61, 59, 60]:

$\dfrac{d\alpha}{dt}+b(t)=2c(t)\alpha+4a(t)\alpha^{2},$ (3.10)
$\dfrac{d\beta}{dt}=(c(t)+4a(t)\alpha(t))\beta,$ (3.11)
$\dfrac{d\gamma}{dt}=a(t)\beta^{2}(t),$ (3.12)
$\dfrac{d\delta}{dt}+2\alpha(t)g(t)=(c(t)+4a(t)\alpha(t))\delta+f(t),$ (3.13)
$\dfrac{d\varepsilon}{dt}=(2a(t)\delta(t)-g(t))\beta(t),$ (3.14)
$\dfrac{d\kappa}{dt}=a(t)\delta^{2}(t)-g(t)\delta(t).$ (3.15)

Considering the substitution

$\alpha=-\dfrac{\mu^{\prime}(t)}{4a(t)\mu(t)}-\dfrac{d(t)}{2a(t)},$ (3.16)

it follows that the Riccati equation (3.10) becomes (characteristic equation)

$\mu^{\prime\prime}-\eta(t)\mu^{\prime}-4\sigma(t)\mu=0,$ (3.17)

with

$\eta(t)=\frac{a^{\prime}}{a}+2c-4d,\hskip
28.45274pt\sigma(t)=ab+cd-d^{2}+\frac{d}{2}\left(\frac{a^{\prime}}{a}-\frac{d^{\prime}}{d}\right).$
(3.18)

Further, if we choose the conditions

$h(t)=a(t)\beta^{2}(t),\quad\quad L_{i}(t)=-b_{i}a(t)\beta^{2}(t),\quad
i=1,2,$

the functions $u(\xi,\tau)$ and $v(\xi,\tau)$ will satisfy the system
(3.3)-(3.4). The Riccati system (3.10)-(3.15) is solved explicitly in terms of
the fundamental solution of the characteristic equation (3.17), as can be
found in the Appendix of this work (see equations (5.1)-(5.1)). ∎

According to Theorem 1, the interaction coefficient $h(t)$ determines how the
dynamics of the solution $\psi$ rely on the function $\varphi$. Clearly, each
coefficient in the equation has a substantial impact on the dynamics of the
solutions, but we will only focus on the diffusion and interaction
coefficients in this work. The examples below demonstrate how the interaction
coefficient might influence the dynamics of solutions when it exhibits two
distinct types of decay.

####  3.1.1. Constant diffusivity and an interaction coefficient with
exponential decay

Consider the general reaction-diffusion system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle\psi_{xx}-\psi-
b_{1}e^{-2t}\psi-x\psi_{x}+e^{-2t}\varphi,$ (3.19) $\displaystyle\varphi_{t}$
$\displaystyle=$ $\displaystyle\varphi_{xx}-\varphi-
b_{2}e^{-2t}\varphi-x\varphi_{x}.$ (3.20)

By solving the Riccati system (3.10)-(3.15) with the initial conditions
$\alpha(0)=\delta(0)=\varepsilon(0)=\kappa(0)=\gamma(0)=0,\ \beta(0)=1,\
\mu(0)=1,$ we get

$\displaystyle\alpha(t)=\delta(t)=\varepsilon(t)=\kappa(t)=0,\quad\beta(t)=e^{-t},\quad\gamma(t)=e^{-t}\sinh
t,\quad\mu(t)=e^{2t}.$

Then, in line with Theorem 1, the system (3.19)-(3.20) admits a solution given
by the equations (3.5)-(3.6) with the functions $u,v$ satisfying the relations
(2.3)-(2.4) (with $a_{1}=1$), i.e.,

$\psi(x,t)=e^{-t}\cos\left(e^{-t}x\right)\left[e^{-(b_{1}+1)e^{-t}\sinh
t}+e^{-(b_{2}+1)e^{-t}\sinh t}\right],$ (3.21)
$\varphi(x,t)=(b_{1}-b_{2})e^{-t}\cos\left(e^{-t}x\right)e^{-(b_{2}+1)e^{-t}\sinh
t}.$ (3.22)

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 1. Solutions for the system (3.19)-(3.20) for the parameters $a_{1}=1,$
$b_{1}=100,$ and $b_{2}=1$. Here, (a) and (c) describe the profiles of the
functions $\psi$ and $\varphi$ respectively. The corresponding contours of
$\psi$ and $\varphi$ are shown in (b) and (d).

The time-evolution of these solutions for $b_{1}=100$ and $b_{2}=1$, are shown
in Figure 1. Therefore we can observe the expected exponential decay of the
solutions.

####  3.1.2. Periodic diffusivity and an interaction coefficient with rational
decay

In this case, we take into account the following system of equations:

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle e^{-2\cos
t}\psi_{xx}+(e^{2\cos
t}x^{2}+1)\psi-\frac{b_{1}e^{2}\psi}{\left[e^{2}+(16-2e^{2})t\right]^{2}}+(2-\sin
t)x\psi_{x}+\frac{e^{2}\varphi}{\left[e^{2}+(16-2e^{2})t\right]^{2}},$ (3.23)
$\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle e^{-2\cos
t}\varphi_{xx}+(e^{2\cos
t}x^{2}+1)\varphi-\frac{b_{2}e^{2}\varphi}{\left[e^{2}+(16-2e^{2})t\right]^{2}}+(2-\sin
t)x\varphi_{x}.$ (3.24)

The solution of the associated Riccati system is

$\displaystyle\alpha(t)=\frac{e^{2\cos
t}\left[4+(8-e^{2})t\right]}{t(-16+2e^{2})-e^{2}},\quad\beta(t)$
$\displaystyle=$ $\displaystyle\frac{e^{1+\cos
t}}{t(16-2e^{2})+e^{2}},\quad\gamma(t)=\frac{t}{t(16-2e^{2})+e^{2}},$

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 2. Solutions for the system (3.23)-(3.24) for the parameters $a_{1}=1,$
$b_{1}=25,$ and $b_{2}=1$. Here, (a) and (c) describe the profiles of the
functions $\psi$ and $\varphi$ respectively. The corresponding contours of
$\psi$ and $\varphi$ are shown in (b) and (d).

and

$\delta(t)\ =\ \varepsilon(t)=\kappa(t)=0,\quad\quad\mu(t)=1+t(16e^{-2}-2).$

Here, we have assumed the initial conditions $\alpha(0)=-4,\ \beta(0)=1,\
\mu(0)=1,\ \delta(0)=\varepsilon(0)=\kappa(0)=\gamma(0)=0.$ Therefore, the
system (3.23)-(3.24) admits the explicit solutions

$\psi(x,t)=\frac{\cos\left(\frac{e^{1+\cos
t}}{t(16-2e^{2})+e^{2}}x\right)}{\sqrt{1+t(16e^{-2}-2)}}\exp\left(\frac{e^{2\cos
t}\left[4+(8-e^{2})t\right]}{t(-16+2e^{2})-e^{2}}x^{2}\right)\left[e^{-(b_{1}+1)t/t(16-2e^{2})+e^{2}}+e^{-(b_{2}+1)t/t(16-2e^{2})+e^{2}}\right],$
(3.25) $\varphi(x,t)=\frac{(b_{1}-b_{2})\cos\left(\frac{e^{1+\cos
t}}{t(16-2e^{2})+e^{2}}x\right)}{\sqrt{1+t(16e^{-2}-2)}}e^{-(b_{2}+1)t/t(16-2e^{2})+e^{2}}.$
(3.26)

The solutions $\psi,\ \varphi$ are shown in Figure 2.

The next result allows the construction of exponential-type solutions for the
generalized linear reaction-diffusion equation. As we will see, the
construction of these solutions requires the introduction of a _modified
Riccati system_ which differs from (3.10)-(3.15) because of the equations for
the new functions $\kappa_{1}$ and $\kappa_{2}$.

######  Theorem 2 (Exponential-type Solutions for the Linear Reaction-
Diffusion System).

The variable coefficient coupled reaction-diffusion system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-(b\left(t\right)x^{2}-d(t)-xf(t))\psi-(g(t)-c\left(t\right)x)\psi_{x}+h(t)\varphi,$
(3.27) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-(b\left(t\right)x^{2}-d(t)-xf(t))\varphi-(g(t)-c\left(t\right)x)\varphi_{x}$
(3.28)

admits the explicit solutions

$\psi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\beta(t)xy+\gamma(t)y^{2}+\delta(t)x+\varepsilon(t)y+\kappa_{1}(t)}$
(3.29)

and

$\varphi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\beta(t)xy+\gamma(t)y^{2}+\delta(t)x+\varepsilon(t)y+\kappa_{1}(t)+\kappa_{2}(t)},$
(3.30)

with $y$ as a real parameter and the functions $\alpha(t),\ \beta(t),\
\gamma(t),\ \delta(t),\ \varepsilon(t),\ \kappa_{1}(t),\ \kappa_{2}(t),\
\mu(t),$ satisfying the modified Riccati system (3.34)-(3.40).

###### Proof.

Let
$S_{i}=\alpha(t)x^{2}+\beta(t)xy+\gamma(t)y^{2}+\delta(t)x+\varepsilon(t)y+\kappa_{1}(t)+(i-1)\kappa_{2}(t)$
for $i=1,2$. In these terms, the derivatives of the function $\psi$ are as
follows (similarly, derivatives of $\phi$ are obtained by changing $S_{1}$ for
$S_{2}$ and adding the term $\kappa_{2}^{\prime}$, in the temporal
derivative):

$\psi_{x}=\mu^{-1/2}e^{S_{1}}\left[2\alpha x+\beta y+\delta\right],$ (3.31)
$\psi_{xx}=\mu^{-1/2}e^{S_{1}}\left[(2\alpha x+\beta
y+\delta)^{2}+2\alpha\right],$ (3.32)
$\psi_{t}=\mu^{-1/2}e^{S_{1}}\left[-\frac{1}{2}\mu^{-1}\mu^{\prime}+\alpha^{\prime}x^{2}+\beta^{\prime}xy+\gamma^{\prime}y^{2}+\delta^{\prime}x+\varepsilon^{\prime}y+\kappa_{1}^{\prime}\right].$
(3.33)

Now, substituting (3.31)-(3.33) into the equations (3.27)-(3.28) we get the
_modified Riccati system_ :

$\dfrac{d\alpha}{dt}+b(t)=2c(t)\alpha+4a(t)\alpha^{2},$ (3.34)
$\dfrac{d\beta}{dt}=(c(t)+4a(t)\alpha(t))\beta,$ (3.35)
$\dfrac{d\gamma}{dt}=a(t)\beta^{2}(t),$ (3.36)
$\dfrac{d\delta}{dt}+2\alpha(t)g(t)=(c(t)+4a(t)\alpha(t))\delta+f(t),$ (3.37)
$\dfrac{d\varepsilon}{dt}=(2a(t)\delta(t)-g(t))\beta(t),$ (3.38)
$\dfrac{d\kappa_{1}}{dt}=a(t)\delta^{2}(t)-g(t)\delta(t)+h(t)e^{\kappa_{2}(t)},$
(3.39) $\dfrac{d\kappa_{2}}{dt}=-h(t)e^{\kappa_{2}(t)}.$ (3.40)

Now, as usual, the substitution

$\alpha=-\dfrac{\mu^{\prime}(t)}{4a(t)\mu(t)}-\dfrac{d(t)}{2a(t)}$ (3.41)

transforms the Riccati equation (3.34) into the characteristic equation
(3.17). In this sense, using equations (5.1)-(5.1) we can solve (3.34)-(3.38),
but equations (3.39)-(3.40) must be solved separately. In fact, the solution
of this Riccati system can be found in the Appendix of this work. ∎

Unlike Theorem 1, this last result enables us to create a family of
exponential-type solutions for the system (3.27)-(3.28) in which both
solutions, $\psi$ and $\varphi$, are directly dependent on the interaction
coefficient $h(t)$. Next, an example of Theorem 2 is shown.

####  3.1.3. Constant diffusivity and a periodic interaction coefficient

The following system with variable coefficients

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\psi_{xx}+2\left(x^{2}+\frac{1}{4}-3x\sin(3t)-x\cos(3t)\right)\psi+\left(2x-\cos(3t)\right)\psi_{x}+3\sin(6t)e^{\sin^{2}(3t)}\varphi,$
(3.42) $\displaystyle\varphi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varphi_{xx}+2\left(x^{2}+\frac{1}{4}-3x\sin(3t)-x\cos(3t)\right)\varphi+\left(2x-\cos(3t)\right)\varphi_{x}$
(3.43)

has the solutions

$\psi(x,t)=\exp\left[-\left(x-\cos(3t)\right)^{2}-\frac{t}{2}\right],$ (3.44)
$\varphi(x,t)=\exp\left[-\left(x-\cos(3t)\right)^{2}-\sin^{2}(3t)-\frac{t}{2}\right].$
(3.45)

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 3. Solutions for the system (3.42)-(3.43). The profiles of the
functions $\psi$ and $\varphi$ are shown in (a) and (c). In the contours of
$\psi$ and $\varphi$, Figures (b) and (d), the bending dynamics are clearly
observed.

The periodicity of the interaction coefficient causes an interesting dynamic:
The central axis of the solution bends left and right over time and eventually
disappears due to time exponential decay. The aforementioned dynamics can be
seen in Figure 3.

###  3.2. Solutions for the generalized Lotka-Volterra system

In the present section, we introduce a general diffusive Lotka-Volterra system
(two and three components) and provide a mechanism for producing solutions by
means of the classical models (2.5)-(2.6) and (2.11)-(2.13).

######  Theorem 3 (Generalized Diffusive Lotka-Volterra System).

The variable coefficient coupled reaction-diffusion system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{1}(t)-xf(t))\psi-(g(t)-c\left(t\right)x)\psi_{x}+(h_{1}(x,t)\psi+r_{1}(x,t)\varphi)\psi,$
(3.46) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{2}(t)-xf(t))\varphi-(g(t)-c\left(t\right)x)\varphi_{x}+(h_{2}(x,t)\psi+r_{2}(x,t)\varphi)\varphi,$
(3.47)

can be transformed into the classical Lotka-Volterra system

$u_{\tau}=u_{\xi\xi}+u(a_{1}-b_{1}u-c_{1}v),$ (3.48)
$v_{\tau}=v_{\xi\xi}+v(a_{2}-b_{2}u-c_{2}v).$ (3.49)

###### Proof.

We will assume solutions for (3.46)-(3.47) in the form

$\psi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}u(\xi,\tau),\qquad\xi=\beta(t)x+\varepsilon(t),\qquad\tau=\gamma(t),$
(3.50)

and

$\varphi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}v(\xi,\tau).\qquad$
(3.51)

After substituting (3.50)-(3.51) in (3.46)-(3.47), we obtain again the Riccati
system (3.10)-(3.15), i.e.,

$\dfrac{d\alpha}{dt}+b(t)=2c(t)\alpha+4a(t)\alpha^{2},$ (3.52)
$\dfrac{d\beta}{dt}=(c(t)+4a(t)\alpha(t))\beta,$ (3.53)
$\dfrac{d\gamma}{dt}=a(t)\beta^{2}(t),$ (3.54)
$\dfrac{d\delta}{dt}+2\alpha(t)g(t)=(c(t)+4a(t)\alpha(t))\delta+f(t),$ (3.55)
$\dfrac{d\varepsilon}{dt}=(2a(t)\delta(t)-g(t))\beta(t),$ (3.56)
$\dfrac{d\kappa}{dt}=a(t)\delta^{2}(t)-g(t)\delta(t).$ (3.57)

Here, the characteristic equation (3.17) can be obtained by the usual
substitution (3.41). Additionally, if we assume the conditions

$L_{i}(t)=a_{i}a(t)\beta^{2}(t),\quad\quad
h_{i}(x,t)=-b_{i}a(t)\beta^{2}(t)\mu^{1/2}(t)e^{-(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},$
$r_{i}(x,t)=-c_{i}a(t)\beta^{2}(t)\mu^{1/2}(t)e^{-(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},\quad
i=1,2,$

the functions $u(\xi,\tau)$ and $v(\xi,\tau)$ will satisfy the system
(3.48)-(3.49). ∎

In simple terms, the integrability of the variable coefficients Lotka-Volterra
system is consequences of the Riccati system and the exponential structure of
the interaction coefficients $h_{i}(x,t)$ and $r_{i}(x,t)$. The subsequent
examples illustrate the preceding result.

####  3.2.1. Constant diffusivity and tanh-type interaction coefficients

In this example we demonstrate the existence of travelling wave-type solutions
for the general model. Indeed, consider the system

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\psi_{xx}+\frac{1}{2}\left(\tanh
t+2a_{1}-1\right)\psi-\left(\frac{b_{1}e^{-\kappa(0)}}{2\sqrt{1+\tanh
t}}\psi+\frac{c_{1}e^{-\kappa(0)}}{2\sqrt{1+\tanh t}}\varphi\right)\psi,$
(3.58) $\displaystyle\varphi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varphi_{xx}+\frac{1}{2}\left(\tanh
t+2a_{2}-1\right)\varphi-\left(\frac{b_{2}e^{-\kappa(0)}}{2\sqrt{1+\tanh
t}}\psi+\frac{c_{2}e^{-\kappa(0)}}{2\sqrt{1+\tanh t}}\varphi\right)\varphi.$
(3.59)

The solution of the Riccati system associated with this problem is

$\displaystyle\alpha(t)=\delta(t)=0,\quad\beta(t)=\sqrt{2},\quad\gamma(t)=t,\quad\varepsilon(t)=\varepsilon(0),\quad\kappa(t)=\kappa(0),\quad\mu(t)=1+\tanh
t.$

Then, a solution for the system (3.58)-(3.59) has the following traveling wave
form:

$\psi(x,t)=\frac{Ae^{\kappa(0)}}{4B\sqrt{1+\tanh
t}}\left[1-\tanh\left(\frac{\sqrt{A}}{24}(\sqrt{2}x+\varepsilon(0))-\frac{5A}{12}t\right)\right]^{2},$
(3.60) $\varphi(x,t)=\frac{e^{\kappa(0)}}{\sqrt{1+\tanh
t}}\left\\{\nu_{0}+\nu_{1}\frac{A}{4B}\left[1-\tanh\left(\frac{\sqrt{A}}{24}(\sqrt{2}x+\varepsilon(0))-\frac{5A}{12}t\right)\right]^{2}\right\\}.$
(3.61)

These solutions exhibit the asymptotic behaviors:

$\left(\psi,\varphi\right)\rightarrow\left(\frac{a_{1}e^{\kappa(0)}}{\sqrt{2}b_{1}},0\right),\quad\mbox{as}\quad
t\rightarrow\infty,$ (3.62)

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 4. Figures (a) and (c) describe the traveling wave solutions for the
system (3.58)-(3.59) with $a_{1}=c_{2}=c_{1}=2,$ $a_{2}=1,$
$b_{1}=b_{2}=\sqrt{2}$, $\nu_{0}=\frac{1}{2},$ $\nu_{1}=-\frac{\sqrt{2}}{4}$,
$A=1,$ $B=\frac{\sqrt{2}}{2}$, $\kappa(0)=1$, and $\varepsilon(0)=2$. The
contours of these solutions are shown in (b) and (d).

or

$\left(\psi,\varphi\right)\rightarrow\left(\frac{a_{1}(\hat{C}-1)e^{\kappa(0)}}{\sqrt{2}b_{2}(\hat{C}-\hat{B})},\frac{a_{1}(1-\hat{B})e^{\kappa(0)}}{\sqrt{2}c_{2}(\hat{C}-\hat{B})}\right),\quad\mbox{as}\quad
t\rightarrow\infty.$ (3.63)

The constants $\hat{B},\hat{C},$ are the same as those established in Section
2. The profiles of these solutions are displayed in Figure 4. A _plausible
interpretation_ of these solutions would be to depict two types of scenarios
in a competition between two species. In the first scenario, one species
dominates the other until it is extinct, whereas in the second situation, both
species coexist over time. This is clear evidence that some aspects of the
classical Lotka-Volterra equations are applicable to the general system.

####  3.2.2. Constant diffusivity and interaction coefficients with periodic
central axes

Let’s assume the system of variable coefficients equations

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\psi_{xx}+\left(2+\frac{1}{2}a_{1}e^{-4t}\right)\psi+\left(2x-\sin
t-\frac{1}{2}\cos
t\right)\psi_{x}-\frac{1}{2}e^{-4t}\exp\left[2\left(x-\frac{1}{2}\sin
t\right)^{2}\right]\left(b_{1}\psi+c_{1}\varphi\right)\psi,$ (3.64)
$\displaystyle\varphi_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varphi_{xx}+\left(2+\frac{1}{2}a_{2}e^{-4t}\right)\varphi+\left(2x-\sin
t-\frac{1}{2}\cos
t\right)\varphi_{x}-\frac{1}{2}e^{-4t}\exp\left[2\left(x-\frac{1}{2}\sin
t\right)^{2}\right]\left(b_{2}\psi+c_{2}\varphi\right)\varphi.$ (3.65)

In this case, we have the functions

$\displaystyle\alpha(t)=-2,\quad\quad\beta(t)=e^{-2t},\quad\quad\gamma(t)=\frac{1}{8}-\frac{e^{-4t}}{8},\quad\quad\delta(t)=2\sin
t,$ $\displaystyle\ \varepsilon(t)=-\frac{1}{2}e^{-2t}\sin
t,\quad\quad\kappa(t)=-\frac{1}{2}\sin^{2}t,\quad\quad\mu(t)=1,$

and the system admits a solution given by the following equations

$\psi(x,t)=\exp\left[-2\left(x-\frac{1}{2}\sin
t\right)^{2}\right]\left[1-\tanh\left(\frac{\sqrt{A}}{24}e^{-2t}(x-\frac{1}{2}\sin
t)-\frac{5A}{96}(1-e^{-4t})\right)\right]^{2},$ (3.66)
$\varphi(x,t)=\exp\left[-2\left(x-\frac{1}{2}\sin
t\right)^{2}\right]\left\\{\nu_{0}+\nu_{1}\frac{A}{4B}\left[1-\tanh\left(\frac{\sqrt{A}}{24}e^{-2t}(x-\frac{1}{2}\sin
t)-\frac{5A}{96}(1-e^{-4t})\right)\right]^{2}\right\\}.$ (3.67)

As can be seen in Figure 5, solutions share an interesting bending dynamic
controlled by the function $\sin t$. Due to the similar behavior of these
solutions and with the intention of providing the reader with visual proof of
the difference between them, we report in Figure 5(d) the profile of
$|\psi-\varphi|$.

The next proposition establishes a natural generalization of Theorem 3,
corresponding to the three-component Lotka-Volterra equations.

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Profile of the difference $|\varphi-\psi|$.

Figure 5. Solutions for the system (3.64)-(3.65) for the parameters
$a_{1}=c_{2}=c_{1}=2,$ $a_{2}=1,$ $b_{1}=b_{2}=\sqrt{2}$,
$\nu_{0}=\frac{1}{2},$ $\nu_{1}=-\frac{\sqrt{2}}{4}$, $A=1,$
$B=\frac{\sqrt{2}}{2}$.

######  Proposition 1 (Generalized Three-Component Lotka-Volterra System).

The general three-component Lotka-Volterra equations

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{1}(t)-xf(t))\psi-(g(t)-c\left(t\right)x)\psi_{x}+(h_{1}\psi+r_{1}\varphi+s_{1}\phi)\psi,$
(3.68) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{2}(t)-xf(t))\varphi-(g(t)-c\left(t\right)x)\varphi_{x}+(h_{2}\psi+r_{2}\varphi+s_{2}\phi)\varphi,$
(3.69) $\displaystyle\phi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\phi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{3}(t)-xf(t))\phi-(g(t)-c\left(t\right)x)\phi_{x}+(h_{3}\psi+r_{3}\varphi+s_{3}\phi)\phi$
(3.70)

can be transformed into the system

$u_{\tau}=u_{\xi\xi}+u(a_{1}-b_{1}u-c_{1}v-e_{1}w),$ (3.71)
$v_{\tau}=v_{\xi\xi}+v(a_{2}-b_{2}u-c_{2}v-e_{2}w).$ (3.72)
$w_{\tau}=w_{\xi\xi}+w(a_{3}-b_{3}u-c_{3}v-e_{3}w).$ (3.73)

###### Proof.

By extending the same ideas from Theorem 3, we use the substitutions

$\psi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}u(\xi,\tau),\qquad\xi=\beta(t)x+\varepsilon(t),\qquad\tau=\gamma(t)$
(3.74)

and

$\varphi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}v(\xi,\tau),\qquad\quad\phi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}w(\xi,\tau).$
(3.75)

As a consequence of imposing the integrability conditions

$r_{i}(x,t)=-c_{i}a(t)\beta^{2}(t)\mu^{1/2}(t)e^{-(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},\quad
s_{i}(x,t)=-e_{i}a(t)\beta^{2}(t)\mu^{1/2}(t)e^{-(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},$
$L_{i}(t)=a_{i}a(t)\beta^{2}(t),\quad\quad
h_{i}(x,t)=-b_{i}a(t)\beta^{2}(t)\mu^{1/2}(t)e^{-(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},\quad
i=1,2,3,$

and the coefficients to satisfy the Riccati system (3.52)-(3.57) (and
therefore equations (3.16)-(3.17)), the functions $u(\xi,\tau),\ v(\xi,\tau),\
w(\xi,\tau),$ will satisfy the system (3.71)-(3.73). ∎

####  3.2.3. Constant diffusivity and interaction coefficients with
exponential growth

To exemplify this case, let us consider the system of equations

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\beta^{-2}(0)\psi_{xx}+\left(a_{1}-1-\operatorname{sech}t\tanh
t\right)\psi-\mu^{1/2}(0)e^{t-\operatorname{sech}t-\kappa(0)}\left(\psi+c_{1}\varphi+e_{1}\phi\right)\psi,$
(3.76) $\displaystyle\varphi_{t}$ $\displaystyle=$
$\displaystyle\beta^{-2}(0)\varphi_{xx}+\left(a_{1}-1-\operatorname{sech}t\tanh
t\right)\varphi-\mu^{1/2}(0)e^{t-\operatorname{sech}t-\kappa(0)}\left(b_{2}\psi+\varphi+e_{2}\phi\right)\varphi,$
(3.77) $\displaystyle\phi_{t}$ $\displaystyle=$
$\displaystyle\beta^{-2}(0)\phi_{xx}+\left(a_{1}-1-\operatorname{sech}t\tanh
t\right)\phi-\mu^{1/2}(0)e^{t-\operatorname{sech}t-\kappa(0)}\left(b_{3}\psi+c_{3}\varphi+\phi\right)\phi.$
(3.78)

Here, the Riccati system (3.52)-(3.57) admits the solution

$\displaystyle\alpha(t)=\delta(t)=0,\quad\beta(t)=\beta(0),\quad\gamma(t)=t+\gamma(0),\quad\varepsilon(t)=\varepsilon(0),\
\kappa(t)=\kappa(0),\quad\mu(t)=\mu(0)e^{2t-2\operatorname{sech}t}.$

In these terms, the proposed system (3.76)-(3.78) has solutions in the form

$\psi(x,t)=\left(2+\theta-\frac{a_{1}}{4}\right)\mu^{-1/2}(0)e^{-t-\operatorname{sech}t+\kappa(0)}\left[1-\tanh\left(\beta(0)x+\varepsilon(0)-\theta(t+\gamma(0))\right)\right]^{2},$
(3.79)
$\varphi(x,t)=\frac{a_{1}}{4}\mu^{-1/2}(0)e^{-t-\operatorname{sech}t+\kappa(0)}\left[1+\tanh\left(\beta(0)x+\varepsilon(0)-\theta(t+\gamma(0))\right)\right]^{2},$
(3.80)
$\phi(x,t)=\left(a_{1}-\theta-2\right)\mu^{-1/2}(0)e^{-t-\operatorname{sech}t+\kappa(0)}\left[1-\tanh\left(\beta(0)x+\varepsilon(0)-\theta(t+\gamma(0))\right)\right],$
(3.81)

provided $a_{1},\ b_{2},\ b_{3},\ c_{1},\ c_{3},\ e_{1},\ e_{2},$ satisfy the
same relations established in Section 2.

###  3.3. Solutions for the generalized Gray-Scott model

This subsection is focused on the study of the generalization of the classical
Gray-Scott system. The main result of this part is stated as follows:

######  Theorem 4 (Generalized Gray-Scott Model).

The Gray-Scott model with variable coefficients

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{1}(t)-xf(t))\psi-(g(t)-c\left(t\right)x)\psi_{x}-h_{1}(x,t)\psi\varphi^{2}+M_{1}(x,t),$
(3.82) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-(b\left(t\right)x^{2}-d(t)-L_{1}(t)-xf(t))\varphi-(g(t)-c\left(t\right)x)\varphi_{x}+h_{1}(x,t)\psi\varphi^{2}+M_{2}(x,t),$
(3.83)

can be transformed into the standard Gray-Scott model

$u_{\tau}=u_{\xi\xi}-uv^{2}+b_{1}(1-u),$ (3.84)
$v_{\tau}=v_{\xi\xi}+uv^{2}-b_{1}v+b_{2}.$ (3.85)

###### Proof.

As in the previous results, we use the substitutions

$\psi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}u(\xi,\tau),\qquad\xi=\beta(t)x+\varepsilon(t),\qquad\tau=\gamma(t),$
(3.86)

and

$\varphi(x,t)=\frac{1}{\sqrt{\mu(t)}}e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)}v(\xi,\tau).\qquad$
(3.87)

By assuming the coefficients satisfy the Riccati system (3.52)-(3.57) and
imposing the integrability conditions

$h_{1}(x,t)=a(t)\beta^{2}(t)\mu(t)e^{-2(\alpha(t)x^{2}+\delta(t)x+\kappa(t))},\quad
L_{1}(t)=-b_{1}a(t)\beta^{2}(t),$
$M_{i}(x,t)=b_{i}a(t)\beta^{2}(t)\mu^{-1/2}(t)e^{\alpha(t)x^{2}+\delta(t)x+\kappa(t)},\quad
i=1,2,$

is easy to see that the functions $u(\xi,\tau),\ v(\xi,\tau),$ satisfy the
Gray-Scott model (3.84)-(3.85). ∎

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 6. Solutions for the system (3.88)-(3.89) for the parameters
$b_{1}=\frac{1}{8},\ b_{2}=0$. Here, (b) and (d) show the contours of $\psi$
and $\varphi$ functions, respectively.

####  3.3.1. Constant diffusivity and an interaction coefficient with
exponential growth

If we assume a system as follows

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\psi_{xx}+\left(4x^{2}+2-\sin
t-b_{1}\right)\psi+4x\psi_{x}-e^{2x^{2}-2\cos
t}\psi\varphi^{2}+b_{1}e^{-x^{2}+\cos t},$ (3.88) $\displaystyle\varphi_{t}$
$\displaystyle=$ $\displaystyle\varphi_{xx}+\left(4x^{2}+2-\sin
t-b_{1}\right)\varphi+4x\varphi_{x}+e^{2x^{2}-2\cos t}\psi\varphi^{2}$ (3.89)

then, Theorem 4 allows the construction of solutions in the following closed
form:

$\psi(x,t)=\exp\left(-x^{2}+\cos
t\right)\left\\{\frac{3-\sqrt{1-4b_{1}}}{4}-\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(x-\theta
t)\right)\right\\},$ (3.90) $\varphi(x,t)=\exp\left(-x^{2}+\cos
t\right)\left\\{\frac{1+\sqrt{1-4b_{1}}}{4}+\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(x-\theta
t)\right)\right\\},$ (3.91)

with $\theta=\frac{\sqrt{2}(1-3\sqrt{1-4b_{1}})}{4}.$ The dynamics of these
solutions can be seen in Figure 6.

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 7. Solutions for the system (3.93)-(3.95) for the parameters
$b_{1}=\frac{1}{8},\ b_{2}=0$. The profiles of the functions $\psi$ and
$\varphi$ are shown in (a) and (c). Figures (b) and (d) show the contours of
$\psi$ and $\varphi$.

####  3.3.2. tanh-type diffusivity and an interaction coefficient with
exponential growth

We consider in this example the general Gray-Scott model

$\displaystyle\psi_{t}=\tanh t\left[\psi_{xx}-4x^{2}\psi+\psi-
b_{1}\operatorname{sech}^{8}t\psi-\operatorname{sech}^{6}te^{2(x-\cosh
t)^{2}}\psi\varphi^{2}+b_{1}\operatorname{sech}^{9}te^{-(x-\cosh
t)^{2}}\right]$ (3.92) $\displaystyle\quad\quad\quad\quad+\sinh
t\left(4x\psi-3\psi_{x}\right),$ (3.93) $\displaystyle\varphi_{t}=\tanh
t\left[\varphi_{xx}-4x^{2}\varphi+\varphi-
b_{1}\operatorname{sech}^{8}t\varphi+\operatorname{sech}^{6}te^{2(x-\cosh
t)^{2}}\psi\varphi^{2}+b_{1}\operatorname{sech}^{9}te^{-(x-\cosh
t)^{2}}\right]$ (3.94) $\displaystyle\quad\quad\quad\quad+\sinh
t\left(4x\varphi-3\varphi_{x}\right).$ (3.95)

By solving the system (3.52)-(3.57), we get

$\displaystyle\alpha(t)=-1,\quad\beta(t)=\operatorname{sech}^{4}t,\quad\gamma(t)=-\frac{\operatorname{sech}^{8}t}{8},\quad\delta(t)=2\cosh
t,$
$\displaystyle\varepsilon(t)=-\frac{1}{3}\operatorname{sech}^{3}t,\quad\quad\kappa(t)=-\cosh^{2}t,\quad\quad\mu(t)=\cosh^{2}t.$

Then, the system (3.93)-(3.95) admits a solution given by the equations

$\psi(x,t)=\operatorname{sech}te^{-(x-\cosh
t)^{2}}\left\\{\frac{3-\sqrt{1-4b_{1}}}{4}-\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(\xi-\theta\tau)\right)\right\\},$
(3.96) $\varphi(x,t)=\operatorname{sech}te^{-(x-\cosh
t)^{2}}\left\\{\frac{1+\sqrt{1-4b_{1}}}{4}+\frac{\sqrt{2+2\sqrt{1-4b_{1}}-4b_{1}}}{4}\tanh\left(\frac{\sqrt{1+\sqrt{1-4b_{1}}-2b_{1}}}{4}(\xi-\theta\tau)\right)\right\\},$
(3.97)

with $\theta=\frac{\sqrt{2}(1-3\sqrt{1-4b_{1}})}{4},$
$\xi=x\operatorname{sech}^{4}t-\frac{1}{3}\operatorname{sech}^{3}t$, and
$\tau=-\frac{1}{8}\operatorname{sech}^{8}t.$ In this situation, we produce
solutions with central axes bent to the right and decreasing in amplitude over
time, see Figure 7. This fascinating phenomenon results of the modulation of
the hyperbolic cosine and the exponential decay induced by the interaction
coefficient.

###  3.4. Solutions for the generalized Burgers system

This section concludes with the introduction of a novel coupled Burgers system
that generalizes on previous models reported in the literature. In this
context, the simplest version of the Riccati system must be used in the
construction of the explicit solutions for such a model. To be more exact, the
result says the following:

######  Theorem 5 (Generalized Burgers System).

The variable coefficient coupled Burgers system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\psi_{xx}-b_{1}a(t)\psi\psi_{x}-c_{1}a(t)(\psi\varphi)_{x}+c(t)\left(\psi+x\psi_{x}\right)-g(t)\psi_{x},$
(3.98) $\displaystyle\varphi_{t}$ $\displaystyle=$ $\displaystyle
a(t)\varphi_{xx}-b_{2}a(t)\varphi\varphi_{x}-c_{2}a(t)(\psi\varphi)_{x}+c(t)\left(\varphi+x\varphi_{x}\right)-g(t)\varphi_{x}$
(3.99)

can be transformed into the Burgers system

$u_{\tau}=u_{\xi\xi}-b_{1}uu_{\xi}-c_{1}(uv)_{\xi}$ (3.100)
$v_{\tau}=v_{\xi\xi}-b_{2}vv_{\xi}-c_{2}(uv)_{\xi}.$ (3.101)

###### Proof.

Consider the substitutions

$\psi(x,t)=\beta(t)u(\xi,\tau),\qquad\xi=\beta(t)x+\varepsilon(t),\qquad\tau=\gamma(t),$
(3.102)

and

$\varphi(x,t)=\beta(t)v(\xi,\tau).$ (3.103)

Now, computing the first derivatives of $\psi$ (similarly the derivatives of
$\varphi$ are obtained by replacing $u$ by $v$) we have:

$\psi_{x}=\beta^{2}(t)u_{\xi},\quad\quad\psi_{xx}=\beta^{3}(t)u_{\xi\xi},$
(3.104)
$\psi_{t}=\beta^{\prime}(t)u+\beta(t)\left[u_{\xi}\left(\beta^{\prime}(t)x+\varepsilon^{\prime}(t)\right)+u_{\tau}\gamma^{\prime}(t)\right].$
(3.105)

Inserting the last equations (3.104)-(3.105) into the system (3.98)-(3.99),
one obtains the Burgers system (3.100)-(3.101) by imposing the conditions

$\frac{d\beta}{dt}=c(t)\beta(t),\quad\quad\frac{d\gamma}{dt}=a(t)\beta^{2}(t),\quad\quad\frac{d\varepsilon}{dt}=-g(t)\beta(t).$
(3.106)

We point out that the system of ODEs (3.106) is a collapsed version with
$\alpha=\delta=\kappa=b=f=0$ of the classical Riccati system (3.52)-(3.57). ∎

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 8. Solutions for the system (3.107)-(3.108) for the parameters
$b_{1}=2,\ b_{2}=1,\ c_{1}=c_{2}=2$, $A=-\frac{111}{8},$ $B=-1,$ $\beta(0)=1,$
$\gamma(0)=\varepsilon(0)=0$.

In the preceding theorem, the integrability of the Burgers system is obtained
by equalizing the diffusion term and the interaction coefficient, $a(t)$.
Given the nature of the solutions expressed in (3.102)-(3.103), $c(t)$ has an
important influence on its dynamics (via the function $\beta(t)$), as it
affects horizontal translations and amplitudes. The following examples provide
insight into the dynamics of solutions when $c(t)$ has periodic structures.

####  3.4.1. Periodic diffusivity and interaction terms

In the present example we will consider the Burgers system

$\displaystyle\psi_{t}$ $\displaystyle=$ $\displaystyle\beta(0)^{-2}e^{-2\sin
t}\left(\psi_{xx}-b_{1}\psi\psi_{x}-c_{1}(\psi\varphi)_{x}\right)+\cos
t\left(\psi+x\psi_{x}-\psi_{x}\right),$ (3.107) $\displaystyle\varphi_{t}$
$\displaystyle=$ $\displaystyle\beta(0)^{-2}e^{-2\sin
t}\left(\varphi_{xx}-b_{2}\varphi\varphi_{x}-c_{2}(\psi\varphi)_{x}\right)+\cos
t\left(\varphi+x\varphi_{x}-\varphi_{x}\right).$ (3.108)

(a) Profile of the solution $\psi$.

(b) Contour of the solution $\psi$.

(c) Profile of the solution $\varphi$.

(d) Contour of the solution $\varphi$.

Figure 9. Solutions for the system (3.111)-(3.112) for the parameters
$b_{1}=2,\ b_{2}=1,\ c_{1}=c_{2}=2$, $A=\frac{111}{8},$ $B=1,$
$\beta(0)=\gamma(0)=\varepsilon(0)=1$.

The solution of the system (3.106) is given by the functions

$\displaystyle\beta(t)=\beta(0)e^{\sin
t},\quad\quad\gamma(t)=t+\gamma(0),\quad\quad\varepsilon(t)=\varepsilon(0)-\beta(0)\left(e^{\sin
t}-1\right).$

Then, we can construct solutions for the system (3.107)-(3.108) as follows

$\psi=\beta(0)e^{\sin
t}\left\\{B-2A\left(2c_{1}-\frac{1}{4c_{1}c_{2}}-1\right)\tanh\left[A\left(20(\beta(0)e^{\sin
t}(x-1)+\beta(0)+\varepsilon(0))-10-2A(t+\gamma(0))\right)\right]\right\\},$
(3.109) $\varphi=\beta(0)Be^{\sin
t}\left(2c_{1}-\frac{1}{4c_{1}c_{2}}-1\right)\tanh\left[A\left(20(\beta(0)e^{\sin
t}(x-1)+\beta(0)+\varepsilon(0))-10-2A(t+\gamma(0))\right)\right].$ (3.110)

The dynamics of these solutions are shown in Figure 8.

####  3.4.2. Diffusivity and interaction terms with exponential growth

Consider the system of equations

$\displaystyle\psi_{t}$ $\displaystyle=$
$\displaystyle\frac{\gamma(0)e^{t}}{\beta^{2}(0)\cos^{2}(\sin
t)}\left(\psi_{xx}-b_{1}\psi\psi_{x}-c_{1}(\psi\varphi)_{x}\right)-\tan(\sin
t)\cos t\left(\psi+x\psi_{x}\right)-\frac{\varepsilon(0)\sin
t}{\beta(0)\cos(\sin t)}\psi_{x},$ (3.111) $\displaystyle\varphi_{t}$
$\displaystyle=$ $\displaystyle\frac{\gamma(0)e^{t}}{\beta^{2}(0)\cos^{2}(\sin
t)}\left(\varphi_{xx}-b_{2}\varphi\varphi_{x}-c_{2}(\psi\varphi)_{x}\right)-\tan(\sin
t)\cos t\left(\varphi+x\varphi_{x}\right)-\frac{\varepsilon(0)\sin
t}{\beta(0)\cos(\sin t)}\varphi_{x}.$ (3.112)

Such a Burgers system admits the functions

$\displaystyle\beta(t)=\beta(0)\cos(\sin
t),\quad\quad\gamma(t)=\gamma(0)e^{t},\quad\quad\varepsilon(t)=\varepsilon(0)\cos
t.$

Then, solutions can be constructed in the following form:

$\psi=\beta(0)\cos(\sin
t)\left\\{B-2A\left(2c_{1}-\frac{1}{4c_{1}c_{2}}-1\right)\tanh\left[A\left(20(\beta(0)\cos(\sin
t)x+\varepsilon(0)\cos t)-10-2A\gamma(0)e^{t}\right)\right]\right\\},$ (3.113)
$\varphi=\beta(0)\cos(\sin
t)\left(2c_{1}-\frac{1}{4c_{1}c_{2}}-1\right)\tanh\left[A\left(20(\beta(0)\cos(\sin
t)x+\varepsilon(0)\cos t)-10-2A\gamma(0)e^{t}\right)\right].$ (3.114)

Figure 9 shows the profiles of these solutions.

##  4\. Conclusions and Final Remarks

In this paper, we introduced and investigated the explicit solutions of a
coupled reaction-diffusion system and a coupled Burgers-type system with
variable coefficients. The structure of such systems extends the classical
linear reaction-diffusion model, the diffusive Lotka-Volterra system, the
Gray-Scott model, and the Burgers equations. We showed that if the
coefficients fulfill a Riccati system, the general models possess solitary
wave solutions as well as solutions with nontrivial behavior (in relation to
solitary waves), presenting bending properties.

We point out that exact solutions always play a crucial role for any nonlinear
PDE (or system of PDEs) describing real-world processes. In fact, even exact
solutions with questionable applications can be used as test problems to
estimate the accuracy and efficiency of numerical methods. As a result, our
research makes a significant contribution in this aspect because the level of
trust in a numerical method for approximating thus general system of equations
increases if it can accurately replicate the dynamics of more complex
solutions.

On the other hand, it is important to mention that the results reported here
apply to any type of solution of the standard reaction-diffusion and Burgers
systems. Likewise, by appropriately selecting the functions $F_{i}$ in the
system (1.1)-(1.2), it can include new reaction-diffusion models such as the
Brusselator model [38, 40], the isothermal chemical system [19, 38], and the
Noyes-Field model [52]. Therefore, the ideas proposed in this paper can be
applicable to these new models.

######  Acknowledgement 1.

J. M. Escorcia thanks Universidad EAFIT for the financial support provided for
this research (Internal Project No. 12330022023).

##  5\. Appendix: Solutions of the Riccati Systems

In this appendix, we present the solutions of the Riccati systems used in the
construction of the explicit solutions of the various reaction-diffusion
equations.

###  5.1. Solution of the Riccati system (3.10)-(3.17)

A solution of this Riccati system including multiparameters is given by the
following expressions [12, 15, 58, 62]:

$\mu\left(t\right)=-2\mu\left(0\right)\mu_{0}\left(t\right)\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right),$
(5.1)
$\alpha\left(t\right)=\alpha_{0}\left(t\right)-\frac{\beta_{0}^{2}\left(t\right)}{4\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.2)
$\beta\left(t\right)=-\frac{\beta\left(0\right)\beta_{0}\left(t\right)}{2\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.3)
$\gamma\left(t\right)=\gamma\left(0\right)-\frac{\beta^{2}\left(0\right)}{4\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.4)
$\delta\left(t\right)=\delta_{0}\left(t\right)-\frac{\beta_{0}\left(t\right)\left(\delta\left(0\right)+\varepsilon_{0}\left(t\right)\right)}{2\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.5)
$\varepsilon\left(t\right)=\varepsilon\left(0\right)-\frac{\beta\left(0\right)\left(\delta\left(0\right)+\varepsilon_{0}\left(t\right)\right)}{2\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.6)
$\kappa\left(t\right)=\kappa\left(0\right)+\kappa_{0}\left(t\right)-\frac{\left(\delta\left(0\right)+\varepsilon_{0}\left(t\right)\right)^{2}}{4\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)},$
(5.7)

subject to the initial arbitrary conditions $\mu\left(0\right),$
$\alpha\left(0\right),$ $\beta\left(0\right)\neq 0,$ $\gamma(0),$ $\delta(0),$
$\varepsilon(0)$ and $\kappa(0)$. Here, $\alpha_{0}$, $\beta_{0}$,
$\gamma_{0}$, $\delta_{0}$, $\varepsilon_{0}$ and $\kappa_{0}$ are given
explicitly by

$\alpha_{0}\left(t\right)=-\frac{1}{4a\left(t\right)}\frac{\mu_{0}^{\prime}\left(t\right)}{\mu_{0}\left(t\right)}-\frac{d\left(t\right)}{2a\left(t\right)},$
(5.8)
$\beta_{0}\left(t\right)=\frac{W\left(t\right)}{\mu_{0}\left(t\right)},\quad
W\left(t\right)=\exp\left(\int_{0}^{t}\left(c\left(s\right)-2d\left(s\right)\right)\
ds\right),$ (5.9)
$\gamma_{0}\left(t\right)=-\frac{1}{2\mu_{1}\left(0\right)}\frac{\mu_{1}\left(t\right)}{\mu_{0}\left(t\right)}+\frac{d\left(0\right)}{2a\left(0\right)},$
(5.10)
$\delta_{0}\left(t\right)=\frac{W\left(t\right)}{\mu_{0}\left(t\right)}\ \
\int_{0}^{t}\left[\left(f\left(s\right)+\frac{d\left(s\right)}{a\left(s\right)}g\left(s\right)\right)\mu_{0}\left(s\right)+\frac{g\left(s\right)}{2a\left(s\right)}\mu_{0}^{\prime}\left(s\right)\right]\
\frac{ds}{W\left(s\right)},$ (5.11)
$\displaystyle\varepsilon_{0}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{2a\left(t\right)W\left(t\right)}{\mu_{0}^{\prime}\left(t\right)}\delta_{0}\left(t\right)-8\int_{0}^{t}\frac{a\left(s\right)\sigma\left(s\right)W\left(s\right)}{\left(\mu_{0}^{\prime}\left(s\right)\right)^{2}}\left(\mu_{0}\left(s\right)\delta_{0}\left(s\right)\right)\
ds$
$\displaystyle+2\int_{0}^{t}\frac{a\left(s\right)W\left(s\right)}{\mu_{0}^{\prime}\left(s\right)}\left[f\left(s\right)+\frac{d\left(s\right)}{a\left(s\right)}g\left(s\right)\right]\
ds,$ $\displaystyle\kappa_{0}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{a\left(t\right)\mu_{0}\left(t\right)}{\mu_{0}^{\prime}\left(t\right)}\delta_{0}^{2}\left(t\right)-4\int_{0}^{t}\frac{a\left(s\right)\sigma\left(s\right)}{\left(\mu_{0}^{\prime}\left(s\right)\right)^{2}}\left(\mu_{0}\left(s\right)\delta_{0}\left(s\right)\right)^{2}\
ds$
$\displaystyle\quad+2\int_{0}^{t}\frac{a\left(s\right)}{\mu_{0}^{\prime}\left(s\right)}\left(\mu_{0}\left(s\right)\delta_{0}\left(s\right)\right)\left[f\left(s\right)+\frac{d\left(s\right)}{a\left(s\right)}g\left(s\right)\right]\
ds,$

with
$\delta_{0}\left(0\right)=g_{0}\left(0\right)/\left(2a\left(0\right)\right),$
$\varepsilon_{0}\left(0\right)=-\delta_{0}\left(0\right),$
$\kappa_{0}\left(0\right)=0.$ Here $\mu_{0}$ and $\mu_{1}$ represent the
fundamental solution of the characteristic equation subject to the initial
conditions $\mu_{0}(0)=0$, $\mu_{0}^{\prime}(0)=2a(0)\neq 0$ and
$\mu_{1}(0)\neq 0$, $\mu_{1}^{\prime}(0)=0$.

We point out that all the formulas involved in such a solution have been
verified previously in [27].

###  5.2. Solution of the modified Riccati system (3.34)-(3.40)

In order to find the solution of this Riccati system, we just need to solve
the equations for the functions $\kappa_{1}(t)$ and $\kappa_{2}(t)$ because
the rest of the equations have been solved previously. Suppose
$\kappa(t)=\kappa_{1}(t)+\kappa_{2}(t).$ Then, the function $\kappa(t)$
satisfies the equation (3.15) of the earlier Riccati system, i.e.,

$\frac{d\kappa}{dt}=a(t)\delta^{2}(t)-g(t)\delta(t).$

On other hand, the equation for $\kappa_{2}(t)$ can be solved explicitly:

$\kappa_{2}(t)=-\ln\Big{|}\int_{0}^{t}h(s)\ ds+e^{-\kappa_{2}(0)}\Big{|}.$

In these terms, using equation (5.7), we have

$\kappa_{1}(t)=\kappa(t)-\kappa_{2}(t)=\kappa_{1}\left(0\right)+\kappa_{2}\left(0\right)+\kappa_{0}\left(t\right)-\frac{\left(\delta\left(0\right)+\varepsilon_{0}\left(t\right)\right)^{2}}{4\left(\alpha\left(0\right)+\gamma_{0}\left(t\right)\right)}+\ln\Big{|}\int_{0}^{t}h(s)ds+e^{-\kappa_{2}(0)}\Big{|},$

and we solved completely the system (3.34)-(3.40).

## References

  * [1] R. Abazaria and A. Borhanifar.  Numerical study of the solution of the Burgers and coupled Burgers equations by a differential transformation method.  Computers and Mathematics with Applications, 59:2711–2722, 2010.
  * [2] P. Acosta-Humánez and E. Suazo.  Liouvillian propagators, Riccati equation and differential Galois theory.  J. Phys. A: Math. Theor., 46:455203, 2013.
  * [3] M. O. Aibinua and S. Moyo.  Constructing exact solutions to systems of reaction-diffusion equations.  Int. J. Nonlinear Anal. Appl., 14 (1):585–595, 2023.
  * [4] G. Amador, K. Colon, N. Luna, G. Mercado, E. Pereira, and E. Suazo.  On solutions for linear and nonlinear Schrödinger equations with variable coefficients: A computational approach.  Symmetry, 8 (39), 2016.
  * [5] A. H. Askar, Á. Nagy, I. F. Barna, and E. Kovács.  Analytical and numerical results for the diffusion-reaction equation when the reaction coefficient depends on simultaneously the space and time coordinates.  Computation, 11 (127):1–27, 2023.
  * [6] N. F. Britton.  Essential mathematical biology.  Berlin: Springer, 2003.
  * [7] J. F. Cariñena, J. Grabowski, J. de Lucas, and C. Sardón.  Dirac–Lie systems and Schwarzian equations.  Journal of Differential Equations, 257(7):2303–2340, 2014.
  * [8] M. A. J. Chaplain.  Reaction–diffusion prepatterning and its potential role in tumour invasion.  Journal of Biological Systems, 03 (4):929–936, 1995.
  * [9] R. Cherniha and V. Davydovych.  Construction and application of exact solutions of the diffusive Lotka-Volterra system: A review and new results.  Communication in Nonlinear Science and Numerical Simulation, 113:106579, 2022.
  * [10] R. M. Cherniha and V. A. Dutka.  Diffusive Lotka–Volterra system: Lie symmetries and exact and numerical solutions.  Ukrainian Mathematical Journal, 56:1665–1675, 2004.
  * [11] C. Chou, Y. Zhang, R. Zhao, and Q. Nie.  Numerical methods for stiff reaction-diffusion systems.  Discrete and continuous Dynamical system- Series B, 7:515–525, 2007.
  * [12] R. Cordero-Soto, R.M. Lopez, E. Suazo, and S.K. Suslov.  Propagator of a charged particle with a spin in uniform magnetic and perpendicular electric fields.  Lett. Math. Phys., 84 (2-3):159–178, 2008.
  * [13] M. Dehghan, A. Hamidia, and M. Shakourifar.  The solution of coupled Burgers’ equations using Adomian–Pade technique.  Applied Mathematics and Computation, 189:1034–1047, 2007.
  * [14] A. J. Ellery, M. J. Simpson, S. W. McCue, and R. E. Baker.  Simplified approach for calculating moments of action for linear reaction-diffusion equations.  Physical Review E, 88:054102, 2013.
  * [15] J. Escorcia and E. Suazo.  Blow-up results and soliton solutions for a generalized variable coefficient nonlinear Schrödinger equation.  Applied Mathematics and Computation, 301:155–176, 2017.
  * [16] J. M. Escorcia and E. Suazo.  On explicit soliton solutions and blow-up for coupled variable coefficient nonlinear Schrödinger equations.  arXiv preprint arXiv:2310.14440, 2023.
  * [17] R. Field and M. Burger.  Oscillations and travelling waves in chemical systems.  Wiley, New York, 1985.
  * [18] M. Freidlin.  Coupled reaction-diffusion equations.  The Annals of Probability, 19 (1):29–57, 1991.
  * [19] C. M. Garcia-Lopez and J. I. Ramos.  Linearized methods, part II: reaction-diffusion equations’.  Comput. Methods Appl. Mech. Eng., 137:357–378, 1996.
  * [20] Z. Ghaemi, O. Nafiu, and E. Tajkhorshid et al.  A computational spatial whole-cell model for hepatitis B viral infection and drug interactions.  Sci Rep, 13:21392, 2023.
  * [21] P. Gray and S. K. Scott.  Autocatalytic reactions in the isothermal continuous stirred tank reactor.  Chemical Engineering Science, 39 (6):1087–1097, 1984.
  * [22] S. Habib, C. Molina-París, and T. S. Deisboeck.  Complex dynamics of tumors: modeling an emerging brain tumor system with coupled reaction-diffusion equations.  Physica A, 327:501–524, 2003.
  * [23] L. Hung.  Traveling wave solutions of competitive–cooperative Lotka–Volterra systems of three species.  Nonlinear Analysis: Real World Applications, 12:3691–3700, 2011.
  * [24] H. M. Jaradat.  Two-mode coupled Burgers equation: Multiple-kink solutions and other exact solutions.  Alexandria Engineering Journal, 57:2151–2155, 2018.
  * [25] W. Jiang, Z. Lu, and J. Wang.  Uniform patterns formation based on Gray-Scott model for 3D printing.  Computer Physics Communications, 295:108974, 2024.
  * [26] N. Kaur and V. Joshi.  Numerical solution to the Gray-Scott reaction-diffusion equation using hyperbolic B-spline.  Journal of Physics: Conference Series, 2267:012072, 2022.
  * [27] C. Koutschan, E. Suazo, and S. K. Suslov.  Fundamental laser modes in paraxial optics: from computer algebra and simulations to experimental observation.  Appl. Phys. B, 2015.
  * [28] Y. Kuang, J. D. Nagy, and S. E. Eikenberry.  Introduction to mathematical oncology.  Boca Raton: CRC Press, 2016.
  * [29] M. Kumar and S. Pandit.  A composite numerical scheme for the numerical simulation of coupled Burgers’ equation.  Computer Physics Communications, 185:809–817, 2014.
  * [30] Y. Kuramoto.  Chemical oscillations, waves and turbulence.  Springer, New York, 1984.
  * [31] R. T. Liu, S. S. Liaw, and P. K. Maini.  Two-stage Turing model for generating pigment patterns on the leopard and the jaguar.  Physical Review E, 74:011914, 2006.
  * [32] A. J. Lotka.  Undamped oscillations derived from the law of mass action.  J. Am. Chem. Soc., 42:1595–1599, 1920.
  * [33] W. Malfliet and W. Hereman.  The tahn method: I. exact solutions of nonlinear evolution and wave equations.  Phys. Scripta, 54:563–568, 1996.
  * [34] M.E. Marhic.  Oscillating hermite-gaussian wave functions of the harmonic oscillator.  Lett. Nuovo Cim., 22 (8):368–378, 1978.
  * [35] M. Maroua and B. Nabila.  Asymptotic behavior of solution for coupled reaction diffusion system by order m.  Int. J. Anal. Appl., 21 (37):1–12, 2023.
  * [36] B. T. Mbopda, S. Issa, R. Guiem, S. C. Oukouomi Noutchie, and H. P. Ekobena.  Travelling waves of a nonlinear reaction-diffusion model of the hepatitis B virus.  Eur. Phys. J. Plus, 138:971, 2023.
  * [37] R.C. Mittal and G. Arora.  Numerical solution of the coupled viscous Burgers equation.  Commun Nonlinear Sci Numer Simulat, 16:1304–1313, 2011.
  * [38] R.C. Mittal and R. Rohila.  Numerical simulation of reaction-diffusion systems by modified cubic B-spline differential quadrature method.  Chaos, Solitons and Fractals, 92:9–19, 2016.
  * [39] J. Morgan and B. Q. Tang.  Global well-posedness for volume-surface reaction-diffusion systems.  Communications in Contemporary Mathematics, 25 (04):2250002, 2023.
  * [40] G. Nicolis and I. Prigogine.  Exploring complexity.  New York: Freeman, 1989.
  * [41] A. Okubo and S. A. Levin.  Diffusion and ecological problems: Modern perspectives.  Second ed. Berlin: Springer, 2001.
  * [42] A. G. Osborne and M. R. Deinert.  Stability instability and Hopf bifurcation in fission waves.  Cell Reports Physical Science, 2:100588, 2021.
  * [43] C. V. Pao.  Global asymptotic stability of Lotka–Volterra competition systems with diffusion and time delays.  Nonlinear Analysis:Real World Applications, 5:91–104, 2004.
  * [44] C. V. Pao and Y. Wang.  Numerical solutions of a three-competition Lotka–Volterra system.  Applied Mathematics and Computation, 204:423–440, 2008.
  * [45] O. Pashaev and G. Tanoglu.  Vector shock soliton and the Hirota bilinear method.  Chaos, Solitons and Fractals, 26:95–105, 2005.
  * [46] E. Pereira, E. Suazo, and J. Trespalacios.  Riccati-Ermakov systems and explicit solutions for variable coefficient reaction-diffusion equations.  Applied Mathematics and Computation, 329:278–296, 2018.
  * [47] A. D. Polyanin.  Construction of exact solutions in implicit form for pdes: New functional separable solutions of non-linear reaction–diffusion equations with variable coefficients.  International Journal of Non-Linear Mechanics, 111:95–105, 2019.
  * [48] A. D. Polyanin.  Functional separable solutions of nonlinear reaction–diffusion equations with variable coefficients.  Applied Mathematics and Computation, 347:282–292, 2019.
  * [49] A. D. Polyanin and A. I. Zhurov.  Multi-parameter reaction–diffusion systems with quadratic nonlinearity and delays: New exact solutions in elementary functions.  Mathematics, 10:1529, 2022.
  * [50] A.D. Polyanin and A.I. Zhurov.  Separation of variables and exact solutions to nonlinear pdes.  Chapman and Hall/CRC, 2021.
  * [51] M. Rodrigo and M. Mimura.  Exact solutions of a competition-diffusion system.  Hirohima Math J, 30:257–270, 2000.
  * [52] M. Rodrigo and M. Mimura.  Exact solutions of reaction-diffusion systems and nonlinear wave equations.  Japan J. Indust. Appl. Math., 18:657–696, 2001.
  * [53] H. Shoji, Y. Iwasa, A. Mochizuki, and S. Kondo.  Directionality of stripes formed by anisotropic reaction-diffusion models.  J. theor. Biol., 214:549–561, 2006.
  * [54] M. J. Simpson and K. A. Landman.  Analysis of split operator methods applied to reactive transport with monod kinetics.  Advances in Water Resources, 30:2026–2033, 2007.
  * [55] M. J. Simpson, J. A. Sharp, L. C. Morrow, and R. E. Baker.  Exact solutions of coupled multispecies linear reaction–diffusion equations on a uniformly growing domain.  PLoS ONE, 10 (9):e0138894, 2015.
  * [56] B. D. Sleeman and E. Tuma.  Comparison principles for strongly coupled reaction-diffusion equations.  Proceedings of the Royal Society of Edinburgh, 1:209–219, 1987.
  * [57] L. W. Somathilake and J. M. J. J. Peiris.  Global solutions of a strongly coupled reaction-diffusion system with different diffusion coefficients.  Journal of Applied Mathematics, 1:23–36, 2005.
  * [58] E. Suazo.  Fundamental solutions of some evolution equations.  Arizona State University, Sep. 2009 (Phd. dissertation).
  * [59] E. Suazo and S. K. Suslov.  Soliton-like solutions for the nonlinear Schrödinger equation with variable quadratic Hamiltonians.  Journal of Russian Laser Research, 33 (1), 2012.
  * [60] E. Suazo, S. K. Suslov, and J. M. Vega-Guzmán.  The Riccati system and a diffusion-type equation.  Mathematics, 2:96–118, 2014.
  * [61] E. Suazo and S.K. Suslov.  An integral form of the nonlinear Schrödinger equation with variable coefficients.  Progress in Electromagnetics Research Symposium (PIERS-Toyama), pages 1214–1220, 2018.
  * [62] S. K. Suslov.  On integrability of nonautonomous nonlinear Schrödinger equations.  Am. Math. Soc., 140 (9):3067–3082, 2012.
  * [63] L. Tang and S. Chen.  Traveling wave solutions for the diffusive Lotka–Volterra equations with boundary problems.  Applied Mathematics and Computation, 413:126599, 2022.
  * [64] G. Tanoglu.  Hirota method for solving reaction-diffusion equations with generalized nonlinearity.  International Journal of Nonlinear Science, 1 (1):30–36, 2006.
  * [65] A. M. Turing.  The chemical basis of morphogenesis.  Bull. Math. Biol., 52(1):153–197, 1990.
  * [66] V. Volterra.  Variazioni e fluttuazioni del numero d’individui in specie animali conviventi.  Mem. Acad. Lincei., 2:31–113, 1926.
  * [67] K. Wang and W. Wang.  Propagation of HBV with spatial dependence.  Math. Biosci., 210(1):78–95, 2007.
  * [68] A. Wazwaz.  Multiple-front solutions for the Burgers equation and the coupled Burgers equations.  Applied Mathematics and Computation, 190:1198–1206, 2007.
  * [69] Y. Zhanga, P. Liu, and W. Hou.  Modeling of glioma growth using modified reaction-diffusion equation on brain MR images.  Computer Methods and Programs in Biomedicine, 227:107233, 2022.


aainstitutetext: Department of Physics, Northeastern University, Boston, MA
02115, USAbbinstitutetext: Computational Science Initiative, Brookhaven
National Laboratory, Upton, NY 11973, USA

# Holographic Entanglement Distillation from the Surface State Correspondence

Ning Bao a Gün Süer<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study correlations between geometric subfactors living on the Ryu-
Takayanagi surface that bounds the entanglement wedge. Using the surface-state
correspondence and the bit threads program, we are able to calculate mutual
information and conditional mutual information between subfactors. This
enables us to count the shared Bell pairs between subfactors, and we propose
an entanglement distillation procedure over these subsystems via a SWAP gate
protocol. We comment on extending to multipartite entanglement.

##  1 Introduction

In recent years there has been great progress made at the intersection of
quantum information and quantum gravity. One of the major areas of interest
has been the study of entanglement entropies in holographic settings Ryu_2006
. In this context, the Ryu-Takayanagi (RT) formula relates the entanglement
entropy of a boundary subregion $A$ to the area of the minimal surface
$\gamma$ in the bulk homologous to that region

$S(A)=\min_{\partial\gamma=A}\frac{\text{Area}(\gamma)}{4G_{N}}.$ (1)

In addition to its conceptual appeal in generalizing the area law for black
hole entropy to entanglement entropies in conformal field theories, it also
enables purely geometrical proofs of various information-theoretic relations
such as strong subadditivity Headrick_2007  and higher party number
generalizations of entanglement entropy inequalities for holographic states
Bao_2015 .

Since entanglement is a central resource of quantum information processing
systems, the question of available entanglement in a given system is of clear
importance Chitambar_2019 . Given an arbitrarily entangled quantum state, the
process of transforming the initial state to a certain number of Bell pairs
using local operations and classical communication (LOCC) is known as
entanglement distillation Deutsch:1996um . While algorithms for performing
entanglement distillation for pure states are known, for generic mixed states
even the amount of distillable entanglement is not known in general.

In this direction, a more recent entry in the holographic dictionary is the
$E_{W}=E_{P}$ conjecture Nguyen_2018 ; Umemoto_2018 . It states that the
entanglement wedge cross-section $E_{W}$ separating two regions is equal to
the entanglement of purification $E_{P}$, on the basis that the holographic
object $E_{W}$ obeys the same set of inequalities as the information-theoretic
object $E_{P}$. Given a bipartite state $\rho_{AB}$, the entanglement of
purification is defined as

$E_{P}(A:B)=\min\left\\{S(AA^{\prime})\ |\ \rho_{AA^{\prime}BB^{\prime}}\
\text{is pure}\right\\},$ (2)

which measures correlations in a bipartite mixed state. Since it involves
optimizing over the purification, it is very challenging to compute. Therefore
finding a holographic dual to entanglement of purification and calculating it
using geometrical tools is appealing. It has been conjectured that the
entanglement of purification is exactly equal to the area of the entanglement
wedge cross-section which is the minimal surface that splits the entanglement
wedge $r_{AB}$ into two regions

$E_{W}(A:B)=\min\left\\{\text{Area}(\Gamma_{AB})\ |\ \Gamma_{AB}\subset
r_{AB}\text{ splits }r_{AB}\text{ into two regions}\right\\}.$ (3)

In this paper we focus on understanding entanglement distillation in a
holographic setting, using the surface state correspondence Miyaji_2015  and
the bit threads program Freedman_2016 . Specifically, we consider two boundary
region entanglement wedges and calculate mutual information and conditional
mutual information for subfactors living on RT surfaces. The states on the
surfaces in the bulk are related to states on the boundary via an isometry,
which is a local unitary preserving the entanglement properties. By selecting
pairs of qubits across different subsystems and performing SWAP tests, we
present a straightforward entanglement distillation procedure that works at
least asymptotically.

Figure 1: Under the surface-state correspondence, the boundary can be pushed
to a closed surface in the bulk via an isometry. The isometry used is a local
unitary acting only on the subregion being pushed, preventing the change of
entanglement properties.

##  2 Surface State Correspondence

The surface state correspondence is a conjectured property of the AdS/CFT
correspondence. It states that one can define a pure state on any closed
surface not containing topological obstructions in the bulk. This surface is
related to the boundary by continuous deformation, or pushing inwards from the
boundary into the bulk. Such pushing, again in the absence of topological
obstructions, will always result in a closed surface; see Fig. 1. The boundary
state is then related to the state on this closed bulk surface via isometry.
In particular, we are allowed to push only a geometric subregion of the
boundary, while keeping the remainder of the boundary fixed; see Fig. 2. At
the conclusion of this pushing, the new closed surface functions as the new
boundary, for the purposes of homology constraints of the RT formula. Once
this pushing has been performed, a new boundary subregion of this redefined
boundary, here taken to not overlap with the first boundary subregion
selected, can be chosen for pushing.111The non-overlapping condition of this
pushing is a special instance of the tree tensor network style pushing of
Bao_20191 .

There will be a portion of the state which is discarded in the process of
implementing the isometry, e.g. a portion of the state which was defined on a
subfactor of the original boundary Hilbert space but is not supported on the
Hilbert space of the closed curve. This portion of the state will be
associated with the region that the isometry pushed through, and in particular
will not be entangled with either the surface state or the remainder of the
boundary CFT state, due to the purity of the surface state and the fact that
the isometry used was a local unitary acting only on the subregion being
pushed and not its complement, preventing the change of entanglement
properties across the bipartition.

In this work, we will focus on pushing from boundary subregions to their RT
surfaces. By purity of the overall closed surface state, this pushing must
preserve the entanglement entropy of the boundary region. However, because the
RT surface has an area precisely equal to the entanglement entropy of the
boundary region up to normalization, by the bottleneck argument in the
equivalent bit threads formalism Freedman_2016  the log of the Hilbert space
dimension of the portion of the state associated with a RT surface is equal to
its entanglement entropy. This means that the portion of the state living on
an RT surface in the surface state correspondence must be maximally mixed.

We will leverage this argument to derive several strong results on
entanglement distillation for mixed states in holography.

##  3 Two Boundary Region Entanglementment Wedges

Let’s start with a connected entanglement wedge bounded by two disconnected
boundary regions, $A$ and $B$. We restrict our attention to this case because
the case of the disconnected entanglement wedge is trivial for our
considerations. We will first push the complementary boundary region to the RT
surface of $\rho_{AB}$, as for pure CFT states the RT surface of the
complement is the same as that of $AB$. Because surface-state implies the
$E_{P}=E_{W}$ conjecture Umemoto_2018 , we can associate the $A^{\prime}$
factor of the purification to the portion of the RT surface on the same side
of the $E_{W}$ surface as $A$, and the $B^{\prime}$ factor of the purification
to the portion of the RT surface on the same side of the $E_{W}$ surface as
$B$. Note that by the arguments of the previous section,
$\rho_{A^{\prime}B^{\prime}}$ must be maximally mixed, to leading order in
$1/N$.

At this point, we will use the iterative step from the previous section to go
further. Keeping the $A^{\prime}B^{\prime}$ portion of the surface fixed, we
can now also push $\rho_{A}$ and $\rho_{B}$ to their respective RT surfaces;
see Fig. 3. We will call the states here $\rho_{A_{D}}$ and $\rho_{B_{D}}$.
Similarly to above, $\rho_{A_{D}}$ and $\rho_{B_{D}}$ must also be maximally
mixed, to leading order in $1/N$.

Figure 2: One can also consider only pushing a subregion while keeping its
complement fixed. This preserves the entanglement properties across the
bipartition.

###  3.1 Mutual Information

Now we can compute mutual information between pairs of these regions. First,
and most simply, $I(A^{\prime}:B^{\prime})=0$, as the mutual information
between two subfactors of a maximally mixed state is always zero.

Now we turn our attention to $I(A:B)$. Because only the $A_{D}$ and $B_{D}$
subfactors can contribute to this mutual information, we have that

$I(A:B)=I(A_{D}:B_{D})=S(A_{D})+S(B_{D})-S(A_{D}B_{D})=S(A)+S(B)-S(AB).$ (4)

This might seem trivial, but it’s important to note that due to $\rho_{A_{D}}$
and $\rho_{B_{D}}$ being maximally mixed, their entanglement entropies are the
logs of their Hilbert space dimensions. Moreover, because they are maximally
mixed, the only internal purification in $\rho_{A_{D}B_{D}}$ allowed are Bell
pairs formed between qubits in $A_{D}$ and qubits in $B_{D}$. Therefore, this
mutual information counts the number of Bell pairs between qubits in $A_{D}$
and those in $B_{D}$, as it sums the individual entropies which count the
number of qubits in $A_{D}$ and $B_{D}$, and subtracts off those that are not
internally purified. Because the mutual information is an upper bound on the
distillable entanglement between $A_{D}$ and $B_{D}$, this is sufficient to
show that to leading order the holographic distillable entanglement is equal
to $I(A_{D}:B_{D})=I(A:B)$, as these Bell pairs between $A_{D}$ and $B_{D}$
are certainly distillable across the bipartition.

Figure 3: We push the regions $A$ and $B$ to their respective RT surfaces,
$A_{D}$ and $B_{D}$. $A^{\prime}B^{\prime}$ is the entanglement wedge. $E_{W}$
is the minimal entanglement wedge cross-section. It bounds the regions
$A^{\prime}$ and $B^{\prime}$ with $A$ and $B$ respectively.

We can also compute $I(A:A^{\prime})$. By the same logic before,

$I(A:A^{\prime})=S(A_{D})+S(A^{\prime})-S(AA^{\prime}).$ (5)

Here we mix the $A_{D}$ and $A$ Hilbert space factors for convenience. The
subtracted term is clearly the $E_{P}$ quantity.222Again, we are using the
fact that surface state implies $E_{P}=E_{W}$. But, by the same arguments in
considering $I(A:B)$, this quantity computes the amount of distillable
entanglement between $A$ and $A^{\prime}$.

Finally, we come to $I(A:B^{\prime})$. By the same arguments as above, this
computes the amount of distillable entanglement between $A$ and $B^{\prime}$.
However, the expanded equation merits further study:

$I(A:B^{\prime})=S(A_{D})+S(B^{\prime})-S(AB^{\prime}).$ (6)

The first two terms are simply the area of the RT surface of $A$ and the
$B^{\prime}$ portion of the RT surface for $B^{\prime}$, respectively.
However, the subtracted term,
$S(AB^{\prime})=\min(S(A)+S(B^{\prime}),S(A^{\prime})+S(B))$, as these are the
two candidate RT surfaces within the surface state region. The former term
would return zero, which means that wlog either $I(A:B^{\prime})=0$ or
$I(B:A^{\prime})=0$, for all holographic states at leading order.333For
symmetric cases both terms in the minimization would be zero, but generically
one of them will be nonzero. This would also mean that the distillable
entanglement between these systems is also zero, as again the mutual
information upper bounds the distillable entanglement. In the nonzero
instance, $I(A:B^{\prime})=S(A)-S(A^{\prime})-(S(B)-S(B^{\prime}))$.

###  3.2 Entanglement Distillation Procedure

While for generic mixed states entanglement distillation procedures are not
known, for this holographic case it appears quite straightforward, at least
asymptotically. We describe the procedure for $\rho_{AB}$, but it works the
same way for the other choices of subsystems. Select a qubit in $A_{D}$ and a
qubit in $B_{D}$ and consider this two-qubit state, while also selecting two
other qubits, both in $A_{D}$, and consider its two-qubit state. The former
can either be a Bell pair or a maximally mixed state, and the latter can only
be a maximally mixed state.444This is true up to subleading corrections, so a
more precise statement would be that the trace distance between this state and
the maximally mixed state is zero up to $1/N$ corrections. Now simply perform
the SWAP test between these two two-qubit states. With enough copies of the
CFT, one can repeat this for all possible two-qubit states with one qubit from
$A_{D}$ and one from $B_{D}$, and find all of the Bell pairs in worst case
$O(N_{A_{D}}N_{B_{D}})$ trials, and from there on simply read off the Bell
pairs from the remaining infinite number of copies of states, asymptotically
achieving $I(A:B)$.

###  3.3 Conditional Mutual Information

We can also compute conditional mutual information using the same technique,
and we describe these calculations below, with minimal comment.

There are twelve distinct conditional mutual informations:
$I(A:B|A^{\prime})$, $I(A:B|B^{\prime})$, $I(A:B^{\prime}|B)$,
$I(A:B^{\prime}|A^{\prime})$, $I(A:A^{\prime}|B)$,
$I(A:A^{\prime}|B^{\prime})$, $I(A^{\prime}:B|A)$,
$I(A^{\prime}:B|B^{\prime})$, $I(A^{\prime}:B^{\prime}:A)$,
$I(A^{\prime}:B^{\prime}:B)$, $I(B:B^{\prime}|A)$, and
$I(B:B^{\prime}|A^{\prime})$. Wlog, we take
$S(A^{\prime}B)=S(A^{\prime})+S(B)=S(AB^{\prime})$.

In order,

$\displaystyle
I(A:B|A^{\prime})=S(AA^{\prime})+S(A^{\prime}B)-S(A^{\prime})-S(B^{\prime})$
(7) $\displaystyle=E_{W}(A:B)+S(B)-S(B^{\prime}),$ (8)

$\displaystyle
I(A:B|B^{\prime})=S(AB^{\prime})+S(B^{\prime}B)-S(A^{\prime})-S(B^{\prime})$
(9) $\displaystyle=E_{W}(A:B)+S(B)-S(B^{\prime}).$ (10)

The first two are the same.

$\displaystyle I(A:B^{\prime}|B)=S(AB)+S(BB^{\prime})-S(B)-S(A^{\prime})$ (11)
$\displaystyle=E_{W}(A:B)+S(B^{\prime})-S(B),$ (12)

$\displaystyle
I(A:B^{\prime}|A^{\prime})=S(AA^{\prime})+S(A^{\prime}B^{\prime})-S(B)-S(A^{\prime})$
(13) $\displaystyle=E_{W}(A:B)+S(B^{\prime})-S(B).$ (14)

These two are also the same.

$\displaystyle I(A:A^{\prime}|B)=S(AB)+S(A^{\prime}B)-S(B)-S(B^{\prime})$ (15)
$\displaystyle=S(AB)+S(A^{\prime})-S(B^{\prime})$ (16)
$\displaystyle=2S(A^{\prime}).$ (17)

$\displaystyle
I(A:A^{\prime}|B^{\prime})=S(AB^{\prime})+S(A^{\prime}B^{\prime})-S(B^{\prime})-S(B)$
(18) $\displaystyle=S(AB)+S(A^{\prime})-S(B^{\prime})-S(B)$ (19)
$\displaystyle=2S(A^{\prime}).$ (20)

These two are the same.

$\displaystyle I(A^{\prime}:B|A)=S(AA^{\prime})+S(AB)-S(A)-S(B^{\prime})$ (21)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A),$ (22)

$\displaystyle
I(A^{\prime}:B|B^{\prime})=S(A^{\prime}B^{\prime})+S(BB^{\prime})-S(A)-S(B^{\prime})$
(23) $\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (24)

These two are the same.

$\displaystyle
I(A^{\prime}:B^{\prime}|A)=S(AA^{\prime})+S(AB^{\prime})-S(A)-S(B)$ (25)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A),$ (26)

$\displaystyle
I(A^{\prime}:B^{\prime}|B)=S(BB^{\prime})+S(AB^{\prime})-S(A)-S(B)$ (27)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (28)

These two are the same.

$\displaystyle I(B:B^{\prime}|A)=S(AB)+S(AB^{\prime})-S(A)-S(A^{\prime})$ (29)
$\displaystyle=2S(A^{\prime})+S(B^{\prime})+S(B)-S(A),$ (30)

$\displaystyle
I(B:B^{\prime}|A^{\prime})=S(AB)+S(AB^{\prime})-S(A)-S(A^{\prime})$ (31)
$\displaystyle=2S(A^{\prime})+S(B^{\prime})+S(B)-S(A).$ (32)

These two are the same.

It is worth noting that all of these quantities must be non-negative, and
moreover must lower bound the mutual information, thereby giving some lower
bounds for the distillable entanglement.

###  3.4 Tripartite Information

We can also compute tripartite information. There are four tripartite
informations to keep track of: $I(A:B:A^{\prime})$, $I(A:B:B^{\prime})$,
$I(A:A^{\prime}:B^{\prime})$, and $I(B:A^{\prime}:B^{\prime})$. Wlog, we are
taking $S(A^{\prime}B)=S(A^{\prime})+S(B)=S(AB^{\prime})$.

In order,

$\displaystyle
I(A:B:A^{\prime})=S(AB)+S(AA^{\prime})+S(A^{\prime}B)-S(A)-S(B)-S(A^{\prime})-S(B^{\prime})$
(33) $\displaystyle=S(AA^{\prime})+S(A^{\prime}B)-S(A)-S(B)$ (34)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (35)

$\displaystyle
I(A:B:B^{\prime})=S(AB)+S(AB^{\prime})+S(BB^{\prime})-S(A)-S(B)-S(A^{\prime})-S(B^{\prime})$
(36) $\displaystyle=S(BB^{\prime})+S(AB^{\prime})-S(A)-S(B)$ (37)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (38)

$\displaystyle
I(A:A^{\prime}:B^{\prime})=S(AA^{\prime})+S(AB^{\prime})+S(A^{\prime}B^{\prime})-S(A)-S(A^{\prime})-S(B^{\prime})-S(B)$
(39) $\displaystyle=S(AA^{\prime})+S(AB^{\prime})-S(A)-S(B)$ (40)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (41)

$\displaystyle
I(B:A^{\prime}:B^{\prime})=S(BA^{\prime})+S(BB^{\prime})+S(A^{\prime}B^{\prime})-S(A)-S(A^{\prime})-S(B^{\prime})-S(B)$
(42) $\displaystyle=S(BB^{\prime})+S(AB^{\prime})-S(A)-S(B)$ (43)
$\displaystyle=E_{W}(A:B)+S(A^{\prime})-S(A).$ (44)

In the end, all four tripartite contributions are the same as the Markov gap
noted in Hayden_2021 , as expected.

##  4 Three Boundary Region Entanglement Wedges

We can extend the analysis from before to three boundary regions, $A$, $B$,
and $C$, with $A^{\prime}$, $B^{\prime}$, and $C^{\prime}$ defined as in
bao2018holographic . We will again do the pushing to the RT surface of $ABC$,
$A$, $B$, and $C$ in the analogous way.

###  4.1 Mutual Information

We can now compute mutual information in analogy to before. The mutual
information of the primed regions remains zero, as before:
$I(A^{\prime}:B^{\prime}C^{\prime})=0$. The mutual information $I(A:BC)$ is
the same and still computes the distillable entanglement between $A$ and $BC$.

Differences arise when we consider $I(A:B^{\prime})$,

$I(A:B^{\prime})=S(A)+S(B^{\prime})-\min(S(A)+S(B^{\prime}),S(B)+S(A^{\prime})+S(CC^{\prime})).$
(45)

So, if Wlog $I(A:B^{\prime})\neq 0$, this equals
$S(A)-S(A^{\prime})-(S(B)-S(B^{\prime}))-S(CC^{\prime})$, which again must be
nonnegative by dint of being a mutual information, providing a nontrivial
bound. Note that this does not become zero in the symmetric case, as with the
two boundary situations.

There are too many different entropic combinations to consider them all for
the three-party case, but clearly, there is a rich structure here that is
worthy of further exploration.

##  5 Conclusion and Outlook

In this work, we have investigated the bipartite and tripartite entanglement
properties of subfactors on RT surfaces inspired by the recently proposed
$E_{W}=E_{P}$ conjecture. By using the surface state correspondence and the
bit threads program we calculate mutual information and conditional mutual
information across different subfactors living on the boundary of the
entanglement wedge. After estimating the number of shared Bell pairs across
different regions, we employ a SWAP test for entanglement distillation, using
the fact that the states on RT surfaces in the surface state correspondence
are maximally mixed. We leave entanglement properties and distillation in
multipartite systems for future work.

###### Acknowledgements.

We are grateful to Charles Cao and Jonathan Harper for discussions. N.B. is
supported by the Computational Science Initiative at Brookhaven National
Laboratory, Northeastern University, and by the U.S. Department of Energy
QuantISED Quantum Telescope award. G.S. is supported by the Graduate
Assistantship from the Department of Physics, Northeastern University.

## References

  * (1) S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from AdS/CFT_ , _Phys. Rev. Lett._ 96 (2006) 181602 [hep-th/0603001].
  * (2) M. Headrick and T. Takayanagi, _A Holographic proof of the strong subadditivity of entanglement entropy_ , _Phys. Rev. D_ 76 (2007) 106013 [0704.3719].
  * (3) N. Bao, S. Nezami, H. Ooguri, B. Stoica, J. Sully and M. Walter, _The Holographic Entropy Cone_ , _JHEP_ 09 (2015) 130 [1505.07839].
  * (4) E. Chitambar and G. Gour, _Quantum resource theories_ , _Reviews of Modern Physics_ 91 (2019) .
  * (5) D. Deutsch, A. Ekert, R. Jozsa, C. Macchiavello, S. Popescu and A. Sanpera, _Quantum privacy amplification and the security of quantum cryptography over noisy channels_ , _Phys. Rev. Lett._ 77 (1996) 2818 [quant-ph/9604039].
  * (6) P. Nguyen, T. Devakul, M.G. Halbasch, M.P. Zaletel and B. Swingle, _Entanglement of purification: from spin chains to holography_ , _JHEP_ 01 (2018) 098 [1709.07424].
  * (7) T. Takayanagi and K. Umemoto, _Entanglement of purification through holographic duality_ , _Nature Phys._ 14 (2018) 573 [1708.09393].
  * (8) M. Miyaji and T. Takayanagi, _Surface/state correspondence as a generalized holography_ , _Progress of Theoretical and Experimental Physics_ 2015 (2015) 073B03.
  * (9) M. Freedman and M. Headrick, _Bit threads and holographic entanglement_ , _Commun. Math. Phys._ 352 (2017) 407 [1604.00354].
  * (10) N. Bao, G. Penington, J. Sorce and A.C. Wall, _Beyond toy models: distilling tensor networks in full AdS/CFT_ , _Journal of High Energy Physics_ 2019 (2019) .
  * (11) P. Hayden, O. Parrikar and J. Sorce, _The markov gap for geometric reflected entropy_ , _Journal of High Energy Physics_ 2021 (2021) .
  * (12) N. Bao and I.F. Halpern, _Holographic Inequalities and Entanglement of Purification_ , _JHEP_ 03 (2018) 006 [1710.07643].