diff --git a/autolens/aggregator/__init__.py b/autolens/aggregator/__init__.py
index 4ba8ba7d1..ebe51e816 100644
--- a/autolens/aggregator/__init__.py
+++ b/autolens/aggregator/__init__.py
@@ -27,5 +27,6 @@
 from autolens.aggregator.subplot import SubplotFitX1Plane as subplot_fit_x1_plane
 from autolens.aggregator.subplot import SubplotFit as subplot_fit
 from autolens.aggregator.subplot import SubplotFitLog10 as subplot_fit_log10
+from autolens.aggregator.subplot import FITSModelGalaxyImages as fits_model_galaxy_images
 from autolens.aggregator.subplot import FITSTracer as fits_tracer
 from autolens.aggregator.subplot import FITSFit as fits_fit
diff --git a/autolens/aggregator/subplot.py b/autolens/aggregator/subplot.py
index 5ff2a926c..c12aa6ee3 100644
--- a/autolens/aggregator/subplot.py
+++ b/autolens/aggregator/subplot.py
@@ -1,6 +1,15 @@
 from enum import Enum
 
 
+class FITSModelGalaxyImages(Enum):
+    """
+    The HDUs that can be extracted from the fit.fits file.
+    """
+
+    lens_light_image = "GALAXY_0"
+    lensed_source_image = "GALAXY_1"
+
+
 class FITSTracer(Enum):
     """
     The HDUs that can be extracted from the fit.fits file.
diff --git a/autolens/imaging/model/analysis.py b/autolens/imaging/model/analysis.py
index a03528f75..e7442ea04 100644
--- a/autolens/imaging/model/analysis.py
+++ b/autolens/imaging/model/analysis.py
@@ -1,5 +1,4 @@
 import logging
-import numpy as np
 
 import autofit as af
 import autogalaxy as ag
diff --git a/autolens/imaging/plot/fit_imaging_plotters.py b/autolens/imaging/plot/fit_imaging_plotters.py
index 4efb0d6ce..20607715e 100644
--- a/autolens/imaging/plot/fit_imaging_plotters.py
+++ b/autolens/imaging/plot/fit_imaging_plotters.py
@@ -428,7 +428,7 @@ def subplot_fit_x1_plane(self):
         self.mat_plot_2d.cmap.kwargs.pop("vmax")
 
         self.residuals_symmetric_cmap = False
-        self.set_title(label="Lens Light Subtracted Image")
+        self.set_title(label="Lens Light Subtracted")
         self.figures_2d(normalized_residual_map=True)
 
         self.mat_plot_2d.cmap.kwargs["vmin"] = 0.0
@@ -468,7 +468,7 @@ def subplot_fit_log10_x1_plane(self):
         self.mat_plot_2d.cmap.kwargs.pop("vmax")
 
         self.residuals_symmetric_cmap = False
-        self.set_title(label="Lens Light Subtracted Image")
+        self.set_title(label="Lens Light Subtracted")
         self.figures_2d(normalized_residual_map=True)
 
         self.residuals_symmetric_cmap = True
@@ -519,7 +519,7 @@ def subplot_fit(self, plane_index: Optional[int] = None):
 
         self.mat_plot_2d.cmap.kwargs["vmin"] = 0.0
 
-        self.set_title(label="Lens Light Subtracted Image")
+        self.set_title(label="Lens Light Subtracted")
         self.figures_2d_of_planes(
             plane_index=plane_index,
             subtracted_image=True,
@@ -527,7 +527,7 @@ def subplot_fit(self, plane_index: Optional[int] = None):
             remove_critical_caustic=True,
         )
 
-        self.set_title(label="Source Model Image (Image Plane)")
+        self.set_title(label="Source Model Image")
         self.figures_2d_of_planes(
             plane_index=plane_index,
             model_image=True,
@@ -551,7 +551,7 @@ def subplot_fit(self, plane_index: Optional[int] = None):
         self.mat_plot_2d.cmap.kwargs["vmin"] = -1.0
         self.mat_plot_2d.cmap.kwargs["vmax"] = 1.0
 
-        self.set_title(label="Normalized Residual Map (1 sigma)")
+        self.set_title(label=r"Normalized Residual Map $1\sigma$")
         self.figures_2d(normalized_residual_map=True)
         self.set_title(label=None)
 
@@ -624,12 +624,12 @@ def subplot_fit_log10(self, plane_index: Optional[int] = None):
 
         self.mat_plot_2d.cmap.kwargs["vmin"] = 0.0
 
-        self.set_title(label="Lens Light Subtracted Image")
+        self.set_title(label="Lens Light Subtracted")
         self.figures_2d_of_planes(
             plane_index=plane_index, subtracted_image=True, use_source_vmax=True, remove_critical_caustic=True
         )
 
-        self.set_title(label="Source Model Image (Image Plane)")
+        self.set_title(label="Source Model Image")
         self.figures_2d_of_planes(
             plane_index=plane_index, model_image=True, use_source_vmax=True, remove_critical_caustic=True
         )
@@ -652,7 +652,7 @@ def subplot_fit_log10(self, plane_index: Optional[int] = None):
         self.mat_plot_2d.cmap.kwargs["vmin"] = -1.0
         self.mat_plot_2d.cmap.kwargs["vmax"] = 1.0
 
-        self.set_title(label="Normalized Residual Map (1 sigma)")
+        self.set_title(label=r"Normalized Residual Map $1\sigma$")
         self.figures_2d(normalized_residual_map=True)
         self.set_title(label=None)
 
diff --git a/autolens/interferometer/plot/fit_interferometer_plotters.py b/autolens/interferometer/plot/fit_interferometer_plotters.py
index 2ca412515..eff749b5b 100644
--- a/autolens/interferometer/plot/fit_interferometer_plotters.py
+++ b/autolens/interferometer/plot/fit_interferometer_plotters.py
@@ -420,7 +420,7 @@ def subplot_fit(self):
         self.mat_plot_2d.cmap.kwargs["vmin"] = -1.0
         self.mat_plot_2d.cmap.kwargs["vmax"] = 1.0
 
-        self.set_title(label="Normalized Residual Map (1 sigma)")
+        self.set_title(label=r"Normalized Residual Map $1\sigma$")
         self.figures_2d(dirty_normalized_residual_map=True)
         self.set_title(label=None)