Physics Notes
Implementation paths and commands in this guide are relative to
src/pairing/ from the repository root.
This solver implements finite-temperature grand-canonical DQMC for a two-flavor Bose-Hubbard model with onsite pairing on a triangular lattice. The active triangular implementation has one orbital per cell, so \(N_s=L_xL_y\).
Hamiltonian
The paper's main-text model has the physical Hamiltonian
where \(\hat n_{b,i}=\hat b_i^+\hat b_i\), \(\hat n_{c,i}=\hat c_i^+\hat c_i\), \(t>0\), and \(U\ge0\). The trace is grand canonical:
The code retains two density-interaction channels,
The main-text benchmark sets \(U_1=0\) and \(U_2=U\). The paired scan has \(U=1\),
\(\mu=-5\), \(\beta=4\), and \(\Delta=0,0.05,\ldots,0.30\); its four observables
form panels (e–h) of the combined benchmark figure. The code reads \(U_1\),
\(U_2\), \(\mu\), and a real \(\Delta\) from paramC_sets.txt; \(t\) is RT=1 in
src/calc_basic.f90. Both BAFQMC and ED use these same input coefficients.
The interaction is equivalently
on each site, as used in the ED construction and energy estimator.
Manuscript notation and pair phase
For the real pairing coefficient used in the calculations, the implementation uses the flavor phase convention
In these code operators, the same physical pair term is
All implementation derivations below use the code operators. Their hats
and code subscripts are omitted for readability, including in the Nambu
basis, Green-function block table, HS scalar factor, and ED identities.
The positive pair coefficient in those formulas is the code representation
of the paper's negative pair term.
Hopping and density interactions are unchanged by this phase rotation. The four plotted observables \(\rho\), \(-E\), \(S_{\rm SF}(K)\), and \(S_{\rm DW}(K)\) therefore use the stored values directly. The anomalous amplitude changes sign:
The total pair energy is \(E_\Delta=N_s e_\Delta\). The physical energy excludes \(-\mu(\hat N_b+\hat N_c)\) throughout both implementations.
For this \(U_1=0\) model, the condition \(\mu<-3t-|\Delta|\) guarantees a finite trace for every real auxiliary-field configuration, as proved in the SM subsection "Convergence throughout the auxiliary-field domain". All seven paired benchmark points satisfy it. Only the relative-density HS field is coupled to the bosons, so the total-density normal-ordering factor is inactive.
Triangular Lattice Vectors And Bonds
The Hamiltonian's triangular lattice is defined by integer cell topology. The spatial cell index is \((x,y)\), with periodic boundaries, and each site stores three forward bonds:
This is a standard triangular nearest-neighbor graph. The pairing and number-conserving solvers use the same real-space and reciprocal-space convention:
The three forward displacements are \(\mathbf a_1\), \(\mathbf a_2\), and
\(-\mathbf a_1+\mathbf a_2\), all with length \(1\). The kinetic estimator
uses the integer bond table and adds both directions explicitly. In
src/lattice.f90, the zero-based displacement labels are used to populate
with zero-based momentum indices \(m=x-1,n=y-1\), and k_dot_r is the
ordinary dot product. For \(3\times3\), the code-facing \(K\) index
Latt%inv_cell_list(2*Nlx/3+1, Nly/3+1) represents
which is one of the triangular-lattice single-particle band minima for the current hopping sign convention.
Continuous Hubbard-Stratonovich Fields
The Trotter step uses two continuous Gaussian fields \(\phi_{1,i,\tau}\) and \(\phi_{2,i,\tau}\). For the intended signs \(U_1\le0\) and \(U_2\ge0\),
The auxiliary fields live on spatial sites, not Nambu sectors:
Conf%phi_list(ns, i_site, ntau), where ns=1 is the \(U_1\) channel and ns=2 is the \(U_2\) channel.
The paper defines \(\alpha_2\) as real; OperatorHubbard%alpha stores the complex coefficient \(i\alpha_2\) for this channel. At zero coupling the corresponding coefficient vanishes.
Nambu Basis Convention
The active pairing implementation uses the four-component bosonic Nambu basis
with
For the one-body triangular hopping, chemical potential, and onsite pairing piece, the commutator matrix has the block form
Here \(T\) is the normal one-body matrix, including the hopping and chemical
potential terms in the code convention. This is the matrix convention used by
src/non_interact.f90 and by the HS update matrices. It is not the reduced
two-component convention
That reduced convention may be a future optimization, but it is not the current implementation. Do not switch formulas between the two conventions without rederiving the determinant weight, Green matrix, and observable block table.
Although the physical Hamiltonian is Hermitian, the commutator matrix \(A\) in this bosonic Nambu representation is generally non-Hermitian when \(\Delta\neq 0\). The one-body propagator must therefore be built as a general matrix exponential of \(A\), not by a Hermitian diagonalization of \(A\). For a single normal mode with normal coefficient \(h\), the paired Nambu block has eigenvalues
in the stable regime. A Hermitian diagonalization would instead produce the
wrong \(\sqrt{h^2+\Delta^2}\) scale and gives a visibly biased finite-\(\Delta\)
benchmark even in the noninteracting limit. src/non_interact.f90 therefore
uses LAPACK zgeev to diagonalize the general complex matrix before forming
\(\exp(\mp\Delta\tau A)\).
HS Constant In The Nambu Convention
Let \(\lambda_1\) and \(\lambda_2\) denote the actual complex coefficients that multiply the normal-ordered density channels on one space-time site:
In the Fortran code, OperatorHubbard%alpha stores \(i\alpha_2\), so its
product with the real field is \(\lambda_2\). The HS density term satisfies
The final \(-\lambda_1\) is the bosonic commutator c-number. It does not appear as a matrix element in the Nambu propagator, but it must appear in the sampling weight. Therefore a local \(U_1\) update contributes
to ratio_constant. In src/operator_Hubbard.f90 this is implemented as
expalpha_old / expalpha_new for IUflag == 1. The \(U_2\) channel has no
scalar c-number because its two flavor coefficients sum to zero.
Green Matrix And Local Determinant Factor
For one full imaginary-time product, define
The physical Green matrix is \(G\). The quantity propagated and updated in the Fortran code is instead
Despite the historical variable name, Prop%Gr stores \(\tilde G\), not the
physical \(G\). The local Metropolis probability is
where
The determinant factor is the current Nambu determinant contribution. In
src/localU.f90, the local update builds a doubled-sector \(4\times4\) block
\(P\) and evaluates
The exponent \(-1/2\) is part of the present doubled Nambu convention. Do not change it to \(-1\) or \(-1/4\) without a full rederivation and benchmark update.
Program Green-Function Block Table
The raw Prop%Gr matrix is \(\tilde G\). In the sector order
\((b,c,b^+,c^+)\), its blocks represent
src/obser_equal.f90 converts these raw blocks into physical contractions by
undoing the explicit minus signs in the table. For example,
\(G_{b,c}=-\texttt{Prop\%Gr(b,cdag)}\) and
\(G_{b^+,b}=-\texttt{Prop\%Gr(bdag,bdag)}\).
Define
These are the contractions used directly in physical normal-ordered
observables. In particular, cross-flavor products use the two sectors that are
present in the explicit Nambu matrix, \(n_b n_c\). They should not be replaced
by an ad hoc conjugated product when evaluating doubleOcc.
Observable Estimator Formulas
The scalar equal-time estimators are accumulated per configuration and then averaged over observations, MPI ranks, and bins. The density and number estimators are
The kinetic output is per spatial site and includes the hopping coefficient:
The Wick contractions used by D_bb, D_cc, and D_bc are
The occupancy and pairing outputs are
squareOcc preserves the density-product estimator above. In the current
model, conservation of \(N_b-N_c\) removes same-flavor anomalous contractions,
so its physical meaning is half the per-site normal-ordered same-flavor
onsite pair density. local_numsquare includes the full Wick contractions
for \(N_s^{-1}\sum_i\langle n_{b,i}^2+n_{c,i}^2\rangle\) and is compared to ED.
See the observable reference for the current output
inventory and the corresponding distinction for more general paired models.
Delta=0 And Corrected No-Pairing doubleOcc
The no-pairing reference code samples only the b determinant and reconstructs
the c sector by complex conjugation. In
src/number_conserving/src/obser_equal.f90, the relevant definitions are
Therefore
and the corrected no-pairing output named doubleOcc is, configuration by
configuration, the physical cross-flavor estimator
This is the same object as the explicit-Nambu doubleOcc at \(\Delta=0\),
where the anomalous blocks vanish and the c normal sector is the complex
conjugate of the b sector. Therefore Benchmark A compares doubleOcc
directly between the two repositories:
against the corrected no-pairing doubleOcc. Older no-pairing output produced
with Grdoc = dconjg(transpose(Grdo)) - ZKRON should not be used as a
reference, because that bug instead gave \(\texttt{Grdoc}=\texttt{Grupc}\) and
turned doubleOcc into \(\mathrm{Re}[\texttt{Grupc}^2]\).
For density correlations, den_upup_sub11, den_dodo_sub11, and den_updo are the \(k=0\) Fourier components of \(D_{bb}\), \(D_{cc}\), and \(D_{bc}\) after the real-space accumulation by separation \(r_i-r_j\).
ED Finite-Temperature Trace
Because \(H_\Delta\) changes total particle number, finite-\(\Delta\) ED uses the grand-canonical Hamiltonian directly rather than fixed-\(N_b,N_c\) blocks. The general ED script builds a no-symmetry boson_basis_general basis over \(2N_s\) boson modes with local cutoff nmax and optional total cutoff ncut.
For the full eigensystem \(H_\mu|\alpha\rangle=E_\alpha|\alpha\rangle\),
The implementation shifts energies by \(E_0=\min_\alpha E_\alpha\) for numerical stability and reports
The ED identity tests check
Both solvers use the same benchmark energy convention. The paper-facing plotted energy is the total physical/internal energy
It excludes the chemical-potential term used in the grand-canonical sampling Hamiltonian. The DQMC and ED JSON files store the per-site value
so report and plotting scripts must multiply energy_density by \(N_s\)
before drawing \(-E\). The chemical-potential contribution and the diagnostic
grand-canonical energy density remain per-site diagnostics:
DQMC/ED Normalization Conversions
| output file | DQMC estimator | DQMC normalization | ED observable | comparison operation | Benchmark A | Benchmark B |
|---|---|---|---|---|---|---|
num_up |
\(\sum_i \mathrm{Re}\,n_b(i,i)\) | total | \(N_b\) | mean |
yes | yes |
num_do |
\(\sum_i \mathrm{Re}\,n_c(i,i)\) | total | \(N_c\) | mean |
yes | yes |
density_up |
num_up / Nsite |
per site | \(N_b/N_s\) | mean |
yes | yes |
density_do |
num_do / Nsite |
per site | \(N_c/N_s\) | mean |
yes | yes |
density |
density_up + density_do |
per site | \((N_b+N_c)/N_s\) | mean |
yes | yes |
density_total |
alias of density |
per site | density_total |
mean |
yes | yes |
density_site_total |
\(\rho_i=\mathrm{Re}[G_{b^+ b}(i,i)+G_{c^+ c}(i,i)]\) | site vector | IPR postprocessing | fixed-width row | no | yes |
kinetic |
hopping expectation including RT |
per site | ED JSON kinetic_total = total \(H_t\) |
mean_times_nsite |
yes | yes |
doubleOcc |
\(N_s^{-1}\sum_i \mathrm{Re}\,D_{bc}(i,i)\) | per site | \(N_s^{-1}\sum_i n_{b,i}n_{c,i}\) | mean |
yes | yes |
squareOcc |
\(N_s^{-1}\sum_i \mathrm{Re}[n_b(i,i)^2+n_c(i,i)^2]\) | per site | half normal-ordered same-flavor onsite pair density in the current model | mean |
yes | no |
local_numsquare |
\(N_s^{-1}\sum_i \mathrm{Re}[D_{bb}(i,i)+D_{cc}(i,i)]\) | per site | \(N_s^{-1}\sum_i(n_{b,i}^2+n_{c,i}^2)\) | mean |
yes when reference exists | yes |
onsite_n2_up |
\(\sum_i\mathrm{Re}\,D_{bb}(i,i)\) | total | onsite_n2_up |
mean |
yes when reference exists | yes |
onsite_n2_do |
\(\sum_i\mathrm{Re}\,D_{cc}(i,i)\) | total | onsite_n2_do |
mean |
yes when reference exists | yes |
numsquare_up |
\(\sum_{i,j}\mathrm{Re}\,D_{bb}(i,j)\) | total | \(N_b^2\) | mean |
yes | yes |
numsquare_do |
\(\sum_{i,j}\mathrm{Re}\,D_{cc}(i,j)\) | total | \(N_c^2\) | mean |
yes | yes |
pair_equal |
\(N_s^{-1}\sum_i\mathrm{Re}\langle b_i c_i+b_i^+ c_i^+\rangle\) | per site | \(N_s^{-1}\partial F/\partial\Delta\) | mean |
zero check | yes |
interaction_energy_density |
\((U_1+U_2)(M_b^{(2)}+M_c^{(2)})/N_s+2(U_1-U_2)\texttt{doubleOcc}\) | per site | interaction_energy_density |
mean |
no | yes |
pairing_energy_density |
\(\Delta\,\texttt{pair\_equal}\) | per site | pairing_energy_density |
mean |
zero check | yes |
energy_density |
\(e=E/N_s=e_t+e_U+e_\Delta\), excluding \(-\mu N\); paper plots use \(-E=-N_s e\) | per site storage | energy_density |
mean |
no | yes |
chemical_energy_density |
\(-\mu\rho\) | per site | chemical_energy_density |
mean |
diagnostic | diagnostic |
grand_energy_density |
energy_density + chemical_energy_density |
per site | grand_energy_density |
mean |
diagnostic | diagnostic |
sf_K |
\(N_s^{-2}\sum_{ij}e^{iK\cdot(r_i-r_j)}\langle b_i^+ b_j+c_i^+ c_j\rangle\) | complex scalar | S_SF_K |
real/imag block stats | no | yes |
dw_K |
\(N_s^{-2}\sum_{ij}e^{iK\cdot(r_i-r_j)}[D_{bb}(i,j)+D_{cc}(i,j)+D_{bc}(i,j)+D_{cb}(i,j)]\) | complex scalar | S_DW_K |
real/imag block stats | no | yes |
psf_Gamma |
\(N_s^{-2}\sum_{ij}\langle b_i^+ c_i^+ c_j b_j\rangle\) | complex scalar | S_PSF_Gamma |
real/imag block stats | no | yes |
Here
For \(K\)-point observables the pairing code uses the same \(K\) index and
phase convention as the number-conserving solver when both lattice dimensions are
multiples of three. ED entrypoints reject incompatible \(K\) labels, and
campaign analysis must mark sf_K and dw_K unavailable rather than treating
placeholders or zeros as physics data.
At compatible sizes sf_K and dw_K are standard Hermitian structure-factor
channels:
Their exact thermal averages are real and non-negative. DQMC still writes two
columns because the estimator is accumulated as a complex Fourier sum; the
real part is the comparison channel and a statistically significant imaginary
part marks the observable unreliable. The same imaginary-part reliability rule
applies to the \(\Gamma\)-point psf_Gamma.
The IPR computed from density_site_total is a DQMC spatial-uniformity
diagnostic. It is useful for spotting broken site indexing or severe sampling
inhomogeneity, but it should not be used as a strict ED agreement gate in
low-density pilots because the ratio
\(\sum_i \rho_i^2/(\sum_i\rho_i)^2\) has a finite-sample upward bias.
The kinetic naming contract is:
DQMC file `kinetic`: per spatial site.
ED JSON `kinetic_total`: total hopping expectation including the coefficient `t`.
DQMC-vs-ED and DQMC-vs-DQMC comparison operation: `mean_times_nsite`.