Finite-size scaling of the photon-blockade breakdown dissipa-
tive quantum phase transition
A. Vukics
1
, A. Dombi
1
, J. M. Fink
2
, and P. Domokos
1
1
Wigner Research Centre for Physics, H-1525 Budapest, P.O. Box 49., Hungary
2
Institute of Science and Technology Austria, 3400 Klosterneuburg, Austria
May 29, 2019
We prove that the observable telegraph sig-
nal accompanying the bistability in the photon-
blockade-breakdown regime of the driven and
lossy Jaynes–Cummings model is the finite-size
precursor of what in the thermodynamic limit is
a genuine first-order phase transition. We con-
struct a finite-size scaling of the system param-
eters to a well-defined thermodynamic limit, in
which the system remains the same microscopic
system, but the telegraph signal becomes macro-
scopic both in its timescale and intensity. The
existence of such a finite-size scaling completes
and justifies the classification of the photon-
blockade-breakdown effect as a first-order dis-
sipative quantum phase transition.
1 Introduction
A number of recent theoretical proposals and experi-
ments refer to a new type of dissipative phase transi-
tion in a quantum system, that is said to be of first
order. The corresponding systems include cluster for-
mation of Rydberg-state atoms, described by Ising-type
spin models [1, 2, 3, 4, 5, 6], the vacuum-superfluid
transition in a dissipative Hubbard model [7], arrays of
nonlinear cavity systems [8], and the ultimately simple
model system of a single semiconductor cavity polariton
mode [9, 10, 11, 12]. Common to these systems is that
the signature of the first-order transition is the bistable
behaviour of an observable. Typically, the bistability
is manifested by a telegraph-like signal, i.e., randomly
alternating ‘bright’ and ‘dark’ periods of a radiation
mode, associated with the mixture of two phases of the
quantum system. This temporal signal can be readily
monitored by a ‘classical’ detector. For example, this
is the case in the photon-blockade breakdown transi-
tion [13, 14, 15], which has recently been experimentally
demonstrated [16]. However, such a signal can be ob-
served in other scenarios that are not regarded as phase
transitions. Think of, for example, the famous quantum
A. Vukics: vukics.andras@wigner.mta.hu
jump experiment [17] measuring the fluorescence of a
driven V-type three-level atom. The fluorescence signal
also exhibits long bright and dark periods, separated by
single quantum jumps. The ‘dark’ period indicates the
atom being in a long-lived, metastable state, whereas a
spontaneous emission event can take it into the ground
state, where the atom is resonantly driven by a laser on
an electric-dipole-allowed ‘bright’ transition. What is
then the distinctive feature of a first-order phase tran-
sition?
This paper is devoted to complete and justify the in-
terpretation of the photon-blockade-breakdown effect as
a first-order dissipative quantum phase transition. To
this end, one needs to introduce the concepts of ther-
modynamic limit and finite-size scaling for the micro-
scopic system of a driven dissipative Jaynes-Cummings
model. The thermodynamic limit will be defined for this
microscopic system in such a way that the number of
relevant degrees of freedom remains fixed. Instead of
growing the system size, we construct a scaling of the
parameters which keeps the form of the stationary so-
lution of the driven-dissipative system invariant. At the
same time, the proposed finite-size scaling leads to an
increasing robustness of the attractor states associated
with the bistability signal, until these states reach full
stability in the thermodynamic limit. On approaching
this limit, at variance with the fluorescence shelving ex-
periment, no single quantum jump or other microscopic
event can flip the system from one phase to the other.
In the thermodynamic limit, the blinking telegraph-like
signal vanishes completely and the state of the system
is determined by the initial condition, similarly to the
usual hysteresis behaviour in classical critical systems.
If such a finite size scaling is possible – and here we show
that this is the case for the photon-blockade breakdown
effect –, then the bistability that can be observed in
a given experimental realization of the system with its
finite parameters not in the thermodynamic limit, can
be considered the finite-size approximation of what is
a genuine first-order phase transition in the thermody-
namic limit.
The paper is structured as follows. In section 2 we
Accepted in Quantum 2019-05-24, click title to verify 1
arXiv:1809.09737v2 [quant-ph] 27 May 2019
clarify the basic notions of dynamical quantum phase
transitions in dissipative systems, in particular, the con-
cept of a first-order transition. In section 3 we recall the
model and the basic ingredients of the photon-blockade-
breakdown effect, together with the equations governing
the underlying classical phase transition. In section 4 we
present the time evolution of the system in steady state,
which we show to be a telegraph signal. We undertake a
detailed numerical analysis of timescales and the statis-
tics of dwell times in the attractors. In section 5, the
finite-size scaling of the parameters is determined, and
section 6 is devoted to the analysis of the stability of
phases. We shortly comment on the role of atomic de-
cay in section 7 and then conclude in section 8.
2 What is a first-order dissipative quan-
tum phase transition?
Quantum phase transitions (QPTs) are abrupt,
symmetry-breaking changes in the quantum state of a
large quantum system at a critical value of an exter-
nal control parameter. The QPT takes place at zero
temperature and thus it generally refers to a change of
the ground state [18]. The ground states on the two
sides of the critical point have different symmetries and
correspond to different macroscopic observables when
a quantum measurement is performed. The concept of
QPT can be straightforwardly extended to dissipative
(open) quantum systems. Assume that the system is ex-
cited by an external coherent field, e.g., by laser irradi-
ation. The coherent driving competing with the various
damping processes sets the system in a dynamical equi-
librium which is not the ground state, moreover, it may
be far from a thermal state due to energy flow through
the system. Clearly, the steady state will depend on ex-
ternal parameters. Just like the ground state of closed
systems, the stationary state of driven-dissipative quan-
tum systems can also exhibit non-analytic behaviour as
a function of some control parameter [19, 20, 21]. All
this can happen at zero temperature. The critical point
is marked by diverging quantum fluctuations and cor-
responds to a symmetry breaking change of the steady-
state [22, 23]. Note that the presence of dissipation in
the open system amounts to information leakage to the
environment. Therefore, if such a transition occurs in
the quantum system, the different steady states under
continuous measurement are associated with different
measured, ‘classical’ outputs. In this sense the quantum
states are ‘amplified’ to phases with different macro-
scopic observables even in the case of a small system
with only few degrees of freedom.
Besides continuous phase transitions, theoretical and
experimental works pointed out recently that first-order
phase transitions can also have their counterpart in
open quantum systems [13, 16]. First-order phase tran-
sitions are benchmarked by the coexistence of phases.
Instead of a critical point separating two phases, there
is a critical region in which bistability and hysteresis
occur [24]. How can the coexistence of ‘macroscopically
distinct’ states be accommodated into quantum theory?
The state of a driven-dissipative system is expressed by
a density operator. The steady-state density operator,
which obeys a linear algebraic equation, must be a con-
tinuous function of the external control parameter of the
system. However, the density operator can represent a
bimodal distribution, a mixture of two macroscopically
distinct components. There can be then a ‘rift’ sepa-
rating the two ‘phases’, meanwhile the transition of the
density matrix can occur smoothly from one phase to
the other when tuning the control parameter. It is the
weights in the mixture which vary continuously as a
function of the control parameter: one of the compo-
nents vanishes gradually with the simultaneous growing
of the weight of the other. The complete transition can
happen in a finite range of the control parameter, which
is the bistability range.
The starting point of our study is that the above pic-
ture of 1st order QPT is, however, not the full story.
In particular, the bimodal steady-state density opera-
tor does not tell us about the stability of the phases. It
allows only for identifying the two phases and their rel-
ative weights. In the bistability range, the system con-
tinuously monitored in its dissipation channel manifests
a randomly blinking signal between the distinct values
associated with the two phases. The same steady-state
density operator can correspond to different characteris-
tic dwelling times in the blinking signal. This is the goal
of unraveling into quantum trajectories, and performing
a finite-size upscaling in the language of these trajecto-
ries. The steady-state density matrix is kept invariant
under the upscaling of the system parameters, only with
increasing intensity in the bright periods, meanwhile the
stability of the phases increases, that is, the characteris-
tic dwell time in the attractors tends to infinity. In the
following, we will perform this analysis in the case of
the photon-blockade-breakdown phase transition.
3 The photon-blockade-breakdown
phase transition in a nutshell
3.1 The driven-lossy Jaynes–Cummings model
The simplest possible system exhibiting the photon-
blockade-breakdown effect is composed of a two-level
system coupled to a harmonic oscillator. The two-level
system can be an atom or artificial atom, whereas the
Accepted in Quantum 2019-05-24, click title to verify 2
Figure 1: (a) In a quantum system consisting of a finite-level
system strongly coupled to a harmonic oscillator, the lower part
of the spectrum is anharmonic, but there are always higher-lying
harmonic subsets of the spectrum. (b) The Jaynes-Cummings
spectrum is a prototype of this behaviour (this part of the fig-
ure is to scale). Panel (i): the anharmonic part of the spectrum.
An excitation tuned to a 1-photon transition misses the second
rung of the ladder by a significant amount. That tuned to a
2- and 3-photon transition, misses the first rung and the first
and second rung, respectively. Panel (ii): closely harmonic part
of the spectrum. An excitation tuned to the transition between
|6, −i and |7, −i is close to resonance also with the transi-
tion between |7, −i and |8, −i (the detuning is invisible on the
figure’s scale).
oscillator can represent a single lossy mode of the radia-
tion field or a longitudinal mode of a stripline resonator.
We describe this interaction within the driven Jaynes-
Cummings model, i.e., using the electric-dipole coupling
and the rotating-wave approximation (RWA):
H = ω
M
a
†
a + ω
A
σ
†
σ + ig
a
†
σ − σ
†
a
+ iη
a
†
e
−iωt
− a e
iωt
, (1)
where ω
M
is the angular frequency of the mode with
boson operator a, ω
A
is that of the atomic transition
with step-down operator σ, and g is the single-photon
Rabi frequency expressing the coupling strength. The
external coherent drive is characterized by the ampli-
tude η, here taken to be real for simplicity, and the
frequency ω. We note that due to the strong coupling,
our treatment is close to the validity limit of the RWA.
A quantum-to-classical transition in a similar system
without the conventional rotating-wave approximation
has been reported in Refs. [25, 26].
Assuming resonance between the mode and the atom,
i.e. ω
M
= ω
A
, and going into a frame rotating at the
drive frequency, one gets a virtually time-independent
Hamiltonian,
H = −∆
a
†
a + σ
†
σ
+ ig
a
†
σ − σ
†
a
+ iη
a
†
− a
,
(2)
where the detuning ∆ = ω −ω
M
is a tunable parameter
of the drive. The mode is that of a high-finesse resonator
and is subject to loss. Similarly, the two-level system can
have decay through spontaneous emission. These inco-
herent processes can be modelled by Liouvillian terms
in the master equation
˙ρ = −i [ H , ρ ] + κ
2aρa
†
− a
†
aρ − ρa
†
a
+ γ
2σρσ
†
− σ
†
σρ − ρσ
†
σ
. (3)
In the rest of the paper we will consider the case κ γ,
most importantly γ = 0. The mode relaxation parame-
ter κ defines the microscopic timescale of the problem.
3.2 The photon-blockade-breakdown effect
For weak drive strengths η g, the excited eigenstates
of the Hamiltonian are close to the Jaynes-Cummings
dressed states
|n, ±i =
1
√
2
(|g, ni ± |e, n − 1i) (4)
with n = 1, 2, . . . , and the energy levels are
E
n,±
= n ω
M
±
√
ng . (5)
In the strong coupling regime g κ, the level shifts
±
√
n g with respect to the bare frequencies exceed sig-
nificantly the linewidth ∼ κ, so the system cannot be
excited out of the ground state |g, 0i by a near-resonant
driving ∆ ≈ 0. This is the photon blockade effect.
1
Above a certain intensity of the driving, however, the
system can be excited via higher-order (multi-photon)
transition processes into the higher-lying part of the
spectrum. Since one of the constituents is a harmonic
oscillator, in the high-lying part there are harmonic sub-
sets of the spectrum (cf. fig. 1), ladders of equidistant
levels which can host a coherent state with large ampli-
tude and well-defined phase. Following a low-probability
multi-photon transition into this range of the spectrum,
the system gets stabilized into such a semiclassical trap-
ping state by the competition of coherent driving and
decay. Such an excitation, i.e., the breakdown of the pho-
ton blockade takes place in the form of a bistability. In
a finite interval of the drive strength, the steady-state
density operator of the system is the mixture of the
‘dark’ and ‘bright’ phases, i.e., the ground state and a
highly excited coherent state of the oscillator [13, 14].
As a main subject of our study, we will analyse this so-
lution of the density matrix in detail starting from the
next section. However, beforehand, we try to capture
the non-analytic behaviour in a corresponding classical
theory.
1
In the literature, the ‘photon blockade’ sometimes denotes
the effect when the first excited state with a single photon can be
excited resonantly, but further excitations are suppressed due to
off-resonance. This is analogous to the effect of Coulomb blockade
to some extent. Here we use the term in a more general sense,
where the system is blocked in the ground state, and no photons
at all can be transferred to the system.
Accepted in Quantum 2019-05-24, click title to verify 3
3.3 Classical phase diagram
Following Carmichael [13], let us first look at what
the Jaynes and Cummings semiclassical (also known as
‘neoclassical’) equations can tell us. These are obtained
by taking expectation values in Heisenberg equations of
motion (e.g. α = hai) and factorizing the expectations
of operator products to obtain a closed set of nonlinear
equations for the expectation values. Taking γ = 0, the
theory leads to the self-consistent equation (valid for
∆ < 0)
|α|
2
N
scale
=
2η
g
2
×
1 +
∆
κ
−
1
p
∆
2
κ
2
/g
4
+ |α|
2
/N
scale
!
2
−1
, (6)
where we introduced the parameter N
scale
= g
2
/4κ
2
.
This nonlinear equation can afford different solution
sets. The various domains are depicted in the phase di-
agram in fig. 2. Below the lower boundary, the phase is
the photon blockade regime, whereas above the higher
boundary, the system is highly excited. In between,
eq. (6) has multivalued solution indicating bistability.
The coordinates were chosen to be the tunable vari-
ables, i.e., the frequency and amplitude of the driving,
which can serve as control parameters of the phase tran-
sition. Only the ∆ < 0 half-plane is shown, the positive-
detuning part being the same mirrored to the ∆ = 0
axis.
The neoclassical result suggests an appropriate ther-
modynamic limit [13] and a corresponding scaling of
system parameters. On fixing the timescale to the mi-
croscopic one, κ = 1, we see that a characteristic photon
number is expected to scale as ∼ g
2
. Hence, the thermo-
dynamic (infinite-system) limit is the strong coupling
limit g → ∞ (in contrast to previously-reported ther-
modynamic or classical limiting cases of quantum phase
transitions in the Jaynes-Cummings and Rabi models
[27, 28, 29, 30]). Simultaneously, the drive amplitude
must also be scaled up. The first guess, cf. the prefactor
on the right hand side of eq. (6), would be η → ∞ such
that η/g is kept invariant. This is why the quantity η/g
is used for the drive amplitude on the vertical axis of the
phase diagram. With this scaling, the lower boundary of
the bistability phase becomes indeed independent of g:
the three curves coincide perfectly. On the other hand,
the upper boundaries reveal a further dependence on g.
That is, the solution sets of the neoclassical equations
are not invariant under the transformation of g, η → ∞
with η/g = const. Later, in sec. 5 we will identify the
non-trivial scaling of η which leads to a properly defined
thermodynamic limit of the system.
Figure 2: Phase diagram of the neoclassical theory based on
a numerical solution of the transcendental eq. (6) for three
different values of the coupling constant g. The low boundaries
closely overlap, the upper boundaries differ for the three values
of g. As the numerics becomes very unstable when approaching
the far-off-resonance closing point of the bistability region, the
closing of the red curves on the left is inferred by extrapolation
(dashed segments). The common closing point on the right
for all g values is the spontaneous dressed-state polarization
critical point. The cyan star denotes the workpoint chosen in
this paper: it is at this detuning and around this drive strength
that we are going to study the bistable solution for several g
values. Magenta dashed line: the inference (8) for the lower
limit of the bistability region.
Nevertheless, since the lower boundary of the neoclas-
sical bistability domain is invariant, and also the upper
boundary does not vary strongly at ∆ = −5κ for the
g values shown in fig. 2 and used in this paper, in the
plots the drive strength η is given in units of g.
Finally, let us make two side remarks. Firstly,
Guti´errez-J´auregui and Carmichael [33] pointed out
that another possible scaling deducible from the mean-
field steady state equations is keeping ∆/g = const.
This would make that in the phase diagram fig. 2, the
upper limiting curves of the bistability domain would
coincide, but the lower ones would differ for different g
values. So, at this point the two scalings ∆/κ = const.
and ∆/g = const. appear as equivalent alternatives. It
can be argued that our choice is more natural in the
sense that the detuning from a (bare) resonance is mea-
sured against the width of that resonance.
Secondly, there is a critical point at η/g = 1/2 specific
to the resonant driving case ∆ = 0, where the lower and
upper limits of the bistability region converge. It sepa-
rates the solution α = 0 with population inversion in-
creasing from −1/2 to 0 from the one with σ
z
= 0 and
increasing α as the drive strength η is increased fur-
ther. This result is in accordance with that of the full
quantum treatment which can be pursued to an analyt-
ical solution in the resonant case [31]. It shows that the
Accepted in Quantum 2019-05-24, click title to verify 4
quasienergies coalesce in this critical point. This criti-
cal behaviour was identified as the spontaneous dressed-
state polarization by Alsing and Carmichael [32].
3.4 Bistability in the quantum solution and in-
tuitive explanations
Numerical simulations of the full quantum problem de-
fined by eq. (3) confirm qualitatively the phase diagram
based on the neoclassical equations. The existence of a
bistability regime close to resonance |∆| g has been
confirmed [13, 14]. The stationary density-operator so-
lution of the master equation is then a statistical mix-
ture of two states: the ‘dim’ state where the field is close
to the vacuum and the atom is in the ground state, and a
‘bright’ state which consists of a highly excited coherent
state of the field (and a completely saturated atom). On
increasing the drive strength η at a fixed detuning, the
relative weight of the two components is continuously
varied such that the probability of the bright component
goes from zero to 1 in a finite range of η [16]. The steady
state is hence a continuous function of the parameters,
however, there is a ‘rift’ between the two components
of the mixture: these are classically discernible states.
The stability of the bright component for high-enough
drive strength can be understood intuitively as a bal-
ance of coherent drive and spontaneous decay in a
harmonic oscillator. The frequency separating adjacent
dressed states |n + 1, −i and |n, −i is ω
M
−g(
√
n + 1 −
√
n) ≈ ω
M
−g/(2
√
n). In the limit of large photon num-
bers n → N, this tends to a harmonic ladder in which
the coherent drive η and the decay κ compete to create
a coherent state with amplitude
α =
η
κ − i
∆ +
g
2
√
N
=⇒ N =
η
2
κ
2
+
∆ +
g
2
√
N
2
,
(7)
where the self-consistent equation for the photon num-
ber N was obtained by taking the absolute value
squared of the amplitude α. The smallest drive strength
for which this equation can be satisfied is in the case of
‘resonance’, i.e. when the expression in the parentheses
in the denominator vanishes: ∆ = −g/(2
√
N). The self-
consistent solution is then N = (η/κ)
2
, from which the
minimum drive strength follows as
η
min
g
'
κ
2|∆|
. (8)
This law is drawn in the thick dashed magenta line in
the phase diagram in fig. 2, and coincides quite accu-
rately with the lower border of the classical bistability
domain. It is remarkable that the solution of a classical
equation obeys a law extracted from intuitive consid-
eration of quantized energy levels. On further increas-
ing the drive strength, the photon number increases,
thereby leading to some detuning ∆ + g/(2
√
N) 6= 0
in eq. (7). However, there may still be a self-consistent
solution. Because the lower part of the harmonic ladder
is missing (that part of the spectrum being the anhar-
monic photon-blockading part), there is no determinis-
tic evolution into such a self-consistent solution. Nev-
ertheless, the displacement operator corresponding to
coherent driving contains multi-photon excitation pro-
cesses, so the excitation into the near-resonant regime
with a self-consistent solution can take place in a prob-
abilistic manner.
The upper limit of the bistability domain cannot be
determined from an argument as simple as the above for
the lower limit. The reason is that the more we increase
the drive strength η, the more the quasi-energy levels,
i.e. the true eigenvalues of the Hamiltonian (1), differ
from the dressed levels of the η = 0 Jaynes-Cummings
model, since they are getting dressed also by the co-
herent drive [31]. However, the analytical form of the
quasi-energy levels in the finite-drive strength case is
not known for ∆ 6= 0, only for the case of ∆ = 0.
Nevertheless, it is clear that the appearance of an η-
dependence of the energy levels makes that the η/g scal-
ing suggested by both eqs. (6) and (7) is disrupted for
large η values.
4 The telegraph signal
Not only the classical phase diagram, but even the
steady-state density operator solution of the full quan-
tum problem defined by eq. (3) does not describe all the
relevant aspects of our phase transition. In the following
we prove by numerical simulation that the components
of the mixture become robust classical attractors in the
thermodynamic limit. To this end, we need to unravel
the density operator into the time domain, so that we
can extract the dwell timescales, we can show their di-
vergence, and we can determine the relevant exponents.
To this end, we use the quantum trajectories generated
by the Monte-Carlo wavefunction method. In principle,
the ensemble average of many trajectories yields a den-
sity operator evolving in time towards the steady-state
one. However, in the ergodic case (which will be our as-
sumption here), the temporal averaging of the stochas-
tic state vectors along a single long quantum trajectory
yields the same steady-state density operator. There-
fore it makes sense to consider a trajectory as an actual
evolution under continuous measurement with an ideal
photodetector.
Obviously, due to the large difference in the photon
numbers, the components of the mixture correspond to
very distinct output signals. The photon number is con-
tinuously monitored via the photons outcoupled into the
Accepted in Quantum 2019-05-24, click title to verify 5
0
100
200
photon number
γ=0
a)
−2.5
0.0
2.5
field phase
b)
0
5
Mandel-Q
c)
1000 1500 2000 2500 3000 3500 4000
κt
0
100
200
photon number
γ=0.01κ
d)
g/κ
25
35
50
70
85
100
Figure 3: Example trajectories showing the blinking behaviour with different g values (color code is the same throughout). Param-
eters: ∆ = −5κ, η = g/4 (corresponding to the cyan star in fig. 2). The dotted lines in panels a) and d) represent the estimate
g
2
/(2∆
2
) for the photon number, following from eq. (7). Panel b) shows the phase of the field along the green trajectory of panel
a) (g = 50κ), assuming a coherent state. The dotted green line is the field phase expected from eq. (7), that is, the phase of the
complex number
κ −i
∆ + g/(2
√
N)
−1
, substituting the bright-state photon number N = 50 that can be read from panel
a). We see that the coincidence with the simulated field phase in the ON periods is very good. Panel c) displays the Mandel-Q
parameter along the same trajectory, exhibiting a larger nonclassicality for the OFF periods than the ON periods, which we will
discuss in section 6. Throughout this work, the numerics was performed with C++QED: a C++/Python framework for simulating
open quantum dynamics [34, 35, 36].
Accepted in Quantum 2019-05-24, click title to verify 6
ON value a
blink-on rate µ = ( dwell time in OFF period )
−1
blink-off rate λ = ( dwell time in ON period )
−1
filling factor F =
µ
µ+λ
expectation value hXi =
aµ
µ+λ
variance var{X} =
a
2
µλ
(µ+λ)
2
temporal correlation hX(t), X(t
0
)i = e
−(µ+λ) |t−t
0
|
var{X}
characteristic timescale τ = (µ + λ)
−1
Table 1: Characteristics of the telegraph process in the random variable X in the special case when the OFF value is 0.
κ loss channel. The classical distinguishability of large
photocurrent versus dark counts amounts to a projec-
tion of the quantum state into only one of the compo-
nents at a time. This is shown in fig. 3, where the instan-
taneous photon number along quantum trajectories is
plotted for various coupling strengths g. The bright and
dark periods alternate sharply in the form of a telegraph
signal. Detailed analysis of the statistical data shows
that the presented signals are indeed very accurately
(even to the limit of numerical accuracy) described by
a telegraph process.
2
This means that such trajectories
have essentially three parameters: the amplitude of the
bright period and the rates of blink-on and -off, µ and
λ, respectively, since a telegraph process is nothing else
than the composition of two temporal Poisson processes
with exponential waiting-time distribution. Hence, the
waiting time for a blink-on (the inverse of the blink-
on rate µ) equals the dwell time in the dim period, the
same being true for blink-off and the bright period. The
characteristics of the telegraph process are summarized
in table 1.
The trajectories for different coupling constants g in
fig. 3 are generated using the drive strength η/g = 1/4
kept invariant, while fixing the detuning ∆ = −5κ.
Hence all these curves correspond to the single point
denoted with the cyan star in the phase diagram fig. 2
within the bistability domain. The photon number in-
creases with increasing g. The numerics shows that the
bright state has the photon number N ≈ g
2
/(2∆
2
),
the corresponding dotted straight lines fitting nicely on
the noisy numerical record. This N is twice the photon
number at the lower boundary of the bistability range
of the neoclassical theory. Interestingly, this analytical
estimate satisfies the simple classical relation (7) with
2
In particular, what we do is to define a binary signal from
the somewhat noisy trajectories whose model is depicted in fig. 3,
simply by assigning the value 1 to the time instants where the
photon number is higher than half of the temporal average and 0
to the others. On this binary signal X(t) we verify the fulfillment
of the relation var{X} = hXi (1 − hXi), that is characteristic of
the telegraph process (cf. table 1). We find agreement up to 10
−14
precision.
high accuracy, which suggests that the intuitive picture
of near-resonantly driving an equidistant ladder holds
even at such drive strength far above the lower bound-
ary of the bistability domain. The same is suggested by
fig. 3(b), where the field phase is plotted together with
an estimate from the same eq. (7), once more to yield
excellent agreement (cf. the figure caption for details).
Since N
scale
∝ g
2
, the photon number measured in
units of N
scale
proves to be invariant for the different
curves. On the other hand, the dwell times in the attrac-
tor states increase significantly with increasing g. This
reveals that there is a thermodynamic limit in which
the phases become robust, the telegraph signal would
disappear and it would be replaced by a hysteresis-like
behaviour in a genuine first-order phase transition.
Figure 3(d) shows that the addition of a small amount
of atomic decay (γ = 0.01κ) does not lead to a qual-
itative change of the conclusions above. The photon
number of the ON period remains the same as with
γ = 0, however, atomic decay does noticeably decrease
the dwell time in this attractor.
5 Filling factor and the scaling of the
drive strength in the thermodynamic limit
In the bimodal steady-state density operator, the
weights of the components vary through the bistability
domain, in particular, on increasing the drive strength
the ‘dim’ state component vanishes gradually in favour
of the ‘bright’ state. In the time domain, the tele-
graph signal manifests this form of transition through
the variation of the filling factor which is theoretically
F =
µ
µ+λ
. The trajectories in fig. 3 show that the fill-
ing factor varies for the telegraph signals with differ-
ent g. This is confirmed in fig. 4. Panel (a) shows the
monotonous increase of F as a function of η/g for a set
of g values. The curves are shifted with respect to each
other. This dependence is made explicit in panel (b),
where g is varied while keeping η/g at various fixed val-
ues. The filling factor is constant in the range of smaller
Accepted in Quantum 2019-05-24, click title to verify 7
Figure 4: Filling factor of the bright periods (cf. table 1) as a
function of η/g and g/κ without (panels a and b) and with a
small (panels c and d) atomic decay γ = 0.01κ. Color code
is the same in the upper and lower row. Panel a): the dotted
lines indicate linear fits in order to determine η
∗
, the drive
strength leading to half filling. The value of η
∗
is indicated for
the coupling strength g = 70κ.
g values, whereas there is a decrease of F in the range
of larger coupling strengths. For example, the green line
represents a closely constant filling factor around 2/3 up
to g ≈ 50κ, but if g is increased further, the filling factor
drops.
3
It follows then that the scaling of the drive strength
η such that η/g is kept invariant does not preserve
the self-similarity of the telegraph signal. Scaling to
the thermodynamic limit means that the system keeps
a self-similar behaviour, only scaled up in time and
brightness of the ON state. Since the latter two scale
with the single parameter g, there is the parameter η left
at our disposal to ensure self-similarity during the up-
scale, which, in the case of such a simple process as the
telegraph one, cannot mean anything else than keeping
the filling factor constant. As the most obvious choice,
we pick the case when the filling factor is 0.5.
Since the concept of self-similarity is a difficult one in
the present context, let us employ a more pictorial ex-
planation. Upscaling means that we are looking on the
system’s time evolution through a telescope that via the
turn of a single knob (here, increasing g, the single scal-
ing parameter), increases its angle of view both in time
and photon number,
4
but keeps projecting the image
3
The reason why the curves span different ranges in panels (a)
and (b) is that we needed to use different ranges of η for different
g values in order to find the value of η
∗
.
4
The increase in these two dimensions is not necessarily at the
same rate, e.g. here the the photon number increases under the
fairly obvious rule g
2
, while the timescale as g
ν
, where ν is one
of the finite-size scaling exponents that we want to find in this
on the same ocular area. (Hence, it has an increasingly
coarse resolution both in time and photon number.)
This is like the “zoom” functionality on modern camera
objectives. It is very important that we are aiming at a
single-parameter scaling theory, that is, we have to find
a rule for how to change the other parameters of the sys-
tem (here, only η) as a function of the scaling parameter
g during the upscale, i.e. how to “scale η with g”. The
rule that defines such a one-dimensional manifold in the
parameter space as the upscale orbit is: self-similarity.
Self-similarity means that we require the image of the
system’s time evolution on the ocular of the telescope
to remain the same during this procedure. In the case
of such a simple process as the telegraph signal, the fill-
ing factor is the single parameter that determines the
image in such a telescope. Hence, self-similarity means
that we are keeping the filling factor constant, namely
0.5.
Therefore, in order to determine the correct scaling of
η in the thermodynamic limit, we have to find the drive
strength for each g value (denoted η
∗
(g) in the follow-
ing) that leads to half filling. Looking at fig. 4(a), we
observe that an approximate linear interpolation (linear
fit) captures appropriately the behaviour of the F(η)
curves for the different g values in the range of interest.
Hence, we use such an interpolation to determine η
∗
(g).
As an example, in the plot we show for the red curve
(g = 70κ) the corresponding η
∗
. In fig. 5, we plot the
function η
∗
(g), embedded in the bistability domain that
we have calculated similarly to fig. 2. With this numer-
ical result we defined the appropriate finite size scaling.
As it turns out, a linear fit reproduces quite exactly the
behaviour of η
∗
for g & 50κ, which indicates that the
correct scaling of the drive strength in the thermody-
namic limit is η ∼ g
2
.
5
6 Dwell times
6.1 Characteristic timescale of the telegraph sig-
nal
Figure 6 presents the characteristic timescale of the
bistable blinking process. This is defined as the inverse
of the sum of blinking rates µ + λ, which is extracted
from an exponential fit on the numerically calculated
temporal self-correlation of the signal (cf. table 1). We
study.
5
It is conceivable that had we included ∆ in the scaling as
suggested in Ref. [33] (∆/g = const.), then the η/g = const. scal-
ing would have been sufficient to preserve self-similarity. However,
to answer this question rigorously, we would need to repeat the
whole numerical work for that scaling also, that is beyond the
scope of the present study.
Accepted in Quantum 2019-05-24, click title to verify 8
Figure 5: The drive strength η
∗
leading to half filling of the
telegraph signal as a function of g, i.e., the equal mixture of
the two phases in the steady-state density matrix. Solid red
curves represent the phase boundaries of the bistability domain
from the classical theory.
use the unit of κ so that the large difference is mani-
fested in the figure: the characteristic times are orders of
magnitude above the microscopic timescale κ
−1
. On in-
creasing the drive strength η in the presented range, the
average dwell time decreases because of the increase of
the rate µ of blink-on, in conjunction with the increase
of the filling factor. Points represent the values of η
∗
where the telegraph signal has half-filling.
From the point of view of the thermodynamic limit,
the dependence of the characteristic timescale on g is
the most relevant (fig. 6(b)). We show the increasing
timescale over two orders of magnitude of the coupling
constant g for various values of η/g. This g range is
given by computational limitations, nevertheless, it is
enough to demonstrate the power-law scaling of the
increase of the timescale and to determine the expo-
nent. To obey the correct self-similar scaling detailed in
section 5, we have to find the timescale for the drive
strength η
∗
(g) for the different g values. To this end,
we again use a linear fit on the curves of η-dependence,
which captures the behavior quite correctly (cf. dotted
lines in the panel (a), with the η
∗
values depicted with
big dots of the corresponding colour). The timescale
change under finite-size scaling τ (g, η
∗
(g)) is shown in
the thick, solid, red line in fig. 6(b). A linear fit in
the log-log scale leads to the numerical estimate of 2.2
for the finite-size scaling exponent of the characteristic
time. The point g = 100κ was omitted from the fitting
because the dwell times were systematically underesti-
mated due to truncation of the trajectories for the long
but finite simulation time.
Figure 7 is devoted to the numerical analysis of the
rates of blink-on and -off, µ and λ, respectively, which
sheds light on the physical processes leading to switch-
ing between the robust classical attractor states in the
finite-size system. These are calculated from combining
the above-discussed characteristic timescale with the
filling factor. For both rates, the increase of g implies
a reduction, in agreement with the expectedly growing
stability of phases on approaching the thermodynamic
limit.
6.2 Blink-off process
The downward process is initiated by a single photon
loss (detection) event. This is because there is a chance
that the state residing in the ladder |n, −i gets projected
into the ladder |n, +i under such an event [13], since
a |n, ±i =
√
n+
√
n−1
2
|n − 1, ±i +
√
n−
√
n−1
2
|n − 1, ∓i
≈
√
n |n − 1, ±i +
1
4
√
n
|n − 1, ∓i. (9)
Once such a jump occurs, there is a downward cas-
cade of photon escapes, because while our (red-detuned)
drive was closely resonant with the high-lying part of
the |n, −i ladder, resulting in an approximately coher-
ent state in this part, it is off-resonant with the |n, +i
ladder (which would be resonant with the drive blue-
detuned with the same amount), so on this ladder there
is no drive to compensate the photon loss. The passage
downward consists of a quick cascade of many jumps
amounting to an exponential decay.
Although a single ladder-switching quantum jump
initiates the blink-off, however, the likelihood of such
a rare jump vanishes completely in the g, N → ∞ ther-
modynamic limit. The rate of ladder switching scales as
κ/N (from Eq. (9), the total jump rate scales as κN ,
while the probability of a ladder switch to occur within
a jump is ∝ 1/N
2
), and as we saw above, the bright-
state photon number scales as g
2
in the thermodynamic
limit. This would suggest the exponent 2 for the finite-
size scaling of the blink-off timescale, which is verified
in the inset of the right panel of fig. 7 to very good accu-
racy. This downward jump rate λ is largely independent
of η, that only slightly influences the photon number in
the bright phase.
6.3 Blink-on process
On the other hand, the switching from the dim to the
bright phase is induced by the external driving and thus
is sensitive to η, as can be seen in the left panel of fig. 7.
The upward process is suppressed by the off-resonance
of low-lying quasi-energy levels (anharmonic part of the
spectrum). Therefore, it is easy to understand that the
larger the coupling g, the larger the shift of the levels
from resonance and the smaller the blink-on rate. In
Accepted in Quantum 2019-05-24, click title to verify 9
Figure 6: The characteristic timescale of the bistable blinking process given by the inverse of µ+λ. The points in panel a) correspond
to the values of η
∗
for the different g values where half-filling of the telegraph signal is achieved. In panel b), the solid red curve
designates the timescale at η
∗
as a function of g, that is, the timescale in the correct finite-size scaling when the telegraph process
is kept at half filling during the passage to the thermodynamic limit. The dotted line shows a log-log linear fit on the values leaving
out g = 100κ. The exponent resulting from the fit is roughly 2.2. Panels c) and d): timescale with finite atomic decay γ = 0.01κ.
Accepted in Quantum 2019-05-24, click title to verify 10
Figure 7: Waiting time for blink-on (µ
−1
) and blink-off (λ
−1
).
Inset: λ values averaged over η as a function of g, with an
exponential fit with exponent 2.
the dim phase, the state of the bosonic mode is close to
the vacuum, however, there must be a small deviation
from that due to the driving. The state is a non-classical
superposition with positive Mandel-Q parameter. The
Mandel-Q parameter measuring the nonclassicality of
the field state is defined as
Q =
var(a
†
a) − ha
†
ai
ha
†
ai
, (10)
where the averaging can be performed either as a quan-
tum average on the actual state vector of the field to re-
flect the nonclassicality at a given time instant, or also
over time. The Q parameter is zero for a classical field
state (coherent state). In fig. 3(c), Q is calculated as a
time-dependent quantum average, and we observe that
the nonclassicality is stronger in the dim phase than in
the bright one. This is consistent with our picture that
the bright phase consists of an approximately coherent
state on the manifold |n, −i.
The state in the dim period has the property that the
projection of the wavefunction after a photon detection
increases the weight of a high-excitation component. It
can be shown that for pure states of a mode, the positiv-
ity of the Mandel-Q parameter, which is what we have
in the dim phase, is equivalent to the nonclassical situ-
ation that a photon escape from the mode increases its
photon number. While in the dim state there is a neg-
ligible amount of excited photon component generated
by the η driving, triggered by a single quantum jump
(which is very rare on account of the very low dim-state
photon number), it can grow in an exponential runaway
process for subsequent photon detections. So the blink-
on also takes place in the form of a cascade of quantum
jumps.
In fig. 8, we plot the field Q parameter averaged over
time, but only in the dim periods.
6
The overall trend
6
Note that a time average of the Mandel-Q calculated for the
Figure 8: Mandel-Q parameter of the field averaged over time
in the dim periods.
is that the nonclassicality of the field in the dim phase
increases both with increasing η and g, and the depen-
dence flattens out for large g at a value close to 1. Hence,
the dim phase remains nonclassical also in the thermo-
dynamic limit.
Let us try to model the dim state of the field in the
following form:
|Ψi =
p
1 − ε
2
|0i + ε |ϕi, (11)
where ε is a small number and |ϕi is a state orthogonal
to the vacuum state. The photon count rate κha
†
ai =
κ ε
2
N
|ϕi
can be made very small with ε → 0, where
N
ϕ
= hϕ|a
†
a|ϕi is the mean photon number associated
with the component |ϕi superposed on the vacuum. At
the same time, the Mandel-Q parameter of the state
(11) is found to be independent of ε in the limit ε 1:
one finds Q
|Ψi
= Q
|ϕi
+ N
|ϕi
. That is, it depends only
on the properties of the state |ϕi. Note that (i) N
|ϕi
> 1
because the vacuum component is missing from photon
number expansion of |ϕi, and (ii) the Mandel-Q pa-
rameter is always limited below 1 in fig. 8. These two
observations imply that the state |ϕi is a nonclassical
state with negative Mandel-Q factor.
instantaneous pure states of trajectories does not reproduce the
Q of the time-averaged density operator, since the dependence of
Q on the density operator is not linear. The aim of our usage is
to substantiate that the quantum state of the mode in the dim
periods is of the form (11), where a photon escape increases the
photon number of the mode.
Accepted in Quantum 2019-05-24, click title to verify 11
(a) (b) (c)
Figure 9: Photon-escape quantum jumps (blue) and the evolution of instantaneous photon-number expectation values (red) during
blink-on and -off events. Green arrows indicate the photon escapes that trigger the switching events. The time-resolved depictions
of quantum jumps are histograms with bin sizes 0.01/κ for (a) and 0.001/κ for (b,c).
6.4 Cascades of quantum jumps switch between
attractors
A study of the photon-number evolution along a single
trajectory together with sufficiently resolved quantum
jump events confirms the above picture both for the
blink-on and -off events, cf. fig. 9. In panel (a) we see
that a single photon escape (marked with the green ar-
row) triggers a shootup of the photon number from the
dim state, in accordance with the dim state being non-
classical with positive Mandel-Q parameter. However,
in this event the surge is not strong enough to break the
blockade. Panel (b) depicts a successful breakthrough
event where we also see that the buildup of the full
bright-state photon number incurs a probabilistic cas-
cade of quantum jump events. In panel (c) also a single
quantum jump triggers the collapse of the bright state
(green arrow at the sudden drop of photon number),
followed by a normal ringdown of the photon number
with rate κ, involving several further photon escapes.
This proves our claim made in the Introduction that the
telegraph signal observed here differs essentially from
the electron-shelving scheme: though there is a trigger
single-photon escape (quantum jump), the switch be-
tween the dim and bright phases is driven by a cascade
of quantum jumps.
7 The role of atomic decay
The condition γ = 0 is essential in the neoclassical the-
ory, since the derivation of the transcendental eq. (6) re-
lies heavily on the fact that the length of atomic pseudo-
spin is conserved, hσ
x
i
2
+ hσ
y
i
2
+ hσ
z
i
2
= 3/4. Allowing
γ 6= 0 hence leads to qualitatively different behaviour
since this conservation law is broken.
Yet, for very small γ values, the behavior of the full
quantum model does not seem to be qualitatively af-
fected: as we have shown in the course of this study in
passing, the photon-blockade-breakdown effect can be
observed in the case of a finite but small γ.
The effect of an atomic decay in the bright state, i.e.,
on a high-lying coherent state can be assessed similarly
to the above:
σ |n, ±i =
1
√
2
(|n − 1, ±i + |n −1, ∓i) . (12)
This means that in the event of an atomic decay, there
is a 1/2 probability of a ladder switch. Hence, a γ on the
order of κ would wipe out the blinking effect that is the
central theme of this paper, since a blink-on would im-
mediately be followed by a collapse of the bright state.
The blinking effect manifests itself most clearly when
γ = 0, the case that we considered mostly here.
If γ κ, which case we exposed in passing in
figs. 3, 4 and 6, the atomic decay emerges as a compet-
ing timescale for high-enough bright-state photon num-
bers, resulting overall in smaller characteristic times
and filling factors. With respect to the phase transi-
tion, the system is interesting as long as the microscopic
timescale of spontaneous emission is longer and thus is
dominated by the shorter macroscopic timescale τ of
the phase stability.
8 Conclusions
The time evolution of a quantum system undergo-
ing a dissipative, first-order quantum phase transition
has been studied. We considered the photon-blockade-
breakdown phase transition, which takes place in the
driven-dissipative Jaynes-Cummings model with strong
coupling between the two-level system and the har-
monic oscillator. For a certain range of drive strength,
the stationary solution corresponds to a bistability of
classically distinguishable states. By unraveling the sta-
tionary solution into quantum trajectories, we resolved
Accepted in Quantum 2019-05-24, click title to verify 12
the nature of coexistence of phases. We constructed
an appropriate scaling of the system parameters such
that the bistability is manifested by a self-similar tele-
graph signal. The finite-size scaling argument verifies
that the bistability solution develops into a first-order
phase transition in the thermodynamic limit. We calcu-
lated the finite-size scaling exponents numerically. Even
in the thermodynamic limit, the stability of phases orig-
inates from the discrete spectrum of a small quantum
system and the dim phase exhibits nonclassical photon
statistics.
Acknowledgements
On behalf of Project WAL, we thank for the usage of
MTA Cloud (https://cloud.mta.hu/) that significantly
helped us achieving the results published in this paper.
This work was supported by the National Research, De-
velopment and Innovation Office of Hungary (NKFIH)
within the Quantum Technology National Excellence
Program (Project No. 2017-1.2.1-NKP-2017-00001) and
by Grant No. K115624. A. Vukics acknowledges support
from the J´anos Bolyai Research Scholarship of the Hun-
garian Academy of Sciences.
References
[1] Ates C, Olmos B, Garrahan J P and Lesanovsky
I 2012 Phys. Rev. A 85(4) 043620 https://doi.
org/10.1103/PhysRevA.85.043620
[2] Carr C, Ritter R, Wade C G, Adams C S
and Weatherill K J 2013 Phys. Rev. Lett.
111(11) 113901 https://doi.org/10.1103/
PhysRevLett.111.113901
[3] Marcuzzi M, Levi E, Diehl S, Garrahan J P and
Lesanovsky I 2014 Phys. Rev. Lett. 113(21) 210401
https://doi.org/10.1103/PhysRevLett.113.
210401
[4] Malossi N, Valado M M, Scotto S, Huillery P, Pillet
P, Ciampini D, Arimondo E and Morsch O 2014
Phys. Rev. Lett. 113(2) 023006 https://doi.org/
10.1103/PhysRevLett.113.023006
[5] Overbeck V R, Maghrebi M F, Gorshkov A V and
Weimer H 2017 Phys. Rev. A 95(4) 042133 https:
//doi.org/10.1103/PhysRevA.95.042133
[6] Letscher F, Thomas O, Niederpr¨um T, Fleis-
chhauer M and Ott H 2017 Phys. Rev. X 7(2)
021020 https://doi.org/10.1103/PhysRevX.7.
021020
[7] Labouvie R, Santra B, Heun S and Ott H 2016
Phys. Rev. Lett. 116(23) 235302 https://doi.
org/10.1103/PhysRevLett.116.235302
[8] Le Boit´e A, Orso G and Ciuti C 2013 Phys. Rev.
Lett. 110(23) 233601 https://doi.org/10.1103/
PhysRevLett.110.233601
[9] Casteels W and Ciuti C 2017 Phys. Rev. A 95(1)
013812 https://doi.org/10.1103/PhysRevA.
95.013812
[10] Casteels W, Fazio R and Ciuti C 2017 Phys.
Rev. A 95(1) 012128 https://doi.org/10.1103/
PhysRevA.95.012128
[11] Rodriguez S R K, Casteels W, Storme F, Car-
lon Zambon N, Sagnes I, Le Gratiet L, Galopin
E, Lemaˆıtre A, Amo A, Ciuti C and Bloch J 2017
Phys. Rev. Lett. 118(24) 247402 https://doi.
org/10.1103/PhysRevLett.118.247402
[12] Fink T, Schade A, H¨ofling S, Schneider C and
Imamoglu A 2018 Nature Physics 14 365–369
https://doi.org/10.1038/s41567-017-0020-9
[13] Carmichael H J 2015 Phys. Rev. X 5(3) 031028
https://doi.org/10.1103/PhysRevX.5.031028
[14] Dombi, Andr´as, Vukics, Andr´as and Domokos, Pe-
ter 2015 Eur. Phys. J. D 69 60 https://doi.org/
10.1140/epjd/e2015-50861-9
[15] P´alyi A, Struck P R, Rudner M, Flens-
berg K and Burkard G 2012 Phys. Rev.
Lett. 108(20) 206811 https://doi.org/10.1103/
PhysRevLett.108.206811
[16] Fink J M, Dombi A, Vukics A, Wallraff A and
Domokos P 2017 Phys. Rev. X 7(1) 011012 https:
//doi.org/10.1103/PhysRevX.7.011012
[17] Bergquist J C, Hulet R G, Itano W M
and Wineland D J 1986 Phys. Rev. Lett.
57(14) 1699–1702 https://doi.org/10.1103/
PhysRevLett.57.1699
[18] Sachdev S 2011 Quantum Phase Transitions (Cam-
bridge University Press) ISBN 978-0-521-51468-2
[19] Marino J and Diehl S 2016 Phys. Rev. Lett.
116(7) 070407 https://doi.org/10.1103/
PhysRevLett.116.070407
[20] Guti´errez-J´auregui R and Carmichael H J 2018
Phys. Rev. A 98(2) 023804 https://doi.org/10.
1103/PhysRevA.98.023804
[21] Reiter F, Nguyen T L, Home J P and Yelin S F
2018 arXiv preprint arXiv:1807.06026
[22] Nagy D, Szirmai G and Domokos P 2011 Phys.
Rev. A 84(4) 043637 https://doi.org/10.1103/
PhysRevA.84.043637
[23] Brennecke F, Mottl R, Baumann K, Landig R,
Donner T and Esslinger T 2013 Proceedings of the
National Academy of Sciences 110 11763–11767
https://doi.org/10.1073/pnas.1306993110
Accepted in Quantum 2019-05-24, click title to verify 13
[24] Bonifacio R, Gronchi M and Lugiato L A 1978
Phys. Rev. A 18(5) 2266–2279 https://doi.org/
10.1103/PhysRevA.18.2266
[25] Pietik¨ainen I, Danilin S, Kumar K S, Veps¨al¨ainen
A, Golubev D S, Tuorila J and Paraoanu G S 2017
Phys. Rev. B 96(2) 020501 https://doi.org/10.
1103/PhysRevB.96.020501
[26] Pietik¨ainen I, Tuorila J, Golubev D and Paraoanu
G 2019 arXiv preprint arXiv:1901.05655
[27] Hwang M J, Puebla R and Plenio M B 2015 Phys.
Rev. Lett. 115(18) 180404 https://doi.org/10.
1103/PhysRevLett.115.180404
[28] Hwang M J and Plenio M B 2016 Phys. Rev.
Lett. 117(12) 123602 https://doi.org/10.1103/
PhysRevLett.117.123602
[29] Larson J and Irish E K 2017 Journal of Physics A:
Mathematical and Theoretical 50 174002 https:
//doi.org/10.1088/1751-8121/aa65dc
[30] Hwang M J, Rabl P and Plenio M B 2018 Phys.
Rev. A 97(1) 013825 https://doi.org/10.1103/
PhysRevA.97.013825
[31] Alsing P, Guo D S and Carmichael H 1992 Physi-
cal Review A 45 5135 https://doi.org/10.1103/
PhysRevA.45.5135
[32] Alsing P and Carmichael H 1991 Quantum Op-
tics: Journal of the European Optical Society Part
B 3 13 https://doi.org/10.1088/0954-8998/
3/1/003
[33] Guti´errez-J´auregui R and Carmichael H J 2018
Phys. Rev. A 98(2) 023804 https://doi.org/10.
1103/PhysRevA.98.023804
[34] Vukics A and Ritsch H 2007 European Physi-
cal Journal D 44 585–599 https://doi.org/10.
1140/epjd/e2007-00210-x
[35] Vukics A 2012 Computer Physics Communica-
tions 183 1381–1396 https://doi.org/10.1016/
j.cpc.2012.02.004
[36] Sandner R and Vukics A 2014 Computer Physics
Communications 185 2380 – 2382 https://doi.
org/10.1016/j.cpc.2014.04.011
Accepted in Quantum 2019-05-24, click title to verify 14