diff --git a/source/mass_distributions.spherical.beta_profile.F90 b/source/mass_distributions.spherical.beta_profile.F90
index a822de87e4..f0bb7909d0 100644
--- a/source/mass_distributions.spherical.beta_profile.F90
+++ b/source/mass_distributions.spherical.beta_profile.F90
@@ -463,9 +463,9 @@ double precision function betaProfilePotential(self,coordinates,componentType,ma
                &  )
        end if
     end if
-    if (.not.self%isDimensionless())                  &
-         & betaProfilePotential=  &
-         &   betaProfilePotential &
+    if (.not.self%isDimensionless()) &
+         & betaProfilePotential=     &
+         &   betaProfilePotential    &
          &   *gravitationalConstantGalacticus
     return
   end function betaProfilePotential
@@ -476,6 +476,8 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
     !!}
     use :: Hypergeometric_Functions, only : Hypergeometric_2F1
     use :: Numerical_Comparison    , only : Values_Agree
+    use :: Error                   , only : Error_Report
+    use :: Gamma_Functions         , only : Gamma_Function
     implicit none
     class           (massDistributionBetaProfile ), intent(inout)           :: self
     double precision                              , intent(in   )           :: moment
@@ -512,7 +514,7 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
                & +radialMomentTwoThirds(specialCaseMoment,fractionalRadiusMaximum)
        else
           ! General solution.
-          betaProfileDensityRadialMoment=                   &
+          betaProfileDensityRadialMoment=                                         &
                & +fractionalRadiusMaximum**(moment+1.0d0)                         &
                & *Hypergeometric_2F1     (                                        &
                &                          [(moment+1.0d0)/2.0d0,1.5d0*self%beta], &
@@ -522,12 +524,34 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
                & /                         (moment+1.0d0)
        end if
     else
-       betaProfileDensityRadialMoment=0.0d0
+       if (moment >= 3.0d0*self%beta-1.0d0) then
+          betaProfileDensityRadialMoment=0.0d0
+          if (present(isInfinite)) then
+             isInfinite=.true.
+             return
+          else
+             call Error_Report('radial moment is infinite at r=0 for m ≥ 3β-1'//{introspection:location})
+          end if
+       else
+          betaProfileDensityRadialMoment=+Gamma_Function((+3.0d0+moment                )/2.0d0) &
+               &                         *Gamma_Function((-1.0d0-moment+3.0d0*self%beta)/2.0d0) &
+               &                         /Gamma_Function((             +3.0d0*self%beta)/2.0d0) &
+               &                         /               (+1.0d0+moment                ) 
+       end if
     end if
     if (present(radiusMinimum)) then
        fractionalRadiusMinimum=radiusMinimum/self%coreRadius
     else
        fractionalRadiusMinimum=0.0d0
+       if (moment <= -1.0d0) then
+          betaProfileDensityRadialMoment=0.0d0
+          if (present(isInfinite)) then
+             isInfinite=.true.
+             return
+          else
+             call Error_Report('radial moment is infinite at r=0 for m ≤ -1'//{introspection:location})
+          end if
+       end if
     end if
     if (specialCaseMoment /= -huge(0)) then
        ! Special case for 0ᵗʰ, 1ˢᵗ, 2ⁿᵈ, and 3ʳᵈ moments of a β=2/3 distribution.
diff --git a/source/tests.mass_distributions.F90 b/source/tests.mass_distributions.F90
index 95ee862506..31ace3f5cf 100644
--- a/source/tests.mass_distributions.F90
+++ b/source/tests.mass_distributions.F90
@@ -137,6 +137,7 @@ program Test_Mass_Distributions
 
   ! Beta profile.
   call Unit_Tests_Begin_Group("Beta profile")
+  call Unit_Tests_Begin_Group("β = 2/3")
   allocate(massDistributionBetaProfile :: massDistribution_)
   select type (massDistribution_)
   type is (massDistributionBetaProfile)
@@ -160,19 +161,19 @@ program Test_Mass_Distributions
           &      absTol=1.0d-6                                                 &
           &     )
      call Assert(                                                              &
-          &      "Radial moment, m=1"                                        , &
+          &      "Radial moment, m=1, from 0 to 1"                           , &
           &      +massDistribution_%densityRadialMoment  (1.0d0,0.0d0,1.0d0) , &
           &      +0.3465735902799727d0                                       , &
           &      absTol=1.0d-6                                                 &
           &     )
      call Assert(                                                              &
-          &      "Radial moment, m=2"                                        , &
+          &      "Radial moment, m=2, from 0 to 1"                           , &
           &      +massDistribution_%densityRadialMoment  (2.0d0,0.0d0,1.0d0) , &
           &      +0.2146018366025517d0                                       , &
           &      absTol=1.0d-6                                                 &
           &     )
      call Assert(                                                              &
-          &      "Radial moment, m=3"                                        , &
+          &      "Radial moment, m=3, from 0 to 1"                           , &
           &      +massDistribution_%densityRadialMoment  (3.0d0,0.0d0,1.0d0) , &
           &      +0.1534264097200273d0                                       , &
           &      absTol=1.0d-6                                                 &
@@ -184,6 +185,36 @@ program Test_Mass_Distributions
           &      absTol=1.0d-6                                                 &
           &     )
   end select
+  call Unit_Tests_End_Group()
+  deallocate(massDistribution_)
+  call Unit_Tests_Begin_Group("β = 2")
+  allocate(massDistributionBetaProfile :: massDistribution_)
+  select type (massDistribution_)
+  type is (massDistributionBetaProfile)
+     massDistribution_=massDistributionBetaProfile(beta=2.0d0,dimensionless=.true.)
+  end select
+  select type (massDistribution_)
+  class is (massDistributionSpherical)
+     call Assert(                                                              &
+          &      "Radial moment, m=1, unbounded"                             , &
+          &      +massDistribution_%densityRadialMoment  (1.0d0            ) , &
+          &      +1.0d0/4.0d0                                                , &
+          &      absTol=1.0d-6                                                 &
+          &     )
+     call Assert(                                                              &
+          &      "Radial moment, m=2, unbounded"                             , &
+          &      +massDistribution_%densityRadialMoment  (2.0d0            ) , &
+          &      +Pi/16.0d0                                                  , &
+          &      absTol=1.0d-6                                                 &
+          &     )
+     call Assert(                                                              &
+          &      "Radial moment, m=3, unbounded"                             , &
+          &      +massDistribution_%densityRadialMoment  (3.0d0            ) , &
+          &      +1.0d0/4.0d0                                                , &
+          &      absTol=1.0d-6                                                 &
+          &     )
+  end select
+  call Unit_Tests_End_Group()
   ! Ensure that a dimensionful profile produces the correct mass inside of its outer radius.
   deallocate(massDistribution_)
   allocate(massDistributionBetaProfile :: massDistribution_)