# Containment strategies and statistical measures for the control of Bovine
Viral Diarrhea spread in livestock trade networks

Jason Bassett<EMAIL_ADDRESS>Pascal Blunk  Jörn Gethmann  Franz J.
Conraths  Philipp Hövel  Institut für Theoretische Physik, Technische
Universität Berlin, Hardenbergstrasse 36, 10623 Berlin, Germany  calimoto
GmbH, Babelsberger Straße 12, 14473 Potsdam, Germany  Institute of
Epidemiology, Friedrich–Loeffler–Institut, Südufer 10, 17493 Greifswald -
Insel Riems, Germany  School of Mathematical Science, University College Cork,
Cork T12 XF64, Ireland  U.S. Army Corps of Engineers Research and Development
Center, Concord Massachusetts, U.S.A.

###### Abstract

Assessing the risk of epidemic spread on networks and developing strategies
for its containment is of tremendous practical importance, both due to direct
effects in public health and its impact on economies. In this work we present
the numerical results of a stochastic, event-driven, hierarchical agent-based
model designed to reproduce the infectious dynamics of the cattle disease
called Bovine Viral Diarrhea (BVD), for which the corresponding movements’
network is the main route of spreading. For the farm-node dynamics, the model
takes into account a vast number of breeding, infectious and animal movement
mechanisms via a susceptible-infected-recovered (SIR) type of dynamics with an
additional permanently infectious class. The interaction between the farms is
described by a supply and demand farm manager mechanism governing the network
structure and dynamics. We discuss the disease and breeding dynamics, study
numerous mitigation strategies of present and past government regulations
taking Germany as a case-study country and perform a sensitivity analysis on
key parameters. We argue that the model, subject to calibration, has universal
predictive potential, can be extended to diseases beyond BVD and demonstrate
that appropriate measures can indeed lead to eradication regarding BVD. We
further present the time-series’ results of the model and conduct a
statistical analysis of and among the different mitigation strategies.

##  1 Introduction

Disease containment in the form of some intervention strategy is the main goal
of any biosecurity program. In the last decades efforts on the behalf of
governments have been implemented at country levels for the purpose of disease
containment in livestock, especially for the European case, with the advent of
integrated union policies. A major such disease that has been afflicting the
livestock industry is Bovine Viral Diarrhea (BVD) [1, 2, 3, 4].

BVD is a disease caused by the virus of the pestivirus family. The main
symptom of an ailing animal is diarrhea but general fatigue, fever and mucosal
secretions are among the accompanying inflictions. Contagion can occur from
all sorts of bodily fluids (including sneezing droplets) as well as from
contaminated material or agricultural equipment. The symptoms persist for one
to two weeks before the animal is recovered, provided complications, leading
to its death, do not occur. The main peculiarity of BVD is its persistence
when affecting early-stage embryos of carrying cows. If the embryos in
question have not yet developed fully their immune system it is highly likely
that they will carry the virus for the rest of their lives and be thus an
enormous source of infection, provided they are born. Respectively, cows
undergoing the disease are likely to be led to abortions [5, 6].

The impact of BVD on the agricultural industry, both beef and dairy, is
immense due to losses in animals as well as due to the reduced milk yield.
This has led to numerous studies exploring the impact of BVD, its containment
practices and potential to the economy for different countries [1, 2, 7, 8].
Such studies have been focused on the in-herd dynamics in an aggregated
fashion regarding the animal population [9], have explored the empirical
network effect of cattle movements on BVD spread with simple dynamics [10],
have combined the two aforementioned elements in metapopulation models
generically [11] and in a multiscale, data-driven fashion (Italy) [12] for
predicting the spatial spread of BVD investigating several intervention
strategies, and have formulated expert systems to test present and potential
policies based on data-driven simulations (Ireland) [13]. Furthermore, there
have been some instances of economic assessment of policies against BVD based
on data [14], one based on the model presented in [15] and its results
scrutinized in this paper [16].

Putting our model to the test on real data, we use the animal trade network of
Germany as a testbed. The current regulations on the containment of BVD are
based on the mandate issued over ten years ago and implemented in 2011. It has
focused on the detection and removal of persistently infected (PI) animals by
testing calves with an antigen (virus detection) test soon after their birth
[17]. An amendment was introduced in 2016 envisioning a quarantine of 40 days
on farms suspected of PI animals through testing [18]. The results have been
impressive in comparison to historical data, although they did not lead to
eradication [19].

With the current work we aim to employ the refined approach of an agent-based
model for the combined dynamics of animal movements and in-herd infections as
opposed to aggregate methods [20]. This sort of modeling allows for a great
level of detail to be included in the system. We additionally include a
simulation plan of scenarios encompassing a wide range of intervention
strategies emulating common biosecurity policies such as testing and
vaccination. From that point we give an overview and interpretation of the
resulting time-series for the global population as well as histograms for the
infected animal populations at the farm level. Finally, we compare the
statistical similarity of the time-series of the system running under no
containment strategy with all the rest of the scenarios in pairs, utilizing
the non-parametric Mann-Whitney U test.

##  2 Materials and methods

##### The model

We start with a short overview of the model that we developed. All the related
details following the ODD (Overview, Design concepts and Details) protocol
[21] can be found in a technical work of ours on the subject [15].

For the simulation of BVD spread through the German cattle movement network we
formulated an agent-based, stochastic, event-driven and hierarchical code in
C++. We decided for the model to be agent-based (the animals being agents)
with the aim to capture a great level of detail as reported in [20].
Furthermore, we introduced hierarchy in the sense that actions could
correspond to different aggregates of agents (animals), thus enabling actions
for a group of animals (herd), for a node in the network (the premises) and
for the systems as a whole (for introducing intervention strategies).
Stochasticity was something we also assumed necessary for both breeding and
infectious features, so as to adequately model their complex, real world
fluctuations. Moreover, we built the model in an event-driven fashion [22] in
order to benefit from the trade-off between following the trajectory of an
agent in continuous time and observing it only when a relevant, previously
scheduled event is executed [23].

Regarding the dynamics of the model, those were designated in two separate
processes, namely to that of the in-farm (premise) dynamics and to that of the
animal flow (movement) between the farms. For the former we employed the
infectious SIR scheme with a permanently infectious class found in Viet et al.
[24] (see figure 1), while for the latter a rule that every premise should
demand or offer animals only when its population lacks or exceeds a certain
quota respectively. In the case of the infectious, in-farm dynamics we assumed
well-mixed conditions, i.e. no spatial structure and equal probabilities for
all the susceptible animals to be infected. Accordingly, following standard
SIR dynamics we defined a stochastic, instantaneous transmission rate for the
infections in question, that dictated the inflow of infections in the pool of
susceptible animals for each farm in a Markovian manner [15, 23]. For the
latter dynamics, essentially the network dynamics, we supplied the farms,
assuming solely dairy farms due to their high significance in the disease
spread in Germany [3, 12], with the simplest supply and demand rule

(a) Vertical states.

(b) Combined horizontal and vertical states.

Figure 1: Schematics (state charts) of the vertical (through births) 1(a) and
combined with horizontal (through contact) 1(b) infectious states SIRP. The
dashed lines denote vertical transmission paths (i.e. from mother to calf),
while the solid horizontal ones. Note that a path to the demise of the animal
is also taken into account with the letter “D”. Additionally, the letters “V”
and “M” take into account temporary immunization effects induced by
vaccination and protection through maternal antibodies respectively.

$\textrm{surplus/deficit}=\left|N_{\textrm{quota}}-N_{\textrm{instant}}\right|$
(1)

i.e. the requirement that the number of animals offered or demanded (no
special criterion applied for their selection) should be determined by the
difference of an animal quota for each particular farm (specified upon
initialization) $N_{\textrm{quota}}$ with the number of animals in the farm at
the time of the movement $N_{\textrm{instant}}$. A surplus occurs when
$N_{\textrm{quota}}>N_{\textrm{instant}}$ and a deficit in the opposite case.
To moreover ascertain that the system is not dissipative nor inflationary
(i.e. that relation (1) for each farm is satisfied by the rest of the farms of
the system), we include a source and a drain farm, which physically correspond
to imports from outside the system (country) and slaughterhouses respectively.
The precise details for both the infectious and the network dynamics can be
found in our technical, ODD report [15].

In addition, the model includes system-level intervention strategies including
animal testing for their infectious status, non-targeted vaccination
(immunization) as well as individual farm isolation (quarantine) for
explicitly specified periods of time. Last but not least, we mention the
sensitivity analysis on control parameters such as infectious transmission
rate, vaccine and infectious identification (test) efficacy in [15]. The
related code is publicly available on GitHub at the URL
https://github.com/Yperidis/bvd_agent_based_model and follows the “Overview,
Design and Details” protocol [21].

##### Simulation setup

For setting up and initializing the model’s simulations we used expert opinion
and the literature [3, 9, 10, 12, 25, 26, 27]. The exact details for the
selection of all the parameters and the distributions employed can be found in
the model’s technical description [15]. The same holds for benchmarking with
previous works of the literature.

Here we only mention the minimal relevant settings, which can also be
referenced in the simulation scheduling as summarized in tables 1 and 2 [28].

  * 1.

The simulation was set to run for 20,000 time steps (the choice of the
simulation parameters makes the time step a calendar day as explained in
[15]).

  * 2.

The units used for all the simulation’s parameters mapped each time step to a
calendar day. Thus the entire runs accounted to approximately 55 years.

  * 3.

The involved tests’ (antigen and antibody) sensitivity was set to 99.8%, while
their specificity was assumed to be 100%.

  * 4.

The probability that a vaccination’s effect would successfully immunize the
animal in question was set to 99.85%.

  * 5.

The infectious transmission coefficients, as defined in [24], were set to
$\beta_{\textrm{TI}}=0.03$ and $\beta_{\textrm{PI}}=0.5$.

  * 6.

Only farms with populations over 10 animals have been included in the
simulation. Smaller farms are deemed to have a negligent contribution to the
spread of BVD in the network.

  * 7.

The periodicity with which all the farms were queried as to their population
status was set to 7 days.

For the terms old and new regulation as used in table 1, a detailed
explanation can be found in [28] and the original text of the regulations can
be found at [17] and [18].

For the initialization of the farms’ population we employed the farm size
distribution of all of Germany scaled down to the size of the federal state of
Thuringia. The scaling was performed due to computational resource limitations
and we comment on that in the “Discussion” section. Regarding the
initialization of the four different states of the population in each farm we
first divided the farms into two classes: those having PI animals and those
devoid of them. Then we set the first class to be a total of 2% of all the
farms. Furthermore, we set a 2% probability for the animals entering the
system in every instance of an import from the source farm to be PI, to
account for a minimal BVD spread breaching the biosurveillance border system
(this feature can be annulled and we make special reference to it in the
supplementary material). Finally, we set the states of the animals for the
population of each farm to be distributed according to the following
probabilities for the two classes: (S, TI, R,
PI)${}_{\textrm{PI}}$=$(0.46,0.06,0.46,0.02)$ and (S, TI, R,
PI)${}_{\textrm{PI-free}}=(0.79,0.005,0.205,0)$, where S stands for
“susceptible”, TI for “transiently infected” and R for “recovered” (lifelong
immunity). According to expert opinion from Dr. Gethmann of the Friedrich-
Loeffler Institute, the aforementioned figures in a normalized population of
farms and corresponding animals depicted the status of BVD in Germany upon the
commencement of the nationwide biosecurity measures in 2011.

Strategy | Description
---|---
1 | No control (baseline)
2 | Old regulation
3 | New regulation
4 | New regulation and vaccination
5a | New regulation and YCW with a semesterly frequency
5b | New regulation and YCW with an annual frequency
6a | New regulation, vaccination and YCW with a semesterly frequency
6b | New regulation, vaccination and YCW with an annual frequency
7 | YCW with a semesterly frequency
8 | Vaccination
9 | Vaccination and YCW with a semesterly frequency
Table 1: The different strategies comprising the scenarios presented in table
2. YCW stands for the “young calf window” testing protocol. Figure 2: The
scenario scheduling plan. STR stands for “strategy” as outlined in table 1. A
color block denotes the effect of the same strategy throughout the different
starting days.

##### Data sets

As implied in the previous paragraph, the data-driven part of the model
appeared upon the initialization with the input of a farm size distribution
from a csv formatted file. This file consisted of two columns of integer
numbers, one standing for farm number and the other for animal count. Each row
linked the number of farms with the number of animals they corresponded to.

For the scaling of the German farm size distribution ($\approx 12^{6}$ animals
in $\approx 157,000$ farms) to the size of the federal state of Thuringia
($\approx 350,000$ animals in $\approx 1,600$ farms) we simply multiplied all
the entries of the farm column with a fraction of the total farm count for
which the total animal count would be that of Thuringia ($\approx 350,000$).
This number turned out to be $\approx 0.0279$.

##### Identification of infected animals

The objective of any biosurveillance policy is to correctly detect infected
animals. In the case of BVD the related policies focus on the detection of the
PI animals either through an antigen or an antibody test [6]. Essentially the
difference of the two consist of the direct detection of the virus, which is
permanent in the PI animals, and the appearance of BVD virus antibodies in an
animal indicating its infection (therefore transient) at some past point
respectively.

Within the code we set two different modes of testing, accounting for both
antigen and antibody tests.

  * 1.

Antigen test – Regular blood test This is a test that identifies infected
(both transiently and persistently) animals. It is scheduled to be performed
at the birth of an animal and is not repeated throughout the animal’s
lifetime. If an animal is identified as infected it is sent to the
slaughterhouse within half a day. We envision additionally a challenge to a
positive test outcome as a rare event, that leads to a second antigen test
whose result determines irrevocable the fate of the animal in question. Due to
the transient nature of the infection in otherwise healthy animals, the
removal policy in question is targeted in its vast majority to the PI animals.

  * 2.

Antibody test – Young calf window (YCW) policy This is a semesterly or annual
test that identifies recovered animals similar to policies in Ireland and
Scotland [13, 29, 30]. It is performed at the specified times at the farm
level, for a random selection of animals as summarized in table 2 with a
methodology based on a hypergeomteric distribution sampling as outlined in
[31]. In case of a positive result an antigen test (regular blood test) is
scheduled for every animal in the farm that has not been already tested. Note
that, as in the antigen test case, in this manner only the PI animals are
targeted.

N |  $\leq$ 10 |  $\leq$ 20 |  $\leq$ 40 |  $\leq$ 80 |  $\leq$ 160 |  $>$ 180
---|---|---|---|---|---|---
n | 8 | 10 | 12 | 13 | 13 | 14
Table 2: Sample sizes given a farm population according to the young calf
window method. Population sizes (N) and the corresponding samples (n) needed
to be tested negative so that the prevalence of the infected animals in the
population will not exceed 20% with a confidence of 95%.

##### Hypothesis testing – Mann-Whitney U test

Having any two sets of series of data with no specified assumption as to the
form of their distribution (e.g. normal), one can use a non-parametric test to
query about them following the same distribution (null hypothesis or $H_{0}$
valid) or not (alternative hypothesis or $H_{1}$ valid). Such a statistical
test assuming independent samples is the Mann-Whitney U test [32, 33].

The whole reasoning behind the Mann-Whitney U test lies in sorting the parts
of the two independent (unpaired) samples in question in ascending order and
measuring their absolute differences, i.e. the shift of their underlying
distributions. The details of the implementation can be found in the original
work [32] as well as in textbooks [33]. In our case any two time-series of the
PI fraction of animals are the pairs of independent samples (because each one
is initiated with a different seed in the stochastic simulation). Since our
baseline scenario is No 1 we compare it pairwise to all the rest of the
scenarios, i.e. we perform the test 12 times for the 13 scenarios. Finally, we
set the significance level for the test to be 0.05 and run the wilcox.test
function in the R language to produce a p-value. A p-value less than 0.05 will
mean that the null hypothesis $H_{0}$ should be rejected otherwise it should
be accepted.

##  3 Results

##### Disease extinction

Having allowed the simulation to run as outlined in the plan 2 we obtained 13
different time-series for the four states of the system (susceptible,
transiently infected, PI and recovered) as exhibited in figures 5 and 6. Two
questions that immediately arise are firstly, whether the PI class of animals
is eradicated in any of the tested intervention strategies beyond a certain
point in time, and secondly how fast (i.e. at which time) this eradication
prevails.

A rather inconspicuous behavior that the system exhibited concerning the
spread of BVD through the PI animals was their introduction through the source
farm, taken that we had set a probability of 2% for the incoming animals to be
PI. All the animal introductions from the source farm for the given farm
distribution with which the simulation was initialized ended at the 1,884th
time step (10,159 introductions in total) and respectively the last PI animal
entering the system was observed at the 561st time step. The former numerical
figure demonstrates that for the given initialization the system quickly
settles on a self-sustained trade to satisfy relation (1) upon every
management period (week), while the latter signifies the end of the BVD spread
amplification from external triggering. To understand if the external PI
amplification intensifies the epidemics we also ran the simulation for a 0%
probability of PI introductions upon animals entering the system through the
source farm. In that case the introductions lasted until the 1,891st time step
(10,369 of total introductions), indicating again a state in time beyond which
the system trade is self-sustained. In this latter case of no PI introductions
in the system, the epidemic peak was slightly lowered compared to that of PI
animals being introduced in the system with a 2% probability (see the
description of scenario 1 further on). Yet the PI prevalence remained the same
in both cases by the 10,000th time step, rendering the effect of the PI
animals’ introduction negligible for the scenarios of table 2, which will be
the focus of any further analysis. For a minimal example of how the PI
prevalence is affected by PI introductions and their absence throughout the
simulation reference the supplementary material (“Effect of PIs from the
Source Farm”).

To address the question of PI eradication and its rate during the simulation,
we draw the histograms of the PI population fraction for all the farms in
ascending order of attributed ID numbers and for all the scenarios of plan 2,
as seen in figure 3. Technically speaking, the ID attribution to the farms
given from the farm size distribution is performed in a way such that the
ascension of IDs corresponds to larger farms in terms of population. Thus, the
aforementioned histogram also holds the answer to the farm size threshold
above which there is a non-zero PI prevalence in scenarios where eradication
is not achieved respectively.

The first direct result we notice by inspection of figure 3 is that PI
eradication is achieved for scenarios 5, 8, 9 and 13 (see again tables 1 and
2). What all these scenarios have in common throughout the duration of their
intervention strategies (i.e. beyond the 10,000th time step) is the
periodical, non-targeted vaccination of animals combined with some removal
through surveillance (testing) strategy, hinting that eradication cannot be
achieved without both strategies.

Regarding the question of speed towards the eradication of the PI animals,
scenario 8 holds the best records. This is intuitively expected as scenario 8
combines the new regulation with vaccination and the shortest (semesterly)
antibody test scheduling per farm. Similarly, a little longer (annual)
antibody test scheduling leads to the second best eradication time. What is
not so obviously expected is the difference between scenario 13 and 5, which
indicates that the YCW strategy when combined with vaccination (strategy 13)
achieves faster PI eradication than its counterpart of the new regulation
combined with vaccination (strategy 5). This is possibly a consequence of the
stronger effect of the YCW semesterly testing when compared to the average
animal arrival times, which in turns determines the time of the antigen test
in the new regulation. Note that the eradication time lags between the
rankings translate to roughly a year.

|  | Eradication time
---|---|---
Scenario ranking | 8 | 14,075
9 | 14,335
13 | 14,800
5 | 15,235
Table 3: Record times of PI eradication for the pertinent scenarios.

Finally, the threshold farm size for PI animals to exist at the end of the
simulation varies depending on the scenario we are looking at. In all the
relevant scenarios, it is intuitive to expect the PI prevalence shifted
towards the large farms of the input distribution, so that the SIR-sort of
dynamics can work better (well-mixed conditions). The sole scenario however
where the threshold is shifted to mid-size farms is the first, where the
system has been left to evolve freely. This is the case where farms with a
population of above 340 animals contribute to the persistence of the PI
animals throughout the simulation runs.

Figure 3: PI bar chart of the system’s farms at the final time step for
scenarios 1 to 13. Eradication is achieved for scenarios 5 (new regulation and
vaccination), 8 (vaccination and YCW with a semesterly period), 9 (vaccination
and YCW with an annual period) and 13 (new regulation, vaccination and YCW
with a semesterly period). See table 3.

##### Feature analysis of time-series

As far as time-series are concerned, there are two prominent, distinct
components we can discern. Firstly, the infectious states’ evolution and
secondly the population in absolute numbers. For the former we allow the
discussion to revolve around the normalized population.

Regarding the global evolution of the infectious states in figures 5, 6 (the
latter focuses only on the infectious states) we should make a number of
remarks for the interpretation of each scenario separately.

  * 1.

Scenario 1 (baseline)

After a well-understood [34, 35] initial epidemic transient peak ($\approx$
after 1,000 days) representing the outbreak from the initial infectious seeds
within the farm level (see the initial conditions in the Materials and Methods
section) and from in-between farms (through movements) the system’s infectious
state fractions tend to an equilibrium state, with a PI prevalence of around
1.1%. Correspondingly, the PI fraction distribution along all farms in figure
3 is the most pronounced of all the scenarios, which is logical as there is no
counter measure. Up to the 10,000th time step the behavior of all the
scenarios is identical.

  * 2.

Scenarios 2, 3

Beyond the 10,000th time step the recovered (R) fraction drops significantly
and the respective susceptible (S) increases. What happens is that due to the
old regulations’ effect the PI animals are being removed relatively fast
(faster in the new regulation than in the old) and therefore the naive
population, i.e. a population of S, is reinstated. When scenario 2 and 3 are
compared graphically on the same scale there is no significant difference
between them. However, the new regulation’s effect in scenario 3 (partly a
shorter time between two antigen tests and mostly the probation of a farm upon
a PI detection) is more pronounced towards the eradication of PI animals.
Taking into account the respective PI population fractions along the different
farms in figure 3, despite the fact that for scenario 3 we witness a higher PI
fraction in some farms than in scenario 2, still the PI fraction distribution
is narrower in scenario 3 than in scenario 2 supporting our previous
assessment that the new regulation should be more effective than the old one
towards the aim of the PI eradication. The higher PI prevalence in the farm
distribution of figure 3 for scenario 3 is also a result of the input farm-
size distribution’s heterogeneities. Overall, scenario 3 depicts the current
situation in Germany. The time frame between the 10,000th day and the 12,006th
day corresponds roughly to the five and a half years (the half year difference
with 2016 refers to a transition period from one mandate to the next) of the
old regulation’s effect in Germany.

  * 3.

Scenario 4

Beyond the 10,000th and up to the 12,373d time step we observe as expected the
same behavior as in scenario 3. After that point a behavior similar to that of
scenario 1 reappears due to the lift of the intervention strategies. The final
PI prevalence nevertheless appears to be slightly less than that of scenario
1, a fact which is reflected in the distribution width of the corresponding
scenario of figure 3 when compared to the width of the distribution for
scenario 1. This could be mainly attributed to the fact that at the point of
the simulation where the intervention strategies are lifted the population has
settled on a fixed point (see figure 4), which is robust to the variations we
implemented in the simulation plan (not presented due to minimal, if any,
observed changes when compared to the baseline scenario in figure 4). Thus the
animals’ movements could lead to fewer infections towards the steady state.

  * 4.

Scenario 5

The behavior for this scenario is similar to scenario 3 with the essential
difference starting from the 12,373d step where vaccination is added to the
new regulation. This action promotes the R population and conversely reduces
the S leading the curves to cross for a third time in figure 5 before they
settle on a steady state. The fact that not all of the animals become immune
at the steady state has to do with the interplay between the vaccination and
the breeding dynamics, which sustain a non-zero susceptible population. As
revealed in figure 3 this is one of the scenarios which leads to a PI
extinction at the end of the simulation.

  * 5.

Scenarios 6, 7

In scenarios 6 and 7 the effect of scenario 3 has nearly eradicated the PI
animals therefore not exhibiting any dramatic difference from the YCW’s
implementation from the 12,373th day onwards, either with a periodicity of a
semester (scenario 6 with strategy 5a) or of a year (scenario 7 with strategy
5b). The two scenarios’ behavior is similar to that of 3, with 6 (semesterly
YCW testing) leading to a higher naive population reestablishment than 7
(annual YCW testing). This last observation is easier to see in the respective
distribution of figure 3: the distribution of scenario 6 is narrower than that
of scenario 7.

  * 6.

Scenarios 8, 9

Similarly to scenarios 6 and 7, in scenarios 8 and 9 the effect of scenario 3
has nearly depleted the PI population by the 12,006th day to see any relevant
effect afterwards. However, the inclusion of vaccination leads to a behavior
similar to that of scenario 5 for that last segment of the simulation. The two
scenarios, like 5, also lead to a final state PI eradication, albeit at faster
rates than scenario 5.

  * 7.

Scenario 10

In scenario 10 the YCW strategy beyond the 12,006th day induces some
periodicity in the PI population. This is an interesting effect which is
attributed to the periodicity of the YCW test. The (declining) peaks’ width of
the PI population’s fraction in figure 5 signify the semesterly periodicity of
the YCW test, which would remove PI animals generated in the time frame of no
testing. This periodicity can also be seen by the declining ripples on the S
and R curves in the corresponding plot of figure 5. This last effect would be
induced by the infectious transmissions of the recurring PI animals. As
suggested in figure 6, the respective distribution of figure 3 indicates that
at the final state there is still some non negligent PI occurrence in the
farms. Although in global population terms this PI prevalence seems to be
rather low (less than 0.5%) it is rather alarming that PI animals are
relatively broadly distributed along the system’s farms, introducing a risk of
recurring future infections.

  * 8.

Scenario 11

Scenario 11 beyond the 12,373d day is again similar in behavior to scenario 5
concerning the S and R curves in the corresponding plot of figure 5, although
less pronounced than scenario 5. Looking closely at the respective plot of
figure 6 nonetheless we observe a short period where the PI population shoots
upwards and is replenished to around the final state without strategy before
dropping to minimal levels in the following segment of the simulation. This
has to do with the success probabilities of the tests’ accuracy and the
working probabilities of the vaccination on the PI prevalence, namely that the
indiscriminate vaccination scheduling protocol has a slower success rate than
that of the ear tag testing. Especially for the peak observed after the
12,373d day, it is a result of non-targeted vaccination (e.g. vaccination on
recovered pregnant cows that are destined to produce a PI will not change the
advent of the PI) in conjunction with the scheduling mechanism of the
vaccination itself, which is intertwined with the insemination, as the latter
needs to be scheduled at least 42 days after the vaccination. To demonstrate
the inefficiency of non-targeted vaccination in one last manner, although the
PI fraction has become very low, it is rather extraordinary that after 10,000
days (around 27 years) of the vaccination strategy’s effect there is still a
relatively (say to that of scenario 6) broad PI fraction distribution along
the system’s farms.

  * 9.

Scenario 12

Scenario 12 combines the new regulation with indiscriminate vaccination beyond
the 12,373d and up to the 12,738th day following similar results to scenario
11. However, the post 12,738th day behavior, which removes the new
regulation’s effect, exhibits again a similar, yet lower peak in PI as in
scenario 11 before reaching minimal levels of the PI population. This can be
attributed to a prolonged declining effect of PI-removed animals induced from
the prior enforcement of the new regulation. Otherwise, in terms of the final
state of the PI fraction’s distribution along farms in figure 3 the situation
is comparable to that of scenario 11 even with slightly broader and higher PI
fractions. This demonstrates how unpredictably and inefficiently non-targeted
vaccination can affect the PI prevalence of the system: even with a previously
combined effect of vaccination with the new regulation, when vaccination is
the sole intervention strategy it fails to lead to the eradication of the PI
population even in a time-span of 26 years.

  * 10.

Scenario 13

For scenario 13 the behavior is similar to scenario 12 up to the 12,738th day.
Beyond that we observe the first peak of the periodical PI increase (roughly
extending to a semester again, which is the period of the YCW test) as in
scenario 10, but with a smaller value due to the vaccination, which
accompanies the YCW strategy at this stage. Henceforth the PI population is
led to extinction by the vaccination’s effect as can be seen in figure 3,
underlying once again the importance of testing in the enforced strategy
towards PI eradication.

In respect to the population’s evolution in absolute numbers (figure 4), we
have already commented on it through the infectious states. After the
transient upon commencement due to the system’s initial conditions, the
population quickly settles on a fixed point, which is robust on the
perturbations introduced by the different intervention strategies at various
points in time. A striking feature that may be clearer after commenting is the
decaying periodicity of the ripples observed in figure 4. That is due to the
initialized insemination events for the cows and their well-defined life-
cycles, thus leading to the series of peaks, which are averaged collectively
in time.

Furthermore, in figures 7 and 8 we explored the variance of the PI prevalence
in the population (in gray) around a mean value and its distribution density
at $t=13,000$ (i.e. the epidemic outbreak’s peak exhibited in scenarios 4 and
10-13). The variance for each scenario was achieved by choosing 100 different
integer seeds (ranging from 2,333,600,960 to 2,333,601,060) for the stochastic
environment of the population and allowing it to run for each one of them. We
performed those runs for all the 13 scenarios of the simulation plan 5. The
fact that the PI population fraction mean is slightly dropping in time is not
alarming as this is a mere indication that the fixed point is being approached
at a rate so low that is insignificant for any practical purpose (i.e. policy
making). The reason for this slow rate of decline in the PI population may be
due to an underestimation of the transmission rates of BVD, which were taken
from the literature [24]. Note that the times shown in figure 7 are only
beyond the 9,000th time step. That is because the strategies vary only after
the 10,000th time step. Extensive sensitivity analysis on the vaccination and
testing can also be found in [15].

Finally, in figure 9 we see the results of the Mann-Whitman U test for all the
13 scenarios in 13 blocks and grouped in different colors with a temporal
increment (columns) of 1,000 time steps. Any changes start from the 10,000th
time step where strategies change according to table 2. When the test accepts
the null hypothesis for a number of scenarios, they are grouped with the same
color within the column. Then for each increment the figure delineates the
differentiation of the underlying population (i.e. the distribution of the
infectious states) to the next for each block-scenario.

Figure 4: Global population evolution for the baseline scenario 1. The rest of
the scenarios exhibit a statistically similar behaviour and are therefor not
shown. Figure 5: Susceptible (solid, blue line), transiently infected (dashed
red line), recovered (dotted-dashed, green line) and PI (dashed-dotted,
magenta line) fractions of the population for scenarios 1 to 13. Merely a
magnification of figure 4 around the infectious states. Figure 6: Transiently
infected (dashed red line) and PI (dashed-dotted, magenta line) fractions of
the population for scenarios 1 to 13. Figure 7: Variance around the mean PI
prevalence ($\%$) for the 13 different scenarios at $t=13,000$. Figure 8:
Variance distributions of the PI prevalence at t=13000 for the 13 different
scenarios. Figure 9: Similarity of the 13 scenarios of table 2 according to
the Mann-Whitman U test.

##  4 Discussion

Our results are among a series [9, 10, 13, 24] demonstrating the predictive
power and insight that computational models can provide in the efforts of
governance regarding epidemic spread containment in domesticated animal
populations. For the case of BVD, we defined general scenarios with disease
spread dynamics limited within the borders of a given country. These scenarios
can be affected by control measures as well as the import of PI animals, which
cannot be wholly prevented in practice when the laws of supply and demand are
in effect and the country in question has no say over the biosurveillance
regulations or their enforcement for the countries exporting to it.
Significantly, all the results pinpoint to the importance of combined herd
immunization (vaccination) and surveillance actions for achieving PI
eradication, given a high sensitivity of BVD tests. As expected, when combined
with vaccination, different surveillance practices lead to a range of record
times up to eradication, with an antibody test augmentation of the current
regulation topping all the alternatives. Furthermore, if vaccination is to be
combined with biosurveillance programs then a refinement of the non-targeted
practice can be certainly implemented with noticeable economical effects. In
any case, the PI eradication times have to always be weighted against the
foreseeable involved costs in the deployment of each scenario at the country
scale. For the case-study of Germany such an analysis can be found in [16].
Finally, for an optimal policy implementation the cost reduction induced by
antibody test practices (YCW strategy), namely the bulk milk tank or serum
test pool versus the individual blood tests should be taken into account.

Regarding the establishment of BVD in a farm population our results hinted on
a critical population size of the farms leading to a persistent epidemic
without intervention. In the cases where intervention still led to an endemic
final state the farm sizes contributing to the epidemics were shifted to the
largest farms of the systems. This is a consequence of the SIR dynamics taking
place at the farm level: a well-mixed, large population with endemic
conditions (infectious versus recovery rate ratio) including demographic
changes can exhibit an endemic equilibrium point [35]. Nevertheless, the
simulation being stochastic, one can expect fluctuations of the farm size
threshold for endemicity.

In terms of extensions of the current work, the influence of the scaling on
the predictions of the system should be investigated. Computationally the
inclusion of more farms requires exponentially more resources in memory even
for small farms, while increasing the number of animals in very few farms
depletes memory resources much slower. Scaling the system to less farms and
population would give a better understanding of the finite-size effects (lower
bound) on the epidemic dynamics. Another parameter that can be further
inquired is the sensitivity of the infection spread among farms to the farm
size distribution given initially, i.e. to the heterogeneity of the
distribution. Moreover, assuming the cattle movement rules are reasonably
similar in most E.U. countries due to unified mandates and similar practices,
the code employed in this work could be applied to extract results for at
least any country within the union. Additionally, the inclusion of subtleties
such as weighted or biased trading or dynamics at the node level would delve
to an ever greater level of detail, attempting to capture features such as the
preferential trade (attachment) or the observed declining trend of farms in
Germany in the last decade (according to the federal statistical bureau of
Germany) respectively. Last but not least and from a theoretical point of
view, the computation of a threshold quantity from parameters of the
simulation for the epidemic outbreak could be computed in a metapopulation
manner as in [36]. From there onwards various vaccination strategies for herd
immunization can be encompassed to alter the epidemic thresholds as in [37].

## Acknowledgments

J. Bassett and P. Hövel acknowledge the support of the German research
association within the framework of the collaborative research center 910.

## References

  * [1] I. Greiser-Wilke, B. Grummer, V. Moennig, Bovine viral diarrhoea eradication and control programmes in europe, Biologicals 31 (2) (2003) 113–118.
  * [2] K. Ståhl, S. Alenius, BVDV control and eradication in Europe – an update, Jpn. J. Vet. Res. 60 (2012) S31–S39.
  * [3] J. Gethmann, T. Homeier, M. Holsteg, H. Schirrmeier, M. Saßerath, B. Hoffmann, M. Beer, F. J. Conraths, BVD-2 outbreak leads to high losses in cattle farms in Western Germany, Heliyon 1 (1) (2015) e00019.
  * [4] V. Richter, E. Kattwinkel, C. L. Firth, T. Marschik, M. Dangelmaier, M. Trauffler, W. Obritzhauser, W. Baumgartner, A. Käsbohrer, B. Pinior, Mapping the global prevalence of Bovine Viral Diarrhoea virus infection and its associated mitigation programme, Vet. Rec. (2019) vetrec–2019.
  * [5] A. L. E. Lindberg, Bovine viral diarrhoea virus infections and its control. a review, Veterinary Quarterly 25 (1) (2003) 1–16, pMID: 12670010.  arXiv:http://dx.doi.org/10.1080/01652176.2003.9695140, doi:10.1080/01652176.2003.9695140.
  * [6] S. R. Lanyon, F. I. Hill, M. P. Reichel, J. Brownlie, Bovine viral diarrhoea: pathogenesis and diagnosis, Vet. J. 199 (2) (2014) 201–209.
  * [7] I. M. G. A. Santman-Berends, M. H. Mars, L. Van Duijn, G. van Schaik, Evaluation of the epidemiological and economic consequences of control scenarios for bovine viral diarrhea virus in dairy herds, J. Dairy Sci. 98 (11) (2015) 7699–7716.
  * [8] B. Pinior, C. L. Firth, V. Richter, K. Lebl, M. Trauffler, M. Dzieciol, S. E. Hutter, J. Burgstaller, W. Obritzhauser, P. Winter, et al., A systematic review of financial and economic assessments of bovine viral diarrhea virus (BVDV) prevention and mitigation activities worldwide, Prev. Vet. Med. 137 (2017) 77–92.
  * [9] A. Damman, A. F. Viet, S. Arnoux, M. C. Guerrier-Chatellet, E. Petit, P. Ezanno, Modelling the spread of bovine viral diarrhea virus (BVDV) in a beef cattle herd and its impact on herd productivity, Vet. Res. 46 (1) (2015) 12\.
  * [10] M. R. Tinsley, F. I. Lewis, F. Brülisauer, Network modeling of BVD transmission, Vet. Res. 43 (1).  doi:10.1186/1297-9716-43-11.
  * [11] A. Courcoul, P. Ezanno, Modelling the spread of bovine viral diarrhoea virus (bvdv) in a managed metapopulation of cattle herds, Vet. Microbiol. 142 (1-2) (2010) 119–128.
  * [12] B. Iotti, E. Valdano, L. Savini, L. Candeloro, A. Giovannini, S. Rosati, V. Colizza, M. Giacobini, Farm productive realities and the dynamics of bovine viral diarrhea (bvd) transmission, bioRxiv (2017) 230045.
  * [13] H. H. Thulke, M. Lange, J. A. Tratalos, T. A. Clegg, G. McGrath, L. O’Grady, P. O’Sullivan, M. L. Doherty, D. A. Graham, S. J. More, Eradicating, reviewing irish programme data and model predictions to support prospective decision making, Prev. Vet. Med. 150 (2017) 151.
  * [14] T. Marschik, W. Obritzhauser, P. Wagner, V. Richter, M. Mayerhofer, C. Egger-Danner, A. Käsbohrer, B. Pinior, A cost-benefit analysis and the potential trade effects of the bovine viral diarrhoea eradication programme in Styria, Austria, The Vet. J. 231 (2018) 19–29.
  * [15] J. Bassett, P. Blunk, T. Isele, J. Gethmann, P. Hövel, An agent-based model for bovine viral diarrhea, arXiv preprint arXiv:1812.06964.
  * [16] J. Gethmann, C. Probst, J. Bassett, P. Blunk, P. Hövel, F. J. Conraths, An epidemiological and economic simulation model to evaluate strategies for the control of bovine virus diarrhea in germany, Frontiers in Vet. Sci. 6 (2019) 406\.
  * [17] F. Ministry of Food, Agriculture, Verordnung zum Schutz der Rinder vor einer Infektion mit dem Bovinen Virusdiarrhoe-Virus (2008).
  * [18] F. Ministry of Food, Agriculture, Zweite Verordnung zur Änderung der BVDV-Verordnung (2016).
  * [19] K. Wernike, J. Gethmann, H. Schirrmeier, R. Schröder, F. J. Conraths, M. Beer, Six years (2011–2016) of mandatory nationwide bovine viral diarrhea control in Germany – a success story, Pathogens 6 (4) (2017) 50.
  * [20] T. Bosse, S. W. Jaffry, G. F. Siddiqui, J. Treur, Comparative analysis of agent-based and population-based modelling in epidemics and economics, Multiagent and Grid Systems 8 (3) (2012) 223–255.
  * [21] V. Grimm, U. Berger, D. L. DeAngelis, J. G. Polhill, J. Giske, S. F. Railsback, The ODD protocol: a review and first update, Ecol. Model. 221 (23) (2010) 2760–2768.
  * [22] G. Fishman, Discrete-event simulation: modeling, programming, and analysis, Springer Science & Business Media, 2013.
  * [23] C. L. Vestergaard, M. Génois, Temporal Gillespie Algorithm: Fast Simulation of Contagion Processes on Time-Varying Networks, PLOS Comput. Biol. 11 (10) (2015) e1004579.  doi:10.1371/journal.pcbi.1004579.
  * [24] A. F. Viet, C. Fourichon, H. Seegers, C. Jacob, C. Guihenneuc-Jouyaux, A model of the spread of the bovine viral-diarrhoea virus within a dairy herd, Prev. Vet. Med. 63 (3) (2004) 211–236.
  * [25] P. Ezanno, C. Fourichon, A. F. Viet, H. Seegers, Sensitivity analysis to identify key-parameters in modelling the spread of bovine viral diarrhoea virus in a dairy herd, Prev. Vet. Med. 80 (1) (2007) 49–64.
  * [26] P. Hoscheit, S. Geeraert, G. Beaunée, H. Monod, C. A. Gilligan, J. A. N. Filipe, E. Vergu, M. Moslonka-Lefebvre, Dynamical network models for cattle trade: towards economy-based epidemic risk assessment, J. Complex Net. (2016) cnw026.
  * [27] L. Bioglio, M. Génois, C. L. Vestergaard, C. Poletto, A. Barrat, V. Colizza, Recalibrating disease parameters for increasing realism in modeling epidemics in closed settings, BMC Infect. Dis. 16 (1) (2016) 676.
  * [28] J. Bassett, Computational and analytical approaches towards epidemic spread containment of temporal animal trade networks, Ph.D. thesis, Technische Universität Berlin (2019).
  * [29] F. Brülisauer, F. I. Lewis, A. G. Ganser, I. J. McKendrick, G. J. Gunn, The prevalence of bovine viral diarrhoea virus infection in beef suckler herds in Scotland, Vet. J. 186 (2) (2010) 226–231.
  * [30] J. A. Tratalos, D. A. Graham, S. J. More, Patterns of calving and young stock movement in Ireland and their implications for BVD serosurveillance, Prev. Vet. Med. 142 (2017) 30–38.
  * [31] F. J. Conraths, J. Gethmann, Epidemiologische Untersuchungen in Tierpopulationen, Tech. rep., FLI, Greifswald - Insel Riems (August 2015).
  * [32] H. B. Mann, D. R. Whitney, On a test of whether one of two random variables is stochastically larger than the other, Ann. Math. Statist. 18 (1) (1947) 50–60.  doi:10.1214/aoms/1177730491.
URL https://doi.org/10.1214/aoms/1177730491

  * [33] D. Wackerly, W. Mendenhall, R. L. Scheaffer, Mathematical statistics with applications, Cengage Learning, 2014.
  * [34] C. D. Keeling, T. P. Whorf, The implications of network structure for epidemic dynamics, Theoretical population biology.
  * [35] P. Rohani, Modeling infectious diseases in humans and animals, Princeton University Press, Princeton, 2008.
  * [36] V. Colizza, A. Vespignani, Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations., J. Theor. Biol. 251 (3) (2008) 450–67.  doi:10.1016/j.jtbi.2007.11.028.
  * [37] T. House, M. J. Keeling, Epidemic prediction and control in clustered populations, J. Theor. Biol. 272 (1) (2011) 1–7.

## Supplementary Material

### Single Farm Dynamics and PI Effect from Source Farm

#### Single Farm Dynamics

We develop some intuition on a minimal example of how the system should behave
in terms of infectious states and population in figures 10(a), 10(b) and 10(c)
respectively. The simulation runs for the same settings as in the main text
without any intervention strategy, as calibrated for the case of Germany in
[15] and with no PI animals originating from the source farm. The population
of the single farm is set to 1,000 animals to achieve well-mixed conditions
for the epidemic dynamics and to diminish finite-size effects. Since there is
an inherent supply and demand mechanism to equalize the imbalance caused by
demographic factors in the system, the minimum working setup requires the
source and the drain farms as well. We notice therefore that the population
fluctuates around the initial value of 1,000.

Moreover, the peak of infection is much more pronounced when only a single
farm is considered as can be seen in figure 10(b). This effect is typical of
SIR dynamics and depends on the involved infectious and recovery rates [35].
The spatial structure may then shift the outbreak peak in value and time
depending on the exact heterogeneities (farm size distribution) and
connectivity (supply and demand) of the farm-nodes involved in the simulation.

(a) All states’ evolution.

(b) Infectious states’ evolution.

(c) Demographic evolution.

Figure 10: All the S, I, R and P states in a single graph 10(a) for the single
farm with no PI inputs and only the infectious I and P in 10(b). Moreover, we
see in 10(c) the evolution of the same farm’s population for no intervention
strategy with animal inputs and outputs (source and drain farm) to satisfy its
demand that its population be constant.

#### Effect of PIs from the Source Farm

For the single farm with 2% probability of introduction of PI animals through
the source farm, the demographic time-series is virtually identical to the one
for no PI animals entering the system through the source farm. The settings
remain always as those in the main text without any intervention strategy.

(a) All states’ evolution.

(b) Infectious states’ evolution.

Figure 11: The respective evolution plots of figure 10 but for a probability
of a PI introduction from the source farm of 2%.

What stands out in figures 11(a) and 11(b) in respect to their counterparts
10(a) and 10(b) (the demographics are statistically identical in both cases
and are therefore skipped in figure 11) of no PI animals’ introduction from
the source farm is the epidemic outbreak peak value and the effect this has on
the immune, R class of animals. The PI peak is over twice times higher than in
the former case and therefore the R class of animals remains consistently at
high levels, without much fluctuation. This is a consequence of the increased
secondary, transient transmissions induced to susceptible animals from the
case of introduced PI animals in the system. These behaviors are both
intuitively expected as the number of secondary, transient infections should
be mainly proportional to the number of PI animals in the system, them being
the main contributor of the epidemic spread for both the well-mixed (local)
and spatially distributed (network) dynamics.

### Sensitivity Analysis

In the following we make a sensitivity analysis on four factors dominating
over the epidemic, testing and vaccination dynamics, as well as on the
stochastic fluctuations of the simulations overall.

  1. 1.

Transmission Rates

The catalytic factor affecting the final state of the PI animals is shown to
be the PI transmission coefficient $\beta_{PI}$ as demonstrated in table 4.
This result is somewhat contrary to that of [25] mainly because the between-
herd (animal group in the authors’ case) dynamics were not considered in our
simulation as opposed to that work. To be precise, despite the form of the
infectious dynamics being nearly identical (see [15]), the authors of [25]
considered only a five animal herd group, therefore introducing pronounced
finite-size effects in their results. Furthermore, in the simulation of our
work, due to the effect of vertical transmissions, i.e. PI animals latently
appearing in the system after a pregnant cow’s infection, the effect of the PI
transmission coefficient $\beta_{\text{PI}}$ dominates over
$\beta_{\text{TI}}$. In addition, the sensitivity analysis of [25] on the PI
animals’ death rates is not directly comparable to our work, as we presume a
uniform random distribution for that purpose.

PIfinal No | $\beta_{\text{TI}}$ | $\beta_{\text{PI}}$
---|---|---
11 | 0.01 | 0.1
210 | 0.01 | 0.5
278 | 0.01 | 0.8
12 | 0.03 | 0.1
305 | 0.03 | 0.5
12 | 0.05 | 0.1
317 | 0.05 | 0.8
21 | 0.1 | 0.1
Table 4: Final number of persistently infected animals as a function of
$\beta_{\text{TI}}$ and $\beta_{\text{PI}}$.

  2. 2.

Ear Tag and Retesting Laps

For the tests’ effect on the PI prevalence their sensitivity was examined by
varying the test’s accuracy probability as seen in figure 12. Not
surprisingly, as the sensitivity probability tends to unity the PI prevalence
tends to eradication as all the infected animals are identified and
subsequently sent to the drain farm (slaughterhouse).

Figure 12: Scenario 3 (reference the simulation plan from the main text) for different test accuracy probabilities ranging from 0 to 1, as outlined in table 5. The effect of the ear tag test protocol is enforced from day 10,000 onwards. Probability | Style (Position)
---|---
0 | Solid, blue (left)
0.1 | Dashed, red (left)
0.2 | Dotted-dashed, green (left)
0.3 | Dashed, dotted, magenta (left)
0.8 | Solid, blue (right)
0.98 | Dashed, red (right)
0.99 | Dotted-dashed, green (right)
0.998 | Dashed, dotted, magenta (right)
1 | Dotted, black (right)
Table 5: Probabilities of the test’s sensitivity for scenario 3 (reference the
simulation plan from the main text).

  3. 3.

Vaccination Success Probability

For the vaccination’s effect on the PI prevalence the vaccination’s working
probability was varied from zero to one as demonstrated in figure 13. As with
the tests’ sensitivity success probability, the increase of certainty for the
vaccine’s working probability leads to eventual extinction of the PI
population. However, the vaccination’s effect on the PI population is slower
than that of the tests’ sensitivity success and requires near certainty values
to lead to extinction. This could be attributed to the fact that
indiscriminate vaccination of animals regardless of their infectious status
reaches its intended targets slower than tests of all animals, which aim to
identify and remove the infected ones.

Figure 13: Strategy 8 (reference the simulation plan from the main text) for different vaccination working probabilities ranging from 0 to 1 as outlined in table 6. Probability | Style (Position)
---|---
0 | Solid, blue (left)
0.1 | Dashed, red (left)
0.2 | Dotted-dashed, green (left)
0.3 | Dashed, dotted, magenta (left)
0.8 | Solid, blue (right)
0.9 | Red, dashed (right)
0.985 (default) | Dotted-dashed, green (right)
1 | Dashed-dotted, magenta (right)
Table 6: Vaccination working probabilities for the sensitivity of strategy 8
(reference the simulation plan from the main text).

  4. 4.

Variance

For completeness and reliability of the stochastic nature of the simulation a
sensitivity analysis was performed on the PI prevalence percentage as a
function of the seeds of the random number generator as displayed in figure
14. The differences proved to be of the order of 0.001% implying that the
results are robust to stochastic fluctuations.

Figure 14: Scenario 1 (reference the simulation plan from the main text) for different seeds as outlined in table 7. Seed No | Style (Position)
---|---
2,333,600,960 | Solid, blue (left)
2,333,600,970 | Dashed, red (left)
2,333,601,960 | Dotted-dashed, green (left)
2,333,620,963 | Dashed, dotted, magenta (left)
2,333,710,962 | Dotted, black (left)
2,333,970,967 | Solid, blue (right)
2,333,650,932 | Dashed, red (right)
4,333,600,860 | Dotted-dashed, green (right)
3,323,601,969 | Dashed, dotted, magenta (right)
3,323,601,970 | Dotted, black (right)
Table 7: Seeds for the variance of scenario 1 (baseline).


# Dissipation and fluctuations in elongated bosonic Josephson junctions

F. Binanti  Dipartimento di Fisica e Astronomia ’Galileo Galilei’, Universita
di Padova, via Marzolo 8, 35131 Padova, Italy  K. Furutani  Dipartimento di
Fisica e Astronomia ’Galileo Galilei’, Universita di Padova, via Marzolo 8,
35131 Padova, Italy  INFN - Sezione di Padova, via Marzolo 8, 35131 Padova,
Italy  L. Salasnich  Dipartimento di Fisica e Astronomia ’Galileo Galilei’,
Università di Padova, via Marzolo 8, 35131 Padova, Italy  INFN - Sezione di
Padova, via Marzolo 8, 35131 Padova, Italy  Padua Quantum Technologies
Research Center, Universita di Padova,
via Gradenigo 6/b, 35131 Padova, Italy  CNR-INO, via Nello Carrara 1, 50019
Sesto Fiorentino, Italy

###### Abstract

We investigate the dynamics of bosonic atoms in elongated Josephson junctions.
We find that these systems are characterized by an intrinsic coupling between
the Josephson mode of macroscopic quantum tunneling and the sound modes. This
coupling of Josephson and sound modes gives rise to a damped and stochastic
Langevin dynamics for the Josephson degree of freedom. From a microscopic
Lagrangian, we deduce and investigate the damping coefficient and the
stochastic noise, which includes thermal and quantum fluctuations. Finally, we
study the time evolution of relative-phase and population-imbalance
fluctuations of the Josephson mode and their oscillating thermalization to
equilibrium.

###### pacs:

03.75.Lm; 03.75.-b; 05.40.-a

##  I Introduction

In the last forty years dissipative quantum systems under the effect of random
noise have been extensively investigated both theoretically and experimentally
caldeira ; ford ; schmidt ; bulgadaev ; dalla . An important tool in the study
of dissipative systems is the Caldeira-Leggett model caldeira , where a bath
of harmonic oscillators is coupled to the main system. This bath leads to
damping and noise in the main system and the dynamics is governed by the
generalized Langevin equation caldeira ; ford . One can find such dynamics in
many different physical systems, for instance in superconducting Josephson
circuits koch80 ; koch82 .

Recently, the dynamics of ultracold atoms in a double-well potential (bosonic
Josephson junction) has attracted marked attention minguzzi ; pigneur ; tononi
. Ref. minguzzi  has revealed the presence of intrinsic coupling between the
Josephson mode and sound modes leading to damping of oscillations of
population imbalance in a one-dimensional bosonic Josephson junction. The
peculiarity of this system is the intrinsic coupling between the Josephson
mode and the bath which is distinct from a system in contact with an external
bath like the Caldeira-Leggett model caldeira . Namely, the one-dimensional
Bose Josephson junction exhibits quantum Langevin dynamics in the Josephson
mode without any external reservoir. Furthermore, the typical approach of a
Caldeira-Leggett model consists in the introduction of an external bath in a
phenomenological way, as well as for the interaction between the system and
the environment. Here instead we manage to pull out an intrinsic phonon bath
(interacting with the Josephson modes) straightly from the model.

In this paper we start from a model of two one-dimensional quasi-condensates
with a Josephson coupling in the head-to-tail configuration discussed in Refs.
minguzzi ; tononi . We treat phase and number density as classical fields to
obtain the effective Lagrangian of relative phase. By adopting this different
procedure with respect to Ref. minguzzi , we show that this system involves
the intrinsic coupling between the Josephson mode and other infinite numbers
of sound modes. Consequently, the dynamics of the Josephson mode is described
by the quantum Langevin equation with damping and noise that includes thermal
and quantum fluctuations. We find an analytic formula for the damping
coefficient moreover, we derive analytic expressions for the fluctuations of
the Josephson relative phase and Josephson population imbalance in the high-
temperature regime, extending the results of Ref. minguzzi . Finally, we
analyze numerically the fluctuations of the Josephson mode in the general
case, in particular in the low-temperature regime where quantum fluctuations
play a crucial role. Remarkably, in this low-temperature regime, an
ultraviolet cutoff for the sound modes is needed. In this way, we determine
the thermalization to the equilibrium of the Josephson fluctuations.

##  II Bosons with Josephson tunnelling

We start from the following Lagrangian density which consists of two weakly-
interacting Bose-Einstein quasi-condensates ($j=1,2$) made of atoms with mass
$m$ in one dimension

$\displaystyle\mathscr{L}$
$\displaystyle=\sum_{j=1}^{2}\left[i\hbar\psi_{j}^{*}\partial_{t}\psi_{j}-\frac{\hbar^{2}}{2m}\absolutevalue{\partial_{x}\psi_{j}}^{2}-\frac{g}{2}\absolutevalue{\psi_{j}}^{4}\right]$
(1)
$\displaystyle+\frac{J(x)}{2}\left[\psi_{1}^{*}\psi_{2}+\psi_{2}^{*}\psi_{1}\right],$

where $\hbar$ is the Planck constant. Here $g$ is the strength of the inter-
atomic potential in the contact-interaction approximation and $J(x)$ is the
space dependent tunneling-energy coupling. The complex field $\psi_{j}(x,t)$
of the $j$-th quasi-condensate can be rewritten by means of the Madelung
representation

$\psi_{j}(x,t)=\sqrt{\rho_{j}(x,t)}\,e^{i\phi_{j}(x,t)},$ (2)

where $\rho_{j}(x,t)=|\psi_{j}(x,t)|^{2}$ is its atomic density. Substituting
the expression of Eq. (2) into the Lagrangian in Eq. (1), we obtain

$\displaystyle\mathscr{L}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{2}\Bigg{[}\frac{i\hbar}{2}\partial_{t}\rho_{j}-\hbar\rho_{j}\partial_{t}\phi_{j}-\frac{\hbar^{2}}{2m}\Bigg{[}\frac{1}{4\rho_{j}}(\partial_{x}\rho_{j})^{2}$
(3) $\displaystyle+$
$\displaystyle\rho_{j}(\partial_{x}\phi_{j})^{2}\Bigg{]}-\frac{g}{2}\rho_{j}^{2}\Bigg{]}+J(x)\sqrt{\rho_{1}\rho_{2}}\cos(\phi_{1}-\phi_{2}).$

A compact description of the system is reached when we introduce the relative
phase and the population imbalance

$\displaystyle\phi$ $\displaystyle=$ $\displaystyle\phi_{1}-\phi_{2},$ (4)
$\displaystyle\zeta$ $\displaystyle=$
$\displaystyle\frac{\rho_{1}-\rho_{2}}{2\bar{\rho}},$ (5)

with $\bar{\rho}=(\rho_{1}+\rho_{2})/2$ the average atomic density.
Introducing also a total phase $\bar{\phi}=\phi_{1}+\phi_{2}$, one can express
the Lagrangian in Eq. (3) in terms of these quantities. In addition, let us
work in a canonical ensemble (no particle leaves or enters the system) for
which $\partial_{t}\bar{\rho}=0$ holds, and assuming it also for the total
phase $\partial_{t}\bar{\phi}=0$, we obtain a new Lagrangian density

$\displaystyle\mathscr{L}$ $\displaystyle=$
$\displaystyle-\hbar\bar{\rho}\zeta\dot{\phi}-\frac{\hbar^{2}\bar{\rho}}{4m}(\partial_{x}\phi)^{2}-g(\bar{\rho}^{2}+\bar{\rho}^{2}\zeta^{2})$
(6) $\displaystyle+$ $\displaystyle
J(x)\bar{\rho}\sqrt{1-\zeta^{2}}\cos(\phi),$

in which we also assumed space homogeneity
$\partial_{x}\bar{\rho}=\partial_{x}\bar{\phi}=0$ and neglected space
variations of the population imbalance $\partial_{x}\zeta=0$ (see also, for
instance, Refs. tononi ; essler ). It is reasonable (also with the
experimental setup) to work with small values of the population imbalance,
such that $\absolutevalue{\zeta}\ll 1$. In this case we obtain the equation of
motion

$\zeta(x,t)=-\frac{\hbar\dot{\phi}(x,t)}{2g\bar{\rho}+J(x)\cos(\phi)}.$ (7)

Setting the Josephson regime $2g\bar{\rho}\gg J(x)$, we can simplify the
expression of Eq. (7) and the Lagrangian density only for the relative phase
can be written as

$\mathscr{L}=\frac{\hbar}{4g}\dot{\phi}^{2}-\frac{\hbar^{2}\bar{\rho}}{4m}(\partial_{x}\phi)^{2}+J(x)\bar{\rho}\cos(\phi).$
(8)

In the rest of the paper, we shall adopt the head-to-tail configuration, as
depicted in Fig. 1, where

$J(x)=J_{0}\,\delta(x)\;.$ (9)

with $J_{0}$ a constant tunneling coupling, $x\in[0,L]$ and $L$ the length of
the two elongated Bose-Einstein quasi-condensates. This head-to-tail
configuration is equivalent to a system where one quasi-condensate is confined
in the region $[-L,0]$ while the other quasi-condensate is confined in the
region $[0,L]$, and the tunneling barrier is located at $x=0$. This is the
same configuration considered in Ref. minguzzi . As we shall see, the Dirac
delta function gives rise to the so-called boundary sine-Gordon model fendley
. More generally, the idea of approximating a realistic finite-range tunneling
energy $J(x)$ with a zero-range one $J_{0}\delta(x)$ implies that
$J_{0}=\int_{-L}^{L}J(x)dx={\tilde{J}}(k=0)$, where ${\tilde{J}}(k)$ is the
Fourier transform of $J(x)$. This approximation is justified if the
characteristic range of $J(x)$, localized around $x=0$, is much smaller than
the characteristic lengths of the problem under consideration. Within this
approximation one can obtain the tunneling coupling given by Eq. (9) starting
from the Lagrangian density of a single bosonic field $\psi(x,t)$, with
$x\in[-L,L]$, setting $\psi(x,t)=\psi_{1}(x,t)+\psi_{2}(x,t)$ where
$\psi_{1}(x,t)$ is mainly localized on one side of a narrow potential barrier
$U(x)$ with maximum at $x=0$ and $\psi_{2}(x,t)$ is mainly localized on the
other side of $U(x)$. In this case one gets $U(x)|\psi(x,t)|^{2}\simeq-
J_{0}\delta(x)[\psi_{1}^{*}(x,t)\psi_{2}(x,t)+\psi_{1}(x,t)\psi_{2}^{*}(x,t)]$,
with $U(x)\simeq-J_{0}\delta(x)$. Alternatively, the delta-function tunneling
term of Eq. (9) can be intepreted as due to a local internal Rabi coupling.
Experimentally one can induce a Rabi coupling in a small region of a Bose-Bose
mixture.

Figure 1: Head-to-tail configuration of the system under investigation. The
tunneling of bosons occurs with strength $J_{0}$ only at $x=0$. Here $J_{0}$
has the dimensions of energy times length.

It is important to stress that the case of a uniform tunneling coupling, i.e.
$J(x)=J_{\mathrm{unif}}$, which is relevant for the experiment of Pigneur et
al. pigneur , was theoretically considered by Bouchoule bouchoule  and also by
Grisins and Mazets grisins .

##  III Quasi-particle description

We will now introduce a quasi-particle description for the phase, based on the
following mode expansion

$\phi(x,t)=\frac{1}{\sqrt{L}}\sum_{n=0}^{+\infty}q_{n}(t)\Phi_{n}(x),$ (10)

where $q_{n}(t)$ are coordinates, and $\Phi(x)$ are real eigenfunctions
satisfying
$-\hbar^{2}/\left(2m\right)\partial_{x}^{2}\Phi_{n}(x)=\epsilon_{n}\Phi_{n}(x)$
where $\epsilon_{n}=\hbar^{2}k_{n}^{2}/\left(2m\right)$ and $k_{n}=\pi n/L$,
and constituting an orthonormal basis
$\int_{0}^{L}\Phi_{n}(x)\Phi_{m}(x)\,dx=\delta_{n,m}$. This mode expansion in
Eq. (10) leads to

$\displaystyle{\cal L}=\int_{0}^{L}\mathscr{L}\,dx$ $\displaystyle=$
$\displaystyle\frac{M}{2}\sum_{n}\dot{q}_{n}^{2}-\frac{M}{2}\sum_{n}\omega_{n}^{2}q_{n}^{2}$
(11) $\displaystyle+$ $\displaystyle
J_{0}\bar{\rho}\cos{\left(\frac{1}{L}\sum_{n}q_{n}\right)}.$

Here we have defined the effective mass

$M=\frac{\hbar^{2}}{2gL},$ (12)

and the mode dispersion

$\omega_{n}=c_{s}\,k_{n},$ (13)

where

$c_{s}=\sqrt{\frac{\bar{\rho}g}{m}},$ (14)

is the speed of sound, as known from the description of superfluids, and that
$J_{0}$ and $g$ have the dimension of energy times a unit of length.

We can distinguish between the Josephson mode we are interested in $(n=0)$ and
a phonon bath (we mention phonons due to the acoustic spectrum above) made of
infinite contributions, that affect the evolution of the oscillations in the
Josephson mode. Separating the Josephson mode from bath modes, one obtains

$\displaystyle{\cal L}$ $\displaystyle=$
$\displaystyle\frac{M}{2}\dot{q}_{0}(t)^{2}+\frac{M}{2}\sum_{n=1}^{+\infty}\dot{q}_{n}(t)^{2}-\frac{M}{2}\sum_{n=1}^{+\infty}\omega_{n}^{2}q_{n}(t)^{2}$
(15) $\displaystyle+$ $\displaystyle
J_{0}\bar{\rho}\cos{\left[\frac{q_{0}(t)}{L}+\frac{1}{L}\sum_{n=1}^{+\infty}q_{n}(t)\right]}.$

This is the Lagrangian describing the Josephson mode coupled to a bath
composed of independent harmonic oscillators.

##  IV Damped dynamics

The Legendre transformation
$\mathcal{H}=\sum_{n=0}^{+\infty}\dot{q}_{n}(t)p_{n}(t)-{\cal L}$ with
$p_{n}(t)=M\dot{q}_{n}(t)$ results in the following Hamiltonian

$\mathcal{H}=\sum_{n=0}^{+\infty}\left[\frac{p_{n}^{2}}{2M}+\frac{M\omega_{n}^{2}}{2}q_{n}^{2}\right]-J_{0}\bar{\rho}\cos(\frac{1}{L}\sum_{n=0}^{+\infty}q_{n}(t)).$
(16)

Now we perform a canonical transformation with new coordinates and momenta

$\displaystyle Q_{0}(t)=q_{0}(t)+\sum_{n=1}^{+\infty}q_{n}(t),$ (17)
$\displaystyle Q_{n}(t)=q_{n}(t)\hskip 41.25648pt\mbox{for}\hskip
8.53581ptn\neq 0,$ (18) $\displaystyle P_{0}(t)=p_{0}(t),$ (19) $\displaystyle
P_{n}(t)=p_{n}(t)-p_{0}(t)\hskip 8.53581pt\mbox{for}\hskip 8.53581ptn\neq 0.$
(20)

Notice that is this is the analogue of the unitary transformation used in Ref.
minguzzi  for the quantum operators. In this way, we obtain the transformed
Hamiltonian

$\mathcal{H}=\frac{P_{0}^{2}}{2M}+\sum_{n=1}^{+\infty}\left[\frac{(P_{0}+P_{n})^{2}}{2M}+\frac{1}{2}M\omega_{n}^{2}Q_{n}^{2}\right]-J_{0}\bar{\rho}\cos(\frac{Q_{0}}{L}).$
(21)

It is worth noting that the canonical transformation introduced an intrinsic
coupling in the harmonic oscillators of the bath, between the zero-mode and
the excited modes. The Hamiltonian of Eq. (21) corresponds to the velocity-
coupling model in which the coupling is through the momentum ford . For the
sake of completeness, we stress that in the case of two parallel tubes with
uniform tunneling energy, i.e. with $J(x)=J_{\mathrm{unif}}$, only the
inclusion of anharmonic terms gives rise to a coupling between the Josephson
mode and the bath of elementary excitations.

Hamilton equations provided by the Hamiltonian in Eq. (21) are

$\displaystyle\dot{Q}_{0}(t)=\frac{P_{0}(t)}{M}+\sum_{n=1}^{+\infty}\frac{P_{0}(t)+P_{n}(t)}{M},$
(22) $\displaystyle\dot{Q}_{n}(t)=\frac{P_{0}(t)+P_{n}(t)}{M},$ (23)
$\displaystyle\dot{P}_{0}(t)=-\frac{J_{0}\bar{\rho}}{L}\sin(\frac{Q_{0}(t)}{L}),$
(24) $\displaystyle\dot{P}_{n}(t)=-M\omega_{n}^{2}Q_{n}(t).$ (25)

These equations of motion lead to

$\displaystyle Q_{n}(t)$ $\displaystyle=$
$\displaystyle\cos(\omega_{n}t)Q_{n}(0)+\frac{\sin(\omega_{n}t)}{\omega_{n}}\dot{Q}_{n}(0)$
(26) $\displaystyle-$
$\displaystyle\frac{J_{0}\bar{\rho}}{ML\omega_{n}^{2}}\biggl{[}\sin(\frac{Q_{0}(t)}{L})-\cos(\omega_{n}t)\sin(\frac{Q_{0}(0)}{L})$
$\displaystyle-$
$\displaystyle\int_{0}^{t}dt^{\prime}\cos{\left[\omega_{n}(t-t^{\prime})\right]}\cos(\frac{Q_{0}(t^{\prime})}{L})\frac{\dot{Q}_{0}(t^{\prime})}{L}\biggl{]}.$

Details for the derivation of Eq. (26) are summarized in Appendix A. One can
easily show that the Josephson mode is related to the relative phase at $x=0$
as

$\phi_{0}(t)=\phi(x=0,t)=\frac{Q_{0}(t)}{L}.$ (27)

Eventually, we reach the equation of motion for a damped harmonic oscillator
with respect to the relative phase

$\ddot{\phi}_{0}(t)+\int_{0}^{t}dt^{\prime}\gamma(t-t^{\prime})\dot{\phi}_{0}(t^{\prime})+\Omega_{0}^{2}\sin(\phi_{0}(t))=\xi_{\phi}(t),$
(28)

in which the damping is ruled by the damping kernel

$\gamma(t-t^{\prime})=\Omega_{0}^{2}\sum_{n=1}^{+\infty}\cos(\omega_{n}(t-t^{\prime}))\cos(\phi_{0}(t^{\prime})),$
(29)

proportional to the Josephson frequency defined as

$\Omega_{0}=\sqrt{\frac{\bar{\rho}J_{0}}{ML^{2}}}.$ (30)

Now let us focus on the case of a small relative phase
$\absolutevalue{\phi_{0}(t)}\ll 1$. In this case
$\cos(\phi_{0}(t^{\prime}))\sim 1$ holds. The expression of Eq. (29) can be
calculated explicitly by moving to the continuum limit as

$\gamma(t-t^{\prime})=\gamma_{0}\,\delta(t-t^{\prime}),$ (31)

where

$\gamma_{0}=\frac{\bar{\rho}J_{0}}{MLc_{s}},$ (32)

is the damping constant that rules the relaxation in time of the phase. Using
the damping constant, one can write Eq. (28) as

$\ddot{\phi}_{0}(t)+\gamma_{0}\dot{\phi}_{0}(t)+\Omega_{0}^{2}\phi_{0}(t)=\xi_{\phi}(t).$
(33)

This is the main equation in the paper. On the left-hand side of this
equation, it involves the Josephson frequency $\Omega_{0}$, given by Eq. (30),
and the damping constant, given by Eq. (32) with the speed of sound $c_{s}$ of
Eq. (14).

The right-hand side of Eq. (33) is a noise term denoted by $\xi_{\phi}(t)$,
composed of the infinite contributions of the bath modes and dependent only on
the initial conditions of coordinates and momenta

$\begin{split}\xi_{\phi}(t)=&-\sum_{n=1}^{+\infty}\left[\omega_{n}^{2}\cos(\omega_{n}t)\frac{Q_{n}(0)}{L}+\omega_{n}\sin(\omega_{n}t)\frac{\dot{Q}_{n}(0)}{L}\right]\\\
&-\gamma_{0}\delta(t)\sin{\left(\phi_{0}(0)\right)}.\end{split}$ (34)

The last term has the nature of a transient, and in order to avoid its
complications, we will get rid of the transient term by setting the initial
condition $\phi_{0}(0)=0$.

###  IV.1 Deterministic dynamics of the Josephson mode

As an homogeneous solution of Eq. (33) without the noise, one gets

$\phi_{0}(t)=e^{-\gamma_{Q}\Omega_{0}t}\frac{\dot{\phi}(0)}{\Omega_{0}\gamma_{J}}\sin(\gamma_{J}\Omega_{0}t),$
(35)

where we defined

$\gamma_{Q}={\gamma_{0}\over 2\Omega_{0}},$ (36)

and

$\gamma_{J}=\sqrt{1-\gamma_{Q}^{2}}.$ (37)

We can confirm that the dimensionless damping constant $\gamma_{Q}$
corresponds to $\sqrt{E_{J}/K}$ in Ref. minguzzi  in the weak-coupling limit
where $K\sim 1/\sqrt{g}$ denotes the Luttinger parameter in Luttinger liquid
theory minguzzi . This indicates that our current analysis under the quasi-
particle description recovers the analysis within Luttinger liquid theory in
the weak-coupling limit $K\sim 1/\sqrt{g}$. In the following, only the
underdamped regime $\gamma_{Q}<1$ is analyzed for simplicity. It reflects the
most similar behaviour with respect to the experimental observations of the
relative phase and population imbalance pigneur . The other two regimes,
overdamped and critical damped, are respectively obtained in the case of
$\gamma_{Q}>1$ and $\gamma_{Q}=1$.

Figure 2: Josephson mode dynamics without the effect of the noise
($\xi_{\phi}(t)=0$). Time evolution of the relative phase $\phi(t)$ (panel
(a)) and of the the zero-mode $\zeta_{0}(t)$ of the population imbalance
(panel (b)). Here we consider the underdamped case with $\gamma_{Q}<1$.
Initial conditions $\phi_{0}(0)=0$ and $\dot{\phi}_{0}(0)/\Omega_{0}=1$;
$\zeta_{0}(0)=0.2$ and $\dot{\zeta}_{0}(0)/\Omega_{0}=0$.

In panel (a) of Fig. 2, we plot the dynamics of the relative phase
$\phi_{0}(t)$, without the effect of noise $\xi_{\phi}(t)$ and setting the
initial condition $\dot{\phi}_{0}(0)/\Omega_{0}=1$ for simplicity. In panel
(b) of Fig. 2, we plot instead the corresponding population imbalance
$\zeta_{0}(t)$, with initial conditions $\zeta_{0}(0)=0.2$ and
$\dot{\zeta}_{0}(0)/\Omega_{0}=0$. The initial conditions for
$\dot{\zeta}(0)/\Omega_{0}$ is dependent on the choice of $\phi_{0}(0)$
through Eq. (24). On the other hand, the choice of $\zeta_{0}(0)$ is necessary
in order to start from a configuration out of equilibrium, in which the number
of atoms in the two wells is not balanced. It illustrates quite intuitive
behaviours. As one increases the damping coefficient $\gamma_{0}$, the
amplitude od oscillation is suppressed and both of $\phi_{0}(t)$ and
$\zeta_{0}(t)$ vanish for $t\to\infty$.

The results in Fig. 2 is, however, missing the relevant role of the noise,
which produces fluctuations around the mean-field behaviour shown in the
figure. In general, these fluctuations have both quantum and thermal
components. We investigate these effects of the noise in the next Sections.

###  IV.2 Dynamics of the Josephson mode in the presence of noise

We refer to $\xi_{\phi}(t)$ as a quasi-stochastic noise due to the analogy of
Eq. (33) with a generalized Langevin equation. The stochastic nature of this
term, however, lies only in the initial conditions $Q_{n}(0)$ and
$\dot{Q}_{n}(0)$, then $\xi_{\phi}(t)$ is not a random quantity in time.
Thanks to this property, we do not need special tools like stochastic calculus
in order to integrate Eq. (33) in time, but we will find the particular
solution using the method of variation of parameters. The solution takes the
form of

$\phi_{0}(t)=e^{-\gamma_{Q}\Omega_{0}t}\frac{\sin(\gamma_{J}\Omega_{0}t)}{\gamma_{J}}+\int_{0}^{t}dt^{\prime}\chi(t-t^{\prime})\xi_{\phi}(t^{\prime}),$
(38)

where we have defined the following retarded function

$\chi(t-t^{\prime})=\frac{2}{\omega_{D}}e^{-\gamma_{Q}\Omega_{0}(t-t^{\prime})}\sin{\left[\frac{\omega_{D}}{2}(t-t^{\prime})\right]}\theta(t-t^{\prime}),$
(39)

with $\theta(t)$ the Heaviside step function and the oscillation frequency

$\omega_{D}=\sqrt{4\Omega_{0}^{2}-\gamma_{0}^{2}}.$ (40)

##  V The population imbalance

Using Eq. (7) in the Josephson regime, one can make a quasi-particle
description also for the population imbalance $\zeta(x,t)$ as

$\displaystyle\zeta(x,t)$
$\displaystyle=\frac{\sqrt{L}}{\hbar\bar{\rho}}\sum_{n=0}^{+\infty}p_{n}(t)\Phi_{n}(x)$
(41)
$\displaystyle=\frac{P_{0}(t)}{\hbar\bar{\rho}}+\frac{\sqrt{L}}{\hbar\bar{\rho}}\sum_{n=1}^{+\infty}\left[P_{0}(t)+P_{n}(t)\right]\Phi_{n}(x).$

We can identify the Josephson mode as the first term

$\zeta_{0}(t)=\frac{P_{0}(t)}{\hbar\bar{\rho}}.$ (42)

Equations of motion for the momenta $P_{0}(t)$ and $P_{n}(t)$ read

$\displaystyle\ddot{P}_{0}(t)+\Omega_{0}^{2}P_{0}(t)=\Omega_{0}^{2}\sum_{n=1}^{+\infty}\frac{\ddot{P}_{n}(t)}{\omega_{n}^{2}},$
(43)
$\displaystyle\ddot{P}_{n}(t)+\omega_{n}^{2}P_{n}(t)=-\omega_{n}^{2}P_{0}(t).$
(44)

They give the damped behaviour of the Josephson mode

$\ddot{\zeta}_{0}(t)+\gamma_{0}\dot{\zeta}_{0}(t)+\Omega_{0}^{2}\zeta_{0}(t)=\xi_{\zeta}(t).$
(45)

Here we have defined the noise term for the population imbalance as

$\begin{split}\xi_{\zeta}(t)=-\frac{\Omega_{0}^{2}}{\hbar\bar{\rho}}\sum_{n=1}^{+\infty}\left[\cos(\omega_{n}t)P_{n}(0)+\frac{\sin(\omega_{n}t)}{\omega_{n}}\dot{P}_{n}(0)\right],\end{split}$
(46)

where we have set the initial condition $P_{0}(0)=0$. The meaning of this
choice will be stressed in the next section. The time evolution of
$\zeta_{0}(t)$, without the effect of the noise $\xi_{\zeta}(t)$, is displayed
in Fig. 2(b).

##  VI Quantum and thermal properties of the noise

Interesting quantities regarding the noise, which allow us to investigate the
quantum and thermal fluctuations of $\phi_{0}(t)$ due to the presence of
$\xi_{\phi}(t)$, are its bath average $\expectationvalue{\xi_{\phi}(t)}$ and,
in particular, its correlation function
$\expectationvalue{\xi_{\phi}(t)\xi_{\phi}(t^{\prime})}$. We will perform a
calculation already done in Ref. ingold  in the context of a Caldeira-Leggett
model caldeira . In order to evaluate the average over the environment, we
have to identify the bath Hamiltonian. From Eq. (21), denoting the bath part
by $\mathcal{H}_{\rm B}$, we have

$\mathcal{H}_{\rm
B}=\sum_{n=1}^{+\infty}\left[\frac{(P_{0}+P_{n})^{2}}{2M}+\frac{M\omega_{n}^{2}}{2}Q_{n}^{2}\right].$
(47)

As we previously pointed out, this Hamiltonian is composed of infinite
harmonic oscillators which are intrinsically coupled to the system. However,
it is reasonable to assume that the size of the phonon bath is huge such that
the system does not affect it. For this reason one can impose $P_{0}(0)=0$,
namely the system is completely decoupled from the bath at $t=0$. Hence, in
this model, we deal with a bath of independent harmonic oscillators ford  at
initial time $t=0$ as

$\mathcal{H}_{\rm
B}=\sum_{n=1}^{+\infty}\left[\frac{P_{n}^{2}}{2M}+\frac{M\omega_{n}^{2}}{2}Q_{n}^{2}\right]$
(48)

which is used to evaluate the ensemble average.

In order to evaluate ensemble averages, we consider $Q_{n}(0)$ and $P_{n}(0)$
to be quantum operators. By means of the Hamiltonian in Eq. (48), we can adopt
the annihilation and creation operators $a_{n}$ and $a_{n}^{\dagger}$ as

$\displaystyle
Q_{n}(0)=\sqrt{\frac{\hbar}{2M\omega_{n}}}(a_{n}+a_{n}^{\dagger}),$ (49)
$\displaystyle P_{n}(0)=-i\sqrt{\frac{\hbar
M\omega_{n}}{2}}(a_{n}-a_{n}^{\dagger}),$ (50)

and then evaluate the average, in the case of $\xi_{\phi}(t)$ for instance, as

$\expectationvalue{\xi_{\phi}(t)}=\frac{\Tr\left[e^{-\beta\mathcal{H}_{\rm
B}}\xi_{\phi}(t)\right]}{\Tr\left[e^{-\beta\mathcal{H}_{\rm B}}\right]},$ (51)

where $\beta=1/(k_{\rm B}T)$ with $k_{\rm B}$ the Boltzmann constant and $T$
the absolute temperature of the bath of oscillators.

One readily obtains the vanishing averages for both of Eqs. (34) and (46)

$\expectationvalue{\xi_{\phi}(t)}=\expectationvalue{\xi_{\zeta}(t)}=0,$ (52)

while we find the two-point correlation functions as

$\displaystyle\expectationvalue{\xi_{\phi}(t)\xi_{\phi}(t^{\prime})}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{+\infty}\frac{\hbar\omega_{n}^{3}}{2ML^{2}}\Bigg{[}\coth(\frac{\beta\hbar\omega_{n}}{2})\cos{\left[\omega_{n}(t-t^{\prime})\right]}$
(53) $\displaystyle-$ $\displaystyle
i\sin{\left[\omega_{n}(t-t^{\prime})\right]}\Bigg{]},$
$\displaystyle\expectationvalue{\xi_{\zeta}(t)\xi_{\zeta}(t^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{M\Omega_{0}^{4}}{\hbar^{2}\bar{\rho}^{2}}\sum_{n=1}^{+\infty}\frac{\hbar\omega_{n}}{2}\Bigg{[}\coth(\frac{\beta\hbar\omega_{n}}{2})\cos{\left[\omega_{n}(t-t^{\prime})\right]}$
(54) $\displaystyle-$ $\displaystyle
i\sin{\left[\omega_{n}(t-t^{\prime})\right]}\Bigg{]}.$

In the high-temperature limit $k_{\rm B}T\gg\hbar\omega_{n}$, one can
approximate

$\coth(\frac{\beta\hbar\omega_{n}}{2})\simeq\frac{2}{\beta\hbar\omega_{n}},$
(55)

resulting in the following correlation functions

$\displaystyle\expectationvalue{\xi_{\phi}(t)\xi_{\phi}(t^{\prime})}=-\frac{\gamma_{0}}{M\Omega_{0}^{2}}k_{\rm
B}T\frac{d^{2}}{dt^{2}}\delta(t-t^{\prime}),$ (56)
$\displaystyle\expectationvalue{\xi_{\zeta}(t)\xi_{\zeta}(t^{\prime})}=\frac{\gamma_{0}M\Omega_{0}^{2}}{\hbar^{2}\bar{\rho}^{2}}k_{\rm
B}T\delta(t-t^{\prime}).$ (57)

It is worth noting the fact that $\xi_{\zeta}(t)$ is a delta-correlated noise
at high temperature likewise the classical fluctuation-dissipation relation,
while $\xi_{\phi}(t)$ is slightly different from a white noise.

In the case of a generic temperature, one has to take into account that Eqs.
(53) and (54) are complex valued quantities. In particular it is convenient to
separate the real part from the purely imaginary one, by rewriting the
correlators as

$\expectationvalue{\xi_{\phi}(t)\xi_{\phi}(t^{\prime})}=\frac{1}{2}\expectationvalue{\anticommutator{\xi_{\phi}(t)}{\xi_{\phi}(t^{\prime})}}+\frac{1}{2}\expectationvalue{\commutator{\xi_{\phi}(t)}{\xi_{\phi}(t^{\prime})}}.$
(58)

The first term, which involves the anti-commutator
$\anticommutator{\cdot}{\cdot}$ is related to the real part of Eq. (53), while
the second one, which involves the commutator $\commutator{\cdot}{\cdot}$,
describes the purely imaginary term of the noise. The anti-symmetric parts
are, however, found to result in no contribution to the variance of phase or
population imbalance because they are odd functions of $\omega$. The Fourier
transformed correlators of the symmetric part are given by

$\displaystyle\int_{-\infty}^{+\infty}dt\expectationvalue{\anticommutator{\xi_{\phi}(t)}{\xi_{\phi}(0)}}e^{-i\omega
t}=\frac{2\gamma_{0}\hbar\omega^{3}}{ML^{2}\Omega_{0}^{2}}\coth(\frac{\beta\hbar\omega}{2}),$
(59)
$\displaystyle\int_{-\infty}^{+\infty}dt\expectationvalue{\anticommutator{\xi_{\zeta}(t)}{\xi_{\zeta}(0)}}e^{-i\omega
t}=\frac{2M\Omega_{0}^{2}\gamma_{0}\hbar\omega}{\hbar^{2}\bar{\rho}^{2}}\coth(\frac{\beta\hbar\omega}{2}).$
(60)

##  VII Fluctuations

In order to understand the effect of the noise in the dynamics of
$\phi_{0}(t)$ and $\zeta_{0}(t)$, we focus on the calculation of the variance
$\Delta\phi_{0}(t)^{2}=\expectationvalue{\phi_{0}(t)^{2}}-\expectationvalue{\phi_{0}(t)}^{2}$
and
$\Delta\zeta_{0}^{2}=\expectationvalue{\zeta_{0}(t)^{2}}-\expectationvalue{\zeta_{0}(t)}^{2}$.
It is obvious from Eq. (38) that they will depend only on the correlators of
the noise.

Figure 3: Relative phase variance $\Delta\phi_{0}(t)^{2}$ (solid curves) and
population imbalance variance $\Delta\zeta_{0}(t)^{2}$ (dashed curves) of the
Josephson mode as a function of time $t$. The curves correspond to three
different underdamped regimes at high temperature, where $k_{\rm
B}T/(\hbar\Omega_{0})\gg 1$. The variances are nomalized by their asymptonic
values, given by Eqs. (63) and (64). The parameter $\gamma_{Q}$ is defined in
Eq. (36). The two upper solid and dashed curves correspond to
$\gamma_{Q}=1/3$, the two intermediate ones to $\gamma_{Q}=1/5$, and the two
lower ones to $\gamma_{Q}=1/10$.

###  VII.1 High temperature

Let us start from the relative phase fluctuations at high temperature with
Eqs. (56) and (57). Eqs. (38), (56), and (57) give

$\displaystyle\Delta\phi_{0}(t)^{2}$ $\displaystyle=$
$\displaystyle\Delta\phi_{0}(\infty)^{2}\,\Bigg{[}1+e^{-2\gamma_{Q}\Omega_{0}t}\biggl{[}\frac{\gamma_{Q}^{2}-\gamma_{J}^{2}}{\gamma_{J}^{2}}$
(61) $\displaystyle-$
$\displaystyle{\gamma_{Q}^{2}\over\gamma_{J}^{2}}\cos(2\gamma_{J}\Omega_{0}t)-{\gamma_{Q}\over\gamma_{J}}\sin(2\gamma_{J}\Omega_{0}t)\biggl{]}\Bigg{]},$
$\displaystyle\Delta\zeta_{0}(t)^{2}$ $\displaystyle=$
$\displaystyle\Delta\zeta_{0}(\infty)^{2}\Bigg{[}1-e^{-2\gamma_{Q}\Omega_{0}t}\biggl{[}\frac{1}{\gamma_{J}^{2}}$
(62) $\displaystyle-$
$\displaystyle{\gamma_{Q}^{2}\over\gamma_{J}^{2}}\cos(2\gamma_{J}\Omega_{0}t)+{\gamma_{Q}\over\gamma_{J}}\sin(2\gamma_{J}\Omega_{0}t)\biggl{]}\Bigg{]},$

where the asymptotic values of these variances for $t\to+\infty$ are given by

$\displaystyle\Delta\phi_{0}(\infty)$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega_{0}}\sqrt{\frac{k_{\rm B}T}{2ML^{2}}},$ (63)
$\displaystyle\Delta\zeta_{0}(\infty)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{Mk_{\rm B}T}{2\hbar^{2}\bar{\rho}^{2}}}.$ (64)

Notice that Eq. (62) is equivalent to the one in Ref. minguzzi .

By using Eqs. (61) and (62) we plot the time evolution of the fluctuations in
the relative phase and population imbalance in Fig. 3. They exhibit strong
similarities in the behaviour of these two quantities, which are ruled by the
damping parameter $\gamma_{Q}$. One can see that curves tend to overlap as
much as one lowers $\gamma_{Q}$.

###  VII.2 Generic temperature

Let us focus on the generic temperature case to see effects of quantum
fluctuations in addition to thermal ones. We already mentioned the fact that
the noise correlators in the most general case are complex valued quantities.
Then we use Eq. (58) for both
$\expectationvalue{\xi_{\phi}(t)\xi_{\phi}(t^{\prime})}$ and
$\expectationvalue{\xi_{\zeta}(t)\xi_{\zeta}(t^{\prime})}$ and we obtain

$\displaystyle\Delta\phi_{0}(t)^{2}=\frac{1}{2}\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}\expectationvalue{\anticommutator{\xi_{\phi}(s)}{\xi_{\phi}(s^{\prime})}}\chi(t-s)\chi(t-s^{\prime}),$
(65)
$\displaystyle\Delta\zeta_{0}(t)^{2}=\frac{1}{2}\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}\expectationvalue{\anticommutator{\xi_{\zeta}(s)}{\xi_{\zeta}(s^{\prime})}}\chi(t-s)\chi(t-s^{\prime}).$
(66)

Eq. (59) provides

$\displaystyle\Delta\phi_{0}(t)^{2}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{0}}{M\Omega_{0}^{2}}\int_{-\infty}^{+\infty}d\omega\hbar\omega^{3}\coth(\frac{\beta\hbar\omega}{2})$
(67) $\displaystyle\times$
$\displaystyle\int_{0}^{t}dt_{1}\chi(t_{1})e^{-i\omega
t_{1}}\int_{0}^{t}dt_{2}\chi(t_{2})e^{i\omega t_{2}},$

where we changed the time variables $t_{1}=t-s$ and $t_{2}=t-s^{\prime}$.
Introducing a new adimensional frequency $\tilde{\omega}=\omega/\Omega_{0}$,
we obtain an expression for the variance of the phase

$\Delta\phi_{0}(t)^{2}=\Gamma\frac{\gamma_{Q}}{\gamma_{J}^{2}}\int_{-\infty}^{+\infty}d\tilde{\omega}\tilde{\omega}^{3}\coth(\frac{\beta\hbar\Omega_{0}\tilde{\omega}}{2})G_{\tilde{\omega}}(t)G_{-\tilde{\omega}}(t),$
(68)

where

$\Gamma={2\hbar\over M\Omega_{0}},$ (69)

Figure 4: (a) Time evolution of the relative phase variance
$\Delta\phi_{0}(t)^{2}$ of the Josephson mode. (b) Time evolution of the
population imbalance variance $\Delta\zeta_{0}(t)^{2}$ of the Josephson mode.
Here we set three different values for $\gamma_{Q}$ and $k_{\rm
B}T/(\hbar\Omega_{0})=10^{-2}$ (low-temperature regime). The normalization
factors $\Gamma$ and $\eta$ are defined in Eqs. (69) and (72). Note that the
results are crucially dependent on an ultraviolet cutoff
${\tilde{\omega}}_{\rm max}$. Here we use ${\tilde{\omega}}_{\rm max}=10^{4}$
(see the main text for details).

and

$\displaystyle G_{\tilde{\omega}}(t)$ $\displaystyle=$
$\displaystyle\frac{e^{-(i\tilde{\omega}+\gamma_{Q})\Omega_{0}t}}{(i\tilde{\omega}+\gamma_{Q})^{2}-\gamma_{J}^{2}}\Big{[}\gamma_{J}e^{(i\tilde{\omega}+\gamma_{Q})(\Omega_{0}t)}$
(70) $\displaystyle-$
$\displaystyle\gamma_{J}\cos(\gamma_{J}\Omega_{0}t)-(i\tilde{\omega}+\gamma_{Q})\sin(\gamma_{J}\Omega_{0}t)\Big{]}.$

What we have just seen can be identically translated into the calculation of
Eq. (66) as

$\Delta\zeta_{0}(t)^{2}=\eta\frac{\gamma_{Q}}{\gamma_{J}^{2}}\int_{-\infty}^{+\infty}d\tilde{\omega}\tilde{\omega}\coth(\frac{\beta\hbar\Omega_{0}\tilde{\omega}}{2})G_{\tilde{\omega}}(t)G_{-\tilde{\omega}}(t),$
(71)

where

$\eta={M\Omega_{0}\over\hbar\bar{\rho}^{2}}.$ (72)

Both of Eqs. (68) and (71) are computed numerically, but it is worth noting
that they have an ultraviolet divergence. For this reason it is necessary to
set a cutoff that, taking into account the phonon dispersion relation we found
previously, can be seen as ${\tilde{\omega}}_{\rm max}=c_{s}\pi
N/(\Omega_{0}L)$. Using typical parameters like $N=10^{3}$, $L=10^{-6}{\rm
m}$, $\Omega_{0}=10^{2}{\rm s^{-1}}$, $c_{s}=10^{-3}{\rm ms^{-1}}$ pigneur ,
we find ${\tilde{\omega}}_{\rm max}=10^{4}$.

In Fig. 4, we report the time evolution of the fluctuations for three values
of $\gamma_{Q}$ in the underdamped regime ($\gamma_{Q}<1)$. Panel (a) in Fig.
4 clearly shows that the variance $\Delta\phi_{0}(t)^{2}$, related to the
relative phase of the Josephson mode, oscillates until it reaches an
asymptotic value in the long-time limit, as shown in the inset. The huge
initial oscillations are crucially dependent on the cutoff
$\tilde{\omega}_{\rm max}$, and the asymptotic value strongly depends on the
choice of $\gamma_{Q}$. It is worth noting that decreasing the value of
$\gamma_{Q}$ leads to a lower asymptotic value. The variance
$\Delta\zeta_{0}(t)^{2}$ of the population imbalance of the Josephson mode in
panel (b) of Fig. 4 also oscillates in time, reaching asymptotically a finite
value. Here the asymptotic value also strongly depends on the damping
parameter $\gamma_{Q}$ and decreasing $\gamma_{Q}$ leads to a larger
asymptotic value.

Figure 5: Variance $\Delta\zeta_{0}(t)^{2}$ of the Josephson population
imbalance as a function of time $t$ for three values of the temperature $T$ of
the bath of phonons. The damping coefficient $\gamma_{Q}$ is set to be
$\gamma_{Q}=1/10$. $\Omega_{0}$ is the Josephson frequency of Eq. (30) while
$\eta$ is given by Eq. (72).

For the sake of completeness, in Fig. 5, we plot $\Delta\zeta_{0}(t)^{2}$ for
three values of the rescaled temperature $k_{\rm B}T/(\hbar\Omega_{0})$ of the
phonon bath. The figure shows that the asymptotic value
$\Delta\zeta_{0}(\infty)^{2}$ grows by increasing the temperature. We have
verified that the same happens with $\Delta\phi_{0}(\infty)^{2}$. These
results are fully consistent with the ones obtained in Ref. grabert  for the
quadratic fluctuations of the damped harmonic oscillator in thermal
equilibrium.

##  VIII Conclusions

We have found that in elongated Josephson junctions the damped dynamics of the
Josephson mode is ruled by the damping constant $\gamma_{0}$, given by Eq.
(32). Moreover, the Josephson dynamics is strictly dependent on the
interaction between the Josephson mode and the quantum-thermal bath of
phonons. We have studied the phase fluctuation $\Delta\phi_{0}(t)$ and the
population-imbalance fluctuation $\Delta\zeta_{0}(t)$ of the Josephson mode.
In our work, as well as in Ref. minguzzi , both noise and damping have the
same origin: the intrinsic coupling to the phonon bath. The damping of
correlations is caused by the friction while nonzero variances are ascribed to
the noise. In the high-temperature regime, where the thermal fluctuations
dominate, we have derived analytic expressions for the Josephson fluctuations.
The generic temperature case has also been taken into account. Here quantum
fluctuations play an important role and we rely on numerical calculations to
determine the time evolution of the Josephson fluctuations, which exhibit
their thermalization to constant values after a transient characterized by
oscillating dynamics.

## Acknowledgments

The authors acknowledge Anna Minguzzi and Juan Polo for useful
e-clarifications. LS thanks Andrea Tononi and Flavio Toigo for enlightening
discussions. KF is supported by a PhD fellowship of the Fondazione Cariparo.

##  Appendix A Derivation of Eq. (26)

In this section, we briefly explain how to obtain the solution of the bath
coordinates in Eq. (26). Taking the derivative of the first two equations of
motion in Eqs. (22) and (23) with respect to time, we reach two differential
equations for the Josephson mode and the excited modes as

$\displaystyle\ddot{Q}_{n}(t)+\omega_{n}^{2}Q_{n}(t)=-\frac{J_{0}\bar{\rho}}{ML}\sin(\frac{Q_{0}(t)}{L}),$
(73)
$\displaystyle\ddot{Q}_{0}(t)+\frac{J_{0}\bar{\rho}}{ML}\sin(\frac{Q_{0}(t)}{L})=\sum_{n=1}^{+\infty}\ddot{Q}_{n}(t).$
(74)

In this way we can start by solving the first equation, finding $Q_{n}(t)$,
and then we exploit what we found in order to study the second one. Eq. (73)
can be solved by taking the Laplace transformation on each member, and we end
up with

$\begin{split}Q_{n}(t)&=\cos(\omega_{n}t)Q_{n}(0)+\frac{\sin(\omega_{n}t)}{\omega_{n}}\dot{Q}_{n}(0)+\\\
&-\frac{J_{0}\bar{\rho}}{ML\omega_{n}}\int_{0}^{t}dt^{\prime}\sin{\left[\omega_{n}(t-t^{\prime})\right]}\sin(\frac{Q_{0}(t^{\prime})}{L}).\end{split}$
(75)

Integrating by parts the last term, we obtain Eq. (26).

## References

  * (1) A. O. Caldeira and A. J. Leggett, Phys. Rev. Lett. 46, 211 (1981).
  * (2) G. W. Ford, J. T. Lewis, and R. F. O’Connell, Phys. Rev. A 37, 4419 (1988).
  * (3) A. Schmidt, Phys. Rev. Lett. 51, 1506 (1983).
  * (4) S. A. Bulgadaev, JETP Lett. 39, 315 (1984).
  * (5) E. D. Torre, E. Demler, T. Giamarchi, and E. Altman, Nat. Phys. 6, 806 (2010).
  * (6) R. H. Koch, D. J. Van Harlingen, and J. Clarke, Phys. Rev. Lett. 45, 2132 (1980).
  * (7) R. H. Koch, D. J. Van Harlingen, and J. Clarke, Phys. Rev. B 26, 74 (1982).
  * (8) M. Pigneur, T. Berrada, M. Bonneau, T. Schumm, E. Demler, and J. Schmiedmayer, Phys. Rev. Lett. 120, 173601 (2017).
  * (9) J. Polo, V. Ahufinger, F. W. J. Hekking, and A. Minguzzi, Phys. Rev. Lett. 2, 090404 (2018).
  * (10) A. Tononi, F. Toigo, S. Wimberger, A. Cappellaro, and L. Salasnich, New J. Phys. 22, 073020 (2020).
  * (11) Y. D van Nieuwkerk and F. H. L. Essler, SciPost Phys. 9, 025 (2020).
  * (12) P. Fendley, F. Lesage, and H. Saleur, J. Stat. Phys. 85, 211 (1996).
  * (13) I. Bouchoule, Eur. Phys. J. D 35, 147 (2005).
  * (14) P. Grisins and I. E. Mazets, Phys. Rev. A 87, 013629 (2013).
  * (15) J. Javanainen, Phys. Rev. Lett. 57, 3164 (1986).
  * (16) B. D. Josephson, Phys. Lett. 1, 251 (1962).
  * (17) A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 (1997).
  * (18) M. Albiez, R. Gati, J. Folling, S. Hunsmann, M. Cristiani, and M.K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
  * (19) H. A. Gersch and G. C. Knollman, Phys. Rev. 129, 959 (1963).
  * (20) G. Ingold, A. Buchleitner, and K. Hornberger, Coherent Evolution in Noisy Environments (Springer, 2002).
  * (21) H. Grabert, U, Weiss, and P. Talkner, Z. Phys. B 55, 87 (1984).