diff --git a/yt/fields/species_fields.py b/yt/fields/species_fields.py
index acd2ccacc09..b8ee8bbcfe9 100644
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -44,7 +44,7 @@ def _create_number_density_func(ftype, species):
 
     def _number_density(field, data):
         weight = formula.weight  # This is in AMU
-        weight *= data.ds.units.physical_constants.amu_cgs
+        weight *= data.ds.quan(1.0, "amu").in_cgs()
         return data[ftype, f"{species}_density"] / weight
 
     return _number_density
@@ -252,7 +252,7 @@ def add_nuclei_density_fields(registry, ftype):
 def _default_nuclei_density(field, data):
     ftype = field.name[0]
     element = field.name[1][: field.name[1].find("_")]
-    amu_cgs = data.ds.units.physical_constants.amu_cgs
+    amu_cgs = data.ds.quan(1.0, "amu").in_cgs()
     if element == "El":
         # This is for determining the electron number density.
         # If we got here, this assumes full ionization!
@@ -273,7 +273,7 @@ def _nuclei_density(field, data):
         return (
             data[(ftype, nuclei_mass_field)]
             / ChemicalFormula(element).weight
-            / data.ds.units.physical_constants.amu_cgs
+            / data.ds.quan(1.0, "amu").in_cgs()
         )
     metal_field = f"{element}_metallicity"
     if (ftype, metal_field) in data.ds.field_info:
@@ -281,7 +281,7 @@ def _nuclei_density(field, data):
             data[ftype, "density"]
             * data[(ftype, metal_field)]
             / ChemicalFormula(element).weight
-            / data.ds.units.physical_constants.amu_cgs
+            / data.ds.quan(1.0, "amu").in_cgs()
         )
 
     field_data = np.zeros_like(
diff --git a/yt/fields/tests/test_species_fields.py b/yt/fields/tests/test_species_fields.py
index 16360b64c51..165cb785f83 100644
--- a/yt/fields/tests/test_species_fields.py
+++ b/yt/fields/tests/test_species_fields.py
@@ -26,7 +26,7 @@ def test_default_species_fields():
     # Test fully ionized case
     dsi = fake_random_ds(32, default_species_fields="ionized")
     spi = dsi.sphere("c", (0.2, "unitary"))
-    amu_cgs = dsi.units.physical_constants.amu_cgs
+    amu_cgs = dsi.quan(1.0, "amu").in_cgs()
 
     mueinv = 1.0 * _primordial_mass_fraction["H"] / ChemicalFormula("H").weight
     mueinv *= spi["index", "ones"]
@@ -56,7 +56,7 @@ def test_default_species_fields():
 
     dsn = fake_random_ds(32, default_species_fields="neutral")
     spn = dsn.sphere("c", (0.2, "unitary"))
-    amu_cgs = dsn.units.physical_constants.amu_cgs
+    amu_cgs = dsn.quan(1.0, "amu").in_cgs()
 
     assert ("gas", "El_number_density") not in ds.derived_field_list
 
diff --git a/yt/frontends/arepo/fields.py b/yt/frontends/arepo/fields.py
index 8c5917b59a1..17d426d1ab1 100644
--- a/yt/frontends/arepo/fields.py
+++ b/yt/frontends/arepo/fields.py
@@ -132,7 +132,7 @@ def _h_p1_fraction(field, data):
             # try first to use the H_fraction, but otherwise we assume the
             # cosmic value for hydrogen to generate the H_number_density
             if (ptype, "NeutralHydrogenAbundance") not in self.field_list:
-                m_u = self.ds.units.physical_constants.amu_cgs
+                m_u = self.ds.quan(1.0, "amu").in_cgs()
                 A_H = ChemicalFormula("H").weight
                 if (ptype, "GFM_Metals_00") in self.field_list: