# Approximate symmetries of the Harry Dym equation

Mehdi Nadjafikhah, Parastoo Kabi-Nejad

###### Abstract.

In this paper, we derive the first order approximate symmetries for the Harry
Dym equation by the method of approximate transformation groups proposed by
Baikov, Gaszizov and Ibragimov [1], [2]. Moreover, we investigate the
structure of the Lie algebra of symmetries of the perturbed Harry Dym
equation. We compute the one-dimensional optimal system of subalgeras as well
as point out some approximately differential invarints with respect to the
generators of Lie algebra and optimal system.

##  1\. Introduction

The following nonlinear partial differential equation

(1) $\displaystyle u_{t}=-\frac{1}{2}u^{3}u_{xxx}.$

is known as the Harry Dym equation [7]. This equation was obtained by Harry
Dym and Martin Kruskal as an evolution equation solvable by a spectral problem
based on the string equation instead of Schrödinger equation. This result was
reported in [9] and rediscovered independently in [15], [16]. The Harry Dym
equation shares many of the properties typical of the soliton equations. It is
a completely integrable equation which can be solved by inverse scattering
transformation [3], [17], [18]. It has a bi-Hamiltonian structure and an
infinite number of conservation laws and infinitely many symmetries [10],
[11].

In this paper, we analyze the perturbed Harry Dym equation

(2) $\displaystyle u_{t}+\frac{1}{2}u^{3}u_{xxx}+\varepsilon u_{x}=0,$

where $\varepsilon$ is a small parameter, with a method which was first
introduced by Baikov, Gazizov and Ibragimov [1], [2]. This method which is
known as ”approximate symmetry” is a combination of Lie group theory and
perturbations. There is a second method which is also known as ”approximate
symmetry” due to Fushchich anf Shtelen [6] and later followed by Euler et al
[4], [5]. For a comparison of these two methods, we refer the interested
reader to the papers [12], [14]. our paper is organized as follows: In section
2, we present some definitions and theorems in the theory of approximate
symmetry. In section 3, we obtain the approximate symmetry of the perturbed
Harry Dym equation. In section 4, we discuss on the structure of its Lie
algebra. In section 5, we construct the one-dimensional optimal system of
subalgebras. In section 6, we compute some approximately differential
invarints with respect to the generators of Lie algebra and optimal system. In
section 7, we summarize our results.

##  2\. Notations and Definitions

In this section, we will provide the background definitions and results in
approximate symmetry that will be used along this paper. Much of it is stated
as in [8]. If a function$f(x,\varepsilon)$ satisfies the condition

(3) $\displaystyle\lim\dfrac{f(x,\varepsilon)}{\varepsilon^{p}}=0,$

it is written $f(x,\varepsilon)=o(\varepsilon^{p})$ and f is said to be
oforder less than $\varepsilon^{p}$. If

(4) $\displaystyle f(x,\varepsilon)-g(x,\varepsilon)=o(\varepsilon^{p}),$

the functions $f$ and $g$are said to be approximately equal (with an error
$o(\varepsilon^{p})$) and written as

(5) $\displaystyle f(x,\varepsilon)=g(x,\varepsilon)+o(\varepsilon^{p}),$

or or, briefly $f\approx g$ when there is no ambiguity. The approximate
equality defines an equivalence relation, and we join functions into
equivalence classes by letting $f(x,\varepsilon)$ and $g(x,\varepsilon)$ to be
members of the same class if and only if $f\approx g$. Given a function
$f(x,\varepsilon)$, let

(6) $\displaystyle f_{o}(x)+\varepsilon
f_{l}(x)+\cdots+\varepsilon^{p}f_{p}(x)$

be the approximating polynomial of degree p in $\varepsilon$ obtained via the
Taylor series expansion of $f(x,\varepsilon)$ in powers of $\varepsilon$ about
$\varepsilon=0$. Then any function $g\approx f$ (in particular, the function
$f$ itself) has the form

(7) $\displaystyle g(x,\varepsilon)=f_{o}(x)+\varepsilon
f_{l}(x)+\cdots+\varepsilon^{p}f_{p}(x)+o(\varepsilon^{p}).$

Consequently the expression (6) is called a canonical representative of the
equivalence class of functions containing $f$· Thus, the equivalence class of
functions $g(x,\varepsilon)f(x,\varepsilon)$ is determined by the ordered set
of $p+1$ functions $f_{0}(x),f_{l}(x),\cdots,f_{p}(x)$. In the theory of
approximate transformation groups, one considers ordered sets of smooth
vector-functions depending on $x$’s and a group parameter $a$:

(8) $\displaystyle f_{0}(x,a),f_{l}(x,a),\cdots,f_{p}(x,a),$

with coordinates

(9) $\displaystyle
f_{0}^{i}(x,a),f_{1}^{i}(x,a),...,f_{p}^{i}(x,a),\,\,i=1,...,n.$

Let us define the one-parameter family G of approximate transformations

(10) $\displaystyle\bar{x}^{i}\approx f_{0}^{i}(x,a)+\varepsilon
f_{1}^{i}(x,a)+\cdots+\varepsilon^{p}f_{p}^{i}(x,a),\,\,i=1,\cdots,n$

of points $x=(x^{1},\cdots,x^{n})\in\mathbb{R}^{n}$ into points
$\bar{x}=(\bar{x}^{1},\cdots,\bar{x}^{n})\in\mathbb{R}^{n}$ as the class of
invertible transformations

(11) $\displaystyle\bar{x}=f(x,a,\epsilon),$

with vector-functions $f=(f^{1},\cdots,f^{n})$ such that

(12) $\displaystyle f^{i}(x,a,\epsilon)\approx f_{0}^{i}(x,a)+\epsilon
f_{1}^{i}(x,a)+\cdots+\varepsilon^{p}f_{p}^{i}~{}(x,a),\,\,i=1,\cdots,n.$

Here $a$ is a real parameter, and the following condition is imposed:

(13) $\displaystyle f(x,0,\epsilon)\approx x.$

#### Definition

The set of transformations (10) is called a one-parameter approximate
transformation group if

(14) $\displaystyle f(f(x,a,\varepsilon),b,\epsilon)\approx
f(x,a+b,\varepsilon)$

for all transformations (11).

#### Definition

Let $G$ be a one-parameter approximate transformation group:

(15) $\displaystyle\bar{z}^{i}\approx f(z,a,\varepsilon)\equiv
f_{0}^{i}(z,a)+\varepsilon f_{1}^{i}(z,a),\,\,i=1,\cdots,N.$

An approximate equation

(16) $\displaystyle F(z,\varepsilon)\equiv F_{0}(z)+\varepsilon
F_{1}(z)\approx 0$

is said to be approximately invariant with respect to $G$, or admits $G$ if

(17) $\displaystyle F(\bar{z},\varepsilon)\approx
F(f(z,a,\varepsilon),\varepsilon)=o(\varepsilon)$

whenever $z=(z^{l},\cdots,z^{N})$ satisfies Eq. (16). If
$z=(x,u,u_{(1)},\cdots,u_{(k)})$ then (16) becomes an approximate differential
equation of order $k$, and $G$ is an approximate symmetry group of the
differential equation.

#### Theorem

Eq. (16) is approximately invariant under the approximate transformation group
(15) with the generator

(18) $\displaystyle X=X_{0}+\varepsilon
X_{1}\equiv\xi_{0}^{i}(z)\frac{\partial}{\partial
z^{i}}+\varepsilon\xi_{1}^{i}\frac{\partial}{\partial z^{i}},$

if and only if

(19) $\displaystyle[X^{(k)}F(z,\varepsilon)]_{F\approx 0}=o(\varepsilon),$

or

(20)
$\displaystyle[X_{0}^{(k)}F_{0}(z)+\varepsilon(X_{1}^{(k)}F_{0}(z)+X_{0}^{(k)}F_{1}(z))]_{(2.5)}=o(\varepsilon).$

where $X^{(k)}$ is the prolongation of $X$ of order $k$. The operator (18)
satisfying Eq. (20) is called an infinitesimal approximate symmetry of, or an
approximate operator admitted by Eq. (16). Accordingly, Eq. (20) is termed the
determining equation for approximate symmetries.

#### Theorem

If Eq. (16) admits an approximate tramformation group with the generator
$X=X_{0}+\varepsilon X_{1}$, where $X_{0}\neq 0$, then the operator

(21) $\displaystyle X_{0}=\xi_{0}^{i}(z)\frac{\partial}{\partial z^{i}}$

is an exact symmetry ofthe equation

(22) $\displaystyle F_{0}(z)=0.$

#### Definition

Eqs. (22) and (16) are termed an unperturbed equation and a perturbed
equation, respectively. Under the conditions of Theorem 2.3, the operator
$X_{0}$ is called a stable symmetry of the unperturbed equation (22). The
corresponding approximate symmetry generator $X=X_{0}+\varepsilon X_{1}$ for
the perturbed equation (16) is called a deformation of the infinitesimal
symmetry $X_{0}$ of Eq. (22) caused by the perturbation $\varepsilon
F_{1}(z)$. In particular, if the most general symmetry Lie algebra of Eq. (22)
is stable, we say that the perturbed equation (16) inherits the symmetries of
the unperturbed equation.

##  3\. Approximate symmetries of the perturbed Harry Dym equation

Consider the perturbed Harry Dym equation

$\displaystyle u_{t}+\frac{1}{2}u^{3}u_{xxx}+\varepsilon u_{x}=0.$

By applying the method of approximate transformation groups, we provide the
infinitesimal approximate symmetries (18) for the perturbed Harry Dym equation
(2).

###  3.1. Exact symmetries

Let us consider the approximate group generators in the form

(23) $\displaystyle X=X_{0}+\varepsilon
X_{1}=(\xi_{0}+\varepsilon\xi_{1})\frac{\partial}{\partial
x}+(\tau_{0}+\varepsilon\tau_{1})\frac{\partial}{\partial
t}+(\phi_{0}+\varepsilon\phi_{1})\frac{\partial}{\partial u}$

where $r\xi_{i},\tau_{i}$ and $\phi_{i}$ for $i=0,1$ are unknown functions
of$x,t$ and $u$. Solving the determining equation

(24) $\displaystyle
X_{0}^{(3)}(u_{t}-\frac{1}{2}u^{3}u_{xxx})\mid_{u_{t}-\frac{1}{2}u^{3}u_{xxx}=0}=0,$

for the exact symmetries $X_{0}$ of the unperturbed equation, we obtain

$\displaystyle\xi_{0}$ $\displaystyle=$
$\displaystyle(A_{1}+A_{2}x+\frac{A_{3}}{2}x^{2}),$ (25)
$\displaystyle\tau_{0}$ $\displaystyle=$ $\displaystyle(A_{4}+3A_{5}t),$
$\displaystyle\phi_{0}$ $\displaystyle=$
$\displaystyle(A_{2}-\frac{1}{3A_{5}}+xA_{3})u,$

where $A_{1},\cdots,A_{5}$ are arbitrary constants. Hence,

(26) $\displaystyle
X_{0}=(A_{1}+A_{2}x+\frac{A_{3}}{2}x^{2})\frac{\partial}{\partial
x}+(A_{4}+3A_{5}t)\frac{\partial}{\partial
t}+(A_{2}-\frac{1}{3A_{5}}+xA_{3})u)\frac{\partial}{\partial u}$

Therefore, the unperturbed Harry Dym equation, admits the five-dimensional Lie
algebra with the basis

(32) $\displaystyle\begin{array}[]{l}\displaystyle
X_{0}^{1}=\frac{\partial}{\partial x},\\\\[5.69054pt] \displaystyle
X_{0}^{2}=\frac{\partial}{\partial t},\\\\[5.69054pt] \displaystyle
X_{0}^{3}=x\frac{\partial}{\partial x}+u\frac{\partial}{\partial
u},\end{array}\qquad\begin{array}[]{l}\displaystyle
X_{0}^{4}=3t\frac{\partial}{\partial t}-u\frac{\partial}{\partial
u},\\\\[5.69054pt] \displaystyle X_{0}^{5}=2x^{2}\frac{\partial}{\partial
x}+xu\frac{\partial}{\partial u}.\end{array}$

###  3.2. Approximate symmetries

At first, we need to determine the auxiliary finction $H$ by vitue of Eqs.
(19), (20) and (16), i.e., by the equation

(33) $\displaystyle H=\frac{1}{\varepsilon}[X_{0}^{(k)}(F_{0}(z)+\varepsilon
F_{1}(z))\mid_{F_{0}(z)+\varepsilon F_{1}(z)=0}]$

Substituting the expression (26) of the generator $X_{0}$ into Eq. (33) we
obtain the auxiliary function

(34) $\displaystyle H=u_{x}(A_{5}-A_{2})+A_{3}(u-xu_{x})$

Now, calculate the operators $X_{1}$ by solving the inhomogeneous determining
equation for deformations:

(35) $\displaystyle X_{1}^{(k)}F_{0}(z)\mid_{F_{0}(z)=0}+H=0.$

So, the above determinig equation for this equation is written as

(36) $\displaystyle
X_{1}^{(3)}(u_{t}+\frac{1}{2}u^{3}u_{xxx})\mid_{u_{t}+\frac{1}{2}u^{3}u_{xxx}=0}+u_{x}(A_{5}-A_{2})+A_{3}(u-xu_{x})=0,$

solving the determining equation yields,

$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle(A_{5}-A_{2})t-A_{3}xt+C_{4}x-C_{5}+\frac{C_{3}}{2}x^{2},$ (37)
$\displaystyle\tau_{1}$ $\displaystyle=$ $\displaystyle(C_{1}t+C_{2}),$
$\displaystyle\phi_{1}$ $\displaystyle=$
$\displaystyle(-A_{3}t+C_{4}+C_{3}x+\frac{C_{1}}{3})u,$

where $C_{1},\cdots,C_{5}$ are arbitrary constants.

Thus, we derive the following approximate symmetries of the perturvbed Harry
Dym equation:

(48)
$\displaystyle\begin{array}[]{l}\displaystyle\mathbf{v_{1}}=\frac{\partial}{\partial
x},\\\\[5.69054pt] \displaystyle\mathbf{v}_{2}=\frac{\partial}{\partial
t},\\\\[5.69054pt] \displaystyle\mathbf{v}_{3}=x\frac{\partial}{\partial
x}+u\frac{\partial}{\partial u},\\\\[5.69054pt]
\displaystyle\mathbf{v}_{4}=3t\frac{\partial}{\partial
t}-u\frac{\partial}{\partial
u},\end{array}\qquad\begin{array}[]{l}\displaystyle\mathbf{v}_{5}=2x^{2}\frac{\partial}{\partial
x}+xu\frac{\partial}{\partial u},\\\\[5.69054pt]
\displaystyle\mathbf{v}_{6}=\varepsilon\frac{\partial}{\partial
x},\\\\[5.69054pt]
\displaystyle\mathbf{v}_{7}=\varepsilon\frac{\partial}{\partial
t},\\\\[5.69054pt]
\displaystyle\mathbf{v}_{8}=\varepsilon(x\frac{\partial}{\partial
x}+u\frac{\partial}{\partial
u}),\end{array}\qquad\begin{array}[]{l}\displaystyle\mathbf{v}_{9}=\varepsilon(3t\frac{\partial}{\partial
t}-u\frac{\partial}{\partial u}),\\\\[5.69054pt]
\displaystyle\mathbf{v}_{10}=\varepsilon(2x^{2}\frac{\partial}{\partial
x}+xu\frac{\partial}{\partial u}).\end{array}$

The following table of commutators, evaluated in the first-order of precision,
shows that the operators (48) span an ten-dimensional approximate Lie algebra
, and hence generate an ten-parameter approximate transformations group.

Table 1. Approximate Commutators of approximate symmetry of perturbed Harry
Dym equation
$\displaystyle\begin{array}[]{c|cccccccccc}&\mathbf{v_{1}}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}&\mathbf{v}_{5}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\hline\cr&&&&&&\\\\[-8.53581pt]
\mathbf{v_{1}}&0&0&\mathbf{v_{1}}&0&2\mathbf{v}_{3}&0&0&\mathbf{v}_{6}&0&2\mathbf{v}_{8}\\\
\mathbf{v}_{2}&0&0&0&12\mathbf{v}_{2}&0&0&0&0&3\mathbf{v}_{7}&0\\\
\mathbf{v}_{3}&-\mathbf{v}_{1}&0&0&0&\mathbf{v}_{5}&-\mathbf{v}_{6}&0&0&0&\mathbf{v}_{10}\\\
\mathbf{v}_{4}&0&-12\mathbf{v}_{2}&0&0&0&0&-3\mathbf{v}_{7}&0&0&0\\\
\mathbf{v}_{5}&-2\mathbf{v}_{3}&0&-\mathbf{v}_{5}&0&0&-2\mathbf{v}_{8}&0&-\mathbf{v}_{10}&0&0\\\
\mathbf{v}_{6}&0&0&\mathbf{v}_{6}&0&2\mathbf{v}_{8}&0&0&0&0&0\\\
\mathbf{v}_{7}&0&0&0&3\mathbf{v}_{7}&0&0&0&0&0&0\\\
\mathbf{v}_{8}&-\mathbf{v}_{6}&0&0&0&\mathbf{v}_{10}&0&0&0&0&0\\\
\mathbf{v}_{9}&0&-3\mathbf{v}_{7}&0&0&0&0&0&0&0&0\\\
\mathbf{v}_{10}&-2\mathbf{v}_{8}&0&-\mathbf{v}_{10}&0&0&0&0&0&0&0\end{array}$

#### Remark.

Equations (48) show that all symmetries (32) of Eq. (1) are stable. Hence, the
perturbed equation (2) inherits the symmetries of the unperturbed equation
(1).

##  4\. The structure of the Lie algebra of symmetries

In this section, we determine the structure of the Lie algebra of symmetries
of the perturbed Harry Dym equation. The Lie algebra $\mathfrak{g}$ is non-
solvable, since

$\displaystyle\mathfrak{g}^{(1)}$ $\displaystyle=$
$\displaystyle[\mathfrak{g},\mathfrak{g}]=\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{5},\mathbf{v}_{6},\mathbf{v}_{7},\mathbf{v}_{8},\mathbf{v}_{10}\\}$
(50) $\displaystyle\mathfrak{g}^{(2)}$ $\displaystyle=$
$\displaystyle[\mathfrak{g}^{(1)},\mathfrak{g}^{(1)}]=\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{1},\mathbf{v}_{3},\mathbf{v}_{5},\mathbf{v}_{6},\mathbf{v}_{8},\mathbf{v}_{10}\\}$
$\displaystyle\mathfrak{g}^{(3)}$ $\displaystyle=$
$\displaystyle[\mathfrak{g}^{(2)},\mathfrak{g}^{(2)}]=\mathfrak{g}^{(2)}$

The Lie algebra $\mathfrak{g}$ admits a Levi decomposition as the following
semi-direct product $\mathfrak{g}=r\propto s$, where

(51) $\displaystyle
r=\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{2},\mathbf{v}_{4},\mathbf{v}_{6},\mathbf{v}_{7},\mathbf{v}_{8},\mathbf{v}_{9},\mathbf{v}_{10}\\}$

is the radical of $\mathfrak{g}$ (the largest solvable ideal contained in
$\mathfrak{g}$) and

(52) $\displaystyle
s=\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{1},\mathbf{v}_{3},\mathbf{v}_{5}\\}$

is a semi-simple subalgebra of $\mathfrak{g}$.

The radical $r$ is solvable with the following chain of ideals

(53) $\displaystyle r^{(1)}\supset r^{(2)}\supset r^{(3)}=\\{0\\},$

where

$\displaystyle r^{(1)}$ $\displaystyle=$
$\displaystyle\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{2},\mathbf{v}_{4},\mathbf{v}_{6},\mathbf{v}_{7},\mathbf{v}_{8},\mathbf{v}_{9},\mathbf{v}_{10}\\},$
(54) $\displaystyle r^{(2)}$ $\displaystyle=$
$\displaystyle\mathrm{Span}_{\mathbb{R}}\\{\mathbf{v}_{2},\mathbf{v}_{7}\\}.$

The semi-simple subalgebra $s$ of $\mathfrak{g}$ is isomorphic to the Lie
algebra $A_{3,8}$ of the classification of three dimensional Lie algebras in
[13], by the following isomorphism.

(55)
$\displaystyle\mathfrak{T}:\\{\mathbf{v}_{1},\mathbf{v}_{3},\mathbf{v}_{5}\\}\rightarrow\\{\mathbf{v}_{1},-\mathbf{v}_{2},-\mathbf{v}_{3}\\}$

##  5\. Optimal system for perturbed Harry Dym equation

#### Definition

Let $G$ be a Lie group. An optimal system of $s$-parameter subgroups is a list
of conjugacy inequivalent $s$-parameter subgroups with the property that any
other subgroup is conjugate to precisely one subgroup in the list. Similarly,
a list of $s$-parameter subalgebras forms an optimal system if every
$s$-parameter subalgebra of $\mathfrak{g}$ is equivalent to a unique member of
the list under some element of the adjoint representation:
$\tilde{\mathfrak{h}}=\mathrm{Ad}g(\mathfrak{h})),\,g\in G$.

#### Proposition

Let $H$ and $\tilde{H}$ be connected, $s$-dimensional Lie subgroups of the Lie
group $G$ with corresponding Lie subalgebras $\mathfrak{h}$ and
$\tilde{\mathfrak{h}}$ of the Lie algebra $\mathfrak{g}$ of $G$. Then
$\tilde{H}=gHg^{-1}$ are conjugate subgroups if and only if
$\tilde{\mathfrak{h}}=\mathrm{Ad}g(\mathfrak{h})$are conjugate
subalgebras.(Proposition 3.7 of [11])

Actually, the proposition says that the problem of finding an optimal system
of subgroups is equivalent to that of finding an optimal system of
subalgebras. For one-dimensional subalgebras, this classification problem is
essentially the same as the problem of classifying the orbits of the adjoint
representation, since each one-dimensional subalgebra is determined by a
nonzero vector in $\mathfrak{g}$. To compute the adjoint representation we use
the Lie series

(56)
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\mu\mathbf{v}_{i}))\mathbf{v}_{j}=\mathbf{v}_{j}-\mu[\mathbf{v}_{i},\mathbf{v}_{j}]+\frac{\mu^{2}}{2}[[\mathbf{v}_{i},[\mathbf{v}_{i},\mathbf{v}_{j}]]-\cdots.$

where $[\mathbf{v}_{i},\mathbf{v}_{j}],\,i,j=1,\cdots,10$ is the commutator
for the Lie algebra and $\mu$ is a parameter. In this manner, we construct the
table with the $(i,j)$-th entry indicating
$\mathrm{Ad}(\mathrm{exp}(\mu\mathbf{v}_{i}))\mathbf{v}_{j}$.

Table 2. Adjoint representation of approximate symmetry of the perturbed Harry
Dym equation

$\displaystyle\begin{array}[]{c|cccccccccc}\mathrm{Ad}&\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}&\mathbf{v}_{5}\\\
\hline\cr&&&&&&\\\\[-8.53581pt]
\mathbf{v}_{1}&\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}-\mu\mathbf{v}_{1}&\mathbf{v}_{4}&\mathbf{v}_{5}-2\mu\mathbf{v}_{3}+\mu^{2}\mathbf{v}_{1}\\\
\mathbf{v}_{2}&\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}-12\mu\mathbf{v}_{2}&\mathbf{v}_{5}\\\
\mathbf{v}_{3}&e^{\mu}\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}&e^{-\mu}\mathbf{v}_{5}\\\
\mathbf{v}_{4}&\mathbf{v}_{1}&e^{12\mu}\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}&\mathbf{v}_{5}\\\
\mathbf{v}_{5}&\mathbf{v}_{1}+2\mu\mathbf{v}_{3}+\mu^{2}\mathbf{v}_{5}&\mathbf{v}_{2}&\mathbf{v}_{3}+\mu\mathbf{v}_{5}&\mathbf{v}_{4}&\mathbf{v}_{5}\\\
\mathbf{v}_{6}&\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}-\mu\mathbf{v}_{6}&\mathbf{v}_{4}&\mathbf{v}_{5}-2\mu\mathbf{v}_{8}\\\
\mathbf{v}_{7}&\mathbf{v}_{1}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}-3\mu\mathbf{v}_{7}&\mathbf{v}_{5}\\\
\mathbf{v}_{8}&\mathbf{v}_{1}+\mu\mathbf{v}_{6}&\mathbf{v}_{2}&\mathbf{v}_{3}&\mathbf{v}_{4}&\mathbf{v}_{5}-\mu\mathbf{v}_{10}\\\
\mathbf{v}_{9}&\mathbf{v}_{1}&\mathbf{v}_{2}+3\mu\mathbf{v}_{7}&\mathbf{v}_{3}&\mathbf{v}_{4}&\mathbf{v}_{5}\\\
\mathbf{v}_{10}&\mathbf{v}_{1}+2\mu\mathbf{v}_{8}&\mathbf{v}_{2}&\mathbf{v}_{3}+\mu\mathbf{v}_{10}&\mathbf{v}_{4}&\mathbf{v}_{5}\end{array}$

$\displaystyle\begin{array}[]{c|cccccccccc}\mathrm{Ad}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\hline\cr&&&&&&\\\\[-8.53581pt]
\mathbf{v}_{1}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}-\mu\mathbf{v}_{6}&\mathbf{v}_{9}&\mathbf{v}_{10}-2\mu\mathbf{v}_{8}+\mu^{2}\mathbf{v}_{6}\\\
\mathbf{v}_{2}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}-3\mu\mathbf{v}_{7}&\mathbf{v}_{10}\\\
\mathbf{v}_{3}&e^{\mu}\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&e^{-\mu}\mathbf{v}_{10}\\\
\mathbf{v}_{4}&\mathbf{v}_{6}&e^{3\mu}\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{5}&\mathbf{v}_{6}+2\mu\mathbf{v}_{8}+\mu^{2}\mathbf{v}_{10}&\mathbf{v}_{7}&\mathbf{v}_{8}+\mu+\mathbf{v}_{10}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{6}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{7}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{8}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{9}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\\\
\mathbf{v}_{10}&\mathbf{v}_{6}&\mathbf{v}_{7}&\mathbf{v}_{8}&\mathbf{v}_{9}&\mathbf{v}_{10}\end{array}$

#### Theorem

An optimal system of one-dimensional approximate Lie algebras of the perturbed
Harry Dym equation is provided by

Table 3. One-dimensional optimal system of the perturbed Harry Dym equation
$\displaystyle\begin{array}[]{l}\mathbf{v}^{1}=\mathbf{v}_{8},\\\
\mathbf{v}^{2}=\mathbf{v}_{7}+a\mathbf{v}_{8},\\\
\mathbf{v}^{3}=\mathbf{v}_{6}+\mathbf{v}_{8},\end{array}\qquad\begin{array}[]{l}\mathbf{v}^{4}=\mathbf{v}_{6}-\mathbf{v}_{7}+\mathbf{v}_{8},\\\
\mathbf{v}^{5}=\mathbf{v}_{6}+\mathbf{v}_{7}+\mathbf{v}_{8},\\\
\mathbf{v}^{6}=\mathbf{v}_{2}+a\mathbf{v}_{8},\end{array}\qquad\begin{array}[]{l}\mathbf{v}^{7}=\mathbf{v}_{2}-\mathbf{v}_{6}+a\mathbf{v}_{8},\\\
\mathbf{v}^{8}=\mathbf{v}_{2}+\mathbf{v}_{6}+a\mathbf{v}_{8},\\\
\mathbf{v}^{9}=\mathbf{v}_{1}+a\mathbf{v}_{2}+b\mathbf{v}_{7}\end{array}$
$\displaystyle\mathbf{v}^{10}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+b\mathbf{v}_{2}+\mathbf{v}_{5}+c\mathbf{v}_{6}+d\mathbf{v}_{7},$
$\displaystyle\mathbf{v}^{11}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+b\mathbf{v}_{2}+\mathbf{v}_{3}+c\mathbf{v}_{5}+d\mathbf{v}_{7}+e\mathbf{v}_{8},$
$\displaystyle\mathbf{v}^{12}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+b\mathbf{v}_{3}+\mathbf{v}_{4}+c\mathbf{v}_{5}+d\mathbf{v}_{6}+e\mathbf{v}_{8},$
$\displaystyle\mathbf{v}^{13}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+b\mathbf{v}_{3}+c\mathbf{v}_{4}+d\mathbf{v}_{5}+e\mathbf{v}_{6}+f\mathbf{v}_{8}+\mathbf{v}_{9},$
$\displaystyle\mathbf{v}^{14}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}-\mathbf{v}_{2}+b\mathbf{v}_{3}+c\mathbf{v}_{4}+d\mathbf{v}_{5}+e\mathbf{v}_{6}+f\mathbf{v}_{8}+\mathbf{v}_{9},$
$\displaystyle\mathbf{v}^{15}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+\mathbf{v}_{2}+b\mathbf{v}_{3}+c\mathbf{v}_{4}+d\mathbf{v}_{5}+e\mathbf{v}_{6}+f\mathbf{v}_{8}+\mathbf{v}_{9},$
$\displaystyle\mathbf{v}^{16}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!a\mathbf{v}_{1}+\mathbf{v}_{2}+b\mathbf{v}_{3}+c\mathbf{v}_{4}+d\mathbf{v}_{5}+e\mathbf{v}_{6}+f\mathbf{v}_{8}+\mathbf{v}_{9}.$

Proof. Consider the approximate symmetry algebra $\mathfrak{g}$ of the
unperturbed Harry Dym equation, whose adjoint representation was determined in
the table 2. Given a nonzero vector

$\displaystyle\mathbf{v}=\sum_{i=1}^{10}a_{i}\mathbf{v_{i}},$

our task is to simplify as many of the coefficients ai as possible through
judicious applications of adjoint maps to $\mathbf{v}$.

Suppose first that $a_{10}\neq 0$. Scaling $\mathbf{v}$ if necessary, we can
assume that $a_{10}=1$. Referring to table (3.24), if we act on such a
$\mathbf{v}$ by

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\frac{a_{8}}{2}\mathbf{v}_{8}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}^{{}^{\prime}}\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}+a_{3}^{{}^{\prime}}\mathbf{v}_{3}+a_{4}\mathbf{v}_{4}+a_{5}\mathbf{v}_{5}+a_{6}^{{}^{\prime}}\mathbf{v}_{6}+a_{7}\mathbf{v}_{7}+a_{9}\mathbf{v}_{9}+\mathbf{v}_{10},$

we can make the coefficient of $a_{8}$ vanish. The remaining one-dimensional
subalgebras are spanned by vectors of the above form with $a_{10}=0$. If
$a_{9}\neq 0$, we scale to make $a_{9}=1$, and then act on $\mathbf{v}$ to
cancel the coefficient of $a_{7}$ as follows:

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\frac{a_{7}}{3}\mathbf{v}_{2}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}^{{}^{\prime}}\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}+a_{3}^{{}^{\prime}}\mathbf{v}_{3}+a_{4}\mathbf{v}_{4}+a_{5}\mathbf{v}_{5}+a_{6}\mathbf{v}_{6}+a_{8}\mathbf{v}_{8}+\mathbf{v}_{9}.$

We can further act on $\mathbf{v}^{\prime}$ by the group generated by
$\mathbf{v}_{4}$; this has the net effect of scaling the coefficients of
$\mathbf{v}_{2}$:

$\displaystyle\mathbf{v}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\mu\mathbf{v}_{4}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}\mathbf{v}_{1}+e^{12\mu}a_{2}\mathbf{v}_{2}+a_{3}\mathbf{v}_{3}+a_{4}\mathbf{v}_{4}+a_{5}\mathbf{v}_{5}+a_{6}\mathbf{v}_{6}+a_{8}\mathbf{v}_{8}+\mathbf{v}_{9}.$

So, depending on the sign of $a_{2}$, we can make the coefficient of
$\mathbf{v}_{2}$ either $+1,-1$ or $0$. If $a_{10}=a_{9}=0$ and $a_{4}\neq 0$,
we scale to make $a_{4}=1$. So, the non-zero vector $\mathbf{v}$ is equivalent
to $\mathbf{v}^{\prime}$ under adjoint maps:

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\frac{a_{7}}{3}\mathbf{v}_{7}))\circ\mathrm{Ad}(\mathrm{exp}(\frac{a_{2}}{12}\mathbf{v}_{2}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}\mathbf{v}_{1}+a_{3}\mathbf{v}_{3}+\mathbf{v}_{4}+a_{5}\mathbf{v}_{5}+a_{6}\mathbf{v}_{6}+a_{8}\mathbf{v}_{8}$

If $a_{10}=a_{9}=a_{4}=0$ and $a_{3}\neq 0$, by scaling $\mathbf{v}$, we can
assume that $a_{3}=1$. Referring to the table, if we act on such a
$\mathbf{v}$ by the following adjoint map, we can arrange that the
coefficients of $a_{6}$ vanish.

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(a_{6}\mathbf{v}_{6}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}+\mathbf{v}_{3}+a_{5}\mathbf{v}_{5}+a_{6}\mathbf{v}_{6}+a_{7}\mathbf{v}_{7}+a_{8}^{{}^{\prime}}\mathbf{v}_{8}.$

If $a_{10}=a_{9}=a_{4}=a_{3}=0$ and $a_{5}\neq 0$, we scale to make $a_{5}=1$.
Thus, $\mathbf{v}$ is equivalent to $\mathbf{v}^{\prime}$ under the adjoint
representations:

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\frac{a_{8}}{2}\mathbf{v}_{6}))\mathbf{v}$
$\displaystyle=$ $\displaystyle
a_{1}\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}+\mathbf{v}_{5}+a_{6}\mathbf{v}_{6}+a_{7}\mathbf{v}_{7}.$

If $a_{10}=a_{9}=a_{4}=a_{3}=a_{5}=0$ and $a_{1}\neq 0$, we scale to make
$a_{1}=1$. So, we can make the coefficients of $a_{6},a_{8}$ zero by using the
following adjoint maps:

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(-\frac{a_{8}}{2}\mathbf{v}_{10}))\circ\mathrm{Ad}(\mathrm{exp}(-a_{6}\mathbf{v}_{8}))\mathbf{v}$
$\displaystyle=$
$\displaystyle\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}+a_{7}\mathbf{v}_{7}.$

If $a_{10}=a_{9}=a_{4}=a_{3}=a_{5}=a_{1}=0$ and $a_{2}\neq 0$, by scaling
$\mathbf{v}$, we can assume that $a_{2}=1$. Therefore, we can arrange that the
coefficients of $a_{7}$ vanish by simplifying the non-zero vector $\mathbf{v}$
as follows:

$\displaystyle\mathbf{v}^{\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(-\frac{a_{7}}{3}\mathbf{v}_{9}))\mathbf{v}$
$\displaystyle=$
$\displaystyle\mathbf{v}_{2}+a_{6}\mathbf{v}_{6}+a_{8}\mathbf{v}_{8},$

We can further act on $\mathbf{v}^{\prime}$ by the group generated by
$\mathbf{v}_{3}$;

$\displaystyle\mathbf{v}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\mathrm{Ad}(\mathrm{exp}(\mu\mathbf{v}_{3})\mathbf{v}$
$\displaystyle=$
$\displaystyle\mathbf{v}_{2}+e^{\mu}a_{6}\mathbf{v}_{6}+a_{8}\mathbf{v}_{8}.$

So, depending on the sign of $a_{6}$, we can make the coefficient of
$\mathbf{v}_{6}$ either $+1,-1$ or $0$. If
$a_{10}=a_{9}=a_{4}=a_{3}=a_{5}=a_{1}=a_{2}=0$ and $a_{6}\neq 0$, by scaling
$\mathbf{v}$, we can assume that $a_{6}=1$. We can act on such a $\mathbf{v}$
by the group generated by $\mathbf{v}_{4}$; So, depending on the sign of
$a_{7}$, we can make the coefficient of $\mathbf{v}_{7}$ either $+1,-1$ or
$0$. The remaining cases $a_{10}=a_{9}=a_{4}=a_{3}=a_{5}=a_{2}=a_{1}=a_{3}=0$
and $a_{7}\neq 0$, no further simplifications are possible. The last remaining
cace occurs when
$a_{10}=a_{9}=a_{4}=a_{3}=a_{5}=a_{1}=a_{2}=a_{4}=a_{6}=a_{7}=0$ and
$a_{8}\neq 0$, for which our earlier simplifications were unnecessary. Since,
the only remaining vectors are the multiples of $\mathbf{v}_{8}$, on which the
adjoint representation acts trivially.

##  6\. Approximately differential invariants for the perturbed Harry Dym
equation

In this section, we compute some approximately differential invariants of the
perturbed Harry Dym equation with respect to the optimal system. Consider the
operator $\mathbf{v}^{2}$. To determine the independent invariants $I$, we
need to solve the first order partial differential equation

(60) $\displaystyle(\varepsilon\frac{\partial}{\partial t}+a\varepsilon
x\frac{\partial}{\partial x}+a\varepsilon u\frac{\partial}{\partial
u})(I(x,t,u))=0,$

that is

(61) $\displaystyle\varepsilon\frac{\partial I}{\partial t}+a\varepsilon
x\frac{\partial I}{\partial x}+a\varepsilon u\frac{\partial I}{\partial u}=0,$

which is a first order homogeneous PDE. The solution can be found by
integrating the corresponding characteristic system of ordinary differential
equation, which is

(62) $\displaystyle\frac{dx}{a\varepsilon
x}=\frac{dt}{\varepsilon}=\frac{du}{a\varepsilon u}$

Hence, the independent approximately differntial invariants are as follows:

(63) $\displaystyle y=\frac{u}{x},\quad v=\frac{\ln x-at}{a}.$

In this manner, we investigate some independent approximately differential
invariants with respect to the optimal system which are listed in Table 4.

Table 4. Approximately differential invariants for the perturbed Harry Dym equation Operator |  | Approximate Differential Invariants
---|---|---
$\displaystyle\mathbf{v}_{1}$ |  | $\displaystyle t,u$
$\displaystyle\mathbf{v}_{2}$ |  | $\displaystyle x,u$
$\displaystyle\mathbf{v}_{3}$ |  | $\displaystyle t,\frac{u}{x}$
$\displaystyle\mathbf{v}_{4}$ |  | $\displaystyle x,ut^{1/3}$
$\displaystyle\mathbf{v}_{5}$ |  | $\displaystyle t,\frac{u}{x^{2}}$
$\displaystyle\mathbf{v}_{7}+a\mathbf{v}_{8}$ |  | $\displaystyle-\frac{\ln x}{a}+t,\frac{u}{x}$
$\displaystyle\mathbf{v}_{6}+\mathbf{v}_{8}$ |  | $\displaystyle t,\frac{u}{x+1}$
$\displaystyle\mathbf{v}_{6}-\mathbf{v}_{7}+\mathbf{v}_{8}$ |  | $\displaystyle\ln(x+1)+t,\frac{u}{x+1}$
$\displaystyle\mathbf{v}_{6}+\mathbf{v}_{7}+\mathbf{v}_{8}$ |  | $\displaystyle-\ln(x+1)+t,\frac{u}{x+1}$
$\displaystyle\mathbf{v}_{2}+a\mathbf{v}_{8}$ |  | $\displaystyle-\frac{\ln x}{a\varepsilon}+t,\frac{u}{x}$
$\displaystyle\mathbf{v}_{2}-\mathbf{v}_{6}+a\mathbf{v}_{8}$ |  | $\displaystyle-\frac{\ln(ax-1)}{a\varepsilon}+t,\frac{u}{ax-1}$
$\displaystyle\mathbf{v}_{2}+\mathbf{v}_{6}+a\mathbf{v}_{8}$ |  | $\displaystyle-\frac{\ln(ax+1)}{a\varepsilon}+t,\frac{u}{ax+1}$
$\displaystyle\mathbf{v}_{1}+a\mathbf{v}_{2}+b\mathbf{v}_{7}$ |  | $\displaystyle-b\varepsilon x-a\varepsilon+t,u$
$\displaystyle a\mathbf{v}_{1}+b\mathbf{v}_{2}+\mathbf{v}_{5}+c\mathbf{v}_{6}+d\mathbf{v}_{7}$ |  | $\displaystyle\frac{-d\varepsilon-b}{\sqrt{c\varepsilon+a}}\arctan\Big{(}\frac{x}{\sqrt{c\varepsilon+a}}\Big{)}+t,\frac{u}{x^{2}+c\varepsilon+a}$

##  7\. Conclusions

In this paper, we investigate the approximate symmetry of the perturbed Harry
Dym equation and discuss on the structure of its Lie algebra. Moreover, we
compute optimal system of one-dimensional approximate Lie algebras of the
perturbed Harry Dym equation and derive some approximately differential
invarints with respect to the generators of Lie algebra and optimal system.

## References

  * [1] V.A. Baikov, R.K. Gazizov and N.H. Ibragimov, Approximate symmetries of equations with a small parameter, Math. Sb. 136 (1988), 435-450 (English Transl. in Math . USSR Sb. 64 (1989), 427-441).
  * [2] V.A. Baikov, R.K. Gazizov and N.H Ibragimov, Approximate transformation groups and deformations of symmetry Lie algebras, Chapter 2, in CRC Handbook of Lie Group Analysis of Differential Equation, Vol. 3, Editor N.H. Ibragimov, Boca Raton, Florida, CRC Press, 1996.
  * [3] F. Calogero and A. Degasperis, Spectral Transform and Solitons 1, North Holland, Amsterdam, 1982.
  * [4] N. Euler, M.W. Shulga, W.H. Steeb, Approximate symmetries and approximate solutions for a multidimensional LandauGinzburg equation, J. Phys. A: Math. Gen. 25 (1992), 1095-1103.
  * [5] M. Euler, N. Euler, A. K¨ohler, On the construction of approximate solutions for a multidimensional nonlinear heat equation, J. Phys. A: Math. Gen. 27 (1994), 2083-2092.
  * [6] W.I. Fushchich, W.H. Shtelen, On approximate symmetry and approximate solutions of the non-linear wave equation with a small parameter, J. Phys. A: Math. Gen. 22 (1989), 887-890.
  * [7] W. Hereman, P. P. Banerjee and M. R. Chatterjee, Derivation and implidt solution of the Harry Dym equation and its connections with the Korteweg–de Vries equation, J. Phys. A, 22:241-255, 1989.
  * [8] N.H. Ibragimov and V.F. Kovalev, Approximate and Renormgroup Symmetries. Beijing (P.R.China): Higher Education Press, 2009. In Series: Nonlinear Physical Science, Ed. Albert C.J Luo and N.H Ibragimov.
  * [9] M. D. Kruskal, Lecture Notes in Physics, vol. 38, 310-354, Springer, Berlin, 1975.
  * [10] F. Magri, 1978, A simple model of the integrable Hamiltonian equation, J. Math. Phys. 19 1156-1162.
  * [11] P.J. Olver, Application of Lie Groups to Differential Equations, 2nd ed., Springer-Verlag, New York, 1993.
  * [12] M. Pakdemirli, M. Y¨ur¨usoy, ˙I.T. Dolap¸ci, Comparison of approximate symmetry methods for differential equations, Acta. Appl. Math. 80 (2004) 243-271.
  * [13] Patera, J., Sharp, R. T., and Winternitz, P., Invariants of real low dimensional Lie algebras, JMP vol 17, No 6, June 1976, 966–994.
  * [14] R.J. Wiltshire, Two approaches to the calculation of approximate symmetry exemplified using a system of advection-diffusion equations, J. Comput. Appl. Math., 197 (2006) 287-301.
  * [15] P. C. Sabatier, On some spectral problems and isospectral evolutions connected with the classical string problem II, Lett. Nuovo Cimento 26:477-483, 1979.
  * [16] L. Yi-Shen, Evolution Equations Associated with the Eigenvalue Problem Based on the Equation $\varphi_{xx}=[\upsilon(x)-\kappa^{2}\rho^{2}(x)]\varphi$, Lett.Nuovo Cimento 70 B:1-12, 1982.
  * [17] M. Wadati, Y. H. Ichikawa and T. Shimizu, Cusp Soliton of a New Integrable Nonlinear Evolution Equation, Prog. Theor. Phys. 64:1959-1967, 1980.
  * [18] M. Wadati, K. Konno and Y. H. Ichikawa, New Integrable Nonlinear Evolution Equations, J. Phys. Soc. Japan 47: 1698-1700, 1979.

Mehdi Nadjafikhah

Iran University of Science and Technology, Narmak-16, Tehran, Iran
E-mail<EMAIL_ADDRESS>
Parastoo Kabi-Nejad

Iran University of Science and Technology, Narmak-16, Tehran, Iran
E-mail<EMAIL_ADDRESS>


UMPA, Ecole Normale Supérieure de Lyon
Email<EMAIL_ADDRESS>

# Velocity Syzygies and Bounding Syzygy Moments in the Planar Three-Body
Problem

Alexei Tsygvintsev

###### Abstract.

We consider the Newtonian planar three-body problem, defining a syzygy
(velocity syzygy) as a configuration where the positions (velocities) of the
three bodies become collinear. We demonstrate that if the total energy is
negative, every collision-free solution has an infinite number of velocity
syzygies. Specifically, the velocities of the three bodies become parallel
within every interval of time containing three consecutive syzygies. Using
comparison theory for matrix Riccati equations, we derive new upper and lower
bounds on the moments when syzygies occur.

###### Key words and phrases:

dynamical systems, celestial mechanics, three-body problem, syzygies

##  1\. Introduction and Preliminary Results

The Newtonian planar three-body problem has been a subject of extensive study
in celestial mechanics. It involves predicting the motion of three bodies
moving under their mutual gravitational attraction. Understanding syzygies,
where the positions or velocities of the three bodies become collinear,
provides significant insight into the dynamics of the system.

Let $P_{1}$, $P_{2}$, and $P_{3}$ be three points in the plane with strictly
positive masses $m_{1},m_{2},m_{3}$ and Euclidean coordinates
$(x_{i},y_{i})\in\mathbb{R}^{2}$, $i=1,2,3$. The Newtonian three-body problem
[2] can be formulated as follows:

$\ddot{z}_{1}=m_{2}\frac{z_{21}}{|z_{21}|^{3}}-m_{3}\frac{z_{13}}{|z_{13}|^{3}},\quad\ddot{z}_{2}=m_{3}\frac{z_{32}}{|z_{32}|^{3}}-m_{1}\frac{z_{21}}{|z_{21}|^{3}},\quad\ddot{z}_{3}=m_{1}\frac{z_{13}}{|z_{13}|^{3}}-m_{2}\frac{z_{32}}{|z_{32}|^{3}}\,,$
(1.1)

where $z_{k}=x_{k}+iy_{k}\in\mathbb{C}$, $k=1,2,3$ and $z_{kl}=z_{k}-z_{l}$.

Given that the total linear momentum (which is conserved) is zero, one can
always set

$\sum\limits_{k}m_{k}\dot{z}_{k}=0,\quad\sum\limits_{k}m_{k}z_{k}=0\,,$ (1.2)

by placing the center of mass of the bodies at the origin.

Let $t\mapsto z_{i}(t)$, $i=1,2,3$ be any collision-free solution of equations
(1.1) defined for $t\in I=[0,a)$, $a>0$, and determined by a set of initial
conditions $(z_{i}(0),\dot{z}_{i}(0))$, $i=1,2,3$.

######  Definition 1.1.

We say that the three bodies $P_{1},P_{2},P_{3}$ form a syzygy (velocity
syzygy) at the moment $t_{0}\in I$ if the complex triplet
$(z_{1},z_{2},z_{3})(t_{0})$ ($(\dot{z}_{1},\dot{z}_{2},\dot{z}_{3})(t_{0})$)
is collinear.

Our first result concerns the existence of velocity syzygies.

######  Theorem 1.1.

Let $t\mapsto(z_{1}(t),z_{2}(t),z_{3}(t))$, $t\in[0,+\infty)$, be a zero
angular momentum, collision-free solution to the three-body problem (1.1) with
negative energy. Then, it has an infinite number of velocity syzygies.

###### Proof.

After introducing the new variables $w_{i}=m_{i}z_{i}$, $i=1,2,3$, the
relations (1.2) yield

$\sum\limits_{i}w_{i}=\sum\limits_{i}\dot{w}_{i}=0.$ (1.3)

Writing

$w_{k}=X_{k}+iY_{k},\quad X_{k}=m_{k}x_{k},\quad Y_{k}=m_{k}y_{k},\quad
k=1,2,3\,,$ (1.4)

and using (1.3), one derives from the equations (1.1) the following $2\times
2$ matrix equation:

$\ddot{X}=AX,\quad X=\begin{pmatrix}X_{1}&Y_{1}\\\
X_{2}&Y_{2}\end{pmatrix},\quad
A=\begin{pmatrix}-m_{2}\rho_{3}-m_{13}\rho_{2}&m_{1}\rho_{32}\\\
m_{2}\rho_{31}&-m_{1}\rho_{3}-m_{32}\rho_{1}\end{pmatrix}\,,$ (1.5)

where

$\rho_{1}=\frac{1}{|z_{32}|^{3}},\quad\rho_{2}=\frac{1}{|z_{13}|^{3}},\quad\rho_{3}=\frac{1}{|z_{21}|^{3}},\quad
m_{ij}=m_{i}+m_{j},\quad\rho_{ij}=\rho_{i}-\rho_{j}\,.$ (1.6)

We define the determinants

$\Delta_{1}(t)=\left|\begin{array}[]{ll}X_{1}&Y_{1}\\\
X_{2}&Y_{2}\end{array}\right|(t),\quad\Delta_{2}(t)=\left|\begin{array}[]{ll}\dot{X}_{1}&\dot{Y}_{1}\\\
\dot{X}_{2}&\dot{Y}_{2}\end{array}\right|(t),\quad t\in[0,+\infty)\,.$ (1.7)

It is sufficient to show that $\Delta_{2}$ has an infinite number of zeros for
$t\geq 0$. As demonstrated by Montgomery [1], and later independently by the
author in [4, 5], the solution $t\mapsto(z_{1}(t),z_{2}(t),z_{3}(t))$, for
$t\in[0,+\infty)$, has an infinite number of consecutive syzygies. That is,
there exist $t_{i}\geq 0$, $i=1,2,\dots$, with $t_{i+1}>t_{i}$ such that
$\Delta_{1}(t_{i})=0$. As shown in [5, p. 6833], in the zero angular momentum
case:

$\dot{\Delta}_{1}(t)^{2}-4\Delta_{1}(t)\Delta_{2}(t)\geq 0,\quad\forall t\geq
0\,,$ (1.8)

where equality occurs only if the positions $z_{i}(t)$ and velocities
$\dot{z}_{i}(t)$ for $i=1,2,3$ are parallel.

Let $t_{i}<t_{i+1}<t_{i+2}$ be three consecutive syzygy moments. By Rolle’s
theorem, there exist $\tau\in(t_{i},t_{i+1})$ and $\eta\in(t_{i+1},t_{i+2})$
such that $\tau<\eta$ and $\dot{\Delta}_{1}(\tau)=\dot{\Delta}_{1}(\eta)=0$.
Without loss of generality, we can assume that $\Delta_{1}(\tau)>0$ and
$\Delta_{1}(\eta)<0$. Indeed, if for some $i\geq 0$,
$\Delta_{1}(t_{i})=\dot{\Delta}_{1}(t_{i})=0$, then according to (1.8), the
corresponding solution is a straight-line one (with all positions and
velocities lying on the same line) and will result in a triple collision since
the energy is negative.

Thus, by the Intermediate Value Theorem and (1.8), $\Delta_{2}$ has a zero in
the interval $[\tau,\eta]$. The proof is complete. ∎

##  2\. Bounds on Syzygy Moments in the Zero Angular Momentum Case

It is well-established that in the case of negative energy and zero angular
momentum, every collision-free solution in the interval $[0,+\infty)$ has an
infinite number of syzygies [1, 4, 5]. In the author’s previous works [4, 5],
an upper bound was determined for the moment when the very first syzygy
occurs, expressed as a function of the energy value and the masses alone. In
this section, we refine our results by providing both upper and lower bounds
that depend on the initial positions and velocities of the bodies, resulting
in significantly more precise estimates.

First, we present a preliminary result from the comparison theory for matrix
Riccati equations, developed by Eschenburg and Heintze in 1990 [6], which will
be utilised later in the proof of our main result.

Let $E$ be a finite-dimensional real vector space equipped with an inner
product $(\cdot,\cdot)$, and let $S(E)$ denote the space of self-adjoint
linear endomorphisms of $E$. Consider the Riccati differential equation with a
smooth coefficient curve $R:(0,t_{0})\to S(E)$:

$\dot{B}+B^{2}+R=0,$ (2.1)

where $B:(0,t_{0})\to S(E)$ is the solution.

Given two smooth coefficient curves $R_{1},R_{2}:(0,t_{0})\to S(E)$ with
$R_{1}\geq R_{2}$ (i.e., $R_{1}-R_{2}$ is positive semidefinite), we can
compare the solutions $B_{1}$ and $B_{2}$ of the Riccati equations:

$\dot{B}_{1}+B_{1}^{2}+R_{1}=0,$ (2.2)

and

$\dot{B}_{2}+B_{2}^{2}+R_{2}=0,$ (2.3)

subject to appropriate initial conditions.

######  Theorem 2.1 ([6]).

Let $R_{1},R_{2}:(0,t_{0})\to S(E)$ be smooth coefficient curves with
$R_{1}\geq R_{2}$. For $j=1,2$, let $B_{j}:(0,t_{j})\to S(E)$ be the solution
to the Riccati equation corresponding to $R_{j}$, with maximal
$t_{j}\in(0,+\infty]$. If the difference $U:=B_{2}-B_{1}$ has a continuous
extension to $0$ with $U(0)\geq 0$, then $t_{1}\leq t_{2}$ and $B_{1}\leq
B_{2}$ on $(0,t_{1})$.

We now assume that at the instant $t=0$, the configuration of the three-body
problem is different from a syzygy, i.e., the positions of the three bodies
are not parallel. In this case, the matrix $C_{0}=\dot{X}_{0}X_{0}^{-1}$,
where $X_{0}=X(0)$ and $\dot{X}_{0}=\dot{X}(0)$, is well-defined since
$\det(X_{0})\neq 0$. In the zero angular momentum case, as shown in [5], both
eigenvalues of $C_{0}$ are real. We denote the spectrum of $C_{0}$ as:

$\text{Spec}(C_{0})=\\{\lambda_{1},\lambda_{2}\\}\subset\mathbb{R}\,,$ (2.4)

and define

$\pi_{s}=\min\\{\lambda_{1},\lambda_{2}\\}\,.$ (2.5)

Finally, we suppose that during the motion, all mutual distances of the bodies
are bounded above and below by the constants $\alpha,\beta>0$, with
$\alpha<\beta$, i.e.,

$\alpha\leq|z_{ij}(t)|\leq\beta,\quad\forall\,i<j\,,$ (2.6)

and we define

$\theta_{\alpha}=\frac{\sqrt{M}}{\alpha^{3/2}},\quad\theta_{\beta}=\frac{\sqrt{M}}{\beta^{3/2}},\quad
M=m_{1}+m_{2}+m_{3}\,.$ (2.7)

######  Theorem 2.2.

Let $T_{s}>0$ be the moment in time when the first syzygy occurs. Then the
following bounds hold:

$\frac{1}{\theta_{\alpha}}\mathrm{arccot}\left(-\frac{\pi_{s}}{\theta_{\alpha}}\right)\leq
T_{s}\leq\frac{1}{\theta_{\beta}}\mathrm{arccot}\left(-\frac{\pi_{s}}{\theta_{\beta}}\right)\,.$
(2.8)

###### Proof.

Let $t\mapsto\phi_{i}(t)$, $i=1,2,3$, be arbitrary positive smooth functions
on $[0,+\infty)$. We consider a linear system of second order

$\ddot{Z}=\mathcal{A}Z,\quad Z=\begin{pmatrix}Z_{11}&Z_{12}\\\
Z_{21}&Z_{22}\end{pmatrix},\quad\mathcal{A}(\phi_{1},\phi_{2},\phi_{3})=\begin{pmatrix}-m_{2}\phi_{3}-m_{13}\phi_{2}&m_{1}\phi_{32}\\\
m_{2}\phi_{31}&-m_{1}\phi_{3}-m_{32}\phi_{1}\end{pmatrix}\,,$ (2.9)

where $\phi_{ij}=\phi_{i}-\phi_{j}$, $m_{ij}=m_{i}+m_{j}$.

If $(Z_{11},Z_{12})=(X_{1},Y_{1})$, $(Z_{21},Z_{22})=(X_{2},Y_{2})$, and
$\phi_{i}=\rho_{i}$ for $i=1,2,3$, then the equations (2.9) and (1.5)
coincide.

It is a straightforward computation to verify that for any solution $Z$ of the
system (2.9), the following identity holds:

$\frac{1}{m_{1}}\begin{vmatrix}Z_{11}&Z_{12}\\\
\dot{Z}_{11}&\dot{Z}_{12}\end{vmatrix}+\frac{1}{m_{2}}\begin{vmatrix}Z_{21}&Z_{22}\\\
\dot{Z}_{21}&\dot{Z}_{22}\end{vmatrix}+\frac{1}{m_{3}}\begin{vmatrix}Z_{11}+Z_{21}&Z_{12}+Z_{22}\\\
\dot{Z}_{11}+\dot{Z}_{21}&\dot{Z}_{12}+\dot{Z}_{22}\end{vmatrix}=k,\quad
k\in\mathbb{R}\,,$ (2.10)

which is an analog of the angular momentum conservation law in the three-body
problem written in the form (1.5).

The matrix $B=\dot{Z}Z^{-1}$ is a solution of the following Riccati equation:

$\dot{B}+B^{2}+R=0,\quad R=-\mathcal{A}\,,$ (2.11)

an equation utilized in the study of the three-body problem in [4].

We introduce the matrices

$I_{2}=\left(\begin{array}[]{rr}1&0\\\
0&1\end{array}\right),\quad\tilde{A}_{1}=\left(\begin{array}[]{rr}\frac{m_{32}}{2}&0\\\
-m_{2}&-\frac{m_{32}}{2}\end{array}\right),\quad\tilde{A}_{2}=\left(\begin{array}[]{rr}-\frac{m_{13}}{2}&-m_{1}\\\
0&\frac{m_{13}}{2}\end{array}\right)\,.$ (2.12)

Then, as shown in [5, p. 6833], $B$ can be presented, if $k=0$, in the
following form:

$B=\frac{\dot{\delta}}{2\delta}I_{2}+\frac{b}{m_{2}\delta}\tilde{A}_{1}-\frac{a}{m_{1}\delta}\tilde{A}_{2}\,,$
(2.13)

where

$\beta=\frac{1}{2}\left(\frac{m_{3}}{m_{1}}+1\right),\quad\gamma=\frac{1}{2}\left(\frac{m_{3}}{m_{2}}+1\right)\,,$
(2.14)

and

$\delta=\begin{vmatrix}Z_{11}&Z_{12}\\\ {Z}_{21}&{Z}_{22}\end{vmatrix},\quad
a=\begin{vmatrix}Z_{11}&Z_{12}\\\ \dot{Z}_{11}&\dot{Z}_{12}\end{vmatrix},\quad
b=\begin{vmatrix}Z_{21}&Z_{22}\\\ \dot{Z}_{21}&\dot{Z}_{22}\end{vmatrix}\,.$
(2.15)

In order to apply Theorem 2.1 to equation (2.11), the solution $B$ and the
matrix $R$ should be symmetric. This can be achieved by the linear
transformation $\tilde{B}=P^{-1}BP$ with the invertible matrix $P$ defined by:

$P=\left(\begin{array}[]{cc}-\frac{m_{1}}{m_{13}}&\frac{1}{m_{13}}\sqrt{\frac{m_{1}m_{3}M}{m_{2}}}\\\
1&0\end{array}\right),\quad M=m_{1}+m_{2}+m_{3}\,.$ (2.16)

Indeed, it is easy to check that the matrices $P^{-1}\tilde{A}_{i}P$, $i=1,2$,
and $\tilde{R}=P^{-1}RP$ are symmetric. Since $\tilde{B}$ is a linear
combination of $I_{2}$ and $P^{-1}\tilde{A}_{i}P$, it is also symmetric. Thus,
equation (2.11) becomes:

$\dot{\tilde{B}}+\tilde{B}^{2}+\tilde{R}=0\,,$ (2.17)

and Theorem 2.1 can be applied.

######  Lemma 2.1.

The matrix
$\tilde{\mathcal{A}}(\phi_{1},\phi_{2},\phi_{3})=P^{-1}\mathcal{A}(\phi_{1},\phi_{2},\phi_{3})P$
is negative semidefinite.

###### Proof.

Since $\tilde{\mathcal{A}}$ is symmetric, it is sufficient to show that the
eigenvalues of $\mathcal{A}$ are negative. We have:

$\det(\mathcal{A})=M(m_{3}\phi_{1}\phi_{2}+m_{2}\phi_{1}\phi_{3}+m_{1}\phi_{2}\phi_{3})\geq
0,\quad\text{Tr}(\mathcal{A})=-\left(m_{32}\phi_{1}+m_{13}\phi_{2}+m_{21}\phi_{3}\right)\leq
0\,,$ (2.18)

since $\phi_{i}\geq 0$ for $i=1,2,3$.

Thus, $\tilde{\mathcal{A}}$ is negative semidefinite. ∎

Let us consider two Riccati equations:

$\dot{\tilde{B}}_{1}+\tilde{B}_{1}^{2}+\tilde{R}_{1}=0,\quad\tilde{R}_{1}=-P^{-1}\mathcal{A}(\rho_{1},\rho_{2},\rho_{3})P\,,$
(2.19)

and

$\dot{\tilde{B}}_{2}+\tilde{B}_{2}^{2}+\tilde{R}_{2}=0,\quad\tilde{R}_{2}=-P^{-1}\mathcal{A}\left(\frac{1}{\beta^{3}},\frac{1}{\beta^{3}},\frac{1}{\beta^{3}}\right)P\,$
(2.20)

subject to the same initial conditions

$\tilde{B}_{1}(0)=\tilde{B}_{2}(0)=P^{-1}\dot{X}_{0}X_{0}^{-1}P\,.$ (2.21)

By the linearity of $\mathcal{A}$ as a function of $\phi_{i}$ ($i=1,2,3$), we
have

$\tilde{R}_{2}-\tilde{R}_{1}=P^{-1}\mathcal{A}\left(\rho_{1}-\frac{1}{\beta^{3}},\rho_{2}-\frac{1}{\beta^{3}},\rho_{3}-\frac{1}{\beta^{3}}\right)P,$

and $\rho_{i}-\frac{1}{\beta^{3}}\geq 0$ for $i=1,2,3$ by the definition of
$\beta$.

According to Lemma 2.1, it follows that $\tilde{R}_{1}\geq\tilde{R}_{2}$.

Equation (2.20) can be easily solved by setting $\tilde{R}_{2}=\dot{Y}Y^{-1}$.
This leads to the equivalent equation:

$\ddot{Y}=-\theta_{\beta}^{2}Y,\quad\theta_{\beta}=\frac{\sqrt{M}}{\beta^{3/2}}\,.$
(2.22)

The solution to the Cauchy problem defined by (2.22) and the initial
conditions $Y(0)=Y_{0}=P^{-1}X_{0}$,
$\dot{Y}(0)=\dot{{Y}}_{0}=P^{-1}\dot{X}_{0}$ is:

$Y(t)=\cos(\theta_{\beta}t)Y_{0}+\frac{1}{\theta_{\beta}}\sin(\theta_{\beta}t)\dot{{Y}}_{0}\,,\quad
t\geq 0\,.$ (2.23)

Therefore, the corresponding maximal solution of (2.20) is:

$\tilde{B}_{2}(t)=\dot{Y}(t)Y(t)^{-1},\quad t\in[0,t_{2})\,,$ (2.24)

where $t_{2}$ is the first positive zero of $t\mapsto\det(Y(t))$.

One has

$\det(Y(t_{2}))=0\,\Longleftrightarrow\,\det(\dot{X}_{0}X_{0}^{-1}+\theta_{\beta}\cot(\theta_{\beta}t_{2})I_{2})=0\,,$
(2.25)

since $X_{0}$ and $P$ are invertible matrices. As a consequence, we obtain:

$-\theta_{\beta}\cot(\theta_{\beta}t_{2})\in\text{Spec}(C_{0}),\quad
C_{0}=\dot{X}_{0}X_{0}^{-1}\,,$ (2.26)

and

$t_{2}=\frac{1}{\theta_{\beta}}\mathrm{arccot}\left(-\frac{\pi_{s}}{\theta_{\beta}}\right)\,.$
(2.27)

Let $T_{s}=t_{1}>0$ be the first syzygy moment for the solution of the three-
body problem $t\mapsto z_{i}(t)$, $i=1,2,3$. Then, using a similar argument,
we show that the maximal solution of the Cauchy problem defined by (2.19) and
(2.21) is defined in the interval $[0,t_{1})$. Therefore, according to Theorem
2.1, $T_{s}\leq t_{2}$, and the upper bound in (2.8) is proven.

To prove the lower bound, consider the Riccati equations:

$\dot{\tilde{B}}_{1}+\tilde{B}_{1}^{2}+\tilde{R}_{1}=0,\quad\tilde{R}_{1}=-P^{-1}\mathcal{A}\left(\frac{1}{\alpha^{3}},\frac{1}{\alpha^{3}},\frac{1}{\alpha^{3}}\right)P\,,$
(2.28)

and

$\dot{\tilde{B}}_{2}+\tilde{B}_{2}^{2}+\tilde{R}_{2}=0,\quad\tilde{R}_{2}=-P^{-1}\mathcal{A}(\rho_{1},\rho_{2},\rho_{3})P\,,$
(2.29)

where $\tilde{R}_{1}\geq\tilde{R}_{2}$. Applying the same arguments as before,
the proof of Theorem 2.2 is complete. ∎

######  Remark 2.1.

The uniform bounds $|z_{ij}|\geq\alpha$, $\forall\,i<j$ have a very natural
astrophysical interpretation: such motion corresponds to the collision-free
movement of three rigid planets, each having the same radius $R=\alpha/2$. As
seen from the proof of Theorem 2.2, the uniform bounds (2.6) and the absence
of collisions are required only for the period of time preceding the first
syzygy.

##  3\. Conclusion: Numerical Validation and Open Questions

To numerically verify our findings, we will consider the figure-eight periodic
solution with equal masses $m_{1}=m_{2}=m_{3}=1$ and zero angular momentum, as
described in [3]. The initial positions and velocities of the bodies are set
as follows:

$z_{1}(0)=1.08075-i0.0126893,\quad z_{2}(0)=-0.570154+i0.350807,\quad
z_{3}(0)=-z_{1}(0)-z_{2}(0),$ (3.1)

and

$\dot{z}_{1}(0)=0.0193421+i0.467219,\quad\dot{z}_{2}(0)=1.0852-i0.174718,\quad\dot{z}_{3}(0)=-\dot{z}_{1}(0)-\dot{z}_{2}(0)\,.$
(3.2)

Numerically, we can find that the first syzygy occurs at $T_{s}=0.55431$ (see
Figure 1). The constants $\alpha$ and $\beta$ are estimated to be
$\alpha=0.690526$ and $\beta=2$. Using formulas (2.7), we determine
$\theta_{\alpha}=3.01849$ and $\theta_{\beta}=0.612372$.

Using the initial conditions (3.1) and (3.2), we compute:

$X_{0}=\begin{pmatrix}0.734528&1.35841\\\
0.755791&-0.470708\end{pmatrix},\quad\text{Spec}(X_{0})=\\{1.31082,-1.047\\},\quad\pi_{s}=-1.047,$
(3.3)

and the bounds (2.8) give us:

$0.409781\leq T_{s}\leq 0.864231,$ (3.4)

which are quite satisfactory.

We would like to highlight several open questions. It would be interesting to
find bounds analogous to (2.8) for the velocity syzygy moments using similar
ideas based on the comparison of matrix Riccati equations. It is important to
note that the estimates provided by our Theorem 2.2 are valid only if the
initial configuration at $t=0$ is not a syzygy. However, with some effort and
by employing Theorem 2.1, one could establish bounds similar to (2.8) even in
the case where the initial configuration is a syzygy.

Finally, generalising to the case of non-zero angular momentum is an
intriguing and important challenge.

Akcnowledgments

I would like to express my special gratitude to Richard Montgomery for useful
suggestions and valuable remarks and to Jean-Claude Sikorav for attracting my
attention to work [6].

Figure 1. The figure shows the trajectory of the figure-eight solution, as
described in [3]. The initial positions of the bodies at $t=0$ are represented
by filled points. Additionally, the configuration of the first syzygy, which
occurs at $t=0.55431$, is indicated by unfilled points.

## References

  * [1] R. Montgomery, ”The zero angular momentum, three-body problem: All but one solution has syzygies,” Ergod. Th. Dynam. Sys., vol. 27, pp. 311-340, 2007.
  * [2] A. Wintner, The analytical foundations of Celestial Mechanics, Princeton University Press, 1941.
  * [3] A. Chenciner and R. Montgomery, ”A remarkable periodic solution of the three-body problem in the case of equal masses,” Annals of Mathematics, vol. 152, pp. 881-901, 2000.
  * [4] A. Tsygvintsev, ”On the existence of generalised syzygies in the planar three-body problem,” C.R. Acad. Sci. Paris, vol. 361, pp. 331-335, 2023.
  * [5] A. Tsygvintsev, ”On some collinear configurations in the planar three-body problem,” Nonlinearity, vol. 36, pp. 6827-6838, 2023.
  * [6] J.H. Eschenburg and E. Heintze, ”Comparison theory for Riccati equations,” Manuscripta Math, vol. 68, pp. 209-214, 1990.


††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work††thanks: These authors contributed
equally to this work

# Sudden death of entanglement with Hamiltonian ensemble assisted by auxiliary
qubits

Congwei Lu  Department of Physics, Applied Optics Beijing Area Major
Laboratory, Beijing Normal University, Beijing 100875, China  Wanting He
Department of Physics, Applied Optics Beijing Area Major Laboratory, Beijing
Normal University, Beijing 100875, China  Jun Wang  Department of Physics,
Applied Optics Beijing Area Major Laboratory, Beijing Normal University,
Beijing 100875, China  Haibo Wang  Department of Physics, Applied Optics
Beijing Area Major Laboratory, Beijing Normal University, Beijing 100875,
China  Qing Ai<EMAIL_ADDRESS>Department of Physics, Applied Optics
Beijing Area Major Laboratory, Beijing Normal University, Beijing 100875,
China

###### Abstract

In this paper, we theoretically propose a method to simulate the longitudinal
relaxation of a single qubit by coupling it to an auxiliary qubit. In order to
mimic the finite-temperature relaxation, we utilize the Hamiltonian-ensemble
approach [Kropf, Gneiting, and Buchleitner, Phys. Rev. X 6, 031023 (2016)] and
in each realization the auxiliary qubit possesses a random level spacing. The
longitudinal relaxation arises as a consequence of the ensemble average and
the interaction between the working qubit and the auxiliary qubit.
Furthermore, we apply this approach to investigate the influence of the
longitudinal relaxation and the transverse relaxation on the entanglement
dynamics of two qubits. It is discovered that the sudden death of the
entanglement will occur as long as the longitudinal relaxation is present. The
transverse relaxation assists the longitudinal relaxation and thus accelerates
the finite-time disentanglement.

##  I Introduction

In an open quantum system, the interaction between the system and the
environment will lead to the exchange of information and energy between them.
As a result, the dynamic behavior of the system is quite different from that
of an isolated system Leggett  _et al._ (1987); Breuer and Petruccione (2002);
Zhang  _et al._ (2012); Breuer  _et al._ (2016); de Vega and Alonso (2017);
Chen  _et al._ (2022). Generally, there are two types of relaxations, i.e.,
transverse relaxation and longitudinal relaxation. The latter will result in
population transfer and decay of the off-diagonal terms of the density matrix,
while the former will only decrease the coherence. These two behaviors play a
crucial role in quantum information processing.

Recently, it was proposed that the quantum dynamics of an open quantum system
can be simulated by the ensemble-averaged state of many random isolated
systems Buluta and Nori (2009); Georgescu  _et al._ (2014); Kropf  _et al._
(2016); Gneiting and Nori (2017); Chen  _et al._ (2018); Wang  _et al._
(2018); Zhang  _et al._ (2021). The process of ensemble averaging over each
random realization will result in averaging all random phases, thus inducing
the loss of phase information, i.e., dephasing. However, due to the classical
property of noise, most of the previous quantum simulation approaches can only
simulate the longitudinal relaxation at the high-temperature limit Soare  _et
al._ (2014a); Wang  _et al._ (2018); Zhang  _et al._ (2021); Soare  _et al._
(2014b); Zhen  _et al._ (2016). The quantum simulation of the longitudinal
relaxation at finite temperature is rarely studied. Since the dissipative
environment can cause finite-time disentanglement Yu and Eberly (2004);
Almeida  _et al._ (2007); Yu and Eberly (2009), enhanced relaxation at avoided
level-crossing Yang  _et al._ (2020); Onizhuk  _et al._ (2021); Head-Marsden
_et al._ (2021), and optical non-reciprocity by detailed balance Yao and Ai
(2022), it may be interesting to in depth understand the influence of
transverse relaxation and longitudinal relaxation on the sudden death of
entanglement and its dynamic simulation.

In this paper, in order to simulate the longitudinal relaxation at finite
temperature, we introduce an auxiliary qubit which interacts with the working
qubit. In order to mimic the longitudinal relaxation, we effectively prepare a
large number of systems including the auxiliary qubit and the working qubit.
They evolve from the same initial state but the auxiliary qubit possesses a
random level spacing. By averaging over the different realizations, we can
effectively simulate the longitudinal relaxation. Our analytical results
demonstrate that this approach can well simulate the finite-temperature
longitudinal relaxation with the dissipation rate linearly dependent on time.
Here, the dissipation rate scales linearly with the variance of the level
spacing of the auxiliary qubit and the interaction strength between the
auxiliary qubit and the working qubit. The initial state of the auxiliary
qubit and the interaction strength jointly determine the distribution of the
steady state. We further investigate the effects of longitudinal and
transverse relaxation on the entanglement of two-qubit. We let the first
working qubit interact with the auxiliary qubit to mimic the longitudinal
relaxation, and apply a random field on the second working qubit to simulate
the transverse relaxation. These two working qubits are initialized in the
maximum-entangled state and the concurrence is utilized to characterize the
dynamics of entanglement. Our simulations show that due to the longitudinal
relaxation, the sudden death of entanglement happens at a finite time. In this
case, the transverse relaxation of the second working qubit will accelerate
the finite-time disentanglement of the two qubits. However, if the
longitudinal relaxation is absent, the entanglement will go to zero when the
time approaches infinity.

The rest of the paper is structured as follows. In Sec. II, we introduce the
quantum-simulation approach by Hamiltonian ensemble. In Sec.III, we simulate
the longitudinal relaxation of a single qubit at a finite temperature. The
effects of the noise fluctuation and interaction strength on entanglement
sudden death are investigated in Sec. IV. Finally, we conclude our main
discoveries in Sec. V.

Figure 1: Schematic diagram for simulating the longitudinal relaxation of a
single qubit by the Hamiltonian-ensemble approach. Each realization is
composed of a working qubit A and an auxiliary qubit B. They evolve
independently from the same initial state $\rho(0)$. The ensemble-averaged
state $\overline{\operatorname{Tr}_{B}\rho_{\varepsilon}(t)}$ is then obtained
by averaging over all reduced density matrix
$\operatorname{Tr}_{B}\rho_{\varepsilon_{i}}(t)$ ($i=1,2,\cdots,N$) of the
working qubit.

##  II Hamiltonian-Ensemble Approach

First of all, we shall give a brief introduction to the quantum-simulation
approach by an ensemble of Hamiltonians Kropf  _et al._ (2016); Chen  _et al._
(2018); Wang  _et al._ (2018); Zhang  _et al._ (2021), as schematically
illustrated in Fig. 1. A general open quantum system can be characterized by a
total Hamiltonian $\hat{H}_{T}=\hat{H}_{S}+\hat{H}_{E}+\hat{H}_{I}$ Breuer and
Petruccione (2002), where $\hat{H}_{S}$ is the system Hamiltonian,
$\hat{H}_{E}$ is the environment Hamiltonian, and $\hat{H}_{I}$ represents
their interaction. The time evolution of the open system can be described as
$\rho_{T}(t)=\hat{U}\rho_{T}(0)\hat{U}^{\dagger}$, with
$\hat{U}=\exp(-i\hat{H}_{T}t/\hbar)$. Thus, the density matrix of the system
can be obtained by partially tracing over the environmental degrees of
freedom, i.e., $\rho_{S}(t)=\operatorname{Tr}_{E}[\rho_{T}(t)]$. To simulate
the open quantum dynamics, we utilize the Hamiltonian ensemble

$\displaystyle\\{(\hat{H}_{\varepsilon},p_{\varepsilon})\\},$ (1)

where the subscript $\varepsilon$ denotes each realization in the ensemble.
The single realization Hamiltonian $\hat{H}_{\varepsilon}$ occurring with
probability $p_{\varepsilon}$ reads

$\displaystyle\hat{H}_{\varepsilon}=\hat{H}_{S}+\hat{H}_{E}^{\varepsilon}+\hat{V}_{\varepsilon}.$
(2)

$\hat{H}_{E}^{\varepsilon}$ and $\hat{V}_{\varepsilon}$ are utilized to
simulate the environment and its interaction with system. We suppose that each
realization begins from the same initial state
$\rho_{\varepsilon}(0)=\rho(0)$. The corresponding evolution at time $t$ is
given by
$\rho_{\varepsilon}(t)=\hat{U}_{\varepsilon}\rho(0)\hat{U}_{\varepsilon}^{\dagger}$,
with $\hat{U}_{\varepsilon}=\exp(-i\hat{H}_{\varepsilon}t/\hbar)$. Finally, we
trace over the environmental degree of freedom in each realization and then
average over all realizations, i.e.,

$\displaystyle\langle\rho(t)\rangle$
$\displaystyle=\overline{\operatorname{Tr}_{E}\rho_{\varepsilon}(t)}=\int
d\varepsilon p_{\varepsilon}\operatorname{Tr}_{E}\rho_{\varepsilon}(t).$ (3)

Hereafter, all ensemble-averaged quantities will be marked with a bar. In the
next section, as an example, we utilize the ensemble-averaged quantum dynamics
of the state $\langle\rho(t)\rangle$ to simulate the longitudinal relaxation
behavior of a single qubit.

##  III Longitudinal Relaxation of A Single Qubit

In the previous investigations, due to the classical property of the noise,
the quantum simulation approach can only simulate the longitudinal relaxation
at the high-temperature limit Wang  _et al._ (2018); Zhang  _et al._ (2021);
Chen  _et al._ (2022). In this paper, we introduce an auxiliary qubit, which
interacts with the working qubit as schematically illustrated in Fig. 1, in
order to simulate the longitudinal relaxation at finite temperatures. In the
Hamiltonian ensemble, the Hamiltonian of a realization reads

$\displaystyle\hat{H}_{\varepsilon}=\frac{\omega_{0}}{2}(\omega_{A}\sigma_{z}^{A}+\varepsilon\sigma_{z}^{B})+f(\omega_{0}\varepsilon)\sigma_{+}^{A}\sigma_{-}^{B}+\textrm{h.c.},$
(4)

where we set $\hbar=1$, $\sigma_{z}$ is the Pauli matrix, $\sigma_{\pm}$ are
the raising and lowering operators and $\omega_{0}$ is the unit for frequency.
Hereafter, we set $\omega_{0}=1$ in the following simulations.
$f(\varepsilon)=f^{\ast}(\varepsilon)$ is the coupling strength between the
working qubit and the auxiliary qubit. For simplicity, the two-qubit product
states are relabeled as

$\begin{split}|1\rangle_{AB}=|+-\rangle_{AB},~{}|2\rangle_{AB}=|-+\rangle_{AB},\\\
|3\rangle_{AB}=|++\rangle_{AB},~{}|4\rangle_{AB}=|--\rangle_{AB},\end{split}$
(5)

where $|\pm\pm\rangle_{AB}\equiv|\pm\rangle_{A}\otimes|\pm\rangle_{B}$ denote
the eigenstates of Pauli operator $\sigma_{z}^{A}\otimes\sigma_{z}^{B}$. Since
$\sigma_{z}^{A}\otimes\sigma_{z}^{B}$ is the conserved quantity, i.e.,
$[\sigma_{z}^{A}\otimes\sigma_{z}^{B},~{}\hat{H}_{\varepsilon}]=0$, we can
rewrite $\hat{H}_{\varepsilon}$ as a block-diagonal matrix in the basis listed
in Eq. (5),

$\displaystyle\hat{H}_{\varepsilon}=\left(\begin{array}[]{cc}\hat{H}_{\varepsilon}^{-}&0\\\
0&\hat{H}_{\varepsilon}^{+}\end{array}\right),$ (8)

where

$\begin{split}\hat{H}_{\varepsilon}^{-}&=\frac{1}{2}(\omega_{A}-\varepsilon)\sigma_{z}+f(\varepsilon)\sigma_{x},\\\
\hat{H}_{\varepsilon}^{+}&=\frac{1}{2}(\omega_{A}+\varepsilon)\sigma_{z}.\end{split}$
(9)

The total evolution operator $\exp(-i\hat{H}_{\varepsilon}t)$ can be
represented as

$\displaystyle\hat{U}_{\varepsilon}=\left(\begin{array}[]{cc}\hat{U}_{\varepsilon}^{-}&0\\\
0&\hat{U}_{\varepsilon}^{+}\end{array}\right),$ (12)

where

$\begin{split}\hat{U}_{\varepsilon}^{-}&=\cos(Et)-i\frac{\sin(Et)}{E}\left[\frac{\omega_{A}-\varepsilon}{2}\sigma_{z}+f(\varepsilon)\sigma_{x}\right],\\\
\hat{U}_{\varepsilon}^{+}&=\cos\frac{(\omega_{A}+\varepsilon)t}{2}-i\sin\frac{(\omega_{A}+\varepsilon)t}{2}\sigma_{z},\end{split}$
(13)

with

$\displaystyle
E=\sqrt{\frac{1}{4}(\omega_{A}-\varepsilon)^{2}+f(\varepsilon)^{2}}.$ (14)

Figure 2: The longitudinal relaxation against $x_{B}$,
$\overline{\varepsilon^{2}}$ and $\alpha$ for $\omega_{A}=0$. The population
of the subsystem A at $|+\rangle$ under ensemble average (a) for
$x_{B}=1,~{}0.8,~{}0.6$, when $\overline{\varepsilon^{2}}=1$ and $\alpha=5$,
(b) for $\overline{\varepsilon^{2}}=10,1,0.5$ when $x_{B}=0.8$ and $\alpha=5$,
(c) for $\alpha=3,1,0.6$, when $x_{B}=0.8$ and $\overline{\varepsilon^{2}}=1$.
The analytic results are shown as solid lines, while the dotted lines are
generated by $N=8000$ random samples with $\overline{\varepsilon}=0$.

We assume that the initial state is a product state
$\ket{\psi(0)}=\ket{\phi(0)}_{A}\otimes\ket{\varphi(0)}_{B}$, where the
working qubit $A$ is in the state $\ket{\phi(0)}_{A}=\ket{-}_{A}$, and the
auxiliary qubit $B$ is in a superposition state
$\ket{\varphi(0)}_{B}=x_{B}\ket{+}_{B}+y_{B}\ket{-}_{B}$. The reduced density
matrix of the working qubit $A$ can be obtained by partially tracing over $B$,

$\displaystyle\rho_{A}(t)=\operatorname{Tr}_{B}(\hat{U}_{\varepsilon}\ket{\psi(0)}\bra{\psi(0)}\hat{U}_{\varepsilon}^{\dagger}).$
(15)

For simplicity, we assume that the coupling strength is proportional to the
detuning between the two-qubit, that is,

$\displaystyle
f(\varepsilon)=\sqrt{\alpha^{2}-\frac{1}{4}}(\varepsilon-\omega_{A}),$ (16)

where $\alpha\geq 1/2$. Specifically, the matrix elements of $\rho_{A}$ read

$\begin{split}\rho_{A}^{++}(t)&\equiv\bra{+}\rho_{A}(t)\ket{+}\\\
&=\frac{c^{2}}{2}|x_{B}|^{2}\left\\{1-\cos[2\alpha(\varepsilon-\omega_{A})t]\right\\},\\\
\rho_{A}^{+-}(t)&\equiv\bra{+}\rho_{A}(t)\ket{-}\\\
&=-icx_{B}y_{B}^{*}e^{-\frac{i}{2}(\varepsilon+\omega_{A})t}\sin\left[\alpha(\varepsilon-\omega_{A})t\right],\end{split}$
(17)

where $c=\sqrt{4\alpha^{2}-1}/(2\alpha)$. Here, we assume that $\varepsilon$
is subject to a Gaussian distribution and the ensemble-averaged state
$\langle\rho(t)\rangle$ defined in Eq. (3) can be written as

$\displaystyle\langle\rho(t)\rangle_{++}\equiv$
$\displaystyle\overline{\rho_{A}^{++}(t)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}c^{2}|x_{B}|^{2}\left[1-\cos(2\alpha\omega_{A}t)e^{-2\alpha^{2}\overline{\varepsilon^{2}}t^{2}}\right],$
$\displaystyle\langle\rho(t)\rangle_{+-}\equiv$
$\displaystyle\overline{\rho_{A}^{+-}(t)}$ (18) $\displaystyle=$
$\displaystyle\frac{1}{2}cx_{B}y_{B}^{*}\left[e^{-\frac{1}{2}(\alpha+\frac{1}{2})^{2}\overline{\varepsilon^{2}}t^{2}}e^{i(\alpha-\frac{1}{2})\omega_{A}t}\right.$
$\displaystyle\left.-e^{-\frac{1}{2}(\alpha-\frac{1}{2})^{2}\overline{\varepsilon^{2}}t^{2}}e^{-i(\alpha+\frac{1}{2})\omega_{A}t}\right].$

Here, we have utilized the $2n$-order moment identity of a Gaussian
distribution with mean zero, i.e.,
$\overline{\varepsilon^{2n}}=\overline{\varepsilon^{2}}(2n)!/(2^{n}n!)$
Goodman (2015).

This result can be considered as the thermalization of the working qubit A in
a thermal bath, which is in a thermal equilibrium at temperature $T$. When the
system A reaches the thermal equilibrium, its probability at the excited state
is $P_{+}=e^{-\beta\Delta}(1+e^{-\beta\Delta})^{-1}$ Breuer and Petruccione
(2002), with $\beta=1/k_{B}T$, and $k_{B}$ being the Boltzmann constant, where
$\Delta$ is the energy-level difference between the two levels. The
dissipation rate $\Gamma(t)=2\alpha^{2}\overline{\varepsilon^{2}}t$ is linear
with respect to time $t$, which can be utilized to improve the quantum
metrology and thus achieve Zeno limit Chin  _et al._ (2012); Matsuzaki  _et
al._ (2011); Long  _et al._ (2022). To simulate this steady-state distribution
at arbitrary temperature $T$, we can effectively tune the coupling strength,
i.e., $\alpha$, and the initial state of the auxiliary qubit B, i.e., $x_{B}$
and $y_{B}$, to fulfill that
$P_{+}=\overline{\rho_{A}^{++}(t\rightarrow\infty)}$. Notice that for any
temperature $T$, this formula can always be satisfied because $0\leq c<1$ and
$0\leq x_{B}\leq 1$, and thus
$0\leq\overline{\rho_{A}^{++}(t\rightarrow\infty)}<1/2$. Here, as a
demonstration, we simulate a process of a single qubit, initialized in the
ground state, relaxation to the equilibrium state at a finite temperature $T$
through interaction with the heat reservoir. In this simulation, this process
can be controlled by the initial state of the auxiliary qubit, the properties
of the noise characterized by $\overline{\varepsilon^{2}}$ and the coupling
strength characterized by $\alpha$. The behaviors of the longitudinal
relaxation against the parameters, i.e., $x_{B}$, $\overline{\varepsilon^{2}}$
and $\alpha$, are plotted in Fig. 2.

In Fig. 2(a), we leave $\overline{\varepsilon^{2}}$ and $\alpha$ unchanged and
only vary $x_{B}$. We find that $x_{B}$ does not change the relaxation time
but the steady-state population. In contrast, we can observe in Fig. 2(b), the
relaxation time will decrease with the increase of the noise variance
$\overline{\varepsilon^{2}}$, which does not change the steady-state
population. Interestingly, the coupling strength between the subsystem A and
the auxiliary qubit B determines both the steady-state population and the
relaxation time. As shown in Fig. 2(c), the greater the coupling strength is,
the shorter the relaxation time becomes while the larger the steady-state
population reaches. It is worth noting that in a real Markovian open quantum
system, the steady-state population is only subject to the temperature of the
environment and the energy-level difference of the system Breuer and
Petruccione (2002), and thus is independent of the coupling strength. However,
in our model, the steady-state population of the subsystem A is also
determined by the coupling strength. Similar discovery has also been observed
in Ref. Dong  _et al._ (2007).

In order to verify our numerical simulation, the analytical and numerical
results of full elements of the density matrix are compared in Fig. 3. We find
that when $\omega_{A}\neq 0$, both the longitudinal and transverse relaxation
behaviors demonstrate an oscillatory decay, but the decay of the transverse
relaxation is slower. As shown in Fig. 3(b), the interaction between the
subsystem A and the auxiliary qubit B will induce the coherence between the
ground state and the excited state, since the auxiliary qubit is initially in
a superposition. However, when the steady state is reached, the coherence of
the subsystem A disappears and thus becomes a mixed state. Moreover, since the
numerical results agree with the analytical results, our numerical simulations
are reliable. To summarize, in this section, we utilize an auxiliary qubit to
effectively simulate the longitudinal relaxation process of a single qubit at
arbitrary temperature. It is found that the initial state $x_{B}$, frequency
variance $\overline{\varepsilon^{2}}$ of the auxiliary qubit, and the
interaction strength between the working qubit and the auxiliary qubit
together determine the relaxation time and steady-state population of the
longitudinal relaxation.

Figure 3:  Comparison of analytical and numerical results for the longitudinal
relaxation of a single qubit. In the numerical calculation, we select $N=5000$
random samples $\left\\{\varepsilon\right\\}$ of Gaussian distribution with
variance 0.6 and expectation 0. The quantum dynamics of (a) the excited-state
population $\langle\rho(t)\rangle_{++}$, (b) the real part of the coherence
$\operatorname{Re}[\langle\rho(t)\rangle_{+-}]$, when $\omega_{A}=4$,
$x_{B}=0.9$, $\overline{\varepsilon^{2}}=0.6$, and $\alpha=1$.

##  IV Finite-Time Disentanglement

In this section, we simulate the quantum dynamics of two-qubit disentanglement
and investigate the effects of longitudinal and transverse relaxation on the
disentanglement behavior. The system includes two working qubits, i.e., qubit
$A_{1}$ and $B_{1}$, where the former interacts with an auxiliary qubit, i.e.,
qubit $A_{2}$, as schematically demonstrated in Fig. 4. Thus, the total
Hamiltonian can be written in two parts as
$\hat{H}_{\varepsilon}=\hat{H}_{\varepsilon}^{A}+\hat{H}_{\varepsilon}^{B}$,
where

$\displaystyle\hat{H}_{\varepsilon}^{A}=$
$\displaystyle\frac{\omega_{0}}{2}(\omega_{A}\sigma_{z}^{A_{1}}+\varepsilon_{A}\sigma_{z}^{A_{2}})+f(\omega_{0}\varepsilon_{A})\sigma_{+}^{A_{1}}\sigma_{-}^{A_{2}}+\textrm{h.c.},$
$\displaystyle\hat{H}_{\varepsilon}^{B}=$
$\displaystyle\frac{\omega_{0}}{2}(\omega_{B}\sigma_{z}^{B_{1}}+\varepsilon_{B}\sigma_{z}^{B_{1}}).$
(19)

The composite system composed of $A_{1}$ and $B_{1}$ is initialized in the
maximum-entangled state, i.e.,
$\ket{\psi(0)}_{A_{1}B_{1}}=(\ket{++}_{A_{1}B_{1}}+\ket{--}_{A_{1}B_{1}})/\sqrt{2}$.
We let $A_{1}$ interact with an auxiliary qubit $A_{2}$ to mimic the
longitudinal relaxation, as depicted in Sec. III, where $f(\varepsilon_{A})$
is the coupling strength between $A_{1}$ and $A_{2}$. For qubit $B_{1}$, we
apply a random energy-level spacing described by $\varepsilon_{B}$ to simulate
the transverse relaxation Zhang  _et al._ (2021); Wang  _et al._ (2018); Kropf
_et al._ (2016).

Figure 4: Schematic illustration of simulating sudden death of entanglement in
a two-qubit system. Qubit $A_{1}$ and $B_{1}$ are initialized in the maximum-
entangled state and have no interaction with each other. The random energy
level-spacing characterized by $\varepsilon_{B}$ is used to simulate the
transverse relaxation of $B_{1}$. We simulate the longitudinal noise of
$A_{1}$ through the interaction between $A_{1}$ and the auxiliary qubit
$A_{2}$. Figure 5: The critical disentanglement time $t_{c}$ against $\alpha$
and $\overline{\varepsilon_{A}^{2}}$ for $x=0.2$,
$\overline{\varepsilon_{B}^{2}}=0$ and $\omega_{A}=0$.

Before the ensemble average, we first of all solve the quantum dynamics of
each realization. The initial state of the three qubits reads

$\displaystyle\rho(0)=x\ket{\psi_{1}(0)}\bra{\psi_{1}(0)}+y\ket{\psi_{2}(0)}\bra{\psi_{2}(0)},$
(20)

where

$\begin{split}\ket{\psi_{1}(0)}&=\frac{1}{\sqrt{2}}\ket{+}_{A_{2}}\otimes(\ket{++}_{A_{1}B_{1}}+\ket{--}_{A_{1}B_{1}}),\\\
\ket{\psi_{2}(0)}&=\frac{1}{\sqrt{2}}\ket{-}_{A_{2}}\otimes(\ket{++}_{A_{1}B_{1}}+\ket{--}_{A_{1}B_{1}}).\end{split}$
(21)

In other words, the two working qubits $A_{1}$ and $B_{1}$ are in a maximum-
entangled state, while the auxiliary qubit $A_{2}$ is in a mixed state. At
time $t$, $\ket{\psi_{1}(0)}$ evolves into the state

$\displaystyle\ket{\psi_{1}(t)}$ $\displaystyle=$
$\displaystyle\eta_{1}\ket{+++}_{A_{2}A_{1}B_{1}}+\eta_{2}\ket{-+-}_{A_{2}A_{1}B_{1}}$
(22) $\displaystyle+\eta_{3}\ket{+--}_{A_{2}A_{1}B_{1}},$

where
$\eta_{1}=\exp\left[-i(\omega_{A}+\varepsilon_{A}+\omega_{B}+\varepsilon_{B})t/2\right]/\sqrt{2}$,
because $\ket{+++}_{A_{2}A_{1}B_{1}}$ is the eigenstate of
$\hat{H}_{\varepsilon}$. In the invariant subspace spanned by the basis
$\\{\ket{-+-}_{A_{2}A_{1}B_{1}},\ket{+--}_{A_{2}A_{1}B_{1}}\\}$, the effective
Hamiltonian can be simplified as

$\displaystyle\hat{H}_{\varepsilon}=-\frac{1}{2}(\omega_{B}+\varepsilon_{B})I+\frac{1}{2}(\omega_{A}-\varepsilon_{A})\sigma_{z}+f(\varepsilon_{A})\sigma_{x}.$
(23)

And thus the evolution operator
$\hat{U}_{\varepsilon}=\exp(-i\hat{H}_{\varepsilon}t)$ reads

$\displaystyle\hat{U}_{\varepsilon}$
$\displaystyle=\left\\{\cos(E_{A}t)+\frac{i\sin(E_{A}t)}{E_{A}}\left[\frac{\varepsilon_{A}-\omega_{A}}{2}\sigma_{z}-f(\varepsilon_{A})\sigma_{x}\right]\right\\}$
$\displaystyle\times e^{i(\varepsilon_{B}+\omega_{B})t/2},$ (24)

where

$\displaystyle
E_{A}=\sqrt{\frac{1}{4}(\omega_{A}-\varepsilon_{A})^{2}+f(\varepsilon_{A})^{2}}.$
(25)

Since
$\left[\eta_{2},\eta_{3}\right]^{T}=\hat{U}_{\varepsilon}\left[0,1/\sqrt{2}\right]^{T}$,
the three coefficients of $\ket{\psi_{1}(t)}$ are explicitly given as

$\begin{split}\eta_{1}\\!\\!&=\\!\\!\frac{1}{\sqrt{2}}e^{-\frac{i}{2}(\omega_{A}+\varepsilon_{A}+\omega_{B}+\varepsilon_{B})t},\\\
\eta_{2}\\!\\!&=\\!\\!\frac{-i}{\sqrt{2}}e^{\frac{i}{2}(\varepsilon_{B}+\omega_{B})t}f(\varepsilon_{A})\frac{\sin
E_{A}t}{E_{A}},\\\
\eta_{3}\\!\\!&=\\!\\!\frac{e^{\frac{i}{2}(\varepsilon_{B}+\omega_{B})t}}{\sqrt{2}}\left[\cos
E_{A}t+\frac{i\sin
E_{A}t}{2E_{A}}(\omega_{A}-\varepsilon_{A})\right].\end{split}$ (26)

Suppose $\ket{\psi_{2}(t)}$ can be expanded as

$\displaystyle\ket{\psi_{2}(t)}$ $\displaystyle=$
$\displaystyle\xi_{1}\ket{+++}_{A_{2}A_{1}B_{1}}+\xi_{2}\ket{-+-}_{A_{2}A_{1}B_{1}}$
(27) $\displaystyle+\xi_{3}\ket{+--}_{A_{2}A_{1}B_{1}}.$

Following the above steps, we can obtain

$\begin{split}\xi_{1}\\!\\!&=\\!\\!\frac{1}{\sqrt{2}}e^{\frac{i}{2}(\omega_{A}+\varepsilon_{A}+\omega_{B}+\varepsilon_{B})t},\\\
\xi_{2}\\!\\!&=\\!\\!\frac{-i}{\sqrt{2}}e^{-\frac{i}{2}(\varepsilon_{B}+\omega_{B})t}f(\varepsilon_{A})\frac{\sin
E_{A}t}{E_{A}},\\\
\xi_{3}\\!\\!&=\\!\\!\frac{e^{-\frac{i}{2}(\varepsilon_{B}+\omega_{B})t}}{\sqrt{2}}\left[\cos
E_{A}t-\frac{i\sin
E_{A}t}{2E_{A}}(\omega_{A}-\varepsilon_{A})\right].\end{split}$ (28)

Because we have solved the quantum dynamics of the three qubits, by tracing
over the auxiliary qubit $A_{2}$, we can obtain the reduced density matrix of
the two working qubits $\rho_{A_{1}B_{1}}(t)=\operatorname{Tr}_{A_{2}}\rho(t)$
in the basis
$\\{\ket{+-}_{A_{1}B_{1}},\ket{++}_{A_{1}B_{1}},\ket{--}_{A_{1}B_{1}},\ket{-+}_{A_{1}B_{1}}\\}$
as

$\displaystyle\rho_{A_{1}B_{1}}(t)=\left(\begin{array}[]{cccc}a(t)&0&0&0\\\
0&b(t)&z(t)&0\\\ 0&z^{*}(t)&c(t)&0\\\ 0&0&0&d(t)\end{array}\right).$ (33)

As in the Sec. III, we define the coupling strength as
$f(\varepsilon_{A})=(\varepsilon_{A}-\omega_{A})\sqrt{\alpha^{2}-1/4}$, where
$\alpha\geq 1/2$. The non-vanishing matrix elements are given by

$\displaystyle a(t)$ $\displaystyle=$
$\displaystyle\frac{xc^{2}}{2}\gamma(t)^{2},$ $\displaystyle b(t)$
$\displaystyle=$
$\displaystyle\frac{x}{2}+\frac{y}{2}\left[1-c^{2}\gamma(t)^{2}\right],$
$\displaystyle c(t)$ $\displaystyle=$
$\displaystyle\frac{y}{2}+\frac{x}{2}\left[1-c^{2}\gamma(t)^{2}\right],$ (34)
$\displaystyle z(t)$ $\displaystyle=$
$\displaystyle\frac{x+y}{2}\zeta(t)\left[1-2\gamma(t)^{2}-\frac{i\gamma(t)}{2\alpha}\right],$
$\displaystyle d(t)$ $\displaystyle=$
$\displaystyle\frac{yc^{2}}{2}\gamma(t)^{2},$

where
$\zeta(t)=\exp\left[\frac{-i}{2}(\omega_{A}+\varepsilon_{A}+2\omega_{B}+2\varepsilon_{B})t\right]$,
$\gamma(t)=\sin(\alpha(\omega_{A}-\varepsilon_{A})t)$ and
$c=\sqrt{4\alpha^{2}-1}/(2\alpha)$. In our numerical simulation, we assume
that $\varepsilon_{A}$ and $\varepsilon_{B}$ are subject to independent
Gaussian distributions. After the ensemble average, the non-vanishing matrix
elements of Eq. (33) are explicitly given as

$\begin{split}\overline{a(t)}\\!\\!&=\\!\\!\frac{xc^{2}}{4}\left[1-\cos(2\alpha\omega_{A}t)e^{-2\alpha^{2}\overline{\varepsilon_{A}^{2}}t^{2}}\right],\\\
\overline{b(t)}\\!\\!&=\\!\\!\frac{x}{2}+\frac{y}{2}\left[1-\frac{c^{2}}{2}(1-\cos(2\alpha\omega_{A}t)e^{-2\alpha^{2}\overline{\varepsilon_{A}^{2}}t^{2}})\right],\\\
\overline{c(t)}\\!\\!&=\\!\\!\frac{y}{2}+\frac{x}{2}\left[1-\frac{c^{2}}{2}(1-\cos(2\alpha\omega_{A}t)e^{-2\alpha^{2}\overline{\varepsilon_{A}^{2}}t^{2}})\right],\\\
\overline{z(t)}\\!\\!&=\\!\\!\frac{1}{4}e^{-\frac{1}{2}\overline{\varepsilon_{B}^{2}}t^{2}}e^{-\frac{i}{2}(\omega_{A}+2\omega_{B})t}[e^{i\alpha\omega_{A}t}e^{-\frac{1}{2}(\alpha+\frac{1}{2})^{2}\overline{\varepsilon_{A}^{2}}t^{2}}\\\
\\!\\!&\times\\!\\!(1-\frac{1}{2\alpha})+e^{-i\alpha\omega_{A}t}e^{-\frac{1}{2}(\alpha-\frac{1}{2})^{2}\overline{\varepsilon_{A}^{2}}t^{2}}(1+\frac{1}{2\alpha})],\\\
\overline{d(t)}\\!\\!&=\\!\\!\frac{yc^{2}}{4}\left[1-\cos(2\alpha\omega_{A}t)e^{-2\alpha^{2}\overline{\varepsilon_{A}^{2}}t^{2}}\right].\end{split}$
(35)

Figure 6:  The concurrence $C(t)$ as a function of time $t$ in the presence of
both longitudinal and transverse relaxation for
$\overline{\varepsilon_{B}^{2}}=0,~{}0.5,~{}2$ when $x=0.2$,
$\overline{\varepsilon_{A}^{2}}=0.5$, $\alpha=1$, and $\omega_{A}=0$. Notice
that the transverse relaxation is turned off when
$\overline{\varepsilon_{B}^{2}}=0$. The analytic results are shown as solid
lines. The dotted lines are generated by $N=300$ random samples with
$\overline{\varepsilon_{A}}=\overline{\varepsilon_{B}}=0$.

To investigate the disentanglement behavior of the composite system composed
of $A_{1}$ and $B_{1}$ under both longitudinal and transverse relaxation, we
utilize the concurrence Wootters (1998) to characterize the entanglement
property
$C(\overline{\rho_{A_{1}B_{1}}})=\text{max}(0,\sqrt{\kappa_{1}}-\sqrt{\kappa_{2}}-\sqrt{\kappa_{3}}-\sqrt{\kappa_{4}})$,
where $\kappa_{i}$’s are the eigenvalues of the matrix $\mathcal{G}$ in
decreasing order

$\displaystyle\mathcal{G}\equiv\overline{\rho_{A_{1}B_{1}}}\left(\sigma_{y}^{A}\otimes\sigma_{y}^{B}\right)\overline{\rho_{A_{1}B_{1}}}^{*}\left(\sigma_{y}^{A}\otimes\sigma_{y}^{B}\right),$
(36)

where $\sigma_{y}^{\alpha}$ ($\alpha=A,B$) are the Pauli operators. When
$C=1$, the two working qubits are in the maximum-entangled state, while $C=0$,
they are disentangled with each other. Here, we can simplify the concurrence
as

$\displaystyle
C(\overline{\rho_{A_{1}B_{1}}})=2\text{max}\left\\{0,|\overline{z}|-\sqrt{\overline{a}\overline{d}}\right\\}.$
(37)

At the beginning, $A_{1}$ and $B_{1}$ are initialized in the maximum-entangled
state with $|\overline{z}|=1/2$ and $\sqrt{\overline{a}\overline{d}}=0$. In
the following, we will show two categories of disentanglement in our
simulation. In the first category, the entanglement tends to vanish only when
the time approaches infinite. In the second category, the entanglement decays
exactly to zero at a critical disentanglement time $t_{c}$ and remains zero
thereafter.

Figure 7: The concurrence $C(t)$ as a function of time $t$ in the presence of
both longitudinal and transverse relaxation for $\omega_{A}=0,~{}3,~{}6$, when
$x=0.2$, $\overline{\varepsilon_{A}^{2}}=0.5$,
$\overline{\varepsilon_{B}^{2}}=0.5$, and $\alpha=1$. The analytic results are
shown as solid lines. The dotted lines are generated by $N=300$ random samples
with $\overline{\varepsilon_{A}}=\overline{\varepsilon_{B}}=0$.

We first investigate the influence of the longitudinal relaxation on the
entanglement properties when there is no transverse relaxation, i.e.,
$\overline{\varepsilon_{B}^{2}}=0$, in the system. In this case, the
disentanglement behavior of the system is dominated by the coupling strength
$\alpha$ and the noise fluctuation $\overline{\varepsilon_{A}^{2}}$.
Obviously, when $\alpha=1/2$, i.e., the subsystem $A_{1}$ does not interact
with the auxiliary qubit $A_{2}$, the system will always be in the maximum-
entangled state. When $\alpha>1/2$, the non-vanishing noise fluctuation
$\overline{\varepsilon_{A}^{2}}$ will determine whether the entanglement of
the system can disappear at a finite time. When $\alpha>1/2$ and
$\overline{\varepsilon_{A}^{2}}=0$, i.e., the longitudinal relaxation is
turned off, $|\overline{z}|$ and $\sqrt{\overline{a}\overline{d}}$ can be
written as

$\begin{split}|\overline{z}|&=\frac{\sqrt{2}}{4}\sqrt{1+\frac{1}{4\alpha^{2}}+\left(1-\frac{1}{4\alpha^{2}}\right)\cos(2\alpha\omega_{A}t)},\\\
\sqrt{\overline{a}\overline{d}}&=\sqrt{xy}\frac{4\alpha^{2}-1}{16\alpha^{2}}\left[1-\cos(2\alpha\omega_{A}t)\right].\end{split}$
(38)

There always exist $t=n\pi/(\alpha\omega_{A})$ with $n\in\mathcal{Z}$, such
that $|\overline{z}|=1/2$ and $\sqrt{\overline{a}\overline{d}}=0$, and thus
the entanglement of the system will not disappear persistently. However, if we
turn on the longitudinal relaxation, i.e., $\alpha>1/2$ and
$\overline{\varepsilon_{A}^{2}}>0$, $|\overline{z}|$ tends to zero and
$\sqrt{\overline{a}\overline{d}}=\sqrt{xy}(4\alpha^{2}-1)/(16\alpha^{2})$ when
the time approaches infinity. Thus, as long as $xy>0$, there exists a finite
$t_{c}$, making the entanglement disappear after time $t_{c}$. The critical
disentanglement time $t_{c}$ against $\overline{\varepsilon_{A}^{2}}$ and
$\alpha$ without transverse relaxation is shown in Fig. 5. We find that when
$\alpha>1/2$ and $\overline{\varepsilon_{A}^{2}}>0$, the entanglement will
disappear at a finite time and $t_{c}$ decays monotonically and rapidly as
$\alpha$ and $\overline{\varepsilon_{A}^{2}}$ increase.

Now we consider the case when there is only transverse relaxation, i.e.,
$\overline{\varepsilon_{B}^{2}}>0$ and $\overline{\varepsilon_{A}^{2}}=0$. The
evolution of $|\overline{z}|$ and $\sqrt{\overline{a}\overline{d}}$ can be
written as

$\displaystyle|\overline{z}|$
$\displaystyle=\frac{\sqrt{2}}{4}e^{-\frac{1}{2}\overline{\varepsilon_{B}^{2}}t^{2}}\sqrt{1+\frac{1}{4\alpha^{2}}+\left(1-\frac{1}{4\alpha^{2}}\right)\cos(2\alpha\omega_{A}t)},$
$\displaystyle\sqrt{\overline{a}\overline{d}}$
$\displaystyle=\sqrt{xy}\frac{4\alpha^{2}-1}{16\alpha^{2}}\left[1-\cos(2\alpha\omega_{A}t)\right].$
(39)

Obviously, $|\overline{z}|$ tends to zero when the time approaches infinite.
No matter how long the time passes, there always exist
$2\alpha\omega_{A}t=2k\pi$ with $k\in\mathcal{Z}$ so that
$\sqrt{\overline{a}\overline{d}}$ vanishes. Thus, the concurrence will not
constantly remain zero after a finite time but will go to zero at infinite
time.

We can conclude that when the longitudinal relaxation exists, i.e.,
$\alpha>1/2$ and $\overline{\varepsilon_{A}^{2}}>0$, entanglement will
disappear at a finite time. The larger the transverse noise fluctuates, the
faster the entanglement disappears. And the oscillation frequency of $A_{1}$
hardly affects the time of disentanglement. The dynamics of concurrence
against $\overline{\varepsilon_{B}^{2}}$ and $\omega_{A}$ in this case are
shown in Figs. 6 and 7, respectively. From Fig. 6, we find that even if
$\overline{\varepsilon_{B}^{2}}=0$, the concurrence will remain zero after the
critical disentanglement time and decay faster with the increase of
$\overline{\varepsilon_{B}^{2}}$. From Fig. 7, we find that when the noise
fluctuation and interaction strength are kept constant, the higher the
frequency of subsystem $A_{1}$ is, the faster the entanglement $C$ oscillates,
but the critical disentanglement time is almost the same. When the system is
immune to the longitudinal relaxation, i.e.,
$\overline{\varepsilon_{A}^{2}}=0$ or $\alpha=1/2$, the entanglement will not
die at a finite time but will disappear at infinite time due to the transverse
relaxation.

##  V conclusion

We have utilized the Hamiltonian-ensemble approach assisted by an auxiliary
qubit to simulate the longitudinal relaxation of a single qubit in open
quantum system. Concretely, many sets of auxiliary qubits with random level
spacings are used to simulate the environmental noise. They are prepared in
the same initial state and interact with the working qubit. The theoretical
results show that the simulated dynamics of the working qubit can be described
by a real thermalization process. We can simulate the equilibrium-state
distribution at any temperature, which is determined by the initial state of
the auxiliary qubit, noise fluctuation and interaction strength. Furthermore,
we simulate the dynamics of two-qubit entanglement initialized in the maximum-
entangled state. We let the first qubit interact with the auxiliary qubit and
relax longitudinally, and let the second qubit relax transversely. We find
that if there is not longitudinal relaxation on the first qubit but transverse
relaxation on the second qubit, the entanglement of the system will disappear
when the time approaches infinity. However, if the longitudinal relaxation
exists, no matter whether the transverse relaxation is present or not, the
entanglement of the two-qubit will disappear after a finite time. And the
larger the noise fluctuation is, the faster the entanglement will decay.

###### Acknowledgements.

Q.A. thanks the financial support from Beijing Natural Science Foundation
under Grant No. 1202017 and the National Natural Science Foundation of China
under Grant Nos. 11674033, 11505007, and Beijing Normal University under Grant
No. 2022129. H.B.W. thanks the financial support from National Natural Science
Foundation of China under Grant No. 61675028 and National Natural Science
Foundation of China under Grant No. 12274037.

## References

  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Rev. Mod. Phys. 59, 1 (1987).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Presss, 2002).
  * Zhang  _et al._ (2012) W. M. Zhang, P. Y. Lo, H. N. Xiong, M. W. Y. Tu, and F. Nori, “General non-Markovian dynamics of open quantum systems,” Phys. Rev. Lett. 109, 170402 (2012).
  * Breuer  _et al._ (2016) H. P. Breuer, E. M. Laine, J. Piilo, and B. Vacchini, “Colloquium: Non-Markovian dynamics in open quantum systems,” Rev. Mod. Phys. 88, 021002 (2016).
  * de Vega and Alonso (2017) I. de Vega and D. Alonso, “Dynamics of non-Markovian open quantum systems,” Rev. Mod. Phys. 89, 015001 (2017).
  * Chen  _et al._ (2022) X. Y. Chen, N. N. Zhang, W. T. He, X. Y. Kong, F. G. Deng, Q. Ai, and G. L. Long, “Global correlation and local information flows in controllable non-Markovian open quantum dynamics,” npj Quantum Inf. 8, 22 (2022).
  * Buluta and Nori (2009) I. Buluta and F. Nori, “Quantum simulators,” Science 326, 108 (2009).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, “Quantum simulation,” Rev. Mod. Phys. 86, 153 (2014).
  * Kropf  _et al._ (2016) C. M. Kropf, C. Gneiting, and A. Buchleitner, “Effective dynamics of disordered quantum systems,” Phys. Rev. X 6, 031023 (2016).
  * Gneiting and Nori (2017) C. Gneiting and F. Nori, “Disorder-induced dephasing in backscattering-free quantum transport,” Phys. Rev. Lett. 119, 176802 (2017).
  * Chen  _et al._ (2018) H. B. Chen, C. Gneiting, P. Y. Lo, Y. N. Chen, and F. Nori, “Simulating open quantum systems with Hamiltonian ensembles and the nonclassicality of the dynamics,” Phys. Rev. Lett. 120, 030403 (2018).
  * Wang  _et al._ (2018) B. X. Wang, M. J. Tao, Q. Ai, T. Xin, N. Lambert, D. Ruan, Y. C. Cheng, F. Nori, F. G. Deng, and G. L. Long, “Efficient quantum simulation of photosynthetic light harvesting,” npj Quantum Inf. 4, 52 (2018).
  * Zhang  _et al._ (2021) N. N. Zhang, M. J. Tao, W. T. He, X. Y. Chen, X. Y. Kong, F. G. Deng, N. Lambert, and Q. Ai, “Efficient quantum simulation of open quantum dynamics at various Hamiltonians and spectral densities,” Front. Phys 16, 51502 (2021).
  * Soare  _et al._ (2014a) A. Soare, H. Ball, D. Hayes, J. Sastrawan, M. C. Jarratt, J. J. McLoughlin, X. Zhen, T. J. Green, and M. J. Biercuk, “Experimental noise filtering by quantum control,” Nat. Phys. 10, 825 (2014a).
  * Soare  _et al._ (2014b) A. Soare, H. Ball, D. Hayes, X. Zhen, M. C. Jarratt, J. Sastrawan, H. Uys, and M. J. Biercuk, “Experimental bath engineering for quantitative studies of quantum control,” Phys. Rev. A 89, 042329 (2014b).
  * Zhen  _et al._ (2016) X. L. Zhen, F. H. Zhang, G. R. Feng, H. Li, and G. L. Long, “Optimal experimental dynamical decoupling of both longitudinal and transverse relaxations,” Phys. Rev. A 93, 022304 (2016).
  * Yu and Eberly (2004) T. Yu and J. H. Eberly, “Finite-time disentanglement via spontaneous emission,” Phys. Rev. Lett. 93, 140404 (2004).
  * Almeida  _et al._ (2007) M. P. Almeida, F. de Melo, M. Hor-Meyll, A. Salles, S. P. Walborn, P. H. Souto Ribeiro, and L. Davidovich, “Environment-induced sudden death of entanglement,” Science 316, 579 (2007).
  * Yu and Eberly (2009) T. Yu and J. H. Eberly, “Sudden death of entanglement,” Science 323, 598 (2009).
  * Yang  _et al._ (2020) Z.-S. Yang, Y.-X. Wang, M.-J. Tao, W. Yang, M. Zhang, Q. Ai, and F.-G. Deng, “Longitudinal relaxation of a nitrogen-vacancy center in a spin bath by generalized cluster-correlation expansion method,” Ann. Phys. (N.Y.) 413, 168063 (2020).
  * Onizhuk  _et al._ (2021) M. Onizhuk, K. C. Miao, J. P. Blanton, H. Ma, C. P. Anderson, A. Bourassa, D. D. Awschalom, and G. Galli, “Probing the coherence of solid-state qubits at avoided crossings,” PRX Quantum 2, 010311 (2021).
  * Head-Marsden  _et al._ (2021) K. Head-Marsden, J. Flick, C. J. Ciccarino, and P. Narang, “Quantum information and algorithms for correlated quantum matter,” Chem. Rev. 121, 3061 (2021).
  * Yao and Ai (2022) Y.-X. Yao and Q. Ai, “Optical non-reciprocity in coupled resonators inspired by photosynthetic energy transfer,” arXiv:2208.05841 (2022).
  * Goodman (2015) J. W. Goodman,  _Statistical Optics_ (Wiley, Hoboken, NJ, 2015).
  * Chin  _et al._ (2012) A. W. Chin, S. F. Huelga, and M. B. Plenio, “Quantum metrology in non-Markovian environments,” Phys. Rev. Lett. 109, 233601 (2012).
  * Matsuzaki  _et al._ (2011) Y. Matsuzaki, S. C. Benjamin, and J. Fitzsimons, “Magnetic field sensing beyond the standard quantum limit under the effect of decoherence,” Phys. Rev. A 84, 012103 (2011).
  * Long  _et al._ (2022) X. Y. Long, W. T. He, N. N. Zhang, K. Tang, Z. D Lin, H. F. Liu, X. F. Nie, G. R. Feng, J. Li, T. Xin, Q. Ai, and D. W. Lu, “Entanglement-enhanced quantum metrology in colored noise by quantum Zeno effect,” Phys. Rev. Lett. 129, 070502 (2022).
  * Dong  _et al._ (2007) H. Dong, S. Yang, X. F. Liu, and C. P. Sun, “Quantum thermalization with couplings,” Phys. Rev. A 76, 044104 (2007).
  * Wootters (1998) W. K. Wootters, “Entanglement of formation of an arbitrary state of two qubits,” Phys. Rev. Lett. 80, 2245 (1998).