diff --git a/optika/_tests/test_systems.py b/optika/_tests/test_systems.py
index 0ea04c04..14093dc5 100644
--- a/optika/_tests/test_systems.py
+++ b/optika/_tests/test_systems.py
@@ -50,7 +50,7 @@ def test_image(
         assert isinstance(result.outputs, na.AbstractScalar)
         assert np.all(result.inputs.wavelength > 0 * u.nm)
         assert na.unit_normalized(result.inputs.position).is_equivalent(u.mm)
-        assert result.outputs.sum() > (0 * u.electron)
+        assert result.outputs.sum() != (0 * u.electron)
 
 
 class AbstractTestAbstractSequentialSystem(
diff --git a/optika/sensors/_materials/_materials.py b/optika/sensors/_materials/_materials.py
index cea766dc..9b56eb9e 100644
--- a/optika/sensors/_materials/_materials.py
+++ b/optika/sensors/_materials/_materials.py
@@ -836,7 +836,7 @@ def signal(
     thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant,
     cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface,
     temperature: u.Quantity | na.ScalarArray = 300 * u.K,
-    method: Literal["exact", "approx"] = "exact",
+    method: Literal["exact", "approx", "expected"] = "exact",
     shape_random: None | dict[str, int] = None,
 ) -> na.AbstractScalar:
     r"""
@@ -875,8 +875,12 @@ def signal(
         The temperature of the light-sensitive silicon layer.
     method
         The method used to generate random samples of measured electrons.
-        The exact method is more accurate for low numbers of photons,
-        but suffers from poor performance for high numbers of photons.
+        The `exact` method simulates every photon so it is accurate for low
+        photon counts but slow for high photon counts.
+        The `approx` method is much faster, but is only accurate if the number
+        of photons is high.
+        The `expected` method does not add any noise to the signal and just
+        returns the expected number of electrons.
     shape_random
         Additional shape used to specify the number of samples to draw.
 
@@ -930,6 +934,18 @@ def signal(
     if absorption is None:
         absorption = optika.chemicals.Chemical("Si").absorption(wavelength)
 
+    if method == "expected":
+        iqy = quantum_yield_ideal(
+            wavelength=wavelength,
+            temperature=temperature,
+        )
+        cce = charge_collection_efficiency(
+            absorption=absorption,
+            thickness_implant=thickness_implant,
+            cce_backsurface=cce_backsurface,
+        )
+        return iqy * absorbance * cce * photons_expected.to(u.ph)
+
     photons_absorbed_expected = absorbance * photons_expected.to(u.ph)
     photons_absorbed = na.random.poisson(
         lam=photons_absorbed_expected,
@@ -967,6 +983,7 @@ def signal(
         self,
         rays: optika.rays.RayVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
+        noise: bool = True,
     ) -> optika.rays.RayVectorArray:
         """
         Given a set of incident rays, compute the number of electrons
@@ -980,6 +997,8 @@ def signal(
             either be in units of photons or energy.
         normal
             The vector perpendicular to the surface of the sensor.
+        noise
+            Whether to add noise to the result.
         """
 
     @abc.abstractmethod
@@ -1032,14 +1051,15 @@ class IdealImagingSensorMaterial(
     AbstractImagingSensorMaterial,
 ):
     """
-    An idealized sensor material with a quantum efficiency of unity
-    and no charge diffusion.
+    An idealized sensor material with a quantum efficiency of unity,
+    no charge diffusion, and no noise model.
     """
 
     def signal(
         self,
         rays: optika.rays.RayVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
+        noise: bool = False,
     ) -> optika.rays.RayVectorArray:
 
         intensity = rays.intensity
@@ -1047,6 +1067,10 @@ def signal(
             h = astropy.constants.h
             c = astropy.constants.c
             intensity = intensity / (h * c / rays.wavelength) * u.photon
+
+        if noise:
+            intensity = na.random.poisson(intensity.to(u.ph)).astype(int)
+
         electrons = intensity * u.electron / u.photon
         electrons = electrons.to(u.electron)
 
@@ -1373,6 +1397,7 @@ def signal(
         self,
         rays: optika.rays.RayVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
+        noise: bool = True,
     ) -> optika.rays.RayVectorArray:
 
         intensity = rays.intensity
@@ -1383,6 +1408,11 @@ def signal(
             c = astropy.constants.c
             intensity = intensity / (h * c / wavelength) * u.photon
 
+        if noise:
+            method = "exact"
+        else:
+            method = "expected"
+
         electrons = signal(
             photons_expected=intensity,
             wavelength=wavelength,
@@ -1391,6 +1421,7 @@ def signal(
             thickness_implant=self.thickness_implant,
             cce_backsurface=self.cce_backsurface,
             temperature=self.temperature,
+            method=method,
         )
 
         result = dataclasses.replace(rays, intensity=electrons)
diff --git a/optika/sensors/_materials/_materials_test.py b/optika/sensors/_materials/_materials_test.py
index 9ffc3c83..88304e03 100644
--- a/optika/sensors/_materials/_materials_test.py
+++ b/optika/sensors/_materials/_materials_test.py
@@ -284,13 +284,18 @@ class AbstractTestAbstractImagingSensorMaterial(
             na.Cartesian3dVectorArray(0, 0, -1),
         ],
     )
+    @pytest.mark.parametrize(
+        argnames="noise",
+        argvalues=[True, False],
+    )
     def test_signal(
         self,
         a: optika.sensors.AbstractImagingSensorMaterial,
         rays: optika.rays.RayVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
+        noise: bool,
     ):
-        result = a.signal(rays, normal)
+        result = a.signal(rays, normal, noise=noise)
         assert isinstance(result, optika.rays.RayVectorArray)
         assert np.all(result.intensity >= 0 * u.electron)
 
diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py
index 1cf2b1ee..2792802f 100644
--- a/optika/sensors/_sensors.py
+++ b/optika/sensors/_sensors.py
@@ -89,6 +89,7 @@ def readout(
         timedelta: None | u.Quantity | na.AbstractScalar = None,
         axis: None | str | Sequence[str] = None,
         where: bool | na.AbstractScalar = True,
+        noise: bool = True,
     ) -> na.FunctionArray[
         na.SpectralPositionalVectorArray,
         na.AbstractScalar,
@@ -113,6 +114,8 @@ def readout(
         where
             A boolean mask used to indicate which photons should be considered
             when calculating the signal measured by the sensor.
+        noise
+            Whether to add noise to the result
         """
         if timedelta is None:
             timedelta = self.timedelta_exposure
@@ -131,6 +134,7 @@ def readout(
         rays = self.material.signal(
             rays=rays,
             normal=normal,
+            noise=noise,
         )
 
         rays = dataclasses.replace(
diff --git a/optika/systems.py b/optika/systems.py
index fe445774..cc885219 100644
--- a/optika/systems.py
+++ b/optika/systems.py
@@ -841,6 +841,7 @@ def image(
         axis_wavelength: None | str = None,
         axis_field: None | tuple[str, str] = None,
         axis_pupil: None | tuple[str, str] = None,
+        noise: bool = True,
         **kwargs,
     ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]:
         """
@@ -875,6 +876,8 @@ def image(
             If :obj:`None`,
             ``set(pupil.shape) - set(self.shape) - {axis_wavelength,} - set(axis_field)``,
             should have exactly two elements.
+        noise
+            Whether to add noise to the result.
         kwargs
             Additional keyword arguments used by subclass implementations
             of this method.
@@ -946,6 +949,7 @@ def image(
             rays=rayfunction.outputs,
             wavelength=wavelength,
             axis=(axis_wavelength,) + axis_field + axis_pupil,
+            noise=noise,
         )
 
     def plot(