diff --git a/final_plots.py b/final_plots.py
index bfbf17a18f3917ac42dd3ce6f1575b4fb40a4695..8d87f84205283a3a5184b7dd61099fc3b37edacc 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -28,6 +28,10 @@ from gen_pulse_data import fourier_coef_gauss_symmetric
 from plot_pulses import load_all_pulses_full, integrate_ft
 from kondo import solveTV0_scalar
 
+# color maps with better sampling
+viridis = plt.cm.viridis.resampled(0x10000)
+seismic = plt.cm.seismic.resampled(0x10000)
+
 # In this program all energies are given in units of the RTRG Kondo
 # temperature Tkrg, which is an integration constant of the E-flow RG
 # equations. The more conventional definition of the Kondo temperature is
@@ -36,6 +40,7 @@ from kondo import solveTV0_scalar
 #TK_VOLTAGE = 3.44249 # for rtol=1e-10, atol=1e-12, voltage_branches=10
 TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE = 3.44334 # with correction
+TK_VOLTAGE_O3 = 3.307473 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE_O3P = 3.4593 # with correction
 TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
@@ -207,74 +212,6 @@ def export_omega5_pgfplots(filename="figdata/omega5_interp.dat", dc_steps=3, ac_
     export_matrix_pgfplots(filename, vdc, vac, gdc, gac, idc, iac, phase, header="vdc vac gdc gac idc iac phase")
 
 
-def export_omega5():
-    "deprecated!"
-    omega = 16.5372
-    dm = DataManager()
-    # Full overview
-    data = dm.list(
-            omega = omega,
-            vdc = None,
-            vac = None,
-            method = "mu",
-            d = 1e9,
-            xL = 0.5,
-            solver_tol_rel = 1e-8,
-            solver_tol_abs = 1e-10,
-            voltage_branches = 4,
-            truncation_order = 3,
-            include_Ga = True,
-            solve_integral_exactly = False,
-            )
-    bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
-            | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
-            | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
-            | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
-            | DataManager.SOLVER_FLAGS["deleted"]
-    data = data.loc[data.solver_flags & bad_flags == 0]
-    vdc, vac, gdc, idc, iac, phase = filter_grid_data(
-            data,
-            omega = omega,
-            vac_min = 0,
-            vac_max = 165.372,
-            vac_num = 51,
-            vdc_min = 0,
-            vdc_max = 165.372,
-            vdc_num = 51,
-            )
-    with open("figdata/vdc_vac_omega16.5372.dat", "w") as file:
-        file.write("vac vdc gdc idc iac")
-        for i in range(51):
-            file.write("\n")
-            np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
-    ## Lines at constant Vac
-    #for i in range(11):
-    #    with open(f"figdata/vdc_vac{i}_omega16.5372.dat", "w") as file:
-    #        file.write("vac vdc gdc idc iac\n")
-    #        np.savetxt(file, np.array([vac[5*i]/omega, vdc[5*i]/omega, np.pi*gdc[5*i], idc[5*i], iac[5*i]]).T)
-    ## Lines at constant Vdc
-    #for i in range(11):
-    #    with open(f"figdata/vac_vdc{i}_omega16.5372.dat", "w") as file:
-    #        file.write("vac vdc gdc idc iac\n")
-    #        np.savetxt(file, np.array([vac[:,5*i]/omega, vdc[:,5*i]/omega, np.pi*gdc[:,5*i], idc[:,5*i], iac[:,5*i]]).T)
-    # Zoom to primary Kondo peak
-    vdc, vac, gdc, idc, iac, phase = filter_grid_data(
-            data,
-            omega = omega,
-            vac_min = 0,
-            vac_max = 16.5372,
-            vac_num = 21,
-            vdc_min = 0,
-            vdc_max = 16.5372,
-            vdc_num = 21,
-            )
-    with open("figdata/vdc_vac_omega16.5372_zoom.dat", "w") as file:
-        file.write("vac vdc gdc idc iac")
-        for i in range(21):
-            file.write("\n")
-            np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
-
-
 def export_interp(
         filename,
         fixed_parameter,
@@ -294,6 +231,9 @@ def export_interp(
         o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+        o3r_scale_x = TK_VOLTAGE_O3/TK_VOLTAGE,
+        o3r_scale_y = TK_VOLTAGE_O3/TK_VOLTAGE,
+        o3r_scale_fixed = TK_VOLTAGE_O3/TK_VOLTAGE,
         **parameters
         ):
     """
@@ -331,6 +271,7 @@ def export_interp(
     order:
         o2    2nd order RG equations
         o3    3rd order RG equations without Ga and with approximated integral
+        o3r   3rd order RG equations corrected Tk scale
         o3a   3rd order RG equations with approximated integral
         o3p   3rd order RG equations with full integral
 
@@ -365,12 +306,14 @@ def export_interp(
     required_flags = dict(
             o2 = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             o3 = 0,
+            o3r = 0,
             o3a = DataManager.SOLVER_FLAGS["include_Ga"],
             o3p = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             )
     bad_flags = dict(
             o2 = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+            o3r = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
             o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
             )
@@ -379,8 +322,13 @@ def export_interp(
     o2_data = o2_data.loc[o2_data.solver_flags & (global_bad_flags | bad_flags["o2"] | required_flags["o2"]) == required_flags["o2"]]
     o2_data = o2_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
 
+    o3r_data = dm.list(**{fixed_parameter:fixed_value*o3r_scale_fixed}, **parameters)
+    o3r_data = o3r_data.loc[o3r_data.solver_flags & (global_bad_flags | bad_flags["o3r"] | required_flags["o3r"]) == required_flags["o3r"]]
+    o3r_data = o3r_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
+    all_data = dict(o2=o2_data, o3r=o3r_data, o3a=data, o3p=data, o3=data)
+
     def selection(order, method, padding=False, **kwargs):
-        thedata = o2_data if order == "o2" else data
+        thedata = all_data[order]
         selected = (thedata.solver_flags & required_flags[order] == required_flags[order]) \
                 & (thedata.solver_flags & bad_flags[order] == 0) \
                 & (thedata.method == method)
@@ -425,6 +373,8 @@ def export_interp(
         xy_data = (x_func(reduced_data), y_func(reduced_data))
         if order == "o2":
             xy_data = (xy_data[0]/o2_scale_x, xy_data[1]/o2_scale_y)
+        elif order == "o3r":
+            xy_data = (xy_data[0]/o3r_scale_x, xy_data[1]/o3r_scale_y)
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
@@ -487,14 +437,17 @@ def export_interp(
             **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
             **gen_data("o2", "mu"),
             **gen_data("o3", "mu"),
+            **gen_data("o3r", "mu"),
             **gen_data("o3a", "mu"),
             **gen_data("o3p", "mu", integral_method=-1),
             **gen_data("o2", "J", False),
             **gen_data("o3", "J", False),
+            **gen_data("o3r", "J", False),
             **gen_data("o3a", "J", False),
             **gen_data("o3p", "J", False, integral_method=-1),
             **gen_data("o2", "J", True),
             **gen_data("o3", "J", True),
+            **gen_data("o3r", "J", True),
             **gen_data("o3a", "J", True),
             **gen_data("o3p", "J", True, integral_method=-1),
             )
@@ -516,9 +469,6 @@ def export_interp_relative_o3a(
         kyorder = 2,
         spline_s = 2e-4,
         special_selection = None,
-        o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
-        o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
-        o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         **parameters
         ):
     """
@@ -554,7 +504,6 @@ def export_interp_relative_o3a(
         gamma Γ(E=0) (not an observable)
 
     order:
-        o2    2nd order RG equations
         o3    3rd order RG equations without Ga and with approximated integral
         o3a   3rd order RG equations with approximated integral
         o3p   3rd order RG equations with full integral
@@ -876,9 +825,9 @@ def prepare_RGconvergence():
     dc_indices = np.array([0, 9, 21, 30, 90, 165])
     ac_indices = np.array([32, 52, 80, 124, 172])
 
-    real_cmap = lambda x: plt.cm.seismic((x+1)/2)[...,:3]
-    imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
-    imwrite(f"figdata/colorbar_seismic_vertical.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
+    real_cmap = lambda x: seismic((x+1)/2)[...,:3]
+    imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
+    imwrite(f"figdata/colorbar_seismic_vertical.png", np.array(0xffff*seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
     def export_img(name, array):
         imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
 
@@ -888,9 +837,9 @@ def prepare_RGconvergence():
 
     file = open("tikz/colorscales.tex", "w")
 
-    # O2
     omega = 16.5372
-    omega_o2 = 48.695054
+    #omega_o2 = 48.695054
+    #omega_o3r = 15.8884
     vdc = data["vdc"][0] / omega
     vac = data["vac"][:,0] / omega
 
@@ -900,7 +849,7 @@ def prepare_RGconvergence():
 \\newcommand\\extentymin{{{1.5*vac[0]-0.5*vac[1]}}}
 \\newcommand\\extentymax{{{1.5*vac[-1]-0.5*vac[-2]}}}""", file=file)
 
-    for method, u, ref in (("o2","mu","o3p"), ("o3","mu","o3p"), ("o3a","mu","o3p"), ("o3a_ispline","mu","o3a"), ("o3a","J","o3a"), ("o3a_p","J","o3a")):
+    for method, u, ref in (("o2","mu","o3p"), ("o3","mu","o3p"), ("o3a","mu","o3p"), ("o3a_ispline","mu","o3a"), ("o3a","J","o3a"), ("o3a_p","J","o3a"), ("o3r","mu","o3p")):
         gdiff = data[f"gdc_{u}_{method}"] - data[f"gdc_mu_{ref}"]
         if u != "mu":
             method += "_" + u
@@ -915,7 +864,7 @@ def prepare_RGconvergence():
         print(f"Color scale for {method} relative: {normalize}")
         print(f"\\newcommand\\scale{tex_name}rel{{{normalize}}}% convergence_{method}_relative.png", file=file)
 
-    for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1)):
+    for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1),("o3","mu",1),("o3","J",1),("o3_p","J",1)):
         gdiff = sign*data_rel[f"gdc_{u}_{method}"]
         if u != "mu":
             method += "_" + u
@@ -1752,7 +1701,7 @@ def export_convergence_pgfplots(simplified_initial_conditions = False):
 
     suffix = "_bad" if simplified_initial_conditions else "_good"
     def export_img(name, array):
-        imwrite(f"figdata/{name}{suffix}.png", np.array(0xffff*plt.cm.viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
+        imwrite(f"figdata/{name}{suffix}.png", np.array(0xffff*viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
     idc_max = np.nanmax(idc_fit_a)
     iac_max = np.nanmax(iac_fit_a)
     i_max = max(idc_max, 2*iac_max)
@@ -1805,7 +1754,7 @@ def export_floquet_matrices_pgfplots():
 
     scale = lambda x: np.log10(np.abs(x))/9 + 1
     def export_img(name, array):
-        imwrite(f"figdata/{name}.png", np.array(0xffff*plt.cm.viridis(scale(array))[...,:3], dtype=np.uint16), format="PNG-FI")
+        imwrite(f"figdata/{name}.png", np.array(0xffff*viridis(scale(array))[...,:3], dtype=np.uint16), format="PNG-FI")
 
     z_ref_j = z_matrices[("J", True, True)]
     z_ref_mu = z_matrices[("mu", True, False)]
@@ -1921,7 +1870,7 @@ def export_Ga_max_pgfplots():
     grid_g2_max = griddata((vdc, vac), g2_max, (vdc_mesh, vac_mesh), method="cubic")
 
     def export_img(name, array):
-        imwrite(f"figdata/{name}.png", np.array(0xffff*plt.cm.viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
+        imwrite(f"figdata/{name}.png", np.array(0xffff*viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
 
     print(f"Extent: xmin={extent[0]/omega:.6g}, xmax={extent[1]/omega:,.6g}, ymin={extent[2]/omega:.6g}, ymax={extent[3]/omega:,.6g}")
     norm = np.nanmax(np.abs(grid_ga_max))
diff --git a/plot.py b/plot.py
index 10a0447b91fd0118f0c5757e26583052232fc5c0..96a0c4520c22ad226eebfea953c3c1c602077608 100644
--- a/plot.py
+++ b/plot.py
@@ -44,6 +44,7 @@ TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE = 3.44334 # with correction
 TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE_O3P = 3.4593 # with correction
+#TK_VOLTAGE_O3 = 3.44428 # for rtol=1e-8, atol=1e-10, voltage_branches=4
 TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
 #TK_VOLTAGE_ORDER2 = 10.13754 # with correction
 
@@ -101,7 +102,7 @@ def main():
             help="Truncation order of RG equations")
     parser.add_argument("--integral_method", metavar="int", type=int, default=-15,
             help="method for solving frequency integral (-1 for exact solution)")
-    parser.add_argument("--include_Ga", metavar="bool", type=bool, default=False,
+    parser.add_argument("--include_Ga", metavar="bool", type=int, default=1,
             help="include Ga in RG equations")
     parser.add_argument("--good_flags", metavar="int", type=int, default=0,
             help = "Require these flags")
@@ -1366,8 +1367,10 @@ def convergence_overview(dm, **parameters):
     vdc_max = 82.686
     vac_max = 82.686
     exponent = 3
-    if parameters["good_flags"] & 8:
+    if parameters["good_flags"] & DataManager.SOLVER_FLAGS["simplified_initial_conditions"]:
         exponent = 2
+    elif parameters["good_flags"] & DataManager.SOLVER_FLAGS["improved_initial_conditions"]:
+        exponent = 4
     d_num = 5
     log10d_min = 5
     log10d_max = 9
@@ -1403,8 +1406,8 @@ def convergence_overview(dm, **parameters):
     ax3.plot(logd_inv3, -gdc_fit_b.reshape((1,-1))*logd_inv3.reshape((d_num,1)), color="gray", linewidth=0.5)
     ax3.plot(logd_inv3, gdc_fit_a.reshape((1,-1)) - gdc.reshape((d_num,-1)))
 
-    #fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True)
-    fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True)
+    fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True)
+    #fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True)
     extent = (-0.5*vdc_max/(vdc_num-1), vdc_max*(1+0.5/(vdc_num-1)), -0.5*vac_max/(vac_num-1), vac_max*(1+0.5/(vac_num-1)))
     ax1.set_title(r"$I_\mathrm{avg}(\Lambda_0=\infty)$")
     ax2.set_title(r"$\frac{|I_\mathrm{avg}(\Lambda_0=10^9 T_K') - I_\mathrm{avg}(\Lambda_0=\infty)|}{I_\mathrm{avg}(\Lambda_0=10^9T_K')}$")
@@ -1427,12 +1430,12 @@ def convergence_overview(dm, **parameters):
     fig.colorbar(img, ax=ax5)
     img = ax6.imshow(-iac_diff, extent=extent, origin="lower")
     fig.colorbar(img, ax=ax6)
-    #img = ax7.imshow(gdc_fit_a, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax7)
-    #img = ax8.imshow(gdc_fit_b*logd_inv3[4]/gdc_fit_a, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax8)
-    #img = ax9.imshow(gdc_diff, extent=extent, origin="lower")
-    #fig.colorbar(img, ax=ax9)
+    img = ax7.imshow(gdc_fit_a, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax7)
+    img = ax8.imshow(gdc_fit_b*logd_inv3[4]/gdc_fit_a, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax8)
+    img = ax9.imshow(gdc_diff, extent=extent, origin="lower")
+    fig.colorbar(img, ax=ax9)
 
     #cmap = plt.cm.seismic
     #gdc_diff = (merged.dc_conductance_1 - merged.dc_conductance_2)/merged.dc_conductance_2
@@ -1833,8 +1836,8 @@ def RGeq_comparison_contours2(
         **trashparams):
     omega = 16.5372
     kwargs = dict(linewidths=0.8)
-    colors = dict(o2="#ff0000", o3="#6000ff", o3a="#00f000", o3p="#000000", ref="#808080")
-    zorder = dict(o2=1, o3=2, o3a=3, o3p=4)
+    colors = dict(o2="#ff0000", o3="#6000ff", o3r="#ff00ff", o3a="#00f000", o3p="#000000", ref="#808080")
+    zorder = dict(o2=1, o3=2, o3r=5, o3a=3, o3p=4)
     approximations = list(zorder.keys())
     data = np.load(filename)
     vdc = data["vdc"] / omega
@@ -1842,8 +1845,11 @@ def RGeq_comparison_contours2(
     levels = (1/np.linspace(1, 12, 12)[::-1]) / prefactor
     fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True, sharey=True)
     for approx in approximations:
-        contour = ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
-        save_contour_coordinates(contour, f"figdata/contour_{observable}_mu_{approx}.dat", "vdc", "vac", "gdc")
+        try:
+            contour = ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+            save_contour_coordinates(contour, f"figdata/contour_{observable}_mu_{approx}.dat", "vdc", "vac", "gdc")
+        except KeyError:
+            settings.logger.warning(f"No data available for {observable}_mu_{approx}")
         try:
             contour = ax.contour(vdc, vac, data[f"{observable}_J_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
             save_contour_coordinates(contour, f"figdata/contour_{observable}_J_{approx}.dat", "vdc", "vac", "gdc")
diff --git a/tikz/contour_wrapper.tex b/tikz/contour_wrapper.tex
deleted file mode 100644
index fd51447292cb401e7b74d4b6c8425c94a76b268f..0000000000000000000000000000000000000000
--- a/tikz/contour_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{contour}
-\end{document}
diff --git a/tikz/convergence-o2-rel.tex b/tikz/convergence-o2-rel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..0cee967eb042322037b14af578ba0d9495bd83cd
--- /dev/null
+++ b/tikz/convergence-o2-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o2}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o2_relative.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {leading order, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o2_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleolrel,
+      xmax = \scaleolrel,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.05, 0, 0.05},
+      xticklabels = {$-5\%$, $0$, $5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleolrel, xmax=\scaleolrel] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3-rel.tex b/tikz/convergence-o3-rel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c8f5fdcae9682784fc82f743f1a45ed18fa13cca
--- /dev/null
+++ b/tikz/convergence-o3-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o3}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o3_relative.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {without $G^a$, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleoxrel,
+      xmax = \scaleoxrel,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.05, 0, 0.05},
+      xticklabels = {$-5\%$, $0$, $5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxrel, xmax=\scaleoxrel] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3a-rel.tex b/tikz/convergence-o3a-rel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..f2217d4f502c4b54223e9e95cdd2c70d87111527
--- /dev/null
+++ b/tikz/convergence-o3a-rel.tex
@@ -0,0 +1,158 @@
+\tikzsetnextfilename{convergence_o3a}%
+\begin{tikzpicture}
+  \pgfplotsset{set layers = axis on top}
+  \begin{axis}[
+      name = img,
+      anchor = north west,
+      width = 4cm,
+      height = 4cm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0, ymax = 10,
+      xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\vac~(\Omega)$},
+      every axis plot post/.append style = {
+        line width=0.2pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \addplot graphics[xmin=\extentxmin, xmax=\extentxmax, ymin=\extentymin, ymax=\extentymax] {../figdata/convergence_o3p_relative_fromdiff.png};
+    \node at (2,9.2) {(b)};
+    % vertical lines
+    \addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
+    \addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
+    \addplot[no marks, black] coordinates {(0,4.0) (10,4.0)};
+    \addplot[no marks, black] coordinates {(0,6.2) (10,6.2)};
+    \addplot[no marks, black] coordinates {(0,8.6) (10,8.6)};
+    \node[anchor=west] at (7.0,1.6) {\small △};
+    \node[anchor=west] at (7.5,2.6) {\small ▽};
+    \node[anchor=west] at (8.0,4.0) {\small □};
+    \node[anchor=west] at (8.5,6.2) {\small ○};
+    \node[anchor=west] at (9.0,8.6) {\small ◊};
+    % horizontal lines
+    \addplot[no marks, black] coordinates {(0,0) (0,10)};
+    \addplot[no marks, black] coordinates {(0.3,0) (0.3,10)};
+    \addplot[no marks, black] coordinates {(0.7,0) (0.7,10)};
+    \addplot[no marks, black] coordinates {(1,0) (1,10)};
+    \addplot[no marks, black] coordinates {(3,0) (3,10)};
+    \addplot[no marks, black] coordinates {(5.5,0) (5.5,10)};
+    \node[anchor=south] at (0,8.5) {\small ▶};
+    \node[anchor=south] at (0.3,7) {\small ◀};
+    \node[anchor=south] at (0.7,5) {\small ■};
+    \node[anchor=south] at (1,4.0) {\large •};
+    \node[anchor=south] at (3,0.2) {\small ♦};
+    \node[anchor=south] at (5.5,0.2) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      at = (img.north),
+      name = ac,
+      anchor = south,
+      width = 4cm,
+      height = 3cm,
+      yshift = 2.5mm,
+      scale only axis,
+      xmin = 0, xmax = 10,
+      ymin = 0.08, ymax = 0.88,
+      %xlabel = {$\vdc~(\Omega)$},
+      ylabel = {$\gdc~(2e^2/h)$},
+      %xticklabel pos = top,
+      xticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+      legend entries = {approx.~integral, full RG equations},
+      legend cell align = left,
+      legend style = {
+        legend pos = outer north east,
+        draw = none,
+        fill = none,
+      },
+    ]
+    \node at (2,0.79) {(a)};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3a_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc4.00}+0.24}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc6.20}+0.12}] {../figdata/convergence_mu_o3p_ac.dat};
+    \addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc8.60}}] {../figdata/convergence_mu_o3p_ac.dat};
+    \node[anchor=west] at (7.0,0.61) {\small △};
+    \node[anchor=west] at (7.5,0.48) {\small ▽};
+    \node[anchor=west] at (8.0,0.365) {\small □};
+    \node[anchor=west] at (8.5,0.251) {\small ○};
+    \node[anchor=west] at (9.0,0.145) {\small ◊};
+  \end{axis}
+  \begin{axis}[
+      name = dc,
+      at = (img.east),
+      anchor = west,
+      width = 3cm,
+      height = 4cm,
+      xshift = 2.5mm,
+      scale only axis,
+      xmin = 0.05, xmax = 1.62,
+      ymin = 0, ymax = 10,
+      xlabel = {$\gdc~(2e^2/h)$},
+      %ylabel = {$\vac~(\Omega)$},
+      %yticklabel pos = right,
+      yticklabel = \empty,
+      every axis plot post/.append style = {
+        line width=0.5pt,
+        line cap=round,
+        line join=round
+      },
+      clip = false,
+    ]
+    \node at (1.3,9.2) {(c)};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, red] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc1.00}+0.24}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc3.00}+0.12}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \addplot[no marks, densely dotted] table [x expr={\thisrow{gdc5.50}}, y=vac] {../figdata/convergence_mu_o3p_dc.dat};
+    \node[anchor=south] at (0.72,10) {\small ▶};
+    \node[anchor=south] at (0.59,10) {\small ◀};
+    \node[anchor=south] at (0.47,10) {\small ■};
+    \node[anchor=south] at (0.35,10) {\large •};
+    \node[anchor=south] at (0.23,10) {\small ♦};
+    \node[anchor=south] at (0.11,10) {\small ★};
+  \end{axis}
+  \begin{axis}[
+      name = colorbar,
+      at = (dc.north west),
+      xmin = -\scaleoxprelX,
+      xmax = \scaleoxprelX,
+      xshift = 5mm,
+      yshift = 11mm,
+      anchor = south west,
+      enlargelimits = false,
+      axis on top,
+      height = .75em,
+      width = 2.5cm,
+      scale only axis,
+      tickpos = bottom,
+      tick align = outside,
+      ytick = \empty,
+      xtick = {-0.005, 0, 0.005},
+      xticklabels = {$-0.5\%$, $0$, $0.5\%$},
+      scaled ticks = false,
+      xlabel = {relative deviation},
+      xlabel shift = -14.6mm,
+    ]
+    \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxprelX, xmax=\scaleoxprelX] {../figdata/colorbar_seismic.png};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o3a_ispline.tex b/tikz/convergence-o3a_ispline.tex
index 53542c15facd74545fce77478a8811fbf7c1e9b9..93e34fc4b044a39547b26ec96e2d0d05af705e61 100644
--- a/tikz/convergence-o3a_ispline.tex
+++ b/tikz/convergence-o3a_ispline.tex
@@ -153,7 +153,7 @@
       xticklabels = {$-5\%$, $0$, $5\%$},
       scaled ticks = false,
       xlabel = {relative deviation},
-      xlabel shift = -14mm,
+      xlabel shift = -14.6mm,
     ]
     \addplot graphics[ymin=0, ymax=1, xmin=-\scaleoxaxisplinerel, xmax=\scaleoxaxisplinerel] {../figdata/colorbar_seismic.png};
   \end{axis}
diff --git a/tikz/g_overview_vdc_vac_wrapper.tex b/tikz/g_overview_vdc_vac_wrapper.tex
deleted file mode 100644
index c8580d48db6f2ced7f4f5300955956d5e7ec7f7a..0000000000000000000000000000000000000000
--- a/tikz/g_overview_vdc_vac_wrapper.tex
+++ /dev/null
@@ -1,13 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{dc}}
-\newcommand\vac{V_\mathrm{ac}}
-\newcommand\gdc{G_\mathrm{dc}}
-\newcommand\gac{G_\mathrm{ac}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\begin{document}
-\input{g_overview_vdc_vac}
-\end{document}
diff --git a/tikz/omega5-3d_wrapper.tex b/tikz/omega5-3d_wrapper.tex
deleted file mode 100644
index 6eeb7c4a107c93d22a63e19c32a1144be7076f34..0000000000000000000000000000000000000000
--- a/tikz/omega5-3d_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{omega5-3d}
-\end{document}
diff --git a/tikz/overview_wrapper.tex b/tikz/overview_wrapper.tex
deleted file mode 100644
index 4efc8fc54b9621be2cd63bf7375fbd85fccb1fc2..0000000000000000000000000000000000000000
--- a/tikz/overview_wrapper.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{standalone}
-\newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\vdc{V_\mathrm{avg}}
-\newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
-\newcommand\gac{G_\mathrm{osc}}
-\usepackage{pgfplots}
-\pgfplotsset{compat=1.18}
-\pgfplotsset{surf shading/precision=pdf}
-\usepackage{fontspec}
-\defaultfontfeatures{Ligatures=TeX}
-\setsansfont[
-  UprightFont = *-Regular,
-  BoldFont = *-Bold,
-  ItalicFont = *-Italic,
-  %BoldItalicFont = LinBiolinum_K, % TODO: Check this font, use a better one
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSans}
-\setmainfont[
-  %UprightFont = *Display-Regular,
-  UprightFont = *-Regular,
-  BoldFont = *-Semibold,
-  ItalicFont = *-Italic,
-  BoldItalicFont = *-SemiboldItalic,
-  SmallCapsFont = *-Regular,
-  SmallCapsFeatures = {Letters=SmallCaps},
-  %SlantedFont = *-Italic,
-  SlantedFont = *-Regular,
-  SlantedFeatures = {FakeSlant=0.15},
-]{LibertinusSerif}
-%\setmonofont[
-%  UprightFont = *-Regular
-%]{LibertinusMono}
-\setmonofont[
-  UprightFont = *-Regular,
-  BoldFont = *-Medium,
-  ItalicFont = *-It,
-  BoldItalicFont = *-MediumIt,
-]{SourceCodePro}
-
-\usepackage{mathtools}
-%\usepackage{amsmath}
-\usepackage{amssymb}
-%\usepackage{amsfonts}
-\usepackage[math-style=ISO, bold-style=TeX, partial=upright, nabla=upright]{unicode-math}
-\setmathfont{Libertinus Math}
-\setmathfont{STIXTwoMath-Regular.otf}[range=cal]
-\setmathfont{STIXTwoMath-Regular.otf}[range={\mathscr},StylisticSet=1]
-\setmathfont{STIXTwoMath-Regular.otf}[range={"02022, "02218, "02200, "0223C, "02248}] % Slightly larger \bullet and \circ; larger \forall; different \sim, \approx
-%\setmathfont{Latin Modern Math}[range={"027F8-"027FA}] % \iff, \implies, \impliedby
-%\setmathfont{STIX Two Math}[range={"027E8, "027E9}] % \langle, \rangle
-\setmathfont{STIXTwoMath-Regular.otf}[range={"0221A-"0221C, "023B7}] % root (sqrt, cbrt, \dots) symbols
-\begin{document}
-\input{overview}
-\end{document}