diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index ea49d2e..285f065 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -34,6 +34,9 @@ jobs:
           - python-version: "3.12"
             numpy-version: "1.26"
             os: ubuntu-latest
+          - python-version: "3.12"
+            numpy-version: "2.0"
+            os: ubuntu-latest
     steps:
       - uses: actions/checkout@v2
 
@@ -128,4 +131,4 @@ jobs:
 
       - name: Check type information
         run: |
-          mypy quantities
\ No newline at end of file
+          mypy quantities
diff --git a/quantities/__init__.py b/quantities/__init__.py
index d315fbb..d534d58 100644
--- a/quantities/__init__.py
+++ b/quantities/__init__.py
@@ -265,6 +265,8 @@
 
 """
 
+class QuantitiesDeprecationWarning(DeprecationWarning):
+    pass
 
 from ._version import __version__
 
diff --git a/quantities/dimensionality.py b/quantities/dimensionality.py
index c8d6002..22207de 100644
--- a/quantities/dimensionality.py
+++ b/quantities/dimensionality.py
@@ -9,6 +9,8 @@
 from .registry import unit_registry
 from .decorators import memoize
 
+_np_version = tuple(map(int, np.__version__.split('.')))
+
 def assert_isinstance(obj, types):
     try:
         assert isinstance(obj, types)
@@ -329,10 +331,11 @@ def _d_copy(q1, out=None):
 
 def _d_clip(a1, a2, a3, q):
     return q.dimensionality
-try:
+
+if _np_version < (2, 0, 0):
     p_dict[np.core.umath.clip] = _d_clip
-except AttributeError:
-    pass # For compatibility with Numpy < 1.17 when clip wasn't a ufunc yet
+else:
+    p_dict[np.clip] = _d_clip
 
 def _d_sqrt(q1, out=None):
     return q1._dimensionality**0.5
diff --git a/quantities/quantity.py b/quantities/quantity.py
index 612992e..29b6df3 100644
--- a/quantities/quantity.py
+++ b/quantities/quantity.py
@@ -3,10 +3,11 @@
 
 import copy
 from functools import wraps
+import warnings
 
 import numpy as np
 
-from . import markup
+from . import markup, QuantitiesDeprecationWarning
 from .dimensionality import Dimensionality, p_dict
 from .registry import unit_registry
 from .decorators import with_doc
@@ -114,15 +115,19 @@ class Quantity(np.ndarray):
     # TODO: what is an appropriate value?
     __array_priority__ = 21
 
-    def __new__(cls, data, units='', dtype=None, copy=True):
+    def __new__(cls, data, units='', dtype=None, copy=None):
+        if copy is not None:
+            warnings.warn(("The 'copy' argument in Quantity is deprecated and will be removed in the future. "
+                           "The argument has no effect since quantities-0.16.0 (to aid numpy-2.0 support)."),
+                           QuantitiesDeprecationWarning, stacklevel=2)
         if isinstance(data, Quantity):
             if units:
                 data = data.rescale(units)
             if isinstance(data, unit_registry['UnitQuantity']):
                 return 1*data
-            return np.array(data, dtype=dtype, copy=copy, subok=True).view(cls)
+            return np.asanyarray(data, dtype=dtype).view(cls)
 
-        ret = np.array(data, dtype=dtype, copy=copy).view(cls)
+        ret = np.asarray(data, dtype=dtype).view(cls)
         ret._dimensionality.update(validate_dimensionality(units))
         return ret
 
@@ -210,8 +215,8 @@ def rescale(self, units=None, dtype=None):
             dtype = self.dtype
         if self.dimensionality == to_dims:
             return self.astype(dtype)
-        to_u = Quantity(1.0, to_dims)
-        from_u = Quantity(1.0, self.dimensionality)
+        to_u = Quantity(1.0, to_dims, dtype=dtype)
+        from_u = Quantity(1.0, self.dimensionality, dtype=dtype)
         try:
             cf = get_conversion_factor(from_u, to_u)
         except AssertionError:
@@ -219,6 +224,8 @@ def rescale(self, units=None, dtype=None):
                 'Unable to convert between units of "%s" and "%s"'
                 %(from_u._dimensionality, to_u._dimensionality)
             )
+        if np.dtype(dtype).kind in 'fc':
+            cf = np.array(cf, dtype=dtype)
         new_magnitude = cf*self.magnitude
         dtype = np.result_type(dtype, new_magnitude)
         return Quantity(new_magnitude, to_u, dtype=dtype)
@@ -272,7 +279,7 @@ def __array_prepare__(self, obj, context=None):
         uf, objs, huh = context
         if uf.__name__.startswith('is'):
             return obj
-        #print self, obj, res, uf, objs
+
         try:
             res._dimensionality = p_dict[uf](*objs)
         except KeyError:
@@ -283,11 +290,21 @@ def __array_prepare__(self, obj, context=None):
                 )
         return res
 
-    def __array_wrap__(self, obj, context=None):
-        if not isinstance(obj, Quantity):
-            # backwards compatibility with numpy-1.3
-            obj = self.__array_prepare__(obj, context)
-        return obj
+    def __array_wrap__(self, obj, context=None, return_scalar=False):
+        _np_version = tuple(map(int, np.__version__.split('.')))
+        # For NumPy < 2.0 we do old behavior
+        if _np_version < (2, 0, 0):
+            if not isinstance(obj, Quantity):
+                return self.__array_prepare__(obj, context)
+            else:
+                return obj
+        # For NumPy > 2.0 we either do the prepare or the wrap
+        else:
+            if not isinstance(obj, Quantity):
+                return self.__array_prepare__(obj, context)
+            else:
+                return super().__array_wrap__(obj, context, return_scalar)
+
 
     @with_doc(np.ndarray.__add__)
     @scale_other_units
@@ -476,7 +493,7 @@ def sum(self, axis=None, dtype=None, out=None):
         ret = self.magnitude.sum(axis, dtype, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -487,8 +504,7 @@ def nansum(self, axis=None, dtype=None, out=None):
         import numpy as np
         return Quantity(
             np.nansum(self.magnitude, axis, dtype, out),
-            self.dimensionality,
-            copy=False
+            self.dimensionality
         )
 
     @with_doc(np.ndarray.fill)
@@ -523,7 +539,7 @@ def argsort(self, axis=-1, kind='quick', order=None):
     @with_doc(np.ndarray.searchsorted)
     def searchsorted(self,values, side='left'):
         if not isinstance (values, Quantity):
-            values = Quantity(values, copy=False)
+            values = Quantity(values)
 
         if values._dimensionality != self._dimensionality:
             raise ValueError("values does not have the same units as self")
@@ -539,7 +555,7 @@ def max(self, axis=None, out=None):
         ret = self.magnitude.max(axis, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -553,8 +569,7 @@ def argmax(self, axis=None, out=None):
     def nanmax(self, axis=None, out=None):
         return Quantity(
             np.nanmax(self.magnitude),
-            self.dimensionality,
-            copy=False
+            self.dimensionality
         )
 
     @with_doc(np.ndarray.min)
@@ -562,7 +577,7 @@ def min(self, axis=None, out=None):
         ret = self.magnitude.min(axis, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -572,8 +587,7 @@ def min(self, axis=None, out=None):
     def nanmin(self, axis=None, out=None):
         return Quantity(
             np.nanmin(self.magnitude),
-            self.dimensionality,
-            copy=False
+            self.dimensionality
         )
 
     @with_doc(np.ndarray.argmin)
@@ -590,10 +604,10 @@ def nanargmax(self,axis=None, out=None):
 
     @with_doc(np.ndarray.ptp)
     def ptp(self, axis=None, out=None):
-        ret = self.magnitude.ptp(axis, None if out is None else out.magnitude)
+        ret = np.ptp(self.magnitude, axis, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -620,7 +634,7 @@ def clip(self, min=None, max=None, out=None):
         )
         dim = self.dimensionality
         if out is None:
-            return Quantity(clipped, dim, copy=False)
+            return Quantity(clipped, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -631,7 +645,7 @@ def round(self, decimals=0, out=None):
         ret = self.magnitude.round(decimals, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -642,7 +656,7 @@ def trace(self, offset=0, axis1=0, axis2=1, dtype=None, out=None):
         ret = self.magnitude.trace(offset, axis1, axis2, dtype, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -652,8 +666,7 @@ def trace(self, offset=0, axis1=0, axis2=1, dtype=None, out=None):
     def squeeze(self, axis=None):
         return Quantity(
             self.magnitude.squeeze(axis),
-            self.dimensionality,
-            copy=False
+            self.dimensionality
         )
 
     @with_doc(np.ndarray.mean)
@@ -661,7 +674,7 @@ def mean(self, axis=None, dtype=None, out=None):
         ret = self.magnitude.mean(axis, dtype, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -672,15 +685,14 @@ def nanmean(self, axis=None, dtype=None, out=None):
         import numpy as np
         return Quantity(
             np.nanmean(self.magnitude, axis, dtype, out),
-            self.dimensionality,
-            copy=False)
+            self.dimensionality)
 
     @with_doc(np.ndarray.var)
     def var(self, axis=None, dtype=None, out=None, ddof=0):
         ret = self.magnitude.var(axis, dtype, out, ddof)
         dim = self._dimensionality**2
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -691,7 +703,7 @@ def std(self, axis=None, dtype=None, out=None, ddof=0):
         ret = self.magnitude.std(axis, dtype, out, ddof)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -701,8 +713,7 @@ def std(self, axis=None, dtype=None, out=None, ddof=0):
     def nanstd(self, axis=None, dtype=None, out=None, ddof=0):
         return Quantity(
             np.nanstd(self.magnitude, axis, dtype, out, ddof),
-            self._dimensionality,
-            copy=False
+            self._dimensionality
         )
 
     @with_doc(np.ndarray.prod)
@@ -715,7 +726,7 @@ def prod(self, axis=None, dtype=None, out=None):
         ret = self.magnitude.prod(axis, dtype, None if out is None else out.magnitude)
         dim = self._dimensionality**power
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -726,7 +737,7 @@ def cumsum(self, axis=None, dtype=None, out=None):
         ret = self.magnitude.cumsum(axis, dtype, None if out is None else out.magnitude)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if not isinstance(out, Quantity):
             raise TypeError("out parameter must be a Quantity")
         out._dimensionality = dim
@@ -743,7 +754,7 @@ def cumprod(self, axis=None, dtype=None, out=None):
         ret = self.magnitude.cumprod(axis, dtype, out)
         dim = self.dimensionality
         if out is None:
-            return Quantity(ret, dim, copy=False)
+            return Quantity(ret, dim)
         if isinstance(out, Quantity):
             out._dimensionality = dim
         return out
diff --git a/quantities/tests/common.py b/quantities/tests/common.py
index fb5d974..ed28c83 100644
--- a/quantities/tests/common.py
+++ b/quantities/tests/common.py
@@ -19,7 +19,11 @@ def assertQuantityEqual(self, q1, q2, msg=None, delta=None):
         Make sure q1 and q2 are the same quantities to within the given
         precision.
         """
-        delta = 1e-5 if delta is None else delta
+        if delta is None:
+            # NumPy 2 introduced float16, so we base tolerance on machine epsilon
+            delta1 = np.finfo(q1.dtype).eps if isinstance(q1, np.ndarray) and q1.dtype.kind in 'fc' else 1e-15
+            delta2 = np.finfo(q2.dtype).eps if isinstance(q2, np.ndarray) and q2.dtype.kind in 'fc' else 1e-15
+            delta = max(delta1, delta2)**0.3
         msg = '' if msg is None else ' (%s)' % msg
 
         q1 = Quantity(q1)
@@ -28,6 +32,7 @@ def assertQuantityEqual(self, q1, q2, msg=None, delta=None):
             raise self.failureException(
                 f"Shape mismatch ({q1.shape} vs {q2.shape}){msg}"
                 )
+
         if not np.all(np.abs(q1.magnitude - q2.magnitude) < delta):
             raise self.failureException(
                 "Magnitudes differ by more than %g (%s vs %s)%s"
diff --git a/quantities/tests/test_arithmetic.py b/quantities/tests/test_arithmetic.py
index 7f955b4..e7e5a05 100644
--- a/quantities/tests/test_arithmetic.py
+++ b/quantities/tests/test_arithmetic.py
@@ -48,32 +48,9 @@ def check(f, *args, **kwargs):
     return (new, )
 
 
-class iter_dtypes:
-
-    def __init__(self):
-        self._i = 1
-        self._typeDict = np.sctypeDict.copy()
-        self._typeDict[17] = int
-        self._typeDict[18] = long
-        self._typeDict[19] = float
-        self._typeDict[20] = complex
-
-    def __iter__(self):
-        return self
-
-    def __next__(self):
-        if self._i > 20:
-            raise StopIteration
-
-        i = self._i
-        self._i += 1
-        return self._typeDict[i]
-
-    def next(self):
-        return self.__next__()
-
 def get_dtypes():
-    return list(iter_dtypes())
+    numeric_dtypes = 'iufc'  # https://numpy.org/doc/stable/reference/generated/numpy.dtype.kind.html
+    return [v for v in np.sctypeDict.values() if np.dtype(v).kind in numeric_dtypes] + [int, long, float, complex]
 
 
 class iter_types:
diff --git a/quantities/tests/test_methods.py b/quantities/tests/test_methods.py
index 8989074..ba184a3 100644
--- a/quantities/tests/test_methods.py
+++ b/quantities/tests/test_methods.py
@@ -1,4 +1,5 @@
-from .. import units as pq
+import warnings
+from .. import QuantitiesDeprecationWarning, units as pq
 from .common import TestCase
 import numpy as np
 
@@ -119,7 +120,10 @@ def methodWithOut(self, name, result, q=None, *args, **kw):
         )
         if isinstance(result, Quantity):
             # deliberately using an incompatible unit
-            out = Quantity(np.empty_like(result.magnitude), pq.s, copy=False)
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", category=QuantitiesDeprecationWarning)
+                out = Quantity(np.empty_like(result.magnitude), pq.s, copy=False)
+            # we can drop 'copy=False' above once the deprecation of the arg has expired.
         else:
             out = np.empty_like(result)
         ret = getattr(q.copy(), name)(*args, out=out, **kw)
diff --git a/quantities/tests/test_umath.py b/quantities/tests/test_umath.py
index d9caa16..4c113d1 100644
--- a/quantities/tests/test_umath.py
+++ b/quantities/tests/test_umath.py
@@ -3,6 +3,8 @@
 from .. import units as pq
 from .common import TestCase, unittest
 
+_np_version = tuple(map(int, np.__version__.split('.')))
+
 
 class TestUmath(TestCase):
 
@@ -54,7 +56,13 @@ def test_cross(self):
         self.assertQuantityEqual(np.cross(a,b), [-15,-2,39]*pq.kPa*pq.m**2)
 
     def test_trapz(self):
-        self.assertQuantityEqual(np.trapz(self.q, dx = 1*pq.m), 7.5 * pq.J*pq.m)
+        # np.trapz is deprecated, np.trapezoid is recommend
+        if _np_version < (2, 0, 0):
+            self.assertQuantityEqual(np.trapz(self.q, dx = 1*pq.m), 7.5 * pq.J*pq.m)
+
+    def test_trapezoid(self):
+        if _np_version >= (2, 0, 0):
+            self.assertQuantityEqual(np.trapezoid(self.q, dx = 1*pq.m), 7.5 * pq.J*pq.m)
 
     def test_sinh(self):
         q = [1, 2, 3, 4, 6] * pq.radian
diff --git a/quantities/uncertainquantity.py b/quantities/uncertainquantity.py
index ad6d4c8..319a6cc 100644
--- a/quantities/uncertainquantity.py
+++ b/quantities/uncertainquantity.py
@@ -2,8 +2,9 @@
 """
 
 import numpy as np
+import warnings
 
-from . import markup
+from . import markup, QuantitiesDeprecationWarning
 from .quantity import Quantity, scale_other_units
 from .registry import unit_registry
 from .decorators import with_doc
@@ -13,15 +14,20 @@ class UncertainQuantity(Quantity):
     # TODO: what is an appropriate value?
     __array_priority__ = 22
 
-    def __new__(cls, data, units='', uncertainty=None, dtype='d', copy=True):
-        ret = Quantity.__new__(cls, data, units, dtype, copy)
+    def __new__(cls, data, units='', uncertainty=None, dtype='d', copy=None):
+        if copy is not None:
+            warnings.warn(("The 'copy' argument in UncertainQuantity is deprecated and will be removed in the future. "
+                           "The argument has no effect since quantities-0.16.0 (to aid numpy-2.0 support)."),
+                           QuantitiesDeprecationWarning, stacklevel=2)
+        ret = Quantity.__new__(cls, data, units, dtype)
         # _uncertainty initialized to be dimensionless by __array_finalize__:
         ret._uncertainty._dimensionality = ret._dimensionality
 
         if uncertainty is not None:
             ret.uncertainty = Quantity(uncertainty, ret._dimensionality)
         elif isinstance(data, UncertainQuantity):
-            if copy or ret._dimensionality != uncertainty._dimensionality:
+            is_copy = id(data) == id(ret)
+            if is_copy or ret._dimensionality != uncertainty._dimensionality:
                 uncertainty = data.uncertainty.rescale(ret.units)
             ret.uncertainty = uncertainty
 
@@ -50,7 +56,7 @@ def uncertainty(self):
     @uncertainty.setter
     def uncertainty(self, uncertainty):
         if not isinstance(uncertainty, Quantity):
-            uncertainty = Quantity(uncertainty, copy=False)
+            uncertainty = Quantity(uncertainty)
         try:
             assert self.shape == uncertainty.shape
         except AssertionError:
@@ -78,7 +84,6 @@ def __array_finalize__(self, obj):
             Quantity(
                 np.zeros(self.shape, self.dtype),
                 self._dimensionality,
-                copy=False
             )
         )
 
@@ -87,7 +92,7 @@ def __array_finalize__(self, obj):
     def __add__(self, other):
         res = super().__add__(other)
         u = (self.uncertainty**2+other.uncertainty**2)**0.5
-        return UncertainQuantity(res, uncertainty=u, copy=False)
+        return UncertainQuantity(res, uncertainty=u)
 
     @with_doc(Quantity.__radd__, use_header=False)
     @scale_other_units
@@ -99,13 +104,13 @@ def __radd__(self, other):
     def __sub__(self, other):
         res = super().__sub__(other)
         u = (self.uncertainty**2+other.uncertainty**2)**0.5
-        return UncertainQuantity(res, uncertainty=u, copy=False)
+        return UncertainQuantity(res, uncertainty=u)
 
     @with_doc(Quantity.__rsub__, use_header=False)
     @scale_other_units
     def __rsub__(self, other):
         if not isinstance(other, UncertainQuantity):
-            other = UncertainQuantity(other, copy=False)
+            other = UncertainQuantity(other)
 
         return UncertainQuantity.__sub__(other, self)
 
@@ -118,7 +123,7 @@ def __mul__(self, other):
             ru = (sru**2+oru**2)**0.5
             u = res.view(Quantity) * ru
         except AttributeError:
-            other = np.array(other, copy=False, subok=True)
+            other = np.asanyarray(other)
             u = (self.uncertainty**2*other**2)**0.5
 
         res._uncertainty = u
@@ -140,7 +145,7 @@ def __truediv__(self, other):
             ru = (sru**2+oru**2)**0.5
             u = res.view(Quantity) * ru
         except AttributeError:
-            other = np.array(other, copy=False, subok=True)
+            other = np.asanyarray(other)
             u = (self.uncertainty**2/other**2)**0.5
 
         res._uncertainty = u
@@ -150,7 +155,7 @@ def __truediv__(self, other):
     def __rtruediv__(self, other):
         temp = UncertainQuantity(
             1/self.magnitude, self.dimensionality**-1,
-            self.relative_uncertainty/self.magnitude, copy=False
+            self.relative_uncertainty/self.magnitude
         )
         return other * temp
 
@@ -165,8 +170,7 @@ def __getitem__(self, key):
         return UncertainQuantity(
             self.magnitude[key],
             self._dimensionality,
-            self.uncertainty[key],
-            copy=False
+            self.uncertainty[key]
         )
 
     @with_doc(Quantity.__repr__, use_header=False)
@@ -198,8 +202,7 @@ def sum(self, axis=None, dtype=None, out=None):
         return UncertainQuantity(
             self.magnitude.sum(axis, dtype, out),
             self.dimensionality,
-            (np.sum(self.uncertainty.magnitude**2, axis))**0.5,
-            copy=False
+            (np.sum(self.uncertainty.magnitude**2, axis))**0.5
         )
 
     @with_doc(np.nansum)
@@ -207,8 +210,7 @@ def nansum(self, axis=None, dtype=None, out=None):
         return UncertainQuantity(
             np.nansum(self.magnitude, axis, dtype, out),
             self.dimensionality,
-            (np.nansum(self.uncertainty.magnitude**2, axis))**0.5,
-            copy=False
+            (np.nansum(self.uncertainty.magnitude**2, axis))**0.5
         )
 
     @with_doc(np.ndarray.mean)
@@ -216,8 +218,8 @@ def mean(self, axis=None, dtype=None, out=None):
         return UncertainQuantity(
             self.magnitude.mean(axis, dtype, out),
             self.dimensionality,
-            ((1.0/np.size(self,axis))**2 * np.sum(self.uncertainty.magnitude**2, axis))**0.5,
-            copy=False)
+            ((1.0/np.size(self,axis))**2 * np.sum(self.uncertainty.magnitude**2, axis))**0.5
+        )
 
     @with_doc(np.nanmean)
     def nanmean(self, axis=None, dtype=None, out=None):
@@ -225,8 +227,8 @@ def nanmean(self, axis=None, dtype=None, out=None):
         return UncertainQuantity(
             np.nanmean(self.magnitude, axis, dtype, out),
             self.dimensionality,
-            ((1.0/size)**2 * np.nansum(np.nan_to_num(self.uncertainty.magnitude)**2, axis))**0.5,
-            copy=False)
+            ((1.0/size)**2 * np.nansum(np.nan_to_num(self.uncertainty.magnitude)**2, axis))**0.5
+        )
 
     @with_doc(np.sqrt)
     def sqrt(self, out=None):