Fix statements that contradict the packages, and clear defects

Facts, checked against sparse-ir 2.1.4 (pylibsparseir 0.9.0) and
SparseIR.jl 2.1.5:

- The default τ sampling points lie in (0, β) (use_positive_taus=True),
  so reversing the array maps τ to β − τ. The migration guide and the
  FLEX_py comment said the v2 default was [−β/2, β/2], and the
  migration guide recommended -gtau*gtau[::-1], which is wrong under
  the default. With use_positive_taus=False the points are folded into
  (−β/2, β/2], and for an odd number of points β/2 is a point but −β/2
  is not, so reversing the array is not τ → −τ in general.
- FLEX_py and FLEX_jl ran with use_positive_taus=False, and their
  temperature sweeps have an odd basis size (L = 43 and 39), so
  -grit*grit[::-1] paired the wrong times. On main, λ_d becomes −10.5
  at T = 0.03 in Python, and the Julia sweep aborts there with
  "U*max(chi0) >= 1". Use the default instead: the sweep then gives
  λ_d from 0.53 to 0.95, and the T = 0.1 run is unchanged. The False
  branch stays as a demonstration but asserts symmetric points. Also
  fix the FLEX_jl comment G(-tau,-r) = G(beta-tau,r) (sign) and
  "beta_min" (β_max = 1/T_min sets Λ = 400).
- DLR poles are the roots of V_L, not the extrema of V'_{L-1}
  (DLR_py, DLR_jl).
- Condition numbers: second_order_perturbation measures 163 (τ) and 466
  (Matsubara), not 180 and 250. sparse_sampling_demo's Matsubara value
  is 213, about four times the τ value, not "slightly larger".
- GW plots against matsb.wn / matsf.wn, which hold the reduced
  frequencies n, so label the x axes n instead of ω.

Defects:

- second_order_perturbation: Σ(τ) = U²G²(τ)G(β − τ) on the default
  points, i.e. srt = U*U*grt*grt*grt[::-1,:]. The old minus sign gave an
  acausal self-energy: Im Σ(iπ/β) = +5.72e-4 at k = Γ and positive at
  every k. It is now −5.72e-4 and negative at every k.
- transformation_py plotted gtau1 instead of gtau2, and set the x label
  twice instead of the y label. transformation_jl hard-coded a
  developer's libsparseir path and loaded Revise.
- DLR_py / DLR_jl: the semicircle √(2/π)√(1−ω²) integrates to 1.25; use
  (2/π)√(1−ω²).
- FLEX_jl docstring F = |G|²Δ (the code uses −). The FLEX_py χ₀ title
  was not a raw string, so "\n" became a newline. FLEX_py c25 passed
  maxiter/U_maxiter and ignored the it_max/U_it_max it defines. The
  FLEX_jl χ₀ plot was titled χ.
- Slips: the Matsubara fit in sparse_sampling_demo was commented "Fit
  G(τ) on the sampling times". The Eliashberg Δ(iν^F) panel had a
  bosonic axis label, and C/T was labeled C(T). SpM dropped the minus
  in G(τ_i) = −Σ U S ρ and in (G_l + S_l ρ_l)². analytic_continuation
  mixed the indices m and i. TPSC wrote χ⁰ for χ₀. BdG had
  lim_{τ→0+} for η → 0⁺ and ê(j) for e⃗(j). second_order_perturbation
  said "three integers" for a 2D point. transformation_py said V_l is
  dense "around τ = 0, β": the roots of V_l cluster around ω = 0.
  giv_tensor is the Matsubara data (was gtau_tensor).
- DMFT_IPT_py linked to second_order_perturbation_py.md, which does
  not exist; link the notebook.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

Author: shinaoka

src/BgG_jl.ipynb

@@ -77,11 +77,11 @@
     "The regarded and advanced Green's functions are, respectively, defined as \n",
     "\n",
     "$$\n",
-    "\\boldsymbol{G}^{\\rm R}(\\omega) = \\lim_{\\tau \\rightarrow 0+} \\boldsymbol{G}(\\omega + \\mathrm{i} \\eta)\n",
+    "\\boldsymbol{G}^{\\rm R}(\\omega) = \\lim_{\\eta \\rightarrow 0^+} \\boldsymbol{G}(\\omega + \\mathrm{i} \\eta)\n",
     "$$\n",
     "\n",
     "$$\n",
-    "\\boldsymbol{G}^{\\rm A}(\\omega) = \\lim_{\\tau \\rightarrow 0+} \\boldsymbol{G}(\\omega - \\mathrm{i} \\eta)\n",
+    "\\boldsymbol{G}^{\\rm A}(\\omega) = \\lim_{\\eta \\rightarrow 0^+} \\boldsymbol{G}(\\omega - \\mathrm{i} \\eta)\n",
     "$$\n",
     "\n",
     "In order to obtain physical quantities (e.g., density of states) from Green's functions, we introduce the following useful $2N$-component unit-vectors $\\vec{e}(i)$ and $\\vec{h}(i)$ $(1 \\le i \\le N)$, which are, respectively, defined as\n",
@@ -140,7 +140,7 @@
     "the superconducting mean-field is expressed as \n",
     "\n",
     "$$\n",
-    "\\langle c_i c_j \\rangle = T \\sum_{\\nu} \\hat{e}(j)^T \\vec{x}(i,\\nu).\n",
+    "\\langle c_i c_j \\rangle = T \\sum_{\\nu} \\vec{e}(j)^T \\vec{x}(i,\\nu).\n",
     "$$\n",
     "\n",
     "For example, s-wave superconducting order parameter is defined as \n",
@@ -445,7 +445,7 @@
     "If one wants to use the Matsubara frequencies, one has to introduce the cutoff of the frequencies. Then the mean-fields are expressed as \n",
     "\n",
     "$$\n",
-    "\\langle c_i c_j \\rangle = T \\sum_{|\\nu| < \\nu_c} \\hat{e}(j)^T \\vec{x}(i,\\nu).\n",
+    "\\langle c_i c_j \\rangle = T \\sum_{|\\nu| < \\nu_c} \\vec{e}(j)^T \\vec{x}(i,\\nu).\n",
     "$$\n",
     "\n",
     "Here, $\\nu_c$ is a cutoff value (`ωc` in the code). \n",

src/DLR_jl.ipynb

@@ -16,7 +16,7 @@
     "    \\rho(\\omega) = \\sum_{p=1}^L c_p \\delta(\\omega - \\bar{\\omega}_p),\n",
     "$$\n",
     "\n",
-    "where sampling frequencies $\\{\\bar{\\omega}_1, \\cdots, \\bar{\\omega}_{L}\\}$ are chosen to the extrema of $V'_{L-1}(\\omega)$ [2].\n",
+    "where the poles $\\{\\bar{\\omega}_1, \\cdots, \\bar{\\omega}_{L}\\}$ are the roots of $V_L(\\omega)$, the first basis function beyond the $L$ functions kept in the basis [2].\n",
     "In the original paper [1], a more and different systematic way of choosing poles was used.\n",
     "We can transform data between $\\rho_l$ and $c_p$ through the relation\n",
     "\n",
@@ -45,7 +45,7 @@
     "We consider fermions and the semi circular DOS\n",
     "\n",
     "$$\n",
-    "\\rho(\\omega) = \\sqrt{\\frac{2}{\\pi}} \\sqrt{1-\\omega^2}.\n",
+    "\\rho(\\omega) = \\frac{2}{\\pi} \\sqrt{1-\\omega^2}.\n",
     "$$\n",
     "\n",
     "The corresponding Green's function is given by\n",
@@ -116,7 +116,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rho(omega) = sqrt(1-omega^2)/sqrt(0.5*π)\n",
+    "rho(omega) = (2/π) * sqrt(1-omega^2)\n",
     "\n",
     "omega = LinRange(-wmax, wmax, 1000)\n",
     "plot(omega, rho.(omega), xlabel=latexstring(\"\\\\omega\"), ylabel=latexstring(\"\\\\rho(\\\\omega)\"))"

src/DLR_py.ipynb

@@ -10,14 +10,14 @@
     "## Theory\n",
     "\n",
     "We explain how to transform expansion coefficients in IR to the discrete Lehmann representation (DLR) [1].\n",
-    "In `sparse-ir`, we choose poles based on the exrema of $V_l(\\omega)$ [2] (In the original paper [1], a more systematic way is proposed).\n",
+    "In `sparse-ir`, we choose the poles based on the roots of the IR basis function $V_L(\\omega)$ [2] (In the original paper [1], a more systematic way is proposed).\n",
     "We model the spectral function as\n",
     "\n",
     "$$\n",
     "    \\rho(\\omega) = \\sum_{p=1}^L c_p \\delta(\\omega - \\bar{\\omega}_p),\n",
     "$$\n",
     "\n",
-    "where sampling frequencies $\\{\\bar{\\omega}_1, \\cdots, \\bar{\\omega}_{L}\\}$ are chosen to the extrema of $V'_{L-1}(\\omega)$.\n",
+    "where the poles $\\{\\bar{\\omega}_1, \\cdots, \\bar{\\omega}_{L}\\}$ are the roots of $V_L(\\omega)$, the first basis function beyond the $L$ functions kept in the basis.\n",
     "This choice is heuristic but allows us a numerically stable transform between $\\rho_l$ and $c_p$ through the relation\n",
     "\n",
     "$$\n",
@@ -45,7 +45,7 @@
     "We consider fermions and the semi circular DOS\n",
     "\n",
     "$$\n",
-    "\\rho(\\omega) = \\sqrt{\\frac{2}{\\pi}} \\sqrt{1-\\omega^2}.\n",
+    "\\rho(\\omega) = \\frac{2}{\\pi} \\sqrt{1-\\omega^2}.\n",
     "$$\n",
     "\n",
     "The corresponding Green's function is given by\n",
@@ -119,7 +119,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rho = lambda omega: np.sqrt(1-omega**2)/np.sqrt(0.5*np.pi)\n",
+    "rho = lambda omega: (2/np.pi) * np.sqrt(1-omega**2)\n",
     "\n",
     "omega = np.linspace(-wmax, wmax, 1000)\n",
     "plt.xlabel(r'$\\omega$')\n",

src/DMFT_IPT_py.ipynb

@@ -44,7 +44,7 @@
     "\n",
     "\n",
     "### IPT Solver\n",
-    "The IPT is an inexpensive way of solving the impurity problem (instead of applying e.g. Quantum Monte Carlo (QMC) or Exact Diagonalization (ED) methods). The self-energy is approximated by the second-order perturbation theory (SOPT/GF2, see other [tutorial](second_order_perturbation_py.md))\n",
+    "The IPT is an inexpensive way of solving the impurity problem (instead of applying e.g. Quantum Monte Carlo (QMC) or Exact Diagonalization (ED) methods). The self-energy is approximated by the second-order perturbation theory (SOPT/GF2, see other [tutorial](second_order_perturbation_py.ipynb))\n",
     "\n",
     "$$\n",
     "    \\Sigma(\\tau) = U^2\\mathcal{G}^3(\\tau).\n",

src/FLEX_jl.ipynb

@@ -401,16 +401,20 @@
     "function ckio_calc(solver::FLEXSolver)\n",
     "    \"\"\" Calculate irreducible susciptibility chi0(iv,q) \"\"\"\n",
     "    # chi0(tau, r) = -G(tau, r) * G(-tau, -r)\n",
-    "    # For positive taus [0, beta]: G(-tau, -r) = G(beta-tau, r)\n",
+    "    # For positive taus in (0, beta): G(-tau, -r) = -G(beta-tau, -r) = -G(beta-tau, r) (inversion symmetry),\n",
+    "    #   and G(beta-tau, r) is obtained by reversing the tau array\n",
     "    #   => chi0(tau, r) = G(tau, r) * G(beta-tau, r)\n",
-    "    # For symmetric taus [-beta/2, beta/2]: G(-tau, r) is obtained by reversing the tau array\n",
+    "    # For taus in (-beta/2, beta/2]: G(-tau, r) is obtained by reversing the tau array\n",
+    "    #   only for an even number of points (for an odd number, beta/2 is a point but -beta/2 is not)\n",
     "    #   => chi0(tau, r) = -G(tau, r) * G(-tau, r)\n",
     "    grit_reversed = reverse(solver.grit, dims=1)\n",
     "    if solver.mesh.use_positive_taus\n",
     "        # G(tau) * G(beta - tau)\n",
     "        crit = solver.grit .* grit_reversed\n",
     "    else\n",
     "        # - G(tau) * G(- tau)\n",
+    "        tau = solver.mesh.IR_basis_set.smpl_tau_f.sampling_points\n",
+    "        @assert reverse(tau) ≈ -tau\n",
     "        crit = -solver.grit .* grit_reversed\n",
     "    end\n",
     "\n",
@@ -497,7 +501,7 @@
    "outputs": [],
    "source": [
     "# initialize calculation\n",
-    "use_positive_taus = false\n",
+    "use_positive_taus = true\n",
     "\n",
     "IR_basis_set = FiniteTempBasisSet(beta, Float64(wmax), IR_tol; use_positive_taus=use_positive_taus)\n",
     "mesh = Mesh(nk1, nk2, IR_basis_set; use_positive_taus=use_positive_taus)\n",
@@ -541,7 +545,7 @@
    "source": [
     "# plot 2D k-dependence of lowest Matsubara frequency of e.g. chi0\n",
     "heatmap(myx, myy, real.(solver.ckio[solver.mesh.iw0_b, :, :]),\n",
-    "    title=latexstring(\"\\\\mathrm{Re}\\\\,\\\\chi(k,i\\\\nu_0)\"), xlabel=latexstring(\"k_x/\\\\pi\"), ylabel=latexstring(\"k_y/\\\\pi\"),\n",
+    "    title=latexstring(\"\\\\mathrm{Re}\\\\,\\\\chi_0(k,\\\\mathrm{i}\\\\nu=0)\"), xlabel=latexstring(\"k_x/\\\\pi\"), ylabel=latexstring(\"k_y/\\\\pi\"),\n",
     "    c=:viridis,\n",
     "    xlim=(0, 2), ylim=(0, 2), aspect_ratio=1.0, size=(370, 300))"
    ]
@@ -681,7 +685,7 @@
     "end\n",
     "\n",
     "function frit_calc(gap_solver::LinearizedGapSolver)\n",
-    "    \"\"\" Calculate (linearized) anomalous Green function F = |G|^2 * delta for evaluating the gap equation \"\"\"\n",
+    "    \"\"\" Calculate (linearized) anomalous Green function F = -|G|^2 * delta for evaluating the gap equation \"\"\"\n",
     "    fkio = -gap_solver.gkio .* conj(gap_solver.gkio) .* gap_solver.delta\n",
     "\n",
     "    # Fourier transform\n",
@@ -793,7 +797,7 @@
     "# numerical parameters\n",
     "nk1, nk2 = 64, 64    # k-mesh sufficiently dense!\n",
     "nk = nk1 * nk2\n",
-    "IR_Lambda = 10^3      # dimensionless IR parameter >= w_max * beta_min = 400\n",
+    "IR_Lambda = 10^3      # dimensionless IR parameter >= wmax * beta_max = 400\n",
     "IR_tol = 1e-7      # accuary for l-cutoff of IR basis functions\n",
     "sfc_tol = 1e-4      # accuracy for self-consistent iteration\n",
     "it_max = 30        # maximal number of iterations in self-consistent cycle\n",
@@ -802,7 +806,7 @@
     "\n",
     "# initialize first IR basis set (no recalculation afterwrds)\n",
     "beta_init = 1.0 / T_values[1]\n",
-    "use_positive_taus = false\n",
+    "use_positive_taus = true\n",
     "IR_basis_set = FiniteTempBasisSet(beta_init, IR_Lambda / beta_init, IR_tol; use_positive_taus=use_positive_taus)\n",
     "\n",
     "# set initial self_energy - will be set to previous calculation step afterwards\n",

src/FLEX_py.ipynb

@@ -94,11 +94,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# If true, we use sampling points in [0, beta].\n",
-    "# If false, we use sampling points in [-beta/2, beta/2].\n",
-    "# sparse-ir version 1 defaults to use_positive_taus=True, \n",
-    "# but sparse-ir version 2 defaults to use_positive_taus=False.\n",
-    "use_positive_taus = False"
+    "# If true, we use sampling points in (0, beta); reversing the array maps tau to beta - tau.\n",
+    "# If false, we use sampling points in (-beta/2, beta/2]; reversing the array maps tau to -tau\n",
+    "# only if the number of points is even (see ckio_calc below).\n",
+    "# sparse-ir defaults to use_positive_taus=True, which reproduces the sampling points of version 1.\n",
+    "use_positive_taus = True"
    ]
   },
   {
@@ -355,6 +355,10 @@
     "            ckio = self.grit * self.grit[::-1, :]\n",
     "        else:\n",
     "            # - G(tau) * G(- tau)\n",
+    "            # Reversing maps tau to -tau only for an even number of points: for an odd number,\n",
+    "            # the unfolded set (-beta/2, beta/2] contains beta/2 but not -beta/2.\n",
+    "            tau = self.mesh.IR_basis_set.smpl_tau_f.sampling_points\n",
+    "            assert np.allclose(tau[::-1], -tau)\n",
     "            ckio = - self.grit * self.grit[::-1, :]\n",
     "\n",
     "        # Fourier transform\n",
@@ -482,7 +486,7 @@
     "ax.set_ylabel('$k_y/\\pi$')\n",
     "ax.set_ylim([0,2])\n",
     "ax.set_aspect('equal')\n",
-    "ax.set_title('Re $\\chi_0(k,i\\nu_0)$')\n",
+    "ax.set_title(r'Re $\\chi_0(k,\\mathrm{i}\\nu^\\mathrm{B}=0)$')\n",
     "plt.colorbar()\n",
     "plt.show()"
    ]
@@ -745,7 +749,7 @@
     "# numerical parameters\n",
     "nk1, nk2  = 64, 64    # k-mesh sufficiently dense!\n",
     "nk        = nk1*nk2\n",
-    "IR_Lambda = 10**3     # dimensionless IR parameter >= w_max * beta_min = 400\n",
+    "IR_Lambda = 10**3     # dimensionless IR parameter >= wmax * beta_max = 400\n",
     "IR_tol    = 1e-8      # desired accuary for l-cutoff of IR basis functions\n",
     "sfc_tol   = 1e-4      # desired accuracy for self-consistent iteration\n",
     "it_max    = 30        # maximal number of iterations in self-consistent cycle\n",
@@ -775,12 +779,12 @@
     "    \n",
     "    # calculate FLEX loop\n",
     "    solver = FLEXSolver(mesh, U, n, sigma_init=sigma_init, sfc_tol=sfc_tol, \n",
-    "                        maxiter=maxiter, U_maxiter=U_maxiter, mix=mix, verbose=False)\n",
+    "                        maxiter=it_max, U_maxiter=U_it_max, mix=mix, verbose=False)\n",
     "    solver.solve()\n",
     "    sigma_init = solver.sigma\n",
     "\n",
     "    # calculate linearized gap equation\n",
-    "    gap_solver = LinearizedGapSolver(solver, maxiter=maxiter, sfc_tol=sfc_tol, verbose=False)\n",
+    "    gap_solver = LinearizedGapSolver(solver, maxiter=it_max, sfc_tol=sfc_tol, verbose=False)\n",
     "    gap_solver.solve()\n",
     "    \n",
     "    # save data for plotting\n",

src/GW_py.ipynb

@@ -273,7 +273,7 @@
     "pl.plot(matsb.wn,P_iw_b.imag,label='imaginary')\n",
     "pl.plot(matsb.wn,P_iw_b.real,label='real')\n",
     "pl.title(r'$P(\\mathrm{i}\\nu^\\mathrm{B})$')\n",
-    "pl.xlabel(r'$\\omega$')\n",
+    "pl.xlabel(r'$n$')\n",
     "pl.legend()"
    ]
   },
@@ -321,7 +321,7 @@
     "pl.plot(matsb.wn,W_iw_b.imag,label='imaginary')\n",
     "pl.plot(matsb.wn,W_iw_b.real,label='real')\n",
     "pl.title(r'$W(\\mathrm{i}\\nu^\\mathrm{B})$')\n",
-    "pl.xlabel(r'$\\omega$')\n",
+    "pl.xlabel(r'$n$')\n",
     "pl.legend()"
    ]
   },
@@ -495,10 +495,10 @@
     "fig, ax = pl.subplots(1,2, figsize=(12,5) )\n",
     "ax[0].plot(matsf.wn,E_iw_f_U.imag)\n",
     "ax[0].set_title(r'Imaginary part of $\\Sigma(\\mathrm{i}\\nu^\\mathrm{F})$')\n",
-    "ax[0].set_xlabel(r'$\\omega$')\n",
+    "ax[0].set_xlabel(r'$n$')\n",
     "ax[1].plot(matsf.wn,E_iw_f_U.real)\n",
     "ax[1].set_title(r'Real part of $\\Sigma(\\mathrm{i}\\nu^\\mathrm{F})$')\n",
-    "ax[1].set_xlabel(r'$\\omega$')"
+    "ax[1].set_xlabel(r'$n$')"
    ]
   },
   {
@@ -536,7 +536,7 @@
     "pl.plot(matsf.wn,G_iw_f.imag,label='imaginary')\n",
     "pl.plot(matsf.wn,G_iw_f.real,label='real')\n",
     "pl.title(r'$G(\\mathrm{i}\\nu^\\mathrm{F})$')\n",
-    "pl.xlabel(r'$\\omega$')\n",
+    "pl.xlabel(r'$n$')\n",
     "pl.legend()"
    ]
   },
@@ -550,10 +550,12 @@
     "    pl.plot(matsf.wn,E_iw_f_arr[i].imag,label='current')\n",
     "    pl.title(r'$\\Sigma(\\mathrm{i}\\nu^\\mathrm{F})$ in the current and previous iteration')\n",
     "    pl.plot(matsf.wn,E_iw_f_arr[i-1].imag,'.',label='previous')\n",
+    "    pl.xlabel(r'$n$')\n",
     "    pl.legend()\n",
     "    pl.show()\n",
     "    pl.plot(matsf.wn,np.abs(E_iw_f_arr[i]-E_iw_f_arr[i-1]))\n",
     "    pl.title(r'$|\\Sigma_i(\\mathrm{i}\\nu^\\mathrm{F})-\\Sigma_{i-1}(\\mathrm{i}\\nu^\\mathrm{F})|$' )\n",
+    "    pl.xlabel(r'$n$')\n",
     "    "
    ]
   },

src/TPSC_py.ipynb

@@ -86,7 +86,7 @@
     "When implementing the TPSC, a few points need to be treated carefully which we adress in the following:\n",
     "\n",
     "* The constant Hartree term $V_{\\mathrm{H}} = U$ in the interaction $V$ and respective self-energy term $\\Sigma_H = U n_{\\sigma} = U\\frac{n}{2}$ can be absorbed into the definition of the chemical potential $\\mu$. Otherwise we would have to treat this term separately, since the IR basis cannot model the $\\delta$-peak in frequency space well/compactly. In this case, evaluate the term analytically and add it after the Fourier transformation step.\n",
-    "* An upper bound for the renormalized spin vertex $U_{\\mathrm{sp}}$ exists. Since the denominator spin susceptibility $\\chi_{\\mathrm{sp}}$ should not diverge, the upper bound is given by the RPA critical interaction value $U_{\\mathrm{crit}} = 1/\\mathrm{max}\\{\\chi^0\\}$. Mathematically, the function $f(U_{\\mathrm{sp}}) = 2\\sum \\chi_{\\mathrm{sp}}(U_{\\mathrm{sp}}) - n + \\frac{U_{\\mathrm{sp}}}{2U}n^2$, from which $U_{\\mathrm{sp}}$ is determined, turns unstable for $U_{\\mathrm{sp}} \\geq U_{\\mathrm{crit}}$ (try plotting $f(U_{\\mathrm{sp}})$!). At this point, TPSC is not applicable and, e.g., the temperature $T$ is too low or the (unrenormalized) interaction $U$ too large.\n",
+    "* An upper bound for the renormalized spin vertex $U_{\\mathrm{sp}}$ exists. Since the denominator spin susceptibility $\\chi_{\\mathrm{sp}}$ should not diverge, the upper bound is given by the RPA critical interaction value $U_{\\mathrm{crit}} = 1/\\mathrm{max}\\{\\chi_0\\}$. Mathematically, the function $f(U_{\\mathrm{sp}}) = 2\\sum \\chi_{\\mathrm{sp}}(U_{\\mathrm{sp}}) - n + \\frac{U_{\\mathrm{sp}}}{2U}n^2$, from which $U_{\\mathrm{sp}}$ is determined, turns unstable for $U_{\\mathrm{sp}} \\geq U_{\\mathrm{crit}}$ (try plotting $f(U_{\\mathrm{sp}})$!). At this point, TPSC is not applicable and, e.g., the temperature $T$ is too low or the (unrenormalized) interaction $U$ too large.\n",
     "* An internal accuracy check $\\frac{1}{2}\\mathrm{Tr}(\\Sigma G) = U \\langle n_{\\uparrow} n_{\\downarrow}\\rangle$ can be employed to test the validity of TPSC (not done here)."
    ]
   },

src/analytic_continuation_py.ipynb

@@ -345,12 +345,12 @@
     "grid $\\{\\omega_1,\\ldots,\\omega_M\\}$ of frequencies and a function $f(x)$ and expand:\n",
     "\n",
     "$$\n",
-    "    \\rho(\\omega) = \\sum_{m=1}^M a_i f(\\omega - \\omega_i),\n",
+    "    \\rho(\\omega) = \\sum_{m=1}^M a_m f(\\omega - \\omega_m),\n",
     "$$\n",
     "\n",
-    "where $a_i$ are now the expansion coefficients. (More advanced methods also optimize over $\\omega_m$ and/or add some shape parameter of $f$ to the optimization parameters.)\n",
+    "where $a_m$ are now the expansion coefficients. (More advanced methods also optimize over $\\omega_m$ and/or add some shape parameter of $f$ to the optimization parameters.)\n",
     "\n",
-    "It is useful to use some probability distribution as $f(x)$, as this allows one to translate non-negativity and norm of $\\rho(\\omega)$ to non-negativity and norm of $a_i$. Since one can only observe \"broadened\" spectra in experiment for any given temperature, a natural choice is the Lorentz (Cauchy) distribution:\n",
+    "It is useful to use some probability distribution as $f(x)$, as this allows one to translate non-negativity and norm of $\\rho(\\omega)$ to non-negativity and norm of $a_m$. Since one can only observe \"broadened\" spectra in experiment for any given temperature, a natural choice is the Lorentz (Cauchy) distribution:\n",
     "\n",
     "$$\n",
     "    f(\\omega) = \\frac1\\pi \\frac\\eta{\\omega^2 + \\eta^2},\n",

src/eliashberg_holstein_py.ipynb

@@ -598,8 +598,8 @@
     "    plt.show()\n",
     "\n",
     "    plt.plot(iv_f.imag, res['delta_iv'].real, marker=\"x\")\n",
-    "    plt.xlabel(r\"$\\nu_m$\")\n",
-    "    plt.ylabel(r\"Re $\\Delta(\\mathrm{i}\\nu_m)$\")\n",
+    "    plt.xlabel(r\"$\\nu^\\mathrm{F}$\")\n",
+    "    plt.ylabel(r\"Re $\\Delta(\\mathrm{i}\\nu^\\mathrm{F})$\")\n",
     "    plt.xlim([-0.8, 0.8])\n",
     "    plt.xticks([-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8])\n",
     "    plt.ylim([-0.05, 0.3])\n",
@@ -745,7 +745,7 @@
     "plt.plot(temps, specific_heat/temps, marker=\"o\")\n",
     "plt.ylim([0, 1200])\n",
     "plt.xlabel(r\"$T$\")\n",
-    "plt.ylabel(r\"$C(T)$\")\n",
+    "plt.ylabel(r\"$C(T)/T$\")\n",
     "plt.show()"
    ]
   },