# Measuring trajectories of environmental noise

Piotr Szańkowski<EMAIL_ADDRESS>Institute of Physics, Polish
Academy of Sciences, al. Lotników 32/46, PL 02-668 Warsaw, Poland

###### Abstract

In classical mechanics, a natural way to simplify a many-body problem is to
“replace” some of the elements of the composite system with surrogate force
fields. In the realm of quantum mechanics, however, such a description is
rarely compatible with the formalism of the theory. Nevertheless, the quantum
version of external field models—the so-called noise representations—can be
employed in certain circumstances. The mathematics behind these models
indicate that the appearing fields typically exhibit random fluctuations,
hence, the name noise field is more apt. In principle, measuring a classical
force field is a trivial task; all that is needed is a probe equipped with an
accelerometer that can be moved around while taking the measure of forces that
affect it. Unfortunately, an analogous method cannot work in the quantum case.
As indicated by the theory, the result of any measurement performed on the
quantum system affected by a noise field always appears as an average over
fluctuations, therefore, it is impossible to observe the noise fluctuations by
measuring the state of the probe. Here, we demonstrate that this limitation
can be circumvented and that it is possible to expose the noise field
fluctuations to direct observation. This allows one to sample the noise
trajectories (stochastic realizations), store them, and later use them to
simulate the dynamics of open quantum systems affected by the noise.

##  I Introduction

In recent years, we are witnessing a rapid pace of advancement in the field of
quantum technologies. This progress is spearheaded by a considerable effort
towards the design of nanoscale devices that employ genuine quantum effects to
perform certain tasks with efficiency that, at least theoretically, surpasses
the limits imposed by classical physics [1]; some examples include computing
[2], metrology [3, 4, 5], and communication [6, 7, 8]. Arguably, the most
significant obstacle standing between the theoretical concept of a quantum
device and its practical implementation is the phenomenon of decoherence that
results in the erosion of subtle, yet indispensable, quantum features. This
hindrance is unavoidable in virtually any realistic setting because of the
inevitable influence of the environment [9]. To understand, predict, and
hopefully counteract the destructive effects of decoherence, a systematic
description of the dynamics of a system open to the environment is required
[10]. An alluring branch of the theory of open systems is the noise
representation of the system–environment interaction, where the environment is
treated as a driving force and the system’s role is reduced to a passive
recipient of this driving. Formally, such representation is implemented by
replacing the full system–environment Hamiltonian with the system-only
operator where the environmental degrees of freedom are “replaced” by an
external noise—the so-called surrogate field111 Note that, traditionally, the
surrogate field is referred to in the literature as the classical noise. Here,
we adopt the stance advocated for in Ref. [11] to phase out such a
nomenclature. A number of examples shown in Ref. [11] demonstrate that the use
of the descriptor “classical” in this context is often inadequate and, in some
cases, can even be misleading; e.g., the observation that the existence of the
noise representation does not preclude the formation of system–environment
entanglement, is a particularly compelling argument.  $\xi$,

$\displaystyle\hat{H}_{SE}=\hat{H}_{S}\otimes\hat{\mathds{1}}+\hat{V}_{S}\otimes\hat{V}_{E}+\hat{\mathds{1}}\otimes\hat{H}_{E}$
$\displaystyle\phantom{H_{SE}=}\xrightarrow{\begin{subarray}{c}\text{noise}\\\
\text{representation}\end{subarray}}\hat{H}_{t}[\xi]=\hat{H}_{S}+\xi(t)\hat{V}_{S}.$
(1)

Here, $\hat{H}_{S(E)}$ are the free Hamiltonians of system (environment) and
$\hat{V}_{S(E)}$ are the system (environment) sides of the coupling. When
noise representation is valid, the system dynamics can be simulated with the
averaged unitary evolution driven by the surrogate field $\xi$,

$\displaystyle\hat{\rho}_{S}(t)$ $\displaystyle=$
$\displaystyle\mathrm{tr}_{E}\left(e^{-it\hat{H}_{SE}}\hat{\rho}_{S}\otimes\hat{\rho}_{E}\,e^{it\hat{H}_{SE}}\right)$
(2) $\displaystyle=$
$\displaystyle\overline{\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi]}.$

Here, $\mathrm{tr}_{E}(\ldots)$ is the partial trace over the $E$ subspace,
$\overline{(\ldots)}$ is the average over realizations (trajectories) of the
stochastic process $\xi$, and each individual realization is governed by
simple (relative to the full $SE$ problem) unitary evolution operator
conditioned by the trajectory $\xi(\tau)$,

$\hat{U}_{t}[\xi]=\mathcal{T}e^{-i\int_{0}^{t}\hat{H}_{\tau}[\xi]d\tau},$ (3)

where $\mathcal{T}$ indicates the time-ordering operation. Recent developments
[11] have provided an explanation of the physical origins of the noise
representation and have quantified the sufficient conditions for the
environment to facilitate the surrogate field that represents it. These
explanations have been built around formal analogy between stochastic average
and the partial trace operation. On the one hand, the average can be expressed
in the terms of functional integral

$\overline{\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi]}=\int\mathcal{D}\xi\,\mathcal{P}[\xi]\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi],$
(4)

where the non-negative functional $\mathcal{P}[\xi]$ is the probability
distribution for trajectories $\xi(t)$ (real-valued functions of time) of
stochastic process $\xi$. On the other hand, as demonstrated in [11], it is
always possible to rewrite the partial trace as

$\displaystyle\mathrm{tr}_{E}\left(e^{-it\hat{H}_{SE}}\hat{\rho}_{S}\otimes\hat{\rho}_{E}e^{it\hat{H}_{SE}}\right)$
$\displaystyle\phantom{\mathrm{tr}_{E}(}=\int\mathcal{D}\xi\int\mathcal{D}\zeta\,\mathcal{Q}[\xi,\zeta]\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}_{t}^{\dagger}[\zeta],$
(5)

where the functional $\mathcal{Q}[\xi,\zeta]$ is complex-valued and, hence, it
can be classified as a quasi-probability distribution for the real-valued
trajectories of the two-component “quantum process” $(\xi(\tau),\zeta(\tau))$.
By comparing (4) and (5), one concludes that the noise representation is
effected when the quasi-probability can be treated as a proper probability
distribution, i.e., when

$\mathcal{Q}[\xi,\zeta]=\delta(\zeta-\xi)\mathcal{P}[\xi]\geqslant 0.$ (6)

Note that $\mathcal{Q}$ is not directly related to the influence functional
used in the description of open-system dynamics within Feynman’s path integral
formalism [12, 13, 14]. In more technical terms, the paths in Feynman’s
formalism and the processes $\xi,\zeta$ in (5) are fundamentally different
entities: the former refer to dynamical variables associated with the
Lagrangians assigned to quantum systems $S$ and $E$, while the latter can be
interpreted as trajectories traced through the domain of the eigenvalues of
the $E$ side of the coupling $\hat{V}_{E}$. More importantly, the influence
functional contains the state $\hat{\rho}_{E}$, the Hamiltonian $\hat{H}_{E}$,
and the $SE$ interaction (thus, making it dependent on $S$), while
$\mathcal{Q}$ leaves out the interaction term (it is incorporated into
$\hat{U}_{t}$, instead); this is a crucial difference.

The key feature of the formulation (5) is that the properties of quasi-
probability $\mathcal{Q}$ are determined exclusively by the environment, i.e.,
$\mathcal{Q}$ is independent of the open system $S$ [11]. Therefore, the
validity of the noise representation (i.e., whether $\mathcal{Q}$ is a proper
probability distribution) is a feature of the environment, and the form of the
resultant surrogate field is established only by the dynamical laws governing
the environmental degrees of freedom (i.e., the triplet of operators
$\hat{H}_{E}$, $\hat{\rho}_{E}$, and $\hat{V}_{E}$). It also means that the
validity of noise representation is not contingent on the effects of
interaction between $E$ and any particular $S$, i.e., it does not matter
whether $S$ exerts back-action onto $E$, or if $S$ entangles with $E$, etc.
[11]. Consequently, when the noise representation exists for a given
environment, it is an objective representation—the dynamics of any system $S$
defined by the choice of $\hat{H}_{S}$, $\hat{\rho}_{S}$, and $\hat{V}_{S}$,
that is coupled to $E$ via a given environment-side operator $\hat{V}_{E}$, is
simulable with the same surrogate field that effectively replaces
$\hat{V}_{E}$. This is the main reason why noise representations are so useful
whenever they can be implemented: once one determines the objective surrogate
field generated by the environment, describing the dynamics of any system open
to this environment is significantly simplified.

From the point of view of practical implementation, the main challenge in
executing the noise representation simulation is the average operation. The
functional integral form (4), as well as the probability distribution
functional $\mathcal{P}[\xi]$, are merely formal theoretical tools that are
never suitable for actual calculations. The solution that is practical
exploits the equivalence between the average over the probability distribution
and the average over the ensemble of samples drawn from this distribution,

$\displaystyle\overline{\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi]}$
$\displaystyle=$
$\displaystyle\int\mathcal{D}\xi\,\mathcal{P}[\xi]\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi]$
(7) $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}\hat{U}_{t}[\xi^{(j)}]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi^{(j)}],$

where the trajectories $\xi^{(j)}(\tau)$ belong to the set
$\mathcal{E}_{N}=\\{\xi^{(j)}(\tau)\\}_{j=1}^{N}$—the ensemble of sample
trajectories of stochastic process $\xi$ drawn at random from the distribution
$\mathcal{P}[\xi]$. Of course, assembling the infinite number of samples is an
impossibility, and so, in practice, one resorts to the approximation

$\overline{\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}_{t}^{\dagger}[\xi]}\approx\frac{1}{N_{e}}\sum_{j=1}^{N_{e}}\hat{U}_{t}[\xi^{(j)}]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi^{(j)}],$
(8)

where $N_{e}$ is the number of samples in the ensemble $\mathcal{E}_{N_{e}}$
of finite size; the larger the number of samples, the more accurate the
approximation. Therefore, to find the average, instead of performing the
abstract functional integral, one takes the ensemble of sample trajectories
$\mathcal{E}_{N_{e}}$ and calculates the unitary evolution
$\hat{\rho}_{S}^{(j)}(t)=\hat{U}_{t}[\xi^{(j)}]\hat{\rho}_{S}\hat{U}^{\dagger}_{t}[\xi^{(j)}]$
for every trajectory $\xi^{(j)}(\tau)$ (e.g., by numerically integrating the
corresponding von Neumann equation). The results are then stored in the
ensemble of density matrices $\\{\hat{\rho}_{S}^{(j)}(t)\\}_{j=1}^{N_{e}}$.
The individual members of this ensemble are meaningless and their only purpose
is to calculate the physically meaningful average that approximates the system
state
$(1/N_{e})\sum_{j=1}^{N_{e}}\hat{\rho}_{S}^{(j)}(t)\approx\hat{\rho}_{S}(t)$.

Clearly, for this procedure to work, first, one needs to populate the ensemble
of trajectories. According to Ref. [11], in principle, the trajectories can be
drawn from so-called joint probability distribution
$p_{E}^{(k)}(\xi_{k}t_{k};\ldots;\xi_{2}t_{2};\xi_{1}t_{1})$ that describes
the probability for the surrogate field $\xi$ to have value $\xi_{1}$ at the
initial time $t_{1}$, followed by $\xi_{2}$ at $t_{2}>t_{1}$, followed by
$\xi_{3}$ at $t_{3}>t_{2}$, etc. The joint probability distributions form an
infinite family $\\{p_{E}^{(k)}\\}_{k=1}^{\infty}$ that fully defines the
stochastic process $\xi$, and it can be used as an alternative to the
description in terms of the distribution functional $\mathcal{P}[\xi]$. The
fundamental advantage one gets for using such a mode of description is that,
unlike the abstract functional, joint distributions are standard functions
given by closed analytical formulas [11],

$\displaystyle
p^{(k)}_{E}(\xi_{k}t_{k};\ldots;\xi_{1}t_{1})=\sum_{\begin{subarray}{c}n_{1}:\\\
\xi_{1}=v_{n_{1}}\end{subarray}}\cdots\sum_{\begin{subarray}{c}n_{k}:\\\
\xi_{k}=v_{n_{k}}\end{subarray}}\langle
n_{1}|\hat{\rho}_{E}(t_{1})|n_{1}\rangle$
$\displaystyle\phantom{p^{(k)}_{E}(\xi_{k}t_{k};\ldots;\xi_{1})}\times\prod_{l=1}^{k-1}|\langle
n_{l+1}|e^{-i(t_{l+1}-t_{l})\hat{H}_{E}}|n_{l}\rangle|^{2},$ (9)

where $\hat{\rho}_{E}(t)=e^{-it\hat{H}_{E}}\hat{\rho}_{E}\,e^{it\hat{H}_{E}}$
and $\\{v_{n}\\}_{n}$, $\\{|n\rangle\\}_{n}$ are (possible degenerate)
eigenvalues and eigenstates of the coupling operator $\hat{V}_{E}$ that is to
be replaced in the simulation by the surrogate field $\xi$. However,
calculating any one of those distributions is almost as difficult as solving
the original system–environment problem (one needs to know $\hat{H}_{E}$,
$\hat{V}_{E}$, and $\hat{\rho}_{E}$, and be able to diagonalize the
Hamiltonian). This is not entirely surprising; after all, the surrogate field,
when it exists, is an exact representation of the environmental degrees of
freedom. This difficulty could be circumvented if one was able to draw the
sample trajectories from some other, ideally more accessible, source. In Sec.
II, we demonstrate that the sequential measurements of the environmental
observable $\hat{V}_{E}$ can be used as such a source of the surrogate field
trajectories. Therefore, for the environments that facilitate noise
representation, one can setup the ensemble of trajectories—and thus enable
future simulations of open-system dynamics—by directly measuring the sample
trajectories instead of drawing them from a calculated probability
distribution.

##  II Trajectory sampling: the proof of concept

Here, we demonstrate that, assuming a valid noise representation, the result
of a sequence of projective measurements [15] of observable $\hat{V}_{E}$
performed directly on the environment $E$, can be considered as a sample
trajectory of the surrogate field $\xi$. To achieve this goal, we make
extensive use of some technical results obtained previously in Ref. [11], that
we cite without redundant derivation.

First, we define the spectral decomposition of the observable of interest,

$\displaystyle\hat{V}_{E}$ $\displaystyle=$
$\displaystyle\sum_{n}v_{n}|n\rangle\langle n|$ (10) $\displaystyle=$
$\displaystyle\sum_{\xi\in\Omega}\xi\sum_{n:v_{n}=\xi}|n\rangle\langle
n|\equiv\sum_{\xi\in\Omega}\xi\,\hat{\Pi}(\xi),$

where $\Omega$ is the set of all its unique eigenvalues (the spectrum). The
operator $\hat{V}_{E}$ above is the same as the environmental side of the
coupling in Eq. (1) that is effectively superseded by the surrogate field
$\xi$ when the noise representation is valid.

Next, in order to showcase how exactly the assumption of a valid noise
representation comes into play, we recall the joint quasi-probability
decomposition of the state of arbitrary system $S$ open to environment $E$.
Following the results of Ref. [11], the density matrix of $S$ in the
interaction picture,
$\hat{\rho}_{S}^{I}(t)=e^{it\hat{H}_{S}}\hat{\rho}_{S}(t)e^{-it\hat{H}_{S}}$,
can be written as

$\displaystyle\hat{\rho}^{I}_{S}(t)=e^{it\hat{H}_{S}}\mathrm{tr}_{E}\left(e^{-it\hat{H}_{SE}}\hat{\rho}_{S}\otimes\hat{\rho}_{E}\,e^{it\hat{H}_{SE}}\right)e^{-it\hat{H}_{S}}$
$\displaystyle=\sum_{k=0}^{\infty}(-i)^{k}\int_{0}^{t}\\!\\!dt_{k}\int_{0}^{t_{k}}\\!\\!dt_{k-1}\cdots\int_{0}^{t_{2}}\\!\\!dt_{1}\\!\sum_{\xi_{1},\zeta_{1}\in\Omega}\\!\cdots\\!\sum_{\xi_{k},\zeta_{k}\in\Omega}$
$\displaystyle\phantom{\sum}\times
q_{E}^{(k)}(\xi_{k}\zeta_{k}t_{k};\ldots;\xi_{1}\zeta_{1}t_{1})\left(\prod_{l=k}^{1}\mathscr{W}(\xi_{l}\zeta_{l}t_{l})\right)\hat{\rho}_{S},$
(11)

where the super operators acting on the initial state are given by the
$S$-side of the $SE$ coupling, $\mathscr{W}(\xi\zeta
t)\hat{A}=\xi\hat{V}_{S}(t)\hat{A}-\zeta\hat{A}\hat{V}_{S}(t)$ with
$\hat{V}_{S}(t)=e^{it\hat{H}_{S}}\hat{V}_{S}e^{-it\hat{H}_{S}}$, and the
symbol $\prod_{l=k}^{1}\mathscr{A}_{l}$ indicates the product where super-
operators are arranged in the decreasing order of the index,
$\mathscr{A}_{k}\mathscr{A}_{k-1}\cdots\mathscr{A}_{1}$. The members of the
infinite family of nonpositively defined functions
$\\{q_{E}^{(k)}\\}_{k=1}^{\infty}$—the joint quasi-probability
distributions—are determined exclusively by the environment and their explicit
form reads

$\displaystyle
q_{E}^{(k)}(\xi_{k}\zeta_{k}t_{k};\ldots;\xi_{1}\zeta_{1}t_{1})$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left(\prod_{l=k}^{1}e^{it_{l}\hat{H}_{E}}\hat{\Pi}(\xi_{l})e^{-it_{l}\hat{H}_{E}}\right)\hat{\rho}_{E}\left(\prod_{l^{\prime}=1}^{k}e^{it_{l^{\prime}}\hat{H}_{E}}\hat{\Pi}(\zeta_{l^{\prime}})e^{-it_{l^{\prime}}\hat{H}_{E}}\right)\right]$
(12) $\displaystyle=$
$\displaystyle\delta_{\xi_{k},\zeta_{k}}\mathrm{tr}\left[\hat{\Pi}(\xi_{k})\left(\prod_{l=k-1}^{1}e^{i(t_{l}-t_{l+1})\hat{H}_{E}}\hat{\Pi}(\xi_{l})\right)\hat{\rho}_{E}(t_{1})\left(\prod_{l^{\prime}=1}^{k-1}\hat{\Pi}(\zeta_{l^{\prime}})e^{i(t_{l^{\prime}+1}-t_{l^{\prime}})\hat{H}_{E}}\right)\right],$

where the symbol $\prod_{l=k}^{1}\hat{A}_{l}$
($\prod_{l^{\prime}=1}^{k}\hat{A}_{l^{\prime}}$) indicates the product where
operators are arranged in the decreasing (increasing) order of the index
$\hat{A}_{k}\hat{A}_{k-1}\cdots\hat{A}_{1}$ ($\hat{A}_{1}\cdots\hat{A}_{k}$).

The joint quasi-probability decomposition (11) is a practical alternative to
the functional-integral description (5); essentially,
$\\{q_{E}^{(k)}\\}_{k=1}^{\infty}$ is to $\mathcal{Q}[\xi,\zeta]$ as
$\\{p_{E}^{(k)}\\}_{k=1}^{\infty}$ is to $\mathcal{P}[\xi]$ [see Eq. (9)].
When translated into the language of joint distributions, the sufficient
condition for valid noise representation
$\mathcal{Q}[\xi,\zeta]\approx\delta(\zeta-\xi)\mathcal{P}[\xi]$ becomes

$\displaystyle
q_{E}^{(k)}(\xi_{k}\zeta_{k}t_{k};\ldots;\xi_{1}\zeta_{1}t_{1})\approx\left(\prod_{l=1}^{k}\delta_{\xi_{l},\zeta_{l}}\\!\right)\\!p_{E}^{(k)}(\xi_{k}t_{k};\ldots;\xi_{1}t_{1}),$
(13)

for all $k=1,2,\ldots$. When the above condition is satisfied, then

$\displaystyle\hat{\rho}^{I}_{S}(t)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}(-i)^{k}\int_{0}^{t}dt_{k}\cdots\int_{0}^{t_{2}}dt_{1}\sum_{\xi_{1}\ldots\xi_{k}\in\Omega}$
(14) $\displaystyle\phantom{\sum}\times
p_{E}^{(k)}(\xi_{k}t_{k};\ldots;\xi_{1}t_{1})\left(\prod_{l=k}^{1}\mathscr{W}(\xi_{l}\xi_{l}t_{l})\right)\hat{\rho}_{S}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}(-i)^{k}\int_{0}^{t}dt_{k}\cdots\int_{0}^{t_{2}}dt_{1}\,\overline{\xi(t_{k})\cdots\xi(t_{1})}$
$\displaystyle\phantom{\sum}\times[\hat{V}_{S}(t_{k}),\cdots[\hat{V}_{S}(t_{2}),[\hat{V}_{S}(t_{1}),\hat{\rho}_{S}]\;]\cdots\;]$
$\displaystyle=$ $\displaystyle
e^{it\hat{H}_{S}}\overline{\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}_{t}^{\dagger}[\xi]}e^{-it\hat{H}_{S}},$

where we have used the relation
$\mathscr{W}(\xi_{l}\xi_{l}t_{l})\hat{A}=\xi_{l}[\hat{V}_{S}(t_{l}),\hat{A}]$,
and the moment of stochastic process $\xi$ reads

$\displaystyle\overline{\xi(t_{k})\cdots\xi(t_{1})}$ $\displaystyle=$
$\displaystyle\sum_{\xi_{1},\ldots,\xi_{k}\in\Omega}\left(\prod_{l=1}^{k}\xi_{l}\right)p_{E}^{(k)}(\xi_{k}t_{k};\ldots;\xi_{1}t_{1}),$

hence the interpretation that $p_{E}^{(k)}$ is the joint distribution that
describes the probability of the trajectory passing through value $\xi_{1}$ at
$t_{1}$ followed by $\xi_{2}$ at $t_{2}$, etc.

Of course, whether the relation (13) holds depends only on the properties of
$E$ (hence, the objectivity of the noise representation), and thus, this
equality can be invoked in any context—its utility is not limited to the
quasi-probability decomposition of density matrices. We use this to our
advantage as we now switch to a different setup where, instead of open-system
dynamics, we examine the sequence of measurements performed directly on $E$.

The environment is initialized in the state $\hat{\rho}_{E}$ and its free
evolution is governed by the Hamiltonian $\hat{H}_{E}$. Assume we can perform
an ideal projective measurement of physical quantity described by the
observable $\hat{V}_{E}$ characterized by the spectral decomposition (10).
Then, according to the Born rule, the probability of obtaining a particular
eigenvalue $\xi_{t_{1}}$ as a result of the measurement carried out after
duration $t_{1}$ has passed is given by
$P(\xi_{t_{1}})=\mathrm{tr}[\hat{\Pi}(\xi_{t_{1}})\hat{\rho}_{E}(t_{1})]$,
which, when compared with Eq. (12), equals

$P(\xi_{t_{1}})=q^{(1)}_{E}(\xi_{t_{1}}\xi_{t_{1}}t_{1}).$ (16)

In accordance with the collapse postulate, assuming that the duration of the
measurement process is effectively instantaneous in comparison with the time
scale set by $\hat{H}_{E}$, the a posteriori state of the environment
continues to evolve as

$\displaystyle\hat{\rho}(t|\xi_{t_{1}})$ $\displaystyle=$
$\displaystyle\frac{e^{-i(t-t_{1})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{1}})\hat{\rho}_{E}(t_{1})\hat{\Pi}(\xi_{t_{1}})e^{i(t-t_{1})\hat{H}_{E}}}{q_{E}^{(1)}(\xi_{t_{1}}\xi_{t_{1}}t_{1})}.$

When the projective measurement is performed again on the same system after
certain period, then the (now conditional) probability of obtaining result
$\xi_{t_{2}}$ at time $t_{2}>t_{1}$ (provided that $\xi_{t_{1}}$ was measured
previously) is given by

$\displaystyle P(\xi_{t_{2}}|\xi_{t_{1}})\\!$ $\displaystyle=$
$\displaystyle\mathrm{tr}\\!\left[\hat{\Pi}(\xi_{t_{2}})\hat{\rho}(t_{2}|\xi_{t_{1}})\right]\\!=\\!\frac{\mathrm{tr}\\!\left[\hat{\Pi}(\xi_{t_{2}})e^{i(t_{1}-t_{2})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{1}})\hat{\rho}_{E}(t_{1})\hat{\Pi}(\xi_{t_{1}})e^{i(t_{2}-t_{1})\hat{H}_{E}}\right]}{q_{E}^{(1)}(\xi_{t_{1}}\xi_{t_{1}}t_{1})}\\!=\\!\frac{q^{(2)}_{E}(\xi_{t_{2}}\xi_{t_{2}}t_{2};\xi_{t_{1}}\xi_{t_{1}}t_{1})}{q_{E}^{(1)}(\xi_{t_{1}}\xi_{t_{1}}t_{1})}.$
(18)

Therefore, the probability of measuring the sequence of results
$(\xi_{t_{1}},\xi_{t_{2}})$ (i.e., $\xi_{t_{1}}$ at $t_{1}$ followed by
$\xi_{t_{2}}$ at $t_{2}$) reads

$\displaystyle P[(\xi_{t_{1}},\xi_{t_{2}})]$ $\displaystyle=$ $\displaystyle
P(\xi_{t_{2}}|\xi_{t_{1}})P(\xi_{t_{1}})=q_{E}^{(2)}(\xi_{t_{2}}\xi_{t_{2}}t_{2};\xi_{t_{1}}\xi_{t_{1}}t_{1}),$
(19)

and the a posteriori state after the latest measurement is

$\displaystyle\hat{\rho}(t|\xi_{t_{1}},\xi_{t_{2}})$ $\displaystyle=$
$\displaystyle\frac{1}{P(\xi_{t_{2}}|\xi_{t_{1}})}e^{-i(t-t_{2})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{2}})\hat{\rho}(t_{2}|\xi_{t_{1}})\hat{\Pi}(\xi_{t_{2}})e^{i(t-t_{2})\hat{H}_{E}}$
(20) $\displaystyle=$
$\displaystyle\frac{e^{-i(t-t_{2})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{2}})e^{i(t_{1}-t_{2})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{1}})\hat{\rho}_{E}(t_{1})\hat{\Pi}(\xi_{t_{1}})e^{i(t_{2}-t_{1})\hat{H}_{E}}\hat{\Pi}(\xi_{t_{2}})e^{i(t-t_{2})\hat{H}_{E}}}{q_{E}^{(2)}(\xi_{t_{2}}\xi_{t_{2}}t_{2};\xi_{t_{1}}\xi_{t_{1}}t_{1})}.$

At this point it is clear that when this measurement procedure is iterated
over consecutive time steps $t_{k}>\cdots>t_{2}>t_{1}$, the probability to
obtain the sequence of results $(\xi_{t_{1}},\xi_{t_{2}},\ldots,\xi_{t_{k}})$
will be given by

$P[(\xi_{t_{1}},\ldots,\xi_{t_{k}})]=q_{E}^{(k)}(\xi_{t_{k}}\xi_{t_{k}}t_{k};\ldots;\xi_{t_{1}}\xi_{t_{1}}t_{1}).$
(21)

That is, the probability of measuring the sequence of length $k$ equals the
“diagonal” part of the quasi-probability $q_{E}^{(k)}$. Even though, from a
technical point of view, the result (21) was straightforward to arrive at, it
is by no means trivial. We have found here a direct link between the dynamics
of a system open to $E$ and the process of sequential measurement performed on
$E$, without any involvement of open systems. This kind of connection, where
the consequences of quantum measurements are quantified with elements of
formalism normally employed in open-system theory, has been made previously,
albeit in a remote context. The quantum trajectory theory [16] in question,
describes the state of a (closed) system that is being continuously monitored,
i.e., the system that undergoes the sequence of measurements performed with
frequency approaching infinity—or, equivalently, infinitesimal delay between
consecutive measurements. It is shown there that when the measurement results
are ignored by the observer, it is possible to parametrize the time evolution
of the a posteriori state of the monitored system with a dynamical equation
that is formally identical to the Lindblad form of the master equation that is
regularly used to describe open-system dynamics.

The connection with surrogate field trajectories is established when we invoke
the noise representation validity (13),

$P[(\xi_{t_{1}},\ldots,\xi_{t_{k}})]\approx
p_{E}^{(k)}(\xi_{t_{k}}t_{k};\ldots;\xi_{t_{1}}t_{1}),$ (22)

so that probability distribution of the result sequence
$(\xi_{t_{1}},\ldots,\xi_{t_{k}})$ is the same as the distribution of
surrogate field trajectory $\xi(t)$ (a continuous function of time) spanned on
a discrete time grid, such that
$\xi_{t_{1}}=\xi(t_{1}),\xi_{t_{2}}=\xi(t_{2}),\ldots,\xi_{t_{k}}=\xi(t_{k})$.

The density of the time grid is, of course, determined by the frequency at
which the projective measurements are performed. Therefore, to sample the
whole continuity of the trajectories’ course, the measurement frequency has to
be increased to infinity, in effect, realizing a scheme for continuous
monitoring of $\hat{V}_{E}$. However, if the surrogate trajectories are being
collected for the purposes of simulating the dynamics of open systems,
sampling on a discrete time grid is sufficient. On the one hand, if the time
grid is fine enough, then the functional integral over continuous trajectories
can be approximated by the ensemble average over discrete samples,

$\displaystyle\hat{\rho}_{S}(t)$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\xi\,\mathcal{P}[\xi]\hat{U}_{t}[\xi]\hat{\rho}_{S}\hat{U}_{t}^{\dagger}[\xi]\approx\int\mathcal{D}\xi\,\mathcal{P}[\xi]\left(\prod_{l=k}^{1}e^{-i(t_{l+1}-t_{l})(\hat{H}_{S}+\xi(t_{l})\hat{V}_{S})}\right)\hat{\rho}_{S}\left(\prod_{l^{\prime}={1}}^{k}e^{i(t_{l^{\prime}+1}-t_{l^{\prime}})(\hat{H}_{S}+\xi(t_{l^{\prime}})\hat{V}_{S})}\right)$
(23) $\displaystyle=$
$\displaystyle\sum_{\xi_{k},\ldots,\xi_{1}\in\Omega}p_{E}^{(k)}(\xi_{k}t_{k};\ldots;\xi_{1}t_{1})\left(\prod_{l=k}^{1}e^{-i(t_{l+1}-t_{l})(\hat{H}_{S}+\xi_{l}\hat{V}_{S})}\right)\hat{\rho}_{S}\left(\prod_{l^{\prime}=1}^{k}e^{i(t_{l^{\prime}+1}-t_{l^{\prime}})(\hat{H}_{S}+\xi_{l^{\prime}}\hat{V}_{S})}\right)$
$\displaystyle\approx$
$\displaystyle\frac{1}{N_{e}}\sum_{(\xi_{t_{1}},\ldots,\xi_{t_{k}})\in\mathcal{E}_{N_{e}}}\left(\prod_{l=k}^{1}\exp\left[-i(t_{l+1}-t_{l})(\hat{H}_{S}+\xi_{t_{l}}\hat{V}_{S})\right]\right)\hat{\rho}_{S}\left(\prod_{l^{\prime}=1}^{k}\exp\left[i(t_{l^{\prime}+1}-t_{l^{\prime}})(\hat{H}_{S}+\xi_{t_{l^{\prime}}}\hat{V}_{S})\right]\right),$

where $t_{k+1}=t$ and $\prod_{l=k}^{1}\hat{A}_{l}$
($\prod_{l^{\prime}=1}^{k}\hat{A}_{l^{\prime}}$) is understood as an ordered
product $\hat{A}_{k}\cdots\hat{A}_{1}$ ($\hat{A}_{1}\cdots\hat{A}_{k}$). On
the other hand, given the ensemble of sample trajectories, the noise
simulation is carried out by first solving the appropriate von Neumann
equation for every member of the ensemble with the surrogate trajectories
incorporated as an external time-dependent field, and then calculating the
ensemble average over so-obtained set of density matrices. To solve the
equations of motion in practice, one would employ some numerical method (e.g.,
the Runge-Kutta method) and, of course, any such method works in discrete time
steps—in other words, the numerical solution only requires trajectories to be
specified on a time grid. In such a case, the frequency of the
measurements—and thus, the density of the grid—would simply set the upper
bound for the accuracy of the numerical solutions.

We stress that the ensemble of trajectories sampled with the measurement
scheme described above is fit for purpose only when the noise representation
is valid [i.e., when (13) and (6) are true]. Of course, it is always possible
(at least in principle) to perform the sequential measurement to gather
results $(\xi_{t_{1}},\ldots,\xi_{t_{k}})$, even when
$\mathcal{Q}[\xi,\zeta]\neq\delta(\zeta-\xi)\mathcal{P}[\xi]$. The issue is
that, unless the noise representation is valid, there is no equivalence
between the average over the ensemble of samples and the functional “average”
over quasi-probability distribution $\mathcal{Q}$; such equivalence exists
only for proper probability distributions. In other words, measuring the
sample trajectories for the explicit purpose of simulating the dynamics of
systems coupled to $E$ through $\hat{V}_{E}$ can work only when such a
simulation can be performed in the first place (i.e., when the noise
representation is valid). However, even when (13) does not hold (or it is not
clear whether noise representation is valid) but one insists on treating the
measured sequences as sample trajectories, there might still be some utility
in employing them in a feigned simulation. First, such a simulation could be
viewed as a type of approximation to the dynamics of an open system. Second,
the feigned simulation could be used as a form of noise representation
witness: if one is able to detect that a simulation failed to predict some
measurable aspects of the dynamics, then it would constitute a proof that
noise representation is not valid. The potential value in the listed use cases
might warrant future investigation but, presently, their analysis lies beyond
the scope of this work.

##  III Illustration of trajectory sampling in numerical experiment

Here, we use an example of a numerical experiment to showcase how the noise
trajectory sampling and subsequent simulation of open-system dynamics could be
implemented in practice.

First, we perform the numerical simulation of sequential measurements of the
observable that can be “replaced” by the surrogate field; this step is meant
to represent the trajectory sampling stage. We take the model of the
environment $E$ that can be divided into two subspaces: the two-level system
$q$, and the large thermal bath $B$,

$\hat{H}_{E}=\hat{H}_{q}\otimes\hat{\mathds{1}}_{B}+\hat{V}_{qB}+\hat{\mathds{1}}_{q}\otimes\hat{H}_{B}.$
(24)

We chose to measure the observable that operates only on the subspace $q$,

$\displaystyle\hat{V}_{E}$ $\displaystyle=$
$\displaystyle\frac{\omega}{2}\hat{\sigma}_{z}\otimes\hat{\mathds{1}}_{B}=\left(\sum_{\xi=\pm\omega/2}\xi|\mathrm{sign}(\xi)\rangle\langle\mathrm{sign}(\xi)|\right)\otimes\hat{\mathds{1}}_{B}$
(25) $\displaystyle=$
$\displaystyle\left(\sum_{\xi=\pm\omega/2}\xi\,\hat{\Pi}(\xi)\right)\otimes\hat{\mathds{1}}_{B},$

where $\hat{\sigma}_{z}=(|{+}\rangle\langle{+}|-|{-}\rangle\langle{-}|)/2$.
Therefore, the spectrum of $\hat{V}_{E}$ is
$\Omega=\\{+\omega/2,{-}\omega/2\\}$ and the corresponding eigenstates are
$\\{|{+}\rangle,|{-}\rangle\\}$. We assume that the relationship between $q$
and $B$ is such that the dynamics of subsystem $q$,
$\hat{\rho}_{q}(t)=\mathrm{tr}_{B}(e^{-it\hat{H}_{E}}\hat{\rho}_{q}\otimes\hat{\rho}_{B}e^{it\hat{H}_{E}})$,
is accurately described by the master equation

$\frac{d}{dt}\hat{\rho}_{q}(t)=\mathscr{L}\hat{\rho}_{q}(t)=-\frac{\gamma}{2}[\hat{\sigma}_{x},[\hat{\sigma}_{x},\hat{\rho}_{q}(t)]],$
(26)

so that $\hat{\rho}_{q}(t)=e^{(t-s)\mathscr{L}}\hat{\rho}_{q}(s)$ for any
$t>s>0$. For simplicity we also assume that the state of $q$ is stationary,
$\mathscr{L}\hat{\rho}_{q}=0$, which implies that
$\hat{\rho}_{q}=\hat{\mathds{1}}/2$. It has been demonstrated in [11] that
such an environment facilitates valid noise representation where the
observable $\hat{V}_{E}$ is replaced by the surrogate field $\xi$ in a form of
random telegraph noise—a stochastic process that switches between
$\pm\omega/2$ at rate $\gamma$.

For this model, the probability of getting result $\xi_{t_{1}}$ in the
measurement of $\hat{V}_{E}$ simplifies as follows

$\displaystyle P(\xi_{t_{1}})$ $\displaystyle=$
$\displaystyle\mathrm{tr}\left[\hat{\Pi}(\xi_{t_{1}})\otimes\hat{\mathds{1}}_{B}\,e^{-it_{1}\hat{H}_{E}}\hat{\rho}_{q}\otimes\hat{\rho}_{B}\,e^{it_{1}\hat{H}_{E}}\right]$
(27) $\displaystyle=$
$\displaystyle\mathrm{tr}_{q}\left[\hat{\Pi}(\xi_{t_{1}})e^{t_{1}\mathscr{L}}\hat{\rho}_{q}\right],$

so that $B$ degrees of freedom are eliminated from the problem. Therefore, it
is enough to consider the a posteriori state reduced to subspace $q$,

$\hat{\rho}_{q}(t|\xi_{t_{1}})=e^{(t-t_{1})\mathscr{L}}\frac{\hat{\Pi}(\xi_{t_{1}})\hat{\rho}_{q}(t_{1})\hat{\Pi}(\xi_{t_{1}})}{P(\xi_{t_{1}})}.$
(28)

In general, the a posteriori state after $l$ measurements is given by the
recurrence relation,

$\displaystyle\hat{\rho}_{q}(t|\xi_{t_{1}},\ldots,\xi_{t_{l}})$
$\displaystyle\phantom{\rho_{q}}=e^{(t-t_{l})\mathscr{L}}\frac{\hat{\Pi}(\xi_{t_{l}})\hat{\rho}_{q}(t_{l}|\xi_{t_{1}},\ldots,\xi_{t_{l-1}})\hat{\Pi}(\xi_{t_{l}})}{P(\xi_{t_{l}}|\xi_{t_{1}},\ldots,\xi_{t_{l-1}})},$
(29)

with the conditional probability of getting the result $\xi_{t_{l}}$ that
reads

$\displaystyle
P(\xi_{t_{l}}|\xi_{t_{1}},\ldots,\xi_{t_{l-1}})=\mathrm{tr}_{q}[\hat{\Pi}(\xi_{t_{l}})\hat{\rho}_{q}(t_{l}|\xi_{t_{1}},\ldots,\xi_{t_{l-1}})].$
(30)

The numerical simulation of the sequential measurement is implemented in the
following way. We choose the interval between consecutive measurements
$\Delta=0.2\times\gamma^{-1}$ and the noise strength $\omega=2\gamma$
($\gamma^{-1}$ is used here as the unit of time). Then, we start by drawing
the first value $\xi_{t_{1}}^{(1)}$ at $t_{1}=0\times\Delta$ from calculated
probability distribution $P(\xi_{t_{1}}^{(1)})$ (recall, the initial state of
$E$ is stationary). Once the value of the first result is set, we then
calculate the a posteriori state $\hat{\rho}_{q}(\Delta|\xi^{(1)}_{t_{1}})$
according to Eq. (29), which, in this case, can be done analytically. Using
this state we calculate the probability distribution
$p(\xi^{(1)}_{t_{2}})=P(\xi^{(1)}_{t_{2}}|\xi_{t_{1}}^{(1)})$ according to
formula (30) (again, can be done exactly) and we draw from it the next result
$\xi_{t_{2}}$ at $t_{2}=\Delta$. These steps are repeated $48$ more times (for
the total of $k=50$ measurements) and we obtain the sequence
$(\xi^{(1)}_{t_{1}},\xi^{(1)}_{t_{2}},\ldots,\xi^{(1)}_{t_{50}})$—the first
sample trajectory. The whole procedure is repeated $1000$ times and each
obtained sample is stored in the ensemble
$\mathcal{E}_{1000}=\\{(\xi_{t_{1}}^{(j)},\ldots,\xi_{t_{50}}^{(j)})\\}_{j=1}^{1000}$;
the obtained samples are plotted in Fig. 1.

Figure 1: First five sample trajectories $\xi^{(j)}$ of the surrogate field
$\xi$ measured and stored in the ensemble $\mathcal{E}_{N_{e}}$ in the
numerical experiment. The trajectories switch at random between two values
$\pm\omega/2$. In this example, the projective measurements were carried out
every $0.2\,\gamma^{-1}$.

Now that we have populated the ensemble of trajectory samples for our model
environment, we can make use of it to calculate the noise representation
simulation of the open-system dynamics of any $S$ coupled to $E$ via
$\hat{V}_{E}$. For this demonstration, we have chosen $S$ as a qubit
interacting with the environment via so-called pure dephasing coupling,

$\hat{H}_{S}=0;\ \hat{V}_{S}=\frac{1}{2}\hat{\sigma}_{z};\
\hat{\rho}_{S}=\frac{1}{2}\left(\begin{array}[]{cc}1&1\\\ 1&1\\\
\end{array}\right),$ (31)

so that

$\hat{H}_{SE}=\frac{1}{2}\hat{\sigma}_{z}\otimes\hat{V}_{E}+\hat{\mathds{1}}\otimes\hat{H}_{E}\to\hat{H}_{t}[\xi]=\frac{1}{2}\xi(t)\hat{\sigma}_{z}.$
(32)

Pure dephasing implies that only the off-diagonal elements of the density
matrix are evolving,
$\langle{\pm}|\hat{\rho}_{S}(t)|{\pm}\rangle=\mathrm{constant}$. The choice of
$S$ was dictated by the fact that the dynamics can be solved exactly and the
state of the system is given by the analytical formula [17, 18], which can be
easily compared with the ensemble average.

We employ the trajectory ensemble to perform the simulation by solving the von
Neumann equation for each measured sample,

$\displaystyle\hat{\rho}_{S}^{(j)}(0)=\frac{1}{2}\left(\begin{array}[]{cc}1&1\\\
1&1\\\ \end{array}\right),$ (35)
$\displaystyle\frac{d}{dt}\hat{\rho}_{S}^{(j)}(t)=-\frac{i}{2}\xi^{(j)}(t)[\hat{\sigma}_{z},\hat{\rho}_{S}^{(j)}(t)],$
(36)

using Runge-Kutta method with time step $h=2\Delta$ and substituting the
trajectory values with the corresponding measured results,
$\xi^{(j)}(l\Delta)=\xi_{t_{l}}^{(j)}$. (Note that the Runge-Kutta algorithm
requires trajectories evaluated at $t=nh=2n\Delta$ and
$t=nh+h/2=(2n+1)\Delta$, hence the trajectories where sampled every
$\Delta=h/2$.) Obtained density matrices are collected in the ensemble
$\\{\hat{\rho}_{S}^{(j)}(t)\\}_{j=1}^{1000}$; the results of averages
calculated with first $10$, $100$ and full $1000$ elements are shown in Fig.
2. The accuracy of the noise representation simulation improves with the size
of the ensemble. Overall, the ensemble average exhibits fluctuations around
the exact value with rms that scales as $1/\sqrt{N_{e}}$—the shot noise
scaling due to averaging with a finite number of sample trajectories. However,
the error is not uniform in $t$: at short times scales ($t\sim\gamma^{-1}$,
$\omega^{2}t\lesssim\gamma$) the error is small even for modest ensemble
sizes, and it grows to maximal value at long times that coincide with the
significant decay of the coherence
$\langle{+}|\hat{\rho}_{S}(t)|{-}\rangle\sim 0$.

Figure 2: The simulated evolution of the system’s coherence
$\langle{+}|\hat{\rho}_{S}(t)|{-}\rangle$. The plot compares the exact result
$(1/2)e^{-\gamma t}\left[\cosh(\mu\gamma t)+\mu^{-1}\sinh(\mu\gamma t)\right]$
with $\mu=\sqrt{1-\omega^{2}/(4\gamma^{2})}$ (solid line) and ensemble
averages
$\hat{\rho}_{S}(t)\approx(1/N_{e})\sum_{j=1}^{N_{e}}\hat{\rho}_{S}^{(j)}(t)$
obtained for three values of ensemble size $N_{e}$: $10$ (green diamonds),
$100$ (red squares), and $1000$ (blue circles). The trajectory-wise density
matrices $\hat{\rho}_{S}^{(j)}(t)$ were obtained by numerical integration of
corresponding von Neumann equation (36).

##  IV Trajectory sampling and noise spectroscopy: a comparison

In essence, the trajectory sampling is a way to characterize the surrogate
field through direct measurement—a scheme for “measuring” the environmental
noise. However, the concept of “measuring” the noise is not new; the most
popular technique of this sort is the dynamical-decoupling-based noise
spectroscopy [19, 20, 21, 22, 23, 24, 25, 26]. It is instructive to compare
and contrast the two approaches. To this end, we first give a brief rundown of
the principle of operation for the spectroscopy.

Noise spectroscopy performs its “measurement” indirectly by utilizing a simple
qubit probe ($p$) brought in contact with the environment through pure
dephasing,

$\displaystyle\hat{H}_{pE}(t)=\frac{1}{2}f_{\mathrm{ctr}}(t)\hat{\sigma}_{z}\otimes\hat{V}_{E}+\hat{\mathds{1}}\otimes\hat{H}_{E}$
$\displaystyle\phantom{H}\to\hat{H}_{t}[\xi]=\frac{1}{2}f_{\mathrm{ctr}}(t)\xi(t)\hat{\sigma}_{z},$
(37)

where we assumed that $E$ has a valid noise representation with surrogate
field $\xi$. The real-valued filter function $f_{\mathrm{ctr}}(t)$
encapsulates the effects of the control protocol exerted over the qubit.
Typically, the employed control scheme has the form of a sequence of precisely
timed pulses that causes effectively instantaneous $\pi$ rotations of the
qubit’s Bloch vector (spin flips)—an approach inspired by the dynamical
decoupling techniques. The control is the key element of the method; usually,
the control is designed to exert a periodic modulation with a well-defined
frequency $\omega_{\mathrm{ctr}}$, so that $f_{\mathrm{ctr}}(t)$ can act as
narrow pass-band frequency filter. The pure dephasing coupling is chosen
because it leads to very simple qubit evolution where the filter action is
easily manipulated and the rate of qubit’s dephasing (the decay of off-
diagonal elements of the density matrix) is most directly related to noise-
characterizing quantities,

$\displaystyle\langle{+}|\hat{\rho}_{p}(t)|{-}\rangle$ $\displaystyle\propto$
$\displaystyle\exp\left[\sum_{k=1}^{\infty}\frac{(-i)^{k}}{k!}\int_{0}^{t}dt_{1}\cdots
dt_{k}\right.$ (38)
$\displaystyle\left.\times\left(\prod_{l=1}^{k}f_{\mathrm{ctr}}(t_{l})\right)\overline{\overline{\xi(t_{1})\cdots\xi(t_{k})}}\right].$

Here, $\overline{\overline{\xi(t_{1})\cdots\xi(t_{k})}}$ is the cumulant of
order $k$ of stochastic process $\xi$ [27]. Formally, the set of cumulants of
all orders (in general, there are infinitely many orders) fully defines the
stochastic process. The cumulant of a given order $k$ can always be expressed
as a combination of moments $\overline{\xi(t_{1})\cdots\xi(t_{l})}$ [see Eq.
(LABEL:eq:moment)] of orders $l\leqslant k$, even though the explicit formulas
get progressively more complex; e.g., the second cumulant, called the
autocorrelation function, is given by

$C(t_{1},t_{2})=\overline{\overline{\xi(t_{1})\xi(t_{2})}}=\overline{\xi(t_{1})\xi(t_{2})}-\overline{\xi(t_{1})}\,\overline{\xi(t_{2})}.$
(39)

The stated objective of noise spectroscopy is to characterize the noise by
recovering its cumulants. We illustrate how this approach works in practice
with a simple example: Consider the case when the qubit’s response (38) is
well described within Gaussian approximation where the contributions from
cumulants of order $>2$ are negligible. In addition, assume the noise is a
stationary process so that $C(t_{1},t_{2})=C(t_{1}-t_{2})$. When the filter
frequency $\omega_{\mathrm{ctr}}$ of the applied control sequence satisfies
$\omega_{\mathrm{ctr}}\gg 2\pi\tau_{c}^{-1}$ where $\tau_{c}$ is the noise
correlation time [the range of $C(t)$] and the total duration is set to
$t=2\pi n\omega_{\mathrm{ctr}}^{-1}$ where the integer $n\gg 1$ so that
$t\gg\tau_{c}$, then the dephasing rate is given by [23]

$\langle{+}|\hat{\rho}_{p}(t)|{-}\rangle\approx\langle{+}|\hat{\rho}_{p}(0)|{-}\rangle\exp\left[-\frac{4t}{\pi^{2}}S\left(\omega_{\mathrm{ctr}}\right)\right],$
(40)

where $S(\omega)=\int_{-\infty}^{\infty}e^{-i\omega s}C(s)ds$ is the noise
power spectral density—the Fourier transform of the autocorrelation function.
Therefore, by measuring the decay rate for a wide range of filter frequencies
$\omega_{\mathrm{ctr}}$ one performs a tomographic reconstruction of spectral
density (hence, the name noise spectroscopy). Since $S(\omega)$ contains
essentially the same information as $C(t)$, the final step of reverting back
to autocorrelation function through inverse Fourier transform is usually
skipped at this point.

Typically, the noise spectroscopy is utilized to extract information that
fills the gaps in the assumed model of the environment (e.g., positions of
spectral lines) [28, 29] and only rarely can it be used to characterize the
surrogate field with enough detail that would allow for the obtained
information to be reused for the purpose of simulating the dynamics of other
open systems. The reason is that, even though the full set of cumulants
defines the process completely, converting cumulants into an ensemble of
trajectories is an excessively difficult process (e.g., using the multivariate
version of Gram–Charlier A series) that is fraught with significant problems
with convergence and accuracy. Naturally, the quality of such a conversion
procedure is directly tied to the number of available cumulants; this turns
out to be the major drawback because, in practice, the noise spectroscopy
methods can feasibly reconstruct only a few first cumulants [30, 31]. In fact,
the standard designs of spectroscopy protocols do not work without Gaussian
approximation.

It is clear that a method for populating the trajectory ensemble through
direct measurement solves the issue of the conversion procedure by simply
circumventing it entirely. Moreover, the ensemble assembled in such a way not
only enables the simulation of open-system dynamics, but it also allows us to
calculate any quantity characterizing the noise itself, including moments,
power spectral density, and cumulants of arbitrary order (Gaussian
approximation is unnecessary); e.g., the autocorrelation function can be
obtained with the following formula:

$\displaystyle C(t_{1},t_{2})$ $\displaystyle\approx$
$\displaystyle\frac{1}{N_{e}}\sum_{j=1}^{N_{e}}\xi^{(j)}(t_{1})\xi^{(j)}(t_{2})$
(41)
$\displaystyle-\frac{1}{N_{e}^{2}}\sum_{j,j^{\prime}=1}^{N_{e}}\xi^{(j)}(t_{1})\xi^{(j^{\prime})}(t_{2}),$

where $\xi^{(j)}\in\mathcal{E}_{N_{e}}$. Therefore, the trajectory sampling
method could be considered as a next evolutionary step of noise spectroscopy.

##  V Conclusions

We have shown that, for the environment $E$ that supports the noise
representation $\xi$ of observable $\hat{V}_{E}$ [i.e., when, for an arbitrary
system $S$ coupled to $E$ via $\hat{V}_{E}$, the dynamics of
$\hat{\rho}_{S}(t)$ can be simulated with stochastic one-body Hamiltonian
$\hat{H}_{t}[\xi]$, see Eq. (1)], the result of sequential projective
measurement of $\hat{V}_{E}$ has a probability distribution identical to
trajectories of surrogate field $\xi$. Therefore, even though it is impossible
to observe fluctuation of the noise field with a probe coupled to $E$ [if $S$
is a probe, then $\hat{\rho}_{S}(t)$ is averaged over fluctuations of $\xi$],
we have demonstrated that, effectively, a single noise trajectory can be
exposed to direct observation through sequential measurement performed on $E$.

If noise trajectories can be measured, then they can be sampled and stored for
later use. One use case is to characterize the noise itself with quantities
such as moments, power spectral densities, or cumulants, all of which can be
easily calculated with sample trajectories [e.g., see Eq. (41)]. Thus,
trajectory sampling can be viewed as an upgrade to standard noise spectroscopy
that circumvents typical limitations, most notably, the reliance on Gaussian
approximation.

Here, we have focused on the utility of trajectory sampling in simulating the
dynamics of systems open to $E$. When the noise representation is valid, the
reduced state of open system is found by averaging the conditional unitary
evolution over surrogate field trajectories [see, Eq. (4)]. However, computing
the necessary probability distribution from first principles is prohibitively
difficult and requires the full knowledge of $E$. Moreover, even when the
distribution is supplied, averaging over it is technically impractical. All
these difficulties are circumvented with trajectory sampling: (i) the average
over distribution can be approximated by the average over finite number of
samples (the ensemble average); (ii) the samples are obtained from direct
measurement, and thus, the knowledge about environmental dynamics is not
necessary.

There are two main reasons why trajectory sampling was not considered before.
First, so far, there was no compelling argument against the view that the
trajectories are nothing more than artifacts of the noise representation, and
as such, they would be impossible to observe or measure. Indeed, the
observable quantities obtained through surrogate field simulation always
involve the average over trajectories [see Eq. (2)]. Therefore, there is no
measurement that could be performed on the system state $\hat{\rho}_{S}(t)$
that would be able to somehow skip the averaging and expose a single
trajectory of the noise. Moreover, given that the noise is, by its very
definition, only a surrogate for actual quantum system $E$, there seems to be
no reason to think of $\xi$ as anything more than a mathematical construct
used in an intermediate step of calculations that only at the terminus give
physically meaningful result. The second reason is that even if trajectories
could be measured somehow, it was unclear whether the noise that generated
them could be repurposed to simulate the dynamics of systems other than the
one used to carry out the measurement in the first place (i.e., the issue of
surrogate field’s objectivity). Presently, we know that the surrogate field is
an objective representation [11], and, thus, the latter objection can be put
to rest. The motivation for assembling the trajectory ensemble is sound, as
the ensemble can be utilized for setting up simulations of any system exposed
to the environment through coupling $\hat{V}_{E}$. The former issue, however,
is more subtle. Here, we have not settle the question whether noise
trajectories are “real” or “artificial.” In fact, we have not even engaged
with this question. What we have shown is that the result of a properly set up
sequential measurement of observable $\hat{V}_{E}$ has a probability
distribution identical to that of the hypothetical noise trajectory. On the
one hand, the fact that trajectories are identically distributed as the
measured sequences allows one to substitute the latter for the former (and
vice versa) when performing the ensemble average. On the other hand, the
equality of probability distributions does not imply that the noise trajectory
and the measured sequence are the same entity, and thus, the trajectory cannot
simply “inherit” tangibility from the measured sequence. After all, physical
setups required to perform the sequence of measurements differ substantially
from the setup where the object of investigation is the system state
$\hat{\rho}_{S}(t)$—in particular, the sequential measurement does not involve
system $S$.

## Acknowledgments

I would like to thank J. Krzywda and Ł. Cywiński for many insightful
discussions and for reviewing drafts of this paper.

## References

  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, and F. Jelezko, The quantum technologies roadmap: a european community view, New Journal of Physics 20, 080201 (2018).
  * Preskill [2018] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nature photonics 5, 222 (2011).
  * Pezze  _et al._ [2018] L. Pezze, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Reviews of Modern Physics 90, 035005 (2018).
  * Demkowicz-Dobrzanski [2020] R. Demkowicz-Dobrzanski, Quantum metrology in the era of quantum information, Bulletin of the American Physical Society 65, R07.00013 (2020).
  * Gisin and Thew [2007] N. Gisin and R. Thew, Quantum communication, Nature photonics 1, 165 (2007).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Northup and Blatt [2014] T. E. Northup and R. Blatt, Quantum information transfer using photons, Nature photonics 8, 356 (2014).
  * Schlosshauer [2007] M. A. Schlosshauer,  _Decoherence: and the quantum-to-classical transition_ (Springer Science & Business Media, 2007).
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Szańkowski and Cywiński [2020] P. Szańkowski and Ł. Cywiński, Noise representations of open system dynamics, Scientific Reports 10, 22189 (2020), 2003.09688 .
  * Feynman and Vernon [1963] R. Feynman and F. Vernon, The theory of a general quantum system interacting with a linear dissipative system, Annals of Physics 24, 118 (1963).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Path integral approach to quantum brownian motion, Physica A: Statistical Mechanics and its Applications 121, 587 (1983).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * Peres [2006] A. Peres,  _Quantum theory: concepts and methods (Vol. 57)_ (Springer Science & Business Media, 2006).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge university press, 2009).
  * Ramon [2012] G. Ramon, Dynamical decoupling of a singlet-triplet qubit afflicted by a charge fluctuator, Physical Review B 86, 125317 (2012).
  * Szańkowski  _et al._ [2013] P. Szańkowski, M. Trippenbach, and Y. B. Band, Spin decoherence due to fluctuating fields, Physical Review E 87, 052112 (2013).
  * Paz-Silva  _et al._ [2016] G. A. Paz-Silva, S.-W. Lee, T. J. Green, and L. Viola, Dynamical decoupling sequences for multi-qubit dephasing suppression and long-time quantum memory, New Journal of Physics 18, 073020 (2016).
  * Szańkowski  _et al._ [2017] P. Szańkowski, G. Ramon, J. Krzywda, D. Kwiatkowski, and Ł. Cywiński, Environmental noise spectroscopy with qubits subjected to dynamical decoupling, Journal of Physics Condensed Matter 29, 333001 (2017).
  * Paz-Silva  _et al._ [2017] G. A. Paz-Silva, L. M. Norris, and L. Viola, Multiqubit spectroscopy of Gaussian quantum noise, Physical Review A 95, 022121 (2017).
  * Ferrie  _et al._ [2018] C. Ferrie, C. Granade, G. Paz-Silva, and H. M. Wiseman, Bayesian quantum noise spectroscopy, New Journal of Physics 20, 123005 (2018).
  * Szańkowski and Cywiński [2018] P. Szańkowski and Ł. Cywiński, Accuracy of dynamical-decoupling-based spectroscopy of Gaussian noise, Physical Review A 97, 032101 (2018).
  * Norris  _et al._ [2018] L. M. Norris, D. Lucarelli, V. M. Frey, S. Mavadia, M. J. Biercuk, and L. Viola, Optimally band-limited spectroscopy of control noise using a qubit sensor, Physical Review A 98, 032315 (2018).
  * Krzywda  _et al._ [2019] J. Krzywda, P. Szańkowski, and Ł. Cywiński, The dynamical-decoupling-based spatiotemporal noise spectroscopy, New Journal of Physics 21, 043034 (2019).
  * Szańkowski [2019] P. Szańkowski, Transition between continuous and discrete spectra in dynamical-decoupling noise spectroscopy, Physical Review A 100, 052115 (2019).
  * Van Kampen [2011] N. Van Kampen,  _Stochastic Processes in Physics and Chemistry_ (North-Holland Personal Library (Elsevier Science), 2011).
  * Lovchinsky  _et al._ [2016] I. Lovchinsky, A. Sushkov, E. Urbach, N. P. de Leon, S. Choi, K. De Greve, R. Evans, R. Gertner, E. Bersin, C. Müller, _et al._ , Nuclear magnetic resonance detection and spectroscopy of single proteins using quantum logic, Science 351, 836 (2016).
  * Krzywda  _et al._ [2017] J. Krzywda, Ł. Cywiński, and P. Szańkowski, Localization of a magnetic moment using a two-qubit probe, Physical Review A 96, 042108 (2017).
  * Norris  _et al._ [2016] L. M. Norris, G. A. Paz-Silva, and L. Viola, Qubit Noise Spectroscopy for Non-Gaussian Dephasing Environments, Physical Review Letters 116, 150503 (2016).
  * Ramon [2019] G. Ramon, Trispectrum reconstruction of non-gaussian noise, Physical Review B 100, 161302(R) (2019).