Kramers–Wannier Duality and Random-Bond Ising Model

Abstract

We present a new combinatorial approach to the Ising model incorporating arbitrary bond weights on planar graphs. In contrast to existing methodologies, the exact free energy is expressed as the determinant of a set of ordered and disordered operators defined on a planar graph and the corresponding dual graph, respectively, thereby explicitly demonstrating the Kramers–Wannier duality. The implications of our derived formula for the Random-Bond Ising Model are further elucidated.

1. Introduction

The well-established duality between order and disorder phases observed in the two-dimensional Ising model was initially exploited by Kramers and Wannier to pinpoint its criticality [1], predating the celebrated solution of Onsager’s free energy [2,3,4]. Furthermore, Kadanoff and Ceva illustrate that the correlation function can be derived by contemplating the disorder operator defined on the dual graph [5], providing substantial insights into the underlying physics [6]. However, the standard methodologies used for the computation of free energy—either algebraically or combinatorially [7]—exhibit the Kramers–Wannier (KW) duality only in the final form after extensive calculations. Despite a tremendous volume of work over the last century devoted to the exact solution of the Ising model, a formula for calculating its free energy with manifest KW duality is still absent in the literature. This gap limits the utility of the exact free energy in broader contexts. For example, in the case of a Random-Bond Ising Model (RBIM) employed for understanding spin glasses, deriving an explicit free energy remains challenging.

It is worth mentioning that the combinatorial approach, pioneered by Kac and Ward [8], provides an alternative pathway to derive Onsager’s free energy. In particular, for a planar graph $G=(V,E)$ with $n=\left|V\right|$ and $m=\left|E\right|$, Kac and Ward introduced a determinant formula,

$$\begin{array}{c}\hfill {\zeta }_F{(G,u)}^{-1}:=det(I_{2m}-u{T}_{KW})=Z_{\mathrm{ising}}^2{(1-u^2)}^m,\end{array}$$

(1)

where $u:=tanh\left(\beta J\right)$ represents the coupling constant, and the Kac–Ward (KW) operator $T_{KW}$ is a $2m\times 2m$ matrix, further elucidated in the subsequent discussion. The original derivation of the Kac–Ward formula (1) was based on heuristic arguments and lacked logical completeness. To address this deficiency, Feynman, in unpublished work cited by Harary [9,10], postulated a path-integral form for the Kac–Ward determinant as

$${\zeta }_F{(G,u)}^{-1}=\prod _{\left[p\right]}\left(1-{(-1)}^{w\left(p\right)}{u}^{l\left(p\right)}\right),$$

(2)

which allows to prove (1) directly. This elegant identity draws an analogy with Riemann’s zeta function, where the Euler product is over all prime cycles, with $l\left(p\right)$ and $w\left(p\right)$ representing the length and winding number of the prime cycles p, respectively. The subscript F highlights the fermionic character of the Ising model, which assigns a negative weight to odd numbers of windings. This formula was later formally proved by Sherman [11,12,13] and Burgoyne [14]. A variant of the combinatorial formulation, mapping the Ising model to the dimer model using Pfaffians, was developed by Green and Hurst [15] and later expanded by others [16,17,18,19]. This approach essentially corresponds to a skew-symmetric version of Equation (1) via a similarity transformation. More recently, a resurgence of the combinatorial approach [20,21,22] focuses on the discrete version of the conformal invariance of the critical Ising model on planar graphs [23].

The combinatorial approach can be seamlessly generalized to accommodate an arbitrary set of bond weights $\mathbf{u}:=\left\{u_e\right|e\in E\}$, thereby presenting a robust numerical tool for probing RBIM. However, Equation (1) reveals little physical insight. For example, the KW duality indicates that ${\zeta }_F$ should remain invariant (up to some prefactor) under the transformation $u\to u^{\ast }:=(1-u)/(1+u)$ on the dual graph $G^{\ast }$, i.e.,

$${\zeta }_F(G,u)\sim {\zeta }_F(G^{\ast },u^{\ast }).$$

Regrettably, the manifestation of the KW duality only emerges following the resolution of the determinant, a process that is impractical for numerous disordered systems. This hidden symmetry within Equation (1) remains far from obvious and poses a significant challenge. It was only recently proven that Equation (1) indeed satisfies KW duality under general conditions [24,25]. Consequently, despite its elegance, the combinatorial approach is primarily utilized as a numerical tool. To the best knowledge of the author, no explicit free energy formula demonstrating manifest KW duality for an arbitrary planar graph and weight set has yet been identified.

In this paper, we propose a new free energy formula for the Ising model with arbitrary bond weights on planar graphs. Contrasting with existing methodologies, our formula is expressed as the determinant of the summation of local ordered and disordered operators, each defined on vertices V and dual vertices $V^{\ast }$, thereby explicitly exhibiting the KW duality. In addition, it establishes a tangible connection with nonlocal ordered and disordered operators, offering insights into the nature of duality. We elucidate the implications of our formula in the context of RBIM.

2. Ihara Zeta Function

To hint at the existence of the manifest dual formula, we initiate our discussion with a warm-up exercise by considering the bosonic counterpart of Equation (2),

$$\begin{array}{c}\hfill {\zeta }_B{\left(u\right)}^{-1}=\prod _{\left[p\right]}\left(1-u^{l\left(p\right)}\right).\end{array}$$

(3)

This definition, known as the Ihara zeta function [26,27], serves as a p-adic analogue of the Selberg zeta function that counts the number of closed geodesics on a hyperbolic surface. Analogous to Equation (1), we express

$${\zeta }_B{\left(u\right)}^{-1}=det(I_{2m}-uT),$$

(4)

where T denotes the $2m\times 2m$ Hashimoto’s edge adjacency operator [28] in analogy to the KW operator $K_{KW}$. Specifically, T applies on the space of oriented edges $\{E,\overline{E}\}$, where an edge $\overline{e}\in \overline{E}$ denotes the directional inverse of a corresponding edge $e\in E$. The matrix $T_{e^{\prime },e}=1$ only if the oriented edge e follows $e^{\prime }$ backtracklessly, meaning that the terminal vertex of e is the starting vertex of $e^{\prime }$ and $e^{\prime }\ne \overline{e}$. The equivalence of Equations (3) and (4) can be demonstrated directly by applying the logarithm and aligning the power expansion term by term.

For a regular graph of degree $q+1$, the Ihara zeta function displays self-duality under the transformation $u\to q/u$, which mirrors Riemann’s functional equation. However, in similarity to its fermionic counterpart, Equation (4) does not explicitly reveal this duality. Intriguingly, a second formula exists, as proposed in Ihara’s original paper [26],

$$\begin{array}{c}\hfill {\zeta }_B{\left(u\right)}^{-1}={(1-u^2)}^{m-n}det(I_n-uA+u^{2}Q),\end{array}$$

(5)

where A is the adjacency matrix, and Q is the diagonal matrix with the degree diminished by one. Assuming the graph is $(q+1)$-regular, that is, $Q=q{I}_n$, it becomes evident that

$${\zeta }_B^{-1}(G,u)\sim {\zeta }_B^{-1}(G,q/u).$$

Consequently, Equation (5) presents a manifestly dual formula for the bosonic zeta function.

The derivation of Equation (5) from Equation (4) provides illuminating insights. Here, we present a streamlined approach based on the original proof by Bass [29,30,31]. We introduce the matrix

$$S:=T+J,$$

where $J_{e^{\prime },e}={\delta }_{e^{\prime },\overline{e}}$ and $S_{e^{\prime },e}$ enumerate all successors $e^{\prime }$ following e, including its inverse $\overline{e}$. A key observation arises from the factorization of the matrix

$$S=Y^{t}X,$$

where X and Y are $n\times 2m$ matrices with $X_{v,e}=1$ or $Y_{v,e}=1$ if v is the starting vertex or the terminal vertex of the orient edge e, respectively. Leveraging this factorization, we obtain

$$\begin{array}{cc}\hfill det(& I_m-uT)=det\left((I_{2m}+uJ)-uS\right)\hfill \\ \hfill & ={(1-u^2)}^{m}det\left(I_n-uX{(I_{2m}+uJ)}^{-1}{Y}^t\right)\hfill \\ \hfill & ={(1-u^2)}^{m-n}det\left(I_n-uA+u^{2}Q\right),\hfill \end{array}$$

where the second line ensues from the generalized matrix determinant lemma and $det(I_{2m}+uJ)={(1-u^2)}^m$. The third line employs the identity

$${(I_{2m}+uJ)}^{-1}={(1-u^2)}^{-1}(I_{2m}-uJ),$$

while noting $A=X{Y}^t$ and $Q=XJ{Y}^t-I_n$. This completes the proof of Equation (5). The critical component of this proof involves the use of the generalized matrix determinant lemma, predicated on the factorability of S, which can be reinterpreted as an index theorem over a chain complex [32].

Consider a more generalized setup involving an arbitrary set of weights $\mathbf{u}:=\left\{u_e\right|e\in E\}$ assigned to each edge. By employing a similar approach, we obtain

$$\begin{array}{c}\hfill {\zeta }_B{\left(\mathbf{u}\right)}^{-1}=\prod _{e\in E}(1-u_e^2)det\left(I_n-\tilde{A}\left(\mathbf{u}\right)+\tilde{D}\left(\mathbf{u}\right)\right),\end{array}$$

(6)

where $\tilde{A}$ represents the weighted adjacency matrix defined as

$${\tilde{A}}_{v,v^{\prime }}:={\displaystyle \frac{u_{v{v}^{\prime }}}{1-u_{v{v}^{\prime }}^2}},$$

and $\tilde{D}$ denotes the weighted degree defined as

$$D_{v,v}=\sum _{(v,v^{\prime })\in E}{\displaystyle \frac{u_{v{v}^{\prime }}^2}{1-u_{v{v}^{\prime }}^2}}.$$

Specifically, for the $\pm J$ disorders, i.e., $u_e=u{\tau }_e$ with ${\tau }_e=\pm 1$, Equation (6) simplifies to

$$\begin{array}{c}\hfill {\zeta }_B{\left(u\right)}^{-1}={(1-u^2)}^{m-n}det(I_n-u{A}^{\prime }+u^{2}Q),\end{array}$$

where $A^{\prime }$ includes entries of 0 and $\pm 1$ to account for bond disorders. It becomes evident that $A^{\prime }+I$ serves as the adjacency matrix of the percolation model, thereby mapping the Ihara zeta function with $\pm J$ disorder to a percolation problem.

3. Manifestly Dual Formula

Consider the fermionic case represented by Equation (2), where the winding number is suitably defined only after immersion into a surface with a given spin structure. This marks a notable distinction from its bosonic counterpart in Equation (3) and creates a host of technical challenges when applying a similar approach to the one used for the bosonic zeta function. For ease of discussion, our examination is confined to a planar graph G embedded in a plane. However, extending this discussion to surfaces with a higher genus is straightforward.

To explicitly reveal the KW duality, we embed the dual graph $G^{\ast }$ over G. In this arrangement, each vertex of G is located inside a face of $G^{\ast }$, and vice versa; each edge in G intersects with the corresponding edge in $G^{\ast }$. For technical convenience, we require these intersections to be perpendicular. Note that the embedding need not be isoradial in general. Figure 1a illustrates such an embedding, where red and blue colors represent G and $G^{\ast }$, respectively.

A crucial element in our depiction involves the quadrilaterals (gray domain in Figure 1a). Each quadrilateral is formed by two neighboring edges in both G and $G^{\ast }$, along with a vertex pair $(v,v^{\ast })$, where $v\in V$ and $v^{\ast }\in V^{\ast }$. We denote the angles associated with the left and right edges of v and $v^{\ast }$ as ${\theta }_L$ and ${\theta }_R$, and ${\theta }_L^{\ast }$ and ${\theta }_R^{\ast }$, respectively. These angles satisfy the relations

$${\theta }_L+{\theta }_R^{\ast }={\theta }_R+{\theta }_L^{\ast }=\pi /2,$$

as depicted in Figure 1a. The collection of $2m$ quadrilaterals, which collectively tile the entire plane, plays a critical role in our new formulation.

The Kac–Ward operator $T_{KW}$ in Equation (1) is defined similarly to Hashimoto’s edge adjacency operator T. However, we must introduce a phase change between two consecutive edges to account for the fermionic nature. Specifically, we impose a gauge transformation

$$\begin{array}{c}\hfill {\left(T_{KW}\right)}_{e^{\prime },e}=e^{i\alpha (e^{\prime },e)/2}=i{e}^{-i\beta (\overline{e},e^{\prime })/2}\end{array}$$

if edge $e^{\prime }$ follows e without backtracking, where $\alpha (e^{\prime },e)$ is the exterior angle from e to $e^{\prime }$, and $\beta (\overline{e},e^{\prime })$ is the interior angle between $e^{\prime }$ and its inverse $\overline{e}$ (Figure 2a). This condition ensures that the summation of the half exterior angles contributes to a total $\pi$ phase change over a cycle, effectively capturing the fermionic sign in Equation (2).

Following a similar approach as that applied to the bosonic zeta function, we introduce the gauged successor operator

$$S^{\prime }:=T_{KW}-iJ,$$

appending an additional element between edge e and its inverse $\overline{e}$ with weight $i{e}^{-i\pi }=-i$. However, $S^{\prime }$ no longer exhibits factorability, preventing directly applying the matrix determinant lemma in this scenario. To tackle this issue, we follow the strategy developed by Cimasoni [24,25], introducing operator Q between neighboring edges e and $e^{\prime }$ that share a common starting vertex v. Specifically, for each quadrilateral q, the operator Q maps its right edge e to the left edge $e^{\prime }$ while inducing a phase shift

$$Q_{e,e^{\prime }}:=e^{i({\theta }_L+{\theta }_R)/2}.$$

As illustrated in Figure 2b,

$$Q_{e^{\prime },e^{″}}{e}^{-i\beta (e^{″},e)/2}=e^{-i\beta (e^{\prime },e)/2}$$

if $e^{\prime }\ne \overline{e}$. However, the fermionic nature possesses a non-trivial monodromy, resulting in a branch cut after a rotation of $2\pi$. This observation leads to the discontinuity

$$\begin{array}{c}\hfill S^{\prime }-Q{S}^{\prime }=-2iJ.\end{array}$$

(7)

As the operator Q acts on the edges attached to the quadrilaterals, it facilitates a natural factorization

$$Q=L^{t}R,$$

where L and R are $2m\times 2m$ matrices associating each quadrilateral q with its left and right edges, respectively. Moreover, we assign weights $L_{q,e}:=e^{i{\theta }_L\left(q\right)/2}$ and $R_{q,e}:=e^{i{\theta }_R\left(q\right)/2}$. Building on this factorization, we find

$$\begin{array}{cc}\hfill & det(I-Q)det(I-u{T}_{KW})=det(I+iuJ-Q(I-iuJ))\hfill \\ \hfill & ={(1+u^2)}^{m}det(I-R(I-iuJ){(I+iuJ)}^{-1}{L}^t)\hfill \\ \hfill & ={(1+u^2)}^{m}det\left(I-{\displaystyle \frac{1-u^2}{1+u^2}}R{L}^t-{\displaystyle \frac{2u}{1+u^2}}R{e}^{-i\pi /2}J{L}^t\right),\hfill \end{array}$$

where the first equality follows from Equation (7), and the second stems from the generalized matrix determinant lemma and $det(I+iuJ)={(1+u^2)}^m$. We then introduce the discrete curl operator

$$d:=R{L}^t,$$

acting on the space of quadrilaterals, with elements

$$\begin{array}{c}\hfill d_{q^{\prime },q}=e^{i({\theta }_L\left(q\right)+{\theta }_R\left(q^{\prime }\right))/2},\end{array}$$

when the quadrilaterals $q^{\prime }$ and q share a common edge with $q^{\prime }$ positioned counterclockwise next to q. This definition suggests that d can be decomposed into a set of operators $d_v$, each acting on the quadrilaterals around vertex v. Therefore, $d={\sum }_{v\in V}{d}_v$, as depicted in Figure 1b. Similarly, the dual operator $d_{\ast }={\sum }_{v^{\ast }\in V^{\ast }}{d}_{v^{\ast }}$ is the summation of the operators $d_{v^{\ast }}$ around the dual vertex $v^{\ast }$. Observe that

$${\left(d_{\ast }^{†}\right)}_{q^{\prime },q}=e^{-i(\pi /2-({\theta }_L\left(q\right)+{\theta }_R\left(q^{\prime }\right))/2}$$

holds if q possesses a left edge e and $q^{\prime }$ a right edge $\overline{e}$ (see Figure 1). This yields

$$d_{\ast }^{†}=R{e}^{-i\pi /2}J{L}^t,$$

which applies to quadrilaterals around dual vertices clockwise. Taken together, we obtain

$$\begin{array}{c}\hfill {\zeta }_F^{-1}=2^{-n}{(1+u^2)}^{m}det\left(I_{2m}-{\displaystyle \frac{1-u^2}{1+u^2}}d-{\displaystyle \frac{1-{u^{\ast }}^2}{1+{u^{\ast }}^2}}{d}_{\ast }^{†}\right),\end{array}$$

(8)

where we employ the identities $det(I-W)={(1-(-1))}^n=2^n$ and

$${\displaystyle \frac{1-{u^{\ast }}^2}{1+{u^{\ast }}^2}}={\displaystyle \frac{2u}{1+u^2}}.$$

Drawing parallels with Equation (6), we can generalize Equation (8) to incorporate a set of bond weights $\mathbf{u}$ as follows:

$$\begin{array}{c}\hfill {\zeta }_F^{-1}(G,\mathbf{u})=2^{-n}\prod _{e\in E}(1+u_e^2)det\left(I_{2m}-D\left(\mathbf{u}\right)-D_{\ast }^{†}\left({\mathbf{u}}^{\ast }\right)\right),\end{array}$$

(9)

where

$$D:=R{\displaystyle \frac{1-{\mathbf{u}}^2}{1+{\mathbf{u}}^2}}{L}^t,$$

and $D_{\ast }$ represent the weighted curl operators around the vertices in V and $V^{\ast }$, respectively. These operators permit a natural decomposition

$$\begin{array}{c}\hfill D\left(\mathbf{u}\right)=\sum _v{D}_v\left(\mathbf{u}\right),D\left({\mathbf{u}}^{\ast }\right)=\sum _{v^{\ast }}{D}_{v^{\ast }}\left({\mathbf{u}}^{\ast }\right).\end{array}$$

(10)

The determinant present in Equation (9) is manifestly symmetric under the dual transformation, thereby reinstating the KW duality for an arbitrary set of bond weights [24,25]

$$2^{-\left|V\right|}\prod _{e\in E}(1+u_e){\zeta }_F(G,\mathbf{u})=2^{-|V^{\ast }|}\prod _{e\in E}(1+u_e^{\ast }){\zeta }_F(G^{\ast },{\mathbf{u}}^{\ast }),$$

(11)

where we applied the identity

$${\displaystyle \frac{1+u_e^2}{1+{u_e^{\ast }}^2}}={\displaystyle \frac{1+u_e}{1+u_e^{\ast }}}.$$

4. Order and Disorder Operators

Equation (9) suggests that the fermionic zeta function

$${\zeta }_F^{-1}=det(I-H\left(\mathbf{u}\right))$$

constitutes the characteristic polynomial of a non-Hermitian Hamiltonian:

$$\begin{array}{c}\hfill H:=D\left(\mathbf{u}\right)+D^{†}\left({\mathbf{u}}^{\ast }\right)=\sum _{v\in V}{D}_v\left(\mathbf{u}\right)+\sum _{v^{\ast }\in V^{\ast }}{D}_{v^{\ast }}^{†}\left({\mathbf{u}}^{\ast }\right).\end{array}$$

Assuming the $\pm J$ disorders of $u_e=u{\tau }_e$, where ${\tau }_e=\pm 1$, the equation simplifies to:

$$\begin{array}{c}\hfill H={\displaystyle \frac{1-u^2}{1+u^2}}\sum _{v\in V}{d}_v+{\displaystyle \frac{2u}{1+u^2}}\sum _{v^{\ast }\in V^{\ast }}{\tilde{d}}_{v^{\ast }}^{†},\end{array}$$

(12)

where

$${\left({\tilde{d}}_{v^{\ast }}\right)}_{q,q^{\prime }}=e^{i({\theta }_L\left(q\right)+{\theta }_R\left(q^{\prime }\right))/2}{\tau }_e,$$

applicable for two quadrilaterals sharing the dual vertex $v^{\ast }$ and edge e. Evidently, only the dual curl operator $D_{\ast }$ represents the disorder, while the curl operator D remains unaffected by the disorder. Thus, we interpret $D_v$ and $D_{v^{\ast }}$ as local order and disorder operators, respectively, echoing the nonlocal disorder operator introduced by Kananoff and Ceva [5].

To establish a connection explicitly, consider a defect that changes a line $\mathcal{L}$ of ferromagnetic bonds to antiferromagnetic, i.e., ${\tau }_e=-1$ only for $e\in \mathcal{L}$. The corresponding correlation function of nonlocal disorder operators involves a shift in the free energy:

$$\begin{array}{c}\hfill \Delta F=-{\displaystyle \frac{kT}{2}}lndet(I-H^{\prime }G),\end{array}$$

(13)

where $G={(I-H_0)}^{-1}$ represents Green’s function of the ferromagnetic Hamiltonian

$$H_0:={\displaystyle \frac{1-u^2}{1+u^2}}d-{\displaystyle \frac{1-{u^{\ast }}^2}{1+{u^{\ast }}^2}}{d}_{\ast }^{†},$$

and the defect operator

$$H^{\prime }:={\displaystyle \frac{2u}{1+u^2}}\sum _{v^{\ast }\in \mathcal{L}}({\tilde{d}}_{v^{\ast }}-d_{v^{\ast }}).$$

It becomes apparent that nonlocal disorder operators correspond to a line integral over the local disorder operator $D_{v^{\ast }}$. Consequently, the KW duality presented in Equation (9) reflects an exact interchange of local order and disorder operators under the duality transformation.

5. Implication to RBIM

We now turn to the implications of our new formula for the RBIM. For the sake of technical convenience, our discussion will primarily focus on $\pm J$ disorder on a square lattice, i.e.,

$$P\left(\tau \right)=(1-p)\delta (\tau -1)+p\delta (\tau +1).$$

A straightforward generalization applies to arbitrary disorder. As we demonstrated earlier, the free energy of $\pm J$ bond disorder is dictated by the spectrum of Equation (12), where only the disorder operator ${\tilde{d}}_{\ast }$ accounts for the bond disorders ${\tau }_e$.

In the absence of disorder ($p=1$), it is straightforward to diagonalize H via the Fourier transform, which consequently recovers Onsager’s free energy. Considering a scenario where p is close to 1, we can apply the Dyson series expansion to Equation (13). This approach enables us to determine the critical coupling $u_c\left(p\right)$ as a series expansion of disorder probability $1-p$. At the leading order, we find

$$u_c\left(p\right)=(\sqrt{2}-1)(1+2\sqrt{2}(1-p))+O\left({(1-p)}^2\right).$$

(14)

This result aligns with the findings initially obtained using the replica trick [33].

We now turn our attention to the zero-temperature limit $\beta \to \infty$ of Equation (12). In this limit, $(1-u^2)/(1+u^2)\approx 2e^{-2\beta J}$ and $2u/(1+u^2)\approx 1$. Consequently, in this scenario, disorder operators dominate the spectrum. Direct computation yields

$$\begin{array}{c}\hfill det(I-d_{v^{\ast }}^{†})=1+W_{v^{\ast }},\end{array}$$

(15)

where the frustration $W_{v^{\ast }}:={\prod }_{e\in P\left(v^{\ast }\right)}{\tau }_e=\pm 1$ is defined as the product of edge disorders around the plaquette of the corresponding dual vertices $v^{\ast }$. It is clear that when the plaquette is frustrated, i.e., $W_{v^{\ast }}=-1$, the determinant acquires a correction from the order operator d with an order at most $e^{-2\beta J}$. On average, there is a ${\displaystyle \frac{1-{(2p-1)}^4}{2}}$ chance of $W_{v^{\ast }}=-1$, which results in a lower bound of the ground state energy density:

$$e/J\ge -2+{\displaystyle \frac{1-{(2p-1)}^4}{2}}.$$

(16)

This finding is consistent with the results in Refs. [34,35] obtained using geometric approaches.

6. Conclusions

In conclusion, we have unveiled a novel combinatorial approach to Ising models with arbitrary bond weights. In contrast to previous methods, our new formulation distinctly manifests the KW duality via order and disorder operators. We have presented preliminary implications for RBIM at the leading order, and our findings are consistent with results derived from alternative approaches. However, our method has the distinct advantage of seamlessly integrating with the standard framework of perturbative techniques, thereby simplifying the extension to higher-order terms. This also opens up the potential to employ non-perturbative methodologies for a more nuanced understanding of the phase diagram of RBIM. Additionally, our approach can be directly applied to other planar graphs, such as triangular and hexagonal lattices. On the other hand, it has been suggested that the RBIM may exhibit a disorder duality based on the replica argument, potentially localizing the multicritical point [36]. The exactness of this duality and its connection to our method remains unclear. We aim to explore these questions in future research.

Our methodology can also be readily generalized to anyonic statistics by considering a non-half-integer phase shift, a topic we plan to discuss elsewhere. Further, a non-Abelian generalization seems feasible. These generalizations have close ties with parafermionic models [37]. Moreover, given that the Ihara zeta function can generalize to high-dimensional objects [38], it is enticing to contemplate a similar higher-dimensional generalization for its fermionic counterpart. Such an extension may hold promise for a solution to the 3D Ising model.