# Bounds and Algorithms for Frameproof Codes and Related Combinatorial
Structures

Marco Dalai1, Stefano Della Fiore12, Adele A. Rescigno2 and Ugo Vaccaro2
1Department of Information Engineering, University of Brescia, Brescia (BS),
Italy
Emails: {marco.dalai<EMAIL_ADDRESS>2Department of Computer
Science, University of Salerno, Fisciano (SA), Italy
Emails: {sdellafiore, arescigno<EMAIL_ADDRESS>
###### Abstract

In this paper, we study upper bounds on the minimum length of frameproof codes
introduced by Boneh and Shaw [3] to protect copyrighted materials. A $q$-ary
$(k,n)$-frameproof code of length $t$ is a $t\times n$ matrix having entries
in $\\{0,1,\ldots,q-1\\}$ and with the property that for any column
$\mathbf{c}$ and any other $k$ columns, there exists a row where the symbols
of the $k$ columns are _all_ different from the corresponding symbol (in the
same row) of the column $\mathbf{c}$. In this paper, we show the existence of
$q$-ary $(k,n)$-frameproof codes of length $t=O(\frac{k^{2}}{q}\log n)$ for
$q\leq k$, using the Lovász Local Lemma, and of length
$t=O(\frac{k}{\log(q/k)}\log(n/k))$ for $q>k$ using the expurgation method.
Remarkably, for the practical case of $q\leq k$ our findings give codes whose
length almost matches the lower bound $\Omega(\frac{k^{2}}{q\log k}\log n)$ on
the length of _any_ $q$-ary $(k,n)$-frameproof code and, more importantly,
allow us to derive an algorithm of complexity $O(tn^{2})$ for the construction
of such codes.

###### Index Terms:

Frameproof codes, strongly-selective codes, superimposed codes, Lovász Local
Lemma.

##  I Introduction

Frameproof codes were proposed by Boneh and Shaw [3] as a tool to protect
copyrighted materials from unauthorized use. The rationale is the following.
When a distributor wants to sell $n$ copies of a digital product, she first
chooses $t$ fixed locations in the digital product. Successively, for each
product copy, she associates to each chosen location a $q$-ary symbol. Such a
collection of chosen symbols and locations in each product copy is known as a
_fingerprint_ , which can be seen as a codeword of length $t$ over an alphabet
of size $q$. The clients do not know the locations and symbols stored in the
data, which means they cannot remove or modify them. However, a group of
malicious clients could collude, share and compare their copies. In this
scenario, the malicious clients could in principle recover the locations and
symbols of the fingerprints in order to create an illegal copy, that is, a
copy of the digital product whose fingerprint equals that of an innocent
client outside the malicious coalition. To prevent this situation, frameproof
codes can be used. Informally, a set of fingerprints (codewords) is called a
$(k,n)$-frameproof code if any coalition of at most $k$ clients cannot create
any other codeword outside their coalition. Essentially, this means that any
codeword possesses a piece of information (in some coordinate) not included in
the same coordinate of any collection of at most $k$ other codewords (a more
formal definition is provided later on).

The main problems in the area consist of both providing bounds on the lengths
of frameproof codes (this is motivated by the practical consideration that the
length of the fingerprint to insert in the digital data represents an obvious
overhead that needs to be minimized) _and_ computationally efficient
procedures to construct frameproof codes of length close to the theoretical
optimum. The seminal paper [3] spurred an interesting line of research to
investigate the two problems mentioned above. Due to space constraints, it is
impossible here to summarize the many results in the area, and we refer the
reader to the papers [1, 2, 4, 5, 10, 15, 16, 17, 18, 19, 20] and references
therein quoted for the relevant literature in the area.

###  I-A Our Results

In this paper we show the existence of $q$-ary $(k,n)$-frameproof codes of
length $t=O(\frac{k^{2}}{q}\log n)$ for $q\leq k$, using the Lovász Local
Lemma, and of length $t=O(\frac{k}{\log(q/k)}\log(n/k))$ for $q>k$, using the
expurgation method. We refer the reader to [8] and [11] for similar
applications of these methods which are used to provide probabilistic
constructions of separating codes and identifiable parent property (IPP)
codes. We will discuss how our results improve on known results present in the
literature after having proved them. Moreover, for the practically important
case of $q\leq k$ (as motivated in [15]) we provide a $O(tn^{2})$ randomized
algorithm to construct codes of length $t=O(\frac{k^{2}}{q}\log n)$, that
almost matches the lower bound $\Omega(\frac{k^{2}}{q\log k}\log n)$ on the
length of _any_ $q$-ary $(k,n)$-frameproof code. To the best of our knowledge,
this is the first polynomial time, in the code size $n$, algorithm with such a
performance. We note that to explicitly construct $(k,n)$-frameproof codes of
length $t$ any algorithm requires $\Omega(tn)$ time.

##  II Preliminaries

Throughout the paper, the logarithms without subscripts are in base two, and
we denote with $\ln(\cdot)$ the natural logarithm. Given integers $a<b$, we
denote with $[a,b]$ the set ${\\{a,a+1,\ldots,b\\}}$. We start by introducing
the combinatorial objects we study in this paper.

######  Definition II.1.

Let $k$, $n$, $q\geq 2$ be positive integers, $n>k$. A $q$-ary
$(k,n)$-_frameproof_ code is a $t\times n$ matrix $M$ with entries in
$[0,q-1]$ such that for any column $\mathbf{c}$ and any other $k$ columns of
$M$ we have that, there exists a row ${i\in[1,t]}$ where the symbols of the
$k$ columns are all different from the symbol in the $i$-th row of the column
$\mathbf{c}$. The number of rows $t$ of $M$ is called the length of the
$q$-ary $(k,n)$-frameproof code.

In order to provide upper bounds on the minimum length of $q$-ary
$(k,n)$-frameproof codes we need to recall a strictly related class of codes,
named $q$-ary $(k,n)$-strongly selective codes, studied in [7] by De Bonis and
Vaccaro.

######  Definition II.2.

Let $k$, $n$, and $q\geq 2$ be positive integers, $n\geq k$. A $q$-ary
$(k,n)$-_strongly selective_ code is a $t\times n$ matrix $M$ with entries in
$[0,q-1]$ such that for any $k$-tuple of the columns of $M$ and for any column
$\mathbf{c}$ of the given $k$-tuple, there exists a row ${i\in[1,t]}$ such
that $\mathbf{c}$ has an entry $s\in[1,q-1]$ in row $i$ whereas the entries in
row $i$ of all the remaining $k-1$ columns of the $k$-tuple belong to
$[0,q-1]\setminus\\{s\\}$. The number of rows $t$ of $M$ is called the length
of the $q$-ary $(k,n)$-strongly selective code.

Given an integer $q\geq 2$ and a $q$-ary vector ${\bf c}\in[0,q-1]^{t}$, we
denote with $w({\bf c})$ the number of nonzero components of ${\bf c}$.

######  Definition II.3.

A $q$-ary $(k,w,n)$-_strongly selective_ code is a $q$-ary $(k,n)$-strongly
selective code with the additional constraint that each column $\mathbf{c}$
has $w(\mathbf{c})=w$.

There is a strong relation between $q$-ary $(k,n)$-frameproof codes and
$q$-ary $(k+1,n)$-strongly selective codes. In fact, if we denote the minimum
length of $q$-ary $(k,n)$-frameproof codes by $t_{\text{FP}}(q,k,n)$ and that
of $q$-ary $(k,n)$-strongly selective codes by $t_{\text{SS}}(q,k,n)$, then
the following lemma holds.

######  Lemma II.4.

$t_{\text{SS}}(q,k+1,n)\geq
t_{\text{FP}}(q,k,n)\geq\frac{1}{2}t_{\text{SS}}(q,k+1,n)\,.$

###### Proof.

The upper bound on $t_{\text{FP}}(q,k,n)$ easily follows by noticing that a
$q$-ary $(k+1,n)$-strongly selective code of length $t$ is a $q$-ary
$(k,n)$-frameproof code of length $t$ by definition of such codes. Conversely,
a $q$-ary $(k,n)$-frameproof code $\mathcal{C}$ is not necessarily a $q$-ary
$(k+1,n)$-strongly selective code since we can have a column $\mathbf{c}$ and
other $k$ columns in $\mathcal{C}$ where the only rows there exists are those
in which the symbols of the $k$ columns do not contain the symbol of the
column $\mathbf{c}$ and where this symbol is equal to $0$. Therefore, we are
satisfying the property for frameproof codes but not the one for strongly
selective codes. However, the minimum length $t_{\text{SS}}(q,k+1,n)$ of
$q$-ary $(k+1,n)$-strongly selective codes is at most twice the minimum length
$t_{\text{FP}}(q,k,n)$ of $q$-ary $(k,n)$-frameproof codes. Indeed, if one is
given a $q$-ary $(k,n)$-frameproof code $\mathcal{C}$ then one can build a
$(k+1,n)$-strongly selective code twice as long as $\mathcal{C}$ by taking the
union of the rows of $\mathcal{C}$ and those of the complementary code
$\overline{\mathcal{C}}$ obtained by replacing each entry $s$ in $\mathcal{C}$
by $q-1-s$. Hence we obtain the lower bound on $t_{\text{FP}}(q,k,n)$ shown in
the statement of the lemma. ∎

Thanks to Lemma II.4 we can concentrate on finding bounds on the minimum
length $t_{\text{SS}}(q,k,n)$ and then obtain, indirectly, bounds on
$t_{\text{FP}}(q,k,n)$.

##  III A New Randomized Algorithm for frameproof codes via Lovász Local Lemma

In this section, we provide new upper bounds on the minimum length
$t_{\text{SS}}(q,k,n)$ of $q$-ary $(k,n)$-strongly selective codes. We will
provide two different bounds, one that is derived using the Lovász Local Lemma
while the other one using the expurgation method.

In [7], it has been proved the existence of $q$-ary $(k,n)$-strongly selective
of length

$t=O\left(\frac{k^{2}}{v}\log(n/k)\right),$ (1)

where $v=q-1$ for $q\leq k$ and $v=k$ for $q>k$. Their probabilistic
contruction produces a randomized algorithm of complexity $\Theta(n^{k})$ to
generate $q$-ary $(k,n)$-strongly selective codes of length of the same order
as the one shown in equation (1). The algorithm is clearly impractical already
for small values of $k$. Here we provide new bounds on $t_{\text{SS}}(q,k,n)$
that are of the same order as the ones provided by [7] but, crucially, this is
accompanied with a randomized construction algorithm of complexity
$O(tn^{2})$. This means that we can construct a $q$-ary $(k,n)$-strongly
selective codes of length as in (1) in time polynomial in $n$ and $k$.

The key idea is to use the following matrices (introduced in [12]) where the
constraints involve only pairs of columns.

######  Definition III.1.

Let $n,w,\lambda$ be positive integers. A $q$-ary $t\times n$ matrix $M$, with
entries in $[0,q-1]$, is a $(\lambda,w,n)$-matrix if the following properties
hold true:
1) each column ${\bf c}$ of $M$ has $w({\bf c})=w$,
2) any pair of columns $\mathbf{c},\mathbf{d}$ of $M$ have at most $\lambda$
nonzero symbols in common, that is, there are at most $\lambda$ rows where
columns $\mathbf{c}$ and $\mathbf{d}$ have both the same entry $s\in[1,q-1]$.

These matrices are related to $q$-ary $(k,w,n)$-strongly selective codes by
the following result.

######  Lemma III.2.

If $M$ is a $q$-ary $t\times n$ $(\lambda,w,n)$-matrix with
$\lambda=\left\lfloor(w-1)/(k-1)\right\rfloor$ then $M$ is a $q$-ary
$(k,w,n)$-strongly selective code of length $t$.

###### Proof.

Let $\mathbf{c}$ be an arbitrary column of $M$ and let $A$ be the set of row
indices in which column $\mathbf{c}$ has nonzero elements. Therefore $|A|=w$.
Let $K$ be a set of arbitrary $k-1$ columns of $M$ where $\mathbf{c}\not\in
K$. Let us denote by $M(A,K)$ the $w\times(k-1)$ submatrix of $M$ constructed
by first selecting the $k-1$ columns in $K$ and then selecting the $w$ rows
whose indices belong to $A$. Since $M$ is a $(\lambda,w,n)$-matrix we have
that the number of nonzero elements that column $\mathbf{c}$ share (in the
same rows) with any column in $K$ is at most
$\lambda=\left\lfloor\frac{w-1}{k-1}\right\rfloor$. Hence, the total number of
nonzero elements in $M(A,K)$ is at most

$(k-1)\lambda\leq(k-1)\left\lfloor\frac{w-1}{k-1}\right\rfloor\leq w-1\,.$

Considering that $M(A,K)$ has $|A|=w$ we have that at least one row in
$M(A,K)$ contains only zero elements. Then the lemma follows. ∎

Now, we want to give a good upper bound on the minimum length of $q$-ary
$(\lambda,w,n)$-matrices that will provide us an upper bound on the minimum
length of $q$-ary $(k,w,n)$-strongly selective codes. We first need to recall
the following well-known facts on binomial coefficients, where for positive
integers $c\leq b\leq a$ we have

$\left(\frac{a}{b}\right)^{b}\leq\binom{a}{b}\leq\frac{a^{b}}{b!}\leq\left(\frac{ea}{b}\right)^{b}\,,$
(2) $\binom{a}{b}\binom{b}{c}=\binom{a}{c}\binom{a-c}{b-c}\,,$ (3)

and the following useful inequality, proved in [9], for positive integers
$c\leq a\leq b$

$\binom{a}{c}\Big{/}\binom{b}{c}\leq\left(\frac{a-\frac{c-1}{2}}{b-\frac{c-1}{2}}\right)^{c}\,.$
(4)

Our main tool is the celebrated algorithmic version of the Lovász Local Lemma
for the symmetric case as given in [13].

######  Lemma III.3.

[13] Let $E_{1},E_{2},\ldots,E_{m}$ be events in a probability space, where
each event $E_{i}$ is mutually independent of the set of all the other events
$E_{j}$ except for at most $D$, and $\Pr(E_{i})\leq P$ for all $1\leq i\leq
m$. If $ePD\leq 1$, then $\Pr(\cap_{i=1}^{m}\overline{E_{i}})>0$. Moreover, a
configuration avoiding all events $E_{i}$ can be effectively found by using an
average number of resampling at most $m/D$.

We are now ready to state our main lemma.

######  Lemma III.4.

There exists a $q$-ary $t\times n$ $(\lambda,w,n)$-matrix with

$\displaystyle t$ $\displaystyle=\max\Bigg{\\{}\left\lceil
2w-(\lambda+1)\right\rceil,$
$\displaystyle\left\lceil\frac{\lambda}{2}+\frac{1}{q-1}\left(\frac{ew}{\lambda+1}\left(w-\frac{\lambda}{2}\right)(e(2n-4))^{\frac{1}{\lambda+1}}\right)\right\rceil\Bigg{\\}}.$
(5)

###### Proof.

Let $M$ be a random $t\times n$ $q$-ary matrix, $t\geq 2w-(\lambda+1)$, where
each column $\mathbf{c}$ is picked uniformly at random among the set of all
distinct $q$-ary vectors ${\bf c}$ of length $t$ such that $w({\bf c})=w$. It
is easy to see that the number of such vectors is equal to
$\binom{t}{w}(q-1)^{w}$, and therefore

$\Pr({\bf c})=\left(\binom{t}{w}(q-1)^{w}\right)^{-1}\,.$

Let $i,j\in[1,n],i\neq j$ and let us consider the event $\overline{E}_{i,j}$
that there exists _at most_ $\lambda$ rows such that both the $i$-th column
and the $j$-th column of $M$ have the same nonzero symbol in _each_ of these
rows. We evaluate the probability of the complementary “ _bad_ ” event
${E}_{i,j}$. Hence $E_{i,j}$ is the event that the $i$-th and $j$-th columns
${\bf c}_{i}$ and ${\bf c}_{j}$ have identical non-zero elements in at least
$\lambda+1$ coordinates. We bound $\Pr(E_{i,j})$ by conditioning on the event
that ${\bf c}_{i}$ is equal to $c$.

For a subset $S\subset[1,t]$ of coordinates, let $E_{i,j}^{S}$ be the event
that in each coordinate of $S$ the $i$-th and $j$-th column have identical
non-zero elements. Finally, let $A$ be the set of coordinates where ${\bf
c}_{i}$ is non-zero. Note that for $S\in\binom{A}{\lambda+1}$, i.e. $S$ is a
subset of $A$ of size $\lambda+1$, we have

$\Pr(E_{i,j}^{S}|{\bf
c}_{i}=c)=\frac{\binom{t-(\lambda+1)}{w-(\lambda+1)}(q-1)^{w-(\lambda+1)}}{\binom{t}{w}(q-1)^{w}}\,.$

Then

$\displaystyle\Pr(E_{i,j}|{\bf c}_{i}=c)$
$\displaystyle\leq\sum_{S\in\binom{A}{\lambda+1}}\Pr(E_{i,j}^{S}|{\bf
c}_{i}=c)$
$\displaystyle=\binom{w}{\lambda+1}\frac{{\binom{t-(\lambda+1)}{w-(\lambda+1)}}(q-1)^{w-(\lambda+1)}}{\binom{t}{w}(q-1)^{w}}.$
(6)

Since the right-hand side of (6) does not depend on the fixed column $c$, it
also holds unconditionally. Hence

$\Pr(E_{i,j})\leq\binom{w}{\lambda+1}\frac{{\binom{t-(\lambda+1)}{w-(\lambda+1)}}(q-1)^{w-(\lambda+1)}}{\binom{t}{w}(q-1)^{w}}.$
(7)

Therefore, by (7) we have

$\displaystyle\Pr({E}_{i,j})$
$\displaystyle\leq\binom{w}{\lambda+1}{\binom{t-(\lambda+1)}{w-(\lambda+1)}}\Big{/}\binom{t}{w}(q-1)^{\lambda+1}$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\binom{w}{\lambda+1}\binom{w}{\lambda+1}\Big{/}\binom{t}{\lambda+1}(q-1)^{\lambda+1}$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\frac{1}{(q-1)^{\lambda+1}}\binom{w}{\lambda+1}\left(\frac{w-\frac{\lambda}{2}}{t-\frac{\lambda}{2}}\right)^{\lambda+1}$
$\displaystyle\stackrel{{\scriptstyle(iii)}}{{\leq}}\frac{1}{(q-1)^{\lambda+1}}\left(\frac{ew}{\lambda+1}\right)^{\lambda+1}\left(\frac{w-\frac{\lambda}{2}}{t-\frac{\lambda}{2}}\right)^{\lambda+1}=P\,,$
(8)

where $(i)$ holds due to equality (3), $(ii)$ is true due to inequality (4),
and finally $(iii)$ holds thanks to inequalities (2).

The number of events ${E}_{i,j}$ is equal to $n(n-1)/2$. Let us fix an event
$E_{i,j}$, the number of events from which $E_{i,j}$ can be dependent is equal
to $D=2n-4$. Hence, according to Lemma III.3, if $P$ (as defined in (8)) and
$D=2n-4$ satisfy $ePD\leq 1$, then the probability that none of the “bad”
events $E_{i,j}$ occurs is strictly positive. One can see that by setting $t$
as in the second term of the maximum in (5) one indeed obtains

$ePD=\frac{e(2n-4)}{(q-1)^{\lambda+1}}\left(\frac{ew}{\lambda+1}\right)^{\lambda+1}\left(\frac{w-\frac{\lambda}{2}}{t-\frac{\lambda}{2}}\right)^{\lambda+1}\leq
1.$

Hence, from Lemma III.3 one can construct a $q$-ary $(\lambda,w,n)$-matrix $M$
whose number of rows $t$ satisfies equality (5) since we also need to consider
the initial assumption $t\geq 2w-(\lambda+1)$ so that all the computation
carried out in this lemma are meaningful. ∎

Now, thanks to Lemmas III.2 and III.4, we can prove the following result.

######  Theorem III.5.

There exists a randomized algorithm to construct a $q$-ary $(k,w,n)$-strongly
selective code with length

$\displaystyle t$ $\displaystyle\leq
1+\max\Bigg{\\{}2w-\frac{w-1}{k-1},\>\>\frac{w-1}{2(k-1)}+$
$\displaystyle\frac{ew(k-1)}{(q-1)(w-1)}\left(w-\frac{w-1}{2(k-1)}+\frac{1}{2}\right)(e(2n-4))^{\frac{k-1}{w-1}}\Bigg{\\}}.$
(9)

The algorithm requires, on average, time $O(tn^{2})$ to construct the code.

###### Proof.

The upper bound on $t$ shown in the statement of the lemma is derived by
substituting the value of $\lambda=\lfloor(w-1)/(k-1)\rfloor$ of Lemma III.2
into equation (5) of Lemma III.4, and by using the inequalities
$\frac{w-1}{k-1}-1\leq\left\lfloor\frac{w-1}{k-1}\right\rfloor\leq\frac{w-1}{k-1}$.
The time complexity $O(tn^{2})$ comes from Lemma III.3 by noticing that
$m/D=n(n-1)/(4n-8)\leq n/3$. Moreover, the algorithm that one obtains from
[13] requires to randomly generate a matrix, checking if the $\Theta(n^{2})$
events $\overline{E}_{i,j}$ are satisfied, and resample _only on non-
satisfied_ events. This means that we need to check if the $i$-th column and
the $j$-th column of the matrix have at most
$\left\lfloor\frac{w-1}{k-1}\right\rfloor$ nonzero elements in common, that
can be done with at most $O(t)$ operations, and resample only over non-
satisfied events. Then we need to check only the events that involve columns
that have been resampled. Altogether, by Lemma III.3 this procedure requires
$O(tn^{2}+n\cdot m/D\cdot t)=O(tn^{2})$ elementary operations. ∎

######  Remark III.6.

Before optimizing (9) over the parameter $w$, we would like to stress that
Theorem III.5 is not a mere technical intermediate result, but might be
important in several practical scenarios. For instance, it has been shown in
[7] that $q$-ary $(k,n)$-strongly selective codes can be used to solve
important communication problems arising in multiple-channel wireless
networks. In the scenario considered in [7], one has a set of $n$
uncoordinated stations, attempting transmission over a set of $q-1$
independent channels. Transmission is successful if and only if no two
stations attempt to transmit over the same channel at the same time instant.
The idea in [7] is the following: Each station is assigned a distinct codeword
of a $q$-ary $(k,n)$-strongly selective code, and such a codeword corresponds
to the transmission schedule of the associated station. The presence of a
symbol $s\in[0,q-1]$ in the $i$-th coordinate of a given codeword ${\bf c}$
naturally translates as the ”instruction”, to the station possessing codeword
${\bf c}$, to stay ”silent” in the $i$-th step of the communication protocol
if $s=0$, and to transmit over the $j$-th channel in the $i$-th step if
$s=j\neq 0$. Under the hypothesis that at any given time at most $k$ stations
are ”active”, the authors of [7] proved that $q$-ary $(k,n)$-strongly
selective codes naturally correspond to conflict resolution protocols in
multiple-channel wireless networks. Now, in many situations, it is important
not only to minimize the length of the protocol (i.e., the number of time
instants before all stations transmit successfully) but it is important also
to restrict the number of attempted transmissions by each station (to save
energy, for example, see [9] for more). Therefore $q$-ary $(k,w,n)$-strongly
selective codes, where each codeword has $w$ non-zero components, could be
useful in these instances. We remark that the techniques of [7] are not able
to deal with this hard constraint on the number of attempted transmissions by
each station, nor suggest efficient algorithms for the construction of $q$-ary
$(k,n)$-strongly selective codes, as our technique does.

We optimize $w$ in equation (9) to get a randomized algorithm for
(unconstrained) $q$-ary $(k,n)$-strongly selective codes.

######  Corollary III.7.

Let $w=\lceil 1+(k-1)\ln(2en)\rceil$, then the algorithm described in Theorem
III.5 constructs a $q$-ary $(k,w,n)$-strongly selective code that is, clearly,
a $q$-ary $(k,n)$-strongly selective code of length $t$ upper bounded as

$\displaystyle t$
$\displaystyle\leq\max\Bigg{\\{}2(k-1)\ln(2en)-\ln(n),\>\>\frac{\ln(n)}{2}+$
$\displaystyle\frac{e^{2}(k-1)^{2}}{q-1}\ln(2en)+\frac{7e^{2}(k-1)}{2(q-1)}\Bigg{\\}}+O(1).$

###### Proof.

Fixing $w$ as in the statement of the corollary and using the following
inequalities

$1+(k-1)\ln(2en)\leq\lceil 1+(k-1)\ln(2en)\rceil\leq 2+(k-1)\ln(2en)$

by Theorem III.5, considering only the second term of the maximum, we have
that

$\displaystyle t$ $\displaystyle\leq
1+\frac{1+(k-1)\ln(2en)}{2(k-1)}+\frac{e}{(q-1)\ln(2en)}$
$\displaystyle(2+(k-1)\ln(2en))\left(\frac{3}{2}+\left(k-1\right)\ln(2en)\right)(2en)^{\frac{1}{\ln(2en)}}$
$\displaystyle\leq\frac{\ln(n)}{2}+\frac{e^{2}(k-1)^{2}}{q-1}\ln(2en)+\frac{7e^{2}(k-1)}{2(q-1)}+O(1)\,,$

since $k\geq 2$, $\ln(2en)\geq 2$ and $(2en)^{\frac{1}{\ln(2en)}}=e$. Then,
the corollary follows since we also need to consider the first term of the
maximum of Theorem III.5. ∎

Therefore due to Lemma II.4 and Corollary III.7 we obtain the main result of
this paper.

######  Theorem III.8.

There exists a randomized algorithm to construct a $q$-ary $(k,n)$-frameproof
code with length

$\displaystyle t\leq\max$
$\displaystyle\Bigg{\\{}2k\ln(2en)-\ln(n),\>\>\frac{\ln(n)}{2}+$
$\displaystyle\frac{e^{2}k^{2}}{q-1}\ln(2en)+\frac{7e^{2}k}{2(q-1)}\Bigg{\\}}+O(1).$

The algorithm requires, on average, time $O(tn^{2})$ to construct the code.

To properly judge the value of Theorem III.8, we recall the following result
that provides a lower bound on the length of any $q$-ary $(k,n)$-frameproof
code.

######  Theorem III.9.

[6, 7, 15] Given positive integers $q,k$, and $n$, with $q\geq 2$ and $2\leq
k\leq\sqrt{n}$, the minimum length of any $q$-ary $(k,n)$-frameproof code
satisfies

$t_{\text{FP}}(q,k,n)=\Omega\left(\frac{k^{2}}{q\log
k}\log\frac{n}{k}\right).$ (10)

Therefore, one can see that the construction method provided by Theorem III.8,
besides being quite efficient, produces codes of almost optimal length.

We can also prove the following result.

######  Theorem III.10.

Let $M$ be a $q$-ary $(k,n)$-frameproof code of _minimum_ length $t$. Then

$\left\lceil\frac{n}{q-1}\right\rceil\geq
t\geq\left\lceil\frac{1}{q}\min\left\\{n,\frac{15+\sqrt{33}}{24}k^{2}\right\\}\right\rceil\,.$
(11)

###### Proof.

We can define a map from $M$ to a binary $(k,n)$-strongly selective code of
length $t^{\prime}$ by mapping each symbol $i\in[0,q-1]$ into
$\mathbf{e}_{i+1}$, where $\mathbf{e}_{i}$ is the binary column vector that
has $1$ in the $i$-th component and $0$ elsewhere. Therefore $t^{\prime}=qt$.
Now, the right-hand side inequality of (11) follows since
${t^{\prime}\geq\min\left\\{n,\frac{15+\sqrt{33}}{24}k^{2}\right\\}}$ by [14,
Theorem 2]. To prove the left-hand side inequality of (11), simply observe
that by taking any $n$ columns from the side-by-side concatenation of $(q-1)$
many $\lceil n/(q-1)\rceil\times\lceil n/(q-1)\rceil$ diagonal matrices, where
the $i$-th matrix has symbol $i\in[1,q-1]$ on its diagonal, and 0 elsewhere,
one gets a $q$-ary $(k,n)$-frameproof code of length $\lceil n/(q-1)\rceil$. ∎

##  IV New improved upper bounds on $t_{\text{FP}}(q,k,n)$ via Expurgation
Method

In this Section, we will provide an existential upper bound on the minimum
length of $q$-ary $(k,n)$-frameproof codes that improves the best results
known in the literature. We first recall such known results.

######  Theorem IV.1.

[19] There exists a $q$-ary $(k,n)$-frameproof code of length

$t\leq-k\ln\left(n\frac{k!}{k!-1}\right)\Big{/}\ln\left(1-\left(1-\frac{1}{q}\right)^{k}\right)\,.$

We note that in the original paper [19], this bound is stated in terms of
bounds on the length of separating hash families. However, it is well known
that a separating hash family of type $(t,n,q,\\{1,k\\})$ is equivalent to a
$q$-ary $(k,n)$-frameproof code of length $t$. In [15], the authors provide
the following bound that improves the one of Theorem IV.1 whenever $q\leq k$
and $k$ is sufficiently large.

######  Theorem IV.2.

[15] If $q\leq k$, then there exists a $q$-ary $(k,n)$-frameproof code of
length

$t\leq\frac{-k\ln n-(k+1)\ln
2}{\ln\left[1-\left(1-\frac{q-1}{k+1}\right)\left(\frac{q-1}{k+1}\right)^{k}-\frac{q-1}{k+1}\left(1-\frac{1}{k+1}\right)^{k}\right]}\,.$

Our result employs the so-called expurgation method or deletion method, the
same technique of [19], but with a more careful optimization of the
parameters. Therefore, we are able to improve the bound of Theorem IV.1 for
every value of the parameters $k$, $q$, and $n$. Here we state the following
theorem.

######  Theorem IV.3.

There exists a $(k,n)$-frameproof code of length $t$, where $t$ is the minimum
integer such that the following inequality holds

$(k+1)\binom{n\left(1+\frac{1}{k}\right)}{k}(1-p_{q,k})^{t}\leq 1\,,$

where

$p_{q,k}:=\begin{cases}\left(1-\frac{1}{q}\right)^{k}&\text{for }q>k,\\\
\left(1-\frac{q-1}{k+1}\right)\left(\frac{q-1}{k+1}\right)^{k}+\frac{q-1}{k+1}\left(1-\frac{1}{k+1}\right)^{k}&\text{otherwise.}\end{cases}$

###### Proof.

Let $M$ be a $q$-ary $t\times(n+\ell)$ matrix, where each element is picked
i.i.d. at random in the set $[0,q-1]$ with distribution
$\mathbf{\mu}=(\mu_{0},\mu_{1},\ldots,\mu_{q-1})$ that will be fixed later.
For a given index $i\in[1,n]$ and a set of column-indices $B$, $|B|=k-1$,
$i\not\in B$, let $E_{i,B}$ be the event such that for every row in which
$\mathbf{c}_{i}$ (the $i$-th column) has a symbol $s$, there exists an index
$j\in B$ such that $\mathbf{c}_{j}$ has symbol $s$ in that same row. Therefore
the probability of each event $E_{i,B}$ can be upper bounded as
$\Pr(E_{i,B})\leq\left(1-\sum_{i=0}^{q-1}\mu_{i}\left(1-\mu_{i}\right)^{k}\right)^{t}$.
The number of such events is equal to $(n+\ell)\binom{n+\ell-1}{k}$. Now, let
$X$ be the random variable that represents the number of events $E_{i,B}$ that
are satisfied. Hence, taking $\mu_{i}=1/q$ for every $i\in[0,q-1]$ (uniform
distribution) when $q>k$ and $\mu_{i}=1/(k+1)$ for every $i\in[1,q-1]$,
$\mu_{0}=1-(q-1)/(k+1)$ when $q\leq k$, we obtain

$\mathbb{E}[X]\leq(n+\ell)\binom{n+\ell-1}{k}\left(1-p_{q,k}\right)^{t}.$

We note that if $\mathbb{E}[X]<\ell+1$ then there exists at most $\ell$ “
_bad_ ” events $E_{i,B}$ that are satisfied. Then, for each of these events
$E_{i,B}$ we remove one column with index in $\\{i\\}\cup B$. Hence we are
left with a $q$-ary matrix with $t$ rows and at least $n$ columns that satisfy
the frameproof property. Therefore we obtain a $(k,n)$-frameproof code with
length $t$. Thus the theorem follows taking $\ell=\lfloor n/k\rfloor$. ∎

######  Corollary IV.4.

Using the inequalities in equation (2), we have that, from Theorem IV.3, the
length of $(k,n)$-frameproof codes is upper bounded as follows

$t\leq\frac{-k\ln\left(n\frac{k+1}{k}\right)-\ln\left(\frac{k+1}{k!}\right)}{\ln\left(1-p_{q,k}\right)}\,,$

where $p_{q,k}$ is the same quantity defined in Theorem IV.3.

######  Theorem IV.5.

The bound of Corollary IV.4 improves the one of Theorem IV.1 for every $n\geq
k\geq 2$ and $q>k$.

###### Proof.

We need to prove the following inequality

$\frac{-k\ln\left(n\frac{k+1}{k}\right)-\ln\left(\frac{k+1}{k!}\right)}{\ln\left(1-p_{q,k}\right)}<\frac{-k\ln\left(n\frac{k!}{k!-1}\right)}{\ln\left(1-\left(1-\frac{1}{q}\right)^{k}\right)}\,.$
(12)

Since for $q>k$, $p_{q,k}=\left(1-1/q\right)^{k}$, we can rearrange and
simplify the terms in (12) to obtain the following inequality.

$\left(\frac{k+1}{k}\right)^{k}\frac{k+1}{k!}<\left(\frac{k!}{k!-1}\right)^{k}\,.$
(13)

Now, since $\left(1+\frac{1}{k}\right)^{k}\leq e$ and $k!\geq
e\left(\frac{k}{e}\right)^{k}$, the left-hand side of (13) is smaller than
$\left(\frac{e}{k}\right)^{k}(k+1)$. To prove inequality (13), it suffices to
show that

$\left(\frac{e}{k}\right)^{k}(k+1)<1\,,$ (14)

since the right-hand side of (13) is greater than $1$ for every $k$.

It can be seen that the left-hand side of (14) is a decreasing function in $k$
for $k\geq 2$. The first integer $k$ for which inequality (14) holds is $k=5$.
Then, the theorem follows since inequality (13) also holds for ${k=2,3,4}$ by
direct computation. ∎

######  Theorem IV.6.

The bound of Corollary IV.4 improves the one of Theorem IV.2 for every $n\geq
k\geq 2$ and $q\leq k$.

###### Proof.

Clearly, we need to prove the following inequality

$\frac{-k\ln\left(n\frac{k+1}{k}\right)-\ln\left(\frac{k+1}{k!}\right)}{\ln\left(1-p_{q,k}\right)}\\\
<\frac{-k\ln n-(k+1)\ln
2}{\ln\left[1-\left(1-\frac{q-1}{k+1}\right)\left(\frac{q-1}{k+1}\right)^{k}-\frac{q-1}{k+1}\left(1-\frac{1}{k+1}\right)^{k}\right]}\,.$
(15)

For $q\leq k$, the denominators in (15) are equal by definition of $p_{q,k}$.
Therefore, we can rearrange and simplify the terms to obtain the following
inequality.

$\left(\frac{k+1}{k}\right)^{k}\frac{k+1}{k!}<2^{k+1}\,.$ (16)

Proceeding as in the proof of Theorem IV.5, the left-hand side of (16) is a
decreasing function in $k$ for $k\geq 2$ while the right-hand side of (16) is
increasing in $k$. By inspection, it is easy to see that inequality (16) holds
even for $k=2$. ∎

## Acknowledgements

The work of A. A. Rescigno and U. Vaccaro was partially supported by project
SERICS (PE00000014) under the NRRP MUR program funded by the EU–NGEU.

## References

  * [1] A. Barg, G. R. Blakley, G. A. Kabatiansky, Digital fingerprinting codes: Problem statements, constructions, identification of traitors. IEEE Transactions on Information Theory 49(4), pp. 852–865, 2003.
  * [2] S.R. Blackburn, Frameproof codes. Siam J. Discrete Mathematics, 13(3), pp. 499–510, 2003.
  * [3] D. Boneh and J. Shaw, Collusion-secure fingerprinting for digital data. IEEE Transactions on Information Theory 44(5), pp.1897–1905, 1998.
  * [4] Y. Chee and X. Zhang, Improved constructions of frameproof codes, IEEE Transactions on Information Theory 58(8), pp. 5449–5453, 2012.
  * [5] M. Cheng and Y. Miao, On anti-collusion codes and detection algorithms for multimedia fingerprinting, IEEE Transactions on Information Theory, 57(7), pp. 4843-4851, 2011,
  * [6] A. De Bonis and U. Vaccaro, A New Kind of Selectors and Their Applications to Conflict Resolution in Wireless Multichannels Networks. In: Chrobak, M., Fernández Anta, A., Gasieniec, L., Klasing, R. (eds) Algorithms for Sensor Systems. ALGOSENSORS 2016. Lecture Notes in Computer Science, vol 10050. Springer, 2016.
  * [7] A. De Bonis and U. Vaccaro, A new kind of selectors and their applications to conflict resolution in wireless multichannels networks, Theoretical Computer Science 806, pp. 219–235, 2020.
  * [8] M. Fernandez and J. Livieratos, Algorithmic aspects on the construction of separating codes, Analysis of Experimental Algorithms: Special Event, Springer International Publishing, 2019.
  * [9] L. Gargano, A. A. Rescigno and U. Vaccaro, Low-weight superimposed codes and related combinatorial structures: Bounds and applications, Theoretical Computer Science 806, pp. 655–672, 2020.
  * [10] C. Guo, D. R. Stinson and T. Van Trung, ”On tight bounds for binary frameproof codes”, Design, Codes, and Cryptography 77(2), pp. 301-319, 2015.
  * [11] H. D. Hollmann, J. H. Van Lint, J. P. Linnartz and L. M. Tolhuizen, On codes with the identifiable parent property, Journal of Combinatorial Theory, Series A 82(2), pp. 121–133, 1998.
  * [12] W. Kautz and R. Singleton, Nonrandom binary superimposed codes, IEEE Transactions on Information Theory 10, pp. 363–377, 1964.
  * [13] R.A. Moser and G. Tardos, A constructive proof of the general Lovász local lemma, Journal of the ACM 57, pp. 1–15, 2010.
  * [14] C. Shangguan and G. Ge, New Bounds on the Number of Tests for Disjunct Matrices, IEEE Transactions on Information Theory 61(12), pp. 7518-7521, 2016.
  * [15] C. Shangguan, X. Wang, G. Ge and Y. Miao, New bounds for frameproof codes. IEEE Transactions on Information Theory 63, pp. 7247–7252, 2017.
  * [16] J. N. Staddon, D. R. Stinson and R. Wei, ”Combinatorial properties of frameproof and traceability codes”, IEEE Transactions on Information Theory 47(3), pp. 1042-1049, 2001.
  * [17] D. R. Stinson and R. Wei, ”Combinatorial properties and constructions of traceability schemes and frameproof codes”, Siam J. Discrete Mathematics 11(1), pp. 41–53, 1998.
  * [18] D. R. Stinson, R. Wei and K. Chen, On generalized separating hash families. Journal of Combinatorial Theory, Series A 115, pp. 105-120, 2008.
  * [19] D. R. Stinson and G. M. Zaverucha, Some improved bounds for secure frameproof codes and related separating hash families, IEEE Transactions on Information Theory 54(6), pp. 2508–2514, 2008.
  * [20] T. Van Trung, ”A tight bound for frameproof codes viewed in terms of separating hash families”, Design, Codes, and Cryptography 72(3), pp. 713–718, 2014.


# Auxiliary open quantum system for the Floquet quantum master equation

Fei Liu<EMAIL_ADDRESS>School of Physics, Beihang University, Beijing
100191, China

###### Abstract

By directly using the probability formulas of quantum trajectories, we
construct an auxiliary open quantum system for a periodically driven open
quantum system whose dynamics is governed by the Floquet quantum master
equation. This auxiliary system can generate a quantum trajectory ensemble
that is consistent with the canonical quantum trajectory ensemble. We find
that, at a long time limit, though the Lindblad operators are modified, the
coherent dynamics of the auxiliary system is the same as that of the original
system. A periodically driven two-level quantum system is used to illustrate
this construction.

##  I Introduction

In the past decade, microcanonical and canonical trajectory ensembles of
stochastic systems have attracted considerable interest Evans (2004a);
Garrahan  _et al._ (2007); Maes and Netočný (2008); Maes  _et al._ (2008);
Garrahan  _et al._ (2009); Garrahan and Lesanovsky (2010); Jack and Sollich
(2010); Nemoto and Sasa (2011); Chetrite and Touchette (2013, 2014); Derrida
and Sadhu (2019). The underlying reason is that the rare events or rare large
fluctuations generated by these trajectory ensembles are crucial for
understanding and quantizing many physical contexts, e.g., sheared fluids
Evans (2004a, 2010, b), dynamic phase transitions in glasses Garrahan  _et
al._ (2007, 2009) and open quantum systems Garrahan and Lesanovsky (2010);
Garrahan  _et al._ (2011); Carollo  _et al._ (2018), and fluctuation theorems
in nonequilibrium processes Lebowitz and Spohn (1999); Kurchan (2000); Maes
(1999); Gallavotti and Cohen (1995); Jarzynski (1997); Crooks (2000); Seifert
(2005).

Generally, microcanonical trajectory ensembles can be regarded as a set of
trajectories in which time-integrated observables defined on trajectories are
constrained by given values. In contrast, in a canonical trajectory ensemble,
constraints are imposed on the trajectory ensemble-averaged values of
observables. Although these two ensembles seem to be very distinct, at long
time limits and under certain conditions, they have been rigorously proven to
be equivalent by using large deviation theory Chetrite and Touchette (2013,
2014); Touchette (2008). Importantly, auxiliary stochastic systems, including
continuous diffusion Chetrite and Touchette (2013, 2014) and discrete jump
stochastic systems Jack and Sollich (2010); Garrahan  _et al._ (2009, 2007),
have been explicitly constructed whose trajectories consist of the
trajectories of canonical or microcanonical trajectory ensembles.

Recently, Carollo et al Carollo  _et al._ (2018) extended the results of
previous research on classical stochastic systems to open quantum systems. The
dynamics of open quantum systems obeys the time-independent quantum master
equations. It has been well established that the quantum master equations can
be unraveled into quantum trajectories and, in particular, that these
trajectories have a classic probability interpretation Breuer and Petruccione
(2002); Carmichael (1993); Wiseman and Milburn (2010); Liu (2018). Therefore,
this extension is not surprising. Inspired by their work, in this paper, we
investigate the construction of an auxiliary quantum system for a periodically
driven open quantum system Kosloff (2013); Szczygielski  _et al._ (2013);
Gasparinetti  _et al._ (2014); Cuetara  _et al._ (2015); Gelbwaser-Klimovsky
_et al._ (2015); Liu and Su (2020). Different from previous stationary
equations, the dynamics of these open quantum systems is governed by the
Floquet quantum master equation Blümel  _et al._ (1991); Kohler  _et al._
(1997); Breuer and Petruccione (1997); Alicki  _et al._ (2006). Our physical
situation of interest is very distinct from that of Carollo et al. Carollo
_et al._ (2018); moreover, our theory is directly based on the notion of
quantum trajectories and thus explicitly uses the probability formulas of
quantum trajectories. Compared with the previous method that employs abstract
operator calculations, our method is more similar to that of Chetrite and
Touchette on classical stochastic systems Chetrite and Touchette (2014).

The remainder of this paper is organized as follows. In Sec. (II), we review
the Floquet quantum master equation and the unraveling of its quantum
trajectories; the essential notations are defined therein. In Sec. (III), we
define the canonical quantum trajectory ensemble. In Sec. (IV), an auxiliary
quantum system that can generate a canonical quantum trajectory ensemble is
constructed for the finite time case and long time limit case. In Sec. (V), a
two-level quantum system is used to illustrate the procedure of constructing
the auxiliary system in the long time limit case. Section (VI) concludes the
paper.

##  II Floquet quantum master equation and trajectory unraveling

Assume that a periodically modulated quantum system interacts with a heat bath
and that the inverse temperature of the heat bath is $\beta$. The Hamiltonian
of the quantum system is $H(t)$ and $H\left(t+2\pi/\Omega\right)=H(t)$, where
$\Omega$ is the driving frequency. According to the Floquet theorem Zeldovich
(1967); Shirley (1965), the Hamiltonian satisfies an eigenvalue equation:

$\displaystyle{\cal H}(t)|u_{n}(t)\rangle=\epsilon_{n}|u_{n}(t)\rangle,$ (1)

where ${\cal H}(t)=H(t)-i\hbar\partial_{t}$ is the Floquet Hamiltonian and
$\epsilon_{n}$ and $|u_{n}(t)\rangle$ ($n$$=$$1,\cdots,N$) are the quasi-
energy and Floquet bases, respectively. Because the quasi-energy
$\epsilon_{n}-q\hbar\Omega$ with basis $\exp(-iq\Omega t)|u_{n}(t)\rangle$ is
also the solution of Eq. (1), where $q$ is an arbitrary integer, we restrict
the quasi-energies in a zone of size $\hbar\Omega$. Under the assumption of a
weak system-bath coupling condition and appropriate time-scale separation, the
evolution of the reduced density matrix of the quantum system $\rho(t)$ can be
described by the Floquet quantum master equation Grifoni and Hänggi (1998);
Breuer and Petruccione (1997); Alicki  _et al._ (2006):

$\displaystyle\partial_{t}\rho(t)={\cal
L}(t)[\rho(t)]=-\frac{i}{\hbar}[H(t),\rho(t)]+D(t)[\rho(t)].$ (2)

The $D(t)$ term in the generator ${\cal L}(t)$ represents dissipation and
dephasing due to the interaction between the system and the heat bath and is
expressed as

$\displaystyle D(t)[\rho]=\sum_{\omega}r(\omega)\left[A(\omega,t)\rho
A^{\dagger}(\omega,t)-\frac{1}{2}\left\\{A^{\dagger}(\omega,t)A(\omega,t),\rho\right\\}\right],$

where the summation is performed with respect to all possible Bohr frequencies
$\omega$, which equal $(\epsilon_{n}-\epsilon_{m})/\hbar+q\Omega$, and $q$ is
a certain integer. The Bohr frequencies may be positive or negative but always
appear in pairs. In the same equation, $A(\omega,t)$ and
$A^{\dagger}(\omega,t)$ are the Lindblad operators and are related by
$A^{\dagger}(\omega,t)=A(-\omega,t)$ due to the Hermitian characteristics of
the interaction Hamiltonian. Note that the Lindblad operators are also the
eigenoperators of the Floquet Hamiltonian Breuer and Petruccione (1997):

$\displaystyle[{\cal H}(t),A(\omega,t)]=-\hbar\omega A(\omega,t),\hskip
5.69046pt[{\cal H}(t),A^{\dagger}(\omega,t)]=\hbar\omega
A^{\dagger}(\omega,t).$ (3)

Because the heat bath is always in the thermal state with the inverse
temperature $\beta$, the Fourier transformation $r(\omega)$ of the correlation
function of the heat bath satisfies the Kubo-Martin-Schwinger (KMS) condition
Breuer  _et al._ (2000): $r(-\omega)=r(\omega)\exp(-\beta\hbar\omega)$.

Eq. (2) can be unraveled into the dynamics of individual quantum systems
Carmichael (1993); Plenio and Knight (1998); Breuer  _et al._ (2000); Wiseman
and Milburn (2010). The evolution of each system is alternately composed of
continuous processes and discrete random jumps. Assume that the jumps occur at
time $t_{i}$ with Bohr frequency $\omega_{i}$, where $i=1,\cdots,M$ and $M$
denotes the total number of jumps. When the evolution ends at time $T$, a
quantum trajectory is generated and is denoted as
$\overrightarrow{\omega}_{M}=\\{\omega_{1},\cdots,\omega_{M}\\}$. Further
assuming that the density matrixes of these individual quantum systems at time
$t$ ($\leq T$) are $\widetilde{\rho}(\overrightarrow{\omega}_{M},t)$, the
solution of Eq. (2) is equal to a quantum trajectory ensemble average:

$\displaystyle\rho(t)$ $\displaystyle=$
$\displaystyle\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{M=0}^{\infty}{\cal
D}(t)\hskip
2.84544ptp(\overrightarrow{\omega}_{M},t){\widetilde{\rho}}(\overrightarrow{\omega}_{M},t),$
(12)

where the sum-integral symbol means that we sum with respect to all possible
quantum trajectories and integrate over all possible times Liu and Xi (2016);
Liu (2018); Liu and Su (2020),

$\displaystyle{\widetilde{\rho}}(\overrightarrow{\omega}_{M},t)=\frac{1}{p(\overrightarrow{\omega}_{M},t)}G_{t,t_{M}}\left(J_{t_{M}}\left(G_{t_{M},t_{M-1}}\left(\cdots
J_{t_{1}}\left(G_{t_{1},0}\left(\rho_{0}\right)\right)\cdots\right)\right)\right),$
(13)

where $\rho_{0}$ is the density matrix at initial time $0$ and the denominator
$p(\overrightarrow{\omega}_{M},t)$ is the probability distribution of
observing the quantum trajectory $\overrightarrow{\omega}_{M}$, which is
simply equal to the trace of the numerator. In Eq. (13), the symbols $G$ and
$J$ denote superoperators, and their action regions are explicitly indicated
by the round brackets. The superoperator $G_{t_{i},t_{i-1}}={\cal
T}_{-}\exp[\int_{t_{i-1}}^{t_{i}}d\tau{L}_{0}(\tau)]$ and its generator is

$\displaystyle L_{0}(t)[O]$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}[H(t),O]-\frac{1}{2}\sum
r(\omega)\left\\{A^{\dagger}(\omega,t)A(\omega,t),O\right\\},$ (14)

where ${\cal T}_{-}$ is the time-ordering operator. $J_{t_{i}}(O)$ is a
shorthand expression of the jump superoperator
$r(\omega_{i})A(\omega_{i},t_{i})OA^{\dagger}(\omega_{i},t_{i})$.

The quantum trajectory ensemble that we are interested in is as follows. At
initial time $0$, the individual quantum systems are in one of the Floquet
pure states $|u_{m}(0)\rangle\langle u_{m}(0)|$, which is randomly selected
based on a probability distribution $p_{m}$, $m=1,\cdots,N$. That is, the
quantum ensemble is initially in a mixed state:
$\rho_{0}=\sum_{m=1}^{N}p_{m}|u_{m}(0)\rangle\langle u_{m}(0)|$. The quantum
system evolves, and a quantum trajectory $\overrightarrow{\omega}_{M}$ is
recorded. At the end time $T$, the system is measured in the Floquet basis,
and a pure state $|u_{n}(T)\rangle\langle u_{n}(T)|$ is obtained. It is not
difficult to argue that, if the joint probability distribution of observing
the quantum trajectory with the special initial and terminal bases is
$p_{m,n}(\overrightarrow{\omega}_{M},T)=p_{m}p_{m|n}(\overrightarrow{\omega}_{M},T)$,
the conditional probability distribution equals Liu and Xi (2016); Liu (2018)

$\displaystyle p_{m|n}(\overrightarrow{\omega}_{M},T)=\langle
u_{n}(T)|G_{T,t_{M}}\left(J_{t_{M}}\left(G_{t_{M},t_{M-1}}\left(\cdots
J_{t_{1}}\left(G_{t_{1},0}\left(|u_{m}(0)\rangle\langle
u_{m}(0)|\right)\right)\cdots\right)\right)\right)|u_{n}(T)\rangle.$ (15)

##  III Canonical quantum trajectory ensemble

In this paper, we investigate canonical quantum trajectory ensembles Garrahan
_et al._ (2009); Chetrite and Touchette (2013, 2014). To define such an
ensemble, we choose a time-integrated observable of interest, namely, the
stochastic heat current along quantum trajectories. According to the
interpretation of quantum measurements Breuer and Petruccione (1997, 2002);
Liu (2018), the occurrence of jump events along a quantum trajectory
$\overrightarrow{\omega}_{M}=\\{\omega_{1},\cdots,\omega_{M}\\}$ indicates
that the quantum system exchanges quanta $\hbar\omega_{i}$ with the heat bath.
From a thermodynamic point of view, these quanta represent the discrete heat
released to the environment Breuer (2003); Roeck and Maes (2006); Crooks
(2008); Horowitz (2012); Hekking and Pekola (2013); Manzano  _et al._ (2015);
Liu and Xi (2016). Hence, given a quantum trajectory
$\overrightarrow{\omega}_{M}$ with duration $T$, the heat current is equal to
$j=Q_{T}/T$, and the total heat production is

$\displaystyle
Q_{T}(\overrightarrow{\omega}_{M})=\sum_{i=1}^{M}\hbar\omega_{i}.$ (16)

Note that our object of interest is different from the average rate of quantum
jumps, which was the focus of Carollo et al. Carollo  _et al._ (2018).
Nevertheless, the extension of our formulas to the latter case is direct.
Given the observable, in a canonical quantum trajectory ensemble, the
probability distribution of the quantum trajectory
$\overrightarrow{\omega}_{M}$ with the initial Floquet basis
$|u_{m}(0)\rangle$ and terminal Floquet basis $|u_{n}(T)\rangle$ is Garrahan
_et al._ (2009); Chetrite and Touchette (2013, 2014)

$\displaystyle
p_{m,n}({\chi},\overrightarrow{\omega}_{M},T)=\frac{1}{\Phi_{T}(\chi)}p_{m,n}(\overrightarrow{\omega}_{M},T)e^{\chi
Q_{T}(\overrightarrow{\omega}_{M})},$ (17)

where the normalized factor $\Phi_{T}(\chi)$ is the moment generating
function,

$\displaystyle\Phi_{T}(\chi)$ $\displaystyle=$
$\displaystyle\sum_{m,n=1}^{N}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{M=0}^{\infty}{\cal
D}(t)\hskip 2.84544ptp_{m,n}(\overrightarrow{\omega}_{M},T)e^{\chi
Q_{T}(\overrightarrow{\omega})}.$ (26)

##  IV Auxiliary open quantum system

###  IV.1 Finite time case

Now, we are in a position to construct an auxiliary open quantum system whose
quantum trajectory ensemble is consistent with Eq. (17). Assume the
Hamiltonian and Lindblad operators of the auxiliary system are $H^{\prime}(t)$
and $A^{\prime}(\omega,t)$, respectively. In the remainder of this paper, we
denote quantities in the auxiliary system with a prime symbol unless otherwise
stated. Inspired by Carollo et al. Carollo  _et al._ (2018), we set

$\displaystyle
A^{\prime}(\omega,t)=l_{\chi}(t)A(\omega,t)l_{\chi}^{-1}(t)e^{\chi\hbar\omega/2},$
(27)

where $l_{\chi}(t)$ is an invertible Hermitian operator to be solved 111This
result is easy to obtain if one considers the consistency between the
classical jump master equation and the diagonal part of the quantum master
equation (2) represented in the Floquet basis Garrahan  _et al._ (2009); Jack
and Sollich (2010). . We preliminarily require that the operator depends on
the parameter $\chi$. Obviously, Eq. (27) implies that the jump superoperator
of the auxiliary system is related to that of the original system as follows:

$\displaystyle
J^{\prime}_{t_{i}}(O)=e^{\chi\hbar\omega_{i}}l_{t_{i}}J_{t_{i}}(l^{-1}_{t_{i}}Ol^{-1}_{t_{i}})l_{t_{i}}.$
(28)

For simplicity, we denote $l_{t_{i}}=l_{\chi}(t_{i})$ and
$l_{t_{i}}^{-1}=l_{\chi}^{-1}(t_{i})$. Substituting Eq. (28) into
$p^{\prime}_{m,n}(\overrightarrow{\omega}_{M},T)$, we obtain the probability
distribution of observing the quantum trajectory $\overrightarrow{\omega}_{M}$
with the initial basis $|u_{m}(0)\rangle$ and terminal basis
$|u_{n}(T)\rangle$ in the auxiliary open quantum system:

$\displaystyle p^{\prime}_{m,n}(\overrightarrow{\omega}_{M},T)$
$\displaystyle=$ $\displaystyle e^{\chi Q(\overrightarrow{\omega}_{M})}\langle
u_{n}(T)|l_{T}l^{-1}_{T}G^{\prime}_{T,t_{M}}\left(l_{t_{M}}J_{t_{M}}\left(l_{t_{M}}^{-1}G^{\prime}_{t_{M},t_{M-1}}\left(l_{t_{M-1}}\right.\right.\right.$
(29) $\displaystyle\left.\left.\left.\cdots
J_{t_{1}}\left(l_{t_{1}}^{-1}G^{\prime}_{t_{1},0}\left(l_{0}l^{-1}_{0}|u_{m}(0)\rangle\langle
u_{m}(0)|l^{-1}_{0}l_{0}\right)l_{t_{1}}^{-1}\right)\cdots
l_{t_{M-1}}\right)l_{t_{M}}^{-1}\right)l_{t_{M}}\right)l^{-1}_{T}l_{T}|u_{n}(T)\rangle
p_{m}.$

Compared with Eq. (17), we find that, if

$\displaystyle
l^{-1}_{t_{i}}G^{\prime}_{t_{i},t_{i-1}}(l_{t_{i-1}}Ol_{t_{i-1}})l^{-1}_{t_{i}}=e^{-\Lambda_{T}(\chi)(t_{i}-t_{i-1})}G_{t_{i},t_{i-1}}(O),$
(30) $\displaystyle\hskip
5.69046ptl_{\chi}(0)|u_{m}(0)\rangle=|u_{m}(0)\rangle,\hskip 5.69046pt{\rm
and}\hskip 5.69046ptl_{\chi}(T)|u_{n}(T)\rangle=|u_{n}(T)\rangle,$ (31)

where

$\displaystyle\Lambda_{T}(\chi)=\frac{1}{T}\ln\Phi_{T}(\chi)$ (32)

is the scaled cumulant generating function Touchette (2008), the auxiliary
system will generate the same quantum trajectory ensemble as the canonical
one. Simple derivations further show that Eq. (30) is also equivalent to the
following conditions:

$\displaystyle H^{\prime}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}l_{\chi}(t)H(t)l_{\chi}^{-1}(t)+\frac{i\hbar}{2}\partial_{t}l_{\chi}(t)l_{\chi}^{-1}(t)+\frac{i\hbar}{4}l_{\chi}^{-1}(t)\left[\sum_{\omega}r(\omega)A^{\dagger}(\omega,t)A(\omega,t)\right]l_{\chi}(t)+H.c.,$
(33)

where $H.c.$ denotes Hermitian conjugation and $l_{\chi}(t)$ must satisfy the
operator equation given by

$\displaystyle\partial_{t}l_{\chi}^{2}(t)+{\cal
L}_{\chi}^{*}(t)[l_{\chi}^{2}(t)]=\Lambda_{T}(\chi)l_{\chi}^{2}(t),$ (34)

where

$\displaystyle{\cal
L}_{\chi}^{*}(t)(O)=\frac{i}{\hbar}[H(t),O]+\sum_{\omega}r(\omega)\left[e^{\chi\hbar\omega}A^{\dagger}(\omega,t)OA(\omega,t)-\frac{1}{2}\left\\{A^{\dagger}(\omega,t)A(\omega,t),O\right\\}\right].$
(35)

The $\chi$-dependence of the operator $l_{\chi}(t)$ is explicit. Using Eqs.
(33) and (27), we can rewrite the generator of the quantum master equation of
the auxiliary quantum system as

$\displaystyle{\cal L}^{\prime}(t)(O)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}[H^{\prime}(t),O]+\sum_{\omega}r(\omega)\left[A^{\prime}(\omega,t)OA^{\prime{\dagger}}(\omega,t)-\frac{1}{2}\left\\{A^{\prime{\dagger}}(\omega,t)A^{\prime}(\omega,t),O\right\\}\right]$
(36) $\displaystyle=$ $\displaystyle l_{\chi}(t){\cal
L}_{\chi}(t)\left[l_{\chi}^{-1}(t)Ol_{\chi}^{-1}(t)\right]l_{\chi}(t)-Ol_{\chi}^{-1}(t){\cal
L}^{*}_{\chi}(t)\left[l_{\chi}^{2}(t)\right]l_{\chi}^{-1}(t)-O(\partial_{t}l_{\chi}(t))l_{\chi}^{-1}(t)+(\partial_{t}l_{\chi}(t))l_{\chi}^{-1}(t)O$
$\displaystyle=$ $\displaystyle l_{\chi}(t){\cal
L}_{\chi}(t)\left[l_{\chi}^{-1}(t)Ol_{\chi}^{-1}(t)\right]l_{\chi}(t)-\Lambda_{T}(\chi)O+Ol_{\chi}^{-1}(t)(\partial_{t}l_{\chi}(t))+(\partial_{t}l_{\chi}(t))l_{\chi}^{-1}(t)O,$

where ${\cal L}_{\chi}(t)$ is the dual of ${\cal L}^{*}_{\chi}(t)$ and is
defined as ${\rm Tr}[O_{1}{\cal L}_{\chi}(t)(O_{2})]={\rm Tr}[{\cal
L}^{*}_{\chi}(t)(O_{1})O_{2}]$, and the third equation is a consequence of
applying Eq. (34). This is a time-dependent quantum Doob transform. We note
that the generator of the quantum master equation is slightly distinct from
that obtained by Carollo et al. Carollo  _et al._ (2018), in which the second
term about the scaled generating function $\Lambda_{T}(\chi)$ is absent.

###  IV.2 Long time limit case

Although the previous arguments are generally valid and the equations are
formally correct, it is not very clear whether Eq. (34) always achieves a
solution that can satisfy both the initial and the terminal conditions of Eq.
(31). In addition, these results do not depend on any special properties of
Floquet open quantum systems. Now we focus on the long time limit case. An
interesting feature of these open quantum systems is that, at a long time
limit, Eq. (2) evolves into a periodic limit cycle
$\rho(T)=\rho(T+2\pi/\Omega)$ Breuer  _et al._ (2000); Szczygielski  _et al._
(2013); Grifoni and Hänggi (1998) and

$\displaystyle\lim_{T\rightarrow\infty}\rho(T)=\sum_{n=1}^{N}P_{n}|u_{n}(T)\rangle\langle
u_{n}(T)|,$ (37)

where $P_{n}$, $n=1,\cdots,N$, is the stationary probability distribution,
which is uniquely determined by the generator of Eq. (2). Importantly, under
the same situation, the stochastic heat (16) obeys the large deviation
principle Touchette (2008); Gasparinetti  _et al._ (2014); Cuetara  _et al._
(2015); Liu and Su (2020): $p(j)\asymp\exp[-TI(j)]$, where the symbol $\asymp$
denotes an asymptotically exponential approximation and $I(j)$ is the rate
function. The scaled cumulant generating function is

$\displaystyle\lim_{T\rightarrow\infty}\Lambda_{T}(\chi)=\Lambda(\chi),$ (38)

where the time-independent $\Lambda(\chi)$ is the maximum eigenvalue of the
generator ${\cal L}_{\chi}(t)$. The rate function and scaled cumulant
generating function are related by the Legendre transform,

$\displaystyle I(j)=\max_{\chi}\\{j\chi-\Lambda(\chi)\\}.$ (39)

Appendix A explains Eq. (37) and the large deviation principle. Hence, we
introduce the logarithmic equivalence between the canonical quantum trajectory
ensemble and the trajectory ensemble generated by the auxiliary quantum system
as that in classical stochastic systems Garrahan  _et al._ (2009); Jack and
Sollich (2010); Chetrite and Touchette (2013, 2014):

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\ln\frac{p^{\prime}_{m,n}(\overrightarrow{\omega}_{M},T)}{p_{m,n}({\chi},\overrightarrow{\omega}_{M},T)}\rightarrow
0.$ (40)

Correspondingly, the conditions (31) are relaxed to

$\displaystyle l_{\chi}(0)|u_{m}(0)\rangle\propto|u_{m}(0)\rangle\hskip
5.69046pt{\rm and}\hskip
5.69046ptl_{\chi}(T)|u_{n}(T)\rangle\propto|u_{n}(T)\rangle.$ (41)

Then, the relationship between the moment generating functions of the
auxiliary and original quantum systems can be intuitively obtained as

$\displaystyle{\Phi}^{\prime}_{T}(\chi^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{m,n=1}^{N}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{M=0}^{\infty}{\cal
D}t\hskip
2.84544pte^{\chi^{\prime}Q_{T}(\overrightarrow{\omega})}p^{\prime}_{m,n}(\overrightarrow{\omega}_{M},T)$
(50) $\displaystyle=$
$\displaystyle\frac{1}{\Phi_{T}(\chi)}\sum_{m,n=1}^{N}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{M=0}^{\infty}{\cal
D}t\hskip
2.84544pte^{(\chi^{\prime}+\chi)Q_{T}(\overrightarrow{\omega})}\frac{p^{\prime}_{m,n}(\overrightarrow{\omega}_{M},T)}{p_{m,n}(\chi,\overrightarrow{\omega}_{M},T)}p_{m,n}(\overrightarrow{\omega}_{M},T).$
(59)

When the duration $T$ tends to infinity and we let the scaled cumulant
generating function of the auxiliary quantum system be
$\Lambda^{\prime}(\chi^{\prime})$, because of the logarithmic equivalence
(40), we immediately have

$\displaystyle{\Lambda}^{\prime}(\chi^{\prime})=\Lambda(\chi^{\prime}+\chi)-\Lambda(\chi).$
(60)

This agrees with the crucial result obtained by Carolla et al. for the time-
independent quantum master equations Carollo  _et al._ (2018) (see Eq. (10)
therein). Due to the Legendre transform (39) and Eq. (60), the typical heat
current of the auxiliary quantum system
$\partial_{\chi^{\prime}}\Lambda^{\prime}|_{\chi^{\prime}=0}$, which is equal
to the average heat current $\langle j^{\prime}\rangle$, is equal to
$\partial_{\chi}\Lambda(\chi)$. Note that the latter is also the average heat
current of the canonical quantum trajectory ensemble. Because of the
equivalence of the canonical and microcanonical trajectory ensembles, if
$\chi\neq 0$, the typical heat current of the auxiliary quantum system becomes
the atypical heat current of the original quantum system conditioned on this
large deviation Garrahan  _et al._ (2009); Jack and Sollich (2010); Chetrite
and Touchette (2013).

Eqs. (40)-(60) depend on the operator $l_{\chi}(t)$, which satisfies both Eqs.
(34) and (41). Here, we argue the existence of this operator. We conjecture
that the solution of Eq. (34) is diagonal in the Floquet basis, and its
elements are time-independent:

$\displaystyle
l_{\chi}^{2}(t)=\sum_{n=1}^{N}l^{2}_{n}(\chi)|u_{n}(t)\rangle\langle
u_{n}(t)|.$ (61)

Apparently, this solution satisfies Eq. (41). In addition, the operator is
periodic. Substituting the solution into Eq. (34) and writing the equation in
the Floquet basis, we obtain a matrix equation:

$\displaystyle\widetilde{\bf
R}(\chi)\overrightarrow{l^{2}}(\chi)=\Lambda(\chi)\overrightarrow{l^{2}}(\chi),$
(62)

where the column vector $\overrightarrow{l^{2}}(\chi)$ is the transpose of the
row vector $(l^{2}_{1}(\chi),\cdots,l^{2}_{N}(\chi))$ and the tilde symbol
represents the transpose; the definition of the matrix ${\bf R}(\chi)$ is
presented in Appendix A. Eq. (62) is simply the left eigenvector
$\overrightarrow{l^{2}}(\chi)$ of the matrix ${\bf R}(\chi)$ with the maximum
eigenvalue $\Lambda(\chi)$. Their positivity is ensured by the definition of
the moment generating function and the Perron-Frobenius theorem. Appendix A
gives some details.

The diagonal structure of $l_{\chi}^{2}(t)$ leads to an interesting
consequence. First, the formally complicated Hamiltonian $H^{\prime}(t)$ of
the auxiliary quantum system is simply equal to the Hamiltonian $H(t)$ of the
original quantum system. The proof includes two steps. In the first step, Eq.
(61) implies that the sum of the first two terms and their Hermitian
conjugations in Eq. (33) is

$\displaystyle\frac{1}{2}l_{\chi}(t)H(t)l_{\chi}^{-1}(t)+\frac{i\hbar}{2}\partial_{t}l_{\chi}(t)l_{\chi}^{-1}(t)+H.c$
(63) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n=1}^{N}(i\hbar\partial_{t}|u_{n}(t)\rangle+\epsilon_{n}|u_{n}(t)\rangle)\langle
u_{n}(t)|+H.c$ $\displaystyle=$ $\displaystyle H(t).$

For this purpose, we use Eq. (1) and the completeness of the Floquet basis.
The second step makes use of the microscopic expressions of the Lindblad
operators. Assume the interaction Hamiltonian between the quantum system and
the heat bath to be $A\otimes B$, where $A$ and $B$ are the Hermitian
operators of the system and heat bath, respectively. Then, the Lindblad
operator is Breuer and Petruccione (1997)

$\displaystyle A(\omega,t)$ $\displaystyle=$
$\displaystyle\sum_{m,n,q}\delta_{\omega,\epsilon_{n}-\epsilon_{m}+q\Omega}\langle\langle
u_{m}|A|u_{n}\rangle\rangle_{q}|u_{m}(t)\rangle\langle u_{n}(t)|e^{-iq\Omega
t},$ (64)

where $\delta$ is the Kronecker symbol and the time-independent coefficient
$\langle\langle u_{m}|A|u_{n}\rangle\rangle_{q}$ is the $q$-th harmonic of the
transition amplitude $\langle u_{m}(t)|A|u_{n}(t)\rangle$, i.e.,

$\displaystyle\langle\langle
u_{m}|A|u_{n}\rangle\rangle_{q}=\frac{\Omega}{2\pi}\int_{0}^{2\pi/\Omega}\langle
u_{m}(t)|A|u_{n}(t)\rangle e^{iq\Omega t}.$ (65)

Using this expression and performing a straightforward argument, we find that
the square bracket term in Eq. (33) is diagonal and especially real in the
Floquet basis. Because the operator $l_{\chi}(t)$ is also diagonal, we
conclude that the sum of the third term and its Hermitian conjugation is
exactly zero. In addition, we can also verify that the Lindblad operators
$A^{\prime}(\omega,t)$ and $A^{\prime{\dagger}}(\omega,t)$ of the auxiliary
quantum system are still the eigenoperators of the Floquet Hamiltonian ${\cal
H}(t)$.

Before closing this section, we want to point out that the scaled cumulant
generating function $\Lambda^{\prime}(\chi^{\prime})$ of the auxiliary quantum
system satisfies a matrix equation analogous to Eq. (62):

$\displaystyle\widetilde{\bf
R^{\prime}}(\chi^{\prime})\overrightarrow{l^{\prime
2}}(\chi^{\prime})=\Lambda^{\prime}(\chi^{\prime})\overrightarrow{l^{\prime
2}}(\chi^{\prime}),$ (66)

where the matrix $\widetilde{\bf R^{\prime}}(\chi^{\prime})$ is determined by
the generator ${\cal L}^{\prime}_{\chi^{\prime}}(t)$. Although we have
obtained Eq. (60), it shall be interesting to see how the same conclusion is
achieved and what the relation is between $\overrightarrow{l^{\prime
2}}(\chi^{\prime})$ and $\overrightarrow{l^{2}}(\chi)$ from an eigenvalue
matrix equation perspective. Some discussion of this topic is presented in
Appendix B.

##  V Two-level quantum system

In this section, we use a two-level quantum system driven by a periodic
external field Breuer and Petruccione (1997); Szczygielski  _et al._ (2013);
Langemeyer and Holthaus (2014); Gasparinetti  _et al._ (2014); Cuetara  _et
al._ (2015); Liu and Su (2020) to illustrate the construction of the auxiliary
quantum system at a long time limit. The Hamiltonian is given as

$\displaystyle
H(t)=\frac{\hbar{\omega_{0}}}{2}\sigma_{z}+\frac{\hbar\Omega_{R}}{2}\left(\sigma_{+}e^{-i\Omega
t}+\sigma_{-}e^{i\Omega t}\right),$ (67)

where $\omega_{0}$ is the transition frequency of the bare system,
$\Omega_{R}$ is the Rabi frequency, and $\Omega$ is the frequency of the
external field. The Floquet basis and quasi-energy are

$\displaystyle|u_{\pm}(t)\rangle=\frac{1}{\sqrt{2\Omega^{\prime}}}\left(\begin{array}[]{c}\pm\sqrt{\Omega^{\prime}\pm\delta}\\\
e^{i\Omega t}\sqrt{\Omega^{\prime}\mp\delta},\end{array}\right),$ (70)

and $\epsilon_{\pm}=\hbar(\Omega\pm\Omega^{\prime})/2$, respectively, where
$\Omega^{\prime}=\sqrt{\delta^{2}+\Omega_{R}^{2}}$ and the detuning parameter
$\delta=\omega_{0}-\Omega$. Here, we additionally set
$\Omega>\Omega^{\prime}$. We assume that the coupling between the quantum
system and heat bath is $\sigma_{x}$-coupling. There are six Lindblad
operators: three of them with Bohr frequencies $\Omega$,
$(\Omega-\Omega^{\prime})$, and $(\Omega+\Omega^{\prime})$ are

$\displaystyle A(\Omega,t)$ $\displaystyle=$
$\displaystyle\frac{\Omega_{R}}{2\Omega^{\prime}}\left(|u_{+}(t)\rangle\langle
u_{+}(t)|-|u_{-}(t)\rangle\langle u_{-}(t)|\right)e^{-i\Omega t},$
$\displaystyle A(\Omega-\Omega^{\prime},t)$ $\displaystyle=$
$\displaystyle\left(\frac{\delta-\Omega^{\prime}}{2\Omega^{\prime}}\right)|u_{+}(t)\rangle\langle
u_{-}(t)|e^{-i\Omega t},$ (71) $\displaystyle A(\Omega+\Omega^{\prime},t)$
$\displaystyle=$
$\displaystyle\left(\frac{\delta+\Omega^{\prime}}{2\Omega^{\prime}}\right)|u_{-}(t)\rangle\langle
u_{+}(t)|e^{-i\Omega t};$

the other three Lindblad operators with Bohr frequencies $-\Omega$,
$-(\Omega-\Omega^{\prime})$, and $-(\Omega+\Omega^{\prime})$ are the adjoint
operators of Eq. (V). Through a simple derivation, we obtain the matrix
elements of ${\textbf{R}}(\chi)$:

$\displaystyle{R}_{11}(\chi)$ $\displaystyle=$
$\displaystyle(e^{\chi\hbar\Omega}-1)\Gamma_{+\Omega}+(e^{-\chi\hbar\Omega}-1)\Gamma_{-\Omega}-\Gamma_{-(\Omega-\Omega^{\prime})}-\Gamma_{+(\Omega+\Omega^{\prime})},$
$\displaystyle{R}_{22}(\chi)$ $\displaystyle=$
$\displaystyle(e^{\chi\hbar\Omega}-1)\Gamma_{+\Omega}+(e^{-\chi\hbar\Omega}-1)\Gamma_{-\Omega}-\Gamma_{+(\Omega-\Omega^{\prime})}-\Gamma_{-(\Omega+\Omega^{\prime})},$
$\displaystyle{R}_{12}(\chi)$ $\displaystyle=$ $\displaystyle
e^{\chi\hbar(\Omega-\Omega^{\prime})}\Gamma_{+(\Omega-\Omega^{\prime})}+e^{-\chi\hbar(\Omega+\Omega^{\prime})}\Gamma_{-(\Omega+\Omega^{\prime})},$
$\displaystyle{R}_{21}(\chi)$ $\displaystyle=$ $\displaystyle
e^{-\chi\hbar(\Omega-\Omega^{\prime})}\Gamma_{-(\Omega-\Omega^{\prime})}+e^{\chi\hbar(\Omega+\Omega^{\prime})}\Gamma_{+(\Omega+\Omega^{\prime})},$
(72)

where the coefficients are

$\displaystyle\Gamma_{\pm\Omega}$ $\displaystyle=$
$\displaystyle\left(\frac{\Omega_{R}}{2\Omega^{\prime}}\right)^{2}r(\pm\Omega),$
$\displaystyle\Gamma_{\pm(\Omega-\Omega^{\prime})}$ $\displaystyle=$
$\displaystyle\left(\frac{\delta-\Omega^{\prime}}{2\Omega^{\prime}}\right)^{2}r(\pm(\Omega-\Omega^{\prime})),$
(73) $\displaystyle\Gamma_{\pm(\Omega+\Omega^{\prime})}$ $\displaystyle=$
$\displaystyle\left(\frac{\delta+\Omega^{\prime}}{2\Omega^{\prime}}\right)^{2}r(\pm(\Omega+\Omega^{\prime})).$

Because the above is a simple $2\times 2$ matrix, we can easily write its
maximum eigenvalue and the corresponding left eigenvector:
$\Lambda(\chi)=[{R}_{11}(\chi)+{R}_{22}(\chi)+B(\chi)]/2$ and

$\displaystyle(l^{2}_{1}(\chi),l^{2}_{2}(\chi))=\left(\frac{{R}_{11}(\chi)-{R}_{22}(\chi)+B(\chi)}{2{R}_{12}(\chi)},1\right),$
(74)

where
$B(\chi)=\sqrt{[{R}_{11}(\chi)-{R}_{22}(\chi)]^{2}+4{R}_{21}(\chi){R}_{12}(\chi)}$.
Based on these results, we can directly construct the auxiliary open quantum
system. For instance, three of the Lindblad operators are

$\displaystyle A^{\prime}(\Omega,t)$ $\displaystyle=$ $\displaystyle
e^{\chi\hbar\Omega/2}A(\Omega,t),$ (75) $\displaystyle
A^{\prime}(\Omega-\Omega^{\prime},t)$ $\displaystyle=$ $\displaystyle
e^{\chi\hbar(\Omega-\Omega^{\prime})/2}\frac{l_{1}(\chi)}{l_{2}(\chi)}A(\Omega-\Omega^{\prime},t),$
(76) $\displaystyle A^{\prime}(\Omega+\Omega^{\prime},t)$ $\displaystyle=$
$\displaystyle
e^{\chi\hbar(\Omega+\Omega^{\prime})/2}\frac{l_{2}(\chi)}{l_{1}(\chi)}A(\Omega+\Omega^{\prime},t).$
(77)

We must emphasize that the other three Lindblad operators
$A^{\prime}(\omega,t)$ with $\omega=-\Omega$, $-(\Omega-\Omega^{\prime})$, and
$-(\Omega+\Omega^{\prime})$ are equal to the adjoints of Eq. (75)- (77) and
that $\chi$ therein is replaced by $-\chi$. For the long time limit case, the
Hamiltonian of the auxiliary quantum system is simply Eq. (67).

Fig. 1 shows the rate functions of the original quantum system and auxiliary
quantum system. These functions are solved by applying the Legendre transform
(39) to the scaled cumulant generating functions (the curves) and by
simulating the quantum trajectories (the symbols). We can clearly see that, by
modulating $\chi$, the typical current of the auxiliary system, i.e., the
value at the minimum of its rate function, is located at different large
deviations from the original system. To intuitively express the behaviors of
the auxiliary quantum system, we also plot several segments of the quantum
trajectories from the original and auxiliary quantum systems in Fig. 2. Figure
(1) reminds us that, by choosing a sufficiently negative $\chi$, the typical
current can become negative. Then, along the quantum trajectories, one would
observe that heat is more frequently absorbed from than released into the heat
bath. Figure (2)(a) confirms this expectation. When we examine the algorithm
of simulating the quantum trajectories Liu and Xi (2016), we note that a
negative $\chi$ value exponentially inhibits the rates of quantum jumps with
heat release and exponentially increases the rates of quantum jumps with heat
absorption (data not shown). These $\chi$-modulated rates do not satisfy the
KMS condition. This observation simply explains the relation between panels
(a) and (b). If $\chi$ is positive, the opposite will happen; see Fig. (2)(c).

Let us comment on the typical negative heat current to finish this work.
According to the first law of thermodynamics, a positive power is output from
the heat bath. Therefore, the auxiliary quantum system seems to violate the
second law of thermodynamics. However, previous studies have confirmed that
Floquet open quantum systems strictly obey the law of thermodynamics
Szczygielski  _et al._ (2013); Cuetara  _et al._ (2015); Gasparinetti  _et
al._ (2014); Langemeyer and Holthaus (2014); Liu and Su (2020). In fact, there
is no contradiction because the auxiliary quantum system is not physical in
the sense that ${A^{\prime}}^{\dagger}(\omega,t)\neq A^{\prime}(-\omega,t)$;
see their definitions in Eq. (27).

Figure 1: Large deviation functions $I(j)$ of the original (solid curve) and
auxiliary two-level quantum systems with $\chi=-0.4$ (dashed curve) and $0.4$
(dotted curve). The functions are solved by applying the Legendre transform
(39) to Eq. (60). The symbols are the approximated large deviation functions
and are calculated by $-\ln p(j)/T$, where the probability distribution $p(j)$
is collected by simulating the quantum trajectories. The durations for the
cross and hollow symbols are 100 and 1000, respectively. The Fourier
transforms of the correlation functions are set to be $r(\omega)$$=$${\cal
A}|{\omega}|^{3}{\cal N}_{k}(\omega)$ for $\omega$$<$$0$; otherwise,
$r(\omega)$$=$${\cal A}|{\omega}|^{3}[{\cal N}_{k}(\omega)+1]$, where ${\cal
N}(\omega)$$=$$1/[\exp(\beta|\omega|)-1]$ and the coefficient $A$ is related
to the coupling strength between the system and the heat bath Breuer  _et al._
(2000). The parameters used are $\omega_{0}$$=$$1$, $\Omega_{R}$$=$$0.8$,
$\Omega$$=$$1.1$, ${\cal A}$$=$$1$, and $\beta=1/3$.  Figure 2: Segments of
the simulated quantum trajectories, where the $\chi$ values in panels (a),
(b), and (c) are $1$, $0.$, and $-1$, respectively, and the time intervals are
50 after the system is in a steady state. The other parameters are the same as
those used in Fig. (1). The arrow in the middle panel denotes the time
orientation. The short lines denote time points at which discrete quantum
jumps occur. The upper and bottom lines represent the released and absorbed
heat, respectively.

##  VI Conclusion

In this paper, we investigate the construction of an auxiliary open quantum
system that can generate a canonical quantum trajectory ensemble. The dynamics
of the original quantum system is governed by the Floquet quantum master
equation. The probability formula of the quantum trajectory is explicitly
used. We see that the theory on quantum trajectory ensembles is highly
analogous to the theory on classical trajectory ensembles. Since most of the
studies in the literature are concerned about trajectory ensembles conditioned
on simple time-integrated observables, it shall be interesting to explore
complex quantities, e.g., stochastic efficiency.

###### Acknowledgements.

This work was supported by the National Science Foundation of China under
Grants No. 12075016 and No. 11575016.

##  Appendix A Explanation of Eqs (37) and (38)

The moment generating function (26) can be rewritten as
$\Phi_{T}(\chi)=\sum_{m,n=1}^{N}p_{m}\phi_{m|n}(\chi,T)$, where

$\displaystyle\phi_{m|n}(\chi,T)=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{M=0}^{\infty}{\cal
D}(t)\hskip 2.84544ptp_{m|n}({\chi},\overrightarrow{\omega}_{M},T)e^{\chi
Q_{T}(\overrightarrow{\omega})}=\left.\langle
u_{n}(t)|\hat{\rho}(t)|u_{n}(t)\rangle\right|_{t=T}.$ (86)

Obviously, $\phi_{m|n}(\chi,T)$ is also the moment generating function of the
stochastic heat current for the special quantum trajectory ensemble, in which
the initial and terminal Floquet bases are constrained at $|u_{m}(0)\rangle$
and $|u_{n}(T)\rangle$, respectively. In the second equation of Eq. (86), the
operator $\hat{\rho}(t)$ ($0\leq t\leq T$) satisfies

$\displaystyle\partial_{t}\hat{\rho}(t)$ $\displaystyle=$ $\displaystyle{\cal
L}_{\chi}(t)[\hat{\rho}(t)],$ (87)

and its initial condition is $|u_{m}(0)\rangle\langle u_{m}(0)|$. To derive
this result, we use the quantum trajectory probability formula of Eq. (12) and
the ensemble average of Eq. (15) Liu and Xi (2016); Liu (2018); Liu and Su
(2020). Eq. (87) has been named the modified quantum master equation Esposito
_et al._ (2009); Garrahan and Lesanovsky (2010); Gasparinetti  _et al._
(2014); Cuetara  _et al._ (2015) or tilted quantum master equation Carollo
_et al._ (2018), where another technique, counting field statistics, is used
in the derivation. If we express the evolution equation in the Floquet basis,
we have

$\displaystyle\frac{d}{dt}\overrightarrow{\phi}_{m}(\chi,t)=\textbf{R}(\chi)\overrightarrow{\phi}_{m}(\chi,t),$
(88)

where $\overrightarrow{\phi}_{m}$ is the transpose of vector
$(\phi_{m|1},\cdots,\phi_{m|N})$ and its initial condition is
$\phi_{m|n}(0)=\delta_{m,n}$, $n=1,\cdots,N$. The non-diagonal matrix element
of ${\bf R}(\chi)$ is

$\displaystyle{R}_{ij}(\chi)=\sum_{\omega\neq
q\Omega}e^{\chi\hbar\omega}r(\omega)|\langle
u_{i}(t)|A(\omega,t)|u_{j}(t)\rangle|^{2},$ (89)

($i\neq j$), and the diagonal elements are

$\displaystyle{R}_{ii}(\chi)=\sum_{\omega=q\Omega}e^{\chi\hbar\omega}r(\omega)|\langle
u_{i}(t)|A(\omega,t)|u_{i}(t)\rangle|^{2}-\sum_{\omega}r(\omega)\sum_{j^{\prime}}|\langle
u_{j^{\prime}}(t)|A(\omega,t)|u_{i}(t)\rangle|^{2}.$ (90)

$i,j,j^{\prime}=1,\cdots N$. We emphasize that ${\bf R}(\chi)$ is a constant
matrix. Then, we formally express the moment generating functions as

$\displaystyle\phi_{m|n}(\chi,T)=\left(e^{T{\bf R}(\chi)}\right)_{mn}.$ (91)

Because these moment generating functions are always positive, according to
the Perron-Frobenius theorem, we immediately arrive at the conclusion that the
matrix ${\bf R}(\chi)$ has a unique maximum positive eigenvalue with a
positive left eigenvector, i.e., Eq. (62). When the duration $T$ tends to
infinity, Eq. (38) is proven. There is a special case, $\chi=0$. Under this
situation, the maximum left eigenvector is a unit vector, $(1,\cdots,1)$, and
the maximum eigenvalue trivially equals $0$. Note that the corresponding right
eigenvector is the stationary probability distribution,
$(P_{1},\cdots,P_{N})$. Hence, Eq. (37) is proven as well.

##  Appendix B $\Lambda^{\prime}(\chi^{\prime})$ of the auxiliary open quantum
system

The moment generating function $\Phi_{T}^{\prime}(\chi^{\prime})$ of the
auxiliary open quantum system can be solved as that of the original quantum
system: $\Phi_{T}^{\prime}(\chi^{\prime})={\rm Tr}[\hat{\rho}^{\prime}(T)]$
and $\partial_{t}\hat{\rho}^{\prime}(t)={\cal
L^{\prime}}_{\chi^{\prime}}(t)[\hat{\rho}^{\prime}(t)]$, where the generator
of the tilted quantum master equation is

$\displaystyle{\cal L}^{\prime}_{\chi^{\prime}}(t)(O)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}[H(t),O]+\sum_{\omega}r(\omega)\left[e^{\chi^{\prime}\hbar\omega}A^{\prime}(\omega,t)OA^{\prime{\dagger}}(\omega,t)-\frac{1}{2}\left\\{A^{\prime{\dagger}}(\omega,t)A^{\prime}(\omega,t),O\right\\}\right].$
(92)

Then, we construct the matrix ${\textbf{R}}^{\prime}(\chi^{\prime})$ and its
elements as follows. The non-diagonal elements are

$\displaystyle{R}^{\prime}_{ij}(\chi^{\prime})$ $\displaystyle=$
$\displaystyle\frac{l_{i}^{2}(\chi)}{l_{j}^{2}(\chi)}\sum_{\omega\neq
q\Omega}e^{(\chi^{\prime}+\chi)\hbar\omega}r(\omega)|\langle
u_{i}(t)|A(\omega,t)|u_{j}(t)\rangle|^{2}$ (93) $\displaystyle=$
$\displaystyle\frac{l_{i}^{2}(\chi)}{l_{j}^{2}(\chi)}{R}_{ij}(\chi^{\prime}+\chi),$

($i\neq j$), and the diagonal elements are

$\displaystyle{R}^{\prime}_{ii}(\chi^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\omega=q\Omega}e^{(\chi^{\prime}+\chi)\hbar\omega}r(\omega)|\langle
u_{i}(t)|A(\omega,t)|u_{i}(t)\rangle|^{2}-\sum_{\omega}e^{\chi\hbar\omega}r(\omega)\sum_{j^{\prime}}\left(\frac{l_{j^{\prime}}}{l_{i}}\right)^{2}|\langle
u_{j^{\prime}}(t)|A(\omega,t)|u_{i}(t)\rangle|^{2},$ (94) $\displaystyle=$
$\displaystyle{R}_{ii}(\chi^{\prime}+\chi)-\Lambda(\chi).$

Eq. (62) is used here. We immediately find that the vector
$\overrightarrow{l^{\prime 2}}(\chi^{\prime})=(l^{\prime
2}_{1}(\chi^{\prime}),\cdots,l^{\prime 2}_{N}(\chi^{\prime}))$ with elements

$\displaystyle l^{\prime
2}_{n}(\chi^{\prime})=\frac{l_{n}^{2}(\chi^{\prime}+\chi)}{l_{n}^{2}(\chi)},$
(95)

$n=1,\cdots,N$, is the eigenvector of Eq. (66), while the eigenvalue is Eq.
(60).

We shall mention that $\Lambda^{\prime}(\chi^{\prime})$ can be alternatively
obtained by abstract operator calculations as Carollo et al. Carollo  _et al._
(2018) did. First, Eq. (66) has an equivalent operator expression:

$\displaystyle\partial_{t}l^{\prime 2}_{\chi^{\prime}}(t)+{\cal
L}^{\prime*}_{\chi^{\prime}}(t)[l^{\prime
2}_{\chi^{\prime}}(t)]=\Lambda^{\prime}(\chi^{\prime})l_{\chi^{\prime}}^{\prime
2}(t),$ (96)

where

$\displaystyle{\cal L^{\prime}}_{\chi^{\prime}}^{*}(t)(O)$ $\displaystyle=$
$\displaystyle\frac{i}{\hbar}[H(t),O]+\sum_{\omega}r(\omega)\left[e^{\chi^{\prime}\hbar\omega}A^{\prime{\dagger}}(\omega,t)OA^{\prime}(\omega,t)-\frac{1}{2}\left\\{A^{\prime{\dagger}}(\omega,t)A^{\prime}(\omega,t),O\right\\}\right].$
(97) $\displaystyle=$ $\displaystyle l^{-1}_{\chi}(t){\cal
L}^{*}_{\chi^{\prime}+\chi}(t)\left[l_{\chi}(t)Ol_{\chi}(t)\right]l^{-1}_{\chi}(t)-\Lambda(\chi)O+l_{\chi}^{-1}(t)(\partial_{t}l_{\chi}(t))O+O(\partial_{t}l_{\chi}(t))l_{\chi}^{-1}(t).$

Then, it is direct to verify that the operator

$\displaystyle l^{\prime 2}_{\chi^{\prime}}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}l^{\prime
2}_{n}(\chi^{\prime})|u_{n}(t)\rangle\langle
u_{n}(t)|=l^{-1}_{\chi}(t)l^{2}_{\chi^{\prime}+\chi}(t)l^{-1}_{\chi}(t)$ (98)

satisfies Eq. (96) and that $\Lambda^{\prime}(\chi^{\prime})$ is simply Eq.
(60).

## References

  * Evans (2004a) R. M. L. Evans, Phys. Rev. Lett. 92, 150601 (2004a).
  * Garrahan  _et al._ (2007) J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, Phys. Rev. Lett. 98, 195702 (2007).
  * Maes and Netočný (2008) C. Maes and K. Netočný, EPL (Europhysics Letters) 82, 30003 (2008).
  * Maes  _et al._ (2008) C. Maes, K. Netočný, and B. Wynants, Physica A: Statistical Mechanics and its Applications 387, 2675 (2008).
  * Garrahan  _et al._ (2009) J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, J. Phys. A: Math. Theor. 42, 075007 (2009).
  * Garrahan and Lesanovsky (2010) J. P. Garrahan and I. Lesanovsky, Phys. Rev. Lett. 104, 160601 (2010).
  * Jack and Sollich (2010) R. L. Jack and P. Sollich, Progr. Theor. Phys. Suppl. ment 184, 304 (2010).
  * Nemoto and Sasa (2011) T. Nemoto and S.-i. Sasa, Phys. Rev. E 84, 061113 (2011).
  * Chetrite and Touchette (2013) R. Chetrite and H. Touchette, Phys. Rev. Lett. 111, 120601 (2013).
  * Chetrite and Touchette (2014) R. Chetrite and H. Touchette, Annales Henri Poincaré 16, 2005 (2014).
  * Derrida and Sadhu (2019) B. Derrida and T. Sadhu, J. Stat. Phys. 176, 773 (2019).
  * Evans (2010) R. M. Evans, Contemporary Physics 51, 413 (2010).
  * Evans (2004b) R. M. L. Evans, J. Phys. A: Math. Gen 38, 293 (2004b).
  * Garrahan  _et al._ (2011) J. P. Garrahan, A. D. Armour, and I. Lesanovsky, Phys. Rev. E 84, 021115 (2011).
  * Carollo  _et al._ (2018) F. Carollo, J. P. Garrahan, I. Lesanovsky, and C. Pérez-Espigares, Phys. Rev. A 98 (2018), 10.1103/physreva.98.010103.
  * Lebowitz and Spohn (1999) J. L. Lebowitz and H. Spohn, J. Stat. Phys. 95, 333 (1999).
  * Kurchan (2000) J. Kurchan, arXiv preprint cond-mat/0007360 (2000).
  * Maes (1999) C. Maes, J. Stat. Phys. 95, 367 (1999).
  * Gallavotti and Cohen (1995) G. Gallavotti and E. G. D. Cohen, Phys. Rev. Lett. 74, 2694 (1995).
  * Jarzynski (1997) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  * Crooks (2000) G. E. Crooks, Phys. Rev.E 61, 2361 (2000).
  * Seifert (2005) U. Seifert, Phys. Rev. Lett. 95, 040602 (2005).
  * Touchette (2008) H. Touchette, Phys. Rep. 478, 1 (2008).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford university press, 2002).
  * Carmichael (1993) H. Carmichael,  _An open systems approach to Quantum Optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_ , Vol. 18 (Springer, 1993).
  * Wiseman and Milburn (2010) H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge University Press, 2010).
  * Liu (2018) F. Liu, Prog. Phys. 38, 1 (2018).
  * Kosloff (2013) R. Kosloff, Entropy 15, 2100 (2013).
  * Szczygielski  _et al._ (2013) K. Szczygielski, D. Gelbwaser-Klimovsky, and R. Alicki, Phys. Rev. E 87, 012120 (2013).
  * Gasparinetti  _et al._ (2014) S. Gasparinetti, P. Solinas, A. Braggio, and M. Sassetti, New J. Phys. 16, 115001 (2014).
  * Cuetara  _et al._ (2015) G. B. Cuetara, A. Engel, and M. Esposito, New J. Phys. 17, 055002 (2015).
  * Gelbwaser-Klimovsky  _et al._ (2015) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, in  _Advances In Atomic, Molecular, and Optical Physics_ (Elsevier, 2015) pp. 329–407.
  * Liu and Su (2020) F. Liu and S. Su, Phys. Rev. E 101, 062144 (2020).
  * Blümel  _et al._ (1991) R. Blümel, A.Buchleitner, R. Graham, L. Sirko, U. Smilansky, and H. Walther, Phys Rev A 44, 4521 (1991).
  * Kohler  _et al._ (1997) S. Kohler, T. Dittrich, and P. Hänggi, Phys. Rev. E 55, 300 (1997).
  * Breuer and Petruccione (1997) H.-P. Breuer and F. Petruccione, Phys. Rev. A 55, 3101 (1997).
  * Alicki  _et al._ (2006) R. Alicki, D. A. Lidar, and P. Zanardi, Phys. Rev. A 73, 052311 (2006).
  * Zeldovich (1967) Y. B. Zeldovich, Sov. Phys. JETP 24, 1006 (1967).
  * Shirley (1965) J. H. Shirley, Phys. Rev. 138, B979 (1965).
  * Grifoni and Hänggi (1998) M. Grifoni and P. Hänggi, Phys. Rep. 304, 229 (1998).
  * Breuer  _et al._ (2000) H. P. Breuer, W. Huber, and F. Petruccione, Phys. Rev. E 61, 4883 (2000).
  * Plenio and Knight (1998) M. Plenio and P. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Liu and Xi (2016) F. Liu and J. Xi, Phys. Rev. E 94, 062133 (2016).
  * Breuer (2003) H.-P. Breuer, Phys. Rev. A 68, 032105 (2003).
  * Roeck and Maes (2006) W. D. Roeck and C. Maes, Rev. Math. Phys. 18, 619 (2006).
  * Crooks (2008) G. E. Crooks, Phys. Rev. A 77, 034101 (2008).
  * Horowitz (2012) J. M. Horowitz, Phys. Rev. E 85, 031110 (2012).
  * Hekking and Pekola (2013) F. W. J. Hekking and J. P. Pekola, Phys. Rev. Lett. 111, 093602 (2013).
  * Manzano  _et al._ (2015) G. Manzano, J. M. Horowitz, and J. M. R. Parrondo, Phys. Rev. E 92, 032129 (2015).
  * Note (1) This result is easy to obtain if one considers the consistency between the classical jump master equation and the diagonal part of the quantum master equation (2) represented in the Floquet basis Garrahan  _et al._ (2009); Jack and Sollich (2010).
  * Langemeyer and Holthaus (2014) M. Langemeyer and M. Holthaus, Phys. Rev. E 89, 012101 (2014).
  * Esposito  _et al._ (2009) M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod. Phys. 81, 1665 (2009).


# qnm: A Python package for calculating Kerr quasinormal modes, separation
constants, and spherical-spheroidal mixing coefficients

Leo C. Stein

††margin:  DOI:
https://doi.org/10.21105/joss.01683https://doi.org10.21105/joss.01683 Software
• https://github.com/openjournals/joss-
reviews/issues/1683https://doi.orgReview •
https://github.com/duetosymmetry/qnmhttps://doi.orgRepository •
https://doi.org/10.5281/zenodo.3459790https://doi.orgArchive Submitted: 21
August 2019
Published: 01 October 2019 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC-BY).

## Background

Black holes can be characterized from far away by their spectroscopic
gravitational-wave “fingerprints,” in analogy to electromagnetic spectroscopy
of atoms, ions, and molecules. The idea of using the quasi-normal modes (QNMs)
of black holes (BHs) for gravitational-wave (GW) spectroscopy was first made
explicit by Detweiler (1980). QNMs of rotating Kerr BHs in general relativity
(GR) depend only on the mass and spin of the BH. Thus GWs containing QNMs can
be used to infer the remnant BH properties in a binary merger, or as a test of
GR by checking the consistency between the inspiral and ringdown portions of a
GW signal (Abbott and others 2016; Isi et al. 2019).

For a review of QNMs see Berti, Cardoso, and Starinets (2009). A Kerr BH’s
QNMs are the homogeneous (source-free) solutions to the Teukolsky equation
(Teukolsky 1973) subject to certain physical conditions. The Teukolsky
equation can apply to different physical fields based on their spin-weight
$s$; for gravitational perturbations, we are interested in $s=-2$ (describing
the Newman-Penrose scalar $\psi_{4}$). The physical conditions for a QNM are
quasi-periodicity in time, of the form $\propto e^{-i\omega t}$ with complex
$\omega$; conditions of regularity, and that the solution has waves that are
only going down the horizon and out at spatial infinity. Separating the
radial/angular Teukolsky equations and imposing these conditions gives an
eigenvalue problem where the frequency $\omega$ and separation constant $A$
must be found simultaneously. This eigenvalue problem has a countably
infinite, discrete spectrum labeled by angular harmonic numbers $(\ell,m)$
with $\ell\geq 2$ (or $\ell\geq|s|$ for fields of other spin weight),
$-\ell\leq m\leq+\ell$, and overtone number $n\geq 0$.

There are several analytic techniques, e.g. one presented by Dolan and
Ottewill (2009), to approximate the desired complex frequency and separation
constant $(\omega_{\ell,m,n}(a),A_{\ell,m,n}(a))$ as a function of spin
parameter $0\leq a<M$ (we follow the convention of using units where the total
mass is $M=1$). These analytic techniques are useful as starting guesses
before applying the numerical method of Leaver (1985) for root-polishing.
Leaver’s method uses Frobenius expansions of the radial and angular Teukolsky
equations to find 3-term recurrence relations that must be satisfied at a
complex frequency $\omega$ and separation constant $A$. The recurrence
relations are made numerically stable to find so-called minimal solutions by
being turned into infinite continued fractions. In Leaver’s approach, there
are thus two “error” functions $E_{r}(\omega,A)$ and $E_{a}(\omega,A)$ (each
depending on $a,\ell,m,n$) which are given as infinite continued fractions,
and the goal is to find a pair of complex numbers $(\omega,A)$ which are
simultaneous roots of both functions. This is typically accomplished by
complex root-polishing, alternating between the radial and angular continued
fractions.

A refinement of this method was put forth by Cook and Zalutskiy (2014) (see
also Appendix A of Hughes (2000)). Instead of solving the angular Teukolsky
equation “from the endpoint” using Leaver’s approach, one can use a spectral
expansion with a good choice of basis functions. The solutions to the angular
problem are the _spin-weighted spheroidal harmonics_ , and the appropriate
spectral basis are the spin-weighted _spherical_ harmonics. This expansion is
written as (spheroidal on the left, sphericals on the right):

${}_{s}Y_{\ell
m}(\theta,\phi;a\omega)={\sum_{\ell^{\prime}=\ell_{\min}(s,m)}^{\ell_{\max}}}C_{\ell^{\prime}\ell
m}(a\omega)\ {}_{s}Y_{\ell^{\prime}m}(\theta,\phi)\,,$

where $\ell_{\min}=\max(|m|,|s|)$, and the coefficients $C_{\ell^{\prime}\ell
m}(a\omega)$ are called the spherical-spheroidal mixing coefficients (we
follow the conventions of Cook and Zalutskiy (2014), but compare Berti and
Klein (2014)). When recast in this spectral form, the angular equation becomes
very easy to solve via standard matrix eigenvector routines, see Cook and
Zalutskiy (2014) for details. If one picks values for $(s,\ell,m,a,\omega)$,
then the separation constant $A(a\omega)$ is returned as an eigenvalue, and a
vector of mixing coefficients $C_{\ell^{\prime}\ell m}(a\omega)$ are returned
as an eigenvector. From this new point of view there is now only one error
function to root-polish, $E_{r}(\omega)=E_{r}(\omega,A(\omega))$ where the
angular separation constant is found from the matrix method at any value of
$\omega$. Polishing roots of $E_{r}$ proceeds via any standard 2-dimensional
root-finding or optimization method.

The main advantage of the spectral approach is rapid convergence, and getting
the spherical-spheroidal mixing coefficients “for free” since they are found
in the process of solving the spectral angular eigenvalue problem.

## Summary

qnm is an open-source Python package for computing the Kerr QNM frequencies,
angular separation constants, and spherical-spheroidal mixing coefficients,
for given values of $(\ell,m,n)$ and spin $a$. There are several QNM codes
available, but some (London 2017) implement either analytic fitting formulae
(which only exist for a range of $s,\ell,m,n$) or interpolation from tabulated
data (so the user can not root-polish); others (Berti 2010) are in proprietary
languages such as Mathematica. We are not aware of any packages that provide
spherical-spheroidal mixing coefficients, which are necessary for multi-mode
ringdown GW modeling.

The qnm package includes a Leaver solver with the Cook-Zalutskiy spectral
approach to the angular sector, thus providing mixing coefficients. We also
include a caching mechanism to avoid repeating calculations. When the user
wants to solve at a new value of $a$, the cached data is used to interpolate a
good initial guess for root-polishing. We provide a large cache of low
$\ell,m,n$ modes so the user can start interpolating right away, and this
precomputed cache can be downloaded and installed with a single function call.
We have adapted the core algorithms so that numba (Lam, Pitrou, and Seibert
2015) can just-in-time compile them to optimized, machine-speed code. We rely
on numpy (Walt, Colbert, and Varoquaux 2011) for common operations such as
solving the angular eigenvalue problem, and we rely on scipy (Jones et al.
2001) for two-dimensional root-polishing, and interpolating from the cache
before root-polishing.

This package should enable researchers to perform ringdown modeling of
gravitational-wave data in Python, without having to interpolate into
precomputed tables or write their own Leaver solver. The author and
collaborators are already using this package for multiple active research
projects. By creating a self-documented, open-source code, we hope to
alleviate the high frequency of re-implemenation of Leaver’s method, and
instead focus efforts on making a single robust, fast, high-precision, and
easy-to-use code for the whole community. In the future, this code can be
extended to incorporate new features (like special handling of algebraically
special modes) or to apply to more general BH solutions (e.g., solving for
QNMs of Kerr-Newman or Kerr-de Sitter).

Development of qnm is hosted on
https://github.com/duetosymmetry/qnmhttps://doi.orgGitHub and distributed
through https://pypi.org/project/qnm/https://doi.orgPyPI; it can be installed
with the single command pip install qnm. Documentation is automatically built
on https://qnm.readthedocs.io/https://doi.orgRead the Docs, and can be
accessed interactively via Python docstrings. Automated testing is run on
https://travis-ci.org/duetosymmetry/qnmhttps://doi.orgTravis CI. The qnm
package is part of the https://bhptoolkit.org/https://doi.orgBlack Hole
Perturbation Theory Toolkit.

## Acknowledgements

We acknowledge E Berti and GB Cook for helpful correspondence and for making
testing data available. We further acknowledge M Giesler, I Hawke, Duncan
Macleod, L Magaña Zertuche, Matt Pitkin, and V Varma for
contributions/testing/feedback/suggestions.

## References

Abbott, B. P., and others. 2016. “Tests of general relativity with GW150914.”
_Phys. Rev. Lett._ 116 (22): 221101.
https://doi.org/10.1103/PhysRevLett.116.221101.

Berti, Emanuele. 2010. “Berti’s ringdown web page.”
https://pages.jh.edu/~eberti2/ringdown/.

Berti, Emanuele, Vitor Cardoso, and Andrei O. Starinets. 2009. “Quasinormal
modes of black holes and black branes.” _Class. Quant. Grav._ 26: 163001.
https://doi.org/10.1088/0264-9381/26/16/163001.

Berti, Emanuele, and Antoine Klein. 2014. “Mixing of spherical and spheroidal
modes in perturbed Kerr black holes.” _Phys. Rev._ D90 (6): 064012.
https://doi.org/10.1103/PhysRevD.90.064012.

Cook, Gregory B., and Maxim Zalutskiy. 2014. “Gravitational perturbations of
the Kerr geometry: High-accuracy study.” _Phys. Rev._ D90 (12): 124021.
https://doi.org/10.1103/PhysRevD.90.124021.

Detweiler, Steven L. 1980. “Black holes and gravitational waves. III. The
resonant frequencies of rotating holes.” _Astrophys. J._ 239: 292–95.
https://doi.org/10.1086/158109.

Dolan, Sam R., and Adrian C. Ottewill. 2009. “On an Expansion Method for Black
Hole Quasinormal Modes and Regge Poles.” _Class. Quant. Grav._ 26: 225003.
https://doi.org/10.1088/0264-9381/26/22/225003.

Hughes, Scott A. 2000. “The Evolution of circular, nonequatorial orbits of
Kerr black holes due to gravitational wave emission.” _Phys. Rev._ D61 (8):
084004. https://doi.org/10.1103/PhysRevD.61.084004.

Isi, Maximiliano, Matthew Giesler, Will M. Farr, Mark A. Scheel, and Saul A.
Teukolsky. 2019. “Testing the no-hair theorem with GW150914.” _Phys. Rev.
Lett._ 123 (11): 111102. https://doi.org/10.1103/PhysRevLett.123.111102.

Jones, Eric, Travis Oliphant, Pearu Peterson, and others. 2001. “SciPy: Open
Source Scientific Tools for Python.” http://www.scipy.org.

Lam, Siu Kwan, Antoine Pitrou, and Stanley Seibert. 2015. “Numba: A LLVM-based
Python JIT Compiler.” In _Proceedings of the Second Workshop on the Llvm
Compiler Infrastructure in Hpc_ , 7:1–7:6. LLVM ’15. New York, NY, USA: ACM.
https://doi.org/10.1145/2833157.2833162.

Leaver, E. W. 1985. “An Analytic representation for the quasi normal modes of
Kerr black holes.” _Proc. Roy. Soc. Lond._ A402: 285–98.
https://doi.org/10.1098/rspa.1985.0119.

London, Lionel. 2017. “London’s QNM codes @ 455f2d5.”
https://github.com/llondon6/kerr_public.

Teukolsky, Saul A. 1973. “Perturbations of a rotating black hole. 1.
Fundamental equations for gravitational electromagnetic and neutrino field
perturbations.” _Astrophys. J._ 185: 635–47. https://doi.org/10.1086/152444.

Walt, Stéfan van der, S. Chris Colbert, and Gaël Varoquaux. 2011. “The NumPy
Array: A Structure for Efficient Numerical Computation.” _Computing in Science
& Engineering_ 13 (2): 22–30. https://doi.org/10.1109/MCSE.2011.37.