diff --git a/data_management.py b/data_management.py
index 06a1f3b98a6013e0f24d89ed3906c94da2bc8802..e28ad52822cb4b64c510667aafa262886e66750d 100644
--- a/data_management.py
+++ b/data_management.py
@@ -72,7 +72,6 @@ def reduce_h5f(filename):
                     h5f.create_array(f"/data/{name}", var, arr)
 
 
-
 def random_string(length : int):
     """
     Generate random strings of alphanumerical characters with given length.
diff --git a/final_plots.py b/final_plots.py
index 800b708a14b807345bd57386978642f38ba3d98c..4cf2016e340c33dda3515f325cada11adddda73c 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -6,6 +6,7 @@
 Kondo FRTRG, generate high-quality plots for publication
 """
 
+import os
 import scipy.constants as sc
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -1044,6 +1045,106 @@ def export_pulse_charge_pgfplots(omega=1.5):
                comments = "")
 
 
+def export_harmonic_modes(
+        #omega = 3.4425351,
+        omega = 68.850702,
+        #omega = 16.5372,
+        n_phase_crossings = 10,
+        ):
+    label = "-omega20"
+    vac_omega = np.array([5, 10, 20, 40, 80])
+    dm = DataManager()
+    data = dm.list(
+                omega = omega,
+                d = 1e9,
+                vdc = 0,
+                method = "J",
+                voltage_branches = 0,
+                include_Ga = True,
+                integral_method = -15,
+                solver_tol_rel = 1e-8,
+                solver_tol_abs = 1e-10,
+                bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+                        | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
+                        | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
+                        | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
+                        | DataManager.SOLVER_FLAGS["deleted"],
+                )
+    current = []
+    t = np.linspace(-0.5, 0.5, 201)
+    current = []
+    phase = []
+    valid_vac_omega = []
+    phase_crossing_times = []
+    phase_crossing_currents = []
+    for vac in vac_omega:
+        sel = np.abs(data.vac-vac*omega) < 1e-6
+        if sel.sum() != 1:
+            continue
+        row = data[sel].iloc[-1]
+        filename = os.path.join(row.dirname, row.basename)
+        kondo, = KondoImport.read_from_h5(filename, row.hash)
+        nmax = row.nmax
+        coef = kondo.gammaL[nmax::-1,nmax]
+        iarr = 2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0]
+        iarr /= TK_VOLTAGE
+        current.append(iarr)
+        pc_times = np.arcsin(2*np.pi/vac*np.arange(n_phase_crossings)-1)/(2*np.pi)
+        pc_currents = 2*sum((c*np.exp(2j*np.pi*pc_times*n)).real for n,c in enumerate(coef)) - coef.real[0]
+        pc_currents /= TK_VOLTAGE
+        phase_crossing_times.append(pc_times)
+        phase_crossing_currents.append(pc_currents)
+        phase.append(vac/(2*np.pi)*(1+np.sin(2*np.pi*t)))
+        valid_vac_omega.append(vac)
+        coef /= coef[1]
+        np.savetxt(
+                f"figdata/harmonic_modes_vac{vac:.3g}{label}.dat",
+                np.array([np.arange(1, coef.size, 2), np.abs(coef[1::2])]).T,
+                header = "n i",
+                comments = "",
+                fmt = "%.6g")
+    np.savetxt(
+            f"figdata/harmonic_current{label}.dat",
+            np.array([t, *current, *phase]).T,
+            header = "t " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"p{vac:.3g}" for vac in valid_vac_omega),
+            comments = "",
+            fmt = "%.6g")
+    np.savetxt(
+            f"figdata/phase_crossings{label}.dat",
+            np.array([np.arange(n_phase_crossings), *phase_crossing_times, *phase_crossing_currents]).T,
+            header = "phase " + " ".join(f"t{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega),
+            comments = "",
+            fmt = "%.6g")
+
+    # adiabatic result
+    n = np.arange(1, 101, 2)
+    adiabatic_modes = []
+
+    data = dm.list(vac=0, omega=0, include_Ga=True, d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, nmax=0, voltage_branches=4, xL=0.5, truncation_order=3, integral_method=-15)
+    data = data.sort_values("vdc")
+    vdc = np.concatenate((-data.vdc[:0:-1], data.vdc))
+    idc = np.concatenate((data.dc_current[:0:-1], data.dc_current))
+    interp = interp1d(vdc, idc, kind="cubic")
+    for i, vac in enumerate(omega*vac_omega):
+        ni, nierr = quad_vec((lambda t: np.sin(n*t) * interp(vac*np.sin(t))), 0, np.pi, limit=200)
+        adiabatic_modes.append(np.abs(ni/ni[0]))
+    np.savetxt(
+            f"figdata/harmonic_modes_adiabatic{label}.dat",
+            np.array([n, *adiabatic_modes]).T,
+            header = "n " + " ".join(f"mode{vac:.0f}" for vac in omega*vac_omega),
+            comments = "",
+            fmt = "%.6g")
+
+    fig, (ax1, ax2) = plt.subplots(nrows=2)
+    for i, (c, p, pt, pc) in enumerate(zip(current, phase, phase_crossing_times, phase_crossing_currents)):
+        color = f"C{i+1}"
+        ax1.plot(t, c, color=color)
+        ax1.plot(pt, pc, "x", color=color)
+        ax2.plot(p, c, color=color)
+        ax2.plot(np.arange(n_phase_crossings), pc, "x", color=color)
+    plt.show()
+
+
 
 if __name__ == "__main__":
     from sys import argv
diff --git a/tikz/harmonic_current.tex b/tikz/harmonic_current.tex
new file mode 100644
index 0000000000000000000000000000000000000000..9f8c552c8dc0c7fa52ef7b1ee815ec581756cbea
--- /dev/null
+++ b/tikz/harmonic_current.tex
@@ -0,0 +1,73 @@
+\definecolor{green}{rgb}{0,.67,0}%
+%\tikzsetnextfilename{harmonic_current}%
+\begin{tikzpicture}
+  \begin{axis}[
+      name = current,
+      anchor = north,
+      width = 5cm,
+      height = 4cm,
+      ylabel = {$I(t)~(e/T)$},
+      xlabel = {$t~(T)$},
+      xmin = -0.5, xmax = 0.5,
+      ymin = -3, ymax = 3,
+      scaled ticks = false,
+      scale only axis,
+    ]
+    \addplot[magenta, no marks] table [x=t, y=i5] {../figdata/harmonic_current-omega1.dat};
+    \addplot[blue, no marks] table [x=t, y=i10] {../figdata/harmonic_current-omega1.dat};
+    \addplot[green, no marks] table [x=t, y=i20] {../figdata/harmonic_current-omega1.dat};
+    \addplot[orange, no marks] table [x=t, y=i40] {../figdata/harmonic_current-omega1.dat};
+    \addplot[red, no marks] table [x=t, y=i80] {../figdata/harmonic_current-omega1.dat};
+    %
+    \addplot[magenta, only marks, mark=|] table [x=t5, y=i5] {../figdata/phase_crossings-omega1.dat};
+    \addplot[blue, only marks, mark=|] table [x=t10, y=i10] {../figdata/phase_crossings-omega1.dat};
+    \addplot[green, only marks, mark=|] table [x=t20, y=i20] {../figdata/phase_crossings-omega1.dat};
+    \addplot[orange, only marks, mark=|] table [x=t40, y=i40] {../figdata/phase_crossings-omega1.dat};
+    \addplot[red, only marks, mark=|] table [x=t80, y=i80] {../figdata/phase_crossings-omega1.dat};
+  \end{axis}
+  \begin{axis}[
+      name = phase,
+      at = (current.east),
+      anchor = west,
+      width = 2.5cm,
+      height = 4cm,
+      %xlabel = {$\varphi$},
+      xlabel = {$\int_{t^*}^tds\, V(s)$},
+      xmin = -0.2, xmax = 7,
+      ymin = -3, ymax = 3,
+      yticklabels = {},
+      xtick = {0, 2, 4, 6},
+      xticklabels = {$0$, $4\pi$, $8\pi$, $12\pi$},
+      minor xtick = {1, 3, 5},
+      scaled ticks = false,
+      scale only axis,
+    ]
+    \addplot[magenta] table [x=p5, y=i5] {../figdata/harmonic_current-omega1.dat};
+    \addplot[blue] table [x=p10, y=i10] {../figdata/harmonic_current-omega1.dat};
+    \addplot[green] table [x=p20, y=i20] {../figdata/harmonic_current-omega1.dat};
+    \addplot[orange] table [x=p40, y=i40] {../figdata/harmonic_current-omega1.dat};
+    \addplot[red] table [x=p80, y=i80] {../figdata/harmonic_current-omega1.dat};
+    %
+    \addplot[magenta, only marks, mark=|] table [x=phase, y=i5] {../figdata/phase_crossings-omega1.dat};
+    \addplot[blue, only marks, mark=|] table [x=phase, y=i10] {../figdata/phase_crossings-omega1.dat};
+    \addplot[green, only marks, mark=|] table [x=phase, y=i20] {../figdata/phase_crossings-omega1.dat};
+    \addplot[orange, only marks, mark=|] table [x=phase, y=i40] {../figdata/phase_crossings-omega1.dat};
+    \addplot[red, only marks, mark=|] table [x=phase, y=i80] {../figdata/phase_crossings-omega1.dat};
+  \end{axis}
+  \begin{axis}[
+      at = (phase.east),
+      xshift = 2cm,
+      anchor = west,
+      width = 4cm,
+      height = 3.2cm,
+      xlabel = {$t~(T)$},
+      ylabel = {$I(t)~(e/T)$},
+      xmin = -0.5, xmax = 0.5,
+      ymin = -0.6, ymax = 0.6,
+      scaled ticks = false,
+      scale only axis,
+    ]
+    \addplot[blue] table [x=t, y=i16] {../figdata/rlm_current.dat};
+    \node[black] at (axis description cs:0.5,0.2) {RLM};
+  \end{axis}
+\end{tikzpicture}%
diff --git a/tikz/harmonic_modes.tex b/tikz/harmonic_modes.tex
new file mode 100644
index 0000000000000000000000000000000000000000..22ea67ac5fba1957c37f6b219d38372520341da5
--- /dev/null
+++ b/tikz/harmonic_modes.tex
@@ -0,0 +1,61 @@
+\definecolor{green}{rgb}{0,.67,0}%
+%\tikzsetnextfilename{harmonic_modes}%
+\begin{tikzpicture}
+  \begin{axis}[
+      name = omega1,
+      width = 12cm,
+      height = 4.5cm,
+      xmode = log,
+      ymode = log,
+      ymin = 1e-6, ymax = 0.2,
+      xmin = 2.5,  xmax = 100,
+      xtick = {5,10,20,40,80},
+      minor xtick = {2,3,4,6,7,8,9,30,50,60,70,90},
+      xticklabels = {$5$, $10$, $20$, $40$, $80$},
+      %log x ticks with fixed point,
+      legend style = {
+        %font=\footnotesize,
+        at = {(axis description cs:0.03, 0.03)},
+        anchor = south west,
+        cells = {anchor=west},
+        draw = none,
+        fill = none,
+      },
+      xlabel = {Fourier mode $n$},
+      ylabel = {$|I_n/I_1|$},
+      %label style = {font=\small, anchor=center},
+      %every tick label/.append style = {font=\small},
+      %xlabel shift = 1ex,
+      %ylabel shift = 1.1ex,
+      scale only axis,
+    ]
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode17] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode34] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode69] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode138] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    \addplot[black, mark size=3pt, only marks, mark=x, forget plot] table [x=n, y=mode275] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %
+    \addplot[magenta, mark size=3pt, only marks, mark=square] table [x=n, y=i] {../figdata/harmonic_modes_vac5-omega1.dat};
+    \addlegendentry{$\vac=\hphantom{1}5\Omega$};
+    \addplot[blue, mark size=3pt, only marks, mark=diamond] table [x=n, y=i] {../figdata/harmonic_modes_vac10-omega1.dat};
+    \addlegendentry{$\vac=10\Omega$};
+    \addplot[green, mark size=3pt, only marks, mark=triangle] table [x=n, y=i] {../figdata/harmonic_modes_vac20-omega1.dat};
+    \addlegendentry{$\vac=20\Omega$};
+    \addplot[orange, mark size=3pt, only marks, mark=asterisk] table [x=n, y=i] {../figdata/harmonic_modes_vac40-omega1.dat};
+    \addlegendentry{$\vac=40\Omega$};
+    \addplot[red, mark size=3pt, only marks, mark=+] table [x=n, y=i] {../figdata/harmonic_modes_vac80-omega1.dat};
+    \addlegendentry{$\vac=80\Omega$};
+    \addplot[black, mark size=3pt, only marks, mark=x] coordinates {(6,10)}; % stupid dummy entry for the legend
+    \addlegendentry{$\vac=80\Omega$, adiabatic};
+    \addplot[black, mark size=1.5pt, only marks, mark=x, forget plot] table [x=n, y=mode5508] {../figdata/harmonic_modes_adiabatic-omega20.dat};
+    %
+    \addplot[magenta, mark size=1.5pt, only marks, mark=square] table [x=n, y=i] {../figdata/harmonic_modes_vac5-omega20.dat};
+    \addplot[blue, mark size=1.5pt, only marks, mark=diamond] table [x=n, y=i] {../figdata/harmonic_modes_vac10-omega20.dat};
+    \addplot[green, mark size=1.5pt, only marks, mark=triangle] table [x=n, y=i] {../figdata/harmonic_modes_vac20-omega20.dat};
+    \addplot[orange, mark size=1.5pt, only marks, mark=asterisk] table [x=n, y=i] {../figdata/harmonic_modes_vac40-omega20.dat};
+    \addplot[red, mark size=1.5pt, only marks, mark=+] table [x=n, y=i] {../figdata/harmonic_modes_vac80-omega20.dat};
+    %\addplot[black, mark size=1.5pt, only marks, mark=x] coordinates {(6,10)}; % stupid dummy entry for the legend
+    %\node[black,anchor=north east] at (axis description cs:0.96,0.92) {$\Omega=20\tkv,~\vdc=0$};
+    %\node[black] at (axis description cs:0.11,0.9) {(a)};
+  \end{axis}
+\end{tikzpicture}
diff --git a/tikz/harmonic_modes_separate.tex b/tikz/harmonic_modes_separate.tex
new file mode 100644
index 0000000000000000000000000000000000000000..be1c135969903e0201552c2d5a3434c2e2594a20
--- /dev/null
+++ b/tikz/harmonic_modes_separate.tex
@@ -0,0 +1,85 @@
+\definecolor{darkgreen}{rgb}{0,0.67,0}
+%\tikzsetnextfilename{harmonic_modes}%
+\begin{tikzpicture}
+  \begin{axis}[
+      name = omega1,
+      width = 6.5cm,
+      height = 4cm,
+      xmode = log,
+      ymode = log,
+      ymin = 1e-6, ymax = 0.2,
+      xmin = 2.5,  xmax = 100,
+      xtick = {5,10,20,40,80},
+      minor xtick = {2,3,4,6,7,8,9,30,50,60,70,90},
+      xticklabels = {$5$, $10$, $20$, $40$, $80$},
+      %log x ticks with fixed point,
+      legend style = {
+        %font=\footnotesize,
+        at = {(axis description cs:0.03, 0.03)},
+        anchor = south west,
+        cells = {anchor=west},
+        draw = none,
+        fill = none,
+      },
+      xlabel = {Fourier mode $n$},
+      ylabel = {$|I_n/I_1|$},
+      %label style = {font=\small, anchor=center},
+      %every tick label/.append style = {font=\small},
+      %xlabel shift = 1ex,
+      %ylabel shift = 1.1ex,
+      scale only axis,
+    ]
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode17] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode34] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode69] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %\addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode138] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    \addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode275] {../figdata/harmonic_modes_adiabatic-omega1.dat};
+    %
+    \addplot[magenta, only marks, mark=square] table [x=n, y=i] {../figdata/harmonic_modes_vac5-omega1.dat};
+    \addlegendentry{$\vac=\hphantom{1}5\Omega$};
+    \addplot[blue, only marks, mark=diamond] table [x=n, y=i] {../figdata/harmonic_modes_vac10-omega1.dat};
+    \addlegendentry{$\vac=10\Omega$};
+    \addplot[darkgreen, only marks, mark=triangle] table [x=n, y=i] {../figdata/harmonic_modes_vac20-omega1.dat};
+    \addlegendentry{$\vac=20\Omega$};
+    \addplot[orange, only marks, mark=asterisk] table [x=n, y=i] {../figdata/harmonic_modes_vac40-omega1.dat};
+    \addlegendentry{$\vac=40\Omega$};
+    \addplot[red, only marks, mark=+] table [x=n, y=i] {../figdata/harmonic_modes_vac80-omega1.dat};
+    \addlegendentry{$\vac=80\Omega$};
+    \addplot[black, only marks, mark=x] coordinates {(6,10)}; % stupid dummy entry for the legend
+    \addlegendentry{$\vac=80\Omega$, adiabatic};
+    \node[black,anchor=north east] at (axis description cs:0.96,0.92) {$\Omega=\tkv,~\vdc=0$};
+    %\node[black] at (axis description cs:0.11,0.9) {(a)};
+  \end{axis}
+  \begin{axis}[
+      at = (omega1.east),
+      anchor = west,
+      width = 6.5cm,
+      height = 4cm,
+      xmode = log,
+      ymode = log,
+      ymin = 1e-6, ymax = 0.2,
+      xmin = 2.5,  xmax = 100,
+      xtick = {5,10,20,40,80},
+      minor xtick = {2,3,4,6,7,8,9,30,50,60,70,90},
+      xticklabels = {$5$, $10$, $20$, $40$, $80$},
+      xlabel = {Fourier mode $n$},
+      %ylabel = {$|I_n/I_1|$},
+      yticklabels = \empty,
+      %label style = {font=\small, anchor=center},
+      %every tick label/.append style = {font=\small},
+      %xlabel shift = 1ex,
+      %ylabel shift = 1.1ex,
+      scale only axis,
+    ]
+    \addplot[black, only marks, mark=x, forget plot] table [x=n, y=mode5508] {../figdata/harmonic_modes_adiabatic-omega20.dat};
+    %
+    \addplot[magenta, only marks, mark=square] table [x=n, y=i] {../figdata/harmonic_modes_vac5-omega20.dat};
+    \addplot[blue, only marks, mark=diamond] table [x=n, y=i] {../figdata/harmonic_modes_vac10-omega20.dat};
+    \addplot[darkgreen, only marks, mark=triangle] table [x=n, y=i] {../figdata/harmonic_modes_vac20-omega20.dat};
+    \addplot[orange, only marks, mark=asterisk] table [x=n, y=i] {../figdata/harmonic_modes_vac40-omega20.dat};
+    \addplot[red, only marks, mark=+] table [x=n, y=i] {../figdata/harmonic_modes_vac80-omega20.dat};
+    \addplot[black, only marks, mark=x] coordinates {(6,10)}; % stupid dummy entry for the legend
+    \node[black,anchor=north east] at (axis description cs:0.96,0.92) {$\Omega=20\tkv,~\vdc=0$};
+    %\node[black] at (axis description cs:0.11,0.9) {(a)};
+  \end{axis}
+\end{tikzpicture}