diff --git a/fiasco/ions.py b/fiasco/ions.py
index 8a77225a..f904d72d 100644
--- a/fiasco/ions.py
+++ b/fiasco/ions.py
@@ -2010,7 +2010,6 @@ def two_photon(self,
         """
         wavelength = np.atleast_1d(wavelength)
         electron_density = np.atleast_1d(electron_density)
-
         final_shape = self.temperature.shape + electron_density.shape + wavelength.shape
         couple_density_to_temperature = kwargs.setdefault('couple_density_to_temperature', False)
         if couple_density_to_temperature:
@@ -2028,8 +2027,9 @@ def two_photon(self,
         else:
             return u.Quantity(np.zeros(final_shape),  'erg cm^3 s^-1 Angstrom^-1')
         level_index = np.where(self.levels.label==label)
+        if level_index[0].shape != (1,):
+            raise ValueError(f'Multiple levels {self.levels.level[level_index]} for label {label}.')
         rest_wavelength = self.levels.energy[level_index].to('AA', equivalencies=u.equivalencies.spectral())
-
         # NOTE: Explicitly setting the boundary condition type here to match the behavior of the
         # IDL spline interpolation functions. See https://github.com/wtbarnes/fiasco/pull/297 for
         # additional details.
@@ -2037,18 +2037,18 @@ def two_photon(self,
         x_new = (rest_wavelength / wavelength).decompose().to_value(u.dimensionless_unscaled)
         psi_interp = cubic_spline(x_new)
         psi_interp = np.where(x_new>1.0, 0.0, psi_interp)
-
         energy_dist = (A_ji * rest_wavelength * psi_interp) / (psi_norm * wavelength**3)
         # NOTE: There are known issues when including the proton rates here for some ions so these
         # are excluded by default. See https://github.com/wtbarnes/fiasco/pull/260#issuecomment-1955770878
         # for more details.
         kwargs.setdefault('include_protons', False)
         level_population = self.level_populations(electron_density, **kwargs)
-        level_population = level_population[..., self.levels[level_index].level-1]
-
+        # NOTE: This is necessary as the level population grid is not necessarily the same as self.levels.level
+        # NOTE: Ensure idx is a scalar so that level_population is always 2D
+        idx = np.where(np.arange(1, self.n_levels+1, 1)==self.levels.level[level_index])[0][0]
+        level_population = level_population[..., idx]
         if couple_density_to_temperature:
             electron_density = electron_density[:, np.newaxis]
         level_density = level_population / electron_density
         matrix = np.outer(level_density, energy_dist).reshape(final_shape)
-
         return const.h * const.c * matrix
diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py
index cdde3d8e..2824a7c5 100644
--- a/fiasco/tests/test_ion.py
+++ b/fiasco/tests/test_ion.py
@@ -467,13 +467,14 @@ def test_free_bound_radiative_loss(ion, h1):
     assert u.isclose(ion.free_bound_radiative_loss()[20], 1.977311794139327e-23 * u.erg * u.cm**3 / u.s)
 
 
-# The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl
-@pytest.mark.requires_dbase_version('>= 8','<= 8.0.7')
 @pytest.mark.parametrize(('ion_fixture_name', 'include_protons', 'expected'), [
     ('c4', False, 0.0 * u.cm**3 * u.erg / u.Angstrom / u.s),
-    ('c5', False, 4.08382025e-25 * u.cm**3 * u.erg / u.Angstrom / u.s),
+    # The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl
+    pytest.param(
+        'c5', False, 4.08382025e-25 * u.cm**3 * u.erg / u.Angstrom / u.s,
+        marks=pytest.mark.requires_dbase_version('>= 8', '<= 8.0.7')
+    ),
     ('c6', False, 8.23471268e-26 * u.cm**3 * u.erg / u.Angstrom / u.s),
-    ('c6', True,  6.77829033e-29 * u.cm**3 * u.erg / u.Angstrom / u.s),
 ])
 def test_two_photon(request, ion_fixture_name, include_protons, expected):
     ion = request.getfixturevalue(ion_fixture_name)
@@ -482,6 +483,18 @@ def test_two_photon(request, ion_fixture_name, include_protons, expected):
     assert u.allclose(emission[30, 0, 0], expected)
 
 
+@pytest.mark.parametrize(('density', 'couple_density_to_temperature', 'expected_shape'), [
+    (1e9*u.cm**(-3), False, (1,)),
+    ([1e9,]*u.cm**(-3), False, (1,)),
+    ([1e8, 1e9, 1e10]*u.cm**(-3), False, (3,)),
+    (1e9*u.cm**(-3)*np.ones((100,)), True, (1,)),
+])
+def test_two_photon_density_shape(c5, density, couple_density_to_temperature, expected_shape):
+    wavelength = np.linspace(20, 100, 5) * u.AA
+    emission = c5.two_photon(wavelength, density, couple_density_to_temperature=couple_density_to_temperature)
+    assert emission.shape == c5.temperature.shape+expected_shape+wavelength.shape
+
+
 def test_free_bound_no_recombining(h1):
     # This is test the case where there is no data available for the recombining
     # ion (H 2)