# Exploiting Vagueness for Multi-Agent Consensus

Michael Crosscombe  Jonathan Lawry
Department of Engineering Mathematics
University of Bristol
BS8 1UB  United Kingdom
<EMAIL_ADDRESS>$\cdot$<EMAIL_ADDRESS>
###### Abstract

A framework for consensus modelling is introduced using Kleene’s three valued
logic as a means to express vagueness in agents’ beliefs. Explicitly
borderline cases are inherent to propositions involving vague concepts where
sentences of a propositional language may be _absolutely true_ , _absolutely
false_ or _borderline_. By exploiting these intermediate truth values, we can
allow agents to adopt a more vague interpretation of underlying concepts in
order to weaken their beliefs and reduce the levels of inconsistency, so as to
achieve consensus. We consider a consensus combination operation which results
in agents adopting the borderline truth value as a shared viewpoint if they
are in direct conflict. Simulation experiments are presented which show that
applying this operator to agents chosen at random (subject to a consistency
threshold) from a population, with initially diverse opinions, results in
convergence to a smaller set of more precise shared beliefs. Furthermore, if
the choice of agents for combination is dependent on the payoff of their
beliefs, this acting as a proxy for performance or usefulness, then the system
converges to beliefs which, on average, have higher payoff.

Keywords: Agent-Based Modelling $\cdot$ Many-Valued Logics $\cdot$ Belief
Aggregation $\cdot$ Consensus

##  1 Introduction

Reaching a consensus by agreeing a shared viewpoint or position is a
fundamental part of many multi-agent decision making and negotiation
scenarios. In this paper we argue that by exploiting vagueness in the form of
explicitly _borderline_ cases we can define an operator for belief combination
which not only allows a population of agents to reach consensus but also
results in them adopting, on average, a more _useful_ set of beliefs. The
basic intuition underlying this operator is that conflicting agents can agree
to allocate borderline truth values to propositions about which they hold
inconsistent beliefs. For example, two individuals, one of which believes that
‘Cameron is an effective prime minister’ whilst the other believes that
‘Cameron is ineffective’, may agree, in some circumstances, to adopt the
shared view that ‘Cameron is borderline effective/ineffective’.

Of course, beliefs about the world do not exist in isolation but inform and
influence our decisions and actions. From this perspective, some sets of
beliefs are more positive or useful than others, resulting in better long term
performance, perhaps by making the individuals concerned richer, happier or
just better able to survive. More generally, in a multi-agent context,
different beliefs result in different actions, collecting different payoffs.
In this paper we present simulation studies which show that implementing our
proposed operator across a population of agents, initially holding diverse
beliefs, results in convergence to a smaller subset of more precise shared
opinions. Furthermore, under the assumption that better performing agents,
i.e. those with higher payoff, are more likely to interact to reach consensus,
we show that the range of beliefs obtained at steady state are on average
better, i.e have higher payoff, than the agents’ initial beliefs. The
formalism adopted here is that of Kleene’s three valued logic and the operator
investigated has been proposed for single propositions in [10] and extended to
multi-propositional languages in [6].

An outline of the paper is as follows: Section 2 gives a brief overview of
consensus modelling. Section 3 introduces Kleene logic and the three valued
consensus combination operator. Section 4 describes simulation experiments in
which agents are selected at random to form a consensus provided that they are
sufficiently consistent with one another. In section 5 we introduce a payoff
function for beliefs, so that the payoff of a particular set of beliefs acts
as a proxy for the performance of an agent holding those beliefs. We then
adapt the experiments described in section 4 so that the probability of an
agent being selected for consensus is proportional to their payoff. Finally,
in section 6 we give some discussions and conclusions.

##  2 Background and Related Work

A number of models for consensus have been proposed in the literature which
have influenced the development of the framework described in this paper. [3]
introduced a model for reaching a consensus involving a weighted, global
updating of beliefs, iterating until an agreement is reached. In DeGroot’s
model, agents assign a weight distribution to the population before forming a
new opinion. By applying their assigned weights to the other agents’ beliefs,
an agent can control the influence that others have on their own beliefs.

As an alternative to DeGroot’s model, the Bounded Confidence (BC) model
introduced in [5] provides agents with a confidence measure. An agent
quantifies their level of confidence in their own opinions and are then able
to limit their interactions to those agents who possess similar beliefs if
they are highly confident (small bounds), or extend the range of possible
interactions if the agents possess low confidence (large bounds). In this
model agents do not a priori assign weights to the beliefs of others, but
instead determine such weightings based on similarity and on their own
confidence levels. This is similar in essence to the inconsistency threshold
that we introduce in section 3, but applied on an individual basis.

The Relative Agreement (RA) model [2] then extends the Bounded Confidence
model to allow agents to assign weights to the beliefs of others by
quantifying the extent of the overlap of their respective confidence bounds.
By having agents declare a confidence interval for their beliefs, the model
then restricts interactions to those pairs of agents with overlapping
intervals. Consequently, agents are only required to assess their own beliefs
and are not required to make explicit judgements about those of other agents.
[2] also moved to a model of pair-wise interactions to better capture social
interactions of individuals, the latter being a setting in which group-wide
updates to beliefs are unintuitive in that they do not reflect typical social
behaviour.

A fundamental difference between our approach and the above models is that we
use Kleene’s three valued logic to represent beliefs in a propositional logic
setting, rather than identify opinions with real values or intervals. [10]
have shown that through use of a three-state model for networked consensus of
complete graphs, nodes converge to a consensus much faster and with greater
accuracy when compared to a restrictive binary model. In the sequel we extend
this approach to a more general setting involving larger languages and
incorporating a measure of payoff for beliefs.

##  3 A Three Valued Consensus Model

In this section we introduce Kleene’s three valued logic [4] as a model of
explicitly borderline cases resulting from the inherent vagueness of
propositions. We adopt a propositional logic setting as follows: Let ${\cal
L}$ be a finite language of propositional logic with connectives $\wedge$,
$\vee$ and $\neg$, and propositional variables ${\cal
P}=\\{p_{1},\ldots,p_{n}\\}$. Also, let $S{\cal L}$ denote the sentences of
${\cal L}$ generated by recursively applying the connectives to the
propositional variables in the usual manner. A Kleene valuation then allocates
truth values $0$ (false), $\frac{1}{2}$ (borderline) and $1$ (true) to the
sentences of ${\cal L}$ as follows:

######  Definition 1.

Kleene Valuations

A Kleene valuation $v$ on $\mathcal{L}$ is a function
$v:S\mathcal{L}\rightarrow\\{0,\frac{1}{2},1\\}$ such that
$\forall\theta,\varphi\in S\mathcal{L}$ the following hold:

  * •

$v(\neg\theta)=1-v(\theta)$

  * •

$v(\theta\land\varphi)=$ min$(v(\theta),v(\varphi))$

  * •

$v(\theta\lor\varphi)=$ max$(v(\theta),v(\varphi))$

The truth table for Kleene valuations are shown in table 1.

$\neg$ | $1$ | $0$
---|---|---
| $\frac{1}{2}$ | $\frac{1}{2}$
| $0$ | $1$

$\wedge$ | $1$ | $\frac{1}{2}$ | $0$
---|---|---|---
$1$ | $1$ | $\frac{1}{2}$ | $0$
$\frac{1}{2}$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$
$0$ | $0$ | $0$ | $0$

$\vee$ | $1$ | $\frac{1}{2}$ | $0$
---|---|---|---
$1$ | $1$ | $1$ | $1$
$\frac{1}{2}$ | $1$ | $\frac{1}{2}$ | $\frac{1}{2}$
$0$ | $1$ | $\frac{1}{2}$ | $0$

Table 1: Kleene truth tables.

It is sometimes convenient to represent a Kleene valuation $v$ by its
associated _orthopair_ $(P,N)$ [6], where $P=\\{p_{i}\in{\cal
P}:v(p_{i})=1\\}$ and $N=\\{p_{i}\in{\cal P}:v(p_{i})=0\\}$. Notice that
$P\cap N=\emptyset$ and that $(P\cup N)^{c}$ corresponds to the set of
borderline propositional variables.

Kleene valuations have been proposed as a suitable formalism in which to
capture explicitly borderline cases as resulting from inherent flexibility in
the definition of vague concepts in natural language [8, 7]. For example,
consider the proposition ‘Ethel is short’. For the concept short, we might
identify a lower height threshold $\underline{h}$ below which any height is
classed as being _absolutely short_ , and similarly there may be an upper
threshold $\overline{h}$ above which any height is _absolutely not short_. If
Ethel’s height lay between $\underline{h}$ and $\overline{h}$ then this would
result in a borderline truth value for the statement ‘Ethel is short’.

It is important to note that the middle truth value $\frac{1}{2}$ is not
intended to represent epistemic uncertainty, but rather explicitly borderline
cases resulting from the inherent vagueness of natural language propositions.
Hence, if we say that the statement ‘Ethel is short’ is borderline true/false
we are not saying that the truth or falsity of this proposition is unknown.
Instead we are indicating that Ethel’s height is a borderline case of the
predicate short. In order to emphasise the difference between the epistemic
and the borderline interpretation of three valued logic it is helpful to think
in terms of conditioning. For instance, if we learn that it is _unknown_
whether or not Ethel is short, then this provides us with no new information
about her height. In contrast, learning that Ethel is _borderline short_ does
provide us with new information about Ethel’s height, namely that it lies on
the borderline between short and not short. A more comprehensive discussion of
these issues can be found in [1]. A consequence of using this interpretation
of the middle truth value is that in the current paper we only model consensus
for sets of propositions which admit borderline cases. In other words, our
approach can be used for propositions such as ‘Ethel is short’ but not, for
example, for the proposition ‘Ethel is strictly less that 1.4 metres tall’.

The following three valued consensus operator was described in detail in [6]:

######  Definition 2.

Consensus Operator

Let $v_{1}$ and $v_{2}$ be Kleene valuations on ${\cal L}$ with associated
orthopairs $(P_{1},N_{1})$ and $(P_{2},N_{2})$. Then the consensus $v_{1}\odot
v_{2}$ is the Kleene valuation with the orthopair

$\displaystyle((P_{1}\scalebox{0.5}[1.0]{$-$}N_{2})\cup(P_{2}\scalebox{0.5}[1.0]{$-$}N_{1}),(N_{1}\scalebox{0.5}[1.0]{$-$}P_{2})\cup(N_{2}\scalebox{0.5}[1.0]{$-$}P_{1}))$

The corresponding truth table for this operator is shown in table 3. From this
we can see that the operator preserves the non-borderline truth values $0$ or
$1$ except in the case of a direct conflict i.e. when one agent has truth
value $1$ and the other $0$. In this case both agents adopt the middle truth
value $\frac{1}{2}$. Alternatively, from definition 2 we can think of $\odot$
as an operator which initially weakens both opinions so as to remove direct
inconsistencies, before then combining them.

$\odot$ | $1$ | $\frac{1}{2}$ | $0$
---|---|---|---
$1$ | $1$ | $1$ | $\frac{1}{2}$
$\frac{1}{2}$ | $1$ | $\frac{1}{2}$ | $0$
$0$ | $\frac{1}{2}$ | $0$ | $0$
Table 2: Truth table for the consensus operator.

$I$ | $1$ | $\frac{1}{2}$ | $0$
---|---|---|---
$1$ | $0$ | $0$ | $1$
$\frac{1}{2}$ | $0$ | $0$ | $0$
$0$ | $1$ | $0$ | $0$
Table 3: Inconsistency truth table.

We now introduce two measures that will be used throughout the subsequent
simulation experiments.

######  Definition 3.

A Measure of Vagueness

Let $v$ be a Kleene valuation on ${\cal L}$ with orthopair $(P,N)$ and $n$
propositional variables. Then we measure the vagueness of $v$ by the
proportion of propositional variables which it classifies as being borderline.
That is:

$\displaystyle V(v)=\frac{\left|(P\cup N)^{c}\right|}{n}$

######  Definition 4.

Inconsistency Measure

Let $v_{1}$ and $v_{2}$ be Kleene valuations on ${\cal L}$ with corresponding
orthopairs $(P_{1},N_{1})$ and $(P_{2},N_{2})$. Then we define the
inconsistency measure of $v_{1}$ and $v_{2}$ to be the proportion of
propositional variables which are in direct conflict between the two
valuations i.e. $v_{1}(p_{i})\neq\frac{1}{2}$, $v_{2}(p_{i})\neq\frac{1}{2}$
and $v_{1}(p_{i})=1-v_{2}(p_{i})$. That is:

$\displaystyle I(v_{1},v_{2})=\frac{\left|(P_{1}\cap
N_{2})\right|+\left|(P_{2}\cap N_{1})\right|}{n}$

Table 3 shows the inconsistency truth table of two valuations for a
propositional variable, highlighting the cases where two valuations are
inconsistent, and consistent otherwise. We can see that there is a probability
of $\frac{2}{9}$ that two valuations will be inconsistent for each
propositional variable in the language. In the sequel we will propose a
threshold $\gamma\in[0,1]$ on inconsistency so that valuations $v_{1}$ and
$v_{2}$ can be combined only if $I(v_{1},v_{2})\leq\gamma$.

##  4 Simulation Experiments based on Random Selection of Agents

Figure 1:  Average vagueness after $50,000$ for varying inconsistency
thresholds $\gamma$ and language sizes.

Figure 2:  Number of distinct valuations after $50,000$ iterations for varying
inconsistency thresholds $\gamma$ and language sizes.

We introduce simulation experiments in order to investigate the convergence
properties of the three valued logic operator when implemented across a multi-
agent system. The experimental set up is loosely based on those proposed in
[2] and [9], although our representation of opinions is quite different with
beliefs taking the form of Kleene valuations on ${\cal L}$, rather than
vectors of bounded real numbers.

We will consider two distinct initialisations of the beliefs of a population
of agents. The _random three valued initialisation_ allocates the truth values
$0$, $\frac{1}{2}$ and $1$ to each agent and each propositional variable at
random i.e. with probability $\frac{1}{3}$ for each truth value. In contrast,
the _random Boolean initialisation_ only allocates the binary truth values $0$
and $1$, each with a probability of $\frac{1}{2}$. This latter initialisation
will be required in section 5 in order to directly compare the proposed three
valued combination operator with a similar two valued operator. In this
section we will use the random three valued initialisation in order to
investigate the extent to which the three valued operator results in
convergence to a shared set of opinions across the population of agents.

We set a fixed maximum number of $50,000$ iterations111In preliminary
experiments we found that $50,000$ was an upper bound on the number of
iterations required for the system to reach steady state across a range of
parameter settings.. At each time step a pair of agents are selected at random
from the population. An inconsistency threshold value $\gamma\in[0,1]$ is set,
so that for any pair of agents with respective valuations $v_{1}$ and $v_{2}$,
if $I(v_{1},v_{2})\leq\gamma$ then both agents replace their beliefs with the
consensus valuation $v_{1}\odot v_{2}$, while if $I(v_{1},v_{2})>\gamma$ then
no combination is performed and both agents retain their original beliefs. For
$\gamma=1$ we obtain what is equivalent to the totally connected graph model
described in [10], in which any pair of agents can combine their beliefs,
whilst taking $\gamma=0$ corresponds to the most conservative scenario in
which only absolutely consistent beliefs can be combined. The parameters for
the simulation experiments are then as follows:

  * •

Population size: $100$

  * •

Language size i.e $|{\cal P}|=n$: $5,10,50,100$

  * •

Initial beliefs: Random three valued.

  * •

Inconsistency threshold: $\gamma\in[0,1]$.

Figures 2 and 2 show the results for the experiments after $50,000$
iterations. In each case the plots show mean values with error bars
representing standard deviation across $100$ independent runs of the
simulation. Figure 2 shows the average vagueness determined by taking the mean
value of $V(v)$ (definition 3) across the population. Note that for a random
three valued initialisation of beliefs we expect a mean vagueness value of
$\frac{1}{3}$ at the start of the simulation. As the threshold $\gamma$
increases then the average vagueness decreases to zero, so that for
$\gamma\geq 0.3$ we are left with almost entirely crisp (i.e. Boolean)
opinions. In general the more conservative the combination rules (i.e.
requiring higher levels of consistency) then the more it is that vague beliefs
are maintained in the population. Figure 2 shows the number of distinct
valuations (i.e. different opinions) remaining in the population after
$50,000$ iterations. Again this decreases with $\gamma$ and for $\gamma>0.4$
agents have on average converged to a single shared belief. This is consistent
with the analytical results presented in [10] for the single propositional,
$\gamma=1$ case.

##  5 Simulation Experiments Incorporating a Payoff Model

In this section we extend the simulation framework described in section 4 to
allow for different payoffs for different beliefs. As outlined in section 1,
payoff is introduced as a proxy for performance, and is motivated by the
intuition that different beliefs result in different actions which then, over
time, lead to different levels of performance. Here we adopt an abstract
simplification of this process in which each Kleene valuation is allocated a
real valued payoff. Then, instead of being selected at random for combination,
an agent is picked from the population according to a probability which is
proportionate to the payoff value of their beliefs. The idea, then, is that
agents with better or more useful opinions will be more successful and
furthermore, it will be these successful agents who will be most likely to
need to reach a consensus between them.

Here the underlying intuition is that, in real systems it is the most
successful agents, with the highest payoff values, who are most likely to find
themselves in conflict with one another, and who will most benefit from
reaching an agreement. We adopt a simple summative payoff model in which each
propositional variable $p_{i}$ is allocated a value in the range $[-1,1]$,
denoted $f(p_{i})$, and the payoff for a valuation with orthopair $(P,N)$ is
then calculated as follows:

$\displaystyle f(P,N)=\sum_{p_{i}\in P}f(p_{i})-\sum_{p_{i}\in N}f(p_{i})$

Another perspective on this type of payoff function is as follows: For each
propositional variable $p_{i}$, a truth value of $1$ results in a payoff
$f(p_{i})$ (which can be either positive or negative), a truth value of $0$
results in the opposite signed payoff $\scalebox{0.5}[1.0]{$-$}f(p_{i})$, and
a borderline truth value $\frac{1}{2}$ results in a neutral payoff of $0$. The
payoff value for a Kleene valuation $v$ is then simply taken to be the sum of
the payoffs for each propositional variable under the truth values allocated
by $v$.

Figure 3:  Average payoff after $50,000$ iterations for varying inconsistency
thresholds $\gamma$, shown as a percentage of the maximal payoff.

Figure 4:  Number of distinct valuations after $50,000$ iterations for varying
inconsistency thresholds $\gamma$.

Based on payoff values we define a probability distribution over the agents in
the population according to which the probability that an agent with beliefs
$(P,N)$ is selected for possible consensus combination is proportional to
$f(P,N)+n$. At each iteration a pair of agents are selected at random
according to this distribution. For each such pair the inconsistency measure
(definition 4) is evaluated and either both the valuations are replaced with
the consensus valuation, or both are left unchanged, depending on the
threshold $\gamma$ as in section 4. The parameters for the simulation
experiments are as follows:

  * •

Population size: 100

  * •

Language size: 5

  * •

Initial beliefs: Random Boolean.

  * •

Inconsistency threshold: $\gamma\in[0,1]$

binary operator | $0$ | $1$
---|---|---
$0$ | $0$ | $0:\frac{1}{2},1:\frac{1}{2}$
$1$ | $0:\frac{1}{2},1:\frac{1}{2}$ | $1$
Table 4: The truth table for the stochastic Boolean consensus operator

Notice that here we are initialising the beliefs as random Boolean valuations
(see section 4)222As a result of this Boolean initialisation, a language size
of $5$ now produces a total of $2^{5}$ ($32$) possible valuations, as opposed
to $3^{5}$ ($243$) possible valuations.. This allows us to make a direct
comparison between the performance of the three valued combination operator
and a similar two valued operator. For the latter we assume that only binary
truth values are available to represent an agent’s beliefs. In this context,
in order for two agents with conflicting truth values for $p_{i}$ (i.e. one
$0$ and the other $1$) to reach consensus, we propose that they simply agree
to pick one of the truth values at random e.g. by tossing a fair coin. Table 4
gives the truth table for the operator in which directly conflicting truth
values leads to a stochastic outcome.

The focus on simulations with $5$ propositional variables is intended to
increase the number of opinions relative to the size of the population, in
order to achieve a good distribution of valuations. For example, a language
size of $5$ allows for $32$ possible Boolean valuations. With a population of
$100$ agents, it is therefore very likely that each opinion will occur at
least once. In comparison, a language size of $10$ produces $1,024$ possible
Boolean valuations which severely decreases the probability of an opinion
being present in a population of the same size.

Figures 4, 4 and 5 show the results for simulation experiments with agent
selection based on payoff. The results shown are mean values with error bars
taken over $100$ independent runs of the simulation. Figure 4 shows the
average population payoff after $50,000$ iterations given as a percentage of
the maximal possible payoff value i.e. the payoff for the valuation $(P,N)$
where $P=\\{p_{i}:f(p_{i})>0\\}$ and $N=\\{p_{i}:f(p_{i})<0\\}$. For both the
binary and the three valued operators we show results for simulations in which
agents are selected according to payoff (three-valued, Boolean) and at random
as in section 4 (three-valued random, Boolean random). We see that for all
values of $\gamma$, the three valued operator with payoff based selection
outperforms all of the other approaches. For the former we can also see that
average payoff increases with $\gamma$. In contrast, for the other approaches,
including the payoff operator with payoff based selection, the mean of the
average population payoff remains close to $0$ after $50,000$ iterations.
Figure 4 shows the mean number of distinct valuations across the population of
agents after $50,000$ iterations. All four versions of the operators converge
on a small set of shared beliefs for sufficiently large $\gamma$. For
$\gamma\geq 0.4$ the mean number of distinct valuations is less than $5$ while
for $\gamma\geq 0.8$ it is $1$. Figure 5 shows a trajectory of how the number
of distinct valuations varies with each iteration when $\gamma=0.7$. We can
see that both the three-valued models converge quickly (in just over $2000$)
iterations while the Boolean models require considerably longer to converge
(over $20,000$ iterations).

Figure 5: Trajectory showing the number of distinct valuations plotted against
iterations for $\gamma=0.7$.

##  6 Conclusions

In this paper we have explored the use of Kleene’s three valued logic as a
framework in which to model multi-agent consensus formation. We have proposed
a three valued combination operator, the intuition behind which is that
conflicting binary truth values are replaced with a borderline (middle) truth
value. A number of simulation experiments have been presented employing this
operator. These can be divided into two main categories. For the first type of
experiments, agents are selected at random from the population and form a
consensus valuation providing that the level of inconsistency of their
respective opinions is below a threshold parameter $\gamma$. Otherwise they do
not form a consensus and instead retain their current opinions. For these
experiments we found that there is convergence to a smaller subset of shared
opinions across the population. For higher $\gamma$ values there is
convergence on average to a single shared opinion and furthermore this opinion
is crisp i.e. it admits no borderlines. For intermediate values of $\gamma$
the system convergences to a small set of opinions which to some extent remain
vague.

In the second type of experiments a payoff function over beliefs is
introduced, and agents are selected for possible combination with probability
proportional to the payoff value of their current beliefs. Here we compare the
three value operator with a similar stochastic Boolean operator. We find that
the three valued operator with payoff based agent selection results in
convergence to a smaller shared set of beliefs with significantly higher
average payoff than that of the initial population. The Boolean operator does
not perform well in this context and does not result in a significant increase
in average payoff, which instead remains close to $0$ after $50,000$
iterations.

The results of the payoff based experiments show how a three valued model for
consensus provides a number of improvements over a traditional Boolean model.
Firstly, we have shown that the introduction of Kleene valuations to capture
the inherent vagueness of propositions does not, in the long run, lead to the
mass adoption of borderline truth values as a result of conflict occurring in
the population. Instead, we have seen how vagueness is reduced at lower
$\gamma$ values, and at higher $\gamma$ values the population converges
towards completely crisp opinions on average, admitting no borderline cases.
In addition to this, we can see that the introduction of a payoff based model
drives consensus towards those valuations which result in higher payoff on
average. By selecting pairs of agents based on their perceived success, we can
achieve an increase to overall payoff in a small number of iterations,
compared to no significant increase in payoff for the Boolean model.
Therefore, we have shown that the three valued approach incorporating a payoff
model can drive convergence across the population towards more successful
opinions.

We suggest that the experiments presented in this paper show the potential of
using three valued logic in consensus modelling. There is also significant
scope to extend the research presented in several new directions. For example,
the above studies concern consensus defined at the level of propositional
variables. However, in many cases agents will be most concerned to reach
agreement about a relevant set of compound statements. For example, they may
need to reach agreement about a particular set of conditional statements, or
equivalences. Hence, an important question is that of how best to extend our
proposed consensus model so as to be applicable to compound logical
expressions. Another significant question concerns uncertainty. Suppose that
in addition to vagueness agents also quantify their uncertainty about beliefs.
[6] propose an extension of the three valued framework in which agents’
beliefs are represented by a probability distribution over Kleene valuations.
Ongoing research concerns the design of simulation studies in which to
evaluate the convergence and payoff based performance of this extended model.
Finally, it would be interesting to consider extensions to the operator which
allows for consensus between groups rather than just pairs of agents.

## Acknowledgements

This research is partially funded by an EPSRC PhD studentship as part of a
doctoral training partnership (grant number EP/L504919/1).

All underlying data is included in full within this paper.

## References

  * [1] D. Ciucci, D. Dubois, and J. Lawry.  Borderline vs. unknown: comparing three-valued representations of imperfect information.  In International Journal of Approximate Reasoning. 55. pp 1866-1889. Elsevier, 2014.
  * [2] G. Deffuant, F. Amblard, G. Weisbuch, and T. Faure.  How can extremism prevail? a study based on the relative agreement interaction model.  In Journal of Artificial Societies and Social Simulation. 5(4). JASSS, 2002.
  * [3] M. H. DeGroot.  Reaching a consensus.  In Journal of the American Statistical Association. 69(345). pp 118-121. JSTOR, 1974.
  * [4] S. C. Kleene.  Introduction to Metamathematics.  North-Holland, The Netherlands, 1st edition, 1952.
  * [5] U. Krause.  A discrete nonlinear and non-autonomous model of consensus formation.  In Communications in Difference Equations: Proceedings of the Fourth International Conference on Difference Equations, pp 227-237. Gordon and Breach, 1998.
  * [6] J. Lawry and D. Dubois.  A bipolar framework for combining beliefs about vague propositions.  In Proceedings of the Thirteenth International Conference on the Principles of Knowledge Representation and Reasoning. pp 530-540. AAAI, 2012\.
  * [7] J. Lawry and I. González-Rodríguez.  A bipolar model of assertability and belief.  In International Journal of Approximate Reasoning. 52(1). pp 76-91. Elsevier, 2011.
  * [8] J. Lawry and Y. Tang.  On truth-gaps, bipolar belief and the assertability of vague propositions.  In Artificial Intelligence. 191. pp 20-41. Elsevier, 2012.
  * [9] M. Meadows and D. Cliff.  Reexamining the relative agreement model of opinion dynamics.  In Journal of Artificial Societies and Social Simulation. 15(4). JASSS, 2012.
  * [10] E. Perron, D. Vasudevan, and M. Vojnovic.  Using three states for binary consensus on complete graphs.  In Proceedings of IEEE Infocom. pp 2527-2535. IEEE Communications Society, 2009.


††thanks: A footnote to the article title

# Wave function realization of a thermal collision model

Uriel Shafir<EMAIL_ADDRESS>The Institute of Chemistry, The
Hebrew University of Jerusalem, Jerusalem 9190401, Israel  Ronnie Kosloff
<EMAIL_ADDRESS>The Institute of Chemistry, The Hebrew University of
Jerusalem, Jerusalem 9190401, Israel

###### Abstract

An efficient algorithm to simulate dynamics of open quantum system is
presented. The method describes the dynamics by unraveling stochastic wave
functions converging to a density operator description. The stochastic
techniques are based on the quantum collision model. Modeling systems dynamics
by wave functions and modeling the interaction with the environment with a
collision sequence reduces the complexity scale significantly. The algorithm
developed, can be implemented on quantum computers. We introduce stochastic
methods that exploit statistical characters of the model, as Markovianity,
Brownian motion and binary distribution. The central limit theorem is employed
to study the convergence of distributions of stochastic dynamics of pure
quantum states represented by wave vectors. By averaging a sample of functions
in the distribution we prove and demonstrate the convergence of the dynamics
to the mixed quantum state described by a density operator.

††preprint: APS/123-QED

##  I introduction

In reality every quantum system is open, while an isolated system is an
exception. Therefore the main setback in simulating and modeling real life
quantum system is the high cost in computation. To analyse the cost, we will
first describe the generic description of an open quantum system, dive into
description of the dynamics and observe the computation problem. After
reviewing the current methods we will demonstrate a scheme able to lower the
computation complexity.

An open quantum system is generically described by the following: the systems’
Hamiltonian and systems density matrix(DM) -$\hat{H}_{s}$, $\hat{\rho}_{s}$
the measuring apparatus - M and the environment Hamiltonian $\hat{H}_{B}$.
Contemporary examples are the IBM and Google quantum machines, the nascent
steps toward quantum computing [1]. If these devices are left alone, the
quantum system -$\hat{\rho}_{s}$, will reach thermal equilibrium with the
extremely cold surrounding temperature $T_{B}$. The device is assembled from
quantum circuits, where the computation output is measured by a measurement
apparatus M. Another example is the NV center in diamond [2]. The systems is
constructed from a Nitrogen impurity adjecent to a negative vacancy in
diamond. The primary quantum system Hamiltonian-$H_{s}$ is a spin triplet. The
neighboring environment is composed of spins such as other nitrogen atoms or
carbon isotopes with nuclear spin (${}^{13}C$) in addition to the fluctuations
of the vacuum phonons. A measurement apparatus is being coupled to the NV
center able to measure changes in population; the measurement apparatus- M
[3].

Open system dynamics address a system interacting with the environment from
the system’s perspective. Different approaches have been employed to construct
reduced descriptions in terms of the system observables. Starting with Bloch a
dynamical derivation based on the weak system bath coupling has led to the
quantum Master equation [4, 5]. An alternative mathematical formulation
employing quantum dynamical semigroups has led to a general structure termed
the Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) equation [6, 7]. Davis has
connected the perturbation derivation to the general structure [8, 9]. Non
Markovian formulations including memory effects have also been suggested [10].

Our mission is to develop a simulation algorithm of an open quantum system. In
this approach the system is viewed from a thermodynamic perspective. Recent
studies by Kosloff and Dann [11] have paved the ground to the conditions that
create consistency between Markovian dynamics and thermodynamic principles in
open quantum systems. The consistency conditions are formulated by a set of
axioms which can be applied to the present study.

The developed quantum simulation should be applicable on both classical and
quantum processors. This means that the method should be based on a wave
functions formulation. On classical computers the wave function description is
computationally preferable to a density operator description.

When simulating quantum dynamics, the computational resources of memory scale
at least polynomially with the size of the Hilbert space. In addition every
degree of freedom in the system increases the size of the Hilbert space
exponentially. For example in a system of spin particles the size of the
Hilbert space is $2^{n}$, where n is the number of particles. Describing the
system by density matrix squares the memory resources of the computation.
Solving for the dynamics the scaling of the number of operations is the
product of the size of the Hilbert space multiplied by the product of
propagation time and energy range. The energy range is also doubled in a
density operator description [12]. Another advantage in a wavefunction
formulation is a pedagogical one. Describing the dynamics of single pure state
gives an intuitive sense of the process.

The starting point to describe the dynamics is a global approach which
includes the system and environment. Typically one constructs a global
Hamiltonian $\hat{H}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}$. Solving for the
combined system is computationally prohibited since it consisted of enormous
amounts of degrees of freedom. To overcome the almost impossible obstacle the
whole setup is divided between system and environment. It is customary to
diagonalize the bath to orthogonal modes either harmonic or composed of an
ensemble of spins [13]. The next step is to obtain an effective reduced
equations of motion for the system where the bath enters implicitly. Assuming
an initial uncorrelated system and bath
$\hat{\rho}=\rho_{S}\otimes\hat{\rho}_{B}$ results in a completely positive
trace preserving dynamical map (CPTP) describing the system propagation
$\hat{\rho}_{S}(t)=\Lambda_{t}\hat{\rho}_{S}(0)$ [14]. Imposing in addition
the condition of Markovianity, Gorini, Kossakowski, Sudarshan and Lindblad
obtained the general form of the Master equation [6, 7]. The GKLS master
equation is describing the dynamics of the open system under the assumption
that the bath is in equilibrium and not effected by the interaction with the
system. The GKLS has become one of the cornerstones of the theory of open
quantum systems. Solving the Master equation is a difficult computational
problem. The state is described by density matrices resulting in a
computational scaling of at least $O(N^{2})$ where $N$ is the size of Hilbert
space. Also, in the GKLS equation we address the bath through its ladder
operators. Finding them is equivalent to diagonalizing the bath- a very
expensive computational task.

To reduce the computational complexity a wavefunction method is desirable. The
algorithm involves a stochastic unraveling of wavefunctions. Stochastic
approaches are currently in use in many fields of quantum dynamics such as
thermal averaging [15, 16] and electronic structure methods [17]. A stochastic
approach has been suggested by Percival and Gisin [18] for unraveling the GKLS
equation. In their approach the GKLS equation was transformed to stochastic
differential equations. This unraveling procedure is non unique. This has the
benefit of the freedom to use the most mathematically convenient choice. The
drawback is that the stochastic wavefunction is not associated with a physical
description. An additional problem of this method is that the dynamics are
formulated by a non linear differential equation. This increase the difficulty
in finding a solution approaching a computational scaling of the size of
Hilbert space squared. A different unraveling approach was developed by Katz,
Torrontegui and Kosloff [19, 20]. The method partitions the environment to a
primary and secondary bath. The primary bath termed surrogate Hamiltonian is
weakly coupled to the system [21, 22]. Stochastically the spins of the primary
bath are refreshed from the secondary bath.

We now concentrate on a method of unraveling wave functions using stochastic
variables modelled on the Collision Model (CM). The quantum CM, first appeared
in 1948 in a paper by Karplus and Schwinger [23] followed by work of J.
Rau[24]. In the 60s, and later on in the 80s, CMs appeared in works on weak
measurements by C. M. Caves and G. J. Milburn [25, 26]. In recent years CM
became more and more popular due to its simplicity and its ability to be
consistent with thermodynamics and the low computational effort in the
description of the bath. [27]. The collision model, models the bath as if it
is composed of sub units- ancillas, with which the system interact with. Thus
there is one basic assumption in the general CM:

  1. 1.

The interaction of the system with the bath is described as an interaction
between the system with a single ancilla from the Environment.

We develop a simple collision model with two additional assumptions regarding
the bath:

  1. 1.

Ancillas do not interact with each other.

  2. 2.

Ancillas are initially uncorrelated.

This collision model can yield a good description of real physical example of
a system in a diluted gas and also the NV center.

This general description of CM requires to address two major complexity
problems:

  1. 1.

In order to describe the interaction of the system with the ancilla we will
have to solve the dynamics of the interaction according to some physical
model. This will require to solve the time dependent Schrödinger or Liouville
equation, a very hard task in general and very expensive for big systems. Even
simpler models to transform only heat between the system and its environment
such as resonance between modes of the system to the ancilla particle result
in diagonalization of the system’s Hamiltonian, which for a system with more
then 10 particles is already computationally expensive.

  2. 2.

The size of the density matrix of a multi particle system grows exponentially
with the number of particles and result in huge matrices very fast.

A significant simplification of the model is obtained by adding the assumption
that the interaction between the system and the bath ancilla is much faster
than the free dynamics of the the system. This enables us to treat the
dissipation and interaction with the bath separately from the internal systems
dynamics.

The complexity is further reduced by an implicit treatment of the bath. We
choose to describe the primary system as a system of coupled spins.
This representation is employed for three main reasons:

  1. 1.

It is computationally inexpensive.

  2. 2.

A system composed of qubits is universal and therefore can simulate other
physical systems.

  3. 3.

Such a system can be implemented on a quantum computer.

The collision model creates a Completely Positive Trace Preserving map. This
process leads to a fixed point of the map. A crucial tool for the simulation
is a partial trace algorithm implemented in a wave functions formalism.

We have developed a stochastic algorithm implemented in wave functions
formulation. The setup is based on the partial trace implementation. In
addition we obtain an intuitive description of the process as an average of
partial measurements of the system. The methods presented in this paper can be
implemented in any case of collision model where the environment is composed
of qubits as in [28], and [29]

$\Ket{\psi}$$Utr_{b}\\{S_{p}\Ket{\psi}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}}\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}\Ket{\psi_{3_{2}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p}\Ket{\psi_{3_{1}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{3_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$$Utr_{b}\\{S_{p}\Ket{\psi}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}}\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}\Ket{\psi_{2_{2}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p}\Ket{\psi_{2_{1}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$$Utr_{b}\\{S_{p}\Ket{\psi}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}}\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}\Ket{\psi_{1_{2}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p}\Ket{\psi_{1_{1}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$
Figure 1:  The unraveling tree: The time evolution of the systems wave
function $\Ket{\psi}$. Three different interactions with an ancilla
$\Ket{\beta_{\theta}}$. Each interaction spawns two wave functions with
different wights. This process is repeated with each ancilla interaction. The
Gaussian distribution exhibits that asymptotically the process obeys the
central limit theorem.

$\Ket{\psi}$$Utr_{b}\\{S_{p}\Ket{\psi}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}}\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{1_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(1-p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{1}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(p)}Utr_{b}\\{S_{p}(\Ket{\psi_{2_{2}}}\otimes{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Ket{\beta_{\theta}}})\\}$
Figure 2:  The reduced unraveling tree: Interaction with a thermal qubit is in
agreement with the definitions of Wiener process thus the tree can be reduced
to a binary tree form.

##  II Implementation

###  II.1 Setup description

The framework of the derivation assumes a unitary evolution generated by the
total Hamiltonian:

$\hat{H}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}~{}~{},$ (1)

composed of the system Hamiltonian $\hat{H}_{S}$ environment Hamiltonian
$\hat{H}_{B}$ and interaction $\hat{H}_{SB}$. We assume $\hbar=1$ all through
the paper.

###  II.2 Representation of the spinor

The model studied is composed of a system of qubits. For such a system the
wave function has dimensions of $1\times 2^{N}$ when N is the number of
qubits. The density matrix representing the system has dimensions of
$2^{N}\times 2^{N}$ . To expand the wave function, we choose a local expansion
constructed by the basis of individual components. The natural choice for
qubits is to construct the wave function in the computational basis as a
linear tensor product of the computation base spanning each qubit space. Each
qubit is represented in the computational basis where
$\Ket{0}=\Ket{\downarrow}$ and $\Ket{1}=\Ket{\uparrow}$. The complete basis
that we will represent the system in is

$\\{\Pi_{i=1}^{N}\otimes\Ket{\delta_{i}}\\}$ (2)

where $\delta_{i}=0,1$

$\\{\Pi_{i=1}^{N}\otimes\Ket{\delta_{i}}\\}=\\{\Ket{00...0},\Ket{00...1},...,\Ket{01...0},\Ket{01...1},\Ket{10...0},\Ket{10...1},\Ket{11...0},\Ket{11...1}\\}$
(3)

The basis is ordered in a raising order of the binary basis.

###  II.3 System and Baths

The multi-qubit systems Hamiltonian $\hat{H_{S}}$ is described as:

$\hat{H_{S}}=\sum_{k}\hat{H_{k}}+\sum_{i,j}\epsilon_{i,j}(\hat{\sigma}_{+_{i,j}}+\hat{\sigma}_{-{i,j}})~{}~{}~{},$
(4)

$\hat{H_{k}}$ is the free k’th particle Hamiltonian, $\epsilon_{i,j}$ is the
interaction coefficient between the i and j qubits and is the interaction
between the i and j particles in the system.
Where N is the number of particles and

$\hat{\sigma}_{+{i,j}}=\hat{I}_{2^{i-1}}\otimes\Ket{0}\bra{1}\otimes\hat{I}_{2^{j-i-1}}\otimes\Ket{1}\bra{0}\otimes\hat{I}_{2^{N-j}}$
$\hat{\sigma}_{-{i,j}}=\hat{I}_{2^{i-1}}\otimes\Ket{1}\bra{0}\otimes\hat{I}_{2^{j-i-1}}\otimes\Ket{0}\bra{1}\otimes\hat{I}_{2^{N-j}}$

We represent the state of the system as the density operator $\hat{\rho}_{S}$.
The density matrix of the bath $\rho_{B}$ is composed of uncorrelated ancilla
qubits.

$\hat{\rho}_{B}=\hat{\rho}_{b_{1}}\otimes\hat{\rho}_{b_{2}}\otimes...\otimes\hat{\rho}_{b_{n}}=\Pi_{i}\otimes\hat{\rho}_{b_{i}}=\Pi_{i}\otimes\frac{e^{-\beta\hat{H}_{b_{i}}}}{Z}$
(5)

$\hat{H}_{b}$ is the Hamiltonian of the individual ancilla qubit and $\beta$
is the inverse temperature times the Boltzman factor $\beta=\frac{1}{k_{B}T}$.
The fact that the bath is uncorrelated to the system initially is with concent
with Kosloff and Dann postulate 2 [11].

####  II.3.1 Observables

An observable $\langle\mathcal{O}\rangle$ is defined as
$tr\\{\hat{\mathcal{O}}\hat{\rho}\\}$. For a pure state described by a
wavefunction $\bra{\psi}\hat{\mathcal{O}}\Ket{\psi}$. In a pure state, it is
equivalent to measurement in the wave function formalism. Let
$\\{\Ket{\psi_{i}}\\}$ be an orthonormal basis with
$\Ket{\psi_{k}}=\Ket{\psi}$ then

$tr\\{\hat{\mathcal{O}}\hat{\rho}\\}=\sum_{i}^{n}\bra{\psi_{i}}\hat{\mathcal{O}}\Ket{\psi}\bra{\psi}\Ket{\psi_{i}}=\bra{\psi}\hat{\mathcal{O}}\Ket{\psi}$
(6)

###  II.4  System ancilla interaction

The interaction between the system and the environment is represented as a
repeated interaction between the system and a subsystem of the environment.
Typically a thermal qubit. A general unitary interaction is employed.
Therefore the interaction can be expressed by its generator, the interaction
Hamiltonian:

$\hat{U}_{int}=e^{-i\hat{H}_{int}\theta}$ (7)

the phase angle $\theta$ has units of time. Since we can write an exponent as
a polynomial sum of $H_{int}$ we get

$[\hat{U}_{int},\hat{H}_{int}]=0$ (8)

###  II.5 Dynamics

In the collision model we assume that the uncorrelated thermal qubit is
employed only once. After interaction the swapped qubit is traced out. It is
therefore assumed that the state of the bath is unchanged. This assumption is
in accordance of an uncorrelated infinitely large bath and imposes that the
systems state is completely dependent on its previous state. The last remark
is a definition of Markovianity and is in accordance with postulate $4$ of
Kosloff and Dann [11]

$\Lambda_{t}=\Lambda_{t-s}\Lambda_{s}$ (9)

The dynamical map propagates density operators. A reduced map generated by a
global Hamiltonian Eq. (1), from an initial uncorrelated state defines a Kraus
map [14]. Such a map $\Lambda$ is a Completely Positive Trace Preserving
(CPTP) map on the system [11]. The generator of the dynamics is defined as:

$\mathcal{L}=\lim_{t\to 0}\frac{\Lambda(t)-\hat{\mathcal{I}}}{dt}$ (10)

Under the assumption that the collision period is much shorter than the
interval between collisions we can write the generator as follow:

$\mathcal{L}=-i[\hat{H},\bullet]+\gamma(tr_{b}\\{\hat{U}_{int}\bullet\otimes\hat{\rho}_{b}\hat{U}_{int}^{\dagger}\\}-\hat{\mathcal{I}}\bullet)~{}~{}~{},$
(11)

where $\gamma$ is the collision rate, eq. (11). The reduced description will
have the form

$\frac{d}{dt}\hat{\rho}_{S}=\mathcal{L}(\hat{\rho}_{S})=-i[\hat{H}_{S},\hat{\rho}_{S}]+\gamma(tr_{b}\\{\hat{U}_{int}\hat{\rho}_{S}\otimes\hat{\rho}_{b}\hat{S}_{p}^{\dagger}\\}-\hat{\rho}_{S})$
(12)

This structure has the Poissonian GKLS form [6]. Assuming the swap is
instantaneous we can write the integrated form of the dynamics as a sequence
of collision events, where $\gamma$ is determined by the average propagation
time dt A single collision event can be described by the super operator
$\mathcal{M}$ acting on the density operator $\rho_{s}$

$\mathcal{M}(\hat{H}_{S},\gamma,\theta,\beta)\hat{\rho}_{S}=(\hat{U}(\hat{H}_{S},dt))tr_{b}\\{\hat{U}_{int}(\theta)\hat{\rho}_{S}\otimes\hat{\rho}_{b}\hat{U}_{int}(\theta)^{\dagger})\\}(\hat{U}^{\dagger}(\hat{H}_{S},dt))$
(13)

And $k$ consecutive collisions

$\prod_{i=0}^{k}\mathcal{M}_{i}(\hat{H}_{S},\gamma,\theta,\beta)$ (14)

###  II.6  Unraveling of the density operator

The density operator $\hat{\rho}$ completely describes the state of the
quantum system. Any observable is determined by the relation
$\langle\mathcal{O}\rangle=tr\\{\hat{\rho}\mathcal{O}\\}$. The density
operator was introduced by von Neuman to describe statistical phenomena in
quantum mechanics [30]. It was observed that a pure state in an entangled
system described by a wave function is reduced to a mixed state when observing
the state of a subsystem.

The statistical character of the density operator is reflected by the
unraveling to an average of outer product of wave functions:

$\hat{\rho}=\sum_{k}p_{k}|\psi_{k}\rangle\langle\psi_{k}|$ (15)

where $\\{\Ket{\psi_{k}}\\}$ is the set of unraveling wavefunctions (not
necessarily orthogonal). The unraveling set $\\{\ket{\psi}\\}$ is not unique
which allows freedom which we will exploit. A straight forward unraveling is
obtained by diagonalizing the density operator.

###  II.7 Stochastic unravelling

The present study employs a stochastic unravelling scheme based on the flowing
lemma: let $\theta$ be a random phase wavefunction then

$\lim_{N\to\infty}\frac{1}{N}\sum_{j,k=1}^{N}e^{i(\theta_{j}-\theta_{k})}=\delta_{j,k}$
(16)

Let $\ket{\psi_{\theta}}$ be a wave function composed of an equal
superposition of an arbitrary orthonormal basis $\\{\Ket{n}\\}$ of size N with
random phase $\theta$

$\Ket{\psi_{\theta}}=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}e^{i\theta_{n}}\Ket{n}$
(17)

The identity operator $\hat{I}_{N}$ can be resolved by an infinite sum of
random wave functions.

$\lim_{K\to\infty}\frac{N}{K}\sum_{i=1}^{K}\Ket{\psi_{\theta_{i}}}\bra{\psi_{\theta_{i}}}=\hat{I}_{N}$
(18)

Since every density matrix can be diagonalized $\hat{\rho}$ can be decomposed
to:

$\hat{\rho}=\sum_{n=1}^{N}p_{i}\Ket{n}\bra{n}$ (19)

This allow to unravel the density operator with stochastic wave functions. We
will prove, that for every basis that an arbitrary $\hat{\rho}$ is diagonal in
we can create a set of stochastic wave functions that unravel $\hat{\rho}$.

$\Ket{\phi}=\sqrt{N}\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}\ket{\psi_{\theta}}$
(20)

$\Ket{\phi}$ is a normalized wave function:

$\braket{\phi}=\sqrt{N}\sum_{n=1}^{N}\sqrt{p_{n}}\bra{\psi_{\theta}}\ket{n}\bra{n}\sqrt{N}\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}\ket{\psi_{\theta}}=1$
(21)

The density operator $\hat{\rho}$ converges to an average over the outer
product of $\Ket{\phi}$ under the condition that it is a stochastic wave
function.

$\hat{\rho}=\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}=\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}~{}\hat{I}~{}\sum_{n=1}^{N}\sqrt{p_{n}}\ket{n}\bra{n}=$
(22)
$\sum_{n=0}^{N-1}\sqrt{p_{n}}\ket{n}\bra{n}\lim_{K\to\infty}\frac{N}{K}\sum_{i=1}^{K}\Ket{\psi_{\theta_{i}}}\bra{\psi_{\theta_{i}}}\sum_{n=0}^{N-1}\sqrt{p_{n}}\ket{n}\bra{n}=$
(23)
$\lim_{K\to\infty}\sum_{i=1}^{K}\frac{1}{K}\sqrt{N}\sum_{n=0}^{N-1}\sqrt{p_{n}}\ket{n}\bra{n}\Ket{\psi_{\theta_{i}}}\sqrt{N}\sum_{n=0}^{N-1}\sqrt{p_{n}}\bra{\psi_{\theta_{i}}}\ket{n}\bra{n}=$
(24)
$\lim_{K\to\infty}\frac{1}{K}\sum_{i=1}^{K}\Ket{\phi_{i}}\bra{\phi_{i}}\;\;\;\;\;\square$

A single wavefunction is a pure state and therefore cannot describe a
statistical distribution. In particular a thermal state is never pure. To
overcome this issue the stochastic unraveling method is employed to represent
the thermal ancilla qubits that collide with the system in accordance with Eq.
(20):

$\hat{\rho}=\sum_{j=0}^{N-1}p_{j}\Ket{\psi_{j}}\bra{\psi_{j}}$ (25)

$\Ket{\beta_{j}}=\sum_{i=0}^{1}\sqrt{N\frac{e^{-\beta\omega_{i}}}{Z}}\Ket{\omega_{i}}\bra{\omega_{i}}\Ket{\psi_{\theta_{j}}}=\sum_{i=0}^{1}\frac{\sqrt{N}e^{-\frac{\beta}{2}\omega_{i}}}{\sqrt{Z}}\Ket{\omega_{i}}\bra{\omega_{i}}\Ket{\psi_{\theta_{j}}}$
(26)
$\displaystyle\hat{\rho}_{b}=\lim_{K\to\infty}\frac{1}{K}\sum_{j=1}^{K}\ket{\beta_{j}}\bra{\beta_{j}}=\frac{e^{-\beta\hat{H}}}{Z}$
(27)

An ensamble average of interactions between the system and many
$\Ket{\beta_{j}}$ will converge to an interaction between the system and a
thermal qubit

$\lim_{K\rightarrow\infty}\frac{1}{K}\sum_{j=1}^{K}\hat{U}_{int}\Ket{\psi}\otimes\Ket{\beta_{j}}\bra{\psi}\otimes\bra{\beta_{j}}\hat{U}_{int}^{\dagger}=$
(28)
$\hat{U}_{int}\Ket{\psi}\bra{\psi}\otimes\Big{(}\lim_{K\rightarrow\infty}\frac{1}{K}\sum_{j=1}^{n}\Ket{\beta_{j}}\bra{\beta_{j}}\Big{)}\hat{U}_{int}^{\dagger}=\hat{U}_{int}\hat{\rho}_{s}\otimes\hat{\rho}_{b}\hat{U}_{int}^{\dagger}$

This method enables to describe a swap interaction with a thermal qubit in the
language of wave functions. In order to completely restore Eq. (14) we need to
translate the partial trace operation into wave function terminology.

###  II.8  Partial trace in a wave function description

Partial trace is an essential operation in obtaining the state of a subsystem
from a composite state. When the subsystem is entangled with its complementary
system the partial trace operation will lead to a mixed state.

To represent such a state with a wave function, stochastic unraveling will be
employed. The algorithm is designed to include in the wave functions, the
correct probabilities such that an average of their outer product will
reproduce the reduced subsystem state.

Even though the algorithm presented can be generalized for the tracing out of
any number of qubits, in this paper we present the algorithm of tracing out a
single qubit. Specifically we assume a system of $n$ spins and we trace out
the $k_{th}$ spin. $for\;\;b\in[0,2^{k-1}]\;\;,\;\;a\in[0,2^{n-k}]$ and
$i\in[0,1]$

We define $\Ket{\phi}^{i}$

$\Ket{\phi}^{i}_{b*2^{n-k}+a}=\frac{\Ket{\psi}_{(2b+i)*2^{n-k}+a}}{N_{i}}$
(29)

where N is the normalization factor and also the square root of the classical
probability of this state.
$\Ket{\phi}^{i}$ is not a random choice, it is embedding the physical meaning
of the measurement of the environment.
$\Ket{\phi}^{0}$ is a normalized vector of all the elements in $\Ket{\psi}$
condition on the state of traced out particle $\Ket{0}$.
$\Ket{\phi}^{1}$ is a normalized vector of all the elements in $\Ket{\psi}$
condition on the state of the traced out particle $\Ket{1}$.
The main result is

$tr_{k}\\{\Ket{\psi}\bra{\psi}\\}=\sum_{i=0}^{1}N_{i}^{2}\Ket{\phi}^{i}\bra{\phi}^{i}$
(30)

The proof is described in appendix (A.4).

The result above, beside complexity reduction also underlines the measurement
postulate in quantum mechanics and the equivalence of the partial trace with
partial measurement. The partial trace is a sum of the system’s possible
states after the traced out particle has ”collapsed” in to its possible states
with the adequate probability.

###  II.9  Stochastic Partial trace

The main result of section [II.8] is Eq. (30), giving an unravelling of the
mixed state. The description of the process requires us to operate on each
wave function of the unraveled tree separately as in figure 1. As will be
described in the next section, the probabilistic nature of the mixed state
enable to employ a Monte Carlo algorithm to randomly chose only one of the
functions. we have developed a generalized Monte Carlo algorithm for tracing
out more then one qubit which is not presented here.

$tr_{b_{x_{r}}}{\Ket{\psi}}=\begin{cases}\Ket{\phi}^{0}\quad\text{if}\quad
x_{r}<N_{0}^{2}\\\ \Ket{\phi}^{1}\quad\text{if}\quad x_{r}\geq
N_{0}^{2}\end{cases}$ (31)

The Monte Carlo algorithm induces the probability $N_{0}^{2}\quad is\quad
tr_{b_{x_{r}}}{\Ket{\psi}}=\Ket{\phi}^{0}$ and the complementary probability
$N_{1}^{2}\quad is\quad tr_{b_{x_{r}}}{\Ket{\psi}}=\Ket{\phi}^{1}$.

###  II.10  Describing the branching process.

The branching process is a sequence of free dynamics following a unitary
collision Cf. Sec. II.5. This process is described in the language of wave
functions, allowing an efficient algorithm by wave functions, that its
unravelling is converging to the density operator dynamics representation. By
combining the stochastic unravelling and partial trace presented in (20) and
(29) we can construct Eq. (14) that represent the consecutive collisions of
the density matrix with a thermal particle :

In details, every wave function $\Ket{\psi}$ will undergo three consecutive
operations:

  1. 1.

Interaction with a thermal wave function-
$\hat{U}_{int}(\Ket{\psi}\otimes\Ket{\beta})$.

  2. 2.

Stochastic partial trace -
${tr_{N+1_{x_{r}}}}\\{\hat{U}_{int}\Ket{\psi}\otimes\Ket{\beta}\\}$.

  3. 3.

Free dynamic of the system-
${\hat{U}}tr_{N+1_{x_{r}}}\\{\hat{U}_{int}\Ket{\psi}\otimes\Ket{\beta}\\}$.

To accurately restore Eq. (32), $\Ket{\psi}$ will accumulate $k$ thermal wave
functions $\Ket{\beta}$ $(k\rightarrow\infty)$. The partial trace will yield
two different wave functions with different probabilities, thus we will have
to average the outer product of all the outcomes with the correct weighs. This
procedure corresponds to a single collision.

$\lim_{K\to\infty}\frac{1}{K}\sum_{j=1}^{K}\int_{x_{r}=0}^{1}\hat{U}tr_{{N+1}_{x_{r}}}\\{\hat{U}_{int}\Ket{\psi}\otimes\Ket{\beta_{j}}\\}tr_{{N+1}_{x_{r}}}\\{\bra{\psi}\otimes\bra{\beta_{j}}\hat{U}_{int}^{\dagger}\\}\hat{U}^{\dagger}dx_{r}$

$\lim_{K\to\infty}\hat{U}tr_{N+1}\\{\hat{U}_{int}\Ket{\psi}\bra{\psi}\otimes\frac{1}{K}\sum_{j=1}^{K}\Ket{\beta_{j}}\bra{\beta_{j}}\hat{U}_{int}\\}\hat{U}^{\dagger}$

$\hat{U}(tr_{N+1}\\{\hat{U}_{int}(\hat{\rho}_{s}\otimes\hat{\rho}_{B})\hat{U}_{int}^{\dagger}\\})\hat{U}^{\dagger}$
(32)

For repeated collisions this process will recur for every outcome as shown in
figure (1). Eq. (14) is restored by consecutively employing Eq. (32) n times.
For a mathematical description we will define the super operator

$\hat{\mathcal{G}}(\psi,dt,\phi,x_{r},\beta)=\hat{U}tr_{{N+1}_{x_{r}}}\\{\hat{U}_{int}\Ket{\psi}\otimes\Ket{\beta_{j}}\\}tr_{{N+1}_{x_{r}}}$
(33)

$\lim_{K\to\infty}\frac{1}{K}\sum_{j=1}^{K}\int_{x_{r}=0}^{1}\hat{\mathcal{G}}^{n}(\psi,dt,\phi,x_{r},\beta)\hat{\mathcal{G}}^{\dagger
n}(\psi,dt,\phi,x_{r},\beta)dx_{r}=\hat{\rho}_{n}$ (34)

###  II.11  Stochastic convergence

This process is extremely computationally expensive, since it grows
exponentially with each collision as $O(2^{nk})$ where n is the number of
collisions and k is the number of thermal wave functions interacting with each
possible state of the system. This resolve in a branching tree illustrated in
Fig. 1.

A solution for this problem is to exploit the stochastic nature of the process
in three ways.

  1. 1.

Reduction to a binary tree due to a Wiener process. If we look at the path of
a single wave function (single branch in the tree in Fig. 1) we observe that
the systems wave function interacts in every collision with a stochastic
thermal wave function. The process satisfies the conditions of a Weiner
process - the process has a fixed initial condition and the stochastic part of
the bath particle in every collision has a mean 0 and a variance $\sigma^{2}$.
Therefore for sufficiently long process, where each collision
$U_{int}(\Ket{\psi_{k}}\bra{\psi_{k}}\otimes\rho_{b})U_{int}^{\dagger}$ is
represented as the average of the outer product of K
$U_{int}(\Ket{\psi_{k}}\otimes\Ket{\beta_{j}})$ and thus splited into K
branches, can be represented by only one $\Ket{\beta_{j}}$. For many collision
it will, follow the Weiner process, undergo a Brownian motion and the process
will converge to a consecutive interaction with a thermal qubit. As a result
the unravelling of the tree in figure 2 will converge to the unravelling of
the larger tree in figure 1.

As can be seen in figure 2 this method will result in a properly weighed
sample of the binomial distribution around the most probable state.

  2. 2.

The nature of tracing out a single bath qubit results in a mixed state
composed of two pure states eq. (30). Using the property of Weiner process,
computation of all possibilities with the correct weights will be resolved in
a binary tree with changing probability weights. We used a Monte Carlo
stochastic partial trace algorithm sec. II.8 in order to stochastically chose
one of the two wave functions constructing the mixed state imposed by the
partial trace in every step.

  3. 3.

Based on the central Limit Theorem, the average of a sequence of independent
and identically distributed random variables drawn from a distribution of
expected value given by $\mu$ and finite variance given by $\sigma^{2}$ will
converge in probability to a normal distribution. The multidimensional Central
Limit Theorem generalize the theory and state that a random vector (satisfying
the vector space axioms) will converge in probability to a multi-variant
Gaussian. Mathematically

$\sqrt{n}(\hat{X_{n}}-\mu)\rightarrow^{d}\mathcal{N}(0,\Sigma)$ (35)

where $\Sigma$ is the covariance matrix.

Employing Eq. (34) the average of the outer product of all wave function
possibilities with the adequate probabilities represented by the tree in Fig.
1 converges to the density matrix $\rho_{n}$ satisfying the collision model in
sec. II.5. Thus a sample of the outer product of $K$ identical wave functions
undergoing n consecutive collisions by eq. (34):
$\mathcal{G}^{n}(\psi,dt,\phi,x_{r},\beta){\mathcal{G}^{n}}^{\dagger}(\psi,dt,\phi,x_{r},\beta)$
is a sequence of independent and random variables drawn from a distribution of
expected value $\rho_{n}$ and a finite variance and thus

$\lim_{K\to\infty}\frac{1}{K}\sum_{i=1}^{K}\mathcal{G}^{n}(\psi,dt,\phi,x_{r},\beta){\mathcal{G}^{n}}^{\dagger}(\psi,dt,\phi,x_{r},\beta)=\hat{\rho}_{n}$
(36)

Moreover with n- number of collisions increasing we expect convergence in
probability

$\frac{1}{K}\sum_{i=1}^{K}\mathcal{G}^{n}(\psi,dt,\phi,x_{r},\beta){\mathcal{G}^{n}}^{\dagger}(\psi,dt,\phi,x_{r},\beta)-\hat{\rho}_{n}\rightarrow^{d}\mathcal{N}(0,\frac{\Sigma}{K})$
(37)

Central limit theorem redundant mathematical use of Brownian motion. Yet, we
have found it useful for physical intuition of the process.

##  III Results

###  III.1 Convergence

To demonstrate the approach in accordance with sec. II.11, we will study a
specific example of a unitary interaction. A partial swap between the last
particle in a system of qubits and an uncorrelated thermal qubit as
interaction is specifically chosen. This type of interaction has been
addressed in the collision model review by Ciccarello, Lorenzo, Giovannetti
and Palma [31].

In a two qubit system the swap algorithm becomes the swap gate:

$\hat{S}=\Ket{0}\bra{0}\otimes\Ket{0}\bra{0}+\Ket{1}\bra{1}\otimes\Ket{1}\bra{1}+\Ket{1}\bra{0}\otimes\Ket{0}\bra{1}+\Ket{0}\bra{1}\otimes\Ket{1}\bra{0}=\begin{pmatrix}1&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\end{pmatrix}$ (38)

In a bigger system composed of $N$ qubits, an operation that swaps between the
$i~{}th$ and $j~{}th$ qubits the swap algorithm looks similar:

let

$\prod_{k=0}^{i-1}\otimes
I_{2^{k}}\otimes\Ket{\delta}\otimes\prod_{k=i+1}^{N}\otimes
I_{2^{k}}=\Ket{\delta_{i}}$ (39) $\prod_{k=0}^{i-1}\otimes
I_{2^{k}}\otimes\bra{\delta}\otimes\prod_{k=i+1}^{N}\otimes
I_{2^{k}}=\bra{\delta_{i}}$ (40)

A swap operation between particles i and j can be described as

$\hat{S}_{i,j}=\Ket{0_{i}}\bra{0_{i}}\Ket{0_{j}}\bra{0_{j}}+\Ket{1_{i}}\bra{1_{i}}\Ket{1_{j}}\bra{1_{j}}+\Ket{0_{i}}\bra{1_{i}}\Ket{1_{j}}\bra{0_{j}}+\Ket{1_{i}}\bra{0_{i}}\Ket{0_{j}}\bra{1_{j}}$
(41)

Swap is a unitary operation - see appendix (A.1)

$\hat{S}\hat{S}^{\dagger}=\hat{S}_{k,m}\hat{S}_{k,m}^{\dagger}=\hat{I}_{2^{n}}$
(42)

$\hat{S}$ is unitary and real, therefor $\hat{S}=\hat{S}^{\dagger}$.

The unitarity of the swap operator makes it a preferable candidate to simulate
interaction. It has the advantage that it can be realized on a quantum
computer.

The partial swap is defined as:

$\hat{S}_{p}=\cos(\theta)\hat{I}_{2^{n}}+i\sin(\theta)\hat{S}$ (43)

The partial swap, like the full swap, is also a unitary operation
$\hat{S}_{p}\hat{S}_{p}^{\dagger}=\hat{I}$\- see appendix (A.2). Unlike the
full swap, the partial swap induces a spectrum of interaction strengths and
maintains correlation between all particles in the system after reduced
description of the setup as we will see in Sec. (II.4). For this reason we
chose it as a preferable candidate for the interaction between the ancilla and
the system.
Since $S_{p}$ is unitary. by definition we can write

$\hat{S}_{p}=e^{-i\hat{H}_{int}\theta}$ (44)

$\theta$ can be either a phase, angle of interaction or time of interaction.
Since we can write an exponent as a polynomial sum of $\hat{H}_{int}$ we get

$[\hat{S}_{p},\hat{H}_{int}]=0$ (45)

if we call the unravelling of K processes

$\frac{1}{K}\sum^{K}\mathcal{G}^{n}{\mathcal{G}^{n}}^{\dagger}=\Theta_{n}(K)$
(46)

We expect the distance

$|\Theta_{n}(K)-\rho_{n}|$ (47)

To converge in probability to $\mathcal{N}(0,\frac{\Sigma}{K})$\- a normal
distribution around zero. In order to examine the convergence rate of our
model, we chose to observe the variance of the distribution of this distance
and expect it to converge as $\frac{1}{K}$. The distance function we chose was
the variance between each element of the unravelling and the density matrix
$\rho_{s}$ of the system under going the same dynamic in the density matrix
form- Eq. (14). $N^{2}$ being the dimension of the system, K- the number of
realizations and n- the number of collisions.

$D(\mathcal{G},\rho,n,N,K)=\frac{1}{N^{2}}\sum_{i,j=1}^{N^{2}}|\rho_{n_{i,j}}-\Theta(K)_{n_{i,j}}|^{2}$
(48)

Figure 3:  The distance function defined by a normalized sum of the absolute
value squared of the euclidean distance between same elements of $\rho_{n}$
and $\Theta_{n}$. That is taken as a normalization of
$\mathcal{D}(\mathcal{G},\rho,n,N)=\frac{1}{N^{2}}\sum_{i=1}^{N^{2}}|\rho_{n_{i,j}}-\Theta_{n_{i,j}}|^{2}$.
We expect that $\cal{D}$ will converge as $\frac{1}{K}$. Since we expect the
variance to converge as $\frac{1}{K}$, we expect the variance of each element
also to converge as $\frac{1}{K}$ and therefore that $\mathcal{D}$, that is a
sum of variance element of each matrix element will converge as $\frac{1}{K}$
as well. The graph exhibits the ${\cal D}$ function for 3 systems undergoing
600 collisions. We can see a decay to zero as $\frac{1}{K}$ for a system
composed of 5,8 and 10 particles as expected

According to the multidimensional Central Limit Theorem we expect a
convergence of the covariance matrix as $\frac{1}{K}$, and thus the variance
of every element $|\Theta_{i,j}-\rho_{i,j}|$ to $\sigma_{i,j}^{2}$. Since each
element of $\Theta$ has the same dependency on the stochastic variable we
expect $\sigma_{i,j}=\sigma_{k,l}$. Thus we expect D to approximate the
average of all element variances.

$D(\mathcal{G},\rho,n,N,K)=\frac{1}{N^{2}}\sum_{i,j=1}^{N^{2}}|\rho_{n_{i,j}}-\Theta_{n_{i,j}}(K)|^{2}\approx\sum_{i,j}^{N^{2}}\frac{1}{K}\sum_{k}^{K}|\rho_{n_{i,j}}-\Theta_{n_{i,j_{k}}}|^{2}$
(49)

From the central limit theorem we expect that the sum of all element variances
will decay as $\frac{1}{K}$ with K the number of realization

$\sum_{i,j}^{N^{2}}\frac{1}{K}\sum_{k}^{K}|\rho_{n_{i,j_{k}}}-\Theta_{n_{i,j_{k}}}|^{2}\sim\frac{1}{K}$
(50)

thus

$D(\mathcal{G},\rho,n,N,K)\sim\frac{1}{K}$ (51)

As we can seen in figure 3 we observe that ${\cal D}$ indeed converges as
expected.

##  IV  Conclusions

In this paper we developed the basic tools for unraveling open system dynamics
with wavefunctions based on a collision model. The stochastic averaging of
these wavefunctions converges for all expectation values, equivalent to the
density operator formalism. The basic algorithm can be divided into three
steps:

  1. 1.

Implementing a Monte Carlo stochastic partial trace algorithm. Sec. II.8,II.9
.

  2. 2.

Restoring a mixed state by unraveling of stochastic wave functions. This was
achieved for the bath particle colliding with the system and for the system
itself after averaging over many stochastic system wave functions undergoing
the dynamics. Sec. (II.7)

  3. 3.

Using statistical properties such as the central limit theorem we observed
convergence to the density matrix described dynamics. The convergence was
achieved with a relatively small number of realizations. This property results
in high computation efficiency.

An illustration of the algorithm is carried out in Sec. (III.1).

The algorithm developed contains a nonlinear component. In the implementation
of the Monte Carlo algorithm in Sec. (II.9), one of the wave functions
composing the mixed state is selected by partial trace. The probability to
choose each of the states, is calculated from the normalization factor of one
of the states. This requires to first calculate one possible outcome. This
outcome might not be used in the next step. If the selected wavefunction is
the one that was not calculated, the other wave function has to be
recalculated. Thus, the non-linearity results from renormalizing the
wavefunction and the possibility of computing an additional wavefunctions if
it was absent in the Monte Carlo lottery.

The modeling method addresses a major problem of the cost of simulating open
quantum system. The wavefunction representation reduces the memory requirement
and as the system becomes larger the speed of convergence to the full
simulation also increases. We therefore expect a reduction in computational
cost of up to a factor of $N$ where $N$ is the size of Hilbert space. As a
result the boundaries of possible simulations are stretched.

In addition the wavefunction method allows additional insight in the physical
process taking place from the viewpoint of a single event.

Finally the simulation can be implemented on quantum computers- Due to the
fact that all operations are unitary and the description of the setup is by
wavefunctions.

###### Acknowledgements.

We thank Christiane Koch, Gil Katz, and Florian Habecker for sharing their
insight. Work supported by the Israel Science Foundation (Grants No. 510/17
and 526/21).

##  Appendix A Appendixes

###  A.1 Proof of Unitarity of $\hat{S}$

$\hat{S}_{i,j}\hat{S}_{i,j}^{\dagger}=(\Ket{0_{i}}\Ket{0_{j}}\bra{0_{j}}\bra{0_{i}}+\Ket{0_{i}}\Ket{1_{j}}\bra{1_{j}}\bra{0_{i}}+\Ket{1_{i}}\Ket{0_{j}}\bra{0_{j}}\bra{1_{i}}+\Ket{1_{i}}\Ket{1_{j}}\bra{1_{j}}\bra{1_{i}}($
(52)
$(\Ket{0_{j}}\Ket{0_{i}}\bra{0_{i}}\bra{0_{j}}+\Ket{0_{j}}\Ket{1_{i}}\bra{1_{i}}\bra{0_{j}}+\Ket{1_{j}}\Ket{0_{i}}\bra{0_{i}}\bra{1_{j}}+\Ket{1_{j}}\Ket{1_{i}}\bra{1_{i}}\bra{1_{j}})=$
$\Ket{0_{i}}\Ket{0_{j}}\bra{0_{j}}\bra{0_{i}}\Ket{0_{j}}\Ket{0_{i}}\bra{0_{i}}\bra{0_{j}}+0+0+0+$
$0+0+\Ket{0_{i}}\Ket{1_{j}}\Ket{1_{j}}\Ket{0_{i}}\bra{0_{i}}\bra{1_{j}}+0+$
$0+\Ket{1_{i}}\Ket{0_{j}}\bra{0_{j}}\bra{1_{i}}\Ket{0_{j}}\Ket{1_{i}}\bra{1_{i}}\bra{0_{j}}+0+0+$
$0+\Ket{1_{i}}\Ket{1_{j}}\bra{1_{j}}\bra{1_{i}}\Ket{1_{j}}\Ket{1_{i}}\bra{1_{i}}\bra{1_{j}})=$
$\Ket{0_{i}}\Ket{0_{j}}\bra{0_{i}}\bra{0_{j}}+\Ket{0_{i}}\Ket{1_{j}}\bra{0_{i}}\bra{1_{j}}+\Ket{1_{i}}\Ket{0_{j}}\bra{1_{i}}\bra{0_{j}}+\Ket{1_{i}}\Ket{1_{j}}\bra{1_{i}}\bra{1_{j}}=$
$(\Ket{0_{i}}\bra{0_{i}}+\Ket{1_{i}}\bra{1_{i}})(\Ket{0_{j}}\bra{0_{j}}+\Ket{1_{j}}\bra{1_{j}})=\hat{I}_{2^{N}}$

###  A.2 Proof that the partial swap is unitary.

$\hat{S}_{p}=\cos(\theta)\hat{I}_{2^{n}}+i\sin(\theta)\hat{S}$ (53)
$\hat{S}_{p}\hat{S}_{p}^{\dagger}=(cos(\theta)\hat{I}_{2^{n}}+i\sin(\theta)\hat{S})(\cos(\theta)\hat{\hat{I}}_{2^{n}}+i\sin(\theta)\hat{S})^{\dagger}=\\\
\cos^{2}(\theta)\hat{I}_{2^{n}}*\hat{I}_{2^{n}}+\sin^{2}(\theta)\hat{S}\hat{S}^{\dagger}+i\sin(\theta)\cos(\theta)\hat{I}_{2^{n}}\hat{S}^{\dagger}-i\sin(\theta)\cos(\theta)\hat{I}_{2^{n}}=\\\
\hat{I}_{2^{n}}(\cos^{2}(\theta)+\sin^{2}(\theta))+\hat{S}(i\sin(\theta)\cos(\theta)-i\sin(\theta)\cos(\theta))=\hat{I}_{2^{n}}+0=\hat{I}_{2^{n}}$
(54)

###  A.3  Tracing out the n’th qubit tensor thermal qubit

$tr_{n+1}\\{\hat{S}_{n,n+1}(\Ket{\psi}\bra{\psi}\otimes\ket{\beta}\bra{\beta})\hat{S}_{n,n+1}^{\dagger}\\}=tr_{n}\\{\ket{\psi}\bra{\psi}\\}\otimes\ket{\beta}\bra{\beta}$
(55)

$\hat{I}_{2^{n-1}}\otimes\Ket{0}=\Ket{\mathcal{O}}\;\;\;\hat{I}_{2^{n-1}}\otimes\Ket{1}=\Ket{\mathcal{I}}$
(56)

$\hat{S}_{n,n+1}=\hat{I}_{2^{n-1}}\otimes(\Big{(}\sum_{i=0}^{1}\Ket{i}\bra{i}\otimes\Ket{i}\bra{i}+\sum_{k\neq
j=0}^{1}\Ket{k}\bra{j}\otimes\Ket{j}\bra{k}\Big{)}$ (57)
$tr_{n+1}\\{\hat{S}_{n,n+1}(\Ket{\psi}\bra{\psi}\otimes\ket{\beta}\bra{\beta})\hat{S}_{n,n+1}^{\dagger}\\}=$
(58)
$\sum_{l=0}^{1}I_{2^{n}}\otimes\bra{l}\Big{(}I_{2^{n-1}}\otimes\Big{(}\sum_{i=0}^{1}\Ket{i}\bra{i}\otimes\Ket{i}\bra{i}+\sum_{k\neq
j=0}^{1}\Ket{k}\bra{j}\otimes\Ket{j}\bra{k}\Big{)}\ket{\psi}\bra{\psi}\otimes\ket{\beta}\bra{\beta}\\\
\Big{(}I_{2^{n-1}}\otimes\Big{(}\sum_{i=0}^{1}\Ket{i}\bra{i}\otimes\Ket{i}\bra{i}+\sum_{k\neq
j=0}^{1}\Ket{k}\bra{j}\otimes\Ket{j}\bra{k}\Big{)}=\\\ $
$I_{2^{n}}\otimes\bra{0}\Big{(}\Ket{\mathcal{O}}\bra{\mathcal{O}}\otimes\Ket{0}\bra{0}+\Ket{\mathcal{I}}\bra{\mathcal{O}}\otimes\Ket{0}\bra{1}\Big{)}\ket{\psi}\bra{\psi}\otimes\\\
\ket{\beta}\bra{\beta}\Big{(}\Ket{\mathcal{O}}\bra{\mathcal{O}}\otimes\Ket{0}\bra{0}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\otimes\Ket{1}\bra{0}\Big{)}I_{2^{n}}\otimes\ket{0}\\\
+I_{2^{n}}\otimes\bra{1}\Big{(}\Ket{\mathcal{I}}\bra{\mathcal{I}}\otimes\Ket{I}\bra{I}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\otimes\Ket{1}\bra{0}\Big{)}\ket{\psi}\bra{\psi}\otimes\\\
\ket{\beta}\bra{\beta}\Big{(}\Ket{\mathcal{I}}\bra{\mathcal{I}}\otimes\Ket{1}\bra{1}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\otimes\Ket{1}\bra{0}\Big{)}I_{2^{n}}\otimes\ket{1}=$
$\Ket{\mathcal{O}}\bra{\mathcal{O}}\ket{\psi}\braket{0}{\beta}\braket{\beta}{0}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{O}}+\Ket{\mathcal{O}}\bra{\mathcal{O}}\ket{\psi}\braket{0}{\beta}\braket{\beta}{1}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{I}}+\\\
\Ket{\mathcal{I}}\bra{\mathcal{O}}\ket{\psi}\braket{1}{\beta}\braket{\beta}{0}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{O}}+\Ket{\mathcal{I}}\bra{\mathcal{O}}\ket{\psi}\braket{1}{\beta}\braket{\beta}{1}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{I}}+\\\
\Ket{\mathcal{I}}\bra{\mathcal{I}}\ket{\psi}\braket{1}{\beta}\braket{\beta}{1}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{I}}+\Ket{\mathcal{I}}\bra{\mathcal{I}}\ket{\psi}\braket{1}{\beta}\braket{\beta}{0}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{O}}+\\\
\Ket{\mathcal{O}}\bra{\mathcal{I}}\ket{\psi}\braket{0}{\beta}\braket{\beta}{1}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{I}}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\ket{\psi}\braket{0}{\beta}\braket{\beta}{1}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{O}}$

$\braket{0}{\beta}\braket{\beta}{0}\Big{(}\Ket{\mathcal{O}}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{O}}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{O}}\Big{)}+\\\
\braket{1}{\beta}\braket{\beta}{1}\Big{(}\Ket{\mathcal{I}}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{I}}+\Ket{\mathcal{I}}\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{I}}\Big{)}+\\\
\braket{1}{\beta}\braket{\beta}{0}\Big{(}\Ket{\mathcal{I}}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{O}}+\Ket{\mathcal{I}}\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{O}}\Big{)}+\\\
\braket{0}{\beta}\braket{\beta}{1}\Big{(}\Ket{\mathcal{O}}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}\bra{\mathcal{I}}+\Ket{\mathcal{O}}\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\bra{\mathcal{I}}\Big{)}$
$\braket{0}{\beta}\braket{\beta}{0}\Ket{\mathcal{O}}\Big{(}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}+\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\Big{)}\bra{\mathcal{O}}+\\\
\braket{1}{\beta}\braket{\beta}{1}\Ket{\mathcal{I}}\Big{(}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}+\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\Big{)}\bra{\mathcal{I}}+\\\
\braket{0}{\beta}\braket{\beta}{1}\Ket{\mathcal{O}}\Big{(}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}+\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\Big{)}\bra{\mathcal{I}}+\\\
\braket{1}{\beta}\braket{\beta}{0}\Ket{\mathcal{I}}\Big{(}\bra{\mathcal{O}}\ket{\psi}\bra{\psi}\Ket{\mathcal{O}}+\bra{\mathcal{I}}\ket{\psi}\bra{\psi}\Ket{\mathcal{I}}\Big{)}\bra{\mathcal{I}}=$

$\braket{0}{\beta}\braket{\beta}{0}I_{2^{n-1}}\otimes\Ket{0}tr_{n}\\{\Ket{\psi}\bra{\psi}\\}I_{2^{n-1}}\otimes\bra{0}+\\\
\braket{1}{\beta}\braket{\beta}{1}I_{2^{n-1}}\otimes\Ket{1}tr_{n}\\{\Ket{\psi}\bra{\psi}\\}I_{2^{n-1}}\otimes\bra{1}+\\\
\braket{0}{\beta}\braket{\beta}{1}I_{2^{n-1}}\otimes\Ket{0}tr_{n}\\{\Ket{\psi}\bra{\psi}\\}I_{2^{n-1}}\otimes\bra{1}+\\\
\braket{1}{\beta}\braket{\beta}{0}I_{2^{n-1}}\otimes\Ket{1}tr_{n}\\{\Ket{\psi}\bra{\psi}\\}I_{2^{n-1}}\otimes\bra{0}=$
$tr_{n}\\{\Ket{\psi}\bra{\psi}\\}\Big{(}\braket{0}{\beta}\braket{\beta}{0}\Ket{0}\bra{0}+\braket{1}{\beta}\braket{\beta}{1}\Ket{1}\bra{1}+\\\
\braket{0}{\beta}\braket{\beta}{1}\Ket{0}\bra{1}+\braket{1}{\beta}\braket{\beta}{0}\Ket{1}\bra{0}\Big{)}=\\\
tr_{n}\\{\ket{\psi}\bra{\psi}\\}\otimes\ket{\beta}\bra{\beta}\;\;\;\;\;\square$

###  A.4  Partial trace is sum of system when bath is measured in different
states

$tr_{k}\\{\Ket{\psi}\bra{\psi}\\}=\sum_{i=0}^{1}I_{2^{k-1}}\otimes\bra{i}\otimes
I_{2^{N-k}}\Ket{\psi}\bra{\psi}I_{2^{k-1}}\otimes\Ket{i}\otimes I_{2^{N-k}}$
(59) $I_{2^{k-1}}\otimes\Ket{i}\otimes I_{2^{N-k}}=P^{i}$ (60)
$I_{2^{k-1}}\otimes\bra{i}\otimes I_{2^{N-k}}=(P^{i})^{\dagger}$ (61)
$tr_{k}\\{\Ket{\psi}\bra{\psi}\\}=\sum_{i=0}^{1}{P^{i}}^{\dagger}\Ket{\psi}\bra{\psi}P^{i}$
(62)

$for\;\;b\in[0,2^{k-1}-1]\;\;and\;\;a\in[0,2^{N-k}-1]$

$\hat{P}^{i}_{[l,m]}=\begin{cases}1\;\;if\;\;l=(2b+i)*2^{N-k}+a,m=b*2^{N-k}+a\\\
0\;\;else\end{cases}$ (63)
$\hat{P^{i}}^{\dagger}_{[l,m]}=\begin{cases}1\;\;if\;\;l=b*2^{N-k}+a,m=(2b+i)*2^{N-k}+a\\\
0\;\;else\end{cases}$ (64)

$l=r*2^{N-k}+x\;\;\;\;\;\;j=y*2^{N-k}+z\;\;\;for\;\;y,r\in[0,k+1]\;\;and\;\;x,z\in[0,2^{N-k}]$

$tr_{k}\\{\Ket{\psi}\bra{\psi}\\}_{l,j}=\sum_{i,d,e}\hat{P^{i}}^{\dagger}_{l,d}(\Ket{\psi}\bra{\psi})_{d,e}\hat{P}^{i}_{e,j}=$
(65)
$\sum_{i=0}^{1}\hat{P}^{i^{\dagger}}_{[r*2^{N-k}+x,(2r+i)*2^{N-k}+x]}(\Ket{\psi}\bra{\psi})_{[(2r+i)*2^{N-k}+x,(2y+i)*2^{N-k}+z]}\hat{P}^{i}_{[2(y+i)*2^{N-k}+z,y*2^{N-k}+z]}=$
(66)
$(\sum_{i=0}^{1}(\Ket{\psi}\bra{\psi})_{[(2r+i)*2^{N-k}+x,(2y+i)*2^{N-k}+z]})=\sum_{i=0}^{1}N_{i}\Ket{\phi}^{i}_{r*2^{N-k}+x}N_{i}\bra{\phi}^{i}_{y*2^{N-k}+z}=$
(67)
$\sum_{i=0}^{1}N_{i}^{2}\Ket{\phi}^{i}\bra{\phi}^{i}_{[r*2^{N-k}+x,y*2^{N-k}+z]}=\sum_{i=0}^{1}N_{i}^{2}\Ket{\phi}^{l}\bra{\phi}^{i}_{l,j}\;\;\;\;\;\square$
(68)

## References

  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Doherty  _et al._ [2012] M. Doherty, F. Dolde, H. Fedder, F. Jelezko, J. Wrachtrup, N. Manson, and L. Hollenberg, Physical Review B 85, 205203 (2012).
  * Doherty  _et al._ [2013] M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. Hollenberg, Physics Reports 528, 1 (2013).
  * Wangsness and Bloch [1953] R. K. Wangsness and F. Bloch, Physical Review 89, 728 (1953).
  * Redfield [1957] A. G. Redfield, IBM Journal of Research and Development 1, 19 (1957).
  * Lindblad [1976] G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Davies [1974] E. B. Davies, Communications in mathematical Physics 39, 91 (1974).
  * Alicki and Kosloff [2018] R. Alicki and R. Kosloff, in  _Thermodynamics in the Quantum Regime_ (Springer, 2018) pp. 1–33.
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Dann and Kosloff [2021] R. Dann and R. Kosloff, Physical Review Research 3, 023006 (2021).
  * Berman and Kosloff [1991] M. Berman and R. Kosloff, Computer physics communications 63, 1 (1991).
  * Prokof’ev and Stamp [2000] N. Prokof’ev and P. Stamp, Reports on Progress in Physics 63, 669 (2000).
  * Kraus [1974] K. Kraus, in  _Foundations of quantum mechanics and ordered linear spaces_ (Springer, 1974) pp. 206–229.
  * Gelman and Kosloff [2003] D. Gelman and R. Kosloff, Chemical physics letters 381, 129 (2003).
  * Ezra  _et al._ [2021] B. Ezra, R. Kosloff, and S. Kallush, arXiv preprint arXiv:2110.12458 (2021).
  * Baer  _et al._ [2013] R. Baer, D. Neuhauser, and E. Rabani, Physical review letters 111, 106402 (2013).
  * Gisin and Percival [1992] N. Gisin and I. C. Percival, Journal of Physics A: Mathematical and General 25, 5677 (1992).
  * Katz  _et al._ [2008] G. Katz, D. Gelman, M. A. Ratner, and R. Kosloff, The Journal of chemical physics 129, 034108 (2008).
  * Torrontegui and Kosloff [2016] E. Torrontegui and R. Kosloff, New Journal of Physics 18, 093001 (2016).
  * Baer and Kosloff [1997] R. Baer and R. Kosloff, The Journal of chemical physics 106, 8862 (1997).
  * Koch  _et al._ [2003] C. P. Koch, T. Klüner, H.-J. Freund, and R. Kosloff, The Journal of chemical physics 119, 1750 (2003).
  * Karplus and Schwinger [1948] R. Karplus and J. Schwinger, Phys. Rev. 73, 1020 (1948).
  * Rau [1963] J. Rau, Phys. Rev. 129, 1880 (1963).
  * Caves [1986] C. M. Caves, Phys. Rev. D 33, 1643 (1986).
  * Caves and Milburn [1987] C. M. Caves and G. J. Milburn, Phys. Rev. A 36, 5543 (1987).
  * Kosloff [2019] R. Kosloff, The Journal of chemical physics 150, 204105 (2019).
  * Li  _et al._ [2022] L. Li, Z.-X. Man, and Y.-J. Xia, Entropy 24, 644 (2022).
  * O’Connor  _et al._ [2021] E. O’Connor, B. Vacchini, and S. Campbell, Entropy 23, 1634 (2021).
  * Von Neumann [2018] J. Von Neumann,  _Mathematical foundations of quantum mechanics_ (Princeton university press, 2018).
  * Ciccarello  _et al._ [2022] F. Ciccarello, S. Lorenzo, V. Giovannetti, and G. M. Palma, Physics Reports 954, 1 (2022).