Bose-Fermi Mapping in Hubbard Models at Imaginary Chemical Potential and Phase-Induced Fermionization

1. Introduction

The Hubbard model is one of the most interesting models in condensed matter physics, because it captures the complex competition between kinetic energy and local interactions. Usually the two specific kinds of the model, the attractive Fermi-Hubbard model and the repulsive Bose-Hubbard model, are treated as two distinct theories since they describe fundamentally different physics: the first helps us to describe pairing and superfluidity, while the other captures Mott insulation and Bose-Einstein condensation [1,2,3,4,5,6,7]. However, despite their different natures, both models exhibit rich crossover phenomena driven by temperature, interaction strength, and external fields. In this work we are trying to answer a fundamental question: whether these two systems, operating under different statistics, can be connected through a common mathematical framework?

Since the extension of fermion-boson mapping [8,9,10,11,12] from zero to finite temperature is always a recurrent theme in high energy and condensed matter physics, we demonstrate an exact duality between the two models, at finite temperature and in the presence of an imaginary chemical potential $\mu =i\theta$. Although this potential is not directly observable in lattice theory, it corresponds to a background $U\left(1\right)$ gauge field creating a holonomy along the thermal circle, analogous to Aharonov-Bohm phases [13]. This temporal gauge field picture connects our analysis to a variety of phenomena where particles acquire phases from moving around closed loops and helps us to connect statistical mechanics to topological phase effects.

At the heart of this work, as a key insight, is the relation between fermionic and bosonic Matsubara frequencies, which leads to an exact duality of the corresponding partition functions. Consequently, this duality ensures that the thermal kernels of the two models are related by a simple shift and a main result is the emergence of universal thermal windows at specific angles at $\varphi =\pi /3$ and $5\pi /3$ for the bosonic model. This observation allows us to ask whether these boundaries map exactly the corresponding boundaries from the attractive Fermi-Hubbard model, and indeed this happens in a remarkable and symmetric way.

While at the Tonks–Girardeau limit [14,15,16,17], where bosons fermionize due to an infinitely strong repulsive interaction, our results suggest that a similar fermion-like behaviour can emerge at interaction of finite strength when the system is at an imaginary chemical potential. Having this idea in mind, we see that fermionization is not driven by a hard-core constraint on the wavefunction [18], but by a redistribution of spectral weight. In this picture, the imaginary chemical potential suppresses low-energy occupation and enhances higher-energy modes within a finite window $\varphi \in [\pi /3,5\pi /3]$ ($\pi /3$ or ${\varphi }_{\mathrm{flux}}=1/6$ is special according to [19,20,21] since this value appear in lattice bosonic systems with real magnetic flux, where it controls strongly correlated and topological states and as a lattice phase geometry, relative to the boundaries), producing an effective exclusion-like behaviour at the level of thermal occupations. This reveals a regime where, despite the lack of a formal Bose-Fermi mapping at the operator level [22], the thermal occupations mimic fermionic behaviour through an interplay of interaction and entropic pressure. On the other hand, dimensionality enters only through the density of states $\rho \left(ϵ\right)$, allowing the framework to be extended to various lattice geometries in the future, including 1D chains, 2D square lattices (with van Hove singularities), and 2D honeycomb lattices (with Dirac cones).

This paper is organized as follows. In Section 2, we present the repulsive Bose-Hubbard model with an imaginary chemical potential and derive its effective action using a large-N expansion. Section 3 derives the gap equation, number equation, and the fermionization-like window of the bosonic model. Section 4 examines the phase structure through the thermal kernel and the extrema of number $N_b$ and section 5 reviews the main results for the attractive Fermi-Hubbard model at imaginary chemical potential and large-N. Section 6 establishes the duality between the attractive Fermi and repulsive Bose models, including the mapping of thermal kernels. In Section 7 we derive a universal statistical transmutation framework for the two models and we conclude in Section 8 with a summary and conclusions. In Appendix we have some useful calculations.

2. Repulsive Bose-Hubbard Model and Large- $\mathit{N}$ Formulation

We consider the spinless Bose-Hubbard model on a lattice at inverse temperature $\beta$, with a fixed total boson number $N_b$, and $n_i=b_i^{†}{b}_i$, so the canonical partition function is [9,10,13,23]

$$Z_{N_b}=\mathrm{Tr}\left[\delta (\widehat{N}-N_b)e^{-\beta H}\right],\widehat{N}=\sum _i{n}_i.$$

(1)

Introducing the Fourier representation of the projector,

$$\delta (\widehat{N}-N_b)={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }d\theta e^{i\theta (\widehat{N}-N_b)},$$

(2)

the canonical partition function can be written as

$$Z_{N_b}={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }d\theta e^{-i\theta N_b}Z\left(\theta \right),$$

(3)

where $Z\left(\theta \right)$ is a grand-canonical partition function with an imaginary chemical potential $\mu =i\theta$, for $\theta$ to be dimensionless. As we will see in next Section, the analysis is performed in the $\theta$-ensemble, where the boson number is not fixed but determined dynamically via

$$N_b\left(\theta \right)={\displaystyle \frac{\partial lnZ\left(\theta \right)}{\partial \left(i\theta \right)}}.$$

(4)

2.1. Particle-Channel Formulation and Large-N Effective Action

To describe the Bose–Fermi crossover, the repulsive Bose-Hubbard interaction must be decoupled in the particle channel, since the formation of composite fermions cannot be captured by a simple mean-field decoupling. We introduce a large-N generalization by extending the bosons to N flavors $\alpha =1,\dots ,N$, which provides a controlled saddle-point limit. The momentum version of Hamiltonian is [24,25]

$$H=\sum _{\mathbf{k},\alpha }{ϵ}_{\mathbf{k}}{b}_{\mathbf{k}\alpha }^{†}{b}_{\mathbf{k}\alpha }+{\displaystyle \frac{U}{2N}}\sum _i\sum _{\alpha ,\beta }{b}_{i\alpha }^{†}{b}_{i\alpha }^{†}{b}_{i\beta }{b}_{i\beta },$$

(5)

where $U>0$ is the repulsive interaction. The factor $1/N$ ensures that the interaction energy scales extensively in the large-N limit. The Euclidean action with imaginary chemical potential then takes the form

$$S={\int }_0^{\beta }d\tau \left[\sum _{\mathbf{k},\alpha }{b}_{\mathbf{k}\alpha }^{†}({\partial }_{\tau }+{ϵ}_{\mathbf{k}}-i\theta )b_{\mathbf{k}\alpha }+{\displaystyle \frac{U}{2N}}\sum _{i,\alpha ,\beta }{b}_{i\alpha }^{†}{b}_{i\alpha }^{†}{b}_{i\beta }{b}_{i\beta }\right].$$

(6)

The quartic interaction is decoupled via a Hubbard–Stratonovich transformation in the particle channel by introducing a complex bosonic field ${\Phi }_i\left(\tau \right)$ (which will become the fermionic order parameter):

$$exp\left[-{\displaystyle \frac{U}{2N}}\int d\tau b^{†}{b}^{†}bb\right]=\int \mathcal{D}\Phi \mathcal{D}{\Phi }^{*}exp\left[-\int d\tau \left({\displaystyle \frac{{N|\Phi |}^2}{2U}}+{\Phi }^{*}bb+\Phi b^{†}{b}^{†}\right)\right].$$

(7)

After this transformation, the bosonic action becomes quadratic in the original boson fields. Introducing the Nambu spinor for bosons,

$${\mathbf{Y}}_{\mathbf{k}\alpha }=\left(\begin{array}{c}{b}_{\mathbf{k}\alpha }\\ b_{-\mathbf{k}\alpha }^{†}\end{array}\right),$$

(8)

the inverse bosonic Green’s function is

$${\mathcal{G}}^{-1}(i{\omega }_n,\mathbf{k})=\left(\begin{array}{cc}(i{\omega }_n+i\theta )-{ϵ}_{\mathbf{k}}& -\Phi \\ {\Phi }^{*}& -(i{\omega }_n+i\theta )-{ϵ}_{\mathbf{k}}\end{array}\right),$$

(9)

where ${\omega }_n=2\pi n/\beta$ are the bosonic Matsubara frequencies. Its determinant yields the quasiparticle spectrum:

$$E_{\mathbf{k}}=\sqrt{{ϵ}_{\mathbf{k}}^2+{|\Phi |}^2},$$

(10)

which governs the Bose–Fermi crossover physics. The physical excitation spectrum is obtained from the poles of the Green’s function. Since the bosonic action is quadratic, the fields can be integrated out and the effective action is

$$S_{\mathrm{eff}}[\Phi ,\theta ]={\int }_0^{\beta }d\tau {\displaystyle \frac{{N|\Phi |}^2}{2U}}+N\mathrm{Tr}ln{\mathcal{G}}^{-1}[\Phi ,\theta ],$$

(11)

where the trace is taken over lattice sites, imaginary time (or Matsubara frequencies), Nambu indices, and internal degrees of freedom. We note the + sign in front of the trace, which arises from the bosonic functional integral (determinant in denominator).

The static saddle-point configuration is ${\Phi }_i\left(\tau \right)=\Phi$ and we may say that it is interpreted as an auxiliary field encoding local two-particle correlations $\Phi \sim \langle b_i{b}_i\rangle$ and also it does not correspond to a conventional order parameter for repulsive bosons. In contrast to the attractive fermionic case, where the analogous field describes Cooper pairing, here $\Phi$ measures the tendency of bosons to occupy the same site. Repulsive interactions suppress such configurations, so that the fermionization regime is characterized by a reduction of $\Phi$. Essentially, we could argue that $\Phi$ serves as an indicator of the degree of bosonic coherence versus fermion-like avoidance, rather than a true symmetry-breaking order parameter.

The trace over Nambu space yields

$$det{\mathcal{G}}^{-1}(i{\omega }_n,\mathbf{k})=-{(i{\omega }_n+i\theta )}^2+\left({ϵ}_{\mathbf{k}}^2+{|\Phi |}^2\right),$$

(12)

The functional trace can therefore be written as

$$\mathrm{Tr}ln{\mathcal{G}}^{-1}=\sum _{\mathbf{k}}\sum _{n}ln\left[-{(i{\omega }_n+i\theta )}^2+E_{\mathbf{k}}^2\right].$$

(13)

Performing the Matsubara frequency sum using contour-integration gives

$$\sum _{n}ln\left[-{(i{\omega }_n+i\theta )}^2+E_{\mathbf{k}}^2\right]=\beta E_{\mathbf{k}}+ln\left(1-e^{-\beta (E_{\mathbf{k}}+i\theta )}\right)+ln\left(1-e^{-\beta (E_{\mathbf{k}}-i\theta )}\right)+\mathrm{constant}.$$

(14)

The constant (independent of $\Phi$ and $\theta$) can be absorbed into the normalization of the partition function. Substituting this result back into the effective action yields

$$S_{\mathrm{eff}}[\Phi ,\theta ]=N\left[{\displaystyle \frac{{|\Phi |}^2}{2U}}+\sum _{\mathbf{k}}\left(E_{\mathbf{k}}+\left[ln\left(1-e^{-\beta (E_{\mathbf{k}}+i\theta )}\right)+ln\left(1-e^{-\beta (E_{\mathbf{k}}-i\theta )}\right)\right]\right)\right].$$

(15)

In the thermodynamic limit, the lattice momentum sum is replaced by an integral over the first Brillouin zone,

$$\sum _{\mathbf{k}}⟶{\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}.$$

(16)

The effective action density then becomes

$$S_{\mathrm{eff}}=N\left[{\displaystyle \frac{{|\Phi |}^2}{2U}}+{\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}\left(E_{\mathbf{k}}+\left[ln\left(1-e^{-\beta (E_{\mathbf{k}}+i\theta )}\right)+ln\left(1-e^{-\beta (E_{\mathbf{k}}-i\theta )}\right)\right]\right)\right].$$

(17)

3. Gap Equations for the Repulsive Bose-Hubbard Model

The first gap equation of the model for the order parameter $\Phi$ is

$${\displaystyle \frac{\partial S_{\mathrm{eff}}}{\partial {\Phi }^{*}}}=0$$

(18)

and yields the equation

$${\displaystyle \frac{1}{U}}={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}{\displaystyle \frac{1}{E_{\mathbf{k}}}}\left[1+f_B(E_{\mathbf{k}}+i\theta )+f_B(E_{\mathbf{k}}-i\theta )\right],$$

(19)

where $f_B\left(x\right)={\displaystyle \frac{1}{e^{\beta x}-1}}$ is the Bose-Einstein distribution function. Note the plus signs and the Bose-Einstein statistics, which differ from the fermionic case. At zero temperature, $f_B\left(x\right)\to 0$, so this reduces to

$${\displaystyle \frac{1}{U}}={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}{\displaystyle \frac{1}{E_{\mathbf{k}}}}.$$

(20)

For $\theta =0$, we obtain the standard gap equation for the repulsive Bose-Hubbard model

$${\displaystyle \frac{1}{U}}={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}{\displaystyle \frac{1+2f_B\left(E_{\mathbf{k}}\right)}{E_{\mathbf{k}}}}.$$

(21)

The canonical constraint is enforced by the saddle-point condition with respect to $\theta$:

$${\displaystyle \frac{\partial S_{\mathrm{eff}}}{\partial \theta }}=0.$$

(22)

This yields the number equation 1

$$N_b={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}\left[f_B(E_{\mathbf{k}}-i\theta )-f_B(E_{\mathbf{k}}+i\theta )\right],$$

(23)

where $N_b$ is the total boson number. In the bosonic case it quantifies particle–hole asymmetry induced by the imaginary chemical potential and reflects fermionization-like regime, or in other words how much the system prefers particle-like versus hole-like excitations, analogous to the fermionic case. For $D=d-1=1$ we will find that the fermionization-like window is where this asymmetry is maximal, corresponding to the point where bosons most closely mimic fermionic behaviour, since fermions have a natural particle-hole asymmetry due to the Pauli principle.

• Fermionization-like window and spectral redistribution.

Studying the phase structure of our theory in next section, we will find that at the special angles $\varphi =\pi /3$ and $\varphi =5\pi /3$, the repulsive Bose–Hubbard model at imaginary chemical potential enters a regime of maximal spectral redistribution. Here the competition between repulsive interactions and phase frustration suppresses Bose condensation and enhances entropy. This results in an effective fermion-like occupation profile, characterized by a flattening of the energy distribution and a reduced tendency for bosons to macroscopically occupy low-energy states. Although there is no true Pauli principle in our case, the thermal kernel induces an emergent exclusion-like behaviour that mimics fermionic statistics at the level of occupation numbers. The two angles correspond to opposite occupations: $\varphi =\pi /3$ enhances particle-like excitations, whereas $\varphi =5\pi /3$ produces a sign-reversed, hole-like response, reflecting a dual structure controlled by the imaginary chemical potential. This situation is analogous to Tonks-Girardeau (TG) gas, where bosons "fermionize" because they cannot pass through each other when they reach the limit of infinite repulsion energy $U\to \infty$. So we may say that our model suggests a different route to the same destination: while TG bosons require infinite energy to fermionize, our system achieves this state at finite interaction, driven by the ’statistical pressure’ of the imaginary chemical potential. On the other hand, we may say that the point $\varphi =\pi$ is not a maximally bosonic or fermionic configuration, but rather a (particle-like)–(hole-like) symmetric point where the number equation vanishes identically.

The gap equation may be split into two parts,

$$I_1={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}{\displaystyle \frac{1}{E_{\mathbf{k}}}}$$

(24)

and

$$I_2={\int }_{\mathrm{BZ}}{\displaystyle \frac{d^{d-1}k}{{\left(2\pi \right)}^{d-1}}}{\displaystyle \frac{1}{E_{\mathbf{k}}}}\left[f_B(E_{\mathbf{k}}+i\theta )+f_B(E_{\mathbf{k}}-i\theta )\right].$$

(25)

Note the plus sign in front of $I_2$ (compared to the minus sign in the fermionic case). For 1D chains at finite temperature, we take simply

$$I_1={\int }_{\mathrm{BZ}}{\displaystyle \frac{dk}{2\pi }}{\displaystyle \frac{1}{E_{\mathbf{k}}}}$$

(26)

and

$$I_2={\int }_{\mathrm{BZ}}{\displaystyle \frac{dk}{2\pi }}{\displaystyle \frac{1}{E_{\mathbf{k}}}}\left[f_B(E_{\mathbf{k}}+i\theta )+f_B(E_{\mathbf{k}}-i\theta )\right].$$

(27)

If we define the vacuum (critical) coupling by

$${\displaystyle \frac{1}{U_c}}\equiv I_1(\Phi =0)={\int }_{\mathrm{BZ}}{\displaystyle \frac{dk}{2\pi }}{\displaystyle \frac{1}{|{ϵ}_k|}},$$

(28)

and subtract this value, the renormalized gap equation becomes

$${\displaystyle \frac{1}{U}}-{\displaystyle \frac{1}{U_c}}={\int }_{\mathrm{BZ}}{\displaystyle \frac{dk}{2\pi }}\left({\displaystyle \frac{1}{\sqrt{{ϵ}_k^2+{|\Phi |}^2}}}-{\displaystyle \frac{1}{|{ϵ}_k|}}\right)+I_2(\Phi ,\theta ).$$

(29)

The sign of ${\displaystyle \frac{1}{U}}-{\displaystyle \frac{1}{U_c}}$ controls the crossover: positive corresponds to the Bose regime (weak repulsion), negative to the Fermi regime (strong repulsion, fermionization), and ${\displaystyle \frac{1}{U}}={\displaystyle \frac{1}{U_c}}$ marks the crossover (fermionization) point.

4. Phase Structure: Bose, Fermionization, and Crossover Regimes

4.1. General Setup for Bosons

Using the identity for Bose-Einstein distributions [26],

$$f_B(E_k+i\theta )+f_B(E_k-i\theta )=-1-{\displaystyle \frac{sinh\left(\beta E_k\right)}{cosh\left(\beta E_k\right)-cos\left(\beta \theta \right)}},$$

(30)

and $\varphi \equiv \beta \theta$ as a dimensionless variable, the bosonic gap equation becomes:

$${\delta }_u\equiv {\displaystyle \frac{1}{U}}-{\displaystyle \frac{1}{U_c}}={\int }_{\mathrm{BZ}}{\displaystyle \frac{dk}{2\pi }}\left[{\displaystyle \frac{sinh\left(\beta E_k\right)}{E_k\left[cosh\left(\beta E_k\right)-cos\varphi \right]}}-{\displaystyle \frac{1}{|{ϵ}_k|}}\right].$$

(31)

Note the minus sign in the denominator ($cosh\left(\beta E_k\right)-cos\varphi$) compared to the fermionic case. For a lattice dispersion ${ϵ}_k=-2tcosk$, the band minimum occurs at $k=0$ (or $k=\pi$ depending on the sign). Near the band minimum $k=0$, we expand:

$${ϵ}_k\approx -2t\left(1-{\displaystyle \frac{k^2}{2}}\right)=-2t+t{k}^2,$$

(32)

so the low-energy dispersion is quadratic:

$${ϵ}_k\approx \mathrm{constant}+{\displaystyle \frac{k^2}{2m^{*}}},$$

(33)

where $m^{*}=1/\left(2t\right)$ is the effective mass. For simplicity, we measure energies from the band bottom, so ${ϵ}_k\approx t{k}^2$. If we split the integral into low-energy contributions near $k=0$ (within a cutoff $\left|k\right|<\Lambda$) and the remainder $C_{\mathrm{UV}}$ (higher momentum part inside BZ), we obtain the low-energy integral:

$${\delta }_u=\underset{{\displaystyle {\mathcal{F}}_{\mathrm{low}}(\Phi ,T,\varphi ;\Lambda )}}{\underbrace{{\int }_0^{\Lambda }{\displaystyle \frac{dk}{\pi }}\left[{\displaystyle \frac{sinh\left(\beta E_k\right)}{E_k\left[cosh\left(\beta E_k\right)-cos\varphi \right]}}-{\displaystyle \frac{1}{t{k}^2}}\right]}}+C_{\mathrm{UV}}(\Lambda ).$$

(34)

The critical coupling is defined as:

$${\displaystyle \frac{1}{U_c}}={\int }_0^{\Lambda }{\displaystyle \frac{dk}{\pi }}{\displaystyle \frac{1}{t{k}^2}}.$$

(35)

Since $C_{\mathrm{UV}}(\Lambda )$ comes from the rest of the Brillouin zone, it is independent of $\Phi$ and T in the low-energy limit ($\Lambda \ll 1$, low T and $\Phi$). Essentially $C_{\mathrm{UV}}$ encodes the contribution of high energy modes within the lattice band where the phase dependence and the structure of thermal windows are governed almost entirely by the low energy part ${\mathcal{F}}_{\mathrm{low}}$. So we may absorb it into a redefinition of the parameter ${\delta }_u$.

4.2. Bosonic Thermal Kernel

The bosonic thermal kernel is defined as:

$$g_B(x,\varphi )\equiv {\displaystyle \frac{sinhx}{coshx-cos\varphi }},x=\beta E_k.$$

(36)

Using the key identity

$$g_B(x,\varphi )-1={\displaystyle \frac{e^{-x}-cos\varphi }{coshx-cos\varphi }},$$

(37)

we obtain:

$$g_B(x,\varphi )⋛1⟺e^{-x}⋛cos\varphi .$$

(38)

Since $e^{-x}\le 1$ for $x\ge 0$, this inequality has interesting structure:

• If $cos\varphi =1$ ($\varphi =0,2\pi$), then $g_B(x,\varphi )>1$ for all $x>0$
• If $cos\varphi <1$, there exists a threshold $x_{*}=-ln(cos\varphi )$ where $g_B=1$

At the special angle $\varphi =\pi /3$ (and similarly at $5\pi /3$), where $cos\varphi =1/2$, the kernel becomes:

$$g_B(x,\pi /3)={\displaystyle \frac{sinhx}{coshx-\frac{1}{2}}},$$

(39)

and satisfies

$$g_B(x,\pi /3)-1={\displaystyle \frac{e^{-x}-\frac{1}{2}}{coshx-\frac{1}{2}}}.$$

(40)

This expression changes sign at $x_{*}=ln2$, exactly as in the fermionic case. Consequently, the integrand of the bosonic gap equation exhibits a nontrivial structure: low-energy modes ($x<ln2$) contribute with $g_B<1$, while high-energy modes ($x>ln2$) contribute with $g_B>1$. At the critical point ${\delta }_u=0$, these contributions balance exactly after inclusion of the high-energy term $C_{\mathrm{UV}}$. Therefore, the angles $\varphi =\pi /3$ and $5\pi /3$ mark the points where the bosonic thermal kernel produces a sign-changing contribution across energy scales. This leads to a particularly sensitive balance in the gap equation at the critical coupling, corresponding to the fermionization-like window where bosons are under an effective exclusion principle at the level of thermodynamics, in order to minimize free energy via entropy gain and interaction avoidance:

$$\begin{array}{|c|}\hline \varphi \in \left[{\displaystyle \frac{\pi }{3}},{\displaystyle \frac{5\pi }{3}}\right]\\ \hline\end{array},$$

(41)

Notably, $ln2$ appears as entropy values of this order in [27] as a similar result. In the low-energy limit, the thermal kernel can be approximated as

$$g(ϵ,\varphi )\approx {\displaystyle \frac{1}{(1-cos\varphi )+\frac{{\left(\beta ϵ\right)}^2}{2}}},$$

(42)

so that at small $ϵ$,

$$g(ϵ,\varphi )-1\approx {\displaystyle \frac{cos\varphi }{1-cos\varphi }}.$$

(43)

This shows that the sign of the low-energy contribution is determined by $cos\varphi$. Since the kernel is peaked at low energies, the integral $\int dϵ\rho \left(ϵ\right)\left[g\right(ϵ,\varphi )-1]$ is dominated by this region, implying that the sign of the integral is approximately controlled by the sign of $g(ϵ,\varphi )-1$. Deviations from this simple picture arise from finite-energy contributions, which shift the precise location of the phase boundaries away from the condition $cos\varphi =0$.

• Bosonic regimes.

For the repulsive Bose–Hubbard model at imaginary chemical potential, the regimes are controlled by the parameter ${\delta }_u=1/U-1/U_c$.

For ${\delta }_u<0$ (strong repulsion), we see that clearly bosons exhibit a fermionization tendency: the system actively suppresses low-energy states, pushing spectral weight toward higher energies. This is not a true transition to fermionic statistics, but rather an interaction-driven flattening of the occupation profile.

On the other hand for ${\delta }_u>0$ (weak repulsion), the system remains in a Bose-like regime, where occupation of low-energy states is enhanced and the phase twist induces only mild distortions of the distribution.

At the critical point ${\delta }_u=0$, the system undergoes a crossover between these different regimes. The effect is strongest at the special angles $\varphi =\pi /3,5\pi /3$, where the thermal kernel defines the fermionization window.

However, it is important to say that we should be careful to distinguish this large-N formulation, based on a particle-channel decoupling, from the classic Mott-insulator-to-superfluid transition, which is associated with the condensation of the single-particle field $\langle b\rangle$. Instead, the present approach focuses on a distinct crossover driven by repulsive interactions and imaginary chemical potential, leading to an effective fermionization of bosonic excitations.

In this sense, the regimes identified here should be interpreted as “Bose-like” and “fermionized” regimes, rather than true superfluid or Mott insulating phases. The strong-repulsion regime (${\delta }_u<0$) shares qualitative features with Mott physics, such as suppression of low-energy occupation, but does not exhibit a true gap associated with integer filling.

At the particle–hole-like symmetric point $\varphi =\pi$, the kernel reduces to

$${\displaystyle \frac{1}{cosh\left(\beta E\right)+1}},$$

which is identical to the fermionic thermal kernel. This establishes a precise thermodynamic correspondence, although the microscopic degrees of freedom remain bosonic. Unlike the BCS–BEC crossover, where composite fermions create bound pairs that behave as bosons, our mechanism does not rely on pairing. Instead, the imaginary chemical potential modifies the thermal kernel, leading to fermion-like occupation constraints for elementary bosons without changing their underlying statistics.

On the other hand, as we have already mentioned, if we interpret $C_{\mathrm{UV}}$ as a renormalization shift of the interaction parameter, we may write

$${\delta }_u^{\mathrm{eff}}={\delta }_u-C_{\mathrm{UV}}={\mathcal{F}}_{\mathrm{low}}(\Phi ,T,\varphi ),$$

(44)

where ${\mathcal{F}}_{\mathrm{low}}$ encodes the low-energy contribution weighted by the thermal kernel $g(ϵ,\varphi )$.

• ${\delta }_u-C_{\mathrm{UV}}>0$: This requires ${\mathcal{F}}_{\mathrm{low}}>0$, which occurs when low-energy modes are effectively enhanced. This corresponds to $g(ϵ,\varphi )>1$ in the infrared, leading to bosonic behaviour.
• ${\delta }_u-C_{\mathrm{UV}}<0$: This requires ${\mathcal{F}}_{\mathrm{low}}<0$, which arises when low-energy modes are suppressed relative to higher-energy contributions. This corresponds to $g(ϵ,\varphi )<1$ in the infrared, inducing an effective exclusion principle that leads to fermionization Figure 1.

Importantly, fermionization does not occur simply because $g<1$, but because $g(ϵ,\varphi )$ creates a competition between low energy suppression and high energy enhancement. This competition is realized inside the thermal window $\varphi \in [\pi /3,5\pi /3]$, where the kernel changes its relative weighting across energy scales. Outside this window, $g(ϵ,\varphi )>1$ for all energies, and the system reduces to a conventional bosonic regime without fermionization effects. We see the phase structure of the model in Figure 1. Notably, we see that lowering $T/t$ means thermal fluctuations are reduced. Inside the thermal window, the low-energy modes dominate the gap equation. These low-energy modes push ${\delta }_u$ negative, favoring $U>U_c$ (strong repulsion), which is the fermionization-like regime. In contrast, at higher temperatures, thermal excitations push ${\delta }_u$ positive, favoring Bose-like behaviour.

Interestingly, the condition $e^{-x}={\displaystyle \frac{1}{2}}$ plays a role analogous to the appearance of $ln2$ in the $C{P}^{N-1}$ model at imaginary chemical potential [9,10], but in our case, while in the continuous model this configuration appears as the free energy density, in our formulation it emerges through the structure of the thermal kernel.

4.3. Extrema of $N_b\left(\varphi \right)$ at $\varphi =\pi /3,5\pi /3$

We analyze the extrema and curvature of the boson number $N_b\left(\varphi \right)$ for the repulsive Bose–Hubbard model at imaginary chemical potential. Using the Bose distribution identity,

$$f_B(E-i\varphi )-f_B(E+i\varphi )={\displaystyle \frac{sin\varphi }{cosh\left(\beta E\right)-cos\varphi }},$$

(45)

the boson number at $\Phi =0$ reads

$$N_b(\varphi ,T)={\int }_{-2t}^{2t}dϵ\rho \left(ϵ\right){\displaystyle \frac{sin\varphi }{cosh\left(\beta \right|ϵ\left|\right)-cos\varphi }}.$$

(46)

Assuming band symmetry,

$$\rho \left(ϵ\right)=\rho (-ϵ),$$

(47)

the integrand is an even function of $ϵ$. If we define

$$I_B\left(\varphi \right)={\int }_{-2t}^{2t}dϵ\rho \left(ϵ\right){\displaystyle \frac{1}{cosh\left(\beta \right|ϵ\left|\right)-cos\varphi }},$$

(48)

and

$$J_B\left(\varphi \right)={\int }_{-2t}^{2t}dϵ\rho \left(ϵ\right){\displaystyle \frac{1}{{\left(cosh\left(\beta \right|ϵ\left|\right)-cos\varphi \right)}^2}}$$

(49)

then the first derivative is

$${\displaystyle \frac{d{N}_b}{d\varphi }}=cos\varphi I_B\left(\varphi \right)+{sin}^2\varphi J_B\left(\varphi \right).$$

(50)

Stationary points are determined by

$$cos\varphi I_B\left(\varphi \right)+{sin}^2\varphi J_B\left(\varphi \right)=0.$$

(51)

At the special angles $\varphi ={\displaystyle \frac{\pi }{3}},{\displaystyle \frac{5\pi }{3}}$, we have $cos\varphi ={\displaystyle \frac{1}{2}},{sin}^2\varphi ={\displaystyle \frac{3}{4}}$ so Equation (51) yields

$${\displaystyle \frac{1}{2}}{I}_B+{\displaystyle \frac{3}{4}}{J}_B=0⟹J_B=-{\displaystyle \frac{2}{3}}{I}_B.$$

(52)

Differentiating Equation (50) gives

$${\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}=-sin\varphi I_B+cos\varphi I_B^{\prime }+2sin\varphi cos\varphi J_B+{sin}^2\varphi J_B^{\prime }.$$

(53)

and using the identity

$$I_B^{\prime }\left(\varphi \right)=-sin\varphi J_B\left(\varphi \right),$$

(54)

we obtain

$${\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}=-sin\varphi I_B+sin\varphi cos\varphi J_B+{sin}^2\varphi J_B^{\prime }.$$

(55)

Using $sin\varphi ={\displaystyle \frac{\sqrt{3}}{2}},cos\varphi ={\displaystyle \frac{1}{2}}$ Equation (55) gives

$${\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}{|}_{\varphi =\pi /3}=-{\displaystyle \frac{\sqrt{3}}{2}}{I}_B+{\displaystyle \frac{\sqrt{3}}{4}}{J}_B+{\displaystyle \frac{3}{4}}{J}_B^{\prime }.$$

(56)

From Equation (52), we obtain

$$\begin{array}{|c|}\hline {\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}{|}_{\varphi =\pi /3}=-{\displaystyle \frac{2\sqrt{3}}{3}}{I}_B+{\displaystyle \frac{3}{4}}{J}_B^{\prime }.\\ \hline\end{array}$$

(57)

The derivative $J_B^{\prime }\left(\varphi \right)$ is given by

$$J_B^{\prime }\left(\varphi \right)=-2sin\varphi {\int }_{-2t}^{2t}dϵ\rho \left(ϵ\right){\displaystyle \frac{1}{{\left(cosh\left(\beta \right|ϵ\left|\right)-cos\varphi \right)}^3}}.$$

(58)

At $\varphi =\pi /3$, $sin\varphi >0$, and the integrand is positive, so

$$J_B^{\prime }(\pi /3)<0$$

(59)

and since $I_B>0\to {\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}{|}_{\varphi =\pi /3}<0$.

At $\varphi =5\pi /3$ we have $sin\varphi =-{\displaystyle \frac{\sqrt{3}}{2}},cos\varphi ={\displaystyle \frac{1}{2}}$ leading to

$$\begin{array}{|c|}\hline {\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}{|}_{\varphi =5\pi /3}={\displaystyle \frac{2\sqrt{3}}{3}}{I}_B+{\displaystyle \frac{3}{4}}{J}_B^{\prime }.\\ \hline\end{array}$$

(60)

In this case, since $sin\varphi <0$, one finds

$$J_B^{\prime }(5\pi /3)>0$$

(61)

and since $I_B>0\to {\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}}{|}_{\varphi =5\pi /3}>0$. So:

• At $\varphi =\pi /3$: ${\displaystyle \frac{d{N}_b}{d\varphi }=0}$ and ${\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}<0}$, therefore $N_b$ has a local maximum. This corresponds to the maximal particle-hole asymmetry with particle-like excitations dominating in fermionization regime.
• At $\varphi =5\pi /3$: ${\displaystyle \frac{d{N}_b}{d\varphi }=0}$ and ${\displaystyle \frac{d^2{N}_b}{d{\varphi }^2}>0}$, therefore $N_b$ has a local minimum. This corresponds to the maximal (negative) particle-hole asymmetry with hole-like excitations dominating in fermionization regime.
  This situation is corresponding to the fermionic case. At the edges of thermal window we have particle-like (maximum) at $2\pi /3$ and hole-like (minimum) at $4\pi /3$. On the other hand for bosons, at the edges of thermal window we have particle-like (maximum) at $\pi /3$ and hole-like (minimum) at $5\pi /3$.

5. Notes on Attractive Fermi-Hubbard Model at Large- $\mathit{N}$ Formulation

In this section it is important to recall the physics of the main results from the attractive Fermi–Hubbard model [28] with an imaginary chemical potential $\mu =i\theta$, which acts as a phase twist in imaginary time and modifies the Matsubara frequencies as $i{\omega }_n\to i{\omega }_n+i\theta$.

The inverse fermionic Green’s function in Nambu space takes the form

$${\mathcal{G}}^{-1}(i{\omega }_n,\mathbf{k})=\left(\begin{array}{cc}(i{\omega }_n+i\theta )-{ϵ}_{\mathbf{k}}& \Delta \\ {\Delta }^{*}& (i{\omega }_n+i\theta )+{ϵ}_{\mathbf{k}}\end{array}\right),$$

(62)

which preserves particle–hole symmetry. The quasiparticle spectrum is determined by

$$E_{\mathbf{k}}=\sqrt{{ϵ}_{\mathbf{k}}^2+{|\Delta |}^2},$$

(63)

and the determinant of the Green’s function is

$$det{\mathcal{G}}^{-1}={(i{\omega }_n+i\theta )}^2-E_{\mathbf{k}}^2.$$

(64)

The thermodynamic potential is given by

$$\Omega ={\displaystyle \frac{{|\Delta |}^2}{2U}}-\int {\displaystyle \frac{dk}{2\pi }}{E}_k-{\displaystyle \frac{1}{\beta }}\int {\displaystyle \frac{dk}{2\pi }}ln\left[1+2e^{-\beta E_k}cos\varphi +e^{-2\beta E_k}\right],$$

(65)

with $\varphi =\beta \theta$. The gap equation becomes

$${\displaystyle \frac{1}{U}}=\int {\displaystyle \frac{dk}{2\pi }}{\displaystyle \frac{1}{E_k}}\left[1-f(E_k+i\theta )-f(E_k-i\theta )\right],$$

(66)

where $f\left(x\right)=1/(e^{\beta x}+1)$ is the Fermi–Dirac distribution. The number equation reads

$$N_f=\int {\displaystyle \frac{dk}{2\pi }}\left[f(E_k-i\theta )-f(E_k+i\theta )\right],$$

(67)

which measures the imbalance between particle-like and hole-like excitations induced by the imaginary chemical potential. The thermodynamics is governed by the fermionic thermal kernel

$$g_F(x,\varphi )={\displaystyle \frac{sinhx}{coshx+cos\varphi }},x=\beta E_k.$$

(68)

The kernel satisfies

$$g_F(x,\varphi )-1=-{\displaystyle \frac{e^{-x}+cos\varphi }{coshx+cos\varphi }},$$

(69)

which determines the sign structure of contributions to the gap equation. The special angles $\varphi =2\pi /3$ and $\varphi =4\pi /3$ (where $cos\varphi =-1/2$) define a universal thermal window in which the kernel changes sign at $x_{*}=ln2$.

Inside this thermal window, spectral weight is redistributed between low- and high-energy modes, leading to enhanced pairing correlations and driving the BCS–BEC crossover. The crossover is controlled by the parameter

$${\delta }_u={\displaystyle \frac{1}{U}}-{\displaystyle \frac{1}{U_c}},$$

(70)

where ${\delta }_u<0$ corresponds to the BEC regime (strong attraction) and ${\delta }_u>0$ to the BCS regime (weak attraction). The imaginary chemical potential acts as a statistical regulator, continuously tuning the system between different pairing regimes by modifying the entropy contribution to the free energy. At the critical point ${\delta }_u=0$, the system exhibits universal behaviour characterized by a thermal window which plays a central role in determining the phase structure.

6. Mapping Between Attractive Fermi and Repulsive Bose Hubbard Models at Imaginary Chemical Potential

6.1. Formal Mapping via Matsubara Frequencies

The partition function for fermions with imaginary chemical potential $i\theta$ is

$${\mathcal{Z}}_F\left(\theta \right)=\prod _{\mathbf{k},n}\left[{(i{\omega }_n^F+i\theta )}^2-E_{\mathbf{k}}^2\right],{\omega }_n^F={\displaystyle \frac{(2n+1)\pi }{\beta }}.$$

(71)

For bosons, the Matsubara frequencies are ${\omega }_n^B=2n\pi /\beta$. Using the relation ${\omega }_n^B={\omega }_n^F+\pi /\beta$, we obtain

$$i{\omega }_n^B+i\theta =i{\omega }_n^F+i(\theta +\pi ).$$

(72)

Consequently,

$$\begin{array}{|c|}\hline {\mathcal{Z}}_B\left(\theta \right)={\mathcal{Z}}_F(\theta +\pi )\\ \hline\end{array}.$$

(73)

Thus, the repulsive Bose-Hubbard model at imaginary chemical potential $\theta$ is exactly dual to the attractive Fermi-Hubbard model at $\theta +\pi$. Within the fermionization-like window, the free energy $F=U-TS$ is minimized by interaction avoidance since the phase twist increases effective density pressure. Then bosons are rearranged into higher energy levels, mimicking Pauli exclusion to avoid the repulsive interaction cost U. Also entropic gain makes the fermionized distribution occupies a larger spectral range than a standard Bose-Einstein condensate, lowering F.

6.2. Mapping Between Fermionic and Bosonic Thermal Kernels

The fermionic thermal kernel derived from the attractive Hubbard model is

$$g_F(x,\varphi )={\displaystyle \frac{sinhx}{coshx+cos\varphi }},x=\beta E_k,\varphi =\beta \theta .$$

(74)

For the repulsive Bose-Hubbard model, the corresponding bosonic thermal kernel is

$$g_B(x,\varphi )={\displaystyle \frac{sinhx}{coshx-cos\varphi }}.$$

(75)

The two kernels are related by a simple shift of the imaginary chemical potential:

$$\begin{array}{|c|}\hline g_B(x,\varphi )=g_F(x,\varphi +\pi )\\ \hline\end{array}.$$

(76)

Indeed, using $cos(\varphi +\pi )=-cos\varphi$, we obtain

$$g_F(x,\varphi +\pi )={\displaystyle \frac{sinhx}{coshx+cos(\varphi +\pi )}}={\displaystyle \frac{sinhx}{coshx-cos\varphi }}=g_B(x,\varphi ).$$

(77)

This mapping reflects the underlying duality between the attractive Fermi and repulsive Bose Hubbard models at imaginary chemical potential, with the shift $\varphi \to \varphi +\pi$ exchanging the roles of the two statistics.

7. Universal Statistical Transmutation Framework at Imaginary Chemical Potential

We consider a generic interacting lattice system at finite temperature with an imaginary chemical potential $\mu =i\theta$. The grand-canonical partition function can be written as

$$\begin{array}{|c|}\hline \mathcal{Z}\left(\varphi \right)={\displaystyle \prod _{\mathbf{k},n}}{\left[{(i{\omega }_n+i\theta )}^2-E_{\mathbf{k}}^2\right]}^{\sigma },\\ \hline\end{array}$$

(78)

where $\sigma =+1$ for fermions and $\sigma =-1$ for bosons. Taking the logarithm and performing the Matsubara sum yields the universal free energy density:

$$\begin{array}{|c|}\hline {\Omega }_{\sigma }\left(\varphi \right)={\displaystyle \frac{{\gamma }_{\sigma }^2}{2U}}-\int {\displaystyle \frac{d^{d}k}{{\left(2\pi \right)}^d}}\left[E_{\mathbf{k}}+{\displaystyle \frac{\sigma }{\beta }}ln\left(1+\sigma 2e^{-\beta E_{\mathbf{k}}}cos\varphi +e^{-2\beta E_{\mathbf{k}}}\right)\right].\\ \hline\end{array}$$

(79)

where ${\gamma }_{+1}\equiv \Delta$ and ${\gamma }_{-1}\equiv \Phi$ and this expression unifies fermionic and bosonic statistics through the parameter $\sigma$. Differentiating the free energy with respect to E yields the universal occupation kernel:

$$\begin{array}{|c|}\hline g_{\sigma }(x,\varphi )={\displaystyle \frac{sinhx}{coshx+\sigma cos\varphi }},x=\beta E.\\ \hline\end{array}$$

(80)

This kernel governs all thermodynamic observables.

• Fermions ($\sigma =+1$):

$$g_F(x,\varphi )={\displaystyle \frac{sinhx}{coshx+cos\varphi }}$$

(81)

• Bosons ($\sigma =-1$):

$$g_B(x,\varphi )={\displaystyle \frac{sinhx}{coshx-cos\varphi }}$$

(82)

and so the free energy densities obey

$$\begin{array}{|c|}\hline {\Omega }_B\left(\varphi \right)={\Omega }_F(\varphi +\pi ).\\ \hline\end{array}$$

(83)

This establishes a precise duality:

$$\begin{array}{c}\hfill \begin{array}{|c|}\hline \mathit{Bosonic}\mathit{thermal}\mathit{kernel}\mathit{at}\mathit{phase}\varphi \mathit{is}\mathit{equivalent}\mathit{to}\mathit{fermionic}\mathit{thermal}\mathit{kernel}\mathit{at}\mathit{phase}\varphi +\pi .\\ \hline\end{array}\end{array}$$

The kernel satisfies

$$g_{\sigma }(x,\varphi )-1=-{\displaystyle \frac{e^{-x}+\sigma cos\varphi }{coshx+\sigma cos\varphi }}.$$

(84)

The sign change occurs at

$$e^{-x_{*}}=-\sigma cos\varphi .$$

(85)

Thus:

• Fermions ($\sigma =+1$): critical angle $\varphi =2\pi /3$
• Bosons ($\sigma =-1$): critical angle $\varphi =\pi /3$

This reveals a universal ${\mathbb{Z}}_3$ structure:

$$\begin{array}{|c|}\hline {\varphi }_c^{\left(B\right)}={\varphi }_c^{\left(F\right)}+\pi .\\ \hline\end{array}$$

(86)

The order parameter $\Delta$ (fermions) or $\Phi$ (bosons) satisfies:

$$\begin{array}{|c|}\hline {\displaystyle \frac{1}{U}}-{\displaystyle \frac{1}{U_c}}=\int dϵ\rho \left(ϵ\right)\left[{\displaystyle \frac{g_{\sigma }(\beta E,\varphi )}{E}}-{\displaystyle \frac{1}{ϵ}}\right].\\ \hline\end{array}$$

(87)

This single equation describes the BCS–BEC crossover for fermions, the fermionization-like crossover for bosons and dimensional effects via $\rho \left(ϵ\right)$. Interestingly, after this map between the two theories, we may say that imaginary chemical potential which shifts Matsubara frequencies, modifies occupation statistics, and creates effective exclusion (fermionization). The free energy then becomes:

$$F=U-TS\left(\varphi \right),$$

(88)

where $\varphi$ controls the entropy functional.

On the other hand, if one wants to find the critical thermal window at unitarity then using the Landau theory of phase transitions by expanding the thermal potential near a critical point, then the phase boundary is determined by:

$$a_2(T,\varphi )=0$$

(89)

with

$$a_2={\delta }_u-\int dϵ\rho \left(ϵ\right){\displaystyle \frac{g_{\sigma }(\beta ϵ,\varphi )-1}{ϵ}}.$$

(90)

arising from analytical calculation in Appendix (A). It is shown in Table 1 the duality picture where BCS–BEC crossover and bosonic fermionization arise from the same underlying mechanism.

8. Discussion

We have established an exact duality between the attractive Fermi-Hubbard model and the repulsive Bose-Hubbard model at finite temperature in the presence of an imaginary chemical potential $\mu =i\theta$. Using the relation ${\omega }_n^B={\omega }_n^F+\pi /\beta$ between bosonic and fermionic Matsubara frequencies, we find ${\mathcal{Z}}_B\left(\theta \right)={\mathcal{Z}}_F(\theta +\pi )$, which maps the repulsive Bose-Hubbard model onto the attractive Fermi-Hubbard model shifted by $\pi$. Consequently, the thermal kernels satisfy $g_B(x,\varphi )=g_F(x,\varphi +\pi )$ with $\varphi =\beta \theta$, and the unified kernel $g_{\sigma }(x,\varphi )=sinhx/(coshx+\sigma cos\varphi )$ ($\sigma =+1$ for fermions, $\sigma =-1$ for bosons) governs both models. The condition $g_{\sigma }(x,\varphi )=1$ yields $x_{*}=ln2$ when $cos\varphi =-\sigma /2$, which interestingly defines the universal thermal window boundaries: $\varphi =2\pi /3,4\pi /3$ for fermions and $\varphi =\pi /3,5\pi /3$ for bosons. Inside these windows, spectral weight is redistributed between lower- and higher-energy modes.

For fermions inside the window, ${\delta }_u=1/U-1/U_c<0$ (strong attraction) favors BEC-like pairing, while ${\delta }_u>0$ (weak attraction) favors BCS-like behaviour. For bosons inside the fermionization window, ${\delta }_u<0$ (strong repulsion) favors Fermi-like spreading, while ${\delta }_u>0$ (weak repulsion) favors Bose-like behaviour. In both cases, lowering $T/t$ enhances the special behaviour inside the window, as the entropy term $-TS$ becomes less dominant and the interaction energy governs the free energy $F=U-TS$. The imaginary chemical potential thus acts as a temporal gauge field that continuously interpolates between Bose and Fermi statistics, with dimensionality entering only through the density of states $\rho \left(ϵ\right)$.

We think that our results offer a new window into the physics of fermionization and Bose-Fermi mapping at finite temperature since this duality unifies BCS-BEC crossover in fermions and fermionization in bosons within a single framework. Notably, like the fermionic case, where the extrema of fermion number $N_f\left(\varphi \right)$ are fixed by particle–hole symmetry, the extrema of boson number $N_b\left(\varphi \right)$ are also universal. This is a result which makes the curvature exhibit a dynamic selection of particle- or hole-like behaviour at the window boundaries which they correspond to fixed maxima/minima.

Although boson–fermion mappings at imaginary chemical potential have been explored in field-theoretic contexts [9,10], the interplay between imaginary chemical potential and fermionization in lattice bosonic and fermionic systems remains largely unexplored. Our findings imply that they could be extended in several new directions and present a remarkable progress in our understanding of duality identities, like an overall Fermi-Boson Hubbard model coupled to a gauge field or it would be interesting to examine our theory with experimental results arising from the appearance of synthetic gauge fields or Floquet engineering [29]. This could be an opened window to ultracold atoms like [30] where anyon statistical angle is the closest experimental analogue to our $\varphi$ from imaginary chemical potential. Also a Bose-Fermi Hubbard map for general d dimensions would be an interesting concept along with the fermionization-like windows as a generalization to the $1D$ results or examine if the edges of thermal windows are at different angles determined by the zeros or extrema of Clausen functions. Finally a possible generalization of this work to disorder Hubbard models like [31] could possibly give interesting results.