ilcssworkingpaper

# Modeling the 15 Tile Puzzle Through the Lens of Group Theory

Viren Khandal

###### Abstract

In this piece, we examine one variant of the infamous 15 Tile Puzzle and
develop a mathematical backing behind why it is unsolvable. Using concepts of
permutations, bijectivity, and cycle transpositions, we not only prove how to
model this puzzle as a group, but also how to determine if a certain
configuration of the puzzle is reachable.

December 10, 2021  The 15 tile puzzle has been a prevalent brain teaser in the
realm of mathematics for several decades. A new variant, invented in the
mid-1870s by Sam Lloyd, was known to drive the world crazy as many deemed it
impossible to solve. In fact, even though there was a $1000 prize was rewarded
to any solvers, no one claimed the prize, stating that the puzzle was
unsolvable.

The below Figure 1 shows an example of the starting position in Lloyd’s 15
tile puzzle. The board consists of a 4 by 4 grid, in which each cell but one
contains a tile, represented by number ranging from 1 to 15. The goal of the
puzzle can be seen as shifting the tiles in an orderly fashion to achieve the
goal state, which is displayed in Figure 2. Tiles can be moved vertically and
horizontally into the empty space in order to reach the aforementioned goal
state.

Figure 1: An example of Sam Lloyd’s 15 slide puzzle.

Figure 2: An example of the solved 15 slide puzzle.

## Discussion

In this section, we discuss the connection of the 15 Tile Puzzle to some of
the concepts covered in the MATH 113: Introduction to Abstract Algebra course
at the University of California, Berkeley.

Note: In this works, the words $mapping$ and $permutation$ are used
interchangeably.

We can model the placement of tiles on the 4x4 grid using a mapping C such
that:

$C:\\{1,2,...,16\\}\xrightarrow{}\\{1,2,...,16\\}$ (1)

Here, C maps each position on the grid labelled 1 through 16 to a tile
labelled 1 to 16, where 16 is the empty tile. Using this mapping C, we can
rewrite a configuration of the below shown board, $A$, as a permutation in two
line and cycle notation:

Figure 3: A: An example of a randomly configured board.

$\displaystyle A=\begin{pmatrix}{1\atop 3}\;{2\atop 2}\;{3\atop 13}\;{4\atop
9}\;{5\atop 6}\;{6\atop 7}\;{7\atop 12}\;{8\atop 5}\;{9\atop 10}\;{10\atop
11}\;{11\atop 8}\;{12\atop 4}\;{13\atop 15}\;{14\atop 14}\;{15\atop
1}\;{16\atop 16}\end{pmatrix}$ (2) $\displaystyle
A=(1\;3\;13\;15)\;(2)\;(4\;9\;10\;11\;8\;5\;6\;7\;12)\;(14)\;(16)$

Further, we can model the movement of tiles as a mapping $M_{1}$ such that:

$M_{1}:\\{1,2,...,16\\}\xrightarrow{}\\{1,2,...,16\\}$ (3)

Here, $M_{1}$ maps each tile labelled 1 through 16, where 16 is the empty
tile, to a position on the grid labelled 1 through 16. The set of all such
mappings, $M$, is a group under composition of moves111Trappa, Peter (2004)
_Permutations and the 15-Puzzle_ , University of Utah, Mathematics..

Using some configuration of the board, $C$, and a permutation of moves, $M$,
we can calculate the configuration of the board, $C^{\prime}$, after $M$ is
performed on $C$:

$C^{\prime}=M*C$ (4)

Starting with some configuration C, we can aggregate a series of moves
$M_{1},M_{2},...,M_{n}$ and formulate a new configuration C’ using:

$C^{\prime}=M_{n}*...*M_{2}*M_{1}*C$ (5)

One key detail to note is that moving a tile into the position occupied by the
empty square has the effect of moving the empty square one spot along to the
former location of the tile that was moved. Thus far, we have formulated both
M and C as permutations in $S_{16}$. Phrased in terms of material from MATH
113, the moves performed on a certain configuration is equivalent to defining
it as a bijection, which is an element of $S_{16}$, on the configuration.

######  Claim 0.1.

Performing a move on a certain configuration of the 15 Tile Puzzle is a
bijection on the set of configurations and thus the elements of M are also
elements of $S_{16}$.

###### Proof.

To prove a certain operation is bijective, we must first prove it is, both,
injective (one-to-one) and surjective (onto).

We begin by proving the injectivity of the operation. Assume $C_{1}$ and
$C_{3}$ to be configurations of the puzzle such that $C_{1}$ $\neq$ $C_{3}$
and let $\varphi_{M_{1}}$ be a function that performs move $M_{1}\in M$, such
that $\varphi_{M_{1}}(C_{1})$ = $\varphi_{M_{1}}(C_{3})$. In order for
$\varphi_{M_{1}}$ to be injective, we must show that $\varphi_{M_{1}}(C_{1})$
= $\varphi_{M_{1}}(C_{3})$ implies $C_{1}$ = $C_{3}$. We know there exists a
$M_{1}^{-1}$ because $M$ is a group, as previously defined.

$\displaystyle\varphi_{M_{1}}(C_{1})=M_{1}*C_{1}=M_{1}*C_{3}=\varphi_{M_{1}}(C_{3})$
(6) $\displaystyle M_{1}^{-1}*M_{1}*C_{1}=M_{1}^{-1}*M_{1}*C_{3}$

Therefore by (6), it must be that $C_{1}$ = $C_{3}$ and so $\varphi_{M_{1}}$
is injective, or one-to-one.

Now that we’ve shown this mapping is well-defined and injective, we can
proceed to proving this operation is surjective, or onto. The proof of
surjectivity is trivial because we have already considered $M$ to be a group,
meaning that there exists some move $M_{1}^{-1}$ that can be performed on
$\varphi_{M_{1}}(C_{1})$ to reach $C_{1}$; and so, $\varphi_{M_{1}}$ is
surjective.
∎

We will now transition into explaining why Lloyd’s puzzle is truly unsolvable.

######  Claim 0.2.

Lloyd’s 15 Tile Puzzle is unsolvable

###### Proof.

If we simplify the moves needed to reach from our starting configuration (Fig
1) to our goal state (Fig 2) can be simply written as $(14\;15)$, since we
just need to swap the 14 and 15 tile.

We can use our implementation of the mapping $M$ (3) to rewrite $(14\;15)$:

$(14\;15)=M_{n}*...*M_{2}*M_{1}$ (7)

Because the empty tile stays in the same location in the goal state as the
ending state, we can claim that tile 16 moves up and down an equal number of
times and left and right an equal number of times. This means that the number
of transpositions on the right side of (7) is even. However, the left side of
our equation has an odd number of transpositions.

This means we have a contradiction in that we are equating an even number of
transpositions (left side of the equation) to an even number of transpositions
(right side of the equation). This contradictory proof shows that it is
impossible to reach the goal state (Fig 2) from our starting state (Fig 1).
∎

As such, we can conclude, through rigorous mathematical backing, that Sam
Lloyd’s variation of the 15 Tile Puzzle was indeed unsolvable. By initially
developing a mappings/permutations of configurations and moves on the 4x4 grid
to eventually grounding the moves as a series of transpositions, we have
developed the mathematical intuition behind Lloyd’s Puzzle and why it is
unsolvable.

Further connections to abstract algebra topics can be found in other variants
of the 15 Tile Puzzle, in which we start with a random configuration and we
can use similar strategies to not only determine if the puzzle is solvable
given that current state, but also develop a unique trajectory or sequence of
moves that will allow us to go from our random configuration (i.e. Fig 3) to
the goal state (Fig 2).

## Acknowledgement

We thank Professor Christopher Ryba and the Mathematics Department of the
University of California, Berkeley for providing the adequate mathematical
background on these topics to allow for this deeper investigation.

## References

  * [1] Trappa, Peter (2004) _Permutations and the 15-Puzzle_ , University of Utah, Mathematics.
  * [2] Conrad, Keith _The 15-Puzzle (And Rubik’S Cube)_ , University of Connecticut, Mathematics.
  * [3] Heinold, Bryan (2012) _Rubik’s Cubes_ , Mount St. Mary’s University, Mathematics.


# Rise and fall of entanglement between two qubits in a non-Markovian bath

Sayan Roy  Theoretische Physik, Universität des Saarlandes, D-66123
Saarbrücken, Germany  Christian Otto  Theoretische Physik, Universität des
Saarlandes, D-66123 Saarbrücken, Germany  Raphaël Menu  Theoretische Physik,
Universität des Saarlandes, D-66123 Saarbrücken, Germany  Giovanna Morigi
Theoretische Physik, Universität des Saarlandes, D-66123 Saarbrücken, Germany

###### Abstract

We analyse the dynamics of quantum correlations between two qubits coupled to
a linear chain of oscillators. The chain mediates interactions between the
qubits and acts as a non-Markovian reservoir. The the model is amenable to an
analytical solution when the initial state of the chain is Gaussian. We study
the dynamics of the qubits concurrence starting from a separable state and
assuming that the chain spectrum is gapped and the chain is initially in a
thermal state. We identify three relevant regimes that depend on the strength
of the qubit-chain coupling in relation to the spectral gap. These are (i) the
weak coupling regime, where the qubits are entangled at the asymptotics; (ii)
the strong coupling regime, where the concurrence can exhibit collapses
followed by revivals with exponentially attenuated amplitude; and (iii) the
thermal damping regime, where the concurrence rapidly vanishes due to the
chain’s thermal excitations. In all cases, if entanglement is generated, this
occurs after a finite time has elapsed. This time scale depends exponentially
on the qubits distance and is determined by the spectral properties of the
chain. Entanglement irreversible decay, on the other hand, is due to the
dissipative effect induced by the coupling with the chain and is controlled by
the coupling strength between the chain and qubits. This study unravels the
basic mechanisms leading to entanglement in a non-Markovian bath and allows to
identify the key resources for realising quantum coherent dynamics of open
systems.

##  I Introduction

The coupling to the surrounding environment is commonly considered the cause
of the fragility of quantum superpositions and entanglement. This fragility is
a central challenge for quantum technologies Preskill (2018), that is usually
addressed by trying to isolate the quantum system and to actively correct the
detrimental effects of the environment on the system’s dynamics. This requires
an accurate knowledge of the noise induced by the external environment. Within
a microscopic theory, the environment is described by a second, large physical
system interacting with the system of interest Gallis and Fleming (1990);
Balian (1991). These interactions establish entanglement between the system
and the environment degrees of freedom, that in turn result in loss of
coherence, namely, decoherence Zurek (2009); Mintert  _et al._ (2005); Konrad
_et al._ (2007). Nevertheless, there are counterexamples where the coupling
with an external bath can even lead to entanglement between the system’s
constituents Plenio  _et al._ (1999); Beige  _et al._ (2000); Braun (2002);
Benatti  _et al._ (2003); Cattaneo  _et al._ (2021); Horovitz and Henkel
(2021). The concepts underlying this environment-induced quantum coherence
have been used to design protocols for quantum state preparation and computing
using dissipation Poyatos  _et al._ (1996); Pielawa  _et al._ (2007); Kraus
_et al._ (2008); Verstraete  _et al._ (2009); Menu  _et al._ (2022). This
progress has put forward the need to systematically understand what are the
features of an environment that are resources, and which are instead
detrimental to quantum coherent dynamics. An important resource are
decoherence free subspaces, namely, subspaces of the system’s Hilbert space
that are effectively decoupled from the environment because of destructive
interference Kempe  _et al._ (2001); Lidar and Birgitta Whaley (2003). This
concept has been recently extended to dynamical symmetries of the master
equation, allowing for the existence of stable limit cycles Buča  _et al._
(2019); Albert (2019). The backflow of information from environment to the
system, characterizing non-Markovianity Breuer  _et al._ (2016); Breuer and
Petruccione (2002), has been shown to be a further important resource Rivas
_et al._ (2010); Benatti  _et al._ (2018). For specific master equations and
settings entanglement generation in a non-Markovian environment can be faster
than for a Markovian one Yu and Eberly (2010); Benatti  _et al._ (2018).
Moreover, non-Markovian baths can distribute entanglement between distant
nodes Plenio and Semião (2005); Campos Venuti  _et al._ (2006); Giampaolo and
Illuminati (2010); Pielawa  _et al._ (2010); Wolf  _et al._ (2011); Fogarty
_et al._ (2013); Kajari  _et al._ (2012); Taketani  _et al._ (2014); Nicacio
and Semião (2016, 2016); Horovitz and Henkel (2021), providing interesting
perspectives for realizing quantum communication in non-unitary channels
Bylicka  _et al._ (2014). In general, classifying the role of the individual
features and understanding their interplay would open the perspective to
design robust and scalable quantum coherent dynamics in noisy environment.

In this work we perform a detailed characterization of a non-Markovian bath in
terms of its capability to establish entanglement between two qubits. The non-
Markovian enviroment is a chain of oscillators, each qubit couples to one of
the oscillators as illustrated in Fig. 1 but does not directly couple to the
other qubit. Previous studies reported that this configuration may support the
onset of entanglement between the qubits Braun (2002); Plenio and Semião
(2005); Wendenbaum  _et al._ (2020), which can even survive at long times
Braun  _et al._ (2001). The model is amenable to an analytical solution of the
propagator Braun  _et al._ (2001); Wendenbaum  _et al._ (2020). This permits
us to shed light into the individual processes of the bath-induced dynamics,
which are otherwise difficult to simulate in a non-Markovian environment Prior
_et al._ (2010). In detail, we provide a systematic analysis of the processes
that lead to the generation of entanglement and of the ones that cause its
decay. We determine the dependence of the corresponding time scales on the
strength of the coupling between qubits and chain, on the properties of the
chain, and on the distance between the qubits.

Figure 1: (color online) Entanglement generation between two qubits induced by
the coupling with a linear chain. The qubits do not mutually interact but are
coupled to one oscillator of a linear chain by means of an optomechanical
interaction. The chain is a non-Markovian bath that mediates interactions and
establishes entanglement between the qubits. The characteristic time scales of
entanglement are determined as a function of the spectral gap of the chain, of
the strength $\gamma$ of the qubit-oscillator coupling, and of the effective
distance $D$ between the qubits, here quantified by the spatial distance
between the chain’s oscillators to which the qubits couple.

This paper is organised as follows. In Sec. II we introduce the model at the
basis of this study and review some of its basic properties. In Sec. III we
analyse the dynamics of correlations between the qubits as a function of
distance. In Sec. IV we characterize the behavior of entanglement as a
function of the coupling strength between qubit and chain. We then determine
the time scales of entanglement generation and decay and their dependence on
the physical parameters, including the spatial distance between the qubits.
The conclusions are drawn in Sec. V. The appendices provide results
complementing the discussion in Sec. III and IV.

##  II Two qubits coupled to a chain of oscillators

The open quantum system is composed by two qubits, here denoted by $a$ and
$b$, that do not mutually interact but couple to the vibrational modes of a
linear chain. The chain is composed of $N$ oscillators, the qubits couple to
the chain at the oscillators labelled by $\ell_{a}$ and $\ell_{b}$,
respectively, according to the Hamiltonian:

$H_{\mathrm{int}}=-\hbar\bar{\gamma}\left(\sigma^{a}_{x}x_{\ell_{a}}+\sigma^{b}_{x}x_{\ell_{b}}\right)\,,$
(1)

where $\sigma^{a}_{x}$ and $\sigma^{b}_{x}$ are the Pauli operators along $x$
for qubits $a$ and $b$ respectively, and $x_{\ell}$ is the (dimensionless)
position operator of the oscillator at site $\ell=1,\ldots,N$ . The parameter
$\bar{\gamma}$ is a positive coupling constant, determining the strength of
qubit-chain coupling. The chain dynamics is governed by the Hamiltonian (with
periodic boundary conditions)

$\displaystyle
H_{\mathrm{c}}=\frac{\hbar\nu}{2}\sum_{\ell}(x_{\ell}^{2}+p_{\ell}^{2})-\frac{\hbar\bar{J}}{2}\sum_{\ell}x_{\ell}x_{\ell+1}\,,$
(2)

where $p_{\ell}$ is the dimensionless canonically conjugated operators to
$x_{\ell}$, while the positive frequencies $\nu$ and $\bar{J}$ respectively
scale the energy of the oscillators and the coupling between neighboring
oscillators. The chain’s stability is warranted by the inequality
$\bar{J}\leq\nu$. The ratio $\bar{J}/\nu$, in particular, controls whether the
chain’s spectrum is gapped ($\bar{J}<\nu$) or gapless ($\bar{J}=\nu$). The
system’s Hamiltonian is here neglected and the total Hamiltonian is thus given
by $H=H_{\mathrm{c}}+H_{\mathrm{int}}$: The Hamiltonian commutes with the
operators $\sigma^{j}_{x}$ and the chain hence acts as a (non-Markovian)
dephasing bath in the eigenbasis of the operators $\sigma^{j}_{x}$.

In the following, we rescale the energy by $\hbar\nu$ and the time by
$\nu^{-1}$. The rescaled oscillator-oscillator and qubit-oscillator couplings
are $J_{0}=\bar{J}/\nu$ and $\gamma=\bar{\gamma}/(\sqrt{2}\nu)$, respectively.
We further introduce the distance $D$ between the qubits in terms of the
distance between the oscillators to which they couple,

$D=|\ell_{b}-\ell_{a}|\,.$

This quantity corresponds to the number of oscillators separating oscillator
$\ell_{a}$ from oscillator $\ell_{b}$. For later convenience, we will use the
definition $\ell_{a}=N-l+1$ and $\ell_{b}=l$, with $l=1,\ldots,\frac{N}{2}$
for $N$ even ($l=1,\ldots,\frac{N+1}{2}$ for $N$ odd). According to this
notation, the distance between the qubits can be rewritten as $D=N+1-2l$.

###  II.1 Symmetries

In order to determine the propagator, we make use of the symmetry of the total
Hamiltonian by reflection about the centers of the chain segments that
connects the two qubits. Given the periodic boundary conditions, there are two
such points and we choose the one separating the shortest segment. This
symmetry defines the separated subspaces containing respectively the symmetric
and antisymmetric states under this reflection. In order to single out this
property we introduce the chain’s symmetric and antisymmetric coordinates

$\displaystyle
x_{N-2\ell+1}^{S}=\frac{1}{\sqrt{2}}\left(x_{N-\ell+1}+x_{\ell}\right)\,,$ (3)
$\displaystyle
x_{N-2\ell+1}^{A}=\frac{1}{\sqrt{2}}\left(x_{N-\ell+1}-x_{\ell}\right)\,,$ (4)

and the respective canonically-conjugated momenta, $p_{N-2\ell+1}^{S}$ and
$p_{N-2\ell+1}^{A}$. The chain Hamiltonian can be decomposed into the sum of
the Hamiltonian for the symmetric and the antisymmetric modes,
$H_{\mathrm{c}}=H_{\mathrm{c}}^{S}+H_{\mathrm{c}}^{A}$, which mutually
commute: $[H_{\mathrm{c}}^{S},H_{\mathrm{c}}^{A}]=0$. Using these coordinates,
the interaction Hamiltonian, Eq. (1), takes the form

$H_{\mathrm{int}}^{(D)}=-\sqrt{2}\gamma(\sigma_{x}^{a}+\sigma_{x}^{b})x_{D}^{S}-\sqrt{2}\gamma(\sigma_{x}^{a}-\sigma_{x}^{b})x^{A}_{D}\,,$
(5)

for $D>0$. For the given Hamiltonian, the operators $\sigma_{x}^{a,b}$ are
constants of motion. Therefore, a convenient qubits basis for analysing the
dynamics is the basis of eigenstates
$\mathinner{|{ij}\rangle}\equiv\mathinner{|{i}\rangle}_{a}\otimes\mathinner{|{j}\rangle}_{b}$,
with
$\sigma_{x}^{a,b}\mathinner{|{\pm}\rangle}_{a,b}=\pm\mathinner{|{\pm}\rangle}_{a,b}$.
In particular, the two-dimensional subspace spanned by the basis vectors
$\\{\mathinner{|{b^{S}}\rangle}\\}:\\{\mathinner{|{++}\rangle}_{x},\mathinner{|{--}\rangle}_{x}\\}$
couples to the symmetric chain, while the antisymmetric chain couples to the
two-dimensional subspace spanned by the basis vectors
$\\{\mathinner{|{b^{A}}\rangle}\\}:\\{\mathinner{|{+-}\rangle}_{x},\mathinner{|{-+}\rangle}_{x}\\}$.
In this basis the coupling does not modify the occupation of the eigenstates,
but affect the evolution of superpositions. The states
$\\{\mathinner{|{b^{S}}\rangle}\\}$ are eigenstates of the operator
$S_{x}^{S}=(\sigma_{x}^{a}+\sigma_{x}^{b})/2$ at the eigenvalues $b^{S}=\pm 1$
and belong to the kernel of operator
$S_{x}^{A}=(\sigma_{x}^{a}-\sigma_{x}^{b})/2$. Vice versa, the states
$\\{\mathinner{|{b^{A}}\rangle}\\}$ are eigenstates of the operator
$S_{x}^{A}$ at the eigenvalues $b^{S}=\pm 1$ and belong to the kernel of
operator $S_{x}^{S}$. This shows that in general no decoherence-free subspace
exists for $D>0$. The case $D=0$ is special. In this case

$H_{\mathrm{int}}^{(0)}=-\gamma(\sigma_{x}^{a}+\sigma_{x}^{b})x_{0}\,,$ (6)

and there is a decoherence-free subspace of the qubits Hilbert space, which
consists of the kernel of operator $\sigma_{x}^{a}+\sigma_{x}^{b}$.

###  II.2 Chain’s normal modes

Before proceeding, it is useful to introduce the normal mode coordinates
$\tilde{x}_{n}^{S}$, $\tilde{p}_{n}^{S}$ and $\tilde{x}_{n}^{A}$,
$\tilde{p}_{n}^{A}$, with normal mode frequency
$\omega_{n}^{j=S,A}=\sqrt{1-J_{0}\cos k_{n}^{j}}$ and wave numbers
$k_{n}\in[-\pi,\pi)$. As we will consider finite chains, one shall distinguish
between an even or odd numbers of oscillators. In detail: for $N$ even, then
the Brillouin zones for the symmetric and antisymmetric modes are
$k_{n}^{S}=2n\pi/N$ and $k_{n}^{A}=k_{n}^{S}+2\pi/N$ with
$n=0,\ldots,\frac{N}{2}-1$. For $N$ odd, then $n=0,\ldots,\frac{N-1}{2}$ for
$k_{n}^{S}$ and $n=0,\ldots,\frac{N-3}{2}$ for $k_{n}^{A}$.

In the normal mode representation, the symmetric and antisymmetric part of the
chain’s Hamiltonian read

$\displaystyle\tilde{H}^{j=S,A}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n}\left(\tilde{p}_{n}^{j2}+\omega_{n}^{j2}\tilde{x}_{n}^{j2}\right)\,.$
(7)

We write the interaction Hamiltonian $H_{i}$ using the chain’s normal modes:

$H_{\mathrm{int}}=-\sum_{n}\left(\tilde{\gamma}_{n}^{S}\tilde{x}_{n}^{S}(\sigma_{x}^{a}+\sigma_{x}^{b})/2+\tilde{\gamma}_{n}^{A}\tilde{x}_{n}^{A}(\sigma_{x}^{a}-\sigma_{x}^{b})/2\right)\,.$
(8)

The coupling constants now depend on the normal modes and read:

$\displaystyle\tilde{\gamma}^{S}_{n}=\frac{2\sqrt{2}\gamma}{\sqrt{N}}\cos\left(\frac{k_{n}^{S}D}{2}\right)\,,$
$\displaystyle\tilde{\gamma}^{A}_{n}=\frac{2\sqrt{2}\gamma}{\sqrt{N}}\sin\left(\frac{k_{n}^{A}D}{2}\right)\,.$
(9)

The coupling to the mode at lowest frequency ($n=0$) is
$\tilde{\gamma}_{0}^{S}=\frac{2\gamma}{\sqrt{N}}$, and the coupling to the
mode at largest frequency ($n=N/2$) takes the form
$\tilde{\gamma}^{A}_{N/2}|_{\text{even}}=\frac{2\gamma}{\sqrt{N}}(-1)^{\frac{D+1}{2}}$
and
$\tilde{\gamma}^{A}_{(N-1)/2}|_{\text{odd}}=\frac{2\gamma}{\sqrt{N}}(-1)^{\frac{D}{2}}$
(the latter distinction is a finite-size effect). Interestingly, the distance
introduces a characteristic wave number $k_{D}=2\pi/D$ which modulates the
coupling strength to a given mode.

We will analyse the dynamics as a function of the coupling strength $\gamma$,
of the distance $D$, and of the elastic strength $J_{0}$. The latter controls
the gap of the chain’s spectrum,

$\omega_{0}=\sqrt{1-J_{0}}\,.$

In this work we choose $J_{0}<1$, thus the spectrum is gapped. The coupling
$J_{0}$ also determines the bandwidth
$\Delta\omega=\sqrt{1+J_{0}}-\sqrt{1-J_{0}}$, as well as the velocity with
which information propagates along the chain. For $\gamma=0$ this is
quantified by the Lieb-Robinson bound, which gives the maximal velocity
$v_{\text{LR}}$ with which information can propagate through a non-
relativistic quantum system Lieb and Robinson (1972):

$v_{\text{LR}}=2\max_{k\in BZ}|\nabla_{k}\,\omega_{k}|,$ (10)

where $BZ$ here indicates the Brillouin zone. Note that $v_{\rm LR}$ vanishes
for $J_{0}\to 0$ and monotonously increases in the interval $J_{0}\in[0,1]$.

###  II.3 Propagator

Let $\rho_{d}(t)$ be the density matrix of the two qubits at time $t$. It is
obtained by tracing out the chain’s degrees of freedom from the density matrix
$\chi(t)$ of the total system, $\rho_{d}(t)=\text{Tr}_{c}[\chi(t)]$, where
$\text{Tr}_{c}$ denotes the trace over the chain’s degrees of freedom. In this
work we assume that at $t=0$ the initial density matrix $\chi(0)$ is a
separable state of qubits and chain, $\chi(0)=\rho_{d}(0)\otimes R_{\beta}$.
The propagator $\Lambda_{t}$, connecting the state of the qubits at time $t$
with the initial state $\rho_{d}(0)$, is defined as

$\rho_{d}(t)=\text{Tr}_{c}[{U}(t)\rho_{d}(0)\otimes
R_{\beta}{U}^{\dagger}(t)]\equiv\Lambda_{t}\rho_{d}(0)\,.$ (11)

Here, $U(t)=\exp(-{\rm i}(H_{\mathrm{int}}+H_{\mathrm{c}})t/\hbar)$ is the
evolution operator. An analytical form for the propagator is found when the
chain is in a Gaussian state. In what follows the chain is initially prepared
in a thermal state with inverse temperature $\beta$, namely,
$R_{\beta}=\exp(-\beta H_{c})/Z$, with $Z$ being the partition function. The
propagator is conveniently expressed in the basis of eigenstates
$\\{\mathinner{|{b_{i}}\rangle}\\}:\\{\mathinner{|{++}\rangle},\mathinner{|{--}\rangle},\mathinner{|{+-}\rangle},\mathinner{|{-+}\rangle}\\}$.
Recalling that the symmetric (antisymmetric) subspace is
$\\{\mathinner{|{b^{S}}\rangle}\\}:\\{\mathinner{|{++}\rangle}_{x},\mathinner{|{--}\rangle}_{x}\\}$
(respectively
$\\{\mathinner{|{b^{A}}\rangle}\\}:\\{\mathinner{|{+-}\rangle}_{x},\mathinner{|{-+}\rangle}_{x}\\}$),
the matrix elements coupling states in the same subspace evolve according to
Braun  _et al._ (2001); Wendenbaum  _et al._ (2020)

$\displaystyle\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(t)\mathinner{|{b_{j}^{\alpha}}\rangle}=\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(0)\mathinner{|{b_{j}^{\alpha}}\rangle}$
(12)
$\displaystyle\times\exp\left(-f^{\alpha}(t)\left(b_{i}^{\alpha}-b_{j}^{\alpha}\right)^{2}\right)\,,$

where $\alpha=S,A$ and

$\displaystyle f^{\alpha}(t)$
$\displaystyle=\frac{1}{2}\sum_{n}\frac{\tilde{\gamma}_{n}^{\alpha\,2}}{\omega_{n}^{\alpha\,3}}\left(2\bar{n}(\omega_{n}^{\alpha})+1\right)\left(1-\cos\left(\omega_{n}^{\alpha}t\right)\right)\,,$
(13)

with the mean thermal occupation of the $\alpha$-th symmetric (antisymmetric)
mode of the chain being given by the Bose-Einstein statistics
$\bar{n}(\omega)=1/(e^{\beta\omega}-1)$. The dynamics leads to damping of the
off-diagonal elements within the same subspace, while the diagonal elements
are constant in time. We will denote this term by damping term (or also
attenuation).

The time evolution of the matrix elements between the symmetric and the
antisymmetric subspaces takes instead the form:

$\displaystyle\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(t)\mathinner{|{b_{j}^{\beta}}\rangle}_{\beta\neq\alpha}=\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(0)\mathinner{|{b_{j}^{\beta}}\rangle}_{\beta\neq\alpha}$
(14)
$\displaystyle\times\exp\left(-f^{\alpha}(t)-f^{\beta}(t)+i(\varphi^{\alpha}(t)-\varphi^{\beta}(t))\right)\,,$

where the time-dependent phases read Braun  _et al._ (2001); Wendenbaum  _et
al._ (2020)

$\displaystyle\varphi^{\alpha}(t)$
$\displaystyle=\frac{1}{2}\sum_{n}\left(\frac{\tilde{\gamma}_{n}^{\alpha\,2}}{\omega_{n}^{\alpha\,2}}t-\frac{\tilde{\gamma}_{n}^{\alpha\,2}}{\omega_{n}^{\alpha\,3}}\sin\left(\omega_{n}^{\alpha}t\right)\,\right)\,.$
(15)

The phase is the sum of two contributions: (i) a contribution linear in time
and (ii) a multi-chromatic, oscillating contribution. The term (i) gives rise
to a periodic oscillation at frequency
$\Omega_{0}=\Omega_{0}^{S}-\Omega_{0}^{A}$ with

$\Omega_{0}^{\alpha}=\sum_{n}\frac{\tilde{\gamma}_{n}^{\alpha\,2}}{\omega_{n}^{\alpha\,2}}\,.$
(16)

The oscillations can be associated to an effective Hamiltonian dynamics
similarly to the Lamb-shift Hamiltonian of quantum electrodynamics Bethe
(1947); Buchheit and Morigi (2016). Due to this analogy and to its collective
nature (see the following section), we will denote this frequency by
collective Lamb shift Gross and Haroche (1982); Scully (2009); Konovalov and
Morigi (2020). For a single qubit the multichromatic oscillations are
responsible for decoherence Braun  _et al._ (2001). For this reason we will
refer to this term as decoherence or also dephasing term.

###  II.4 Discussion

The form of the propagator allows us to identify some relevant time scales. We
first observe that real and imaginary parts exhibit multichromatic
oscillations at the frequency of the chain’s spectrum. At the time $t\ll
t_{\rm max}\equiv 1/\omega_{\rm max}$ one finds the scaling behavior $f\propto
t^{2}$ while $\varphi\propto t^{3}$, which is independent of the chain
spectrum Strunz  _et al._ (2003); Braun  _et al._ (2001). During the
evolution, as time goes by the influence of the chain’s normal modes at
decreasing frequency start to be important. When $\Delta\omega\gg\omega_{0}$,
one can identify a time scale separation at which the low-frequency modes of
the reservoir spectrum become relevant in determining the system evolution
Gualdi and Koch (2013). In this work this is not the case, since the chain
spectrum is gapped, $\omega_{0}>0$ and the largest frequency of the normal
modes, $\omega_{\rm max}=\omega_{0}+\Delta\omega$, is chosen to be of the same
order as $\omega_{0}$. For our dynamics, an important time scale is $t_{\rm
chain}\sim\pi/\Delta\omega$, at which the oscillations at the smallest and at
the largest chain frequencies are out of phase by $\pi$. At this time scale
the decoherence and damping terms start to oscillate. There is a further
important time scale to consider in our analysis. This is the time scale over
which the finite size of the chain becomes important and is essentially the
Poincaré time Balian (1991). We estimate the Poincaré time using the Lieb-
Robinson bound, $t_{P}\sim N/(2v_{\rm LR})$. For $N\gg 1$ and times $t\ll
t_{P}$ the dynamics is independent of the specific chain’s size.

Let us now discuss the real and imaginary parts of the propagator separately.
The dissipative component scales with $\bar{n}(\omega)$ and thus increases
with the initial temperature of the chain. Its value as a function of time is
bound from above: $f\leq 16(2\bar{n}(\omega_{0})+1)\gamma^{2}/\omega_{0}^{3}$.
Hence, for weak couplings and moderate temperatures the off-diagonal elements
of the density matrix are attenuated with respect to the initial value, while
for large couplings they vanish after a time $t\gg t_{\rm max}$. The function
$\varphi$, instead, is independent of the temperature. Interestingly, for a
single qubit it vanishes identically Braun  _et al._ (2001). For two qubits,
this term vanishes when the qubits are either in an eigenstate $|ij\rangle$ of
$\sigma_{x}^{a,b}$ or in a maximally entangled state. This also implies that
the frequency $\Omega_{0}$ is different from zero only in the presence of a
second qubit. Note that the frequency $\Omega_{0}$ corresponds to an
effective, coherent interaction between the qubits that emerges because of the
coupling with the bath.

##  III Bath-induced dynamics

The physical problem here considered is an example of bath-induced dynamics,
where the reservoir is intrinsically non-Markovian. Notably, for a single
qubit the dynamics is solely dissipative, while coherent oscillations and
decoherence are induced only in the presence of a second qubit. In order to
identify the salient physical regimes, it is useful to first consider three
asymptotic, well-defined limits, taking care that the symmetry of the two-
qubit dynamics is preserved. (i) In the trivial limit $\gamma=0$ the qubits
are isolated and undergo no dynamics. (ii) When $\overline{J}\to 0$ and at
finite $\gamma$, the oscillators of the chain become decoupled, and the
dynamics reduces to a qubit coupled to a single oscillator at frequency
$\omega_{0}=1$. In this case, correlations can be established between the
qubits only when they couple to the same oscillator, namely, for $D=0$. (iii)
Finally, in the limit $\gamma\to\infty$ the dynamics reduces to an exchange of
excitations between qubit and the oscillator to which it couples. The
interaction Hamiltonian tends to freeze the oscillator in an eigenstate of the
position operator and to suppress propagation of excitations along the chain.
Outside of these asymptotic cases, the dynamics is characterized by
oscillations at frequency $\Omega_{0}$ that are dephased by the second term in
the function $\varphi$ and damped by the attenuation function.

###  III.1 Weak and strong coupling regimes

Some salient regimes can be identified for a gapped chain spectrum
($\overline{J}<1$), which is the situation we analyse in this work. For
$\gamma\ll\omega_{0}^{3/2}$, the coupling to the chain gives rise to the
coherent interaction $\Omega_{0}$ while decoherence and damping (at moderate
temperatures) are small corrections that limit the visibility of the
oscillations. We will denote this regime by weak-coupling regime. In the
opposite limit, the strong coupling regime $\gamma\gg\omega_{0}^{3/2}$,
damping and decoherence become dominant at sufficiently long times, which we
identify with $t\gg t_{\rm max}$. There is a third regime, where the
temperature is such that $\gamma\ll\omega_{0}^{3/2}$ but
$\gamma\sqrt{\bar{n}}\gtrsim\omega_{0}^{3/2}$: in this case one observes
decaying oscillations due to damping, while dephasing is negligible. The
oscillations vanish at the asymptotics when
$\gamma\sqrt{\bar{n}}\gg\omega_{0}^{3/2}$. We will denote this regime by
thermal-damping regime. These considerations are useful for understanding the
dynamics of observables and correlations.

###  III.2 Dynamics of the correlation functions

For a single qubit the coupling with the chain induces decoherence and
dissipation, but no dynamics Strunz  _et al._ (2003). Interestingly, the
coupling of a second qubit gives rise to observable effects on the first
qubit. In this section we discuss these dynamics by analysing the single-
particle expectation values $\langle\sigma_{\alpha}^{a,b}\rangle$ and the two-
particle correlation function,

$g_{\alpha,\beta}(t,t^{\prime})=\langle\sigma_{\alpha}^{a}(t)\sigma_{\beta}^{b}(t^{\prime})\rangle-\langle\sigma_{\alpha}^{a}(t)\rangle\langle\sigma_{\beta}^{b}(t^{\prime})\rangle\,,$
(17)

where $\langle\cdot\rangle={\rm Tr}\\{\cdot\rho_{d}(0)\\}$. The qubits are
initially aligned along $y$, namely, they are prepared in the pure states
$\mathinner{|{\psi}\rangle}_{j=a,b}=(\mathinner{|{+}\rangle}_{x}+\mathinner{|{-}\rangle}_{x})/\sqrt{2}$.
Moreover, we assume that the chain is at low temperatures, such that the mean
occupations of the normal modes are $\bar{n}(\omega)\ll 1$.

Figure 2: Dynamics of the observable $\langle\sigma_{y}^{a}(t)\rangle$ when
the qubits are at distance $D=0$ (upper row) and at distance $D=10$ (lower
row). The coupling strength is $\gamma=0.04$ in the left subplots, (a) and
(c), and $\gamma=0.3$ in the right subplots (b) and (d). The chain consists of
$N=10^{4}$ spins with $\omega_{0}=0.45$ and inverse temperature
$\beta=10^{5}$. The qubits are initially prepared in the pure state
$\mathinner{|{\psi}\rangle}=\mathinner{|{\psi}\rangle}_{a}\otimes\mathinner{|{\psi}\rangle}_{b}$
with
$\mathinner{|{\psi}\rangle}_{j=a,b}=(\mathinner{|{+}\rangle}_{x}+\mathinner{|{-}\rangle}_{x})/\sqrt{2}$,
corresponding to eigenstates of $\sigma_{y}^{a,b}$ at eigenvalue 1.

Figure 2 shows the dynamics of the $y$-component of one of the qubits in the
weak (left) and in the strong (right) coupling regime, when the second qubit
is at distance $D=0$ (upper row) and $D>0$ (lower row). We first discuss the
case $D=0$, where there is a decoherence-free subspace. In the weak coupling
regime, Fig. 2(a), the qubit undergoes slow oscillations at the frequency
$\Omega_{0}$, whose visibility is close to unity. We emphasize that the qubit
oscillation is a cooperative effect determined by the collective Lamb shift.

The strong coupling regime is shown in Fig. 2(b): The oscillations become
multichromatic, the characteristic frequencies are the normal modes of the
chain which are excited during the corresponding time scale, the corresponding
amplitude is quickly damped to the asymptotic value.

When instead the qubits are at finite distance, here $D=10$, the dynamics at
weak coupling is effectively frozen out, see Fig. 2(c), while at strong
coupling the expectation value first quickly decays to a non-vanishing value,
and then undergoes a slow dynamics. This behavior for relatively large
distances is understood by inspecting the dependence on $D$ of the parameter
$\bar{\gamma}_{n}$, Eq. (II.2), which scales the Fourier components of the
propagator. The damping of the symmetric and antisymmetric modes is multiplied
by the factor ${\cos^{2}(kD/2)}$ and ${\sin^{2}(kD/2)}$, respectively, while
the Fourier components of the imaginary part, $\varphi^{S}-\varphi^{A}$, are
scaled by the factor ${\cos(kD)}$. In particular,

$\Omega_{0}(D)\simeq 4\gamma^{2}/N\sum_{n}\cos(k_{n}D)/\omega_{n}^{2}\,.$ (18)

As $D$ increases the value of $\Omega_{0}$ decreases and for the spectrum
$\omega_{n}$ here considered, then $\Omega_{0}\to 0$ for $D\to\infty$.
Therefore, for $D\gg 1$ the oscillation frequency is such that
$|\Omega_{0}|\ll 1$, and in the weak coupling limit the dynamics is
essentially frozen over the considered time scale. In the strong coupling
limit, Fig. 2(d), the expectation value quickly decreases due to damping,
while the slower dynamics is due to the slow oscillation at frequency
$\Omega_{0}$.

Figure 3 displays the correlations between the two qubits, corresponding to
the subplots of Fig. 2. The correlations are here evaluated at equal time, and
we choose to show the representative case $\alpha=\beta=y$. Starting from a
separable state, we observe that the correlations grow with $t^{3}$ at very
short times. This demonstrates that correlations are first established by the
imaginary component of the propagator. After this transient, in subplots (a),
(c) and (d) they exhibit a slow oscillation at the characteristic frequency
$\Omega_{0}$, about which we observe fast oscillations at the normal mode
frequencies of the chain (which are not visible in the weak coupling regime,
subplot (a), because the amplitude is very small). The strong coupling regime
for $D=0$, Fig. 3(b), shall be discussed apart. In this case damping and
dephasing are dominant and the correlations quickly reach a stationary state
about which fast oscillations occur. We refer the reader to Appendix A for the
analytic expression of all two-point correlations at equal and different
times. These results allow us to identify mechanisms, which determine the
dynamics of entanglement. Interestingly, they also show that the measurement
of single and two-body correlations provide an insightful probe of the
surrounding, non-Markovian environment.

Figure 3:  Qubit-qubit correlations, (17) for $\alpha=\beta=y$ and at equal
time, $g_{yy}(t,t)$. The parameters of the subplots are the same as in Fig. 2.

##  IV Entanglement

Having shown that the chain establishes correlations between the qubits, we
verify when these correlations are non-classical. Several measures have been
discussed in the literature for studying entanglement between qubits Horodecki
_et al._ (2009). For our purpose, the concurrence is a useful quantity, since
for two qubits it also quantifies entanglement. The concurrence is analysed as
a function of time and is defined as Wootters (1998):

$\displaystyle
C(\rho(t))=\max\\{0,\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4}\\},$ (19)

where $\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}$ are the eigenvalues,
in decreasing order, of the matrix
$\omega(t)=\sqrt{\sqrt{\rho_{d}(t)}\tilde{\rho}_{d}(t)\sqrt{\rho_{d}(t)}}$
with
$\tilde{\rho}_{d}(t)=(\sigma_{y}\otimes\sigma_{y})\rho_{d}^{*}(t)(\sigma_{y}\otimes\sigma_{y})$
and the complex conjugate is taken in the eigenbasis of the Pauli matrix
$\sigma_{z}$. The concurrence being a monotonous function of entanglement, it
provides information about the growth or decay of quantum correlations between
the qubits. The definition we use, moreover, applies also to mixed states.

In order to analyse the capability of the chain to mediate entanglement, we
assume that the qubits are initially uncorrelated. Their initial density
matrix is $\rho_{d}(0)=\rho_{a}(0)\otimes\rho_{b}(0)$, with $\rho_{j=a,b}(0)$
the density matrix of qubit $a,b$ and
$\rho_{j}(0)=\mathinner{|{\Psi}\rangle}_{j}\mathinner{\langle{\Psi}|}$ is a
pure state. We parametrize the initial state in the Bloch-sphere
representation

$\mathinner{|{\Psi}\rangle}_{j}=\cos\theta_{j}\mathinner{|{+}\rangle}_{x}+e^{i\phi_{j}}\sin\theta_{j}\mathinner{|{-}\rangle}_{x},$
(20)

with $\theta_{j},\phi_{j}$ real and $\mathinner{|{+}\rangle}_{x}$ and
$\mathinner{|{-}\rangle}_{x}$ the eigenvectors of $\sigma_{x}$ Pauli matrix
with eigenvalues $\pm 1$ respectively.

We first focus on the characteristic time scales of entanglement for the case
in which there exists a decoherence-free subspace for the qubits, namely, the
qubits couple to the same oscillator. We then analyse the characteristic
length scales by studying the situation where the qubits couple to different
oscillators as a function as the distance $D$ between the oscillators.

###  IV.1 Time scales of entanglement

When the qubits couple to the same oscillator, $D=0$, the antisymmetric
subspace is decoupled from the chain and is a decoherence-free subspace. Deep
into the weak-coupling regime, in leading order the dynamics is coherent and
the concurrence is given by $C(\rho)=\sqrt{2(1-\text{Tr}\\{\rho_{R}^{2}\\})}$
where $\rho_{R}$ is the reduced density matrix of one qubit obtained after
tracing out the Hilbert space of the other one Horodecki  _et al._ (2009).
When the qubits are initially prepared in the initial state (20), then the
concurrence takes the form

$C(\rho(t))\simeq|\sin(2\theta_{a})\sin(2\theta_{b})\sin(\Omega_{0}t)|\,,$
(21)

and periodically oscillates at the frequency of the collective Lamb shift
$\Omega_{0}$. The periodic oscillation characterizes also the dynamics of the
correlations in the weak-coupling regime and indicates that the collective
Lamb-shift term of the propagator is responsible for the appearance of
entanglement. The initial state of the qubits determines the maximum reached
by the concurrence: The concurrence varies between 0 and $C_{\rm
max}=|\sin(2\theta_{a})\sin(2\theta_{b})|$ with period $\pi/\Omega_{0}$ and
reaches the maximal amplitude $C_{\rm max}=1$ for initial states with Bloch
angles $\theta_{a},\theta_{b}=(2n+1)\pi/4$, with $n\in\mathbb{Z}$, namely,
when the qubits are aligned along $y$ or $z$.

We note that Eq. (21) is valid after a transient time, and specifically for
$t\gtrsim t_{\rm max}$. For $t\ll t_{\rm max}$, instead, the term
$\Omega_{0}t$ is of the same order as the dephasing term. This is the regime
where decoherence exhibits the universal behavior $\varphi\propto t^{3}$ and
where we do not expect to find entanglement because of dephasing. This
argument shows that there is a finite time scale at which entanglement is
generated between the qubits. We will extensively characterize it in the next
subsection.

We now discuss the behavior of the concurrence for different values of the
coupling strength $\gamma$, ranging from the weak- to the strong-coupling
regime. In what follows we assume that the qubits are initially prepared in
the same eigenstate of $\sigma_{y}$. Unless otherwise stated, the chain is
initially at very low temperatures. Figure 4(a) displays the time evolution of
the concurrence deep in the weak coupling regime, for
$\gamma\ll\omega_{0}^{3/2}$. The concurrence displays a periodic oscillation
as in Eq. (21), with slightly reduced visibility due to damping. There are
also fast oscillations about the slow envelope, that become visible only after
zooming in and are due to the dephasing term of the propagator. For
$\gamma\sim\omega_{0}^{3/2}$, Fig. 4(b), the concurrence exhibits a
multichromatic oscillation whose mean amplitude decreases with $\gamma$. At
even larger values, deep in the strong coupling regime, the concurrence
features collapses and revivals with quickly decreasing amplitude, see Fig.
4(c). The decay of the signal is due to the attenuation, the collapses and
revivals are due to dephasing and rephasing of the function $\varphi$. The
rephasing time is mostly determined by the beating of $\Omega_{0}$ and
$\omega_{0}$. The dephasing time scales with $t_{\rm chain}$ and is determined
by the finite bandwidth $\Delta\omega$ of the chain’s spectrum.

Figure 4: Concurrence $C$ for different values of $\gamma$ when the qubits are
coupled to the same oscillator, $D=0$. Subplot (a) shows $C$ for three
different values of $\gamma$ deep in the weak coupling regime, subplot (b)
shows the concurrence for $\gamma\simeq\omega_{0}^{3/2}$, where the effects of
damping and decoherence starts becoming relevant on the time scale of the
oscillations, and subplot (c) displays the concurrence for
$\gamma\gg\omega_{0}^{3/2}$, in the regime dominated by damping and
decoherence. The qubits are initially prepared in the same pure state (20)
with $\theta_{a}=\theta_{b}=\pi/4$ and $\phi_{a}=\phi_{b}=0$. The chain is
composed by $N=10^{4}$ oscillators with $\omega_{0}=0.45$ and inverse
temperature $\beta=10^{5}$.

In order to study the entanglement at long times we use the time average
concurrence $\bar{C}$, which we define as

$\bar{C}=\frac{1}{t_{\rm end}}\int_{0}^{t_{\rm end}}C(t^{\prime})dt^{\prime},$
(22)

where $t_{\rm end}$ is the integration time. In order to capture a
sufficiently large number of oscillations, the dynamics is evolved over times
$t_{\rm end}$ which are generally longer than the Poincaré time $t_{P}$. We
show in Appendix B that the behavior found at $t_{\rm end}$ is in qualitative
agreement with the one found for $t<t_{P}$. The results we are going to
discuss, hence, give a reliable indication of the be havior in the
thermodynamic limit at finite $\gamma$. The behavior in the region at
$\gamma\to 0$, instead, is determined by finite-size effects, as we will argue
below.

Figure 5(a) displays the entanglement phase diagram as a function of the
coupling strength $\gamma$ and of the characteristic frequency $\omega_{0}$ of
the chain. The diagram qualitatively correspond to the steady state ($t_{\rm
end}\to\infty$) with the exception of the small stripe at $\gamma\to 0$. In
this region the average concurrence vanishes because the time scale, at which
entanglement is generated, is larger than the integration time $t_{\rm end}$
(In Appendix B we show that the size of this region shrinks as $t_{\rm end}$
increases). Outside of this region, the concurrence is maximum deep in the
weak coupling regime, while it decays as $\gamma$ is increased. More
specifically, the contour lines correspond to good approximation to the
constant values $\gamma^{2}/\omega_{0}^{3}$: The average concurrence
monotonically decreases to zero as the ratio $\gamma^{2}/\omega_{0}^{3}$
grows, and with it the role of dissipation and dephasing on the dynamics.

Figure 5(b) shows the effect of the initial chain temperature, here given by
the mean occupation $n_{0}\equiv n(\omega_{0})$, on the asymptotic behavior of
entanglement. At low temperatures, $n_{0}<1$, the average concurrence is
almost independent of the temperature: it decays to zero as $\gamma$ reaches
the strong coupling regime. By increasing $n_{0}$ above unity the system
enters the thermal damping regime: entanglement vanishes at lower values of
$\gamma$. The contour lines, in particular, follow the functional behavior
$\gamma\sqrt{n_{0}}=$ constant, and show that damping is majorly responsible
for the disappearance of the average concurrence.

Figure 5: (color online) Color plot of the average concurrence $\bar{C}$, Eq.
(22) for $D=0$. Subplot (a) shows $\bar{C}$ as a function of $\gamma$ and
$\omega_{0}$ for $\beta=10^{5}$. Subplot (b) displays $\bar{C}$ as a function
of $\gamma$ and $n_{0}\equiv n(\omega_{0})$, the average occupation of the
mode at frequency $\omega_{0}$, for $\omega_{0}=0.45$. The plots have been
evaluated for a chain with $N=10^{4}$ oscillators and integrating over the
time $t_{\rm end}=10^{6}$. The qubits are initially in the same initial state
with Bloch angles $\theta_{a}=\theta_{b}=\pi/4$ and $\phi_{a}=\phi_{b}=0$.

###  IV.2 Entanglement as a function of the distance

Figure 6: Concurrence as a function of time (a) for different values of
$\gamma$ when the qubits are at distance $D=10$ and (b) as a function of $D$
for $\gamma=0.04$. The chain is composed by $N=10^{4}$ oscillators with
$\omega_{0}=0.45$ and $\beta=10^{5}$. The qubits are initially in a separable,
pure state with Bloch angles $\theta_{a}=\theta_{b}=\pi/4$ and
$\phi_{a}=\phi_{b}=0$. It can be observed that there is a finite time
$t_{\text{gen}}$ after which entanglement is generated.

We now analyse the features of entanglement as a function of the distance $D$.
Figure 6(a) displays the time evolution of the concurrence for $D=10$ and for
two different values of $\gamma$. Because of the relatively large distance,
the effective coupling strength of the function $\varphi$ is reduced, so that
the dynamics determined by $\varphi$ is the one corresponding to the weak
coupling regime: the concurrence oscillates at the frequency $\Omega_{0}$. Due
to the different scaling with $D$ the effect of the attenuation becomes
relevant at $\gamma=0.25$ and significantly reduces the maximal value that the
concurrence can attain. At larger values of $\gamma$ we do not observe
entanglement.

Whenever we find entanglement, we observe that at short time the concurrence
starts to grow after a finite time $t_{\rm gen}$ has elapsed, while for
$t<t_{\rm gen}$ there is no entanglement between the qubits. Figure 6(b)
compares the concurrence at relatively short times and for different distances
$D$, showing that $t_{\rm gen}$ increases monotonically with $D$. This result
suggests that entaglement propagates at a finite velocity.

Figure 7: Relation between the distance $D$ and the delay time $t_{\rm gen}$
after which the concurrence grows. $t_{\rm gen}$ is extracted from the
numerical data and corresponds to the time for which the concurrence starts to
grow continuously. The data have been evaluated for $\gamma=0.04$ (orange
triangle) and $\gamma=0.1$ (blue dot). The other parameters are the same as in
Fig. 6. The red stars report the values predicted by the estimate of Eq. (23).
The green diamonds give the propagation time of the light cone with the Lieb-
Robinson bound velocity $v_{LR}$ and are plotted for comparison. The threshold
was set such that $C(t_{\rm gen})=0.001$, corresponding to the maximal value
that the concurrence reaches in the universal regime, where the decoherence
function in the propagator scales with $t^{3}$. In this regime, in particular,
$C$ has no monotonous behavior but fluctuates between 0 and 0.001.

Figure 7 shows $t_{\rm gen}$ as a function of $D$. The time scale of
entanglement generation is relatively independent of $\gamma$ and scales
exponentially with the distance. We compare $t_{\rm gen}$ with the time scale
$t_{LR}=D/v_{LR}$ that an excitation needs for covering the distance $D$ at
the velocity of the Lieb-Robinson bound, Eq. (10). The comparison shows that
$t_{\rm gen}$ is consistent with the Lieb-Robinson bound, but cannot be
related to cone-light propagation across the chain. The dependence of $t_{\rm
gen}$ on $D$ can be understood by considering that entanglement is generated
by the collective Lamb shift. Therefore, it is generated at times where the
coherent term in the imaginary part of the propagator, $\Omega_{0}(D)t$,
exceeds the dephasing term. For short times dephasing is dominant while for
$t\gtrsim t_{\rm max}$ the dephasing term starts to oscillate with the upper
bound $4\gamma^{2}/\omega_{0}^{3}$. We identify the time scale for
entanglement generation using an equation that overestimates its value,
$\Omega_{0}t_{\rm gen}=4\gamma^{2}/\omega_{0}^{3}$. It gives

$\displaystyle t^{\rm
est}_{\mathrm{gen}}\sim\frac{2\pi/\omega_{0}^{3}}{\int_{0}^{2\pi}{\rm
d}k\frac{\cos(kD)}{\omega(k)^{2}}}\,,$ (23)

where we have taken the continuum limit of Eq. (18),
$\Omega_{0}(D)\simeq\frac{4\gamma^{2}}{2\pi}\int_{0}^{2\pi}{\rm
d}k\frac{\cos(kD)}{\omega(k)^{2}}$. The resulting expression is independent of
$\gamma$ and is reported in Fig. 7: $t^{\rm est}_{\mathrm{gen}}$ overestimates
the time extracted from the numerical data but has the same functional
dependence on the distance $D$. These considerations also clarify why we do
not observe entanglement deep in the strong coupling regime for $D>0$. In this
case, in fact, the density matrix decays to a statistical mixture at a faster
time scale than $t_{\rm gen}$, or, in other words, the lifetime of
entanglement is shorter than the time scale at which it can be generated.

Figure 8(a) displays the average concurrence as a function of the coupling
strength $\gamma$ and of the frequency $\omega_{0}$ when the qubits are at
distance $D=1$. By comparing with Fig. 5(a) for $D=0$, we observe the same
qualitative behavior with two salient differences: (i) for the same value of
$\gamma$ and $\omega_{0}$ the concurrence is generally smaller, which we
attribute to the fact that $D>0$ leads to an effectively reduced coupling
strength. Moreover, (ii) the concurrence decreases to zero for
$\omega_{0}\gtrsim 0.9$. This latter behavior is due to the fact that the band
width becomes very small, $\Delta\omega\lesssim 0.1$, the band becomes
flatter, and the collective Lamb shift, that is essential for generating
entanglement, vanishes as
$\Omega_{0}\sim\int_{0}^{2\pi}{\mathrm{d}k\cos(2kD)}=0$.

The average concurrence as a function of $D$ and $\gamma$ is shown in Fig.
8(b). Deep in the weak coupling regime we observe that entanglement is
generated when $\gamma$ exceeds a threshold value that depends on $D$, and
more specifically, it increases monotonically with the distance. This behavior
is a consequence of the scaling of $t_{\rm gen}$ with the distance $D$: In the
region where the concurrence vanishes, $t_{\rm gen}$ is larger than the total
evolution time (Simulations performed for smaller $t_{\rm end}$ show, in fact,
that this threshold moves towards larger values of $\gamma$, see Appendix B).
Entanglement vanishes again at large $\gamma$, when the coupling strength
reaches the strong coupling regime and damping suppresses quantum coherence.

Figure 8: (color online) Color plot of the average Concurrence $\bar{C}$ (a)
as a function of $\gamma$ and $\omega_{0}$ when the qubit are at distance
$D=1$ and (b) as a function of the distance $D$ and of $\gamma$ for
$J_{0}=0.8$ ($\omega_{0}=0.45$). The chain is composed by $N=10^{4}$ and is
initially at inverse temperature $\beta=10^{5}$. The two qubits are in a
separable, pure state, Eq. (20) with Bloch angles
$\theta_{a}=\theta_{b}=\pi/4$ and $\phi_{a}=\phi_{b}=0$. The concurrence was
averaged over $t_{\rm end}=10^{6}$.

###  IV.3 Discussion

The analysis so far is based on the properties of the propagator under the
assumption that the symmetric and the antisymmetric subspaces are decoupled.
We discuss now the case when this property, emerging from a symmetry of the
Hamiltonian, is no longer fulfilled. We now assume that the qubits couple with
different coupling strengths, denoted by $\gamma_{a}$ and $\gamma_{b}$, to the
chain oscillator(s). In this case the propagator for the symmetric and
antisymmetric subspaces, separately, is the same as for the symmetric case,
with now $\gamma=(\gamma_{a}+\gamma_{b})/2$. Additionally, the two subspaces
are now coupled with one another with strength
$\Delta\gamma=(\gamma_{a}-\gamma_{b})/2$. In order to assess the effect of the
asymmetry, we compare different cases at constant $\gamma$ and varying
$\Delta\gamma$.
In Appendix C we provide some details of our study. The results show that,
when the qubits couple to the same oscillator, then the analysis for
$\Delta\gamma=0$ essentially applies also for $\Delta\gamma>0$ as long as
$\Delta\gamma\ll\gamma$. In detail, in the weak coupling regime the
entanglement is reduced by an amount proportional to $\Delta\gamma/\gamma$,
while in the strong coupling regime the effects of the asymmetry
$\Delta\gamma$ are irrelevant, since the corresponding time scale is longer
than the lifetime of entanglement. When the qubits couple to different
oscillators, the asymmetry can be neglected as long as the corresponding time
scale is longer than the time scale $t_{\rm gen}$ at which entanglement is
generated, and which depends solely on the chain spectral properties.

##  V Conclusions

In this work we have analysed the dynamics of entanglement that is generated
between two qubits coupled to a non-Markovian bath, here modelled by a chain
of oscillators with a gapped spectrum. Using the exact solution for the qubits
propagator we could unravel the processes that lead to entanglement between
the two qubits and to its decay. Entanglement is generated for a certain class
of initial, separable states of the qbits by coherent, Hamiltonian processes,
which are reminiscent of the collective Lamb shift of dipolar systems.
Dephasing mechanism lead to collapse and revival of the concurrence, while
damping tends to suppress entanglement. Their interplay is controlled by the
coupling strength of the qubits to the chain and by the elastic constant
coupling the oscillators: Large coupling strengths tend to suppress
entanglement, while instead a large elastic constant tends to favour it. We
emphasize that, in all situations here discussed, the environment is initially
in a thermal state. Its capability to generate correlations requires that the
initial temperature is below an upper bound, which our model allows to
determine.

Interestingly, entanglement is generated after a finite time has been elapsed.
This behavior is due to the interplay between the effective coherent dynamics
and the dephasing mechanism, the corresponding time scale is determined by the
spectral properties of the chain’s normal modes: entanglement is generated on
time scales where the propagator does not exhibit the characteristic universal
scaling with time. The scaling of the velocity, with which the qubits become
entangled, is exponential with the qubit distance. This is consistent with the
Lieb-Robinson bound, but cannot be related to light-cone propagation in the
chain. This behaviour seems at odds with the picture developed by Calabrese
and Cardy on the post-quench dynamics of entanglement as a consequence of
ballistically propagating quasi-particles Calabrese and Cardy (2009). Indeed,
entanglement is generated by the effect of the collective Lamb shift,
introducing an additional characteristic frequency scale in the problem.
Moreover, the geometry of the coupling between chain and qubit introduces a
characteristic wavelength scaling with the qubits distance $D$ that has a
similar effect as the free spectral range of a Fabry-Perot resonator Berman
(1994); Wolf  _et al._ (2011).

Our findings provide important guidelines for designing quantum steering
protocols Roy  _et al._ (2020) which go beyond the Markovian paradigm of
projective measurements Puente  _et al._ . The configuration here considered
could be implemented in quantum optics experiment: the environment would be
mimicked by coupled microcavities Tanese  _et al._ (2013); Wachter  _et al._
(2019), by an optomechanical array Eichenfield  _et al._ (2009); Ludwig and
Marquardt (2013), or by the transverse modes of an ion chain Porras and Cirac
(2004); Serafini  _et al._ (2009). For the dynamics here discussed the qubits
frequencies shall be the smallest frequency scale of the problem. Entanglement
between the qubits could be detected using the protocol discussed in Ref.
Serafini  _et al._ (2009); Taketani  _et al._ (2014) or the measures discussed
in Ref. Cattaneo  _et al._ (2021).

The dynamics of entanglement shall change substantially for a gapless
spectrum, where the weight of damping and dephasing terms in general increases
with time. Future works shall address how the velocity of entanglement
propagation behaves in the presence of disorder Burrell and Osborne (2007),
and how it is modified when the oscillators interact with long-range
interactions, where the bound on the propagation of information follows a
different scaling Tran  _et al._ (2020).

To conclude, we have analysed the entanglement generated between two qubits by
a non-Markovian quantum channel and characterised the mechanisms that
determine its onset, stability and/or decay. Our analysis can be extended to
determine quantum information scrambling in the chain, by considering a qubit
chain where information is mediated by the non-Markovian environment and the
entanglement entropy in order to analyse the dynamics of correlations Harrow
_et al._ (2021). .

## Acknowledgments

The authors are grateful to E. Kajari, E. King, D. Karevski, F. Marquardt, B.
Taketani, and P. Wendenbaum for discussions. This work was funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Project-ID
429529648 TRR 306 QuCoLiMa(”Quantum Cooperativity of Light and Matter”) and by
the Bundesministerium für Bildung und Forschung (BMBF, Project ”NiQ: Noise in
Quantum Algorithms”).

##  Appendix A Correlation functions

Below we report the expectation value of the Pauli operators and of the
correlation functions at equal time and at different times. These can be
straightforwardly evaluated using the propagator. We compute the correlations
with the reduced density matrix in the tensor product basis
$\\{\mathinner{|{++}\rangle}_{x},\mathinner{|{+-}\rangle}_{x},\mathinner{|{-+}\rangle}_{x},\mathinner{|{--}\rangle}_{x}\\}$
(labelled from $1$ to $4$ respectively for notation purpose).

The time evolution for the expectation values of the Pauli matrices are:

$\displaystyle\langle\sigma_{x}^{a}\rangle$ $\displaystyle=$
$\displaystyle\rho_{11}(0)+\rho_{22}(0)-\rho_{33}(0)-\rho_{44}(0)$
$\displaystyle\langle\sigma_{y}^{a}\rangle$ $\displaystyle=$ $\displaystyle
2\,\text{Re}[\rho_{13}(0)e^{i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle+$ $\displaystyle
2\,\text{Re}[\rho_{24}(0)e^{-i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle\langle\sigma_{z}^{a}\rangle$ $\displaystyle=$ $\displaystyle
2\,\text{Im}[\rho_{13}(0)e^{i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
(24) $\displaystyle+$ $\displaystyle
2\,\text{Im}[\rho_{24}(0)e^{-i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$

The correlation functions at equal time read:

$\displaystyle\langle\sigma_{x}^{a}\sigma_{x}^{b}\rangle$ $\displaystyle=$
$\displaystyle\rho_{11}(0)-\rho_{22}(0)-\rho_{33}(0)+\rho_{44}(0)$
$\displaystyle\langle\sigma_{y}^{a}\sigma_{y}^{b}\rangle$ $\displaystyle=$
$\displaystyle
2\,\text{Re}[\rho_{14}(0)]\,e^{-4f^{S}(t)}+2\,\text{Re}[\rho_{23}(0)]\,e^{-4f^{A}(t)}$
$\displaystyle\langle\sigma_{z}^{a}\sigma_{z}^{b}\rangle$ $\displaystyle=$
$\displaystyle-2\,\text{Re}[\rho_{14}(0)]\,e^{-4f^{S}(t)}+2\,\text{Re}[\rho_{23}(0)]\,e^{-4f^{A}(t)}$
$\displaystyle\langle\sigma_{x}^{a}\sigma_{y}^{b}\rangle$ $\displaystyle=$
$\displaystyle
2\,\text{Re}[\rho_{12}(0)e^{i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle-$ $\displaystyle
2\,\text{Re}[\rho_{34}(0)e^{-i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle\langle\sigma_{x}^{a}\sigma_{z}^{b}\rangle$ $\displaystyle=$
$\displaystyle
2\,\text{Im}[\rho_{12}(0)e^{i(\phi^{S}(t)-\phi^{A}(t))}]\,e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle-$ $\displaystyle
2\,\text{Im}[\rho_{34}(0)e^{-i(\phi^{S}(t)-\phi^{A}(t))}]e^{-f^{S}(t)-f^{A}(t)}$
$\displaystyle\langle\sigma_{y}^{a}\sigma_{z}^{b}\rangle$ $\displaystyle=$
$\displaystyle
2\,\text{Im}[\rho_{14}(0)]\,e^{-4f^{S}(t)}-2\,\text{Im}[\rho_{23}(0)]\,e^{-4f^{A}(t)}$

Figure 9: (color online) Same as Fig. 5(a) but for $t_{\rm end}=10^{3}$.
Subplot (b) shows a zoom into the region at small $\gamma$. Subplot (c) shows
the zoom in the same region for the propagation time of Fig. 5(a).

The correlation function at different times can be computed by first evolving
the density matrix $(\,\rho_{t}=U(t)\rho_{0}U^{\dagger}(t)\,)$ for time $t$,
and then modifying it accordingly. We denote the modified density matrix as
$\tilde{\rho}^{\alpha}(t)=\rho_{t}(\,\sigma_{\alpha}^{a}\otimes\mathds{1}\,)$
and let $t^{\prime}=t+\tau$, then the correlation functions have the form:

$\displaystyle\langle\sigma_{x}^{a}(t)\sigma_{x}^{b}(t^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\tilde{\rho}^{x}_{11}(t)-\tilde{\rho}^{x}_{22}(t)+\tilde{\rho}^{x}_{33}(t)-\tilde{\rho}^{x}_{44}(t)$
$\displaystyle\langle\sigma_{y}^{a}(t)\sigma_{y}^{b}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\,\text{Re}[\tilde{\rho}^{y}_{12}(t)e^{i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle+$ $\displaystyle
2\text{Re}[\tilde{\rho}^{y}_{34}(t)e^{-i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle\langle\sigma_{x}^{a}(t)\sigma_{y}^{b}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\,\text{Re}[\tilde{\rho}^{x}_{12}(t)e^{i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle+$ $\displaystyle
2\text{Re}[\tilde{\rho}^{x}_{34}(t)e^{-i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle\langle\sigma_{z}^{a}(t)\sigma_{z}^{b}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\,\text{Im}[\tilde{\rho}^{z}_{12}(t)e^{i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle+$ $\displaystyle
2\text{Im}[\tilde{\rho}^{z}_{34}(t)e^{-i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle\langle\sigma_{x}^{a}(t)\sigma_{z}^{b}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\,\text{Im}[\tilde{\rho}^{x}_{12}(t)e^{i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle+$ $\displaystyle
2\text{Im}[\tilde{\rho}^{x}_{34}(t)e^{-i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle\langle\sigma_{y}^{a}(t)\sigma_{z}^{b}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\,\text{Im}[\tilde{\rho}^{y}_{12}(t)e^{i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$
$\displaystyle+$ $\displaystyle
2\text{Im}[\tilde{\rho}^{y}_{34}(t)e^{-i(\phi^{S}(\tau)-\phi^{A}(\tau))}]e^{-f^{S}(\tau)-f^{A}(\tau)}$

##  Appendix B Average concurrence for $t_{\rm end}=10^{3}$

We report the average concurrence for $t_{\rm end}=10^{3}$, which for the
considered chain size is below the Poincaré time. Figure 9(a) shall be
compared with Fig. 5(a): apart for larger fluctuations, the qualitative
behavior is the same. Subplot (b) and (c) show the zoom in the region at small
$\gamma$ for $t_{\rm end}=10^{3}$ and $t_{\rm end}=10^{6}$, respectively: For
longer integration times entanglement is generated at smaller values of
$\gamma$.

Figure 10 is calculated for the same parameter of Fig. 8(b) but integration
time $t_{\rm end}=10^{3}$: entanglement is established till distance $D\sim 6$
due to the exponential scaling with the distance of the velocity with which it
propagates along the chain.

Figure 10: (color online) Same as Fig. 8(b) but for $t_{\rm end}=10^{3}$.

##  Appendix C Effect of Asymmetric Coupling Strength in the Interaction
Hamiltonian

In the main text we have used the uniform coupling strength for the two qubits
in the interaction Hamiltonian (1), however one may also ask how the results
modify for the case of asymmetric coupling strength. In order to analyze the
role of asymmetry, we introduce the new interaction Hamiltonian

$H_{\text{int}}=-\hbar\left(\bar{\gamma}_{a}\sigma_{x}^{a}x_{l_{a}}+\bar{\gamma}_{b}\sigma_{x}^{b}x_{l_{b}}\right).$
(27)

Following the same treatment as in the main text, the interaction Hamiltonian
can be written in the normal coordinates as

$\displaystyle H_{\text{int}}=-\sum_{n}\Big{(}$
$\displaystyle\frac{\tilde{\gamma}^{S}_{a,n}+\tilde{\gamma}^{S}_{b,n}}{2}\,\tilde{x}_{n}^{S}(\sigma_{x}^{a}+\sigma_{b}^{x})/2$
(28) $\displaystyle+$
$\displaystyle\frac{\tilde{\gamma}^{A}_{a,n}+\tilde{\gamma}^{A}_{b,n}}{2}\,\tilde{x}_{n}^{A}(\sigma_{x}^{a}-\sigma_{b}^{x})/2$
$\displaystyle+$
$\displaystyle\frac{\tilde{\gamma}^{A}_{a,n}-\tilde{\gamma}^{A}_{b,n}}{2}\,\tilde{x}_{n}^{A}(\sigma_{x}^{a}+\sigma_{b}^{x})/2$
$\displaystyle+$
$\displaystyle\frac{\tilde{\gamma}^{S}_{a,n}-\tilde{\gamma}^{S}_{b,n}}{2}\,\tilde{x}_{n}^{S}(\sigma_{x}^{a}-\sigma_{b}^{x})/2\Big{)},$

where the coupling constant reads as (II.2) with $\gamma$ replaced with
$\gamma_{j}$, with $j=a,b$. One can easily recover the case of same coupling
strength (8) by using the condition $\gamma_{a}=\gamma_{b}=\gamma$. In the
interaction Hamiltonian, the first two terms is analogous to the case of
uniform coupling strength. For understanding the role of asymmetry, the last
two terms containing the difference between coupling strength have to be
analyzed. For the same we keep the sums
$\gamma_{+}^{S,A}\equiv\tilde{\gamma}^{S,A}_{a,n}+\tilde{\gamma}^{S,A}_{b,n}$
constant and vary the differences
$\gamma_{-}^{S,A}\equiv\tilde{\gamma}^{S,A}_{a,n}-\tilde{\gamma}^{S,A}_{b,n}$.
Following the same recipe as discussed in Braun  _et al._ (2001); Wendenbaum
_et al._ (2020), the reduced density matrix elements in the same subspace
evolves as

$\displaystyle\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(t)\mathinner{|{b_{j}^{\alpha}}\rangle}=\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(0)\mathinner{|{b_{j}^{\alpha}}\rangle}$
(29)
$\displaystyle\times\exp\left(-[f_{+}^{\alpha}(t)+f_{-}^{\beta}(t)]\left(b_{i}^{\alpha}-b_{j}^{\alpha}\right)^{2}\right)\,,$

where $\alpha,\beta=S,A$, $\alpha\neq\beta$ and

$\displaystyle f_{\pm}^{\alpha}(t)$
$\displaystyle=\frac{1}{2}\sum_{n}\frac{\tilde{\gamma}_{n,\pm}^{\alpha\,2}}{\omega_{n}^{\alpha\,3}}\left(2\bar{n}(\omega_{n}^{\alpha})+1\right)\left(1-\cos\left(\omega_{n}^{\alpha}t\right)\right)\,,$
(30)

with the mean thermal occupation of the $\alpha$-th symmetric (antisymmetric)
mode of the chain being given by the Bose-Einstein statistics
$\bar{n}(\omega)=1/(e^{\beta\omega}-1)$. The dynamics leads to additional
damping (with respect to uniform coupling strength) of the off-diagonal
elements within the same subspace which scales linearly with the differences
$\gamma_{-}^{S,A}$ , while the diagonal elements are still constant in time.

Figure 11: Concurrence $C$ for different values of difference in coupling
strength $\Delta\gamma$ while keeping the sum $\gamma_{a}+\gamma_{b})/2$ fixed
for the case when qubits are coupled to the same oscillator. Subplot (a) shows
the dynamics of concurrence in the weak coupling regime,
$(\gamma_{a}+\gamma_{b})/2\ll\omega_{0}^{3/2}$ . Here we chose
$(\gamma_{a}+\gamma_{b})/2=0.04$. The plot for the concurrence dynamics in
strong coupling regime ($(\gamma_{a}+\gamma_{b})/2\gg\omega_{0}^{3/2}$) is
shown in subplot (b) with $(\gamma_{a}+\gamma_{b})/2=0.6$. The qubits are
initially prepared in the same pure state (20) with
$\theta_{a}=\theta_{b}=\pi/4$ and $\phi_{a}=\phi_{b}=0$. The chain is composed
by $N=10^{4}$ oscillators with $\omega_{0}=0.45$ and inverse temperature
$\beta=10^{5}$.

The time evolution of the matrix elements between the symmetric and the
antisymmetric subspaces takes instead the form:

$\displaystyle\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(t)\mathinner{|{b_{j}^{\beta}}\rangle}_{\beta\neq\alpha}$
$\displaystyle=$
$\displaystyle\mathinner{\langle{b_{i}^{\alpha}}|}\rho_{d}(0)\mathinner{|{b_{j}^{\beta}}\rangle}_{\beta\neq\alpha}$
(31) $\displaystyle\times\exp\Bigg{(}$ $\displaystyle-$ $\displaystyle
f_{+}^{\alpha}(t)-f_{-}^{\alpha}(t)-f_{+}^{\beta}(t)-f_{-}^{\beta}(t)$
$\displaystyle-$ $\displaystyle(-1)^{i+j}[F^{\alpha}+F^{\beta}]$
$\displaystyle+$ $\displaystyle
i(\varphi_{+}^{\alpha}(t)+\varphi_{-}^{\alpha}(t)-\varphi_{+}^{\beta}(t))-\varphi_{-}^{\beta}(t)\Bigg{)}\,,$

where the time-dependent phases read

$\displaystyle\varphi_{\pm}^{\alpha}(t)$
$\displaystyle=\frac{1}{2}\sum_{n}\left(\frac{\tilde{\gamma}_{n,\pm}^{\alpha\,2}}{\omega_{n}^{\alpha\,2}}t-\frac{\tilde{\gamma}_{n,\pm}^{\alpha\,2}}{\omega_{n}^{\alpha\,3}}\sin\left(\omega_{n}^{\alpha}t\right)\,\right)\,.$
(32)

and the additional decay term

$\displaystyle F^{\alpha}(t)$
$\displaystyle=\frac{1}{2}\sum_{n}\frac{\tilde{\gamma}_{n,+}^{\alpha}\tilde{\gamma}_{n,-}^{\alpha}}{\omega_{n}^{\alpha\,3}}\left(2\bar{n}(\omega_{n}^{\alpha})+1\right)\left(1-\cos\left(\omega_{n}^{\alpha}t\right)\right)\,,$
(33)

Figure 12: Concurrence as a function of time for different distances with
constant different coupling strength. Here, $\gamma_{a}=0.05$,
$\gamma_{b}=0.03$. The qubits are initially prepared in the same pure state
(20) with $\theta_{a}=\theta_{b}=\pi/4$ and $\phi_{a}=\phi_{b}=0$. The chain
is composed by $N=10^{4}$ oscillators with $\omega_{0}=0.45$ and inverse
temperature $\beta=10^{5}$.

The same site dynamics can be recovered by putting all the antisymmetric
coupling vector to 0. In contrast to the uniform strength case, here the
antisymmetric subspace is not protected from decoherence due to difference in
coupling strength. The damping in the antisymmetric subspace decays linearly
with the difference.

One can identify the same coupling regimes as before with a modified coupling
strength. The weak coupling regime is when
$(\gamma_{a}+\gamma_{b})/2\ll\omega_{0}^{3/2}$. This regime is mainly governed
by coherent oscillations with negligible decoherence and damping. The strong
coupling regime is the opposite limit when both decoherence and damping
effects become dominant that is
$(\gamma_{a}+\gamma_{b})/2\gg\omega_{0}^{3/2}$. In this framework, we plot
below the dynamics of entanglement by keeping the sum
$(\gamma_{a}+\gamma_{b})/2$ fixed and vary the difference
$\Delta\gamma=|\gamma_{a}-\gamma_{b}|/2$.

For the case of $D=0$, we observe in Fig. 11(a), that in the weak coupling
regime the difference in coupling strength benefits the concurrence to grow
quicker and this is due to the additional frequency in the collective Lamb
shift $\Omega_{0}$ which scales with $\Delta\gamma$. For the case of strong
coupling regime Fig. 11(b) , we observe additional decay of concurrence
maxima, but it retains the collapse and revival type of dynamics.

As a function of distance, the entanglement generation time have the same
behavior even with the difference in coupling strength which is clearly
visible in Fig. 12.

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Gallis and Fleming (1990) M. R. Gallis and G. N. Fleming, Phys. Rev. A 42, 38 (1990).
  * Balian (1991) R. Balian,  _From Microphysics to Macrophysics: Methods and Applications of Statistical Physics_, From Microphysics to Macrophysics: Methods and Applications of Statistical Physics No. v. 2 (Springer-Verlag, 1991).
  * Zurek (2009) W. H. Zurek, Nature Physics 5, 181 (2009).
  * Mintert  _et al._ (2005) F. Mintert, A. R. Carvalho, M. Kuś, and A. Buchleitner, Physics Reports 415, 207 (2005).
  * Konrad  _et al._ (2007) T. Konrad, F. de Melo, M. Tiersch, C. Kasztelan, A. Aragão, and A. Buchleitner, Nat. Phys. 4, 99 (2007).
  * Plenio  _et al._ (1999) M. B. Plenio, S. F. Huelga, A. Beige, and P. L. Knight, Phys. Rev. A 59, 2468 (1999).
  * Beige  _et al._ (2000) A. Beige, D. Braun, B. Tregenna, and P. L. Knight, Phys. Rev. Lett. 85, 1762 (2000).
  * Braun (2002) D. Braun, Phys. Rev. Lett. 89, 277901 (2002).
  * Benatti  _et al._ (2003) F. Benatti, R. Floreanini, and M. Piani, Phys. Rev. Lett. 91, 070402 (2003).
  * Cattaneo  _et al._ (2021) M. Cattaneo, G. L. Giorgi, S. Maniscalco, G. S. Paraoanu, and R. Zambrini, Annalen der Physik 533, 2100038 (2021).
  * Horovitz and Henkel (2021) B. Horovitz and C. Henkel, Phys. Rev. B 104, L081405 (2021).
  * Poyatos  _et al._ (1996) J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 77, 4728 (1996).
  * Pielawa  _et al._ (2007) S. Pielawa, G. Morigi, D. Vitali, and L. Davidovich, Phys. Rev. Lett. 98, 240401 (2007).
  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nature Physics 5, 633 (2009).
  * Menu  _et al._ (2022) R. Menu, J. Langbehn, C. P. Koch, and G. Morigi, Phys. Rev. Res. 4, 033005 (2022).
  * Kempe  _et al._ (2001) J. Kempe, D. Bacon, D. A. Lidar, and K. B. Whaley, Phys. Rev. A 63, 042307 (2001).
  * Lidar and Birgitta Whaley (2003) D. A. Lidar and K. Birgitta Whaley, “Decoherence-free subspaces and subsystems,” in  _Irreversible Quantum Dynamics_, edited by F. Benatti and R. Floreanini (Springer Berlin Heidelberg, Berlin, Heidelberg, 2003) pp. 83–120.
  * Buča  _et al._ (2019) B. Buča, J. Tindall, and D. Jaksch, Nature Communications 10, 1730 (2019).
  * Albert (2019) V. V. Albert, Quantum 3, 151 (2019).
  * Breuer  _et al._ (2016) H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Great Clarendon Street, 2002).
  * Rivas  _et al._ (2010) A. Rivas, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050403 (2010).
  * Benatti  _et al._ (2018) F. Benatti, L. Ferialdi, and S. Marcantoni, Journal of Physics A: Mathematical and Theoretical 52, 035305 (2018).
  * Yu and Eberly (2010) T. Yu and J. Eberly, Optics Communications 283, 676 (2010), quo vadis Quantum Optics?
  * Plenio and Semião (2005) M. B. Plenio and F. L. Semião, New Journal of Physics 7, 73 (2005).
  * Campos Venuti  _et al._ (2006) L. Campos Venuti, C. Degli Esposti Boschi, and M. Roncaglia, Phys. Rev. Lett. 96, 247206 (2006).
  * Giampaolo and Illuminati (2010) S. M. Giampaolo and F. Illuminati, New Journal of Physics 12, 025019 (2010).
  * Pielawa  _et al._ (2010) S. Pielawa, L. Davidovich, D. Vitali, and G. Morigi, Phys. Rev. A 81, 043802 (2010).
  * Wolf  _et al._ (2011) A. Wolf, G. De Chiara, E. Kajari, E. Lutz, and G. Morigi, EPL (Europhysics Letters) 95, 60008 (2011).
  * Fogarty  _et al._ (2013) T. Fogarty, E. Kajari, B. G. Taketani, A. Wolf, T. Busch, and G. Morigi, Phys. Rev. A 87, 050304 (2013).
  * Kajari  _et al._ (2012) E. Kajari, A. Wolf, E. Lutz, and G. Morigi, Phys. Rev. A 85, 042318 (2012).
  * Taketani  _et al._ (2014) B. G. Taketani, T. Fogarty, E. Kajari, T. Busch, and G. Morigi, Phys. Rev. A 90, 012312 (2014).
  * Nicacio and Semião (2016) F. Nicacio and F. L. Semião, Phys. Rev. A 94, 012327 (2016).
  * Nicacio and Semião (2016) F. Nicacio and F. L. Semião, Journal of Physics A: Mathematical and Theoretical 49, 375303 (2016).
  * Bylicka  _et al._ (2014) B. Bylicka, D. Chruściński, and S. Maniscalco, Scientific Reports 4, 5720 (2014).
  * Wendenbaum  _et al._ (2020) P. Wendenbaum, B. G. Taketani, E. Kajari, G. Morigi, and D. Karevski, The European Physical Journal Plus 135, 182 (2020).
  * Braun  _et al._ (2001) D. Braun, F. Haake, and W. T. Strunz, Phys. Rev. Lett. 86, 2913 (2001).
  * Prior  _et al._ (2010) J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
  * Lieb and Robinson (1972) E. Lieb and D. Robinson, Communications in Mathematical Physics 28, 251 (1972).
  * Bethe (1947) H. A. Bethe, Phys. Rev. 72, 339 (1947).
  * Buchheit and Morigi (2016) A. A. Buchheit and G. Morigi, Phys. Rev. A 94, 042111 (2016).
  * Gross and Haroche (1982) M. Gross and S. Haroche, Physics Reports 93, 301 (1982).
  * Scully (2009) M. O. Scully, Phys. Rev. Lett. 102, 143601 (2009).
  * Konovalov and Morigi (2020) A. Konovalov and G. Morigi, Phys. Rev. A 102, 013724 (2020).
  * Strunz  _et al._ (2003) W. T. Strunz, F. Haake, and D. Braun, Phys. Rev. A 67, 022101 (2003).
  * Gualdi and Koch (2013) G. Gualdi and C. P. Koch, Phys. Rev. A 88, 022122 (2013).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Wootters (1998) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * Calabrese and Cardy (2009) P. Calabrese and J. Cardy, Journal of Physics A: Mathematical and Theoretical 42, 504005 (2009).
  * Berman (1994) P. R. Berman,  (1994).
  * Roy  _et al._ (2020) S. Roy, J. T. Chalker, I. V. Gornyi, and Y. Gefen, Phys. Rev. Res. 2, 033347 (2020).
  * (54) D. A. Puente, F. Motzoi, T. Calarco, G. Morigi, and M. Rizzi, arXiv:2305.08641 .
  * Tanese  _et al._ (2013) D. Tanese, H. Flayac, D. Solnyshkov, A. Amo, A. Lemaitre, E. Galopin, R. Braive, P. Senellart, I. Sagnes, G. Malpuech, and J. Bloch, Nature Communications 4, 1749 (2013).
  * Wachter  _et al._ (2019) G. Wachter, S. Kuhn, S. Minniberger, C. Salter, P. Asenbaum, J. Millen, M. Schneider, J. Schalko, U. Schmid, A. Felgner, D. Hüser, M. Arndt, and M. Trupke, Light: Science & Applications 8, 37 (2019).
  * Eichenfield  _et al._ (2009) M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O. Painter, Nature 462, 78 (2009).
  * Ludwig and Marquardt (2013) M. Ludwig and F. Marquardt, Phys. Rev. Lett. 111, 073603 (2013).
  * Porras and Cirac (2004) D. Porras and J. I. Cirac, Phys. Rev. Lett. 93, 263602 (2004).
  * Serafini  _et al._ (2009) A. Serafini, A. Retzker, and M. B. Plenio, New Journal of Physics 11, 023007 (2009).
  * Burrell and Osborne (2007) C. K. Burrell and T. J. Osborne, Phys. Rev. Lett. 99, 167201 (2007).
  * Tran  _et al._ (2020) M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo, A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov, and A. Lucas, Phys. Rev. X 10, 031009 (2020).
  * Harrow  _et al._ (2021) A. W. Harrow, L. Kong, Z.-W. Liu, S. Mehraban, and P. W. Shor, PRX Quantum 2, 020339 (2021).