ABSTRACT (original article): It is commonly believed that there are only two types of particle exchange statistics in quantum mechanics, fermions and bosons, with the exception of anyons in two dimensions. In principle, a second exception known as parastatistics, which extends outside two dimensions, has been considered but was believed to be physically equivalent to fermions and bosons. Here we show that non-trivial parastatistics inequivalent to either fermions or bosons can exist in physical systems. These new types of identical particle obey generalized exclusion principles, leading to exotic free-particle thermodynamics distinct from any system of free fermions and bosons. We formulate our theory by developing a second quantization of paraparticles that naturally includes exactly solvable non-interacting theories and incorporates physical constraints such as locality. We then construct a family of exactly solvable quantum spin models in one and two dimensions, in which free paraparticles emerge as quasiparticle excitations, and their exchange statistics can be physically observed and are notably distinct from fermions and bosons. This demonstrates the possibility of a new type of quasiparticle in condensed matter systems and—more speculatively—the potential for previously unconsidered types of elementary particle. CITATION (original article): Wang, Z., Hazzard, K.R.A. Particle exchange statistics beyond fermions and bosons. Nature 637, 314–318 (2025). https://doi.org/10.1038/s41586-024-08262-7​
GitHub: https://github.com/lagrenge94/Mathematica-codes-for-parastatistics

R-matrices, exclusion statistics, and the 1D spin model

We first construct all the nontrivial
R
-matrices that appear in this paper [Exs. 3 and 4 in Tab. I, and the set-theoretical R-matrix in Eqs. (29,31)] and verify that they satisfy the YBE (5). Then we compute the exclusion statistics and single mode partition function
z
R
(x)
for each R-matrix, using the method in Sec. S2 C. Next we construct the local operators of the 1D solvable spin model for each R-matrix, and verify that they satisfy the defining relations in Eq. (S21). The local operators
m

±

x
a
,
±

y
a

a=1
and
±

Τ
ab
are constructed in the file SpinModelOperators.wl, using the explicit matrix representation in Eq. (S24) and Eq. (S28), respectively. Note that
m

+

x
a

a=1
is stored as a
m×d×d
tensor, where
d
is the on-site Hilbert space dimension of the 1D spin model, and similarly for
m

-

x
a

a=1
,
m

+

y
a

a=1
,
m

-

y
a

a=1
.

+

Τ
ab
1<=a,b<=m
is stored as a
m×m×d×d
tensor, for example,
Τ
p
[[a,b]]
in the code gives
+

Τ
ab
, and similarly for

-

Τ
ab
1<=a,b<=m
.
Find all needed packages here: https://github.com/lagrenge94/Mathematica-codes-for-parastatistics
In[]:=
SetDirectory[NotebookDirectory[]];​​<<RMatTools`;<<TensorTool`;<<SpinModelOperators`;

Check YBE for R-matrices, and compute exclusion statistics

The set-theoretical R-matrix

Import the set - theoretical r - matrix from file, and print it
rsetth=Import["rsetth.mx"];SetthPrint[rsetth,4]//MatrixForm
Out[]//MatrixForm=
43
12
24
31
21
34
42
13
14
41
33
22
32
23
11
44
Check that rsetth satisfies the set - theoretical YBE
CheckSetthYBE[rsetth,4]
Out[]=
True
Convert the set - theoretical r - matrix to an R - matrix, and check YBE
Rsetth=-SetthToMat[rsetth,4];CheckYBE[Rsetth]
Out[]=
True
Compute the exclusion statistics
{
d
n
|0≤n≤ncutoff}
of the set - theoretical R - matrix, the last argument is ncutoff
In[]:=
ncutoff=5;ExclusionStatistics[Rsetth,ncutoff]
Out[]=
{1,4,6,4,1,0}
Compute the single mode partition function
z
R
(x)
PartitionFunction[Rsetth,ncutoff,x]//Factor
Out[]=
4
(1+x)

R-matrix in Ex .3 of Tab.I

In[]:=
R3[m_]:=-Id[
2
m
];CheckYBE[R3[4]]
Out[]=
True
In[]:=
m=4;ncutoff=3;ExclusionStatistics[R3[m],ncutoff]
Out[]=
{1,4,0,0}

R-matrix in Ex .4 of Tab.I

Construct the R - matrix in Ex .4 of Tab . I, and check YBE
σz={{1,0},{0,-1}};σx={{0,1},{1,0}};σy={{0,-},{,0}};​​λTensor[m_]:=Module{λ0=IdentityMatrix[m]},​​IfMod[m,2]0,​​λ0[[3;;m,3;;m]]=KroneckerProduct[IdentityMatrix[m/2-1],σz];,​​λ0[[1,1]]=-1;​​λ0[[2;;3,2;;3]]=MatrixExpLog
5
+3
2
σy;​​Ifm≥5,λ0[[4;;m,4;;m]]=KroneckerProductIdentityMatrix
m-1
2
-1,σz;;​​;​​λ0​​;​​R4[m_]:=Module[{c=IdentityMatrix[m],λ=λTensor[m]},KP[λ//Flatten,c//Flatten]-IdentityMatrix[m*m]//Simplify];CheckYBE[R4[3]]
Out[]=
True
In[]:=
m=3;ncutoff=4;ExclusionStatistics[R4[m],ncutoff]
Out[]=
{1,3,1,0,0}

Construct the local operators of the 1D spin model and check their algebra

The set-theoretical R-matrix

In[]:=
R=Rsetth;nmax=4;{Yp,Ym,Xp,Xm,Tp,Tm}=Simplify@ParaSpinOperators[R,nmax];​​CheckXYAlg[R,Yp,Ym,Xp,Xm]
Out[]=
True

R-matrix in Ex .3 of Tab.I

In[]:=
m=2;R=-Id[m*m];nmax=1;{Yp,Ym,Xp,Xm,Tp,Tm}=Simplify@ParaSpinOperators[R,nmax];​​CheckXYAlg[R,Yp,Ym,Xp,Xm]
Out[]=
True
In[]:=
m=3;R=-Id[m*m];nmax=1;{Yp,Ym,Xp,Xm,Tp,Tm}=Simplify@ParaSpinOperators[R,nmax];​​CheckXYAlg[R,Yp,Ym,Xp,Xm]
Out[]=
True

R-matrix in Ex .4 of Tab.I

In[]:=
m=3;R=R4[m];nmax=2;{Yp,Ym,Xp,Xm,Tp,Tm}=Simplify@ParaSpinOperators[R,nmax];​​CheckXYAlg[R,Yp,Ym,Xp,Xm]
Out[]=
True
In[]:=
m=5;R=R4[m];nmax=2;{Yp,Ym,Xp,Xm,Tp,Tm}=Simplify@ParaSpinOperators[R,nmax];​​CheckXYAlg[R,Yp,Ym,Xp,Xm]
Out[]=
True

The 2D solvable spin model

We first load the Hopf algebra
ℋ
64
from file, and verify that it satisfies the Hopf algebra axioms. The toolbox for finite dimensional Hopf algebras is stored in HopfAlgebraTools.wl. Note that in HopfAlgebraTools.wl we verify the Hopf algebra axioms step by step, for the related definitions such as prebialgebras and weak-Hopf algebras, see Ref. [3]. Then we verify that the set-theoretical R-matrix is indeed constructed from the universal R-matrix of
ℋ
64
(see, e.g., Sec.8.1.2 of Ref. [2] for the general construction), via the relation
R=[(ρ⊗ρ)ℛ]X
,
where ℛ is the universal
R
-matrix and X is the swap operator.
Then we load the tensors
u
L
,
u
R
,
v
L
,
v
R
of the 2D solvable spin model, and verify they satisfy all the CRs in Fig. S5.
Here
u
is given as
+
u
represented as a
mQ×mQ
matrix (where
m=4,Q=64
) as defined in Eq. (S31):
for
u=
u
L
,
u
R
,
v
L
,
v
R
. The function VerifyCRuv[R, S, u, v] verifies that u; v satisfies the
(R,S)CR
, by directly checking the first line of Eq. (S33), and the function VerifyCG[R, X, uL, uR, vL, vR] verifies all the CRs in Fig. S5.Next we use the tensors
u
L
,
u
R
,
v
L
,
v
R
to construct

H
1
on a
3×3
lattice, as shown in Fig. M1, using the definition in Eq. (32) in Methods. More precisely, we construct the matrix representation of

H
1
in the subspace where all black sites are in the state
|0〉
, as described in Fact. S5.2. Since there are only two white sites in the 3 lattice, the dimension of this subspace is
2
64
=4096
. We then go on to verify the claim in Fact. S5.2 that

H
1
has a unique ground state on this
3×3
lattice.The next part of the code verifies the properties of the MPO JW strings [defined in Eq. (34) in Methods] as shown in Fig. M1. Specifically, we verify that the action of the JW string on the ground state is path independent, as claimed in Fact. S5.2, and that the JW string connecting upper left and lower right corners acts on the ground state as a delta function

W
bc
|G〉=
δ
bc
|G〉
, a property we used with only a partial proof in Sec. S5 D.At last, we verify the relation of

H
1
to the definition of Kitaev’s quantum double model in the literature [1], as we explain below Fact. S5.2. More precisely, we verify that
′

A

,
′

B
p
that appears in the definition of Kitaev’s quantum double model
′

H
1
=-
∑

′

A

-
∑
p
′

B
p
in Ref. [1] are indeed projection operators to the highest eigenstates of

A

,

B
p
, and that

H
1
and
′

H
1
have the same ground state.
In[]:=
SetDirectory[NotebookDirectory[]];​​<<HopfAlgebraTools`;<<RMatTools`;<<TensorTool`;

The Hopf algebra
H
64
and construction of the set-theoretical R-matrix from the universal R-matrix of
H
64

Load the Hopf algebra H64 from data file
mul : multiplication; Δ : comultiplication; η : unit; ϵ : counit; S : antipode; R : universal R - matrix; Rb : inverse of R
HAName="H64";​​{mul,Δ,η,ϵ,S,R,Rb}=LoadHA[HAName];​​(*PrintoutthenamesofalltheavailablerepresentationsofH64*)RepNames
Verify the Hopf algebra axioms for H64, may take half a minute
Construct the R - matrix from Hopf algebra H64, and check the YBE
Load the set - theoretical r - matrix from file
Check that R4m is indeed the matrix form of the set - theoretical R (with a minus sign in front)

Load tensors uL,uR,vL,vR and verify their algebra

Import the tensors uL, uR, vL, vR from data file, each one is given as a m dA * m dA matrix .
Verify that u, v satisfy the (R, S) commutation relation
Verify that uL, uR, vL, vR satisfy the diagram of CRs in Fig . S5

The ground state properties of the quantum double model on a 3*3 lattice

Hamiltonian construction

Dimension of the parastatistical index m = 4
Dimension of the Hopf algebra H64, dA = 64
for a JW string W, take partial trace in the m - dimensional auxiliary dimension, i . e ., trace out the parastatistical index
for a JW string W, swap the tensor factor spaces of the two white sites
for w = uL, uR, vL, vR, construct the corresponding matrix acting on the m*dHil dimensional space
Construct the loop terms Av and Bp
Construct the quantum double Hamiltonian
Verify that Av, Bp commute
Check that both Av and Bp are Hermitian
Find the largest and smallest eigenvalues of Av, Bp

Obtaining the ground state using exact diagonalization

Using exact diagonalization to find the ground state*
Define the GS and verify its energy, and that it is frustration free

Verify the properties of the MPO JW strings

Define all the JW strings shown in Fig . M1 . This may take several minutes
Verify that the JW strings W1 and W2 have the same action on the state Ψ
Verify that JW1 | G >  JW2 | G >, i . e ., path independence of the JW string
Verify that JW3, JW4, JW5, and JW6 acts on GS as a delta function

Check agreement with definition in literature

Check that AvKQDM is a projection operator
Check that AvKQDM is indeed the projection operator to the highest eigenstate of Av
Check that BpKQDM is a projection operator
Check that BpKQDM is indeed the projection operator to the highest eigenstate of Av
Define KQDM Hamiltonian, and use exact diagonalization to find the ground state
Define the GS and verify its energy, and that
Verify that H1 and HKQDM have the same ground state

References

[1] O. Buerschaper, J. M. Mombelli, M. Christandl, and M. Aguado. A hierarchy of topological tensor network states. J. Math. Phys., 54(1):12201, Jan 2013.
[2] A. Klimyk and K. Schmüdgen. Quantum groups and their representations. Springer-Verlag, 1997.
[3] D. Molnar, Andras and de Alarcón, Alberto Ruiz and Garre-Rubio, José and Schuch, Norbert and Cirac, J Ignacio and Pérez-García. Matrix product operator algebras I: representations of weak Hopf algebras and projected entangled pair states. arXiv Prepr. arXiv2204.05940, 2022.

CITE THIS NOTEBOOK

Particle exchange statistics beyond fermions and bosons​
by Zhiyuan Wang & Kaden R. A. Hazzard
Wolfram Community, STAFF PICKS, January 26, 2025
​https://community.wolfram.com/groups/-/m/t/3365104