diff --git a/etl/prs/prs_analytical_representations.py b/etl/prs/prs_analytical_representations.py
index 09e51956148c141385082634e7d13835aeff7089..942e1e64a6ebf65569eed275530cac50e8c548ad 100644
--- a/etl/prs/prs_analytical_representations.py
+++ b/etl/prs/prs_analytical_representations.py
@@ -69,7 +69,11 @@ def main(ratios_to_fit: Union[List[str], None] = None):
         plt.scatter(data_clean[_column_to_fit], data_clean['avg_nh2'], marker='+', alpha=0.1, color='red')
         for components in x0[ratio_string]:
             plt.plot(x_reg, approx(x_reg, *components), color='cyan')
-        plt.show()
+        plt.savefig(os.path.join('..',
+        'publications',
+        '6373bb408e4040043398e495',
+        'referee',
+        f'analytical_expressions_comparison_{ratio_string}.png'))
 
 
 if __name__ == '__main__':
diff --git a/etl/prs/prs_poc_figures.py b/etl/prs/prs_poc_figures.py
index 40993c8c0c9828777dcc66261e693e60b28c8100..8f3a8f8f55f60c76cca5b5146a56bac3f7c3e5f5 100644
--- a/etl/prs/prs_poc_figures.py
+++ b/etl/prs/prs_poc_figures.py
@@ -86,7 +86,7 @@ def make_volume_densities_comparison_figure(df_full: pd.DataFrame):
                  df_merge['best_fit'],
                  yerr=uncertainty_array,
                  fmt='none')
-    plt.plot([2e4, 8e5], [2e4, 8e6], color='red')
+    plt.plot([2e4, 2e6], [2e4, 2e6], color='red')
     plt.loglog()
     plt.xlabel('<$n$(H$_{2, dust}$)> [cm$^{-3}$]')
     plt.ylabel('<$n$(H$_{2, SAK}$)> [cm$^{-3}$]')
@@ -98,10 +98,10 @@ def make_volume_densities_comparison_figure(df_full: pd.DataFrame):
 def main():
     df_full = get_poc_results(line_fit_filename='ch3oh_data_top35.csv')
     df_full.rename(columns={'class_phys': 'Class'}, inplace=True)
+    make_volume_densities_comparison_figure(df_full=df_full)
     df_full = df_full[(df_full['mass'] < 10000) & (df_full['mass'] > 300)]
     make_ecdf_figure(df_full=df_full)
     make_violin_figure(df_full=df_full)
-    make_volume_densities_comparison_figure(df_full=df_full)
 
 
 if __name__ == '__main__':
diff --git a/etl/prs/prs_poc_latex_table.py b/etl/prs/prs_poc_latex_table.py
index b0a57490016dba15eeb7c0af62b1ae4b8bf64b3c..30dcdf04156a7b93cce8b7063b23f1e6a74e4b5f 100644
--- a/etl/prs/prs_poc_latex_table.py
+++ b/etl/prs/prs_poc_latex_table.py
@@ -61,21 +61,21 @@ remap_columns = OrderedDict([
     ('hpd_interval', ('67\%  HPD', density_format)),
     ('mass', ('Mass', '$[10^2 \mathrm{M_\odot}]$')),
     ('distance', ('Distance', '[kpc]')),
-    ('tpeak_86', ('$T_{MB, 86}$', '[K]')),
+    ('tpeak_86', ('$T_{MB,(2_{-1}-1_{-1})}$', '[K]')),
     ('linewidth_86', ('$FWHM_{(2_K-1_K)}$', fwhm_units)),
-    ('tpeak_87', ('$T_{MB, 87}$', '[K]')),
+    ('tpeak_87', ('$T_{MB,(2_{0}-1_{0})}$', '[K]')),
     ('linewidth_87', ('$FWHM_{87}$', '[km s$^{-1}]$')),
-    ('tpeak_88', ('$T_{MB, 88}$', '[K]')),
+    ('tpeak_88', ('$T_{MB,(2_{1}-1_{1})}$', '[K]')),
     ('linewidth_88', ('$FWHM_{88}$', '[km s$^{-1}]$')),
     # ('rms_noise_86', ('$\sigma_{96.7GHz}$', '[K]')),
-    ('tpeak_256', ('$T_{MB, 256}$', '[K]')),
+    ('tpeak_256', ('$T_{MB,(5_{0}-4_{0})}$', '[K]')),
     ('linewidth_256', ('$FWHM_{256}$', fwhm_units)),
-    ('tpeak_257', ('$T_{MB, 257}$', '[K]')),
+    ('tpeak_257', ('$T_{MB,(5_{-1}-4_{-1})}$', '[K]')),
     ('linewidth_257', ('$FWHM_{(5_K-4_K)}$', fwhm_units)),
     # ('rms_noise_256', ('$\sigma_{241.7GHz}$', '[K]')),
-    ('tpeak_380', ('$T_{MB, 380}$', '[K]')),
+    ('tpeak_380', ('$T_{MB,(7_{0}-6_{0})}$', '[K]')),
     ('linewidth_380', ('$FWHM_{380}$', fwhm_units)),
-    ('tpeak_381', ('$T_{MB, 381}$', '[K]')),
+    ('tpeak_381', ('$T_{MB,(7_{-1}-6_{-1})}$', '[K]')),
     ('linewidth_381', ('$FWHM_{(7_K-6_K)}$', fwhm_units)),
     # ('rms_noise_380', ('$\sigma_{338.3GHz}$', '[K]')),
 ])
@@ -98,14 +98,14 @@ with open(os.path.join('prs', 'output', 'poc_tables', 'fit_results.tex'), 'w') a
     outfile.write(latex_table)
 
 table_cols = [
-    [('Source', ''), ('$T_{MB, 86}$', '[K]'), ('$T_{MB, 87}$', '[K]'), ('$T_{MB, 88}$', '[K]'), ('$FWHM_{(2_K-1_K)}$',
+    [('Source', ''), ('$T_{MB,(2_{-1}-1_{-1})}$', '[K]'), ('$T_{MB,(2_{0}-1_{0})}$', '[K]'), ('$T_{MB,(2_{1}-1_{1})}$', '[K]'), ('$FWHM_{(2_K-1_K)}$',
                                                                                                  fwhm_units)],
-    [('Source', ''), ('$T_{MB, 256}$', '[K]'), ('$T_{MB, 257}$', '[K]'), ('$FWHM_{(5_K-4_K)}$', fwhm_units), ('$T_{MB, 380}$', '[K]'), ('$T_{MB, 381}$', '[K]'), ('$FWHM_{(7_K-6_K)}$',
+    [('Source', ''), ('$T_{MB,(5_{0}-4_{0})}$', '[K]'), ('$T_{MB,(5_{-1}-4_{-1})}$', '[K]'), ('$FWHM_{(5_K-4_K)}$', fwhm_units), ('$T_{MB,(7_{0}-6_{0})}$', '[K]'), ('$T_{MB,(7_{-1}-6_{-1})}$', '[K]'), ('$FWHM_{(7_K-6_K)}$',
                                                                                                                                                          fwhm_units), ]
 ]
 captions = [
-    'Line properties of the lines in the $(2_K-1_K)$ methanol band. Only one FWHM is listed because the fit is performed forcing all lines to have the same width. A non-detection is indicated with three dots.',
-    'Line properties of the lines in the $(5_K-4_K)$ and $(7_K-6_K)$ methanol bands. Only one FWHM is listed per band because the fit is performed forcing all lines to have the same width. A non-detection is indicated with three dots, while missing data are indicated with \'N/A\'.',
+    "Line properties of the lines in the $(2_K-1_K)$ methanol band. Only one FWHM is listed because the fit is performed forcing all lines to have the same width. The main-beam temperature of the lines is indicated as $T_{MB,J_K-J',_K'}$. A non-detection is indicated with three dots.",
+    "Line properties of the lines in the $(5_K-4_K)$ and $(7_K-6_K)$ methanol bands. Only one FWHM is listed per band because the fit is performed forcing all lines to have the same width. The main-beam temperature of the lines is indicated as $T_{MB,J_K-J',_K'}$. A non-detection is indicated with three dots, while missing data are indicated with \'N/A\'.",
 ]
 labels = [
     'tab:poc_lines_3mm',