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

$c_i, c_i^\dagger$0

By vectorizing $c_i, c_i^\dagger$1 and introducing "right-acting" fermions $c_i, c_i^\dagger$2, the Lindbladian is represented as

$c_i, c_i^\dagger$3

After conjugating with the twist operator $c_i, c_i^\dagger$4 (flipping the sign of $c_i, c_i^\dagger$5 on $c_i, c_i^\dagger$6) and relabeling, the Lindbladian generator $c_i, c_i^\dagger$7 becomes the non-Hermitian Hubbard Hamiltonian $c_i, c_i^\dagger$8. 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 $c_i, c_i^\dagger$9 is approximated as an average over random unitaries: $\{c_i, c_j^\dagger\} = \delta_{ij}$0 where $\{c_i, c_j^\dagger\} = \delta_{ij}$1 with $\{c_i, c_j^\dagger\} = \delta_{ij}$2 (i.i.d., probability $\{c_i, c_j^\dagger\} = \delta_{ij}$3), and

$\{c_i, c_j^\dagger\} = \delta_{ij}$4

Iterating $\{c_i, c_j^\dagger\} = \delta_{ij}$5 steps of size $\{c_i, c_j^\dagger\} = \delta_{ij}$6 and averaging over $\{c_i, c_j^\dagger\} = \delta_{ij}$7 Monte Carlo samples reconstructs the time evolution up to an error $\{c_i, c_j^\dagger\} = \delta_{ij}$8. 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 $\{c_i, c_j^\dagger\} = \delta_{ij}$9
2a	For each sample $$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}.$$0, draw random signs $$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}.$$1 for each $$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}.$$2rotter step and site
2b	Initialize $$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}.$$3
2c	For $$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}.$$4, apply $$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}.$$5 as above
2d	Extract amplitude $$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}.$$6 by tracing against desired configuration
3	Report mean amplitude as estimator for $$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}.$$7

Each unitary step is simulated in $$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}.$$8 (or better), where $$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}.$$9 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$0, the sampled operators are unitary, ensuring controlled variance and efficient convergence. The number of samples needed for accuracy $U$1 with confidence $U$2 is

$U$3

yielding $U$4. Each sample costs $U$5 time, with $U$6 to control Trotter error. Memory cost is only $U$7.

For generic complex parameters—i.e., if $U$8 or $U$9 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 $i\gamma$0. The required sample count obeys

$i\gamma$1

rendering the algorithm intractable for large system size, time, or nonzero $i\gamma$2 or $i\gamma$3. 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 $i\gamma$4 at a time. Constructing the full wavefunction across all $i\gamma$5 configurations remains exponentially costly.
• Parameter restrictions: Efficient simulation requires precisely imaginary on-site interaction $i\gamma$6 and real hopping. Any deviation (e.g., nonzero $i\gamma$7 or $i\gamma$8) induces exponential sample complexity.
• Trotter-sampling tradeoff: The Trotter error $i\gamma$9 and sampling cost $$\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},$$0 must be balanced. In practice, $$\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},$$1, giving polynomial scaling with $$\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},$$2 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).