diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py
index 7a8f3e899..e80dbf3ac 100644
--- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py
+++ b/python/lsst/pipe/tasks/computeExposureSummaryStats.py
@@ -950,7 +950,8 @@ def compute_magnitude_limit(
     # Total background counts derived from Eq. 45 of LSE-40
     background = (skyBg/gain + sigma_inst**2) * neff
     # Source counts to achieve the desired SNR (from quadratic formula)
-    source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain) + snr**2 * background)
+    # Note typo in gain exponent in Eq. 45 of LSE-40 (DM-53234)
+    source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain**2) + snr**2 * background)
 
     # Convert to a magnitude using the zeropoint.
     # Note: Zeropoint currently includes exposure time
diff --git a/tests/test_computeExposureSummaryStats.py b/tests/test_computeExposureSummaryStats.py
index 8d10e6c49..976abacf8 100644
--- a/tests/test_computeExposureSummaryStats.py
+++ b/tests/test_computeExposureSummaryStats.py
@@ -36,6 +36,15 @@
 from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsTask
 from lsst.pipe.tasks.computeExposureSummaryStats import compute_magnitude_limit
 
+# Values from syseng_throughputs notebook assuming 30s exposure
+# consisting of 2x15s snaps each with readnoise of 9e-
+fwhm_eff_fid = {'g': 0.87, 'r': 0.83, 'i': 0.80}
+skycounts_fid = {'g': 463.634122, 'r': 988.626863, 'i': 1588.280513}
+zeropoint_fid = {'g': 28.508375, 'r': 28.360838, 'i': 28.171396}
+readnoise_fid = {'g': 12.73, 'r': 12.73, 'i': 12.73}
+# Output magnitude limit from syseng_throughputs notebook
+m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10}
+
 
 class ComputeExposureSummaryTestCase(lsst.utils.tests.TestCase):
 
@@ -167,17 +176,11 @@ def testComputeMagnitudeLimit(self):
         """Test the magnitude limit calculation using fiducials from SMTN-002
         and syseng_throughputs."""
 
-        # Values from syseng_throughputs notebook assuming 30s exposure
-        # consisting of 2x15s snaps each with readnoise of 9e-
-        fwhm_eff_fid = {'g': 0.87, 'r': 0.83, 'i': 0.80}
-        skycounts_fid = {'g': 463.634122, 'r': 988.626863, 'i': 1588.280513}
-        zeropoint_fid = {'g': 28.508375, 'r': 28.360838, 'i': 28.171396}
-        readnoise_fid = {'g': 12.73, 'r': 12.73, 'i': 12.73}
         # Assumed values from SMTN-002
         snr = 5
         gain = 1.0
-        # Output magnitude limit from syseng_throughputs notebook
-        m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10}
+        # For testing non-unity gain
+        gain_test = 2.0
 
         for band in ['g', 'r', 'i']:
             # Translate into DM quantities
@@ -190,6 +193,12 @@ def testComputeMagnitudeLimit(self):
             m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, readNoise, gain, snr)
             self.assertFloatsAlmostEqual(m5, m5_fid[band], atol=1e-2)
 
+            # Changing the gain consistently should not change M5
+            m5_gain = compute_magnitude_limit(psfArea, skyBg/gain_test,
+                                              zeroPoint - 2.5*np.log10(gain_test),
+                                              readNoise/gain_test, gain_test, snr)
+            self.assertFloatsAlmostEqual(m5_gain, m5, atol=1e-6)
+
         # Check that input NaN lead to output NaN
         nan = float('nan')
         m5 = compute_magnitude_limit(nan, skyBg, zeroPoint, readNoise, gain, snr)