# Exploring the early Bitcoin Mesoscale

Nicolò Vallarano Blockchain & Distributed Ledger Technologies Group,
Department of Informatics, University of Zurich, Switzerland
UZH Blockchain Center, University of Zurich, Switzerland Tiziano Squartini IMT
School for Advanced Studies, Piazza San Francesco 19, 55100 - Lucca (Italy)
Institute for Advanced Study (IAS), University of Amsterdam, Oude Turfmarkt
145, 1012 GC - Amsterdam (The Netherlands) Claudio J. Tessone Blockchain &
Distributed Ledger Technologies Group, Department of Informatics, University
of Zurich, Switzerland
UZH Blockchain Center, University of Zurich, Switzerland

###### Abstract

The open availability of the entire history of the Bitcoin transactions opens
up the possibility to study this system at an unprecedented level of detail.
This contribution is devoted to the analysis of the _mesoscale_ structural
properties of the Bitcoin User Network (BUN), across its entire history (i.e.
from 2009 to 2017). What emerges from our analysis is that the BUN is
characterized by a _core-periphery_ structure a deeper analysis of which
reveals a certain degree of _bow-tieness_ (i.e. the presence of a Strongly-
Connected Component, an IN- and an OUT-component together with some tendrils
attached to the IN-component). Interestingly, the evolution of the BUN
structural organization experiences fluctuations that seem to be correlated
with the presence of _bubbles_ , i.e. periods of price surge and decline
observed throughout the entire Bitcoin history: our results, thus, further
confirm the interplay between structural quantities and price movements
observed in previous analyses.

###### keywords:

Bitcoin. Networks. Motifs. Transactions.

##  1 Introduction

Introduced in 2008 by Satoshi Nakamoto with the release of a white paper [1],
Bitcoin is the first and most widely adopted _cryptocurrency_. Loosely
speaking, it consists of a decentralised peer-to-peer network to which users
connect to exchange native tokens (i.e. the _bitcoins_). After having been
validated by the so-called _miners_ \- according to the consensus rules that
are part of the Bitcoin protocol [2, 3] \- each transaction is included in a
replicated database, i.e. the _blockchain_. The cryptographical protocols
Bitcoin rests upon aim at preventing the possibility for the same digital
token to be spent more than once, in absence of a central, third party that
guarantees the validity of the transactions themselves [1, 4]: remarkably, the
transaction-verification mechanism Bitcoin relies on allows its entire
transaction history to be openly accessible - a feature that, in turn, allows
it to be analyzable in the preferred representation.

The Bitcoin structural properties have been only recently started to be
investigated: in [5], the authors consider the Bitcoin user network at the
_macroscale_ , in order to check for its small-worldness; in [6], the authors
explore the evolution of the _local_ properties of four different
representations of Bitcoin (i.e. both the _user_ and the _address_ network, at
both the _daily_ and _weekly_ time scale) with the aim of investigating their
relationship with the price movements; in [7], the authors analyze the so-
called Bitcoin Lightning Network (BLN), highlighting the increasingly
centralized character of such a network structure.

With the present paper, we aim at contributing to this stream of research by
studying the structure of the Bitcoin User Network (BUN) at the _mesoscale_.

##  2 Data

###  2.1 Data Subsection

As said before, Bitcoin relies on a decentralised public ledger, the
blockchain, that records all transactions of bitcoins among users. A
transaction is nothing else than a set of input and output addresses: the
output addresses that are said ‘unspent’, i.e. not yet recorded on the ledger
as input addresses, can be claimed, and therefore spent, only by the owner of
the corresponding cryptographic key. This is the reason why one speaks of
_pseudonimity_ : an observer of the blockchain can see all unspent addresses
but cannot link them to the actual owners. Techniques exist, however, to infer
the identity of users (hereby, this term will indicate ‘groups of addresses’):
they rest upon the the so-called _heuristics_ , i.e. sets of rules taking
advantage of the implementation of the Bitcoin protocol.

Bubble | Start | End | Days
---|---|---|---
1 | 25 May 2012 | 18 Aug 2012 | 84
2 | 3 Jan 2013 | 11 Apr 2013 | 98
3 | 7 Oct 2013 | 23 Nov 2013 | 47
4 | 31 Mar 2017 | 18 Dec 2017 | 155
Table 1: The four bitcoin bubbles as detected in [8].

#### Bitcoin Address Network (BAN).

The Bitcoin Address Network is the simplest network that can be constructed
from the blockchain records: it is a directed, weighted graph whose nodes
represent addresses. The direction and the weight of links are provided by the
input-output relationships defining the transactions recorded on the
blockchain. The only free parameter is represented by the temporal window
chosen for the data aggregation. In this paper, we chose to aggregate these
data at a daily time scale, i.e. the shortest scale that still guarantees that
the resulting network is connected [1].

#### Bitcoin User Network (BUN).

Since the same owner may control several addresses, one can derive a network
of users whose nodes are _clusters of addresses_. These clusters are derived
by implementing different heuristics: let us now provide a brief description
of the ones that have been employed here and that have been have derived from
the state-of-the-art literature [9, 10, 11, 12].

The first heuristics is the so-called _multi-input heuristics_ : it is based
on the assumption that two (or more) addresses that are part of the input of
the same transaction are controlled by the same user. The key idea behind this
heuristics is that the private keys of all addresses must be accessible to the
creator of a transaction, in order to produce it. It is generally believed to
be the safest one for clustering addresses.

The second heuristics is the so-called _change-address identification
heuristics_ and is based upon the observation that transaction outputs must be
fully spent upon re-utilisation; hence, the transaction creator usually
controls also one of the output addresses. More specifically, we assume that
if an output address is new and the amount transferred to it is lower than all
the inputs, then it must belong to the input user.

Whenever the Bitcoin User Network is mentioned in the paper, we refer to the
representation obtained by clustering the nodes of the BAN according to a
combination of the two heuristics above. Naturally, users can employ different
wallets that are not necessarily linked together by transactions: as a
consequence, the user network we obtain should not be considered as a perfect
representation of the actual network of users but, rather, an attempt enabling
us to group addresses while minimising the presence of false positives. Once
the addresses are grouped into ‘users’, we build the network as follows: any
two nodes $i$ and $j$ are connected via a directed edge from $i$ to $j$ if at
least one transaction from one of the addresses defining $i$ to one of the
addresses defining $j$ occurs at the considered time scale.

In the present work, when not otherwise stated, we consider the Bitcoin User
Network at a weekly time scale.

#### Price bubbles data

The valuation of cryptocurrencies is an emerging field of studies, and there’s
no unanimous consensus on what the bitcoin real value should be nor on the
actual definition to detect bitcoin bubbles. In this paper we rely on the
bubbles identified in [8]. The method to identify bubbles developed in [8]
uses a generalized version of Metcalfe’s law to determine bitcoin fundamental
value, which is shown to be heavily exceeded, on at least four occasions,
listed in table 1.

##  3 Results

#### Connected Components

(a) CC size (b) CC numerosity (c) LCC ratio

Figure 1: The evolution of connected components(CC) over time. Panel (a) shows
the percentage of nodes in the weak and strong large connected component.
Panel (b) depicts the number of connected components over time. Panel (c)
shows the ratio between the largest connected component(LCC) and the seocnd
largest connected component

As a first inspection we observe the transaction networks connected components
evolution over time.

Aside from the initial phase until fall 2010, we can detect the emergence of
both the huge Weakly-Connected Component(WCC) and the Strongly-Connected
Component(SCC) in two different points in time. Figure 1 panel (c) shows the
ratio between the largest connected component and the second largest one. As
one can see both at the strong ad weak level the huge connected component
emerge after 2012.

In panel (a) we observe that after the initial phase two clear different
trends emerge at the strong and weak level: the huge weak connected component
size is stable around the $80\%$ of the total number of nodes, while the huge
strong connected component oscillates between the $10\%$ and the $40\%$ of the
nodes; in the specific we observe a long plateau from 2014 to 2016 where the
huge strong connected component size is around $40\%$ of the total nodes and
then gets back to the stable value of $10\%$.

At last, panel (b) shows the number of connected components: it depicts a
situation where the large majority of connected components is composed of a
very small number of nodes.

It seems clear that while a large majority of nodes is connected, the paths
drawn by transactions are usually one way, implying the reduced dimension of
the strongly connected component. Then nodes out of the main cluster tend to
be isolated, thus explaining the large number of connected components.

#### Bitcoin Disassortativity

Figure 2: Evolution of the four directed variants of Newman’s assortativity
coefficient, revealing the (weakly) disassortative character of our BUNs.
Moreover, since $r_{dir}^{out-in}$ is ‘asymptotically’ zero, one can conclude
that $e_{jk}\simeq q_{j}^{out}q_{k}^{in}$ (and analogously for the other
indices).

Let us now consider the assortativity of Bitcoin Transactions Networks. A
network is said to be assortative when nodes with large degree tend to connect
to each other, as opposed to disassortative networks, where nodes with large
degree connects to nodes with low degree. Following [13], on undirected
networks one defines the assortativity coefficient $r_{und}$ as:

$r_{und}=\frac{\sum_{j,k}jk(e_{jk}-q_{j}q_{k})}{\sigma^{2}_{q}}$ (1)

where the sum runs over the ’excess degrees’ of a node: imagine to reach a
vertex following a specific edge, the ’excess degree’ of the vertex is the
degree of the vertex minus that specific edge you followed. $q_{k}$ is the
‘excess degree’ probability distribution, reading

$q_{k}\propto p_{k+1}$ (2)

(with $p_{k+1}$ being the plain degree distribution), $\sigma^{2}_{q}$ is its
standard deviation and $e_{jk}$ is the fraction of edges in the network
connecting nodes of degree $j$ with nodes of degree $k$. Naturally,
$\sum_{j}e_{jk}=q_{k}$.

When considering directed networks [14], instead, four variants of the
aforementioned Pearson coefficient can be calculated, i.e. the ones accounting
for the correlation between out-degrees and out-degrees, out-degrees and in-
degrees, in-degrees and out-degrees, in-degrees and in-degrees. For example,
one of thee variants reads

$r_{dir}^{out-
in}=\frac{\sum_{j,k}jk(e_{jk}-q_{j}^{out}q_{k}^{in})}{\sigma_{q^{out}}\sigma_{q^{in}}}$
(3)

where $e_{jk}$ now represents the percentage of edges starting from nodes
whose out-degree is $j$ and ending on nodes whose in-degree is $k$. Naturally,
it also holds true that $\sum_{j}e_{jk}=q_{k}^{in}$.

Plotting the evolution of the aforementioned coefficients on our BUNs shows
their weakly disassortative nature (see figure 2). In particular, since
$r_{dir}^{out-in}$ is ‘asymptotically’ zero, one can conclude that
$e_{jk}\simeq q_{j}^{out}q_{k}^{in}$ \- and analogously for the other indices
of direct assortativity (the years until 2011 can be considered as a
‘transient’ period where the Bitcoin ecosystem was still of reduced
dimensions, hence sensitive to even small structural changes).

#### Bitcoin core-periphery structure.

Our first result concerns the Bitcoin _core-periphery-ness_. In order to
analyse it, we have run a recently-proposed method [15] based on the
multivariate extension of the _surprise_ score function.

Figure 3: Evolution of the percentage of nodes belonging to the core portion
of the Bitcoin User Network.  Figure 4: Evolution of the temporal z-score
$z_{X}^{(t)}=\frac{X^{(t)}-\overline{X}}{s_{X}}$ for the percentage of core
nodes (the mean $\overline{X}$ and the standard deviation $s_{X}$ have been
computed over a sample of values covering the six months before time $t$).
Peaks are clearly visible in correspondence of the bubbles identified by the
shaded areas, thus indicating the existence of periods during which the price
and the structural quantities of interest co-evolve. Figure 5: Evolution of
the temporal z-score $z_{X}^{(t)}=\frac{X^{(t)}-\overline{X}}{s_{X}}$ for the
number of periphery links.

The application of the multivariate surprise to the partition induced by the
bow-tie structure reveals that it indeed induces a significant core-periphery
structure (the surprise score function is steadily below the threshold of
$5\%$), the core being the SCC and the periphery being composed by all the
other nodes. Figure 3 shows the evolution of the percentage of nodes composing
the core and the periphery of our BUNs. As expected from the results
concerning the SCC, the periphery contains the vast majority of nodes
throughout the entire Bitcoin history.

As fig. 3 shows, there seem to be three different phases: the first one
coincides with the biennium 2012-2014, during which the core portion of the
BUN steadily rises until it reaches the 40% of the network; afterwards, during
the biennium 2014-2016, it remains quite constant; then, during the last two
years covered by our dataset (i.e. 2016-2018), the core portion of the BUN
shrinks and the percentage of nodes belonging to it goes back to the pre-2012
values. In order to gain insight into the correlations between the evolution
of purely topological quantities and the Bitcoin price, let us plot the trend
of the _temporal z-score_ , defined as

$z_{X}^{(t)}=\frac{X^{(t)}-\overline{X}}{s_{X}}$ (4)

for a generic quantity $X$, where the mean
$\overline{X}=\sum_{t}\frac{X^{t}}{T}$ and the standard deviation
$s_{X}=\sqrt{\overline{X^{2}}-\overline{X}^{2}}$ have been computed over a
sample of values covering the six months before time $t$.

As fig. 4 and 5 show, the calculation of the temporal z-score for the
percentage of core nodes/periphery links reveals the presence of peaks in
correspondence of the first three bubbles (identified by the shaded areas),
thus indicating the existence of periods during which the price and the
structural quantities of interest co-evolve. In particular, while the number
of links within the core and the periphery rises significantly, with respect
to the previous temporal interval, periods in-between the bubbles are,
instead, characterized by a decrease of the statistical significance of the
same quantities.

#### Bitcoin bow-tie structure.

The core-periphery structure of the BUN is characterized by an additional
structure, known as _bow-tie_ structure. The definition of _bow-tieness_ rests
upon the concept of _reachability_ : we say that $j$ is reachable from $i$ if
there exist a path from $i$ to $j$. A directed graph is said to be strongly
connected if any two nodes are mutually reachable. Mutual reachability is an
equivalence relation on the vertices of a graph, the equivalence classes being
the strongly connected components of the graph itself. The bow-tie
decomposition of a graph, hence, consists of the following sets of nodes [16]:

  * •

Strongly-connected Component:$\text{SCC}\equiv S$.
Each node within the Strongly-Connected Component(SCC) can be reached by any
other node within it. This means that a directed path exists connecting each
node with each other node;

  * •

In-Component $\equiv\\{i\in V\setminus S\>|\>\text{$S$ is reachable from
$i$}\\}$;

  * •

Out-Component $\equiv\\{i\in V\setminus S\>|\>\text{$i$ is reachable from
$S$}\\}$;

  * •

Tubes $\equiv\\{i\in V\setminus
S\>\cup\>\text{IN}\cup\>\text{OUT}\>|\>\text{$i$ is reachable from IN and OUT
is reachable from $i$}\\}$;

  * •

In-Tendrils $\equiv\\{i\in V\setminus S\>|\>\text{$i$ is reachable from IN and
OUT is not reachable from $i$}\\}$;

  * •

Out-Tendrils $\equiv\\{i\in V\setminus S\>|\>\text{$i$ is not reachable from
IN and OUT is reachable from $i$}\\}$;

  * •

Others: $\equiv\\{i\in V\setminus
S\>\cup\>\text{IN}\>\cup\>\text{OUT}\>\cup\>\text{TUBES}\>\cup\>\text{IN-
TENDRILS}\>\cup\>\text{OUT-TENDRILS}\\}$.

Figure 6: Temporal evolution of the percentage of nodes constituting the bow-
tie components listed above. While the SCC coincides with the core portion of
the BUN, the periphery gathers all other node subsets. The green vertical line
indicates the Mt.Gox breakdown.

Generally speaking, a large SCC, incorporating the vast majority of nodes,
starts emerging in 2012, ‘stabilizes’ around mid-2013 and persists until 2016.
More specifically, during the biennium 2012-2013 the SCC steadily rises until
it reaches $\simeq 30\%$ of the network size; afterwards, during the biennium
2014-2015, it remains quite constant; then, during the last two years covered
by our data set (i.e. 2016-2018), it shrinks and the percentage of nodes
belonging to it goes back to the pre-2012 values. While in the biennium
2014-2016 the percentage of nodes constituting the SCC is larger than the
percentage of nodes belonging to the other subsets, since 2016 this is no
longer true: in fact, while both the SCC and the OUT-component shrink, the IN-
component becomes the dominant portion of the network.

Different results have been reported in [17]: however, this may be due to the
different data collection and data mining processes implemented there.

#### Dyadic motifs.

Let us now consider the structure of dyadic motifs, meaning the patterns
involving two nodes. Althought the number of dyadic motifs is very limited and
of easy definition, these may be of interest to study the evolution of
economic networks [18, 19].

On a directed network $G(N,E)$, given two nodes $i,j\in E$ you can observe
three disjoint possibilities:

  * •

a reciprocated dyad: $(i,j),(j,i)\in E$, meaning a connection exists both from
$i$ to $j$ and vice versa. We denote the total number of reciprocated dyads in
the network as $L^{\leftrightarrow}$

  * •

a non reciprocated dyad: $(i,j)\in E$ or $(j,i)\in E$, meaning only one of the
two possible links exists. We denote the total number of non reciprocated
dyads as $L^{\rightarrow}$

  * •

an empty dyad$(i,j),(j,i)\notin E$, meaning no link exists between the two
nodes. The total number of empty dyads is denoted by
$L^{\centernot{\leftrightarrow}}$

Dyadic motif $m$ | Abundance $N_{m}$
---|---
$L^{\nleftrightarrow}$ | $\sum_{i\neq j}(1-a_{ij})(1-a_{ji})$
$L^{\rightarrow}$ | $\sum_{i\neq j}a_{ij}(1-a_{ji})$
$L^{\leftrightarrow}$ | $\sum_{i\neq j}a_{ij}a_{ji}$
Table 2: Definition of the abundance of dyadic motifs.

Let us now study dyadic motifs by adopting a different approach with respect
to the one employed so far. Instead of using the temporal $z$-score to spot
‘temporal’ outliers, i.e. values that are statistically significant with
respect to a time average, let us consider an index that points out quantities
not compatible with a given null model. In order to do so, we will employ the
Directed Binary Configuration Model as null model, from the Exponential Random
Graph Model(ERGM) family[18, 20]. In this framework, a $z$-score of the kind

$z[X]=\frac{X(\mathbf{A}^{*})-\langle X\rangle}{\sigma[X]}$ (5)

remains naturally defined, where $X(\mathbf{A}^{*})$ is the empirical value of
the quantity of interest (i.e. observed on the original network
$\mathbf{A}^{*}$), $\langle X\rangle$ and $\sigma[X]$ are, respectively, its
expectation value and its standard deviation, both computed on the ensemble
induced by the DBCM. The interpretation of this $z$-score is the following
one: values such that $z[X]>+3$ signal that the empirical value is
significantly larger than expected while values such that $z[X]<-3$) signal
that the empirical value is significantly smaller than expected. In both cases
one may conclude that the empirical value $X(\mathbf{A}^{*})$ is not
compatible with the specific model and something else is required to fully
account for it. On the other hand, if $-3\leq z[X]\leq+3$, there is no
evidence of a significant deviation from the expected value and one may
conclude that $X(\mathbf{A}^{*})$ is completely explained by the constraints
defining the model at hand.

This exercise gives us an idea of how much the weekly BTC transaction network
differ from the randomized network derived from its degree distribution. The
results are shown in 7. The first thing that stands out are the out-of-scale
values of the reciprocated dyads z-score, which makes totale sense:the BUNs
(and in general the networks obtained from BTC transactions) are very sparse
network, were the large majority of nodes have very few connections; it makes
sense that each bilateral connection would be very odd in comparison to a
randomized network with the same number of nodes and links.

The interpretation of the behavior of empty and single dyads is analogous: by
chance, a larger-than-observed number of non-reciprocated dyads are created,
whence their over-representation within the DBCM ensemble and the negative
$z$-score recovered by our analysis. In order to understand why this implies
that the DBCM tends to create less-than-observed empty dyads, let us imagine
to ‘destroy’ a reciprocal dyad, by decoupling the two paired links: in order
to create ‘more’ single dyads, one of the two links mus be redirected towards
a previously disconnected node; upon doing so, a reciprocal dyad disappears,
as well as an empty dyad, while two single dyads are created.

Figure 7: Left panels: evolution of the empirical value of the dyads (from top
to bottom; empty, non-reciprocated and reciprocated dyads). Right panels:
evolution of the dyads $z$-score, computed over the ensemble induced by the
Directed Binary Configuration Model. Dashed, gray lines signal the values of
$\pm 2$ and $\pm 3$. Points are coloured according to the value of the log-
return of the Bitcoin price in USD, in that week. Shaded areas indicate
periods during which the price grows.

#### Centrality and centralization

We computed four different centrality measures: degree centrality, closeness
centrality, betweenness centrality and eigenvector centrality.

The normalized closeness centrality is the average length of the shortest path
between a node and all the others. The more central a node is, the closer it
is to all the other nodes. Betweenness centrality denotes the number of times
a node is on the shortest path connecting two other nodes. It was introduced
to measure the control a node had on the communication layer between other
nodes on the network Eigenvector centrality measures the influence of a node
over the network. Each node is assigned a relative score, based on the idea
that being connected to high-scored nodes increases a node score.

We condensate the information given by centralities for each network datapoint
by considering two different aggregate indices derived from each centrality
distribution: the Gini index and the centralisation index.

The Gini coefficient attempts to measure the unevenness of a distribution of a
certain quantity111Usually, it is employed to measure the unevenness of the
income distribution[21]: given a set of values $\\{c_{i}\\}_{i=1}^{N}$, the
Gini index is defined as

$G_{c}=\frac{\sum_{i=1}^{N}\sum_{j=1}^{N}|c_{i}-c_{j}|}{2N\sum_{i=1}^{N}{c_{i}}}$
(6)

and assumes values between $0$ and $1$; while a Gini index of $0$ indicates
perfect evenness (e.g. everyone has exactly the same income), a Gini index of
$1$ indicates perfect unevenness (e.g. a population whose entire income is
concentrated in the hands of a single individual).

Applying the Gini coefficient to the degrees of our BUNs aims at shedding
light on the (un)evenness of the nodes degree centralization: while a value
close to $0$ would depict an ecosystem where all actors have exactly the same
number of interactions with each other, a value close to $1$ would indicate
that there are nodes participating to the vast majority of transactions.

Centralization indices are global measures intended to measure the centrality
of the entire network (instead of providing a rank of its nodes). In
mathematical terms, the centralization reads

$C_{c}=\frac{\sum_{i=1}^{N}(c^{*}-c_{i})}{\max\left\\{\sum_{i=1}^{N}(c^{*}-c_{i})\right\\}}$
(7)

where $c^{*}=\max\\{c_{i}\\}_{i=1}^{N}$ represents the empirical, maximum
value of the chosen centrality measure (i.e. computed on the network under
consideration) and the denominator is calculated over a benchmark graph,
defined as the one providing the maximum attainable value of the quantity
$\sum_{i=1}^{N}(c^{*}-c_{i})$.

The benchmark graph for which is nothing else than a star graph (with the same
number of nodes of the network under inspection). For each centrality we
should compute the centralization index by computing the star-graph
centralization maximum value.

$C_{k}=\frac{\sum_{i=1}^{N}(k^{*}-k_{i})}{(N-1)(N-2)}$ (8)

and the degree centralization would reveal to us if (and, in case, ‘how much’)
Bitcoin has become similar to a star graph at a certain point during its
history.

(a) Gini Indices (b) Centralisation Indices

Figure 8: Temporal evolution of Gini coefficient and centralisation indices
for the four centrality measures under examination

From the results in figure 8 appears that, after a period of growth lasted
until mid-2013, during which it reached values as large as $0.75$, the Gini
coefficient has decreased and is now steadily around the value of $0.5$.
Overall, we would like to stress that $0.5$ is not a small value: in fact, it
describes an ecosystem where the $50\%$ of connections are incident to the
$1\%$ of nodes. It is also interesting to notice the big leap down of the Gini
coefficient in 2013: during that year, Mt. Gox (which managed $\simeq 70\%$ of
transactions at the time [22]) started the down-ward spiral which eventually
led to its bankruptcy in 2014: USD withdrawals halting, financial
investigations and expensive lawsuits weakened the trading website ability to
stay on the market. The final blow was the public discovery of a huge theft of
around $750.000$ bitcoins, which went on undetected for years.

The huge decrease of the Gini coefficient may be, thus, related to the Mt. Gox
‘loss of prominence’ in the Bitcoin ecosystem. On the other hand, bubble
periods seem to have little correlation with the evolution of the Gini
coefficient.

The evolution of the Gini coefficient may lead us to imagine that the Bitcoin
ecosystem has become similar to a very centralized structure, pretty much
similar to a star graph, at some point during its history. In order to answer
this question, we have computed the so-called centralization index at the
weekly time scale (from [7]).

While during the initial phases of its life, Bitcoin was indeed quite similar
to a star graph, figure 8 reveals that the degree-centralization has quickly
stabilized around very small values. Overall, we may, thus, conclude that
Bitcoin is not evolving towards a star-like structure, where a single central
node participates to all transactions. However, the large value of the Gini
coefficient let us suspect that there may be several hubs: hence, the
unrealistic picture of a star-like structure may be replaced by the more
realistic one depicting several ‘locally star-like’ structures. The centers of
these structures are ‘local hubs’, i.e. vertices with a large number of
connections, that are crossed by a large percentage of paths and that are
connected among them. In order to justify the apparent contradiction, in the
appendix we show a toy model which partially reproduces the results.

#### Small-world

(a) APL (b) Clustering

Figure 9: Panel (a): scatter plot of the logarithm of the number of nodes
versus the average path length $L$. Each marker is a Weekly BUN between 2011
and 2014. On the right, the time series of the BUNs clustering coefficient
versus the clustering coefficient of the Erdos Renyi random graph generated
with the same number of nodes and edges of the same week BUN counterpart.

In order to explore the small-word properties of Bitcoin Transaction Networks,
we compute the Average Path Length (APL) $L$ and the global clustering
coefficient $c$ evolution of the BUNs in our dataset[23]. Because of
computational constraints, the analysis is limited to the first half of the
period under examination(until 2014). Please note that the measures are taken
on the giant connected component. The results are displayed in figure 9: in
the period under exam, both the measures seems to be stable. The clustering
coefficient floats around very small values, in the order of $10^{-2}$, while
the APL is stable around $10$ and $15$, with rare exceptions.

From figure 9 (a) it is clear that a linear relation between $L$ and $\log{N}$
does hold. This means that at least the close-connected behaviour of small-
world network is present in Bitcoin Transaction Networks. Keeping this in
mind, we can’t ignore that the BUNs we analysed have a very low clustering
nature, meaning that it’s hard to observe two neighbours nodes(meaning users)
of a node connected among each others.

In order to put the clustering values in perspective we plot the BUNs
coefficient versus the Erdos-Renyi randomized clustering coefficient. In
figure 9 (b) we can see that while the magnitude of the clustering is small
both on the original and on the Erdos-Renyi randomization, we notice that
constantly over time the original clustering is at least one order of
magnitude larger than the Erdos-Renyi: this means that while there are not
many triangles, Bitcoin Transactions Networks have more triangles than an
equivalent random graph.

##  4 Discussion

In this paper we analyse the impact of structural properties of the Bitcoin
transaction network on the generation and crash of bubbles in the exchange
with respect to fiat currencies. Specifically, we examine network features
such as heterogeneity of the degree distributions and frequency of
connectivity patterns (i.e., motifs). We find significant changes in these
properties during the period of price bubbles. A more detailed analysis
unveils that, during the first bubble, the frequency of motifs indicating the
relationship hubs have with new, low-degree users changes significantly; this
suggests that hubs take an even more important role in becoming liquidity
providers. These results are confirmed in the second bubble: There, by
analysing the heterogeneity of the in-degree, out-degree, and total degree
distributions, we find that there is a significant widening (narrowing) of the
out-degree (in-degree) distributions, whereas the total degree does not change
its distribution significantly. By performing additional analyses on the two
largest hubs, we find that these structural changes - similar to what is
observed during the first bubble - is likely to be caused by the
_centralising_ role hubs take on as liquidity providers.

Although we find that measures can explain well some price bubbles but not
others, these results highlight that tracking properties of hubs in the
transaction network is key for understanding the underlying mechanisms of a
bubble. Moreover, at least in the first three Bitcoin bubbles, the behaviour
of hubs significantly increased the systemic risk of the Bitcoin economy,
eventually leading to systemic failure and sudden price crashes.

These results also suggest that Bitcoin bubbles are difficult to forecast, but
can be prevented, or at least alleviated, by introducing policies that aim at
reducing the importance of large hubs in the network. In future work, we plan
to extend our analysis by introducing new structural measures and by covering
all the bubbles that happened to date.

These structural changes suggest the presence of two hubs that centralise the
market by selling bitcoins to most of the traders that enter the market during
the bubble, resulting in a significant increase of the systemic risk. Indeed,
if only a few hubs account for most of the transactions in the network, if at
any point in time one of them fails, the whole network may crash. This is
exactly what happened on April 10th 2013, when Mt Gox, the major Bitcoin
exchange, broke under the high trading volume, triggering the burst of the
bubble.

## References

  * [1] Nakamoto, S. “Bitcoin: A peer-to-peer electronic cash system.” (2008).
  * [2] Halaburda, H., Sarvary, M. _Beyond bitcoin: The Economics of Digital Currencies_. Palgrave Macmillan (2016).
  * [3] Glaser, F. “Pervasive decentralisation of digital infrastructures: a framework for blockchain enabled system and use case analysis.” In _Proceedings of the Hawaii International Conference on System Sciences 2017 (HICSS-50)_ (2017) .
  * [4] Antonopoulos, A. M. _Mastering Bitcoin: Programming the Open Blockchain_. O’Reilly Media, Inc. (2017).
  * [5] Javarone, M. A., Wright, C. S. “From Bitcoin to Bitcoin Cash: a network analysis.” _arXiv:1804.02350v2_ .
  * [6] Bovet, A., _et al._ “The evolving liaisons between the transaction networks of Bitcoin and its price dynamics.” _arXiv:1907.03577_ .
  * [7] Lin, J.-H., Primicerio, K., Squartini, T., Decker, C., Tessone, C. J. “Lightning Network: a second path towards centralisation of the Bitcoin economy.” _arXiv:2002.02819_ .
  * [8] Wheatley, S., Sornette, D., Huber, T., Reppen, M., Gantner, R. N. “Are Bitcoin bubbles predictable? Combining a generalized Metcalfe’s Law and the Log-Periodic Power Law Singularity model.” _Royal Society Open Science_ 6.6 180538 (2019) 10.1098/rsos.180538 https://royalsocietypublishing.org/doi/pdf/10.1098/rsos.180538 URL https://royalsocietypublishing.org/doi/abs/10.1098/rsos.180538.
  * [9] Androulaki, E., Karame, G. O., Roeschlin, M., Scherer, T., Capkun, S. “Evaluating user privacy in bitcoin.” In _International Conference on Financial Cryptography and Data Security_ Springer 34–51 (2013) .
  * [10] Tasca, P., Hayes, A., Liu, S. “The evolution of the bitcoin economy: extracting and analyzing the network of payment relationships.” _J. Risk Financ._ 19.2 94–126 (2018).
  * [11] Harrigan, M., Fretter, C. “The unreasonable effectiveness of address clustering.” In _Ubiquitous Intelligence & Computing, Advanced and Trusted Computing, Scalable Computing and Communications, Cloud and Big Data Computing, Internet of People, and Smart World Congress (UIC/ATC/ScalCom/CBDCom/IoP/SmartWorld), 2016 Intl IEEE Conferences_ IEEE 368–373 (2016) .
  * [12] Ron, D., Shamir, A. “Quantitative analysis of the full bitcoin transaction graph.” In _International Conference on Financial Cryptography and Data Security_ Springer 6–24 (2013) .
  * [13] Newman, M. E. “Mixing patterns in networks.” _Physical Review E_ 67.2 026126 (2003).
  * [14] Noldus, R., Van Mieghem, P. “Assortativity in complex networks.” _Journal of Complex Networks_ 3.4 507–542 (2015).
  * [15] de Jeude, J. v. L., Caldarelli, G., Squartini, T. “Detecting core-periphery structures by surprise.” _Europhysics Letters_ 125.6 68001 (2019).
  * [16] de Jeude, J. v. L., Di Clemente, R., Caldarelli, G., Saracco, F., Squartini, T. “Reconstructing mesoscale network structures.” _Complexity_ 2019.5120581.
  * [17] Maesa, D. D. F., Marino, A., Ricci, L. “The bow tie structure of the bitcoin users graph.” _Applied Network Science_ 4.1 1–22 (2019).
  * [18] Squartini, T., van Lelyveld, I., Garlaschelli, D. “Early-warning signals of topological collapse in interbank networks.” _Sci. Rep._ 3.1 3357 (2013).
  * [19] Squartini, T., Garlaschelli, D. “Stationarity, non-stationarity and early warning signals in economic networks.” _Journal of Complex Networks_ 3.1 1–21 (2015).
  * [20] Garlaschelli, D., Loffredo, M. I. “Maximum likelihood: Extracting unbiased information from complex networks.” _Physical Review E_ 78.1 015101 (2008).
  * [21] Dixon, P. M., Weiner, J., Mitchell-Olds, T., Woodley, R. “Bootstrapping the Gini coefficient of inequality.” _Ecology_ 68.5 1548–1551 (1987).
  * [22] Decker, C., Wattenhofer, R. “Bitcoin Transaction Malleability and MtGox.” In M. Kutyłowski, J. Vaidya (Eds.), _Computer Security - ESORICS 2014_ Cham: Springer International Publishing 313–326 (2014) .
  * [23] Albert, R., Barabási, A.-L. “Statistical mechanics of complex networks.” _Rev. Mod. Phys._ 74.1 47 (2002).
  * [24] Aldecoa, R., Marin, I. “Surprise maximization reveals the community structure of complex networks.” _Scientific Reports_ 3.1060.
  * [25] Aldecoa, R., Marin, I. “Exploring the limits of community detection strategies in complex networks.” _Scientific Reports_ 3.2216.
  * [26] Nicolini, C., Bifone, A. “Exploring the limits of community detection strategies in complex networks.” _Scientific Reports_ 6.19250.
  * [27] Freeman, L. C. “Centrality in social networks conceptual clarification.” _Social networks_ 1.3 215–239 (1978).

##  Appendix A Surprise

Originally proposed to detect communities [24, 25, 26], the surprise reads

$\mathscr{S}=\sum_{l\geq
l_{\bullet}^{*}}^{\min\\{L,V_{\bullet}\\}}\frac{\binom{V_{\bullet}}{l}\binom{V_{\circ}}{L-l}}{\binom{V}{L}}$
(9)

where $V$ is the volume of the network (coinciding with the total number of
node pairs, i.e. $V=N(N-1)$ in the directed case), $V_{\bullet}$ is the total
number of intracluster pairs (i.e. the number of node pairs _within_ the
individuated communities), $L$ is the total number of links and $l^{*}$ is the
observed number of intracluster links (i.e. _within_ the individuated
communities). In other words, $\mathscr{S}$ is the p-value of an
hypergeometric distribution, aimed at testing the statistical significance of
a given partition. Such a distribution, in fact, describes the probability of
observing $l$ successes in $L$ draws (without replacement) from a finite
population of size $V$ that contains exactly $V_{\bullet}$ objects with the
desired feature (in our case, being an intracluster pair), each draw being
interpreted as either a ‘success’ or a ‘failure’: the lower this probability,
the ‘better’ the individuated partition (i.e. the more surprising its
observation).

When carrying out a community detection exercise, links are understood as
belonging to two different categories, i.e. the ones _within_ the clusters and
the ones _between_ the clusters. The surprise-based formalism, however, can be
extended to detect _bimodular_ structures, a terminology that is intended to
describe either _bipartite_ or _core-periphery_ structures. In this case,
three different ‘species’ of links are needed (e.g. the core, the periphery
and the core-periphery ones): for this reason, we need to consider the
multinomial version of the surprise [15], reading

$\mathscr{S}_{\parallel}=\sum_{i\geq l_{\bullet}^{*}}\sum_{j\geq
l_{\circ}^{*}}\frac{\binom{V_{\bullet}}{i}\binom{V_{\circ}}{j}\binom{V-(V_{\bullet}+V_{\circ})}{L-(i+j)}}{\binom{V}{L}}$
(10)

(the presence of three different binomial coefficients allows three different
kinds of links to be accounted for). From a technical point of view,
$\mathscr{S}_{\parallel}$ is the p-value of a multivariate hypergeometric
distribution, describing the probability of $i+j$ successes in $L$ draws
(without replacement), from a finite population of size $V$ that contains
exactly $V_{\bullet}$ objects with a first specific feature and $V_{\circ}$
objects with a second specific feature, wherein each draw is either a
‘success’ or a ‘failure’: analogously to the univariate case,
$i+j\in[l_{\bullet}^{*}+l_{\circ}^{*},\min\\{L,V_{\bullet}+V_{\circ}]$.

According to the interpretation proposed in [15], revealing the core-periphery
structure minimizing the surprise means individuating the partition that is
least likely to be explained by invoking the Random Graph Model (RGM) with
respect to the Stochastic Block Model (SBM). This, in turn, means that the BUN
is indeed characterized by subgraphs with very different link densities, an
evidence that cannot be explained by a model defined by just one global
parameter (i.e. the one characterizing the RGM).

##  Appendix B Reciprocity

The ratio of reciprocated dyads over the total number of links is called
_reciprocity_ , and reads:

$r=\frac{\sum_{i}\sum_{j(\neq i)}a_{ij}a_{ji}}{\sum_{i}\sum_{j(\neq
i)}a_{ij}}\equiv\frac{L^{\leftrightarrow}}{L}$ (11)

The evolution of the BUN reciprocity structure is plotted in fig. 10: the
first observation concerns the absolute value of $r$ which amounts at few
percents. This means that, overall, the reciprocity of the BUN is quite low
(although the picture changes if we consider the value of the reciprocity
within the SCC). Interestingly, however, the evolution of the reciprocity
shows some peaks of activity in correspondence of the bubbles; the relevance
of this kind of information is confirmed upon calculating the reciprocity
temporal z-score defined in eq. 4: it rises significantly in correspondence of
the bubbles and decreases during the inter-bubble periods.

Figure 10: Evolution of the empirical reciprocity $r$ and of its temporal
z-score. Although the absolute value of the reciprocity is, overall, quite
low, its evolution is characterized by peaks whose statistical relevance is
confirmed upon calculating the temporal z-score. Analogously to other
topological quantities, the BUN reciprocity significantly rises in
correspondence of the bubbles and decreases during the inter-bubble periods.

##  Appendix C Centralization formulas

The benchmark graph for which is nothing else than a star graph (with the same
number of nodes of the network under inspection): hence, here we report the
mathematical formulas from [7] for the centralization indices:

$C_{k}=\frac{\sum_{i=1}^{N}(k^{*}-k_{i})}{(N-1)(N-2)}$ (12)

where $\\{k_{i}\\}_{i}$ is the degree centralization of each node.

$C_{c^{C}}=\frac{\sum_{i=1}^{N}(c^{*}-c^{c}_{i})}{(N-1)(N-2)/(2n-3)}$ (13)

Where $\\{c_{i}^{c}\\}_{i}$ is the normalized closeness centrality of each
node.

$C_{b^{C}}=\frac{\sum_{i=1}^{N}(b^{*}-b^{c}_{i})}{(N-1)^{2}(N-2)/2}$ (14)

Where $\\{b_{i}^{c}\\}_{i}$ is the normalized betweenness centrality of each
node. [27].

$C_{e^{C}}=\frac{\sum_{i=1}^{N}(e^{*}-e^{c}_{i})}{(\sqrt{N-1}-1)(N-1)/(\sqrt{N-1}+N-1)}$
(15)

Where $\\{b_{i}^{c}\\}_{i}$ is the normalized eigenvector centrality of each
node.

Figure 11: Upper row, from left to right: the Gini coefficient of degree
centrality, closeness centrality, betweenness centrality and eigenvector
centrality versus the number of nodes for each data point. Bottom row, form
left to right: centralisation index of degree centrality, closeness
centrality, betweenness centrality and eigenvector centrality.

#### Toy model

A toy model can help understanding the two apparently contradictory results
provided by the Gini coefficient and the degree-centralization. Imagine
$N_{h}$ hubs connected between them and $N_{l}$ leaves connected to each of
them; hence, total number of nodes is $N=N_{h}(N_{l}+1)$, the degree of each
hub reads $k_{h}=(N_{h}-1)+N_{l}$, the degree of each leave reads $k_{l}=1$
and

$\displaystyle G_{k}$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{N}\sum_{j=1}^{N}|k_{i}-k_{j}|}{2N\sum_{i=1}^{N}{k_{i}}}=\frac{2[(N_{h}-1)+N_{l}-1]N_{h}^{2}N_{l}}{2N\left[N_{h}(N_{h}-1)+2N_{h}N_{l}\right]}$
$\displaystyle=$
$\displaystyle\frac{[(N_{h}-1)+N_{l}-1]N_{h}^{2}N_{l}}{N_{h}(N_{l}+1)\left[N_{h}(N_{h}-1)+2N_{h}N_{l}\right]}$
$\displaystyle=$
$\displaystyle\frac{(N_{h}+N_{l}-2)N_{l}}{(N_{l}+1)\left[(N_{h}-1)+2N_{l}\right]}\simeq\frac{N_{h}+N_{l}}{N_{h}+2N_{l}};$

now, $G_{k}\simeq 2/3$ if $N_{l}=N_{h}$, i.e. if each hub is linked to a
number of ‘leaves’ that matches the total number of hubs and $G_{k}\rightarrow
1/2$ as $N_{l}\rightarrow+\infty$, i.e. if the number of leaves per hub
becomes ‘very large’. In this setting, we have that

$C_{k}=\frac{\sum_{i=1}^{N}(k^{*}-k_{i})}{(N-1)(N-2)}=\frac{((N_{h}-1)+N_{l}-1)*N_{h}N_{l}}{(N_{h}(N_{l}+1)-1)(N_{h}(N_{l}+1)-2)}\simeq\frac{N_{h}+N_{l}}{N_{h}N_{l}}$
(17)

which amounts at $C_{k}\simeq 0.02$ if we set $N_{h}=N_{l}=100$. Hence,
opportunely tuning the parameters of our model we can recover core-periphery
structures for which a large Gini coefficient co-exists with a small degree-
centralization.


# Reconciling quantum and classical spectral theories of ultrastrong coupling:
Role of cavity bath coupling and gauge corrections

Stephen Hughes<EMAIL_ADDRESS>Department of Physics, Engineering Physics
and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada  Chris Gustin
<EMAIL_ADDRESS>Edward L. Ginzton Laboratory, Stanford University,
Stanford, California 94305, USA<EMAIL_ADDRESS>Franco Nori
<EMAIL_ADDRESS>Theoretical Quantum Physics Laboratory, Cluster for Pioneering
Research, RIKEN, Wako-shi, Saitama 351-0198, Japan  Center for Quantum
Computing, RIKEN, RIKEN, Wako-shi, Saitama 351-0198, Japan  Physics
Department, The University of Michigan, Ann Arbor, Michigan 48109-1040, USA

(September 3, 2024)

###### Abstract

Focusing on the widely adopted Hopfield model with cavity dissipation, we show
how the linear spectrum of an ultrastrongly coupled cavity and a dipole can be
described either classically or quantum mechanically, but only when the
quantum model includes (i) corrections to maintain gauge invariance, and (ii)
a specific type of cavity bath coupling. We also show the impact of this bath
model on the quantum Rabi model, which has no classical analogue in
ultrastrong coupling.

Introduction.—Strong coupling between a single cavity mode and a dipole or
two-level system (TLS) [1, 2, 3] can be well explained quantum mechanically or
classically (or semiclassically) [4, 5, 6]. The characteristic signature of
strong coupling is a splitting in the emitted spectrum by $2g$, where $g$ is
the dipole-cavity coupling rate, which exceeds any losses in the system, e.g.,
$g^{2}>\kappa^{2}/16$ [3] (or more strictly $g^{2}>\kappa^{2}/8$ [7]), with
$\kappa$ the cavity decay rate. Quantum mechanically, this is referred to as
vacuum Rabi splitting, or in classical physics as normal-mode splitting.

Strong coupling has been observed in atoms [2], molecules [8, 9], quantum dots
[10, 11], and circuit QED [12, 13], and is often considered a prerequisite for
exploring unique quantum effects when one moves beyond a weak excitation
approximation or linear response [14, 15]. In a quantum description of cavity-
TLS coupling, multi-photon effects manifest in an anharmonic response [16, 17,
18, 19], which is not captured by the physics of two coupled classical
harmonic oscillators (HOs). Nevertheless, a classical description of the
emitted spectrum under weak excitation is an adequate description of the
system, and one recovers a perfect quantum to classical correspondence of the
light-matter system, albeit with a different interpretation. Indeed, the
quantum interpretation of spontaneous emission can be described in terms of
radiation reaction, vacuum fluctuations, or a mixture of both these effects
[20, 21].

Quantum and classical descriptions of certain light-matter coupling have lead
to interesting interpretations and insights, including the difference between
quantum and classical oscillations in phase qubits [22], classical pseudo-Rabi
oscillations in flux qubits [22], and vacuum Rabi splitting as a manifestation
of linear-dispersion theory [4].

Beyond strong coupling, recent interest in cavity-QED has turned to
ultrastrong coupling (USC) [23, 24, 25, 26, 27, 28, 29, 30], where one cannot
invoke a rotating wave approximation (RWA), typically when $g\geq
0.1\omega_{0}$, where $\omega_{0}$ is the dipole resonance frequency. The
regime of USC presents some fascinating uniquely quantum concepts such as
virtual photons in the ground state [31, 29, 30, 28]. Squeezed vacuum states,
also with no classical analogue, occur in both bosonic and TLS emitter systems
in the USC, which are described using the Hopfield model and the quantum Rabi
model (QRM), respectively [29]. Exciton and many-emitter Dicke systems also
take on the form of the Hopfield model, such as cavity coupling to 2D electron
systems including Landau levels in THz cavities [32]. In the USC regime, these
systems exhibit spectral signatures that reflect the nature of the quantum
Hopfield model. On the other hand, the classical theory of coupled oscillators
(including collectively coupled TLSs in the dilute thermodynamic limit) does
not require a RWA, and it might be expected that the USC regime should also
have a quantum to classical correspondence.

In probing the Hopfield regime, the optical spectrum is typically measured,
and most features are argued as being quantum mechanical in nature, e.g.,
stemming from diamagnetic coupling terms. Yet, the spectral locations of the
upper and lower polaritons can be matched with coupled classical HO theory
[23, 33, 34]. In the presence of dissipation, however, a quantum-classical
correspondence is not known. Although there exists some quantum Langevin
approaches for simplified geometries [23], these are not suitable for
calculating higher-order quantum correlation functions and arbitrary
geometries. Furthermore, the physics of a cavity-coupled TLS is often said to
reduce to cavity-HO physics if one neglects saturation effects [35, 36].
However, this is not true in the USC regime, even with linear response, since
multiple photon states emerge already in excited states, and saturation
effects are unavoidable in the USC regime due to the virtual excitation of the
TLS, including the ground state.

The USC regime presents additional challenges for quantum field models,
including: (i) gauge corrections because of a truncated Hilbert space [27, 37,
38], and (ii) the specific form of the system-bath interaction for the cavity
mode matters [39, 40]. Since a classically-coupled HO model (expressed
entirely in terms of classical electromagnetic fields) has no issues with
gauge invariance, it is essential to seek out if and when such a quantum-
classical correspondence can be made. This is not just motivated by
fundamental physics reasons, but is practically important since many of the
emerging USC experiments require some sort of modelling with classical Maxwell
solvers [41, 33].

In this work, we first show that, for a lossless system, the spectral poles
(resonances) of the Hopfield model precisely overlap a classical HO solution,
and these deviate from the QRM as soon as one enters the USC regime. We then
introduce a spectral theory of the dissipative Hopfield model, using a gauge-
invariant master equation theory expressed in the multi-polar gauge, and show
how it is possible to identify a specific form of the system-bath coupling
that matches the classical solution.

We choose a common and established classical theory, based on a normal-mode
expansion of the cavity Green function with phenomenological decay; a more
rigorous approach could use quasinormal modes [42, 43]. Finally, we study the
impact of this model on both the dissipative Hopfield model and the
dissipative QRM, and show the striking differences between these two models
for different $\eta\equiv g/\omega_{c}$. Usually $\eta>0.1$ is the criterion
for the USC regime. A simple schematic of lossy cavity-QED systems are shown
in Fig. 1.

Figure 1: Schematic of two systems in dissipative cavity-QED that can realize
USC, including (a) an atom inside a cavity, which has a decay rate $\kappa$,
and (b) a planar cavity coupled to a collective emitter system. The emitters
can be treated as a bosonic (Hopfield model) or as a TLS (QRM).

Theory.—We first consider the interaction between a bosonic cavity mode, with
creation (annihilation) operator $a^{\dagger},a$, and a bosonic dipole, with
creation (annihilation) operator $b^{\dagger},b$. Neglecting bath losses for
now, in the multi-polar gauge, and with the dipole approximation, the Hopfield
model can be written as ($\hbar=1$)

$H_{\rm
Hop}=\omega_{c}a^{\dagger}a+\omega_{0}b^{\dagger}b+ig(a^{\dagger}-a)(b+b^{\dagger})+D(b+b^{\dagger})^{2},$
(1)

where $\omega_{c}$ ($\omega_{0}$) is the cavity (atom) transition frequency,
and $D=\eta g$ is the diamagnetic amplitude [44]. Note that a naive truncation
of the multipolar gauge Hamiltonian would render this Hopfield $D$ term
infinite, and one must account for the gauge invariance of the truncated
single-mode model to obtain physical and correct results [38, 45].

Using a Bogoliubov transformation, we rewrite this Hamiltonian as $H_{\rm
Hop}=\omega_{c}a^{\dagger}a+\tilde{\omega}_{0}\tilde{b}^{\dagger}\tilde{b}+i\tilde{g}(a^{\dagger}-a)(\tilde{b}+\tilde{b}^{\dagger})+D$,
where $\tilde{\omega}_{0}=\omega_{0}(1+4\eta^{2})^{0.5}$ and
$\tilde{g}^{2}=g^{2}/(1+4\eta^{2})^{0.5}$. Diagonalization yields two
polariton poles [46]:
$\omega^{2}_{\pm}=\frac{1}{2}\left[\tilde{\omega}_{0}^{2}+\omega_{c}^{2}\pm\sqrt{(\tilde{\omega}_{0}^{2}-\omega_{c}^{2})^{2}+16\tilde{g}^{2}\tilde{\omega}_{0}\omega_{c}}\right]$.
Assuming on resonance conditions ($\omega_{0}=\omega_{c}$), then

$\omega_{\pm}=\omega_{0}\sqrt{1+2\eta^{2}\pm 2\eta(1+\eta^{2})^{1/2}}.$ (2)

To lowest order in the counter rotating-wave effects, i.e., to order
$\eta^{2}$ (Bloch-Siegert regime), we obtain $\omega_{\pm}|_{\rm
BS}=\omega_{0}(1+\eta^{2}/2)\pm g$. If one neglects the diamagnetic term, then
$\omega_{\pm}^{0}=\omega_{0}(1\pm 2\eta)^{0.5}$, which is problematic for
$\eta\geq 0.5$.

We can also compare this solution to the QRM, with

$H_{\rm
QRM}=\omega_{c}a^{\dagger}a+\omega_{0}\sigma^{+}\sigma^{-}+ig(a^{\dagger}-a)(\sigma^{+}+\sigma^{-}),$
(3)

where $\sigma^{+}$ ($\sigma^{-}$) is the creation (annihilation) operator for
the TLS, which has important saturation effects. In this case, we have no
diamagnetic term to consider as the relevant TLS term has no effect on energy
differences with the TLS operators. For the QRM, we have an infinite set of
anharmonic eigenenergies, which modifies the Jaynes-Cummings ladder states
because of counter-rotating wave terms. Considering again a Bloch-Siegert
regime (order $\eta^{2}$ coupling) [30, 40], we obtain the poles of the lowest
order polaritons $\omega_{\pm}|_{\rm BS}^{\rm QRM}=\omega_{0}\pm
g(1+\eta^{2}/4)^{0.5}\approx\omega_{0}\pm g$. These QRM pole resonances differ
from the Hopfield model in the Bloch-Siegert regime, even with linear
response.

Linear spectral shifts beyond those in the Jaynes-Cummings model are often
termed vacuum Bloch-Siegert shifts [47], but below we quantify why there is
nothing uniquely quantum about such resonance shifts in a Hopfield model. This
is in contrast to the QRM, which becomes uniquely quantum in nature in the USC
regime.

In classical electromagnetism, the bare polarizability volume of an oscillator
is $\alpha(\omega)={A_{0}\omega_{0}}/{(\omega_{0}^{2}-\omega^{2})}$, with
$A_{0}={2d^{2}}/{\epsilon_{0}}$ and $d$ the dipole moment. Considering the
emitter position ${\bf r}_{0}$, the photonic Green function, under a single
mode expansion (and assuming scalar fields) is [48]

$G_{c}({\bf r}_{0},{\bf
r}_{0},\omega)=\frac{A_{c}\omega^{2}}{\omega_{c}^{2}-\omega^{2}},$ (4)

where $A_{c}=1/V_{\rm eff}\epsilon_{b}$, with $V_{\rm eff}$ the effective mode
volume and $\epsilon_{b}$ the background dielectric constant. Embedding the HO
dipole in the cavity, we obtain the exact polarizability:

$\alpha(\omega)=\frac{A_{0}\omega_{0}}{\omega_{0}^{2}-\omega^{2}-(\omega_{0}/\omega_{c})4g^{2}\omega^{2}/(\omega_{c}^{2}-\omega^{2})},$
(5)

where $4g^{2}=d^{2}\omega_{c}/(2\epsilon_{0}V_{\rm eff}\epsilon_{b})$.
Considering the on-resonance case again, the classical poles are
$\omega_{\pm}^{G}=\omega_{\pm}$ [as in Eq. (2)], so it is identical to the
solution of the lossless Hopfield model. This is consistent with experimental
results on ultrastrongly-coupled molecular vibrational dipoles in IR cavities,
where the same classical-quantum correspondence was observed in the oscillator
frequencies [49]. In a quantum picture, the blueshift is caused by the $P^{2}$
term (or $A^{2}$ term in the Coulomb gauge). However, in a classical picture,
this blueshift is caused from the poles of the cavity-renormalized
polarizability; thus, there is nothing uniquely quantum about this spectral
blueshift. The blueshift is caused by counter rotating-wave effects, in both
models.

This correspondence with the poles of the quantum Hopfield model and classical
electromagnetism is partly known [33, 34], yet sometimes not recovered in
quantum models. Moreover, in a linear optical material, the classical Green
function of the hybrid system must have poles in the upper complex plane [50],
even with linear gain [51]. The classical pole correspondence does not mean
that there are no unique quantum effects in the Hopfield model, since the
ground and excited states are squeezed states [31, 52]. Indeed, for
$\omega_{c}=\omega_{0}$, the quantum ground state, $\ket{0_{+}0_{-}}$ has
energy:
$\omega_{0,0}{=}(\omega_{+}{+}\omega_{-})/2{-}\omega_{0}=\omega_{0}(1+\eta^{2})^{0.5}{-}\omega_{0}>0$.
The ground state contains virtual photons and is an entangled state [53].
Despite this, there appears to be no unique quantum effects that affect the
polariton eigenfrequencies.

What is less well known is to what degree the predicted optical spectra agree
(or not) between the quantum and classical coupled mode theories in the USC
regime, and how to describe such a regime with open-system master equations.
Since all cavity-QED systems have dissipation and input-output channels, it is
essential to model them as open quantum systems. Within the RWA, the vacuum
Rabi doublets are well described classically or quantum mechanically [4], for
both boson and TLS emitters. In the USC regime, things are much more subtle,
and the quantum models have technical problems related to how to properly
include dissipation as well as gauge correction terms (caused by material and
cavity mode truncation).

In a classical light-matter model, we can include a heuristic cavity decay
rate, $\kappa$, in the cavity Green function, and derive the the cavity-
emitted spectrum as

$S^{\rm Class}=F({\bf
R})\left|\frac{E_{0}\,g^{2}\omega^{2}}{(\omega^{2}-\omega_{c}^{2}-i\omega\kappa)(\omega^{2}-\omega_{0}^{2})-4g^{2}\omega^{2}}\right|^{2},$
(6)

where $E_{0}$ is the excitation field strength and $F({\bf R})$ is a geometric
factor. The solution is non-Markovian, causal, and contains no RWA [48]. Also,
this phenomenological approach to dissipation ensures a symmetric spectrum
outside of the USC regime for a resonant cavity and TLS.

In a quantum picture, to include cavity dissipation, we use an open-system
approach [54, 40], at the level of a generalized master equation (GME) [55,
40, 56, 57],

$\frac{\rm{d}}{\rm{dt}}\rho=-\frac{i}{\hbar}[H_{\rm
S},\rho]+\mathcal{L}_{\rm{G}}\rho+\frac{P_{c}}{2}\mathcal{D}[X^{-}]\rho,$ (7)

where $P_{c}$ is an incoherent pump term, with ${\cal D}[O]\rho=2O\rho
O^{\dagger}-\rho O^{\dagger}O-O^{\dagger}O\rho$, and the cavity dissipator is

$\displaystyle\mathcal{L}_{\rm{G}}\rho$
$\displaystyle=\frac{1}{2}\sum\limits_{\omega,\omega^{\prime}>0}\Gamma_{c}(\omega)[X^{+}(\omega)\rho
X^{-}(\omega^{\prime})-X^{-}(\omega^{\prime})X^{+}(\omega)\rho]$
$\displaystyle+\Gamma_{c}(\omega^{\prime})[X^{+}(\omega)\rho
X^{-}(\omega^{\prime})-\rho X^{-}(\omega^{\prime})X^{+}(\omega)].$ (8)

The dressed-state operators, $X^{\pm}$, are defined from

$X^{+}(\omega)=\braket{j}{\Pi_{c}}{k}\ket{j}\bra{k},$ (9)

where $\omega=\omega_{k}-\omega_{j}>0$, $X^{-}=(X^{+})^{\dagger}$, and
$\Pi_{c}$ is a cavity operator linear in the photon creation and destruction
operators. We neglect atom/emitter decay channels, since these are typically
negligible. The cavity decay rates are obtained from $\Gamma_{c}(\omega)=2\pi
J_{c}(\omega)$, where $J_{c}(\omega)$ is the spectral bath function. Below, we
use $\Gamma_{c}=\kappa$, and show that this is sufficient to recover the
classical spectral form if the appropriate $\Pi_{c}$ operator can be
identified.

Figure 2:  Hopfield GME results versus the classical model for $\eta=0.5$, and
$\kappa=0.05g$, with three different bath models, (a) $\Pi_{c}=P$, (b)
$\Pi_{c}=Q$, and (c) $\Pi_{c}=(P+Q)/\sqrt{2}$ ($\Pi_{c}=(P-Q)/\sqrt{2}$ gives
identical results). Gauge-corrected results use $X^{\pm}_{\rm GC}$ and primed
cavity operators for $\Pi_{c}$. Only model (c), with gauge corrections (‘GC’),
overlaps with the classical solution.  Figure 3: Dissipative (GME) Hopfield
model (a-c) and QRM (d-f) versus the classical solution, for three values of
$\eta$. We only present the $\Pi_{c}=(P+Q)/\sqrt{2}$ bath modeland show GME
results, with and without gauge corrections. Once again, only the Hopfield
model with gauge corrections overlaps with the classical solution. At higher
values of $\eta$, the QRM also shows multiple resonances, and these are also
substantially different with gauge corrections. The QRM fails to recover the
classical solution in the USC regime.

The precise form of $\Pi_{c}$ matters in the USC regime [40]. For example, one
could choose $\Pi_{c}=i(a^{\dagger}-a)\equiv P$, or
$\Pi_{c}=a^{\dagger}+a\equiv Q$, and obtain significantly different
predictions, or any linear combination of these two. This is not the case with
a RWA. Furthermore, in the USC regime, there a gauge ambiguity for the
electric field operator [37], because $P$ represents the Coulomb gauge
electric field, and we are using a system Hamiltonian in the dipole gauge. For
a restricted TLS subspace, this ambiguity is corrected through the
transformation $a^{\prime}\rightarrow{\cal U}a{\cal
U}^{\dagger}=a+i\eta\sigma_{x}$ [58], where ${\cal
U}=\exp(-i\eta(a+a^{\dagger})\sigma_{x})$ is the projected unitary operator
[37, 59, 38], with $\sigma_{x}=b+b^{\dagger}$ (Hopfield model) or
$\sigma_{x}=\sigma^{+}+\sigma^{-}$ (QRM). Thus, one must use $a^{\prime}$ and
$a^{\prime\dagger}$ in the computation of the dissipators to ensure gauge
invariance. The fact that the $\Pi_{c}$ operator should consist only of
bosonic $a$ and $a^{\dagger}$ operators in the Coulomb gauge is a consequence
of photon loss being associated only with electromagnetic degrees of freedom,
and not the TLS [38].

In the USC regime, the system has transition operators $\ket{j}\bra{k}$ which
cause transitions between the dressed eigenstates of the system
$\\{\ket{j},\ket{k}\\}$. For the cavity mode operator, these transitions are
obtained from the dressed operators $X^{\pm}$, and again these must be gauge
corrected. To make the notation clearer, we can use $X^{\pm}_{\rm GC}$ to
indicate that we are applying $\Pi_{c}$ operators with gauge corrections.
Thus, the cavity-emitted quantum spectrum is obtained from

$S^{\rm QM}\propto\text{Re}\left[\int_{0}^{\infty}d\tau
e^{i\omega\tau}\int_{0}^{\infty}\left\langle X_{\rm GC}^{-}(t)X_{\rm
GC}^{+}(t+\tau)\right\rangle dt\right],$ (10)

and calculations without gauge corrections simply use $X^{\pm}$, i.e., without
primed cavity operators in the computation of the dissipators and cavity-mode
obervables. We will show both solutions to better highlight the role of these
gauge corrections, and also show how they are required to recover classical
correspondence. In all GME calculations below, we use $P_{c}\ll\kappa$, to
ensure weak excitation, and the numerical results are carried out using Python
and QuTiP [60, 61].

It is important to stress that our gauge-corrected results are necessary to
ensure gauge invariance. For example, we could also use a Hopfield model in
the Coulomb gauge, where $H_{\rm Hop}^{\rm
GC}=\omega_{c}a^{\dagger}a+\omega_{0}b^{\dagger}b+ig\frac{\omega_{0}}{\omega_{c}}(b^{\dagger}-b)(a+a^{\dagger})+D(a+a^{\dagger})^{2}$,
and then we could use unprimed operators for the cavity mode, where indicated
above, specifically in $X^{\pm}$ and $P$. Note also that $D$ must be the same
in both gauges to ensure gauge invariance, which has been proven also for a
dilute Dicke model [34]. For the QRM, the gauge-corrected system Hamiltonian
in the Coulomb gauge is [37], ${H}_{\rm QRM}^{\rm
CG}=\omega_{c}a^{\dagger}a+\frac{\omega_{0}}{2}\\{\sigma_{z}\cos(2\eta(a+a^{\dagger}))+\sigma_{y}\sin(2\eta(a+a^{\dagger}))\\}$,
which contains field operators to all orders. For simplicity, we use only the
dipole gauge below, but we have checked that all results below are identical
in the Coulomb gauge.

Computed spectra.—Figure 2, shows the classical and quantum solutions for the
dissipative Hopfield model, with three types of bath coupling models: (a)
$\Pi_{c}=P$, (b) $\Pi_{c}=Q$, and (c) $\Pi_{c}=(P\pm Q)/\sqrt{2}$, where
primed indices are used for the gauge-corrected models. We first choose
$\eta=0.5$ here, which is well into the USC regime, with $\kappa=0.05g$. As
recognized, we find very good agreement with the classical solution only when
$1/\sqrt{2}(P\pm Q)$ and only with gauge corrections. To the best of our
knowledge, this is the first time that such a solution and correspondence has
been made, and these results also demonstrate the significant problem with
choosing an arbitrary system-bath interaction form in the USC regime. It
should be noted that such a correspondence does not necessarily indicate that
this is the _correct_ choice of dissipation model. Rather, this result allows
for unambiguous connection and comparison between quantum and classical
heuristic models of dissipation, and is further evidence of the limitations of
purely phenomenological approaches, when treating losses in open quantum
systems.

Next, we look at the role of this bath coupling model, for different $\eta$,
using $\Pi_{c}=(P\pm Q)/\sqrt{2}$, for both the Hopfield model and the QRM.
The spectral calculations are shown in Fig. 3, along with the classical
solution. We see that the Hopfield model and QRM differ substantially in all
USC regimes, and the QRM takes on multiple resonances when $\eta$ is
sufficiently large, even for weak excitation, as well as pronounced spectral
asymmetries. Moreover, we find that the dissipative Hopfield model, with
$\Pi_{c}=(P+Q)/\sqrt{2}$, agrees very well with the classical two oscillator
model at all coupling regimes shown. This is clearly not the case for the QRM,
and such substantial differences should be easy to identify in experiments.

Discussions and summary.—We have shown how the optical spectra for a
dissipative Hopfield model in USC can be described quantum mechanically or
classically. To achieve correspondence with the classical dissipative result,
quantum models must properly respect gauge invariance and implement the
appropriate bath coupling operator. Without such a correspondence, any open-
system master equations in this regime with ad hoc system-bath interactions
are ambiguous and can predict wildly differing spectra.

We have also clarified how the dissipative Hopfield model and QRM
substantially differ, even for weak excitation, at all USC regimes, including
the perturbative Bloch-Siegert regime. Thus, only the Hopfield model yields a
classical correspondence under linear response, and this correspondence only
occurs with a careful treatment of quantum dissipation and gauge corrections.
While we used a normal mode expansion with heuristic broadening, this form is
well established for high $Q$ cavities outside the USC regime, and future work
could improve such models using classical and quantized quasinormal mode
theories [42, 43]. One possible clue as to the significance of the
$\Pi_{c}\propto P\pm Q$ coupling can be seen by noting that the classical
phenomenological loss model is phase-insensitive. In the quantum loss model, a
choice of $P\pm Q$ for the quadrature coupling to the bath is the only choice
which gives equal coupling magnitude to each quadrature (i.e., is phase-
insensitive in magnitude).

Broadly, these findings are important for a wide class of light-matter systems
now emerging to study the USC regime, including lossy Landau systems and
metallic systems [32, 33]. Apart from showing a direct classical
correspondence for dissipative modes, our results can be used to guide open-
system quantum models that are needed when observations are uniquely quantum
in nature, e.g., in the QRM for any excitation including coherent excitation,
and the Hopfield model excited with non-classical fields.

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC), the National Research Council of Canada (NRC), the
Canadian Foundation for Innovation (CFI), and Queen’s University, Canada. S.H.
acknowledges the Japan Society for the Promotion of Science (JSPS) for funding
support through an Invitational Fellowship. F.N. is supported in part by
Nippon Telegraph and Telephone Corporation (NTT) Research, the Japan Science
and Technology Agency (JST) [via the Quantum Leap Flagship Program (Q-LEAP),
and the Moonshot R&D Grant No. JPMJMS2061], the Asian Office of Aerospace
Research and Development (AOARD) (via Grant No. FA2386-20-1-4069), and the
Foundational Questions Institute Fund (FQXi) via Grant No. FQXi-IAF19-06. We
thank Jun Kono and Hideo Mabuchi for useful discussions.

## References

  * Agarwal [1984] G. S. Agarwal, Vacuum-field rabi splittings in microwave absorption by rydberg atoms in a cavity, Phys. Rev. Lett. 53, 1732 (1984).
  * Boca  _et al._ [2004] A. Boca, R. Miller, K. M. Birnbaum, A. D. Boozer, J. McKeever, and H. J. Kimble, Observation of the vacuum Rabi spectrum for one trapped atom, Phys. Rev. Lett. 93, 233603 (2004).
  * Andreani  _et al._ [1999] L. C. Andreani, G. Panzarini, and J.-M. Gérard, Strong-coupling regime for quantum boxes in pillar microcavities: Theory, Phys. Rev. B 60, 13276 (1999).
  * Zhu  _et al._ [1990] Y. Zhu, D. J. Gauthier, S. E. Morin, Q. Wu, H. J. Carmichael, and T. W. Mossberg, Vacuum rabi splitting as a feature of linear-dispersion theory: Analysis and experimental observations, Phys. Rev. Lett. 64, 2499 (1990).
  * Khitrova  _et al._ [1999] G. Khitrova, H. M. Gibbs, F. Jahnke, M. Kira, and S. W. Koch, Nonlinear optics of normal-mode-coupling semiconductor microcavities, Rev. Mod. Phys. 71, 1591 (1999).
  * Hughes and Kamada [2004] S. Hughes and H. Kamada, Single-quantum-dot strong coupling in a semiconductor photonic crystal nanocavity side coupled to a waveguide, Phys. Rev. B 70, 195313 (2004).
  * Cazé  _et al._ [2013] A. Cazé, R. Pierrat, and R. Carminati, Strong coupling to two-dimensional Anderson localized modes, Phys. Rev. Lett. 111, 053901 (2013).
  * Herrera and Spano [2016] F. Herrera and F. C. Spano, Cavity-Controlled Chemistry in Molecular Ensembles, Physical Review Letters 116, 238301 (2016).
  * Flick  _et al._ [2017] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry, Proceedings of the National Academy of Sciences 114, 3026 (2017).
  * Reithmaier  _et al._ [2004] J. P. Reithmaier, G. Sęk, A. Löffler, C. Hofmann, and S. o. Kuhn, Strong coupling in a single quantum dot–semiconductor microcavity system, Nature 432, 197 (2004).
  * Yoshie  _et al._ [2004] T. Yoshie, A. Scherer, J. Hendrickson, _et al._ , Vacuum Rabi splitting with a single quantum dot in a photonic crystal nanocavity, Nature 432, 200 (2004).
  * Gu  _et al._ [2017] X. Gu, A. F. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, Microwave photonics with superconducting quantum circuits, Physics Reports 718-719, 1 (2017).
  * Cao  _et al._ [2010] X. Cao, J. Q. You, H. Zheng, A. Kofman, and F. Nori, Dynamics and quantum Zeno effect for a qubit in either a low-or high-frequency bath beyond the rotating-wave approximation, Physical Review A 82, 022119 (2010).
  * Schuster  _et al._ [2008] I. Schuster, A. Kubanek, A. Fuhrmanek, T. Puppe, P. W. H. Pinkse, K. Murr, and G. Rempe, Nonlinear spectroscopy of photons bound to one atom, Nature Physics 4, 382 (2008).
  * Bose  _et al._ [2014] R. Bose, T. Cai, K. R. Choudhury, G. S. Solomon, and E. Waks, All-optical coherent control of vacuum Rabi oscillations, Nature Photonics 8, 858 (2014).
  * Fink  _et al._ [2008] J. M. Fink, M. Göppl, M. Baur, R. Bianchetti, P. J. Leek, A. Blais, and A. Wallraff, Climbing the Jaynes–Cummings ladder and observing its nonlinearity in a cavity QED system, Nature 454, 315 (2008).
  * Bishop  _et al._ [2008] L. S. Bishop, J. M. Chow, J. Koch, _et al._ , Nonlinear response of the vacuum Rabi resonance, Nature Physics 5, 105 (2008).
  * Illes  _et al._ [2015] E. Illes, C. Roy, and S. Hughes, Spectral multiphoton effects and quantum anharmonicities in dissipative cavity-QED systems via off-resonant coherent excitation, Optica 2, 689 (2015).
  * Hamsen  _et al._ [2017] C. Hamsen, K. N. Tolazzi, T. Wilk, and G. Rempe, Two-Photon Blockade in an Atom-Driven Cavity QED System, Physical Review Letters 118, 133604 (2017).
  * Milonni and Smith [1975] P. W. Milonni and W. A. Smith, Radiation reaction and vacuum fluctuations in spontaneous emission, Phys. Rev. A 11, 814 (1975).
  * Milonni [1976] P. Milonni, Semiclassical and quantum-electrodynamical approaches in nonrelativistic radiation theory, Physics Reports 25, 1 (1976).
  * Omelyanchouk  _et al._ [2008] A. N. Omelyanchouk, S. N. Shevchenko, _et al._ , Pseudo-Rabi oscillations in superconducting flux qubits in the classical regime, Phys. Rev. B 78, 054512 (2008).
  * Ciuti and Carusotto [2006] C. Ciuti and I. Carusotto, Input-output theory of cavities in the ultrastrong coupling regime: The case of time-independent cavity parameters, Phys. Rev. A 74, 033811 (2006).
  * Anappara  _et al._ [2009] A. A. Anappara, S. De Liberato, A. Tredicucci, _et al._ , Signatures of the ultrastrong light-matter coupling regime, Physical Review B 79, 201303(R) (2009).
  * Beaudoin  _et al._ [2011] F. Beaudoin, J. M. Gambetta, and A. Blais, Dissipation and ultrastrong coupling in circuit QED, Physical Review A 84, 043832 (2011).
  * Zueco  _et al._ [2009] D. Zueco, G. M. Reuther, S. Kohler, and P. Hänggi, Qubit-oscillator dynamics in the dispersive regime: Analytical theory beyond the rotating-wave approximation, Phys. Rev. A 80, 033846 (2009).
  * De Bernardis  _et al._ [2018] D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Breakdown of gauge invariance in ultrastrong-coupling cavity QED, Physical Review A 98, 053819 (2018).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nature Reviews Physics 1, 19 (2019).
  * Forn-Díaz  _et al._ [2019] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong coupling regimes of light-matter interaction, Reviews of Modern Physics 91, 025005 (2019).
  * Boité [2020] A. L. Boité, Theoretical methods for ultrastrong light–matter interactions, Advanced Quantum Technologies 3, 1900140 (2020).
  * Liberato  _et al._ [2007] S. D. Liberato, C. Ciuti, and I. Carusotto, Quantum vacuum radiation spectra from a semiconductor microcavity with a time-modulated vacuum rabi frequency, Phys. Rev. Lett. 98, 103602 (2007).
  * Zhang  _et al._ [2016] Q. Zhang, M. Lou, X. Li, J. L. Reno, W. Pan, J. D. Watson, M. J. Manfra, and J. Kono, Collective non-perturbative coupling of 2d electrons with high-quality-factor terahertz cavity photons, Nature Physics 12, 1005 (2016).
  * Rajabali  _et al._ [2021] S. Rajabali, E. Cortese, M. Beck, S. De Liberato, J. Faist, and G. Scalari, Polaritonic nonlocality in light–matter interaction, Nature Photonics 15, 690 (2021).
  * Garziano  _et al._ [2020] L. Garziano, A. Settineri, O. Di Stefano, _et al._ , Gauge invariance of the Dicke and Hopfield models, Phys. Rev. A 102, 023718 (2020).
  * Bouchet and Carminati [2019] D. Bouchet and R. Carminati, Quantum dipole emitters in structured environments: a scattering approach: tutorial, Journal of the Optical Society of America A 36, 186 (2019).
  * Jørgensen and Wubs [2022] M. A. Jørgensen and M. Wubs, Quantifying the breakdown of the rotating-wave approximation in single-photon superradiance, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 195401 (2022).
  * Di Stefano  _et al._ [2019] O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics, Nature Physics 15, 803 (2019).
  * Gustin  _et al._ [2023] C. Gustin, S. Franke, and S. Hughes, Gauge-invariant theory of truncated quantum light-matter interactions in arbitrary media, Phys. Rev. A 107, 013722 (2023).
  * Bamba and Ogawa [2014] M. Bamba and T. Ogawa, Recipe for the Hamiltonian of system-environment coupling applicable to the ultrastrong-light-matter-interaction regime, Phys. Rev. A 89, 023817 (2014).
  * Salmon  _et al._ [2022] W. Salmon, C. Gustin, A. Settineri, O. D. Stefano, D. Zueco, S. Savasta, F. Nori, and S. Hughes, Gauge-independent emission spectra and quantum correlations in the ultrastrong coupling regime of open system cavity-QED, Nanophotonics 11, 1573 (2022).
  * Baranov  _et al._ [2020] D. G. Baranov, B. Munkhbat, E. Zhukova, A. Bisht, A. Canales, B. Rousseaux, G. Johansson, T. J. Antosiewicz, and T. Shegai, Ultrastrong coupling between nanoparticle plasmons and cavity photons at ambient conditions, Nature Communications 11 (2020).
  * Franke  _et al._ [2019] S. Franke, S. Hughes, M. K. Dezfouli, _et al._ , Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics, Phys. Rev. Lett. 122, 213901 (2019).
  * Franke  _et al._ [2020] S. Franke, M. Richter, J. Ren, _et al._ , Quantized quasinormal-mode description of nonlinear cavity-QED effects from coupled resonators with a Fano-like resonance, Phys. Rev. Research 2, 033456 (2020).
  * Savasta  _et al._ [2020] S. Savasta, O. D. Stefano, and F. Nori, Thomas–Reiche–Kuhn (TRK) sum rule for interacting photons, Nanophotonics 10, 465 (2020).
  * Taylor  _et al._ [2022] M. A. D. Taylor, A. Mandal, and P. Huo, Resolving ambiguities of the mode truncation in cavity quantum electrodynamics, Opt. Lett. 47, 1446 (2022).
  * Emary and Brandes [2003] C. Emary and T. Brandes, Quantum Chaos Triggered by Precursors of a Quantum Phase Transition: The Dicke Model, Phys. Rev. Lett. 90, 044101 (2003).
  * Li  _et al._ [2018] X. Li, M. Bamba, Q. Zhang, S. Fallahi, G. C. Gardner, W. Gao, M. Lou, K. Yoshioka, M. J. Manfra, and J. Kono, Vacuum Bloch–Siegert shift in Landau polaritons with ultra-high cooperativity, Nature Photonics 12, 324 (2018).
  * SI [2023] (2023), see Supplemental Material at [URL will be inserted by publisher] for more detailed information, including further details on the classical scatteering theory and spectrum, the role of the cavity bath spectrum function, and influence of cavity pump rate, where we also reference papers [62, 63, 64, 65, 66].
  * George  _et al._ [2016] J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura, C. Genet, and T. W. Ebbesen, Multiple Rabi Splittings under Ultrastrong Vibrational Coupling, Phys. Rev. Lett. 117, 153601 (2016).
  * Bamba and Ogawa [2013] M. Bamba and T. Ogawa, System-environment coupling derived by Maxwell’s boundary conditions from the weak to the ultrastrong light-matter-coupling regime, Phys. Rev. A 88, 013814 (2013).
  * Franke  _et al._ [2021] S. Franke, J. Ren, M. Richter, A. Knorr, and S. Hughes, Fermi’s golden rule for spontaneous emission in absorptive and amplifying media, Phys. Rev. Lett. 127, 013602 (2021).
  * Makihara  _et al._ [2021] T. Makihara, K. Hayashida, G. Timothy Noe II, _et al._ , Ultrastrong magnon–magnon coupling dominated by antiresonant interactions, Nature Communications 12 (2021).
  * Miao and Agarwal [2023] Q. Miao and G. S. Agarwal, Polaritonic ultrastrong coupling: Quantum entanglement in ground state (2023), arXiv:2304.00680 .
  * Carmichael [2013] H. J. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ (Springer Science & Business Media, 2013).
  * Settineri  _et al._ [2018] A. Settineri, V. Macrí, A. Ridolfo, _et al._ , Dissipation and thermal noise in hybrid quantum systems in the ultrastrong-coupling regime, Physical Review A 98, 053834 (2018).
  * Mercurio  _et al._ [2022] A. Mercurio, V. Macrì, C. Gustin, S. Hughes, S. Savasta, and F. Nori, Regimes of cavity QED under incoherent excitation: From weak to deep strong coupling, Phys. Rev. Res. 4, 023048 (2022).
  * Akbari  _et al._ [2023] K. Akbari, W. Salmon, F. Nori, and S. Hughes, Generalized Dicke model and gauge-invariant master equations for two atoms in ultrastrongly-coupled cavity quantum electrodynamics, Phys. Rev. Res. 5, 033002 (2023).
  * Settineri  _et al._ [2021] A. Settineri, O. Di Stefano, D. Zueco, S. Hughes, S. Savasta, and F. Nori, Gauge freedom, quantum measurements, and time-dependent interactions in cavity QED, Phys. Rev. Res. 3, 023079 (2021).
  * Savasta  _et al._ [2021] S. Savasta, O. Di Stefano, A. Settineri, D. Zueco, S. Hughes, and F. Nori, Gauge principle and gauge invariance in two-level systems, Phys. Rev. A 103, 053703 (2021).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013).
  * de Vries  _et al._ [1998] P. de Vries, D. V. van Coevorden, and A. Lagendijk, Point scatterers for classical waves, Rev. Mod. Phys. 70, 447 (1998).
  * Van Vlack [2012] C. P. Van Vlack,  _Dyadic Green Functions and their Applications_ , Ph.D. thesis, Queen’s University, Canada (2012).
  * Kristensen  _et al._ [2011] P. T. Kristensen, J. Mørk, P. Lodahl, and S. Hughes, Decay dynamics of radiatively coupled quantum dots in photonic crystal slabs, Phys. Rev. B 83, 075305 (2011).
  * Laucht  _et al._ [2009] A. Laucht, N. Hauke, J. M. Villas-Bôas, F. Hofbauer, G. Böhm, M. Kaniber, and J. J. Finley, Dephasing of exciton polaritons in photoexcited ingaas quantum dots in gaas nanocavities, Phys. Rev. Lett. 103, 087405 (2009).
  * Yao  _et al._ [2010] P. Yao, P. K. Pathak, E. Illes, _et al._ , Nonlinear photoluminescence spectra from a quantum-dot–cavity system: Interplay of pump-induced stimulated emission and anharmonic cavity QED, Phys. Rev. B 81, 033309 (2010).