# Structured quantum collision models: generating coherence with thermal
resources

Stefano Cusumano<EMAIL_ADDRESS>International Centre for Theory of
Quantum Technologies, University of Gdańsk, Wita Stwosza 63, 80-308 Gdańsk,
Poland  Università degli studi di Napoli “Federico II”, Dipartimento di Fisica
“E. Pancini”, Via Cupa Cinthia 21, 80126 Napoli, Italy  Gabriele De Chiara
Centre for Quantum Materials and Technology, School of Mathematics and
Physics, Queen’s University Belfast, Belfast BT7 1NN, United Kingdom

###### Abstract

Quantum collision models normally consist of a system interacting with a set
of ancillary units representing the environment. While these ancillary systems
are usually assumed to be either two level systems (TLS) or harmonic
oscillators, in this work we move further and represent each ancillary system
as a structured system, i.e., a system made out of two or more subsystems. We
show how this scenario modifies the kind of master equation that one can
obtain for the evolution of the open systems. Moreover, we are able to
consider a situation where the ancilla state is thermal yet has some
coherence. This allows the generation of coherence in the steady state of the
open system and, thanks to the simplicity of the collision model, this allows
us to better understand the thermodynamic cost of creating coherence in a
system. Specifically, we show that letting the system interact with the
coherent degrees of freedom requires a work cost, leading to the natural
fulfillment of the first and second law of thermodynamics without the
necessity of ad hoc formulations.

##  I Introduction

Open quantum systems are those systems that do not evolve unitarily, i.e.
according to the Schrödinger equation, due to the presence of a usually very
large external system, dubbed the environment, which interacts with the system
influencing its dynamical evolution blum ; breuer .

The main problem when studying such systems is the very large number of
degrees of freedom of the environment, which make it practically impossible to
compute the exact dynamics of the joint system. In order to overcome this
problem, several techniques have been developed in order to take into account
the effects of the environment, both analytical and numerical. Among the
analytical methods we mention master equations methods alicki , stochastic
differential equations gardiner  and path integrals techniques, while among
numerical methods one has the density matrix renormalization group schwollock
; rossini , matrix product states cirac ; orus  and other efficient truncation
techniques to efficiently simulate the dynamics of the open system tamascelli1
; tamascelli2 .

Among these methods, growing attention has been given in recent years to
collision models coll_mod_rev ; review_cusumano ; campbell , which proved to
be useful also for efficient quantum simulations cattaneo . In these models
one depicts the environment as a collection of smaller units, called ancillae,
which interact piecewise with the system giving rise to an open system
dynamics. The main strength of collision models is their versatility in
describing many different physical situations by acting only on few parameters
of the model, like the state of the ancillae, the interaction time and the
interaction Hamiltonian between the ancillae and the system or allowing the
ancillae to interact more than once with the system. Acting on these few knobs
it has been possible to examine different regimes of open system dynamics,
ranging from Markovian gisin ; scarani ; ziman  to non-Markovian dynamics
ciccarello1 ; altamirano ; ciccarello3 ; ciccarello4 , and many different
physical systems, ranging from optical systems cusumano1 ; cusumano2 ;
ciccarello2  to spin chains dechiara2 . Collision models have been used to
study different phenomena such as dissipative dynamics, bath engineering,
quantum darwinism, the difference between local and global Lindblad equations
cattaneo2 .

A very important field where collision models have found application is
quantum thermodynamics strasberg ; cusumano3 ; molitor ; rodrigues ; barra ;
piccione . While in this context one is interested in the influence of quantum
effects on work production and heat dissipation, collision models have proved
useful in investigating the efficiency of quantum heat engines and the role of
quantum coherence in thermodynamic phenomena dechiara ; dechiara3 .

While in most scenarios the environment is represented as a collection of
qubits and independent harmonic oscillators, in recent years some attention
has been given to the case of structured environments, that is, environments
which are not a simple collection of harmonic oscillators, but rather
structured quantum systems, for instance by allowing the harmonic oscillators
to have interaction among them or by considering lattice systems as
environments breuer2 ; strunz ; lukin ; keeling ; vacchini ; manatuly ; dag .
Inspired by this, we introduce a collision model in which the environment is
depicted as a collection of ancillae which are intrinsically composite
systems, that is, for instance, two or more interacting qubits. The crucial
difference with respect to other works in the literature filip ; li ;
purkayastha ; purkayastha2 , is the presence of interactions between the
components of an ancillary system, while in the cited works such interactions
are absent.

These structured ancillae may give rise to driving terms in the master
equation, providing a connection with the field of open quantum systems with
structured environments. We will also use this model to introduce a situation
where, in spite of the fact that the ancillae are in a thermal state, yet they
possess coherence in the energy basis. This represents an advancement with
respect to previous works, where coherence in the ancillae was assumed a
priori, thus hindering the work cost of accessing that coherence and leading
to violations of the second law of thermodynamics and ad hoc formulations in
order to preserve it. Another novelty with respect to previous literature is
that the generation of steady-state coherence is not engineered through a
particular form of the interaction between the system and the bath, but rather
it is due to the presence in the interaction Hamiltonian of operators that are
not eigenoperators of the free Hamiltonian of the ancillae. This also gives
rise explicitly to unitary driving terms in the master equation of the open
system, making it clear from where the coherence stems from.

This can be used to generate quantum coherence in the system’s steady state,
using solely thermal resources. Through this model, we are able to clarify the
role of coherence in thermodynamics, highlighting how the interaction of the
system with an ancilla with coherence requires work in order to be activated.
This makes our work consistent with previous works on steady state coherence
generation filip2  and coherence broadcasting lostaglio ; spekkens ; faist ,
and once more confirms that there is no free lunch.

The manuscript is organized as follows. In Sec. II we introduce our model and
derive the most generic master equation for the case where the ancillae are
structured instead of being simple TLSs or harmonic oscillators, highlighting
the main differences between the two scenarios. Then in Sec. III we show some
features of our model through some examples. In Sec. IV we use our model to
explore the role of coherence in work production and heat dissipation. Finally
in Sec. V we draw our conclusions and give some outlook for future works.

##  II Collision models with structured ancillae

In what follows we want to analyze the dynamics of a quantum system $S$
interacting with a structured environment $\mathcal{E}$ by using a collision
model. Here by structured environment we mean an environment composed of
ancillae which are composite subsystems. As an example, one can think of each
ancilla as being composed by two interacting TLSs or harmonic oscillators.

To represent the collision model, we assume to have an Hamiltonian of the
form:

$\displaystyle\hat{H}_{tot}=\hat{H}_{S}+\hat{H}_{\cal E}+\hat{H}_{S\cal E}$
(1)

where the Hamiltonians $\hat{H}_{\cal E},\hat{H}_{S\cal E}$ are of the form:

$\displaystyle\hat{H}_{\cal E}=\sum_{i}\hat{H}_{E_{i}}\quad\hat{H}_{S\cal
E}=\sum_{i}\hat{H}_{SE_{i}}$ (2)

The Hamiltonians $\\{\hat{H}_{E_{i}}\\}$ are the free Hamiltonians of the
ancillae, while $\hat{H}_{SE_{i}}$ are the interaction Hamiltonians between
the system $S$ and the ancillae. In order to compute the master equation
describing the dynamics of $S$, we start our analysis from the free
Hamiltonian $\hat{H}_{E_{i}}$ of the ancillae. Note that, as we are dealing
with structured ancillae, each free Hamiltonian $\hat{H}_{E_{i}}$ will be in
general a many-body operator. Thus, we start our analysis by studying the most
generic $\hat{H}_{E_{i}}$ in the next subsection.

###  II.1 Structured ancillae

Figure 1: In panel (a) one can see an example of structured ancilla made out
of three subsystems each interacting with each other. Panels (b)-(d):
pictorial representation of the Markovian collision model. In panel (b) the
system interacts with the first ancilla, whose state is factorized from the
one of all the other ancillas. After the system has interacted with the first
ancilla, the latter is traced out and the system interacts with the second
ancilla, as in panel (c). The dynamics, first interacting with an ancilla then
tracing it out, as in panel (d).

In our model, we want to write the Hamiltonian of each ancilla in the form:

$\displaystyle\hat{H}_{E_{i}}$ $\displaystyle=$
$\displaystyle\hat{H}_{E_{i}}^{\text{free}}+\hat{H}_{E_{i}}^{\text{int}},$ (3)
$\displaystyle\hat{H}_{E_{i}}^{\text{free}}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{K}\hbar\omega_{k}\hat{C}_{E_{i}k}^{\dagger}\hat{C}_{E_{i}k},$
(4) $\displaystyle\hat{H}_{E_{i}}^{\text{int}}$ $\displaystyle=$
$\displaystyle\sum_{j\neq
k=1}^{K}\kappa_{jk}\hat{D}_{E_{i}j}\hat{D}_{E_{i}k}^{\dagger}+h.c.$ (5)

Here the operator $\hat{H}_{E_{i}}^{\text{free}}$ includes the free
Hamiltonians of each subsystem composing the ancilla, while the operator
$\hat{H}_{E_{i}}^{\text{int}}$ contains all the interaction terms between the
subsystems of the ancilla.

Note that we are not making any assumption over the nature of the subsystems
of the ancillae, as the operators $\hat{C}_{k}$ can account for both the
bosonic and fermionic case by imposing respectively the canonical commutation
relations (CCR) or the canonical anti-commutation relations (CAR), that is:

$\displaystyle\commutator{\hat{C}_{E_{i}k}}{\hat{C}_{E_{i}j}^{\dagger}}=\delta_{jk},\;\;\commutator{\hat{C}_{E_{i}k}}{\hat{C}_{E_{i}j}}=\commutator{\hat{C}_{E_{i}k}^{\dagger}}{\hat{C}_{E_{i}j}^{\dagger}}=0,$

for the case of bosons and

$\displaystyle\anticommutator{\hat{C}_{E_{i}k}}{\hat{C}_{E_{i}j}^{\dagger}}=\delta_{jk},\;\;\anticommutator{\hat{C}_{E_{i}k}}{\hat{C}_{E_{i}j}}=\anticommutator{\hat{C}_{E_{i}k}^{\dagger}}{\hat{C}_{E_{i}j}^{\dagger}}=0.$

for the case of fermions, where $\commutator{\cdot}{\cdot}$ is the commutator
between the operators, while $\anticommutator{\cdot}{\cdot}$ stands for the
anti-commutator.

Let us now consider the trivial case of unstructured ancillas, that is to say
when $\kappa_{jk}=0\;\forall j,k$. In this case we can easily define the
energy eigenbasis of the ancilla as given by product states of each subsystems
energy eigenstates:

$\displaystyle\ket{n_{1},n_{2},\cdots,n_{K}}=\ket{\vec{n}}=\bigotimes_{k=1}^{K}\ket{n_{k}}.$
(8)

where $\ket{n_{k}}$ is the n-th energy eigenstate of the k-th subsystem.

However, as soon as interactions are present, i.e. for $\kappa_{jk}\neq 0$,
the Hamiltonian $\hat{H}_{E_{i}}$ is not diagonal anymore in the basis
$\\{\ket{\vec{n}}\\}$. This implies that a thermal state of the form
$\exp[-\beta\hat{H}_{E_{i}}]/Z$ is non-diagonal in the basis
$\\{\ket{\vec{n}}\\}$ as well, so that if the interaction between the ancilla
$E_{i}$ and the system $S$ $\hat{H}_{SE_{i}}$ is defined in terms of
eigenoperators of the Hamiltonian $\hat{H}_{E_{i}}^{\text{free}}$, a unitary
driving term may appear in the master equation for $S$.

To continue setting up the notation, we define as $\\{\ket{E^{\prime}_{j}}\\}$
the eigenbasis of the total Hamiltonian $\hat{H}_{E_{i}}$, with corresponding
eigenvalues $E^{\prime}_{j}$. The eigenbasis of $\hat{H}_{E_{i}}$ and
$\hat{H}_{E_{i}}^{\text{free}}$ are related via a unitary transformation
$\hat{V}^{E_{i}}$:

$\displaystyle\ket{E^{\prime}_{k}}=\sum_{\vec{n}}V^{E_{i}}_{k\vec{n}}\ket{\vec{n}}.$
(9)

Practically, the columns of the matrix $\hat{U}^{E_{i}}$ contain the
coefficients of the expansion of the state $\ket{E^{\prime}_{k}}$ in terms of
the non-interacting basis $\\{\ket{\vec{n}}\\}$.

Once we have completed our analysis of the Hamiltonian $\hat{H}_{E_{i}}$ of
the ancillae, we can move to the interaction between $S$ and the ancillae.

###  II.2 System-ancilla interaction and Markovian dynamics

We assume a system-ancilla Hamiltonian of the form:

$\displaystyle\hat{H}_{SE_{i}}=\sum_{j,k}\hat{A}_{S}^{(k,j)}\otimes\hat{B}_{E_{i}k}^{(j)},$
(10)

where the only assumption is that the Hamiltonian $\hat{H}_{SE_{i}}$ is
Hermitian.

We also assume the initial states of $S$ and the ancillae to be factorized:

$\displaystyle\hat{R}_{S\cal
E}(0)=\hat{\rho}_{S}(0)\bigotimes_{i=1}\hat{\eta}_{E_{i}}.$ (11)

More generally, since after the interaction the ancilla is discarded, we
indicate the joint state of $S$ and $\mathcal{E}$ after the collision as:

$\displaystyle\hat{R}_{S\cal
E}(i)=\hat{\rho}_{S}(i)\bigotimes_{j=i+1}\hat{\eta}_{E_{j}}.$ (12)

The joint state of $S$ and $\mathcal{E}$ evolves unitarily according to:

$\displaystyle\hat{R}_{S\cal E}(i)\rightarrow\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger},$ (13)

where the unitary operator $\hat{U}_{SE_{i}}$ is defined as:

$\displaystyle\hat{U}_{SE_{i}}=\exp\left[-\frac{i}{\hbar}\delta
t\left(\hat{H}_{S}+\hat{H}_{E_{i}}+\hat{H}_{SE_{i}}\right)\right].$ (14)

Here $\delta t$ is the collision time during which $S$ and $E_{i}$ interact.

In order to connect the state of the system $S$ before and after its
interaction with an ancilla, we can expand the unitary operator in Eq. (14) in
powers of $\delta t$, obtaining:

$\displaystyle\hat{U}_{SE_{i}}$ $\displaystyle\simeq$
$\displaystyle\hat{1}-\frac{i}{\hbar}\delta
t(\hat{H}_{S}+\hat{H}_{E_{i}}+\hat{H}_{SE_{i}})$ (15)
$\displaystyle-\frac{\delta
t^{2}}{2\hbar^{2}}(\hat{H}_{S}+\hat{H}_{E_{i}}+\hat{H}_{SE_{i}})^{2},$

where $\hat{1}$ is the identity operator. The dynamics now proceeds in
discrete steps. First, the system $S$ interacts with the first ancilla
$E_{1}$. After the interaction the degrees of freedom of $E_{1}$ are traced
away, leaving us with the state of the system after the first collision
$\hat{\rho}_{S}(1)$. Then the system interacts with the second ancilla
$E_{2}$, which is then traced away to leave us with the state
$\hat{\rho}_{S}(2)$ and so on. After the $i$-th collision one has:

$\displaystyle\hat{\rho}_{S}(i+1)=\Tr_{E_{i+1}}\left[\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger}\right].$ (16)

To obtain a master equation for $\hat{\rho}_{S}$, we can use the power
expansion in Eq. (15) and substitute it into Eq. (16). Retaining only terms up
to $\delta t^{2}$ one obtains:

$\displaystyle\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger}=$ (17)
$\displaystyle\mathcal{U}_{0}\left(\hat{R}_{S\cal
E}(i)\right)+\mathcal{U}_{1}\left(\hat{R}_{S\cal
E}(i)\right)+\mathcal{U}_{2}\left(\hat{R}_{S\cal E}(i)\right),$

where the superoperators $\mathcal{U}_{0},\mathcal{U}_{1},\mathcal{U}_{2}$ are
the zero, first and second order contributions to the dynamics respectively.
The explicit form of these contributions reads:

$\displaystyle\mathcal{U}_{0}\left(\hat{R}_{S\cal E}(i)\right)=\hat{R}_{S\cal
E}(i),$ (18) $\displaystyle\mathcal{U}_{1}(\hat{R}_{S\cal E}(i))=$
$\displaystyle\qquad\qquad-\frac{i}{\hbar}\delta
t\commutator{\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}}}{\hat{R}_{S\cal
E}(i)},$ (19) $\displaystyle\mathcal{U}_{2}(\hat{R}_{S\cal E}(i))=\frac{\delta
t^{2}}{2\hbar^{2}}$
$\displaystyle\Big{(}2(\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}})\hat{R}_{S\cal
E}(i)(\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}})$
$\displaystyle-\anticommutator{(\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}})^{2}}{\hat{R}_{S\cal
E}(i)}\Big{)}.$ (20)

It is immediate to observe that the zeroth order contribution corresponds
simply to the identity superoperator, while the first order contribution has
the form of a unitary contribution to the dynamics. Finally, the second order
contribution is the one describing dissipative dynamics.

After computing all these contributions one can trace out the degrees of
freedom of $E_{i+1}$, obtaining:

$\displaystyle\hat{\rho}_{S}(i+1)=$
$\displaystyle\mathcal{L}_{0}(\hat{\rho}_{S}(i))+\mathcal{L}_{1}(\hat{\rho}_{S}(i))+\mathcal{L}_{2}(\hat{\rho}_{S}(i)),$
(21)

where once again the superoperators
$\mathcal{L}_{0},\mathcal{L}_{1},\mathcal{L}_{2}$ represent the zero, first
and second order contribution to the dynamics respectively. Using the
expression of the interaction Hamiltonian in Eq. (10), we can write explicitly
these contributions as:

$\displaystyle\mathcal{L}_{0}(\hat{\rho}_{S}(i))=\hat{\rho}_{S}(i),$ (22)
$\displaystyle\mathcal{L}_{1}(\hat{\rho}_{S}(i))=-\frac{i}{\hbar}\delta
t\commutator{\hat{H}_{S}+\sum_{k,j}g_{k}^{(j)}\hat{A}_{S}^{(k,j)}}{\hat{\rho}_{S}(i)},$
(23) $\displaystyle\mathcal{L}_{2}(\hat{\rho}_{S}(i))=\frac{\delta
t^{2}}{2\hbar^{2}}\Bigg{[}$
$\displaystyle\sum_{k,k^{\prime},j,j^{\prime}}g_{kk^{\prime}}^{(j,j^{\prime})}\left(2\hat{A}_{S}^{(k,j)}\hat{\rho}_{S}(i)\hat{A}_{S}^{(k^{\prime},j^{\prime}){\dagger}}-\anticommutator{\hat{A}_{S}^{(k^{\prime},j^{\prime})}\hat{A}_{S}^{(k,j)}}{\hat{\rho}_{S}(i)}\right)$
$\displaystyle+\sum_{k,j}g_{k}^{(j)}\left(2\hat{H}_{S}\hat{\rho}_{S}\hat{A}_{S}^{(k,j){\dagger}}-\anticommutator{\hat{A}_{S}^{(k,j){\dagger}}\hat{H}_{S}}{\hat{\rho}_{S}(i)}\right)$
$\displaystyle+\sum_{k,j}g_{k}^{(j)}\left(2\hat{A}_{S}^{(j,k)}\hat{\rho}_{S}(i)\hat{H}_{S}-\anticommutator{\hat{H}_{S}\hat{A}_{S}^{(k,j)}}{\hat{\rho}_{S}(i)}\right)$
$\displaystyle+\sum_{k,j}\commutator{g_{kH_{E}}^{(j)}\hat{A}_{S}^{(k,j)}-g_{kH_{E}}^{(j)*}\hat{A}_{S}^{(k,j){\dagger}}}{\hat{\rho}_{S}(i)}$
$\displaystyle+\left(2\hat{H}_{S}\hat{\rho}_{S}(i)\hat{H}_{S}-\anticommutator{\hat{H}_{S}^{2}}{\hat{\rho}_{S}(i)}\right)\Bigg{]},$
(24)

where the factors $g_{k}^{(j)}$, $g_{k}^{(j,j^{\prime})}$ and
$g_{kH_{E}}^{(j)}$ are:

$\displaystyle g_{k}^{(j)}$ $\displaystyle=$
$\displaystyle\Tr\left[\hat{B}_{E_{i}k}^{(j)}\hat{\eta}_{E_{i}}\right],$ (25)
$\displaystyle g_{kk^{\prime}}^{(j,j^{\prime})}$ $\displaystyle=$
$\displaystyle\Tr\left[\hat{B}_{E_{i}k^{\prime}}^{(j^{\prime}){\dagger}}\hat{B}_{E_{i}k}^{(j)}\hat{\eta}_{E_{i}}\right],$
(26) $\displaystyle g_{kH_{E}}^{(j)}$ $\displaystyle=$
$\displaystyle\Tr\left[\hat{H}_{E_{i+1}}\hat{B}_{E_{i+1}k}^{(j)}\hat{\eta}_{E_{i+1}}\right].$
(27)

As expected, the term $\mathcal{L}_{0}$ is simply the identity channel,
mapping the density matrix $\hat{\rho}_{S}$ onto itself. The term
$\mathcal{L}_{1}$ gives instead a unitary contribution, adding a driving term
to the free Hamiltonian $\hat{H}_{S}$. Finally, the term $\mathcal{L}_{2}$
describes dissipation phenomena on $S$. The effect of a structured environment
can be observed by the presence of the double index $kk^{\prime}$ in the
correlation functions $g_{kk^{\prime}}^{(j,j^{\prime})}$, which signal the
presence of terms stemming from the common interaction of $S$ with two
different subsystems $k,k^{\prime}$ of the ancilla $E_{i}$.

Notice also that we did not make the usual assumption $g_{k}^{(j)}=0$. This is
usually a consequence of assuming the environment to be in a thermal state,
for which first order correlation functions are identically null. However,
since we are now considering a structured environment, even if it is in a
thermal state the terms $g_{k}^{(j)}$ could still be different from zero, due
to the coherence between different components of the ancillae.

Let us comment on the continuous limit $\delta t\rightarrow 0$. For systems
where the first order correlation functions of the environment are null, e.g.
ancillae in a thermal state, one can simply introduce a dependence of
$(\sqrt{\delta t})^{-1}$ in the system-environment interaction Hamiltonian,
giving rise to a finite limit coll_mod_rev ; review_cusumano . However, using
the same trick when the first order correlation functions are non-null, as it
could happen in our model, would lead to a divergent equation. This is
sometimes solved by inserting a dependence on $\sqrt{\delta t}$ in the
coherences of the ancilla state, healing the divergence rodrigues ; dechiara ;
hammam . However, in our model this route is not viable, as there is no way of
inserting such an “ad hoc” dependence. Thus in this work we decided to
consider the master equation in Eq. (21) as a discrete equation without
performing the continuous time limit.

This being said, let us now see two examples of a system interacting with a
structured bath, in order to better understand when the presence of structured
ancillae has effects on the dynamics, and which are those effects.

Figure 2: Panels (a): Energy levels of two non interacting TLS. Panel (b):
Energy levels of two TLS interacting through the Hamiltonian in Eq. (28). It
is immediate to see that the energy levels are changed by the interaction
between the two TLS, and that the corresponding eigenstates are linear
combinations of the non-interacting eigenstates according to the coefficients
$a_{\pm},b_{\pm}$, reported in Eq. (79).

##  III Examples of master equation with structured ancillae

###  III.1 Example 1

As first example we want to consider the case where each ancilla is made out
of two interacting TLS, so that the ancilla Hamiltonian reads:

$\displaystyle\hat{H}_{E_{i}}=\frac{\omega_{1}}{2}\hat{\sigma}_{1}^{z}+\frac{\omega_{2}}{2}\hat{\sigma}_{2}^{z}+\kappa_{12}\left(\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{-}+\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{+}\right),$
(28)

where from now on we set $\hbar=1$ for ease of notation. Here
$\hat{\sigma}_{k}^{z}$ is the usual diagonal Pauli matrix, and
$\hat{\sigma}_{k}^{\pm}=\frac{1}{2}(\hat{\sigma}_{k}^{x}\pm
i\hat{\sigma}_{k}^{y})$. The bare energy levels of the two TLS are shown in
the left panel of Fig. 2. To find the eigenenergies of $\hat{H}_{E_{i}}$ shown
in the right panel of Fig. 2 one needs to diagonalize the Hamiltonian in Eq.
(28). To find the eigenenergies one has to solve the equation:

$\displaystyle\lambda^{4}-\lambda^{2}(x^{2}+y^{2}+\kappa_{12}^{2})+x^{2}(y^{2}+\kappa_{12}^{2})=0$
(29) $\displaystyle x=\frac{\omega_{1}+\omega_{2}}{2},\quad
y=\frac{\omega_{1}-\omega_{2}}{2}$ (30)

whose solutions are:

$\displaystyle E^{\prime}_{1}=-x=-\frac{\omega_{1}+\omega_{2}}{2},$ (31)
$\displaystyle
E^{\prime}_{2}=-\sqrt{y^{2}+\kappa_{12}^{2}}=-\frac{1}{2}\sqrt{(\omega_{1}-\omega_{2})^{2}+4\kappa_{12}^{2}},\quad$
(32) $\displaystyle
E^{\prime}_{3}=\sqrt{y^{2}+\kappa_{12}^{2}}=\frac{1}{2}\sqrt{(\omega_{1}-\omega_{2})^{2}+4\kappa_{12}^{2}},$
(33) $\displaystyle E^{\prime}_{4}=+x=\frac{\omega_{1}+\omega_{2}}{2}.$ (34)

The corresponding eigenvectors are:

$\displaystyle\ket{E^{\prime}_{1}}=\ket{E_{1}},$ (35)
$\displaystyle\ket{E^{\prime}_{2}}=\frac{E_{2}-y}{\sqrt{(E_{2}-y)^{2}+\kappa_{12}^{2}}}\left(\frac{\kappa_{12}}{E_{2}-y}\ket{E_{2}}+\ket{E_{3}}\right)\;\;$
(36)
$\displaystyle\ket{E^{\prime}_{3}}=\frac{E_{3}-y}{\sqrt{(E_{3}-y)^{2}+\kappa_{12}^{2}}}\left(\frac{\kappa_{12}}{E_{3}-y}\ket{E_{2}}+\ket{E_{3}}\right)$
(37) $\displaystyle\ket{E^{\prime}_{4}}=\ket{E_{4}}.$ (38)

Note that the weights in $\ket{E^{\prime}_{2}},\ket{E^{\prime}_{3}}$ depend on
the detuning $y$. If the two TLS have the same frequency, i.e.
$\omega_{1}=\omega_{2}$, then $\ket{E^{\prime}_{2}}$ and
$\ket{E^{\prime}_{3}}$ become a singlet and a triplet state with zero
magnetization respectively.

We now assume that the open system $S$ is also a TLS with free Hamiltonian:

$\displaystyle\hat{H}_{S}=\frac{\omega_{S}}{2}\hat{\sigma}_{S}^{z},$ (39)

interacting with the ancillae through the interaction Hamiltonian:

$\displaystyle\hat{H}_{SE_{i}}=\alpha_{1}(\hat{\sigma}_{S}^{+}\hat{\sigma}_{1}^{-}+\hat{\sigma}_{S}^{-}\hat{\sigma}_{1}^{+})+\alpha_{2}(\hat{\sigma}_{S}^{+}\hat{\sigma}_{2}^{-}+\hat{\sigma}_{S}^{-}\hat{\sigma}_{2}^{+}),$
(40)

where in the right hand side (rhs) of the equation we have dropped for
simplicity the index $i$ and the operators $\hat{\sigma}_{1,2}^{\pm}$ are the
same appearing in Eq. (28). Moreover, we assume the ancillae to be initially
in a thermal state at inverse temperature $\beta$, that is:

$\displaystyle\hat{\eta}_{E_{i}}=\frac{e^{-\beta\hat{H}_{E_{i}}}}{Z},\;Z=\Tr[e^{-\beta\hat{H}_{E_{i}}}].$
(41)

We first look at the possible unitary contributions, i.e. the ones stemming
from the terms proportional to $g_{k}^{(j)}$ or $g_{kH_{E}}^{(j)}$. As shown
in App. A all these terms are zero, and thus there is no unitary contribution
to the dynamics in this case.

Next we look at the structure of the dissipative part of the dynamics, the one
described by the term $\mathcal{L}_{2}(\hat{\rho}_{S})$. To do this, we give a
look at the correlation functions $g_{kk^{\prime}}^{j,j^{\prime}}$ of the
environment, starting with the ones where $k=k^{\prime}$. There are four non-
zero such correlation functions:

$\displaystyle
g_{11}^{(+,-)}=\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{1}^{+}\hat{\eta}_{E_{i}}\right],g_{11}^{(-,+)}=\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{1}^{-}\hat{\eta}_{E_{i}}\right],$
(42) $\displaystyle
g_{22}^{(+,-)}=\Tr\left[\hat{\sigma}_{2}^{-}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right],g_{22}^{(-,+)}=\Tr\left[\hat{\sigma}_{2}^{+}\hat{\sigma}_{2}^{-}\hat{\eta}_{E_{i}}\right],$
(43)

which are computed explicitly in App. A. In addition to this, we also have
cross correlation functions, i.e. those with $k\neq k^{\prime}$. In this case
we have eight possible terms, which are complex conjugates in couples:

$\displaystyle
g_{12}^{(+,-)}=g_{21}^{(+,-)*}=\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]$
(44) $\displaystyle
g_{12}^{(-,+)}=g_{21}^{(-,+)*}=\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{-}\hat{\eta}_{E_{i}}\right]$
(45) $\displaystyle
g_{12}^{(-,-)}=g_{21}^{(-,-)*}=\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]$
(46) $\displaystyle
g_{12}^{(+,+)}=g_{21}^{(+,+)*}=\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-}\hat{\eta}_{E_{i}}\right].$
(47)

As shown in App. A, the only non zero elements are the ones in Eqs.( 44, 45).
Conversely, the terms with $g_{12}^{(+,+)}$ and $g_{12}^{(-,-)}$ are equal to
zero. In the end, we can write the master equation for the dynamics of $S$ in
the interaction picture as follows:

$\displaystyle\hat{\rho}_{S}(i+1)=\hat{\rho}_{S}(i)-i\delta
t\commutator{\hat{H}_{S}}{\hat{\rho}_{S}(i)}$ $\displaystyle+\frac{\delta
t^{2}}{2}\Bigg{[}\sum_{i,j=1,2}g_{ij}^{(+,-)}\mathcal{D}[\hat{\sigma}_{S}^{+}](\hat{\rho}_{S}(i))$
$\displaystyle\sum_{i,j=1,2}g_{ij}^{(-,+)}\mathcal{D}[\hat{\sigma}_{S}^{-}](\hat{\rho}_{S}(i))+\frac{\omega_{S}^{2}}{4}\mathcal{D}[\hat{\sigma}_{S}^{z}](\hat{\rho}_{S}(i))\Bigg{]},$
(48)

where we have defined:

$\displaystyle\mathcal{D}[\hat{A}](\cdots)=2\hat{A}\cdots\hat{A}^{\dagger}-\anticommutator{\hat{A}^{\dagger}\hat{A}}{\cdots},$
(49)
$\displaystyle\mathcal{D}^{\prime}[\hat{A},\hat{B}](\cdots)=2\hat{A}\cdots\hat{B}-\anticommutator{\hat{B}\hat{A}}{\cdots}.$
(50)

This is the usual discrete Markovian master equation for a system interacting
with a thermal bath. The presence of coherence in the thermal bath has no
effect on the dynamics and thus the system will simply reach thermal
equilibrium at the steady state, without the formation of any coherence. One
can in fact compute explicitly the steady state satisfying
$\hat{\rho}_{S}(i+1)=\hat{\rho}_{S}(i)$ to obtain:

$\displaystyle\hat{\rho}_{S}^{steady}=\frac{1}{\sum_{i,j,\pm}g_{ij}^{(\pm,\mp)}}\begin{bmatrix}\sum_{i,j}g_{ij}^{(-,+)}&0\\\
0&\sum_{i,j}g_{ij}^{(+,-)}\end{bmatrix}.$ (51)

With this example we can see that when system $S$ interacts with each qubit of
the ancilla independently, then the presence of coherence in the bath degrees
of freedom has no effects on the dynamics. In the next example we are going to
see how, when the interaction between the system and the ancilla involves
joint degrees of freedom of the ancilla, i.e. degrees of freedom that have
coherence, then the presence of coherence becomes significant for the
dynamics.

###  III.2 Example 2

As a second example we consider the same system we examined in Sec. III.1, but
this time with an interaction Hamiltonian between the system and the ancillae
of the form:

$\displaystyle\hat{H}_{SE_{i}}=\alpha\hat{\sigma}_{S}^{+}\hat{\sigma}_{E_{i}}^{-}+\alpha^{*}\hat{\sigma}_{S}^{-}\hat{\sigma}_{E_{i}}^{+}$
(52)

where the operators $\hat{\sigma}_{E_{i}}^{\pm}$ are defined as:

$\displaystyle\hat{\sigma}_{E_{i}}^{+}=\outerproduct{E_{3}}{E_{2}}\quad\hat{\sigma}_{E_{i}}^{-}=\outerproduct{E_{2}}{E_{3}},$
(53)

and the states $\ket{E_{2}},\ket{E_{3}}$ are the bare levels of the ancilla
Hamiltonian. Thanks to the results of the previous section we can immediately
write the Hamiltonian and dissipative parts of the dynamics. We have:

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=$
$\displaystyle-i\commutator{\left(g_{E_{i}}^{(+)}\right)\hat{\sigma}_{S}^{-}+\left(g_{E_{i}}^{(-)}\right)\hat{\sigma}_{S}^{+}}{\hat{\rho}_{S}}$
(54) $\displaystyle=$
$\displaystyle-i\commutator{\Re{g_{E_{i}}^{(+)}}\hat{\sigma}_{S}^{x}+\Im{g_{E_{i}}^{(+)}}\hat{\sigma}_{S}^{y}}{\hat{\rho}_{S}},$

where the constants $g_{E_{i}}^{(\pm)}$ are:

$\displaystyle
g_{E_{i}}^{(+)}=\Tr\left[\hat{\sigma}_{E_{i}}^{+}\hat{\eta}_{E_{i}}\right]=g_{E_{i}}^{(-)*}.$
(55)

Note that for an unstructured environment these terms would be identically
zero, as they are proportional to the amount of coherence between the levels
$\ket{E_{2}}$ and $\ket{E_{3}}$. In practice, the term
$\mathcal{L}_{1}(\hat{\rho}_{S})$ describe the unitary evolution due to an
effective Hamiltonian of the form:

$\displaystyle\hat{H}_{S}^{\rm
eff}=\Re{g_{E_{i}}^{(+)}}\hat{\sigma}_{S}^{x}+\Im{g_{E_{i}}^{(+)}}\hat{\sigma}_{S}^{y}.$

Figure 3: Plots of the evolution of populations and coherence. In all the
plots $\delta t=0.1$, $\omega_{S}=1$, $\omega_{1}=0.5$, $\omega_{2}=1.5$,
$\kappa_{12}=0.3$, $\alpha=0.1$. The initial state is a thermal state of the
form $\hat{\rho}_{S}(0)=\exp[-\beta_{S}\hat{H}_{S}]/Z_{S}$, $Z_{S}$ being the
partition function and $\beta_{S}=1$. Panel (a): plot of the expectation value
$\langle\hat{\sigma}_{z}\rangle$ as a function of the number of collisions. As
one can see, the population tends to equilibrate towards the steady state
value, which is not simply the thermal population, as the steady state also
has coherence. This observation in particular justifies the behavior of the
population for the case $\beta=\beta_{S}$, which without the creation of
coherence should just be a line. In panels (b) the $\ell_{1}$ norm of
coherence $C_{\ell_{1}}(\hat{\rho}_{S})$ is plotted. One can observe how the
amount of coherence is maximal when the system and the environment have the
same temperature. In panel (c) and (d) the same quantities are plotted, this
time for different values of $\kappa_{12}$, the environmental temperature
being $\beta=\beta_{S}=1$. As one can observe, there is a small difference in
the equilibrium populations, due to the fact that the amount of coherence
created is different. Moreover, one can observe that the stronger the
interaction between the two TLSs composing the ancilla, the greater the amount
of coherence created.

While we leave the explicit calculations of the coefficients to App. B, we can
immediately write down the master equation as:

$\displaystyle\hat{\rho}_{S}(i+1)=\hat{\rho}_{S}(i)-i\delta
t\commutator{\hat{H}_{S}+g_{E_{i}}^{(-)}\hat{\sigma}_{S}^{+}+g_{E_{i}}^{(+)}\hat{\sigma}_{S}^{-}}{\hat{\rho}_{S}(i)}$
$\displaystyle+\frac{\delta
t^{2}}{2}\Bigg{[}\frac{\omega_{S}^{2}}{4}\mathcal{D}[\hat{\sigma}_{S}^{z}](\hat{\rho}_{S}(i))+g_{E_{i}}^{(+,-)}\mathcal{D}[\hat{\sigma}_{S}^{-}](\hat{\rho}_{S}(i))$
$\displaystyle+g_{E_{i}}^{(-,+)}\mathcal{D}[\hat{\sigma}_{S}^{-}](\hat{\rho}_{S}(i))+\frac{g_{E_{i}}^{(+)}\omega_{S}}{2}\mathcal{D}^{\prime}[\hat{\sigma}_{S}^{-},\hat{\sigma}_{S}^{z}](\hat{\rho}_{S}(i))$
$\displaystyle+\frac{g_{E_{i}}^{(+)}\omega_{S}}{2}\mathcal{D}^{\prime}[\hat{\sigma}_{S}^{z},\hat{\sigma}_{S}^{-}](\hat{\rho}_{S}(i))+\frac{g_{E_{i}}^{(-)}\omega_{S}}{2}\mathcal{D}^{\prime}[\hat{\sigma}_{S}^{+},\hat{\sigma}_{S}^{z}](\hat{\rho}_{S}(i))$
$\displaystyle\qquad+\frac{g_{E_{i}}^{(-)}\omega_{S}}{2}\mathcal{D}^{\prime}[\hat{\sigma}_{S}^{z},\hat{\sigma}_{S}^{+}](\hat{\rho}_{S}(i))\Bigg{]}.$
(56)

The dynamics described by this master equation consists in the usual
dissipative contribution and dephasing term driving the system towards the
thermal equilibrium, plus a unitary contribution which creates coherence in
the system. On one hand one can see that the populations, i.e. $\rho_{11}$ and
$\rho_{22}$, evolve according to the usual term due to dissipation, plus
another term due to the presence of the unitary driving. On the other hand,
the coherences $\rho_{12}$ and $\rho_{21}$, besides the usual decay term due
to the presence of the environment, evolves also according to terms due to the
presence of the unitary driving, thus leading to a non-zero value in the
steady state.

In fact it is possible to compute the steady state of system $S$ by imposing
$\hat{\rho}(i+1)=\hat{\rho}(i)$, obtaining analytical expressions for
$\langle\hat{\sigma}_{x}\rangle$, $\langle\hat{\sigma}_{x}\rangle$,
$\langle\hat{\sigma}_{z}\rangle$, which are reported in App. B. What can be
seen from these expressions is that the coherence is non null in the steady
state, as both $\langle\hat{\sigma}_{x}\rangle$ and
$\langle\hat{\sigma}_{y}\rangle$ are non-zero. Moreover, it is worth noticing
that the equilibrium populations will not simply be the thermal populations,
due to the fact that coherence is created in the system $S$.

Figure 4: Detail of the evolution of $\langle\hat{\sigma}_{z}\rangle$. Panel
(a): plot of the evolution of $\langle\hat{\sigma}_{z}\rangle$ for different
values of the environmental temperature $\beta$. The values of the other
parameters are $\delta t=0.1$, $\omega_{S}=1$, $\omega_{1}=0.5$,
$\omega_{2}=1.5$, $\kappa_{12}=0.3$, $\alpha=0.2$. It can be observed that,
while going towards its steady state value, $\langle\hat{\sigma}_{z}\rangle$
oscillates due to the presence of the unitary term. Panel (b): evolution of
$\langle\hat{\sigma}_{z}\rangle$ for different values of $\kappa_{12}$, where
the environmental temperature is $\beta=\beta_{S}=1$. From this plot one can
see that the amplitude of the oscillations of $\langle\hat{\sigma}_{z}\rangle$
is larger for larger values of $\kappa_{12}$.

We now want to analyze with more detail the dynamics of $S$, in particular
checking how the populations and the coherence evolve, that is, we want to
look at the transient dynamics before reaching the steady state. In order to
quantify the amount of coherence in the state $\hat{\rho}_{S}$ we use the
$\ell_{1}$ norm of coherence plenio  $C_{\ell_{1}}(\hat{\rho})$, defined as:

$\displaystyle C_{\ell_{1}}(\hat{\rho})=\sum_{i\neq
j}\left|\hat{\rho}_{ij}\right|,$ (57)

which in the present example is simply evaluated as:

$\displaystyle
C_{\ell_{1}}(\hat{\rho}_{S})=\left|\langle\hat{\sigma}_{x}\rangle+i\langle\hat{\sigma}_{y}\rangle\right|.$
(58)

In Fig. 3 some plots for the populations $\langle\hat{\sigma_{z}}\rangle$ and
$C_{\ell_{1}}(\hat{\rho}_{S})$ are shown. In the upper plots the evolution of
these quantities is shown for different values of the environmental
temperature $\beta$, while in the bottom plots the same quantities are plotted
for different values of $\kappa_{12}$.

Some observations are to be made. First, one can notice that the equilibrium
populations are not the thermal ones, consistently with the fact that
coherence is created in system $S$. Furthermore, one can observe oscillations
in the value of $C_{\ell_{1}}(\hat{\rho}_{S})$. These oscillations are slowly
suppressed as the system reaches its steady state. Finally, one can clearly
observe that the total amount of coherence available to system $S$ is strongly
dependent on the amount of coherence in the ancillae, ultimately depending on
the value of $\kappa_{12}$.

Figure 5: Plots of $C_{\ell_{1}}(\hat{\rho}_{S})$ at the steady state. The
values of the parameters are $\delta t=0.1$, $\omega_{S}=1$, $\omega_{1}=0.5$,
$\omega_{2}=1.5$, $\alpha=0.2$ Panel (a): plot of
$C_{\ell_{1}}(\hat{\rho}_{S})$ as a function of $\kappa_{12}$ for different
values of $\beta$. In this case it can be noticed that in general the amount
of coherence generated in the steady state is larger for lower temperatures.
However the behavior of the coherence for large values of $\kappa_{12}$ shows
that the amount of coherence reaches a steady value independently of the
temperature. Moreover, one can see that for specific values of $\kappa_{12}$
the amount of coherence is larger. Panel (b): plot of
$C_{\ell_{1}}(\hat{\rho}_{S})$ as a function of $\beta$ for different values
of $\kappa_{12}$. On one hand, it can be immediately noticed that, for values
of $\kappa_{12}\lesssim 1$ one has larger values of
$C_{\ell_{1}}(\hat{\rho}_{S})$ for larger values of $\kappa_{12}$. On the
other hand, for values of $\kappa_{12}\gtrsim 1$, in the limit of very large
$\beta$, the coherence reaches always the same steady value. The case
$\kappa_{12}\simeq 1$ represents the border between the two scenarios, where
the coherence, in the limit of large $\beta$ reaches the higher value.

Moreover, in Fig. 4 we have plotted a detail of the evolution of the
populations for various temperatures $\beta$ and different $\kappa_{12}$. It
can be observed, consistently with the presence of the unitary driving term,
that while going to their steady state value, the populations also oscillates.
One can also observe that these oscillations become wider for larger values of
$\kappa_{12}$, as the strength of the unitary driving becomes larger.

In Fig. 5 we have plotted the steady state value of
$C_{\ell_{1}}(\hat{\rho}_{S})$ as a function of $\kappa_{12}$ for several
values of the environmental temperature $\beta$ and as a function of $\beta$
for different values of $\kappa_{12}$. As one can see in Fig. 5(a), on one
hand for all values of $\beta$, as $\kappa_{12}$ grows large the steady state
coherence converges towards an asymptotic value. On the other hand, one can
observe a peak in $C_{\ell_{1}}(\hat{\rho}_{S})$ for specific values of
$\kappa_{12}$ and $\beta$. In Fig. 5(b) one can instead observe a different
behavior as a function of the temperature for different values of
$\kappa_{12}$. In fact, for $\kappa_{12}\lesssim 1$ it can be observed that
the steady state coherence goes to zero as $\beta$ grows large. Conversely,
for values of $\kappa_{12}\gtrsim 1$, as $\beta$ grows large one can observe
that the steady state coherence goes towards a non zero value, the largest
value being obtained for $\kappa_{12}\simeq 1$. In conclusion, the dependence
of $C_{\ell_{1}}(\hat{\rho}_{S})$ on $\kappa_{12}$ and $\beta$ is in general
non trivial, though its main characteristics can be extrapolated. These plots
demonstrate the result we anticipated and one of the main contributions of
this paper: the generation of steady-state coherence using only thermal
resources.

##  IV Thermodynamics

In this section we want to analyze the thermodynamics of a system when
interacting with a bath made out of structured ancillae. We will first define
the relevant thermodynamic quantities, such as internal energy, work and heat.
Then we will make some general considerations based on the general form of the
master equation before moving to a practical example based on the system
described in Sec. III.2. We will show that the creation of coherence in the
steady state requires work to be injected into the system, thus recovering the
first principle of thermodynamics in its standard form and fulfilling the
second principle of thermodynamics.

###  IV.1 Defining the thermodynamic quantities

We define the thermodynamic quantities starting from the internal energy of
system $S$, which is simply:

$\displaystyle U=\Tr\left[\hat{H}_{S}\hat{\rho}_{S}\right],$ (59)

so that at each step the internal energy variation can be defined as:

$\displaystyle\Delta U(i+1)$ $\displaystyle=$
$\displaystyle\Tr\left[\hat{H}_{S}\left(\hat{\rho}_{S}(i+1)-\hat{\rho}_{S}(i)\right)\right]$
(60) $\displaystyle=$
$\displaystyle\Tr\left[\hat{H}_{S}\left(\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}-\hat{R}_{S\cal E}(i)\right)\right].\qquad$

By inserting the expressions in Eqs. ( 19, 20) into Eq. (60) we obtain:

$\displaystyle\Delta U(i+1)=\Tr\left[\hat{H}_{S}\mathcal{U}_{1}(\hat{R}_{S\cal
E}(i))\right]+\Tr\left[\hat{H}_{S}\mathcal{U}_{2}(\hat{R}_{S\cal
E}(i))\right]$ $\displaystyle=-i\delta
t\Tr\left[\commutator{\hat{H}_{S}}{\hat{H}_{SE_{i+1}}}\hat{R}_{S\cal
E}(i)\right]+\frac{\delta
t^{2}}{2}\Tr\left[\commutator{\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}}}{\commutator{\hat{H}_{S}}{\hat{H}_{SE_{i+1}}}}\hat{R}_{S\cal
E}(i)\right].$ (61)

Figure 6: Plot of $\Delta U$, $\delta Q$, $W_{sw}$ for the system in Example 2
and the same parameters value s in Fig. 3. Panel (a), (b), (c): plots of the
three quantities for the initial part of the dynamics. As one can see, as both
the internal energy change and the switching work depend on the coherences in
system $S$, one can observe oscillations in these quantities. Conversely, heat
does not oscillate, as it only depends on the populations of $S$ and the
ancilla matrix elements.

We now define the heat exchanged by the system, assumed to be positive when it
is absorbed by the system. At each interaction between the system and an
ancilla, we thus define the heat exchanged as the opposite of the energy
variation of the ancilla:

$\displaystyle\delta
Q(i+1)=-\Tr\left[\hat{H}_{E_{i}}(\hat{\eta}_{E_{i+1}}^{out}-\hat{\eta}_{E_{i+1}}^{in})\right]\quad$
$\displaystyle=-\Tr\left[\hat{H}_{E_{i+1}}\left(\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger}-\hat{R}_{S\cal E}(i)\right)\right],\qquad$
(62)

where we have defined:

$\displaystyle\hat{\eta}_{E_{i+1}}^{out}=\Tr_{S}\left[\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger}\right].$ (63)

Definition (62) is justified since the initial state of the environment is
thermal. Once again we can insert the expressions in Eqs. ( 19, 20) to obtain:

$\displaystyle\delta
Q(i+1)=-\Tr\left[\hat{H}_{E_{i+1}}\mathcal{U}_{1}(\hat{R}_{S\cal
E}(i))\right]+\Tr\left[\hat{H}_{E_{i+1}}\mathcal{U}_{2}(\hat{R}_{S\cal
E})\right]$ $\displaystyle=i\delta
t\Tr\left[\commutator{\hat{H}_{E_{i+1}}}{\hat{H}_{SE_{i+1}}}\hat{R}_{S\cal
E}(i)\right]-\frac{\delta
t^{2}}{2}\Tr\left[\commutator{\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}}}{\commutator{\hat{H}_{E_{i+1}}}{\hat{H}_{SE_{i+1}}}}\hat{R}_{S\cal
E}(i)\right],$

Moreover, any variation of the system Hamiltonian will contribute to work, so
that we define:

$\displaystyle\delta W=\Tr[\Delta\hat{H}_{S}\hat{\rho}_{S}].$ (64)

As a consequence, we also have to take into consideration the switching work
contribution strasberg , namely the one stemming from turning on and off of
the interaction between the system and the ancilla, which reads:

$\displaystyle W_{sw}(i+1)$ $\displaystyle=$
$\displaystyle\Tr\left[\hat{H}_{SE_{i}}\left(\hat{R}_{S\cal
E}(i)-\hat{U}_{SE_{i+1}}\hat{R}_{S\cal
E}(i)\hat{U}_{SE_{i+1}}^{\dagger}\right)\right],$

Also in this case we can use the expressions in Eqs. ( 19, 20) and obtain:

$\displaystyle
W_{sw}(i+1)=-\Tr\left[\hat{H}_{SE_{i+1}}\mathcal{U}_{1}(\hat{R}_{S\cal
E}(i))\right]-\Tr\left[\hat{H}_{SE_{i+1}}\mathcal{U}_{2}(\hat{R}_{S\cal
E}(i))\right]=i\delta
t\Tr\left[\commutator{\hat{H}_{SE_{i+1}}}{\hat{H}_{S}+\hat{H}_{E_{i+1}}}\hat{R}_{S\cal
E}(i)\right]$ (66) $\displaystyle-\frac{\delta
t^{2}}{2}\Tr\left[\commutator{\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}}}{\commutator{\hat{H}_{SE_{i+1}}}{\hat{H}_{S}}}\hat{R}_{S\cal
E}(i)\right]-\frac{\delta
t^{2}}{2}\Tr\left[\commutator{\hat{H}_{S}+\hat{H}_{E_{i+1}}+\hat{H}_{SE_{i+1}}}{\commutator{\hat{H}_{SE_{i+1}}}{\hat{H}_{E_{i+1}}}}\hat{R}_{S\cal
E}(i)\right].$

It is important to note that this term will be different from zero in case the
interaction Hamiltonian is energy non-conserving. Note also that in our
convention the work is positive when it is absorbed by the system and negative
when it is extracted.

Putting together Eqs.( 61, IV.1, 66) we can observe that the first principle
of thermodynamics is automatically fulfilled, as it reads:

$\displaystyle\Delta U=\delta Q+\delta W+W_{sw},$ (67)

where in absence of external driving $\delta W=0$. The interpretation of Eqs.(
61, IV.1, 66) is straightforward as one can actually perform explicitly the
traces, obtaining the expressions reported in App. C.. At each collision, work
is needed in order to turn on and off the interaction between the system and
the ancilla. Part of this work is then dissipated as heat into the bath, while
another part is used to change the internal energy of the system.

Let us now turn our attention to the entropic quantities. The von Neumann
entropy of a quantum state is defined as:

$\displaystyle
S(\hat{\rho}_{S})=-\Tr\left[\hat{\rho}_{S}\log\hat{\rho}_{S}\right].$ (68)

We thus define the system $S$ entropy variation during one collision as:

$\displaystyle\Delta S_{S}(i+1)=S(\hat{\rho}_{S}(i+1))-S(\hat{\rho}_{S}(i)).$
(69)

The entropy production, describing the total amount of entropy created during
one collision, including the contributions due to the change of the ancilla
state and the one due to the loss of information, can be defined as strasberg
:

$\displaystyle\Sigma(i)=\mathcal{I}(\hat{\rho}_{S}(i);\hat{\eta}_{E_{i}}^{out})+S(\hat{\eta}_{E_{i}}^{out}||\hat{\eta}_{E_{i}}^{in}),$
(70)

where

$\displaystyle\mathcal{I}(\hat{\rho}_{S}(i);\hat{\eta}_{E_{i}}^{out})=$
$\displaystyle=S(\hat{\rho}_{S}(i))+S(\hat{\eta}_{E_{i}}^{out})-S(\hat{U}_{SE_{i}}R_{S\cal
E}(i)\hat{U}_{SE_{i}})$
$\displaystyle=S(\hat{\rho}_{S}(i))+S(\hat{\eta}_{E_{i}}^{out})-S(\hat{\rho}_{S}(i-1))-S(\hat{\eta}_{E_{i}}^{in}),$

is the mutual information between system $S$ and the ancilla $E_{i}$ after
their interaction, while:

$\displaystyle
S(\hat{\eta}_{E_{i}}^{out}||\hat{\eta}_{E_{i}}^{in})=\Tr\left[\hat{\eta}_{E_{i}}^{out}\log\hat{\eta}_{E_{i}}^{out}\right]-\Tr\left[\hat{\eta}_{E_{i}}^{out}\log\hat{\eta}_{E_{i}}^{in}\right],$
(71)

is the relative entropy between the state of the ancilla after the interaction
and the one before the interaction. Summing the two quantities together one
obtains:

$\displaystyle\Sigma(i)=\Delta S_{S}(i)-\beta\,\delta Q(i),$ (72)

where we have exploited the fact that the initial ancilla state
$\hat{\eta}_{E_{i}}^{in}$ is thermal. The second principle of thermodynamics
can thus be expressed as:

$\displaystyle\Sigma(i)\geq 0\quad\forall i.$ (73)

###  IV.2 Thermodynamics of the system in Example 2

Figure 7: Plot of $W_{sw}$ and $\delta Q$ at the steady state. In both plots
$\omega_{S}=1$, $\omega_{1}=0.5$, $\omega_{2}=1.5$, $\alpha=0.2$. In panel (a)
the quantities are plotted as a function of $\kappa_{12}$ for different values
of the environmental temperature. It can be noticed that the work cost for
having coherence in the steady state is always equal to the heat dissipated
into the environment. Moreover, the amount of work is proportional to the
amount of coherence created in system $S$. In panel (b) the same quantities
are plotted as a function of the environmental temperature $\beta$ for
different values of $\kappa_{12}$. Also in this case it can be noticed that
the work cost is higher for lower temperatures, that is, for larger amounts of
coherence.

We now want to apply the formulas derived in the previous section to the
system in Sec. III.2. By substituting the Hamiltonians in Eqs.( 28, 39, 52) we
obtain explicit analytical expressions for the corresponding quantities, which
are reported in App. C.

The first term in the switching work is simply due to the mismatch of the
energy difference $E_{3}-E_{2}$ in the structured ancilla and the system’s
energy separation $\omega_{S}$ and is common to all collision models. The last
two terms arise because of the specific form of the interaction Hamiltonian of
the system with the structured ancillae.

In Fig. 6 we plot $\Delta U$, $\delta Q$ and $W_{sw}$ for the initial part of
the dynamics. One can observe, consistently with Eqs. ( C, C), that since
$\Delta U$ and $W_{sw}$ depend on the coherences in system $S$, both
quantities presents oscillations before reaching their steady state values.
Conversely, as the dissipated heat only depends on
$\langle\hat{\sigma}_{z}\rangle$ and on the ancillae coherences, no
oscillations are observed in this case.

After reaching the steady state, the variation of internal energy $\Delta U$
becomes identically zero. In the plots in Fig. 7 one can observe that in all
cases, at the steady state, a positive non-null amount of work is injected
into the system in order to keep the steady state coherence. This energy is
then exactly dissipated as heat into the thermal bath. This is perfectly
consistent with the first principle of thermodynamics, which after the
reaching of the steady state simply becomes $W_{sw}=-\delta Q$. This is also
consistent with the second law of thermodynamics: no work can be created with
a device operating with only one thermal bath. Moreover, in both plots one can
observe that the amount of work is proportional to the amount of coherence
created in the steady state, i.e. the larger the coherence the larger the work
needed.

Figure 8: Plot of the entropy production $\Sigma(i)$ at each step of the
dynamics. Panel (a): plot of $\Sigma(i)$ for $\kappa_{12}=0.3$ for different
values of $\beta$. One can notice that at the steady state (for large $n$) the
entropy variation is larger for larger values of the coherence in the steady
state. Panel (b): entropy production plotted for $\beta=1$ and for different
values of $\kappa_{12}$. Also in this case one can observe that the entropy
increase is larger for larger values of the coherence in the steady state. In
both cases one can also observe that the entropy variation is always positive,
thus fulfilling the second principle of thermodynamics at every step of the
dynamics.

Furthermore, in Fig. 8 we have plotted the step-by-step entropy production
$\Sigma(i)$. It is immediate to note that this quantity is always positive, in
fulfillment of the second principle of thermodynamics, and that it does never
get to zero, since even at the steady state work has to be injected into the
system and then dissipated as heat in order to retain coherence.

##  V Conclusions

In this work we have introduced a collision model in which the ancillae are
composite systems, rather than simple TLS or harmonic oscillators. We have
seen how to compute the discrete master equation for an open quantum system
interacting with a structured environment using our model. We have shown how
the presence of structured ancillae has consequences on the dynamics of $S$,
most noticeably the possibility of creating coherence in the system in the
steady state even if the environment is in a thermal state. We have seen
practically through a couple of examples involving qubits how the presence of
the structured environment can lead to additional terms in the master
equation, in particular the unitary driving terms leading to the steady state
coherence, and how this is possible only when the system interacts with the
degrees of freedom of the environment that possess coherence.

Furthermore, we have shown how the interaction between the system and the
degrees of freedom of the environment possessing coherence requires a work
cost to be activated, thus leading to the natural fulfillment of the first and
second principle of thermodynamics without invoking any further assumption. We
have also been able to compute analytically all the relevant thermodynamic
quantities, and shown that one needs to keep injecting work into the system in
order to retain the coherence in the steady state.

In future works we plan to use the model presented here to study the
efficiency of heat engines and refrigerators exploiting coherence, and how
including the cost of creating coherence in the system influences quantities
such as the efficiency and the efficiency at maximum power. Finally, the
present model might be useful also in different context than thermodynamics,
for instance for investigating whether the presence of coherence in the
environment can influence non Markovian dynamics or whether more complicated
structured ancillae can be exploited to further engineer the steady state of
the open quantum systems in order to obtain useful resources for quantum
technologies.

## Acknowledgements

SC thanks D. Tamascelli for a useful comment during IQIS 2022 conference. SC
acknowledges support by the Foundation for Polish Science (IRAP project,
ICTQT, Contract No. 2018/MAB/5, cofinanced by the EU within the Smart Growth
Operational Programme). G.D.C. acknowledges the support by the UK EPSRC
EP/S02994X/1, the Royal Society IEC\R2\222003.

## References

  * (1) K. Blum, Density Matrix Theory and Applications, Springer Berlin, Heidelberg (2012).
  * (2) H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, Oxford (2007).
  * (3) R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, Springer Berlin, Heidelberg (2007).
  * (4) C. Gardiner, Quantum Noise, Springer Berlin, Heidelberg (2004).
  * (5) H. Ma, U. Schollwöck, Z. Shuai, Density Matrix Renormalization Group (DMRG), Elsevier, (2021).
  * (6) G. De Chiara, M. Rizzi, D. Rossini and S. Montangero, Density Matrix Renormalization Group for Dummies, Journal of Computational and Theoretical Nanoscience, Vol. 5, 7, 1277-1288 (2008).
  * (7) J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * (8) R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117-158 (2014).
  * (9) D. Tamascelli, A. Smirne, S. F. Huelga and M. B. Plenio, Nonperturbative Treatment of non-Markovian Dynamics of Open Quantum Systems, Phys. Rev. Lett. 120, 030402 (2018).
  * (10) D. Tamascelli, A. Smirne, J. Lim,S. F. Huelga and M. B. Plenio, Efficient Simulation of Finite-Temperature Open Quantum Systems, Phys. Rev. Lett. 123, 090402 (2019).
  * (11) F. Ciccarello, S. Lorenzo, V. Giovannetti, G. M. Palma, Quantum collision models: Open system dynamics from repeated interactions, Phys. Rep. 954, 1 (2022).
  * (12) S. Cusumano, Quantum Collision Models: A Beginner Guide, Entropy 24 (9), 1258 (2022).
  * (13) S. Campbell and B. Vacchini, Collision models in open system dynamics: A versatile tool for deeper insights?, EPL 133, 60001 (2021).
  * (14) M. Cattaneo, G. De Chiara, S. Maniscalco, R. Zambrini and G. L. Giorgi, Collision Models Can Efficiently Simulate Any Multipartite Markovian Quantum Dynamics, Phys. Rev. Lett. 126, 130403 (2021).
  * (15) M. Ziman, P. Štelmachovič, V. Bužek, M. Hillery, V. Scarani and N. Gisin, Diluting quantum information: An analysis of information transfer in system-reservoir interactions, Phys. Rev. A 65, 042105 (2002).
  * (16) V. Scarani, M. Ziman, P. Štelmachovič, N. Gisin and V. Bužek, Thermalizing Quantum Machines: Dissipation and Entanglement, Phys. Rev. Lett. 88, 097905 (2002).
  * (17) M. Ziman and V. Bužek. Open system dynamics of simple collision models, Quantum Dynamics And Information, 199–227, (2010).
  * (18) F. Ciccarello, G. M. Palma and V. Giovannetti, Collision-model-based approach to non-Markovian quantum dynamics, Phys. Rev. A 87 4 (2013).
  * (19) N. Altamirano, P. Corona-Ugalde, R. B. Mann and M. Zych, Unitarity, feedback, interactions—dynamics emergent from repeated measurements, New J. Phys. 19, 013035 (2017).
  * (20) S. Lorenzo, F. Ciccarello and G. M. Palma, Class of exact memory-kernel master equations, Phys. Rev. A 93, 052111 (2016).
  * (21) F. Ciccarello and V. Giovannetti, A quantum non-Markovian collision model: incoherent swap case, Phys. Scr. 2013, 014010 (2013).
  * (22) S. Cusumano, A. Mari and V. Giovannetti, Interferometric quantum cascade systems, Phys. Rev. A 95, 053838 (2017).
  * (23) S. Cusumano, A. Mari and V. Giovannetti, Interferometric modulation of quantum cascade interactions, Phys. Rev. A 97, 053811 (2018).
  * (24) F. Ciccarello, Collision models in quantum optics, Quantum Meas. Quantum Metrol. 4.1 (2017).
  * (25) G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia and M. Antezza, Reconciliation of quantum local master equations with thermodynamics, New J. Phys. 20 113024 (2018).
  * (26) M. Cattaneo, G. L. Giorgi, R. Zambrini and S. Maniscalco, A Brief Journey through Collision Models for Multipartite Open Quantum Dynamics, Open Sys. and Inf. Dyn. 29, 03, 2250015 (2022).
  * (27) P. Strasberg, G. Schaller, T. Brandes, M. Esposito, Quantum and Information Thermodynamics: A Unifying Framework based on Repeated Interactions, Phys. Rev. X 7, 021003 (2017).
  * (28) S. Cusumano, V. Cavina, M. Keck, A. De Pasquale and V. Giovannetti, Entropy production and asymptotic factorization via thermalization: A collisional model approach, Phys. Rev. A 98, 032119 (2018).
  * (29) O. A. D. Molitor and G. T. Landi, Stroboscopic two-stroke quantum heat engines, Phys. Rev. A 102, 042217 (2020).
  * (30) F. L. S. Rodrigues, G. De Chiara, M. Paternostro and G. T. Landi, Thermodynamics of Weakly Coherent Collisional Models, Phys. Rev. Lett. 123, 140601 (2019).
  * (31) F. Barra, The thermodynamic cost of driving quantum systems by their boundaries, Sci. Rep. 5, 14873 (2015).
  * (32) N. Piccione, G. De Chiara and B. Bellomo, Power maximization of two-stroke quantum thermal machines, Phys. Rev. A 103, 032211 (2021).
  * (33) Kenza Hammam, Heather Leitch, Yassine Hassouni, Gabriele De Chiara, Exploiting coherence for quantum thermodynamic advantage, New J. Phys. 24 113053 (2022).
  * (34) G. De Chiara and M. Antezza, Quantum machines powered by correlated baths, Phys. Rev. Research 2, 033315 (2020).
  * (35) H-P. Breuer, Non-Markovian generalization of the Lindblad theory of open quantum systems, Phys. Rev. A 75, 022103 (2007).
  * (36) D. Suess, A. Eisfeld, and W. T. Strunz, Hierarchy of Stochastic Pure States for Open Quantum System Dynamics, Phys. Rev. Lett. 113, 150403 (2014).
  * (37) C-K. Chan, G-D. Lin, S. F. Yelin, and M. D. Lukin, Quantum interference between independent reservoirs in open quantum systems, Phys. Rev. A 89, 042117 (2014).
  * (38) G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett and J. Keeling, Efficient Exploration of Hamiltonian Parameter Space for Optimal Control of Non-Markovian Open Quantum Systems, Phys. Rev. Lett. 126, 200401 (2021).
  * (39) H-P. Breuer, E-M. Laine, J. Piilo and B. Vacchini, Colloquium: Non-Markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * (40) A. Manatuly, W. Niedenzu, R. Román-Ancheyta, B. Çakmak, Ö. E. Müstecaplıoğlu and G. Kurizki, Collectively enhanced thermalization via multiqubit collisions, Phys. Rev. E 99, 042145 (2019).
  * (41) C. B. Dağ, W. Niedenzu, Ö. E. Müstecaplıoğlu and G. Kurizki, Multiatom Quantum Coherences in Micromasers as Fuel for Thermal and Nonthermal Machines, Entropy, 18(7), 244 (2016).
  * (42) R. Román-Ancheyta, M. Kolář, G. Guarnieri and R. Filip, Enhanced steady-state coherence via repeated system-bath interactions, Phys. Rev. A 104, 062209 (2021).
  * (43) Y. Li, X. Li and J. Jin, Quantum nonstationary phenomena of spin systems in collision models, Phys. Rev. A 107, 042205 (2023).
  * (44) A. Purkayastha, G. Guarnieri, S. Campbell, J. Prior and J. Goold, Periodically refreshed baths to simulate open quantum many-body dynamics, Phys. Rev. B 104, 045417 (2021).
  * (45) A. Purkayastha, G. Guarnieri, S. Campbell, J. Prior and J. Goold, Periodically refreshed quantum thermal machines, Quantum 6, 801 (2022).
  * (46) G. Guarnieri, M. Kolář, and R. Filip, Steady-State Coherences by Composite System-Bath Interactions, Phys. Rev. Lett. 121, 070401 (2018); Erratum Phys. Rev. Lett. 127, 129901 (2021).
  * (47) M. Lostaglio and M. P. Müller, Coherence and Asymmetry Cannot be Broadcast, Phys. Rev. Lett. 123, 020403 (2019).
  * (48) I. Marvian and R. W. Spekkens, No-Broadcasting Theorem for Quantum Asymmetry and Coherence and a Trade-off Relation for Approximate Broadcasting, Phys. Rev. Lett. 123, 020404 (2019).
  * (49) Philippe Faist, Jonathan Oppenheim and Renato Renner, Gibbs-preserving maps outperform thermal operations in the quantum regime, New J. Phys. 17, 043003 (2015).
  * (50) T. Baumgratz,M. Cramer and M. B. Plenio, Quantifying Coherence, Phys. Rev. Lett. 113, 140401 (2014).
  * (51) K. Hammam, Y. Hassouni, R. Fazio and G. Manzano, Optimizing autonomous thermal machines powered by energetic coherence, New J. Phys. 23 043024 (2021).

##  Appendix A Calculations of Example 1

In order to compute the constants $g_{k}^{(j)}$ and
$g_{kk^{\prime}}^{(j,j^{\prime})}$ appearing in Eq. (48) we define the non-
interacting energy base of the 2 qubits as:

$\displaystyle\ket{E_{1}}=\ket{\downarrow\downarrow},\;\ket{E_{2}}=\ket{\downarrow\uparrow},$
$\displaystyle\ket{E_{3}}=\ket{\uparrow\downarrow},\;\ket{E_{4}}=\ket{\uparrow\uparrow}.$
(74)

The energy eigenbasis of the Hamiltonian $\hat{H}_{E_{i}}$ can be easily
computed through standard methods, and is defined as:

$\displaystyle\ket{E^{\prime}_{1}}=\ket{E_{1}},$ (75)
$\displaystyle\ket{E^{\prime}_{2}}=\left(a_{-}\ket{E_{2}}+b_{-}\ket{E_{3}}\right),$
(76)
$\displaystyle\ket{E^{\prime}_{3}}=\left(a_{+}\ket{E_{2}}+b_{+}\ket{E_{3}}\right),$
(77) $\displaystyle\ket{E^{\prime}_{4}}=\ket{E_{4}},$ (78)

where the constants $a_{\pm},b_{\pm}$ are worth:

$\displaystyle
a_{-}=\frac{\kappa_{12}}{\sqrt{\kappa_{12}^{2}+(E^{\prime}_{2}-y)^{2}}},\quad
b_{-}=\frac{E^{\prime}_{2}-y}{\sqrt{\kappa_{12}^{2}+(E^{\prime}_{2}-y)^{2}}},$
(79) $\displaystyle
a_{+}=\frac{\kappa_{12}}{\sqrt{\kappa_{12}^{2}+(E^{\prime}_{3}-y)^{2}}},\quad
b_{+}=\frac{E^{\prime}_{3}-y}{\sqrt{\kappa_{12}^{2}+(E^{\prime}_{3}-y)^{2}}}.$

The two energy basis are related via the unitary $\hat{V}$ as:

$\displaystyle\ket{E^{\prime}_{k}}=\sum_{i}V_{ik}\ket{E_{i}},$ (80)

where the unitary $\hat{V}$ is worth:

$\displaystyle\hat{V}=\begin{bmatrix}1&0&0&0\\\ 0&b_{-}&b_{+}&0\\\
0&a_{-}&a_{+}&0\\\ 0&0&0&1\end{bmatrix}.$ (81)

At this point we are ready to compute the constants
$g_{k}^{(j,j^{\prime})},g_{kk^{\prime}}^{(j,j^{\prime})}$. The expression for
the thermal state is:

$\displaystyle\hat{\eta}_{E_{i}}=\frac{e^{-\beta\hat{H}_{E_{i}}}}{Z}$
$\displaystyle=$ $\displaystyle\frac{1}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}\outerproduct{E^{\prime}_{k}}{E^{\prime}_{k}}=\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V_{kj}^{\dagger}\outerproduct{E_{i}}{E_{j}}.$ (82)

To illustrate how the calculation works, we show how to compute $g_{1}^{(+)}$:

$\displaystyle g_{1}^{(+)}$ $\displaystyle=$
$\displaystyle\alpha_{1}\Tr\left[\hat{\sigma}_{1}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{1}}{Z}\Tr\left[\hat{\sigma}_{1}^{+}\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V^{\dagger}_{kj}\outerproduct{E_{i}}{E_{j}}\right]$ (83)
$\displaystyle=$
$\displaystyle\frac{\alpha_{1}}{Z}\Tr\left[\sum_{j,k}e^{-\beta
E^{\prime}_{k}}\left(V_{1k}V^{\dagger}_{kj}\outerproduct{E_{3}}{E_{j}}+V_{2k}V^{\dagger}_{kj}\outerproduct{E_{4}}{E_{j}}\right)\right]$
$\displaystyle=$ $\displaystyle\frac{\alpha_{1}}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}\left(V_{1k}V^{\dagger}_{k3}+V_{2k}V^{\dagger}_{k4}\right)=0.$

and similarly for $g_{1}^{(-)},g_{2}^{(\pm)}$.

As for $g_{11}^{(+,-)}$ we have:

$\displaystyle
g_{11}^{(+,-)}=\alpha_{1}^{2}\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{1}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{1}^{2}}{Z}\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{1}^{+}\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V_{kj}^{\dagger}\outerproduct{E_{i}}{E_{j}}\right]$
$\displaystyle=\frac{\alpha_{1}^{2}}{Z}\Tr\left[(\outerproduct{E_{1}}{E_{1}}+\outerproduct{E_{2}}{E_{2}})\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V_{kj}^{\dagger}\outerproduct{E_{i}}{E_{j}}\right]$
$\displaystyle=\frac{\alpha_{1}^{2}}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}\left(V_{1k}V^{\dagger}_{k1}+V_{2k}V^{\dagger}_{k2}\right)$
$\displaystyle=\frac{\alpha_{1}^{2}}{Z}\left(e^{-\beta
E^{\prime}_{1}}+|b_{-}|^{2}e^{-\beta E^{\prime}_{2}}+|b_{+}|^{2}e^{-\beta
E^{\prime}_{3}}\right).$ (84)

and one can compute in a similar manner:

$\displaystyle g_{11}^{(-,+)}$ $\displaystyle=$
$\displaystyle\alpha_{1}^{2}\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{1}^{-}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{1}^{2}}{Z}\left(|a_{-}|^{2}e^{-\beta
E^{\prime}_{2}}+|a_{+}|^{2}e^{-\beta E^{\prime}_{3}}+e^{-\beta
E^{\prime}_{4}}\right),$ (85) $\displaystyle g_{22}^{(+,-)}$ $\displaystyle=$
$\displaystyle\alpha_{2}^{2}\Tr\left[\hat{\sigma}_{2}^{-}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{2}^{2}}{Z}\left(e^{-\beta
E^{\prime}_{1}}+|a_{-}|^{2}e^{-\beta E^{\prime}_{2}}+|a_{+}|^{2}e^{-\beta
E^{\prime}_{3}}\right),$ (86) $\displaystyle g_{22}^{(-,+)}$ $\displaystyle=$
$\displaystyle\alpha_{2}^{2}\Tr\left[\hat{\sigma}_{2}^{+}\hat{\sigma}_{2}^{-}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{2}^{2}}{Z}\left(|b_{-}|^{2}e^{-\beta
E^{\prime}_{2}}+|b_{+}|^{2}e^{-\beta E^{\prime}_{3}}+e^{-\beta
E^{\prime}_{4}}\right).$ (87)

Finally, we compute the cross term $g_{12}^{(+,-)}$ as:

$\displaystyle g_{12}^{(+,-)}$ $\displaystyle=$ $\displaystyle
g_{21}^{(+,-)*}=\alpha_{1}\alpha_{2}^{*}\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{1}\alpha_{2}}{Z}\Tr\left[\outerproduct{E_{2}}{E_{3}}\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V_{kj}^{\dagger}\outerproduct{E_{i}}{E_{j}}\right]$ (88)
$\displaystyle=$ $\displaystyle\frac{\alpha_{1}\alpha_{2}}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}V_{3k}V_{k2}^{\dagger}=\frac{\alpha_{1}\alpha_{2}}{Z}\left(e^{-\beta
E^{\prime}_{2}}a_{-}b_{-}^{*}+e^{-\beta E^{\prime}_{3}}a_{+}b_{+}^{*}\right).$

Similarly we can compute the other cross terms:

$\displaystyle g_{12}^{(-,+)}$ $\displaystyle=$ $\displaystyle
g_{21}^{(-,+)*}=\alpha_{1}^{*}\alpha_{2}\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha_{1}\alpha_{2}}{Z}\left(e^{-\beta
E^{\prime}_{2}}b_{-}a_{-}^{*}+e^{-\beta E^{\prime}_{3}}b_{+}a_{+}^{*}\right),$
(89) $\displaystyle g_{12}^{(+,+)}$ $\displaystyle=$ $\displaystyle
g_{21}^{(+,+)}=\alpha_{1}\alpha_{2}\Tr\left[\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-}\hat{\eta}_{E_{i}}\right]=0,$
(90) $\displaystyle g_{12}^{(-,-)}$ $\displaystyle=$ $\displaystyle
g_{21}^{(-,-)}=\alpha_{1}\alpha_{2}\Tr\left[\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{+}\hat{\eta}_{E_{i}}\right]=0.$
(91)

##  Appendix B Calculations of Example 2

First we want to compute the factors $g_{E_{i}}^{(\pm)}$:

$\displaystyle g_{E_{i}}^{(+)}$ $\displaystyle=$
$\displaystyle\alpha\Tr\left[\hat{\sigma}_{E_{i}}^{+}\hat{\eta}_{E_{i}}\right]=\frac{\alpha}{Z}\Tr\left[\outerproduct{E_{3}}{E_{2}}\sum_{k}e^{-\beta
E^{\prime}_{k}}\outerproduct{E^{\prime}_{k}}{E^{\prime}_{k}}\right]=\frac{\alpha}{Z}\Tr\left[\outerproduct{E_{3}}{E_{2}}\sum_{i,j,k}V_{ik}V^{\dagger}_{kj}e^{-\beta
E_{i}}\outerproduct{E_{j}}{E_{j}}\right]$ (92) $\displaystyle=$
$\displaystyle\frac{\alpha}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}V_{2k}V^{\dagger}_{k3}=\frac{\alpha}{Z}\left[e^{-\beta
E^{\prime}_{2}}a_{-}b_{-}^{*}+e^{-\beta E^{\prime}_{3}}a_{+}b_{+}^{*}\right],$

and similarly we get:

$\displaystyle g_{E_{i}}^{(-)}=\frac{\alpha^{*}}{Z}\left[e^{-\beta
E^{\prime}_{2}}a_{-}^{*}b_{-}+e^{-\beta
E^{\prime}_{3}}a_{+}^{*}b_{+}\right]=g_{E_{i}}^{(+)*}.$ (93)

Then we move to the calculation of the factors $g_{E_{i}}^{(\pm,\pm)}$,
starting from $g_{E_{i}}^{(-,+)}$:

$\displaystyle
g_{E_{i}}^{(-,+)}=\Tr\left[\hat{\sigma}_{E_{i}}^{+}\hat{\sigma}_{E_{i}}^{-}\hat{\eta}_{E_{i}}\right]=\frac{|\alpha|^{2}}{Z}\Tr\left[\outerproduct{E_{3}}{E_{2}}\outerproduct{E_{2}}{E_{3}}\sum_{k}e^{-\beta
E^{\prime}_{k}}\outerproduct{E^{\prime}_{k}}{E^{\prime}_{k}}\right]$ (94)
$\displaystyle=$
$\displaystyle\frac{|\alpha|^{2}}{Z}\Tr\left[\outerproduct{E_{3}}{E_{3}}\sum_{i,j,k}e^{-\beta
E^{\prime}_{k}}V_{ik}V^{\dagger}_{kj}\outerproduct{E_{i}}{E_{j}}\right]=\frac{|\alpha|^{2}}{Z}\sum_{k}e^{-\beta
E^{\prime}_{k}}V_{3k}V^{\dagger}_{k3}=\frac{|\alpha|^{2}}{Z}\left[e^{-\beta
E^{\prime}_{2}}|a_{-}|^{2}+e^{-\beta E^{\prime}_{3}}|a_{+}|^{2}\right],$

and similarly for $g_{E_{i}}^{(+,-)}$ we find:

$\displaystyle g_{E_{i}}^{(+,-)}=\frac{|\alpha|^{2}}{Z}\left[e^{-\beta
E^{\prime}_{2}}|b_{-}|^{2}+e^{-\beta E^{\prime}_{3}}|b_{+}|^{2}\right].$ (95)

Finally, it is immediate to see that:

$\displaystyle
g_{E_{i}}^{(+,+)}=\Tr\left[\hat{\sigma}_{E_{i}}^{+}\hat{\sigma}_{E_{i}}^{+}\hat{\eta}_{E_{i}}\right]=0,$
(96) $\displaystyle
g_{E_{i}}^{(-,-)}=\Tr\left[\hat{\sigma}_{E_{i}}^{-}\hat{\sigma}_{E_{i}}^{-}\hat{\eta}_{E_{i}}\right]=0.$
(97)

The matrix elements are mapped by the master equation in (56) according to:

$\displaystyle\rho_{11}\rightarrow\rho_{11}+\delta
t^{2}\left(g_{E_{i}}^{(-,+)}\rho_{22}-g_{E_{i}}^{(+,-)}\rho_{11}\right)+\frac{\delta
t^{2}}{2}\omega_{S}\left(g_{E_{i}}^{(+)}\rho_{12}+g_{E_{i}}^{(-)}\rho_{21}\right)+i\delta
t\left(g_{E_{i}}^{(+)}\rho_{12}-g_{E_{i}}^{(-)}\rho_{21}\right),$ (98)
$\displaystyle\rho_{22}\rightarrow\rho_{22}-\delta
t^{2}\left(g_{E_{i}}^{(-,+)}\rho_{22}-g_{E_{i}}^{(+,-)}\rho_{11}\right)-\frac{\delta
t^{2}}{2}\omega_{S}\left(g_{E_{i}}^{(+)}\rho_{12}+g_{E_{i}}^{(-)}\rho_{21}\right)-i\delta
t\left(g_{E_{i}}^{(+)}\rho_{12}-g_{E_{i}}^{(-)}\rho_{21}\right),$ (99)
$\displaystyle\rho_{12}\rightarrow\rho_{12}+\frac{\delta
t^{2}}{2}\omega_{S}g_{E_{i}}^{(-)}(\rho_{11}-\rho_{22})+i\delta
tg_{E_{i}}^{(-)}(\rho_{11}-\rho_{22})-i\delta
t\omega_{S}\rho_{12}-\frac{\delta
t^{2}}{2}\rho_{12}\left(g_{E_{i}}^{(+,-)}+g_{E_{i}}^{(-,+)}+\omega_{S}^{2}\right),\qquad$
(100) $\displaystyle\rho_{21}\rightarrow\rho_{21}+\frac{\delta
t^{2}}{2}\omega_{S}g_{E_{i}}^{(+)}(\rho_{11}-\rho_{22})-i\delta
tg_{E_{i}}^{(+)}(\rho_{11}-\rho_{22})+i\delta
t\omega_{S}\rho_{21}-\frac{\delta
t^{2}}{2}\rho_{21}\left(g_{E_{i}}^{(+,-)}+g_{E_{i}}^{(-,+)}+\omega_{S}^{2}\right).$
(101)

The steady state values of $\langle\hat{\sigma}_{x}\rangle$,
$\langle\hat{\sigma}_{y}\rangle$, $\langle\hat{\sigma}_{z}\rangle$ are:

$\displaystyle\langle\hat{\sigma}_{z}\rangle=\frac{\left(g_{E}^{(-,+)}-g_{E}^{(+,-)}\right)\left[4\omega_{S}^{2}+\delta
t^{2}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\omega_{S}^{2}\right)^{2}\right]}{K},$
(102)
$\displaystyle\langle\hat{\sigma}_{x}\rangle=\frac{\left(g_{E}^{(-,+)}-g_{E}^{(+,-)}\right)2\real\left[g_{E}^{(-)}\left(4\omega_{S}+2i\delta
t\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\right)+\delta
t^{2}\omega_{S}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}+\omega_{S}^{2}\right)\right)\right]}{K},$
(103)
$\displaystyle\langle\hat{\sigma}_{y}\rangle=-\frac{\left(g_{E}^{(-,+)}-g_{E}^{(+,-)}\right)2\imaginary\left[g_{E}^{(-)}\left(4\omega_{S}+2i\delta
t\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\right)+\delta
t^{2}\omega_{S}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}+\omega_{S}^{2}\right)\right)\right]}{K},\qquad\quad$
(104) $\displaystyle\quad
K=8g_{E}^{(-)}g_{E}^{(+)}\left(g_{E}^{(-,+)}+g_{E}^{(-,+)}-\omega_{S}^{2}\right)+4\omega_{S}^{2}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\right)$
$\displaystyle\quad+\delta
t^{2}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\right)\left[\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}\right)^{2}+2\omega_{S}^{2}\left(g_{E}^{(-,+)}+g_{E}^{(+,-)}-g_{E}^{(-)}g_{E}^{(+)}\right)+\omega_{S}^{4}\left(1-2g_{E}^{(-)}g_{E}^{(+)}\right)\right].$
(105)

##  Appendix C Thermodynamic expressions

Explicit expressions for Eqs.( 61, IV.1):

$\displaystyle\Delta U=-i\delta
t\sum_{j,k}g_{k}^{(j)}\Tr\left[\commutator{\hat{H}_{S}}{\hat{A}_{S}^{(k,j)}}\hat{\rho}_{S}(i)\right]\frac{\delta
t^{2}}{2}\sum_{k,j}g_{k}^{(j)}\Tr\left[\commutator{\hat{A}_{S}^{(k,j){\dagger}}}{\hat{H}_{S}}\hat{H}_{S}\hat{\rho}_{S}(i)\right]$
$\displaystyle+\frac{\delta
t^{2}}{2}\sum_{k,j}g_{k}^{(j)}\Tr\left[\commutator{\hat{H}_{S}}{\hat{A}_{S}^{(k,j)}}\hat{\rho}_{S}(i)\hat{H}_{S}\right]+\frac{g\delta
t^{2}}{2}\sum_{k,k^{\prime},j,j^{\prime}}g_{kk^{\prime}}^{(j,j^{\prime})}\Tr\left[\left(2\hat{A}_{S}^{(k,j)}\hat{H}_{S}\hat{A}_{S}^{(k,j){\dagger}}-\anticommutator{\hat{A}_{S}^{(k,j){\dagger}}\hat{A}_{S}^{(k,j)}}{\hat{H}_{S}}\right)\hat{\rho}_{S}(i)\right]$
$\displaystyle+\frac{\delta
t^{2}}{2}\sum_{j,k}\Tr\left[\commutator{\hat{H}_{S}}{g_{kH_{E}}^{(j)}\hat{A}_{S}^{(k,j)}-g_{kH_{E}}^{(j)*}\hat{A}_{S}^{(k,j){\dagger}}}\hat{\rho}_{S}(i)\right],$
(106)

for the variation of internal energy, while the heat exchanged reads:

$\displaystyle\delta Q=i\delta
t\sum_{k,j}(g_{kH_{E}}^{(j)}-g_{kH_{E}}^{(j)*})\Tr\left[\hat{A}_{S}^{(k,j)}\hat{\rho}_{S}(i)\right]-\frac{\delta
t^{2}}{2}\sum_{j,k}\Tr\left[\commutator{\hat{H}_{S}}{g_{kH_{E}}^{(j)}\hat{A}_{S}^{(k,j)}-g_{kH_{E}}^{(j)*}\hat{A}_{S}^{(k,j){\dagger}}}\hat{\rho}_{S}(i)\right]$
$\displaystyle+\frac{\delta
t^{2}}{2}\sum_{j,k}\Tr\left[\hat{A}_{S}^{(j,k)}\hat{\rho}_{S}(i)\right]\Tr\left[(2\hat{H}_{E}\hat{B}_{Ek}^{(j)}\hat{H}_{E}-\anticommutator{\hat{H}_{E}^{2}}{\hat{B}_{Ek}^{(j)}})\hat{\eta}_{E_{i+1}}\right]$
$\displaystyle-\frac{\delta
t^{2}}{2}\sum_{k,k^{\prime},j,j^{\prime}}\Tr\left[\hat{A}_{S}^{(k,j){\dagger}}\hat{A}_{S}^{(k,j)}\hat{\rho}_{S}\right]\Tr\left[\commutator{\hat{H}_{E}}{\hat{B}_{Ek}^{(j)}}\hat{\eta}_{E_{i+1}}\hat{B}_{Ek}^{(j){\dagger}}\right]+\frac{\delta
t^{2}}{2}\sum_{k,k^{\prime},j,j^{\prime}}\Tr\left[\hat{A}_{S}^{(k,j)}\hat{A}_{S}^{(k,j){\dagger}}\hat{\rho}_{S}\right]\Tr\left[\commutator{\hat{H}_{E}}{\hat{B}_{Ek}^{(j)}}\hat{B}_{Ek}^{(j){\dagger}}\hat{\eta}_{E_{i+1}}\right].$

Finally, we will not write explicitly the switching work in Eq. (66), since
one can easily derive it from:

$\displaystyle W_{sw}=\Delta U-\delta Q.$ (108)

The explicit expressions of $\Delta U$, $\delta Q$ and $W_{sw}$ for the system
in Example 2 read:

$\displaystyle\Delta U=$ $\displaystyle-$ $\displaystyle\frac{\delta
t^{2}}{2}\omega_{S}\left[\expectationvalue{\hat{\sigma}_{S}^{z}}\left(g_{E_{i}}^{(+,-)}+g_{E_{i}}^{(-,+)}\right)+\left(g_{E_{i}}^{(+,-)}-g_{E_{i}}^{(-,+)}\right)\right]$
$\displaystyle+$ $\displaystyle\frac{\delta
t^{2}}{2}\omega_{S}^{2}\left(g_{E_{i}}^{(+)}\rho_{12}+g_{E_{i}}^{(-)}\rho_{21}\right)+i\delta
t\omega_{S}\left(g_{E_{i}}^{(+)}\rho_{12}-g_{E_{i}}^{(-)}\rho_{21}\right),$
$\displaystyle\delta Q$ $\displaystyle=$ $\displaystyle-\frac{\delta
t^{2}}{2}(E_{3}-E_{2})\left[\expectationvalue{\hat{\sigma}_{S}^{z}}\left(g_{E_{i}}^{(+,-)}+g_{E_{i}}^{(-,+)}\right)+\left(g_{E_{i}}^{(+,-)}-g_{E_{i}}^{(-,+)}\right)\right]+\frac{\delta
t^{2}}{2}\kappa_{12}|\alpha|^{2}\left(\eta_{23}+\eta_{32}\right),$ (110)
$\displaystyle W_{sw}=$ $\displaystyle\frac{\delta
t^{2}}{2}(E_{3}-E_{2}-\omega_{S})\left[\expectationvalue{\hat{\sigma}_{S}^{z}}\left(g_{E_{i}}^{(+,-)}+g_{E_{i}}^{(-,+)}\right)+\left(g_{E_{i}}^{(+,-)}-g_{E_{i}}^{(-,+)}\right)\right]$
$\displaystyle-$ $\displaystyle\frac{\delta
t^{2}}{2}\kappa_{12}|\alpha|^{2}\left(\eta_{23}+\eta_{32}\right)+\frac{\delta
t^{2}}{2}\omega_{S}^{2}\left(g_{E_{i}}^{(+)}\rho_{12}+g_{E_{i}}^{(-)}\rho_{21}\right)+i\delta
t\omega_{S}t\left(g_{E_{i}}^{(+)}\rho_{12}-g_{E_{i}}^{(-)}\rho_{21}\right).$


# Correlation between the solar magnetic field strength and the millimeter
brightness temperature

I. SKOKIĆ<EMAIL_ADDRESS>D. SUDAR and R. BRAJŠA
Hvar Observatory Faculty of Geodesy  University of Zagreb
Kačićeva 26 HR–10000 Zagreb  Croatia

###### Abstract

Images of the Sun at millimeter wavelengths obtained by ALMA show a
significant correspondence with the magnetograms. In this paper, we
investigate this correspondence by comparing ALMA full-disk solar image taken
at 1.2 mm with a SDO/HMI magnetogram and analyze their correlation. It is
found that chromospheric network and active regions show a positive
correlation where brightness temperature is increasing with the line-of-sight
magnetic field strength, while sunspots have a negative correlation. Quiet Sun
regions do not show any dependence of the brightness temperature with the
magnetic field. Thermal bremsstrahlung is given as the best explanation for
the observed correlations.

###### keywords:

Sun: magnetic fields - Sun: chromosphere - Sun: radio radiation

##  1 Introduction

Solar radiation in the (sub)millimeter wavelength range is mostly generated by
free–free (thermal bremsstrahlung) emission from the interaction of free
electrons with ions within the chromosphere, while in regions with strong
magnetic fields, such as sunspots, a gyroresonant emission coming from the
electrons spiraling along the magnetic field lines also plays an important
role (see, e.g., ???). Theoretical models predict that magnetic field
modulates free–free opacity differently for ordinary and extraordinary wave
modes, which should be measurable as somewhat lower and higher intensity from
the non-magnetic level, when observing in circularly polarized radiation (?).
For gyroresonant emission, however, models predict that very strong magnetic
fields are neccessary to have significant contribution to total radiation in
the (sub)millimeter range, much stronger than those observed (?).

Solar observations with the Atacama Large Millimeter/submillimeter Array
(ALMA) have provided us with a new tool to study the Sun at (sub)millimeter
wavelength range (?). It has already produced new and interesting results
(e.g. ???????). However, ALMA is not yet capable of providing circular
polarization measurements of the Sun.

Various features in solar ALMA images have been identified to have a good
correspondence with known solar structures. ? have reported that, at 1.2 and 3
mm, sunspots and inversion lines are darker, active regions and coronal bright
points are brighter, while coronal holes and prominences have very small
contrast from the quiet Sun background. They also found a very good
correspondence with the Solar Dynamics Observatory (SDO; ?) Atmospheric
Imaging Assembly (AIA; ?) 304 Å and 1600 Å channels, and Helioseismic and
Magnetic Imager (HMI; ?) magnetograms.

The main goal of this paper is to analyze and quantize this good
correspondence of ALMA solar images with HMI magnetograms indicating a close
relationship between the observed millimeter brightness temperature and the
magnetic field strength.

##  2 Data and method

We used a single dish image of the Sun from the ALMA Solar Commissioning and
Science Verification Campaign (CSV), taken on December 18, 2015 at 20:12 UT
(?). This image was selected because it was taken when the Sun was still very
active and it contains examples of many different solar regions such as active
regions, quiet Sun network and inter-network regions, sunspots, etc. Moreover,
this data is well known and has been analyzed by many researchers. The image
was made by fast scanning the Sun using a double circle pattern with a 12 m PM
antenna in the band 6, in the spectral window 3 with a rest frequency of 248
GHz or equivalent wavelength of 1.2 mm. Limb brightening was removed from the
ALMA image by the method presented in ? and the image was scaled to the quiet
Sun level of 5900 K, as suggested by ?.

For comparison with the ALMA image, a line-of-sight (LOS) magnetogram from the
HMI instrument aboard SDO was used. The HMI image was first convolved with a
Gaussian kernel with a standard deviation of
$\sigma=\textrm{FWHM}/(2\sqrt{2\ln 2})\approx\textrm{FWHM}/2.355$, where FWHM
is the Full Width at Half Maximum of the ALMA band 6 beam of the 12 m antenna
(equal to 26.7 arcsec, which is also the effective spatial resolution). After
that, the 4096x4096 pixel HMI image was resized to match 800x800 pixel ALMA
image with an image scale of 3 arcsec per pixel. The ALMA image, the resized
HMI image, and the convolved and resized HMI image are shown in Fig. 1, in the
left, middle and right panels, respectively. ALMA contours of 6100 K were
overlaid on both HMI images for easier comparison with the ALMA image. A good
correspondence can be seen between the bright ALMA regions and the HMI areas
with a large magnetic field intensity.

Figure 1: HMI magnetogram resized (center image) to ALMA image (left) size and
convolved with the ALMA beam (right). ALMA contours of 6100 K were overlaid on
both HMI images. For better contrast, HMI and convolved HMI images were
normalized to $\pm 100$ G and $\pm 20$ G, respectively. Figure 2: Profile of
the line cutting through the center of the Sun. HMI data resized to ALMA image
size is shown in gray, convolved and resized HMI data is shown in blue, while
ALMA intensity is shown in red. Note the different axes for HMI and ALMA data.
Black horizontal line denotes LOS magnetic field strength of zero Gauss and
quiet Sun brightness temperature for ALMA band 6 of 5900 K.

The main idea is to compare the ALMA brightness temperature with the HMI LOS
magnetic field strength and look for possible correlation between the two.
First, the single line profiles from ALMA and HMI, cutting through the center
of the Sun, are compared. Then, the whole images are compared on a pixel-by-
pixel basis.

##  3 Results and discussion

Profiles of the line cutting through the center of the Sun are shown in Fig.
2. The gray line denotes resized only HMI data, the blue line denotes
convolved and resized HMI data, while the red line represents ALMA data. A
horizontal black line of zero Gauss is added to mark the polarity reversal of
the HMI data, corresponding also to the quiet Sun level of 5900 K. It can be
seen that ALMA brightness temperature ($T_{b}$) is correlated pretty well with
the positive, and anti-correlated with the negative values of the convolved
HMI LOS field strength ($B_{LOS}$). This (anti)correlation is not that obvious
when using the non-convolved HMI data, implicating the importance of correctly
preprocessing the data and performing the comparison at the same spatial scale
and resolution.

Figure 3: ALMA 248 GHz brightness temperature compared with HMI LOS magnetic
field $B_{LOS}$ (left). Data set was divided into groups based on $B_{LOS}$
and marked with different colors, and the corresponding regions are shown in
the HMI image on the right. Tbe two wings, left (yellow) and right (green),
come from two sunspots of opposite polarity, as denoted in the right panel.
Figure 4: ALMA brightness temperature and HMI LOS magnetic field strength
comparison as in Fig. 3, but zoomed to a region of $\pm 60$ G (left) and with
a logarithmic scale to show the dependence in the region of low magnetic field
values (right). In the right panel, absolute LOS values are shown. The same
colors are used as in Fig. 3 Colored lines denote binned values of data
subgroups while white/black line is binned over the whole data set. Bins have
equal number of data in each bin.

The comparison of the whole disk $T_{b}$ and $B_{LOS}$ values is shown in the
left panel of Fig. 3. Only pixels inside the solar disk with a cutoff radius
of 300 pixels (900 arcsec) were used. The apparent photospheric radius of the
Sun ($6.96\times 10^{8}$ m) had a value of 975.3 arcsec or 325 pixels.

The result is an interesting but complicated shape that suggests that there
might be several $T_{b}$ – $B_{LOS}$ correlations hidden inside the graph. To
investigate this further, pixels were divided into five subgroups denoting
different regions of the graph, based on their HMI $B_{LOS}$ value:
$B_{LOS}<-120$ G (yellow), $-120<B_{LOS}<-1.2$ G (blue), $-1.2<B_{LOS}<1.2$ G
(black), $1.2<B_{LOS}<120$ G (red), and $B_{LOS}>120$ G (green). Not
surprisingly, these subgroups actually trace different regions on the Sun,
namely, positive (red) and negative (blue) polarity active regions and
chromospheric/magnetic network, positive (green) and negative (yellow)
polarity sunspots, and the remaining pixels come from the quiet Sun/inter-
network regions.

Interesting is a negative correlation between $T_{b}$ and absolute value of
$B_{LOS}$ within both the positive and negative polarity sunspot; higher
absolute values of $B_{LOS}$ correspond to lower $T_{b}$. The correlation is
not perfect, there is some structure most probably coming from systematic
errors like the antenna scanning pattern, which is still present in the ALMA
data as small regular variations of the brightness temperature. The ALMA Solar
Development Team is actively working to correct these systematic errors and
other artifacts still present in many solar ALMA single-dish data.

On the other hand, active regions suggest a positive correlation where higher
absolute $B_{LOS}$ values correspond to higher $T_{b}$. This is visible as a
narrow V shape in the left panel of Fig. 3. This V shape is more obvious in a
zoomed-in image in the left panel of Fig. 4, indicating a clear link between
the $B_{LOS}$ and $T_{b}$ in the active solar regions. However, a small
asymmetry between positive and negative active regions visible in Fig. 3 (red
and blue stripe), is not visible in the zoomed-in Fig. 4.

To check if there is a correlation in quiet Sun regions as well, a logarithmic
scale for $B_{LOS}$ is used in the right panel of Fig. 4, since average
magnetic field strength is low in these regions. They are denoted with black
pixels. It can be seen that these regions do not correlate with $B_{LOS}$,
rather they have more or less a constant value of $T_{b}\approx 5900$ K.

To explore further these interesting relationships, we binned the whole data
set into bins with equal number of data points. The result is shown as white
and black line in the right panel of Fig. 4. Each subgroup was also binned and
shown with the line of the same but darker color, with 1-$\sigma$ values
denoted by crosses at bin centers. Bins just confirm previous conclusions:
regions of $|B_{LOS}|<1.2$ G do not show any correlation, regions of
$|B_{LOS}|$ between 1.2 and 120 G show a positive correlation, while regions
with $|B_{LOS}|$ above $120$ G show a negative correlation. We calculated
correlation coefficients for the five subgroups of the data and got 0.735
(0.736) for positive (negative) polarity active regions, -0.916 (-0.619) for
the negative (positive) polarity sunspot, and 0.074 for the quiet Sun region.

There is a significant scatter in the presented graphs with possible hints of
a substructure (Figs. 3 and 4). One reason, already mentioned, is that ALMA
maps are not yet perfectly calibrated and some systematic errors are still
present in the data, such as variability in the brightness temperature, scan
pattern visibility, "spokes" on the limb, etc. Possible causes include
variable receiver gain, atmospheric conditions, antenna response, and even
dynamic solar conditions, since the effective exposure time was more than 1
min without calibration scans. Moreover, by convolving the HMI data to the
resolution of ALMA, averaged magnetic field values are decreased, as can be
seen in Fig. 2. The analysis only included LOS magnetic field without
information of the field direction which introduces geometrical effects which
are more pronounced near the limb.

Despite all these effects which were not taken into account, the correlation
is definitely present. In the case of sunspots, the correlation is close to
linear, while active regions show pronounced nonlinearity. A similar negative
correlation for sunspots was observed by ? (see their Fig. 8). They used
RATAN-600 radiotelescope at 6–18 GHz and assumed that the radiation is coming
from a gyroresonance layer above sunspots with the brightness temperature
equal to the kinetic temperature of electrons which depends on the magnetic
field strength and the observed frequency.

The apparent small asymmetry between the positive and negative active regions
visible in Fig. 3 (opposite polarity sunspots also show different behavior)
might be the result of different properties of the ordinary and extraordinary
wave modes, coming preferentially from different polarity regions, but it is
highly uncertain.

We examine two possible radiation mechanisms to explain the observed
correlation. One is gyroresonance, where non-relativistic electrons spiral
along magnetic field lines, radiating along the way. The absorption
coefficient is strongest at harmonic/resonant frequencies of the gyrofrequency
of the electrons and quickly decreases elsewhere (?). Also, the contribution
of higher resonant modes is negligible (?). For the ALMA band 6 frequency of
248 GHz, a field of $\approx 30$ kG is needed for gyroresonance emission at
3rd harmonic, which is an order of magnitude higher than the observed field in
sunspots. This leaves us with the free-free emission as the most probable
explanation, where radiation is a result of free electrons interacting with
ions in the plasma. In this case, magnetic field can change the local plasma
properties, as is the case in sunspots where it inhibits convection. In active
regions, the main source of millimeter free–free opacity are plages,
chromospheric regions with higher temperature and density than the surrounding
quiet Sun regions. The free–free millimeter intensity is approximately
proportional to the electron temperature, so an ALMA image practically maps
the temperature variations of the chromosphere.

##  4 Summary and conclusion

We analyzed the correlation between the photopsheric LOS magnetic field and
the ALMA brightness temperature by comparing the HMI magnetogram and the ALMA
single-dish image of the Sun. The analysis was performed on five subgroups of
data based on their magnetic field strength. No correlation was observed for
the quiet Sun regions, strong negative correlation was found for both positive
and negative polarity sunspots, while network and active regions show a
positive correlation. The character of the observed correlation seems to
change at $|B_{LOS}|\approx 1$ G (no correlation to positive one) and $\approx
120$ G (positive to negative). The probable explanation is that observed
millimeter radiation is mostly coming from the thermal bremsstrahlung of solar
plasma, where the magnetic field affects the temperature and density of the
plasma and thereby modulates the radiated intensity. Measured field strengths
are far too weak for the gyroresonance to contribute significantly to the
free-free opacity.

In this comparison between the LOS magnetic field and the brightness
temperature, we have removed the limb brightening profile from the ALMA image,
which is a questionable step because it removed the $T_{b}$ dependency with
the line-of-sight angle. However, we also used a radial limit for the data, so
regions near the limb with a significant brightening would be removed anyway
from the analysis, even without the limb-brightening correction. The effect of
limb brightening on the results will be a topic of further research on the
subject.

## Acknowledgements

This work has been supported by the Croatian Science Foundation under the
project 7549 "Millimeter and submillimeter observations of the solar
chromosphere with ALMA". This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2011.0.00020.SV. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. We acknowledge the use of the ALMA Solar Ephemeris Generator (?) for
preparation of solar ALMA observations. This research used version 2.1.2 (?)
of the SunPy open source software package (?). HMI image is courtesy of
NASA/SDO and the AIA, EVE, and HMI science teams.

## Bibliography

  *   * Alissandrakis, C. E.: 2020, Frontiers in Astronomy and Space Sciences 7, 74.
  * Bastian, T. S., Bárta, M., Brajša, R., Chen, B., Pontieu, B. D., Gary, D. E., Fleishman, G. D., Hales, A. S., Iwai, K., Hudson, H., Kim, S., Kobelski, A., Loukitcheva, M., Shimojo, M., Skokić, I., Wedemeyer, S., White, S. M., and Yan, Y.: 2018, The Messenger 171, 25–30.
  * Brajša, R., Romštajn, I., Wöhl, H., Benz, A. O., Temmer, M., and Roša, D.: 2009, A&A 493(2), 613–621.
  * Brajša, R., Sudar, D., Benz, A. O., Skokić, I., Bárta, M., De Pontieu, B., Kim, S., Kobelski, A., Kuhar, M., Shimojo, M., Wedemeyer, S., White, S., Yagoubov, P., and Yan, Y.: 2018, A&A 613, A17.
  * Iwai, K., Loukitcheva, M., Shimojo, M., Solanki, S. K., and White, S. M.: 2017, ApJ 841(2), L20.
  * Korzhavin, A. N., Opeikina, L. V., and Peterova, N. G.: 2010, Astrophysical Bulletin 65(1), 60–74.
  * Lemen, J. R., Title, A. M., Akin, D. J., Boerner, P. F., Chou, C., Drake, J. F., Duncan, D. W., Edwards, C. G., Friedlaender, F. M., Heyman, G. F., Hurlburt, N. E., Katz, N. L., Kushner, G. D., Levay, M., Lindgren, R. W., Mathur, D. P., McFeaters, E. L., Mitchell, S., Rehse, R. A., Schrijver, C. J., Springer, L. A., Stern, R. A., Tarbell, T. D., Wuelser, J.-P., Wolfson, C. J., Yanari, C., Bookbinder, J. A., Cheimets, P. N., Caldwell, D., Deluca, E. E., Gates, R., Golub, L., Park, S., Podgorski, W. A., Bush, R. I., Scherrer, P. H., Gummin, M. A., Smith, P., Auker, G., Jerram, P., Pool, P., Soufli, R., Windt, D. L., Beardsley, S., Clapp, M., Lang, J., and Waltham, N.: 2012, Sol. Phys. 275(1-2), 17–40.
  * Loukitcheva, M.: 2020, Frontiers in Astronomy and Space Sciences 7, 45.
  * Loukitcheva, M. A., White, S. M., and Solanki, S. K.: 2019, ApJ 877(2), L26.
  * Mumford, S. J., Freij, N., Christe, S. et al.: 2021, Sunpy.
URL: https://doi.org/10.5281/zenodo.4580466

  * Nindos, A.: 2020, Frontiers in Astronomy and Space Sciences 7, 57.
  * Nindos, A., Alissandrakis, C. E., Bastian, T. S., Patsourakos, S., De Pontieu, B., Warren, H., Ayres, T., Hudson, H. S., Shimizu, T., Vial, J. C., Wedemeyer, S., and Yurchyshyn, V.: 2018, A&A 619, L6.
  * Pesnell, W. D., Thompson, B. J., and Chamberlin, P. C.: 2012, Sol. Phys. 275(1-2), 3–15.
  * Rodger, A. S., Labrosse, N., Wedemeyer, S., Szydlarski, M., Simões, P. J. A., and Fletcher, L.: 2019, ApJ 875(2), 163.
  * Schou, J., Scherrer, P. H., Bush, R. I., Wachter, R., Couvidat, S., Rabello-Soares, M. C., Bogart, R. S., Hoeksema, J. T., Liu, Y., Duvall, T. L., Akin, D. J., Allard, B. A., Miles, J. W., Rairden, R., Shine, R. A., Tarbell, T. D., Title, A. M., Wolfson, C. J., Elmore, D. F., Norton, A. A., and Tomczyk, S.: 2012, Sol. Phys. 275(1-2), 229–259.
  * Selhorst, C. L., Simões, P. J. A., Brajša, R., Valio, A., Giménez de Castro, C. G., Costa, J. E. R., Menezes, F., Rozelot, J. P., Hales, A. S., Iwai, K., and White, S.: 2019, ApJ 871(1), 45.
  * Skokić, I., and Brajša, R.: 2019, Rudarsko-geološko-naftni zbornik 34(2), 59–65.
  * Sudar, D., Brajša, R., Skokić, I., and Benz, A. O.: 2019a, Sol. Phys. 294(11), 163.
  * Sudar, D., Skokić, I., and Brajša, R.: 2019b, Central European Astrophysical Bulletin 43, 1–8.
  * The SunPy Community, Barnes, W. T., Bobra, M. G., Christe, S. D., Freij, N., Hayes, L. A., Ireland, J., Mumford, S., Perez-Suarez, D., Ryan, D. F., Shih, A. Y., Chanda, P., Glogowski, K., Hewett, R., Hughitt, V. K., Hill, A., Hiware, K., Inglis, A., Kirk, M. S. F., Konge, S., Mason, J. P., Maloney, S. A., Murray, S. A., Panda, A., Park, J., Pereira, T. M. D., Reardon, K., Savage, S., Sipőcz, B. M., Stansby, D., Jain, Y., Taylor, G., Yadav, T., Rajul, Dang, T. K., and Contributors, S.: 2020, The Astrophysical Journal 890(1), 68.
URL: https://doi.org/10.3847/1538-4357/ab4f7a

  * Wedemeyer, S., Bastian, T., Brajša, R., Hudson, H., Fleishman, G., Loukitcheva, M., Fleck, B., Kontar, E. P., De Pontieu, B., Yagoubov, P., Tiwari, S. K., Soler, R., Black, J. H., Antolin, P., Scullion, E., Gunár, S., Labrosse, N., Ludwig, H. G., Benz, A. O., White, S. M., Hauschildt, P., Doyle, J. G., Nakariakov, V. M., Ayres, T., Heinzel, P., Karlicky, M., Van Doorsselaere, T., Gary, D., Alissandrakis, C. E., Nindos, A., Solanki, S. K., Rouppe van der Voort, L., Shimojo, M., Kato, Y., Zaqarashvili, T., Perez, E., Selhorst, C. L., and Barta, M.: 2016, Space Sci. Rev. 200(1-4), 1–73.
  * White, S. M., Iwai, K., Phillips, N. M., Hills, R. E., Hirota, A., Yagoubov, P., Siringo, G., Shimojo, M., Bastian, T. S., Hales, A. S., Sawada, T., Asayama, S., Sugimoto, M., Marson, R. G., Kawasaki, W., Muller, E., Nakazato, T., Sugimoto, K., Brajša, R., Skokić, I., Bárta, M., Kim, S., Remijan, A. J., de Gregorio, I., Corder, S. A., Hudson, H. S., Loukitcheva, M., Chen, B., De Pontieu, B., Fleishmann, G. D., Gary, D. E., Kobelski, A., Wedemeyer, S., and Yan, Y.: 2017, Sol. Phys. 292(7), 88.