diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index 23d89b06e5..b7cb44a5f4 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -18,6 +18,7 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add ``return_components`` kwarg to :py:func:`pvlib.irradiance.isotropic`. (:issue:`2750`, :pull:`2787`)
 * Ensure all timezones are available in all OSs. (:issue:`2795`, :pull:`2809`)
 
 
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 98bd948286..07d89a1a05 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -599,7 +599,7 @@ def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None):
     return diffuse_irrad
 
 
-def isotropic(surface_tilt, dhi):
+def isotropic(surface_tilt, dhi, return_components=False):
     r'''
     Determine diffuse irradiance from the sky on a tilted surface using
     the isotropic sky model.
@@ -623,10 +623,29 @@ def isotropic(surface_tilt, dhi):
     dhi : numeric
         Diffuse horizontal irradiance, must be >=0. See :term:`dhi`.
 
+    return_components : bool, default ``False``
+        If ``False``, ``poa_sky_diffuse`` is returned.
+        If ``True``, ``diffuse_components`` is returned.
+        For this model, ``return_components`` contains the same data as
+        ``poa_sky_diffuse``, but it is included for consistency with the
+        other sky diffuse models.
+
     Returns
     -------
-    diffuse : numeric
-        The sky diffuse component of the solar radiation. [Wm⁻²]
+    numeric, dict, or DataFrame
+        Return type controlled by ``return_components`` argument.
+        If ``False``, ``poa_sky_diffuse`` is returned.
+        If ``True``, ``diffuse_components`` is returned.
+
+    poa_sky_diffuse : numeric
+        The sky diffuse component of irradiance on a tilted plane. [Wm⁻²]
+
+    diffuse_components : dict (array input) or DataFrame (Series input)
+        Keys/columns are:
+            * poa_sky_diffuse: The sky diffuse component of irradiance on a
+            tilted plane. [Wm⁻²]
+            * poa_isotropic: The portion of sky diffuse irradiance on a tilted
+            plane from the isotropic sky dome. [Wm⁻²]
 
     References
     ----------
@@ -640,9 +659,20 @@ def isotropic(surface_tilt, dhi):
        Energy vol. 201. pp. 8-12
        :doi:`10.1016/j.solener.2020.02.067`
     '''
-    sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5
+    poa_sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5
 
-    return sky_diffuse
+    if return_components:
+        diffuse_components = {
+            'poa_sky_diffuse': poa_sky_diffuse,
+            'poa_isotropic': poa_sky_diffuse
+        }
+
+        if isinstance(poa_sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+
+        return diffuse_components
+    else:
+        return poa_sky_diffuse
 
 
 def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index fbe9906358..ff51b5a306 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -168,6 +168,32 @@ def test_isotropic_series(irrad_data):
     assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4)
 
 
+def test_isotropic_components(irrad_data):
+    keys = ['poa_sky_diffuse', 'poa_isotropic']
+    expected = pd.DataFrame(np.array(
+        [[0, 35.728402, 104.601328, 54.777191],
+         [0, 35.728402, 104.601328, 54.777191]]).T,
+        columns=keys,
+        index=irrad_data.index
+    )
+    # pandas
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'], return_components=True)
+    assert_frame_equal(result, expected, check_less_precise=4)
+    # numpy
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'].to_numpy(), return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key], atol=1e-4)
+    assert isinstance(result, dict)
+    # scalar
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'].to_numpy()[-1], return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key].iloc[-1], atol=1e-4)
+    assert isinstance(result, dict)
+
+
 def test_klucher_series_float():
     # klucher inputs
     surface_tilt, surface_azimuth = 40.0, 180.0