# Assessment of the degree of non-Markovianity of open quantum systems and the
necessary generalization of quantum state distance measures

G. Mouloudakis1,2<EMAIL_ADDRESS>I. Stergou1 P. Lambropoulos1,2
1Department of Physics, University of Crete, P.O. Box 2208, GR-71003
Heraklion, Crete, Greece
2Institute of Electronic Structure and Laser, FORTH, P.O.Box 1527, GR-71110
Heraklion, Greece

###### Abstract

We provide a quantitative evaluation of the degree of non-Markovianity (DNM)
for an XX chain of interacting qubits with one end coupled to a reservoir. The
DNM is assessed in terms of various quantum state distance (QSD) measures and
various cases of non-Markovian spectral densities. Our approach is based on
the construction of the density matrix of the open chain, without the
necessity of a master equation. For the quantification of the DNM we calculate
the dynamics of the QSD measures between the Markovian-damped and various
types of non-Markovian-damped cases. Since in the literature several QSD
measures, appear in forms that imply trace preserving density matrices, we
introduced appropriate modifications so as to render them applicable to open
systems with damped traces. Our results produce remarkable consistency between
the various QSD measures. They also reveal a subtle and potentially useful
interplay between qubit-qubit interaction and non-Markovian damping. Our
calculations have also uncovered a surprisingly dramatic slowing-down of
dissipation by the squared Lorentzian reservoir.

##  I Introduction

The dynamics of open quantum systems coupled to non-Markovian reservoirs is at
the forefront of current interest, as it is relevant to practical as well as
to fundamental problems. In that endeavor, one dimensional many-body systems
such as quantum spin chains arise in many contexts, as they offer valuable and
versatile ground for quantitative studies. Their relevance extends from
quantum information theory to condensed matter physics, owing to their
versatility as basic resources for the implementation of solid-state devices
for quantum computing and quantum communication tasks ref1 . Among these
tasks, faithful quantum state transfer ref2 ; ref3 ; ref4 ; ref5 ; ref6 ; ref7
; ref8 ; ref9  and long-distance entanglement ref10 ; ref11 ; ref12 ; ref13 ;
ref14 ; ref15  have been investigated in detail for a variety of spin chain
configurations, over the last 20 years or so, with research on these fields
continuing at an undiminished pace.

In that general context, the dynamics of an arbitrary length spin XX chain
driven by non-Markovian reservoirs at both ends was studied in a recent paper
by two of us ref16 . Our contribution in that paper was the detailed
quantitative study of reservoirs, with Lorentzian as well as Ohmic spectral
densities. Through a formulation in terms of the wave-function, combined with
Laplace transform, we were able to obtain results for chains of arbitrary
length, i.e. number of sites, and study dynamical as well as state transfer
aspects of the system. Guided by the behaviour of the excitation dynamics, we
assessed the qualitative differences of the system from that driven by a
Markovian reservoir which is known to entail monotonic decay. That assessment,
indirect as it was, can only be considered qualitative, although it did
provide a glimpse of the similarities between these two non-Markovian cases.
Yet, there are a number of technical measures aiming at a quantitative
assessment of the difference between the states of a quantum system, as it
evolves under coupling to different external causes; in this case reservoirs.
The magnitude of the difference between the evolution under a given reservoir
from that under a Markovian one, should therefore provide a quantitative
measure for what, following standard terminology ref17 ; ref18 , we shall
refer to as “degree of non-Markovianity” (DNM). The term degree of
Markovianity (DM) can also be found in the literature, meant to indicate the
similarity to the evolution under a Markovian reservoir. Obviously, a large
DNM corresponds to a small DM. We have adopted the term DNM in this paper, as
it appears to have prevailed in the literature. Since the terms environment
and reservoir appear in the literature as synonyms, we are using them
interchangeably in this paper.

Having examined much of the existing literature on the issue, we have come to
the realization that our scheme provides excellent territory for the
quantitative study of non-Markovianity. The relevant literature is vast,
ranging from general formal considerations ref19 ; ref20 ; ref21 ; ref22 ;
ref23 ; ref24  to examples in small specific systems ref25 ; ref26 ; ref27 ;
ref28 , including methods for the estimation the DNM using machine learning
tools ref29 . The advantage of our system stems from the combination of a
realistic arrangement of an XX chain of qubits, coupled to each other, with
one end coupled to a reservoir. In our previous paper ref16 , we had
considered a chain coupled to reservoirs at both ends. Although the present
quantitative study could as easily address that arrangement, since our
objective is the comparative study of non-Markovian reservoirs, we have chosen
to consider the same system with only one end qubit coupled to a variety of
reservoirs. We calculate first the time evolution of the system under each
reservoir, one of which is Markovian. In the second step, we calculate the
difference of the state of the system evolved under a given non-Markovian
reservoir from the state evolved under the Markovian one. The time evolution
of that difference is what we refer to as the DNM of the given reservoir.

We have chosen to address the issue in terms of perhaps the most direct
notion, namely the quantum state distance (QSD). A variety of measures
appropriate as tools for our task have been proposed over the years ref30 ;
ref31 . Their objective is the quantification of the difference between two
different quantum states of a system, which is of fundamental importance in
quantum information processing ref32 ; ref33 , state transfer being a case in
point ref34 . The most common quantity pertaining to that issue is the
Fidelity (F), for which there are various expressions in the literature ref35
; ref36 ; ref37 . It is meant to characterize the similarity or difference
between the density matrices representing the two states. Obviously it is
related to the notion of distance between two quantum states, which explains
why the formulation of some distance measures involve one of the F
expressions. There is no unique expression or value for the F, which is also
the case for the QSD. In some sense, it may be matter of taste or perhaps
usefulness in a particular context which value of F or QSD is adopted.

Our stated objective then is the quantitative evaluation and calibration of
non-Markovianity of several non-Markovian reservoirs by means of QSD. In view
of the diversity of expressions for the QSD, inevitably an additional
component has been injected into our task, namely, the comparative analysis of
various QSD measures. Specifically, since the value of DNM does depend on the
particular QSD measure, it is important to know whether and to what extend the
classification of various reservoirs in terms of DNM is independent of the QSD
measure employed. And this completes the road map of our present project.

In the next session, we review various QSD measures, their properties and the
modifications necessary in order to accommodate issues pertaining to open
systems. In section III we describe the theoretical approach to our problem in
terms of the time-dependent Schrödinger’s equations of motion in Laplace
space, from which one can obtain even closed form solutions for the
judiciously transformed amplitudes, for chains of arbitrary length, as well as
arbitrary initial conditions, within the single-excitation subspace. Using the
amplitudes of the chain sites, we construct the density matrices necessary for
the calculation of the various QSD measures. In section IV we present the
results of our study while in section V we provide an overview as well as some
concluding remarks related to our work.

##  II Brief summary and necessary modifications of QSD measures

QSD measures need to be formulated in terms of density operators, so as to
allow their applicability to mixed states. In books and papers, more often
than not, expressions for F or QSD are assumed, either tacitly or explicitly,
to involve Hermitian density operators of trace equal to unity. Coupling of an
initially closed system to any reservoir, however, leads to decay. As a
consequence, the reduced density operator of the system becomes non-Hermitian,
while the trace decays usually to zero, as time tends to infinity. This
implies that in expressions involving the density operator, its Hermitian
adjoint, as well as the trace must be retained explicitly, even though
typically at $t=0$ the trace is normalized to unity.

To illustrate the issue, let us consider the most common and straightforward
QSD measure, referred to as the trace distance. It is defined by

$D_{T}\equiv\frac{1}{2}\Tr\absolutevalue{ρ-σ}$ (1)

where $\absolutevalue{A}$, for a matrix $A$, stands for $\sqrt{A^{\dagger}A}$
which stands for the positive definite square root of the matrix under the
radical symbol, with $A^{\dagger}$ being the Hermitian adjoint of $A$.
Applying the above definition to the expression for the trace distance, we
obtain

$D_{T}=\frac{1}{2}\Tr\sqrt{\left(ρ^{\dagger}ρ+σ^{\dagger}σ-ρ^{\dagger}σ-σ^{\dagger}ρ\right)}$
(2)

Now, if both density operators under the radical are Hermitian, the expression
for the trace distance shown in Eqn. (2) reduces to

$D_{T}=\frac{1}{2}\Tr\left(ρ-σ\right)$ (3)

This expression may be valid for closed systems. Incidentally, if the two
density operators are normalized to unity, this expression will always be
zero, which seems like a tautology. However, since for a system coupled to a
reservoir, the density operators become non-Hermitian, the correct expression
for the trace distance is given by Eqn. (2).

Let us consider now another expression for QSD, namely the Hellinger measure
$D_{H}$, defined as

$D_{H}^{2}\equiv\Tr\left(\sqrt{ρ}-\sqrt{σ}\right)^{2}$ (4)

which upon expanding the square becomes

$D_{H}^{2}=\Tr ρ+\Tr σ-2\Tr\left(\sqrt{ρ}\sqrt{σ}\right)$ (5)

where the invariance of the trace under cyclic permutation of the factors has
been used. If the system is closed, with density operators normalized to
unity, Eqn. (5) reduces to
$D_{H}^{2}=2\left[1-\Tr\left(\sqrt{ρ}\sqrt{σ}\right)\right]$, leading to the
expression

$D_{H}=\sqrt{2}\left[1-\Tr\left(\sqrt{ρ}\sqrt{σ}\right)\right]^{1/2}$ (6)

for the Hellinger distance measure, which typically is the expression cited
for this QSD measure. However, in the general case of open systems, following
Eqn. (5), we have:

$D_{H}=\left[\Tr ρ+\Tr σ-2\Tr\left(\sqrt{ρ}\sqrt{σ}\right)\right]^{1/2}$ (7)

When we are dealing with open systems, which implies coupling to reservoirs,
the density operators tend to zero as $t\rightarrow\infty$. In that limit,
Eqn. (7) correctly yields zero for the distance between the two states, as
they both decay. On the other hand, Eqn. (6) would lead to the value
$\sqrt{2}$ which on physical grounds, is at best problematic, as it does not
account for the decay of the traces of the density operators.

There is one more measure that we employ in our calculations. It is known as
the Bures distance measure, usually denoted by $D_{B}$. Since it is typically
defined in terms of one of the expressions for F, we list below the three most
common expressions, labelled for our convenience $F_{1}$, $F_{2}$ and $F_{3}$.
They are:

$F_{1}\left(ρ,σ\right)=\left(\Tr\sqrt{\sqrt{ρ}σ\sqrt{ρ}}\right)^{2}$ (8a)

$F_{2}\left(ρ,σ\right)=\Tr\sqrt{\sqrt{ρ}σ\sqrt{ρ}}=\sqrt{F_{1}\left(ρ,σ\right)}$
(8b)

$F_{3}\left(ρ,σ\right)=\Tr\left(ρσ\right)$ (8c)

The Bures QSD is usually defined as $D_{B}^{2}\left(ρ,σ\right)\equiv
2\left[1-\sqrt{F_{1}\left(ρ,σ\right)}\right]$, which, in view of Eqn. (8b)
above, can be simply written as

$D_{B}^{2}\left(ρ,σ\right)=2\left[1-F_{2}\left(ρ,σ\right)\right]$ (9)

from which we obtain

$D_{B}\left(ρ,σ\right)=\sqrt{2}\left[1-\Tr\sqrt{\sqrt{ρ}σ\sqrt{ρ}}\right]^{1/2}$
(10)

Again, this expression is valid as long as the traces of both density
operators remain equal to one, because if in any distance measure we set
$ρ=σ$, we should obtain zero. This condition is satisfied as long as $\Tr
ρ=\Tr σ=1$ for all times. In the presence of dissipation, however, the traces
decay to zero, in which case for $ρ=σ$ we would obtain the nonphysical value
$\sqrt{2}$. This contradiction is amended if the above expression is modified
as shown in the following equation:

$D_{B}\left(ρ,σ\right)=\sqrt{2}\left[\frac{1}{2}\left(\Tr ρ+\Tr
σ\right)-\Tr\sqrt{\sqrt{ρ}σ\sqrt{ρ}}\right]^{1/2}$ (11)

This is the expression we are using for the Bures QSD in this paper. In the
long time limit and in the presence of dissipation, the modified Bures
distance tends to zero as it should. Part of our investigation addresses the
rate with which the various QSD measures tend to their final value. We are in
addition interested in the consistency of the values of the DNM obtained by
different measures. For whatever it is worth, it seems to us that one could as
well define a QSD measure through an extension of Eqn. (9), by inserting any
of the expressions for F.

##  III Theory

Figure 1: Schematic representation of the system at study and our methodology.
A chain of an arbitrary number of coupled qubits is interacting with an
external reservoir at its one end. The system is studied for two different
cases of external reservoirs, namely for a Markovian and a non-Markovian
reservoir with spectral density $J(ω)$. Through our methodology we calculate
the qubit amplitudes in the single-excitation subspace via which we construct
the density matrix of the chain for each case. The QSD between the two density
matrices is then calculated using various measures, allowing us to assess the
degree of non-Markovianity of the system for different types of spectral
densities $J(ω)$.

As discussed in section II, all expressions for QSD measures involve the
density matrices $ρ$ and $σ$ of the quantum states whose distance is to be
calculated. In general, the temporal evolution of the density matrix of a
system coupled to an external environment is given by the solution of the
master equation, obtained after tracing out the environmental degrees of
freedom from the density matrix of the compound system (system + environment)
ref38 . Despite the ongoing progress in the derivation of closed form
expressions for master equations of quantum systems coupled to non-Markovian
reservoirs ref39 ; ref40 ; ref41 ; ref42 ; ref43 , such expressions are
usually obtained for special cases of environmental spectral densities under
certain approximations and often are too complicated to be handled even
numerically. Consequently, an attempt to quantify the degree of non-
Markovianity in systems coupled to environments with a variety of possible
spectral densities using numerical solutions of the respective master
equations would be pointless -if possible at all- for certain forms of
environmental spectral densities, such as those studied in this work.

To address this difficulty, we have developed a formulation enabling the
calculation of the amplitudes of the wavefunction of the open system, from
which we then construct its density matrix. This formulation has been
described in detail in a previous paper ref16 , where it was used to calculate
aspects of the dynamics of a XX chain boundary-driven by non-Markovian
environments. As demonstrated below, one of the advantages of this approach is
its applicability to any possible form of environmental spectral density, as
well as an arbitrary number of qubits in the chain.

The system under consideration consists of $N$ identical qubits in a one-
dimensional configuration coupled to each other via a nearest-neighbor
coupling $\mathcal{J}$. The $N^{th}$ qubit of the chain is coupled, via a
coupling strength $g$, to an external environment of an arbitrary spectral
density $J(ω)$. Our goal is to assess the degree of non-Markovianity of that
system by means of the QSD between its density matrix and the density matrix
of the system consisting again by $N$ qubits but with the $N^{th}$ qubit
coupled to a Markovian environment. The Markovian environment is known to
induce a shift of the energy of the last qubit, as well as a decay ref38 . We
neglect the shift, as it is of no significance in the issue of dissipation and
denote the decay by $\gamma_{M}$. The problem is studied using several QSD
measures, for various forms of non-Markovian spectral densities $J(ω)$.
Hereafter, we denote the density matrix of the Markovian open system by $σ(t)$
and the density matrix of the non-Markovian open system by $ρ(t)$. A schematic
presentation of the systems under study and our methodology is depicted in
Fig. 1.

The Hamiltonian of the compound system
$\mathcal{H}=\mathcal{H}_{S}+\mathcal{H}_{E}+\mathcal{H}_{I}$, consists of
three parts; namely, the Hamiltonian of the XX chain $\mathcal{H}_{S}$, the
Hamiltonian of the environment which is coupled to $N^{th}$ (last) qubit of
the chain $\mathcal{H}_{E}$, and the chain-environment interaction Hamiltonian
$\mathcal{H}_{I}$, given by the relations $(\hbar=1)$:

$\begin{split}\mathcal{H}_{S}=&\omega_{e}\sum_{i=1}^{N}\ket{e}_{i}\prescript{}{i}{\bra{e}}+\omega_{g}\sum_{i=1}^{N}\ket{g}_{i}\prescript{}{i}{\bra{g}}\\\
&+\sum_{i=1}^{N-1}\frac{\mathcal{J}}{2}\left(\sigma_{i}^{+}\sigma_{i+1}^{-}+\sigma_{i}^{-}\sigma_{i+1}^{+}\right),\end{split}$
(12a)

$\mathcal{H}_{E}=\sum_{\lambda}\omega_{\lambda}{a^{E}_{\lambda}}^{\dagger}a_{\lambda}^{E},$
(12b)

$\mathcal{H}_{I}=\sum_{\lambda}g(\omega_{\lambda})\left(a_{\lambda}^{E}\sigma_{N}^{+}+{a_{\lambda}^{E}}^{\dagger}\sigma_{N}^{-}\right),$
(12c)

where $\omega_{g}$ and $\omega_{e}$ are, respectively, the energies of the
ground and excited state of each spin (all spins are assumed identical),
$\omega_{\lambda}$ is the energy of the $\lambda$-mode of the environment
consisting of an infinite number of bosons, $a_{\lambda}^{E}$ and
${a^{E}_{\lambda}}^{\dagger}$, are the bosonic annihilation and creation
operators of the environment, and
$\sigma_{i}^{+}=\ket{e}_{i}\prescript{}{i}{\bra{g}}$ and
$\sigma_{i}^{-}=\ket{g}_{i}\prescript{}{i}{\bra{e}}$, $i=1,\dots,N$ are the
qubit raising and lowering operators, respectively.

The wavefunction of the compound system in the single-excitation space is
expressed as:

$\ket{\Psi(t)}=\sum_{i=1}^{N}c_{i}(t)\ket{\psi_{i}}+\sum_{\lambda}c_{\lambda}^{E}(t)\ket{\psi_{\lambda}^{E}},$
(13)

where,

$\ket{\psi_{i}}\equiv\ket{g}_{1}\ket{g}_{2}\dots\ket{g}_{i-1}\ket{e}_{i}\ket{g}_{i+1}\dots\ket{g}_{N}\ket{0}_{E},$
(14a) and
$\ket{\psi_{\lambda}^{E}}\equiv\ket{g}_{1}\ket{g}_{2}\dots\ket{g}_{N}\ket{00\dots
01_{\lambda}0\dots 00}_{E}.$ (14b)

###  III.1 Construction of $σ(t)$ (Markovian system)

As already noted above, a Markovian environment (i.e. an environment with
slowly varying spectral density) interacting with a two-level quantum system
(qubit) within the Born approximation, induces a shift as well as a decay in
the energy of the latter. Formally, this result is obtained by tracing the
environmental degrees of freedom out of the density matrix of the compound
system, a procedure found in many quantum optics and open quantum systems
textbooks ref38 ; ref44 . In our problem, the Markovian environment will only
affect the energy of the last qubit ($N^{th}$ qubit) of the chain. Neglecting
the shift and keeping the decay rate $\gamma_{M}$, the time dependent
Schrödinger equation (TDSE) for the chain amplitudes reads:

$\begin{split}i\frac{dc_{i}(t)}{dt}=&\left[\omega_{e}+(N-1)\omega_{g}\right]c_{i}(t)+\frac{\mathcal{J}}{2}\left[c_{i-1}(t)+c_{i+1}(t)\right],\\\
&\quad i=1,\dots,N-1,\end{split}$ (15a)
$i\frac{dc_{N}(t)}{dt}=\left[\omega_{e}+(N-1)\omega_{g}-iγ_{M}\right]c_{N}(t)+\frac{\mathcal{J}}{2}c_{N-1}(t)$
(15b)

where, in order to compress the notation in our equations, we allowed $i$ to
run from 1 to $N-1$ in Eqn. (15a), noting that $c_{i-1}(t)=0$ for $i=1$,
because qubit $1$ is coupled only to qubit $2$. To simplify the set of Eqns.
(15b), it is useful to introduce the amplitude transformations
$c_{i}(t)=e^{-i\left[\omega_{e}+(N-1)\omega_{g}\right]t}\tilde{c}_{i}(t)$,
$i=1,\dots,N-1$ and
$c_{N}(t)=e^{-i\left[\omega_{e}+(N-1)\omega_{g}-i\gamma_{M}\right]t}\tilde{c}_{N}(t)$,
in terms of which our equations become:

$\frac{d\tilde{c}_{i}(t)}{dt}=-i\frac{\mathcal{J}}{2}\left[\tilde{c}_{i-1}(t)+\tilde{c}_{i+1}(t)\right],\quad
i=1,\dots,N-1,$ (16a)
$\frac{d\tilde{c}_{N}(t)}{dt}=-i\frac{\mathcal{J}}{2}e^{\gamma_{M}t}\tilde{c}_{N-1}(t)$
(16b)

where, as above, $\tilde{c}_{i-1}(t)=0$ for $i=1$. The above set of $N$
coupled linear differential equations can be readily solved numerically to
yield the time dependence of the tilde amplitudes of the chain sites. The
density matrix $\sigma(t)$ of the chain can then be constructed in terms of
the site amplitudes as shown the following equation:

$\sigma(t)=\sum_{i,j}^{N}c_{i}(t)c_{j}^{*}(t)\ket{\psi_{i}}\bra{\psi_{j}}$
(17)

###  III.2 Construction of $ρ(t)$ (non-Markovian system)

The construction of $ρ(t)$, i.e. the density matrix of the chain with the last
qubit interacting with a non-Markovian reservoir, is a subtler task because
the memory effects of the reservoir make the equation of the last qubit
amplitude more complicated. In two of our recent papers, we examined the
dynamics of XX chain boundary driven by non-Markovian reservoirs ref16 , as
well as the connection between exceptional points and the Quantum Zeno effect
in a simpler quantum system consisting of two interacting qubits one of which
is coupled to a non-Markovian reservoir ref45 . On the basis of the
formulation developed in those papers, as applied to the present system, we
outline here the basic procedure to obtain the amplitudes of the chain sites.
In order to discriminate the amplitudes from those of the Markovian damped
chain, we denote them by $c^{\prime}_{i}(t)$, $i=1,\dots N$.

Adopting the amplitude transformations
$c_{i}^{\prime}(t)=e^{-i\left[\omega_{e}+(N-1)\omega_{g}\right]t}{\tilde{c}^{\prime}}_{i}(t)$
, $i=1,\dots N$ and
${c^{\prime}_{\lambda}}^{E}(t)=e^{-i\left(N\omega_{g}+\omega_{\lambda}\right)t}{{\tilde{c}}_{\lambda}}^{{}^{\prime}E}(t)$,
the TDSE yields the following equations of motion:

$i\frac{d\tilde{c}^{\prime}_{i}(t)}{dt}=\frac{\mathcal{J}}{2}\left[\tilde{c}^{\prime}_{i-1}(t)+\tilde{c}^{\prime}_{i+1}(t)\right],\quad
i=1,\dots,N-1,$ (18a)

$i\frac{d\tilde{c}^{\prime}_{N}(t)}{dt}=\frac{\mathcal{J}}{2}\tilde{c}^{\prime}_{N-1}(t)+\sum_{\lambda}g(\omega_{\lambda})e^{-i\Delta_{\lambda}t}\tilde{c}_{\lambda}^{{}^{\prime}E}(t),$
(18b)

$i\frac{d\tilde{c}_{\lambda}^{{}^{\prime}E}(t)}{dt}=g(\omega_{\lambda})e^{i\Delta_{\lambda}t}\tilde{c}^{\prime}_{N}(t),$
(18c)

where
$\Delta_{\lambda}\equiv\omega_{\lambda}-(\omega_{e}-\omega_{g})\equiv\omega_{\lambda}-\omega_{eg}$
and $\tilde{c}^{\prime}_{i-1}(t)=0$ for $i=1$.

Formal integration of Eqn. (18c), under the initial condition
$\tilde{c}_{\lambda}^{{}^{\prime}E}(0)=0$, and substitution into Eqn. (18b),
yields:

$\begin{split}\frac{d\tilde{c}^{\prime}_{N}(t)}{dt}=&-i\frac{\mathcal{J}}{2}\tilde{c}^{\prime}_{N-1}(t)\\\
&-\int_{0}^{t}\sum_{\lambda}e^{-i\Delta_{\lambda}\left(t-t^{\prime}\right)}\left[g(\omega_{\lambda})\right]^{2}\tilde{c}^{\prime}_{N}(t^{\prime})dt^{\prime}.\end{split}$
(19)

The summation over the environmental modes is at this point replaced by an
integration which requires the specification of the environment’s spectral
density $J(\omega)$, according to the relation
$\sum_{\lambda}\left[g(\omega_{\lambda})\right]^{2}\rightarrow\int d\omega
J(\omega)$. The resulting set of differential equations then becomes:

$\frac{d\tilde{c}^{\prime}_{i}(t)}{dt}=-i\frac{\mathcal{J}}{2}\left[\tilde{c}^{\prime}_{i-1}(t)+\tilde{c}^{\prime}_{i+1}(t)\right],\quad
i=1,\dots,N-1,$ (20a)

$\frac{d\tilde{c}^{\prime}_{N}(t)}{dt}=-i\frac{\mathcal{J}}{2}\tilde{c}^{\prime}_{N-1}(t)-\int_{0}^{t}R(t-t^{\prime})\tilde{c}^{\prime}_{N}(t^{\prime})dt^{\prime},$
(20b)

where

$R(t)\equiv\int_{0}^{\infty}J(\omega)e^{-i\left(\omega-\omega_{eg}\right)t}d\omega.$
(21)

The second term in Eqn. (20b) reflects the possibility for the excitation to
be transferred from the $N^{th}$ qubit to the non-Markovian reservoir and vice
versa.

Taking the Laplace transform of Eqns. (20b) and following the procedure
presented in ref16  we find that the Laplace transform $F_{1}(s)$ of the tilde
amplitude $\tilde{c}^{\prime}_{1}(t)$ is given by the expression:

$F_{1}(s)=\frac{ikc_{N}(0)-k^{2}\left[s+B(s)\right]A_{1}(s)c_{N-1}(0)+(ik)\sum_{m=1}^{N-2}\Big{\\{}ik\left[s+B(s)\right]A_{N-m}(s)-A_{N-1-m}(s)\Big{\\}}c_{m}(0)}{ik\left[s+B(s)\right]A_{N}(s)-\left[1+s+B(s)\right]A_{N-1}(s)}.$
(22)

where $k\equiv\frac{2}{\mathcal{J}}$, $B(s)$ is the Laplace transform of the
function $R(t)$ and $A_{m}(s)$ is given by the following expression:

$A_{m}(s)\equiv\frac{\left[(iks)+i\sqrt{k^{2}s^{2}+4}\right]^{m}-\left[(iks)-i\sqrt{k^{2}s^{2}+4}\right]^{m}}{2^{m}i\sqrt{k^{2}s^{2}+4}},\quad
m=1,\dots,N.$ (23)

The Laplace transforms $F_{i}(s)$ of the remaining qubit amplitudes are
connected to $F_{1}(s)$ via the relation:

$\begin{split}F_{i}(s)&=A_{i}(s)F_{1}(s)-(ik)\sum_{n=1}^{i-1}A_{i-n}(s)c_{n}(0),\\\
&\quad i=2,\dots,N,\end{split}$ (24)

The effectiveness of this method rests on the fact that it allows the
derivation of closed form expressions for the Laplace transforms of the tilde
amplitudes for an arbitrary number of qubits $N$, an arbitrary spectral
density $J(\omega)$ for the environment that interacts with the $N^{th}$ qubit
of the chain, as well as arbitrary initial conditions. The Laplace inversion
integrals can be readily calculated numerically (or even analytically in some
cases) upon specification of these parameters to yield the tilde amplitudes of
the chain sites, in terms of which we can calculate $c_{i}^{\prime}(t)$ ,
$i=1,\dots N$, and thus construct the density matrix for the chain via the
relation:

$\rho(t)=\sum_{i,j}^{N}c^{\prime}_{i}(t){c_{j}^{\prime}}^{*}(t)\ket{\psi_{i}}\bra{\psi_{j}}$
(25)

Specific types of non-Markovian reservoirs, such as Lorentzian ref46 ; ref47 ,
Lorentzian squared ref48  and Ohmic ref49 ; ref50  are employed in order to
assess the non-Markovianity of the chain in each case. The analytical
derivations of the functions $R(t)$ and $B(s)$ for these types of reservoirs
are given in Appendix A.

##  IV Results $\&$ Discussion

Figure 2: Dynamics of the excitation survival probability of a single qubit
coupled to various types of reservoirs: (a) Markovian reservoir with
$γ_{M}=0.01$, (b) Lorentzian reservoir with $g=1$, $γ=0.03$ and $Δ_{c}=0$, (c)
Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_{c}=0$, (d) Ohmic
reservoir with $g=1$, $S=1.5$, $ω_{c}=8$ and qubit frequency $ω_{eg}=10$.
Figure 3: QSD measures between the Markovian-damped and the non-Markovian-
damped system using various types of non-Markovian reservoirs in the single-
qubit case. The types of reservoirs chosen and their parameters are: Markovian
reservoir with $γ_{M}=0.01$, Lorentzian reservoir with $g=1$, $γ=0.03$ and
$Δ_{c}=0$, Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_{c}=0$,
Ohmic reservoir with $g=1$, $S=1.5$, $ω_{c}=8$ and qubit frequency
$ω_{eg}=10$. (a) Markovian-Lorentzian Trace Distance, (b) Markovian-
Lorentzian-squared Trace Distance, (c) Markovian-Ohmic Trace Distance, (d)
Markovian-Lorentzian Hellinger Distance, (e) Markovian-Lorentzian-squared
Hellinger Distance, (f) Markovian-Ohmic Hellinger Distance, (g) Markovian-
Lorentzian Bures Distance, (h) Markovian-Lorentzian-squared Bures Distance,
(i) Markovian-Ohmic Bures Distance.

The stated objective of this paper is the quantitative assessment of the DNM
of non-Markovian reservoirs in a realistic context provided by a spin chain.
As outlined in the previous sections, for the task we have chosen to employ
measures of QSD between the state of the chain coupled to various reservoirs.
An obviously suitable for the purpose quantity is the QSD between the density
operator of the chain evolved under a Markovian from that under a non-
Markovian reservoir. From a physical standpoint, it stands to reason that the
larger the values of the QSD, irrespective of the measure under consideration,
the less Markovian the system’s evolution will be, which means higher DNM.
Given that reservoirs by definition involve certain parameters, characterising
their properties, in addition to the system-reservoir coupling constant, a
decision has to be made as to the choice of parameters that make the value of
QSD most meaningful. It is rather obvious that a direct comparison between,
for example, an evolution under a Markovian reservoir with damping rate
$\gamma_{M}$ and a Lorentzian one with arbitrarily chosen parameters $g$,
$\gamma$ and $\Delta_{c}$ would be pointless. We suggest that, for the
comparison to be meaningful, a feature common to the effect of all reservoirs
entering the comparison should be adopted. With that in mind, we note that for
all reservoirs the excitation of all qubits will eventually go to zero. Since
initially the excitation is in the first qubit, without claim to uniqueness,
we have chosen the parameters of the various reservoirs such that the half
life of the excitation of the first qubit is approximately the same to its
half life under a Markovian reservoir with a given $\gamma_{M}$.

###  IV.1 Single qubit

We begin our analysis by considering the simplest case of a single qubit
coupled to a variety of non-Markovian reservoirs. These results serve as a
point of reference in the comparative analysis of the interplay between non-
Markovianity and qubit-qubit interaction, as reflected in the behaviour of the
QSD measures as a function of the time, analysed in section IVB.

In Fig. 2 we plot the time dynamics of the single qubit excitation survival
probability for various types of reservoirs. In particular, we examine the
cases of Markovian, Lorentzian, Lorentzian squared and Ohmic reservoirs, with
the parameters of each non-Markovian reservoir chosen as described above. The
Lorentzian and Lorentzian squared spectral densities result in general to fast
oscillations in the single-qubit survival of excitation dynamics, indicating
the exchange of the excitation between the qubit and the environment within
finite times. The Ohmic spectral distribution tends in general to be much
broader than the Lorentzian and Lorentzian squared distributions, leading to a
qubit excitation survival probability that resembles more the respective
survival probability of the excitation of the Markovian-damped qubit.

The QSD measures between the Markovian-damped single-qubit system and the non-
Markovian ones are presented in Fig. 3. In all three of the considered non-
Markovian reservoirs, the coupling strength $g$ between the single-qubit and
the reservoir has been kept the same. As expected, all QSD measures are zero
for $t=0$, as well as in the long time limit, as there is practically no
excitation left in the open system.

We focus first on the comparison of the degrees of non-Markovianity, for
different types of non-Markovian reservoirs, resulting from the same QSD
measure. The trace distance between the Markovian and the Lorentzian damped
system (Fig. 3(a)) exhibits rapid oscillations between two bounds, whose
values change over time. The upper bound is initially at its maximum value and
over time decreases to zero, whereas the lower bound initially increases,
reaches a maximum, after which it decreases merging with the upper bound,
tending eventually to zero for long times. The dynamics of the trace distance
between the Markovian and the Lorentzian-squared damped system (Fig. 3(b))
also exhibit rapid oscillations, following a trend similar to the Markovian-
Lorentzian case, with the exception that the lower bound can be zero within
finite times, indicating that, at those instants, the Lorentzian-squared
damped system resembles the Markovian-damped one. The oscillations in both
Figs. 3(a) and 3(b) are indicative of the fast excitation exchange between the
qubit and the non-Markovian reservoirs. In Fig. 3(c) we plot the trace
distance between the Markovian-damped and the Ohmic-damped qubit state as a
function of time. Clearly, the behaviour of the trace distance dynamics in
this case differs substantially from the dynamics of the Markovian-Lorentzian
and the Markovian-Lorentzian-squared cases, in that it does not exhibit any
oscillations, reaches a maximum after some finite time before decreasing, with
the overall values of the trace distance remaining much smaller than the
respective values in panels (a) and (b) of Fig. 3. We may therefore safely
conclude that, under the prescribed conditions on parameters, the Ohmic-damped
single qubit system displays the lowest DNM.

In addition to the trace distance, we have calculated QSD’s in terms of two
other relevant measures, namely the Hellinger and Bures measures, with the
results shown in Fig. 3. The overall behavior bears significant similarity to
that of the trace distance, with some non-negligible quantitative differences.
The oscillations reflecting the exchange of excitation between qubit and
reservoir are still present in both measures for the Markovian-Lorentzian and
Markovian-Lorentzian-squared cases, whereas they are absent for the Markovian-
Ohmic one. However both of these measures produce generally slower decrease of
the QSD from the respective trace distance measure. For example, in terms of
the trace distance, the value of the QSD between Lorentzian and Markovian in
the vicinity of $t=100$ hovers around $0.15$, whereas the Hellinger and Bures
measures indicate values around a mean of $0.4$. Moreover, the oscillations of
the QSD for the Markovian obtained from these two measures are damped much
more slowly than the oscillations in the trace distance. The Ohmic case seems
to follow the same trend, in that these two measures produce larger values for
the QSD and a much slower evolution towards the expected value of $0$. In
summary, for all three non-Markovian reservoirs in the single qubit case, the
trace distance is found to show faster decay of non-Markovianity than the
Hellinger and Bures measures do. It is worth noting that, despite the non-
trivial formal difference of the expressions for the Hellinger and Bures QSD
(Eqns. (7) and (11)), they lead to practically the same dynamics, indicating
the same DNM.

###  IV.2 Five mutually interacting qubits

Figure 4: Dynamics of the excitation survival probability of the first qubit
of a chain consisting of $N=5$ qubits with qubit-qubit coupling
$\mathcal{J}=1$. The initial excitation is chosen to be on the first qubit
while the last qubit of the chain is coupled to various types of reservoirs:
(a) Markovian reservoir with $γ_{M}=0.01$, (b) Lorentzian reservoir with
$g=1$, $γ=0.03$ and $Δ_{c}=0$, (c) Lorentzian squared reservoir with $g=1$,
$γ=0.3$ and $Δ_{c}=0$, (d) Ohmic reservoir with $g=1$, $S=1.5$, $ω_{c}=8$ and
qubit frequency $ω_{eg}=10$. Figure 5: QSD measures between the Markovian-
damped and the non-Markovian-damped system using various types of non-
Markovian reservoirs for the case of $N=5$ qubits with qubit-qubit coupling
$\mathcal{J}=1$. The types of reservoirs chosen and their parameters are:
Markovian reservoir with $γ_{M}=0.01$, Lorentzian reservoir with $g=1$,
$γ=0.03$ and $Δ_{c}=0$, Lorentzian squared reservoir with $g=1$, $γ=0.3$ and
$Δ_{c}=0$, Ohmic reservoir with $g=1$, $S=1.5$, $ω_{c}=8$ and qubit frequency
$ω_{eg}=10$. (a) Markovian-Lorentzian Trace Distance, (b) Markovian-
Lorentzian-squared Trace Distance, (c) Markovian-Ohmic Trace Distance, (d)
Markovian-Lorentzian Hellinger Distance, (e) Markovian-Lorentzian-squared
Hellinger Distance, (f) Markovian-Ohmic Hellinger Distance, (g) Markovian-
Lorentzian Bures Distance, (h) Markovian-Lorentzian-squared Bures Distance,
(i) Markovian-Ohmic Bures Distance. In the inset of panel (f) we show the
dynamics of the probability of finding the excitation in the Ohmic environment
for the case of $N=5$ qubits. The time window where the probability dynamics
abruptly become oscillatory coincides with the time window where the
Markovian-Ohmic QSD measures follow the same behaviour.

The merit of our approach rests with the possibility to obtain
straightforwardly the Laplace transforms of the qubits tilde state amplitudes
for an arbitrary number of qubits $N$, because $N$ appears in our equations as
a modifiable parameter. Having calibrated the measures of QSD in
characterising non-Markovianity in the simplest set-up of a single qubit, in
this section we explore the role of non-Markovianity in the more realistic
situation of a chain involving qubit-qubit interaction. For a concrete
quantitative analysis, we have chosen five qubits in an XX chain with nearest
neighbor coupling $\mathcal{J}$. The calculation follows the formulation
presented in section III, with the initial excitation in qubit $1$ and the
last qubit $5$ coupled to the reservoir. The question now is how the non-
Markovianity of the reservoir affects the exchange of excitation between the
qubits, as manifested in the evolution of QSD between the Markovian and the
various non-Markovian reservoirs.

In a set-up of this type, the crucial parameter is $\mathcal{J}$ as it
controls the communication between qubits, which typically is the dominant
aspect in applications. In order to remain within a realistic scenario, we
have chosen the Markovian damping rate $\gamma_{M}$ much smaller than
$\mathcal{J}$, keeping the choice of the parameters of the non-Markovian
reservoirs as in the previous section, so that the initial excitation of the
first qubit decays on a time scale approximately the same for all reservoirs.
At first, we calculate the time evolution of the excitation probability of the
first qubit for all of the reservoirs under consideration, with the results
shown in Fig. 4. The glaring difference between the single-qubit dynamics in
Fig. 2 and the $N=5$ dynamics of the first qubit excitation in Fig. 4, is the
appearance of oscillations in the dynamics of the latter, indicating the
exchange of population between the qubits of the chain. In the cases of
Lorentzian and Lorentzian squared reservoirs, such oscillations are
superimposed on the oscillations between the last qubit and the reservoir,
resulting to the characteristic oscillatory dynamics of $P_{1}(t)$ depicted in
Figs. 4(b) and 4(c). For the Markovian case (Fig. 4(a)) the oscillations are
due solely to the coupling between the qubits, because as we have seen in the
case of the single qubit (Fig. 2(a)) the Markovian damping is monotonic.
Again, for the Ohmic-damped case, the dynamics of the first qubit population
(Fig. 4(d)) bears strong similarity to the Markovian-damped chain, with a
rather unexpected exception, in that the oscillations are found to collapse
abruptly after a finite time, a behaviour that has been also reported in
previous works involving qubits coupled to Ohmic reservoirs ref16 ; ref51 .

We consider now the dynamics of the QSD measures exploring the non-
Markovianity of the three reservoirs in the context of a chain of $N=5$
interacting qubits, with the results shown in Fig. 5. First we note that the
time scales of decay for all QSD measures are longer than the respective time
scales in the single-qubit case of Fig. 3. This behaviour was to be expected,
since the excitation is now spread over the whole chain and naturally it takes
more time to be lost in the reservoir. At the same time, the values of the QSD
measures are in general larger than the respective values of Fig. 3,
indicating that the DNM of the open system is affected by the number of qubits
in the chain. The Markovian-Lorentzian and Markovian-Lorentzian-squared QSD
measures are seen to display an overall similar behaviour. We should, however,
keep in mind that the values of the parameters for the Lorentzian and
Lorentzian-squared spectral densities have been chosen according to the
criterion explained in the discussion of Fig. 1. In contrast, an open system
damped by a Lorentzian squared environment would display a much larger DNM
compared to the same system damped by a Lorentzian spectral environment
characterized by the same parameters ($g$, $γ$ and $Δ_{c}$). In fact, under
those conditions, the excitation throughout the chain in the Lorentzian
squared case remains equal to one for a very long time such that it resembles
a practically steady state. All QSD measures we examined (Trace distance,
Hellinger distance and Bures distance) exhibit dynamics following the same
trend, although the exact distance values may differ in general between
different measures. As clearly seen in Fig. 5, the Hellinger and Bures
distances display minor differences in their dynamical behaviour for all of
the reservoirs studied.

A rather interesting phenomenon appears in the dynamical behaviour of the QSD
measures in the case of the Ohmic-damped chain. In contrast to the single-
qubit case, the interplay between qubit-qubit and qubit- Ohmic reservoir
interactions gives rise to an intricate behavior in the dynamics for all of
the Markovian-Ohmic QSD measures considered. The dynamics are initially smooth
for all QSD measures with the distance being zero at $t=0$ and slowly
increasing afterwards. However, there is a time when they suddenly change and
exhibit rapid oscillations. These rapid oscillations appear only in a time
window which is the same for all QSD measures considered. Beyond that time
window, the dynamics are smoothed out with the distance slowly decaying back
to zero. The time window of that sudden modification of the behaviour of the
QSD measures coincides with the time window within which the total probability
of finding the excitation in the Ohmic reservoir, exhibits similar behaviour
(see inset of Fig. 5(f)). Changes in the dynamical behaviour of the system can
be also observed through abrupt changes in the dynamics of the qubit
excitations of the chain, some of which were reported in one of our previous
papers ref16 . In summary, the qubit-qubit interaction is found to enhance the
DNM obtained in all QSD measures. Conversely, this may be interpreted as an
enhancement of the excitation storage in the chain due to the non-Markovian
character of the reservoir.

##  V Summary and Outlook

The work in this paper was undertaken with the dual purpose of evaluating the
DNM of some typical non-Markovian reservoirs, as well as the generalization of
the most common QSD measures employed for the purpose. In order to offer an
analysis connected with a realistic system typical in quantum information
processing, we chose an XX chain of interacting qubits with a reservoir
connected to one end of the chain. QSD measures require the density matrices
of the systems under consideration. Since non-Markovian reservoirs do not in
general lend themselves to descriptions in terms of a master equation, we
employed an approach developed in earlier work ref16  based on the solution of
the Schrödinger equation, in which the amplitudes of the excitation of the
sites as a function of time are obtained by means of Laplace transform. With
the amplitudes in hand we construct the density matrix of the system which is
necessary for the calculation of the desired QSD measure.

The expressions for QSD measures found in most of the literature explicitly or
implicitly assume density matrices of trace unity. However the density matrix
of an open system by definition will eventually decay to zero, which means the
trace cannot remain equal to unity. That led us to a generalization of the
expressions of certain QSD measures so that they are valid for open as well as
closed systems.

For the DNM we have chosen the distance of the quantum state of our system as
evolved under a given non-Markovian reservoir from the state evolved under a
Markovian one. A meaningful comparison required the judicious choice of the
parameters entering the expressions of the spectral densities of the
reservoirs involved in each QSD calculation. Given the flexibility of our
approach as to the number of qubits it can handle, we studied first the case
of a single-qubit coupled to a reservoir, so as to have a frame for the
evaluation of the interplay between qubit-qubit interaction and qubit-
reservoir. The results from the single qubit analysis showed that the
Hellinger and Bures measures indicated slightly higher DNM for all three
reservoirs than the trace distance measure did. On the other hand, all three
measures produced by far much lower DNM for the Ohmic reservoir compared to
the Lorentzian or Lorentzian squared. Moreover, the absence of any noticeable
back-flow of excitation from the Ohmic reservoir to the qubit suggests that
the Ohmic reservoir behaves practically in a Markovian fashion.

The Lorentzian and Lorentzian squared reservoirs, on the basis of all three
measures were found to exhibit essentially the same DNM. This might seem
counter intuitive, because inspection of the relevant formal expressions in
the appendix suggests that the squared Lorentzian is more peaked than the
Lorentzian with the same parameters, which would imply higher DNM. However,
the Lorentzian squared employed in all of our comparative calculations (Figs.
2-5) was not the square of the accompanying Lorentzian, owing to the
adjustment of the parameters explained in section III. If we calculated
evolution of the qubit excitation dynamics and DNM for a Lorentzian and a
Lorentzian squared spectral density characterized by the same parameters ($g$,
$γ$ and $Δ_{c}$), we would indeed find more back-flow and higher DNM for the
Lorentzian squared, which brings up an interesting issue in reservoir
engineering. If in any open system subjected to dissipation through a
Lorentzian reservoir, it were possible to modify the reservoir to the square
of the Lorentzian, dissipation would slow down significantly. Our
consideration of the Lorentzian squared spectral density is not a mere
mathematical exercise. It was inspired by the observation of non-Lorentzian
line shapes ref48 ; ref52 , including the Lorentzian squared, albeit in
different physical context.

The results for the case of five interacting qubits chain bear similarities to
as well as differences from the single-qubit case. The similarity is seen in
the DNM by the three QSD measures which produce more or less similar values,
with the value of the trace distance somewhat lower than the values of the
Hellinger and Bures measures. In other words, the three measures provide
mutually compatible results. However, the values of the DNM for all three
reservoirs for the chain are noticeably higher than those for the single-
qubit. Moreover, for the Ohmic reservoir, the QSD obtained by all three
measures is not only quantitatively but even qualitatively different from that
of the single qubit. The smooth dynamics obtained in the single qubit case, is
in this case interrupted by a region of abrupt oscillations, within a certain
time window in which the dynamics of the probability of finding the excitation
in the Ohmic environment was found to exhibit the same behaviour. However, the
stunning difference between the single and five qubit cases lies in the
retention of population distributed among the interacting qubits and the much
slower decay of the excitation in all of the reservoirs in the latter case,
which is compatible with the consistently higher DNMs reported in our work. It
can then be argued that the exchange of excitation between the qubits slows
down the dissipation which implies that the larger the number of the qubit the
slower the loss into the reservoir. Perhaps the most significant message
conveyed by the comparison of the single-qubit and chain of interacting qubits
results is that the DNM is not independent of the quantum system in the
context of which it is evaluated.

In closing, we believe that the appropriate modifications of the QSD measures,
that extend their use to non-trace-preserving processes, which are inevitable
in open systems, and the detailed analysis of their applicability to a qubit
chain system, have revealed subtle aspects of non-Markovianity. Extension of
our approach and analysis to other open systems can be expected to offer new
insights into the performance of protocols pertaining to quantum information
technology tasks as well as to the statistical mechanics of spin chains. Even
for the single chain coupled to non-Markovian reservoirs, issues such as the
dynamics of entanglement, and the DNM beyond single excitation represent
potentially fruitful territory for exploration.

## Acknowledgments

GM would like to acknowledge the Hellenic Foundation for Research and
Innovation (HFRI) for financially supporting this work under the 3rd Call for
HFRI PhD Fellowships (Fellowship Number: 5525).

##  APPENDIX A Derivations of $R(t)$ and $B(s)$ for special types of non-
Markovian reservoirs

In this appendix we provide the analytical derivations of the functions
$R(t)\equiv\int_{0}^{\infty}J(\omega)e^{-i\left(\omega-\omega_{eg}\right)t}d\omega$
and their Laplace transforms $B(s)$ for various forms of non-Markovian
spectral densities $J(\omega)$. The calculation of $B(s)$ in each case is
necessary in order to invert the Laplace transforms $F_{i}(s)$, $i=1,\dots,N$
and find the expressions of the tilde amplitudes of the chain.

###  A.1 Lorentzian spectral density

For reservoirs characterized by Lorentzian spectral densities, $J(\omega)$ is
given by:

$J(\omega)=\frac{g^{2}}{\pi}\frac{\frac{\gamma}{2}}{\left(\omega-\omega_{c}\right)^{2}+(\frac{\gamma}{2})^{2}},$
(A.1)

where $g$ is the coupling strength between the last qubit of the chain and the
reservoir, while $\gamma$ and $\omega_{c}$ is the width and the peak of the
distribution, respectively.

An approximation that considerably simplifies analytically the calculation of
the function $R(t)$ is the extension of the lower limit of the integration
over frequency from 0 to $-\infty$. Note that such extension is not in general
valid for any spectral density; it is however well justified in the case of a
Lorentzian with positive peak frequency and width such that the distribution
has practically negligible extension to negative frequencies. The necessary
condition for this approximation is $\gamma$ $<<$ $\omega_{c}$. In this case
the frequency integral can be calculated analytically, yielding:

$R(t)=g^{2}e^{-i\Delta_{c}t}e^{-\frac{\gamma}{2}t},$ (A.2)

where $\Delta_{c}\equiv\omega_{c}^{j}-\omega_{eg}$. The Laplace transforms of
$R(t)$ is then given by the expression:

$B(s)=\frac{g^{2}}{s+\frac{\gamma}{2}+i\Delta_{c}}.$ (A.3)

###  A.2 Lorentzian squared spectral density

The Lorentzian squared spectral density is a special case of a broader set of
spectral densities referred as powers of Lorentzian that are given by the
form:

$J(\omega)=g^{2}\mathcal{N}_{n}\frac{\left(\frac{\gamma}{2}\right)^{2n-1}}{\left[\left(\omega-\omega_{c}\right)^{2}+(\frac{\gamma}{2})^{2}\right]^{n}},$
(A.4)

where $\mathcal{N}_{n}$ is an appropriate normalization factor that depends
upon the non-zero integer $n$. For $n=1$ we capture the case of the Lorentzian
spectral density, while for $n=2$ we capture the Lorentzian squared spectral
density given by:

$J(\omega)=\frac{2g^{2}}{\pi}\frac{\left(\frac{\gamma}{2}\right)^{3}}{\left[\left(\omega-\omega_{c}\right)^{2}+(\frac{\gamma}{2})^{2}\right]^{2}},$
(A.5)

Following the same procedure as with the case of the Lorentzian, we find that
the function $R(t)$ is given by:

$R(t)=g^{2}\left(1+\frac{\gamma
t}{2}\right)e^{-i\Delta_{c}t}e^{-\frac{\gamma}{2}t}.$ (A.6)

Notice that Eqn. (A.6) is the same as Eqn. (A.2) with an addition of a
multiplication factor $\left(1+\frac{\gamma t}{2}\right)$. The Laplace
transform of Eqn. (A.6) reads:

$B(s)=g^{2}\frac{\left(s+\gamma+i\Delta_{c}\right)}{\left(s+\frac{\gamma}{2}+i\Delta_{c}\right)^{2}}$
(A.7)

###  A.3 Ohmic spectral density

The Ohmic spectral density is given by the expression:

$J(\omega)=\mathcal{N}g^{2}\omega_{c}\left(\frac{\omega}{\omega_{c}}\right)^{\mathcal{S}}\exp\left(-\frac{\omega}{\omega_{c}}\right),$
(A.8)

where $\mathcal{S}$ is the Ohmic parameter, characterizing whether the
spectrum of the reservoir is sub-Ohmic ($\mathcal{S}<1$), Ohmic ($\mathcal{S}$
= 1) or super-Ohmic ($\mathcal{S}>1$), $\omega_{c}$ is the Ohmic cut-off
frequency and $\mathcal{N}$ is a normalization constant given by the relation
$\mathcal{N}=\frac{1}{\left({\omega_{c}}\right)^{2}\Gamma\left(1+\mathcal{S}\right)}$,
where $\Gamma(z)$ is the gamma function.

Substituting Eqn. (A.8) back to the expression of $R(t)$ and calculation of
the frequency integral yields:

$R(t)=g^{2}e^{i\omega_{eg}t}\left(i\omega_{c}t+1\right)^{-1-\mathcal{S}}$
(A.9)

The Laplace transform of $R(t)$ can then be calculated analytically, yielding:

$B(s)=-g^{2}\frac{i^{1-\mathcal{S}}}{\omega_{c}}e^{-iK(s)}\left[K(s)\right]^{\mathcal{S}}\Gamma\left(-\mathcal{S},-iK(s)\right),$
(A.10)

where $K(s)\equiv\left(s-i\omega_{eg}\right)/\omega_{c}$ and $\Gamma(a,z)$ is
the incomplete gamma function.

## References

  * (1) S. Paganelli, S. Lorenzo, T. J. G. Apollaro, F. Plastina, and G. L. Giorgi, Phys. Rev. A 87, 062309 (2013).
  * (2) J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi, Phys. Rev. Lett. 78, 3221 (1997).
  * (3) M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. 92, 187902, (2004).
  * (4) C. Di Franco, M. Paternostro, and M. S. Kim, Phys. Rev. Lett. 101, 230502 (2008).
  * (5) D. Petrosyan, G. M. Nikolopoulos, and P. Lambropoulos, Phys. Rev. A 81, 042307 (2010).
  * (6) N. Y. Yao, L. Jiang, A. V. Gorshkov, Z.-X. Gong, A. Zhai, L.-M. Duan, and M. D. Lukin, Phys. Rev. Lett. 106, 040505 (2011).
  * (7) A. K. Pavlis, G. M. Nikolopoulos, and P. Lambropoulos, Quantum Inf. Process. 15, 2553–2568 (2016).
  * (8) A. Gratsea, G. M. Nikolopoulos, and P. Lambropoulos, Phys. Rev. A 98, 012304 (2018).
  * (9) M. B. Plenio, J. Hartley and J. Eisert, New. J. Phys. 6, 36 (2004).
  * (10) S. Bose, Phys. Rev. Lett. 91, 207901 (2003).
  * (11) L. Campos Venuti, S. M. Giampaolo, F. Illuminati, and P. Zanardi, Phys. Rev. A 76, 052328 (2007).
  * (12) S. M. Giampaolo and F. Illuminati, New J. Phys. 12, 025019 (2010).
  * (13) S. Sahling, G. Remenyi, C. Paulsen, P. Monceau, V. Saligrama, C. Marin, A. Revcolevschi, L. P. Regnault, S. Raymond, and J. E. Lorenzo, Nature Physics 11, 255–260 (2015).
  * (14) M. P. Estarellas, I. D’Amico, and T. P. Spiller, Phys. Rev. A 95, 042335 (2017).
  * (15) F. Plastina and T. J. G. Apollaro, Phys. Rev. Lett. 99, 177210 (2017).
  * (16) G. Mouloudakis, T. Ilias, and P. Lambropoulos, Phys. Rev. A 105, 012429 (2022).
  * (17) H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009).
  * (18) E.-M. Laine, J. Piilo, and H.-P. Breuer, Phys. Rev. A, 81, 062115 (2010).
  * (19) D. Chruściński, A. Kossakowski, and Á. Rivas, Phys. Rev. A 83, 052128 (2011).
  * (20) S. Wißmann, H.-P. Breuer, and B. Vacchini, Phys. Rev. A 92, 042108 (2015).
  * (21) L. Li, M. J. W. Hall, H. M. Wiseman, Phys. Rep. 759, 1–51 (2018).
  * (22) S. Luo, S. Fu, and H. Song, Phys. Rev. A 86, 044101 (2012).
  * (23) S. Lorenzo, F. Plastina, and M. Paternostro, Phys. Rev. A 88, 020102(R) (2013).
  * (24) H.-P. Breuer, J. Phys. B: At. Mol. Opt. Phys. 45, 154001 (2012).
  * (25) P. Haikka, J. D. Cresser, and S. Maniscalco Phys. Rev. A 83, 012112 (2011).
  * (26) A. Mortezapour, M. A. Borji, D. Park and R. L. Franco, Open Sys. Inf. Dyn. 24, 1740006 (2017).
  * (27) C. Addis, B. Bylicka, D. Chruściński, and S. Maniscalco, Phys. Rev. A 90, 052103 (2014).
  * (28) M. Mannone, R. L. Franco, and G. Compagno, Phys. Scr. 2013, 014047 (2013).
  * (29) F. F. Fanchini, G. Karpat, D. Z. Rossatto, A. Norambuena, and R. Coto, Phys. Rev. A 103, 022425 (2021).
  * (30) Á. Rivas, S. F. Huelga, M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * (31) J. Dajka, J. Łuczka, and Peter Hänggi, Phys. Rev. A 84, 032120 (2011).
  * (32) B. Bylicka, D. Chruściński, and S. Maniscalco, Scientific Reports 4, 5720 (2014).
  * (33) M. Ozawa, Phys. Lett. A 268, 158-160 (2000).
  * (34) F.-H. Ren, Z.-M. Wang, and Y.-J. Gu, Q. Inf. Process. 18, 193 (2019).
  * (35) R. Jozsa, J. Mod. Opt. 41, 2315 (1994).
  * (36) A. Uhlmann, Phys. Rev. A 62, 032307 (2000).
  * (37) P. E. M. F. Mendonça, R. d. J. Napolitano, M. A. Marchiolli, C. J. Foster, and Y.-C. Liang, Phys. Rev. A 78, 052330 (2008).
  * (38) P. Lambropoulos and D. Petrosyan, Fundamentals of Quantum Optics and Quantum Information (Springer Verlag, 2007).
  * (39) F. Shibata, Y. Takahashi, and N. Hashitsume, J. Stat. Phys. 17, 171–187 (1977).
  * (40) H. Z. Shen, M. Qin, X.-M. Xiu, and X. X. Yi Phys. Rev. A 89, 062113 (2014).
  * (41) L. Ferialdi, Phys. Rev. A 95, 020101(R) (2017); Phys. Rev. A 95, 069908 (2017).
  * (42) S. Maniscalco and F. Petruccione, Phys. Rev. A 73, 012111 (2006); Phys. Rev. A 73, 029902 (2006); Phys. Rev. A 75, 059905 (2007).
  * (43) B. Vacchini, Phys. Rev. A 87, 030101(R) (2013).
  * (44) See, for example, H. P. Breuer and F. Petruccione, The theory of open quantum systems, (Oxford University Press, Oxford, 2002).
  * (45) G. Mouloudakis and P. Lambropoulos, Phys. Rev. A Phys. Rev. A 106, 053709 (2022).
  * (46) U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 2001).
  * (47) G. Mouloudakis and P. Lambropoulos, Quantum Inf. Process. 20, 331 (2021).
  * (48) E. R. Whipple, Nucl. Inst. Meth. 180, 241-248 (1981).
  * (49) A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, Anupam Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987); Rev. Mod. Phys. 67, 725 (1995).
  * (50) H.-M. Zou, R. Liu, D. Long, J. Yang, and D. Lin, Phys. Scr. 95, 085105 (2020).
  * (51) M. Abdi and M. B. Plenio, Phys. Rev. A 98, 040303(R) (2018).
  * (52) D. Pellegrino et al., Phys. Rev. Lett. 124, 123902 (2020).


# Conserving elastic turbulence numerically using artificial diffusivity

V. Dzanic  School of Mechanical, Medical, and Process Engineering, Faculty of
Engineering, Queensland University of Technology, QLD 4001, Australia  C. S.
From<EMAIL_ADDRESS>Department of Chemical Engineering,
University of Manchester, Manchester M13 9PL, UK  E. Sauret
<EMAIL_ADDRESS>School of Mechanical, Medical, and Process
Engineering, Faculty of Engineering, Queensland University of Technology, QLD
4001, Australia

###### Abstract

To simulate elastic turbulence, where viscoelasticity dominates, numerical
solvers introduce an artificial stress diffusivity term to handle the steep
polymer stress gradients that ensue. This has recently been shown [Gupta &
Vincenzi, J. Fluid Mech. 870, 405-418 (2019); Dzanic et al., J. Fluid Mech.
937, A31 (2022)] to introduce unphysical artifacts with a detrimental impact
on simulations. In this Letter, we propose that artificial diffusion is
limited to regions where stress gradients are steep instead of seeking the
zero-diffusivity limit. Through the cellular forcing and four-roll mill
problem, we demonstrate that this modified artificial diffusivity is devoid of
unphysical artifacts, allowing all features of elastic turbulence to be
retained. Results are found to conform with direct simulations while only
requiring a fraction of the numerical resolution.

It is well-known that fluid flow with the addition of polymer molecules to a
solvent (i.e., viscoelastic fluids) in the absence of inertial instabilities
$Re=\rho U\ell/\mu\leq 1$ generates an anisotropic elastic stress contribution
that transitions the flow to a novel chaotic regime, known as elastic
turbulence (ET) [1, 2, 3]. ET is purely driven by viscoelastic instabilities,
where the viscous to elastic effects are measured by the Weissenberg number
$Wi=\tau U/\ell\gg 1$, and viscoelastic and inertial effects characterized by
the elasticity number $El\equiv Wi/Re=\tau\mu/(\rho\ell^{2})\gg 1$. Here,
$\tau$ is the longest polymer relaxation time, $\ell$ characteristic length
scale, $\rho$ is the density, $\mu$ is the total viscosity, and $U$ is the
average velocity. This purely elastic instability shares a lot of features
with traditional inertial turbulence, namely, (i) increased flow resistance,
(ii) enhanced mixing, (iii) random flow fluctuations characterized by a
broadband spectrum of spatial and temporal frequencies [3]. In an effort to
better understand the role of viscoelastic instabilities in these novel
chaotic flow regimes, a plethora of numerical studies have been conducted [4,
5, 6, 7, 8, 9]. The majority of which, involve resolving the hydrodynamic
field through the incompressible Navier-Stokes equations,

$\bm{\nabla}\bm{\cdot}\bm{u}=0,~{}~{}\frac{D\mathbf{u}}{Dt}=-\bm{\nabla}P+\mu_{s}\mathbf{\Delta}\mathbf{u}+\bm{\nabla}\bm{\cdot}\bm{\sigma}+\mathbf{F},$
(1)

coupled with the the polymer stress tensor
$\bm{\sigma}=f(r)\frac{\mu_{p}}{\tau}\left(\mathsfbi{C}-\mathsfbi{I}\right)$,
described by a space-time dependent conformation tensor ($\mathsfbi{C}$)
constitutive equation,

$\frac{D\mathsfbi{C}}{Dt}=\mathsfbi{C}\bm{\cdot}\left(\bm{\nabla}\mathbf{u}\right)+\left(\bm{\nabla}\mathbf{u}\right)^{\mathsf{T}}\bm{\cdot}\mathsfbi{C}-\frac{f(r)}{\tau}\left(\mathsfbi{C}-\mathsfbi{I}\right),$
(2)

where the function $f(r)$ allows for various constitutive polymer models to be
used. $\mathsfbi{I}$ is the identity tensor, $\mathbf{u}$ is the velocity
field, $\mathbf{F}$ is the external force, $\mu_{s}$ and $\mu_{p}$ are the
solvent and polymer dynamic viscosity, respectively.

Numerical simulations of turbulent viscoelastic fluid flows is, however, far
from trivial. By definition, the conformation tensor is a symmetric positive-
definite (SPD) tensor $\mathsfbi{C}\gneq 0$, where conserving this property is
important to prevent the rapid growth of Hadamard instabilities [10]. To
overcome stability issues, numerical methods apply decomposition techniques in
Eq. (2) to conserve the SPD properties of $\mathsfbi{C}$ by construction [11,
12]. However, these specialized techniques alone cannot control the steep
polymer stress gradients that ensue at high elastic effects, captured by the
infamous high-$Wi$ number problem, owing to the inherent hyperbolic nature and
lack of numerical regularization terms in Eq. (2). In turn, solvers introduce
an additional global artificial diffusivity (GAD) term in the Laplacian form,
$\mathsfbi{D}=k\mathbf{\Delta}\mathsfbi{C}$, into Eq. (2) [8, 9, 13], where
the level of diffusivity $k$ is characterized with respect to the solvent
kinematic viscosity $\nu_{s}$ by the Schmidt number $Sc=\nu_{s}/k$. GAD
effectively converts the constitutive polymer model to a parabolic form by
smoothing steep polymer stresses over large regions of the flow. Indeed,
polymer-stress diffusion is physically present at $Sc\sim 10^{6}$; however, to
achieve numerical stability requires $Sc\lesssim 10^{3}$ [14, 8]. The
unphysically large $k$ values required to achieve such $Sc$ numbers with GAD
are known to promote laminarization for elasto-inertial turbulence (EIT)
($Re\gg 10$) [15], suppressing the elastic instabilities that are necessary to
promote turbulence. By combining SPD-conserving solvers with high-order,
shock-capturing schemes have enabled simulations of EIT without the need for
GAD [16]. This same approach was recently used to simulate ET at $Sc=\infty$,
where direct comparison against GAD with $Sc=10^{3}$ revealed the dramatic
effect on the large-scale properties of the flow, suppressing chaotic
fluctuations [8]. In the authors’ recent work [17], it was further shown that
the global effect of GAD is so large that a single unit cell in fully-periodic
problems does not conserve unicity, leading to numerical artifacts represented
by qualitative anomalies, thus rendering the problem unphysical. Conserving
unicity allows features of ET to be recovered at $Sc=10^{3}$, with chaotic
fluctuations increasing as $Sc>10^{3}$, whereas ET is suppressed at $Sc\sim
10^{2}$, demonstrating the detrimental effect of GAD. Indeed in the absence of
artificial diffusivity, an exact physical representation can be obtained
directly from Eq. (2), however, for problems where $El\gg 1$, the steep
polymer stress gradients that develop require significant resolutions, and
hence computational costs, to overcome the numerical stability issues that
ensue [7, 8, 17].

In this Letter, a simple yet elegant approach to include artificially
diffusivity while retaining all features of ET is proposed. We demonstrate the
ability to avoid the unphysical artifacts that have plagued the GAD scheme
while still offering the same level of numerical benefits, particularly in
terms of convenience and stability.

Earlier works on turbulent polymer drag reduction have seen success using a
local artificial diffusivity (LAD) scheme in which
$\mathsfbi{D}(\bm{x},t)=k\mathbf{\Delta}\mathsfbi{C}(\bm{x},t)$ is applied in
regions of the computational domain $\bm{x}$ where $\mathsfbi{C}\leq 0$ [18,
19]. However, LAD does not guarantee SPD and is inherently not applicable to
modern, SPD-conserving solvers [11, 12, 9]. Motivated by the LAD scheme and
stability rationale, i.e., the singular purpose of artificial diffusivity is
for numerical regularity, we seek a suitable approach with which $k$ is not
decreased but instead restricted where it is applied to minimize the potential
consequent global effects. To this aim, we propose a modified artificial
diffusivity (MAD) scheme such that
$\mathsfbi{D}(\bm{x},t)=\kappa(\bm{x},t)\mathbf{\Delta}\mathsfbi{C}(\bm{x},t)$
with $\max_{\bm{x}}[\kappa(\bm{x},t)]=k~{}\forall t$, where $\bm{x}$ is the
domain, which is directly applicable to all modern, SPD-conserving, numerical
solvers. To illustrate, an intuitive way to apply $\kappa$ is to consider the
characteristic feature of ET, that the stress gradients are sharp and
localised. By restricting artificial diffusivity within those regions to a
maximum of $k=\nu_{s}/Sc$, we allow for the required numerical regularity to
be retained while minimizing global smearing where regularization is not
needed. Leveraging this idea, MAD considers the Laplacian form,

$\mathsfbi{D}(\bm{x},t)=\frac{k\mathcal{Q}(\bm{x},t)}{\mathcal{Q}_{max}(t)}\mathbf{\Delta}\mathsfbi{C}(\bm{x},t),$
(3)

where $\mathcal{Q}$ is the sum of the magnitude of polymer stress component
gradients, i.e.,

$\mathcal{Q}(\bm{x},t)=\sum_{i,j}^{D}\sqrt{\sum_{q}^{D}\Big{[}\nabla_{q}\mathsfit{C}_{ij}(\bm{x},t)\Big{]}^{2}},$

for a $D$ dimensional domain $\bm{x}$ and
$\mathcal{Q}_{max}(t)=\max_{\bm{x}}[\mathcal{Q}(\bm{x},t)]$ is a normalisation
factor. In turn, in MAD (3) the artificial diffusivity is scaled linearly with
stress gradients to a maximum of $k$, i.e.,
$\lim_{\mathcal{Q}\rightarrow\mathcal{Q}_{max}}\mathsfbi{D}=k\mathbf{\Delta}\mathsfbi{C}$.

Figure 1: Representative snapshots of $\ln($tr$\mathsfbi{C})$ for FRM (a,b)
with $n=2$ and CF (c,d) using GAD (a,c) and MAD (b,d). The red borders in
(a,b) illustrate a single unit cell for FRM. For both cases, simulations with
GAD (a,c) result in a loss of symmetry spawned from artificial diffusivity,
whereas with MAD (b,d) they are largely constrained to the background forcing
symmetry. Figure 2: The $|$tr$\mathsfbi{D}|$ obtained with (a) GAD
$\mathsfbi{D}=k\mathbf{\Delta}\mathsfbi{C}$ and (b) MAD [Eq. (3)] using
$\mathsfbi{C}$ from the simulation of CF with MAD in Fig. 1 (d). With MAD,
$\mathsfbi{D}$ is largely concentrated at steep stress gradients where
$\forall\kappa\gneq 0,~{}Sc=10^{3}\sim 10^{8}$. Figure 3: Time series of the
dimensionless mean kinetic energy $\langle E_{k}\rangle_{V}/E_{k,0}$ for the
(a) FRM and (b) CF using the GAD (blue) and MAD (red) schemes. Note, insets
provided in (a) correspond to the snapshots of $\ln($tr$\mathsfbi{C})$ at
locations A ($t=400T$) and B ($t=2500T$) for the GAD scheme, illustrating the
different transition points. The inset provided in (b) is a zoom-in
corresponding to the late-time dynamics in the dashed box region.

To demonstrate the suitability of the MAD scheme while further demonstrating
the severe limitations of the GAD scheme for simulations of ET, we apply both
MAD and GAD to simulate two popular ET cases, namely, the four-roll mill (FRM)
problem [20, 21, 22] and the cellular forcing scheme (CF) [7, 8]. In recent
investigations of ET, simulations of FRM and CF using GAD were shown to depict
various numerical artifacts that have been investigated extensively in both
[17] and [8], respectively, which presents strict test criteria for MAD. Both
FRM and CF are numerical recreations of popular viscoelastic experimental
cases [23, 24, 25], and are solved in a 2D domain $\bm{x}=[0,~{}n\times
2\pi]^{2}$ with double periodic boundary conditions (PBCs) where a single unit
cell is $\left[0,2\pi\right]^{2}$. Here, $n$ is the level of periodicity, for
which $n>1$ results in $n^{2}$ unit cells to be solved. The experimental
effect of rollers, which create an elongational flow regime, are mimicked
through a constant external forcing, which for the FRM forcing scheme is given
by

$\mathbf{F}(\bm{x})=F_{0}\left(2\sin\left(K{x}\right)\cos\left(K{y}\right),~{}-2\cos\left(K{x}\right)\sin\left(K{y}\right)\right),$
(4)

and for the CF forcing scheme,

$\mathbf{F}(\bm{x})=F_{0}\left(-\sin\left(K{y}\right),~{}\sin\left(K{x}\right)\right),$
(5)

in which the forcing amplitude is $F_{0}=U\nu_{s}K^{2}$ and $K$ is the spatial
frequency (i.e., $\ell=1/K$), resulting in a turnover time $T=\nu_{s}K/F_{0}$.
For FRM (4) $K=1$ and for CF (5) $K=2$. To simulate ET using FRM and CF, a
small perturbation $\bm{\delta}$ is added to the initial conformation tensor
$\mathsfbi{C}=\mathsfbi{I}+\bm{\delta}$, as originally proposed in [21].
Equations (1) and (2), with GAD and MAD (3), are solved using a SPD-conserving
numerical solver comprising of the lattice Boltzmann coupled with a high-order
finite-difference scheme (see [26]), which was applied in our previous
investigation of ET [17]. We realistically capture the main physical behavior
of polymers using the FENE-P constitutive model,
$f(r)=\left(L^{2}-2\right)/\left(L^{2}-r\right)$ in Eq. (2) with
$r=\mathrm{tr}\left(\mathsfbi{C}\right)$ and $L^{2}=2.5\times 10^{3}$, which
imposes a maximum finite extensibility $L>r$ [27]. In all simulations,
dimensionless groups are set following previous numerical investigations of ET
including, $Wi=10$, artificial diffusivity $Sc=10^{3}$ (as done in [21, 22,
8]), and set $Re=1$ (as done in [5, 7, 17]) below the critical value at which
inertial instabilities arise $Re_{c}=\sqrt{2}$ [28]. The polymer
concentration, characterized by $\beta=\nu_{p}/\nu_{s}$, is set to $\beta=0.5$
and $\beta=0.4$ in order to match previous investigations for FRM [21, 22] and
CF[8], respectively.

Based on our previous investigation [17], GAD with $n=1$ will result in
unphysical artifacts arising from insufficient periodicity, characterized by
qualitative anomalies which contaminate the base flow, thus preventing the
accurate simulation of ET. For CF, we set $n=1$ to compare against results
from [8], which were solved on a $N^{2}=1024^{2}$ grid for both $Sc=10^{3}$
(using GAD) and $Sc=\infty$. In addition, this is useful for testing whether
unicity can be retained. For FRM, where the external forcing is larger, we set
$n=2$ to conserve unicity allowing numerical artifacts arising purely from
artificial diffusivity to be studied. All simulations are purposely conducted
with $(n\times N)^{2}=256^{2}$ grid points, where each unit cell has $N^{2}$
grid points, i.e., the resolution $2\pi/N$ for FRM is half of that used for
CF. This is a notable reduction compared to previous studies [7, 8] to
highlight the ability to retain numerical robustness with the proposed
approach of modified artificial diffusivity.

Figure 4: The temporal power spectral density $E(f)$ of the dimensionless mean
kinetic energy fluctuations for the (a) FRM and (b) CF using the GAD (blue)
and MAD (red) schemes. Note, results with MAD follow a power law, $E\propto
f^{-3.44}$ and $E\propto f^{-3.39}$ behavior for the FRM and CF cases,
respectively. Alternatively, with GAD follows a steeper power law, $E\propto
f^{-4.31}$ behavior for the FRM, while failing to follow a power law behavior
for the CF.

In Fig. 1, contour plots of the polymer field $\ln($tr$\mathsfbi{C})$ using
the GAD and MAD are compared for both the FRM [Fig. 1 (a) and (b)] and CF
[Fig. 1 (c) and (d)] cases. The stark differences are immediately clear on a
qualitative-scale. More specifically, it can be seen that FRM with GAD leads
to a single-leading vortex within each unit cell [note, a unit cell is
characterized by four rollers as illustrated by the red borders in Fig. 1 (a)
and (b)], which appears during the onset of ET. Similarly, CF with GAD
eventually transitions from the initial background forcing symmetry into a
state with which dynamics are periodic and then into the double-leading
vortices in Fig. 1 (c), where dynamics become aperiodic. Analogous results are
obtained for FRM at $n=1$ using GAD [21]. Ultimately, these are a product of
the GAD excessively spreading polymer stresses over large regions of the flow,
even within the vortical regions where there should be no polymer stretching,
as illustrated in Fig. 2 (a). In turn, the polymer field eventually
destabilizes and loses the initial forcing symmetry [8]. The initial breakdown
in symmetry during the onset of ET has been shown to be independent of
periodicity [17]. Furthermore, considering the same observations were made for
CF with GAD [8] but at significantly larger resolution ($N^{2}=1024^{2}$)
demonstrates the grid-independence of numerical artifacts arising from
artificial diffusivity and periodicity. On the other hand, numerical
simulations of FRM and CF with MAD are devoid of such numerical artifacts. For
CF with MAD, despite experiencing momentary losses of symmetry, the polymer
stretching is mostly constrained to the background forcing symmetry, in
qualitative agreement with the results obtained using direct simulations at
$Sc=\infty$ in [8] (see animations in the Supplementary Materials). Notably,
the momentary losses of symmetry are expected to be eradicated with $n>1$
[17], nevertheless, we shall see that these have a negligible impact on the
ability to simulate ET. Essentially, MAD overcomes numerical artifacts by
concentrating $\mathsfbi{D}$ only within the steep polymer stress gradients
[Fig. 2 (b)], whereas with GAD $\mathsfbi{D}$ is diffused over all gradients
[Fig. 2 (a)].

The quantitative results further reflect the stark differences observed for
the two artificial diffusivity schemes. The time history response for the
dimensionless mean kinetic energy $\langle E_{k}\rangle_{V}/E_{k,0}$, where
$\langle\cdot\rangle_{V}$ denotes the spatial average and $E_{k,0}=U^{2}/2$,
is shown for both FRM and CF in Fig. 3. Notably, the time series of the
spatially-averaged polymer trace is also provided in the Supplementary
Material and shows analogous results for both FRM and CF. In Fig. 3 (a), FRM
with GAD initially reaches a steady-state, an artifact due to GAD suppressing
intrinsic instabilities, followed by an onset of ET for which there is a rapid
loss of symmetry [inset point A in Fig. 3 (a)], which is characterized by slow
oscillations. In our previous study [17], we showed that a complex interplay
exists between periodicity and artificial diffusivity, whereby the classic FRM
problem, in which $n=1$, is unable to overcome this initial loss of symmetry,
as similarly observed for CF with GAD in Fig. 3 (b). Imposing sufficient
$n$-levels of periodicity, as done here for FRM with $n=2$, in turn, conserves
unicity [e.g., the background forcing symmetry, shown in inset point B in Fig.
3 (a)] and enables features of ET to be recovered, such as chaotic flow
fluctuations. Notably, defining $n>2$ increases the rate at which ET is
recovered after the initial loss of symmetry [17]. Remarkably, with MAD, the
FRM results are devoid of all unphysical artifacts spawned from artificial
diffusivity. When observing the results for CF with GAD in Fig. 3 (b),
multiple instability modes exist. The first mode occurs within the early
stages during an initial partial loss of flow symmetry, which periodically
cycles between two different states, reflected by the fully periodic dynamics.
This is the first observable numerical artifact caused by artificial
diffusivity [8] and PBCs with limited periodicity [17]. Qualitatively, this
instability mode can be appreciated in the animation provided in the
Supplementary Materials. Beyond this state, the flow experiences a complete
loss of symmetry, dominated by two vortices [refer to Fig 1 (c)] where the
dynamics transition into an aperiodic state, as previously reported in [8] and
analogous results can be observed for FRM with $n=1$ [21, 22, 26]. While the
initial loss of symmetry is a consequence of artificial diffusivity, the
periodic and aperiodic states are largely due to the limited periodicity,
where conserving unicity, as done here for FRM with GAD, prevents such
unphysical behavior [17]. On the other hand, CF with MAD [Fig. 3 (b)] almost
instantly transitions into a fully chaotic state, where the heavy fluctuations
perturb the cellular vortices, however, the large-scale structure remains
largely constrained to the background forcing with momentary losses of
symmetry [refer to Fig. 1 (d)]. In comparing the late-time dynamics within the
statistically homogeneous state (approximately $t>500T$, for both FRM and CF),
it is noticeable that for CF with MAD, fluctuations are at a much higher
frequency compared to the slow aperiodic state observed for GAD. Similar
observations were made in comparing GAD against the direct solution (i.e.,
$Sc=\infty$) [8]. Thus, the excessive diffusivity from GAD suppresses the
chaotic fluctuations characteristic of the ET regime, a well-reported artifact
of artificial diffusivity [8, 17]. High-frequency oscillations with MAD occur
almost instantly for both FRM and CF, which is in agreement with the direct
simulation of CF at $Sc=\infty$ [8] and experimental analog of FRM [25].

In Fig. 4, we examine the temporal power-law spectrum of the time signals from
Fig. 3 within the statistically homogeneous regime. A distinctive feature of
ET is a fairly steep power-law spectrum of velocity fluctuations [3, 29]. The
fluctuations for FRM with $n=2$ using GAD [Fig. 4 (a)] follows a steep power-
law with an exponent of $E\propto f^{-4.31}$. This steepness is independent of
periodicity for $n\gneq 1$, [17] and the results here for FRM with MAD confirm
this to be attributed to the numerical artifacts. The aperiodic state observed
for CF with GAD [Fig. 4 (b)] fails to follow any apparent power law, as
originally observed in [8] with significantly greater resolution. This is
attributed to GAD suppressing the high-wavenumber fluctuations of polymer
stresses and is analogously observed for FRM with $n=1$ [17]. Ultimately, with
GAD and $n=1$, numerical simulations of FRM and CF are unable to conserve any
features of ET. On the other hand, for both cases, the fluctuations retained
in simulations with MAD behave as a power law, with exponents similar to
previous experimental and numerical studies of ET, in which the decay rate
varied with the setup, but the exponent was always smaller than $-3$ [1, 2, 5,
30, 3, 29]. The strong retention of ET features using MAD is a result of
restricting artificial diffusivity within the critical regions of the flow
[i.e., steep polymer stress gradients, as shown in Fig. 2], thus avoiding the
consequent global deviation from the exact polymer representation (i.e.,
$Sc=\infty$). The remarkable feature of MAD is the numerical robustness while
retaining all features of ET with an absence of numerical artifacts despite
using 16 times less collocation points (i.e., a quarter resolution) compared
to that required by direct simulations of ET at $Sc=\infty$ [8, 7].

Summarising, a new view on including numerical regularity for simulating ET
through a modified artificial diffusivity (MAD) has been proposed. Its
applicability is demonstrated by applying the MAD to two stringent numerical
experiments. It is shown that with MAD, all characteristic features of ET can
be simulated while overcoming the recently discovered unphysical numerical
artifacts of the global artificial diffusivity used traditionally. The impact
of artificial diffusion with MAD is effectively reduced from a qualitative
scale to a quantitative scale, requiring a numerical grid that is one order of
magnitude smaller compared to that required to numerically simulate ET at
$Sc=\infty$. It is difficult to overlook the apparent trend and analogy that
artificial diffusivity schemes share with traditional turbulence modeling. In
a similar fashion, the MAD scheme is an additional numerical tool that offers
numerical features that could be of paramount importance in progressing
towards more complex and computationally expensive ET cases, notably, 3D
simulations.

The authors acknowledge the High-Performance Computing facilities at QUT.
Prof. E. Sauret is the recipient of an Australian Research Council Future
Fellowship (FT200100446) funded by the Australian Government. V. Dzanic
gratefully acknowledges QUT for support through a Ph.D. scholarship.

## References

  * Groisman and Steinberg [2000] A. Groisman and V. Steinberg, Nature 405 (2000), 10.1038/35011019.
  * Groisman and Steinberg [2004] A. Groisman and V. Steinberg, New J. of Phys. 6 (2004), 10.1088/1367-2630/6/1/029.
  * Steinberg [2021] V. Steinberg, Ann. Rev. of Fluid Mech. 53, 27 (2021).
  * Poole  _et al._ [2007] R. J. Poole, M. A. Alves, and P. J. Oliveira, Phys. Rev. Lett. 99, 164503 (2007).
  * Berti  _et al._ [2008] S. Berti, A. Bistagnino, G. Boffetta, A. Celani, and S. Musacchio, Phys. Rev. E 77, 055306 (2008).
  * Grilli  _et al._ [2013] M. Grilli, A. Vázquez-Quesada, and M. Ellero, Phys. Rev. Lett. 110, 174501 (2013).
  * Plan  _et al._ [2017] E. L. C. V. M. Plan, A. Gupta, D. Vincenzi, and J. D. Gibbon, J. of Fluid Mech. 822, R4 (2017).
  * Gupta and Vincenzi [2019] A. Gupta and D. Vincenzi, J. of Fluid Mech. 870, 405–418 (2019).
  * Alves  _et al._ [2021] M. Alves, P. Oliveira, and F. Pinho, Ann. Rev. of Fluid Mech. 53, 509 (2021).
  * Sureshkumar and Beris [1995] R. Sureshkumar and A. N. Beris, J. of Non-Newtonian Fluid Mech. 60, 53 (1995).
  * Vaithianathan and Collins [2003] T. Vaithianathan and L. R. Collins, J. of Comp. Phys. 187, 1 (2003).
  * Fattal and Kupferman [2005] R. Fattal and R. Kupferman, J. of Non-Newtonian Fluid Mech. 126, 23 (2005).
  * Lee  _et al._ [2021] J. Lee, W. R. Hwang, and k. s. Cho, J. of Non-Newtonian Fluid Mech. 297, 104650 (2021).
  * Vaithianathan  _et al._ [2006] T. Vaithianathan, A. Robert, J. G. Brasseur, and L. R. Collins, J. of Non-Newtonian Fluid Mech. 140, 3 (2006).
  * Sid  _et al._ [2018] S. Sid, V. E. Terrapon, and Y. Dubief, Phys. Rev. Fluids 3, 011301 (2018).
  * Shekar  _et al._ [2020] A. Shekar, R. M. McMullen, B. J. McKeon, and M. D. Graham, J. of Fluid Mech. 897, A3 (2020).
  * Dzanic  _et al._ [2022a] V. Dzanic, C. S. From, and E. Sauret, J. of Fluid Mech. 937, A31 (2022a).
  * Min  _et al._ [2001] T. Min, J. Yoo, and H. Choi, J. of Non-Newtonian Fluid Mech. 100, 27 (2001).
  * Dubief  _et al._ [2005] Y. Dubief, V. Terrapon, C. White, E. Shaqfeh , P. Moin, and S. Lele, Flow Turbul. Combust. 74, 311 (2005).
  * Thomases and Shelley [2007] B. Thomases and M. Shelley, Phys. Fluids 19 (2007), 10.1063/1.2783426.
  * Thomases and Shelley [2009] B. Thomases and M. Shelley, Phys. Rev. Lett. 103, 094501 (2009).
  * Thomases [2011] B. Thomases, J. of Non-Newtonian Fluid Mech. 166, 1221 (2011).
  * Cardoso  _et al._ [1994] O. Cardoso, D. Marteau, and P. Tabeling, Phys. Rev. E 49, 454 (1994).
  * Rothstein  _et al._ [1999] D. M. Rothstein, E. Henry, and J. P. Gollub, Nature 401, 770 (1999).
  * Liu  _et al._ [2012] B. Liu, M. Shelley, and J. Zhang, J. of Non-Newtonian Fluid Mech. 175-176, 38 (2012).
  * Dzanic  _et al._ [2022b] V. Dzanic, C. S. From, and E. Sauret, Computers & Fluids 235, 105280 (2022b).
  * Peterlin [1961] A. Peterlin, Polymer 2, 257 (1961).
  * Gotoh and Yamada [1984] K. Gotoh and M. Yamada, J. of the Phys. Society of Japan 53, 3395 (1984).
  * Steinberg [2019] V. Steinberg, Phys. Rev. Lett. 123, 234501 (2019).
  * Canossi  _et al._ [2020] D. O. Canossi, G. Mompean, and S. Berti, Europhys. Lett. 129, 24002 (2020).