diff --git a/final_plots.py b/final_plots.py
index 8cf83724c587700a0e9dffabfa45e5e5dcdb5c9e..837063880223a25af7c014fd30f9f6a88f0a5619 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -615,8 +615,10 @@ def prepare_RGconvergence():
     vdc = data["vdc"][0] / omega
     vac = data["vac"][:,0] / omega
 
-    for method, ref in (("o2","o3p"), ("o3","o3p"), ("o3a","o3p"), ("o3a_ispline","o3a")):
-        gdiff = data[f"gdc_mu_{method}"] - data[f"gdc_mu_{ref}"]
+    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")):
+        gdiff = data[f"gdc_{u}_{method}"] - data[f"gdc_mu_{ref}"]
+        if u != "mu":
+            method += "_" + u
         normalize = np.nanmax(np.abs(gdiff))
         export_img(f"convergence_{method}_diff", gdiff/normalize)
         print(f"Color scale for {method} difference: {np.pi*normalize}")
@@ -625,14 +627,33 @@ def prepare_RGconvergence():
         export_img(f"convergence_{method}_relative", grel/normalize)
         print(f"Color scale for {method} relative: {normalize}")
 
-    for method  in ("o2", "o3", "o3a", "o3p", "o3a_ispline"):
-        g = np.pi*data[f"gdc_mu_{method}"]
-        with open(f"figdata/convergence_mu_{method}_dc.dat", "w") as file:
+    gdiff = data[f"gdc_mu_o3a"] - data_vb7[f"gdc_mu_o3a"]
+    normalize = np.nanmax(np.abs(gdiff))
+    export_img(f"convergence_o3a_vb7_diff", gdiff/normalize)
+    print(f"Color scale for o3a_vb7 difference: {np.pi*normalize}")
+    grel = gdiff/data_vb7[f"gdc_mu_o3a"]
+    normalize = np.nanmax(np.abs(grel))
+    export_img(f"convergence_o3a_vb7_relative", grel/normalize)
+    print(f"Color scale for o3a_vb7 relative: {normalize}")
+
+    for method, u  in (("o2","mu"), ("o3","mu"), ("o3a","mu"), ("o3p","mu"), ("o3a_ispline","mu")):
+        if u != "mu":
+            method += "_" + u
+        g = np.pi*data[f"gdc_{u}_{method}"]
+        with open(f"figdata/convergence_{u}_{method}_dc.dat", "w") as file:
             file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
             np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
-        with open(f"figdata/convergence_mu_{method}_ac.dat", "w") as file:
+        with open(f"figdata/convergence_{u}_{method}_ac.dat", "w") as file:
             file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
             np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
+    g = np.pi*data_vb7[f"gdc_mu_o3a"]
+    with open(f"figdata/convergence_mu_o3a_vb7_dc.dat", "w") as file:
+        file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
+        np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
+    with open(f"figdata/convergence_mu_o3a_vb7_ac.dat", "w") as file:
+        file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
+        np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
+
     print("DC voltage: " + " ".join("%.6g"%vdc[i] for i in dc_indices))
     print("AC voltage: " + " ".join("%.6g"%vac[i] for i in ac_indices))
 
@@ -1419,6 +1440,77 @@ def export_convergence_pgfplots(simplified_initial_conditions = False):
     print(f"I diff range: 0, {i_diff_max:.6g}")
 
 
+def export_floquet_matrices_pgfplots():
+    dm = DataManager()
+    omega = 16.5371763
+    vdc = 19.84461156
+    vac = 29.76691734
+    vb = 4
+    params = dict(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, voltage_branches=vb, omega=omega, vdc=vdc, vac=vac)
+    kondos = {}
+    g200_matrices = {}
+    g201_matrices = {}
+    z_matrices = {}
+    gamma_matrices = {}
+    deltaGammaL_matrices = {}
+    nmax_ref = 19
+    for kondo in dm.load_all(**params):
+        nmax = kondo.nmax
+        idx = (kondo.method, kondo.nmax > 20, kondo.padding>0)
+        if nmax > nmax_ref:
+            z_matrices[idx] = kondo.z[vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
+            gamma_matrices[idx] = kondo.gamma[vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
+            deltaGammaL_matrices[idx] = kondo.deltaGammaL[nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
+            g200_matrices[idx] = kondo.g2[0,0,vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
+            g201_matrices[idx] = kondo.g2[0,1,vb, nmax-nmax_ref:nmax_ref-nmax, nmax-nmax_ref:nmax_ref-nmax]
+        elif nmax == nmax_ref:
+            z_matrices[idx] = kondo.z[vb]
+            gamma_matrices[idx] = kondo.gamma[vb]
+            deltaGammaL_matrices[idx] = kondo.deltaGammaL
+            g200_matrices[idx] = kondo.g2[0,0,vb]
+            g201_matrices[idx] = kondo.g2[0,1,vb]
+        else:
+            raise ValueError(f"nmax={nmax} should not be smaller than nmax_ref={nmax_ref}")
+        kondos[idx] = kondo
+
+    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")
+
+    z_ref_j = z_matrices[("J", True, True)]
+    z_ref_mu = z_matrices[("mu", True, False)]
+
+    g201_ref_j = g201_matrices[("J", True, True)]
+    g201_ref_mu = g201_matrices[("mu", True, False)]
+
+    export_img("floquet_matrix_z_mu_small_p", z_matrices[("mu",False,True)])
+    export_img("floquet_matrix_z_mu_small", z_matrices[("mu",False,False)])
+    export_img("floquet_matrix_z_mu_large_p", z_matrices[("mu",True,True)])
+    export_img("floquet_matrix_z_mu_large", z_matrices[("mu",True,False)])
+    export_img("floquet_matrix_z_J_small_p", z_matrices[("J",False,True)])
+    export_img("floquet_matrix_z_J_small", z_matrices[("J",False,False)])
+    export_img("floquet_matrix_z_J_large_p", z_matrices[("J",True,True)])
+    export_img("floquet_matrix_z_J_large", z_matrices[("J",True,False)])
+    export_img("floquet_matrix_z_mu_diff_p", z_matrices[("mu",False,True)] - z_ref_mu)
+    export_img("floquet_matrix_z_mu_diff", z_matrices[("mu",False,False)] - z_ref_mu)
+    export_img("floquet_matrix_z_J_diff_p", z_matrices[("J",False,True)] - z_ref_j)
+    export_img("floquet_matrix_z_J_diff", z_matrices[("J",False,False)] - z_ref_j)
+
+    export_img("floquet_matrix_g201_mu_small_p", g201_matrices[("mu",False,True)])
+    export_img("floquet_matrix_g201_mu_small", g201_matrices[("mu",False,False)])
+    export_img("floquet_matrix_g201_mu_large_p", g201_matrices[("mu",True,True)])
+    export_img("floquet_matrix_g201_mu_large", g201_matrices[("mu",True,False)])
+    export_img("floquet_matrix_g201_J_small_p", g201_matrices[("J",False,True)])
+    export_img("floquet_matrix_g201_J_small", g201_matrices[("J",False,False)])
+    export_img("floquet_matrix_g201_J_large_p", g201_matrices[("J",True,True)])
+    export_img("floquet_matrix_g201_J_large", g201_matrices[("J",True,False)])
+    export_img("floquet_matrix_g201_mu_diff_p", g201_matrices[("mu",False,True)] - g201_ref_mu)
+    export_img("floquet_matrix_g201_mu_diff", g201_matrices[("mu",False,False)] - g201_ref_mu)
+    export_img("floquet_matrix_g201_J_diff_p", g201_matrices[("J",False,True)] - g201_ref_j)
+    export_img("floquet_matrix_g201_J_diff", g201_matrices[("J",False,False)] - g201_ref_j)
+
+
+
 
 if __name__ == "__main__":
     from sys import argv
diff --git a/tikz/current_convergence.tex b/tikz/current_convergence.tex
index 46dfa5be58b9d522f24c5c032c64b64e390f94ba..cbbaa54c5d7085d9dbefcb39f98e38fa7d13fe16 100644
--- a/tikz/current_convergence.tex
+++ b/tikz/current_convergence.tex
@@ -89,7 +89,7 @@
     \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1.54775] {../figdata/colorbar_viridis.png};
   \end{axis}
   \begin{axis}[
-      name = colorbar,
+      name = colorbar_diff,
       at = (iac_diff.east),
       anchor = west,
       ymin = 0,
diff --git a/tikz/floquet_matrices.tex b/tikz/floquet_matrices.tex
new file mode 100644
index 0000000000000000000000000000000000000000..f05d5a10db3c078519e559916d6be04a28df091a
--- /dev/null
+++ b/tikz/floquet_matrices.tex
@@ -0,0 +1,313 @@
+\newcommand\matrixsize{24mm}
+\begin{tikzpicture}
+  \begin{axis}[
+      name = z_mu_small,
+      title = {no $U$, no extrap.},
+      title style = {yshift = -2mm},
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      ylabel = {absolute},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_mu_small.png};
+  \end{axis}
+  \begin{axis}[
+      name = z_mu_diff,
+      at = (z_mu_small.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      ylabel = {deviation},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_mu_diff.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = z_mu_small_p,
+      title = {no $U$, extrap.},
+      title style = {yshift = -2mm},
+      at = (z_mu_small.east),
+      xshift = 3mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_mu_small_p.png};
+  \end{axis}
+  \begin{axis}[
+      name = z_mu_diff_p,
+      at = (z_mu_small_p.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_mu_diff_p.png};
+  \end{axis}
+  %
+  \node[xshift=2mm, yshift=8mm] at (z_mu_small_p.north east) {$Z$};
+  %
+  \begin{axis}[
+      name = z_J_small,
+      title = {with $U$, no extrap.},
+      title style = {yshift = -2mm},
+      at = (z_mu_small_p.east),
+      xshift = 8mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_J_small.png};
+  \end{axis}
+  \begin{axis}[
+      name = z_J_diff,
+      at = (z_J_small.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_J_diff.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = z_J_small_p,
+      title = {with $U$, extrap.},
+      title style = {yshift = -2mm},
+      at = (z_J_small.east),
+      xshift = 3mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_J_small_p.png};
+  \end{axis}
+  \begin{axis}[
+      name = z_J_diff_p,
+      at = (z_J_small_p.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_J_diff_p.png};
+  \end{axis}
+  %%%%
+  %%%%
+  %%%%
+  %%%%
+  \begin{axis}[
+      name = g201_mu_small,
+      title = {no $U$, no extrap.},
+      title style = {yshift = -2mm},
+      at = (z_mu_diff.south),
+      yshift = -20mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      ylabel = {absolute},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_mu_small.png};
+  \end{axis}
+  \begin{axis}[
+      name = g201_mu_diff,
+      at = (g201_mu_small.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      ylabel = {deviation},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_mu_diff.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = g201_mu_small_p,
+      title = {no $U$, extrap.},
+      title style = {yshift = -2mm},
+      at = (g201_mu_small.east),
+      xshift = 3mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_mu_small_p.png};
+  \end{axis}
+  \begin{axis}[
+      name = g201_mu_diff_p,
+      at = (g201_mu_small_p.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_mu_diff_p.png};
+  \end{axis}
+  %
+  \node[xshift=2mm, yshift=8mm] at (g201_mu_small_p.north east) {$G^2_{LR}$};
+  %
+  \begin{axis}[
+      name = g201_J_small,
+      title = {with $U$, no extrap.},
+      title style = {yshift = -2mm},
+      at = (g201_mu_small_p.east),
+      xshift = 8mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_J_small.png};
+  \end{axis}
+  \begin{axis}[
+      name = g201_J_diff,
+      at = (g201_J_small.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      ytick = {-5, 0, 5},
+      yticklabels = {$5$, $0$, $-5$},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_J_diff.png};
+  \end{axis}
+  %
+  \begin{axis}[
+      name = g201_J_small_p,
+      title = {with $U$, extrap.},
+      title style = {yshift = -2mm},
+      at = (g201_J_small.east),
+      xshift = 3mm,
+      anchor = west,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      xticklabels = {},
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_J_small_p.png};
+  \end{axis}
+  \begin{axis}[
+      name = g201_J_diff_p,
+      at = (g201_J_small_p.south),
+      yshift = -3mm,
+      anchor = north,
+      width = \matrixsize,
+      height = \matrixsize,
+      xmin = -9.5, xmax = 9.5,
+      ymin = -9.5, ymax = 9.5,
+      scale only axis,
+      axis on top,
+      yticklabels = {},
+    ]
+    \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_J_diff_p.png};
+  \end{axis}
+  %%%%
+  %%%%
+  %%%%
+  %%%%
+  \begin{semilogyaxis}[
+      name = colorbar,
+      at = (z_J_diff_p.south east),
+      anchor = west,
+      ymin = 1e-9,
+      ymax = 1,
+      xshift = 5mm,
+      yshift = -10mm,
+      enlargelimits = false,
+      axis on top,
+      width = .75em,
+      height = 5cm,
+      scale only axis,
+      tickpos = right,
+      tick align = outside,
+      xtick = \empty,
+    ]
+    \addplot graphics[xmin=0, xmax=1, ymin=1e-9, ymax=1] {../figdata/colorbar_viridis.png};
+  \end{semilogyaxis}
+\end{tikzpicture}