@@ -73,6 +73,7 @@
     "#### Notes on practical implementation\n",
     "* We will try to find the Mott transition which is a phase transition of first order. To illustrate this, we will do calculations starting from the metallic and the insulating phase at the end of this tutorial. This will show a hysteresis, signaling the metastable/coexistence region of both solutions. For converged results, this will need many iteration steps.\n",
     "* In the insulating phase, the low-frequency part of the self-energy turns downwards and diverges. Using the above formula for calculating the renormalization factor $Z$ from a linearized self-energy, it would be negative. To circumvent this unphysical notion, $Z$ then is defined as zero in the implementation.\n",
+    "* At half filling with a symmetric density of states the exact self-energy is particle-hole symmetric: $\\Sigma(\\mathrm{i}\\nu)$ is purely imaginary and odd in $\\nu$. The Fourier transform via the IR basis preserves this only up to rounding errors, and the symmetric solution is an *unstable* fixed point of the DMFT iteration with respect to perturbations that break the symmetry: a deviation of $10^{-16}$ grows by roughly 40% per iteration until, after a few hundred iterations, the loop ends up in an unphysical, symmetry-broken state. We therefore project each new self-energy onto its symmetric part, $\\Sigma(\\mathrm{i}\\nu) \\to \\mathrm{i}\\,\\frac{1}{2}[\\mathrm{Im}\\,\\Sigma(\\mathrm{i}\\nu) - \\mathrm{Im}\\,\\Sigma(-\\mathrm{i}\\nu)]$. Without this projection, the converged calculations below show a spurious \"transition\" near $U/t_\\star \\approx 3.5$.\n",
     "* We include a mixing $p<1$ in each iteration step, such that the Green function of step $n+1$ is partially constructed from the old and new Green function as $G^{n+1} = p\\,G^{n+1} + (1-p)\\,G^{n}$. This smoothes too strong oscillations of the solution which is important especially near the phase transition point."
    ]
   },
@@ -149,6 +150,14 @@
     "        # frequency mesh (for plotting Green function/self-energy)\n",
     "        self.iwn_f = 1j * self.IR_basis_set.wn_f * np.pi * T\n",
     "\n",
+    "        # index of -n for every fermionic sampling frequency n (they come in ± pairs)\n",
+    "        wn_f = list(self.IR_basis_set.wn_f)\n",
+    "        self.mirror_f = np.array([wn_f.index(-n) for n in wn_f])\n",
+    "\n",
+    "    def ph_symmetrize(self, obj_wn):\n",
+    "        \"\"\" Particle-hole symmetric part of a fermionic function at half filling: i*Im, odd in iwn \"\"\"\n",
+    "        return 0.5j * (obj_wn.imag - obj_wn.imag[self.mirror_f])\n",
+    "\n",
     "    def smpl_obj(self, statistics):\n",
     "        \"\"\" Return sampling object for given statistic \"\"\"\n",
     "        smpl_tau = {'F': self.IR_basis_set.smpl_tau_f, 'B': self.IR_basis_set.smpl_tau_b}[statistics]\n",
@@ -251,6 +260,8 @@
     "\n",
     "            # Solve impurity problem in IPT\n",
     "            sigma = IPTSolver(self.mesh, self.g_weiss, self.U)\n",
+    "            # Enforce particle-hole symmetry (see the notes above)\n",
+    "            sigma = self.mesh.ph_symmetrize(sigma)\n",
     "            self.sigma = sigma*self.mix + self.sigma*(1-self.mix)\n",
     "\n",
     "            # Set new Green function and Weiss field\n",
@@ -478,7 +489,7 @@
     "    Z = CalcRenormalizationZ(solver)\n",
     "    Z_arr[it] = Z\n",
     "\n",
-    "    if it in [30, 34, 35, 36, 37, 40]:\n",
+    "    if it in [50, 54, 56, 57, 58, 60]:\n",
     "        plt.plot(np.imag(mesh.iwn_f[mesh.iw0_f:]), np.imag(solver.sigma[mesh.iw0_f:]),'-o', label='$U$={:.2f}'.format(U))\n",
     "\n",
     "\n",
@@ -505,7 +516,7 @@
    "id": "43279056-8223-432e-801f-364730395f07",
    "metadata": {},
    "source": [
-    "You can quite well see a jump of $Z$ between $U/t_\\star=3$ and 4 which is accompanied by the low-frequency part of the self-energy diverging instead of going to zero.\n",
+    "You can quite well see a jump of $Z$ between $U/t_\\star=5.7$ and 5.8 which is accompanied by the low-frequency part of the self-energy diverging instead of going to zero.\n",
     "\n",
     "In a last step, we want to map out the metastable coexistence region. In order to do this, we have to start every new $U$ calculation from the previous converged solution. If we start from small $U$, we start from the metallic region, whereas starting from large $U$ corresponds to starting from the insulating solution."
    ]
@@ -558,7 +569,9 @@
    "id": "e61f48f4-f538-4c19-b181-3db4490a7c17",
    "metadata": {},
    "source": [
-    "For comparison, we also plotted the transition line where we always start from the non-interacting Green function $G_0$. Since we do only a few iteration steps, we do not have completetly converged solutions, so that the $G_0$ line has a slightly larger critical $U$ value compared to the line which starts from the metallic solution. Try increasinging the parameter `maxiter` (to 1000 or so) and decreasing the parameter `sfc_tol` (to 1e-14 or so) and look at the results! (The lines should become equal.)\n",
+    "The metallic solution exists up to $U/t_\\star = 5.7$ ($U_{c2}$ lies between 5.7 and 5.8), while the insulating solution survives down to $U/t_\\star = 5.5$ ($U_{c1}$ lies between 5.4 and 5.5). In between, both solutions are stable and which one is obtained depends on the starting point: this is the coexistence region of the first-order Mott transition. For comparison, we also plotted the line where we always start from the non-interacting Green function $G_0$. It coincides with the line starting from the metallic solution, since $G_0$ lies in the basin of the metal wherever the metal exists; a scan from $G_0$ alone therefore only finds $U_{c2}$.\n",
+    "\n",
+    "Try also removing the line `sigma = self.mesh.ph_symmetrize(sigma)` in `DMFT_loop.solve`: rounding errors then drive the loop away from the particle-hole symmetric solution, and the scan shows a spurious \"transition\" near $U/t_\\star \\approx 3.5$.\n",
     "\n",
     "Often the energy unit is not the hopping $t_\\star$ but the half-bandwidth $D$, so we can also plot the results in this convention:"
    ]