diff --git a/CHANGES.rst b/CHANGES.rst
index 320134df..d429dcac 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -8,6 +8,7 @@ New Features
   computing basis function coefficients for the SCF potential from discrete particles.
 
 - Added the Burkert potential as a built-in cpotential.
+- Added a method to generate the Burkert potential with just r0 as an input
 
 Bug fixes
 ---------
diff --git a/gala/potential/potential/builtin/core.py b/gala/potential/potential/builtin/core.py
index 8889c420..3e56820f 100644
--- a/gala/potential/potential/builtin/core.py
+++ b/gala/potential/potential/builtin/core.py
@@ -407,6 +407,24 @@ class BurkertPotential(CPotentialBase):
 
     Wrapper = BurkertWrapper
 
+    
+    @classmethod
+    def from_r0(cls, r0, units=None):
+        r"""
+        from_r0(r0, units=None)
+
+        Initialize a Burkert potential from the core radius, ``r0``.
+        See Equations 4 and 5 of Mori & Burkert.
+
+        Parameters
+        ----------
+        r0 : :class:`~astropy.units.Quantity`, numeric [length]
+            The core radius of the Burkert potential.
+        """
+        a = 0.021572405792749372 * u.Msun / u.pc**3  # converted: 1.46e-24 g/cm**3
+        rho_d0 = a * (r0 / (3.07 * u.kpc))**(-2/3)
+        return cls(rho=rho_d0, r0=r0, units=units)
+
 
 # ============================================================================
 # Flattened, axisymmetric models
diff --git a/gala/potential/potential/tests/test_all_builtin.py b/gala/potential/potential/tests/test_all_builtin.py
index 9ed53e32..8649fe8f 100644
--- a/gala/potential/potential/tests/test_all_builtin.py
+++ b/gala/potential/potential/tests/test_all_builtin.py
@@ -553,3 +553,14 @@ def test_against_sympy(self):
     @pytest.mark.skip(reason="Hessian not implemented for Burkert potential")
     def test_hessian(self):
         pass
+
+    def test_from_r0(self):
+        # Test against values from Zhu+2023
+        pot = p.BurkertPotential.from_r0(r0=11.87 * u.kpc, units=galactic)
+
+        rho = pot.parameters['rho'].to(u.g / u.cm ** 3)
+        rho_check = 5.93e-25 * u.g / u.cm ** 3
+
+        # Check a 1% tolerance on inferred density against published values
+        assert abs(rho - rho_check) / rho_check < 0.01
+