# Linear multistep methods with repeated global Richardson extrapolation

I. Fekete, L. Lóczi<EMAIL_ADDRESS>corresponding author, Department
of Applied Analysis and Computational Mathematics, ELTE Eötvös Loránd
University, Pázmány P. s. 1/c, H-1117 Budapest<EMAIL_ADDRESS>Department of Numerical Analysis, ELTE Eötvös Loránd University, Pázmány P. s.
1/c, H-1117 Budapest, Hungary, and Department of Differential Equations, BME
Budapest University of Technology and Economics

###### Abstract

In this work, we further investigate the application of the well-known
Richardson extrapolation (RE) technique to accelerate the convergence of
sequences resulting from linear multistep methods (LMMs) for numerically
solving initial-value problems of systems of ordinary differential equations.
By extending the ideas of our previous paper, we now utilize some advanced
versions of RE in the form of repeated RE (RRE). Assume that the underlying
LMM—the base method—has order $p$ and RE is applied $\ell$ times. Then we
prove that the accelerated sequence has convergence order $p+\ell$. The
version we present here is global RE (GRE, also known as passive RE), since
the terms of the linear combinations are calculated independently. Thus, the
resulting higher-order LMM-RGRE methods can be implemented in a parallel
fashion and existing LMM codes can directly be used without any modification.
We also investigate how the linear stability properties of the base method
(e.g. $A$\- or $A(\alpha)$-stability) are preserved by the LMM-RGRE methods.

Keywords: linear multistep methods; Richardson extrapolation; Adams methods;
BDF methods; convergence; region of absolute stability
Mathematics Subject Classification (2020): 65L05, 65L06

##  1 Introduction

Let us consider the initial-value problem

$y^{\prime}(t)=f(t,y(t)),\quad y(t_{0})=y_{0},$ (1)

where $m\in\mathbb{N}^{+}$ and
$f:\mathbb{R}\times\mathbb{R}^{m}\to\mathbb{R}^{m}$ is a given sufficiently
smooth function. Suppose we approximate the unique solution $y$ of (1) on an
interval $[t_{0},t_{\text{final}}]$ by applying a $k$-step linear multistep
method (LMM)

$\sum_{j=0}^{k}\alpha_{j}y_{n+j}=\sum_{j=0}^{k}h\beta_{j}f_{n+j}$ (2)

of order $p\geq 1$ on a uniform grid $\\{t_{n}\\}$, where, as usual,
$h:=t_{n+1}-t_{n}>0$ is the step size (or grid length),
$f_{m}:=f(t_{m},y_{m})$, and the numbers $\alpha_{j}\in\mathbb{R}$ and
$\beta_{j}\in\mathbb{R}$ ($j=0,\ldots,k$) are the given method coefficients
with $\alpha_{k}\neq 0$. The LMM (2) generates a sequence $y_{n}(h)$ which is
supposed to approximate the exact solution at $t_{n}$, that is,
$y_{n}(h)\approx y(t_{n})$. In this work, this LMM will also be referred to as
the base method or the underlying method.

Classical Richardson extrapolation (RE) [11, 12] is a technique to accelerate
the convergence of numerical sequences depending on a small parameter, by
eliminating the lowest order error term(s) from the corresponding asymptotic
expansion. When solving (1) numerically, the parameter in RE can be chosen as
the discretization step size $h$. The application of RE to sequences generated
by one-step (e.g., Runge–Kutta) methods is described, for example, in [3, 7].
In [13], global (also known as passive) and local (or active) versions of RE
are implemented with Runge–Kutta sequences. These combined methods find
application in many areas (see, e.g., [4, 15]).

When carrying out global RE (GRE), one considers a suitable linear combination
of two approximations, one generated on a coarser grid and one on a finer
grid, to obtain a better approximation of the solution $y$ of (1). This
extrapolation is called global, because the sequences on the two grids are
computed independently and their linear combination is formed only in the last
step. Taking this idea further, one can consider several approximations on
finer and finer grids to improve convergence even further. Such a procedure is
called repeated global Richardson extrapolation (RGRE). To describe RGRE, let
us consider some nested grids with grid lengths, say, $h$, $h/2$, $h/4$, and
$h/8$, and the corresponding sequences $y_{n}(h)$,
$y_{2n}\left(\frac{h}{2}\right)$, $y_{4n}\left(\frac{h}{4}\right)$, and
$y_{8n}\left(\frac{h}{8}\right)$. Suppose that RE is applied
$\ell\in\mathbb{N}^{+}$ times, and let us denote the corresponding linear
combination of the approximations by $r_{n}^{[\ell]}(h)$ (thus, the $\ell=1$
case corresponds to classical RE). In [14], the authors show the formulae
$r_{n}^{[\ell]}(h)$ for several values of $\ell$, here we just reproduce the
first 3 cases:

$r_{n}^{[1]}(h):=\frac{2^{p}\cdot
y_{2n}\left(\frac{h}{2}\right)-y_{n}(h)}{2^{p}-1}$ (3)

$r_{n}^{[2]}(h):=\frac{2^{2p+1}\cdot y_{4n}\left(\frac{h}{4}\right)-3\cdot
2^{p}\cdot
y_{2n}\left(\frac{h}{2}\right)+y_{n}(h)}{\left(2^{p}-1\right)\left(2^{p+1}-1\right)}$
(4)

$r_{n}^{[3]}(h):=\frac{2^{3p+3}\cdot y_{8n}\left(\frac{h}{8}\right)-7\cdot
2^{2p+1}\cdot y_{4n}\left(\frac{h}{4}\right)+7\cdot 2^{p}\cdot
y_{2n}\left(\frac{h}{2}\right)-y_{n}(h)}{\left(2^{p}-1\right)\left(2^{p+1}-1\right)\left(2^{p+2}-1\right)}.$
(5)

In [14], all the sequences $y_{n}$ are assumed to be generated by Runge–Kutta
methods applied to (1).

######  Remark 1.1.

Note the slight difference in the terminology: in [14], $q$ denotes the number
of RE repetitions, while here $\ell$ denotes the number of RE applications. In
other words, $q=\ell-1$, and classical RE corresponds to $\ell=1$ or $q=0$.

######  Remark 1.2.

In [1], the authors describe multiple RE, which is another advanced version of
Richardson extrapolation besides repeated RE. Based on (3), multiple RE is
defined as

$r_{n}^{\mathrm{multiple}}(h):=\frac{2^{p+1}\cdot
r_{2n}^{[1]}\left(\frac{h}{2}\right)-r_{n}^{[1]}(h)}{2^{p+1}-1}.$

This formula is identical to $r_{n}^{[2]}(h)$ in (4). (The context in [1] is,
however, different, since they deal with local RE applied to Runge–Kutta
methods as underlying methods.)

In our recent work [5], we have shown that the sequence $r_{n}^{[1]}(h)$ also
converges to the solution of (1) if its component sequences $y_{2n}$ and
$y_{n}$ are generated by a LMM (2).

In the present work, we extend these ideas and prove that if the base method
(2) (satisfying some assumptions $None$–$None$ listed in Section 2) is of
order $p$, and global Richardson extrapolation is applied $\ell$ times (for
any $\ell\geq 1$), then the sequences $r_{n}^{[\ell]}(h)$ converge to the
solution of (1) with order of convergence $p+\ell$. We can refer to such
procedure as LMM-RGRE, or, when the number of RE applications is made
explicit, LMM-$\ell$GRE.

The structure of our paper is as follows. Section 1.2 summarizes some
notation. In Section 2, the improved convergence of LMM-RGREs is proved if the
underlying LMM is an Adams–Bashforth, Adams–Moulton or BDF method (see
Corollary 2.7). Linear stability of LMM-RGREs is investigated in Section 3.
Finally, several numerical tests are presented in Section 4 that demonstrate
the expected convergence order.

###  1.1 Conclusions

  * •

To implement a LMM-RGRE method, existing LMM codes can directly be used
without any modification.

  * •

The proof of Theorem 2.2 does not depend on the number of RE applications
(although, for simplicity, we have listed only the cases $\ell=1,2,3$), one
only needs to assume more smoothness on the function $f$ for larger values of
$\ell$.

  * •

Although our recent work [5] can be considered as the $\ell=1$ special case of
Theorem 2.2 of this paper, still, the assumptions of Theorem 2.2 on the
closeness of the initial values are relaxed (see Remark 2.5 below): instead of
${\cal{O}}(h^{p+1})$-closeness, it is sufficient to assume
${\cal{O}}(h^{p})$-closeness for the starting values of the LMM.

  * •

The computational cost of a LMM-RGRE increases for larger values of the RE
repetitions. In general, due to the computations on finer grids, applying GRE
$\ell$ times requires approximately $2^{\ell+1}-1$ as much computation as the
underlying LMM. However, this is compensated by the higher convergence order
$p+\ell$.

  * •

Regarding linear stability, when the underlying LMM is a BDF method, for
example, LMM-RGREs preserve the $A(\alpha)$-stability angles. In particular,
when $\ell=2$ and the base method is a BDF2 method, we obtain a
$4^{\text{th}}$-order $A$-stable method. Another example is the BDF5 method
with $\ell=2$, resulting in a $7^{\text{th}}$-order method with
$A(\alpha)$-stability angle $\approx 51.839^{\circ}$ (see also [5, Table 1]
for the $\ell=1$ case).

  * •

The proof of Theorem 2.2 is based on [7, Section III.9]. Since they establish
the existence of asymptotic expansions of the global error for (strictly
stable) general linear methods (GLMs), one can clearly apply RGRE in that
context as well, and accelerate the convergence of numerical sequences
generated by such GLMs.

  * •

The application of local RE to LMMs is the subject of an ongoing study.

###  1.2 Notation

We assume throughout this work that $0\in\mathbb{N}$. The Kronecker product of
two matrices is denoted by $\otimes$ (for its definition and properties, see,
e.g., [2]). For a (complex) square matrix $A$, let $\mathfrak{M}_{A}$ denote
its minimal polynomial (that is, the unique univariate polynomial with leading
coefficient equal to 1 and having the least degree such that
$\mathfrak{M}_{A}(A)$ is the zero matrix). For a set $S\subset\mathbb{C}$ and
for $a>0$, we define $a\,S:=\\{az:z\in S\\}$. For $2\leq k\in\mathbb{N}$,
AB$k$, AM$k$ and BDF$k$ denote, respectively, the $k$-step Adams–Bashforth,
Adams–Moulton, and BDF methods.

##  2 Convergence analysis

Given the $\alpha$-coefficients of the LMM (2), we define—according to [7,
Section III.4]—the matrix

$A:=\left(\begin{array}[]{ccccc}-\alpha_{k-1}/\alpha_{k}&-\alpha_{k-2}/\alpha_{k}&\ldots&-\alpha_{1}/\alpha_{k}&-\alpha_{0}/\alpha_{k}\\\
1&0&\ldots&0&0\\\ 0&1&\ldots&0&0\\\ 0&0&\ddots&\vdots&\vdots\\\
0&0&\ldots&1&0\\\ \end{array}\right)$ (6)

so that $A\in\mathbb{R}^{k\times k}$ and its lower left block is the
$(k-1)\times(k-1)$ identity matrix.

We assume that

  * $None$

the right-hand side $f$ of the initial-value-problem (1) is sufficiently
smooth;

  * $None$

the LMM (2) has order $p\geq 1$;

  * $None$

for any $h\in[0,h_{0}]$, the $k$ starting values of the sequence

$\\{y_{0}(h),y_{1}(h),\ldots,y_{k-1}(h)\\}$

to initiate the LMM are each ${\cal{O}}(h^{p})$-close to the corresponding
exact solution values

$\\{y(0),y(h),\ldots,y((k-1)h)\\}$

of the initial-value-problem (1);

  * $None$

the eigenvalues of $A$ lie in the closed unit disk of the complex plane, the
only eigenvalue with modulus 1 is 1, and its algebraic multiplicity is also 1.

######  Remark 2.1.

It is known that, due to consistency of the LMM (2), 1 is always an eigenvalue
of $A$.

The following theorem establishes the convergence of sequences (3)–(5). The
grids have grid length $h$, $h/2$, $h/4$, and $h/8$, respectively, so the
fixed grid point $t^{*}:=t_{0}+nh$ is part of all of these grids. (As always
in this context, any implied constant in an ${\cal{O}}$ symbol is independent
of $n$ and $h$.)

######  Theorem 2.2.

Under the above assumptions $None$–$None$ and for any fixed grid point
$t^{*}\in[t_{0},t_{\emph{final}}]$, the sequences (3)–(5) satisfy

$r_{n}^{[1]}(h)-y(t^{*})={\cal{O}}(h^{p+1}),$
$r_{n}^{[2]}(h)-y(t^{*})={\cal{O}}(h^{p+2}),$
$r_{n}^{[3]}(h)-y(t^{*})={\cal{O}}(h^{p+3}).$

###### Proof.

Extrapolation techniques are based on the existence of an asymptotic expansion
of the global error with respect to the small parameter $h$. We will show in
the second half of the proof that [7, Section III.9, Theorem 9.1] is
applicable in the present situation. This theorem then guarantees the
existence of some functions $\mathbf{e}_{p}$,
$\mathbf{e}_{p+1}$,$\mathbf{e}_{p+2}$ and $\mathbf{E}$ such that

$y_{n}(h)-y(t^{*})=\mathbf{e}_{p}(t^{*})\cdot
h^{p}+\mathbf{e}_{p+1}(t^{*})\cdot h^{p+1}+\mathbf{e}_{p+2}(t^{*})\cdot
h^{p+2}+\mathbf{E}(t^{*},h)h^{p+3},$ (7)

where, given any $\delta>0$ small enough, the function $\mathbf{E}$ is
uniformly bounded for any $h\in[0,h_{0}]$ and any
$t^{*}\in[t_{0}+\delta,t_{\text{final}}]$. (Clearly, for fixed $t^{*}$, we can
choose $\delta$ such that $t^{*}\in[t_{0}+\delta,t_{\text{final}}]$.) The
linear combinations in the definition of the sequences (3)–(5) have been set
up such that the first, first two, and the first three terms, respectively, on
the right-hand side of (7) are eliminated. More precisely, to construct the
coefficients of $\eqref{ketszeresRE}$, for example, one considers the
expression

$\gamma_{1}\cdot(y_{n}(h)-y(t^{*}))+\gamma_{2}\cdot\left(y_{2n}\left(\frac{h}{2}\right)-y(t^{*})\right)+\gamma_{3}\cdot\left(y_{4n}\left(\frac{h}{4}\right)-y(t^{*})\right),$

then applies (7) with $h$, $h/2$ and $h/4$, and solves the linear system for
$\gamma_{1,2,3}$ obtained by setting $\gamma_{1}+\gamma_{2}+\gamma_{3}$ equal
to 1, the coefficient of $h^{p}$ equal to 0, and the coefficient of $h^{p+1}$
equal to 0.

To finish the proof, what remains is to show that our assumptions
$None$–$None$ on the LMM imply the assumptions of [7, Section III.9, Theorem
9.1]. Theorem 9.1 is a general theorem about the asymptotic expansion of the
global error of strictly stable general linear methods (GLMs). It is known [7,
Section III.4] that any LMM (2) can be interpreted as a GLM written as a one-
step method in a higher dimensional space as follows:

$Y_{n+1}=(A\otimes I)Y_{n}+h\Phi(t_{n},Y_{n},h)\quad\quad(n\in\mathbb{N}),$
(8)

where $Y_{n}:=(y_{n+k-1},y_{n+k-2},\ldots,y_{n})^{\top}$ ($n\in\mathbb{N}$),
$\Phi$ is the increment function of the numerical method,
$A\in\mathbb{R}^{k\times k}$ is the matrix (6), $I\in\mathbb{R}^{m\times m}$
is the identity matrix (recall that $m$ is the dimension of (1), hence
$A\otimes I\in\mathbb{R}^{m\cdot k\times m\cdot k}$). Then the matrix $S$
appearing in [7, Section III.9, Theorem 9.1] is $S:=A\otimes I$.
Step 1. Our assumption $None$ implies the smoothness assumption A3 of Theorem
9.1.
Step 2. Assumption A2 of Theorem 9.1 holds automatically since in our case the
increment function $\Phi$ does not depend on $n$.
Step 3. Assumptions $None$–$None$ imply the consistency of order $p$ of the
GLM. Indeed, assumption $None$ implies that the local error of the LMM is
${\cal{O}}(h^{p+1})$ in the sense of [7, Section III.2, Definition 2.3]. But
this quantity is just the first component of the vector
$Y(x_{n+1})-\hat{Y}_{n+1}$ in [7, Section III.4, Lemma 4.3] describing the
local error of the one-step reformulation (8) of the LMM; the remaining
components of $Y(x_{n+1})-\hat{Y}_{n+1}$ are $0$. This means that the local
error $d_{n+1}$ ($n\in\mathbb{N}$) of the GLM in the sense of [7, Section
III.8, Definition 8.9] is also ${\cal{O}}(h^{p+1})$. Since
$d_{0}={\cal{O}}(h^{p})$ due to assumption $None$, we get, by using [7,
Section III.8, Lemma 8.11], that the consistency order of the GLM is $p$. (In
(8.15) of [7, Section III.8, Lemma 8.11], the relation with the spectral
projector $E\delta_{p}(x)=0$ now clearly holds, since $\delta_{p}(x)=0$ due to
$d_{n+1}={\cal{O}}(h^{p+1})$.)
Step 4. We finally show that assumption A1 of Theorem 9.1 follows from $None$.
We need to show that the GLM is strictly stable, that is,

  1. (i)

the matrix $S$ is power bounded, in other words,
$\sup_{n\in\mathbb{N}}\|S^{n}\|\in\mathbb{R}$;

  2. (ii)

1 is the only eigenvalue of $S$ with modulus one.

It is known [2, Chapter 7] that the eigenvalues of a Kronecker product
$S=A\otimes I$ are the product of the eigenvalues of its components, hence
(ii) is verified. Finally, to check (i) we recall (see, e.g., [8, Chapter 7,
Theorem 2.2.5]) that a square matrix $S$ is power bounded if and only if

  * (ia)

each zero of its minimal polynomial lies in the closed unit disk, and

  * (ib)

(potential) zeros of its minimal polynomial on the unit circle have
multiplicity 1.

The minimal polynomial divides the characteristic polynomial, so $None$ imples
(ia). We also know from $None$ that $1$ is a simple zero of the characteristic
polynomial of $A$, but from the above properties of the Kronecker product we
see that its multiplicity in the _characteristic_ polynomial of $S$ is $m$.
However, $1$ is also a simple zero of the minimal polynomial of $A$. Hence
Lemma 2.3 below with $M:=A$ implies property (ib). ∎

The following auxiliary result, which we state separately, has been used in
the above proof.

######  Lemma 2.3.

Let $M\in\mathbb{C}^{k\times k}$ be any square matrix, and let
$I\in\mathbb{R}^{m\times m}$ be the identity matrix ($k,m\in\mathbb{N}^{+}$).
Then the minimal polynomials of $M$ and $M\otimes I$ coincide.

###### Proof.

In the proof, we denote identity matrices of appropriate size by the same
symbol $I$. Consider a Jordan canonical form of $M=TJT^{-1}$. Then one easily
checks that $M\otimes I=(T\otimes I)(J\otimes I)(T\otimes I)^{-1}$. We know
that if two matrices are similar, then their minimal polynomials coincide. So
to prove the lemma, it is sufficient to show that
$\mathfrak{M}_{J}=\mathfrak{M}_{J\otimes I}$.

Let $J_{1},\ldots,J_{n}$ denote the blocks of the block diagonal matrix $J$.
Then the blocks of the block diagonal matrix $J\otimes I$ are $J_{1}\otimes
I,\ldots,J_{n}\otimes I$. It is easily seen that

$\mathfrak{M}_{J}=\text{lcm}\left(\mathfrak{M}_{J_{1}},\mathfrak{M}_{J_{2}},\ldots,\mathfrak{M}_{J_{n}}\right),$

where lcm denotes the least common multiple of the given polynomials with
leading coefficient set to $1$. Similarly, we have

$\mathfrak{M}_{J\otimes I}=\text{lcm}\left(\mathfrak{M}_{J_{1}\otimes
I},\ldots,\mathfrak{M}_{J_{n}\otimes I}\right).$

Finally, we claim that $\mathfrak{M}_{J_{j}}=\mathfrak{M}_{J_{j}\otimes I}$
for any block $J_{j}\in\mathbb{R}^{s_{j}\times s_{j}}$ ($j=1,2,\ldots,n$).
Indeed, let $\lambda$ denote the independent variable of the minimal
polynomial, $\lambda_{j}$ the eigenvalue corresponding to $J_{j}$, and $N$ the
appropriate sized (nilpotent) square matrix with $1$s in its (first)
superdiagonal and $0$s everywhere else. Then, for both possible types of
Jordan blocks, we have

  * •

$J_{j}=\lambda_{j}I\Longrightarrow\mathfrak{M}_{J_{j}}=\lambda-\lambda_{j}=\mathfrak{M}_{J_{j}\otimes
I}$,

  * •

$J_{j}=\lambda_{j}I+N\Longrightarrow\mathfrak{M}_{J_{j}}=(\lambda-\lambda_{j})^{s_{j}}=\mathfrak{M}_{J_{j}\otimes
I}$,

completing the proof. ∎

######  Remark 2.4.

Power boundedness of the matrix $S=A\otimes I$ in the above proof of Theorem
2.2 also follows from [7, Section III.4, Lemma 4.4], since now $A$ is power
bounded. Hence our Lemma 2.3 can be considered as an alternative (more
algebraic) proof of this fact.

######  Remark 2.5.

[5, Lemma 2.1] contains the convergence proof of the sequence (3) based on [6,
Section 6.3.4, Theorem 6.3.5]. This is the reason we have assumed in [5, Lemma
2.1] that the $k$ starting values of the sequence
$\\{y_{0}(h),y_{1}(h),\ldots,y_{k-1}(h)\\}$ to initiate the LMM are each
${\cal{O}}(h^{p+1})$-close to the corresponding exact solution values. The
present discussion shows that it is in fact enough to assume the weaker
closeness assumption $None$.

######  Remark 2.6.

In [14, Theorem 4.1], the authors did not explicitly reference a result that
guarantees the existence of an asymptotic expansion of the global error for
Runge–Kutta methods, on which expansion their proof is based. One such
classical theorem for one-step methods is due to Gragg (1964); for a modern
treatment, see, e.g., [7, Section II.8, Theorem 8.1].

Many LMMs satisfy assumptions $None$ and $None$. In particular, all AB and AM
methods satisfy $None$ (since, for any of these methods, the only non-zero
eigenvalue of the corresponding matrix (6) is $1$ with algebraic multiplicity
$1$), and one directly checks that each BDF$k$ method also satisfies $None$.
Therefore, we have the following convergence result.

######  Corollary 2.7.

According to Theorem 2.2, the convergence of the sequences obtained by
applying any of the LMMs

  * •

AB$k$ (with $k\geq 2$ steps),

  * •

AM$k$ (with $k\geq 2$ steps), or

  * •

BDF$k$ (with $2\leq k\leq 6$ steps)

to problem (1) can be accelerated by using formulae (3)–(5).

These combined methods will be referred to as AB$k$-RGRE, AM$k$-RGRE, and
BDF$k$-RGRE in general, and AB$k$-$\ell$GRE, AM$k$-$\ell$GRE, and
BDF$k$-$\ell$GRE in particular (where $\ell$ is the number of RE applications
appearing in $r_{n}^{[\ell]}(h)$ in (3)–(5)).

##  3 Linear stability analysis

In this section, let LMM denote any of the base methods AB$k$ ($k\geq 2$),
AM$k$ ($k\geq 2$), or BDF$k$ ($2\leq k\leq 6$), and let
$\ell\in\mathbb{N}^{+}$ denote the number of RE applications. The region of
absolute stability of the base method is denoted by
$\mathcal{S}_{\text{LMM}}$, while that of the combined LMM-RGRE method will be
denoted by ${\mathcal{S}}_{\text{RGRE}}^{[\ell]}$. Let us apply the LMM-RGRE
method to the scalar linear test equation $y^{\prime}(t)=\lambda y(t)$,
$y(t_{0})=y_{0}$ with some $h>0$ and $\lambda\in\mathbb{C}$. Similarly to [5,
Section 2.3], we define
${\mathcal{S}}_{\text{RGRE}}^{[\ell]}\subset\mathbb{C}$ as

the set of numbers $\mu:=h\lambda$ for which the sequence $n\mapsto
r_{n}^{[\ell]}(h)$ is bounded for any choice of the starting values of any of
its component sequences $n\mapsto y_{n}\left({h}/{2^{j}}\right)$
($j=1,\ldots,\ell$), but excluding the values of $\mu$ for which any of the
leading coefficients $\alpha_{k}-\left({\mu}/{2^{j}}\right)\beta_{k}$
($j=1,\ldots,\ell$) vanishes.

The proofs of the lemmas of [5, Section 2.3] can be extended in a
straightforward way to the present situation, so we only state the results.

######  Lemma 3.1.

We have the inclusions

$\bigcap_{j=0}^{\ell}\big{(}{2^{j}\mathcal{S}}_{\text{\emph{LMM}}}\big{)}\subseteq{\mathcal{S}}_{\text{RGRE}}^{[\ell]}\subseteq{\mathcal{S}}_{\text{\emph{LMM}}}.$

######  Lemma 3.2.

For any $2\leq k\leq 6$ and $\ell\in\mathbb{N}^{+}$, the
$\text{BDF}k$-$\ell$GRE method has the same $A(\alpha)$-stability angle as
that of the underlying $\text{BDF}k$-method.

######  Lemma 3.3.

Assume that ${\mathcal{S}}_{\text{\emph{LMM}}}$ is convex. Then
${\mathcal{S}}_{\text{RGRE}}^{[\ell]}={\mathcal{S}}_{\text{\emph{LMM}}}$.

######  Remark 3.4.

In [5], we have shown that ${\mathcal{S}}_{\text{\emph{LMM}}}$ is convex, for
example, for the AB2 and AM2 methods, but not convex for the AB3 method.

##  4 Numerical tests

In this section, we verify the expected order of convergence of LMM-RGREs on
three benchmark problems: we chose a Dahlquist test problem

$y^{\prime}(t)=-5y(t),\quad y(0)=1,$

a Lotka–Volterra system

$y^{\prime}_{1}(t)=0.1y_{1}(t)-0.3y_{1}(t)y_{2}(t),\quad
y^{\prime}_{2}(t)=0.5(y_{1}(t)-1)y_{2}(t)$

for $t\in[0,62]$ with initial condition $y(0)=(1,1)^{\top}$, and a mildly
stiff van der Pol equation

$y^{\prime}_{1}(t)=y_{2}(t),\quad\quad
y^{\prime}_{2}(t)=2(1-y_{1}^{2}(t))y_{2}(t)-y_{1}(t)$

for $t\in[0,20]$ with initial condition $y(0)=(2,0)^{\top}$.

As base LMMs, we considered the $2^{\text{nd}}$\- and $3^{\text{rd}}$-order
AB, AM, and BDF methods. For starting methods, we chose the $2^{\text{nd}}$\-
and $3^{\text{rd}}$-order Ralston methods, having minimum error bounds [10].
As it is usual, the AM methods were implemented in predictor-corrector style.
For the nonlinear algebraic equations arising in connection with implicit
LMMs, we use MATLAB’s fsolve command. When the goal is to achieve high
convergence order or when we apply RE several times, one should change to a
multi-precision environment due to the double precision limitation of MATLAB
and to the limitations of fsolve. The fine-grid solutions obtained by a
$6^{\text{th}}$-order Runge–Kutta method with $2^{16}$ grid points are used as
a reference solution to measure the global error in maximum norm and to
estimate the order of convergence of LMM-RGREs [9, Appendix A].

Table 1 and Figure 1 illustrate the expected $4^{\text{th}}$-order of
convergence for all tested LMM-RGREs.

Dahlquist problem | Lotka–Volterra system
---|---
$\text{AB}2$-2GRE |  $\text{AM}2$-2GRE |  $\text{BDF}2$-2GRE |  $\text{AB}2$-2GRE |  $\text{AM}2$-2GRE |  $\text{BDF}2$-2GRE
3.9437 | 4.0923 | 4.4613 | 3.9707 | 4.3757 | 3.8136
3.9833 | 4.0469 | 4.2500 | 3.9883 | 4.2391 | 3.9265
3.9942 | 4.0237 | 4.1333 | 3.9951 | 4.1376 | 3.9679
3.9977 | 4.0119 | 4.0342 | 3.9983 | 4.0928 | 3.9908
Table 1: The estimated order of convergence for the Dahlquist problem and the
Lotka–Volterra system for different $4^{\text{th}}$-order LMM-RGREs with
$64,128,\ldots,1024$ and $512,1024,\ldots,8192$ grid points, respectively
Figure 1: Results for the van der Pol equation for LMM-2GREs, number of grid
points versus the global error in maximum norm

Tables 2 and 3, and Figure 2 illustrate the expected $5^{\text{th}}$-order of
convergence for all tested LMM-RGREs.

$\text{AB}3$-2GRE |  $\text{AM}3$-2GRE |  $\text{BDF}3$-2GRE |  $\text{AB}2$-3GRE |  $\text{AM}2$-3GRE |  $\text{BDF}2$-3GRE
---|---|---|---|---|---
5.0121 | 5.0319 | 5.2081 | 4.8145 | 5.0564 | 5.1449
Table 2: The estimated order of convergence for the Dahlquist problem for different $5^{\text{th}}$-order LMM-RGREs with $256$ and $512$ grid points $\text{AB}3$-2GRE |  $\text{AM}3$-2GRE |  $\text{BDF}3$-2GRE |  $\text{AB}2$-3GRE |  $\text{AM}2$-3GRE |  $\text{BDF}2$-3GRE
---|---|---|---|---|---
4.4933 | 4.8252 | 5.4540 | 6.1243 | 4.2922 | 4.0779
5.0397 | 4.9769 | 5.4343 | 5.4060 | 4.8533 | 4.7738
5.0747 | 4.9996 | 5.2692 | 5.2996 | 4.9615 | 4.9942
4.9983 | 4.9509 | 5.2278 | 5.2431 | 4.9864 | 4.9914
Table 3: The estimated order of convergence for the Lotka–Volterra system for
different $5^{\text{th}}$-order LMM-RGREs with $256,512,\ldots,4096$ grid
points Figure 2: Results for the van der Pol equation for
$5^{\text{th}}$-order LMM-RGREs, number of grid points versus the global error
in maximum norm

## Acknowledgement

I. Fekete was supported by the János Bolyai Research Scholarship of the
Hungarian Academy of Sciences. Supported by the ÚNKP-22-5 New National
Excellence Program of the Ministry for Culture and Innovation from the source
of the National Research, Development and Innovation Fund. The research was
supported by the Ministry of Innovation and Technology NRDI Office within the
framework of the Artificial Intelligence National Laboratory Program
RRF-2.3.1-21-2022-00004.

## References

  * [1] T. Bayleyegn, Á. Havasi, Multiple Richardson Extrapolation Applied to Explicit Runge–Kutta Methods, Advances in High Performance Computing. Springer, 902, 262–270, (2020), http://doi.org/10.1007/978-3-030-55347-0_22
  * [2] D. S. Bernstein, Matrix Mathematics, Princeton Univ. Press (2009), https://doi.org/10.1515/9781400833344
  * [3] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd edition, Wiley (2016), https://doi.org/10.1002/9781119121534
  * [4] R. D. Falgout, T. A. Manteuffel, B. O’Neill, J. B. Schroder, Multigrid reduction in time with Richardson extrapolation, Electron. Trans. Numer. Anal., 54, 210–233 (2021), https://doi.org/10.1553/etna_vol54s210
  * [5] I. Fekete, L. Lóczi, Linear multistep methods and global Richardson extrapolation, Appl. Math. Letters 133 (2022), https://doi.org/10.1016/j.aml.2022.108267
  * [6] W. Gautschi, Numerical Analysis, 2nd edition, Birkhäuser (2012), https://doi.org/10.1007/978-0-8176-8259-0
  * [7] E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, 2nd revised edition, Springer-Verlag (1993), https://doi.org/10.1007/978-3-540-78862-1
  * [8] Z. Jackiewicz, General Linear Methods for Ordinary Differential Equations, Wiley (2009), https://doi.org/10.1002/9780470522165
  * [9] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM, Philadelphia (2007), https://doi.org/10.1137/1.9780898717839
  * [10] A. Ralston, Runge–Kutta methods with minimum error bounds, Math. Comp., 16, 431–437 (1962), https://doi.org/10.1090/S0025-5718-1962-0150954-0
  * [11] L. F. Richardson, The approximate arithmetical solution by finite differences of physical problems including differential equations, with an application to the stresses in a masonry dam, Philos. Trans. Roy. Soc. London Ser. A., 210, 307–357 (1911), https://doi.org/10.1098/rsta.1911.0009
  * [12] L. F. Richardson, The deferred approach to the limit, Philos. Trans. Roy. Soc. London Ser. A., 226, 299–361 (1927), https://doi.org/10.1098/rsta.1927.0008
  * [13] Z. Zlatev, I. Dimov, I. Faragó, Á. Havasi, Richardson Extrapolation: Practical Aspects and Applications, Berlin, Boston: De Gruyter (2017), https://doi.org/10.1515/9783110533002
  * [14] Z. Zlatev, I. Dimov, I. Faragó, K. Georgiev, Á. Havasi, Explicit Runge–Kutta Methods Combined with Advanced Versions of the Richardson Extrapolation, Comput. Methods Appl. Math. 20(4) 739–762 (2020), https://doi.org/10.1515/cmam-2019-0016
  * [15] Z. Zlatev, I. Dimov, I. Faragó, Á. Havasi, Efficient implementation of advanced Richardson extrapolation in an atmospheric chemical scheme, J. Math. Chem., 60 (1), 219–238 (2022), https://doi.org/10.1007/s10910-021-01300-z


# Open Quantum Systems in Noninertial Frames

Salman Khan and M. K. Khan
Department of Physics, Quaid-i-Azam University,
Islamabad 45320, Pakistan<EMAIL_ADDRESS>
###### Abstract

We study the effects of decoherence on the entanglement generated by Unruh
effect in noninertial frames by using bit flip, phase damping and depolarizing
channels. It is shown that decoherence strongly influences the initial state
entanglement. The entanglement sudden death can happens irrespective of the
acceleration of the noninertial frame under the action of phase flip and phase
damping channels. It is investigated that an early sudden death happens for
large acceleration under the depolarizing environment. Moreover, the
entanglement increases for a highly decohered phase flip channel.
PACS: 03.65.Ud; 03.65.Yz; 03.67.Mn;04.70.Dy

Keywords: Entanglement; Decoherence; Noninertial frames.

##  1 Introduction

Entanglement is one of the potential sources of quantum theory. It is the key
concept and major resource for quantum communication and computation [1]. In
the last few years, enormous efforts has been made to investigate various
aspects of quantum entanglement and its benefits in a number of setups, such
as teleportation of unknown states [2] , quantum key distribution [3], quantum
cryptography [4] and quantum computation [5, 6]. Recently, the study of
quantum entanglement of various fields has been extended to the relativistic
setup [7, 8, 9, 10, 11, 12] and interesting results about the behavior of
entanglement have been obtained. The study of entanglement in the relativistic
framework is important not only from quantum information perspective but also
to understand deeply the black hole thermodynamics [13, 14] and the black hole
information paradox [15, 16].

The earlier investigations on quantum entanglement in the relativistic
framework is mainly focused by considering isolated quantum systems. In fact,
no quantum system can be completely isolated from its environment and may
results in a non-unitary dynamics of the system. Therefore, it is important to
study the effect of environment on the entanglement in an initial state of a
quantum system during its evolution. The interaction between an environment
and a quantum system leads to the phenomenon of decoherence and it gives rise
to an irreversible transfer of information from the system to the environment
[17, 18, 19].

\put(-350.0,220.0){} |  |
---|---|---
Figure 1: Rindler spacetime diagram: A uniformly accelerated observer Rob (R)
moves on a hyperbola with acceleration $a$ in region $I$ and is causally
disconnected from region $II$.

In this paper we work out the effect of decoherence on the entanglement of
Dirac field in a noninertial system. Alsing et al [7] have shown that the
entanglement between two modes of a free Dirac field is degraded by the Unruh
effect and asymptotically reaches a nonvanishing minimum value in the infinite
acceleration. We investigate that how the loss of entanglement through Unruh
effect is influenced in the presence of decoherence by using a phase flip, a
phase damping and a depolarizing channel in the Kraus operators formalism. The
effect of amplitude damping channel on Dirac field in a noninertial system is
recently studied by Wang and Jing [20]. We consider two observers, Alice and
Rob, that share a maximally entangled initial state of two qubits at a point
in flat Minkowski spacetime. Then Rob moves with a uniform acceleration and
Alice stays stationary. To achieve our goal, we consider two cases. In one
instance we allow only Rob’s qubit to interact with a noisy environment and in
the second instance both qubits of the two observers interact with a noisy
environment. Let the two modes of Minkowski spacetime that correspond to Alice
and Rob are, respectively, given by $|n\rangle_{A}$ and $|n\rangle_{R}$.
Moreover, we assume that the observers are equipped with detectors that are
sensitive only to their respective modes and share the following maximally
entangled initial state

$|\psi\rangle_{A,R}=\frac{1}{\sqrt{2}}\left(|00\rangle_{A,R}+|11\rangle_{A,R}\right),$
(1)

where the first entry in each ket corresponds to Alice and the second entry
corresponds to Rob. From the accelerated Rob’s frame, the Minkowski vacuum
state is found to be a two-mode squeezed state [7],

$|0\rangle_{M}=\cos r|0\rangle_{I}|0\rangle_{II}+\sin
r|1\rangle_{I}|1\rangle_{II},$ (2)

where $\cos r=\left(e^{-2\pi\omega c/a}+1\right)^{-1/2}$. The constant
$\omega$, $c$ and $a$, in the exponential stand, respectively, for Dirac
particle’s frequency, light’s speed in vacuum and Rob’s acceleration. In Eq.
(2) the subscripts $I$ and $II$ of the kets represent the Rindler modes in
region $I$ and $II$, respectively, in the Rindler spacetime diagram (see Fig.
(1)). The excited state in Minkowski spacetime is related to Rindler modes as
follow [7]

$|1\rangle_{M}=|1\rangle_{I}|0\rangle_{II}.$ (3)

In terms of Minkowski modes for Alice and Rindler modes for Rob, the maximally
entangled initial state of Eq. (1) by using Eqs. (2) and (3) becomes

$|\psi\rangle_{A,I,II}=\frac{1}{\sqrt{2}}\left(\cos
r|0\rangle_{A}|0\rangle_{I}|0\rangle_{II}+\sin
r|0\rangle_{A}|1\rangle_{I}|1\rangle_{II}+|1\rangle_{A}|1\rangle_{I}|0\rangle_{II}\right).$
(4)

Since Rob is causally disconnected from region $II$, we must take trace over
all the modes in region $II$. This leaves the following mixed density matrix
between Alice and Rob, that is,

$\displaystyle\rho_{A,I}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\cos^{2}r|00\rangle_{A,I}\langle 00|+\cos
r(|00\rangle_{A,I}\langle 11|+|11\rangle_{A,I}\langle 00|)$ (5)
$\displaystyle\sin^{2}r|01\rangle_{A,I}\langle 01|+|11\rangle_{A,I}\langle
11|].$

Table 1: A single qubit Kraus operators for phase flip channel, phase damping channel and depolarizing channel.  phase flip channel | $E_{o}=\sqrt{1-p}\left(\begin{array}[]{cc}1&0\\\ 0&1\end{array}\right),\qquad E_{1}=\sqrt{p}\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right)$
---|---
phase damping channel | $E_{o}=\left(\begin{array}[]{cc}1&0\\\ 0&\sqrt{1-p}\end{array}\right),\qquad E_{1}=\left(\begin{array}[]{cc}0&0\\\ 0&\sqrt{p}\end{array}\right)$
depolarizing channel | $\begin{array}[]{c}E_{o}=\sqrt{1-p}\left(\begin{array}[]{cc}1&0\\\ 0&1\end{array}\right),\qquad E_{1}=\sqrt{p/3}\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right),\\\ E_{2}=\sqrt{p/3}\left(\begin{array}[]{cc}0&-i\\\ i&0\end{array}\right),\qquad E_{3}=\sqrt{p/3}\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right)\end{array}$

[htb]

##  2 single qubit in a noisy environment

In this section we consider that only the Rob’s qubit is coupled to a noisy
environment. The final density matrix of the system in the Kraus operators
representation becomes

$\rho_{f}=\sum_{i}\left(\sigma_{o}\otimes
E_{i}\right)\rho_{A,I}\left(\sigma_{o}\otimes E_{i}^{{\dagger}}\right),$ (6)

where $\rho_{A,I}$ is the initial density matrix of the system given by Eq.
(5), $\sigma_{o}$ is the single qubit identity matrix and $E_{i}$ are a single
qubit Kraus operators of the channel under consideration. The Kraus operators
of the channels we use are given in Table $1$. The spin-flip matrix of the
final density matrix of Eq. (6) is defined as
$\tilde{\rho}_{f}=\left(\sigma_{2}\otimes\sigma_{2}\right)\rho_{f}\left(\sigma_{2}\otimes\sigma_{2}\right)$,
where $\sigma_{2}$ is the Pauli matrix. The degree of entanglement in the two
qubits mixed state in a noisy environment can be quantified conveniently by
concurrence $C$, which is given as [21, 22]

$C=\max\left\\{0,\sqrt{\lambda_{1}}-\sqrt{\lambda_{2}}-\sqrt{\lambda_{3}}-\sqrt{\lambda_{4}}\right\\}\qquad\lambda_{i}\geq\lambda_{i+1}\geq
0,$ (7)

where $\lambda_{i}$ are the eigenvalues of the matrix
$\rho_{f}\tilde{\rho}_{f}$. The eigenvalues under the action of phase-flip
channel becomes

\put(-350.0,220.0){} |  |
---|---|---
Figure 2: The concurrence $C$ under the action of phase flip channel is
plotted against decoherence parameter $p$ for the case when only Rob’s qubit
is coupled to a noisy environment.

$\displaystyle\lambda_{1}^{\mathrm{PF}}$ $\displaystyle=$
$\displaystyle(1-2p+p^{2})\cos^{2}r,$ $\displaystyle\lambda_{2}^{\mathrm{PF}}$
$\displaystyle=$ $\displaystyle p^{2}\cos^{2}r,$
$\displaystyle\lambda_{3}^{\mathrm{PF}}$ $\displaystyle=$
$\displaystyle\lambda_{4}^{\mathrm{PF}}=0,$ (8)

where the superscript PF corresponds to phase flip channel. Similarly, the
eigenvalues under the action of phase damping and depolarizing channels are,
respectively, given by

$\displaystyle\lambda_{1,2}^{\mathrm{PD}}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(2-p\pm 2\sqrt{1-p})\cos^{2}r,$
$\displaystyle\lambda_{3}^{\mathrm{PD}}$ $\displaystyle=$
$\displaystyle\lambda_{4}^{\mathrm{PD}}=0,$ (9)
$\displaystyle\lambda_{1}^{\mathrm{DP}}$ $\displaystyle=$
$\displaystyle(-1+p)^{2}\cos^{2}r,$ $\displaystyle\lambda_{2}^{\mathrm{DP}}$
$\displaystyle=$
$\displaystyle\lambda_{3}^{\mathrm{DP}}=\lambda_{4}^{\mathrm{DP}}=\frac{1}{9}p^{2}\cos^{2}r,$
(10)

where the superscripts PD and DP stand for phase damping and depolarizing
channels, respectively. In all these equations $p\in\left[0,1\right]$ is the
decoherence parameter. The upper and lower values of $p$ correspond to
undecohered and fully decohered case of the channels, respectively. The
concurrence under the action of every channel reduces to the result of Ref.
[7] when the decoherence parameter $p=0$.

To see how the concurrence and hence the entanglement is influenced by
decoherence parameter $p$ in the presence of Unruh effect, we plot the
concurrence for each channel against $p$ for various values of $r$. In Fig.
(2), the concurrence under the action of phase flip channel is plotted against
$p$. The figure shows that for smaller values of $p$, the entanglement is
strongly acceleration dependent, such that for large values of Rob’s
acceleration (the value of $r$) it gets weakened. However, as $p$ increases
the dependence of entanglement on acceleration decreases and the increasing
value of $p$ causes a rapid loss of entanglement.

\put(-350.0,220.0){} |  |
---|---|---
Figure 3: The concurrence $C$ under the action of phase damping channel is
plotted against decoherence parameter $p$ for the case when only Rob’s qubit
is coupled to a noisy environment.

The entanglement sudden death happens irrespective of the acceleration of
Rob’s frame for a $50\%$ decoherence. Fig. (3) shows the effect of decoherence
on the concurrence under the action of phase damping channel. In this case,
the degradation of entanglement due to decoherence is smaller as compare to
the the degradation in the case of phase flip. The entanglement vanishes for
all values of acceleration only when the channel is fully decohered. The
concurrence under the action of the depolarizing channel is exactly equal to
the one for phase flip channel. Hence it influences the entanglement in a way
exactly similar to the phase flip channel as shown in Fig. ($2$).

##  3 Both qubits in a noisy environment

In this section we consider that both Alice’s and Rob’s qubits are influenced
simultaneously by a noisy environment. The final density matrix in this case
can be written in the Kraus operators formalism as follows

$\rho_{f}=\sum_{k}E_{k}\rho_{A,I}E_{k}^{{\dagger}},$ (11)

where $\rho_{A,I}$ is given by Eq. (5) and $E_{k}$ are the Kraus operators for
a two qubit system, satisfying the completeness relation
$\sum_{k}E_{k}E_{k}=I$ and are constructed from a single qubit Kraus operators
of a channel by taking tensor product of all the possible combinations in the
following way

$E_{k}=\sum_{i,j}E_{i}\otimes E_{j},$ (12)

where $E_{i,j}$ are the single qubit Kraus operators of a channel given in
Table $1$. We consider that both Alice’s and Bob’s qubits are influenced by
the same environment, that is, the decoherence parameter $p$ for both qubits
is same. Proceeding in a similar way like the case of single qubit coupled to
the environment, the eigenvalues of the matrix $\rho_{f}\tilde{\rho}_{f}$
under the action of phase flip channel become

$\displaystyle\lambda_{1}^{\mathrm{PF}}$ $\displaystyle=$
$\displaystyle(1+2(-1+p)p)^{2}\cos^{2}r,$
$\displaystyle\lambda_{2}^{\mathrm{PF}}$ $\displaystyle=$ $\displaystyle
4(-1+p)^{2}p^{2}\cos^{2}r,$ $\displaystyle\lambda_{3}^{\mathrm{PF}}$
$\displaystyle=$ $\displaystyle\lambda_{4}^{\mathrm{PF}}=0,$ (13)

Likewise the eigenvalues for phase damping and depolarizing channels,
respectively, becomes

$\displaystyle\lambda_{1}^{\mathrm{PD}}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(-2+p)^{2}\cos^{2}r,$
$\displaystyle\lambda_{2}^{\mathrm{PD}}$ $\displaystyle=$
$\displaystyle\frac{1}{4}p^{2}\cos^{2}r,$
$\displaystyle\lambda_{3}^{\mathrm{PD}}$ $\displaystyle=$
$\displaystyle\lambda_{4}^{\mathrm{PD}}=0,$ (14)

$\displaystyle\lambda_{1,3}^{\mathrm{DP}}$ $\displaystyle=$
$\displaystyle\frac{1}{1296}[324+p(-3+2p)(387+152p(-3+2p))$
$\displaystyle+4(3-4p)^{2}(9+5p(-3+2p))\cos 2r$
$\displaystyle+(3-4p)^{2}p(-3+2p)\cos 4r\pm 4(3-4p)^{2}\cos r$
$\displaystyle\times\\{3(54+p(-3+2p)(33+8p(-3+2p)))$
$\displaystyle+(3-4p)^{2}(2(9-6p+4p^{2})\cos 2r+p(-3+2p)\cos 4r)\\}^{1/2}],$
$\displaystyle\lambda_{2}^{\mathrm{DP}}$ $\displaystyle=$
$\displaystyle\lambda_{4}^{\mathrm{DP}}=\frac{1}{648}p(-3+2p)(-9+4p$ (15)
$\displaystyle+(-3+4p)\cos 2r)(3+4p+(-3+4p)\cos 2r),$

The $"\pm"$ sign in Eq. (15), correspond to the eigenvalues $\lambda_{1}$, and
$\lambda_{3}$ respectively. It is necessary to point out here that the
concurrence under the action of each channel reduces to the result of Ref. [7]
when we set the decoherence parameter $p=0$.

\put(-350.0,220.0){} |  |
---|---|---
Figure 4: The concurrence $C$ under the action of phase flip channel is
plotted against decoherence parameter $p$ for the case when both qubits are
coupled to a noisy environment.

To see how the entanglement behaves when both the qubits are coupled to the
noisy environment, we plot the concurrence against the decoherence parameter
$p$ for different values of $r$ under the action of each channel separately.
Fig. (4) shows the dependence of concurrence on decoherence parameter $p$
under the action of phase flip channel. The dependence of entanglement on
acceleration of Rob’s frame is obvious in the region of lower values of $p$.
However, this dependence diminishes as $p$ increases and a rapid decrease in
the degree of entanglement develops. At a $50\%$ decoherence level, the
entanglement sudden death occurs irrespective of Rob’s acceleration. It’s
interesting to see that beyond this point onward, the entanglement regrows as
$p$ increases. The dependence of entanglement on acceleration of the Rob’s
frame reemerges and the entanglement reaches to the corresponding undecohered
maximum value for a fully decohered case. The concurrence varies as a
parabolic function of decoherence parameter $p$ with its vertex at $p=0.5$.

\put(-350.0,220.0){} |  |
---|---|---
Figure 5: The concurrence $C$ under the action of phase damping channel is
plotted against decoherence parameter $p$ for the case when both qubits are
coupled to a noisy environment.

The dependence of entanglement on $p$ under the action of phase damping
channel is shown in Fig. (5). In this case the entanglement decreases linearly
as $p$ increases and the dependence on acceleration diminishes. Whatever the
acceleration of Rob’s frame may be, the entanglement sudden death occurs when
the channel is fully decohered. The influence of depolarizing channel on the
entanglement is shown in Fig. (6). Unlike the other two channels, the
depolarizing channel does not diminish the effect of acceleration on the
entanglement as the $p$ increases. However a rapid decrease in entanglement
appears which leads to entanglement sudden death at different values of
decoherence parameter for different acceleration of Rob’s frame. The larger
the acceleration the earlier the entanglement sudden death occurs.

\put(-350.0,220.0){} |  |
---|---|---
Figure 6: The concurrence $C$ under the action of depolarizing channel is
plotted against decoherence parameter $p$ for the case when both qubits are
coupled to a noisy environment.

If we compare the single qubit and the both qubits decohering situations, it
becomes obvious that the entanglement loss is rapid when both the qubits are
coupled to the noisy environment. For example, in the case of bit flip channel
the concurrence behaves as a linear function of $p$ for single qubit
decohering case whereas in the case of both qubits decohering case it varies
as a parabolic function. Nevertheless, the sudden death happens at the same
value of $p$, irrespective of the acceleration, for both cases under the
action of bit flip and phase damping channels. For depolarizing channel,
however, this is not true.

##  4 Conclusion

In conclusion, we have investigated that the entanglement in Dirac fields is
strongly dependent on coupling with a noisy environment. This result is
contrary to the case of an isolated system in which the entanglement of Dirac
fields survives even in the limit of infinite acceleration of Rob’s frame. In
the presence of decoherence, the entanglement rapidly decreases and
entanglement sudden death occurs even for zero acceleration. Under the action
of phase flip channel, the entanglement can regrow when both qubits are
coupled to a noisy environment in the limit of large values of decoherence
parameter. The entanglement disappears, irrespective of the acceleration,
under the action of phase damping channel only when the channel is fully
decohered both for single qubit and the two qubits decohering cases. However,
under the action of depolarizing channel an early sudden death occurs for
larger acceleration when both qubits are coupled to the environment. In
summary, the entanglement generated by Unruh effect in noninertial frame is
strongly influenced by decoherence.

## References

  * [1] The Physics of Quantum Information, D. Bouwmeester, A. Ekert, A. Zeilinger (Springer-Verlag, Berlin, 2000)
  * [2] C. H. Bennett, G. Brassard, C. Crepeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * [3] A. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * [4] C.H. Bennett, G. Brassard, N.D. Mermin, Phys. Rev. Lett. 68, 557 (1992).
  * [5] L.K. Grover, Phys. Rev. Lett. 79, 325 (1997).
  * [6] D.P. DiVincenzo, Science 270, 255 (1995)
  * [7] P.M. Alsing, I. Fuentes-Schuller, R. B.Mann, and T. E. Tessier, Phys. Rev. A 74, 032326 (2006).
  * [8] Yi Ling et al, J. Phys. A: Math. Theor. 40, 9025 (2007).
  * [9] R. M. Gingrich and C. Adami Phys. Rev. Lett. 89, 270402 (2002).
  * [10] Q. Pan and J. Jing, Phys. Rev. A 77, 024302 (2008).
  * [11] I. Fuentes-Schuller and R. B. Mann, Phys. Rev. Lett. 95, 120404 (2005).
  * [12] H. Terashima and M. Ueda, Int. J. Quantum Inf. 1, 93 (2003).
  * [13] L. Bombelli, R. K. Koul, J. Lee, and R. Sorkin, Phys. Rev. D 34, 373 (1986).
  * [14] C. Callen and F. Wilzcek, Phys. Lett. B 333, 55 (1994).
  * [15] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975); Phys. Rev. D 14, 2460 (1976).
  * [16] H. Terashima, Phys. Rev. D 61, 104016 (2000).
  * [17] Zurek W H et al. Phys. Today 44, 36 (1991).
  * [18] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford), 2002; H. Carmichael, An Open Systems Approach to Quantum Optics (Springer, Berlin, 1993).
  * [19] W. H. Zurek, Rev. Mod. Phys. 75, 715 (2003).
  * [20] J. Wang and J. Jing, arxiv: 1005.2865v4 (2010).
  * [21] W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * [22] V. Coffman, J. Kundu, and W. K. Wootters, Phys. Rev. A 61, 052306 (2000).