Introduction
Introduction. —Quantum impurity models describe the dynamics of correlated electronic orbitals in quantum matter and charge transport through mesoscopic devices by treating an interacting region coupled to its environment as an open quantum system [1 ] . The same framework also underpins quantum embedding approaches to correlated materials [2 ] . Accounting for general multiorbital interactions at practical cost, both in and out of equilibrium, would greatly expand the range of phenomena accessible to impurity solvers.
Treating several correlated orbitals remains challenging for established impurity solvers. Continuous-time quantum Monte Carlo (CT-QMC), including interaction-expansion (CT-INT) and hybridization-expansion (CT-HYB), provides numerically exact results but can encounter severe sign problems in models with off-diagonal hybridization and general interactions [10]. Quantics tensor trains accelerate deterministic diagrammatic expansions, but the difficulty of compressing multiorbital integrands has so far restricted this approach to single-orbital models [23 ] [11]. Separable imaginary-time diagram methods reach realistic three-band models but are most applicable to strong-coupling regimes [24 ] [13]. Even in equilibrium, a few correlated orbitals with general interactions can already be prohibitively expensive.
Dynamical and transport properties are even more challenging to obtain. Equilibrium calculations typically compute imaginary-time Green’s functions, and obtaining real-frequency spectra requires ill-posed analytic continuation. Nonequilibrium properties instead call for real-time or Keldysh-contour approaches [19 ] . Many-body perturbation theories (MBPTs) like self-consistent second-Born (2ndB) and G W GW G W [25 ] [15]. Real-time QMC suffers from a dynamical sign problem in long-time simulations [26 ] [17]. Inchworm methods reach longer times but remain extremely costly for general multiorbital systems [20 ] [19]. The fork tensor-product-state solver avoids sign problems, but its efficiency relies on orbital-diagonal hybridization [20]. Tensor-network influence-functional methods compress temporal correlations, but extending current single-impurity [27 ] [21] and two-impurity Anderson implementations [23] to general multiorbital models with off-diagonal hybridizations remains challenging.
These limitations motivate a deterministic approach to general multiorbital impurities with direct access to real-time Green’s functions. In this Letter, we introduce Coupled-Cluster Coupled Lindblad (C 3 L \mathrm{C}^{3}\mathrm{L} C 3 L T T T [12 ] . C 3 L \mathrm{C}^{3}\mathrm{L} C 3 L [30 ] [27], which reliably reproduces the bath hybridization to accuracy ε \varepsilon ε T T T polylog ( T / ε ) \operatorname{polylog}(T/\varepsilon) polylog ( T / ε )
Nevertheless, efficiently computing the steady state and Green’s functions of the resulting interacting Lindbladian is still demanding. For instance, general impurity interactions and Lindbladian couplings would complicate real-time tensor-network calculations [15 ] . A mixed-configuration approximation within the auxiliary master equation approach (AMEA) seeks to address this by averaging conditional single-orbital solutions, which could limit its description of complex interimpurity quantum fluctuations [33 , 34 ] . Coupled-cluster (CC) theory provides an efficient route by explicitly including general interorbital interactions and organizing correlations into a hierarchy of excitations. Its singles-and-doubles truncation (CCSD) is known to be accurate for weak to intermediate correlations in molecules and solids [35 , 36 , 37 , 38 , 39 , 40 ] . Its equation-of-motion (EOM) formulation provides direct access to charged and neutral excitations and Green’s functions [41 , 42 , 43 , 44 , 45 , 46 , 47 ] .
Extending CC and EOM-CC from Hamiltonians to general Lindbladians, beyond the specific settings treated previously [33 , 48 , 49 , 50 ] , requires nontrivial theoretical machinery. We formulate a new parity-resolved fermionic superoperator tensor algebra on a doubled Fock space for steady states and Green’s functions [51 , 52 , 53 ] . A matrix-free Newton–Krylov solver yields the CC steady state, and Arnoldi-based EOM-CC yields real-time and real-frequency response functions without sampling or analytic continuation. A two-orbital Anderson-impurity benchmark that took inchworm QMC tens of thousands of CPU hours [20 ] now takes our prototype C3 ^{3} 3
Problem setup. —We consider n I n_{\mathrm{I}} n I
H ^ I = ∑ i j σ h i j a ^ i σ † a ^ j σ + ∑ i j k l , σ σ ′ V i j k l a ^ i σ † a ^ j σ ′ † a ^ k σ ′ a ^ l σ . (1) \hat{H}_{\mathrm{I}} = \sum_{ij\sigma} h_{ij}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{j\sigma} + \sum_{ijkl,\sigma\sigma'} V_{ijkl}\hat{a}^{\dagger}_{i\sigma}\hat{a}^{\dagger}_{j\sigma'}\hat{a}_{k\sigma'}\hat{a}_{l\sigma}. \tag*{(1)} H ^ I = ij σ ∑ h ij a ^ iσ † a ^ j σ + ij k l , σ σ ′ ∑ V ij k l a ^ iσ † a ^ j σ ′ † a ^ k σ ′ a ^ l σ . (1) For a Gaussian bath, the lesser and greater hybridization functions Δ ≷ ( t ) \Delta^{\gtrless}(t) Δ ≷ ( t ) [51 , 52 ] . The coupled-Lindblad construction reconstructs the bath with accuracy ε \varepsilon ε T T T N p m = polylog ( T / ε ) N_{\mathrm{pm}}=\operatorname{polylog}(T/\varepsilon) N pm = polylog ( T / ε ) ρ s s \rho_{\mathrm{ss}} ρ ss
G i σ , j σ ′ R ( t ) = − i Θ ( t ) Tr { { a ^ i σ ( t ) , a ^ j σ ′ † } ρ s s G^{R}_{i\sigma,j\sigma'}(t)=-i\Theta(t)\operatorname{Tr}\{\{\hat{a}_{i\sigma}(t),\hat{a}^{\dagger}_{j\sigma'}\}\rho_{\mathrm{ss}} G iσ , j σ ′ R ( t ) = − i Θ ( t ) Tr {{ a ^ iσ ( t ) , a ^ j σ ′ † } ρ ss for the Lindblad equation d ρ d t = L [ ρ ] = − i [ H ^ C L , ρ ] + D [ ρ ] \frac{d\rho}{dt}=\mathcal{L}[\rho]=-i[\hat{H}_{\mathrm{CL}},\rho]+\mathcal{D}[\rho] d t d ρ = L [ ρ ] = − i [ H ^ CL , ρ ] + D [ ρ ]
Coupled-cluster steady state and response
With the bath reduced to a few pseudomodes, it remains to compute the steady state and Green’s functions of the interacting Lindbladian. CC achieves polynomial cost with good accuracy at intermediate correlation by parametrizing the correlated steady state in terms of excitation amplitudes t \mathbf{t} t ∣ ρ C C ⟩ = e T ^ ∣ ρ 0 ⟩ |\rho_{\mathrm{CC}}\rangle=e^{\hat{T}}|\rho_{0}\rangle ∣ ρ CC ⟩ = e T ^ ∣ ρ 0 ⟩ ∣ ∙ ⟩ |\bullet\rangle ∣ ∙ ⟩ T ^ \hat{T} T ^ ρ 0 \rho_{0} ρ 0 S 1 , 3 \mathcal{S}_{1,3} S 1 , 3 S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5
C3 ^{3} 3 Figure 1 ) addresses two distinct computational costs by reducing the number of orbitals needed to represent the bath, and using low-scaling solvers for the resulting correlated problem. We note that real-time QMC methods that treat the bath continuously can still entail costs growing exponentially with impurity size n I n_{\mathrm{I}} n I T T T [54 ] . By comparison, C3 ^{3} 3 n I n_{\mathrm{I}} n I T T T
Figure 1. Overview of the C3 ^{3} 3 polylog ( T / ε ) \operatorname{polylog}(T/\varepsilon) polylog ( T / ε )
Single- and two-impurity Anderson models. —We first benchmark systems accessible to ED and time-dependent variational principle (TDVP) calculations, starting with the particle-hole-symmetric single-impurity Anderson model with a semicircular bath. Numerical setups and additional results are given in Supplemental Material Sec. S4. In Figure 2 a (β = 100 \beta=100 β = 100 N p m = 2 N_{\mathrm{pm}}=2 N pm = 2 S 1 , 3 \mathcal{S}_{1,3} S 1 , 3 A ( ω ) = − 1 π Im G ~ R ( ω ) A(\omega)=-\frac{1}{\pi}\operatorname{Im}\tilde{G}^{R}(\omega) A ( ω ) = − π 1 Im G ~ R ( ω ) U = 2 U=2 U = 2 ∣ S ∣ = 5 |S|=5 ∣ S ∣ = 5 U U U U = 10 U=10 U = 10 [55 , 56 ] . At U = 10 U = 10 U = 10 β = 1 \beta= 1 β = 1 Figure 2 b shows that both the hybridization fit and the spectra converge exponentially with the number of pseudomodes, illustrating the compactness of the coupled-Lindblad representation [52 ] [27].
Figure 2. (a) Single-impurity spectra for U = 2 , 4 , 6 , 8 , 10 U = 2, 4, 6, 8, 10 U = 2 , 4 , 6 , 8 , 10 N p m = 2 N_{\mathrm{pm}} = 2 N pm = 2 β = 100 \beta= 100 β = 100 U = 10 U = 10 U = 10 β = 1 \beta= 1 β = 1 N p m = 1 , … , 6 N_{\mathrm{pm}} = 1,\ldots,6 N pm = 1 , … , 6 G + R ( t ) = G 11 R + G 12 R G^{R}_{+}(t) = G^{R}_{11} + G^{R}_{12} G + R ( t ) = G 11 R + G 12 R U = 2 U = 2 U = 2 β = ∞ \beta= \infty β = ∞ M t o t = 48 , 72 M_{\mathrm{tot}} = 48, 72 M tot = 48 , 72
We next benchmark a two-impurity Anderson model against TDVP propagation with a discretized unitary bath [28 ] [24]. Each impurity couples to an independent semicircular bath at β = ∞ \beta= \infty β = ∞ C 3 L C^{3}\mathrm{L} C 3 L G R ( t ) G^{R}(t) G R ( t ) T T T Figure 2 c, the TDVP bonding response G + R G^{R}_{+} G + R t ≈ 2 t \approx2 t ≈ 2 M t o t = 48 M_{\mathrm{tot}} = 48 M tot = 48 M t o t = 72 M_{\mathrm{tot}} = 72 M tot = 72 C 3 L C^{3}\mathrm{L} C 3 L t = 3 t = 3 t = 3 N p m = 4 N_{\mathrm{pm}} = 4 N pm = 4 C 3 L C^{3}\mathrm{L} C 3 L
General interactions in multiorbital impurities. —We compare with imaginary- and real-time QMC solvers for multiorbital problems with general interactions and bath hybridizations (Figure 3 a). We first consider the rotationally invariant two-orbital Kanamori interaction, relevant to two e g e_g e g [22 , 57 ] ,
Figure 3. (a) Schematic of general interorbital interactions and off-diagonal hybridization. (b) Kanamori model at β = 8 \beta= 8 β = 8 β = 64 \beta= 64 β = 64 [22 ] : G M ( τ ) G^{\mathrm{M}}(\tau) G M ( τ ) 3 ^{3} 3 A ( ω ) A(\omega) A ( ω ) 8.32 × 10 − 4 8.32 \times10^{-4} 8.32 × 1 0 − 4 β = 64 \beta= 64 β = 64 β = 0.1 \beta= 0.1 β = 0.1 [20 ] : imaginary-time Green’s functions (left); Im G i j R ( t ) \operatorname{Im}G^{R}_{ij}(t) Im G ij R ( t ) A i j ( ω ) A_{ij}(\omega) A ij ( ω ) Im G ± R \operatorname{Im}G^{R}_{\pm} Im G ± R A ± A_{\pm} A ± ± \pm ± G ± R = G 11 R + G 22 R ± G 12 R ± G 21 R 2 G^{R}_{\pm} = \frac{G^{R}_{11}+G^{R}_{22}\pm G^{R}_{12}\pm G^{R}_{21}}{2} G ± R = 2 G 11 R + G 22 R ± G 12 R ± G 21 R A ± A_{\pm} A ±
H ^ I = U ∑ i n ^ i ↑ n ^ i ↓ + ∑ i ≠ j , σ σ ′ ( U ′ − J H δ σ σ ′ ) n ^ i σ n ^ j σ ′ + J H ∑ i ≠ j ( a ^ i ↑ † a ^ j ↓ † a ^ i ↓ a ^ j ↑ + a ^ i ↑ † a ^ i ↓ † a ^ j ↑ a ^ j ↓ ) . (2) \hat{H}_{\mathrm{I}} = U \sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow} + \sum_{i\ne j,\sigma\sigma'} \left(U' - J^{\mathrm{H}}\delta_{\sigma\sigma'}\right)\hat{n}_{i\sigma}\hat{n}_{j\sigma'} + J^{\mathrm{H}}\sum_{i\ne j}\left(\hat{a}^{\dagger}_{i\uparrow}\hat{a}^{\dagger}_{j\downarrow}\hat{a}_{i\downarrow}\hat{a}_{j\uparrow} + \hat{a}^{\dagger}_{i\uparrow}\hat{a}^{\dagger}_{i\downarrow}\hat{a}_{j\uparrow}\hat{a}_{j\downarrow}\right). \tag*{(2)} H ^ I = U i ∑ n ^ i ↑ n ^ i ↓ + i = j , σ σ ′ ∑ ( U ′ − J H δ σ σ ′ ) n ^ iσ n ^ j σ ′ + J H i = j ∑ ( a ^ i ↑ † a ^ j ↓ † a ^ i ↓ a ^ j ↑ + a ^ i ↑ † a ^ i ↓ † a ^ j ↑ a ^ j ↓ ) . (2) We consider the half-filled model with U = 2 U = 2 U = 2 J H = 0.2 J^{\mathrm{H}} = 0.2 J H = 0.2 U ′ = U − 2 J H U' = U - 2J^{\mathrm{H}} U ′ = U − 2 J H C 3 L C^{3}\mathrm{L} C 3 L [58 ] and 1,200 CPU hours per temperature. Although CT-HYB generally provides numerically exact imaginary-time Green’s functions, Kanamori interactions and off-diagonal bath couplings cause a severe low-temperature sign problem [22 ] [13]. As shown in Figure 3 b, C 3 L C^{3}\mathrm{L} C 3 L S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 β = 8 \beta= 8 β = 8 β = 64 \beta= 64 β = 64 C 3 L C^{3}\mathrm{L} C 3 L
The second two-orbital case tests the solver directly against state-of-the-art real-time inchworm QMC, using the general interacting model of Ref. [20 ] :
H ^ I = ∑ i σ [ ϵ 0 n ^ i σ + V n ^ i σ ( a ^ 1 σ † a ^ 2 σ + h . c . ) ] + U ∑ i n ^ i ↑ n ^ i ↓ + ∑ σ [ V ( a ^ 1 σ † a ^ 2 σ + h . c . ) + V ′ a ^ 1 σ † a ^ 1 σ ˉ † a ^ 2 σ ˉ a ^ 2 σ ] + ν ′ ( a ^ 1 ↑ † a ^ 1 ↓ † a ^ 2 ↑ a ^ 2 ↓ + h . c . ) . (3) \begin{aligned}
\hat{H}_{\mathrm{I}} = \sum_{i\sigma}\left[\epsilon_{0}\hat{n}_{i\sigma} + V\hat{n}_{i\sigma}\left(\hat{a}^{\dagger}_{1\sigma}\hat{a}_{2\sigma} + \mathrm{h.c.}\right)\right] + U\sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow} \\
&+ \sum_{\sigma}\left[V\left(\hat{a}^{\dagger}_{1\sigma}\hat{a}_{2\sigma} + \mathrm{h.c.}\right) + V'\hat{a}^{\dagger}_{1\sigma}\hat{a}^{\dagger}_{1\bar{\sigma}}\hat{a}_{2\bar{\sigma}}\hat{a}_{2\sigma}\right] \\
&+ \nu'\left(\hat{a}^{\dagger}_{1\uparrow}\hat{a}^{\dagger}_{1\downarrow}\hat{a}_{2\uparrow}\hat{a}_{2\downarrow} + \mathrm{h.c.}\right). \tag*{(3)}
\end{aligned} H ^ I = iσ ∑ [ ϵ 0 n ^ iσ + V n ^ iσ ( a ^ 1 σ † a ^ 2 σ + h.c. ) ] + U i ∑ n ^ i ↑ n ^ i ↓ + σ ∑ [ V ( a ^ 1 σ † a ^ 2 σ + h.c. ) + V ′ a ^ 1 σ † a ^ 1 σ ˉ † a ^ 2 σ ˉ a ^ 2 σ ] + ν ′ ( a ^ 1 ↑ † a ^ 1 ↓ † a ^ 2 ↑ a ^ 2 ↓ + h.c. ) . (3) The parameters and the off-diagonal wide-band baths are specified in Supplemental Material Sec. S4. We compare C 3 L C^{3}\mathrm{L} C 3 L G W GW G W Figure 3 c, C 3 L C^{3}\mathrm{L} C 3 L C 3 L C^{3}\mathrm{L} C 3 L G W GW G W C 3 L C^{3}\mathrm{L} C 3 L S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
Resonance-resolved transport in a molecular junction. —A tunneling electron can leave a molecule in different many-body states even when it exits through the same orbital. In STM experiments on perylenetetracarboxylic dianhydride (PTCDA), bias-dependent conductance and light emission reveal selective access to neutral triplet excitons through electron removal from the anion [59 ] . Controlled single-electron spectroscopy further distinguishes transitions between molecular charge states and neutral excitations [60 ] . When an extra electron occupies the lowest unoccupied molecular orbital (LUMO), removing an electron from the highest occupied molecular orbital (HOMO) leaves two unpaired electrons. Their spins can form a neutral triplet or excited singlet, D− → ^{-} \rightarrow − → 1 _1 1 0 → _0 \rightarrow 0 → + ^{+} +
We model these competing removal processes with four spinful frontier orbitals (h2 _2 2 1 _1 1 1 _1 1 2 _2 2 Figure 4 a), with tip bias Φ \Phi Φ
Figure 4. Resonance-resolved molecular transport. (a) Frontier orbitals coupled to tip and substrate reservoirs, with charge-extraction transitions [59 ] . (b) Equilibrium G L U M O M ( τ ) G_{\mathrm{LUMO}}^{\mathrm{M}}(\tau) G LUMO M ( τ ) R R R I t I_t I t Φ = 0 \Phi= 0 Φ = 0 J t \mathcal{J}_t J t G W GW G W C 3 L C^{3}L C 3 L Φ = 0 \Phi= 0 Φ = 0
H ^ I = ∑ i j σ h i j a ^ i σ † a ^ j σ + U ∑ i n ^ i ↑ n ^ i ↓ + V ∑ i < j n ^ i n ^ j − ∑ a , b = 1 2 J h a l b E ( 2 S ^ h a ⋅ S ^ l b + 1 2 n ^ h a n ^ l b ) − ∑ a , b = 1 2 J h a l b P ( a ^ h a ↑ † a ^ h a ↓ † a ^ l b ↓ a ^ l b ↑ + h . c . ) . (4) \begin{aligned}
\hat{H}_{\mathrm{I}} ={}& \sum_{ij\sigma} h_{ij}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{j\sigma} + U\sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow} + V\sum_{i<j}\hat{n}_i\hat{n}_j \\
&- \sum_{a,b=1}^{2} J^{\mathrm{E}}_{h_a l_b}\left(2\hat{\mathbf{S}}_{h_a}\cdot\hat{\mathbf{S}}_{l_b}+\frac{1}{2}\hat{n}_{h_a}\hat{n}_{l_b}\right) \\
&- \sum_{a,b=1}^{2} J^{\mathrm{P}}_{h_a l_b}\left(\hat{a}^{\dagger}_{h_a\uparrow}\hat{a}^{\dagger}_{h_a\downarrow}\hat{a}_{l_b\downarrow}\hat{a}_{l_b\uparrow}+\mathrm{h.c.}\right). \tag*{(4)}
\end{aligned} H ^ I = ij σ ∑ h ij a ^ iσ † a ^ j σ + U i ∑ n ^ i ↑ n ^ i ↓ + V i < j ∑ n ^ i n ^ j − a , b = 1 ∑ 2 J h a l b E ( 2 S ^ h a ⋅ S ^ l b + 2 1 n ^ h a n ^ l b ) − a , b = 1 ∑ 2 J h a l b P ( a ^ h a ↑ † a ^ h a ↓ † a ^ l b ↓ a ^ l b ↑ + h.c. ) . (4) Here n ^ i : = n ^ i ↑ + n ^ i ↓ \hat{n}_i := \hat{n}_{i\uparrow}+\hat{n}_{i\downarrow} n ^ i := n ^ i ↑ + n ^ i ↓ S ^ i \hat{\mathbf{S}}_i S ^ i J E J^{\mathrm{E}} J E J P J^{\mathrm{P}} J P
The additional orbitals allow competing configurations beyond the HOMO–LUMO pair. This effective electronic impurity model is motivated by the experiments, and its parameters are given in Supplementary Material Sec. S4.
We first benchmark the equilibrium response. At Φ = 0 \Phi= 0 Φ = 0 β = 10 \beta= 10 β = 10 3 ^{3} 3 Figure 4 b). Increasing the extraction bias − Φ -\Phi − Φ I t I_t I t Figure 4 c). The driven steady state’s neutral spin response changes markedly as − Φ -\Phi − Φ R R R [61 , 62 ] . For the active orbital pair, R > 1 R > 1 R > 1 R < 1 R < 1 R < 1 R R R [59 ] .
Although CCSD, GW, and 2ndB predict similar total currents (Figure 4 c), the energy distributions of the tip current are different, as revealed by the nonequilibrium energy-resolved current J t ( ω , Φ ) J_t(\omega,\Phi) J t ( ω , Φ ) J t ( ω , Φ ) J_t(\omega,\Phi) J t ( ω , Φ ) A ≶ ( ω , Φ ) A^{\lessgtr}(\omega,\Phi) A ≶ ( ω , Φ ) [63 , 64 ] . In Figure 4 d to Figure 4 f, GW and 2ndB produce broad envelopes, whereas C3 ^{3} 3 Φ = − 0.6 \Phi= -0.6 Φ = − 0.6 S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 C 3 L C^{3}L C 3 L [55 ] , electron-gas satellites [43 ] , and X-ray absorption features [65 ] .
Exchange interactions produce multiplet structure in the removal spectrum, whose resolution requires correlating the particle–hole pairs created by electron removal. To illustrate a possible failure mechanism, we consider a minimal exchange model in Supplemental Material Sec. S5. Although the 2ndB and G W GW G W
Conclusions and outlook. —C 3 L C^{3}L C 3 L
Robust nonlinear solvers [66 ] and Krylov methods tailored to nonnormal operators [67 ] could broaden the regimes accessible to C 3 L C^{3}L C 3 L [19 ] [1]. Connections to diagrammatic many-body theory for Lindbladian Green’s functions remain to be explored [68 , 69 , 70 ] .
Extending C 3 L C^{3}L C 3 L [71 ] and trapped-ion systems [72 , 73 ] . C 3 L C^{3}L C 3 L [74 , 75 , 76 ] .
Data and Code Availability. —All C 3 L C^{3}L C 3 L [77 ] .
Acknowledgments. —This material is based upon work supported by the U.S. Department of Energy, Office of Science, Accelerated Research in Quantum Computing Centers, Quantum Utility through Advanced Computational Quantum Algorithms, grant no. DE-SC0025572 (H.L.), and by the Simons Targeted Grants in Mathematics and Physical Sciences on Moiré Materials Magic (Z.H.). H.L. thanks the Department of Computing and Mathematical Sciences at Caltech, where part of this work was done, for its hospitality. The Flatiron Institute is a division of the Simons Foundation. L.L. is a Simons Investigator in Mathematics. We thank Andre Erpenbeck, Antoine Georges, Emanuel Gull, Andy Millis, Gunhee Park, Xiao-Yu Ouyang, and Avijit Shee for helpful discussions and suggestions.
OpenAI GPT 5.6 Sol and 6 Astra were used to assist with code development and algebraic derivations under the authors’ direction and careful examination. We estimate that completing the same work without this assistance would have taken years. We take full responsibility for the content of this work.
References
End Matter
Coupled-Lindblad bath representation. —In the coupled-Lindblad construction, we introduce a finite set of auxiliary pseudomodes by fitting the greater and lesser hybridization functions as
Δ ≷ ( t ) ≃ ( g ≷ ) † e − i H ≷ t − i Γ ≷ t g ≷ (5) \Delta^{\gtrless}(t) \simeq(\mathbf{g}^{\gtrless})^\dagger e^{-i\mathbf{H}^{\gtrless}t-i\mathbf{\Gamma}^{\gtrless}t}\mathbf{g}^{\gtrless} \tag*{(5)} Δ ≷ ( t ) ≃ ( g ≷ ) † e − i H ≷ t − i Γ ≷ t g ≷ (5) where g ≷ ∈ C N p m ≷ × n 1 \mathbf{g}^{\gtrless} \in\mathbb{C}^{N_{\mathrm{pm}}^{\gtrless}\times n_1} g ≷ ∈ C N pm ≷ × n 1 H ≷ , Γ ≷ ∈ C N p m ≷ × N p m ≷ \mathbf{H}^{\gtrless},\mathbf{\Gamma}^{\gtrless} \in\mathbb{C}^{N_{\mathrm{pm}}^{\gtrless}\times N_{\mathrm{pm}}^{\gtrless}} H ≷ , Γ ≷ ∈ C N pm ≷ × N pm ≷ Γ ≷ ≥ 0 \mathbf{\Gamma}^{\gtrless}\geq0 Γ ≷ ≥ 0 N p m = N p m ≷ N_{\mathrm{pm}}=N_{\mathrm{pm}}^{\gtrless} N pm = N pm ≷
H ^ C L = H ^ I + ∑ k , l , σ ( H k l ≶ p ^ k σ † p ^ l σ + H k l ≷ d ^ k σ † d ^ l σ ) + ∑ i , k , σ [ ( g ≶ ) k i p ^ k σ † a ^ i σ + ( g ≷ ) k i d ^ k σ † a ^ i σ + h . c . ] . (6) \begin{aligned}
\hat{H}_{\mathrm{CL}}=\hat{H}_{I}+\sum_{k,l,\sigma}\left(H_{kl}^{\lessgtr}\hat{p}_{k\sigma}^{\dagger}\hat{p}_{l\sigma}+H_{kl}^{\gtrless}\hat{d}_{k\sigma}^{\dagger}\hat{d}_{l\sigma}\right) \\
&+\sum_{i,k,\sigma}\left[(\mathbf{g}^{\lessgtr})_{ki}\hat{p}_{k\sigma}^{\dagger}\hat{a}_{i\sigma}+(\mathbf{g}^{\gtrless})_{ki}\hat{d}_{k\sigma}^{\dagger}\hat{a}_{i\sigma}+\mathrm{h.c.}\right]. \tag*{(6)}
\end{aligned} H ^ CL = H ^ I + k , l , σ ∑ ( H k l ≶ p ^ k σ † p ^ l σ + H k l ≷ d ^ k σ † d ^ l σ ) + i , k , σ ∑ [ ( g ≶ ) k i p ^ k σ † a ^ iσ + ( g ≷ ) k i d ^ k σ † a ^ iσ + h.c. ] . (6) Here p ^ k σ \hat{p}_{k\sigma} p ^ k σ d ^ k σ \hat{d}_{k\sigma} d ^ k σ
D [ ρ ] = ∑ k , l , σ Γ k l < ( 2 p ^ k σ † ρ p ^ l σ − { p ^ l σ p ^ k σ † , ρ } ) + ∑ k , l , σ Γ k l > ( 2 d ^ l σ ρ d ^ k σ † − { d ^ k σ † d ^ l σ , ρ } ) . (7) \begin{aligned}
\mathcal{D}[\rho]=\sum_{k,l,\sigma}\Gamma_{kl}^{<}\left(2\hat{p}_{k\sigma}^{\dagger}\rho\hat{p}_{l\sigma}-\{\hat{p}_{l\sigma}\hat{p}_{k\sigma}^{\dagger},\rho\}\right) \\
&+\sum_{k,l,\sigma}\Gamma_{kl}^{>}\left(2\hat{d}_{l\sigma}\rho\hat{d}_{k\sigma}^{\dagger}-\{\hat{d}_{k\sigma}^{\dagger}\hat{d}_{l\sigma},\rho\}\right). \tag*{(7)}
\end{aligned} D [ ρ ] = k , l , σ ∑ Γ k l < ( 2 p ^ k σ † ρ p ^ l σ − { p ^ l σ p ^ k σ † , ρ } ) + k , l , σ ∑ Γ k l > ( 2 d ^ l σ ρ d ^ k σ † − { d ^ k σ † d ^ l σ , ρ } ) . (7) Together, the Hamiltonian and dissipator define the Lindblad dynamics d ρ d t = L [ ρ ] = − i [ H ^ C L , ρ ] + D [ ρ ] \frac{d\rho}{dt}=\mathcal{L}[\rho]=-i[\hat{H}_{\mathrm{CL}},\rho]+\mathcal{D}[\rho] d t d ρ = L [ ρ ] = − i [ H ^ CL , ρ ] + D [ ρ ] d d d n = n I + 2 d N p m n=n_I+2dN_{\mathrm{pm}} n = n I + 2 d N pm [79 ] . If the hybridization function is fitted to accuracy ε \varepsilon ε [ 0 , T ] [0,T] [ 0 , T ] ρ I ( t ) \rho_I(t) ρ I ( t ) t ∈ [ 0 , T ] t\in[0,T] t ∈ [ 0 , T ] ε → 0 \varepsilon\to0 ε → 0 [52 ] . The coupled-Lindblad construction requires only N p m = poly ( T / ε ) N_{\mathrm{pm}}=\operatorname{poly}(T/\varepsilon) N pm = poly ( T / ε )
The fitting problem admits a convex formulation that can be solved robustly and efficiently using semidefinite programming (SDP), avoiding the nonconvex optimization required by many existing approaches, which can suffer from numerical instabilities in practice [50 , 80 ] . The procedure is implemented in the publicly available RealTimeBath package on GitHub [81 ] .
Coupled cluster equations. —We vectorize ρ \rho ρ ∣ ρ ⟩ |\rho\rangle ∣ ρ ⟩ K = i L K=iL K = i L L ∣ ρ ⟩ : = ∣ L [ ρ ] ⟩ L|\rho\rangle:=|\mathcal{L}[\rho]\rangle L ∣ ρ ⟩ := ∣ L [ ρ ]⟩ c p c_p c p c ~ p \tilde{c}_p c ~ p ∣ ρ 0 ⟩ |\rho_0\rangle ∣ ρ 0 ⟩ { P κ , P ~ κ } \{P_\kappa,\tilde{P}_\kappa\} { P κ , P ~ κ } K = K 0 + W K=K_0+W K = K 0 + W
K 0 = ∑ κ λ κ P ~ κ P κ satisfies K 0 ∣ ρ 0 ⟩ = 0 , (8) K_0=\sum_{\kappa}\lambda_\kappa\tilde{P}_\kappa P_\kappa\quad\text{satisfies}\quad K_0|\rho_0\rangle=0, \tag*{(8)} K 0 = κ ∑ λ κ P ~ κ P κ satisfies K 0 ∣ ρ 0 ⟩ = 0 , (8) and W W W { P ~ κ , P κ } \{\tilde{P}_{\kappa},P_{\kappa}\} { P ~ κ , P κ } ∣ Φ S ⟩ ⟩ = τ S ∣ ρ 0 ⟩ ⟩ |\Phi_S\rangle\rangle=\tau_S|\rho_0\rangle\rangle ∣ Φ S ⟩⟩ = τ S ∣ ρ 0 ⟩⟩ τ S \tau_S τ S P ~ κ \tilde{P}_{\kappa} P ~ κ κ ∈ S \kappa\in S κ ∈ S Span { ∣ Φ S ⟩ ⟩ : Q ↑ ( S ) = Q ↓ ( S ) = 0 } \operatorname{Span}\{|\Phi_S\rangle\rangle:Q_{\uparrow}(S)=Q_{\downarrow}(S)=0\} Span { ∣ Φ S ⟩⟩ : Q ↑ ( S ) = Q ↓ ( S ) = 0 }
CC represents the correlated steady state ∣ ρ s s ⟩ ⟩ |\rho_{\mathrm{ss}}\rangle\rangle ∣ ρ ss ⟩⟩ ∣ ρ C C ⟩ ⟩ : = e T ∣ ρ 0 ⟩ ⟩ |\rho_{\mathrm{CC}}\rangle\rangle:=e^T|\rho_0\rangle\rangle ∣ ρ CC ⟩⟩ := e T ∣ ρ 0 ⟩⟩ [35 , 82 ] , with
T : = ∑ r ≥ 1 T r , T r : = ∑ ∣ S ∣ = 2 r , Q ↑ / ↓ ( S ) = 0 t S τ S , (9) T:=\sum_{r\geq1}T_r,\qquad T_r:=\sum_{\substack{|S|=2r,Q_{\uparrow/\downarrow}(S)=0}}t_S\tau_S, \tag*{(9)} T := r ≥ 1 ∑ T r , T r := ∣ S ∣ = 2 r , Q ↑ / ↓ ( S ) = 0 ∑ t S τ S , (9) where t : = { t S } t:=\{t_S\} t := { t S } r r r
S ∈ S C C S D : = { ∣ S ∣ = 2 , 4 , Q ↑ / ↓ ( S ) = 0 } , (10) S\in\mathcal{S}_{\mathrm{CCSD}}:=\{|S|=2,4,Q_{\uparrow/\downarrow}(S)=0\}, \tag*{(10)} S ∈ S CCSD := { ∣ S ∣ = 2 , 4 , Q ↑ / ↓ ( S ) = 0 } , (10) requiring O ( n 4 ) \mathcal{O}(n^4) O ( n 4 ) n n n [35 , 82 ] . The CCSD amplitude equations are obtained by projecting K ∣ ρ C C ⟩ ⟩ = 0 K|\rho_{\mathrm{CC}}\rangle\rangle=0 K ∣ ρ CC ⟩⟩ = 0 ⟨ ⟨ Φ ~ S ∣ \langle\langle\tilde{\Phi}_S| ⟨⟨ Φ ~ S ∣
z S ( t ) : = ⟨ ⟨ Φ ~ S ∣ e − T ( t ) K e T ( t ) ∣ ρ 0 ⟩ ⟩ = 0 , ∀ S ∈ S C C S D . (11) z_S(t):=\langle\langle\tilde{\Phi}_S|e^{-T(t)}Ke^{T(t)}|\rho_0\rangle\rangle=0,\quad\forall S\in\mathcal{S}_{\mathrm{CCSD}}. \tag*{(11)} z S ( t ) := ⟨⟨ Φ ~ S ∣ e − T ( t ) K e T ( t ) ∣ ρ 0 ⟩⟩ = 0 , ∀ S ∈ S CCSD . (11) The residual can be evaluated at O ( n 6 ) \mathcal{O}(n^6) O ( n 6 ) T T T ⟨ ⟨ I ∣ ρ C C ⟩ ⟩ = 1 \langle\langle I|\rho_{\mathrm{CC}}\rangle\rangle=1 ⟨⟨ I ∣ ρ CC ⟩⟩ = 1 [82 ] . The residual in Eq. (11) and the Jacobian are evaluated through low-rank tensor contractions, enabling a matrix-free Newton–Krylov solver that improves nonlinear convergence [66 , 83 ] , as detailed in Supplemental Material Sec. S2. We recall the retarded Green’s function, with p = ( i , σ ) p=(i,\sigma) p = ( i , σ )
G p q R ( t ) = − i Θ ( t ) Tr { [ a ^ p ( t ) , a ^ q † ] ρ s s } . (12) G^R_{pq}(t)=-i\Theta(t)\operatorname{Tr}\{[\hat{a}_p(t),\hat{a}^{\dagger}_q]\rho_{\mathrm{ss}}\}. \tag*{(12)} G pq R ( t ) = − i Θ ( t ) Tr {[ a ^ p ( t ) , a ^ q † ] ρ ss } . (12) Since a ^ p \hat{a}_p a ^ p a ^ p ( t ) = e t L † [ a ^ p ] \hat{a}_p(t)=e^{t\mathcal{L}^{\dagger}}[\hat{a}_p] a ^ p ( t ) = e t L † [ a ^ p ] L † \mathcal{L}^{\dagger} L †
L [ ∙ ] = Π ^ L [ Π ^ ∙ ] , Π ^ = e i π N ^ , (13) \mathcal{L}[\bullet]=\hat{\Pi}\mathcal{L}[\hat{\Pi}\bullet],\qquad\hat{\Pi}=e^{i\pi\hat{N}}, \tag*{(13)} L [ ∙ ] = Π ^ L [ Π ^ ∙ ] , Π ^ = e iπ N ^ , (13) where N ^ \hat{N} N ^ [51 , 52 , 53 ] . A Fourier transform of the projected dynamics yields the frequency-domain EOM-CC expression
G ~ p q R ( z ) = ⟨ ⟨ I ∣ c p P S E O M ( z − K ‾ ) − 1 P S E O M P o ( c ˉ q † − c ˉ q ) ∣ ρ 0 ⟩ ⟩ . (14) \widetilde{G}^{R}_{pq}(z)=\langle\langle I|c_pP_{\mathrm{SEOM}}(z-\overline{K})^{-1}P_{\mathrm{SEOM}}P_o(\bar{c}^{\dagger}_q-\bar{c}_q)|\rho_0\rangle\rangle. \tag*{(14)} G pq R ( z ) = ⟨⟨ I ∣ c p P SEOM ( z − K ) − 1 P SEOM P o ( c ˉ q † − c ˉ q ) ∣ ρ 0 ⟩⟩ . (14) In Eq. (14) , P o = ( I − Π ) / 2 P_o=(I-\Pi)/2 P o = ( I − Π ) /2 Π \Pi Π S 1 , 3 \mathcal{S}_{1,3} S 1 , 3 S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5 ∣ S ∣ = 1 , 3 |S|=1,3 ∣ S ∣ = 1 , 3 1 , 3 , 5 1,3,5 1 , 3 , 5 + 1 +1 + 1 K ‾ \overline{K} K K ‾ : = e − T ( t C C S D ) K e T ( t C C S D ) \overline{K}:=e^{-T(t_{\mathrm{CCSD}})}Ke^{T(t_{\mathrm{CCSD}})} K := e − T ( t CCSD ) K e T ( t CCSD ) O ( n 3 ) \mathcal{O}(n^3) O ( n 3 ) O ( n 5 ) \mathcal{O}(n^5) O ( n 5 )
We also consider the lesser and greater Green’s functions, defined for t ≥ 0 t\geq0 t ≥ 0
G p q < ( t ) = i Tr [ a ^ q † a ^ p ( t ) ρ s s ] , G p q > ( t ) = − i Tr [ a ^ p ( t ) a ^ q † ρ s s ] . (15) G^<_{pq}(t)=i\operatorname{Tr}[\hat{a}^{\dagger}_q\hat{a}_p(t)\rho_{\mathrm{ss}}],\qquad G^>_{pq}(t)=-i\operatorname{Tr}[\hat{a}_p(t)\hat{a}^{\dagger}_q\rho_{\mathrm{ss}}]. \tag*{(15)} G pq < ( t ) = i Tr [ a ^ q † a ^ p ( t ) ρ ss ] , G pq > ( t ) = − i Tr [ a ^ p ( t ) a ^ q † ρ ss ] . (15) The CC calculation replaces ρ s s \rho_{\mathrm{ss}} ρ ss ρ C C \rho_{\mathrm{CC}} ρ CC (14) gives the addition and removal contributions of G ~ R ( z ) \widetilde{G}^{R}(z) G R ( z )
G ~ p q R , > ( z ) = ⟨ ⟨ I ∣ c p P S E O M ( z − K ‾ ) − 1 P S E O M P o c ˉ q † ∣ ρ 0 ⟩ ⟩ , \widetilde{G}^{R,>}_{pq}(z)=\langle\langle I|c_pP_{\mathrm{SEOM}}(z-\overline{K})^{-1}P_{\mathrm{SEOM}}P_o\bar{c}^{\dagger}_q|\rho_0\rangle\rangle, G pq R , > ( z ) = ⟨⟨ I ∣ c p P SEOM ( z − K ) − 1 P SEOM P o c ˉ q † ∣ ρ 0 ⟩⟩ , G ~ p q R , < ( z ) = − ⟨ ⟨ I ∣ c p P S E O M ( z − K ‾ ) − 1 P S E O M P o c ˉ q ∣ ρ 0 ⟩ ⟩ . (16) \widetilde{G}^{R,<}_{pq}(z)=-\langle\langle I|c_pP_{\mathrm{SEOM}}(z-\overline{K})^{-1}P_{\mathrm{SEOM}}P_o\bar{c}_q|\rho_0\rangle\rangle. \tag*{(16)} G pq R , < ( z ) = − ⟨⟨ I ∣ c p P SEOM ( z − K ) − 1 P SEOM P o c ˉ q ∣ ρ 0 ⟩⟩ . (16) On the real-frequency axis, the full greater and lesser functions can be evaluated as
G ~ > ( ω ) = G ~ R , > ( ω + i 0 + ) − [ G ~ R , > ( ω + i 0 + ) ] † , \widetilde{G}^{>}(\omega)=\widetilde{G}^{R,>}(\omega+i0^+)-[\widetilde{G}^{R,>}(\omega+i0^+)]^{\dagger}, G > ( ω ) = G R , > ( ω + i 0 + ) − [ G R , > ( ω + i 0 + ) ] † , G ~ < ( ω ) = − G ~ R , < ( ω + i 0 + ) + [ G ~ R , < ( ω + i 0 + ) ] † . (17) \widetilde{G}^{<}(\omega)=-\widetilde{G}^{R,<}(\omega+i0^+)+[\widetilde{G}^{R,<}(\omega+i0^+)]^{\dagger}. \tag*{(17)} G < ( ω ) = − G R , < ( ω + i 0 + ) + [ G R , < ( ω + i 0 + ) ] † . (17) We evaluate the similarity-transformed actions through tensor contractions and solve the resulting nonnormal linear systems using matrix-free Arnoldi-type iterative methods [84 ] . Shift invariance allows the same Krylov basis to be reused across frequencies. The detailed derivation is given in Supplemental Material Sec. S3.
C3 ^3 3
H ^ I = − J E S ^ 1 ⋅ S ^ 2 + ∑ i = 1 , 2 U ( n ^ i ↑ − 1 2 ) ( n ^ i ↓ − 1 2 ) (18) \hat{H}_{I}=-J^{E}\hat{\mathbf{S}}_1\cdot\hat{\mathbf{S}}_2+\sum_{i=1,2}U\left(\hat{n}_{i\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i\downarrow}-\frac{1}{2}\right) \tag*{(18)} H ^ I = − J E S ^ 1 ⋅ S ^ 2 + i = 1 , 2 ∑ U ( n ^ i ↑ − 2 1 ) ( n ^ i ↓ − 2 1 ) (18) at half filling, with U = 4 U=4 U = 4 β = 100 \beta=100 β = 100 J E = 0 J^{E}=0 J E = 0 J E = 2 J^{E}=2 J E = 2 N p m = 6 N_{\mathrm{pm}}=6 N pm = 6 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
Figure 5 shows that C3 ^{3} 3 3 ^{3} 3
Figure 5. Low-temperature Kondo spectra at U = 4 U=4 U = 4 β = 100 \beta=100 β = 100 J E = 0 J^{E}=0 J E = 0 3 ^3 3 S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5 N p m = 6 N_{\mathrm{pm}}=6 N pm = 6
Hartree–Fock reference state and normal modes of the noninteracting Lindbladian
The pseudomode HF calculation solves a self-consistent quasi-free Lindbladian in the full impurity-plus-pseudomode one-particle space. Starting from a 1-RDM γ \gamma γ F [ γ ] F[\gamma] F [ γ ] γ ⋆ \gamma_\star γ ⋆ ρ 0 \rho_0 ρ 0 L 0 \mathcal{L}_0 L 0 { P κ , P ~ κ } \{P_\kappa,\widetilde{P}_\kappa\} { P κ , P κ }
A Lindbladian with a quadratic Hamiltonian, H = ∑ p , q F p q c ^ p † c ^ q H = \sum_{p,q} F_{pq}\hat{c}^{\dagger}_p\hat{c}_q H = ∑ p , q F pq c ^ p † c ^ q (7) ) generates quasi-free dynamics with a Gaussian steady state [85 ] . Its 1-RDM γ q p = Tr ( c ^ p † c ^ q ρ s s ) \gamma_{qp} = \operatorname{Tr}(\hat{c}^{\dagger}_p\hat{c}_q\rho_{\mathrm{ss}}) γ q p = Tr ( c ^ p † c ^ q ρ ss )
− i [ F , γ ] − { Γ < + Γ > , γ } + 2 Γ < = 0. (19) -\mathrm{i}[F,\gamma]-\{\Gamma^{<}+\Gamma^{>},\gamma\}+2\Gamma^{<}=0. \tag*{(19)} − i [ F , γ ] − { Γ < + Γ > , γ } + 2 Γ < = 0. (19) The Hartree–Fock approximation replaces the interaction by a quadratic mean-field Hamiltonian whose coefficients depend on the 1-RDM:
F [ γ ] p q = E p q + ∑ r , s V p q r s γ s r . (20) F[\gamma]_{pq}=E_{pq}+\sum_{r,s}V_{pqrs}\gamma_{sr}. \tag*{(20)} F [ γ ] pq = E pq + r , s ∑ V pq r s γ sr . (20) With this dependence, Eq. (19) becomes a nonlinear equation for γ \gamma γ γ ↦ γ ′ \gamma\mapsto\gamma' γ ↦ γ ′ F [ γ ] F[\gamma] F [ γ ] [86 , 87 ] . Spin-restricted HF calculations impose identical spin blocks and zero spin mixing at every iteration.
The Green’s function of the quasi-free HF Lindbladian can be evaluated analytically. The Hamiltonian associated with the converged noninteracting Lindbladian L 0 \mathcal{L}_0 L 0
H C L , 0 = a ^ † h 0 a ^ : = a ^ † ( F [ γ ⋆ ] ( g < ) † ( g > ) † g < H < 0 g > 0 H > ) a ^ , a ^ : = ( ( c ^ i σ ) i = 1 n I + N p m < + N p m > ) T : = ( ( a ^ i σ ) i = 1 n I , ( p ^ i σ ) i = 1 N p m < , ( d ^ i σ ) i = 1 N p m > ) T . (21) H_{\mathrm{CL},0}=\hat{\mathbf{a}}^\dagger\mathbf{h}_0\hat{\mathbf{a}}:=\hat{\mathbf{a}}^\dagger
\begin{pmatrix}
F[\gamma_\star] & (\mathbf{g}^{<})^\dagger& (\mathbf{g}^{>})^\dagger\\
\mathbf{g}^{<} & H^{<} & 0\\
\mathbf{g}^{>} & 0 & H^{>}
\end{pmatrix}
\hat{\mathbf{a}},\qquad
\hat{\mathbf{a}}:=\left((\hat{c}_{i\sigma})_{i=1}^{n_I+N_{\mathrm{pm}}^{<}+N_{\mathrm{pm}}^{>}}\right)^{\mathsf{T}}
:=\left((\hat{a}_{i\sigma})_{i=1}^{n_I},(\hat{p}_{i\sigma})_{i=1}^{N_{\mathrm{pm}}^{<}},(\hat{d}_{i\sigma})_{i=1}^{N_{\mathrm{pm}}^{>}}\right)^{\mathsf{T}}. \tag*{(21)} H CL , 0 = a ^ † h 0 a ^ := a ^ † F [ γ ⋆ ] g < g > ( g < ) † H < 0 ( g > ) † 0 H > a ^ , a ^ := ( ( c ^ iσ ) i = 1 n I + N pm < + N pm > ) T := ( ( a ^ iσ ) i = 1 n I , ( p ^ iσ ) i = 1 N pm < , ( d ^ iσ ) i = 1 N pm > ) T . (21) For the dissipator, the corresponding block matrices are
Γ < = ( 0 0 0 0 Γ < 0 0 0 0 ) , Γ > = ( 0 0 0 0 0 0 0 0 Γ > ) . (22) \Gamma^{<}=
\begin{pmatrix}
0&0&0\\
0&\Gamma^{<}&0\\
0&0&0
\end{pmatrix},
\qquad
\Gamma^{>}=
\begin{pmatrix}
0&0&0\\
0&0&0\\
0&0&\Gamma^{>}
\end{pmatrix}. \tag*{(22)} Γ < = 0 0 0 0 Γ < 0 0 0 0 , Γ > = 0 0 0 0 0 0 0 0 Γ > . (22) The odd-parity part L 0 \mathcal{L}_0 L 0 Z 2 \mathbb{Z}_2 Z 2
L 0 [ ∙ ] = − i [ H C L , 0 , ∙ ] − ∑ i , j , σ Γ i j < ( 2 p ^ i σ † ∙ p ^ j σ + { p ^ j σ p ^ i σ † , ∙ } ) − ∑ i , j , σ Γ i j > ( 2 d ^ j σ ∙ d ^ i σ † + { d ^ i σ † d ^ j σ , ∙ } ) . (23) \mathcal{L}_0[\bullet]=-\mathrm{i}[H_{\mathrm{CL},0},\bullet]-\sum_{i,j,\sigma}\Gamma^{<}_{ij}\left(2\hat{p}^{\dagger}_{i\sigma}\bullet\hat{p}_{j\sigma}+\{\hat{p}_{j\sigma}\hat{p}^{\dagger}_{i\sigma},\bullet\}\right)-\sum_{i,j,\sigma}\Gamma^{>}_{ij}\left(2\hat{d}_{j\sigma}\bullet\hat{d}^{\dagger}_{i\sigma}+\{\hat{d}^{\dagger}_{i\sigma}\hat{d}_{j\sigma},\bullet\}\right). \tag*{(23)} L 0 [ ∙ ] = − i [ H CL , 0 , ∙ ] − i , j , σ ∑ Γ ij < ( 2 p ^ iσ † ∙ p ^ j σ + { p ^ j σ p ^ iσ † , ∙ } ) − i , j , σ ∑ Γ ij > ( 2 d ^ j σ ∙ d ^ iσ † + { d ^ iσ † d ^ j σ , ∙ } ) . (23) For the noninteracting retarded Green’s function G 0 , p q R ( t ) = − i Θ ( t ) Tr ( { c ^ p ( t ) , c ^ q † ( 0 ) } ρ 0 ) = : − i Θ ( t ) Y ( t ) G^R_{0,pq}(t)=-i\Theta(t)\operatorname{Tr}(\{\hat{c}_p(t),\hat{c}^{\dagger}_q(0)\}\rho_0) =: -i\Theta(t)Y(t) G 0 , pq R ( t ) = − i Θ ( t ) Tr ({ c ^ p ( t ) , c ^ q † ( 0 )} ρ 0 ) =: − i Θ ( t ) Y ( t )
d d t Y ( t ) = [ Tr ( { L 0 † [ c ^ p ( t ) ] , c ^ q † ( 0 ) } ρ 0 ) ] p , q = ( − i H 0 − Γ < − Γ > ) Y ( t ) . (24) \frac{\mathrm{d}}{\mathrm{d}t}Y(t)=\left[\operatorname{Tr}\left(\{\mathcal{L}_0^\dagger[\hat{c}_p(t)],\hat{c}^{\dagger}_q(0)\}\rho_0\right)\right]_{p,q}=(-iH_0-\Gamma^{<}-\Gamma^{>})Y(t). \tag*{(24)} d t d Y ( t ) = [ Tr ( { L 0 † [ c ^ p ( t )] , c ^ q † ( 0 )} ρ 0 ) ] p , q = ( − i H 0 − Γ < − Γ > ) Y ( t ) . (24) Hence G 0 R ( t ) = − i Θ ( t ) e ( − i H 0 − Γ < − Γ > ) t G^R_0(t)=-i\Theta(t)e^{(-iH_0-\Gamma^{<}-\Gamma^{>})t} G 0 R ( t ) = − i Θ ( t ) e ( − i H 0 − Γ < − Γ > ) t
G ~ 0 R ( z ) = ∫ 0 ∞ e i z t G 0 R ( t ) d t = ( z − H 0 + i Γ < + i Γ > ) − 1 = [ z − ( F [ γ ⋆ ] ( g < ) † ( g > ) † g < H < − i Γ < 0 g > 0 H > − i Γ > ) ] − 1 . (25) \widetilde{G}^R_0(z)=\int_0^\infty e^{izt}G^R_0(t)\,\mathrm{d}t=(z-H_0+i\Gamma^{<}+i\Gamma^{>})^{-1}=\left[z-\begin{pmatrix}F[\gamma^\star]&(\mathbf{g}^{<})^\dagger&(\mathbf{g}^{>})^\dagger\\ \mathbf{g}^{<}&H^{<}-i\Gamma^{<}&0\\ \mathbf{g}^{>}&0&H^{>}-i\Gamma^{>}\end{pmatrix}\right]^{-1}. \tag*{(25)} G 0 R ( z ) = ∫ 0 ∞ e i z t G 0 R ( t ) d t = ( z − H 0 + i Γ < + i Γ > ) − 1 = z − F [ γ ⋆ ] g < g > ( g < ) † H < − i Γ < 0 ( g > ) † 0 H > − i Γ > − 1 . (25) The impurity Green’s function is the ( 1 , 1 ) (1,1) ( 1 , 1 ) G ~ 0 R ( z ) \widetilde{G}^R_0(z) G 0 R ( z )
G ~ 0 , 11 R ( z ) = ( z I − F [ γ ⋆ ] − Δ < ( z ) − Δ > ( z ) ) − 1 , Δ ≶ ( z ) = ( g ≶ ) † ( z I − H ≶ + i Γ ≶ ) − 1 g ≶ . (26) \widetilde{G}^R_{0,11}(z)=(zI-F[\gamma^\star]-\Delta^{<}(z)-\Delta^{>}(z))^{-1},\qquad\Delta^{\lessgtr}(z)=(\mathbf{g}^{\lessgtr})^\dagger(zI-H^{\lessgtr}+i\Gamma^{\lessgtr})^{-1}\mathbf{g}^{\lessgtr}. \tag*{(26)} G 0 , 11 R ( z ) = ( z I − F [ γ ⋆ ] − Δ < ( z ) − Δ > ( z ) ) − 1 , Δ ≶ ( z ) = ( g ≶ ) † ( z I − H ≶ + i Γ ≶ ) − 1 g ≶ . (26) We next construct the normal modes that define the HF reference for the post-HF CCSD calculation. The fermionic superoperators { c ^ p , c ^ p † , c ~ p , c ~ p † } \{\hat{c}_p,\hat{c}^{\dagger}_p,\widetilde{c}_p,\widetilde{c}^{\dagger}_p\} { c ^ p , c ^ p † , c p , c p † }
c p : ∣ ρ ⟩ ⟩ ↦ ∣ c ^ p ρ ⟩ ⟩ , c p † : ∣ ρ ⟩ ⟩ ↦ ∣ c ^ p † ρ ⟩ ⟩ , c ~ p : ∣ ρ ⟩ ⟩ ↦ ∣ − i Π ^ c ^ p † ρ Π ^ ⟩ ⟩ , c ~ p † : ∣ ρ ⟩ ⟩ ↦ ∣ − i Π ^ c ^ p ρ Π ^ ⟩ ⟩ . (27) c_p:|\rho\rangle\rangle\mapsto|\hat{c}_p\rho\rangle\rangle,\qquad c_p^\dagger:|\rho\rangle\rangle\mapsto|\hat{c}_p^\dagger\rho\rangle\rangle,\qquad\widetilde{c}_p:|\rho\rangle\rangle\mapsto|-i\hat{\Pi}\hat{c}_p^\dagger\rho\hat{\Pi}\rangle\rangle,\qquad\widetilde{c}_p^\dagger:|\rho\rangle\rangle\mapsto|-i\hat{\Pi}\hat{c}_p\rho\hat{\Pi}\rangle\rangle. \tag*{(27)} c p : ∣ ρ ⟩⟩ ↦ ∣ c ^ p ρ ⟩⟩ , c p † : ∣ ρ ⟩⟩ ↦ ∣ c ^ p † ρ ⟩⟩ , c p : ∣ ρ ⟩⟩ ↦ ∣ − i Π ^ c ^ p † ρ Π ^ ⟩⟩ , c p † : ∣ ρ ⟩⟩ ↦ ∣ − i Π ^ c ^ p ρ Π ^ ⟩⟩ . (27) They satisfy the CAR because ${\hat{\Pi},\hat{c}_p}={\hat{\Pi},\hat{c}^{\dagger}_p}=0.
We partition the Lindbladian into noninteracting and interacting contributions,
K : = i L = K 0 + W , (28) K:=iL=K_0+W, \tag*{(28)} K := i L = K 0 + W , (28) where the noninteracting part K 0 K_0 K 0 K 0 = ∑ σ c σ ‡ h σ c σ K_0=\sum_\sigma c^\ddagger_\sigma h_\sigma c_\sigma K 0 = ∑ σ c σ ‡ h σ c σ c σ = ( c 0 σ , ⋯ , c n − 1 , σ , c ~ 0 σ † , ⋯ , c ~ n − 1 , σ † ) T c_\sigma=(c_{0\sigma},\cdots,c_{n-1,\sigma},\widetilde{c}^{\dagger}_{0\sigma},\cdots,\widetilde{c}^{\dagger}_{n-1,\sigma})^\mathsf{T} c σ = ( c 0 σ , ⋯ , c n − 1 , σ , c 0 σ † , ⋯ , c n − 1 , σ † ) T c σ ‡ = ( c 0 σ † , ⋯ , c n − 1 , σ † , c ~ 0 σ , ⋯ , c ~ n − 1 , σ ) c^\ddagger_\sigma=(c^\dagger_{0\sigma},\cdots,c^\dagger_{n-1,\sigma},\widetilde{c}_{0\sigma},\cdots,\widetilde{c}_{n-1,\sigma}) c σ ‡ = ( c 0 σ † , ⋯ , c n − 1 , σ † , c 0 σ , ⋯ , c n − 1 , σ ) h σ h_\sigma h σ 2 n × 2 n 2n\times2n 2 n × 2 n ζ = 0 , … , 2 n − 1 \zeta=0,\ldots,2n-1 ζ = 0 , … , 2 n − 1 h σ h_\sigma h σ h σ V σ = V σ diag ( ε ξ , σ ) ξ = 0 2 n − 1 h_\sigma V_\sigma=V_\sigma\operatorname{diag}(\varepsilon_{\xi,\sigma})_{\xi=0}^{2n-1} h σ V σ = V σ diag ( ε ξ , σ ) ξ = 0 2 n − 1
P ξ σ ⋆ = ∑ ζ = 0 2 n − 1 ( V σ − 1 ) ξ ζ c ζ σ , P ~ ξ σ ⋆ = ∑ ζ = 0 2 n − 1 ( V σ ) ζ ξ c ζ σ ‡ . (29) P^\star_{\xi\sigma}=\sum_{\zeta=0}^{2n-1}(V_\sigma^{-1})_{\xi\zeta}c_{\zeta\sigma},\qquad\widetilde{P}^\star_{\xi\sigma}=\sum_{\zeta=0}^{2n-1}(V_\sigma)_{\zeta\xi}c^\ddagger_{\zeta\sigma}. \tag*{(29)} P ξ σ ⋆ = ζ = 0 ∑ 2 n − 1 ( V σ − 1 ) ξ ζ c ζ σ , P ξ σ ⋆ = ζ = 0 ∑ 2 n − 1 ( V σ ) ζ ξ c ζ σ ‡ . (29) These modes satisfy the biorthogonal CAR, { P ξ σ ⋆ , P ~ ξ ′ σ ′ ⋆ } = δ ξ ξ ′ δ σ σ ′ \{P^\star_{\xi\sigma},\widetilde{P}^\star_{\xi'\sigma'}\}=\delta_{\xi\xi'}\delta_{\sigma\sigma'} { P ξ σ ⋆ , P ξ ′ σ ′ ⋆ } = δ ξ ξ ′ δ σ σ ′ { P ξ σ ⋆ , P ξ ′ σ ′ ⋆ } = { P ~ ξ σ ⋆ , P ~ ξ ′ σ ′ ⋆ } = 0 \{P^\star_{\xi\sigma},P^\star_{\xi'\sigma'}\}=\{\widetilde{P}^\star_{\xi\sigma},\widetilde{P}^\star_{\xi'\sigma'}\}=0 { P ξ σ ⋆ , P ξ ′ σ ′ ⋆ } = { P ξ σ ⋆ , P ξ ′ σ ′ ⋆ } = 0 P ξ σ P_{\xi\sigma} P ξ σ P ~ ξ σ \widetilde{P}_{\xi\sigma} P ξ σ λ ξ , σ \lambda_{\xi,\sigma} λ ξ , σ q ξ , σ q_{\xi,\sigma} q ξ , σ
( P ξ σ , P ~ ξ σ , λ ξ σ , q ξ σ ) : = { ( P ξ σ ⋆ , P ~ ξ σ ⋆ , ε ξ σ , 1 ) , if Im ε ξ σ < 0 , ( P ~ ξ σ ⋆ , P ξ σ ⋆ , − ε ξ σ , − 1 ) , if Im ε ξ σ > 0. (30) (P_{\xi\sigma},\widetilde{P}_{\xi\sigma},\lambda_{\xi\sigma},q_{\xi\sigma}):=
\begin{cases}
(P^\star_{\xi\sigma},\widetilde{P}^\star_{\xi\sigma},\varepsilon_{\xi\sigma},1),&\text{if }\operatorname{Im}\varepsilon_{\xi\sigma}<0,\\
(\widetilde{P}^\star_{\xi\sigma},P^\star_{\xi\sigma},-\varepsilon_{\xi\sigma},-1),&\text{if }\operatorname{Im}\varepsilon_{\xi\sigma}>0.
\end{cases} \tag*{(30)} ( P ξ σ , P ξ σ , λ ξ σ , q ξ σ ) := { ( P ξ σ ⋆ , P ξ σ ⋆ , ε ξ σ , 1 ) , ( P ξ σ ⋆ , P ξ σ ⋆ , − ε ξ σ , − 1 ) , if Im ε ξ σ < 0 , if Im ε ξ σ > 0. (30) The final modes { P ξ σ , P ~ ξ σ } \{P_{\xi\sigma},\widetilde{P}_{\xi\sigma}\} { P ξ σ , P ξ σ } P ξ σ ⋆ P^\star_{\xi\sigma} P ξ σ ⋆ P ~ ξ σ ⋆ \widetilde{P}^\star_{\xi\sigma} P ξ σ ⋆ K 0 K_0 K 0 K 0 = ∑ κ = 1 4 n λ κ P ~ κ P κ K_0=\sum_{\kappa=1}^{4n}\lambda_\kappa\widetilde{P}_\kappa P_\kappa K 0 = ∑ κ = 1 4 n λ κ P κ P κ κ = ( σ , ξ ) ∈ { ↑ , ↓ } × { 0 , … , 2 n − 1 } \kappa=(\sigma,\xi)\in\{\uparrow,\downarrow\}\times\{0,\ldots,2n-1\} κ = ( σ , ξ ) ∈ { ↑ , ↓ } × { 0 , … , 2 n − 1 } ∣ ρ 0 ⟩ ⟩ |\rho_0\rangle\rangle ∣ ρ 0 ⟩⟩
P κ ∣ ρ 0 ⟩ ⟩ = 0 , ⟨ ⟨ I ∣ P ~ κ = 0 , ⟨ ⟨ I ∣ ρ 0 ⟩ ⟩ = Tr ( ρ 0 ) = 1 , ∀ κ . (31) P_\kappa|\rho_0\rangle\rangle=0,\qquad\langle\langle I|\widetilde{P}_\kappa=0,\qquad\langle\langle I|\rho_0\rangle\rangle=\operatorname{Tr}(\rho_0)=1,\qquad\forall\kappa. \tag*{(31)} P κ ∣ ρ 0 ⟩⟩ = 0 , ⟨⟨ I ∣ P κ = 0 , ⟨⟨ I ∣ ρ 0 ⟩⟩ = Tr ( ρ 0 ) = 1 , ∀ κ . (31) The configurations ∣ Φ S ⟩ ⟩ |\Phi_S\rangle\rangle ∣ Φ S ⟩⟩
∣ Φ S ⟩ ⟩ = τ S ∣ ρ 0 ⟩ ⟩ = ∏ κ ∈ S ↗ P ~ κ ∣ ρ 0 ⟩ ⟩ . (32) |\Phi_S\rangle\rangle=\tau_S|\rho_0\rangle\rangle=\overset{\nearrow}{\prod_{\kappa\in S}}\widetilde{P}_\kappa|\rho_0\rangle\rangle. \tag*{(32)} ∣ Φ S ⟩⟩ = τ S ∣ ρ 0 ⟩⟩ = κ ∈ S ∏ ↗ P κ ∣ ρ 0 ⟩⟩ . (32) where ↗ \nearrow ↗ ( ↑ , + ) ≺ ( ↑ , − ) ≺ ( ↓ , + ) ≺ ( ↓ , − ) (\uparrow,+)\prec(\uparrow,-)\prec(\downarrow,+)\prec(\downarrow,-) ( ↑ , + ) ≺ ( ↑ , − ) ≺ ( ↓ , + ) ≺ ( ↓ , − ) S = ( κ 1 ≺ ⋯ ≺ κ ∣ S ∣ ) S=(\kappa_1\prec\cdots\prec\kappa_{|S|}) S = ( κ 1 ≺ ⋯ ≺ κ ∣ S ∣ ) ξ \xi ξ q ξ σ q_{\xi\sigma} q ξ σ τ S \tau_S τ S P ~ κ \widetilde{P}_\kappa P κ
Due to the U ( 1 ) U(1) U ( 1 ) ρ s s \rho_{\mathrm{ss}} ρ ss [ N ^ σ , ρ s s ] = 0 [\hat{N}_\sigma,\rho_{\mathrm{ss}}]=0 [ N ^ σ , ρ ss ] = 0
( N σ − N ~ σ ) ∣ ρ s s ⟩ ⟩ = 0 , where N σ = ∑ i = 0 n − 1 c i σ † c i σ , N ~ σ = ∑ i = 0 n − 1 c ~ i σ † c ~ i σ . (33) (N_\sigma-\widetilde{N}_\sigma)|\rho_{\mathrm{ss}}\rangle\rangle=0,\qquad\text{where }N_\sigma=\sum_{i=0}^{n-1}c^\dagger_{i\sigma}c_{i\sigma},\qquad\widetilde{N}_\sigma=\sum_{i=0}^{n-1}\widetilde{c}^{\dagger}_{i\sigma}\widetilde{c}_{i\sigma}. \tag*{(33)} ( N σ − N σ ) ∣ ρ ss ⟩⟩ = 0 , where N σ = i = 0 ∑ n − 1 c iσ † c iσ , N σ = i = 0 ∑ n − 1 c iσ † c iσ . (33) Thus the steady state lies in the charge-neutral sector of the doubled Fock space
∣ ρ s s ⟩ ⟩ ∈ Span { ∣ Φ S ⟩ ⟩ : Q ↑ ( S ) = Q ↓ ( S ) = 0 } , (34) |\rho_{ss}\rangle\rangle\in\operatorname{Span}\{|\Phi_S\rangle\rangle: Q_{\uparrow}(S)=Q_{\downarrow}(S)=0\}, \tag*{(34)} ∣ ρ ss ⟩⟩ ∈ Span { ∣ Φ S ⟩⟩ : Q ↑ ( S ) = Q ↓ ( S ) = 0 } , (34) where the total particle–hole charge of a configuration S S S Q ↑ / ↓ ( S ) = ∑ κ ∈ S , σ ( κ ) = ↑ / ↓ q κ Q_{\uparrow/\downarrow}(S)=\sum_{\kappa\in S,\sigma(\kappa)=\uparrow/\downarrow}q_{\kappa} Q ↑ / ↓ ( S ) = ∑ κ ∈ S , σ ( κ ) =↑ / ↓ q κ
The dual configuration is ⟨ ⟨ Φ ~ S ∣ = ⟨ ⟨ I ∣ ∏ κ ∈ S ← P κ \langle\langle\widetilde{\Phi}_S|=\langle\langle I|\prod_{\kappa\in S}^{\leftarrow}P_{\kappa} ⟨⟨ Φ S ∣ = ⟨⟨ I ∣ ∏ κ ∈ S ← P κ (32) . The left and right configurations are biorthogonal ⟨ ⟨ Φ ~ S ∣ Φ S ′ ⟩ ⟩ = δ S S ′ \langle\langle\widetilde{\Phi}_S|\Phi_{S'}\rangle\rangle=\delta_{SS'} ⟨⟨ Φ S ∣ Φ S ′ ⟩⟩ = δ S S ′
CCSD equations for interacting Lindbladians
Overview of the CCSD equations
The exact steady state ∣ ρ s s ⟩ ⟩ |\rho_{ss}\rangle\rangle ∣ ρ ss ⟩⟩ L = − i K L=-iK L = − i K
K ∣ ρ s s ⟩ ⟩ = 0 , ⟨ ⟨ I ∣ ρ s s ⟩ ⟩ = Tr ( ρ s s ) = 1. (35) K|\rho_{ss}\rangle\rangle=0,\qquad\langle\langle I|\rho_{ss}\rangle\rangle=\operatorname{Tr}(\rho_{ss})=1. \tag*{(35)} K ∣ ρ ss ⟩⟩ = 0 , ⟨⟨ I ∣ ρ ss ⟩⟩ = Tr ( ρ ss ) = 1. (35) CCSD uses the exponential ansatz
∣ ρ C C S D ⟩ ⟩ = e T ∣ ρ 0 ⟩ ⟩ , T = T 1 + T 2 , T r = ∑ ∣ S ∣ = 2 r , Q ↑ / ↓ ( S ) = 0 t S τ S . (36) |\rho_{\mathrm{CCSD}}\rangle\rangle=e^T|\rho_0\rangle\rangle,\qquad T=T_1+T_2,\qquad T_r=\sum_{\substack{|S|=2r,Q_{\uparrow/\downarrow}(S)=0}}t_S\tau_S. \tag*{(36)} ∣ ρ CCSD ⟩⟩ = e T ∣ ρ 0 ⟩⟩ , T = T 1 + T 2 , T r = ∣ S ∣ = 2 r , Q ↑ / ↓ ( S ) = 0 ∑ t S τ S . (36) Here T 1 T_1 T 1 T 2 T_2 T 2 q = + 1 q=+1 q = + 1 q = − 1 q=-1 q = − 1 a , b a,b a , b q = + 1 q=+1 q = + 1 i , j i,j i , j q = − 1 q=-1 q = − 1
T 1 = ∑ σ , a i t a i σ P ~ a σ + P ~ i σ − , T 2 = ∑ a i b j t a i , b j ↑ ↓ P ~ a ↑ + P ~ i ↑ − P ~ b ↓ + P ~ j ↓ − + 1 4 ∑ σ a b i j t a b , i j σ σ P ~ a σ + P ~ b σ + P ~ i σ − P ~ j σ − . (37) T_1=\sum_{\sigma,ai}t_{ai}^{\sigma}\widetilde{P}_{a\sigma}^{+}\widetilde{P}_{i\sigma}^{-},\qquad T_2=\sum_{aibj}t_{ai,bj}^{\uparrow\downarrow}\widetilde{P}_{a\uparrow}^{+}\widetilde{P}_{i\uparrow}^{-}\widetilde{P}_{b\downarrow}^{+}\widetilde{P}_{j\downarrow}^{-}+\frac{1}{4}\sum_{\sigma abij}t_{ab,ij}^{\sigma\sigma}\widetilde{P}_{a\sigma}^{+}\widetilde{P}_{b\sigma}^{+}\widetilde{P}_{i\sigma}^{-}\widetilde{P}_{j\sigma}^{-}. \tag*{(37)} T 1 = σ , ai ∑ t ai σ P aσ + P iσ − , T 2 = aibj ∑ t ai , bj ↑↓ P a ↑ + P i ↑ − P b ↓ + P j ↓ − + 4 1 σ abij ∑ t ab , ij σ σ P aσ + P bσ + P iσ − P j σ − . (37) Interchanging two equal-sign creators changes the sign of their product, so the full amplitude tensors satisfy t a b , i j σ σ = − t b a , i j σ σ = − t a b , j i σ σ t_{ab,ij}^{\sigma\sigma}=-t_{ba,ij}^{\sigma\sigma}=-t_{ab,ji}^{\sigma\sigma} t ab , ij σ σ = − t ba , ij σ σ = − t ab , j i σ σ
# T 1 amplitudes = 2 n 2 , # T 2 amplitudes = n 4 + 2 ( n 2 ) 2 , # CCSD amplitudes = 2 n 2 + n 4 + 2 ( n 2 ) 2 . (38) \#T_1\text{ amplitudes}=2n^2,\qquad\#T_2\text{ amplitudes}=n^4+2\binom{n}{2}^{2},\qquad\#\text{CCSD amplitudes}=2n^2+n^4+2\binom{n}{2}^{2}. \tag*{(38)} # T 1 amplitudes = 2 n 2 , # T 2 amplitudes = n 4 + 2 ( 2 n ) 2 , # CCSD amplitudes = 2 n 2 + n 4 + 2 ( 2 n ) 2 . (38) Every term in T T T (31) , ⟨ ⟨ I ∣ T = 0 \langle\langle I|T=0 ⟨⟨ I ∣ T = 0 ⟨ ⟨ I ∣ e T = ⟨ ⟨ I ∣ \langle\langle I|e^T=\langle\langle I| ⟨⟨ I ∣ e T = ⟨⟨ I ∣
Tr ( ρ C C ) = ⟨ ⟨ I ∣ ρ C C ⟩ ⟩ = ⟨ ⟨ I ∣ e T ∣ ρ 0 ⟩ ⟩ = ⟨ ⟨ I ∣ ρ 0 ⟩ ⟩ = 1. (39) \operatorname{Tr}(\rho_{\mathrm{CC}})=\langle\langle I|\rho_{\mathrm{CC}}\rangle\rangle=\langle\langle I|e^T|\rho_0\rangle\rangle=\langle\langle I|\rho_0\rangle\rangle=1. \tag*{(39)} Tr ( ρ CC ) = ⟨⟨ I ∣ ρ CC ⟩⟩ = ⟨⟨ I ∣ e T ∣ ρ 0 ⟩⟩ = ⟨⟨ I ∣ ρ 0 ⟩⟩ = 1. (39) The CC ansatz therefore preserves trace normalization without a separate constraint.
The CCSD amplitude equations follow by projecting the similarity-transformed residual ∣ r ( t C C S D ) ⟩ ⟩ = K ‾ ( t C C S D ) ∣ ρ 0 ⟩ ⟩ |r(t_{\mathrm{CCSD}})\rangle\rangle=\overline{K}(t_{\mathrm{CCSD}})|\rho_0\rangle\rangle ∣ r ( t CCSD )⟩⟩ = K ( t CCSD ) ∣ ρ 0 ⟩⟩ K ‾ ( t C C S D ) = e − T ( t C C S D ) K e T ( t C C S D ) \overline{K}(t_{\mathrm{CCSD}})=e^{-T(t_{\mathrm{CCSD}})}Ke^{T(t_{\mathrm{CCSD}})} K ( t CCSD ) = e − T ( t CCSD ) K e T ( t CCSD )
z S ( t C C S D ) : = ⟨ ⟨ Φ ~ S ∣ r ( t C C S D ) ⟩ ⟩ = ⟨ ⟨ Φ ~ S ∣ K ‾ ( t C C S D ) ∣ ρ 0 ⟩ ⟩ = 0 , ∀ S ∈ S C C S D : = { S : ∣ S ∣ = 2 , 4 , Q ↑ / ↓ ( S ) = 0 } . (40) z_S(t_{\mathrm{CCSD}}):=\langle\langle\widetilde{\Phi}_S|r(t_{\mathrm{CCSD}})\rangle\rangle=\langle\langle\widetilde{\Phi}_S|\overline{K}(t_{\mathrm{CCSD}})|\rho_0\rangle\rangle=0,\qquad\forall S\in\mathcal{S}_{\mathrm{CCSD}}:=\{S:|S|=2,4,\ Q_{\uparrow/\downarrow}(S)=0\}. \tag*{(40)} z S ( t CCSD ) := ⟨⟨ Φ S ∣ r ( t CCSD )⟩⟩ = ⟨⟨ Φ S ∣ K ( t CCSD ) ∣ ρ 0 ⟩⟩ = 0 , ∀ S ∈ S CCSD := { S : ∣ S ∣ = 2 , 4 , Q ↑ / ↓ ( S ) = 0 } . (40) Derivation of the CCSD equations using exterior algebra
The CCSD equations can be derived through several equivalent formalisms. One standard route is to normal-order the Baker–Campbell–Hausdorff (BCH) expansion in terms of particle–hole creation and annihilation operators and evaluate the resulting products via Wick’s theorem [88 , 89 ] . Alternatively, contractions can also be organized as Goldstone coupled-cluster diagrams [89 , 90 ] . In the non-normal Lindbladian setting, however, both algebraic normal ordering and diagrammatic expansions demand additional bookkeeping [91 ] . Because P κ P_{\kappa} P κ P ~ κ \widetilde{P}_{\kappa} P κ ⟨ ⟨ I ∣ \langle\langle I| ⟨⟨ I ∣ ∣ ρ 0 ⟩ |\rho_0\rangle ∣ ρ 0 ⟩
Wick’s theorem remains valid for the superoperator CAR algebra, and the BCH expansion of K ‾ \overline{K} K [89 Sec. 10.1] :
K ‾ = K + [ K , T ] + 1 2 [ [ K , T ] , T ] + 1 6 [ [ [ K , T ] , T ] , T ] + 1 24 [ [ [ [ K , T ] , T ] , T ] , T ] (41) \overline{K}=K+[K,T]+\frac{1}{2}[[K,T],T]+\frac{1}{6}[[[K,T],T],T]+\frac{1}{24}[[[[K,T],T],T],T] \tag*{(41)} K = K + [ K , T ] + 2 1 [[ K , T ] , T ] + 6 1 [[[ K , T ] , T ] , T ] + 24 1 [[[[ K , T ] , T ] , T ] , T ] (41) We first express the physical superoperators in terms of normal modes. For each component ζ \zeta ζ c σ c_{\sigma} c σ c σ † c_{\sigma}^{\dagger} c σ † ξ = 0 , … , 2 n − 1 \xi=0,\ldots,2n-1 ξ = 0 , … , 2 n − 1 ϵ ξ = ϵ ξ σ \epsilon_{\xi}=\epsilon_{\xi\sigma} ϵ ξ = ϵ ξ σ
c ζ σ = ∑ Im ϵ ξ < 0 ( V σ ) ζ ξ P ξ σ + ∑ Im ϵ ξ > 0 ( V σ ) ζ ξ P ~ ξ σ , c ζ σ † = ∑ Im ϵ ξ < 0 ( V σ − 1 ) ξ ζ P ~ ξ σ + ∑ Im ϵ ξ > 0 ( V σ − 1 ) ξ ζ P ξ σ (42) c_{\zeta\sigma}=\sum_{\operatorname{Im}\epsilon_{\xi}<0}(V_{\sigma})_{\zeta\xi}P_{\xi\sigma}+\sum_{\operatorname{Im}\epsilon_{\xi}>0}(V_{\sigma})_{\zeta\xi}\widetilde{P}_{\xi\sigma},\qquad c_{\zeta\sigma}^{\dagger}=\sum_{\operatorname{Im}\epsilon_{\xi}<0}(V_{\sigma}^{-1})_{\xi\zeta}\widetilde{P}_{\xi\sigma}+\sum_{\operatorname{Im}\epsilon_{\xi}>0}(V_{\sigma}^{-1})_{\xi\zeta}P_{\xi\sigma} \tag*{(42)} c ζ σ = Im ϵ ξ < 0 ∑ ( V σ ) ζ ξ P ξ σ + Im ϵ ξ > 0 ∑ ( V σ ) ζ ξ P ξ σ , c ζ σ † = Im ϵ ξ < 0 ∑ ( V σ − 1 ) ξ ζ P ξ σ + Im ϵ ξ > 0 ∑ ( V σ − 1 ) ξ ζ P ξ σ (42) Each of { c p , c p † , c ~ p , c ~ p † } \{c_{p},c_{p}^{\dagger},\widetilde{c}_{p},\widetilde{c}_{p}^{\dagger}\} { c p , c p † , c p , c p † } (42) into W = W ( c p , c p † , c ~ p , c ~ p † ) W=W(c_{p},c_{p}^{\dagger},\widetilde{c}_{p},\widetilde{c}_{p}^{\dagger}) W = W ( c p , c p † , c p , c p † )
W = ∑ α g α F α 1 ⋯ F α ℓ α , F α j = ∑ κ = 1 4 n ( v α j , κ P κ + u α j , κ P ~ κ ) , v α j , κ , u α j , κ ∈ C . (43) W=\sum_{\alpha}g_{\alpha}F_{\alpha1}\cdots F_{\alpha\ell_{\alpha}},\qquad F_{\alpha j}=\sum_{\kappa=1}^{4n}(v_{\alpha j,\kappa}P_{\kappa}+u_{\alpha j,\kappa}\widetilde{P}_{\kappa}),\qquad v_{\alpha j,\kappa},u_{\alpha j,\kappa}\in\mathbb{C}. \tag*{(43)} W = α ∑ g α F α 1 ⋯ F α ℓ α , F α j = κ = 1 ∑ 4 n ( v α j , κ P κ + u α j , κ P κ ) , v α j , κ , u α j , κ ∈ C . (43) Here ℓ α \ell_{\alpha} ℓ α α \alpha α W W W
We introduce the exterior algebra on the doubled Fock space through the identities needed below; a general treatment is given in Ref. [92 ] . Let V = V ↑ ⊕ V ↓ = C 4 n \mathcal{V}=\mathcal{V}_{\uparrow}\oplus\mathcal{V}_{\downarrow}=\mathbb{C}^{4n} V = V ↑ ⊕ V ↓ = C 4 n V σ = C 2 n \mathcal{V}_{\sigma}=\mathbb{C}^{2n} V σ = C 2 n 1 , … , 4 n 1,\ldots,4n 1 , … , 4 n q κ q_{\kappa} q κ a ∈ ⋀ r V a\in\bigwedge^{r}\mathcal{V} a ∈ ⋀ r V
a i π ( 1 ) , … , i π ( r ) = sgn ( π ) a i 1 , … , i r , ∀ π ∈ S r , where i 1 , … , i r ∈ { 1 , … , 4 n } . (44) a_{i_{\pi(1)},\ldots,i_{\pi(r)}}=\operatorname{sgn}(\pi)a_{i_{1},\ldots,i_{r}},\qquad\forall\pi\in\mathcal{S}_{r},\quad\text{where }i_{1},\ldots,i_{r}\in\{1,\ldots,4n\}. \tag*{(44)} a i π ( 1 ) , … , i π ( r ) = sgn ( π ) a i 1 , … , i r , ∀ π ∈ S r , where i 1 , … , i r ∈ { 1 , … , 4 n } . (44) Here S r \mathcal{S}_{r} S r r r r sgn ( π ) \operatorname{sgn}(\pi) sgn ( π ) π \pi π r r r a ∈ ⋀ r V a\in\bigwedge^{r}\mathcal{V} a ∈ ⋀ r V
A r ( a ) = 1 r ! ∑ i 1 , … , i r = 1 4 n a i 1 , … , i r P ~ i 1 ⋯ P ~ i r = ∑ ∣ S ∣ = r a S τ S . (45) A_{r}(a)=\frac{1}{r!}\sum_{i_{1},\ldots,i_{r}=1}^{4n}a_{i_{1},\ldots,i_{r}}\widetilde{P}_{i_{1}}\cdots\widetilde{P}_{i_{r}}=\sum_{\lvert S\rvert=r}a_{S}\tau_{S}. \tag*{(45)} A r ( a ) = r ! 1 i 1 , … , i r = 1 ∑ 4 n a i 1 , … , i r P i 1 ⋯ P i r = ∣ S ∣ = r ∑ a S τ S . (45) The subscript r r r a ∈ ⋀ r V a\in\bigwedge^{r}\mathcal{V} a ∈ ⋀ r V b ∈ ⋀ q V b\in\bigwedge^{q}\mathcal{V} b ∈ ⋀ q V ( r + q ) (r+q) ( r + q ) [92 Sec. 3.4]
( a ∧ b ) i 1 , … , i r + q : = ∑ π ∈ S r + q , π ( 1 ) < ⋯ < π ( r ) , π ( r + 1 ) < ⋯ < π ( r + q ) sgn ( π ) a i π ( 1 ) , … , i π ( r ) b i π ( r + 1 ) , … , i π ( r + q ) . (46) (a\wedge b)_{i_{1},\ldots,i_{r+q}}:=\sum_{\substack{\pi\in\mathcal{S}_{r+q},\,\pi(1)<\cdots<\pi(r),\\ \pi(r+1)<\cdots<\pi(r+q)}}\operatorname{sgn}(\pi)a_{i_{\pi(1)},\ldots,i_{\pi(r)}}b_{i_{\pi(r+1)},\ldots,i_{\pi(r+q)}}. \tag*{(46)} ( a ∧ b ) i 1 , … , i r + q := π ∈ S r + q , π ( 1 ) < ⋯ < π ( r ) , π ( r + 1 ) < ⋯ < π ( r + q ) ∑ sgn ( π ) a i π ( 1 ) , … , i π ( r ) b i π ( r + 1 ) , … , i π ( r + q ) . (46) This definition gives the algebra homomorphism
A r ( a ) A q ( b ) = A r + q ( a ∧ b ) , ∀ a ∈ ⋀ r V , b ∈ ⋀ q V . (47) A_{r}(a)A_{q}(b)=A_{r+q}(a\wedge b),\qquad\forall a\in\bigwedge^{r}\mathcal{V},b\in\bigwedge^{q}\mathcal{V}. \tag*{(47)} A r ( a ) A q ( b ) = A r + q ( a ∧ b ) , ∀ a ∈ ⋀ r V , b ∈ ⋀ q V . (47) For a covector, i.e. linear functional v = ( v 1 , … , v 4 n ) ∈ V ∗ v=(v_{1},\ldots,v_{4n})\in\mathcal{V}^{*} v = ( v 1 , … , v 4 n ) ∈ V ∗ ι v : ⋀ r V → ⋀ r − 1 V \iota_{v}:\bigwedge^{r}\mathcal{V}\to\bigwedge^{r-1}\mathcal{V} ι v : ⋀ r V → ⋀ r − 1 V [92 Sec. 3.6]
( ι v a ) i 1 , … , i r − 1 : = ∑ j = 1 4 n v j a j , i 1 , … , i r − 1 , ∀ a ∈ ⋀ r V . (48) (\iota_{v}a)_{i_{1},\ldots,i_{r-1}}:=\sum_{j=1}^{4n}v_{j}a_{j,i_{1},\ldots,i_{r-1}},\qquad\forall a\in\bigwedge^{r}\mathcal{V}. \tag*{(48)} ( ι v a ) i 1 , … , i r − 1 := j = 1 ∑ 4 n v j a j , i 1 , … , i r − 1 , ∀ a ∈ ⋀ r V . (48) The corresponding operator identities are
Proposition 1. For v ∈ V ∗ v\in\mathcal{V}^{*} v ∈ V ∗ u ∈ V u\in\mathcal{V} u ∈ V
P ~ [ u ] : = ∑ κ = 1 4 n u κ P ~ κ , P [ v ] : = ∑ κ = 1 4 n v κ P κ . (49) \widetilde{P}[u]:=\sum_{\kappa=1}^{4n}u_{\kappa}\widetilde{P}_{\kappa},\qquad P[v]:=\sum_{\kappa=1}^{4n}v_{\kappa}P_{\kappa}. \tag*{(49)} P [ u ] := κ = 1 ∑ 4 n u κ P κ , P [ v ] := κ = 1 ∑ 4 n v κ P κ . (49) Then
P ~ [ u ] A r ( a ) = A r + 1 ( u ∧ a ) , P [ v ] A r ( a ) = A r − 1 ( ι v a ) + ( − 1 ) r A r ( a ) P [ v ] . (50) \widetilde{P}[u]A_r(a)=A_{r+1}(u\wedge a),\qquad P[v]A_r(a)=A_{r-1}(\iota_v a)+(-1)^r A_r(a)P[v]. \tag*{(50)} P [ u ] A r ( a ) = A r + 1 ( u ∧ a ) , P [ v ] A r ( a ) = A r − 1 ( ι v a ) + ( − 1 ) r A r ( a ) P [ v ] . (50) In particular,
P [ v ] A r ( a ) ∣ ρ 0 ⟩ ⟩ = A r − 1 ( ι v a ) ∣ ρ 0 ⟩ ⟩ . (51) P[v]A_r(a)|\rho_0\rangle\rangle=A_{r-1}(\iota_v a)|\rho_0\rangle\rangle. \tag*{(51)} P [ v ] A r ( a ) ∣ ρ 0 ⟩⟩ = A r − 1 ( ι v a ) ∣ ρ 0 ⟩⟩ . (51) These identities follow directly from the definitions of the wedge product and interior contraction. They give the action of the similarity-transformed normal modes on ∣ ρ 0 ⟩ ⟩ |\rho_0\rangle\rangle ∣ ρ 0 ⟩⟩
Proposition 2. Let F α j = ∑ κ = 1 4 n ( v α j , κ P κ + u α j , κ P ~ κ ) = P [ v α j ] + P ~ [ u α j ] F_{\alpha j}=\sum_{\kappa=1}^{4n}(v_{\alpha j,\kappa}P_\kappa+u_{\alpha j,\kappa}\widetilde{P}_\kappa)=P[v_{\alpha j}]+\widetilde{P}[u_{\alpha j}] F α j = ∑ κ = 1 4 n ( v α j , κ P κ + u α j , κ P κ ) = P [ v α j ] + P [ u α j ] (43) ). Denote
w α j ( 1 ) = u α j + ι v α j t ( 2 ) , w α j ( 3 ) = ι v α j t ( 4 ) , (52) w_{\alpha j}^{(1)}=u_{\alpha j}+\iota_{v_{\alpha j}}t^{(2)},\qquad w_{\alpha j}^{(3)}=\iota_{v_{\alpha j}}t^{(4)}, \tag*{(52)} w α j ( 1 ) = u α j + ι v α j t ( 2 ) , w α j ( 3 ) = ι v α j t ( 4 ) , (52) where t ( 2 ) ∈ ⋀ 2 V t^{(2)}\in\bigwedge^2\mathcal{V} t ( 2 ) ∈ ⋀ 2 V t ( 4 ) ∈ ⋀ 4 V t^{(4)}\in\bigwedge^4\mathcal{V} t ( 4 ) ∈ ⋀ 4 V F ‾ α j = e − T F α j e T \overline{F}_{\alpha j}=e^{-T}F_{\alpha j}e^T F α j = e − T F α j e T ∣ ρ 0 ⟩ ⟩ |\rho_0\rangle\rangle ∣ ρ 0 ⟩⟩
F ‾ α j A ( a ) ∣ ρ 0 ⟩ ⟩ = A ( φ ‾ α j [ a ] ) ∣ ρ 0 ⟩ ⟩ , where φ ‾ α j [ a ] = ι v α j a + ( w α j ( 1 ) + w α j ( 3 ) ) ∧ a . (53) \overline{F}_{\alpha j}A(a)|\rho_0\rangle\rangle=A(\overline{\varphi}_{\alpha j}[a])|\rho_0\rangle\rangle,\qquad\text{where }\overline{\varphi}_{\alpha j}[a]=\iota_{v_{\alpha j}}a+(w_{\alpha j}^{(1)}+w_{\alpha j}^{(3)})\wedge a. \tag*{(53)} F α j A ( a ) ∣ ρ 0 ⟩⟩ = A ( φ α j [ a ]) ∣ ρ 0 ⟩⟩ , where φ α j [ a ] = ι v α j a + ( w α j ( 1 ) + w α j ( 3 ) ) ∧ a . (53) Here A ( a ) = ∑ r A r ( a ( r ) ) A(a)=\sum_r A_r(a^{(r)}) A ( a ) = ∑ r A r ( a ( r ) ) a = ∑ r a ( r ) a=\sum_r a^{(r)} a = ∑ r a ( r ) a ( r ) ∈ ⋀ r V a^{(r)}\in\bigwedge^r\mathcal{V} a ( r ) ∈ ⋀ r V
Proof. We write the cluster operator as T = A 2 ( t ( 2 ) ) + A 4 ( t ( 4 ) ) T=A_2(t^{(2)})+A_4(t^{(4)}) T = A 2 ( t ( 2 ) ) + A 4 ( t ( 4 ) ) T T T
[ P [ v α j ] , T ] = A 1 ( ι v α j t ( 2 ) ) + A 3 ( ι v α j t ( 4 ) ) , [ P ~ [ u α j ] , T ] = 0. (54) [P[v_{\alpha j}],T]=A_1(\iota_{v_{\alpha j}}t^{(2)})+A_3(\iota_{v_{\alpha j}}t^{(4)}),\qquad[\widetilde{P}[u_{\alpha j}],T]=0. \tag*{(54)} [ P [ v α j ] , T ] = A 1 ( ι v α j t ( 2 ) ) + A 3 ( ι v α j t ( 4 ) ) , [ P [ u α j ] , T ] = 0. (54) The first commutator [ F α j , T ] [F_{\alpha j},T] [ F α j , T ] T T T (41) gives
F ‾ α j = F α j + [ F α j , T ] = P [ v α j ] + A ( w α j ( 1 ) + w α j ( 3 ) ) . (55) \overline{F}_{\alpha j}=F_{\alpha j}+[F_{\alpha j},T]=P[v_{\alpha j}]+A(w_{\alpha j}^{(1)}+w_{\alpha j}^{(3)}). \tag*{(55)} F α j = F α j + [ F α j , T ] = P [ v α j ] + A ( w α j ( 1 ) + w α j ( 3 ) ) . (55) Applying this identity to A ( a ) ∣ ρ 0 ⟩ ⟩ A(a)|\rho_0\rangle\rangle A ( a ) ∣ ρ 0 ⟩⟩ A A A
F ‾ α j A ( a ) ∣ ρ 0 ⟩ ⟩ = A ( ι v α j a + ( w α j ( 1 ) + w α j ( 3 ) ) ∧ a ) ∣ ρ 0 ⟩ ⟩ , (56) \overline{F}_{\alpha j}A(a)|\rho_0\rangle\rangle=A\left(\iota_{v_{\alpha j}}a+(w_{\alpha j}^{(1)}+w_{\alpha j}^{(3)})\wedge a\right)|\rho_0\rangle\rangle, \tag*{(56)} F α j A ( a ) ∣ ρ 0 ⟩⟩ = A ( ι v α j a + ( w α j ( 1 ) + w α j ( 3 ) ) ∧ a ) ∣ ρ 0 ⟩⟩ , (56) which proves Eq. (53) . □
The quadratic superoperator K 0 = ∑ κ = 1 4 n λ κ P ~ κ P κ K_0=\sum_{\kappa=1}^{4n}\lambda_\kappa\widetilde{P}_\kappa P_\kappa K 0 = ∑ κ = 1 4 n λ κ P κ P κ r r r
K 0 A r ( a ( r ) ) ∣ ρ 0 ⟩ ⟩ = : A r ( D λ a ( r ) ) ∣ ρ 0 ⟩ ⟩ , ( D λ a ( r ) ) i 1 , … , i r = ∑ j = 1 r λ i j ⋅ a i 1 , … , i r ( r ) . (57) K_0A_r(a^{(r)})|\rho_0\rangle\rangle=:A_r(D_\lambda a^{(r)})|\rho_0\rangle\rangle,\qquad(D_\lambda a^{(r)})_{i_1,\ldots,i_r}=\sum_{j=1}^{r}\lambda_{i_j}\cdot a^{(r)}_{i_1,\ldots,i_r}. \tag*{(57)} K 0 A r ( a ( r ) ) ∣ ρ 0 ⟩⟩ =: A r ( D λ a ( r ) ) ∣ ρ 0 ⟩⟩ , ( D λ a ( r ) ) i 1 , … , i r = j = 1 ∑ r λ i j ⋅ a i 1 , … , i r ( r ) . (57) For any ordered interaction term F α 1 ⋯ F α ℓ F_{\alpha_1}\cdots F_{\alpha_\ell} F α 1 ⋯ F α ℓ ϱ j \varrho_j ϱ j j = 1 , … , ℓ , ℓ + 1 j=1,\ldots,\ell,\ell+1 j = 1 , … , ℓ , ℓ + 1
F ‾ α j ⋯ F ‾ α ℓ ∣ ρ 0 ⟩ ⟩ = A ( ϱ α , j ) ∣ ρ 0 ⟩ ⟩ , j = 1 , … , ℓ , ϱ α , ℓ + 1 = 1. (58) \overline{F}_{\alpha_j}\cdots\overline{F}_{\alpha_\ell}|\rho_0\rangle\rangle=A(\varrho_{\alpha,j})|\rho_0\rangle\rangle,\qquad j=1,\ldots,\ell,\quad\varrho_{\alpha,\ell+1}=1. \tag*{(58)} F α j ⋯ F α ℓ ∣ ρ 0 ⟩⟩ = A ( ϱ α , j ) ∣ ρ 0 ⟩⟩ , j = 1 , … , ℓ , ϱ α , ℓ + 1 = 1. (58) By Proposition 2, these tensors satisfy the recursion
ϱ j ( r ) = ι v j ϱ j + 1 ( r + 1 ) + w j ( 1 ) ∧ ϱ j + 1 ( r − 1 ) + w j ( 3 ) ∧ ϱ j + 1 ( r − 3 ) , j = 1 , … , ℓ . (59) \varrho_j^{(r)}=\iota_{v_j}\varrho_{j+1}^{(r+1)}+w_j^{(1)}\wedge\varrho_{j+1}^{(r-1)}+w_j^{(3)}\wedge\varrho_{j+1}^{(r-3)},\qquad j=1,\ldots,\ell. \tag*{(59)} ϱ j ( r ) = ι v j ϱ j + 1 ( r + 1 ) + w j ( 1 ) ∧ ϱ j + 1 ( r − 1 ) + w j ( 3 ) ∧ ϱ j + 1 ( r − 3 ) , j = 1 , … , ℓ . (59) We suppress α \alpha α (59) and the next few equations. Components of negative degree vanish. For a quadratic term F α 1 F α 2 F_{\alpha_1}F_{\alpha_2} F α 1 F α 2 ℓ = 2 \ell=2 ℓ = 2 (59) gives
ϱ 2 ( 1 ) = w 2 ( 1 ) , ϱ 2 ( 3 ) = w 2 ( 3 ) , ϱ 1 ( 2 ) = ι v 1 w 2 ( 3 ) + w 1 ( 1 ) ∧ w 2 ( 1 ) , ϱ 1 ( 4 ) = w 1 ( 1 ) ∧ w 2 ( 3 ) + w 1 ( 3 ) ∧ w 2 ( 1 ) . (60) \varrho_2^{(1)}=w_2^{(1)},\qquad\varrho_2^{(3)}=w_2^{(3)},\qquad\varrho_1^{(2)}=\iota_{v_1}w_2^{(3)}+w_1^{(1)}\wedge w_2^{(1)},\qquad\varrho_1^{(4)}=w_1^{(1)}\wedge w_2^{(3)}+w_1^{(3)}\wedge w_2^{(1)}. \tag*{(60)} ϱ 2 ( 1 ) = w 2 ( 1 ) , ϱ 2 ( 3 ) = w 2 ( 3 ) , ϱ 1 ( 2 ) = ι v 1 w 2 ( 3 ) + w 1 ( 1 ) ∧ w 2 ( 1 ) , ϱ 1 ( 4 ) = w 1 ( 1 ) ∧ w 2 ( 3 ) + w 1 ( 3 ) ∧ w 2 ( 1 ) . (60) For a quartic interaction term F α 1 F α 2 F α 3 F α 4 F_{\alpha_1}F_{\alpha_2}F_{\alpha_3}F_{\alpha_4} F α 1 F α 2 F α 3 F α 4 ℓ = 4 \ell=4 ℓ = 4
ϱ 4 ( 1 ) = w 4 ( 1 ) , ϱ 4 ( 3 ) = w 4 ( 3 ) , ϱ 3 ( 0 ) = ι v 3 ϱ 4 ( 1 ) , ϱ 3 ( 2 ) = ι v 3 ϱ 4 ( 3 ) + w 3 ( 1 ) ∧ ϱ 4 ( 1 ) , ϱ 4 ( 4 ) = w 3 ( 1 ) ∧ ϱ 3 ( 3 ) + w 3 ( 3 ) ∧ ϱ 3 ( 1 ) , ϱ 3 ( 6 ) = w 3 ( 3 ) ∧ ϱ 3 ( 3 ) , ϱ 2 ( 1 ) = ι v 2 ϱ 3 ( 2 ) + w 2 ( 1 ) ∧ ϱ 3 ( 0 ) , ϱ 3 ( 3 ) = ι v 2 ϱ 3 ( 4 ) + w 2 ( 1 ) ∧ ϱ 3 ( 2 ) + w 2 ( 3 ) ∧ ϱ 3 ( 0 ) , ϱ 2 ( 5 ) = ι v 2 ϱ 3 ( 6 ) + w 2 ( 1 ) ∧ ϱ 3 ( 4 ) + w 2 ( 3 ) ∧ ϱ 3 ( 2 ) , ϱ 1 ( 2 ) = ι v 1 ϱ 2 ( 3 ) + w 1 ( 1 ) ∧ ϱ 2 ( 1 ) , ϱ 1 ( 4 ) = ι v 1 ϱ 2 ( 5 ) + w 1 ( 1 ) ∧ ϱ 2 ( 3 ) + w 1 ( 3 ) ∧ ϱ 2 ( 1 ) . (61) \begin{aligned}
\varrho_4^{(1)} &= w_4^{(1)}, & \varrho_4^{(3)} &= w_4^{(3)}, & \varrho_3^{(0)} &= \iota_{v_3}\varrho_4^{(1)}, & \varrho_3^{(2)} &= \iota_{v_3}\varrho_4^{(3)} + w_3^{(1)} \wedge\varrho_4^{(1)}, & \varrho_4^{(4)} &= w_3^{(1)} \wedge\varrho_3^{(3)} + w_3^{(3)} \wedge\varrho_3^{(1)}, \\
\varrho_3^{(6)} &= w_3^{(3)} \wedge\varrho_3^{(3)}, & \varrho_2^{(1)} &= \iota_{v_2}\varrho_3^{(2)} + w_2^{(1)} \wedge\varrho_3^{(0)}, & \varrho_3^{(3)} &= \iota_{v_2}\varrho_3^{(4)} + w_2^{(1)} \wedge\varrho_3^{(2)} + w_2^{(3)} \wedge\varrho_3^{(0)}, \\
\varrho_2^{(5)} &= \iota_{v_2}\varrho_3^{(6)} + w_2^{(1)} \wedge\varrho_3^{(4)} + w_2^{(3)} \wedge\varrho_3^{(2)}, & \varrho_1^{(2)} &= \iota_{v_1}\varrho_2^{(3)} + w_1^{(1)} \wedge\varrho_2^{(1)}, & \varrho_1^{(4)} &= \iota_{v_1}\varrho_2^{(5)} + w_1^{(1)} \wedge\varrho_2^{(3)} + w_1^{(3)} \wedge\varrho_2^{(1)}. \tag*{(61)}
\end{aligned} ϱ 4 ( 1 ) ϱ 3 ( 6 ) ϱ 2 ( 5 ) = w 4 ( 1 ) , = w 3 ( 3 ) ∧ ϱ 3 ( 3 ) , = ι v 2 ϱ 3 ( 6 ) + w 2 ( 1 ) ∧ ϱ 3 ( 4 ) + w 2 ( 3 ) ∧ ϱ 3 ( 2 ) , ϱ 4 ( 3 ) ϱ 2 ( 1 ) ϱ 1 ( 2 ) = w 4 ( 3 ) , = ι v 2 ϱ 3 ( 2 ) + w 2 ( 1 ) ∧ ϱ 3 ( 0 ) , = ι v 1 ϱ 2 ( 3 ) + w 1 ( 1 ) ∧ ϱ 2 ( 1 ) , ϱ 3 ( 0 ) ϱ 3 ( 3 ) ϱ 1 ( 4 ) = ι v 3 ϱ 4 ( 1 ) , = ι v 2 ϱ 3 ( 4 ) + w 2 ( 1 ) ∧ ϱ 3 ( 2 ) + w 2 ( 3 ) ∧ ϱ 3 ( 0 ) , = ι v 1 ϱ 2 ( 5 ) + w 1 ( 1 ) ∧ ϱ 2 ( 3 ) + w 1 ( 3 ) ∧ ϱ 2 ( 1 ) . ϱ 3 ( 2 ) = ι v 3 ϱ 4 ( 3 ) + w 3 ( 1 ) ∧ ϱ 4 ( 1 ) , ϱ 4 ( 4 ) = w 3 ( 1 ) ∧ ϱ 3 ( 3 ) + w 3 ( 3 ) ∧ ϱ 3 ( 1 ) , (61) Only ϱ 1 ( 2 ) \varrho_1^{(2)} ϱ 1 ( 2 ) ϱ 1 ( 4 ) \varrho_1^{(4)} ϱ 1 ( 4 )
∣ r ( t C C S D ) ⟩ ⟩ = A ( r ) ∣ ρ 0 ⟩ ⟩ , r ( 2 ) = D λ t ( 2 ) + ∑ α g α ϱ α , 1 ( 2 ) , r ( 4 ) = D λ t ( 4 ) + ∑ α g α ϱ α , 1 ( 4 ) . (62) \lvert r(t_{\mathrm{CCSD}})\rangle\rangle= A(r)\lvert\rho_0\rangle\rangle,\qquad r^{(2)} = D_{\lambda}t^{(2)} + \sum_{\alpha}g_{\alpha}\varrho_{\alpha,1}^{(2)},\qquad r^{(4)} = D_{\lambda}t^{(4)} + \sum_{\alpha}g_{\alpha}\varrho_{\alpha,1}^{(4)}. \tag*{(62)} ∣ r ( t CCSD )⟩⟩ = A ( r ) ∣ ρ 0 ⟩⟩ , r ( 2 ) = D λ t ( 2 ) + α ∑ g α ϱ α , 1 ( 2 ) , r ( 4 ) = D λ t ( 4 ) + α ∑ g α ϱ α , 1 ( 4 ) . (62) Projection onto { ⟨ ⟨ Φ ~ S ∣ : S ∈ S C C S D } \{\langle\langle\widetilde{\Phi}_S\lvert:S\in\mathcal{S}_{\mathrm{CCSD}}\} {⟨⟨ Φ S ∣ : S ∈ S CCSD }
( D λ t ( 2 ) + ∑ α g α ϱ α , 1 ( 2 ) ) ∣ Q = 0 = 0 , ( D λ t ( 4 ) + ∑ α g α ϱ α , 1 ( 4 ) ) ∣ Q = 0 = 0. (63) \left.\left(D_{\lambda}t^{(2)}+\sum_{\alpha}g_{\alpha}\varrho_{\alpha,1}^{(2)}\right)\right|_{Q=0}=0,\qquad\left.\left(D_{\lambda}t^{(4)}+\sum_{\alpha}g_{\alpha}\varrho_{\alpha,1}^{(4)}\right)\right|_{Q=0}=0. \tag*{(63)} ( D λ t ( 2 ) + α ∑ g α ϱ α , 1 ( 2 ) ) Q = 0 = 0 , ( D λ t ( 4 ) + α ∑ g α ϱ α , 1 ( 4 ) ) Q = 0 = 0. (63) This gives the exterior-algebra formulation of the CCSD equations.
Several optimizations could be applied to reduce the cost of these recursions. For instance, only contraction paths that can reach a charge-neutral degree-2 or degree-4 tensor need to be evaluated. An intermediate can carry nonzero charge, but it can be discarded if the remaining factors cannot return it to Q = 0 Q=0 Q = 0 a a a b b b ι v 1 ι v 2 ( a ∧ b ) \iota_{v_1}\iota_{v_2}(a\wedge b) ι v 1 ι v 2 ( a ∧ b ) [92 Sec. 3.6]
ι v 1 ι v 2 ( a ∧ b ) = ( ι v 1 ι v 2 a ) ∧ b + ( ι v 2 a ) ∧ ( ι v 1 b ) − ( ι v 1 a ) ∧ ( ι v 2 b ) + a ∧ ( ι v 1 ι v 2 b ) , (64) \iota_{v_1}\iota_{v_2}(a\wedge b)=(\iota_{v_1}\iota_{v_2}a)\wedge b+(\iota_{v_2}a)\wedge(\iota_{v_1}b)-(\iota_{v_1}a)\wedge(\iota_{v_2}b)+a\wedge(\iota_{v_1}\iota_{v_2}b), \tag*{(64)} ι v 1 ι v 2 ( a ∧ b ) = ( ι v 1 ι v 2 a ) ∧ b + ( ι v 2 a ) ∧ ( ι v 1 b ) − ( ι v 1 a ) ∧ ( ι v 2 b ) + a ∧ ( ι v 1 ι v 2 b ) , (64) so only tensors of degree at most four need to be formed. Moreover, paths that share an ordered prefix generate identical intermediate tensors, which need to be evaluated only once and can be reused across multiple paths. The optimized contractions can be represented as array operations and compiled into GPU kernels for efficient evaluation. The same contraction graph provides the Jacobian–vector products required by the Newton–Krylov solver described in Sec. S2 C.
Solving the CCSD equations with Newton–Krylov iterations
We solve the nonlinear system Eq. (S45) with a Newton–Krylov method, which has been shown to improve the convergence behavior of Hamiltonian CC equations [95 ] . By adapting the forcing term at each Newton step, the Eisenstat–Walker strategy [96 ] dynamically sets the Krylov tolerance to prevent oversolving. The required Jacobian–vector product (JVP) acts on the independent charge-neutral components of the degree-2 and degree-4 amplitude tensors. In practice, the implementation packs these complex amplitudes into a real vector. The derivation below only uses complex notation for clarity.
Differentiating Eq. (S34) with respect to the CCSD amplitudes gives
δ w α j ( 1 ) = ι v α j δ t ( 2 ) , δ w α j ( 3 ) = ι v α j δ t ( 4 ) . (65) \delta w_{\alpha j}^{(1)}=\iota_{v_{\alpha j}}\delta t^{(2)},\qquad\delta w_{\alpha j}^{(3)}=\iota_{v_{\alpha j}}\delta t^{(4)}. \tag*{(65)} δ w α j ( 1 ) = ι v α j δ t ( 2 ) , δ w α j ( 3 ) = ι v α j δ t ( 4 ) . (65) Applying the product rule to Eq. (S41) gives
δ ϱ j ( r ) = ι v j δ ϱ j + 1 ( r + 1 ) + δ w j ( 1 ) ∧ ϱ j + 1 ( r − 1 ) + w j ( 1 ) ∧ δ ϱ j + 1 ( r − 1 ) + δ w j ( 3 ) ∧ ϱ j + 1 ( r − 3 ) + w j ( 3 ) ∧ δ ϱ j + 1 ( r − 3 ) , j = 1 , … , ℓ , δ ϱ ℓ + 1 ( 0 ) = 0. (66) \delta\varrho_j^{(r)}=\iota_{v_j}\delta\varrho_{j+1}^{(r+1)}+\delta w_j^{(1)}\wedge\varrho_{j+1}^{(r-1)}+w_j^{(1)}\wedge\delta\varrho_{j+1}^{(r-1)}+\delta w_j^{(3)}\wedge\varrho_{j+1}^{(r-3)}+w_j^{(3)}\wedge\delta\varrho_{j+1}^{(r-3)},\quad j=1,\ldots,\ell,\quad\delta\varrho_{\ell+1}^{(0)}=0. \tag*{(66)} δ ϱ j ( r ) = ι v j δ ϱ j + 1 ( r + 1 ) + δ w j ( 1 ) ∧ ϱ j + 1 ( r − 1 ) + w j ( 1 ) ∧ δ ϱ j + 1 ( r − 1 ) + δ w j ( 3 ) ∧ ϱ j + 1 ( r − 3 ) + w j ( 3 ) ∧ δ ϱ j + 1 ( r − 3 ) , j = 1 , … , ℓ , δ ϱ ℓ + 1 ( 0 ) = 0. (66) The resulting JVP is
D z ( t C C S D ) ( δ t ) = ( D λ δ t + ∑ α g α δ ϱ α , 1 ) ∣ Q = 0 , ∣ S ∣ = 2 , 4 , z ( t C C S D ) = [ z S ( t C C S D ) ] S ∈ S C C S D . (67) \mathcal{D}_{\mathbf{z}}(t_{\mathrm{CCSD}})(\delta t)=\left.\left(D_{\lambda}\delta t+\sum_{\alpha}g_{\alpha}\delta\varrho_{\alpha,1}\right)\right|_{Q=0,\lvert S\rvert=2,4},\qquad\mathbf{z}(t_{\mathrm{CCSD}})=[z_S(t_{\mathrm{CCSD}})]_{S\in\mathcal{S}_{\mathrm{CCSD}}}. \tag*{(67)} D z ( t CCSD ) ( δ t ) = ( D λ δ t + α ∑ g α δ ϱ α , 1 ) Q = 0 , ∣ S ∣ = 2 , 4 , z ( t CCSD ) = [ z S ( t CCSD ) ] S ∈ S CCSD . (67) The Newton vector contains the same independent charge-neutral degree-2 and degree-4 components as t C C S D t_{\mathrm{CCSD}} t CCSD
D z [ t C C S D k ] ( δ t k ) = − z ( t C C S D k ) . (68) \mathcal{D}_{\mathbf{z}}[t_{\mathrm{CCSD}}^k](\delta t^k)=-\mathbf{z}(t_{\mathrm{CCSD}}^k). \tag*{(68)} D z [ t CCSD k ] ( δ t k ) = − z ( t CCSD k ) . (68) We left-precondition Newton–GMRES with D λ − 1 D_{\lambda}^{-1} D λ − 1 ( D λ ) S S = ∑ κ ∈ S λ κ (D_{\lambda})_{SS}=\sum_{\kappa\in S}\lambda_{\kappa} ( D λ ) S S = ∑ κ ∈ S λ κ
∥ G [ t C C S D k ] ( δ t k ) + z ( t C C S D k ) ∥ ≤ η k ∥ z ( t C C S D k ) ∥ , 0 ≤ η k < 1 , t C C S D k + 1 = t C C S D k + s k δ t k . (69) \left\|\mathcal{G}\left[t_{\mathrm{CCSD}}^{k}\right](\delta t^{k})+z\left(t_{\mathrm{CCSD}}^{k}\right)\right\|\leq\eta_{k}\left\|z\left(t_{\mathrm{CCSD}}^{k}\right)\right\|,\qquad0\leq\eta_{k}<1,\qquad t_{\mathrm{CCSD}}^{k+1}=t_{\mathrm{CCSD}}^{k}+s_{k}\delta t^{k}. \tag*{(69)} G [ t CCSD k ] ( δ t k ) + z ( t CCSD k ) ≤ η k z ( t CCSD k ) , 0 ≤ η k < 1 , t CCSD k + 1 = t CCSD k + s k δ t k . (69) The adaptive forcing term allows a loose inner tolerance when the nonlinear residual is large and tightens the GMRES solve near convergence. This avoids oversolving the early Newton steps while preserving rapid local convergence. A backtracking line search with an Armijo condition determines s k s_{k} s k [99 ] .
EOM-CC equations for Green’s functions
We derive the EOM-CC equations for the Lindbladian Green’s function [89 , 100 ] . For impurity sites i , j ∈ I i,j\in\mathcal{I} i , j ∈ I
G i j , σ R ( t ) = − i Θ ( t ) Tr [ ρ C C { c ^ i σ ( t ) , c ^ j σ † ( 0 ) } ] , c ^ i σ ( t ) = e t L † [ c ^ i σ ] (70) G_{ij,\sigma}^{R}(t)=-i\Theta(t)\operatorname{Tr}\left[\rho_{\mathrm{CC}}\left\{\hat{c}_{i\sigma}(t),\hat{c}_{j\sigma}^{\dagger}(0)\right\}\right],\qquad\hat{c}_{i\sigma}(t)=e^{t\mathcal{L}^{\dagger}}[\hat{c}_{i\sigma}] \tag*{(70)} G ij , σ R ( t ) = − i Θ ( t ) Tr [ ρ CC { c ^ iσ ( t ) , c ^ j σ † ( 0 ) } ] , c ^ iσ ( t ) = e t L † [ c ^ iσ ] (70) where the odd-parity Z 2 \mathbb{Z}_{2} Z 2 L ‾ [ ∙ ] : = Π ^ L [ Π ^ ⋅ ∙ ] \underline{\mathcal{L}}[\bullet]:=\hat{\Pi}\mathcal{L}[\hat{\Pi}\cdot\bullet] L [ ∙ ] := Π ^ L [ Π ^ ⋅ ∙ ]
L ‾ [ ∙ ] = − i [ H ^ C L , ∙ ] − ∑ i , j , σ Γ i j p ( 2 p ^ i σ † ∙ p ^ j σ + { p ^ j σ p ^ i σ † , ∙ } ) − ∑ i , j , σ Γ i j h ( 2 d ^ j σ ∙ d ^ i σ † + { d ^ i σ † d ^ j σ , ∙ } ) . (71) \underline{\mathcal{L}}[\bullet]=-i[\hat{H}_{\mathrm{CL}},\bullet]-\sum_{i,j,\sigma}\Gamma_{ij}^{p}\left(2\hat{p}_{i\sigma}^{\dagger}\bullet\hat{p}_{j\sigma}+\{\hat{p}_{j\sigma}\hat{p}_{i\sigma}^{\dagger},\bullet\}\right)-\sum_{i,j,\sigma}\Gamma_{ij}^{h}\left(2\hat{d}_{j\sigma}\bullet\hat{d}_{i\sigma}^{\dagger}+\{\hat{d}_{i\sigma}^{\dagger}\hat{d}_{j\sigma},\bullet\}\right). \tag*{(71)} L [ ∙ ] = − i [ H ^ CL , ∙ ] − i , j , σ ∑ Γ ij p ( 2 p ^ iσ † ∙ p ^ j σ + { p ^ j σ p ^ iσ † , ∙ } ) − i , j , σ ∑ Γ ij h ( 2 d ^ j σ ∙ d ^ iσ † + { d ^ iσ † d ^ j σ , ∙ } ) . (71) Therefore,
G i j , σ R ( t ) = − i Θ ( t ) Tr [ ρ C C { e t L † [ c ^ i σ ] , c ^ j σ † } ] = − i Θ ( t ) Tr [ c ^ i σ e t L ( c ^ j σ † ρ C C + ρ C C c ^ j σ † ) ] = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ i σ e − i t K ‾ P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ C C ⟩ ⟩ = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ i σ e − i t K ‾ P o ( c ^ j σ † − i c ~ j σ ) e T [ t C C ] ∣ ρ 0 ⟩ ⟩ = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ i σ e T [ t C C ] e − i t K ‾ P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ 0 ⟩ ⟩ . (72) \begin{aligned}
G_{ij,\sigma}^{R}(t)&=-i\Theta(t)\operatorname{Tr}\left[\rho_{\mathrm{CC}}\left\{e^{t\mathcal{L}^{\dagger}}[\hat{c}_{i\sigma}],\hat{c}_{j\sigma}^{\dagger}\right\}\right]\\
&=-i\Theta(t)\operatorname{Tr}\left[\hat{c}_{i\sigma}e^{t\mathcal{L}}\left(\hat{c}_{j\sigma}^{\dagger}\rho_{\mathrm{CC}}+\rho_{\mathrm{CC}}\hat{c}_{j\sigma}^{\dagger}\right)\right]\\
&=-i\Theta(t)\langle\!\langle I|\hat{c}_{i\sigma}e^{-it\overline{K}P_{\mathrm{o}}}\left(\hat{c}_{j\sigma}^{\dagger}-i\tilde{c}_{j\sigma}\right)|\rho_{\mathrm{CC}}\rangle\!\rangle\\
&=-i\Theta(t)\langle\!\langle I|\hat{c}_{i\sigma}e^{-it\overline{K}P_{\mathrm{o}}}\left(\hat{c}_{j\sigma}^{\dagger}-i\tilde{c}_{j\sigma}\right)e^{T[t_{\mathrm{CC}}]}|\rho_{0}\rangle\!\rangle\\
&=-i\Theta(t)\langle\!\langle I|\hat{c}_{i\sigma}e^{T[t_{\mathrm{CC}}]}e^{-it\overline{K}P_{\mathrm{o}}}\left(\hat{c}_{j\sigma}^{\dagger}-i\tilde{c}_{j\sigma}\right)|\rho_{0}\rangle\!\rangle. \tag*{(72)}
\end{aligned} G ij , σ R ( t ) = − i Θ ( t ) Tr [ ρ CC { e t L † [ c ^ iσ ] , c ^ j σ † } ] = − i Θ ( t ) Tr [ c ^ iσ e t L ( c ^ j σ † ρ CC + ρ CC c ^ j σ † ) ] = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ iσ e − i t K P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ CC ⟩ ⟩ = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ iσ e − i t K P o ( c ^ j σ † − i c ~ j σ ) e T [ t CC ] ∣ ρ 0 ⟩ ⟩ = − i Θ ( t ) ⟨ ⟨ I ∣ c ^ iσ e T [ t CC ] e − i t K P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ 0 ⟩ ⟩ . (72) Here P o = 1 2 ( I − ( − 1 ) N + N ~ ) P_{\mathrm{o}}=\frac{1}{2}(I-(-1)^{N+\tilde{N}}) P o = 2 1 ( I − ( − 1 ) N + N ~ ) L † \mathcal{L}^{\dagger} L † L \mathcal{L} L T T T P o P_{\mathrm{o}} P o
EOM-CC restricts the response to an odd-parity subspace with a specified excitation-rank cutoff. For the retarded Green’s function, we use the charge-1 sector with maximum rank 3 or 5:
S R ( σ ) = { S : ∣ S ∣ ∈ R , Q σ ( S ) = 1 , Q σ ˉ ( S ) = 0 } , R = { 1 , 3 } or { 1 , 3 , 5 } . (73) \mathcal{S}_{\mathcal{R}}^{(\sigma)}=\{S:|S|\in\mathcal{R},Q_{\sigma}(S)=1,Q_{\bar{\sigma}}(S)=0\},\qquad\mathcal{R}=\{1,3\}\ \text{or}\ \{1,3,5\}. \tag*{(73)} S R ( σ ) = { S : ∣ S ∣ ∈ R , Q σ ( S ) = 1 , Q σ ˉ ( S ) = 0 } , R = { 1 , 3 } or { 1 , 3 , 5 } . (73) Let P R ( σ ) P_{\mathcal{R}}^{(\sigma)} P R ( σ ) { ∣ Φ S ⟩ ⟩ : S ∈ S R ( σ ) } \{|\Phi_{S}\rangle\!\rangle:S\in\mathcal{S}_{\mathcal{R}}^{(\sigma)}\} { ∣ Φ S ⟩ ⟩ : S ∈ S R ( σ ) }
P R ( σ ) = ∑ S ∈ S R ( σ ) ∣ Φ S ⟩ ⟩ ⟨ ⟨ Φ ~ S ∣ . (74) P_{\mathcal{R}}^{(\sigma)}=\sum_{S\in\mathcal{S}_{\mathcal{R}}^{(\sigma)}}|\Phi_{S}\rangle\!\rangle\langle\!\langle\widetilde{\Phi}_{S}|. \tag*{(74)} P R ( σ ) = S ∈ S R ( σ ) ∑ ∣ Φ S ⟩ ⟩ ⟨ ⟨ Φ S ∣. (74) The EOM amplitudes x j σ = ∑ r ∈ R x j σ ( r ) x_{j\sigma}=\sum_{r\in\mathcal{R}}x_{j\sigma}^{(r)} x j σ = ∑ r ∈ R x j σ ( r ) x j σ ( r ) ∈ R V x_{j\sigma}^{(r)}\in\mathbb{R}^{\mathcal{V}} x j σ ( r ) ∈ R V A ( x j σ ) ∣ ρ 0 ⟩ ⟩ = ∑ r ∈ R A r ( x j σ ( r ) ) ∣ ρ 0 ⟩ ⟩ A(x_{j\sigma})|\rho_{0}\rangle\!\rangle=\sum_{r\in\mathcal{R}}A_{r}(x_{j\sigma}^{(r)})|\rho_{0}\rangle\!\rangle A ( x j σ ) ∣ ρ 0 ⟩ ⟩ = ∑ r ∈ R A r ( x j σ ( r ) ) ∣ ρ 0 ⟩ ⟩ e − i t K ‾ P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ 0 ⟩ ⟩ e^{-it\overline{K}P_{\mathrm{o}}}(\hat{c}_{j\sigma}^{\dagger}-i\tilde{c}_{j\sigma})|\rho_{0}\rangle\!\rangle e − i t K P o ( c ^ j σ † − i c ~ j σ ) ∣ ρ 0 ⟩ ⟩ A ( x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ A(x_{j\sigma}(t))|\rho_{0}\rangle\!\rangle A ( x j σ ( t )) ∣ ρ 0 ⟩ ⟩ { ⟨ ⟨ Φ ~ S ∣ : S ∈ S R ( σ ) } \{\langle\!\langle\widetilde{\Phi}_{S}|:S\in\mathcal{S}_{\mathcal{R}}^{(\sigma)}\} {⟨ ⟨ Φ S ∣ : S ∈ S R ( σ ) }
P R ( σ ) [ i A ( ∂ t x j σ ( t ) ) − K ‾ A ( x j σ ( t ) ) ] ∣ ρ 0 ⟩ ⟩ = 0. (75) P_{\mathcal{R}}^{(\sigma)}\left[iA(\partial_{t}x_{j\sigma}(t))-\overline{K}A(x_{j\sigma}(t))\right]|\rho_{0}\rangle\!\rangle=0. \tag*{(75)} P R ( σ ) [ i A ( ∂ t x j σ ( t )) − K A ( x j σ ( t )) ] ∣ ρ 0 ⟩ ⟩ = 0. (75) Since A ( x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ A(x_{j\sigma}(t))|\rho_{0}\rangle\!\rangle A ( x j σ ( t )) ∣ ρ 0 ⟩ ⟩
A ( ∂ t x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ = − i P R ( σ ) K ‾ A ( x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ = − i P R ( σ ) [ K ‾ , A ( x j σ ( t ) ) ] ∣ ρ 0 ⟩ ⟩ − i P R ( σ ) A ( x j σ ( t ) ) K ‾ ∣ ρ 0 ⟩ ⟩ . (76) A(\partial_{t}x_{j\sigma}(t))|\rho_{0}\rangle\!\rangle=-iP_{\mathcal{R}}^{(\sigma)}\overline{K}A(x_{j\sigma}(t))|\rho_{0}\rangle\!\rangle=-iP_{\mathcal{R}}^{(\sigma)}[\overline{K},A(x_{j\sigma}(t))]|\rho_{0}\rangle\!\rangle-iP_{\mathcal{R}}^{(\sigma)}A(x_{j\sigma}(t))\overline{K}|\rho_{0}\rangle\!\rangle. \tag*{(76)} A ( ∂ t x j σ ( t )) ∣ ρ 0 ⟩ ⟩ = − i P R ( σ ) K A ( x j σ ( t )) ∣ ρ 0 ⟩ ⟩ = − i P R ( σ ) [ K , A ( x j σ ( t ))] ∣ ρ 0 ⟩ ⟩ − i P R ( σ ) A ( x j σ ( t )) K ∣ ρ 0 ⟩ ⟩ . (76) The last term contains the static residual K ‾ ∣ ρ 0 ⟩ ⟩ \overline{K}|\rho_{0}\rangle\!\rangle K ∣ ρ 0 ⟩ ⟩
Multiplication by the pure-creator polynomial A ( x ) A(x) A ( x ) S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
A ( ∂ t x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ = − i P S R ( σ ) [ K ‾ , A ( x j σ ( t ) ) ] ∣ ρ 0 ⟩ ⟩ , A ( x j σ ( 0 ) ) ∣ ρ 0 ⟩ ⟩ = P S R ( σ ) ( c ‾ j σ † − i c ‾ j σ ) ∣ ρ 0 ⟩ ⟩ . (77) A\left(\partial_t x_{j\sigma}(t)\right)|\rho_0\rangle\rangle=-iP_{\mathcal{S}_{R}^{(\sigma)}}[\overline{K},A(x_{j\sigma}(t))]|\rho_0\rangle\rangle,\qquad A(x_{j\sigma}(0))|\rho_0\rangle\rangle=P_{\mathcal{S}_{R}^{(\sigma)}}(\overline{c}_{j\sigma}^{\dagger}-i\overline{c}_{j\sigma})|\rho_0\rangle\rangle. \tag*{(77)} A ( ∂ t x j σ ( t ) ) ∣ ρ 0 ⟩⟩ = − i P S R ( σ ) [ K , A ( x j σ ( t ))] ∣ ρ 0 ⟩⟩ , A ( x j σ ( 0 )) ∣ ρ 0 ⟩⟩ = P S R ( σ ) ( c j σ † − i c j σ ) ∣ ρ 0 ⟩⟩ . (77) Define the linear map A R , σ A_{R,\sigma} A R , σ b j σ R b_{j\sigma}^{R} b j σ R
[ A R , σ x ] S : = ⟨ ⟨ Φ ~ S ∣ [ K ‾ , A ( x ) ] ∣ ρ 0 ⟩ ⟩ , [ b j σ R ] S : = ⟨ ⟨ Φ ~ S ∣ ( c ‾ j σ † − i c ‾ j σ ) ∣ ρ 0 ⟩ ⟩ , S ∈ S R ( σ ) . (78) [A_{R,\sigma}x]_{S}:=\langle\langle\widetilde{\Phi}_{S}|[\overline{K},A(x)]|\rho_0\rangle\rangle,\qquad[b_{j\sigma}^{R}]_{S}:=\langle\langle\widetilde{\Phi}_{S}|(\overline{c}_{j\sigma}^{\dagger}-i\overline{c}_{j\sigma})|\rho_0\rangle\rangle,\qquad S\in\mathcal{S}_{R}^{(\sigma)}. \tag*{(78)} [ A R , σ x ] S := ⟨⟨ Φ S ∣ [ K , A ( x )] ∣ ρ 0 ⟩⟩ , [ b j σ R ] S := ⟨⟨ Φ S ∣ ( c j σ † − i c j σ ) ∣ ρ 0 ⟩⟩ , S ∈ S R ( σ ) . (78) Projecting Eq. (77) onto ⟨ ⟨ Φ ~ S ∣ \langle\langle\widetilde{\Phi}_{S}| ⟨⟨ Φ S ∣
∂ t x j σ ( t ) = − i A R , σ ( x j σ ( t ) ) , x j σ ( 0 ) = b j σ R , thus x j σ ( t ) = e − i t A R , σ b j σ R . (79) \partial_t x_{j\sigma}(t)=-iA_{R,\sigma}(x_{j\sigma}(t)),\qquad x_{j\sigma}(0)=b_{j\sigma}^{R},\qquad\text{thus}\quad x_{j\sigma}(t)=e^{-itA_{R,\sigma}}b_{j\sigma}^{R}. \tag*{(79)} ∂ t x j σ ( t ) = − i A R , σ ( x j σ ( t )) , x j σ ( 0 ) = b j σ R , thus x j σ ( t ) = e − i t A R , σ b j σ R . (79) Reconstruction of the Green’s function also requires the row vectors
[ s i σ R ] S : = ⟨ ⟨ I ∣ c ‾ i σ ∣ Φ S ⟩ ⟩ , S ∈ S R ( σ ) , i ∈ I . (80) [s_{i\sigma}^{R}]_{S}:=\langle\langle I|\overline{c}_{i\sigma}|\Phi_{S}\rangle\rangle,\qquad S\in\mathcal{S}_{R}^{(\sigma)},\quad i\in\mathcal{I}. \tag*{(80)} [ s iσ R ] S := ⟨⟨ I ∣ c iσ ∣ Φ S ⟩⟩ , S ∈ S R ( σ ) , i ∈ I . (80) The EOM-CC retarded Green’s function is then
G ~ i j , σ R , R ( t ) : = − i Θ ( t ) ⟨ ⟨ I ∣ c i σ e T A ( x j σ ( t ) ) ∣ ρ 0 ⟩ ⟩ = − i Θ ( t ) ∑ S ∈ S R ( σ ) ⟨ ⟨ I ∣ c i σ e T ∣ Φ S ⟩ ⟩ [ x j σ ( t ) ] S = − i Θ ( t ) ∑ S ∈ S R ( σ ) ⟨ ⟨ I ∣ c i σ ∣ Φ S ⟩ ⟩ [ x j σ ( t ) ] S = − i Θ ( t ) ( s i σ R ) T e − i t A R , σ b j σ R . (81) \begin{aligned}
\widetilde{G}_{ij,\sigma}^{R,R}(t):=-i\Theta(t)\langle\langle I|c_{i\sigma}e^{T}A(x_{j\sigma}(t))|\rho_0\rangle\rangle=-i\Theta(t)\sum_{S\in\mathcal{S}_{R}^{(\sigma)}}\langle\langle I|c_{i\sigma}e^{T}|\Phi_{S}\rangle\rangle[x_{j\sigma}(t)]_{S} \\
&=-i\Theta(t)\sum_{S\in\mathcal{S}_{R}^{(\sigma)}}\langle\langle I|c_{i\sigma}|\Phi_{S}\rangle\rangle[x_{j\sigma}(t)]_{S}=-i\Theta(t)(s_{i\sigma}^{R})^{\mathsf{T}}e^{-itA_{R,\sigma}}b_{j\sigma}^{R}. \tag*{(81)}
\end{aligned} G ij , σ R , R ( t ) := − i Θ ( t ) ⟨⟨ I ∣ c iσ e T A ( x j σ ( t )) ∣ ρ 0 ⟩⟩ = − i Θ ( t ) S ∈ S R ( σ ) ∑ ⟨⟨ I ∣ c iσ e T ∣ Φ S ⟩⟩ [ x j σ ( t ) ] S = − i Θ ( t ) S ∈ S R ( σ ) ∑ ⟨⟨ I ∣ c iσ ∣ Φ S ⟩⟩ [ x j σ ( t ) ] S = − i Θ ( t ) ( s iσ R ) T e − i t A R , σ b j σ R . (81) To obtain the penultimate equality, we apply Proposition 2 to the linear superoperator c ‾ i σ \overline{c}_{i\sigma} c iσ c ‾ i σ = c i σ + [ c i σ , T ] \overline{c}_{i\sigma}=c_{i\sigma}+[c_{i\sigma},T] c iσ = c iσ + [ c iσ , T ] ⟨ ⟨ I ∣ \langle\langle I| ⟨⟨ I ∣ ⟨ ⟨ I ∣ e T = ⟨ ⟨ I ∣ \langle\langle I|e^{T}=\langle\langle I| ⟨⟨ I ∣ e T = ⟨⟨ I ∣
⟨ ⟨ I ∣ c i σ e T = ⟨ ⟨ I ∣ c ‾ i σ = ⟨ ⟨ I ∣ c i σ . (82) \langle\langle I|c_{i\sigma}e^{T}=\langle\langle I|\overline{c}_{i\sigma}=\langle\langle I|c_{i\sigma}. \tag*{(82)} ⟨⟨ I ∣ c iσ e T = ⟨⟨ I ∣ c iσ = ⟨⟨ I ∣ c iσ . (82) A Fourier transform of Eq. (81) gives the frequency-domain EOM-CC Green’s function
G ~ i j , σ R , R ( z ) = ∫ − ∞ + ∞ G i j , σ R , R ( t ) e i z t d t = ( s i σ R ) T ( z I − A R , σ ) − 1 b j σ R , Im z > max Im spec ( A R , σ ) . (83) \widetilde{G}_{ij,\sigma}^{R,R}(z)=\int_{-\infty}^{+\infty}G_{ij,\sigma}^{R,R}(t)e^{izt}\,\mathrm{d}t=(s_{i\sigma}^{R})^{\mathsf{T}}(zI-A_{R,\sigma})^{-1}b_{j\sigma}^{R},\qquad\operatorname{Im}z>\max\operatorname{Im}\operatorname{spec}(A_{R,\sigma}). \tag*{(83)} G ij , σ R , R ( z ) = ∫ − ∞ + ∞ G ij , σ R , R ( t ) e i z t d t = ( s iσ R ) T ( z I − A R , σ ) − 1 b j σ R , Im z > max Im spec ( A R , σ ) . (83) Real-frequency values are evaluated from the resolvent with the stated broadening. We evaluate A R , σ x j σ A_{R,\sigma}x_{j\sigma} A R , σ x j σ A R , σ A_{R,\sigma} A R , σ F ‾ α 1 ⋯ F ‾ α ℓ \overline{F}_{\alpha_1}\cdots\overline{F}_{\alpha_\ell} F α 1 ⋯ F α ℓ ς α , j [ x ] \varsigma_{\alpha,j}[x] ς α , j [ x ] j = 1 , … , ℓ + 1 j=1,\ldots,\ell+1 j = 1 , … , ℓ + 1
F ‾ α j ⋯ F ‾ α ℓ A ( x ) ∣ ρ 0 ⟩ ⟩ = A ( ς α , j [ x ] ) ∣ ρ 0 ⟩ ⟩ , j = 1 , … , ℓ , ς α , ℓ + 1 [ x ] = x . (84) \overline{F}_{\alpha_j}\cdots\overline{F}_{\alpha_\ell}A(x)|\rho_0\rangle\rangle=A(\varsigma_{\alpha,j}[x])|\rho_0\rangle\rangle,\qquad j=1,\ldots,\ell,\qquad\varsigma_{\alpha,\ell+1}[x]=x. \tag*{(84)} F α j ⋯ F α ℓ A ( x ) ∣ ρ 0 ⟩⟩ = A ( ς α , j [ x ]) ∣ ρ 0 ⟩⟩ , j = 1 , … , ℓ , ς α , ℓ + 1 [ x ] = x . (84) The corresponding recursion is
ς j ( r ) [ x ] = ι v j ς j + 1 ( r + 1 ) [ x ] + w j ( 1 ) ∧ ς j + 1 ( r − 1 ) [ x ] + w j ( 3 ) ∧ ς j + 1 ( r − 3 ) [ x ] , j = 1 , … , ℓ . (85) \varsigma_{j}^{(r)}[x]=\iota_{v_j}\varsigma_{j+1}^{(r+1)}[x]+w_{j}^{(1)}\wedge\varsigma_{j+1}^{(r-1)}[x]+w_{j}^{(3)}\wedge\varsigma_{j+1}^{(r-3)}[x],\qquad j=1,\ldots,\ell. \tag*{(85)} ς j ( r ) [ x ] = ι v j ς j + 1 ( r + 1 ) [ x ] + w j ( 1 ) ∧ ς j + 1 ( r − 1 ) [ x ] + w j ( 3 ) ∧ ς j + 1 ( r − 3 ) [ x ] , j = 1 , … , ℓ . (85) This recursion has the form of Eq. (59) , but is initialized with the odd EOM tensor ς ℓ + 1 [ x ] = x \varsigma_{\ell+1}[x]=x ς ℓ + 1 [ x ] = x ∣ ρ 0 ⟩ ⟩ |\rho_0\rangle\rangle ∣ ρ 0 ⟩⟩
[ F ‾ α 1 ⋯ F ‾ α ℓ , A ( x ) ] ∣ ρ 0 ⟩ ⟩ = A ( ς α , 1 [ x ] ) ∣ ρ 0 ⟩ ⟩ − A ( x ) F ‾ α 1 ⋯ F ‾ α ℓ ∣ ρ 0 ⟩ ⟩ = A ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ∣ ρ 0 ⟩ ⟩ . (86) [\overline{F}_{\alpha_1}\cdots\overline{F}_{\alpha_\ell},A(x)]|\rho_0\rangle\rangle=A(\varsigma_{\alpha,1}[x])|\rho_0\rangle\rangle-A(x)\overline{F}_{\alpha_1}\cdots\overline{F}_{\alpha_\ell}|\rho_0\rangle\rangle=A(\varsigma_{\alpha,1}[x]-x\wedge\varrho_{\alpha,1})|\rho_0\rangle\rangle. \tag*{(86)} [ F α 1 ⋯ F α ℓ , A ( x )] ∣ ρ 0 ⟩⟩ = A ( ς α , 1 [ x ]) ∣ ρ 0 ⟩⟩ − A ( x ) F α 1 ⋯ F α ℓ ∣ ρ 0 ⟩⟩ = A ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ∣ ρ 0 ⟩⟩ . (86) Therefore, for K = ∑ κ = 1 4 n λ κ P ~ κ P κ + ∑ α g α F ‾ α 1 ⋯ F ‾ α ℓ K=\sum_{\kappa=1}^{4n}\lambda_{\kappa}\widetilde{P}_{\kappa}P_{\kappa}+\sum_{\alpha}g_{\alpha}\overline{F}_{\alpha_1}\cdots\overline{F}_{\alpha_\ell} K = ∑ κ = 1 4 n λ κ P κ P κ + ∑ α g α F α 1 ⋯ F α ℓ
[ K ‾ 0 , A ( x ) ] ∣ ρ 0 ⟩ ⟩ = [ K ‾ 0 + A ( D λ t ) , A ( x ) ] ∣ ρ 0 ⟩ ⟩ = A ( D λ x + ( D λ t ) ∧ x − x ∧ ( D λ t ) ) ∣ ρ 0 ⟩ ⟩ = A ( D λ x ) ∣ ρ 0 ⟩ ⟩ . (87) [\overline{K}_0,A(x)]|\rho_0\rangle\rangle=[\overline{K}_0+A(D_{\lambda}t),A(x)]|\rho_0\rangle\rangle=A(D_{\lambda}x+(D_{\lambda}t)\wedge x-x\wedge(D_{\lambda}t))|\rho_0\rangle\rangle=A(D_{\lambda}x)|\rho_0\rangle\rangle. \tag*{(87)} [ K 0 , A ( x )] ∣ ρ 0 ⟩⟩ = [ K 0 + A ( D λ t ) , A ( x )] ∣ ρ 0 ⟩⟩ = A ( D λ x + ( D λ t ) ∧ x − x ∧ ( D λ t )) ∣ ρ 0 ⟩⟩ = A ( D λ x ) ∣ ρ 0 ⟩⟩ . (87) Since x x x D λ t D_{\lambda}t D λ t x ∧ ( D λ t ) = ( D λ t ) ∧ x x\wedge(D_{\lambda}t)=(D_{\lambda}t)\wedge x x ∧ ( D λ t ) = ( D λ t ) ∧ x [92 Sec. 3.4] . The projected commutator can therefore be evaluated as
P S R ( σ ) [ K ‾ , A ( x ) ] ∣ ρ 0 ⟩ ⟩ = A [ ( D λ x + ∑ α g α ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ) ∣ S ∈ S R ( σ ) ] ∣ ρ 0 ⟩ ⟩ . (88) P_{\mathcal{S}_{R}^{(\sigma)}}[\overline{K},A(x)]|\rho_0\rangle\rangle=A\left[\left(D_{\lambda}x+\sum_{\alpha}g_{\alpha}\left(\varsigma_{\alpha,1}[x]-x\wedge\varrho_{\alpha,1}\right)\right)\bigg|_{S\in\mathcal{S}_{R}^{(\sigma)}}\right]|\rho_0\rangle\rangle. \tag*{(88)} P S R ( σ ) [ K , A ( x )] ∣ ρ 0 ⟩⟩ = A [ ( D λ x + α ∑ g α ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ) S ∈ S R ( σ ) ] ∣ ρ 0 ⟩⟩ . (88) In coefficient space, Eq. (88) gives the linear EOM map
( A R , σ x ) S ∈ S R ( σ ) = ( D λ x + ∑ α g α ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ) S ∈ S R ( σ ) . (89) \left(A_{R,\sigma}x\right)_{S\in\mathcal{S}_{R}^{(\sigma)}}=\left(D_{\lambda}x+\sum_{\alpha}g_{\alpha}\left(\varsigma_{\alpha,1}[x]-x\wedge\varrho_{\alpha,1}\right)\right)_{S\in\mathcal{S}_{R}^{(\sigma)}}. \tag*{(89)} ( A R , σ x ) S ∈ S R ( σ ) = ( D λ x + α ∑ g α ( ς α , 1 [ x ] − x ∧ ϱ α , 1 ) ) S ∈ S R ( σ ) . (89) These contractions can be optimized using the same strategies as for the CCSD residual. The similarity transformation of the linear odd operator c j σ † − i c ~ j σ c_{j\sigma}^{\dagger}-\mathrm{i}\tilde{c}_{j\sigma} c j σ † − i c ~ j σ T T T ∣ ρ 0 ⟩ ⟩ \lvert\rho_{0}\rangle\rangle ∣ ρ 0 ⟩⟩ ⟨ ⟨ I ∣ c i σ \langle\langle I\lvert c_{i\sigma} ⟨⟨ I ∣ c iσ b j σ R b_{j\sigma}^{R} b j σ R s i σ R s_{i\sigma}^{R} s iσ R
b j σ 1 , 3 , 5 = ( b j σ 1 , 3 0 ) , ( s i σ 1 , 3 , 5 ) T = ( ( s i σ 1 , 3 ) T 0 ) . (90) b_{j\sigma}^{1,3,5}=\begin{pmatrix}b_{j\sigma}^{1,3}\\0\end{pmatrix},\qquad\left(s_{i\sigma}^{1,3,5}\right)^{\mathsf{T}}=\left(\left(s_{i\sigma}^{1,3}\right)^{\mathsf{T}}\ 0\right). \tag*{(90)} b j σ 1 , 3 , 5 = ( b j σ 1 , 3 0 ) , ( s iσ 1 , 3 , 5 ) T = ( ( s iσ 1 , 3 ) T 0 ) . (90) To solve ( z I − A R , σ ) x j σ = b j σ R \left(zI-A_{R,\sigma}\right)x_{j\sigma}=b_{j\sigma}^{R} ( z I − A R , σ ) x j σ = b j σ R m K m_{K} m K K m K ( z I − A R , σ , b j σ R ) \mathcal{K}_{m_{K}}\left(zI-A_{R,\sigma},b_{j\sigma}^{R}\right) K m K ( z I − A R , σ , b j σ R ) K m K ( z I − A R , σ , b j σ R ) = K m K ( A R , σ , b j σ R ) \mathcal{K}_{m_{K}}\left(zI-A_{R,\sigma},b_{j\sigma}^{R}\right)=\mathcal{K}_{m_{K}}\left(A_{R,\sigma},b_{j\sigma}^{R}\right) K m K ( z I − A R , σ , b j σ R ) = K m K ( A R , σ , b j σ R ) z ∈ C z\in\mathbb{C} z ∈ C j σ j\sigma j σ z z z i σ i\sigma iσ G i j , σ R , R ( z ) G_{ij,\sigma}^{R,R}(z) G ij , σ R , R ( z ) ( z I − D λ ) − 1 \left(zI-D_{\lambda}\right)^{-1} ( z I − D λ ) − 1
Numerical setups and additional results
Here we specify the models and numerical settings for the calculations in the main text and present additional results. Unless otherwise stated, bath fits were obtained with the RealTimeBath package [105 ] , using the semidefinite-programming (SDP) construction of Ref. [106 ] . For fermionic baths, the fitted lesser and greater hybridization functions Δ ≶ ( t ) \Delta^{\lessgtr}(t) Δ ≶ ( t ) J ( ω ) = [ J i j ( ω ) ] i , j = 1 n I J(\omega)=[J_{ij}(\omega)]_{i,j=1}^{n_{I}} J ( ω ) = [ J ij ( ω ) ] i , j = 1 n I
Δ ≶ ( t ) = ∫ − ∞ + ∞ J ( ω ) f ≶ ( ω ) e − i ω t d ω . (91) \Delta^{\lessgtr}(t)=\int_{-\infty}^{+\infty}J(\omega)f^{\lessgtr}(\omega)e^{-\mathrm{i}\omega t}\,\mathrm{d}\omega. \tag*{(91)} Δ ≶ ( t ) = ∫ − ∞ + ∞ J ( ω ) f ≶ ( ω ) e − i ω t d ω . (91) Here, f < ( ω ) = f F D β , μ ( ω ) f^{<}(\omega)=f_{\mathrm{FD}}^{\beta,\mu}(\omega) f < ( ω ) = f FD β , μ ( ω ) f > ( ω ) = 1 − f F D β , μ ( ω ) f^{>}(\omega)=1-f_{\mathrm{FD}}^{\beta,\mu}(\omega) f > ( ω ) = 1 − f FD β , μ ( ω ) f F D β , μ ( ω ) = 1 e β ( ω − μ ) + 1 f_{\mathrm{FD}}^{\beta,\mu}(\omega)=\frac{1}{e^{\beta(\omega-\mu)}+1} f FD β , μ ( ω ) = e β ( ω − μ ) + 1 1 β \beta β μ \mu μ r r r (91) is generalized by replacing J ( ω ) f ≶ ( ω ) J(\omega)f^{\lessgtr}(\omega) J ( ω ) f ≶ ( ω ) ∑ r J r ( ω ) f r ≶ ( ω ) \sum_{r}J_{r}(\omega)f_{r}^{\lessgtr}(\omega) ∑ r J r ( ω ) f r ≶ ( ω ) f F D β r , μ r ( ω ) f_{\mathrm{FD}}^{\beta_{r},\mu_{r}}(\omega) f FD β r , μ r ( ω )
Single- and two-impurity Anderson models. —The particle–hole-symmetric single-impurity Anderson model is given by
H ^ I = − ϵ ∑ σ n ^ σ + U n ^ ↑ n ^ ↓ , ϵ = U / 2. (92) \hat{H}_{I}=-\epsilon\sum_{\sigma}\hat{n}_{\sigma}+U\hat{n}_{\uparrow}\hat{n}_{\downarrow},\qquad\epsilon=U/2. \tag*{(92)} H ^ I = − ϵ σ ∑ n ^ σ + U n ^ ↑ n ^ ↓ , ϵ = U /2. (92) The impurity couples to a semicircular bath with spectral density J ( ω ) = 1 π 1 − ω 2 W 2 ∣ 1 ( − W , W ) ( ω ) J(\omega)=\frac{1}{\pi}\sqrt{1-\frac{\omega^{2}}{W^{2}}\vphantom{\big|}}\mathbf{1}_{(-W,W)}(\omega) J ( ω ) = π 1 1 − W 2 ω 2 1 ( − W , W ) ( ω ) W = 10 W=10 W = 10 N p m = 2 N_{\mathrm{pm}}=2 N pm = 2 Figure 2 a) and N p m = 6 N_{\mathrm{pm}}=6 N pm = 6 Figure S2 ) coupled pseudomodes per lesser or greater component. For N p m = 2 N_{\mathrm{pm}}=2 N pm = 2 β = 100 \beta=100 β = 100 3 ^{3} 3 η = 0 \eta=0 η = 0 U = 2 U=2 U = 2 10 10 10 m K m_{K} m K S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5 S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5 U = 10 U=10 U = 10 β = 1 \beta=1 β = 1 N p m = 6 N_{\mathrm{pm}}=6 N pm = 6
Figure S2. The single-impurity Im G R ( t ) \operatorname{Im}G^{R}(t) Im G R ( t ) β = ∞ \beta=\infty β = ∞ U = 2 U=2 U = 2 2.5 2.5 2.5 3 3 3 N p m = 6 N_{\mathrm{pm}}=6 N pm = 6 M t o t = 48 M_{\mathrm{tot}}=48 M tot = 48
The N p m = 2 N_{\mathrm{pm}}=2 N pm = 2 U = 10 U=10 U = 10 β = 100 \beta=100 β = 100 [101 ] . The benchmark with ED therefore assesses solver accuracy for this small pseudomode model, without establishing its low-energy continuum physics quantitatively. Figure S1 compares the spectral functions and self-energies as the impurity level is shifted away from half filling. Across the fillings considered, the S 1 , 3 , 5 \mathcal{S}_{1,3,5} S 1 , 3 , 5 Figure S5 , which compares different pseudomode numbers for the low-temperature two-impurity bath, including the refined N p m = 6 N_{\mathrm{pm}}=6 N pm = 6
Figure S1. Comparison of C³L with S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 A ↑ ( ω ) A_{\uparrow}(\omega) A ↑ ( ω ) Σ ~ ↑ R ( ω ) \widetilde{\Sigma}^{R}_{\uparrow}(\omega) Σ ↑ R ( ω )
Figure S5. Bath representations for the strongly correlated two-impurity benchmark at β = 100 \beta= 100 β = 100 μ = 0 \mu= 0 μ = 0 W = 10 W = 10 W = 10 f F D ( ω ) f_{\mathrm{FD}}(\omega) f FD ( ω ) J f i t o c c ( ω ) J f i t o c c ( ω ) + J f i t e m p t y ( ω ) \frac{J_{\mathrm{fit}}^{\mathrm{occ}}(\omega)}{J_{\mathrm{fit}}^{\mathrm{occ}}(\omega)+J_{\mathrm{fit}}^{\mathrm{empty}}(\omega)} J fit occ ( ω ) + J fit empty ( ω ) J fit occ ( ω ) N p m = 2 , 4 , 6 N_{\mathrm{pm}} = 2, 4, 6 N pm = 2 , 4 , 6 N p m = 6 N_{\mathrm{pm}} = 6 N pm = 6 Figure 5 . The frequency-focused refinement improves Fermi-edge resolution while increasing the time-window relative L 2 L^2 L 2
Coupled-Lindblad pseudomodes can represent bath correlations over longer times with fewer modes [106 , 107 ] . To demonstrate this, we compare N p m = 6 N_{\mathrm{pm}}=6 N pm = 6 M = 48 M=48 M = 48 [108 ] . For a flat spectral density, this coincides with the Gauss-quadrature discretization underlying TEDOPA [109 ] . For the semicircular baths used here, both give polynomially convergent unitary representations. The comparison is performed at β = ∞ \beta=\infty β = ∞ ∣ ψ U ⟩ |\psi_U\rangle ∣ ψ U ⟩ [110 ] .
[ G U R ( t ) ] p q = − i Θ ( t ) [ ⟨ ψ U ∣ a ^ p e − i t ( H ^ U − E U ) a ^ q † ∣ ψ U ⟩ + ⟨ ψ U ∣ a ^ q † e i t ( H ^ U − E U ) a ^ p ∣ ψ U ⟩ ] . (93) \left[G^{R}_{U}(t)\right]_{pq}=-i\Theta(t)\left[\langle\psi_U|\hat{a}_{p}e^{-it(\hat{H}_{U}-E_{U})}\hat{a}^{\dagger}_{q}|\psi_U\rangle+\langle\psi_U|\hat{a}^{\dagger}_{q}e^{it(\hat{H}_{U}-E_{U})}\hat{a}_{p}|\psi_U\rangle\right]. \tag*{(93)} [ G U R ( t ) ] pq = − i Θ ( t ) [ ⟨ ψ U ∣ a ^ p e − i t ( H ^ U − E U ) a ^ q † ∣ ψ U ⟩ + ⟨ ψ U ∣ a ^ q † e i t ( H ^ U − E U ) a ^ p ∣ ψ U ⟩ ] . (93) Here H ^ U \hat{H}_{U} H ^ U E U E_{U} E U p , q p,q p , q U = 2 U=2 U = 2 2.5 2.5 2.5 3 3 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 G U R ( t ) G^{R}_{U}(t) G U R ( t )
For the two-impurity Anderson model
H ^ I = − t 12 ∑ σ ( a ^ 1 σ † a ^ 2 σ + a ^ 2 σ † a ^ 1 σ ) + U ∑ i = 1 , 2 n ^ i ↑ n ^ i ↓ , (94) \hat{H}_{I}=-t_{12}\sum_{\sigma}\left(\hat{a}^{\dagger}_{1\sigma}\hat{a}_{2\sigma}+\hat{a}^{\dagger}_{2\sigma}\hat{a}_{1\sigma}\right)+U\sum_{i=1,2}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}, \tag*{(94)} H ^ I = − t 12 σ ∑ ( a ^ 1 σ † a ^ 2 σ + a ^ 2 σ † a ^ 1 σ ) + U i = 1 , 2 ∑ n ^ i ↑ n ^ i ↓ , (94) we set t 12 = 0.5 t_{12}=0.5 t 12 = 0.5 U = 2 U=2 U = 2 μ = U / 2 \mu=U/2 μ = U /2 β = ∞ \beta=\infty β = ∞ J ( ω ) J(\omega) J ( ω ) W = 10 W=10 W = 10 N p m = 4 N_{\mathrm{pm}}=4 N pm = 4 2 n = 2 ( n I + 2 d N p m ) = 36 2n=2(n_I+2dN_{\mathrm{pm}})=36 2 n = 2 ( n I + 2 d N pm ) = 36 M = 24 M=24 M = 24 36 36 36 52 52 52 76 76 76 Figure S3 . The M t o t = 48 M_{\mathrm{tot}} = 48 M tot = 48 t = 1.5 t = 1.5 t = 1.5 Figure 2 c.
Figure S3. Lesser hybridization function Δ < ( t ) \Delta^{<}(t) Δ < ( t ) β = ∞ \beta= \infty β = ∞ M t o t = 48 , 72 M_{\mathrm{tot}} = 48, 72 M tot = 48 , 72 N p m = 4 N_{\mathrm{pm}} = 4 N pm = 4 ∣ δ Δ < ( t ) ∣ = ∣ Δ f i t < ( t ) − Δ < ( t ) ∣ \lvert\delta\Delta^{<}(t)\rvert= \lvert\Delta_{\mathrm{fit}}^{<}(t) - \Delta^{<}(t)\rvert ∣ δ Δ < ( t )∣ = ∣ Δ fit < ( t ) − Δ < ( t )∣
*Imaginary-time reconstruction.—*At equilibrium, we evaluate the EOM-CC resolvent in Eq. (83) at the Matsubara frequencies { i ω n } \{i\omega_n\} { i ω n }
G p q M ( τ ) = 1 β ∑ n ∈ Z e − i ω n τ G ~ p q R ( i ω n ) , ω n = ( 2 n + 1 ) π β (95) G^{\mathrm{M}}_{pq}(\tau) = \frac{1}{\beta}\sum_{n\in\mathbb{Z}}e^{-i\omega_n\tau}\widetilde{G}^{\mathrm{R}}_{pq}(i\omega_n), \qquad\omega_n = \frac{(2n+1)\pi}{\beta} \tag*{(95)} G pq M ( τ ) = β 1 n ∈ Z ∑ e − i ω n τ G pq R ( i ω n ) , ω n = β ( 2 n + 1 ) π (95) using the standard Matsubara summation; see Refs. [111 Appendix C] and [112 Appendix A] .
*Continuum Dyson reconstruction.—*To remove the pseudomode-fitting error from the impurity response, we extract the C3 ^3 3 Δ R ( z ) \Delta^{R}(z) Δ R ( z ) (112) :
Σ ~ R ( z ) = z − h − Δ f i t R ( z ) − [ G ~ R ( z ) ] − 1 , G ~ c o n t R ( z ) = [ z − h − Δ R ( z ) − Σ ~ R ( z ) ] − 1 , (96) \widetilde{\Sigma}^{R}(z) = z - h - \Delta_{\mathrm{fit}}^{R}(z) - [\widetilde{G}^{R}(z)]^{-1}, \qquad\widetilde{G}_{\mathrm{cont}}^{R}(z) = [z - h - \Delta^{R}(z) - \widetilde{\Sigma}^{R}(z)]^{-1}, \tag*{(96)} Σ R ( z ) = z − h − Δ fit R ( z ) − [ G R ( z ) ] − 1 , G cont R ( z ) = [ z − h − Δ R ( z ) − Σ R ( z ) ] − 1 , (96) where h h h
*Benchmark on the Kanamori model.—*For the interaction in Eq. (2) , we set U = 2 U = 2 U = 2 J H = 0.2 J^{\mathrm{H}} = 0.2 J H = 0.2 μ = ( 3 U − 5 J H ) / 2 = 2.5 \mu= (3U - 5J^{\mathrm{H}})/2 = 2.5 μ = ( 3 U − 5 J H ) /2 = 2.5
J i j ( ω ) = J i j 0 J 0 ( ω ) and J 0 ( ω ) = 1 π 1 − ω 2 W 2 1 ( − W , W ) ( ω ) , (97) J_{ij}(\omega) = J^{0}_{ij}J_0(\omega) \qquad\text{and} \qquad J_0(\omega) = \frac{1}{\pi}\sqrt{1-\frac{\omega^2}{W^2}}\,\mathbf{1}_{(-W,W)}(\omega), \tag*{(97)} J ij ( ω ) = J ij 0 J 0 ( ω ) and J 0 ( ω ) = π 1 1 − W 2 ω 2 1 ( − W , W ) ( ω ) , (97) where the half-bandwidth W = 2 W = 2 W = 2 J i j 0 = δ i j + r ( 1 − δ i j ) J^{0}_{ij} = \delta_{ij} + r(1-\delta_{ij}) J ij 0 = δ ij + r ( 1 − δ ij ) r = 1 r = 1 r = 1 N p m = 5 N_{\mathrm{pm}} = 5 N pm = 5 β = 8 \beta= 8 β = 8 β = 64 \beta= 64 β = 64 r = 1 r = 1 r = 1 J 0 = ( 1 1 1 1 ) J^0 = \begin{pmatrix}1 & 1 \\ 1 & 1\end{pmatrix} J 0 = ( 1 1 1 1 ) 1 1 1 d = 1 d = 1 d = 1 2 n = 2 ( n I + 2 d N p m ) = 24 2n = 2(n_I + 2dN_{\mathrm{pm}}) = 24 2 n = 2 ( n I + 2 d N pm ) = 24 3 ^3 3 (96) . For the spectra, we use the broadening ω ↦ ω + i η \omega\mapsto\omega+ i\eta ω ↦ ω + i η η = 0.05 \eta= 0.05 η = 0.05 m K m_K m K S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 β = 8 \beta= 8 β = 8 β = 64 \beta= 64 β = 64
Figure S4 and Figure S4 show additional results away from half filling for d = 4 , 6 d = 4, 6 d = 4 , 6 μ = 1 \mu= 1 μ = 1 N p m = 5 N_{\mathrm{pm}} = 5 N pm = 5 η = 0.05 \eta= 0.05 η = 0.05 m K ( S 1 , 3 ) = 1024 m_K(S_{1,3}) = 1024 m K ( S 1 , 3 ) = 1024 m K ( S 1 , 3 , 5 ) = 1536 m_K(S_{1,3,5}) = 1536 m K ( S 1 , 3 , 5 ) = 1536 β = 8 \beta= 8 β = 8 m K ( S 1 , 3 ) = 1024 m_K(S_{1,3}) = 1024 m K ( S 1 , 3 ) = 1024 m K ( S 1 , 3 , 5 ) = 2048 m_K(S_{1,3,5}) = 2048 m K ( S 1 , 3 , 5 ) = 2048 β = 64 \beta= 64 β = 64 β = 8 \beta= 8 β = 8 G 00 M G^{\mathrm{M}}_{00} G 00 M G 01 M G^{\mathrm{M}}_{01} G 01 M [117 ] ; the fifth-order strong coupling expansion (SCE) results of Ref. [118 ] are taken from the numerical data accompanying Ref. [120 ] . The inchworm results at both temperatures and the CT-HYB results at β = 64 \beta= 64 β = 64 [116 ] and are for G 00 M G^{\mathrm{M}}_{00} G 00 M
Figure S4. Additional Kanamori results at μ = 1 \mu= 1 μ = 1 β = 8 \beta= 8 β = 8 β = 64 \beta= 64 β = 64 η = 0.05 \eta= 0.05 η = 0.05 [116 ] .
*Benchmark against real-time inchworm QMC.—*For the general interacting Hamiltonian in Eq. (3) [121 ] , we use ϵ 0 = V = V ′ ′ = ν ′ = − 2 \epsilon_0 = V = V^{\prime\prime} = \nu^{\prime} = -2 ϵ 0 = V = V ′′ = ν ′ = − 2 U = ν = 2 U = \nu= 2 U = ν = 2 β = 0.1 \beta= 0.1 β = 0.1
J i j ( ω ) = J i j 0 J 0 ( ω ) and J 0 ( ω ) = 1 π ( 1 + e υ ( ω − ω c ) ) ( 1 + e − υ ( ω + ω c ) ) , (98) J_{ij}(\omega) = J^0_{ij}J_0(\omega) \quad\text{and} \quad J_0(\omega) = \frac{1}{\pi(1+e^{\upsilon(\omega-\omega_c)})(1+e^{-\upsilon(\omega+\omega_c)})}, \tag*{(98)} J ij ( ω ) = J ij 0 J 0 ( ω ) and J 0 ( ω ) = π ( 1 + e υ ( ω − ω c ) ) ( 1 + e − υ ( ω + ω c ) ) 1 , (98) where υ = 10 \upsilon= 10 υ = 10 ω c = 25 \omega_c = 25 ω c = 25 J 11 0 = J 22 0 = 3 / 4 J^0_{11} = J^0_{22} = 3/4 J 11 0 = J 22 0 = 3/4 J 12 0 = J 21 0 = 1 / 4 J^0_{12} = J^0_{21} = 1/4 J 12 0 = J 21 0 = 1/4 N p m = 5 N_{\mathrm{pm}} = 5 N pm = 5 J 0 = ( 3 / 4 1 / 4 1 / 4 3 / 4 ) J^0 = \begin{pmatrix}3/4 & 1/4 \\ 1/4 & 3/4\end{pmatrix} J 0 = ( 3/4 1/4 1/4 3/4 ) d = 2 d = 2 d = 2 2 n = 2 ( n I + 2 d N p m ) = 44 2n = 2(n_I + 2dN_{\mathrm{pm}}) = 44 2 n = 2 ( n I + 2 d N pm ) = 44 A ( ω ) A(\omega) A ( ω ) η = 0 \eta= 0 η = 0 m K m_{\mathrm{K}} m K S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
Exchange-coupled molecular junction. —For the Hamiltonian in Eq. (4) and the transport in Figure 4 , we use four orbitals ( h 2 , h 1 , l 1 , l 2 ) (\mathrm{h2},\mathrm{h1},\mathrm{l1},\mathrm{l2}) ( h2 , h1 , l1 , l2 ) a ^ i σ \hat{a}_{i\sigma} a ^ iσ β = 10 \beta= 10 β = 10 n ^ i σ = a ^ i σ † a ^ i σ \hat{n}_{i\sigma} = \hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\sigma} n ^ iσ = a ^ iσ † a ^ iσ n ^ i = ∑ σ n ^ i σ \hat{n}_i = \sum_{\sigma}\hat{n}_{i\sigma} n ^ i = ∑ σ n ^ iσ S ^ i = ( S ^ i x , S ^ i y , S ^ i z ) \hat{\mathbf{S}}_i = (\hat{S}^x_i,\hat{S}^y_i,\hat{S}^z_i) S ^ i = ( S ^ i x , S ^ i y , S ^ i z ) S ^ i ν = 1 2 ∑ σ σ ′ a ^ i σ † σ σ σ ′ ν a ^ i σ ′ \hat{S}^{\nu}_i = \frac{1}{2}\sum_{\sigma\sigma'}\hat{a}^{\dagger}_{i\sigma}\sigma^{\nu}_{\sigma\sigma'}\hat{a}_{i\sigma'} S ^ i ν = 2 1 ∑ σ σ ′ a ^ iσ † σ σ σ ′ ν a ^ i σ ′ σ ν \sigma^\nu σ ν ( ν = x , y , z ) (\nu= x,y,z) ( ν = x , y , z )
h = ( − 2.95 0.05 0 0 0.05 − 2.55 0 0 0 0 − 0.71 0.04 0 0 0.04 − 0.48 ) , U = 0.8 , V = 0.3 , J E = − J P = ( 0.04 0.18 0.25 0.04 ) . (99) h = \begin{pmatrix}
-2.95 & 0.05 & 0 & 0 \\
0.05 & -2.55 & 0 & 0 \\
0 & 0 & -0.71 & 0.04 \\
0 & 0 & 0.04 & -0.48
\end{pmatrix}, \quad U = 0.8,\quad V = 0.3,\quad J^{\mathrm{E}} = -J^{\mathrm{P}} = \begin{pmatrix}0.04 & 0.18 \\ 0.25 & 0.04\end{pmatrix}. \tag*{(99)} h = − 2.95 0.05 0 0 0.05 − 2.55 0 0 0 0 − 0.71 0.04 0 0 0.04 − 0.48 , U = 0.8 , V = 0.3 , J E = − J P = ( 0.04 0.25 0.18 0.04 ) . (99) The rows of J E , P J^{\mathrm{E,P}} J E , P ( h 2 , h 1 ) (\mathrm{h2},\mathrm{h1}) ( h2 , h1 ) ( l 1 , l 2 ) (\mathrm{l1},\mathrm{l2}) ( l1 , l2 )
J r ( ω ) = J r 0 J 0 ( ω − μ r ) u r u r † ( r ∈ { s , t } ) , J 0 ( ω ) = 1 π 1 − ω 2 W 2 1 ( − W , W ) ( ω ) , W = 4 , J s 0 = 0.10 , J t 0 = 0.01 , (100) J_r(\omega) = J^0_r J_0(\omega-\mu_r)\mathbf{u}_r\mathbf{u}^{\dagger}_r \quad(r \in\{\mathrm{s},\mathrm{t}\}), \quad J_0(\omega) = \frac{1}{\pi}\sqrt{1-\frac{\omega^2}{W^2}}\,\mathbf{1}_{(-W,W)}(\omega), \quad W=4,\quad J^0_{\mathrm{s}}=0.10,\quad J^0_{\mathrm{t}}=0.01, \tag*{(100)} J r ( ω ) = J r 0 J 0 ( ω − μ r ) u r u r † ( r ∈ { s , t }) , J 0 ( ω ) = π 1 1 − W 2 ω 2 1 ( − W , W ) ( ω ) , W = 4 , J s 0 = 0.10 , J t 0 = 0.01 , (100) with u s = ( 1 , 1 , 1 , 1 ) T \mathbf{u}_{\mathrm{s}} = (1,1,1,1)^{\mathsf{T}} u s = ( 1 , 1 , 1 , 1 ) T u t = ( 0.4 , 1 , 0.12 , 0.08 ) T ∣ ( 0.4 , 1 , 0.12 , 0.08 ) T ∣ \mathbf{u}_{\mathrm{t}} = \frac{(0.4,1,0.12,0.08)^{\mathsf{T}}}{\lvert(0.4,1,0.12,0.08)^{\mathsf{T}}\rvert} u t = ∣( 0.4 , 1 , 0.12 , 0.08 ) T ∣ ( 0.4 , 1 , 0.12 , 0.08 ) T μ s = 0 \mu_{\mathrm{s}} = 0 μ s = 0 μ t = Φ \mu_{\mathrm{t}} = \Phi μ t = Φ Δ ≷ = Δ s ≷ + Δ t ≷ \Delta^{\gtrless} = \Delta^{\gtrless}_{\mathrm{s}} + \Delta^{\gtrless}_{\mathrm{t}} Δ ≷ = Δ s ≷ + Δ t ≷ Δ r ≷ ( t ) = ∫ J r ( ω ) f r ≷ ( ω ) e − i ω t d ω \Delta^{\gtrless}_r(t) = \int J_r(\omega)f^{\gtrless}_r(\omega)e^{-i\omega t}\,d\omega Δ r ≷ ( t ) = ∫ J r ( ω ) f r ≷ ( ω ) e − iω t d ω f r < = f F D β , μ r f^{<}_r = f^{\beta,\mu_r}_{\mathrm{FD}} f r < = f FD β , μ r f r > = 1 − f F D β , μ r f^{>}_r = 1-f^{\beta,\mu_r}_{\mathrm{FD}} f r > = 1 − f FD β , μ r Φ < 0 \Phi< 0 Φ < 0 N p m , r = N p m , r > = 4 N_{\mathrm{pm},r} = N^{>}_{\mathrm{pm},r} = 4 N pm , r = N pm , r > = 4 2 n = 40 2n = 40 2 n = 40 H ^ r t ↦ H ^ r t + Φ I \hat{H}^{\mathrm{t}}_r \mapsto\hat{H}^{\mathrm{t}}_r + \Phi I H ^ r t ↦ H ^ r t + Φ I g ^ t ≷ \hat{g}^{\gtrless}_{\mathrm{t}} g ^ t ≷ Γ ^ t ≷ \hat{\Gamma}^{\gtrless}_{\mathrm{t}} Γ ^ t ≷ Φ = − 0.2 , − 0.4 , … , − 1.6 \Phi= -0.2,-0.4,\ldots,-1.6 Φ = − 0.2 , − 0.4 , … , − 1.6 η = 0 \eta= 0 η = 0 m K m_{\mathrm{K}} m K S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
Molecular-orbital basis and equilibrium comparison. —The fixed molecular operators f ^ m σ = ∑ i ψ i m ∗ a ^ i σ \hat{f}_{m\sigma} = \sum_i \psi^*_{im}\hat{a}_{i\sigma} f ^ mσ = ∑ i ψ im ∗ a ^ iσ m m m − - − + + +
H ^ I = ∑ m n σ h ~ m n f ^ m σ † f ^ n σ + ∑ m n r s , σ σ ′ V ~ m n r s f ^ m σ † f ^ n σ ′ † f ^ r σ ′ f ^ s σ . (101) \hat{H}_{\mathrm{I}} = \sum_{mn\sigma}\widetilde{h}_{mn}\hat{f}^{\dagger}_{m\sigma}\hat{f}_{n\sigma} + \sum_{mnrs,\sigma\sigma'}\widetilde{V}_{mnrs}\hat{f}^{\dagger}_{m\sigma}\hat{f}^{\dagger}_{n\sigma'}\hat{f}_{r\sigma'}\hat{f}_{s\sigma}. \tag*{(101)} H ^ I = mnσ ∑ h mn f ^ mσ † f ^ nσ + mn r s , σ σ ′ ∑ V mn r s f ^ mσ † f ^ n σ ′ † f ^ r σ ′ f ^ s σ . (101) where h ~ = Ψ † h Ψ \widetilde{h}=\Psi^\dagger h\Psi h = Ψ † h Ψ V ~ m n r s = ∑ i j k l Ψ i m ∗ Ψ j n ∗ V i j k l Ψ k r Ψ l s \widetilde{V}_{mnrs}=\sum_{ijkl}\Psi^*_{im}\Psi^*_{jn}V_{ijkl}\Psi_{kr}\Psi_{ls} V mn r s = ∑ ij k l Ψ im ∗ Ψ j n ∗ V ij k l Ψ k r Ψ l s V i j k l V_{ijkl} V ij k l (1) . The molecular-basis Green’s functions are G M , M O = Ψ † G M Ψ G^{\mathrm{M,MO}}=\Psi^\dagger G^{\mathrm{M}}\Psi G M , MO = Ψ † G M Ψ G L U M O M = ( G M , M O ) L U M O , L U M O G^{\mathrm{M}}_{\mathrm{LUMO}}=(G^{\mathrm{M,MO}})_{\mathrm{LUMO,LUMO}} G LUMO M = ( G M , MO ) LUMO , LUMO f ^ L U M O , σ \hat{f}_{\mathrm{LUMO},\sigma} f ^ LUMO , σ (95) , with the leading 1 / ( i ω n ) 1/(i\omega_n) 1/ ( i ω n ) [111 Appendix C] and [112 Appendix A] ).
We diagonalize H ^ I \hat{H}_I H ^ I { ψ i } \{\psi_i\} { ψ i } b ^ t σ = ∑ i ( u t ) i a ^ i σ \hat{b}_{t\sigma}=\sum_i(u_t)_i\hat{a}_{i\sigma} b ^ t σ = ∑ i ( u t ) i a ^ iσ M \mathcal{M} M
A i s o , M < ( ω ) = 1 ∣ M ∣ ∑ k ∈ M , f , σ ∣ ⟨ ψ f ∣ b ^ t σ ∣ ψ k ⟩ ∣ 2 δ ( ω − ( E k − E f ) ) . (102) A^{<}_{\mathrm{iso},\mathcal{M}}(\omega)=\frac{1}{|\mathcal{M}|}\sum_{k\in\mathcal{M},f,\sigma}\left|\langle\psi_f|\hat{b}_{t\sigma}|\psi_k\rangle\right|^2\delta\left(\omega-(E_k-E_f)\right). \tag*{(102)} A iso , M < ( ω ) = ∣ M ∣ 1 k ∈ M , f , σ ∑ ⟨ ψ f ∣ b ^ t σ ∣ ψ k ⟩ 2 δ ( ω − ( E k − E f ) ) . (102) Here M \mathcal{M} M f f f Figure 4 d to Figure 4 f have heights given by these weights ∣ ⟨ ψ f ∣ b ^ t σ ∣ ψ k ⟩ ∣ 2 \left|\langle\psi_f|\hat{b}_{t\sigma}|\psi_k\rangle\right|^2 ⟨ ψ f ∣ b ^ t σ ∣ ψ k ⟩ 2
Orbital-transfer strengths. —In the original basis, define a ^ i = ( a ^ i ↑ , a ^ i ↓ ) T \hat{\boldsymbol{a}}_i=(\hat{a}_{i\uparrow},\hat{a}_{i\downarrow})^\mathsf{T} a ^ i = ( a ^ i ↑ , a ^ i ↓ ) T l 1 → h 1 \mathrm{l}_1\to\mathrm{h}_1 l 1 → h 1 [122 , 123 ] ,
B ^ 0 = 1 2 a ^ h 1 † a ^ l 1 , B ^ 1 ν = 1 2 a ^ h 1 † σ ν a ^ l 1 ( ν = x , y , z ) . (103) \hat{B}_0=\frac{1}{\sqrt{2}}\hat{a}^\dagger_{\mathrm{h}_1}\hat{a}_{\mathrm{l}_1},\qquad\hat{B}_1^\nu=\frac{1}{\sqrt{2}}\hat{a}^\dagger_{\mathrm{h}_1}\sigma^\nu\hat{a}_{\mathrm{l}_1}\qquad(\nu=x,y,z). \tag*{(103)} B ^ 0 = 2 1 a ^ h 1 † a ^ l 1 , B ^ 1 ν = 2 1 a ^ h 1 † σ ν a ^ l 1 ( ν = x , y , z ) . (103) Their integrated steady-state strengths are
W 0 = ⟨ B ^ 0 † B ^ 0 ⟩ = 1 2 ⟨ n ^ l 1 ⟩ − 1 4 ⟨ n ^ h 1 n ^ l 1 ⟩ − ⟨ S ^ h 1 ⋅ S ^ l 1 ⟩ , W 1 = ∑ ν = x , y , z ⟨ ( B ^ 1 ν ) † B ^ 1 ν ⟩ = 3 2 ⟨ n ^ l 1 ⟩ − 3 4 ⟨ n ^ h 1 n ^ l 1 ⟩ + ⟨ S ^ h 1 ⋅ S ^ l 1 ⟩ . (104) W_0=\langle\hat{B}_0^\dagger\hat{B}_0\rangle=\frac{1}{2}\langle\hat{n}_{\mathrm{l}_1}\rangle-\frac{1}{4}\langle\hat{n}_{\mathrm{h}_1}\hat{n}_{\mathrm{l}_1}\rangle-\langle\hat{\boldsymbol{S}}_{\mathrm{h}_1}\cdot\hat{\boldsymbol{S}}_{\mathrm{l}_1}\rangle,\qquad W_1=\sum_{\nu=x,y,z}\langle(\hat{B}_1^\nu)^\dagger\hat{B}_1^\nu\rangle=\frac{3}{2}\langle\hat{n}_{\mathrm{l}_1}\rangle-\frac{3}{4}\langle\hat{n}_{\mathrm{h}_1}\hat{n}_{\mathrm{l}_1}\rangle+\langle\hat{\boldsymbol{S}}_{\mathrm{h}_1}\cdot\hat{\boldsymbol{S}}_{\mathrm{l}_1}\rangle. \tag*{(104)} W 0 = ⟨ B ^ 0 † B ^ 0 ⟩ = 2 1 ⟨ n ^ l 1 ⟩ − 4 1 ⟨ n ^ h 1 n ^ l 1 ⟩ − ⟨ S ^ h 1 ⋅ S ^ l 1 ⟩ , W 1 = ν = x , y , z ∑ ⟨( B ^ 1 ν ) † B ^ 1 ν ⟩ = 2 3 ⟨ n ^ l 1 ⟩ − 4 3 ⟨ n ^ h 1 n ^ l 1 ⟩ + ⟨ S ^ h 1 ⋅ S ^ l 1 ⟩ . (104) The ratio is defined as R = W 1 / ( 3 W 0 ) R=W_1/(3W_0) R = W 1 / ( 3 W 0 ) R R R G W GW G W l 1 → h 1 \mathrm{l}_1\to\mathrm{h}_1 l 1 → h 1
For an observable O ^ \hat{O} O ^ O O O O ^ \hat{O} O ^ O S = ⟨ ⟨ I ∣ O ∣ Φ S ⟩ ⟩ O_S=\langle\!\langle I|O|\Phi_S\rangle\!\rangle O S = ⟨ ⟨ I ∣ O ∣ Φ S ⟩ ⟩ ∣ Φ 0 ⟩ ⟩ = ∣ ρ 0 ⟩ ⟩ |\Phi_0\rangle\!\rangle=|\rho_0\rangle\!\rangle ∣ Φ 0 ⟩ ⟩ = ∣ ρ 0 ⟩ ⟩
⟨ O ^ ⟩ C C S D = ⟨ ⟨ I ∣ O e T ∣ ρ 0 ⟩ ⟩ = O 0 + ∑ ∣ S ∣ = 2 O S t S ( 2 ) + ∑ ∣ S ∣ = 4 O S ( t S ( 4 ) + 1 2 t S ( 2 ) ∧ t S ( 2 ) ) . (105) \langle\hat{O}\rangle_{\mathrm{CCSD}}=\langle\!\langle I|Oe^T|\rho_0\rangle\!\rangle=O_0+\sum_{|S|=2}O_S t_S^{(2)}+\sum_{|S|=4}O_S\left(t_S^{(4)}+\frac{1}{2}t_S^{(2)}\wedge t_S^{(2)}\right). \tag*{(105)} ⟨ O ^ ⟩ CCSD = ⟨ ⟨ I ∣ O e T ∣ ρ 0 ⟩ ⟩ = O 0 + ∣ S ∣ = 2 ∑ O S t S ( 2 ) + ∣ S ∣ = 4 ∑ O S ( t S ( 4 ) + 2 1 t S ( 2 ) ∧ t S ( 2 ) ) . (105) Energy-resolved contact currents. —All matrices below use the original basis. The per-spin removal and addition spectra are A < ( ω ) = − i G ~ e < ( ω ) / ( 2 π ) A^{<}(\omega)=-i\widetilde{G}^{<}_{e}(\omega)/(2\pi) A < ( ω ) = − i G e < ( ω ) / ( 2 π ) A > ( ω ) = i G ~ e > ( ω ) / ( 2 π ) A^{>}(\omega)=i\widetilde{G}^{>}_{e}(\omega)/(2\pi) A > ( ω ) = i G e > ( ω ) / ( 2 π ) Ξ r < \Xi_r^{<} Ξ r < Ξ r > \Xi_r^{>} Ξ r > [124 , 125 ] :
j r ( ω ) = 2 Re Tr [ Ξ r < ( ω ) A > ( ω ) − Ξ r > ( ω ) A < ( ω ) ] , I r ( ω ) = − j r ( ω ) , I r = ∫ R I r ( ω ) d ω . (106) j_r(\omega)=2\operatorname{Re}\operatorname{Tr}\left[\Xi_r^{<}(\omega)A^{>}(\omega)-\Xi_r^{>}(\omega)A^{<}(\omega)\right],\qquad\mathcal{I}_r(\omega)=-j_r(\omega),\qquad I_r=\int_{\mathbb{R}}\mathcal{I}_r(\omega)\,d\omega. \tag*{(106)} j r ( ω ) = 2 Re Tr [ Ξ r < ( ω ) A > ( ω ) − Ξ r > ( ω ) A < ( ω ) ] , I r ( ω ) = − j r ( ω ) , I r = ∫ R I r ( ω ) d ω . (106) Here I r \mathcal{I}_r I r I r I_r I r Ξ r ≷ = 2 π f r ≷ J r \Xi_r^{\gtrless}=2\pi f_r^{\gtrless}J_r Ξ r ≷ = 2 π f r ≷ J r H r ≷ H_r^{\gtrless} H r ≷ Γ r ≷ \Gamma_r^{\gtrless} Γ r ≷ g r ≷ g_r^{\gtrless} g r ≷ (6) and (7) :
g r R , ≷ ( ω ) = [ ω I − H r ≷ + i Γ r ≷ ] − 1 , g r A , ≷ = ( g r R , ≷ ) † , Ξ r ≷ = 2 ( g r ≷ ) † g r R , ≷ Ξ r ≷ g r A , ≷ . (107) g_r^{R,\gtrless}(\omega)=\left[\omega I-H_r^{\gtrless}+i\Gamma_r^{\gtrless}\right]^{-1},\qquad g_r^{A,\gtrless}=(g_r^{R,\gtrless})^\dagger,\qquad\Xi_r^{\gtrless}=2(g_r^{\gtrless})^\dagger g_r^{R,\gtrless}\Xi_r^{\gtrless}g_r^{A,\gtrless}. \tag*{(107)} g r R , ≷ ( ω ) = [ ω I − H r ≷ + i Γ r ≷ ] − 1 , g r A , ≷ = ( g r R , ≷ ) † , Ξ r ≷ = 2 ( g r ≷ ) † g r R , ≷ Ξ r ≷ g r A , ≷ . (107) The static extraction current is evaluated independently from the full spin-orbital 1-RDM γ p q = ⟨ c ^ q † c ^ p ⟩ \gamma_{pq}=\langle\hat{c}_q^\dagger\hat{c}_p\rangle γ pq = ⟨ c ^ q † c ^ p ⟩ I I I I r I_r I r r r r
I r = − 2 Im ∑ p ∈ I , q ∈ I r H p q γ q p . (108) I_r=-2\operatorname{Im}\sum_{p\in I,q\in I_r}H_{pq}\gamma_{qp}. \tag*{(108)} I r = − 2 Im p ∈ I , q ∈ I r ∑ H pq γ q p . (108) Here H H H H ^ C L \hat{H}_{\mathrm{CL}} H ^ CL
The combined runtime for computing the CCSD steady states, the imaginary-time Green’s functions and real-frequency spectra on an NVIDIA H200 GPU is 122 s for S 1 , 3 S_{1,3} S 1 , 3 S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5
Strongly correlated two-impurity Anderson model with exchange. —For the benchmark in Figure 5 , we consider J E = 2 J^{E}=2 J E = 2 U = 4 U=4 U = 4 β = 100 \beta=100 β = 100 W = 10 W=10 W = 10
J i j ( ω ) = δ i j π [ 1 + e ( ω − W ) / 0.1 ] [ 1 + e ( − ω − W ) / 0.1 ] ≈ δ i j π Θ ( W − ∣ ω ∣ ) . (109) J_{ij}(\omega)=\frac{\delta_{ij}}{\pi[1+e^{(\omega-W)/0.1}][1+e^{(-\omega-W)/0.1}]}\approx\frac{\delta_{ij}}{\pi}\Theta(W-\lvert\omega\rvert). \tag*{(109)} J ij ( ω ) = π [ 1 + e ( ω − W ) /0.1 ] [ 1 + e ( − ω − W ) /0.1 ] δ ij ≈ π δ ij Θ ( W − ∣ ω ∣) . (109) For each impurity-spin channel, we use N p m = 6 N_{\mathrm{pm}}=6 N pm = 6 t ∈ [ 0 , 100 ] t\in[0,100] t ∈ [ 0 , 100 ] [126 , 127 ] , with max nfev=3000 in the final stage. For one channel, we write J = J i i J=J_{ii} J = J ii
Δ f i t ≶ ( z ) = ( g ≶ ) † ( z I − H ≶ + i Γ ≶ ) − 1 g ≶ , Δ f i t R = Δ f i t < + Δ f i t > , J f i t o c c ( ω ) = − 1 π Im Δ f i t < ( ω ) , J f i t e m p t y ( ω ) = − 1 π Im Δ f i t > ( ω ) . (110) \Delta_{\mathrm{fit}}^{\lessgtr}(z)=(g^{\lessgtr})^{\dagger}(zI-H^{\lessgtr}+i\Gamma^{\lessgtr})^{-1}g^{\lessgtr},\qquad\Delta_{\mathrm{fit}}^{R}=\Delta_{\mathrm{fit}}^{<}+\Delta_{\mathrm{fit}}^{>},\qquad J_{\mathrm{fit}}^{\mathrm{occ}}(\omega)=-\frac{1}{\pi}\operatorname{Im}\Delta_{\mathrm{fit}}^{<}(\omega),\qquad J_{\mathrm{fit}}^{\mathrm{empty}}(\omega)=-\frac{1}{\pi}\operatorname{Im}\Delta_{\mathrm{fit}}^{>}(\omega). \tag*{(110)} Δ fit ≶ ( z ) = ( g ≶ ) † ( z I − H ≶ + i Γ ≶ ) − 1 g ≶ , Δ fit R = Δ fit < + Δ fit > , J fit occ ( ω ) = − π 1 Im Δ fit < ( ω ) , J fit empty ( ω ) = − π 1 Im Δ fit > ( ω ) . (110) The corresponding time-domain hybridization is Δ f i t < ( t ) = ( g < ) † e ( − i H < − Γ < ) t g < \Delta_{\mathrm{fit}}^{<}(t)=(g^{<})^{\dagger}e^{(-iH^{<}-\Gamma^{<})t}g^{<} Δ fit < ( t ) = ( g < ) † e ( − i H < − Γ < ) t g <
min θ 1 ∣ Ω ∣ ∑ ω ∈ Ω w ω [ ∣ Δ f i t R ( ω ) − Δ R ( ω ) ∣ 2 + π 2 ∣ J f i t o c c ( ω ) − J ( ω ) f F D 100 , 0 ( ω ) ∣ 2 ] + χ t ∑ t ∈ T w t ∣ Δ f i t < ( t ) − Δ < ( t ) ∣ 2 ∑ t ∈ T w t ∣ Δ < ( t ) ∣ 2 + χ r ∥ θ − θ 0 ∥ 2 2 . (111) \min_{\theta}\frac{1}{\lvert\Omega\rvert}\sum_{\omega\in\Omega}w_{\omega}\left[\left\lvert\Delta_{\mathrm{fit}}^{R}(\omega)-\Delta^{R}(\omega)\right\rvert^{2}+\pi^{2}\left\lvert J_{\mathrm{fit}}^{\mathrm{occ}}(\omega)-J(\omega)f_{\mathrm{FD}}^{100,0}(\omega)\right\rvert^{2}\right]+\chi_{t}\frac{\sum_{t\in\mathcal{T}}w_{t}\left\lvert\Delta_{\mathrm{fit}}^{<}(t)-\Delta^{<}(t)\right\rvert^{2}}{\sum_{t\in\mathcal{T}}w_{t}\left\lvert\Delta^{<}(t)\right\rvert^{2}}+\chi_{r}\lVert\theta-\theta_{0}\rVert_{2}^{2}. \tag*{(111)} θ min ∣ Ω ∣ 1 ω ∈ Ω ∑ w ω [ Δ fit R ( ω ) − Δ R ( ω ) 2 + π 2 J fit occ ( ω ) − J ( ω ) f FD 100 , 0 ( ω ) 2 ] + χ t ∑ t ∈ T w t ∣ Δ < ( t ) ∣ 2 ∑ t ∈ T w t ∣ Δ fit < ( t ) − Δ < ( t ) ∣ 2 + χ r ∥ θ − θ 0 ∥ 2 2 . (111) Here, the exact retarded hybridization is the Cauchy transform of the bath spectral density. The inverse relation follows from the Plemelj formula:
Δ R ( z ) = ∫ − ∞ ∞ J ( ϵ ) z − ϵ d ϵ ( Im z > 0 ) , J ( ω ) = i 2 π lim η ↓ 0 [ Δ R ( ω + i η ) − Δ R ( ω + i η ) † ] . (112) \Delta^{R}(z)=\int_{-\infty}^{\infty}\frac{J(\epsilon)}{z-\epsilon}\,d\epsilon\quad(\operatorname{Im}z>0),\qquad J(\omega)=\frac{i}{2\pi}\lim_{\eta\downarrow0}\left[\Delta^{R}(\omega+i\eta)-\Delta^{R}(\omega+i\eta)^{\dagger}\right]. \tag*{(112)} Δ R ( z ) = ∫ − ∞ ∞ z − ϵ J ( ϵ ) d ϵ ( Im z > 0 ) , J ( ω ) = 2 π i η ↓ 0 lim [ Δ R ( ω + i η ) − Δ R ( ω + i η ) † ] . (112) We set χ t = 0.04 \chi_{t}=0.04 χ t = 0.04 χ r = 10 − 12 \chi_{r}=10^{-12} χ r = 1 0 − 12 ℓ 2 \ell^{2} ℓ 2 w ω = 1 + 15 e − ( ω / 0.3 ) 2 w_{\omega}=1+15e^{-(\omega/0.3)^{2}} w ω = 1 + 15 e − ( ω /0.3 ) 2 w t = 1 + t / 2 w_{t}=1+t/2 w t = 1 + t /2 θ \theta θ H < H^{<} H < Υ \Upsilon Υ g < g^{<} g < Γ < = Υ Υ † ⪰ 0 \Gamma^{<}=\Upsilon\Upsilon^{\dagger}\succeq0 Γ < = Υ Υ † ⪰ 0 θ 0 \theta_{0} θ 0 H > = − ( H < ) ∗ H^{>} = -(H^{<})^{*} H > = − ( H < ) ∗ Γ > = ( Γ < ) ∗ \Gamma^{>} = (\Gamma^{<})^{*} Γ > = ( Γ < ) ∗ g > = ( g < ) ∗ g^{>}=(g^{<})^{*} g > = ( g < ) ∗ Ω \Omega Ω [ − 12 , 12 ] [-12,12] [ − 12 , 12 ] [ − 2 , 2 ] [-2,2] [ − 2 , 2 ] [ − 0.2 , 0.2 ] [-0.2,0.2] [ − 0.2 , 0.2 ] T = { 0 , 0.025 , … , 2 } ∪ { 0 , 0.2 , … , 200 } \mathcal{T}=\{0,0.025,\ldots,2\}\cup\{0,0.2,\ldots,200\} T = { 0 , 0.025 , … , 2 } ∪ { 0 , 0.2 , … , 200 } S 1 , 3 , 5 S_{1,3,5} S 1 , 3 , 5 η = 0.01 \eta=0.01 η = 0.01 m K = 768 m_{K}=768 m K = 768 ∣ ω ∣ ≤ 0.5 \lvert\omega\rvert\leq0.5 ∣ ω ∣ ≤ 0.5 m K = 256 m_{K}=256 m K = 256 0.1 < ∣ ω ∣ ≤ 0.5 0.1<\lvert\omega\rvert\leq0.5 0.1 < ∣ ω ∣ ≤ 0.5 m K = 768 m_{K}=768 m K = 768 ∣ ω ∣ ≤ 0.1 \lvert\omega\rvert\leq0.1 ∣ ω ∣ ≤ 0.1 Figure 5 compares the total bath spectral density and Fermi-edge resolution for different numbers of pseudomodes.
A simple mechanism for MBPT descriptions of exchange-induced multiplet Green’s functions
We present an isolated two-orbital model with closed-form exact, self-consistent 2ndB and GW Green’s functions. We will observe that the missing physical interaction in those MBPT approximations can lead to qualitatively incorrect descriptions of the exchange-induced multiplet structure. The Hamiltonian and reference state are
H ^ = ϵ n ^ 1 + U ( n ^ 1 ↑ n ^ 1 ↓ + n ^ 2 ↑ n ^ 2 ↓ ) − J E ( 2 S ^ 1 ⋅ S ^ 2 + 1 2 n ^ 1 n ^ 2 ) , ∣ ψ ⟩ = ∣ 2 ↑ ⟩ , ϵ , U , J E > 0. (113) \hat{H}=\epsilon\hat{n}_{1}+U(\hat{n}_{1\uparrow}\hat{n}_{1\downarrow}+\hat{n}_{2\uparrow}\hat{n}_{2\downarrow})-J^{E}\left(2\hat{\mathbf{S}}_{1}\cdot\hat{\mathbf{S}}_{2}+\frac{1}{2}\hat{n}_{1}\hat{n}_{2}\right),\qquad\lvert\psi\rangle=\lvert2\uparrow\rangle,\qquad\epsilon,U,J^{E}>0. \tag*{(113)} H ^ = ϵ n ^ 1 + U ( n ^ 1 ↑ n ^ 1 ↓ + n ^ 2 ↑ n ^ 2 ↓ ) − J E ( 2 S ^ 1 ⋅ S ^ 2 + 2 1 n ^ 1 n ^ 2 ) , ∣ ψ ⟩ = ∣ 2 ↑ ⟩ , ϵ , U , J E > 0. (113) The one-particle reference ∣ ψ ⟩ \lvert\psi\rangle ∣ ψ ⟩ n ^ 1 , n ^ 2 \hat{n}_{1},\hat{n}_{2} n ^ 1 , n ^ 2 S ^ z \hat{S}^{z} S ^ z G i σ , j σ ′ ≶ ( t ) = δ i j δ σ σ ′ G i σ ≶ ( t ) G^{\lessgtr}_{i\sigma,j\sigma'}(t)=\delta_{ij}\delta_{\sigma\sigma'}G^{\lessgtr}_{i\sigma}(t) G iσ , j σ ′ ≶ ( t ) = δ ij δ σ σ ′ G iσ ≶ ( t ) i , j = 1 , 2 i,j=1,2 i , j = 1 , 2 G p > ( t ) = − i ⟨ ψ ∣ a ^ p ( t ) a ^ p † ∣ ψ ⟩ G_{p}^{>}(t)=-i\langle\psi|\hat{a}_{p}(t)\hat{a}_{p}^{\dagger}|\psi\rangle G p > ( t ) = − i ⟨ ψ ∣ a ^ p ( t ) a ^ p † ∣ ψ ⟩ G p < ( t ) = i ⟨ ψ ∣ a ^ p † a ^ p ( t ) ∣ ψ ⟩ G_{p}^{<}(t)=i\langle\psi|\hat{a}_{p}^{\dagger}\hat{a}_{p}(t)|\psi\rangle G p < ( t ) = i ⟨ ψ ∣ a ^ p † a ^ p ( t ) ∣ ψ ⟩ G p R ( t ) = Θ ( t ) [ G p > ( t ) − G p < ( t ) ] G_{p}^{R}(t)=\Theta(t)[G_{p}^{>}(t)-G_{p}^{<}(t)] G p R ( t ) = Θ ( t ) [ G p > ( t ) − G p < ( t )] ∣ 2 ↑ 2 ↓ ⟩ |2 \uparrow2 \downarrow\rangle ∣2 ↑ 2 ↓ ⟩ ∣ 1 ↑ , 2 ↑ ⟩ |1 \uparrow, 2 \uparrow\rangle ∣1 ↑ , 2 ↑ ⟩ 0 0 0 U U U ϵ − J E \epsilon- J^{E} ϵ − J E
Table S96.
For the remaining 1 ↓ 1 \downarrow 1 ↓ 1 ↓ 1 \downarrow 1 ↓ ∣ ψ ⟩ = ∣ 2 ↑ ⟩ |\psi\rangle= |2 \uparrow\rangle ∣ ψ ⟩ = ∣2 ↑ ⟩ ∣ p ⟩ = ∣ 1 ↓ , 2 ↑ ⟩ |p\rangle= |1 \downarrow, 2 \uparrow\rangle ∣ p ⟩ = ∣1 ↓ , 2 ↑ ⟩ ∣ q ⟩ = ∣ 1 ↑ , 2 ↓ ⟩ |q\rangle= |1 \uparrow, 2 \downarrow\rangle ∣ q ⟩ = ∣1 ↑ , 2 ↓ ⟩ ∣ ψ ⟩ → a ^ 1 ↓ † ∣ p ⟩ ← − J E ∣ q ⟩ |\psi\rangle\xrightarrow{\hat{a}^{\dagger}_{1\downarrow}} |p\rangle\xleftarrow{\ \ -J^{E}\ \ } |q\rangle ∣ ψ ⟩ a ^ 1 ↓ † ∣ p ⟩ − J E ∣ q ⟩ ∣ p ⟩ |p\rangle ∣ p ⟩ ∣ q ⟩ |q\rangle ∣ q ⟩ ∣ T 0 ⟩ = ( ∣ p ⟩ + ∣ q ⟩ ) / 2 |T_{0}\rangle= (|p\rangle+ |q\rangle)/\sqrt{2} ∣ T 0 ⟩ = ( ∣ p ⟩ + ∣ q ⟩) / 2 ∣ S ⟩ = ( ∣ p ⟩ − ∣ q ⟩ ) / 2 |S\rangle= (|p\rangle- |q\rangle)/\sqrt{2} ∣ S ⟩ = ( ∣ p ⟩ − ∣ q ⟩) / 2 ϵ − J E \epsilon- J^{E} ϵ − J E ϵ + J E \epsilon+ J^{E} ϵ + J E a ^ 1 ↓ ∣ ψ ⟩ = 0 \hat{a}_{1\downarrow}|\psi\rangle= 0 a ^ 1 ↓ ∣ ψ ⟩ = 0 H ^ ∣ ψ ⟩ = 0 \hat{H}|\psi\rangle= 0 H ^ ∣ ψ ⟩ = 0 a ^ 1 ↓ † ∣ ψ ⟩ = ∣ p ⟩ \hat{a}^{\dagger}_{1\downarrow}|\psi\rangle= |p\rangle a ^ 1 ↓ † ∣ ψ ⟩ = ∣ p ⟩ G 1 ↓ R G^{R}_{1\downarrow} G 1 ↓ R { ∣ p ⟩ , ∣ q ⟩ } \{|p\rangle, |q\rangle\} { ∣ p ⟩ , ∣ q ⟩} ∣ p ⟩ = ( ∣ T 0 ⟩ + ∣ S ⟩ ) / 2 |p\rangle= (|T_{0}\rangle+ |S\rangle)/\sqrt{2} ∣ p ⟩ = ( ∣ T 0 ⟩ + ∣ S ⟩) / 2
G 1 ↓ > ( t ) = − i 2 [ e − i ( ϵ − J E ) t + e − i ( ϵ + J E ) t ] , G 1 ↓ < ( t ) = 0. (114) G^{>}_{1\downarrow}(t)=-\frac{\mathrm{i}}{2}\left[e^{-\mathrm{i}(\epsilon-J^{E})t}+e^{-\mathrm{i}(\epsilon+J^{E})t}\right], \qquad G^{<}_{1\downarrow}(t)=0. \tag*{(114)} G 1 ↓ > ( t ) = − 2 i [ e − i ( ϵ − J E ) t + e − i ( ϵ + J E ) t ] , G 1 ↓ < ( t ) = 0. (114) The Fourier transform of G 1 ↓ R ( t ) G^{R}_{1\downarrow}(t) G 1 ↓ R ( t )
G ~ 1 ↓ R ( z ) = 1 z − ϵ − ( J E ) 2 / ( z − ϵ ) = 1 2 [ 1 z − ϵ + J E + 1 z − ϵ − J E ] . (115) \widetilde{G}^{R}_{1\downarrow}(z)=\frac{1}{z-\epsilon-(J^{E})^{2}/(z-\epsilon)}=\frac{1}{2}\left[\frac{1}{z-\epsilon+J^{E}}+\frac{1}{z-\epsilon-J^{E}}\right]. \tag*{(115)} G 1 ↓ R ( z ) = z − ϵ − ( J E ) 2 / ( z − ϵ ) 1 = 2 1 [ z − ϵ + J E 1 + z − ϵ − J E 1 ] . (115) The 1 ↓ 1 \downarrow 1 ↓ 2 J E 2J^{E} 2 J E U U U G W GW G W
We use the interaction convention
H ^ i n t = 1 2 ∑ i j k l , σ τ C i j , k l a ^ i σ † a ^ k τ † a ^ l τ a ^ j σ , C 11 , 11 = C 22 , 22 = U , C 12 , 21 = C 21 , 12 = J E , (116) \hat{H}_{\mathrm{int}}=\frac{1}{2}\sum_{ijkl,\sigma\tau}C_{ij,kl}\hat{a}^{\dagger}_{i\sigma}\hat{a}^{\dagger}_{k\tau}\hat{a}_{l\tau}\hat{a}_{j\sigma},\qquad C_{11,11}=C_{22,22}=U,\qquad C_{12,21}=C_{21,12}=J^{E}, \tag*{(116)} H ^ int = 2 1 ij k l , σ τ ∑ C ij , k l a ^ iσ † a ^ k τ † a ^ l τ a ^ j σ , C 11 , 11 = C 22 , 22 = U , C 12 , 21 = C 21 , 12 = J E , (116) with all other entries zero. The spin-orbital interaction elements are V ( i σ ) ( k τ ) , ( j σ ′ ) ( l τ ′ ) = C i j , k l δ σ σ ′ δ τ τ ′ V_{(i\sigma)(k\tau),(j\sigma^{\prime})(l\tau^{\prime})}=C_{ij,kl}\delta_{\sigma\sigma^{\prime}}\delta_{\tau\tau^{\prime}} V ( iσ ) ( k τ ) , ( j σ ′ ) ( l τ ′ ) = C ij , k l δ σ σ ′ δ τ τ ′ V α β , γ δ = V α β , γ δ − V α β , δ γ V_{\alpha\beta,\gamma\delta}=V_{\alpha\beta,\gamma\delta}-V_{\alpha\beta,\delta\gamma} V α β , γ δ = V α β , γ δ − V α β , δ γ Σ ~ R ( z ) = Σ H F + Σ ~ c ( z ) \widetilde{\Sigma}^{R}(z)=\Sigma^{\mathrm{HF}}+\widetilde{\Sigma}^{c}(z) Σ R ( z ) = Σ HF + Σ c ( z ) Σ α β H F = ∑ γ V ‾ α γ , β γ ⟨ ψ ∣ n ^ γ ∣ ψ ⟩ = V ‾ α 2 ↑ , β 2 ↑ \Sigma^{\mathrm{HF}}_{\alpha\beta}=\sum_{\gamma}\overline{V}_{\alpha\gamma,\beta\gamma}\langle\psi|\hat{n}_{\gamma}|\psi\rangle=\overline{V}_{\alpha\,2\uparrow,\beta\,2\uparrow} Σ α β HF = ∑ γ V α γ , β γ ⟨ ψ ∣ n ^ γ ∣ ψ ⟩ = V α 2 ↑ , β 2 ↑ ( 1 ↑ , 1 ↓ , 2 ↑ , 2 ↓ ) (1 \uparrow, 1 \downarrow, 2 \uparrow, 2 \downarrow) ( 1 ↑ , 1 ↓ , 2 ↑ , 2 ↓ ) Σ H F = diag ( − J E , 0 , 0 , U ) \Sigma^{\mathrm{HF}}=\operatorname{diag}(-J^{\mathrm{E}},0,0,U) Σ HF = diag ( − J E , 0 , 0 , U ) ϵ 1 σ = ϵ \epsilon_{1\sigma}=\epsilon ϵ 1 σ = ϵ ϵ 2 σ = 0 \epsilon_{2\sigma}=0 ϵ 2 σ = 0 f p = ϵ p + Σ p p H F f_p=\epsilon_p+\Sigma^{\mathrm{HF}}_{pp} f p = ϵ p + Σ pp HF f 1 ↑ = ϵ − J E f_{1\uparrow}=\epsilon-J^{\mathrm{E}} f 1 ↑ = ϵ − J E f 1 ↓ = ϵ f_{1\downarrow}=\epsilon f 1 ↓ = ϵ f 2 ↑ = 0 f_{2\uparrow}=0 f 2 ↑ = 0 f 2 ↓ = U f_{2\downarrow}=U f 2 ↓ = U G ~ p R ( z ) = [ z − f p − Σ ~ p c ( z ) ] − 1 \widetilde{G}^{R}_{p}(z)=[z-f_p-\widetilde{\Sigma}^{c}_{p}(z)]^{-1} G p R ( z ) = [ z − f p − Σ p c ( z ) ] − 1 G ~ 1 ↓ R ( z ) = 1 z − ϵ − Σ ~ 1 ↓ c ( z ) \widetilde{G}^{R}_{1\downarrow}(z)=\frac{1}{z-\epsilon-\widetilde{\Sigma}^{c}_{1\downarrow}(z)} G 1 ↓ R ( z ) = z − ϵ − Σ 1 ↓ c ( z ) 1 Σ ~ c ( z ) \widetilde{\Sigma}^{c}(z) Σ c ( z ) (115) gives the exact correlation self-energy,
Σ ~ 1 ↓ c ( z ) = ( J E ) 2 z − ϵ . (117) \widetilde{\Sigma}^{c}_{1\downarrow}(z)=\frac{(J^{\mathrm{E}})^2}{z-\epsilon}. \tag*{(117)} Σ 1 ↓ c ( z ) = z − ϵ ( J E ) 2 . (117) Second Born. —For diagonal propagators G ≷ G^{\gtrless} G ≷ [128 , 129 ]
Σ α β c , ≷ ( t ) = Σ α β ( 2 ) , ≷ ( t ) = 1 2 ∑ γ μ ν V ‾ α γ , μ ν V ‾ μ ν , β γ G μ ≷ ( t ) G ν ≷ ( t ) G γ ≶ ( − t ) . (118) \Sigma^{c,\gtrless}_{\alpha\beta}(t)=\Sigma^{(2),\gtrless}_{\alpha\beta}(t)=\frac{1}{2}\sum_{\gamma\mu\nu}\overline{V}_{\alpha\gamma,\mu\nu}\overline{V}_{\mu\nu,\beta\gamma}G^{\gtrless}_{\mu}(t)G^{\gtrless}_{\nu}(t)G^{\lessgtr}_{\gamma}(-t). \tag*{(118)} Σ α β c , ≷ ( t ) = Σ α β ( 2 ) , ≷ ( t ) = 2 1 γ μν ∑ V α γ , μν V μν , β γ G μ ≷ ( t ) G ν ≷ ( t ) G γ ≶ ( − t ) . (118) For the greater component, G γ < G^{<}_{\gamma} G γ < γ = 2 ↑ \gamma=2\uparrow γ = 2 ↑ γ = 2 ↑ \gamma=2\uparrow γ = 2 ↑ ( μ , ν ) (\mu,\nu) ( μ , ν ) ( 1 ↑ , 2 ↓ ) (1\uparrow,2\downarrow) ( 1 ↑ , 2 ↓ ) ( 2 ↓ , 1 ↑ ) (2\downarrow,1\uparrow) ( 2 ↓ , 1 ↑ ) − J E -J^{\mathrm{E}} − J E + J E +J^{\mathrm{E}} + J E
Σ 1 ↓ ( 2 ) , > ( t ) = ( J E ) 2 G 1 ↑ > ( t ) G 2 ↓ > ( t ) G 2 ↑ < ( − t ) = − i ( J E ) 2 e − i ( ϵ + U − J E ) t . (119) \Sigma^{(2),>}_{1\downarrow}(t)=(J^{\mathrm{E}})^2G^{>}_{1\uparrow}(t)G^{>}_{2\downarrow}(t)G^{<}_{2\uparrow}(-t)=-i(J^{\mathrm{E}})^2e^{-i(\epsilon+U-J^{\mathrm{E}})t}. \tag*{(119)} Σ 1 ↓ ( 2 ) , > ( t ) = ( J E ) 2 G 1 ↑ > ( t ) G 2 ↓ > ( t ) G 2 ↑ < ( − t ) = − i ( J E ) 2 e − i ( ϵ + U − J E ) t . (119) Every other greater entry vanishes by the same interaction-vertex selection. All lesser entries require ( μ , ν ) = ( 2 ↑ , 2 ↑ ) (\mu,\nu)=(2\uparrow,2\uparrow) ( μ , ν ) = ( 2 ↑ , 2 ↑ ) V ‾ \overline{V} V Σ c ( t ) = Θ ( t ) [ Σ c , > ( t ) − Σ c , < ( t ) ] = Θ ( t ) Σ ( 2 ) , > ( t ) \Sigma^{c}(t)=\Theta(t)[\Sigma^{c,>}(t)-\Sigma^{c,<}(t)]=\Theta(t)\Sigma^{(2),>}(t) Σ c ( t ) = Θ ( t ) [ Σ c , > ( t ) − Σ c , < ( t )] = Θ ( t ) Σ ( 2 ) , > ( t )
Σ ~ 1 ↓ , 2 n d B c ( z ) = ( J E ) 2 z − ( ϵ + U − J E ) . (120) \widetilde{\Sigma}^{c}_{1\downarrow,\mathrm{2ndB}}(z)=\frac{(J^{\mathrm{E}})^2}{z-(\epsilon+U-J^{\mathrm{E}})}. \tag*{(120)} Σ 1 ↓ , 2ndB c ( z ) = z − ( ϵ + U − J E ) ( J E ) 2 . (120) Substituting this closed-form self-energy into Dyson’s equation gives the 2ndB Green’s function,
G ~ 1 ↓ , 2 n d B R ( z ) = 1 z − ϵ − ( J E ) 2 / [ z − ( ϵ + U − J E ) ] . (121) \widetilde{G}^{R}_{1\downarrow,\mathrm{2ndB}}(z)=\frac{1}{z-\epsilon-(J^{\mathrm{E}})^2/[z-(\epsilon+U-J^{\mathrm{E}})]}. \tag*{(121)} G 1 ↓ , 2ndB R ( z ) = z − ϵ − ( J E ) 2 / [ z − ( ϵ + U − J E )] 1 . (121) Both its poles and residues depend on U U U (115) . The physical picture of this failure is clear. The intermediate configuration is ∣ q ⟩ = a ^ 1 ↑ † a ^ 2 ↓ † a ^ 2 ↑ ∣ ψ ⟩ |q\rangle=\hat{a}^{\dagger}_{1\uparrow}\hat{a}^{\dagger}_{2\downarrow}\hat{a}_{2\uparrow}|\psi\rangle ∣ q ⟩ = a ^ 1 ↑ † a ^ 2 ↓ † a ^ 2 ↑ ∣ ψ ⟩ (119) and Figure S6 ). However, the two J E J^{\mathrm{E}} J E ∣ q ⟩ |q\rangle ∣ q ⟩ ∣ p ⟩ |p\rangle ∣ p ⟩ f 1 ↑ + f 2 ↓ − f 2 ↑ = ϵ + U − J E f_{1\uparrow}+f_{2\downarrow}-f_{2\uparrow}=\epsilon+U-J^{\mathrm{E}} f 1 ↑ + f 2 ↓ − f 2 ↑ = ϵ + U − J E (120) . Adding 2 ↓ 2\downarrow 2 ↓ U U U 2 ↑ 2\uparrow 2 ↑ 1 ↑ 1\uparrow 1 ↑ − J E -J^{\mathrm{E}} − J E ∣ q ⟩ |q\rangle ∣ q ⟩ 2 ↑ 2\uparrow 2 ↑ U U U − J E -J^{\mathrm{E}} − J E
Figure S6. Self-energy diagrams for 1 ↓ addition. Solid arrows denote bold propagators, a single wavy line denotes a bare interaction, and the double wavy line denotes the correlation part of the screened interaction, W c = W − C W^{c}=W-C W c = W − C 4 + J E +J^{E} + J E − U -U − U G W GW G W W c W^{c} W c + J E +J^{E} + J E U U U G W GW G W
GW. —The screening approximation uses [128 , 129 ]
Σ i j , σ c , ≷ ( t ) = i ∑ k l G k l , σ ≷ ( t ) W i k , l j ≶ ( t ) . (122) \Sigma^{c,\gtrless}_{ij,\sigma}(t)=i\sum_{kl}G^{\gtrless}_{kl,\sigma}(t)W^{\lessgtr}_{ik,lj}(t). \tag*{(122)} Σ ij , σ c , ≷ ( t ) = i k l ∑ G k l , σ ≷ ( t ) W ik , l j ≶ ( t ) . (122) Here W W W W c W^{c} W c P P P
P ~ i j , k l ≷ ( t ) = − i ∑ σ G j k , σ ≷ ( t ) G l i , σ ≶ ( − t ) , W = C + C ∗ P ∗ W , W c = W − C . (123) \widetilde{P}^{\gtrless}_{ij,kl}(t)=-i\sum_{\sigma}G^{\gtrless}_{jk,\sigma}(t)G^{\lessgtr}_{li,\sigma}(-t),\qquad W=C+C*P*W,\qquad W^{c}=W-C. \tag*{(123)} P ij , k l ≷ ( t ) = − i σ ∑ G j k , σ ≷ ( t ) G l i , σ ≶ ( − t ) , W = C + C ∗ P ∗ W , W c = W − C . (123) Here ∗ * ∗ C C C [129 ] . The screening equation becomes W ~ R = C + C ~ P ~ R W ~ R \widetilde{W}^{R}=C+\widetilde{C}\widetilde{P}^{R}\widetilde{W}^{R} W R = C + C P R W R
( C ~ P ~ R W ~ R ) i j , k l ( ω ) = ∑ a b , c d C ~ i j , a b P ~ a b , c d R ( ω ) W ~ c d , k l R ( ω ) . (124) (\widetilde{C}\widetilde{P}^{R}\widetilde{W}^{R})_{ij,kl}(\omega)=\sum_{ab,cd}\widetilde{C}_{ij,ab}\widetilde{P}^{R}_{ab,cd}(\omega)\widetilde{W}^{R}_{cd,kl}(\omega). \tag*{(124)} ( C P R W R ) ij , k l ( ω ) = ab , c d ∑ C ij , ab P ab , c d R ( ω ) W c d , k l R ( ω ) . (124) There are only two nonzero polarization entries, associated with 2 ↑ ↔ 1 ↑ 2\uparrow\leftrightarrow1\uparrow 2 ↑↔ 1 ↑
P ~ 21 , 12 R ( ω ) = 1 ω − ( ϵ − J E ) , P ~ 12 , 21 R ( ω ) = − 1 ω + ( ϵ − J E ) , Im ω > 0. (125) \widetilde{P}^{R}_{21,12}(\omega)=\frac{1}{\omega-(\epsilon-J^{\mathrm{E}})},\qquad\widetilde{P}^{R}_{12,21}(\omega)=-\frac{1}{\omega+(\epsilon-J^{\mathrm{E}})},\qquad\operatorname{Im}\omega>0. \tag*{(125)} P 21 , 12 R ( ω ) = ω − ( ϵ − J E ) 1 , P 12 , 21 R ( ω ) = − ω + ( ϵ − J E ) 1 , Im ω > 0. (125) Thus, for C = C U + C J E C=C_{U}+C_{J^{E}} C = C U + C J E C U P = P C U = 0 C_{U}P=PC_{U}=0 C U P = P C U = 0 (116) . The two surviving screening equations are
W ~ 12 , 21 R ( ω ) = J E + J E ω − ( ϵ − J E ) W ~ 12 , 21 R ( ω ) , W ~ 21 , 12 R ( ω ) = J E − J E ω + ( ϵ − J E ) W ~ 21 , 12 R ( ω ) . (126) \widetilde{W}^{R}_{12,21}(\omega)=J^{E}+\frac{J^{E}}{\omega-(\epsilon-J^{E})}\widetilde{W}^{R}_{12,21}(\omega),\qquad\widetilde{W}^{R}_{21,12}(\omega)=J^{E}-\frac{J^{E}}{\omega+(\epsilon-J^{E})}\widetilde{W}^{R}_{21,12}(\omega). \tag*{(126)} W 12 , 21 R ( ω ) = J E + ω − ( ϵ − J E ) J E W 12 , 21 R ( ω ) , W 21 , 12 R ( ω ) = J E − ω + ( ϵ − J E ) J E W 21 , 12 R ( ω ) . (126) Subtracting the bare interaction C C C
W ~ 12 , 21 c , R ( ω ) = ( J E ) 2 ω − ϵ , W ~ 21 , 12 c , R ( ω ) = − ( J E ) 2 ω + ϵ , W 12 , 21 c , > ( t ) = − i ( J E ) 2 e − i ϵ t . (127) \widetilde{W}^{c,R}_{12,21}(\omega)=\frac{(J^{E})^{2}}{\omega-\epsilon},\qquad\widetilde{W}^{c,R}_{21,12}(\omega)=-\frac{(J^{E})^{2}}{\omega+\epsilon},\qquad W^{c,>}_{12,21}(t)=-i(J^{E})^{2}e^{-i\epsilon t}. \tag*{(127)} W 12 , 21 c , R ( ω ) = ω − ϵ ( J E ) 2 , W 21 , 12 c , R ( ω ) = − ω + ϵ ( J E ) 2 , W 12 , 21 c , > ( t ) = − i ( J E ) 2 e − i ϵ t . (127) The only nonzero contribution to Eq. (122) closes W c W^{c} W c Figure S6 .
Σ 1 ↓ c , > ( t ) = i G 2 ↓ > ( t ) W 12 , 21 c , > ( t ) = − i ( J E ) 2 e − i ( ϵ + U ) t , Σ 1 ↓ c , < ( t ) = 0 , Σ ~ 1 ↓ , G W c ( z ) = ( J E ) 2 z − ( ϵ + U ) . (128) \Sigma^{c,>}_{1\downarrow}(t)=iG^{>}_{2\downarrow}(t)W^{c,>}_{12,21}(t)=-i(J^{E})^{2}e^{-i(\epsilon+U)t},\qquad\Sigma^{c,<}_{1\downarrow}(t)=0,\qquad\widetilde{\Sigma}^{c}_{1\downarrow,GW}(z)=\frac{(J^{E})^{2}}{z-(\epsilon+U)}. \tag*{(128)} Σ 1 ↓ c , > ( t ) = i G 2 ↓ > ( t ) W 12 , 21 c , > ( t ) = − i ( J E ) 2 e − i ( ϵ + U ) t , Σ 1 ↓ c , < ( t ) = 0 , Σ 1 ↓ , G W c ( z ) = z − ( ϵ + U ) ( J E ) 2 . (128) Similar to Eq. (121) , G W GW G W ϵ + U \epsilon+U ϵ + U ϵ \epsilon ϵ U U U ϵ \epsilon ϵ ∣ 2 ↑ ⟩ |2\uparrow\rangle ∣2 ↑ ⟩ ∣ 1 ↑ ⟩ |1\uparrow\rangle ∣1 ↑ ⟩ G W GW G W U U U J E J^{E} J E U U U ∣ q ⟩ |q\rangle ∣ q ⟩ (117) requires interactions between the 2 ↓ particle and the 2 ↑ hole, as illustrated in Figure S6 .
Consequences for the multiplet spectral features.— We discuss the consequence of the missing qualitative interactions on the spectral features. Since the other three propagators retain their HF poles and updating G 1 ↓ G_{1\downarrow} G 1 ↓
G ~ 1 ↓ R ( z ) = 1 z − ϵ − ( J E ) 2 / ( z − ϵ − a ) = ( z I − ( ϵ − J E J E − J E ϵ + a ) ) 11 − 1 , (129) \widetilde{G}^{R}_{1\downarrow}(z)=\frac{1}{z-\epsilon-(J^{E})^{2}/(z-\epsilon-a)}=\left(zI-\begin{pmatrix}\epsilon-J^{E}&J^{E}\\-J^{E}&\epsilon+a\end{pmatrix}\right)^{-1}_{11}, \tag*{(129)} G 1 ↓ R ( z ) = z − ϵ − ( J E ) 2 / ( z − ϵ − a ) 1 = ( z I − ( ϵ − J E − J E J E ϵ + a ) ) 11 − 1 , (129) where a = 0 a=0 a = 0 a = U − J E a=U-J^{E} a = U − J E a = U a=U a = U G W GW G W U U U G ~ 1 ↓ R \widetilde{G}^{R}_{1\downarrow} G 1 ↓ R
ω ± = ϵ + a ± a 2 + 4 ( J E ) 2 2 , Res ω ± G ~ 1 ↓ R = 1 2 ( 1 ∓ a a 2 + 4 ( J E ) 2 ) . (130) \omega_{\pm}=\epsilon+\frac{a\pm\sqrt{a^{2}+4(J^{E})^{2}}}{2},\qquad\operatorname{Res}_{\omega_{\pm}}\widetilde{G}^{R}_{1\downarrow}=\frac{1}{2}\left(1\mp\frac{a}{\sqrt{a^{2}+4(J^{E})^{2}}}\right). \tag*{(130)} ω ± = ϵ + 2 a ± a 2 + 4 ( J E ) 2 , Res ω ± G 1 ↓ R = 2 1 ( 1 ∓ a 2 + 4 ( J E ) 2 a ) . (130) Figure S7 illustrates the resulting responses at ϵ = 6 \epsilon=6 ϵ = 6 U = 4 U=4 U = 4 J E = 0.5 J^{E}=0.5 J E = 0.5
![Im G ~ 1 ↓ R ( ω + i η ) \widetilde{G}^{R}_{1\downarrow}(\omega+i\eta) G 1 ↓ R ( ω + i η )
Figure S7. Im G ~ 1 ↓ R ( ω + i η ) \widetilde{G}^{R}_{1\downarrow}(\omega+i\eta) G 1 ↓ R ( ω + i η ) G W GW G W ϵ = 6 \epsilon=6 ϵ = 6 U = 4 U=4 U = 4 J E = 0.5 J^{\mathrm{E}}=0.5 J E = 0.5 η = 0.1 \eta=0.1 η = 0.1 ω − ϵ = ± J E \omega-\epsilon=\pm J^{\mathrm{E}} ω − ϵ = ± J E G W GW G W
Recovering the exact self-energy pole requires correlations among the different pairs of the 2 p 1 h 2p1h 2 p 1 h ↑ \uparrow ↑ ↑ \uparrow ↑ ↓ \downarrow ↓ ↑ \uparrow ↑ 2 p 1 h 2p1h 2 p 1 h + J E +J^{\mathrm{E}} + J E − U -U − U ϵ \epsilon ϵ [128 , 130 , 131 ] . Hamiltonian EOM-CCSD accommodates those interactions [89 ] by explicitly including both ∣ p ⟩ \lvert p\rangle ∣ p ⟩ ∣ q ⟩ \lvert q\rangle ∣ q ⟩ 1 p / 2 p 1 h 1p/2p1h 1 p /2 p 1 h
References 
