# Optimized diamond inverted nanocones for enhanced color center to fiber
coupling

Cem Güney Torun  Integrated Quantum Photonics, Institute for Physics,
Humboldt-Universität zu Berlin, Newtonstraße 15, 12489 Berlin, Germany
Philipp-Immanuel Schneider  JCMwave GmbH, Bolivarallee 22, D-14050 Berlin,
Germany  Zuse Institute Berlin (ZIB), Takustraße 7, D-14195 Berlin, Germany
Martin Hammerschmidt  JCMwave GmbH, Bolivarallee 22, D-14050 Berlin, Germany
Zuse Institute Berlin (ZIB), Takustraße 7, D-14195 Berlin, Germany  Sven
Burger  JCMwave GmbH, Bolivarallee 22, D-14050 Berlin, Germany  Zuse Institute
Berlin (ZIB), Takustraße 7, D-14195 Berlin, Germany  Joseph H.D. Munns
Integrated Quantum Photonics, Institute for Physics, Humboldt-Universität zu
Berlin, Newtonstraße 15, 12489 Berlin, Germany  Tim Schröder
<EMAIL_ADDRESS>Integrated Quantum Photonics, Institute for
Physics, Humboldt-Universität zu Berlin, Newtonstraße 15, 12489 Berlin,
Germany  Diamond Nanophotonics, Ferdinand-Braun-Institut, Gustav-Kirchhoff-
Straße 4, 12489 Berlin, Germany

(August 27, 2024)

###### Abstract

Nanostructures can be used for boosting the light outcoupling of color centers
in diamond; however, the fiber coupling performance of these nanostructures is
rarely investigated. Here, we use a finite element method for computing the
emission from color centers in inverted nanocones and the overlap of this
emission with the propagation mode in a single-mode fiber. Using different
figures of merit, the inverted nanocone parameters are optimized to obtain
maximal fiber coupling efficiency, free-space collection efficiency, or rate
enhancement. The optimized inverted nanocone designs show promising results
with 66% fiber coupling or 83% free-space coupling efficiency at the tin-
vacancy center zero-phonon line wavelength of 619 nm. Moreover, when evaluated
for broadband performance, the optimized designs show 55% and 76% for fiber
coupling and free-space efficiencies respectively, for collecting the full
tin-vacancy emission spectrum at room temperature. An analysis of fabrication
insensitivity indicates that these nanostructures are robust against
imperfections. For maximum emission rate into a fiber mode, a design with a
Purcell factor of 2.34 is identified. Finally, possible improvements offered
by a hybrid inverted nanocone, formed by patterning into two different
materials, are investigated, and increases the achievable fiber coupling
efficiency to 71%.

Figure 1:  (a) Qualitative properties of nanostructures relevant for photon
collection seen from analogy with a plane mirror, and nanopillar structures:
Directivity of the emission to the collection (upper) side; divergence
determining emission into a particular solid angle; and mode shape determines
how much of the emission can be collected into the fundamental mode of the
waveguide. (b) Dimensions used to parametrize INC design. (c) Optical
collection schemes for determining a figure of merit. The free space
collection efficiency ($\eta_{\mathrm{FS}}$) is defined as the fraction of
emission that can be collected by the first collection optics. The fiber
coupling efficiency ($\eta_{\mathrm{Fib}}$) is defined as the fraction of
emission coupling to the fiber after a second lens is used to image the
emission onto the fiber facet.

Color centers in diamond are solid-state atom-like systems that are a
promising platforms for quantum information processing [1, 2, 3] and quantum
sensing [4, 5] due to their ability to encode optically addressable qubits.
However, challenges regarding the collection of emission from these centers
remain a limiting factor. The refractive index mismatch between diamond
($n=2.41$) and air ($n=1$) presents a critical angle of $24.5^{\circ}$ for
total internal reflection. This results in only around 4% of the emission to
outcouple from the sample to propagating air modes from an emitter in bulk
diamond[6]. Additionally, the dipole emission pattern of color centers is
distributed into a solid angle of 4$\pi$ and thus, constitutes a mismatch with
the targeted single fiber modes.

In order to overcome these limitations, fabricating optical nanostructures
into diamond is a common approach[7]. A wide range of strategies have been
demonstrated, including applications with microlenses[8, 9], nanopillars[10,
6], bullseye gratings[11, 12], nanocones[13, 14], parabolic reflectors[15],
and recently, inverted nanocones[16, 17]. In order to quantify the efficiency
achieved by these nanostructures, a free-space collection efficiency
($\eta_{\mathrm{FS}}$) figure of merit is generally employed. This value is
determined by the emission into a solid angle that is determined by the
numerical aperture of the collection optics. However, for quantum technology
applications, coupling this emission to a single fiber mode is desirable such
that multiple devices can be connected efficiently, and scalably. Furthermore,
for many quantum information processing algorithms, indistinguishable single
photons require collection into a single-mode. In this case, it is the
mismatch between the mode profile of the emission impinging upon the fiber
facet and the fundamental mode of a single-mode fiber that determines the
achievable fiber coupling efficiency ($\eta_{\mathrm{Fib}}$). While the fiber
coupling properties of ‘in-contact’ systems such as waveguide to tapered
fibers[18, 19] and nanodiamond to tapered fibers[20, 21] have been
investigated, this measure of performance has rarely been considered during
the development of nanostructures that rely upon free-space propagation and
collection into fiber[22, 23].

This work focuses on exploring the fiber-coupling efficiency of an inverted
nanocone (INC) nanostructure, due to promising enhancements of the fiber
coupling efficiency offered in this geometry. This particular nanostructure
design has several appealing advantages compared to other designs. The
contributing factors can be conceptually understood in terms of reflection and
waveguiding, as illustrated in figure 1(a). In this context, reflectivity
captures the degree of directivity of the emission towards the collection
optics ($+z$ direction) instead of radiating into $4\pi$ thus enabling free-
space collection efficiencies exceeding 50%. The complementary effect of
waveguiding can be used to confine emission into a smaller solid angle thereby
reducing the diversion of the emitted light. Furthermore, guiding can also be
exploited in order to manipulate the mode shape, which enables higher fiber
coupling efficiencies through engineering emission profiles similar to
Gaussian beams instead of a dipole profile. As an example of this, a ‘trumpet’
structure[24, 25] in semiconductor-quantum dot systems, which has a similar
form to the INC but larger heights (greater than 10 $\mathrm{\upmu m}$), has
been shown to have a very strong mode-shaping property to obtain such beams
with an $\eta_{\mathrm{FS}}$ of few ten percent without utilizing mirrors
below the sample. For the diamond INC, the reflectivity results from total-
and partial internal reflection at the diamond-air interface at the angled
walls of the inverted nanocone. Waveguiding, similarly to nanopillars, arises
since the INC can act as a mode-shaping device supporting quasi-normal modes
to which the quantum emitter emission can couple.

In this letter the performance of INC structures is investigated by optimizing
the design parameters for three alternative metrics for collection
performance: fiber coupling, free-space collection, and rate enhancement.
These three different metrics offer different designs with 66% fiber coupling
efficiency, 83% free-space collection, and 80% fiber coupling of the emission
scaled to the power emitted by a dipole emitter in bulk.

Finite element method simulations are used to model fiber coupling as it is
well suited for simulating accurately point emitters in complex shaped
geometries. These evaluations are computed using the nanophotonics simulation
package JCMsuite[26]. The fiber coupling simulations consist of computing the
emission that scatters from the INC, propagating this through an imaging
system, and calculating its overlap with the HE11 fiber mode[22]. The
presented results are given with a numerical precision well below 0.01 to
ensure numerical stability during the optimization procedures. Details
regarding the simulations, including convergence and computational times, are
presented in Supplementary 1.

The INC geometry, as illustrated in figure 1(b), is parametrized in terms of
full height ($h$) between the top facet of the inverted nanocone and the
sample surface, the top radius ($r_{\mathrm{t}}$) of the emitting facet of the
INC, and the bottom radius ($r_{\mathrm{b}}$) of the circle connecting the INC
to the sample. The color center is modeled as a point dipole centered on the
INC axis ($z$) and aligned in the $x-z$ plane, at a depth of $d_{z}$ from the
top facet. In subsequent analyses, a lateral displacement, $d_{x}$, of the
dipole within the $x-y$ plane is also considered. An emission wavelength of
619 nm is selected, corresponding to the zero-phonon line (ZPL) transition
wavelength of tin-vacancy (SnV) color centers[27] due to their high emission
into the ZPL ($\sim 0.57$ at 10 K)[28] in combination with their possible
favourable spin coherence properties at temperatures achieved by a standard
helium cryostat ($T^{*}_{2}=540\pm 40\,\mathrm{ns}$ at 2.9 K)[29]. The chosen
fiber parameters are based on a 630HP Thorlabs single-mode fiber, and is
modeled with a $3.5\,\upmu\mathrm{m}$ diameter fiber core with
$n_{\mathrm{core}}=1.4632$ and an infinitely large cladding with
$n_{\mathrm{clad}}=1.4574$ at 619 nm.

In all simulations, an optical imaging system with a numerical aperture (NA)
of 0.9 is used to collect emission from the nanostructure. The fraction of
emission coupled into this solid angle determines $\eta_{\mathrm{FS}}$. For
fiber-coupling efficiency, an additional lens-system is implemented for
further mode matching to the fiber mode. A schematic drawing of the collection
schemes is shown in figure 1(c). A comparative discussion on the efficiency of
a direct coupling scheme without the imaging system is provided in
Supplementary 4.a.

For determining the optimal parameter values for maximal efficiencies, the
Bayesian optimizer of JCMsuite is employed. The method enables a very
efficient global search in a high dimensional parameter space with many local
minima [30]. After performing an optimization (labeled $\Lambda$), the set of
parameters obtained are referred to with an identification label
INC-$\Lambda$.

We first consider the performance of an “ideal” INC that is optimized for
fiber coupling (INC-Fib) by varying over the design parameters in figure 1(b).
For maximal coupling to the INC modes, and thus an optimized
$\eta_{\mathrm{Fib}}$, the dipole orientation is fixed parallel to the sample
surface, aligned along the $x$ axis. The optimization was performed over $h$,
$d_{z}$, and $r_{\mathrm{t}}$ in order to maximize $\eta_{\mathrm{Fib}}$,
where $r_{\mathrm{b}}=1\,\mathrm{nm}$ was set as a constraint. This study
obtained $\eta_{\mathrm{Fib}}=0.66$ with design parameters
$h=635\,\mathrm{nm}$, $r_{\mathrm{t}}=391\,\mathrm{nm}$, and
$d_{z}=507\,\mathrm{nm}$, resulting in a side-wall angle of $58^{\circ}$
between the sample surface and the INC. Further details of the optimization
conditions are provided in Supplementary 2.

As discussed, however, free-space coupling is a more commonplace figure of
merit. For comparison with prior work, the optimization procedure is therefore
repeated to maximize $\eta_{\mathrm{FS}}$ instead, yielding design parameters
for an INC that is optimized for free-space collection (INC-FS). In this case
$\eta_{\mathrm{FS}}=0.83$ with $h=511\,\mathrm{nm}$,
$r_{\mathrm{t}}=848\,\mathrm{nm}$, and $d_{z}=82\,\mathrm{nm}$ resulting in a
side-wall angle of $31^{\circ}$.

While the fiber coupling properties of other diamond nanostructures have not
been investigated, the free-space efficiency, $\eta_{\mathrm{FS}}=0.83$ of the
INC-FS demonstrates significant improvement to simulations showing
$\eta_{\mathrm{FS}}=0.298$ for a microlens[8] (NA:0.9, averaged over 600-800
nm) and $\eta_{\mathrm{FS}}=0.53$ for a nanopillar[6] (NA:0.9, at 637 nm).
Furthermore, it slightly outperforms $\eta_{\mathrm{FS}}=0.752$ of an
optimized bullseye grating[12] (NA:0.95, at 637 nm) and
$\eta_{\mathrm{FS}}=0.77$ of a parabolic reflector[15] (NA:1.3, averaged over
600-800 nm).

From the definitions of the alternative efficiency measures, both
$\eta_{\mathrm{FS}}$ and $\eta_{\mathrm{Fib}}$ necessitate emission into the
solid angle determined by the collection optics. Fiber coupling, however,
places additional requirements on the mode quality of emission. Consequently,
structures optimized for fiber coupling efficiency therefore also enable high
free-space collection efficiency, whereas designs considering only free-space
collection will not necessarily couple well to a fiber mode (see Supplementary
4.b). In the far-field, it is readily seen that INC-Fib produces a near
uniform Gaussian beam profile, figure 2(a), whereas INC-FS exhibits a
multimode profile with significant sidelobes, figure 2(c), due to the larger
top radius.

Since $\eta_{\mathrm{FS}}$ is insensitive to the emitted mode profile, the
convergence to a large top radius of the INC-FS is driven by increasing the
directivity of the emission in $+z$, through increasing the opening angle of
the cone and hence, the fraction of emission fulfilling the total internal
reflection condition (seen in figure 2(d) showing the INC cross-section and
the local field distribution). Conversely, INC-Fib tends towards a larger
$d_{z}$ and reduced $r_{\mathrm{t}}$ which have the effect of increasing the
waveguiding effect of the INC, which together result in a near single-mode
emission profile figure 2(b).

Figure 2:  Far-field emission profiles and field distributions along the
computational domain for a dipole radiating at 619 nm. (a) Far-field emission
profile and (b) field distribution resulting from an INC optimized for fiber
coupling (INC-Fib). (c, d) Corresponding field profiles for an INC optimized
for free-space coupling (INC-FS). The circle in the far-field profiles (a, c)
indicates the solid angle from which emission can be collected, and is used to
determine the free-space coupling efficiency. The outline in the field
distribution sections (b, d) highlights the outer dimensions of the INC and
sample substrate.

Having identified the essential elements for the design of INC structures for
optimized collection efficiencies, the feasibility of fabricating devices in
practice must also be considered. Fabrication processes introduce
uncertainties in the nanostructure shape and dimensions, which in turn impact
the final performance of a device. In the following analysis, a fabrication
precision of $\sim 10\,\mathrm{nm}$ is assumed to be achievable.

To provide an estimate of tolerances during fabrication, a fabrication
insensitivity score $\mathcal{S}(\zeta)$ is obtained from the variation of a
chosen performance metric $\zeta$, as each of the geometrical parameters of a
given INC-$\Lambda$ is individually swept about the target value, over a range
reflecting physically achievable precision, while maintaining other dimensions
constant. Further information regarding how these metrics are obtained are
detailed in Supplementary 3.

Evaluation against this metric yields similarly high values for both
structures considered, with $\mathcal{S}(\eta_{\mathrm{Fib}})=0.92$ for INC-
Fib, and $\mathcal{S}(\eta_{\mathrm{FS}})=0.94$ for INC-FS. This suggests that
the INC structure may offer a repeatable – and therefore scalable – design for
enhancing the achievable collection efficiencies from quantum emitters.

Figure 3:  Comparison of the resulting performance and geometries of INC
optimized for different figures of merit. A full tabulation of results can be
found in table 1. (a) Bar chart representing the performances of the five
optimized INCs. The INC parameters are optimized using the figure of merit
$\zeta$ indicated by the star. Further performance metrics
$\left\langle\eta_{\mathrm{Fib}}\right\rangle_{\lambda}$ (light yellow bar)
and $\mathcal{S}(\zeta)$ (pink bar) are evaluated using $\zeta$. (b)
Geometries of the five INC plotted up to scale with inscribed side-wall
angles.

Depending on the application, the total collectible emission rate rather than
purely the fraction of light coupling into a single-mode fiber may be a more
significant measure of performance. The local density of states is modified by
the INC and determines the absolute fluorescence rates of the color centers
through the Purcell effect. In comparison to bulk diamond, the effective
refractive index of the INC is in general reduced, since the mode is
distributed between diamond and air; hence, unless a cavity effect is achieved
in an INC, its fluorescence emission rate is suppressed with respect to
emitters in bulk. The ratio of emission inside the nanostructure and the bulk
is expressed by the Purcell factor ($F$). The INC-Fib, optimized for fiber
coupling, has a $F=0.39$, whereas the INC-FS achieves a weak cavity effect
with $F=1.32$.

Analogous to brightness enhancement[31] for free space coupling, an
alternative figure of merit may be employed to quantify the fiber-coupled rate
enhancement, i.e. the fraction of fiber coupled light under normalization
against a dipole emitter in bulk diamond. The rate enhancement therefore
provides a means to optimize the total usable count rate. In this case, we
define $\mathcal{R}=F\eta_{\mathrm{Fib}}$ to measure the rate enhancement in
conjunction with fiber coupling.

Optimizing the INC parameters for brightness with fiber coupling (INC-Fib:RE)
yields $\mathcal{R}=0.80$, with $\eta_{\mathrm{Fib}}=0.34$, and $F=2.34$,
which constitutes a factor of $\sim 3$ improvement in rate compared with the
INC-Fib ($\mathcal{R}=0.26$); however, since the rate increase originates from
a cavity enhancement, the performance becomes more sensitive to any variations
in the structure geometry – thus introducing challenges in fabrication – seen
in $\mathcal{S}(\mathcal{R})=0.19$ for INC-Fib:RE, as compared to $\sim 0.9$
for INC-Fib and INC-FS.

A further important consideration is the impact of dipole orientation within
the nanostructure. To obtain a horizontally oriented dipole for SnV color
centers, a (110) terminated surface crystal could be used, in which two of the
four SnV dipole orientations within the lattice
($\left[1\overline{1}\overline{1}\right]$,
$\left[\overline{1}1\overline{1}\right]$) are parallel to the diamond surface.
However, many commercial diamond samples have a (100) terminated surface, for
which all 4 SnV dipole orientations, $\left\langle 111\right\rangle$, are
equivalent, and result in an angle of 35.3∘ between the dipole axis and the
(100) plane.

Optimizing for $\eta_{\mathrm{Fib}}$ with the tilted dipole (TD) orientation
($\left\langle 111\right\rangle$ dipole in $(100)$ crystal) yields the
parameters for INC-Fib:TD, for which $\eta_{\mathrm{Fib}}=0.47$. The
approximately 20% reduction of efficiency with respect to the
$\eta_{\mathrm{Fib}}=0.66$ for INC-Fib arises from the reduced symmetry, which
therefore impacts the far-field mode quality. The optimization attempts to
compensate for this and converges to approximately twice the height of INC-Fib
in order to increase the wave-guiding property of the structure, whilst
retaining similar proportions in terms of the INC opening angle and relative
depth of the emitter. Conversely, if the dipole in INC-Fib:TD is horizontal
instead, then the fiber coupling increases to $\eta_{\mathrm{Fib}}=0.60$, near
the optimized case.

The fabrication insensitivity without changing $d_{x}$ remains reasonably
robust with $\mathcal{S}(\eta_{\mathrm{Fib}})=0.81$, but is more sensitive to
$x$ position due to the off-axis dipole. An analysis of this effect and
further comparative discussions of the optimized INCs are provided in
Supplementary 4.d-e.

Figure 4:  Efficiencies achieved by the optimized inverted nanocones at
different wavelengths and the tin-vacancy color center fluorescence spectrum
measured at room temperature under 520 nm non-resonant excitation.

While many quantum information processing applications require
indistinguishable photons collected from the zero-phonon line, other
applications may make use of the full emission spectrum. We therefore define a
measure of the broadband performance for metric $\zeta$ of a given structure
by averaging over the spectral intensity, $I(\lambda)$, of the emitter:
$\left\langle\zeta\right\rangle_{\lambda}=\int\\!\mathrm{d}\lambda\,\zeta(\lambda)I(\lambda),$
where $I(\lambda)$ is taken from a measurement of an SnV fluorescence spectrum
measured at room temperature, in the region 580-750 nm (see Supplementary 5
for experimental details). Using this measure, INC-Fib achieves a broadband
fiber coupling efficiency
$\left\langle\eta_{\mathrm{Fib}}\right\rangle_{\lambda}=0.55$.

Figure 4 shows $I(\lambda)$, and $\zeta(\lambda)$ for the different INC
designs. In most cases, the performance is only weakly modulated as a function
of the wavelength (see Supplementary 4.f); however, a strong modulation is
observed for INC-Fib:RE. Since the Purcell enhancement necessitates coupling
to quasi-normal modes (QNMs) of the structure, the stronger oscillatory
behavior can be attributed to the presence of QNMs with narrow bandwidths.
Similar modulations are exhibited for the INC-Fib:TD, which has a comparable
geometry to INC-Fib:RE, and may support similar modes, and thus has a stronger
waveguiding effect in contrast to the shorter INCs. Furthermore, the asymmetry
of source within the INC could allow for the possibility of weakly coupling
with more QNMs, and correspondingly has less trivial wavelength dependence of
the final mode-mismatch with the single-mode fiber.

As a final remark, the waveguiding property of an optical structure is not
only a consequence of the geometry but also the material refractive index. A
supplementary approach to engineering the guiding properties is to consider a
hybrid structure incorporating multiple materials. As a simple example of
this, a multi-material INC is considered, in which the top section of the INC
above the emitter is substituted for a higher $n$ material, resulting in
stronger confinement of the light and thereby directing emission towards the
upper half. Here we consider a layer of a conceptual material with $n=3.65$,
deposited onto a bulk diamond sample, and an INC fabricated out of both
materials together. In practice, AlGaAs may achieve this value[32] while GaP
has a similar refractive index as high as $n=3.33$[33] and both materials
could in principle be used for fabrication.

To provide an indicative analysis, an analogous optimization to INC-Fib was
performed with the additional free parameter of the high-$n$ coating
thickness. Optimisation (INC-Fib:Hyb) achieves both an improvement in fiber
coupling to the pure diamond INC-Fib with $\eta_{\mathrm{Fib}}=0.71$, while
also reaching a comparable free-space efficiency to the INC-FS with
$\eta_{\mathrm{FS}}=0.80$. Furthermore, the fabrication insensitivity of fiber
coupling also remains high with $\mathcal{S}(\eta_{\mathrm{Fib}})=0.90$. While
such a hybrid material system has not been fabricated and it is therefore
difficult to predict its feasibility, this design offers an effective and
relatively simple means for further enhancing the collection efficiency from
color centers in diamond compared to pure diamond designs.

In conclusion, the fiber coupling properties of an optical nanostructure, the
INC, is investigated. This detailed numerical investigation focuses
particularly on fiber-coupling as an important – and often neglected – figure
of merit that is requisite for providing indistinguishable single modes for
the generation of indistinguishable photons, as well as a means for coupling
multiple quantum photonic devices in a modular and scalable way. Although
defects in diamond have been considered here, the methods presented are more
generally applicable to other quantum emitters and host materials, which
require free-space collection in their operation.

INC structures providing typically 60% fiber coupling efficiency (70% free-
space coupling efficiency) were found, and indicate promising tolerance to
fabrication errors. A summary of results and comparisons is tabulated in table
1.

INC-$\Lambda$ | $\eta_{\mathrm{fib}}$ | $\eta_{\mathrm{fs}}$ | $\mathcal{R}$ | $\left\langle\zeta\right\rangle_{\lambda}$ | $h$ | $r_{\mathrm{t}}$ | $d_{z}$ | $\mathcal{S}(\zeta)$
---|---|---|---|---|---|---|---|---
|  |  |  | (nm) | (nm) | (nm) |
Fib111Evaluated with $\zeta=\eta_{\mathrm{Fib}}$. | 0.66 | 0.74 | 0.26 | 0.55 | 635 | 391 | 507 | 0.92
FS222Evaluated with $\zeta=\eta_{\mathrm{FS}}$. | 0.22 | 0.83 | 0.29 | 0.76 | 511 | 848 | 82 | 0.94
Fib:RE33footnotemark: 3 | 0.34 | 0.49 | 0.80 | 0.53 | 1 599 | 914 | 1 286 | 0.19
Fib:TD111Evaluated with $\zeta=\eta_{\mathrm{Fib}}$. | 0.47 | 0.63 | 0.25 | 0.44 | 1 375 | 953 | 1 261 | 0.81
Fib:Hyb111Evaluated with $\zeta=\eta_{\mathrm{Fib}}$.444Total INC dimensions, where $h_{\mathrm{diamond}}=855\,\mathrm{nm}$; $h_{\mathrm{coating}}=607\,\mathrm{nm}$. | 0.71 | 0.80 | 0.28 | 0.62 | 1 462 | 624 | 1 295 | 0.90

22footnotetext: Evaluated with $\zeta=\mathcal{R}$.

Table 1:  Performance metrics and corresponding dimensions for INCs optimized
according to different figures of merit (purple text) for a dipole radiating
at 619 nm. After optimization, performance against alternative metrics is
evaluated and given.

These results highlight that while free-space coupling only considers
directivity, fiber coupling additionally requires improved mode-quality of
emission, and is a more practical figure of merit in quantum technology
applications. Optimization of inverted nanocones for fiber coupling generally
finds wide side-wall angles between $\sim 55-70^{\circ}$, with the emitter
near the base of the structure at depths of $\sim 80-90\%$ of the full
structure height. The initial indicative analyses on its tolerance to
fabrication errors show that within reasonable fabricable precision of 10 nm
the performance remains robust. The INC, therefore, presents an appealing
design for the enhancement of photon collection from color centers in diamond,
due to its favorable geometry which simultaneously enables high directivity,
controlled divergence, and single-mode emission, that is provided in a package
promising fabrication simplicity, and is readily integrable into cryogenic
setups.

## SUPPLEMENTARY MATERIAL

See supplementary material for details regarding the simulation steps,
optimizations, fabrication insensitivity analysis, further discussions and
experimental acquisition of the spectrum.

###### Acknowledgements.

The authors would like to thank Tommaso Pregnolato, Gregor Pieplow, and Julian
Bopp for insightful discussions and technical assistance. This research was
supported through the Federal Ministry of Education and Research of Germany
(BMBF, project DiNOQuant, 13N14921) and the European Research Council (ERC
Starting Grant “QUREP”). P.I.S., M.H. and S.B. acknowledge support by BMBF
Forschungscampus MODAL, project number 05M20ZBM. JHDM gratefully acknowledges
support from the Alexander von Humboldt Foundation.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## REFERENCES

## References

  * Yao  _et al._ [2012] N. Y. Yao, L. Jiang, A. V. Gorshkov, P. C. Maurer, G. Giedke, J. I. Cirac, and M. D. Lukin, “Scalable architecture for a room temperature solid-state quantum information processor,” Nature Communications 3, 800 (2012).
  * Childress and Hanson [2013] L. Childress and R. Hanson, “Diamond NV centers for quantum computing and quantum networks,” MRS Bulletin 38, 134–138 (2013).
  * Atatüre  _et al._ [2018] M. Atatüre, D. Englund, N. Vamivakas, S.-Y. Lee, and J. Wrachtrup, “Material platforms for spin-based photonic quantum technologies,” Nature Reviews Materials 3, 38–51 (2018).
  * Degen, Reinhard, and Cappellaro [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, “Quantum sensing,” Reviews of Modern Physics 89, 035002 (2017).
  * Wu  _et al._ [2016] Y. Wu, F. Jelezko, M. B. Plenio, and T. Weil, “Diamond Quantum Devices in Biology,” Angewandte Chemie International Edition 55, 6586–6598 (2016).
  * Hausmann  _et al._ [2010] B. J. Hausmann, M. Khan, Y. Zhang, T. M. Babinec, K. Martinick, M. McCutcheon, P. R. Hemmer, and M. Lončar, “Fabrication of diamond nanowires for quantum information processing applications,” Diamond and Related Materials 19, 621–629 (2010).
  * Schröder  _et al._ [2016] T. Schröder, S. L. Mouradian, J. Zheng, M. E. Trusheim, M. Walsh, E. H. Chen, L. Li, I. Bayn, and D. Englund, “Quantum nanophotonics in diamond (invited),” J. Opt. Soc. Am. B 33, B65–B83 (2016).
  * Hadden  _et al._ [2010] J. P. Hadden, J. P. Harrison, A. C. Stanley-Clarke, L. Marseglia, Y.-L. D. Ho, B. R. Patton, J. L. O’Brien, and J. G. Rarity, “Strongly enhanced photon collection from diamond defect centers under microfabricated integrated solid immersion lenses ,” Applied Physics Letters 97, 241901 (2010).
  * Jamali  _et al._ [2014] M. Jamali, I. Gerhardt, M. Rezai, K. Frenner, H. Fedder, and J. Wrachtrup, “Microscopic diamond solid-immersion-lenses fabricated around single defect centers by focused ion beam milling,” The Review of scientific instruments 85, 123703 (2014).
  * Babinec  _et al._ [2010] T. M. Babinec, B. J. M. Hausmann, M. Khan, Y. Zhang, J. R. Maze, P. R. Hemmer, and M. Loncar, “A diamond nanowire single-photon source,” Nature Nanotechnology 5, 195–199 (2010).
  * Li  _et al._ [2015] L. Li, E. H. Chen, J. Zheng, S. L. Mouradian, F. Dolde, T. Schröder, S. Karaveli, M. L. Markham, D. J. Twitchen, and D. Englund, “Efficient photon collection from a nitrogen vacancy center in a circular bullseye grating,” Nano Letters 15, 1493–1497 (2015).
  * Zheng  _et al._ [2017] J. Zheng, A. C. Liapis, E. H. Chen, C. T. Black, and D. Englund, “Chirped circular dielectric gratings for near-unity collection efficiency from quantum emitters in bulk diamond,” Optics Express 25, 32420 (2017).
  * Momenzadeh  _et al._ [2015] S. A. Momenzadeh, R. J. Stöhr, F. F. de Oliveira, A. Brunner, A. Denisenko, S. Yang, F. Reinhard, and J. Wrachtrup, “Nanoengineered diamond waveguide as a robust bright platform for nanomagnetometry using shallow nitrogen vacancy centers,” Nano Letters 15, 165–169 (2015).
  * Zhu  _et al._ [2020] T. Zhu, Y. Liang, Z. Liu, J. Fu, Y. Wang, G. Shao, D. Zhao, J. Wang, R. Wang, Q. Wei, W. Wang, F. Wen, T. Min, and H. Wang, “Nanocone Structures Enhancing Nitrogen-Vacancy Center Emissions in Diamonds,” Coatings 10, 513 (2020).
  * Wan  _et al._ [2018] N. H. Wan, B. J. Shields, D. Kim, S. Mouradian, B. Lienhard, M. Walsh, H. Bakhru, T. Schröder, and D. Englund, “Efficient Extraction of Light from a Nitrogen-Vacancy Center in a Diamond Parabolic Reflector,” Nano Letters 18, 2787–2793 (2018).
  * Torun and Schröder [2020] C. G. Torun and T. Schröder, “Poster: FE Simulation of a Dipole Emitter Coupled to an Inverted Diamond Nanocone,” (Workshop on Theoretical and Numerical Tools for Nanophotonics, Berlin, 2020).
  * Jeon  _et al._ [2020] S.-W. Jeon, J. Lee, H. Jung, S.-W. Han, Y.-W. Cho, Y.-S. Kim, H.-T. Lim, Y. Kim, M. Niethammer, W. C. Lim, J. Song, S. Onoda, T. Ohshima, R. Reuter, A. Denisenko, J. Wrachtrup, and S.-Y. Lee, “Bright Nitrogen-Vacancy Centers in Diamond Inverted Nanocones,” ACS Photonics 7, 2739–2747 (2020).
  * Patel  _et al._ [2016] R. N. Patel, T. Schröder, N. Wan, L. Li, S. L. Mouradian, E. H. Chen, and D. R. Englund, “Efficient photon coupling from a diamond nitrogen vacancy center by integration with silica fiber,” Light: Science & Applications 5, e16032 (2016).
  * Burek  _et al._ [2017] M. J. Burek, C. Meuwly, R. E. Evans, M. K. Bhaskar, A. Sipahigil, S. Meesala, B. Machielse, D. D. Sukachev, C. T. Nguyen, J. L. Pacheco, E. Bielejec, M. D. Lukin, and M. Lončar, “Fiber-coupled diamond quantum nanophotonic interface,” Physical Review Applied 8, 024026 (2017).
  * Almokhtar  _et al._ [2014] M. Almokhtar, M. Fujiwara, H. Takashima, and S. Takeuchi, “Numerical simulations of nanodiamond nitrogen-vacancy centers coupled with tapered optical fibers as hybrid quantum nanophotonic devices,” Optics Express 22, 20045–20059 (2014).
  * Vorobyov  _et al._ [2016] V. V. Vorobyov, V. V. Soshenko, S. V. Bolshedvorskii, J. Javadzade, N. Lebedev, A. N. Smolyaninov, V. N. Sorokin, and A. V. Akimov, “Coupling of single NV center to adiabatically tapered optical single mode fiber,” The European Physical Journal D 70, 1–7 (2016).
  * Schneider  _et al._ [2018] P.-I. Schneider, N. Srocka, S. Rodt, L. Zschiedrich, S. Reitzenstein, and S. Burger, “Numerical optimization of the extraction efficiency of a quantum-dot based single-photon emitter into a single-mode fiber,” Optics Express 26, 8479–8492 (2018).
  * Rickert  _et al._ [2019] L. Rickert, T. Kupko, S. Rodt, S. Reitzenstein, and T. Heindel, “Optimized designs for telecom-wavelength quantum light sources based on hybrid circular Bragg gratings,” Optics Express 27, 36824–36837 (2019).
  * Munsch  _et al._ [2013] M. Munsch, N. S. Malik, E. Dupuy, A. Delga, J. Bleuse, J.-M. Gérard, J. Claudon, N. Gregersen, and J. Mørk, “Dielectric GaAs antenna ensuring an efficient broadband coupling between an InAs quantum dot and a Gaussian optical beam,” Physical Review Letters 110, 177402 (2013).
  * Stepanov  _et al._ [2015] P. Stepanov, A. Delga, N. Gregersen, E. Peinke, M. Munsch, J. Teissier, J. Mørk, M. Richard, J. Bleuse, J.-M. Gérard, and J. Claudon, “Highly directive and Gaussian far-field emission from “giant” photonic trumpets,” Applied Physics Letters 107, 141106 (2015).
  * Burger  _et al._ [2008] S. Burger, L. Zschiedrich, J. Pomplun, and F. Schmidt, “JCMsuite: An Adaptive FEM Solver for Precise Simulations in Nano-Optics,” Integrated Photonics and Nanophotonics Research and Applications , ITuE4 (2008).
  * Iwasaki  _et al._ [2017] T. Iwasaki, Y. Miyamoto, T. Taniguchi, P. Siyushev, M. H. Metsch, F. Jelezko, and M. Hatano, “Tin-Vacancy Quantum Emitters in Diamond,” Physical Review Letters 119, 253601 (2017).
  * Görlitz  _et al._ [2020] J. Görlitz, D. Herrmann, G. Thiering, P. Fuchs, M. Gandil, T. Iwasaki, T. Taniguchi, M. Kieschnick, J. Meijer, M. Hatano, A. Gali, and C. Becher, “Spectroscopic investigations of negatively charged tin-vacancy centres in diamond,” New Journal of Physics 22, 013048 (2020).
  * Trusheim  _et al._ [2020] M. E. Trusheim, B. Pingault, N. H. Wan, M. Gündoğan, L. de Santis, R. Debroux, D. Gangloff, C. Purser, K. C. Chen, M. Walsh, J. J. Rose, J. N. Becker, B. Lienhard, E. Bersin, I. Paradeisanos, G. Wang, D. Lyzwa, A. R.-P. Montblanch, G. Malladi, H. Bakhru, A. C. Ferrari, I. A. Walmsley, M. Atatüre, and D. Englund, “Transform-Limited Photons From a Coherent Tin-Vacancy Spin in Diamond,” Physical Review Letters 124, 023602 (2020).
  * Schneider  _et al._ [2019] P.-I. Schneider, X. Garcia Santiago, V. Soltwisch, M. Hammerschmidt, S. Burger, and C. Rockstuhl, “Benchmarking Five Global Optimization Approaches for Nano-optical Shape Optimization and Parameter Reconstruction,” ACS Photonics 6, 2726–2733 (2019).
  * Abudayyeh and Rapaport [2017] H. A. Abudayyeh and R. Rapaport, “Quantum emitters coupled to circular nanoantennas for high-brightness quantum light sources,” Quantum Science and Technology 2, 034004 (2017).
  * Aspnes  _et al._ [1986] D. E. Aspnes, S. M. Kelso, R. A. Logan, and R. Bhat, “Optical properties of AlxGa1-xAs,” Journal of Applied Physics 60, 754–767 (1986).
  * Aspnes and Studna [1983] D. E. Aspnes and A. A. Studna, “Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV,” Physical Review B 27, 985–1009 (1983).


# Non-equilibrium Dynamics and Phase Transitions in Potts model and
Interacting Ehrenfest urn model

Chi-Ho<EMAIL_ADDRESS>and Pik-Yin<EMAIL_ADDRESS>aDept. of Physics, National Changhua University of Education, Changhua 500,
Taiwan, R.O.C.  bDept. of Physics and Center for Complex Systems, National
Central University, Chung-Li District, Taoyuan City 320, Taiwan, R.O.C.
cPhysics Division, National Center for Theoretical Sciences, Taipei 10617,
Taiwan, R.O.C.

###### Abstract

We show that the recently proposed interacting Ehrenfest $M$-urn model at
equilibrium[Phys. Rev. E 101 (2020) 012123] can be exactly mapped to a mean-
field $M$-state Potts model. By exploiting this correspondence, we show that
the $M$-state Potts model with $M\geq 3$, with transition rates motivated by
the non-equilibrium urn model, can exhibit rich non-equilibrium spin dynamics
such as non-equilibrium steady states and non-equilibrium periodic states.
Monte Carlo simulations of the 3-state Potts model are performed to
demonstrate explicitly the first-order transitions for the equilibrium and
non-equilibrium steady states, as well as the far-from-equilibrium periodic
states.

###### keywords:

Potts model, Ehrenfest urn model, non-equilibrium phase transition, spin
dynamics

††journal: Journal of LaTeX Templates

##  1 Introduction

The $M$-state Potts model[1] is a generalization of the classic Ising model[2]
to $M$ spin components. It is known that the Potts model is related to several
paradigm systems in statistical physics [3], in particular for phase
transitions in equilibrium statistical mechanics and its critical behavior has
also been shown to be richer than that of the Ising model. For example, it was
shown[4] that the problem of the bond percolation can be formulated in terms
of the $M=1$ Potts model which was extended further to the site percolation
problem[5]. The Potts model also finds its applications in a variety of
systems related to graphs or networks, such as the chromatic number problem of
graph coloring which is closely related to the ground state of the anti-
ferromagnetic Potts model[3, 6]. Other applications include the graph
partition problem[7, 8], the network-community detection[9], and even the
cell-differentiation model in developmental biology[10]. Although the
equilibrium Potts model and its phase transition behavior have been rather
thoroughly investigated, its non-equilibrium dynamics are much less studied as
compared to the non-equilibrium Ising spin dynamics.

On the other hand, the less well-known Ehrenfest urn model[11, 12, 13],
proposed more than a decade before the Ising model aimed to explain or
illustrate the Second Law of thermodynamics, was one of the earliest models
for non-equilibrium statistical physics. The original Ehrenfest two-urn
model[11] is a simple model with no interaction that is tractable and can
illustrate the conceptual foundation of statistical mechanics for the
relaxation towards equilibrium. Recently, the two-urn Ehrenfest model was
extended to include particle interactions inside an urn[14] in which particles
can interact with all other particles inside the same urn, but particles
belonging to different urns do not interact. In addition, a jumping rate
(asymmetric in general) from one urn to another can also be introduced, which
is independent of the particle interaction. The system was shown to exhibit
interesting phase transitions and the Poincaré cycle together with the
relaxation times can be calculated[14]. The interacting Ehrenfest two-urn
model was subsequently extended to $M$ urns ($M\geq 2$) for the equilibrium
case in which detailed balance can be achieved. The equilibrium phase behavior
is rich and can be investigated in detail[15]. Similar to the two-urn
case[14], $N$ particles are distributed into the $M$ urns. Pairwise all-to-all
interaction is introduced only for particles in the same urn and particles in
different urns do not interact. As the inter-particle interaction strength is
varied, phases of different levels of non-uniformity emerge and their
stabilities are calculated analytically. In particular, the coexistence of
locally stable uniform and non-uniform phases connected by first-order
transition occurs[15]. The phase transition threshold and energy barrier can
be derived exactly together with the phase diagram obtained analytically.
Analytic and exact results are derived for the condition for the emergence of
coexisting uniform or non-uniform phases and the associated first-order phase
transition and energy barrier.

Very recently, the $M$-urn model with intra-urn interactions has been extended
to the non-equilibrium scenario by arranging the urns in a ring and
introducing asymmetric clockwise and counter-clockwise jumping rates[16]. It
was demonstrated that the system can exhibit two distinct non-equilibrium
steady states (NESS) of uniform and non-uniform particle distributions with
the associated non-equilibrium thermodynamic laws revealed. In addition, a
first-order non-equilibrium phase transition occurs[16] between these two
NESSs as the inter-particle attraction varies. The phase boundaries, the NESS
particle distributions near the NESSs and the associated particle fluxes,
average urn population fractions, and the relaxational dynamics to the NESSs
were obtained analytically and verified numerically[16].

As for spin systems, non-equilibrium spin dynamics usually result from the
relaxation of a quenched external field or external time-depending driving.
Periodic oscillations in magnetization can be easily obtained with an external
time-dependent oscillatory drive, but autonomous oscillations in interacting
spin systems are hardly reported.

In this paper, we show and demonstrate that the above seemingly independent
Ehrenfest $M$-urn model and $M$-state Potts model are in fact closely related
and can be mapped to each other. Furthermore, the non-equilibrium dynamics of
these systems reveal non-trivial NESSs and non-equilibrium period states
(NEPS), which may serve as a paradigm statistical physics system to
investigate systems with different levels away from equilibrium. Sec. 2 gives
a summary of the $M$-urn model and its phase transition behavior at
equilibrium, together with the non-equilibrium model of the $M$ urns arranged
on a ring. The equivalence of the $M$-urn model with intra-urn interaction at
equilibrium and the mean-field $M$-state Potts model is established in Sec. 3.
The Potts spin dynamics that correspond to the non-equilibrium $M$-urn model
on a ring are derived in Sec. 4. These theoretical results are verified by the
Monte Carlo simulations of the 3-state Potts model in Sec. 5. Sec. 6 gives the
concluding remarks and future outlooks.

##  2 Summary of the equilibrium and non-equilibrium Ehrenfest urn models with
interactions

In this section, we give a summary of the equilibrium and non-equilibrium
Ehrenfest urn models and their associated properties. Details of these models
can be found in Refs. [15] and [16]. For the $M$-urn model, a total of $N$
particles are placed inside $M$ urns and the number of particles in the
$\alpha^{th}$ urn is denoted by $n_{\alpha}$, with $\alpha=0,1,\cdots,M-1$.
Particles from different urns do not interact whereas any pair of particles in
the same urn interact with an energy specified by the urn. Since the total
particle number $N$ is fixed, the constraint that
$n_{0}+n_{2}+\cdots+n_{M-1}=N$ is always satisfied. The energy or Hamiltonian
of the interacting particles in the urns is given by

$\beta{\cal
H}=\frac{1}{2N}\sum_{\alpha=1}^{M}g_{\alpha}n_{\alpha}(n_{\alpha}-1),$ (1)

where $\beta\equiv 1/(k_{B}T)$ is the inverse temperature and $g_{\alpha}$ is
the pair-wise interaction (in the unit of $k_{B}T$) of the particles inside
the $\alpha^{th}$ urn. It is found that for homogeneous coupling with
$g_{\alpha}=g$, an equilibrium first-order phase transition occurs between the
uniform and the first non-uniform states as $g$ varies and the first-order
transition occurs at $g=g_{t}$ which is given by[15]

$\displaystyle g_{t}$ $\displaystyle=$
$\displaystyle-\frac{2(M-1)}{M-2}\ln(M-1).$ (2)

To quantify how non-uniform the state is, one can define

$\Psi=\sqrt{\frac{1}{M(M-1)}\sum_{\alpha\neq\alpha^{\prime}}(x_{\alpha}-x_{\alpha^{\prime}})^{2}}$
(3)

as the non-uniformity of the state, where
$x_{\alpha}\equiv\frac{n_{\alpha}}{N}$ is the population fraction in the
$\alpha^{th}$ urn. $\Psi$ can also serve as an order parameter for the phase
transition: $\Psi=0$ for the uniform (disordered) state and $\Psi>0$ for the
non-uniform (order) state.

For the non-equilibrium urn model, particles undergo transitions among the
urns, in general one can consider $M$ urns are placed on a network with
possible transitions between the urns represented by the edges. The transition
probability of a particle in the $\alpha^{th}$ urn jumps to the
$\alpha^{\prime th}$ urn is denoted by $T_{\alpha\to\alpha^{\prime}}$. As in
Ref. [16], here we focus on urns placed on a one-dimensional ring for
convenience to investigate the cyclic particle fluxes. In addition, a direct
jumping rate is further introduced such that the probability of anticlockwise
(clockwise) direction is $p$ ($q$). The condition $p+q=1$ is imposed which
changes only the time scale. For the 3-urn system, it has been shown[16] that
two distinct NESSs of uniform and non-uniform particle distributions,
separated by a first-order non-equilibrium phase transition, occur as the
interparticle attraction varies, with the presence of a coexistence regime for
smaller values of $|p-q|$.

NESS corresponds to the situation that an urn receives several particles per
unit of time from the upstream neighboring urn and sends out an equal number
of particles per unit of time to the downstream neighboring urn such that the
number of particles in this urn remains basically unchanged (subject to
stochastic fluctuations). If the inter-particle attraction is not so strong,
the particle numbers in each urn are the same on average resulting in a
uniform NESS. On the other hand, when the inter-particle attraction becomes
strong enough, an urn can have more particles and maintain the population on
average because of stronger attraction resulting in a non-uniform NESS.
Remarkably, under some parameter regimes in which the driving is strong and
the particle attraction is also significant, the fast flow of particles from
upstream will be accumulated on average in an urn for some duration before
their particles can be released to the downstream urn to deplete the
population below average. This can lead to a stable periodic oscillation of
the population in an urn and the particle flux in the ring, resulting in the
NEPS[17].

##  3 Equivalence of the Ehrenfest Urns Hamiltonian to the Potts model

The $M$-urn model at equilibrium can be mapped to the mean-field $M$-state
Potts model with the $i^{th}$ spin (particle) denoted by $\sigma_{i}$ which
can take values $\sigma_{i}=0,1,\cdots,M-1$. We shall use indices
$i,j,k,\cdots$ for the spin/particle label and the Greek alphabet
$\alpha=0,1,\cdots,M-1$ for the urn or Potts state label. Note that the
population fraction in the $\alpha^{th}$ urn is given by

$x_{\alpha}\equiv\frac{n_{\alpha}}{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{\sigma_{i}\alpha},$
(4)

which can be viewed as the mean “magnetization” of the Potts spin along the
$\alpha$-direction (see Fig. 1a). The interacting urn Hamiltonian in Eq. (1)
can be expressed as an interacting Potts spin Hamiltonian using Eq. (4):

$\displaystyle\beta{\cal H}$ $\displaystyle=$
$\displaystyle\frac{1}{2N}\sum_{\alpha=0}^{M-1}g_{\alpha}n_{\alpha}^{2}-\frac{1}{2N}\sum_{\alpha=0}^{M-1}g_{\alpha}n_{\alpha}=\frac{N}{2}\sum_{\alpha=0}^{M-1}g_{\alpha}x_{\alpha}^{2}-\frac{1}{2}g_{\alpha}x_{\alpha}$
(5) $\displaystyle=$
$\displaystyle\frac{1}{2N}\sum_{i,j}^{N}\sum_{\alpha=0}^{M-1}g_{\alpha}\delta_{\sigma_{i}\alpha}\delta_{\sigma_{j}\alpha}-\frac{1}{2N}\sum_{i,j}^{N}\sum_{\alpha=0}^{M-1}g_{\alpha}\delta_{\sigma_{i}\alpha}$
$\displaystyle=$
$\displaystyle\frac{1}{2N}\sum_{i,j}^{N}g_{\sigma_{i}}\sum_{\alpha=0}^{M-1}\delta_{\sigma_{i}\alpha}\delta_{\sigma_{j}\alpha}-\frac{1}{2N}\sum_{i,j}^{N}g_{\sigma_{i}}\sum_{\alpha=0}^{M-1}\delta_{\sigma_{i}\alpha}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i<j}^{N}g_{\sigma_{i}}\delta_{\sigma_{i}\sigma_{j}},$

which is a Potts model with heterogeneous couplings that depend on the spin
directions. In general, the energy change due to the flipping of the $k^{th}$
spin is given by

$\beta\Delta{\cal
H}(\sigma_{k}\to\sigma_{k}^{\prime})=g_{\sigma_{k}^{\prime}}x_{\sigma_{k}^{\prime}}-g_{\sigma_{k}}x_{\sigma_{k}}.$
(6)

In particular, for homogeneous coupling with $g_{\alpha}=g$ for all $\alpha$,
Eq. (5) is the mean-field $M$-state Potts model[3, 18]. To describe the phase
transition in the Pott model, one can define the complex order parameter[3] as

$\phi=\sum_{\alpha=1}^{M}x_{\alpha}e^{\frac{2\pi i\alpha}{M}}\equiv\Phi
e^{i\Theta}.$ (7)

A real-value order parameter can also be defined as[3, 18]

$s\equiv\frac{Mx_{\rm max}-1}{M-1}$ (8)

where $x_{\rm max}$ is the maximal of {$\langle
x_{\alpha}\rangle\\}_{\alpha=0,1\cdots,M-1}$.

From the above mapping, the $M$ urns on a ring can be viewed as interacting
Potts spin system with the spin state related by $Z_{M}$ discrete symmetry as
depicted schematically in Fig. 1. The magnetization along the $\alpha^{th}$
direction can be represented by a vector pointing in that direction with
magnitude $x_{\alpha}$ as shown in Fig. 1a. The disordered state refers to the
situation that the magnetization vectors in all $M$ directions are of equal
magnitude which corresponds to the uniform state with equal particle fractions
in each urn in the $M$-urn model. The transition of a particle to a
neighboring urn can be represented by a “tick of the clock” in Fig. 1a in the
spin direction of a Potts spin.

Figure 1: Schematic picture showing the equivalence of the interacting
Ehrenfest $M$-urn model on a ring and the $M$-state Potts spin with discrete
$Z_{M}$ symmetry. (a) $M$ urns (denoted by the open circles) are placed on a
ring. The $M$ open circles also represent the $M$ discrete Potts spin
directions on the unit circle. The non-equilibrium dynamics are introduced by
the direct counter-clockwise and clockwise jumping rates $p$ and $q$
respectively. (b) Schematic pictures showing the average magnetization in each
direction of the 3-state Potts model under uniform/disordered state (left),
non-uniform/ordered state (middle), and non-uniform non-equilibrium steady
state (right).

Throughout this paper, we shall focus on the case of homogeneous ferromagnetic
coupling case with $g_{\alpha}=g<0$ hereafter. The first-order transition
given in Eq. (2) corresponds to the order(first non-uniform
state)-disordered(uniform state) transition in the mean-field $M$-state Potts
model which agrees with the results in Ref. [18]. We note that the non-
uniformity of the first non-uniform state is related to the order parameter of
the Potts model $s$ in Eq. (8) (which characterizes the order-disordered
transition)[3] by $\Psi^{(1)}=\sqrt{\frac{2}{M}}s$.

The various uniform or non-uniform equilibrium, non-equilibrium steady states
discovered in the Ehrenfest urn models in Refs. [15] and [16] can be
conveniently represented graphically in terms of the Potts spin directions as
depicted in Fig. 1b for the case of $M=3$. For the equilibrium case, the
uniform state in the urn model (or the disordered state in the Potts model)
can be represented by the $M$ symmetric magnetization vectors separated by
$360\deg/M$ as depicted by the left panel in Fig. 1b. The non-uniform state
(the order state in Potts model) is signified by one longer magnetization and
$M-1$ shorter ones of equal lengths due to the spontaneous broken $Z_{M}$
symmetry as depicted by the middle panel in Fig. 1b. For the NESSs in the urn
model (and the corresponding Potts model), the spin magnetization picture is
the same as the uniform equilibrium case. On the other hand, for the non-
uniform NESS with $p>q$, the symmetry of the $M-1$ shorter directions is
further broken and has different lengths(see the right panel in Fig. 1b).

##  4  Non-equilibrium model of $M$ urns on a ring and the Potts spin Dynamics

For the urns/Potts model at equilibrium, detailed balance is obeyed and one
possible (microscopic) transition probability for a spin from the state
$\alpha$ to state $\alpha^{\prime}$ is given by the Glauber dynamics [16]:

$T_{\alpha\to\alpha^{\prime}}=\frac{1}{1+e^{\beta\Delta{\cal
H}}}=\frac{1}{1+e^{\beta g(x_{\alpha^{\prime}}-x_{\alpha})}},$ (9)

and other spin updating dynamics is possible as long as the detailed balance
is satisfied. On the other hand, for non-equilibrium spin models, one needs to
specify the dynamical or updating rules for the spins. For example, there are
different types of kinetic Ising models[19] using Glauber or Kawasaki dynamics
in combination with sequential or parallel updating for the spins, which give
different non-equilibrium properties.

The non-equilibrium dynamics in the urn model on a ring described in Sec. 2
violates the detailed balance condition but still can be phrased in terms of
the Potts spins and Hamiltonian in Eq. (5). In particular, the system of $M$
urns arranged in a ring with intra-urn interactions has recently been studied
showing interesting non-equilibrium thermodynamics and phase transition
behavior. Here we will map the non-equilibrium $M$-urn model on a ring in Ref.
[16] to a Potts model with appropriate spin transition dynamics. In this case,
the flow of particles is restricted to neighboring urns on the ring and thus
the only possible transitions are $T_{\alpha\to\alpha\pm 1}$. In terms of the
Potts spins, the net transition rate for spin $i$ to make the transition from
state $\alpha$ to $\alpha+1$, as depicted schematically in Fig. 1a, is given
by

$K^{(i)}_{\alpha\to\alpha+1}=p\delta_{\sigma_{i}\alpha}T_{\alpha\to\alpha+1}-q\delta_{\sigma_{i}\alpha+1}T_{\alpha+1\to\alpha}.$
(10)

Since the urn labels correspond to the Potts spin states, the transition rates
per particle for urns in a 1D ring in Ref. [16] can be expressed in terms of
the Potts spins as

$\displaystyle W_{\sigma_{k}\to\sigma_{k}+n}$ $\displaystyle=$ $\displaystyle
0\quad\hbox{if }n\neq\pm 1,$ $\displaystyle W_{\sigma_{k}\to\sigma_{k}+1}$
$\displaystyle=$
$\displaystyle\frac{p}{N}T_{\sigma_{k}\to\sigma_{k}+1}\sum_{i=1}^{N}\delta_{\sigma_{i}\sigma_{k}}=px_{\sigma_{k}}T_{\sigma_{k}\to\sigma_{k}+1}$
(11) $\displaystyle W_{\sigma_{k}\to\sigma_{k}-1}$ $\displaystyle=$
$\displaystyle\frac{q}{N}T_{\sigma_{k}\to\sigma_{k}-1}\sum_{i=1}^{N}\delta_{\sigma_{i}\sigma_{k}}=qx_{\sigma_{k}}T_{\sigma_{k}\to\sigma_{k}-1},$

where it is understood that the value of the spin state is always under
modulus $M$ to respect the cyclic representation of the Potts states (or the
periodic boundary condition of the urns in a ring). For the special case of
$p=q$, it reduces to the equilibrium case, and the detailed balance is obeyed.
Thus the non-equilibrium urn model on a ring introduced in Ref. [16] can be
viewed as a mean-field $M$-state Potts with special transition rates
restricted to neighboring spin states as given by Eq. (11) that breaks
detailed balance in general and gives rise to a variety of non-equilibrium
states. The transition between the neighboring Potts spin states for spin $i$
can be viewed as a stochastic ‘ticking clock’ in which the spin direction
$\sigma_{i}$ undergoes clockwise or counter-clockwise discrete rotation of an
angular step of $\frac{2\pi}{M}$ (see Fig. 1a). The net transition rate of the
mean magnetization from $\alpha\to\alpha+1$ of the system is then

$\displaystyle K_{\alpha\to\alpha+1}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}K^{(i)}_{\alpha\to\alpha+1}=px_{\alpha}T_{\alpha\to\alpha+1}-qx_{\alpha+1}T_{\alpha+1\to\alpha}$
(12) $\displaystyle=$ $\displaystyle
W_{\alpha\to\alpha+1}-W_{\alpha+1\to\alpha}=\frac{px_{\alpha}e^{gx_{\alpha}}-qx_{\alpha+1}e^{gx_{\alpha}+1}}{e^{gx_{\alpha}}+e^{gx_{\alpha}+1}},$
(13)

which is also the mean particle flux from the $\alpha$ to $\alpha+1$ urns in
the urn model[16].

The NESSs in the urns model correspond to the steady flipping of the Potts
spin state $\alpha\to\alpha+1$ in a cyclic manner with the mean magnetizations
remaining constant in time and that is why the uniform NESS and the uniform
equilibrium state have the same magnetization diagram shown in the left panel
of Fig. 1b. It is clear from the above discussion that the mean magnetizations
are the same for the uniform equilibrium state and uniform NESS, thus to
distinguish the non-equilibrium and equilibrium states in the Potts variables,
one needs to measure the net mean transition rates of the spin directions from
$\alpha\to\alpha+1$ by monitoring the microscopic transitions in the MC
simulations. One can measure the number of transitions from
$\sigma_{\alpha}\to\sigma_{\alpha+1}$ in a given period of time and average
over all the spins as given by Eq. (12). One expects $\langle
K_{\alpha\to\alpha+1}\rangle\simeq 0$ for equilibrium states, whereas
$|\langle K_{\alpha\to\alpha+1}\rangle|>0$ for non-equilibrium states. In
addition, one can also measure the time-dependent transition rates
$K_{\alpha\to\alpha+1}(t)$ from the time-series data of $x_{\alpha}(t)$ using
Eq.(13).

##  5 Monte Carlo Simulation results

To explicitly verify the theoretical results in previous sections, one can
carry out Monte Carlo simulations for the $M$-state Potts system. In such a
simulation, a total of $N$ ($N$ is an integer multiple of $M$) Potts spins are
placed in the system. The simulation algorithm starts by choosing a spin $i$
(whose state is $\sigma_{i}$) at random with uniform probability and the final
state of the spin is also chosen randomly from the allowed final states with
the transition rate given by Eq. (11). In particular, for the urn model on a
ring, the two allowed final states are chosen with probability $p$ for
$\sigma_{i}+1$ and $q=1-p$ for $\sigma_{i}-1$ respectively. Then the trial
spin transition is made with probability given by Eq. (9). If
$p=q=\frac{1}{2}$, then the transition rules satisfy the detailed balance
condition and the system is at equilibrium. In general if $p>q$, there will be
a counter-clockwise flux and a non-equilibrium state can be achieved. After
some sufficiently long transient time, the mean magnetization (or population
fraction of an urn) $x_{\alpha}(t)$, is measured using Eq. (4), and the
complex order parameter ($\phi(t)$) in Eq. (7) of the Potts spin system is
monitored.

On the other hand, due to the equivalence of the mean-field Potts model and
the urn model, one can also simulate directly the $M$-urn system on a network
as was done in Ref. [15] for the equilibrium case. For simulations of the
particles in the urns, the algorithm starts by choosing a particle at random
with uniform probability, say the chosen particle is in the $\alpha^{th}$ urn,
then its jumping probability to another urn $\alpha^{\prime}$ is given by the
weighted transition $p_{\alpha\alpha^{\prime}}T_{\alpha\to\alpha^{\prime}}$,
where $p_{\alpha\alpha^{\prime}}$ is a normalized weighing factor and
$T_{\alpha\to\alpha^{\prime}}$ is given by Eq. (9). For the equilibrium case,
there is no bias in the weighting factor, i.e. unweighted, and one can simply
put $p_{\alpha\alpha^{\prime}}=1$. On the other hand, for the non-equilibrium
model of $M$ urns on a ring, the allowed jumping probabilities are
$pT_{\alpha\to\alpha+1}$ and $(1-p)T_{\alpha\to\alpha-1}$ respectively for the
counter-clockwise and clockwise jumps respectively.

In this paper, we carry out MC simulations of the 3-state Potts model with the
clockwise jumping rate $p$ and transition rates specified by Eq. (11). Time
evolution of the detailed dynamics of the order parameter is monitored and the
probability distributions of the amplitude and phase of the Potts order
parameter are measured. Time is in the unit of Monte Carlo Steps per particle
(MCS/N). One MCS/N means that on average every spin has attempted a flipping.
Most simulations are performed with $N=1500$ spins, which is already large
enough to demonstrate the predicted phases and their transitions (see Figs. 2
and 3).

###  5.1 Equilibrium first-order transition

We first simulate the 3-state and 4-state Potts model at equilibrium (i.e.
with equal jumping rates to other spin states) with the mean-field Hamiltonian
in Eq. (5). Fig. 2a shows the average amplitude of the order parameter given
by Eq. (7) as a function of the coupling $g$ in the unit of $k_{B}T$.
$\langle\Phi\rangle$ shows an abrupt jump as the coupling increases to a
threshold value that agrees well with the theoretical first-order transition
point in Eq. (2), for both the 3-state and 4-state Potts model. The first-
order transition nature is revealed clearly in the two coexisting peaks in the
distribution function of $\Phi$ as shown in Fig. 2b near the transition point
for the 3-state Potts model. Results for smaller ($N=300$) and larger
($N=3000$) system sizes are also shown to confirm that the system size
employed in the simulations is sufficiently large to observe the
characteristic behavior of coexisting uniform and non-uniform equilibrium
phases, and that the barrier between these coexisting phases increases with
$N$. The distributions for two of the magnetizations near the transition
point, which coincide with each other, display three peaks that correspond to
the coexistence of the disordered states($x_{\alpha}\simeq\frac{1}{3}$) and
the ordered state, as shown in Fig. 2c. The state of the system can be
revealed in detail by examining the evolution of the order parameter $\phi(t)$
on its complex plane. Fig. 2d displays the time evolution of $\phi(t)$ for a
duration of 5000 MCS/N for $g=-2.7$ (the disordered state, black) and of
$g=-3.2$ (the three symmetry-related disordered states, red). Each disordered
state is prepared separately with a different initial condition, and the
hopping between these three symmetry-related disordered states will take a
very long time due to the high free energy barrier separating them.

Figure 2: Monte Carlo simulation results of the ferromagnetic Potts model at
equilibrium for $N=1500$. (a) The average order parameter amplitude, $\Phi$,
of the 3-state and 4-state Potts models plotted against the coupling $g$ (in
the unit of $k_{B}T$). The vertical dashed lines are the theoretical first-
order transition values given by Eq. (2). The result with $N=3000$ for the
3-state model is also shown. (b) The measured distribution function of $\Phi$,
for the 3-state Potts model near the first-order transition at $g=-2.77$. The
co-existing phases are clearly revealed from the two peaks. Results for
$N=300$ and $N=3000$ are also displayed to show that the barrier between the
coexisting phases increases with $N$. $10^{6}$ MCS/N are used to obtain good
statistics. (c) The measured distributions of the two directions of the
magnetization, $P(x_{\alpha})$ for the case in (b). (d) Time evolution of the
complex order parameter $\phi$ in its complex plane for a duration of 5000
MCS/N for $g=-2.7$ (black) and the three symmetry-related states of $g=-3.2$
(red) which are prepared separately with different initial conditions. The
dashed lines mark the three Potts spin directions of 0 and
$\pm\frac{2\pi}{3}$.

###  5.2 NESSs and spin state transition rates

Here we consider the 3-state Potts model with asymmetric jumping rates for
$\alpha\to\alpha+1$ (with jumping rate $p$) and $\alpha\to\alpha-1$ (with
jumping rate $q=1-p\neq p$), as depicted in Fig. 1a. Fig. 3a shows the average
amplitude of the order parameter given by Eq. (7) as a function of $g$ (in the
unit of $k_{B}T$). For smaller values of $p$, $\langle\Phi\rangle$ shows an
abrupt jump as the coupling increases at some threshold value, signaling a
first-order transition similar to the equilibrium ($p=\frac{1}{2}$) case. But
for larger values of $p$ (see the $p=0.9$ curve), the sharp jump is replaced
by a gradual increase in $\langle\Phi\rangle$ for $-g$ between 3 and $\sim
3.5$, then followed by a steeper increase, suggesting a more complex non-
equilibrium dynamics. The first-order non-equilibrium phase transition is
revealed from the two peaks of $P(\Phi)$ at $g=-3.07$ for $p=0.7$, as shown in
Fig. 3b. Results for smaller ($N=300$) and larger ($N=3000$) system sizes are
also shown to confirm that the system size employed in the simulations is
sufficiently large to observe the characteristic behavior of coexistence of
the two NESS phases, and that the barrier between them increases with $N$. The
relative heights of the two peaks vary with the system sizes arising from the
systematic shift of the phase transition thresholds resulting from the finite
size effects, which is also present for the equilibrium phase transition as
shown in Fig. 2b. The associated distribution for the phase of the order
parameter is plotted in Fig. 3c, showing that $\phi$ points along the three
Potts directions. For smaller values of $p$, the system is in a NESS and its
nature cannot be revealed by merely measuring the steady-state distributions
such as $P(\Phi)$ and $P(\Theta)$ as they show similar features as the
equilibrium case. As suggested by the result of the 3-urn model in Ref. [16],
the equilibrium state and NESS differ in the non-vanishing mean particle flux
(which corresponds to the microscopic spin direction transition rate for the
Potts spin model). The microscopic transition rates of the spin directions
from $\alpha\to\alpha+1$, $\langle K_{\alpha\to\alpha+1}\rangle$, are measured
in the simulations and the results are shown in Fig. 3d for the equilibrium
and non-equilibrium Potts models. It is clear that $\langle
K_{\alpha\to\alpha+1}\rangle$ vanishes for the equilibrium case while $\langle
K_{\alpha\to\alpha+1}\rangle$ is significantly positive for the cases of
$p=0.7$ and 0.9. In addition, $\langle K_{\alpha\to\alpha+1}\rangle$ shows a
sharp change at some value of $g$, indicating the occurrence of a non-
equilibrium phase transition.

Figure 3: Monte Carlo simulation results of the non-equilibrium ferromagnetic
3-state Potts model with $p>\frac{1}{2}$. $N=1500$. (a) The average amplitude
of the order parameter, $\langle\Phi\rangle$, is plotted against the coupling
$g$ for different values of $p$. System size of $N=3000$ is also shown for the
$p=0.9$ case. (b) The measured distribution function of $\Phi$ for $p=0.7$
near the non-equilibrium first-order transition at $g=-3.07$. The co-existing
phases are clearly revealed from the two peaks. (c) The corresponding
distribution function of $\Theta$ for the case in (b). (d) The measured
microscopic transition rates of the spin directions from $\alpha\to\alpha+1$
for the $p=0.5$ (equilibrium case, bottom line), $p=0.7$ (middle curve), and
$p=0.9$ (upper curve) cases.

As $p$ becomes large, the system is far from equilibrium and the net
transition rate of mean magnetizations from $\alpha\to\alpha+1$ becomes
larger. We first examine the case when the ferromagnetic coupling is not so
strong and the system is still in a disordered state or the uniform NESS (the
fraction of particles in the urns is uniform) as discussed in Ref. [16] in
terms of the urn model. MC simulation results for the uniform NESS in the
3-state Potts model with $p=0.9$ and $g=-2.9$ are shown in Fig. 4. The
disordered nature of the uniform state can be seen in the single peak at a
small value of $\Phi$ in the distribution of the amplitude of the order
parameter shown in Fig. 4a. The phase of the order parameter shows broad
distribution peaks slightly ahead of the three Potts spin directions of 0 and
$\pm\frac{2\pi}{3}$ as shown in Fig. 4b. More details about the uniform NESS
state can be seen from the time traces of $\Phi(t)$ and $\Theta(t)$ shown
respectively in Figs. 4c and 4d. $\Phi(t)$ displays stochastic bursts
accompanied with the oscillatory features in $\Theta(t)$. The trajectory of
$\phi(t)$ in Fig. 4e reveals that the order parameter fluctuates around with
very small magnitude most of the time with occasional stochastic periodic
counter-clockwise excursions. The time traces of the mean magnetization in the
three Potts directions are shown in Fig. 4f, indicating that the system
steadily fluctuates around the disordered (uniform) NESS. The key feature
signifying the uniform NESS nature can be seen in Fig. 4g for the time traces
of the spin transition rates measured using Eq. (13), which shows that all
three net transitions fluctuate about the same finite mean value.

Figure 4: MC simulations for the uniform NESS in the 3-state Potts model with
$p=0.9$ and $g=-2.9$, $N=1500$. (a) The measured distribution function of
$\Phi$. (b) The corresponding distribution function of $\Theta$ for the case
in (a). (c) Time course of the amplitude of the order parameter. $\Phi(t)$,
showing the stochastic bursts of the amplitude. (d) Time courses of the phase
$\Theta(t)$ in (c), showing the oscillatory feature during the amplitude
burst. Time is in the unit of Monte Carlo Steps per spin (MCS/N). (e) Time
evolution of the complex order parameter $\phi$ in its complex plane. (f) Time
evolution of the Potts spin magnetizations $x_{\alpha}(t)$. (g) Time evolution
of the mean transition rates of the Potts spin directions as measured using
Eq. (13).

As the ferromagnetic coupling grows, the NESS undergoes a first-order
transition from a disordered (uniform) NESS to an ordered (nonuniform) NESS.
Fig. 5 shows the MC simulation results for the non-uniform NESS with $p=0.9$
and $g=-3.6$. The distribution of the amplitude of the order parameter in Fig.
5a displays a major peak at $\Phi\lesssim 1$ corresponding to the order (non-
uniform) NESS. The transient states for the transitions between the three
symmetry-related NESSs can be seen from the shoulder at small values of $\Phi$
revealed in the semi-log plot. The corresponding distribution of $\Theta$ is
plotted in Fig. 5b, showing three sharp peaks slightly ahead of the three
Potts directions. The peaks in $P(\Theta)$ are much sharper as compared to the
uniform NESS indicating that the directional fluctuations of Potts spins are
much smaller for the non-uniform NESS due to the strong ferromagnetic
coupling. Details for the non-uniform NESS state can be seen from the time
traces of $\Phi(t)$ and $\Theta(t)$ shown respectively in Figs. 5c and 5d.
$\Phi(t)$ fluctuates about some large mean value with fast stochastic dips
that occur occasionally. $\Theta(t)$ reveals that these stochastic dips in
$\Phi$ are due to the collective global flipping of the spins from
$\alpha\to\alpha+1$ as shown in Fig. 5d. The above picture for the non-uniform
NESS is confirmed by the trajectory of $\phi(t)$ shown in Fig. 5e, indicating
that the order parameter points at one of three symmetry-related Potts
directions with some fluctuations for most of the time, and $\phi$ flips one
step counter-clockwise stochastically. Such a counter-clockwise stochastic
global flipping of the spins does not occur frequently due to the free energy
barrier separating the three symmetric non-uniform NESSs. By comparing with
the non-uniform state at equilibrium in Sec. 5.1 in which such a stochastic
global flipping of the spin was never observed in similar simulation
durations, one can deduce that the free energy barrier separating the three
symmetric non-uniform states is significantly lower for the NESS case. The
time traces of the mean magnetization in the three Potts directions are
displayed in Fig. 5f showing that the system predominantly fluctuates steadily
around the ordered (non-uniform) NESS, with occasional stochastic global
flipping of the spins. The NESS nature of the non-uniform state can be seen in
the time traces of the spin transition rates shown in Fig. 5g. All three net
transitions fluctuate about the same finite mean value which is somewhat lower
than that of the uniform NESS, apart from the occasional stochastic bursts due
to the global flipping of the spins between the three symmetry-related Potts
spin states.

Figure 5:  MC simulations for the non-uniform NESS with $p=0.9$ and $g=-3.6$,
$N=1500$. (a) The measured distribution function of $\Phi$. (b) The
corresponding distribution function of $\Theta$ for the case in (a). (c) Time
courses of the amplitude $\Phi(t)$, and (d) phase $\Theta(t)$ of the order
parameter. Time is in the unit of Monte Carlo Steps per particle (MCS/N). (e)
Time evolution of the complex order parameter $\phi$ in its complex plane. The
amplitude $\Phi$ and phase $\Theta$ of the complex order parameter $\phi$ are
also shown schematically. (f) Time evolution of the Potts spin magnetizations
$x_{\alpha}(t)$. (g) Time evolution of the mean transition rates of the Potts
spin directions as measured using Eq. (13).

###  5.3 NEPS and spin rotation

For sufficiently large values of $p$, a new non-steady non-equilibrium state
is observed in the non-equilibrium Potts model that was not reported in detail
in the corresponding urn model in Ref. [16]. Such a non-steady non-equilibrium
state occurs in some intermediate strength of ferromagnetic coupling and
undergoes periodic oscillation masked by thermal fluctuation, and is termed
non-equilibrium period state (NEPS). Fig. 6 shows the MC simulation results
for the NEPS with $p=0.9$ and $g=-3.5$. The distribution of the amplitude of
the order parameter in Fig. 6a displays a major peak at larger $\Phi$ and a
smaller peak at lower $\Phi$. Given the kink in the $\langle\Phi\rangle$ vs.
$-g$ curve in Fig. 3a, our result suggests there is a first-order phase
transition between the non-uniform NESS and NEPS. The order parameter
$\phi(t)$ undergoes continuous rotation on its complex plane with a non-
uniform angular velocity, slowing down around the three major Potts spin
directions, as depicted in Fig. 6b. The corresponding distribution of $\Theta$
is plotted in Fig. 6b showing three peaks slightly ahead of the three Potts
directions. The peaks in $P(\Theta)$ are less sharp as compared to the non-
uniform NESS (see Fig. 5b). The periodic nature of NEPS cannot be seen merely
from the result of the $P(\Phi)$ and $P(\Theta)$ distributions, but instead
one needs to examine in detail the time traces of $\Phi$ and $\Theta$ which
are displayed in Fig. 6d. $\Phi(t)$ shows oscillations with fluctuating
amplitudes accompanied by small fluctuations in the period as shown in Fig.
6c. The clear periodic oscillations in $\Theta(t)$ reveal the periodic nature
of the NEPS along with the continuous counter-clockwise rotation of the order
parameter. The periodic counter-clockwise rotation with a non-uniform angular
speed being slowed down around the three Potts directions is also confirmed in
the trajectory of $\phi(t)$ shown in Fig. 6e. The time traces of the mean
magnetization in the three Potts directions are displayed in Fig. 6f, showing
periodic phase-locked oscillations between the three magnetizations (or the
phase-locked periodic oscillations of particle fractions in the three urns).
Finally, the periodic dynamics also show up in the time traces of the spin
transition rates shown in Fig. 6g, displaying the same phase-locked periodic
oscillations for the three transition rates (or the phase-locked periodic
oscillations of particle fluxes between three urns).

Figure 6: MC simulations for the NEPS with $p=0.9$ and $g=-3.5$. $N=1500$. (a)
The measured distribution function of $\Phi$. The co-existing phases are
revealed from the two peaks. (b) The corresponding distribution function of
$\Theta$ for the case in (a). (c) Time courses of the amplitude $\Phi(t)$, and
(d) phase $\Theta(t)$ of the order parameter. Time is in the unit of Monte
Carlo Steps per particle (MCS/N). (e) Time evolution of the complex order
parameter $\phi$ in its complex plane. $\phi$ undergoes a continuous counter-
clockwise rotation with non-uniform angular speed in its complex plane. (f)
Time evolution of the Potts spin magnetizations $x_{\alpha}(t)$. (g) Time
evolution of the mean transition rates of the Potts spin directions as
measured using Eq. (13).

##  6 Summary and Outlook

In this paper, we established the equivalence of the multi-urn Ehrenfest model
with intra-run interactions and the Potts spin model for the equilibrium and
non-equilibrium cases. The non-equilibrium dynamics are introduced with bias
jumping rates for the clockwise and counter-clockwise spin transition
directions. Such a bias jumping rate, specified by the parameter $p$ in the
current model, might appear ad hoc at first sight, but as shown in Ref. [16]
it can be interpreted in terms of an effective chemical potential
$\mu=k_{B}T\ln(\frac{p}{q})$ that drives the Potts spin direction from
$\alpha\to\alpha+1$. Such a driving can possibly be realized by a rotating
magnetic field. Similar to the urn model, there are non-trivial non-
equilibrium phase transitions for the non-equilibrium Potts spin model,
reflecting different levels of non-equilibria. Monte Carlo simulation of the
Potts spin model is performed to investigate and confirm that the equilibrium
and non-equilibrium phase transitions observed in the Ehrenfest multi-urn
model are also observed in the equivalent Potts spin model. In particular,
there are interesting non-equilibrium phase transitions between the disordered
and ordered Potts NESSs. In addition, new NEPS exists for the far-from-
equilibrium case as shown in the simulations of the non-equilibrium Potts
model, whose detailed mechanism for its occurrence can be further investigated
theoretically in the framework of the $M$-urn model in terms of bifurcation
theory in nonlinear dynamics.

For higher Potts states with $M\geq 4$, due to the complexity in the
underlying nonlinear dynamics, it is possible that the system could display
very far-from-equilibrium dynamics such as quasi-periodic or chaotic non-
equilibrium states. The $M\geq 4$ cases are under investigation and the
results will be reported in the future. Furthermore, the associated entropy
production in these non-equilibrium states with different degrees of non-
equilibrium is another interesting issue, which can provide valuable insight
and understanding of how the energetics and dissipation are being transported.

## Acknowledgement

This work has been supported by the National Science and Technology Council of
Taiwan under grants Nos. 110-2112-M-008-026-MY3 and 111-2112-M-018-005.

## References

  * [1] R. B. Potts,  Some generalized order-disorder transformations,  Mathematical Proceedings of the Cambridge Philosophical Society, 48 (1952) 106–109. https://doi.org/10.1017/S0305004100027419
  * [2] E. Ising,  Beitrag zur theorie des ferromagnetismus,  Z. Phys., 31 (1925) 253-258. https://doi.org/10.1007/BF02980577
  * [3] F. Y. Wu,  The Potts model,  Rev. Mod. Phys., 54 (1982) 235–268. https://doi.org/10.1103/RevModPhys.54.235
  * [4] P. W. Kasteleyn and C. M. Fortuin,  J. Phys. Soc. Jpn. 26, (Suppl.) (1969) 11. https://api.semanticscholar.org/CorpusID:117993636
  * [5] M. R. Giri, M. J. Stephen, and G. S. Grest,  Spin models and cluster distributions for bond and site percolation models,  Phys. Rev. B 16 (1977) 4971–4977. https://doi.org/10.1103/PhysRevB.16.4971
  * [6] Rodney J. Baxter,  Exactly solved models in statistical mechanics,  Academic Press London, 1982.
  * [7] P.Y. Lai and Y.Y. Goldschmidt,  Application of statistical mechanics to combinatorial optimization problems: The chromatic number problem and q-partitioning of a graph,  J. Stat. Phys. 48 (1987) 513. https://doi.org/10.1007/BF01019685
  * [8] Y.Y. Goldschmidt and P.Y. Lai,  q-partitioning of graphs with finite coordination number,  J. Phys. A 21 (1988) L1043. https://dx.doi.org/10.1088/0305-4470/21/22/001
  * [9] P. Ronhovde and Z. Nussinov,  Local resolution-limit-free Potts model for community detection,  Phys. Rev. E 81 (2010) 046114. https://doi.org/10.1103/PhysRevE.81.046114
  * [10] F. Graner and J. A. Glazier,  Simulation of biological cell sorting using a two-dimensional extended Potts model,  Phys. Rev. Lett. 69 (1992) 2013–2016. https://doi.org/10.1103/PhysRevLett.69.2013
  * [11] T. Ehrenfest, P. Ehrenfest,  Uber zwei bekannte einwände gegen das Boltzmannsche H-theorem,  Z. Phys. 8 (1907) 311.
  * [12] S. Karlin and J. McGregor,  Ehrenfest urn models,  J. App. Prob. 2 (1965) 352–376. https://doi.org/10.2307/3212199
  * [13] F. P. Kelly,  The Ehrenfest model, in Reversibility and Stochastic Networks,  Wiley Chichester, 1979.
  * [14] C.-H. Tseng, Y.-M. Kao, and C.-H. Cheng,  Ehrenfest urn model with interaction,  Phys. Rev. E 96 (2017) 032125. https://doi.org/10.1103/PhysRevE.96.032125
  * [15] C.-H. Cheng, B. Gemao, and P. Y. Lai,  Phase transitions in Ehrenfest urn model with interactions: Coexistence of uniform and nonuniform states,  Phys. Rev. E 101 (2020) 012123. https://doi.org/10.1103/PhysRevE.101.012123
  * [16] C.-H. Cheng and P. Y. Lai,  Nonequilibrium thermodynamics and phase transition of Ehrenfest urns with interactions,  Phys. Rev. Res. 3 (2021) 023134. https://doi.org/10.1103/PhysRevResearch.3.023134
  * [17] C.-H. Cheng and P. Y. Lai,  Theory of non-equilibrium asymptotic state thermodynamics: Interacting Ehrenfest urn ring as an example,  Phys. Rev. E 108 (2023) 064114. https://doi.org/10.1103/PhysRevE.108.064114
  * [18] L Mittag and M J Stephen,  Mean-field theory of the many component Potts model,  J. Phys. A 7 (1974) L109. https://dx.doi.org/10.1088/0305-4470/7/9/003
  * [19] S. Dattagupta and S. Puri,  Kinetic Ising Models,  Springer Berlin Heidelberg, 2004.