# Effect of inter-system coupling on heat transport in a microscopic collision
model

Feng Tian1, Jian Zou1,†, Lei Li2, Hai Li3 and Bin Shao1 1School of Physics,
Beijing Institute of Technology, Beijing 100081, China
2School of Physical Science and Technology, Inner Mongolia University, Hohhot
010021, China
3 School of Information and Electronic Engineering, Shandong Technology and
Business University, Yantai 264000, China<EMAIL_ADDRESS>
###### Abstract

In this paper we consider a bipartite system composed of two subsystems each
coupled to its own thermal environment. Based on a collision model, we mainly
study whether the approximation (i.e., the inter-system interaction is ignored
when modeling the system-environment coupling) is valid or not. We also
address the problem of heat transport unitedly for both conventional energy-
preserving system-environment interactions and non-energy preserving system-
environment interactions. For the former interaction, as the inter-system
interaction strength increases, at first this approximation gets worse as
expected, but then counterintuitively gets better even for a stronger inter-
system coupling. For the latter interaction with asymmetry, this approximation
gets progressively worse. In this case we realize a perfect thermal
rectification, and we can not find apparent rectification effect for the
former interaction. Finally and more importantly, our results show that
whether this approximation is valid or not is closely related to the quantum
correlations between the subsystems, i.e., the weaker the quantum
correlations, the more justified the approximation and vice versa.

  * December 2020

Keywords: inter-system coupling, quantum correlation, heat current, thermal
rectification, collision model

##  1 Introduction

In most practical situations, a quantum system inevitably interacts with its
environment which induces decoherence and dissipation [1]. In this case, its
dynamics is usually described by Gorini-Kossakowski-Lindblad-Sudarshan (GKLS)
quantum master equation with a series of approximations [2]. When deriving a
master equation for an open quantum system, one may always obtain a global
master equation which considers the full system Hamiltonian, i.e., including
the direct coupling between the subsystems [3, 4, 5, 6]. However, such
derivation is very complicated when the system is composed of two or more
interacting subsystems. Hence a local master equation that ignores the direct
interactions between the subsystems is often used as a substitute [3, 4, 5, 6,
7, 8, 9, 10, 11]. By means of these two kinds of master equations, a lot of
efforts has recently been devoted to the heat transport for thermodynamic
systems [7, 8, 9, 10, 11, 12, 13, 14, 15]. However, by comparing the dynamics
resulting from the corresponding master equations with exact numerical
simulations, both approaches may lead to seeming thermodynamic inconsistencies
or just suit to some parameter regimes, as pointed out in references [15, 16,
17, 18, 19, 20, 21, 22]. In reference [3] the local description for two
coupled quantum nodes may predict heat currents from a cold to a hot thermal
reservoir, or the existence of currents even in the absence of a temperature
gradient. The origin of these effects, as discussed in references [5, 23],
lies in the fact that there is an external work cost related to the breaking
of global detailed balance. By including this work cost, this inconsistencies
can be resolved. For weak inter-system coupling, it was shown that the global
approach fails in non-equilibrium situations, whereas the local approach
agrees with the exact solution [19]. In reference [21], due to a failure of
the secular approximation, it was reported that the global master equation
erroneously gave a vanishing heat current through a spin-$\frac{1}{2}$
Heisenberg chain in the presence of a finite-temperature gradient. Recently,
for two coupled qubits interacting with common and separate baths, it was
shown that the global approach with partial secular rather than full secular
approximation always provides the most accurate choice for the master equation
[24]. In reference [25], it was shown that the completely positive version of
the Redfield equation obtained using coarse-grain and an appropriate time-
dependent convex mixture of the local and global solutions gives rise to the
most accurate semigroup approximations of the whole exact system dynamics.
Moreover, the two approaches above, from the viewpoint of the system, often
rely on approximate Markovian master equations derived under the assumptions
of weak system-environment coupling, which would become challenging under
strong coupling.

Furthermore, the manipulation of heat transport in non-equilibrium steady-
state has been identified as one of the crucial studies of quantum
thermodynamics, which gives us an improved understanding of classical
thermodynamics in quantum domain [26, 27, 28, 29, 30, 31, 32, 33, 34]. For
example, heat transport between two bosonic reservoirs was predicted for a
coupled two-state system, and a formula for thermal conductance was derived
based on a rate equation formalism [32]. Bandyopadhyay  _et al_. described a
numerical scheme for exactly simulating the heat transport in a quantum
harmonic chain with self-consistent reservoirs [33]. By means of an effective
harmonic Hamiltonian, a quantum thermal transport through anharmonic systems
was studied within the framework of the nonequilibrium Green’s function method
[34]. Besides, quantum devices such as heat rectifier, thermal memory, and
thermal ratchet, have also become goals of controlling thermal transport in
quantum thermodynamics [35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. It was found
that, by using the quantum master equation, thermal rectification in
anisotropic Heisenberg spin chains could change sign when the external
homogeneous magnetic field was varied [42]. An optimal rectification in the
ultrastrong-coupling regime of two coupled two-level systems was shown [43].
In [44], Jose  _et al_. studied two interacting spin-like systems
characterized by different excitation frequencies, which can be used as a
quantum thermal diode.

Recently, collision model, also called repeated interactions, has drawn
attention for its potential advantage in simulating open quantum system [45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. It was assumed that the environment
consists of a large collection ancillas and the system of interest interacts,
or collides, with an ancilla at each time step. In the framework of collision
model, a continuous-time description in terms of a Lindblad master equation
can be derived in the short-time limit provided some assumptions are made
about the system-ancilla interaction [39, 48]. Recently, for instance, in
reference [56] the system’s dynamics embodied by the stroboscopic map can be
approximated by a Lindblad master equation in a short-time limit. In a similar
way, in reference [5] a local master equations with a Lindblad form for two
coupled harmonic oscillator was derived by using the method of repeated
interactions. Moreover, also from the viewpoint of the system, the
corresponding reduced dynamics can be obtained in many cases without any
approximations [57, 58, 59, 60, 61, 62, 63]. This is because collision model
allows for the possibility to decompose a complicated open dynamics in terms
of discrete elementary processes. It is particularly suited for addressing the
thermodynamics of engineered reservoirs [23, 47, 48, 60, 64, 65]. For example,
under an energy-preserving system-environment interaction within the framework
of collision model, it has been found that the non-monotonic time behaviour of
the heat exchange between the system and environment could serve as an
indicator of non-Markovian behaviors [60]. Besides, based on a microscopic
collision model, heat current between a coupled system can flow from the cold
nonthermal reservoir to the hot one due to the contribution of coherence [64].
And a link between information and thermodynamics, for a multipartite open
quantum system with a finite temperature reservoir, was displayed in term of a
collision model [65].

In this paper, we consider a bipartite system which interacts with its local
environment consisting of a large collection of identical ancillas. Inspired
by the previous work on the local and global master equations, we consider an
approximate situation within the framework of collision model, specifically,
environment acts on subsystem without considering the inter-system
interaction, i.e., ignores the direct interaction between the subsystems when
modeling the system-environment coupling. Note that although this
approximation in this collision model is not strictly a substitute for the
local master equation we described above, we expect that the results derived
from this simple and solvable model can provide a reference for those in more
involved but less tractable models. We mainly study whether this approximation
is valid or not, i.e., whether or when this direct interaction can be ignored.
We systematically examine the heat transport through the system in both weak
and strong system-ancilla coupling. Within the framework of collision model,
we give a more general definition of heat current under non-energy preserving
system-ancilla interactions (there is external work performing on the whole
system-ancilla compound). In the case of conventional energy-preserving
system-ancilla interactions, we find that the results predicted in the case of
ignoring the inter-system interaction is valid even at strong inter-system
interactions. In particular, we realize a perfect rectification in the
asymmetric systems and give a discussion about the mechanism of this
phenomenon. Moreover, we find that whether or not this approximation (i.e.,
ignores the direct interaction between the subsystems) is valid is closely
related to the quantum correlations between the subsystems.

Figure 1: Schematic sketch of a bipartite system $S$ made up of two
interacting subsystems connected to two independent subenvironments. In the
$n$th round of the dynamics, after a free evolution of the whole system,
$S_{1}$ interacts with $E_{n}^{1}$ and next $S_{2}$ interacts with
$E_{n}^{2}$. The system then moves to the $(n+1)$th round and this process is
repeated over and over.

##  2 Model

We consider a bipartite system $S$ consisting of two identical two-level
subsystems $S_{i}$ $({i}=1,2)$, and each subsystem is coupled to its local
thermal subenvironment $E_{i}$ set at temperature $T_{i}$. Here the
subenvironment $E_{i}$ is a sequence of non-interacting ancillas
$(E_{1}^{i},E_{2}^{i},...,E_{n}^{i})$ all in the same initial state
$\eta_{n}^{i}$. Subsystem $S_{i}$ interacts with the connected subenvironment
$E_{i}$ via a series of short subsystem-ancilla interactions. The joint state
of system and environment is initially factorized:

$\rho^{SE}(0)=\rho^{S}(0)(\otimes_{j=1}^{n}\eta_{j}^{1})(\otimes_{j=1}^{n}\eta_{j}^{2})$
(1)

where $\rho^{S}(0)$ and $\otimes_{j=1}^{n}\eta_{j}^{i}$ are the initial states
of the open system and sub-environment $E_{i}$, respectively.

In figure 1 we show a schematic sketch of the collision model considered. It
is composed of a series of repeated rounds. In any one round of the collision
process, first the whole system $S$ undergoes a free evolution lasting a time
interval $\tau$, subsequently $S_{1}$ and $S_{2}$ locally interact with only
one ancilla of its environment, respectively. Then by tracing out the
environment’s degrees of freedom and repeating the above process in the next
round, we can obtain the reduced state of the system in the full time
evolution. Note that the subenvironment is assumed to be large enough, so that
the subsystem never collides twice with the same ancilla. As a consequence, at
each collision round $n$, the subsystem $S_{i}$ collides with a ”fresh”
$E_{n}^{i}$.

We assume throughout this paper that $S$ and each ancilla $E_{n}^{i}$ of the
sub-environment $E_{i}$ are qubits with logical states
$\\{|0\rangle,|1\rangle\\}$. And the corresponding free Hamiltonians for the
subsystem $S_{i}$ $({i}=1,2)$ and the ancilla $E_{n}^{i}$ are
$\hat{H}_{S_{i}}=1/2\omega_{i}\hat{\sigma}_{z}$ and
$\hat{H}_{E_{n}^{i}}=1/2\omega_{0}\hat{\sigma}_{z}$, respectively. Here
$\hat{\sigma}_{z}$ is the usual Pauli operator (we set $\hbar=1$). The free
dynamics of the whole system $S$ is described by the unitary evolution
operator

$\hat{U}_{S_{1},S_{2}}=\exp{[-i(\hat{H}_{0}+\hat{H}_{int}^{S_{1},S_{2}})\tau]},$
(2)

with $\hat{H}_{0}=\hat{H}_{S_{1}}+\hat{H}_{S_{2}}$ and the inter-system
interaction Hamiltonian $\hat{H}_{int}^{S_{1},S_{2}}$. We use unitary operator
$\hat{V}_{S_{i},E_{n}^{i}}$ to model the collisions between the subsystem
$S_{i}$ and ancilla $E_{n}^{i}$ of sub-environment $E_{i}$ (Its exact
definition will be given later in different cases). And it is assumed that all
the subsystem-ancilla collisions have the same duration $\tau$. So the
dynamical map $\hat{\Lambda}_{S_{1},S_{2}}$ that governs the free evolution of
the system $S$ and map $\hat{\Psi}_{S_{1},E_{n}^{i}}$ that governs the system-
subenvironment coupling at the $n$th round can be written as

$\hat{\Lambda}_{S_{1},S_{2}}(\rho)=\hat{U}_{S_{1},S_{2}}\rho\hat{U}_{S_{1},S_{2}}^{\dagger},$
(3)
$\hat{\Psi}_{S_{i},E_{n}^{i}}(\rho)=\hat{V}_{S_{i},E_{n}^{i}}\rho\hat{V}_{S_{i},E_{n}^{i}}^{\dagger},$
(4)

respectively. Following the repeated interaction approach mentioned above, the
joint state of system and ancillas is brought from the $n$th round to the
$(n+1)$th round through the process

$\rho^{S}_{n-1}\otimes\eta^{1}_{n}\otimes\eta^{2}_{n}\rightarrow\rho^{SE}_{n}=\hat{\mathcal{U}}[\rho^{S}_{n-1}\otimes\eta^{1}_{n}\otimes\eta^{2}_{n}]\hat{\mathcal{U}}^{\dagger},$
(5)

where
$\hat{\mathcal{U}}=\hat{V}_{S_{2},E_{n}^{2}}\hat{V}_{S_{1},E_{n}^{1}}\hat{U}_{S_{1},S_{2}}$.
Hence after round $n$, the reduced density matrix of the system $\rho^{S}_{n}$
is

$\rho^{S}_{n}=\mathrm{Tr}_{E_{n}^{1},E_{n}^{2}}[\rho^{SE}_{n}],$ (6)

where $\mathrm{Tr}_{E_{n}^{1},E_{n}^{2}}[\cdot]$ denotes the partial trace
over the two ancillas $E_{n}^{1}$ and $E_{n}^{2}$. Similarly, the reduced
state $\tilde{\eta}_{n}^{1(2)}$ of the $n$th ancilla of subenvironment
$E_{1(2)}$ is

$\tilde{\eta}_{n}^{1(2)}=\mathrm{Tr}_{SE_{n}^{2(1)}}[\rho^{SE}_{n}].$ (7)

Throughout we assume each ancilla of subenvironment $E_{i}$ to be initially in
a thermal state with inverse temperature $\beta_{i}=1/T_{i}$ (we set $k=1$),
namely,

$\eta_{n}^{i}=\frac{1}{Z}\exp{(-\beta_{i}\hat{H}_{E_{n}^{i}})},$ (8)

where $Z=\mathrm{Tr}[\exp{(-\beta_{i}\hat{H}_{E_{n}^{i}})}]$ is the partition
function.

Now we begin to give the subsystem-ancilla evolution operator in every round
of the collision model considered. When ancilla $E_{n}^{i}$ collides with
subsystem $S_{i}$ not ignoring the direct interaction between the two
subsystems, the corresponding unitary time evolution operator can be written
as

$\hat{V}_{S_{i},E_{n}^{i}}=\exp{[-i(\hat{H}_{0}^{\prime}+\hat{H}_{int}^{S_{1},S_{2}}+\hat{H}_{int}^{S_{i},E_{n}^{i}})\tau]},$
(9)

where
$\hat{H}_{0}^{\prime}=\hat{H}_{S_{1}}+\hat{H}_{S_{2}}+\hat{H}_{E_{n}^{i}}$ and
$\hat{H}_{int}^{S_{i},E_{n}^{i}}$ is the interaction Hamiltonian between the
subsystem $S_{i}$ and ancilla $E_{n}^{i}$. However, it is convenient to ignore
the direct interaction between the two subsystems, namely, neglect the
interaction between the subsystems when modeling the system-environment
interaction. In this case, the corresponding unitary time evolution operator
is written as

$\hat{V}_{S_{i},E_{n}^{i}}^{app}=\exp{[-i(\hat{H}_{0}^{\prime\prime}+\hat{H}_{int}^{S_{i},E_{n}^{i}})\tau]},$
(10)

where $\hat{H}_{0}^{\prime\prime}=\hat{H}_{S_{i}}+\hat{H}_{E_{n}^{i}}$.

It is obvious that the difference between equations (9) and (10) resides in
the subsystem-ancilla interaction in each round: equation (9) arises naturally
when modeling the subenvironment-subsystem coupling from a microscopic model
considering the full system Hamiltonian, i.e., not ignore the direct
interaction between its subsystems, while equation (10) ignores this
interaction.

##  3 Symmetric system

In this section we consider a symmetric system, from which to investigate
whether the inter-system interaction can be ignored or not. Among the possible
choices for the interaction between the subsystem $S_{i}$ and their connected
ancilla $E_{n}^{i}$, we choose the conventional energy-preserving interaction
Hamiltonian as

$\hat{H}_{int}^{S_{i},E_{n}^{i}}=\gamma(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y})$
(11)

with subsystem-subenvironment coupling strength $\gamma$. And the interaction
between $S_{1}$ and $S_{2}$ takes the same form

$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y})$
(12)

with coupling strength $\delta$. And we consider that the energy gaps of
subsystems are the same, i.e., $\omega_{1}=\omega_{2}=\omega_{0}$.

###  3.1 Heat current

Now we are going to investigate the heat currents to evaluate the performance
of this approximation. We also consider the heat currents in the case of not
ignoring the inter-system interactions, which serves as a benchmark. In the
case of ignoring the inter-system interactions, for the system characterized
by Hamiltonians of equations (11) and (12), its corresponding unitary system-
environment operator preserves the energy, i.e.,

$[\hat{V}_{S_{i},E_{n}^{i}}^{ign},\hat{H}_{S_{i}}+\hat{H}_{E_{i}}]=0.$ (13)

This implies that all the energy leaving the ancilla enters the system. Hence,
during the $(n+1)$th round, the heat exchange $\triangle{Q}$ between system
$S_{1(2)}$ and ancilla $E_{n}^{1(2)}$ is given by

$\triangle{Q}_{E_{n}^{1(2)}}=\mathrm{Tr}[\hat{H}_{E_{n}^{1(2)}}(\tilde{\eta}_{n}^{1(2)}-\eta_{n}^{1(2)})].$
(14)

Due to $\triangle{Q}_{E_{n}^{1}}=-\triangle{Q}_{E_{n}^{2}}$ for steady state,
the stationary heat current flowing from subenvironment $E_{1}$ to the system
$S_{1}$ can be defined as

$J_{h}=-\triangle{Q}_{E_{n}^{1}}.$ (15)

This definition was often used when equation (13) is satisfied [64, 56, 41,
42, 43, 44]. However, in the case of not ignoring the inter-system
interactions, the corresponding unitary system-environment operator no longer
preserves the energy because of

$[\hat{V}_{S_{i},E_{n}^{i}}^{con},\hat{H}_{S_{1}}+\hat{H}_{S_{2}}+\hat{H}_{int}^{S_{1},S_{2}}+\hat{H}_{E_{i}}]\neq
0.$ (16)

That is, not all the energy leaving the ancilla enters the system. So an
external work is required to turn the system-ancilla coupling on and off [65,
66], which we refer to as the switching work (labeled $W_{SW}$). In each round
of collisions illustrated in figure 1, the corresponding work $W_{SW}$ on the
system can be written as

$W_{SW}=W_{1}+W_{2}.$ (17)

The first term $W_{1}$ is from the sudden on and off switching of
$S_{1}-E_{n}^{1}$ interaction, and it reads

$W_{1}=\mathrm{Tr}[\hat{H}_{int}^{S_{1},E_{n}^{1}}(\rho^{\prime}_{S_{1}S_{2}}\otimes\eta_{n}^{1}-\rho_{S_{1}S_{2}E_{n}^{1}})],$
(18)

where $\rho^{\prime}_{S_{1}S_{2}}$ is the state of system before
$S_{1}-E_{n}^{1}$ collision, and $\rho_{S_{1}S_{2}E_{n}^{1}}$ is the global
state of system $S$ and ancilla $E_{n}^{1}$ after $S_{2}-E_{n}^{1}$ collision.
Similarly, the second term associated with $S_{2}-E_{n}^{2}$ interaction takes
the form

$W_{2}=\mathrm{Tr}[\hat{H}_{int}^{S_{2},E_{n}^{2}}(\rho^{\prime\prime}_{S_{1}S_{2}}\otimes\eta_{n}^{2}-\rho_{S_{1}S_{2}E_{n}^{2}})],$
(19)

where $\rho^{\prime\prime}_{S_{1}S_{2}}$ is the reduced state of system after
$S_{1}-E_{n}^{1}$ collision, and $\rho_{S_{1}S_{2}E_{n}^{2}}$ is the global
state of system $S$ and ancilla $E_{n}^{2}$ after $S_{2}-E_{n}^{2}$ collision.
As a result, equation (15) is no longer valid to calculate the heat current in
the case of not ignoring the inter-system interactions.

Based on this collision model, now let us focus on deriving the expression of
heat current that can be applicable in the cases of ignoring and not ignoring
the inter-system interactions, as well as for weak and strong system-
environment couplings. Physically, heat current is defined as the energy going
through the system. We consider the energy change of the system in a complete
round (i.e., free evolution of $S$, $S_{1}-E_{n}^{1}$ and $S_{2}-E_{n}^{2}$
interactions). The energy change $\triangle{E}_{S}^{1}$ of system $S$ in each
$S_{1}-E_{n}^{1}$ interaction can be written as

$\triangle{E}_{S}^{1}=\mathrm{Tr}[(\hat{H}_{0}+\hat{H}_{int}^{S_{1}S_{2}})(\rho^{\prime\prime}_{S_{1}S_{2}}-\rho^{\prime}_{S_{1}S_{2}})].$
(20)

Similar expressions can be obtained for $S_{2}-E_{n}^{2}$ interaction, and the
corresponding energy change $\triangle{E}_{S}^{2}$ of system $S$ is given by

$\triangle{E}_{S}^{2}=\mathrm{Tr}[(\hat{H}_{0}+\hat{H}_{int}^{S_{1}S_{2}})(\rho^{\prime\prime\prime}_{S_{1}S_{2}}-\rho^{\prime\prime}_{S_{1}S_{2}})],$
(21)

where $\rho^{\prime\prime\prime}_{S_{1}S_{2}}$ is the reduced state of system
$S$ after $S_{2}-E_{n}^{2}$ collision. Since the inter-system dynamics is
unitary, the energy of system is preserved in each free evolution. For steady
state, the state of the system cannot change in a complete round, leading to

$\triangle{E}_{S}^{1}=-\triangle{E}_{S}^{2}.$ (22)

From the above discussion, it is clear that $\triangle{E}_{S}^{1}$ is the
energy going through the system. Thus, the stationary heat current $J_{h}$
flowing from $E_{1}$ to $E_{2}$ can be rewritten as

$J_{h}=\triangle{E}_{S}^{1}.$ (23)

In the case of ignoring the inter-system interactions, if the corresponding
system satisfies the commutation condition equation (13), i.e., all energy
changes in the system can be attributed to energy flowing to or from the
ancilla, it is thus clear that equations (23) reduces to equation (15).

Here we should mention that exchange of energy between two quantum systems,
characterized by particular commutation relation between the local
Hamiltonians and the interaction operator, can always be split into work and
heat [67, 68]. And this formalism was extended to situations where one, or
both, subsystems are coupled to a thermal environment. In the case of not
ignoring the inter-system interaction, according to reference [68] the
stationary heat current and work in each $S_{1}-E_{n}^{1}$ collision are
calculated as

$J_{h}^{\prime}=-i\int_{0}^{\tau}\mathrm{Tr}{[\mathbb{I}\otimes(\hat{H}_{0}+\hat{H}_{int}^{S_{1},S_{2}}),\hat{H}_{int}^{S_{1},E_{n}^{1}}\otimes\mathbb{I}]C_{12}}dt,$
(24)

$W=-i\int_{0}^{\tau}\mathrm{Tr}{[\hat{H}_{S}^{eff},\hat{H}_{S}]\tilde{\eta}_{n}^{1}}dt,$
(25)

where $C_{12}$ represents the correlations between the system $S$ and
$E_{n}^{1}$ after $S_{1}-E_{n}^{1}$ collision, and $\hat{H}_{S}^{eff}$ is the
diagonal part of matrix
$\mathrm{Tr}_{E_{n}^{1}}[\hat{H}_{int}^{S_{1},E_{n}^{1}}(\tilde{\eta}_{n}^{1}\otimes\mathbb{I})]$.
After a lot of numerical calculations, we find that the heat current
$J_{h}^{\prime}$ given by equation (24) is equal to that given by equation
(23). Thus our previous definition of equation (23) is consistent with
reference [68].

In figure 2 we plot $J_{h}$ as a function of the inter-system coupling
strength $\delta$ for various subsystem-ancilla coupling strength $\gamma$ in
the cases of ignoring and not ignoring the inter-system interactions. We fix
the initial state of the system to be $|11\rangle$ and set $T_{1}=5\omega_{0}$
and $T_{2}=\omega_{0}$ for subenvironments $E_{1}$ and $E_{2}$, respectively.
It can be seen from figure 2 that, for fixed $\gamma$, the heat currents in
both cases increase rapidly with the increase of $\delta$. And they eventually
reach the same steady value, before which the heat current in the case of
ignoring the inter-system interaction is always smaller than that of not
ignoring the inter-system interactions. In general, for weak inter-system
coupling, one would expect smaller deviation from not ignoring the inter-
system interaction. As expected, in the limit $\delta\rightarrow 0$, it gives
the same (zero) heat current as that of not ignoring the inter-system
interactions. And for fixed $\gamma$, when $\delta$ increases we observe that
the heat current predicted from the approximation gradually deviates from that
of not ignoring the inter-system interaction, and this deviation quickly
increases from zero to its maximum. However, it might seem counterintuitive
that this deviation gradually decreases with $\delta$ and eventually vanishes
at larger values of $\delta$. That is, for larger $\delta$, heat currents from
the approximation can still be consistent with that of not ignoring the inter-
system interaction, and this approximation consequently may not necessarily
break down.

Figure 2: Steady heat currents $J_{h}$ as a function of $\delta$ for different
$\gamma$, $\gamma=0.2\omega_{0}$ (blue), $\gamma=0.5\omega_{0}$ (red) and
$\gamma=0.8\omega_{0}$ (green). The solid and dashed lines correspond to not
ignoring and ignoring the inter-system interaction, respectively. The inset is
the magnified $J_{h}$ for $\gamma=0.2\omega_{0}$. In both cases the system is
initialized in $|11\rangle$ and each ancilla is initialized in its thermal
state. The plots are obtained for $\omega_{1}=\omega_{2}=\omega_{0}$,
$\omega_{0}\tau=0.1$, $T_{1}=5\omega_{0}$ and $T_{2}=\omega_{0}$.

Moreover, we find that the difference of the heat currents between those of
ignoring and not ignoring the inter-system interaction strongly depends on
$\gamma$, as shown in figure 2. A bigger difference is obtained for stronger
coupling strength $\gamma$. For instance, for $\gamma=0.8\omega_{0}$, the
difference which maintains nonzero values is within a larger region of
$\delta$, compared with $\gamma=0.2\omega_{0}$ and $0.5\omega_{0}$. That is,
when assessing the heat currents the results of ignoring the inter-system
interactions is inconsistent with that in the case of not ignoring the inter-
system interactions. Smaller difference is obtained for weaker coupling
strength $\gamma$. For $\gamma=0.2\omega_{0}$, it can be seen from the inset
of figure 2 that the difference between two cases maintaining nonzero values
is within a very small region below $\delta\sim 0.04\omega_{0}$. In other
words, for smaller values of $\gamma$, the results of ignoring the inter-
system interactions agrees well with that of not ignoring the inter-system
interactions, and consequently it also approximately predicts the correct heat
currents within almost all region of $\delta$. Physically, this can be easily
understood as following: decreasing the system-ancilla coupling strength would
weaken the influence of the inter-system coupling ignored.

###  3.2 Trace distance

Figure 3: Trace distance $\mathcal{D}_{\mathcal{T}}$ between the density
matrices obtained from ignoring and not ignoring the inter-system interactions
against $\delta$ at steady state. The blue line, the red line and the green
line correspond to $\gamma=0.2\omega_{0}$, $\gamma=0.5\omega_{0}$ and
$\gamma=0.8\omega_{0}$, respectively. The other parameters are the same as
those in figure 2. Inset shows $\delta^{\prime}$ for the maximum of
$\mathcal{D}_{\mathcal{T}}$ as a function of $\gamma$.

When investigating the accuracy of the local and global master equations,
steady state is often used as a reference for predicting the results of the
nonequilibrium dynamics [16, 19]. To further assess the validity of this
approximation (i.e., ignore the inter-system interaction), we also consider
the obtained steady states. As a measure of distinguish ability, here we
examine the trace distance between the steady state obtained from the case of
ignoring the inter-system interactions and the steady state obtained from the
case of not ignoring the inter-system interactions:

$\mathcal{D}_{\mathcal{T}}=\frac{1}{2}\parallel\rho^{S}-\rho^{Sapp}\parallel_{1},$
(26)

where $\parallel\cdot\parallel_{1}$ is the trace norm while $\rho^{S}$ and
$\rho^{Sapp}$ are the reduced steady states of the system in the cases of not
ignoring and ignoring the inter-system interactions, respectively. The trace
distance is equal to unity for fully distinguishable states, while it is null
for identical states. We assume that the system $S$ in both cases evolves from
the same initial state $\rho^{S}(0)$. Figure 3 shows the dependence of
$\mathcal{D}_{\mathcal{T}}$ on $\delta$ for different strength $\gamma$. For
fixed $\gamma$, it can be seen that $\mathcal{D}_{\mathcal{T}}$ quickly
increases and then gradually decreases with the increase of $\delta$ (e.g.,
for $\gamma=0.1\omega_{0}$). That is, when $\delta$ increases, the steady
state $\rho_{n}^{Sapp}$ predicted from this approximation gradually deviates
from the steady state $\rho_{n}^{S}$ predicted without the approximation, then
this deviation reaches its maximum and eventually decreases. It is obvious
that this behavior of the steady state from this approximation is indeed
similar to that of heat current. And here we also arrive at a conclusion
similar to that of heat current: even for a large $\delta$, this approximation
can be well justified when predicting the steady state. The inset of figure 3
plots the value of $\delta^{\prime}$ for the maximum of
$\mathcal{D}_{\mathcal{T}}$ as a function of $\gamma$ (i.e., when
$\delta=\delta^{\prime}$, $\mathcal{D}_{\mathcal{T}}$ reaches its maximum for
fixed $\gamma$). It can be seen that $\delta^{\prime}$ nonlinearly increases
with the subenvironment-ancilla interaction strength $\gamma$.

Figure 3 also shows the effect of system-ancilla coupling strength $\gamma$ on
$\mathcal{D}_{\mathcal{T}}$. For $\gamma=0.2\omega_{0}$, it can be seen that
$\mathcal{D}_{\mathcal{T}}$ maintains a larger value only within a smaller
region of $\delta$, and quickly decreases compared to $\gamma=0.5\omega_{0}$
and $\gamma=0.8\omega_{0}$. Again, for weak $\gamma$, steady state predicted
from this approximation agrees well with that of not ignoring the inter-system
interactions even at strong interaction strengths $\delta$. In other words,
this approximation is well justified to describe steady state in this case.

Figure 4: Quantum discord $\mathcal{D}$ of the qubit system as a function of
$\delta$, without approximation. The blue line and red line correspond to
$\gamma=0.2\omega_{0}$ and $\gamma=0.5\omega_{0}$, respectively. Inset a is
the magnified $\mathcal{D}$ for $\gamma=0.2\omega_{0}$. Inset b corresponds to
$\gamma=0.8\omega_{0}$. The remaining parameters are the same as those in
figure 2.

###  3.3 Quantum discord

Why can the approximation (i.e., ignore the inter-system interactions) give a
good estimate of heat current and steady state even for a large $\delta$ ? In
references [70, 71], it was found that the compositeness (two particles behave
like a single particle) is closely related to the quantum correlations between
the constituent particles. It is thus very interesting to consider the quantum
correlations between the bipartite system $S_{1}$ and $S_{2}$, which can be
captured by the discord [69]

$\mathcal{D}(\rho_{S})={\min}_{\Pi_{A}}\\{I(\rho_{S})-I(\Pi_{A}\rho_{S})\\}.$
(27)

Here $\Pi_{A}$ is a set of rank-one POVM projectors on system $S_{1}$, and
$I(\rho_{S})=S(\rho_{S})-S(\rho_{S_{1}})-S(\rho_{S_{2}})$ is the quantum
mutual information associated with the von Neumann entropy. Without
approximation, in figure 4 we plot $\mathcal{D}$ as a function of $\delta$. It
can be seen that, for fixed $\gamma$, $\mathcal{D}$ first increases from zero
to its maximum with the increase of $\delta$, then it gradually decreases,
i.e., the quantum correlations between two subsystems first increases and then
decreases. It is obvious that such behavior is associated with the heat
current or the steady state, i.e., the greater the quantum correlation
$\mathcal{D}$ between two subsystems, the greater the deviation of the heat
current and steady state predicted from the approximation (i.e., ignore the
inter-system interactions), and vice versa. In general, one would think that
whether this approximation is valid or not should depend on inter-system
coupling strength. However, here we find that it depends strongly on quantum
correlations between the two qubits, i.e., the higher the quantum
correlations, the worse the approximation, and vice versa. In other words,
higher correlations means these two subsystems behave more like a composite
and can not be treated separately. As a result, this approximation gets worse.
Only when the correlations becomes smaller, can it get better.

Moreover, for this model, we also consider the other situations by replacing
equations (11) and (12) with any of the following three combinations of
conventional energy-preserving interactions:

$\bullet$
$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z})$
and
$\hat{H}_{int}^{S_{i},E_{n}^{i}}=\gamma(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y})$.

$\bullet$
$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y})$
and
$\hat{H}_{int}^{S_{i},E_{n}^{i}}=\gamma(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z})$.

$\bullet$
$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z})$
and
$\hat{H}_{int}^{S_{i},E_{n}^{i}}=\gamma(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z})$.

For these three combinations, we also investigate the heat current, trace
distance and quantum discord with and without ignoring the inter-system
interactions. The corresponding results are similar to those from equations
(11) and (12). It is noted that, in the case of ignoring the inter-system
interaction, the corresponding unitary operators of all three combinations
above, like equations (11) and (12), also satisfy the condition
$[\hat{V}_{S_{i},E_{n}^{i}}^{ign},\hat{H}_{S_{i}}+\hat{H}_{E_{i}}]=0$, i.e.,
they correspond to energy-preserving system-environment interactions and no
external work exists. In the case of not ignoring the inter-system
interaction, in each round (i.e., free evolution of $S$, $S_{1}-E_{n}^{1}$ and
$S_{2}-E_{n}^{2}$ interactions) there is an external work due to
$[\hat{V}_{S_{i},E_{n}^{i}}^{con},\hat{H}_{S_{1}}+\hat{H}_{S_{2}}+\hat{H}_{int}^{S_{1},S_{2}}+\hat{H}_{E_{i}}]\neq
0$. But after a lot of numerical calculations we find that the external work
in each round is much smaller than the heat flowing through the system in each
round.

##  4 Asymmetric system

In this section, we turn our attention to the systems featured by distinct
asymmetries, such as the qubits with different frequencies or asymmetric
couplings to their reservoirs. Indeed, inspired by the manipulation and
control of the thermal transport in micro-scale, such asymmetric systems were
widely investigated in the thermal diode and thermal transistor [36, 37, 38,
39, 40, 41, 42, 43, 44]. In the following, we will study the heat transport
and, at the same time, study whether the inter-system interaction can be
ignored or not when modeling the system-environment interaction.

###  4.1 Off-resonant interacting qubits

First we consider two qubits are characterized by different energy gaps
$\omega_{1}$ and $\omega_{2}$, i.e., $\omega_{1}\neq\omega_{2}$. The
interaction Hamiltonian between two subsystems $S_{1}$ and $S_{2}$ is given as

$\hat{H}_{int}^{S_{1},S_{2}}=\delta\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}.$
(28)

Such system may arise, for example, when a nonuniform magnetic field is
applied to a pair of interacting spins in the $z$ directions. The coupling to
each ancilla of environment is described by the following Hamiltonian:

$\hat{H}_{int}^{S_{i},E_{n}^{i}}=\gamma\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}.$
(29)

For this interaction, in both cases with and without ignoring the inter-system
interactions, not all the energy leaving the ancilla enters the system due to
$[\hat{V}_{S_{i},E_{n}^{i}}^{ign},\hat{H}_{S_{i}}+\hat{H}_{E_{i}}]\neq 0$ and
$[\hat{V}_{S_{i},E_{n}^{i}}^{con},\hat{H}_{S_{1}}+\hat{H}_{S_{2}}+\hat{H}_{int}^{S_{1},S_{2}}+\hat{H}_{E_{i}}]\neq
0$. Accordingly, after a lot of numerical calculations we find that, different
from the results of section 3 where the external work in each round is much
smaller than the heat flowing through the system in each round, in this case
the external work in each round is almost the same order as the heat flowing
through the system in each round. As a result, equation (15) cannot be used to
calculated the heat current. In particular, when the two qubits are uncoupled,
equation (15) even gives (unphysical result) nonzero heat currents for these
two cases. As mentioned before, this is because an external work is required
to turn the subsystem-ancilla coupling on and off.

Figure 5: $J_{h}$ as a function of $\delta$ with various temperature
difference in the case of not ignoring the inter-system interaction. We use
solid line for $T_{1}=10\omega_{0}$ and different $T_{2}$,
$T_{2}=0.1\omega_{0}$ (green), $4\omega_{0}$ (red), and $8\omega_{0}$ (blue);
dashed line for $T_{2}=10\omega_{0}$ and different $T_{1}$,
$T_{1}=0.1\omega_{0}$ (green), $4\omega_{0}$ (red), and $8\omega_{0}$ (blue).
The two qubits are off-resonant with $1/2\omega_{1}=\omega_{2}=\omega_{0}$.
Here we set $\gamma=0.3\omega_{0}$, and the other parameters are the same as
those in figure 2. Insets a, b and c show the corresponding rectification
factor $R$ for $T_{1}=10\omega_{0}$ and different $T_{2}$,
$T_{2}=8\omega_{0}$, $4\omega_{0}$, and $0.1\omega_{0}$, respectively.

From equation (23), we calculate the stationary heat currents for system
characterized by Hamiltonians of equations (28) and (29), also addressing for
ignoring and not ignoring the inter-system interaction separately. From our
study we find that the heat current in the case of ignoring the inter-system
interaction is always zero even in the presence of a finite temperature
gradient $T_{1}>T_{2}$ or $T_{2}>T_{1}$, regardless of the coupling strength
$\delta$ and $\gamma$. This result deviates from that in the case without the
approximation, in which heat currents increases as $\delta$ increases, so this
approximation gets progressively worse with the increase of $\delta$.

Moreover, for this model we compare the heat currents for various temperature
difference $\triangle{T}$ between two environments, which shows an asymmetric
heat transport of the system (i.e., heat rectification). To quantify this
rectification efficiency, we also give a rectification factor $R$ defined as
follows [42, 43, 44]:

$R=\frac{|J_{h}(\triangle{T})+J_{h}(-\triangle{T})|}{\max(|J_{h}(\triangle{T})|,|J_{h}(-\triangle{T})|)},$
(30)

where $J_{h}(\triangle{T})$ is the forward heat current for $T_{1}>T_{2}$, and
$J_{h}(-\triangle{T})$ is the backward heat current when the temperature
gradient is reversed. It can be seen from figure 5 that in the case
$T_{1}=10\omega_{0}$ and $T_{2}=8\omega_{0}$, its corresponding forward heat
current is greater than the backward heat currents ($T_{2}=10\omega_{0}$ and
$T_{1}=8\omega_{0}$). Hence an asymmetric conduction
($|J_{h}(\triangle{T})|\neq|J_{h}(-\triangle{T})|$) emerges which varies with
the increase of $\delta$. If the temperature gradient raises, such as
$T_{1}=10\omega_{0}$ and $T_{2}=4\omega_{0}$, this asymmetric conduction
increases, i.e., the higher the gradient, the stronger the heat rectification.
Especially, when $T_{1}=10\omega_{0}$ and $T_{2}=0.1\omega_{0}$, its
corresponding backward heat current [cf. green dashed line corresponds to
$T_{1}=0.1\omega_{0}$ and $T_{2}=10\omega_{0}$ in figure 5] is almost zero for
any $\delta$, i.e., an optimal rectification is realized with $R>0.996$.
Insets in figure 5 plot factor $R$ as a function of $\delta$ for various
temperature differences $\triangle{T}$. And it can be seen that $R$ increases
with the increase of $\triangle{T}$.

Figure 6: $J_{h}$ as a function of $\delta$ with various temperature
difference in the case of not ignoring the inter-system interaction. The black
dashed line corresponds to $T_{1}=T_{2}=10\omega_{0}$. We use blue line for
$T_{1}=10\omega_{0}$ and $T_{2}=5\omega_{0}$; red and green lines for
$T_{2}=10\omega_{0}$ and $T_{1}=5\omega_{0}$ and $0.05\omega_{0}$,
respectively. The other parameters are the same as those in figure 5.

The above rectification phenomena can be explained as follows. There are two
ingredients affecting the rectification effect. One is nonzero heat currents
in the absence of temperature gradient and the other is the temperature
difference $\triangle{T}$ between two environments. (i) It can be seen from
figure 6 that the case of not ignoring the inter-system interaction can give
$J_{h}>0$ even when the two environments are at the same temperature [cf.
black dashed line corresponds to $T_{1}=T_{2}=10\omega_{0}$ in figure 6]. It
seems to violate the second law of thermodynamics but justified. This is
because there is an external work which is closely related to the internal
subsystem-ancilla interaction. By taking this into account, the thermodynamic
consistency can be guaranteed. We also find that the heat transport through
the two-qubit system is dramatically affected by this work cost and can even
be inverted, leading to $J_{h}>0$ for $T_{1}<T_{2}$, i.e., heat current from a
cold to a hot bath [cf. dashed lines in figure 5]. Similarly, in reference [5]
heat transport through a chain of harmonic oscillator is investigated by using
a local master equation based on repeated interactions, and all thermodynamic
inconsistencies (such as a heat current from the cold to the hot bath) can be
resolved correctly with the consideration of external work turning the
subsystem-ancilla coupling on and off. (ii) Because there is a non-zero heat
current in the absence of temperature gradient (labeled $J_{0}$), so when
$T_{1}>T_{2}$, it can be seen from figure 6 that the corresponding heat
current ($J_{h}>J_{0})$ increases with the increase of temperature gradient
[cf. blue line in figure 6]. However, when $T_{1}<T_{2}$, as the temperature
gradient raises, the corresponding heat current decreases from positive to
zero [cf. red line in figure 6]. If the temperature difference increases
further (if it is large enough), the corresponding heat current can change
from positive to negative [cf. green line in figure 6].

Figure 7: Trace distance $\mathcal{D}_{\mathcal{T}}$ between the density
matrices obtained from ignoring and not ignoring the inter-system interaction
against $\delta$ for steady state, and the corresponding plot is obtained for
$T_{1}=10\omega_{0}$ and different $T_{2}$, $T_{2}=0.1\omega_{0}$ (green),
$4\omega_{0}$ (red), and $8\omega_{0}$ (blue).

Figure 7 plots the trace distance $\mathcal{D}_{\mathcal{T}}$ between the
density matrices with and without the approximation against $\delta$ for
steady state. As expected, for $\delta=0$, we observe that
$\mathcal{D}_{\mathcal{T}}=0$, which means that this approximation gives the
same steady state as that of not ignoring the interaction. For $\delta\neq 0$,
$\mathcal{D}_{\mathcal{T}}$ increases with the increase of $\delta$.
Therefore, the steady state from the approximation would deviates from that
without the approximation, and this deviation increases as $\delta$ increases.

Now we consider quantum discord of the bipartite system. Without
approximation, in figure 8 we plot $\mathcal{D}$, i.e., quantum correlations
of two qubits, as a function of $\delta$ at steady state. It is clear that
$\mathcal{D}$ increases with the increase of $\delta$. This behavior can also
be closely related to that of figure 5 (or figure 7), i.e., as $\delta$
increases, the heat currents (or the steady state) predicted by ignoring the
inter-system interaction gradually deviates from that predicted by the case of
not ignoring the inter-system interaction. Therefore this further confirms our
conclusion in section 3.3, i.e., the higher the quantum correlations between
two qubits, the greater the deviation from not ignoring the inter-system
interaction.

Figure 8: Quantum discord $\mathcal{D}$ of the qubit system as a function of
$\delta$, without approximation. The plots are obtained for
$T_{1}=10\omega_{0}$ and different $T_{2}$, $T_{2}=0.1\omega_{0}$ (green),
$4\omega_{0}$ (red), and $8\omega_{0}$ (blue). The remaining parameters are
the same as those in figure 5

In addition, we also consider the other situations by replacing equation (28)
with
$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y})$
or
$\hat{H}_{int}^{S_{1},S_{2}}=\delta(\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z})$.
In spite of a quantitative difference, the corresponding results are similar
to those obtained above.

In section 3 we only consider a resonant two-qubit system (Note that system-
environment interactions that we choose in section 3 is different from
equation (29) in this section). For a more comprehensive study of heat
transfer, now we consider two off-resonant qubits under the same interactions
as in section 3. From our numerical calculations we find that the
corresponding results are similar to those obtained in section 3, and are
different from the results in this section above. That is, as $\delta$
increases, heat current predicted in the case of ignoring the inter-system
interactions gradually deviates from that in the case of not ignoring the
inter-system interactions, then this deviation decreases at larger values of
$\delta$. Moreover, for the system-environment interactions in section 3, we
find that although there is also an asymmetric heat conduction in the case of
not ignoring the inter-system interactions, no obvious thermal rectification
emerges. This shows that thermal rectification are not only related to
asymmetry of the system, but also related to the system-environment
interactions. Here it is noted that in this situation, the external work is
much smaller than the heat flowing through the system in each round in both
cases with and without ignoring the inter-system interactions. This is
different from the results above that the external work is almost the same
order as the heat flowing through the system in each round.

Figure 9: $J_{h}$ as a function of inter-system coupling strength $\delta$.
The solid and dashed lines correspond to not ignoring and ignoring the inter-
system interaction, respectively. We use red line for $T_{1}=10\omega_{0}$ and
$T_{2}=\omega_{0}$. Blue line is for $T_{1}=\omega_{0}$ and
$T_{2}=10\omega_{0}$ in the case of not ignoring the inter-system interaction.
Here we set $\gamma=0.5\omega_{0}$. The other parameters are the same as those
in figure 2. Inset shows corresponding rectification factor in the case of not
ignoring the inter-system interaction.

###  4.2 Anisotropic interacting qubits

Connecting two qubits $S_{1}$ and $S_{2}$ by anisotropic interactions is
another possible way to introduce asymmetry which may also exhibit
rectification. Here we consider the case where the inter-system coupling is
taken as an anisotropic exchange interaction [41]

$\hat{H}_{int}^{S_{1},S_{2}}=\delta\hat{\sigma}_{z}\otimes\hat{\sigma}_{x},$
(31)

and the subsystem-ancilla coupling is chosen as equation (29). As mentioned
before, for this coupling an external work in almost the same order as the
heat flowing through the system in each round with and without ignoring the
inter-system interactions. Thus, the stationary heat current $J_{h}$ flowing
out the subenvironment $E_{i}$ should also be obtained by equation (23). In
figure 9, we plot the heat current as a function of $\delta$ for these two
cases. It can be seen from figure 9 that the heat current increases with the
increase of $\delta$. When $T_{1}>T_{2}$, the heat current obtained from the
approximation (i.e., ignore the inter-system interaction) [cf. red dashed line
in figure 9] gradually deviates from the corresponding one obtained from the
case of not ignoring the inter-system interaction [cf. red solid line in
figure 9] with the increase of $\delta$. We also investigate the corresponding
trace distance $\mathcal{D}_{\mathcal{T}}$, and despite quantitative
difference, its behaviors are similar to those of figure 7, more specifically,
as $\delta$ increases, the steady state predicted in the case of ignoring the
inter-system interaction gradually deviates from that predicted in the case of
not ignoring the inter-system interaction. So this approximation gets
progressively worse as $\delta$ increases.

Moreover, in figure 9, when $T_{1}<T_{2}$ we also consider the heat current in
the case of not ignoring the inter-system interaction. It can be seen that in
the case $T_{1}=10\omega_{0}$ and $T_{2}=\omega_{0}$ [cf. red solid line in
figure 9], its corresponding forward heat current is greater than the backward
heat currents ($T_{2}=10\omega_{0}$ and $T_{1}=\omega_{0}$) [cf. blue solid
line in figure 9]. So there is also a heat rectification, and the
corresponding rectification factor is plotted in the inset of figure 9, which
decreases with the increase of $\delta$. It can be seen that an optimal
rectification is almost realized with a high rectification factor $R>0.985$.
Besides, under the same conventional energy-preserving system-environment
interactions as those of section 3, we also compare the forward and backward
heat currents for this asymmetric system. The corresponding result shows that,
despite an asymmetric heat conduction, there is no obvious thermal
rectification. With lots of numerical calculations, we find that as long as
the system-environment interactions is chosen to be like those in section 3,
no matter how we introduce asymmetry into the system, there is no apparent
thermal rectification effect. Therefore, we conclude that the rectification
effect strongly depends on the form of the system-environment interaction. It
might be the reason why most work [41, 42, 43, 44] chose system-bath
interactions like equation (29) to achieve the thermal rectification for
various kinds of asymmetric systems.

Figure 10: Quantum discord $\mathcal{D}$ of the qubit system as a function of
$\delta$, without approximation. The plot is obtained for $T_{1}=10\omega_{0}$
and $T_{2}=\omega_{0}$. The remaining parameters are the same as those in
figure 9

Next we move to the quantum correlation $\mathcal{D}$ between two subsystems.
Without the approximation, in figure 10 we plot $\mathcal{D}$ as a function of
$\delta$ at steady state. It is clear that its behavior can be closely related
to the failure of the approximation. We again confirm the same conclusion: the
weaker the correlation between the two subsystems, the more justified the
approximation (i.e., ignore the inter-system interaction) and vice versa.
Moreover we do a lot of numerical calculations for various circumstances and
all the results confirm our conclusion: the greater the correlation, the worse
the approximation and vice versa.

##  5 Conclusions

In this paper we have considered a bipartite system consisting of two
identical qubits each coupled to its own heat bath. Based on the collision
model, we mainly study whether the approximation (the inter-system interaction
is ignored when modeling the system-environment coupling) is valid or not for
describing the nonequilibrium dynamics.

We have first investigated the stationary heat current and steady state for
symmetric systems characterized by conventional energy-preserving system-
ancilla interactions. In this situation, no external work exists in the case
of ignoring the inter-system interaction, while an external but small work is
required to turn the system-ancilla coupling on and off in the case of not
ignoring the inter-system interaction. Surprisingly, we have found that the
approximation is still valid to describe either heat current or steady state
even at strong inter-system couplings.

Then we have turned to the asymmetric systems in the case of non-energy
preserving system-ancilla interactions. Here an external work related to the
system-ancilla couplings is almost the same order as the heat flowing through
the system in each round. In this case we have found that the heat current or
(steady state) predicted from the approximation gradually deviates from that
without the approximation as $\delta$ increases, namely, this approximation
gets progressively worse with the increase of $\delta$. In particular, without
the approximation we have realized a perfect thermal rectification effect in
this asymmetric systems and further explained why this phenomenon can happen.
Through the analysis of different forms of interaction between system and
environment, we have found that the thermal rectification effect is closely
related to the form of system-ancilla interaction.

We have also considered the quantum discord between the two qubits without the
approximation. For symmetric systems, we have found that quantum correlation
between two subsystems first increases and then decreases as the inter-system
coupling strength increases. And the deviation of heat current (or steady
state) from the approximation also increases and then decreases in this case.
For asymmetric systems, the quantum correlation between two subsystems
increases with the increase of inter-system coupling strength. And the
deviation of heat current (or steady state) from the approximation in this
case increases with the increase of inter-system coupling strength. Although
one would think that whether this approximation is valid or not should depend
on inter-system coupling strength, through a lot of numerical calculations we
have found that it is closely related to the quantum correlations between the
two subsystems. Specifically, the higher the quantum correlations means these
two qubits behave more like a composite and must be treated as a whole, so
this approximation gets worse, and vice versa.

## Acknowledgement

This work is supported by the National Natural Science Foundation of China
(Grant Nos. 11775019 and 11875086), and Special Funds for Theoretical Physics
of the National Natural Science Foundation of China (Grant No. 11947047).

## References

## References

  * [1] Weiss U 1993 _Quantum Dissipative Systems_ (World Scientific, Singapore)
  * [2] Breuer H P and Petruccione F 2006 _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford).
  * [3] Levy A and Kosloff R 2014 Europhys. Lett. 107 20004
  * [4] Werlang T and Valente D 2015 Phys. Rev. E 91 012143
  * [5] De Chiara G, Landi G, Hewgill A, Reid B, Ferraro A, Roncaglia A J and Antezza M 2018 New J. Phys. 20 113024
  * [6] Deçordi G L and Vidiella-Barranco A 2017 Opt. Commun. 387 366
  * [7] Tahir Naseem M, Xuereb A and Müstecaplıoğlu özgür E 2018 Phys. Rev. A 98 052123
  * [8] Čápek V 2002 Eur. Phys. J. B 25 101
  * [9] Novotný T 2002 Europhys. Lett. 59 648
  * [10] Betzholz R, Torres J M and Bienert M 2014 Phys. Rev. A 90 063818
  * [11] Plenio M B, Almeida J and Huelga S F 2013 J. Chem. Phys. 139 235102
  * [12] Atalaya J and Gorelik L Y 2012 Phys. Rev. B 85 245309
  * [13] Mohseni M, Rebentrost P, Lloyd S and Aspuru-Guzik A 2008 J. Chem. Phys. 129 174106
  * [14] Wilson-Rae I, Zoller P and Imamoḡlu A 2004 Phys. Rev. Lett. 92 075507
  * [15] Scali S, Anders J and Correa L A 2020 arXiv:2009.11324v1
  * [16] Gonzlez J O, Correa L A, Nocerino G, Palao J P, Alonso D and Adesso G 2017 Open Syst. Inf. Dyn. 24 1740010
  * [17] Du J Y and Zhang F L 2018 Commun. Theor. Phys. 70 038
  * [18] Manrique P D, Rodri-guez F, Quiroga L and Johnson N F 2015 Adv. Condens. Matter Phys. 19 615727
  * [19] Hofer P P, Perarnau-Llobet M, Miranda L D M, Haack G, Silva R, Brask J B and Brunner N 2017 New J. Phys. 19 123037
  * [20] Trushechkin A S and Volovich I V 2016 Europhys. Lett. 113 30005
  * [21] Wichterich H, Henrich M J, Breuer H P, Gemmer J and Michel M 2007 Phys. Rev. E 76 031115
  * [22] Benatti F, Floreanini R and Memarzadeh L 2020 arXiv:2004.10433v3
  * [23] Barra F 2015 Sci. rep. 5 14873
  * [24] Cattaneo M, Giorgi G L, Maniscalco S and Zambrini R 2019 New J. Phys. 21 113045
  * [25] Farina D, Filippis G D, Cataudella V, Polini M and Giovannetti V 2020 arXiv:2003.11845v1
  * [26] Wang J S, Wang J and Lü J T 2008 Eur. J. Phys. B 62 381
  * [27] Yan Y, Wu C Q and Li B 2009 Phys. Rev.B 79 014207
  * [28] Thingna J, García-Palacios J L and Wang J S 2012 85 195452
  * [29] Michel M, Hess O, Wichterich H and Gemmer J 2008 Phys. Rev. B 77 104303
  * [30] Manzano D, Tiersch M, Asadian A and Briegel H J 2012 Phys. Rev. E 86 061118
  * [31] Li H, Agarwalla B K, Li B and Wang J S 2013 Eur. Phys. J B 86 500
  * [32] Segal D 2014 Phys. Rev. E 90 012148
  * [33] Bandyopadhyay M and Segal D 2011 Phys. Rev. E 84 011151
  * [34] He D, Thingna J, Wang J S and Li B 2016 Phys. Rev. B 94 155411
  * [35] Man Z X, An N B and Xia Y J 2016 Phys. Rev. E 94 042135
  * [36] Guo B Q, Liu T and Yu C S 2018 Phys. Rev. E 98 022118
  * [37] Wang L and Li B 2008 Phys. Rev. Lett. 101 267203
  * [38] Ruokola T, Ojanen T and Jauho A P 2016 Phys. Rev. B 79 144306
  * [39] Landi G T, Novais E, Oliveira M J and Karevski D 2014 Phys. Rev. E 90 042142
  * [40] Pons M, Cuı Y Y, Ruschhaupt A, Simón M A and Muga J G 2017 Europhys. Lett. 119 64001
  * [41] Kargi C, Naseem M T, Opatrný T, Müstecaplıoğlu Ö E and G. Kurizki 2019 Phys. Rev. E 99 042121
  * [42] Zhang L, Yan Y H, Wu C Q, Wang J S and Li B 2016 Phys. Rev. B 80 172301
  * [43] Werlang T, Marchiori M A, Cornelio M F and Valente D 2014 Phys. Rev. E 89 062109
  * [44] Ordonez-Miranda J, Ezzahri Y and Joulain K 2017 Phys. Rev. E 95 022128
  * [45] Rau J 1963 Phys. Rev. 129 1880
  * [46] Bernardes N K, Carvalho A R R, Monken C H and Santos M F 2014 Phys. Rev. A 90 032111
  * [47] Lorenzo S, McCloskey R, Ciccarello F, Paternostro M and Palma G M 2015 Phys. Rev. Lett. 115 120403
  * [48] Karevski D and Platini T 2009 Phys. Rev. Lett. 102 207207
  * [49] Bernardes N K, Carvalho A R R, Monken C H and SantosW M F 2017 Phys. Rev. A 95 032117
  * [50] Benenti G and Palma G M 2007 Phys. Rev. A 75 052110
  * [51] Filippov S N, Piilo J, Maniscalco S and Ziman M 2017 Phys. Rev. A 96 032111
  * [52] Giovannetti V and Palma G M 2012 Phys. Rev. Lett. 108 040401
  * [53] Daryanoosh S, Baragiola B Q, Guf T and Gilchrist A 2018 Phys. Rev. A 98 062104
  * [54] Man Z X, Xia Y J and Franco R L 2019 Phys. Rev. A 99 042106
  * [55] Manatuly A, Niedenzu W, Ancheyta R R, Çakmak B, Müstecaplıoğlu Ö E and Kurizki G 2019 Phys. Rev. E 99 042145
  * [56] Rodrigues F L S, De ChiaraPhys G, Paternostro M and LandiRev G T 2019 Phys. Rev. Lett. 123 140601
  * [57] Çakmak B, Campbell SVacchini B, Müstecaplıoğlu Ö E and Paternostro M 2018 Phys. Rev. A 99 012319
  * [58] Man Z X, Xia Y J and Franco R L 2018 Phys. Rev. A 97 062104
  * [59] McCloskey R and Paternostro M 2014 Phys. Rev. A 89 052120
  * [60] Campbell S, Ciccarello F, Palma G M and Vacchini B 2018 Phys. Rev. A 98 012142
  * [61] Çakmak B, Pezzutto M, Paternostro M and Müstecaplıoğlu Ö E 2017 Phys. Rev. A 96 022109
  * [62] Beyer K, Luoma K and Strunz W T 2017 Phys. Rev. A 97 032113
  * [63] Lorenzo S, Farace A, Ciccarello F, Palma G M and Giovannetti V 2015 Phys. Rev. A 91 022121
  * [64] Li L, Zou J, Li H, Xu B M, Wang Y M, and Shao B 2018 Phys. Rev. E 97 022111
  * [65] Strasberg P, Schaller G, Brandes T and Esposito M 2017 Phys. Rev. X 7 021003
  * [66] Seah S, Nimmrichter S, and Scarani V 2019 Phys. Rev. E 99 042103
  * [67] Weimear H, Henrich M J, Rempp F, Schröder H and Mahler G 2008 Europhys. Lett. 83 30008
  * [68] Hossein-Nejad H, O’Reilly E J, and Olaya-Castro A 2015 New J. Phys. 17 075014
  * [69] Wang C Z, Li C X, Nie L Y and Li J F 2011 J. Phys. B 44 015503
  * [70] Law C K 2005 Phys. Rev. A 71 034306
  * [71] Lasmar Z, Sajna A S, Lee S Y and Kurzyński P 2018 Phys. Rev. A 98 062105


constitute a generalization, and even more, a simple and natural
generalization of Newton’s laws.

What is also clear is that Hilbert discussed GR in a superficial way only,
concentrating on the mathematical structure of the equations and on their
Lagrangian formulation, but probably without understanding in depth their
physical meaning, quite the opposite of Einstein’s approach to this issue.

  2. 2.

Astronomers who made most important contributions to understand the expansion
of our Univers were Vesto Slipher, Henrietta Leavitt and Edwin Hubble.
However, the first person who clearly realized that the Universe is expanding,
was a Belgian priest, mathematician and physicist, named Georges Lemaître, who
published his results in 1927. We should note that, in those years everybody
believed that the universe was static, for very strong physical reasons.
Indeed, as any ordinary physical system that had had more than enough time to
evolve (infinite, in theory, since the Universe was considered to have always
existed), it should have necessarily reached the stationary state. It could
not be otherwise. But Lemaître proved quite the contrary! on the basis of
Einstein’s GR and matching the theory with the astronomical observations of
Slipher and Hubble, in a masterful way, he proved that the Universe was
expanding; and later, he observed that it was not eternal, that it had an
origin.

During many decades (even now, it is still so declared in most places) people
believed that it was the astronomer Edwin Hubble the person who first
discovered the expansion of the Universe. Only recently, without denying at
all Hubble’s important contributions, historians have put things in the right
place. A very detailed account of this thrilling story can be found in a book
recently published by one of the authors eliz21 .

  3. 3.

In his calculations, using his table of distances (obtained in part with the
help of Leavitt’s law), and Slipher’s table of velocities (obtained as optical
Doppler shifts), Hubble got a rather large value for the expansion rate, of
$H_{0}\cong 500$ $\frac{{\mbox{km/s}}}{{\mbox{Mpc}}},$ which is off the
presently accepted value by almost one order of magnitude. The reason is that
it is extremely difficult to measure cosmological distances. On the contrary,
obtaining velocities by means of the Doppler shifts is somehow easier.
However, here there is also the problem of appropriately disentangling the
recession redshift from other contributions to the observed Doppler shift,
coming from the gravitational influence of other massive celestial objects.
Even today, there is still a sharp controversy about the right value of
$H_{0}$ (see, e.g., ekom20  and references therein). Results from astronomical
observations by different groups, each one reporting uncertainties of just 1
or 2% differ by some 5 to 10%, an unpleasant situation termed the ”Hubble
parameter tension” tension1 .

  4. 4.

In 1922, Alexander Friedmann was the first to discover full families of
solutions of the EEs, which he rightly interpreted as corresponding to
expanding and to contracting universes. In 1922 and 1924 Friedman published
two seminal papers in the prestigious German journal Zeitschrift für Physik
friedm12 , of which Einstein was an Editor. In those papers, Friedmann showed
that there were solutions to Einstein’s equations where the universe evolved
in an expanding or contracting way. Recall (see above) that, in this epoch,
everybody believed that our Universe was static. However, Friedmann explicitly
declared in 1924 that, based on some of his solutions, our Universe might well
be expanding.

Einstein was the ’referee’ (this figure did not actually have the same meaning
and function at that time) of Friedmann’s first paper and, after having
studied it, he wrote a letter to the journal saying that Friedmann’s
calculations should not be published since they contained an error. When
Friedmann (indirectly) learned of this opinion, he sent all the details of his
calculation to Einstein, asking him to check them for himself. After some
discussions (described in detail in e.g. eliz21 ) Einstein finally recognized
that Friedmann did no mistake, and hurried to publish another letter
recognizing his own error and saying that Friedmann’s paper should be
published. Unfortunately, in those days it took a long time to see an article
published after it was finished, and Friedmann died in 1925 before the
publication of Einsteins’ retraction.

  5. 5.

During a break at the very famous Solvay meeting of 1927 in Brussels, Einstein
said literarilly to Lemaître (as was reported later by the last): ”Vos calculs
sont corrects, mais votre physique est abominable´´. He was referring to
Lemaître’s paper of the same year 1927 lemai27 , he had handed to Einstein
during a previous conversation. In the paper (in French) Lemaître had obtained
for the first time ever the Hubble law and, moreover, he had interpreted it in
the right way, as being a proof of the expansion of the Universe, of the very
fabric of the cosmos (an interpretation that Hubble never admitted, in his
whole life). Einstein himself did not accept the expansion of the cosmos until
1932, when he was finally convinced by Richard Tolmann and Willem de Sitter
eliz21 .

  6. 6.

In 1956, George Gamow wrote in Scientific American, that Einstein had told
him, long ago, that the idea of the cosmic repulsion associated with the
cosmological constant had been ”the greatest blunder of his life” (“Die
grösste Eselei meines Lebens,” in German). For years to come, this was the
only testimony of such a claim; what led many to question it, because of
Gamow’s well-known imaginative character. Recently, however, it has been
discovered by historians of physics (oraife1  and references therein) that
Einstein made a similar statement on at least two more occasions. Indeed, John
Wheeler wrote in his book, ”Exploring Black Holes: Introduction to General
Relativity”, that he had personally been present when Einstein said the above
words to Gamow, outside the hall of the Institute for Advanced Studies in
Princeton. Moreover, Ralph Alpher also testified once that he had heard
Einstein make such claim.

In addition, it is a proven fact that Einstein never wanted to use the CC
again, not even when someone suggested that it might be interesting to put it
back, so to better adjust the age of the Universe to the results of
observations of the oldest galaxies, which seemed at one point (erroneously),
to clearly exceed the age of the Universe. As Einstein explained in a footnote
in the appendix to the second edition of his book “The meaning of relativity”
94 : “If Hubble’s expansion had been discovered at the time of the creation of
the general theory of relativity, the cosmological constant would never had
been added. It now seems much less justified to add a term like this in the
field equations, since its introduction loses the only justification it
originally had.”

Finally, an important consideration, which very few mention, is the following
eliz21 : in taking this position, Einstein was even more radical than
Friedmann and Lemaître (the defendants of the expanding universe!), since
those always included the CC term in their models for the Universe; even if
such term was not necessary at all in their equations, contrary to the case of
Einstein’s static model, where it was crucial (for an expanding universe
solution, such an additional term plays a secondary role). Anyhow, if there
was no reason for its presence, it should not be put there, under any
circumstance; this is what Einstein said (please, see eliz21  for additional
details).

  7. 7.

Fred Hoyle an English nuclear physicists and astronomer, who formulated the
theory of stellar nucleosynthesis (and, with it, the remarkable fact that we
are all stardust), was one of the authors of the Steady State Theory of the
Universe, an attempt to maintain a static model for the cosmos able to account
for Hubble’s empirical law of expansion eliz21 . Hoyle did not buy Lemaître’s
conclusion that the universe had an origin, much less his hypothesis of the
primeval atom that latter exploded; which he, as a serious nuclear physicists,
understood it lacked any physical rigor. On the BBC radio’s Third Program
broadcast on 28 March 1949, Hoyle explained to the audience that, when
comparing Lemaître’s model (by then improved by Gamow) with his steady state
theory (where a smooth creation of matter had to take place in order to
compensate for Hubble’s expansion and keep the matter density of the Universe
constant), in Lemaître’s model a sudden creation of all the matter in the
universe had to occur at the very beginning of it. And, for this to happen, an
unbelievably huge expansion (a Big Bang) of the fabric of space was absolutely
necessary. But, of course, such phenomenon was fully impossible! and,
therefore, he pronounced these famous words in a very disdainful tone (see a
much more detailed explanation in eliz21 ). Hoyle anticipated the idea of
cosmic inflation very clearly, albeit as an impossible thought, exactly thirty
years before Alan Guth, on an inspired night, could formulate it precisely.

But even if Hoyle had spoken these two words, Big Bang, in a disrespectful
manner –trying on purpose to mock Lemaître’s model (which had, by then,
substantially been improved by George Gamow)– from this moment on everybody,
starting with Gamow himself (a very peculiar character, as is well known)
began to use this term to refer to the origin of the Universe.

  8. 8.

The cosmological horizon problem (aka the homogeneity problem) is a fine-
tuning issue that affects classical Big Bang models of the universe. It arises
due to the impossibility of explaining the homogeneity reported by
astronomical surveys of very distant regions of space –which are causally
disconnected in these Big Bang models– unless one invokes a mechanism that
sets the same initial conditions everywhere with very high accuracy. This
problem was first pointed out by Wolfgang Rindler in 1956. The most commonly
accepted solution is cosmic inflation, as we have discussed here, but an
explanation in terms of a variable speed of light has been proposed, too.

  9. 9.

The flatness problem is another important issue that appeared in the old,
classical Big Bang model of the Universe. It was first mentioned by Robert
Dicke in 1969, in the Jane Lectures he gave for the American Philosophical
Society that year. The total normalized energy density of our present Universe
has been measured to be very close to 1, with very small uncertainty, what
points towards a very flat universe. Any departure from the conditions leading
to this value in the past would had been magnified enormously over cosmic
time. This leads to the conclusion that one would need an unbelievably
accurate fine tuning in the initial conditions of the Universe, with an energy
density that should have been incredibly closer to the critical value at the
very beginning of the Universe. As we have seen, cosmic inflation provides the
solution to this problem on making the Universe extremely flat, to the needed
precision.

##  VII Conclusions

Our main aim in this review was to explain important issues in Modern
Cosmology in a simple and comprehensible way. To start, from Hubble’s law and
by using an ”homogeneous” version of the Einstein-Hilbert action together with
the first law of thermodynamics, we have easily derived the constituent
equations of Cosmology. Then we dealt with their mathematical singularities,
namely the famous Big Bang singularity and some possible future singularities,
such as the Big Rip and others, providing different ways to remove them and
eventually obtain physically sound results.

Next, we have explained about a number of shortcomings that had appeared in
the old Big Bang Cosmology, such as the horizon and the flatness problems,
whose solution was given in terms of the inflationary paradigm introduced by
Alan Guth. We have reviewed inflation in detail and in very understandable
terms. We have focused on the slow-roll regime, explicitly showing its
attractor behavior and its precise relation with the slow-roll parameters. We
have also calculated, step by step, the number of e-folds the Universe must
necessarily expand, in order to overcome all these problems.

Later, we have focused on the study of the current cosmic acceleration via the
introduction of different sorts of dark energy. Specifically, we have first
considered the most simple model, which uses the cosmological constant. Here
we have shown, as a warm up exercise, that the original approach by Einstein
–namely his static model of the Universe– was unstable (in fact, we show in
detail that it corresponds to a saddle point of the model, viewed as a
dynamical system). Then, we have considered a quintessence field, which can
also be used to unify both periods of inflation, under the form of a most
popular theory, named quintessential inflation, which we have discussed in
some detail.

We have also analysed the Universe’s reheating mechanism, which is a very
important and necessary stage after inflation ends; by analysing again a basic
system: the quantum harmonic oscillator with a time dependent frequency.
Within this simple example, we have introduced the standard diagonalization
method, based on the calculation of the time-dependent Bogoliubov
coefficients. As an application, we have derived the bounds of the reheating
temperature for the model of an $\alpha$-attractor in the context of
quintessential inflation.

To finish, in a closing section we have added a few short notes, which provide
updated descriptions of a number of important historical events. A much more
detailed account of them is given in eliz21 , which is a perfect complement to
the present, more technical albeit still very pedagogical review.

## Acknowledgments

JdH is supported by grant MTM2017-84214-C2-1-P funded by
MCIN/AEI/10.13039/501100011033 and by ”ERDF A way of making Europe”. EE is
supported by MICINN (Spain), project PID2019-104397GB-I00, financed by the
State Research Agency, and by the program Unidad de Excelencia María de Maeztu
CEX2020-001058-M. Both authors are also supported in part by the Catalan
Government, AGAUR project 2017-SGR-247.

## References

  * (1) E. Elizalde, The True Story of Modern Cosmology: Origins, Protagonists and Breakthroughs (Springer Verlag, Berlin, 2021); ISBN: 978-3-030-80653-8. E. Elizalde, Zeta functions and the cosmos: A basic brief review, Universe 7, 5 (2021). E. Elizalde, Some issues on the foundations of Modern Cosmology, Gravitation and Quantum Physics, Universe 6, 189 (2020). E. Elizalde, Reasons in favor of a Hubble-Lemaître-Slipher’s (HLS) law, Symmetry 11, 35 (2019).
  * (2) A. Einstein, Kosmologische Betrachtungen zur allgemeinen Relativitätstheorie. Sitz. König. Preuss. Akad. 142-152 (1917). Translation: Cosmological considerations in the general theory of relativity in Lorentz et al., The Principle of Relativity, Dover 1952, 175-188.
  * (3) W. de Sitter, On Einstein’s theory of gravitation and its astronomical consequences, Mon. Not. Roy. Ast. Soc. 78, 3-28 (1917).
  * (4) A. Friedmann, Über die Krümmung des Raumes, Zeitschrift Phys. 10, 377-386 (1922), On the curvature of space (English translation) https://wwwphy.princeton.edu// steinh/ph563/friedmann.pdf A. Friedmann, Über die Möglichkeit einer Welt mit konstanter negativer Krümmung des Raumes, Zeitschrift Phys. 21, 326-332 (1924).
  * (5) G. Lemaître, Un univers homogène de masse constante et de rayon croissant, rendant compte de la vitesse radiale des nébuleuses extra-galactiques, Annales Soc. Sci. Bruxelles A 47 49-59 (1927); (English translation) Gen. Rel. Grav. 45 8, 1635-1646 (2013).
  * (6) G. Lemaître, The Beginning of the World from the Point of View of Quantum Theory, Nature 127, 706 (1931). https://doi.org/10.1038/127706b0
  * (7) M. V. Fischetti, J. B. Hartle, and B. L. Hu, Quantum Effects in the Early Universe. 1. Influence of Trace Anomalies on Homogeneous, Isotropic, Classical Geometries, Phys. Rev. D20 (1979) 1757–1771.
  * (8) A. Ashtekar and P. Singh, Loop Quantum Cosmology: A Status Report, Class. Quant. Grav. 28, 213001 (2011) [arXiv:1108.0893 [gr-qc]].
  * (9) S. Nojiri and S. D. Odintsov, Introduction to modified gravity and gravitational alternative for dark energy, Proceedings, 42nd Karpacz Winter School, Ladek, Poland, February 6-11, 2006, eConf C0602061 (2006) 06; Int. J. Geom. Meth. Mod. Phys. 4, 115-146 (2007). [arXiv:0601.213 [hep-th]]. T. P. Sotiriou and V. Faraoni, $f(R)$ Theories of Gravity, Rev. Mod. Phys. 82, 451-497 (2010) [arXiv:0805.1726 [gr-qc]].
  * (10) A. Guth, The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems, Phys. Rev. D23, (1981) 347-356.
  * (11) A. Linde, A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy and primordial monopole problems, Phys. Lett. B 108, 389 (1982).
  * (12) A. Albrecht and P.J. Steinhardt, Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking, Phys. Rev. Lett. 48, 1220 (1982).
  * (13) G. Chibisov and V. Mukhanov, Galaxy formation and phonons, Mon. Not. Roy. Astron. Soc. 200, 535 (1982).
  * (14) A. A. Starobinsky, Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations, Phys. Lett. B 117, 175 (1982).
  * (15) A. H. Guth and S-Y. Pi, Fluctuations in the New Inflationary Universe, Phys. Rev. Lett. 49, 1110 (1982).
  * (16) J. M. Bardeen, P. J. Steinhardt and M. S. Turner, Spontaneous creation of almost scale-free density perturbations in an inflationary universe, Phys. Rev. D 28, 679 (1983).
  * (17) A. D. Linde, Scalar Field Fluctuations in Expanding Universe and the New Inflationary Universe Scenario, Phys. Lett. B 116, 335 (1982).
  * (18) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594, A20 (2016) [arXiv:1502.02114 [astro-ph.CO]].
  * (19) A. G. Riess et al, Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant, Astron. J.116, 1009 (1998) [arXiv:astro-ph/9805201].
  * (20) S. Perlmutter et al, Measurements of Omega and Lambda from $42$ High-Redshift Supernovae, Astrophys. J. 517, 565 (1999) [arXiv:astro-ph/9812133].
  * (21) E. J. Copeland, M. Sami and S. Tsujikawa, Dynamics of dark energy, Int. J. Mod. Phys. D 15, 1753 (2006) [arXiv:hep-th/0603057].
  * (22) https://cerncourier.com/a/exploring-the-hubble-tension/; https://en.wikipedia.org/wiki/Hubble
  * (23) S. Tsujikawa, Quintessence: A Review, Class. Quant. Grav. 30, 214003 (2013) [arXiv:1304.1961 [gr-qc]].
  * (24) P. J. E. Peebles and A. Vilenkin, Quintessential Inflation, Phys. Rev. D 59, 063505 (1999) [arXiv:astro-ph/9810509].
  * (25) M. Joyce, Electroweak Baryogenesis and the Expansion Rate of the Universe, Phys. Rev. D 55, 1875 (1997) [arXiv:9606223].
  * (26) L. Parker, Phys. Rev. Lett. 21, 562 (1968). ; L. Parker, Phys. Rev. 183, 1057 (1969); L. Parker, Phys. Rev. D 3, 346 (1970).
  * (27) V. M. Folov, S. G. Mamayev and V. M. Mostepanenko, Phys. Lett. A 55, 7 (1976).
  * (28) A. A. Grib, B. A. Levitskii and V. M. Mostepanenko, Theoreticheskaya i Matematicheskaya Fizika 19, 59 (1974).
  * (29) A. A. Grib, S. G. Mamayev and V. M. Mostepanenko, Gen. Rel. Grav. 7, 535 (1976). A. A. Grib, S. G. Mamayev and V. M. Mostepanenko, Soviet Physics Journal 17, 1700 (1974).
  * (30) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * (31) Ya B. Zeldovich and A. A. Starobinsky, JETP Lett. 26, 252 (1977).
  * (32) B. Spokoiny, Deflationary Universe Scenario, Phys. Lett. B 315, 40 (1993) [arXiv:gr-qc/9306008].
  * (33) S. Hashiba and J. Yokoyama, Gravitational reheating through conformally coupled superheavy scalar particles, JCAP 01, 028 (2019) [arXiv:1809.05410 [gr-qc]].
  * (34) J. Haro, W. Yang and S. Pan, Reheating in quintessential inflation via gravitational production of heavy massive particles: A detailed analysis, JCAP 01, 023 (2019) [arXiv:1811.07371 [gr-qc]].
  * (35) N. D. Birrell and P. C. W. Davies, J. Phys. A: Math. Gen. 13, 2109 (1980)
  * (36) D. J. H. Chung, E.W. Kolb and A. Riotto, Superheavy dark matter, Phys. Rev. D 59, 023501 (1998) [arXiv:hep-ph/9802238].
  * (37) D. J. H. Chung, P. Crotty, E. W. Kolb and A. Riotto, On the gravitational production of superheavy dark matter, Phys. Rev. D 64, 043503 (2001) [arXiv:hep-ph/0104100].
  * (38) G. Felder, L. Kofman and A. Linde, Instant Preheating, Phys. Rev. D 59, 123523 (1999) [arXiv:hep-ph/9812289]
  * (39) G. Felder, L. Kofman and A. Linde, Inflation and Preheating in NO models, Phys.Rev. D60, 103505 (1999) [arXiv:hep-ph/9903350].
  * (40) https://inspirehep.net/files/b039be17085960901edb8c913f3640c7; https://inspirehep.net/files/bd0d1b777e1e505ad92dacf5731aac64; https://pweb.cfa.harvard.edu/big-questions/how-do-stars-and-planets-form-and-evolve.
  * (41) Rodrigo S. Gonçalves, Gabriela C. Carvalho, Uendert Andrade, Carlos A. P. Bengaly, Joel C. Carvalho, Jailson Alcaniz, Measuring the cosmic homogeneity scale with SDSS-IV DR16 Quasars, JCAP 03 (2021) 029, http://curious.astro.cornell.edu/about-us/101-the-universe/cosmology-and-the-big-bang/general-questions/574-what-do-homogeneity-and-isotropy-mean-intermediate; Houjun Mo, Frank van den Bosch, Simon White, Galaxy Formation and Evolution (Cambridge U. Press, 2010).
  * (42) C. Sagan, Cosmos, Chap. 8,  https://www.youtube.com/watch?v=2ktRTHCBBjg
  * (43) E. Hubble, A Relation Between Distance and Radial Velocity Among Extra-Galactic Nebulae. Proceedings of the National academy of Sciences of the United States of America 15: 168–173 (1929).
  * (44) J. P. Huchra, The Hubble constant, https://lweb.cfa.harvard.edu/ dfabricant/huchra/hubble/; J. P. Huchra, Estimates of the Hubble Constant, https://lweb.cfa.harvard.edu/ dfabricant/huchra/hubble.plot.dat
  * (45) E. Elizalde, M. Khurshudyan, S. D. Odintsov and R. Myrzakulov, An analysis of the $H_{0}$ tension problem in a universe with a viscous dark fluid, Physical Review D102, 123501 (2020). arXiv:2006.01879
  * (46) S. Nojiri, S. D. Odintsov, and S. Tsujikawa, Properties of singularities in (phantom) dark energy universe, Phys. Rev. D71 (2005) 063004, [arXiv:hep-th/0501025].
  * (47) J. D. Barrow, Sudden future singularities, Class. Quant. Grav. 21 (2004) L79–L82, [arXiv:gr-qc/0403084].
  * (48) S. Nojiri and S. D. Odintsov, Inhomogeneous equation of state of the universe: Phantom era, future singularity and crossing the phantom barrier, Phys. Rev. D72 (2005) 023003, [arXiv:hep-th/0505215].
  * (49) P. Singh, Loop cosmological dynamics and dualities with Randall-Sundrum braneworlds, Phys. Rev. D73, 063508 (2006) [arXiv:0603043].
  * (50) P. Singh, K. Vandersloot and G. Vereshchagin, Non-singular bouncing universes in loopquantum cosmology, Phys. Rev. D74, 043510 (2006) [arXiv:0606032].
  * (51) M. Sami, P. Singh and S. Tsujikawa, Avoidance of future singularities in loop quantum cosmology, Phys. Rev. D74, 043514 (2006) [arXiv:0605113].
  * (52) J. Amorós, J. de Haro and S. D. Odintsov,  Bouncing Loop Quantum Cosmology from $F(T)$ gravity, Physical Review D 87, 104037 (2013) [arXiv:1305.2344].
  * (53) A. De Felice and S. Tsujikawa, $f(R)$ theories, Living Rev. Rel. 13, 3 (2010) [arXiv:1002.4928 [gr-qc]]. S. Nojiri and S. D. Odintsov, Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models, Phys. Rept. 505, 59-144 (2011). [arXiv:1011.0544 [gr-qc]]. S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution, Phys. Rept. 692, 1-104 (2017). [arXiv:1705.11098 [gr-qc]].
  * (54) P. H. Frampton, K. J. Ludwick, and R. J. Scherrer, The Little Rip, Phys. Rev. D84 (2011) 063003, [arXiv:1106.4996 [astro-ph.CO]].
  * (55) P. H. Frampton, K. J. Ludwick, S. Nojiri, S. D. Odintsov, and R. J. Scherrer, Models for Little Rip Dark Energy, Phys. Lett. B708 (2012) 204–211, [arXiv:1108.0067 [hep-th]].
  * (56) P. C. W. Davies, Singularity Avoidance and Quantum Conformal Anomalies, Phys. Lett. B68 (1977) 402–404.
  * (57) R. M. Wald, Trace Anomaly of a Conformally Invariant Quantum Field in Curved Space-Time, Phys. Rev. D17 (1978) 1477–1484
  * (58) T. Azuma and S. Wada, Classification of Spatially Flat Cosmological Solutions in the Presence of the Cosmological Constant and Backreaction of Conformally Invariant Quantum Fields, Prog. Theor. Phys. 75, 845 (1986).
  * (59) A.A. Starobinsky, A new type of isotropic cosmological models without singularity, Phys. Lett. B91, 99 (1980).
  * (60) L. Kofman, A. Linde and A. Starobinsky, Towards the Theory of Reheating After Inflation, Lev Kofman, Andrei Linde, Alexei Starobinsky Phys. Rev. D56 3258-3295 (1997) [arXiv:hep-ph/9704452].
  * (61) Y. Ema, K. Nakayama and Y. Tang, Production of Purely Gravitational Dark Matter, JHEP 09, 135 (2018) [arXiv:1804.07471 [hep-ph]].
  * (62) S. Hashiba and J. Yokoyama, Gravitational particle creation for dark matter and reheating, Phys. Rev. D99, 043008 (2019) [arXiv:1812.10032 [hep-ph]].
  * (63) T. Damour and A. Vilenkin, String Theory and Inflation, Phys. Rev. D 53, 2981 (1996) [arXiv:9503149].
  * (64) M. Giovannini, Gravitational waves constraints on post-inflationary phases stiffer than radiation, Phys. Rev. D 58, 083504 (1998) [arXiv:9806329].
  * (65) A. R. Liddle and L. A. Urena-Lopez, Curvaton reheating: an application to braneworld inflation, Phys. Rev. D68, 043517 (2003) [arXiv:0302054]
  * (66) B. Feng and M. Li, Curvaton Reheating in Non-oscillatory Inflationary Models, Phys. Lett. B564, 169-174 (2003) [arXiv:0212213].
  * (67) A. Agarwal, S. Bekov and K. Myrzakulov, Quintessential Inflation and curvaton reheating, [arXiv:1807.03629].
  * (68) C. Campuzano, S. del Campo and R. Herrera, Extended Curvaton reheating in inflationary models, JCAP 0606, 017 (2006) [arXiv:astro-ph/0606392].
  * (69) J.C. Bueno Sanchez and K. Dimopoulos, Curvaton reheating allows TeV Hubble scale in NO inflation, JCAP 0711, 007 (2007) [arXiv:0707.3967 [hep-ph]].
  * (70) K. Bamba, E. Elizalde, S. D. Odintsov and T. Paul, Inflationary magnetogenesis with reheating phase from higher curvature coupling, JCAP 04, 009 (2021) [arXiv:2012.12742 [gr-qc]]. D. Maity, S. Pal and T. Paul, Effective theory of inflationary magnetogenesis and constraints on reheating, JCAP 05, 045 (2021) [arXiv:2103.02411 [gr-qc]]. Md. R. Haque, D. Maity, T. Paul and L. Sriramkumarl, Decoding the phases of early and late time reheating through imprints on primordial gravitational waves, Phys. Rev. D 104, 063513 (2021) [arXiv:2105.09242 [astro-ph.CO]]. K. Bamba, S. D. Odintsov, T. Paul and D. Maity, Helical magnetogenesis with reheating phase from higher curvature coupling and baryogenesis, arXiv:2107.11524 [gr-qc].
  * (71) J. de Haro and L. Aresté Saló, A review of Quintessential Inflation, Galaxies 73, 9(4) (2021) [arXiv:2108.11144 [gr-qc]].
  * (72) G. F. Giudice, E. W. Kolb and A. Riotto, Largest temperature of the radiation era and its cosmological implications, Phys. Rev. D 64, 023508 (2001) [arXiv:hep-ph/0005123].
  * (73) J. Ellis, D.V. Nanopoulos and S. Sarkar, Nuc. Phys. B 259, 175 (1985).
  * (74) J. Ellis, A. Linde and D. Nanopoulos, Phys. Lett. B 118, 59 (1982).
  * (75) S. Nojiri and S. D. Odintsov, Non-singular modified gravity: the unification of the inflation, dark energy and dark mater, AIP Conf.Proc.1241:1094-1099, (2010) [arXiv:0910.1464 [hep-th]].
  * (76) S. Nojiri, S. D. Odintsov and V.K. Oikonomou, Unifying Inflation with Early and Late-time Dark Energy in $F(R)$ Gravity, Physics of Dark Universe 29(1):100602, (2020) [arXiv:1912.13128 [gr-qc]].
  * (77) T. Harko, F. S.N. Lobo, S. Nojiri and S, D. Odintsov, $f(R,T)$ gravity, Phys. Rev. D 84, 024020 (2011) [arXiv:1104.2669 [gr-qc]].
  * (78) J. Haro, J. Amorós and S. Pan, The Peebles-Vilenkin quintessential inflation model revisited, Eur.Phys.J. C 79 no.6, 505 (2019) [arXiv:1901.00167 [gr-qc]].
  * (79) K. Dimopoulos,Towards a model of Quintessential Inflation, Nucl. Phys. Proc. Suppl 95, 70 (2001) [arXiv:astro-ph/0012298].
  * (80) Md. Wali Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, A class of quintessential inflation models with parameter space consistent with BICEP2, Phys. Rev. D 89, 123513 (2014) [arXiv:1404.1445 [gr-qc]]
  * (81) Md. Wali Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, Unification of inflation and dark energy à la quintessential inflation, Int. J. Mod. Phys. D 24, 1530014 (2015) [arXiv:1410.6100 [gr-qc]].
  * (82) J. de Haro, J. Amorós and S. Pan, Simple inflationary quintessential model II: Power law potentials, Phys. Rev. D 94, 064060 (2016) [arXiv:1607.06726].
  * (83) C. Q. Geng, C. C. Lee, M. Sami, E. N. Saridakis and A. A. Starobinsky, Observational constraints on successful model of quintessential Inflation, JCAP 06, 011 (2017) [arXiv:1705.01329 [gr-qc]].
  * (84) J. Haro, J. Amorós and S. Pan, Scaling solutions in quintessential inflation, Eur. Phys. J. C 80, no.5, 404 (2020) [arXiv:1908.01516 [gr-qc]].
  * (85) J. de Haro, S. Pan and L. Aresté Saló, Understanding gravitational particle production in quintessential inflation, JCAP 06, 056 (2019) [arXiv:1903.01181 [gr-qc]].
  * (86) J. de Haro, J. Amorós and S. Pan, Simple inflationary quintessential model, Phys. Rev. D 93, no.8, 084018 (2016) [arXiv:1601.08175 [gr-qc]].
  * (87) D. Bettoni and J. Rubio, Quintessential inflation: A tale of emergent and broken symmetries, (2021) [arXiv:2112.11948 [astro-ph.CO]].
  * (88) K. Dimopoulos, J. W. F. Valle, Modeling Quintessential Inflation, Astroparticle Physics 18 3, pages 287-306 (2002) [arXiv:astro-ph/0111417]
  * (89) D. Benisty and E.I. Guendelman, Lorentzian Quintessential Inflation, Int.J.Mod.Phys. D29, no. 14, 2042002 (2020) [arXiv:2004.00339 [astro-ph.CO]].
  * (90) D. Benisty and E.I. Guendelman, Quintessential Inflation from Lorentzian Slow Roll, Eur. Phys. J. C 80, 577 (2020) [arXiv:2006.04129 [astro-ph.CO]].
  * (91) L. Aresté Saló, D. Benisty, E. I. Guendelman and J.d. Haro, Quintessential inflation and cosmological seesaw Mechanism: reheating and observational constraints, JCAP 07, 007 (2021) [arXiv:2102.09514 [astro-ph.CO]].
  * (92) Y. Akrami, R. Kallosh, A. Linde and V. Vardanyan, Dark energy, $\alpha$-attractors, and large-scale structure surveys, JCAP 1806, 041 (2018) [arXiv:1712.09693 [hep-th]].
  * (93) K. Dimopoulos and C. Owen, Quintessential Inflation with $\alpha$-attractors, JCAP 1706, 027 (2017) [arXiv:1703.00305 [gr-qc]].
  * (94) J. Haro and S. Pan, Bulk viscous quintessential inflation, Int. J. Mod. Phys. D 27, 1850052 (2018) [arXiv:1512.03033 [gr-qc]].
  * (95) J.C. Bueno Sanchez and K. Dimopoulos, Trapped Quintessential Inflation, Phys. Lett. B642, 294-301 (2006); Erratum-ibid. B647, 526 (2007) [arXiv:hep-th/0605258].
  * (96) K. Dimopoulos, M. Karciauskas and C. Owen, Quintessential inflation with a trap and axionic dark matter, Phys. Rev. D100, 083530 (2019) [arXiv:1907.04676 [hep-ph]].
  * (97) K. Dimopoulos and S. Sánchez López, Quintessential inflation in Palatini $f(R)$ gravity, Phys. Rev. D103, 043533 (2021) [arXiv:2012.06831 [gr-qc]].
  * (98) L. Aresté Saló, D. Benisty, E. I. Guendelman and J.d. Haro, $\alpha$-attractors in Quintessential Inflation motivated by Supergravity, Phys. Rev. D103, 123535 (2021) [arXiv:2103.07892 [astro-ph.CO]].
  * (99) A. A. Grib, S.G. Mamayev and V. M. Mostepanenko, Friedmann Laboratory Publishing for Theoretical Physics, St. Petersburg (1994).
  * (100) L. Aresté Saló and J. de Haro, Gravitational particle production of superheavy massive particles in Quintessential Inflation: A numerical analysis, Phys. Rev. D104, 083544 (2021) [arXiv:2108.10795 [gr-qc]].
  * (101) T. S. Bunch, Adiabatic regularisation for scalar fields with arbitrary coupling to the scalar curvature, J. Phys. A 13, 1297 (1980).
  * (102) C. O’Raifeartaigh, Investigating the legend of Einstein’s “biggest blunder”, Physics Today, 30 Oct 2018; https://physicstoday.scitation.org/do/10.1063/PT.6.3.20181030a/full/.
  * (103) A. Einstein, The meaning of relativity, Four lectures delivered at Princeton University, May, 1921, EBook-No. 36276 (Project Gutenberg, 2011).