Self-organization in many-body systems with short-range interactions: clustering, correlations and topology

In this section we describe how a bosonic fluid with controllable properties, emerges from an one-dimensional (1D) Hubbard chain of strongly interacting particles, where only one particle is allowed per site. An example is shown in figure 1 for N = 6 particles with a nearest-neighbor interaction U in a chain of M sites, described by the Hubbard-like Hamiltonian term [34–36]

$$\begin{equation}{H}_{U}=U\sum _{i=1}^{M-1}{n}_{i}{n}_{i+1}.\end{equation} \tag{ 1 }$$

where ${n}_{i}={c}_{i}^{{\dagger}}{c}_{i}$ is the number operator and ${c}_{i}^{{\dagger}},{c}_{i}$ the creation and annihilation operators for a particle at site i. For simplicity in figure 1 we have drawn only the particles, ignoring the empty sites in the system. Firstly we notice that the particles self-organize at different energies condensing into clusters, resulting in the formation of various microstates at each energy with large degeneracy. Minimal Hamiltonians like equation (1) are useful to reproduce the behavior of a many-body system at the strong interaction limit when the hopping between the particles can be treated as a small perturbation which breaks the degeneracy at each energy. Although the Hamiltonian equation (1) is trivially diagonal in the Fock state basis, non-trivial features emerge as the system lies in a superposition of Fock states with rich microstructures due to the short-range interaction. The number of clusters at each energy remains constant, for example there are always four clusters at energy E = 2U. Every pair of neighboring particles (bond) contributes energy U in the total energy of each microstate. The bonds between two particles act as bosons. As the particles condense, the formed clusters can be considered effectively as states that can be filled by these bosons. The number of available bosons is simply the energy of the system in units of U. The occupation number of each effective state is given by the length of the corresponding cluster. The number of available states for filling can be found, by stacking all particles together, forming a long cluster. This cluster plus the free remaining particles gives the total number of available effective states to be filled with bosons at energy E,

$$\begin{equation}{C}^{1\mathrm{D}}\left(N,E\right)=N-E.\end{equation} \tag{ 2 }$$

By using Boson–Einstein statistics we can calculate the number of microstates at each energy, that is, the degeneracy at energy E. If we ignore the spatial gaps between the clusters, formed by the empty sites in the Hubbard chain, then we have E bosons filling C^1D sites giving the degeneracy

$$\begin{equation}{D}_{P}^{1\mathrm{D}}\left(N,E\right)=\left(\genfrac{}{}{0.0pt}{}{{C}^{1\mathrm{D}}+E-1}{E}\right)=\left(\genfrac{}{}{0.0pt}{}{N-1}{E}\right).\end{equation} \tag{ 3 }$$

If we ignore the length of the clusters, (considering them as particles) then the different ways to distribute them among M sites is [34, 36]

$$\begin{equation}{D}_{H}^{1\mathrm{D}}\left(M,N,E\right)=\left(\genfrac{}{}{0.0pt}{}{M-N+1}{N-E}\right).\end{equation} \tag{ 4 }$$

Then in order to find the total number of microstates at each energy we simply multiply equations (3) with (4) giving the total degeneracy of the system at energy E

$$\begin{equation}{D}^{1\mathrm{D}}\left(M,N,E\right)=\left(\genfrac{}{}{0.0pt}{}{N-1}{E}\right)\left(\genfrac{}{}{0.0pt}{}{M-N+1}{N-E}\right).\end{equation} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

The system described above is an emergent bosonic fluid with E free bosons represented by the corresponding bonds between the original particles. The E free bosons distribute among C^1D clusters/sites. Each cluster contains N_i particles and E_i = N_i − 1 bonds. If the total mass of the system is N and the mass of each cluster is N_i , then we can define an effective mass for each bond (quasi-particle) in a cluster, which would be ${m}_{i}=\frac{{N}_{i}}{\left({N}_{i}-1\right)}{ >}1$, satisfying ${\sum }_{i=1}^{{C}^{1\mathrm{D}}}\left({N}_{i}-1\right){m}_{i}=N$.

The wavefunction of the system at energy E will be a superposition of all the allowed microstates. In the effective boson picture, by ignoring the empty space between the clusters, we will have a wavefunction of the form

$$\begin{equation}\vert {\Psi}\rangle =\frac{1}{\sqrt{{D}_{P}^{1\mathrm{D}}\left(N,E\right)}}\left(\vert E000\cdots \rangle +\vert \left(E-1\right)100\cdots \rangle +\cdots \right),\end{equation} \tag{ 6 }$$

where each number inside the kets represents the occupation number (number of bosons) or length of each effective cluster/state. An equal superposition as in equation (6) is the simplest solution that treats all the Fock states uniformly. Other superpositions can be considered also but are beyond the scope of the current manuscript. We expect that the distribution of the amplitudes does not qualitatively affect the type of correlations emerging [34]. Exchanging two bosons in equation (6) will always lead to an even number of particle permutations in the original wavefunction of the system as shown in figure 1. Therefore the wavefunction of the emergent particles (bonds) will stay symmetric independently of the symmetry/antisymmetry under exchange of the wavefunction of the original particles. Taking account of the fact that we allow only one particle per site in the Hubbard chain equation (1), the original particles can be either hard-core bosons or fermions.

In overall, our analysis shows that an 1D quantum fluid of N fermions/hard-core bosons at energy E with short-range interactions, is equivalent to an 1D bosonic fluid with E free bosons distributed among C^1D = N − E sites. Also it is useful to note that equation (1) can be mapped to an Ising model of spins by the substitution n_i=(σ_i ^z+1)/2, where σ_i ^z=±1. This will result in H_U=∑_in_in_i+1=U∑_i(σ_i ^z+1)(σ_i+1 ^z+1)/4=(U/4)∑_i(σ_i ^zσ_i+1 ^z+σ_i ^z+σ_i+1 ^z+1), i.e. an Ising model with a nearest neighbor interaction term but also on-site magnetic field terms. Note that with open boundaries the magnetic fields on the edges of the chain are half the value in the bulk.

In order to quantify the correlations at a given energy we can use the partition method that is widely applied to calculate the entanglement in quantum many-body systems. We apply this method to calculate the correlations that arise when the system lies in a superposition of the possible microstates at energy E in the form of equation (6). We split the system in two partitions, one containing A clusters and the other one the remaining C^1D − A clusters. Then we can calculate the reduced density matrix of the partition containing A clusters by tracing out the rest of the system ρ_A ≡ tr_B |Ψ⟩⟨Ψ|, where |Ψ⟩ is equation (6). The elements of the reduced density matrix can be calculated via ${\rho }_{A}^{ij}={\sum }_{k\in B}{{\Psi}}_{ik}^{{\ast}}{{\Psi}}_{jk}$, where ${{\Psi}}_{ik}=\frac{1}{\sqrt{{D}_{P}^{1\mathrm{D}}}}$ is the amplitude for each partitioned state |ik⟩, where i(k) is the corresponding microstate inside each partition. We notice that the density matrix elements are zero unless the microstates i and j contain the same number of bonds connecting neighboring particles. Therefore, the density matrix splits in blocks each one corresponding to a fixed number of bonds m inside the partition A. Each block contains ${D}_{P}^{1\mathrm{D}}\left(A,m\right)$ elements coming from all the possible ways to arrange m bonds in A clusters as given by equation (3). In addition, each block is a full matrix with identical elements ${D}_{P}^{1\mathrm{D}}\left({C}^{1\mathrm{D}}-A,E-m\right)/{D}_{P}^{1\mathrm{D}}\left({C}^{1\mathrm{D}},E\right)$. Consequently, each block in the density matrix gives only one eigenvalue

$$\begin{equation}f\left({C}^{1\mathrm{D}},E,A,m\right)={D}_{P}^{1\mathrm{D}}\left(A,m\right)\frac{{D}_{P}^{1\mathrm{D}}\left({C}^{1\mathrm{D}}-A,E-m\right)}{{D}_{P}^{1\mathrm{D}}\left({C}^{1D},E\right)}.\end{equation} \tag{ 7 }$$

Note that the fraction in the above equation, for A = 1, gives the probability of cluster containing m bonds, which is essentially the probability distribution of the length of the clusters

$$\begin{equation}P\left({C}^{1\mathrm{D}},E,m\right)=\frac{{D}_{P}^{1\mathrm{D}}\left({C}^{1\mathrm{D}}-1,E-m\right)}{{D}_{P}^{1\mathrm{D}}\left({C}^{1\mathrm{D}},E\right)}.\end{equation} \tag{ 8 }$$

The von Neumann entanglement entropy S ≡ −tr(ρ_A log(ρ_A )) of the partition A can be calculated, by summing over the eigenvalues of each block of the reduced density matrix

$$\begin{equation}S\left({C}^{1\mathrm{D}},E,A\right)\equiv -\sum _{m=0}^{E}f\left({C}^{1\mathrm{D}},E,A,m\right)\mathrm{log}\enspace f\left({C}^{1\mathrm{D}},E,A,m\right).\end{equation} \tag{ 9 }$$

Using this method we have calculated the bipartite entropy at different energies as a function of the partition size A for N = 10⁴. The results are shown in figure 2. For low energies near the ground state (E = 0), shown in figure 2(a), all the points can be fitted by an expression ${S}_{A}^{1\mathrm{D}}=aA+b{A}^{c}+d\enspace \mathrm{log}\left(A\right)$ with four fitting parameters a, b, c, d represented by the continuous curves. Therefore, at low energies the entanglement entropy of the 1D system follows a mixture of power law and logarithmic scaling with the partition size. Despite its low energy the system does not obey the area law, S ∼ L^D−1 ⇒ S ∼ constant, expected for the ground state of gapped many-body systems. The logarithmic term however is encountered also for the critical phases in the ground state of Ising spin chains [7, 37].

Zoom In Zoom Out Reset image size

On the other hand, we have obtained a pure logarithmic scaling for high energies (E > 100) shown in figure 2(b). All points at this energy range can be fitted sufficiently well by the expression ${S}_{A}^{1\mathrm{D}}=a\enspace \mathrm{log}\left(A\right)+b$ using two fitting parameters a, b. Parameter a reaches gradually the asymptotic value a = 0.493 709 ≈ 1/2 as the energy E = 1000 is reached. A similar universal prefactor that takes values 1/3 for free bosons and 1/6 for free fermions has been shown to be related to conformal field theories effectively describing these systems [7, 37].

2.2.1. Thermodynamic limit N → ∞

In this subsection we derive analytically the scaling behavior of the entropy with partition size for a large number of particles (N → ∞). By setting $\alpha =\frac{N-E}{A}$ equation (7) becomes

$$\begin{equation}{f}^{1\mathrm{D}}\left(N,E,\alpha ,m\right)={D}_{P}^{1\mathrm{D}}\left(\frac{N-E}{\alpha },m\right)\frac{{D}_{P}^{1\mathrm{D}}\left(\left(N-E\right)\left(1-\frac{1}{\alpha }\right),E-m\right)}{{D}_{P}^{1\mathrm{D}}\left(N-E,E\right)}.\end{equation} \tag{ 10 }$$

Then we can take its thermodynamic limit (N → ∞) obtaining,

$$\begin{equation}\underset{N\to \infty }{\mathrm{lim}}\enspace {f}^{1\mathrm{D}}\left(E,\alpha ,m\right)={\left(\alpha -1\right)}^{\left(E-m\right)}{\alpha }^{-E}\left(\genfrac{}{}{0.0pt}{}{E}{m}\right).\end{equation} \tag{ 11 }$$

Then, by setting $p=\frac{1}{\alpha }$ and q = 1 − p equation (11) becomes,

$$\begin{equation}{f}^{1\mathrm{D}}\left(E,p,q,m\right)={p}^{m}{q}^{\left(E-m\right)}\left(\genfrac{}{}{0.0pt}{}{E}{m}\right).\end{equation} \tag{ 12 }$$

In the limit Epq ≫ 1 we obtain the normal (Gaussian) distribution with mean value μ and variance σ,

$$\begin{equation}{f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\approx \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\enspace {\mathrm{e}}^{-\frac{{\left(m-\mu \right)}^{2}}{2{\sigma }^{2}}}\end{equation} \tag{ 13 }$$

$$\begin{equation}\mu =E/\alpha \end{equation} \tag{ 14 }$$

$$\begin{equation}\sigma =\sqrt{\frac{E}{\alpha }\left(1-\frac{1}{\alpha }\right)},\end{equation} \tag{ 15 }$$

which satisfies,

$$\begin{equation}{\int }_{0}^{E}{f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\mathrm{d}m=1.\end{equation} \tag{ 16 }$$

The von Neumann entropy of partition A can be calculated by transforming the summation in equation (7) to an integral as,

$$\begin{equation}S\left(E,\alpha \right)=-\sum _{m}{f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\mathrm{log}\enspace {f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\end{equation} \tag{ 17 }$$

$$\begin{equation}\approx -{\int }_{0}^{E}{f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\mathrm{log}\enspace {f}^{1\mathrm{D}}\left(E,\alpha ,m\right)\mathrm{d}m.\end{equation} \tag{ 18 }$$

After some algebra, in the limit of E ≫ 1, we obtain

$$\begin{align}\hfill S\left(E,\alpha \right)& =-\frac{1}{2}\left(\mathrm{L}\mathrm{o}\mathrm{g}\left[2E\left(1-\frac{1}{\alpha }\right)\frac{1}{\alpha }\right]+\mathrm{L}\mathrm{o}\mathrm{g}\left[\pi \right]\right)\hfill \\ \hfill & -\frac{1}{4}\left(2-\frac{\frac{E}{\alpha }}{\sqrt{E\left(1-\frac{1}{\alpha }\right)\frac{1}{\alpha }}}\enspace {\mathrm{e}}^{-\frac{{\left(\frac{E}{\alpha }\right)}^{2}}{2E\left(1-\frac{1}{\alpha }\right)\frac{1}{\alpha }}}\sqrt{\frac{2}{\pi }}\right)\hfill \end{align} \tag{ 19 }$$

where α > 1. Expanding the above equation, taking also account of N ≫ E, $\alpha =\frac{N-E}{A}\approx \frac{N}{A}$ and N ≫ A we get

$$\begin{equation}S\left(N,E,A\right)=-\frac{1}{2}\enspace \mathrm{log}\left(A\right)+Q\left(N,E\right),\end{equation} \tag{ 20 }$$

with

$$\begin{equation}Q\left(N,E\right)=\frac{1}{2}\left(\mathrm{log}\left(N\right)-\mathrm{log}\left(\pi \right)-\mathrm{log}\left(2E\right)-1\right).\end{equation} \tag{ 21 }$$

A comparison between the analytical result equations (20) and (21) and the exact result using equations (7) and (9) is shown in figure 2(c). In addition from equation (21) we can see that the entanglement grows logarithmically, becoming stronger, as the energy of the system E is increased. This behavior can be seen in figure 2(d).

The prefactor 1/2 in front of the log in equation (20) can be understood mathematically as a consequence of the Gaussian form in equation (13). Physically this prefactor can be related to certain permutational symmetries present in the wavefunction of the system, which is expressed as a superposition of the different microstates allowed at energy E (equation (6)). Exchanging any two cluster lengths or alternatively any two occupation numbers in the effective boson picture, leaves the wavefunction unchanged, i.e., the wavefunction is permutationally invariant under cluster length permutations. The prefactor 1/2 is also present in the logarithmic growth of the entropy for the ground state of our system (E = 0), which contains microstates with single particles as clusters and at least one empty site between them [36]. The prefactor 1/2 for this ground state is also due to partial permutational invariance of the respective wavefunction.

The number of clusters contained in each microstate remains constant at a fixed energy. Therefore, the number given by equation (2) can be used as a topological number to characterize the microstates with the same energy. Loosely speaking the clusters or the spatial gaps between them could be considered as holes in an 1D manifold. Alternatively, each collection of microstates at a fixed energy can be thought as a set of vertex lines whose total length, determined by the energy, remains constant. Each vertex line is a path graph P_n where n is the number of its vertices.

In mathematical terms the microstates at each energy form a set of cofinite topology that has a topological invariant, the number of clusters/line segments. Any line cluster with a finite number of vertices, can be transformed to any other by the homeomorphic processes of subdivision (adding vertices) or smoothing out vertices, which is essentially like increasing or reducing the length of the clusters. Therefore all the microstates at the same energy can be continuously deformed between each other. On the other hand changes between microstates with a different number of clusters, which lie at different energies, require non-homeomorphic processes such us splitting a line segment in two. In conclusion the microstates at different energies have different topologies.

In 2D the particles form more intricate structures than in 1D, due to the additional degrees of freedom present. The system can be described by the following Hubbard Hamiltonian [35]

$$\begin{equation}H=U\sum _{x=1}^{{M}_{x}}\sum _{y=1}^{{M}_{y}}\left({n}_{x,y}{n}_{x+1,y}+{n}_{x,y}{n}_{x,y+1}\right)\end{equation} \tag{ 22 }$$

where ${n}_{x,y}={c}_{x,y}^{{\dagger}}{c}_{x,y}$ is the number operator with ${c}_{x,y}^{{\dagger}},{c}_{x,y}$ being the creation and annihilation operators for spinless particles at site with coordinates x, y, in a square lattice with periodic boundary conditions N + 1 → 1.

The self-organization of the interacting particles can be split in two clustering mechanisms. One is the same as in 1D, resulting in line clusters whose total length is determined by the energy of the system. The other mechanism is related to formation of closed vertex-lines (cycles) that is not possible in 1D. One way to help categorize all the possible structures generated by these two mechanisms, is by using the Euler characteristic χ of the graph/network structures formed by the particles as they condense. The Euler characteristic of a graph can be defined as (handshake definition)

$$\begin{equation}\chi =N-E,\end{equation} \tag{ 23 }$$

where N is the number of vertices and E the number of edges between the vertices in the graph. Each cluster will have its own number of vertices (particles) N_i , edges (bonds between particles) E_i and an Euler number χ_i = N_i − E_i . In the graph mathematical language the structures in 2D consist of either path graphs P_n , grid graphs (square lattices) P_n × P_m with N_i = nm or compositions/combinations of both [38].

If there are C^2D clusters in total at energy E, then the following expression should be satisfied

$$\begin{equation}\chi =\sum _{i=1}^{{C}^{2\mathrm{D}}}{\chi }_{i}.\end{equation} \tag{ 24 }$$

As for the 1D system, the number of clusters, along with the allowed values of χ_i , are both determined by the energy of the system and the total number of particles N. Each microstate can be described by a set of 2C^2D numbers, either, $\left\{{N}_{1},\dots ,{N}_{{C}^{2\mathrm{D}}},{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$, $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}},{N}_{1},\dots ,{N}_{{C}^{2\mathrm{D}}}\right\}$ or $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}},{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$. For our analysis we choose the set $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}},{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$. In the network/graph mathematical language, the problem consists essentially of investigating the different ways E edges can be distributed among N vertices and the resulting topological structures emerging, represented by C^2D clusters. In 1D the only possible values of χ_i are 0 when no cluster is present and 1 for the clusters that consist of either a single particle or a line of particles (vertex-line). In 2D however additional negative values of χ_i are allowed since the clusters can contain closed vertex-lines (cycles), which result in N_i ⩽ E_i . In general in both 1D and 2D, the Euler characteristic will take values in the interval

$$\begin{equation}{\chi }_{i}\in \left[{\chi }_{\mathrm{min}},1\right].\end{equation} \tag{ 25 }$$

Some examples of the two clustering mechanisms for 2D can be seen in figure 3 for various microstates in a system with N = 14, E = 15. We notice that microstates with the same energy, can contain a different number of clusters in them, unlike in 1D. For convenience, in our analysis we consider a fixed number of clusters C^2D at each energy, which counts the real clusters but also the empty clusters that have (χ_i = E_i = 0). The number of vertices N_i in each cluster constrains the corresponding number of edges E_i that can fit in the cluster, by the following equation

$$\begin{equation}{E}_{i}{\leqslant}2\left({N}_{i}-\sqrt{{N}_{i}}\right).\end{equation} \tag{ 26 }$$

The minimum value χ_min for all microstates at energy E, can be obtained by stacking all the E edges/bonds in one of the clusters. The number of vertices in this cluster can be found by using the equality in equation (26) with E_i = E,

$$\begin{equation}{N}_{1}=\frac{1}{2}\left(1+E+\sqrt{1+2E}\right).\end{equation} \tag{ 27 }$$

In case that the above equation does not give an integer value, then N₁ has to be rounded to the next largest integer. The cluster containing N₁ vertices gives the lowest possible value of the Euler characteristic at energy E,

$$\begin{equation}{\chi }_{\mathrm{min}}\left(N,E\right)={N}_{1}-E=\frac{1}{2}\left(1-E+\sqrt{1+2E}\right).\end{equation} \tag{ 28 }$$

This cluster plus the remaining free particles (vertices) gives the total number of clusters C^2D(E, N) available for filling with edges, in analogy to the idea applied in the 1D system,

$$\begin{equation}{C}^{2\mathrm{D}}\left(N,E\right)=N-{N}_{1}+1=N-\frac{1}{2}\left(-1+E+\sqrt{1+2E}\right).\end{equation} \tag{ 29 }$$

This is essentially the maximum number of clusters for each microstate at energy E. Notice that the corresponding number of vertices in 1D, for the cluster that contains all the edges is N₁ = E + 1, which is replaced by equation (27) for the 2D case. Alternatively the number of clusters can be found by C^2D = N − N₁ + 1 = N − (χ_min + E) + 1 = χ − χ_min + 1.

Zoom In Zoom Out Reset image size

In principle, by using equations (23)–(29) we can find all the possible cluster structures in the 2D system at energy E for N particles, characterized by the set of 2C^2D numbers $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}},{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$. We analyze the procedure below.

We have C^2D Euler numbers forming the set {χ₁, ..., χ_C }, which partly characterizes the cluster structures, if their corresponding number of edges characterized by the set $\left\{{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$, is ignored. The integer values of χ_i for each cluster are constrained by equations (25) and (28). The number of all possible configurations of χ can be found by distributing $\left\vert {\chi }_{\mathrm{min}}\right\vert$ bosons among C^2D states, via the binomial

$$\begin{equation}{D}_{x}^{2\mathrm{D}}\left(N,E\right)=\left(\genfrac{}{}{0.0pt}{}{{C}^{2\mathrm{D}}+\left\vert {\chi }_{\mathrm{min}}\right\vert -1}{\left\vert {\chi }_{\mathrm{min}}\right\vert }\right).\end{equation} \tag{ 30 }$$

In addition, each cluster structure for a fixed set $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}}\right\}$ can contain different distributions of edges $\left\{{E}_{1},\dots ,{E}_{{C}^{2\mathrm{D}}}\right\}$, that give the same total χ (equation (24)).

3.1.1. Clusters with χ ≠ 0

A cluster with Euler χ_i ≠ 0 will contain a minimum number of edges ${E}_{i}^{\mathrm{min}}$, which can be derived by substituting N_i = χ_i + E_i in equation (26), giving,

$$\begin{equation}{E}_{i}^{\mathrm{min}}=2\left(1-{\chi }_{i}+\sqrt{1-{\chi }_{i}}\right).\end{equation} \tag{ 31 }$$

If the above equation does not give an integer value, then it has to be rounded to the next largest integer. Notice that the same idea applies for the 1D system, but since the only possible structures are vertex-lines in this case, then we always have χ_min = 1 and E_i ⩾ 0. This means that, in 1D there are no constraints on the arrangement of the edges among the clusters. In contrast, for the 2D system the arrangement of the edges is constrained by equation (31).

3.1.2. Clusters with χ = 0

The minimum number of edges for a cluster with χ_i ≠ 0 is determined by equation (31). Clusters with χ_i = 0 have to be treated separately however, since there are two possibilities for their minimum number of edges. Euler χ_i = 0 represents either an empty cluster with N_i = E_i = χ_i = 0, or a cluster with at least one cycle (closed vertex-line) with N_i = E_i , χ_i = 0 giving a minimum number of edges ${E}_{i}^{\mathrm{min}}=4$.

If there are n clusters with χ_i = 0 then there will be 2ⁿ different distributions of minimum edges, since each cluster with χ_i = 0 can have either zero or four minimum number of edges, as stated above. If there are ${C}_{E=0}^{\chi =0}$ empty clusters and ${C}_{E{\geqslant}4}^{\chi =0}$ clusters with at least one cycle, then the number of χ ≠ 0 clusters is ${C}_{1}^{2\mathrm{D}}={C}^{2\mathrm{D}}-{C}_{E=0}^{\chi =0}-{C}_{E{\geqslant}4}^{\chi =0}$. Then there will be

$$\begin{equation}{E}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=E-\left(\sum _{i=1}^{{C}_{1}^{2\mathrm{D}}}{E}_{i}^{\mathrm{min}}+4{C}_{E{\geqslant}4}^{\chi =0}\right)\end{equation} \tag{ 32 }$$

edges which can be distributed freely among ${C}_{2}^{2\mathrm{D}}={C}^{2\mathrm{D}}-{C}_{E=0}^{\chi =0}$ clusters, which are left after excluding the empty clusters. The minimum number of edges for the χ ≠ 0 clusters ${E}_{i}^{\mathrm{min}}$ is given by equation (31). The number of different ways to distribute the free edges will be given by a summation over 2ⁿ binomials as,

$$\begin{equation}{D}^{2\mathrm{D}}\left(N,E\right)=\sum _{i=1}^{{2}^{n}}\left(\genfrac{}{}{0.0pt}{}{{C}_{2}^{2\mathrm{D}}+{E}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}-1}{{E}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}\right).\end{equation} \tag{ 33 }$$

The sum runs over the different combinations of minimum edges determined by the number of clusters with χ_i = 0. Each binomial in equation (33) corresponds to a different configuration/set of minimum edges which gives a different number for E_free and ${C}_{2}^{2\mathrm{D}}$. For a fixed number of non-empty clusters ${C}_{2}^{2\mathrm{D}}$, the free edges are dangling bonds (vertices connected with one edge) that stack on top of the clusters that contain the minimum number of edges. For example the microstate in figure 3(c) has one free dangling bond/boson which can be rearranged as shown in figure 3(g), conserving the same set of Euler numbers {−3, 1, 1, 0}. Another remark is that when there are not enough edges to fill all the χ = 0 clusters with four edges, then E_free < 0 and the respective binomial in equation (33) is zero.

From our analysis we can see that the interacting particles in 2D condense into a fluid of E bosons distributed among a number of clusters/states determined by N and E, but constrained also by the minimum number of edges/bosons that fit in each cluster. There are two major differences with the 1D system. Firstly, in 1D the number of clusters inside the microstates with the same energy is fixed. Secondly, the edges/bosons in 1D distribute freely among this fixed number of clusters/states without any constraints. Consequently, if we interpret the clusters as composite particles, then their number in 1D is conserved at a fixed energy, in contrast to 2D. A fractionalization mechanism can be easily distinguished by calculating the effective mass of the edges/bonds/quasi-particles. This is ${m}_{i}=\frac{{N}_{i}}{{E}_{i}}$. Since in 2D the cluster structures can have E_i > N_i , which is not possible in 1D, each edge in the cluster has a fractionalized effective mass m_i < 1.

We remark that in principle, the core of the above results could be extended to any random network that can be arranged as a planar graph, one that contains no overlapping edges, formed by the interactions between many particles. In general, the total number of bonds/edges between two particles determines the energy of the system. These edges can be distributed among N particles/vertices giving many different microstates with diverse topological/geometrical cluster structures, similar to the ones that we have demonstrated specifically for the square lattice, which is a grid graph.

The number of cycles (closed vertex lines) in each cluster can be used as a topological invariant to characterize the different cluster configurations, corresponding to different sets $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}}\right\}$. For instance, the cycles are five, four and three for the cluster residing on the left side of figures 3(a)–(c). The number of cycles can be changed for example, by removing one bond which opens/breaks a cycle, creating dangling bonds (vertices connected with one edge). This act corresponds to a non-continuous deformation analogous to the process of pocking a hole through a sphere to create a torus, two topologically distinct geometrical shapes. The Euler number is analogous to the number of cycles in a cluster. This can be understood by starting with one cycle with χ = 0. Every additional cycle in the cluster will change χ by V − E = V − (V + 1) = −1 resulting in

$$\begin{equation}\chi =-\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{s}+1.\end{equation} \tag{ 34 }$$

The above formula is consistent with Euler's formula for planar graphs V − E + F = 2 where F are the faces in the graph. For our case F = cycles + 1 since F counts the faces enclosed by the square lattice cluster structures (cycles), along with the external face in the surface where the clusters are embedded. Since χ = V − E we have χ + cycles + 1 = 2 ⇒ χ = −cycles + 1. This formula is analogous to the respective formula that connects the Euler characteristic and the genus (number of holes) g in topological objects described by differential geometry χ = 2(−g + 1). Geometrical shapes with the same g correspond to the same topology.

Microstates that are described by the same set $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}}\right\}$ have the same topology and the corresponding cluster structures can be continuously deformed between each other, for example the clusters in figures 3(e) and (i). Note that when the cluster contains only one cycle, like the rightmost cluster in figure 3(e), it can be shrunken down to an empty cluster, which is understood topologically by taking a circle and shrinking it's radius to zero. Different sets $\left\{{\chi }_{1},\dots ,{\chi }_{{C}^{2\mathrm{D}}}\right\}$ correspond to cluster structures with different topologies. Notice that the number of clusters inside each microstate at each energy is not constant, unlike the 1D case, and therefore it cannot be used as a topological invariant. In overall, we see that while in the 1D system each sub-Hilbert space of microstates at energy E defines a topology characterized by the number of clusters, in 2D each sub-Hilbert space at energy E contains microstates with different topologies, characterized by the set of Euler numbers describing the corresponding cluster structures.

We remark that equation (34) is valid for any planar graph, one with no overlapping edges, even random graphs formed by interactions between many particles. In this sense the topological properties described above are valid for any many-body network that can be arranged as a planar graph.