Imaginary-Interaction Fermi-Hubbard Models

• Imaginary-interaction Fermi-Hubbard models are non-Hermitian extensions that replace real on-site interactions with a purely imaginary coupling to explore quantum dissipation.
• They map many-body quantum evolution to Lindbladian dephasing, enabling efficient classical sampling via stochastic unitary channel approximations.
• The model provides practical insights for probing non-Hermitian phenomena, benchmarking quantum hardware, and designing advanced simulation algorithms under strict parameter constraints.

An imaginary-interaction Fermi-Hubbard model is a non-Hermitian extension of the traditional Fermi-Hubbard model, in which the on-site interaction term is replaced by a purely imaginary coupling. This variant is defined on bipartite lattices and is distinguished by a duality with dephasing Lindbladian evolution of free fermions. The resulting model enables a classically efficient sampling algorithm for simulating its time dynamics under particular parameter constraints, in contrast to the computational hardness generically associated with interacting fermion models. Such models provide an exact mapping between non-Hermitian many-body quantum evolution and open-system Lindblad dynamics, revealing unexpected algorithmic tractability and new pathways for probing non-Hermitian quantum phenomena (Santos, 18 Jan 2026).

1. Hamiltonian Formulation and Structure

The Fermi-Hubbard model on a bipartite graph $\mathcal G = \mathcal A \cup \mathcal B$ is formulated using spinless fermion operators $c_i, c_i^\dagger$ with canonical anticommutation relations $\{c_i, c_j^\dagger\} = \delta_{ij}$. The standard Hermitian Hamiltonian is

$$H_{\rm Hubb} = \sum_{i\in\mathcal A \atop j\in\mathcal B} h_{ij} (c_i^\dagger c_j + c_j^\dagger c_i) + U \sum_i n_{i\uparrow} n_{i\downarrow}, \qquad n_{i\sigma} = c_{i\sigma}^\dagger c_{i\sigma}.$$

The imaginary-interaction extension substitutes the real on-site interaction $U$ by a purely imaginary $i\gamma$: $$\mathcal H = \sum_{i,j \in \mathcal G \atop \sigma = \uparrow, \downarrow} h_{ij} a_{i,\sigma}^\dagger a_{j,\sigma} + i\gamma\sum_i n_{i\uparrow} n_{i\downarrow} - \frac{i\gamma}{2}\sum_{i,\sigma}n_{i\sigma},$$ where $a_{i,\downarrow} \equiv c_i$ and $a_{i,\uparrow} \equiv \bar c_i$ are two "copies" (interpreted as spin-up and spin-down) of spinless fermions. The hopping is restricted to inter-sublattice links due to the bipartite structure, which is crucial for mapping to the Lindbladian formalism. The coefficient $\gamma>0$ parameterizes the strength of the imaginary interaction (corresponding to a dephasing rate in the dual description).

2. Duality with Lindbladian Dephasing

The mapping between imaginary-interaction Fermi-Hubbard evolution and Lindbladian dephasing employs the thermofield vectorized representation of density matrices for free fermions. The Lindblad generator for dephasing is

$$\mathcal L(\rho) = -i[H_0, \rho] + \gamma \sum_j \left( n_j \rho n_j - \frac{1}{2} \{ n_j, \rho \} \right), \qquad H_0 = \sum_{ij}h_{ij}c_i^\dagger c_j.$$

By vectorizing $\rho$ and introducing "right-acting" fermions $\bar c_j$, the Lindbladian is represented as

$$\mathcal L = -i(H_0 - \bar H_0) + \gamma \sum_j (n_j \bar n_j - \tfrac{1}{2}(n_j + \bar n_j)),\qquad \bar H_0 = \sum_{ij} h_{ij} \bar c_i^\dagger \bar c_j.$$

After conjugating with the twist operator $$\mathcal U_{\mathcal A}(\pi) = \exp(i\pi\sum_{j\in \mathcal A} \bar n_j)$$ (flipping the sign of $\bar c_j$ on $\mathcal A$) and relabeling, the Lindbladian generator $-i\mathcal L$ becomes the non-Hermitian Hubbard Hamiltonian $\mathcal H$. Thus, time evolution under imaginary-interaction Hubbard models is isomorphic to Lindblad master equations for dephasing.

3. Classical Sampling Algorithm: Mixed-Unitary Channels

The duality enables an efficient classical simulation via stochastic sampling of unitary channels, as shown by Wang et al. (Wang et al., 9 Jan 2026). The discrete evolution of the Lindbladian $e^{\delta \mathcal L}$ is approximated as an average over random unitaries: $$e^{\delta \mathcal L} = \mathbb{E}_{\mathbf s}\left[ \mathcal U_{\delta, \mathbf s} \right] + O(\delta^2),$$ where $\mathbf s = (s_1, \dots, s_N)$ with $s_j = \pm 1$ (i.i.d., probability $1/2$), and

$$\mathcal U_{\delta, \mathbf s}(\rho) = e^{-iH_0\delta}\left( \prod_j e^{s_j i \sqrt{\gamma\delta} n_j} \right) \rho \left( \prod_j e^{-s_j i \sqrt{\gamma\delta} n_j} \right) e^{iH_0\delta}.$$

Iterating $R$ steps of size $\delta = t/R$ and averaging over $M$ Monte Carlo samples reconstructs the time evolution up to an error $O(t^2/R)$. The sampling is performed in the free-fermion Gaussian formalism, with polynomial cost per sample. High-level pseudocode for the method is:

Step	Description
1	Initialize $$\rho_0 = |\Psi_\uparrow\rangle\langle\Psi_\downarrow|$$
2a	For each sample $k=1,\dots,M$, draw random signs $s_j^{(r)}$ for each $T$rotter step and site
2b	Initialize $\rho^{(k)} \leftarrow \rho_0$
2c	For $r=1,\dots,R$, apply $\mathcal{U}_{\delta, \mathbf{s}^{(r)}}$ as above
2d	Extract amplitude $\Psi_{\mathbf n, \mathbf m}^{(k)}(t)$ by tracing against desired configuration
3	Report mean amplitude as estimator for $\Psi_{\mathbf n, \mathbf m}(t)$

Each unitary step is simulated in $O(N^3)$ (or better), where $N$ is the number of sites, by Gaussian-fermion techniques (Santos, 18 Jan 2026).

4. Complexity Analysis: Efficiency and Breakdown

For the imaginary-interaction line $(U = i\gamma, J \in \mathbb{R})$, the sampled operators are unitary, ensuring controlled variance and efficient convergence. The number of samples needed for accuracy $\epsilon$ with confidence $1-\delta$ is

$$M \gtrsim \frac{\Delta^2}{2\epsilon^2} \ln\left(\frac{2}{\delta}\right), \qquad \Delta \leq 1,$$

yielding $M = O(\epsilon^{-2}\log(1/\delta))$. Each sample costs $O(R\cdot\mathrm{poly}(N))$ time, with $R \sim t^2/\epsilon_{\mathrm{Trot}}^2$ to control Trotter error. Memory cost is only $O(N^2)$.

For generic complex parameters—i.e., if $\mathrm{Im}(J) \neq 0$ or $U$ has a real part—the simulation involves non-unitary similarity transforms. In this case, the operator norm grows exponentially, causing the variance of the estimator to diverge exponentially in $N t |\mathrm{Im}(J)|$. The required sample count obeys

$M \gtrsim \exp(\alpha N t |\mathrm{Im}(J)|),$

rendering the algorithm intractable for large system size, time, or nonzero $\mathrm{Im}(J)$ or $\mathrm{Re}(U)$. Efficient simulation is strictly limited to the imaginary-interaction, real-hopping axis (Santos, 18 Jan 2026).

5. Applications in Non-Hermitian Quantum Dynamics

The practical importance of imaginary-interaction Fermi-Hubbard models is multifaceted:

• Probing non-Hermitian many-body phenomena: The method provides exact time-domain simulations of systems with exceptional points, PT-broken phases, and other non-Hermitian effects.
• Quantum hardware benchmarking: Enables reference simulation for quantum-computer experiments that implement non-Hermitian or open-system dynamics.
• Analytical continuation and Lefschetz thimble methods: The tractable imaginary-interaction sector can serve as a starting point for path-integral deformation or analytic continuation approaches aimed at solving the conventional (real-interaction) Hubbard model.

A plausible implication is that this duality may facilitate new analytical and numerical schemes for exploring quantum many-body systems with engineered dissipation or complex interactions, provided the model admits a bipartite structure and the interaction lies on the imaginary axis.

6. Limitations and Rigorous Constraints

The efficient sampling algorithm is subject to strict limitations:

• Single-amplitude extraction: The method recovers one amplitude $\Psi_{\mathbf n \mathbf m}(t)$ at a time. Constructing the full wavefunction across all $$\binom{N}{N_\uparrow} \times \binom{N}{N_\downarrow}$$ configurations remains exponentially costly.
• Parameter restrictions: Efficient simulation requires precisely imaginary on-site interaction $(U = i\gamma)$ and real hopping. Any deviation (e.g., nonzero $\mathrm{Im}(J)$ or $\mathrm{Re}(U)$) induces exponential sample complexity.
• Trotter-sampling tradeoff: The Trotter error $O(t^2/R)$ and sampling cost $M$ must be balanced. In practice, $R \sim O(t^2/\epsilon)$, giving polynomial scaling with $t$ for the imaginary-interaction line but not otherwise.

This suggests that while the imaginary-interaction Fermi-Hubbard model is classically tractable in a specific parameter regime, it does not eliminate the exponential complexity generic to interacting fermion models outside this regime (Santos, 18 Jan 2026).