Exploring the sign problem

In this notebook, we use pyALF to investigate the average sign of the square-lattice Hubbard model at different fillings, temperatures, lattice sizes and Hubbard-Stratonovich (HS) decompositions.

7.2.1Hubbard Hamiltonian¶

The Hamiltonian with a chemical potential $\mu$ is written as

\hat{H} = -t\sum_{\langle i,j\rangle,\sigma} \left( \hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma} +\text{H.c.} \right) -\mu\sum_i \hat{n}_i +\hat{H}_U,

(1)

where $\hat{n}_i=\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}$ . Following the notation of the ALF documentation, the interaction can be written in the charge-coupled form,

\hat{H}_U^{\mathrm{charge}} = \frac{U}{2}\sum_i \left(\hat{n}_i-1\right)^2,

(2)

or in the $M_z$ -coupled form,

\hat{H}_U^{M_z} = -\frac{U}{2}\sum_i \left( \hat{n}_{i\uparrow}-\hat{n}_{i\downarrow} \right)^2.

(3)

For the spin- $1/2$ Hubbard model, the two expressions differ only by an additive constant and therefore describe the same physics. In ALF, the charge and $M_z$ representations are selected using Mz=False and Mz=True, respectively. Although the physical Hamiltonian is the same, the two choices lead to different auxiliary-field representations and can therefore have different sign or phase properties.

We investigate how the average reweighting factor changes with chemical potential, inverse temperature, lattice size, and the choice of Hubbard-Stratonovich decomposition.

7.2.2Preparing pyALF¶

We first import the pyALF classes and python libraries required to obtain the ALF source code and set up the simulations and later analysis.

from py_alf import ALF_source, Simulation

import pandas as pd

Next, we specify the location of the ALF source code:


alf_src = ALF_source(
    alf_dir="./ALF",
    branch="master",
)

print("ALF source:", alf_src.alf_dir)

Here, alf_dir points to the directory containing the ALF source code. If this directory does not yet exist, pyALF automatically downloads ALF from its Git repository and stores it at the specified location. We then compile ALF for the Hubbard model.

compile_sim = Simulation(alf_src=alf_src, ham_name='Hubbard',sim_dict={}) #sim_dict is used to choose parameters.
compile_sim.compile()

7.2.3Simulating the Hubbard model¶

We simulate the Hubbard model at $U/t=4$ for inverse temperatures $\beta t=1$ , 2, and 4 on square lattices of size $L\times L$ , with $L=4$ and 6. For every parameter set, we consider both the $M_z$ -coupled and charge-coupled Hubbard-Stratonovich decompositions. These are selected in ALF using Mz=True and Mz=False, respectively.

To investigate the effect of doping, we vary the chemical potential from $\mu/t=0$ to 2.5. At $\mu=0$ , the model is particle-hole symmetric and corresponds to half filling, where the simulations are free of the sign problem. Increasing $\mu$ moves the system away from half filling, where a sign problem can occur.

The imaginary-time discretization is fixed to $\Delta\tau\,t=0.1$ for all simulations.

7.2.3.1Preparing the simulations¶

sims = []

betas = [1.0, 2.0,4.0]
lattice_sizes = [4,6]
hs_decompositions = [True, False]
chemical_potentials = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5]

dtau = 0.1

for L in lattice_sizes:
    for beta in betas:
        for mz in hs_decompositions:
            for mu in chemical_potentials:

                hs_name = "Mz" if mz else "density"

                print(
                    f"Creating simulation: L={L}, beta={beta}, "
                    f"mu={mu}, HS={hs_name}, dtau={dtau}"
                )

                sim = Simulation(
                    alf_src,
                    "Hubbard",
                    {
                        "Lattice_type": "Square",
                        "L1": L,
                        "L2": L,
                        "Checkerboard": False,
                        "Symm": True,
                        "ham_T": 1.0,
                        "ham_U": 4.0,
                        "ham_chem": mu,
                        "ham_Tperp": 0.0,
                        "beta": beta,
                        "Ltau": 0,
                        "NSweep": 100,
                        "NBin": 10,
                        "Dtau": dtau,
                        "Mz": mz,
                    },
                    sim_root="ALF_data/sign_problem",
                    sim_dir=(
                        f"L{L}_beta{beta:g}_mu{mu:g}_{hs_name}"
                    ),
                    machine="GNU",
                )

                sims.append(sim)


print(f"\nCreated {len(sims)} simulations.")

The nested loop constructs one Simulation object for every combination of lattice size, inverse temperature, chemical potential, and HS decomposition. At this stage, the simulations are only configured; no Monte Carlo calculation has been performed yet.

7.2.3.2Running¶

We now run all simulations and analyse the resulting Monte Carlo data directly. On a typical laptop, this should take approximately 8 minutes, although the precise runtime depends on the available hardware.

for sim in sims:
    sim.run()
    sim.analysis()

7.2.4Results¶

During the analysis, pyALF calculates expectation values and their statistical uncertainties from the Monte Carlo bins. The average sign is reported for each observable, including the energy, but it is not dependent on the observable itself. For energy measurements, it is stored in Ener_scal_sign, with its estimated statistical error stored in Ener_scal_sign_err.

7.2.4.1Collecting the results¶

Before plotting the results, we collect the analysed observables from all simulations in a single table. Each row corresponds to one simulation and contains both its parameters and the measured observables.

results = pd.concat(
    [sim.get_obs() for sim in sims],
    ignore_index=True,
)

sign_data = (
    results[
        [
            "l1",
            "beta",
            "ham_chem",
            "mz",
            "Ener_scal_sign",
            "Ener_scal_sign_err",
        ]
    ]
    .rename(
        columns={
            "l1": "L",
            "ham_chem": "mu",
            "Ener_scal_sign": "average_sign",
            "Ener_scal_sign_err": "average_sign_error",
        }
    )
    .sort_values(["L", "beta", "mz", "mu"])
)

sign_data["decomposition"] = sign_data["mz"].map(
    {True: "$M_z$", False: "charge"}
)

sign_data

7.2.4.2Average sign as a function of chemical potential¶

We now plot the average sign as a function of chemical potential for one of the two HS decompositions. Colors distinguish the inverse temperatures, while different line styles distinguish the lattice sizes. The error bars show the statistical uncertainty returned by the pyALF analysis.

# Select the Hubbard-Stratonovich decomposition
selected_mz = False  # True: M_z, False: charge



import matplotlib.pyplot as plt
from matplotlib.lines import Line2D


plot_data = sign_data[
    sign_data["mz"] == selected_mz
].copy()

# Colors distinguish inverse temperatures
colors = {
    1.0: "tab:blue",
    2.0: "tab:red",
    4.0: "tab:green",
}

# Line styles distinguish lattice sizes
line_styles = {
    4: "-",
    6: "--",
}

fig, ax = plt.subplots(figsize=(6, 4))

for (L, beta), data in plot_data.groupby(["L", "beta"]):
    data = data.sort_values("mu")

    ax.errorbar(
        data["mu"],
        data["average_sign"],
        yerr=data["average_sign_error"],
        color=colors[beta],
        linestyle=line_styles[L],
        marker="o",
        linewidth=1.5,
        markersize=5,
        capsize=3,
    )

ax.set_xlabel(r"Chemical potential $\mu/t$")
ax.set_ylabel(r"Average sign $\langle s\rangle_{\mathrm{abs}}$")

if selected_mz:
    ax.set_title(r"$M_z$ decomposition")
else:
    ax.set_title("Charge decomposition")

ax.set_ylim(0, 1.05)
ax.grid(alpha=0.25)

# Handles showing the color assigned to each inverse temperature
beta_handles = [
    Line2D(
        [0],
        [0],
        color=colors[beta],
        linewidth=2,
        label=fr"$\beta t={beta:g}$",
    )
    for beta in sorted(plot_data["beta"].unique())
]

# Handles showing the line style assigned to each lattice size
size_handles = [
    Line2D(
        [0],
        [0],
        color="black",
        linestyle=line_styles[L],
        linewidth=2,
        label=fr"$L={L}$",
    )
    for L in sorted(plot_data["L"].unique())
]

# Single movable legend
legend = ax.legend(
    handles=beta_handles + size_handles,
    ncol=2,
    loc="center right",
    frameon=False,
)

legend.set_draggable(True)

hs_name = "Mz" if selected_mz else "charge"
filename = f"average_sign_{hs_name}.png"

fig.savefig(
    filename,
    dpi=200,
    bbox_inches="tight",
)

print(f"Figure saved as: {filename}")

Average sign of the Hubbard model — Average sign as a function of chemical potential for different inverse temperatures and lattice sizes.

At half filling, corresponding to $\mu=0$ , particle-hole symmetry protects the simulations from the sign problem. Away from half filling, cancellations between configurations with different signs become possible, and the average sign can decrease.

In general, the sign problem becomes more severe with increasing inverse temperature and system size. Its precise behavior also depends on the chosen HS decomposition, even though the charge and $M_z$ formulations describe the same physical spin- $1/2$ Hubbard model. This illustrates that the sign problem is not solely a property of the Hamiltonian, but also depends on its auxiliary-field representation.

Exercise. Repeat the plot for the other HS decomposition. How does the average sign change with chemical potential, temperature, and lattice size? Which decomposition produces the more severe sign problem for the parameters considered here?

7.2 Exploring the sign problem