diff --git a/astropy/cosmology/_utils.py b/astropy/cosmology/_utils.py
index 84617b29e76..cb1f42ebcfd 100644
--- a/astropy/cosmology/_utils.py
+++ b/astropy/cosmology/_utils.py
@@ -9,7 +9,7 @@
 from collections.abc import Callable
 from dataclasses import Field
 from numbers import Number
-from typing import TYPE_CHECKING, Any, TypeVar
+from typing import TYPE_CHECKING, Any, Generic, TypeVar, overload
 
 import numpy as np
 
@@ -20,8 +20,11 @@
 
 if TYPE_CHECKING:
     from astropy.cosmology import Parameter
+    from astropy.cosmology.core import Cosmology
+
 
 _F = TypeVar("_F", bound=Callable[..., Any])
+R = TypeVar("R")
 
 
 def vectorize_redshift_method(func=None, nin=1):
@@ -140,3 +143,55 @@ def _depr_kws(func: _F, /, kws: tuple[str, ...], since: str) -> _F:
     wrapper = _depr_kws_wrap(func, kws, since)
     functools.update_wrapper(wrapper, func)
     return wrapper
+
+
+class CachedInDictPropertyDescriptor(Generic[R]):
+    """Descriptor for a property that is cached in the instance's dictionary.
+
+    Note that this is a non-data descriptor, not NOT a data descriptor, so
+    after the property is accessed and cached with the same key as the
+    property, the instance's dictionary will have the property value
+    directly, and the descriptor will not be called again, until the
+    property is deleted from the instance's dictionary. See
+    https://docs.python.org/3/howto/descriptor.html#descriptor-protocol.
+    """
+
+    # __slots__ = ("fget", "name")  # TODO: when __doc__ is supported by __slots__
+
+    def __init__(self, fget: Callable[[Cosmology], R]) -> None:
+        self.fget = fget
+        self.__doc__ = fget.__doc__
+
+    def __set_name__(self, cosmo_cls: type[Cosmology], name: str) -> None:
+        self.name: str = name
+
+    @overload
+    def __get__(
+        self, cosmo: None, cosmo_cls: Any
+    ) -> CachedInDictPropertyDescriptor[R]: ...
+
+    @overload
+    def __get__(self, cosmo: Cosmology, cosmo_cls: Any) -> R: ...
+
+    def __get__(
+        self, cosmo: Cosmology | None, cosmo_cls: type[Cosmology] | None
+    ) -> R | CachedInDictPropertyDescriptor[R]:
+        # Accessed from the class, return the descriptor itself
+        if cosmo is None:
+            return self
+
+        # If the property is not in the instance's dictionary, calculate and store it.
+        if self.name not in cosmo.__dict__:
+            cosmo.__dict__[self.name] = self.fget(cosmo)
+
+        # Return the property value from the instance's dictionary
+        # This is only called once, thereafter the property value is accessed directly
+        # from the instance's dictionary.
+        return cosmo.__dict__[self.name]
+
+
+def cached_on_dict_property(
+    fget: Callable[[Cosmology], R],
+) -> CachedInDictPropertyDescriptor[R]:
+    """Descriptor for a property that is cached in the instance's dictionary."""
+    return CachedInDictPropertyDescriptor(fget)
diff --git a/astropy/cosmology/flrw/base.py b/astropy/cosmology/flrw/base.py
index 2cd7be91435..dd98dcdbece 100644
--- a/astropy/cosmology/flrw/base.py
+++ b/astropy/cosmology/flrw/base.py
@@ -21,6 +21,7 @@
 import astropy.units as u
 from astropy.cosmology._utils import (
     aszarr,
+    cached_on_dict_property,
     deprecated_keywords,
     vectorize_redshift_method,
 )
@@ -209,37 +210,6 @@ class FLRW(Cosmology, _ScaleFactorMixin):
     )
 
     def __post_init__(self):
-        # Derived quantities:
-        # Dark matter density; matter - baryons, if latter is not None.
-        object.__setattr__(
-            self, "_Odm0", (None if self._Ob0 is None else (self._Om0 - self._Ob0))
-        )
-
-        # 100 km/s/Mpc * h = H0 (so h is dimensionless)
-        object.__setattr__(self, "_h", self._H0.value / 100.0)
-        # Hubble distance
-        object.__setattr__(self, "_hubble_distance", (const.c / self._H0).to(u.Mpc))
-        # H0 in s^-1
-        H0_s = self._H0.value * _H0units_to_invs
-        # Hubble time
-        object.__setattr__(self, "_hubble_time", (_sec_to_Gyr / H0_s) << u.Gyr)
-
-        # Critical density at z=0 (grams per cubic cm)
-        cd0value = _critdens_const * H0_s**2
-        object.__setattr__(self, "_critical_density0", cd0value << u.g / u.cm**3)
-
-        # Compute photon density from Tcmb
-        object.__setattr__(
-            self,
-            "_Ogamma0",
-            _a_B_c2 * self._Tcmb0.value**4 / self._critical_density0.value,
-        )
-
-        # Compute Neutrino temperature:
-        # The constant in front is (4/11)^1/3 -- see any cosmology book for an
-        # explanation -- for example, Weinberg 'Cosmology' p 154 eq (3.1.21).
-        object.__setattr__(self, "_Tnu0", 0.7137658555036082 * self._Tcmb0)
-
         # Compute neutrino parameters:
         if self._m_nu is None:
             nneutrinos = 0
@@ -280,26 +250,12 @@ def __post_init__(self):
         # to do integrals with (perhaps surprisingly! But small python lists
         # are more efficient than small NumPy arrays).
         if self._massivenu:  # (`_massivenu` set in `m_nu`)
-            nu_y = (self._massivenu_mass / (_kB_evK * self._Tnu0)).value
+            nu_y = (self._massivenu_mass / (_kB_evK * self.Tnu0)).value
             nu_y_list = nu_y.tolist()
-            object.__setattr__(self, "_nu_y", nu_y)
-            object.__setattr__(self, "_nu_y_list", nu_y_list)
-            Onu0 = self._Ogamma0 * self.nu_relative_density(0)
         else:
-            # This case is particularly simple, so do it directly The 0.2271...
-            # is 7/8 (4/11)^(4/3) -- the temperature bit ^4 (blackbody energy
-            # density) times 7/8 for FD vs. BE statistics.
-            Onu0 = 0.22710731766 * self._Neff * self._Ogamma0
             nu_y = nu_y_list = None
-            object.__setattr__(self, "_nu_y", nu_y)
-            object.__setattr__(self, "_nu_y_list", nu_y_list)
-
-        object.__setattr__(self, "_Onu0", Onu0)
-
-        # Compute curvature density
-        object.__setattr__(
-            self, "_Ok0", 1.0 - self._Om0 - self._Ode0 - self._Ogamma0 - self._Onu0
-        )
+        object.__setattr__(self, "_nu_y", nu_y)
+        object.__setattr__(self, "_nu_y_list", nu_y_list)
 
         # Subclasses should override this reference if they provide
         #  more efficient scalar versions of inv_efunc.
@@ -358,64 +314,75 @@ def m_nu(self, param, value):
     @property
     def is_flat(self):
         """Return bool; `True` if the cosmology is flat."""
-        return bool((self._Ok0 == 0.0) and (self.Otot0 == 1.0))
+        return bool((self.Ok0 == 0.0) and (self.Otot0 == 1.0))
 
     @property
     def Otot0(self):
         """Omega total; the total density/critical density at z=0."""
-        return self._Om0 + self._Ogamma0 + self._Onu0 + self._Ode0 + self._Ok0
+        return self._Om0 + self.Ogamma0 + self.Onu0 + self._Ode0 + self.Ok0
 
-    @property
+    @cached_on_dict_property
     def Odm0(self):
         """Omega dark matter; dark matter density/critical density at z=0."""
-        return self._Odm0
+        return None if self.Ob0 is None else (self.Om0 - self.Ob0)
 
-    @property
+    @cached_on_dict_property
     def Ok0(self):
         """Omega curvature; the effective curvature density/critical density at z=0."""
-        return self._Ok0
+        return 1.0 - self.Om0 - self.Ode0 - self.Ogamma0 - self.Onu0
 
-    @property
+    @cached_on_dict_property
     def Tnu0(self):
         """Temperature of the neutrino background as |Quantity| at z=0."""
-        return self._Tnu0
+        # The constant in front is (4/11)^1/3 -- see any cosmology book for an
+        # explanation -- for example, Weinberg 'Cosmology' p 154 eq (3.1.21).
+        return 0.7137658555036082 * self.Tcmb0
 
     @property
     def has_massive_nu(self):
         """Does this cosmology have at least one massive neutrino species?"""
-        if self._Tnu0.value == 0:
+        if self.Tnu0.value == 0:
             return False
         return self._massivenu
 
-    @property
+    @cached_on_dict_property
     def h(self):
         """Dimensionless Hubble constant: h = H_0 / 100 [km/sec/Mpc]."""
-        return self._h
+        return self.H0.value / 100.0
 
-    @property
+    @cached_on_dict_property
     def hubble_time(self):
         """Hubble time as `~astropy.units.Quantity`."""
-        return self._hubble_time
+        return (_sec_to_Gyr / (self.H0.value * _H0units_to_invs)) << u.Gyr
 
-    @property
+    @cached_on_dict_property
     def hubble_distance(self):
         """Hubble distance as `~astropy.units.Quantity`."""
-        return self._hubble_distance
+        return (const.c / self.H0).to(u.Mpc)
 
-    @property
+    @cached_on_dict_property
     def critical_density0(self):
         """Critical density as `~astropy.units.Quantity` at z=0."""
-        return self._critical_density0
+        return (
+            _critdens_const * (self.H0.value * _H0units_to_invs) ** 2
+        ) << u.g / u.cm**3
 
-    @property
+    @cached_on_dict_property
     def Ogamma0(self):
         """Omega gamma; the density/critical density of photons at z=0."""
-        return self._Ogamma0
+        # photon density from Tcmb
+        return _a_B_c2 * self.Tcmb0.value**4 / self.critical_density0.value
 
-    @property
+    @cached_on_dict_property
     def Onu0(self):
         """Omega nu; the density/critical density of neutrinos at z=0."""
-        return self._Onu0
+        if self._massivenu:  # (`_massivenu` set in `m_nu`)
+            return self.Ogamma0 * self.nu_relative_density(0)
+        else:
+            # This case is particularly simple, so do it directly The 0.2271...
+            # is 7/8 (4/11)^(4/3) -- the temperature bit ^4 (blackbody energy
+            # density) times 7/8 for FD vs. BE statistics.
+            return 0.22710731766 * self.Neff * self.Ogamma0
 
     # ---------------------------------------------------------------
 
@@ -554,13 +521,13 @@ def Odm(self, z):
         This does not include neutrinos, even if non-relativistic at the
         redshift of interest.
         """
-        if self._Odm0 is None:
+        if self.Odm0 is None:
             raise ValueError(
                 "Baryonic density not set for this cosmology, "
                 "unclear meaning of dark matter density"
             )
         z = aszarr(z)
-        return self._Odm0 * (z + 1.0) ** 3 * self.inv_efunc(z) ** 2
+        return self.Odm0 * (z + 1.0) ** 3 * self.inv_efunc(z) ** 2
 
     @deprecated_keywords("z", since="7.0")
     def Ok(self, z):
@@ -581,9 +548,9 @@ def Ok(self, z):
             Returns `float` if the input is scalar.
         """
         z = aszarr(z)
-        if self._Ok0 == 0:  # Common enough to be worth checking explicitly
+        if self.Ok0 == 0:  # Common enough to be worth checking explicitly
             return np.zeros(z.shape) if hasattr(z, "shape") else 0.0
-        return self._Ok0 * (z + 1.0) ** 2 * self.inv_efunc(z) ** 2
+        return self.Ok0 * (z + 1.0) ** 2 * self.inv_efunc(z) ** 2
 
     @deprecated_keywords("z", since="7.0")
     def Ode(self, z):
@@ -629,7 +596,7 @@ def Ogamma(self, z):
             Returns `float` if the input is scalar.
         """
         z = aszarr(z)
-        return self._Ogamma0 * (z + 1.0) ** 4 * self.inv_efunc(z) ** 2
+        return self.Ogamma0 * (z + 1.0) ** 4 * self.inv_efunc(z) ** 2
 
     @deprecated_keywords("z", since="7.0")
     def Onu(self, z):
@@ -654,7 +621,7 @@ def Onu(self, z):
             Returns `float` if the input is scalar.
         """
         z = aszarr(z)
-        if self._Onu0 == 0:  # Common enough to be worth checking explicitly
+        if self.Onu0 == 0:  # Common enough to be worth checking explicitly
             return np.zeros(z.shape) if hasattr(z, "shape") else 0.0
         return self.Ogamma(z) * self.nu_relative_density(z)
 
@@ -694,7 +661,7 @@ def Tnu(self, z):
         Tnu : `~astropy.units.Quantity` ['temperature']
             The temperature of the cosmic neutrino background in K.
         """
-        return self._Tnu0 * (aszarr(z) + 1.0)
+        return self.Tnu0 * (aszarr(z) + 1.0)
 
     @deprecated_keywords("z", since="7.0")
     def nu_relative_density(self, z):
@@ -866,15 +833,15 @@ def efunc(self, z):
         It is not necessary to override this method, but if de_density_scale
         takes a particularly simple form, it may be advantageous to.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
         return np.sqrt(
-            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0)
+            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0)
             + self._Ode0 * self.de_density_scale(z)
         )
 
@@ -897,15 +864,15 @@ def inv_efunc(self, z):
             Returns `float` if the input is scalar.
         """
         # Avoid the function overhead by repeating code
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
         return (
-            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0)
+            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0)
             + self._Ode0 * self.de_density_scale(z)
         ) ** (-0.5)
 
@@ -1067,7 +1034,7 @@ def _lookback_time(self, z, /):
         t : `~astropy.units.Quantity` ['time']
             Lookback time in Gyr to each input redshift.
         """
-        return self._hubble_time * self._integral_lookback_time(z)
+        return self.hubble_time * self._integral_lookback_time(z)
 
     @vectorize_redshift_method
     def _integral_lookback_time(self, z, /):
@@ -1156,7 +1123,7 @@ def _age(self, z, /):
         t : `~astropy.units.Quantity` ['time']
             The age of the universe in Gyr at each input redshift.
         """
-        return self._hubble_time * self._integral_age(z)
+        return self.hubble_time * self._integral_age(z)
 
     @vectorize_redshift_method
     def _integral_age(self, z, /):
@@ -1198,7 +1165,7 @@ def critical_density(self, z):
         rho : `~astropy.units.Quantity`
             Critical density in g/cm^3 at each input redshift.
         """
-        return self._critical_density0 * (self.efunc(z)) ** 2
+        return self.critical_density0 * (self.efunc(z)) ** 2
 
     @deprecated_keywords("z", since="7.0")
     def comoving_distance(self, z):
@@ -1285,7 +1252,7 @@ def _integral_comoving_distance_z1z2(self, z1, z2, /):
         d : `~astropy.units.Quantity` ['length']
             Comoving distance in Mpc between each input redshift.
         """
-        return self._hubble_distance * self._integral_comoving_distance_z1z2_scalar(z1, z2)  # fmt: skip
+        return self.hubble_distance * self._integral_comoving_distance_z1z2_scalar(z1, z2)  # fmt: skip
 
     @deprecated_keywords("z", since="7.0")
     def comoving_transverse_distance(self, z):
@@ -1340,12 +1307,12 @@ def _comoving_transverse_distance_z1z2(self, z1, z2, /):
         -----
         This quantity is also called the 'proper motion distance' in some texts.
         """
-        Ok0 = self._Ok0
+        Ok0 = self.Ok0
         dc = self._comoving_distance_z1z2(z1, z2)
         if Ok0 == 0:
             return dc
         sqrtOk0 = sqrt(abs(Ok0))
-        dh = self._hubble_distance
+        dh = self.hubble_distance
         if Ok0 > 0:
             return dh / sqrtOk0 * np.sinh(sqrtOk0 * dc.value / dh.value)
         else:
@@ -1518,11 +1485,11 @@ def comoving_volume(self, z):
         V : `~astropy.units.Quantity`
             Comoving volume in :math:`Mpc^3` at each input redshift.
         """
-        Ok0 = self._Ok0
+        Ok0 = self.Ok0
         if Ok0 == 0:
             return 4.0 / 3.0 * pi * self.comoving_distance(z) ** 3
 
-        dh = self._hubble_distance.value  # .value for speed
+        dh = self.hubble_distance.value  # .value for speed
         dm = self.comoving_transverse_distance(z).value
         term1 = 4.0 * pi * dh**3 / (2.0 * Ok0) * u.Mpc**3
         term2 = dm / dh * np.sqrt(1 + Ok0 * (dm / dh) ** 2)
@@ -1558,7 +1525,7 @@ def differential_comoving_volume(self, z):
             input redshift.
         """
         dm = self.comoving_transverse_distance(z)
-        return self._hubble_distance * (dm**2.0) / (self.efunc(z) << u.steradian)
+        return self.hubble_distance * (dm**2.0) / (self.efunc(z) << u.steradian)
 
     @deprecated_keywords("z", since="7.0")
     def kpc_comoving_per_arcmin(self, z):
@@ -1678,9 +1645,9 @@ def __post_init__(self):
         object.__setattr__(self, "_Ode0", 0)
         super().__post_init__()
         # Do some twiddling after the fact to get flatness
-        object.__setattr__(self, "_Ok0", 0.0)
-        object.__setattr__(
-            self, "_Ode0", 1.0 - (self._Om0 + self._Ogamma0 + self._Onu0 + self._Ok0)
+        self.__dict__["Ok0"] = 0.0
+        object.__setattr__(  # managed by a Property
+            self, "_Ode0", 1.0 - (self.Om0 + self.Ogamma0 + self.Onu0 + self.Ok0)
         )
 
     @lazyproperty
@@ -1693,8 +1660,9 @@ def nonflat(self: _FlatFLRWMixinT) -> _FLRWT:
         # Make new instance, respecting args vs kwargs
         inst = self.__nonflatclass__(*ba.args, **ba.kwargs)
         # Because of machine precision, make sure parameters exactly match
-        for n in (*inst._parameters_all, "Ok0"):
+        for n in tuple(inst._parameters_all):
             object.__setattr__(inst, "_" + n, getattr(self, n))
+        inst.__dict__["Ok0"] = self.Ok0
 
         return inst
 
diff --git a/astropy/cosmology/flrw/lambdacdm.py b/astropy/cosmology/flrw/lambdacdm.py
index 61a07ba0368..e7eb66db35c 100644
--- a/astropy/cosmology/flrw/lambdacdm.py
+++ b/astropy/cosmology/flrw/lambdacdm.py
@@ -95,8 +95,8 @@ def __post_init__(self):
         # about what is being done here.
         if self._Tcmb0.value == 0:
             inv_efunc_scalar = scalar_inv_efuncs.lcdm_inv_efunc_norel
-            inv_efunc_scalar_args = (self._Om0, self._Ode0, self._Ok0)
-            if self._Ok0 != 0:
+            inv_efunc_scalar_args = (self._Om0, self._Ode0, self.Ok0)
+            if self.Ok0 != 0:
                 object.__setattr__(
                     self,
                     "_comoving_distance_z1z2",
@@ -107,16 +107,16 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0 + self._Onu0,
+                self.Ok0,
+                self.Ogamma0 + self.Onu0,
             )
         else:
             inv_efunc_scalar = scalar_inv_efuncs.lcdm_inv_efunc
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0,
+                self.Ok0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -124,7 +124,7 @@ def __post_init__(self):
         object.__setattr__(self, "_inv_efunc_scalar", inv_efunc_scalar)
         object.__setattr__(self, "_inv_efunc_scalar_args", inv_efunc_scalar_args)
 
-        if self._Tcmb0.value == 0 and self._Ok0 == 0:
+        if self._Tcmb0.value == 0 and self.Ok0 == 0:
             self._optimize_flat_norad()
 
     def _optimize_flat_norad(self):
@@ -241,10 +241,10 @@ def _elliptic_comoving_distance_z1z2(self, z1, z2, /):
 
         # The analytic solution is not valid for any of Om0, Ode0, Ok0 == 0.
         # Use the explicit integral solution for these cases.
-        if self._Om0 == 0 or self._Ode0 == 0 or self._Ok0 == 0:
+        if self._Om0 == 0 or self._Ode0 == 0 or self.Ok0 == 0:
             return self._integral_comoving_distance_z1z2(z1, z2)
 
-        b = -(27.0 / 2) * self._Om0**2 * self._Ode0 / self._Ok0**3
+        b = -(27.0 / 2) * self._Om0**2 * self._Ode0 / self.Ok0**3
         kappa = b / abs(b)
         if (b < 0) or (2 < b):
 
@@ -260,8 +260,8 @@ def phi_z(Om0, Ok0, kappa, y1, A, z, /):
             g = 1 / sqrt(A)
             k2 = (2 * A + kappa * (1 + 3 * y1)) / (4 * A)
 
-            phi_z1 = phi_z(self._Om0, self._Ok0, kappa, y1, A, z1)
-            phi_z2 = phi_z(self._Om0, self._Ok0, kappa, y1, A, z2)
+            phi_z1 = phi_z(self._Om0, self.Ok0, kappa, y1, A, z1)
+            phi_z2 = phi_z(self._Om0, self.Ok0, kappa, y1, A, z2)
         # Get lower-right 0<b<2 solution in Om0, Ode0 plane.
         # Fot the upper-left 0<b<2 solution the Big Bang didn't happen.
         elif (0 < b) and (b < 2) and self._Om0 > self._Ode0:
@@ -276,12 +276,12 @@ def phi_z(Om0, Ok0, y1, y2, z, /):
             y3 = (1.0 / 3) * (-1 + yb - yc)
             g = 2 / sqrt(y1 - y2)
             k2 = (y1 - y3) / (y1 - y2)
-            phi_z1 = phi_z(self._Om0, self._Ok0, y1, y2, z1)
-            phi_z2 = phi_z(self._Om0, self._Ok0, y1, y2, z2)
+            phi_z1 = phi_z(self._Om0, self.Ok0, y1, y2, z1)
+            phi_z2 = phi_z(self._Om0, self.Ok0, y1, y2, z2)
         else:
             return self._integral_comoving_distance_z1z2(z1, z2)
 
-        prefactor = self._hubble_distance / sqrt(abs(self._Ok0))
+        prefactor = self.hubble_distance / sqrt(abs(self.Ok0))
         return prefactor * g * (ellipkinc(phi_z1, k2) - ellipkinc(phi_z2, k2))
 
     def _dS_comoving_distance_z1z2(self, z1, z2, /):
@@ -314,7 +314,7 @@ def _dS_comoving_distance_z1z2(self, z1, z2, /):
         except ValueError as e:
             raise ValueError("z1 and z2 have different shapes") from e
 
-        return self._hubble_distance * (z2 - z1)
+        return self.hubble_distance * (z2 - z1)
 
     def _EdS_comoving_distance_z1z2(self, z1, z2, /):
         r"""Einstein-de Sitter comoving LoS distance in Mpc between two redshifts.
@@ -347,7 +347,7 @@ def _EdS_comoving_distance_z1z2(self, z1, z2, /):
         except ValueError as e:
             raise ValueError("z1 and z2 have different shapes") from e
 
-        prefactor = 2 * self._hubble_distance
+        prefactor = 2 * self.hubble_distance
         return prefactor * ((z1 + 1.0) ** (-1.0 / 2) - (z2 + 1.0) ** (-1.0 / 2))
 
     def _hypergeometric_comoving_distance_z1z2(self, z1, z2, /):
@@ -388,7 +388,7 @@ def _hypergeometric_comoving_distance_z1z2(self, z1, z2, /):
 
         s = ((1 - self._Om0) / self._Om0) ** (1.0 / 3)
         # Use np.sqrt here to handle negative s (Om0>1).
-        prefactor = self._hubble_distance / np.sqrt(s * self._Om0)
+        prefactor = self.hubble_distance / np.sqrt(s * self._Om0)
         return prefactor * (
             self._T_hypergeometric(s / (z1 + 1.0))
             - self._T_hypergeometric(s / (z2 + 1.0))
@@ -435,7 +435,7 @@ def _dS_age(self, z, /):
             The age of the universe in Gyr at each input redshift.
         """
         t = inf if isinstance(z, Number) else np.full_like(z, inf, dtype=float)
-        return self._hubble_time * t
+        return self.hubble_time * t
 
     def _EdS_age(self, z, /):
         r"""Age of the universe in Gyr at redshift ``z``.
@@ -461,7 +461,7 @@ def _EdS_age(self, z, /):
         .. [1] Thomas, R., & Kantowski, R. (2000). Age-redshift relation for
                standard cosmology. PRD, 62(10), 103507.
         """
-        return (2.0 / 3) * self._hubble_time * (aszarr(z) + 1.0) ** (-1.5)
+        return (2.0 / 3) * self.hubble_time * (aszarr(z) + 1.0) ** (-1.5)
 
     def _flat_age(self, z, /):
         r"""Age of the universe in Gyr at redshift ``z``.
@@ -489,7 +489,7 @@ def _flat_age(self, z, /):
         """
         # Use np.sqrt, np.arcsinh instead of math.sqrt, math.asinh
         # to handle properly the complex numbers for 1 - Om0 < 0
-        prefactor = (2.0 / 3) * self._hubble_time / np.emath.sqrt(1 - self._Om0)
+        prefactor = (2.0 / 3) * self.hubble_time / np.emath.sqrt(1 - self._Om0)
         arg = np.arcsinh(
             np.emath.sqrt((1 / self._Om0 - 1 + 0j) / (aszarr(z) + 1.0) ** 3)
         )
@@ -548,7 +548,7 @@ def _dS_lookback_time(self, z, /):
         t : `~astropy.units.Quantity` ['time']
             Lookback time in Gyr to each input redshift.
         """
-        return self._hubble_time * log(aszarr(z) + 1.0)
+        return self.hubble_time * log(aszarr(z) + 1.0)
 
     def _flat_lookback_time(self, z, /):
         r"""Lookback time in Gyr to redshift ``z``.
@@ -596,14 +596,14 @@ def efunc(self, z):
         """
         # We override this because it takes a particularly simple
         # form for a cosmological constant
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
-        return np.sqrt(zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0) + self._Ode0)
+        return np.sqrt(zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0) + self._Ode0)
 
     @deprecated_keywords("z", since="7.0")
     def inv_efunc(self, z):
@@ -624,14 +624,14 @@ def inv_efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H_z = H_0 / E`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
-        return (zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0) + self._Ode0) ** (
+        return (zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0) + self._Ode0) ** (
             -0.5
         )
 
@@ -706,21 +706,21 @@ def __post_init__(self):
             inv_efunc_scalar_args = (self._Om0, self._Ode0)
             # Repeat the optimization reassignments here because the init
             # of the LambaCDM above didn't actually create a flat cosmology.
-            # That was done through the explicit tweak setting self._Ok0.
+            # That was done through the explicit tweak setting self.Ok0.
             self._optimize_flat_norad()
         elif not self._massivenu:
             inv_efunc_scalar = scalar_inv_efuncs.flcdm_inv_efunc_nomnu
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0 + self._Onu0,
+                self.Ogamma0 + self.Onu0,
             )
         else:
             inv_efunc_scalar = scalar_inv_efuncs.flcdm_inv_efunc
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -749,10 +749,10 @@ def efunc(self, z):
         """
         # We override this because it takes a particularly simple
         # form for a cosmological constant
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
@@ -777,10 +777,10 @@ def inv_efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H_z = H_0 / E`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
         return (zp1**3 * (Or * zp1 + self._Om0) + self._Ode0) ** (-0.5)
diff --git a/astropy/cosmology/flrw/tests/test_base.py b/astropy/cosmology/flrw/tests/test_base.py
index 71ab2882cf7..41182fe680a 100644
--- a/astropy/cosmology/flrw/tests/test_base.py
+++ b/astropy/cosmology/flrw/tests/test_base.py
@@ -18,6 +18,7 @@
     Parameter,
     Planck18,
 )
+from astropy.cosmology._utils import CachedInDictPropertyDescriptor
 from astropy.cosmology.core import _COSMOLOGY_CLASSES, dataclass_decorator
 from astropy.cosmology.flrw.base import _a_B_c2, _critdens_const, _H0units_to_invs, quad
 from astropy.cosmology.parameter._core import MISSING
@@ -554,13 +555,11 @@ def test_init_Tcmb0_zeroing(self, cosmo_cls, ba):
     # Properties
 
     def test_Odm0(self, cosmo_cls, cosmo):
-        """Test property ``Odm0``."""
+        """Test CachedInDictPropertyDescriptor ``Odm0``."""
         # on the class
-        assert isinstance(cosmo_cls.Odm0, property)
-        assert cosmo_cls.Odm0.fset is None  # immutable
+        assert isinstance(cosmo_cls.Odm0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.Odm0 is cosmo._Odm0
         # Odm0 can be None, if Ob0 is None. Otherwise DM = matter - baryons.
         if cosmo.Ob0 is None:
             assert cosmo.Odm0 is None
@@ -568,13 +567,11 @@ def test_Odm0(self, cosmo_cls, cosmo):
             assert np.allclose(cosmo.Odm0, cosmo.Om0 - cosmo.Ob0)
 
     def test_Ok0(self, cosmo_cls, cosmo):
-        """Test property ``Ok0``."""
+        """Test CachedInDictPropertyDescriptor ``Ok0``."""
         # on the class
-        assert isinstance(cosmo_cls.Ok0, property)
-        assert cosmo_cls.Ok0.fset is None  # immutable
+        assert isinstance(cosmo_cls.Ok0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.Ok0 is cosmo._Ok0
         assert np.allclose(
             cosmo.Ok0, 1.0 - (cosmo.Om0 + cosmo.Ode0 + cosmo.Ogamma0 + cosmo.Onu0)
         )
@@ -590,13 +587,11 @@ def test_is_flat(self, cosmo_cls, cosmo):
         assert cosmo.is_flat is bool((cosmo.Ok0 == 0.0) and (cosmo.Otot0 == 1.0))
 
     def test_Tnu0(self, cosmo_cls, cosmo):
-        """Test property ``Tnu0``."""
+        """Test CachedInDictPropertyDescriptor ``Tnu0``."""
         # on the class
-        assert isinstance(cosmo_cls.Tnu0, property)
-        assert cosmo_cls.Tnu0.fset is None  # immutable
+        assert isinstance(cosmo_cls.Tnu0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.Tnu0 is cosmo._Tnu0
         assert cosmo.Tnu0.unit == u.K
         assert u.allclose(cosmo.Tnu0, 0.7137658555036082 * cosmo.Tcmb0, rtol=1e-5)
 
@@ -615,54 +610,44 @@ def test_has_massive_nu(self, cosmo_cls, cosmo):
     def test_h(self, cosmo_cls, cosmo):
         """Test property ``h``."""
         # on the class
-        assert isinstance(cosmo_cls.h, property)
-        assert cosmo_cls.h.fset is None  # immutable
+        assert isinstance(cosmo_cls.h, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.h is cosmo._h
         assert np.allclose(cosmo.h, cosmo.H0.value / 100.0)
 
     def test_hubble_time(self, cosmo_cls, cosmo):
         """Test property ``hubble_time``."""
         # on the class
-        assert isinstance(cosmo_cls.hubble_time, property)
-        assert cosmo_cls.hubble_time.fset is None  # immutable
+        assert isinstance(cosmo_cls.hubble_time, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.hubble_time is cosmo._hubble_time
         assert u.allclose(cosmo.hubble_time, (1 / cosmo.H0) << u.Gyr)
 
     def test_hubble_distance(self, cosmo_cls, cosmo):
-        """Test property ``hubble_distance``."""
+        """Test CachedInDictPropertyDescriptor ``hubble_distance``."""
         # on the class
-        assert isinstance(cosmo_cls.hubble_distance, property)
-        assert cosmo_cls.hubble_distance.fset is None  # immutable
+        assert isinstance(cosmo_cls.hubble_distance, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.hubble_distance is cosmo._hubble_distance
         assert cosmo.hubble_distance == (const.c / cosmo._H0).to(u.Mpc)
 
     def test_critical_density0(self, cosmo_cls, cosmo):
-        """Test property ``critical_density0``."""
+        """Test CachedInDictPropertyDescriptor ``critical_density0``."""
         # on the class
-        assert isinstance(cosmo_cls.critical_density0, property)
-        assert cosmo_cls.critical_density0.fset is None  # immutable
+        assert isinstance(cosmo_cls.critical_density0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.critical_density0 is cosmo._critical_density0
         assert cosmo.critical_density0.unit == u.g / u.cm**3
 
         cd0value = _critdens_const * (cosmo.H0.value * _H0units_to_invs) ** 2
         assert cosmo.critical_density0.value == cd0value
 
     def test_Ogamma0(self, cosmo_cls, cosmo):
-        """Test property ``Ogamma0``."""
+        """Test CachedInDictPropertyDescriptor ``Ogamma0``."""
         # on the class
-        assert isinstance(cosmo_cls.Ogamma0, property)
-        assert cosmo_cls.Ogamma0.fset is None  # immutable
+        assert isinstance(cosmo_cls.Ogamma0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.Ogamma0 is cosmo._Ogamma0
         # Ogamma cor \propto T^4/rhocrit
         expect = _a_B_c2 * cosmo.Tcmb0.value**4 / cosmo.critical_density0.value
         assert np.allclose(cosmo.Ogamma0, expect)
@@ -671,13 +656,11 @@ def test_Ogamma0(self, cosmo_cls, cosmo):
             assert cosmo.Ogamma0 == 0
 
     def test_Onu0(self, cosmo_cls, cosmo):
-        """Test property ``Onu0``."""
+        """Test CachedInDictPropertyDescriptor ``Onu0``."""
         # on the class
-        assert isinstance(cosmo_cls.Onu0, property)
-        assert cosmo_cls.Onu0.fset is None  # immutable
+        assert isinstance(cosmo_cls.Onu0, CachedInDictPropertyDescriptor)
 
         # on the instance
-        assert cosmo.Onu0 is cosmo._Onu0
         # neutrino temperature <= photon temperature since the neutrinos
         # decouple first.
         if cosmo.has_massive_nu:  # Tcmb0 > 0 & has massive
@@ -1044,10 +1027,8 @@ def test_init(self, cosmo_cls):
         super().test_init(cosmo_cls)
 
         cosmo = cosmo_cls(*self.cls_args, **self.cls_kwargs)
-        assert cosmo._Ok0 == 0.0
-        assert cosmo._Ode0 == 1.0 - (
-            cosmo._Om0 + cosmo._Ogamma0 + cosmo._Onu0 + cosmo._Ok0
-        )
+        assert cosmo.Ok0 == 0.0
+        assert cosmo.Ode0 == 1.0 - (cosmo.Om0 + cosmo.Ogamma0 + cosmo.Onu0 + cosmo.Ok0)
 
     def test_Ok0(self, cosmo_cls, cosmo):
         """Test property ``Ok0``."""
@@ -1125,7 +1106,7 @@ def test_is_equivalent(self, cosmo, nonflatcosmo):
             Ode0=1.0 - cosmo.Om0 - cosmo.Ogamma0 - cosmo.Onu0,
             **self.cls_kwargs,
         )
-        object.__setattr__(flat, "_Ok0", 0.0)
+        flat.__dict__["Ok0"] = 0.0
         assert flat.is_equivalent(cosmo)
         assert cosmo.is_equivalent(flat)
 
diff --git a/astropy/cosmology/flrw/w0cdm.py b/astropy/cosmology/flrw/w0cdm.py
index 3106de640a2..66d36b8542c 100644
--- a/astropy/cosmology/flrw/w0cdm.py
+++ b/astropy/cosmology/flrw/w0cdm.py
@@ -91,14 +91,14 @@ def __post_init__(self):
         # about what is being done here.
         if self._Tcmb0.value == 0:
             inv_efunc_scalar = scalar_inv_efuncs.wcdm_inv_efunc_norel
-            inv_efunc_scalar_args = (self._Om0, self._Ode0, self._Ok0, self._w0)
+            inv_efunc_scalar_args = (self._Om0, self._Ode0, self.Ok0, self._w0)
         elif not self._massivenu:
             inv_efunc_scalar = scalar_inv_efuncs.wcdm_inv_efunc_nomnu
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0 + self._Onu0,
+                self.Ok0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
             )
         else:
@@ -106,8 +106,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0,
+                self.Ok0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -190,15 +190,15 @@ def efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H(z) = H_0 E(z)`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
         return sqrt(
-            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0)
+            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0)
             + self._Ode0 * zp1 ** (3.0 * (1.0 + self._w0))
         )
 
@@ -221,15 +221,15 @@ def inv_efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H_z = H_0 / E`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
         return (
-            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self._Ok0)
+            zp1**2 * ((Or * zp1 + self._Om0) * zp1 + self.Ok0)
             + self._Ode0 * zp1 ** (3.0 * (1.0 + self._w0))
         ) ** (-0.5)
 
@@ -312,7 +312,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0 + self._Onu0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
             )
         else:
@@ -320,7 +320,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -348,10 +348,10 @@ def efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H(z) = H_0 E(z)`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
@@ -378,10 +378,10 @@ def inv_efunc(self, z):
             Returns `float` if the input is scalar.
             Defined such that :math:`H(z) = H_0 E(z)`.
         """
-        Or = self._Ogamma0 + (
-            self._Onu0
+        Or = self.Ogamma0 + (
+            self.Onu0
             if not self._massivenu
-            else self._Ogamma0 * self.nu_relative_density(z)
+            else self.Ogamma0 * self.nu_relative_density(z)
         )
         zp1 = aszarr(z) + 1.0  # (converts z [unit] -> z [dimensionless])
 
diff --git a/astropy/cosmology/flrw/w0wacdm.py b/astropy/cosmology/flrw/w0wacdm.py
index b84852e7b64..4b38a10f49f 100644
--- a/astropy/cosmology/flrw/w0wacdm.py
+++ b/astropy/cosmology/flrw/w0wacdm.py
@@ -111,7 +111,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
+                self.Ok0,
                 self._w0,
                 self._wa,
             )
@@ -120,8 +120,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0 + self._Onu0,
+                self.Ok0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
                 self._wa,
             )
@@ -130,8 +130,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0,
+                self.Ok0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -294,7 +294,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0 + self._Onu0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
                 self._wa,
             )
@@ -303,7 +303,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
diff --git a/astropy/cosmology/flrw/w0wzcdm.py b/astropy/cosmology/flrw/w0wzcdm.py
index 4d020bbecfe..7126946215b 100644
--- a/astropy/cosmology/flrw/w0wzcdm.py
+++ b/astropy/cosmology/flrw/w0wzcdm.py
@@ -104,7 +104,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
+                self.Ok0,
                 self._w0,
                 self._wz,
             )
@@ -113,8 +113,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0 + self._Onu0,
+                self.Ok0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
                 self._wz,
             )
@@ -123,8 +123,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0,
+                self.Ok0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -278,7 +278,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0 + self._Onu0,
+                self.Ogamma0 + self.Onu0,
                 self._w0,
                 self._wz,
             )
@@ -287,7 +287,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
diff --git a/astropy/cosmology/flrw/wpwazpcdm.py b/astropy/cosmology/flrw/wpwazpcdm.py
index 96fc1e26ddb..ca6eace25f1 100644
--- a/astropy/cosmology/flrw/wpwazpcdm.py
+++ b/astropy/cosmology/flrw/wpwazpcdm.py
@@ -128,7 +128,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
+                self.Ok0,
                 self._wp,
                 apiv,
                 self._wa,
@@ -138,8 +138,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0 + self._Onu0,
+                self.Ok0,
+                self.Ogamma0 + self.Onu0,
                 self._wp,
                 apiv,
                 self._wa,
@@ -149,8 +149,8 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ok0,
-                self._Ogamma0,
+                self.Ok0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,
@@ -328,7 +328,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0 + self._Onu0,
+                self.Ogamma0 + self.Onu0,
                 self._wp,
                 apiv,
                 self._wa,
@@ -338,7 +338,7 @@ def __post_init__(self):
             inv_efunc_scalar_args = (
                 self._Om0,
                 self._Ode0,
-                self._Ogamma0,
+                self.Ogamma0,
                 self._neff_per_nu,
                 self._nmasslessnu,
                 self._nu_y_list,