diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 1ee0da3dc7..590a5698a7 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -53,6 +53,7 @@ if sys.version_info < (3, 9):
 
 
 import qiskit._accelerate
+import qiskit._numpy_compat
 
 # Globally define compiled submodules. The normal import mechanism will not find compiled submodules
 # in _accelerate because it relies on file paths, but PyO3 generates only one shared library file.
diff --git a/qiskit/_numpy_compat.py b/qiskit/_numpy_compat.py
new file mode 100644
index 0000000000..a6c06671c9
--- /dev/null
+++ b/qiskit/_numpy_compat.py
@@ -0,0 +1,73 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Compatiblity helpers for the Numpy 1.x to 2.0 transition."""
+
+import re
+import typing
+import warnings
+
+import numpy as np
+
+# This version pattern is taken from the pypa packaging project:
+# https://github.com/pypa/packaging/blob/21.3/packaging/version.py#L223-L254 which is dual licensed
+# Apache 2.0 and BSD see the source for the original authors and other details.
+_VERSION_PATTERN = r"""
+    v?
+    (?:
+        (?:(?P<epoch>[0-9]+)!)?                           # epoch
+        (?P<release>[0-9]+(?:\.[0-9]+)*)                  # release segment
+        (?P<pre>                                          # pre-release
+            [-_\.]?
+            (?P<pre_l>(a|b|c|rc|alpha|beta|pre|preview))
+            [-_\.]?
+            (?P<pre_n>[0-9]+)?
+        )?
+        (?P<post>                                         # post release
+            (?:-(?P<post_n1>[0-9]+))
+            |
+            (?:
+                [-_\.]?
+                (?P<post_l>post|rev|r)
+                [-_\.]?
+                (?P<post_n2>[0-9]+)?
+            )
+        )?
+        (?P<dev>                                          # dev release
+            [-_\.]?
+            (?P<dev_l>dev)
+            [-_\.]?
+            (?P<dev_n>[0-9]+)?
+        )?
+    )
+    (?:\+(?P<local>[a-z0-9]+(?:[-_\.][a-z0-9]+)*))?       # local version
+"""
+
+VERSION = np.lib.NumpyVersion(np.__version__)
+VERSION_PARTS: typing.Tuple[int, ...]
+"""The numeric parts of the Numpy release version, e.g. ``(2, 0, 0)``.  Does not include pre- or
+post-release markers (e.g. ``rc1``)."""
+if match := re.fullmatch(_VERSION_PATTERN, np.__version__, flags=re.VERBOSE | re.IGNORECASE):
+    # Assuming Numpy won't ever introduce epochs, and we don't care about pre/post markers.
+    VERSION_PARTS = tuple(int(x) for x in match["release"].split("."))
+else:
+    # Just guess a version.  We know all existing Numpys have good version strings, so the only way
+    # this should trigger is from a new or a dev version.
+    warnings.warn(
+        f"Unrecognized version string for Numpy: '{np.__version__}'.  Assuming Numpy 2.0.",
+        RuntimeWarning,
+    )
+    VERSION_PARTS = (2, 0, 0)
+
+COPY_ONLY_IF_NEEDED = None if VERSION_PARTS >= (2, 0, 0) else False
+"""The sentinel value given to ``np.array`` and ``np.ndarray.astype`` (etc) to indicate that a copy
+should be made only if required."""
diff --git a/qiskit/circuit/__init__.py b/qiskit/circuit/__init__.py
index 2e8603af9f..9fbefb4c5d 100644
--- a/qiskit/circuit/__init__.py
+++ b/qiskit/circuit/__init__.py
@@ -833,7 +833,7 @@ method <https://numpy.org/devdocs/user/basics.interoperability.html#the-array-me
 ``__array__``.  This is used by :meth:`Gate.to_matrix`, and has the signature:
 
 .. currentmodule:: None
-.. py:method:: __array__(dtype=None)
+.. py:method:: __array__(dtype=None, copy=None)
 
     Return a Numpy array representing the gate.  This can use the gate's :attr:`~Instruction.params`
     field, and may assume that these are numeric values (assuming the subclass expects that) and not
@@ -875,7 +875,9 @@ are unitary, so should be a :class:`Gate`::
             # Also we have an efficient representation of power.
             return RXZGate(exponent * self.params[0])
 
-        def __array__(self, dtype=None):
+        def __array__(self, dtype=None, copy=None):
+            if copy is False:
+                raise ValueError("unable to avoid copy while creating an array as requested")
             cos = math.cos(0.5 * self.params[0])
             isin = 1j * math.sin(0.5 * self.params[0])
             return np.array([
@@ -1340,6 +1342,7 @@ In both these cases, the matrix form of :class:`.CCXGate` in ``ctrl_state = 1``
 """
 
 from .exceptions import CircuitError
+from . import _utils
 from .quantumcircuit import QuantumCircuit
 from .classicalregister import ClassicalRegister, Clbit
 from .quantumregister import QuantumRegister, Qubit, AncillaRegister, AncillaQubit
diff --git a/qiskit/circuit/_utils.py b/qiskit/circuit/_utils.py
index cfde85bad8..86a058e885 100644
--- a/qiskit/circuit/_utils.py
+++ b/qiskit/circuit/_utils.py
@@ -15,6 +15,8 @@ This module contains utility functions for circuits.
 
 import math
 import numpy
+
+from qiskit import _numpy_compat
 from qiskit.exceptions import QiskitError
 from qiskit.circuit.exceptions import CircuitError
 from .parametervector import ParameterVectorElement
@@ -117,8 +119,9 @@ def with_gate_array(base_array):
     nonwritable = numpy.array(base_array, dtype=numpy.complex128)
     nonwritable.setflags(write=False)
 
-    def __array__(_self, dtype=None):
-        return numpy.asarray(nonwritable, dtype=dtype)
+    def __array__(_self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = nonwritable.dtype if dtype is None else dtype
+        return numpy.array(nonwritable, dtype=dtype, copy=copy)
 
     def decorator(cls):
         if hasattr(cls, "__array__"):
@@ -149,15 +152,21 @@ def with_controlled_gate_array(base_array, num_ctrl_qubits, cached_states=None):
     if cached_states is None:
         nonwritables = [matrix_for_control_state(state) for state in range(2**num_ctrl_qubits)]
 
-        def __array__(self, dtype=None):
-            return numpy.asarray(nonwritables[self.ctrl_state], dtype=dtype)
+        def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+            arr = nonwritables[self.ctrl_state]
+            dtype = arr.dtype if dtype is None else dtype
+            return numpy.array(arr, dtype=dtype, copy=copy)
 
     else:
         nonwritables = {state: matrix_for_control_state(state) for state in cached_states}
 
-        def __array__(self, dtype=None):
-            if (out := nonwritables.get(self.ctrl_state)) is not None:
-                return numpy.asarray(out, dtype=dtype)
+        def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+            if (arr := nonwritables.get(self.ctrl_state)) is not None:
+                dtype = arr.dtype if dtype is None else dtype
+                return numpy.array(arr, dtype=dtype, copy=copy)
+
+            if copy is False and copy is not _numpy_compat.COPY_ONLY_IF_NEEDED:
+                raise ValueError("could not produce matrix without calculation")
             return numpy.asarray(
                 _compute_control_matrix(base, num_ctrl_qubits, self.ctrl_state), dtype=dtype
             )
diff --git a/qiskit/circuit/delay.py b/qiskit/circuit/delay.py
index 16d84d15cb..a333125a5a 100644
--- a/qiskit/circuit/delay.py
+++ b/qiskit/circuit/delay.py
@@ -17,9 +17,11 @@ import numpy as np
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.gate import Gate
+from qiskit.circuit import _utils
 from qiskit.circuit.parameterexpression import ParameterExpression
 
 
+@_utils.with_gate_array(np.eye(2, dtype=complex))
 class Delay(Instruction):
     """Do nothing and just delay/wait/idle for a specified duration."""
 
@@ -53,10 +55,6 @@ class Delay(Instruction):
         """Set the duration of this delay."""
         self.params = [duration]
 
-    def __array__(self, dtype=None):
-        """Return the identity matrix."""
-        return np.array([[1, 0], [0, 1]], dtype=dtype)
-
     def to_matrix(self) -> np.ndarray:
         """Return a Numpy.array for the unitary matrix. This has been
         added to enable simulation without making delay a full Gate type.
diff --git a/qiskit/circuit/library/generalized_gates/pauli.py b/qiskit/circuit/library/generalized_gates/pauli.py
index e8b063a75e..01bbd09c19 100644
--- a/qiskit/circuit/library/generalized_gates/pauli.py
+++ b/qiskit/circuit/library/generalized_gates/pauli.py
@@ -63,13 +63,13 @@ class PauliGate(Gate):
         r"""Return inverted pauli gate (itself)."""
         return PauliGate(self.params[0])  # self-inverse
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the pauli gate.
         i.e. tensor product of the paulis"""
         # pylint: disable=cyclic-import
         from qiskit.quantum_info.operators import Pauli
 
-        return Pauli(self.params[0]).__array__(dtype=dtype)
+        return Pauli(self.params[0]).__array__(dtype=dtype, copy=copy)
 
     def validate_parameter(self, parameter):
         if isinstance(parameter, str):
diff --git a/qiskit/circuit/library/generalized_gates/permutation.py b/qiskit/circuit/library/generalized_gates/permutation.py
index 8888344c78..776c69d94f 100644
--- a/qiskit/circuit/library/generalized_gates/permutation.py
+++ b/qiskit/circuit/library/generalized_gates/permutation.py
@@ -147,8 +147,11 @@ class PermutationGate(Gate):
 
         super().__init__(name="permutation", num_qubits=num_qubits, params=[pattern])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the Permutation gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
+
         nq = len(self.pattern)
         mat = np.zeros((2**nq, 2**nq), dtype=dtype)
 
diff --git a/qiskit/circuit/library/generalized_gates/unitary.py b/qiskit/circuit/library/generalized_gates/unitary.py
index 6180411422..1fd36e52e0 100644
--- a/qiskit/circuit/library/generalized_gates/unitary.py
+++ b/qiskit/circuit/library/generalized_gates/unitary.py
@@ -18,6 +18,7 @@ import math
 import typing
 import numpy
 
+from qiskit import _numpy_compat
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.controlledgate import ControlledGate
 from qiskit.circuit.annotated_operation import AnnotatedOperation, ControlModifier
@@ -118,10 +119,10 @@ class UnitaryGate(Gate):
             return False
         return matrix_equal(self.params[0], other.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
         """Return matrix for the unitary."""
-        # pylint: disable=unused-argument
-        return self.params[0]
+        dtype = self.params[0].dtype if dtype is None else dtype
+        return numpy.array(self.params[0], dtype=dtype, copy=copy)
 
     def inverse(self, annotated: bool = False):
         """Return the adjoint of the unitary."""
diff --git a/qiskit/circuit/library/hamiltonian_gate.py b/qiskit/circuit/library/hamiltonian_gate.py
index a87504a97b..2997d01ed4 100644
--- a/qiskit/circuit/library/hamiltonian_gate.py
+++ b/qiskit/circuit/library/hamiltonian_gate.py
@@ -21,6 +21,7 @@ import typing
 from numbers import Number
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.quantumregister import QuantumRegister
@@ -92,18 +93,22 @@ class HamiltonianGate(Gate):
         times_eq = self.params[1] == other.params[1]
         return operators_eq and times_eq
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return matrix for the unitary."""
-        # pylint: disable=unused-argument
         import scipy.linalg
 
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         try:
-            return scipy.linalg.expm(-1j * self.params[0] * float(self.params[1]))
+            time = float(self.params[1])
         except TypeError as ex:
             raise TypeError(
                 "Unable to generate Unitary matrix for "
                 "unbound t parameter {}".format(self.params[1])
             ) from ex
+        arr = scipy.linalg.expm(-1j * self.params[0] * time)
+        dtype = complex if dtype is None else dtype
+        return np.array(arr, dtype=dtype, copy=_numpy_compat.COPY_ONLY_IF_NEEDED)
 
     def inverse(self, annotated: bool = False):
         """Return the adjoint of the unitary."""
diff --git a/qiskit/circuit/library/standard_gates/global_phase.py b/qiskit/circuit/library/standard_gates/global_phase.py
index 50576bf17a..ccd758e472 100644
--- a/qiskit/circuit/library/standard_gates/global_phase.py
+++ b/qiskit/circuit/library/standard_gates/global_phase.py
@@ -69,10 +69,12 @@ class GlobalPhaseGate(Gate):
         """
         return GlobalPhaseGate(-self.params[0])
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the global_phase gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta = self.params[0]
-        return numpy.array([[numpy.exp(1j * theta)]], dtype=dtype)
+        return numpy.array([[numpy.exp(1j * theta)]], dtype=dtype or complex)
 
     def __eq__(self, other):
         if isinstance(other, GlobalPhaseGate):
diff --git a/qiskit/circuit/library/standard_gates/p.py b/qiskit/circuit/library/standard_gates/p.py
index 179be025bc..6de0307dc7 100644
--- a/qiskit/circuit/library/standard_gates/p.py
+++ b/qiskit/circuit/library/standard_gates/p.py
@@ -140,8 +140,10 @@ class PhaseGate(Gate):
         """
         return PhaseGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the Phase gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         lam = float(self.params[0])
         return numpy.array([[1, 0], [0, exp(1j * lam)]], dtype=dtype)
 
@@ -279,8 +281,10 @@ class CPhaseGate(ControlledGate):
         r"""Return inverted CPhase gate (:math:`CPhase(\lambda)^{\dagger} = CPhase(-\lambda)`)"""
         return CPhaseGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CPhase gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         eith = exp(1j * float(self.params[0]))
         if self.ctrl_state:
             return numpy.array(
diff --git a/qiskit/circuit/library/standard_gates/r.py b/qiskit/circuit/library/standard_gates/r.py
index 3fc537baef..9d4905e278 100644
--- a/qiskit/circuit/library/standard_gates/r.py
+++ b/qiskit/circuit/library/standard_gates/r.py
@@ -93,8 +93,10 @@ class RGate(Gate):
         """
         return RGate(-self.params[0], self.params[1])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the R gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, phi = float(self.params[0]), float(self.params[1])
         cos = math.cos(theta / 2)
         sin = math.sin(theta / 2)
diff --git a/qiskit/circuit/library/standard_gates/rx.py b/qiskit/circuit/library/standard_gates/rx.py
index 3483d5ebc9..eaa73cf87c 100644
--- a/qiskit/circuit/library/standard_gates/rx.py
+++ b/qiskit/circuit/library/standard_gates/rx.py
@@ -120,8 +120,10 @@ class RXGate(Gate):
         """
         return RXGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RX gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         cos = math.cos(self.params[0] / 2)
         sin = math.sin(self.params[0] / 2)
         return numpy.array([[cos, -1j * sin], [-1j * sin, cos]], dtype=dtype)
@@ -263,8 +265,10 @@ class CRXGate(ControlledGate):
         """
         return CRXGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CRX gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         half_theta = float(self.params[0]) / 2
         cos = math.cos(half_theta)
         isin = 1j * math.sin(half_theta)
diff --git a/qiskit/circuit/library/standard_gates/rxx.py b/qiskit/circuit/library/standard_gates/rxx.py
index 03e9d22dcc..c4e35e53d5 100644
--- a/qiskit/circuit/library/standard_gates/rxx.py
+++ b/qiskit/circuit/library/standard_gates/rxx.py
@@ -122,8 +122,10 @@ class RXXGate(Gate):
         """
         return RXXGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the RXX gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta2 = float(self.params[0]) / 2
         cos = math.cos(theta2)
         isin = 1j * math.sin(theta2)
diff --git a/qiskit/circuit/library/standard_gates/ry.py b/qiskit/circuit/library/standard_gates/ry.py
index b902887ee0..633a518bca 100644
--- a/qiskit/circuit/library/standard_gates/ry.py
+++ b/qiskit/circuit/library/standard_gates/ry.py
@@ -119,8 +119,10 @@ class RYGate(Gate):
         """
         return RYGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RY gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         cos = math.cos(self.params[0] / 2)
         sin = math.sin(self.params[0] / 2)
         return numpy.array([[cos, -sin], [sin, cos]], dtype=dtype)
@@ -258,8 +260,10 @@ class CRYGate(ControlledGate):
         ."""
         return CRYGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CRY gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         half_theta = float(self.params[0]) / 2
         cos = math.cos(half_theta)
         sin = math.sin(half_theta)
diff --git a/qiskit/circuit/library/standard_gates/ryy.py b/qiskit/circuit/library/standard_gates/ryy.py
index 50ce9b0c4f..98847b7b21 100644
--- a/qiskit/circuit/library/standard_gates/ryy.py
+++ b/qiskit/circuit/library/standard_gates/ryy.py
@@ -122,8 +122,10 @@ class RYYGate(Gate):
         """
         return RYYGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RYY gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta = float(self.params[0])
         cos = math.cos(theta / 2)
         isin = 1j * math.sin(theta / 2)
diff --git a/qiskit/circuit/library/standard_gates/rz.py b/qiskit/circuit/library/standard_gates/rz.py
index c7311b4a6e..3040f95683 100644
--- a/qiskit/circuit/library/standard_gates/rz.py
+++ b/qiskit/circuit/library/standard_gates/rz.py
@@ -130,10 +130,12 @@ class RZGate(Gate):
         """
         return RZGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RZ gate."""
         import numpy as np
 
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         ilam2 = 0.5j * float(self.params[0])
         return np.array([[exp(-ilam2), 0], [0, exp(ilam2)]], dtype=dtype)
 
@@ -276,10 +278,12 @@ class CRZGate(ControlledGate):
         """
         return CRZGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CRZ gate."""
         import numpy
 
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         arg = 1j * float(self.params[0]) / 2
         if self.ctrl_state:
             return numpy.array(
diff --git a/qiskit/circuit/library/standard_gates/rzx.py b/qiskit/circuit/library/standard_gates/rzx.py
index d59676663d..1f930ab422 100644
--- a/qiskit/circuit/library/standard_gates/rzx.py
+++ b/qiskit/circuit/library/standard_gates/rzx.py
@@ -166,10 +166,12 @@ class RZXGate(Gate):
         """
         return RZXGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RZX gate."""
         import numpy
 
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         half_theta = float(self.params[0]) / 2
         cos = math.cos(half_theta)
         isin = 1j * math.sin(half_theta)
diff --git a/qiskit/circuit/library/standard_gates/rzz.py b/qiskit/circuit/library/standard_gates/rzz.py
index 3a00fb7b73..5ca974764d 100644
--- a/qiskit/circuit/library/standard_gates/rzz.py
+++ b/qiskit/circuit/library/standard_gates/rzz.py
@@ -130,10 +130,12 @@ class RZZGate(Gate):
         """
         return RZZGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RZZ gate."""
         import numpy
 
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         itheta2 = 1j * float(self.params[0]) / 2
         return numpy.array(
             [
diff --git a/qiskit/circuit/library/standard_gates/u.py b/qiskit/circuit/library/standard_gates/u.py
index 81b48536f2..3d63189885 100644
--- a/qiskit/circuit/library/standard_gates/u.py
+++ b/qiskit/circuit/library/standard_gates/u.py
@@ -12,7 +12,7 @@
 
 """Two-pulse single-qubit gate."""
 import cmath
-import copy
+import copy as _copy
 import math
 from cmath import exp
 from typing import Optional, Union
@@ -136,8 +136,10 @@ class UGate(Gate):
             )
         return gate
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the U gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, phi, lam = (float(param) for param in self.params)
         cos = math.cos(theta / 2)
         sin = math.sin(theta / 2)
@@ -146,7 +148,7 @@ class UGate(Gate):
                 [cos, -exp(1j * lam) * sin],
                 [exp(1j * phi) * sin, exp(1j * (phi + lam)) * cos],
             ],
-            dtype=dtype,
+            dtype=dtype or complex,
         )
 
     def __eq__(self, other):
@@ -337,8 +339,10 @@ class CUGate(ControlledGate):
             ctrl_state=self.ctrl_state,
         )
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CU gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, phi, lam, gamma = (float(param) for param in self.params)
         cos = math.cos(theta / 2)
         sin = math.sin(theta / 2)
@@ -372,5 +376,5 @@ class CUGate(ControlledGate):
         # assuming that `params` will be a view onto the base gate's `_params`.
         memo = memo if memo is not None else {}
         out = super().__deepcopy__(memo)
-        out._params = copy.deepcopy(out._params, memo)
+        out._params = _copy.deepcopy(out._params, memo)
         return out
diff --git a/qiskit/circuit/library/standard_gates/u1.py b/qiskit/circuit/library/standard_gates/u1.py
index b92bea51a2..1d59cabae1 100644
--- a/qiskit/circuit/library/standard_gates/u1.py
+++ b/qiskit/circuit/library/standard_gates/u1.py
@@ -160,8 +160,10 @@ class U1Gate(Gate):
         """
         return U1Gate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the U1 gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         lam = float(self.params[0])
         return numpy.array([[1, 0], [0, numpy.exp(1j * lam)]], dtype=dtype)
 
@@ -304,8 +306,10 @@ class CU1Gate(ControlledGate):
         """
         return CU1Gate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CU1 gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         eith = exp(1j * float(self.params[0]))
         if self.ctrl_state:
             return numpy.array(
diff --git a/qiskit/circuit/library/standard_gates/u2.py b/qiskit/circuit/library/standard_gates/u2.py
index 021a38f4da..c8e4de96ef 100644
--- a/qiskit/circuit/library/standard_gates/u2.py
+++ b/qiskit/circuit/library/standard_gates/u2.py
@@ -127,8 +127,10 @@ class U2Gate(Gate):
         """
         return U2Gate(-self.params[1] - pi, -self.params[0] + pi)
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the U2 gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         isqrt2 = 1 / sqrt(2)
         phi, lam = self.params
         phi, lam = float(phi), float(lam)
@@ -137,5 +139,5 @@ class U2Gate(Gate):
                 [isqrt2, -exp(1j * lam) * isqrt2],
                 [exp(1j * phi) * isqrt2, exp(1j * (phi + lam)) * isqrt2],
             ],
-            dtype=dtype,
+            dtype=dtype or complex,
         )
diff --git a/qiskit/circuit/library/standard_gates/u3.py b/qiskit/circuit/library/standard_gates/u3.py
index c92a48ab52..62c1e33b96 100644
--- a/qiskit/circuit/library/standard_gates/u3.py
+++ b/qiskit/circuit/library/standard_gates/u3.py
@@ -149,8 +149,10 @@ class U3Gate(Gate):
         qc.u(self.params[0], self.params[1], self.params[2], 0)
         self.definition = qc
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the U3 gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, phi, lam = self.params
         theta, phi, lam = float(theta), float(phi), float(lam)
         cos = math.cos(theta / 2)
@@ -160,7 +162,7 @@ class U3Gate(Gate):
                 [cos, -exp(1j * lam) * sin],
                 [exp(1j * phi) * sin, exp(1j * (phi + lam)) * cos],
             ],
-            dtype=dtype,
+            dtype=dtype or complex,
         )
 
 
@@ -305,8 +307,10 @@ class CU3Gate(ControlledGate):
             -self.params[0], -self.params[2], -self.params[1], ctrl_state=self.ctrl_state
         )
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CU3 gate."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, phi, lam = self.params
         theta, phi, lam = float(theta), float(phi), float(lam)
         cos = math.cos(theta / 2)
@@ -319,7 +323,7 @@ class CU3Gate(ControlledGate):
                     [0, 0, 1, 0],
                     [0, exp(1j * phi) * sin, 0, exp(1j * (phi + lam)) * cos],
                 ],
-                dtype=dtype,
+                dtype=dtype or complex,
             )
         else:
             return numpy.array(
@@ -329,7 +333,7 @@ class CU3Gate(ControlledGate):
                     [exp(1j * phi) * sin, 0, exp(1j * (phi + lam)) * cos, 0],
                     [0, 0, 0, 1],
                 ],
-                dtype=dtype,
+                dtype=dtype or complex,
             )
 
 
diff --git a/qiskit/circuit/library/standard_gates/xx_minus_yy.py b/qiskit/circuit/library/standard_gates/xx_minus_yy.py
index 387a23ad05..4bf4ab80ec 100644
--- a/qiskit/circuit/library/standard_gates/xx_minus_yy.py
+++ b/qiskit/circuit/library/standard_gates/xx_minus_yy.py
@@ -169,8 +169,10 @@ class XXMinusYYGate(Gate):
         theta, beta = self.params
         return XXMinusYYGate(-theta, beta)
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Gate matrix."""
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         theta, beta = self.params
         cos = math.cos(theta / 2)
         sin = math.sin(theta / 2)
diff --git a/qiskit/circuit/library/standard_gates/xx_plus_yy.py b/qiskit/circuit/library/standard_gates/xx_plus_yy.py
index b69ba49de3..a7b62175f2 100644
--- a/qiskit/circuit/library/standard_gates/xx_plus_yy.py
+++ b/qiskit/circuit/library/standard_gates/xx_plus_yy.py
@@ -15,6 +15,9 @@ import math
 from cmath import exp
 from math import pi
 from typing import Optional
+
+import numpy
+
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.quantumregister import QuantumRegister
 from qiskit.circuit.parameterexpression import ParameterValueType
@@ -167,10 +170,10 @@ class XXPlusYYGate(Gate):
         """
         return XXPlusYYGate(-self.params[0], self.params[1])
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the XX+YY gate."""
-        import numpy
-
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
         half_theta = float(self.params[0]) / 2
         beta = float(self.params[1])
         cos = math.cos(half_theta)
diff --git a/qiskit/primitives/backend_sampler_v2.py b/qiskit/primitives/backend_sampler_v2.py
index 6a990465c5..51d1ded150 100644
--- a/qiskit/primitives/backend_sampler_v2.py
+++ b/qiskit/primitives/backend_sampler_v2.py
@@ -268,7 +268,7 @@ def _memory_array(results: list[list[str]], num_bytes: int) -> NDArray[np.uint8]
             # no measure in a circuit
             data = np.zeros((len(memory), num_bytes), dtype=np.uint8)
         lst.append(data)
-    ary = np.array(lst, copy=False)
+    ary = np.asarray(lst)
     return np.unpackbits(ary, axis=-1, bitorder="big")
 
 
diff --git a/qiskit/primitives/containers/observables_array.py b/qiskit/primitives/containers/observables_array.py
index 0d0322dc6a..21c415d758 100644
--- a/qiskit/primitives/containers/observables_array.py
+++ b/qiskit/primitives/containers/observables_array.py
@@ -101,10 +101,10 @@ class ObservablesArray(ShapedMixin):
         """Convert to a nested list"""
         return self._array.tolist()
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Convert to an Numpy.ndarray"""
         if dtype is None or dtype == object:
-            return self._array
+            return self._array.copy() if copy else self._array
         raise ValueError("Type must be 'None' or 'object'")
 
     @overload
diff --git a/qiskit/primitives/containers/shape.py b/qiskit/primitives/containers/shape.py
index 952916cd67..6d893f46c1 100644
--- a/qiskit/primitives/containers/shape.py
+++ b/qiskit/primitives/containers/shape.py
@@ -85,7 +85,7 @@ def array_coerce(arr: ArrayLike | Shaped) -> NDArray | Shaped:
     """
     if isinstance(arr, Shaped):
         return arr
-    return np.array(arr, copy=False)
+    return np.asarray(arr)
 
 
 def _flatten_to_ints(arg: ShapeInput) -> Iterable[int]:
diff --git a/qiskit/quantum_info/operators/channel/chi.py b/qiskit/quantum_info/operators/channel/chi.py
index 7cd5fce525..ee0ddaa453 100644
--- a/qiskit/quantum_info/operators/channel/chi.py
+++ b/qiskit/quantum_info/operators/channel/chi.py
@@ -16,10 +16,11 @@ Chi-matrix representation of a Quantum Channel.
 """
 
 from __future__ import annotations
-import copy
+import copy as _copy
 import math
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
@@ -131,10 +132,9 @@ class Chi(QuantumChannel):
             raise QiskitError("Input is not an n-qubit Chi matrix.")
         super().__init__(chi_mat, num_qubits=num_qubits)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     @property
     def _bipartite_shape(self):
@@ -181,7 +181,7 @@ class Chi(QuantumChannel):
 
     @classmethod
     def _tensor(cls, a, b):
-        ret = copy.copy(a)
+        ret = _copy.copy(a)
         ret._op_shape = a._op_shape.tensor(b._op_shape)
         ret._data = np.kron(a._data, b.data)
         return ret
diff --git a/qiskit/quantum_info/operators/channel/choi.py b/qiskit/quantum_info/operators/channel/choi.py
index afd8e4fca6..1c9579e456 100644
--- a/qiskit/quantum_info/operators/channel/choi.py
+++ b/qiskit/quantum_info/operators/channel/choi.py
@@ -16,10 +16,11 @@ Choi-matrix representation of a Quantum Channel.
 """
 
 from __future__ import annotations
-import copy
+import copy as _copy
 import math
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
@@ -134,10 +135,9 @@ class Choi(QuantumChannel):
             choi_mat = _to_choi(rep, data._data, input_dim, output_dim)
         super().__init__(choi_mat, op_shape=op_shape)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     @property
     def _bipartite_shape(self):
@@ -152,12 +152,12 @@ class Choi(QuantumChannel):
     # ---------------------------------------------------------------------
 
     def conjugate(self):
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.conj(self._data)
         return ret
 
     def transpose(self):
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._op_shape = self._op_shape.transpose()
         # Make bipartite matrix
         d_in, d_out = self.dim
@@ -206,7 +206,7 @@ class Choi(QuantumChannel):
 
     @classmethod
     def _tensor(cls, a, b):
-        ret = copy.copy(a)
+        ret = _copy.copy(a)
         ret._op_shape = a._op_shape.tensor(b._op_shape)
         ret._data = _bipartite_tensor(
             a._data, b.data, shape1=a._bipartite_shape, shape2=b._bipartite_shape
diff --git a/qiskit/quantum_info/operators/channel/ptm.py b/qiskit/quantum_info/operators/channel/ptm.py
index 1bdf1b5ef2..84db071121 100644
--- a/qiskit/quantum_info/operators/channel/ptm.py
+++ b/qiskit/quantum_info/operators/channel/ptm.py
@@ -16,10 +16,11 @@ Pauli Transfer Matrix (PTM) representation of a Quantum Channel.
 """
 
 from __future__ import annotations
-import copy
+import copy as _copy
 import math
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
@@ -133,10 +134,9 @@ class PTM(QuantumChannel):
             raise QiskitError("Input is not an n-qubit Pauli transfer matrix.")
         super().__init__(ptm, num_qubits=num_qubits)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     @property
     def _bipartite_shape(self):
@@ -194,7 +194,7 @@ class PTM(QuantumChannel):
 
     @classmethod
     def _tensor(cls, a, b):
-        ret = copy.copy(a)
+        ret = _copy.copy(a)
         ret._op_shape = a._op_shape.tensor(b._op_shape)
         ret._data = np.kron(a._data, b.data)
         return ret
diff --git a/qiskit/quantum_info/operators/channel/superop.py b/qiskit/quantum_info/operators/channel/superop.py
index 0d7116ef50..19867696ec 100644
--- a/qiskit/quantum_info/operators/channel/superop.py
+++ b/qiskit/quantum_info/operators/channel/superop.py
@@ -15,12 +15,13 @@ Superoperator representation of a Quantum Channel."""
 
 from __future__ import annotations
 
-import copy
+import copy as _copy
 import math
 from typing import TYPE_CHECKING
 
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.exceptions import QiskitError
@@ -127,10 +128,9 @@ class SuperOp(QuantumChannel):
         # Initialize QuantumChannel
         super().__init__(super_mat, op_shape=op_shape)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     @property
     def _tensor_shape(self):
@@ -149,18 +149,18 @@ class SuperOp(QuantumChannel):
     # ---------------------------------------------------------------------
 
     def conjugate(self):
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.conj(self._data)
         return ret
 
     def transpose(self):
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.transpose(self._data)
         ret._op_shape = self._op_shape.transpose()
         return ret
 
     def adjoint(self):
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.conj(np.transpose(self._data))
         ret._op_shape = self._op_shape.transpose()
         return ret
@@ -177,7 +177,7 @@ class SuperOp(QuantumChannel):
 
     @classmethod
     def _tensor(cls, a, b):
-        ret = copy.copy(a)
+        ret = _copy.copy(a)
         ret._op_shape = a._op_shape.tensor(b._op_shape)
         ret._data = _bipartite_tensor(
             a._data, b.data, shape1=a._bipartite_shape, shape2=b._bipartite_shape
diff --git a/qiskit/quantum_info/operators/dihedral/dihedral.py b/qiskit/quantum_info/operators/dihedral/dihedral.py
index af15e0ed3a..4f49879063 100644
--- a/qiskit/quantum_info/operators/dihedral/dihedral.py
+++ b/qiskit/quantum_info/operators/dihedral/dihedral.py
@@ -357,10 +357,11 @@ class CNOTDihedral(BaseOperator, AdjointMixin):
         _append_circuit(elem, circuit)
         return elem
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.to_matrix(), dtype=dtype)
-        return self.to_matrix()
+    def __array__(self, dtype=None, copy=None):
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
+        arr = self.to_matrix()
+        return arr if dtype is None else arr.astype(dtype, copy=False)
 
     def to_matrix(self):
         """Convert operator to Numpy matrix."""
diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py
index 90d40b93cc..d119a38124 100644
--- a/qiskit/quantum_info/operators/operator.py
+++ b/qiskit/quantum_info/operators/operator.py
@@ -16,13 +16,14 @@ Matrix Operator class.
 
 from __future__ import annotations
 
-import copy
+import copy as _copy
 import re
 from numbers import Number
 from typing import TYPE_CHECKING
 
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.library.standard_gates import HGate, IGate, SGate, TGate, XGate, YGate, ZGate
 from qiskit.circuit.operation import Operation
@@ -117,10 +118,9 @@ class Operator(LinearOp):
             shape=self._data.shape,
         )
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     def __repr__(self):
         prefix = "Operator("
@@ -447,13 +447,13 @@ class Operator(LinearOp):
 
     def conjugate(self):
         # Make a shallow copy and update array
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.conj(self._data)
         return ret
 
     def transpose(self):
         # Make a shallow copy and update array
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.transpose(self._data)
         ret._op_shape = self._op_shape.transpose()
         return ret
@@ -523,7 +523,7 @@ class Operator(LinearOp):
         """
         if self.input_dims() != self.output_dims():
             raise QiskitError("Can only power with input_dims = output_dims.")
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         if isinstance(n, int):
             ret._data = np.linalg.matrix_power(self.data, n)
         else:
@@ -550,7 +550,7 @@ class Operator(LinearOp):
 
     @classmethod
     def _tensor(cls, a, b):
-        ret = copy.copy(a)
+        ret = _copy.copy(a)
         ret._op_shape = a._op_shape.tensor(b._op_shape)
         ret._data = np.kron(a.data, b.data)
         return ret
@@ -585,7 +585,7 @@ class Operator(LinearOp):
         self._op_shape._validate_add(other._op_shape, qargs)
         other = ScalarOp._pad_with_identity(self, other, qargs)
 
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = self.data + other.data
         return ret
 
@@ -603,7 +603,7 @@ class Operator(LinearOp):
         """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = other * self._data
         return ret
 
@@ -643,7 +643,7 @@ class Operator(LinearOp):
         Returns:
             Operator: the operator with reversed subsystem order.
         """
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         axes = tuple(range(self._op_shape._num_qargs_l - 1, -1, -1))
         axes = axes + tuple(len(axes) + i for i in axes)
         ret._data = np.reshape(
diff --git a/qiskit/quantum_info/operators/scalar_op.py b/qiskit/quantum_info/operators/scalar_op.py
index d856a39c2a..38f3619373 100644
--- a/qiskit/quantum_info/operators/scalar_op.py
+++ b/qiskit/quantum_info/operators/scalar_op.py
@@ -15,7 +15,7 @@ ScalarOp class
 """
 
 from __future__ import annotations
-import copy
+import copy as _copy
 from numbers import Number
 import numpy as np
 
@@ -52,10 +52,11 @@ class ScalarOp(LinearOp):
         self._coeff = coeff
         super().__init__(input_dims=dims, output_dims=dims)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.to_matrix(), dtype=dtype)
-        return self.to_matrix()
+    def __array__(self, dtype=None, copy=None):
+        if copy is False:
+            raise ValueError("could not produce matrix without calculation")
+        arr = self.to_matrix()
+        return arr if dtype is None else arr.astype(dtype, copy=False)
 
     def __repr__(self):
         return f"ScalarOp({self.input_dims()}, coeff={self.coeff})"
@@ -104,7 +105,7 @@ class ScalarOp(LinearOp):
         # If other is also an ScalarOp we only need to
         # update the coefficient and dimensions
         if isinstance(other, ScalarOp):
-            ret = copy.copy(self)
+            ret = _copy.copy(self)
             ret._coeff = self.coeff * other.coeff
             ret._op_shape = new_shape
             return ret
@@ -112,7 +113,7 @@ class ScalarOp(LinearOp):
         # If we are composing on the full system we return the
         # other operator with reshaped dimensions
         if qargs is None:
-            ret = copy.copy(other)
+            ret = _copy.copy(other)
             ret._op_shape = new_shape
             # Other operator might not support scalar multiplication
             # so we treat the identity as a special case to avoid a
@@ -148,7 +149,7 @@ class ScalarOp(LinearOp):
             other = Operator(other)
 
         if isinstance(other, ScalarOp):
-            ret = copy.copy(self)
+            ret = _copy.copy(self)
             ret._coeff = self.coeff * other.coeff
             ret._op_shape = self._op_shape.tensor(other._op_shape)
             return ret
@@ -160,7 +161,7 @@ class ScalarOp(LinearOp):
             other = Operator(other)
 
         if isinstance(other, ScalarOp):
-            ret = copy.copy(self)
+            ret = _copy.copy(self)
             ret._coeff = self.coeff * other.coeff
             ret._op_shape = self._op_shape.expand(other._op_shape)
             return ret
diff --git a/qiskit/quantum_info/operators/symplectic/clifford.py b/qiskit/quantum_info/operators/symplectic/clifford.py
index f0d30cf59f..47ef8b629d 100644
--- a/qiskit/quantum_info/operators/symplectic/clifford.py
+++ b/qiskit/quantum_info/operators/symplectic/clifford.py
@@ -122,10 +122,11 @@ class Clifford(BaseOperator, AdjointMixin, Operation):
     _COMPOSE_PHASE_LOOKUP = None
     _COMPOSE_1Q_LOOKUP = None
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.to_matrix(), dtype=dtype)
-        return self.to_matrix()
+    def __array__(self, dtype=None, copy=None):
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
+        arr = self.to_matrix()
+        return arr if dtype is None else arr.astype(dtype, copy=False)
 
     def __init__(self, data, validate=True, copy=True):
         """Initialize an operator object."""
@@ -164,8 +165,17 @@ class Clifford(BaseOperator, AdjointMixin, Operation):
 
         # Initialize StabilizerTable directly from the data
         else:
-            if isinstance(data, (list, np.ndarray)) and np.asarray(data, dtype=bool).ndim == 2:
-                data = np.array(data, dtype=bool, copy=copy)
+            if (
+                isinstance(data, (list, np.ndarray))
+                and (data_asarray := np.asarray(data, dtype=bool)).ndim == 2
+            ):
+                # This little dance is to avoid Numpy 1/2 incompatiblities between the availability
+                # and meaning of the 'copy' argument in 'array' and 'asarray', when the input needs
+                # its dtype converting.  'asarray' prefers to return 'self' if possible in both.
+                if copy and np.may_share_memory(data, data_asarray):
+                    data = data_asarray.copy()
+                else:
+                    data = data_asarray
                 if data.shape[0] == data.shape[1]:
                     self.tableau = self._stack_table_phase(
                         data, np.zeros(data.shape[0], dtype=bool)
diff --git a/qiskit/quantum_info/operators/symplectic/pauli.py b/qiskit/quantum_info/operators/symplectic/pauli.py
index 1bdff0cf8f..e1bcfa29eb 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli.py
@@ -222,10 +222,11 @@ class Pauli(BasePauli):
             return front + "..."
         return self.to_label()
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.to_matrix(), dtype=dtype)
-        return self.to_matrix()
+    def __array__(self, dtype=None, copy=None):
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
+        arr = self.to_matrix()
+        return arr if dtype is None else arr.astype(dtype, copy=False)
 
     @classmethod
     def set_truncation(cls, val: int):
diff --git a/qiskit/quantum_info/operators/symplectic/pauli_list.py b/qiskit/quantum_info/operators/symplectic/pauli_list.py
index 1bb9ae3d1a..3d348d2363 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli_list.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli_list.py
@@ -148,14 +148,15 @@ class PauliList(BasePauli, LinearMixin, GroupMixin):
         """Return settings."""
         return {"data": self.to_labels()}
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Convert to numpy array"""
-        # pylint: disable=unused-argument
+        if copy is False:
+            raise ValueError("cannot provide a matrix without calculation")
         shape = (len(self),) + 2 * (2**self.num_qubits,)
         ret = np.zeros(shape, dtype=complex)
         for i, mat in enumerate(self.matrix_iter()):
             ret[i] = mat
-        return ret
+        return ret if dtype is None else ret.astype(dtype, copy=False)
 
     @staticmethod
     def _from_paulis(data):
diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 0d29cd098d..dc445509b0 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -142,7 +142,12 @@ class SparsePauliOp(LinearOp):
         if coeffs is None:
             coeffs = np.ones(pauli_list.size, dtype=complex)
         else:
-            coeffs = np.array(coeffs, copy=copy, dtype=dtype)
+            coeffs_asarray = np.asarray(coeffs, dtype=dtype)
+            coeffs = (
+                coeffs_asarray.copy()
+                if copy and np.may_share_memory(coeffs, coeffs_asarray)
+                else coeffs_asarray
+            )
 
         if ignore_pauli_phase:
             # Fast path used in copy operations, where the phase of the PauliList is already known
@@ -166,10 +171,11 @@ class SparsePauliOp(LinearOp):
         # Initialize LinearOp
         super().__init__(num_qubits=self._pauli_list.num_qubits)
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.to_matrix(), dtype=dtype)
-        return self.to_matrix()
+    def __array__(self, dtype=None, copy=None):
+        if copy is False:
+            raise ValueError("unable to avoid copy while creating an array as requested")
+        arr = self.to_matrix()
+        return arr if dtype is None else arr.astype(dtype, copy=False)
 
     def __repr__(self):
         prefix = "SparsePauliOp("
diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index 07cc656857..1c66d8bcf5 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -15,10 +15,11 @@ DensityMatrix quantum state class.
 """
 
 from __future__ import annotations
-import copy
+import copy as _copy
 from numbers import Number
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
@@ -110,10 +111,9 @@ class DensityMatrix(QuantumState, TolerancesMixin):
             raise QiskitError("Invalid DensityMatrix input: not a square matrix.")
         super().__init__(op_shape=OpShape.auto(shape=self._data.shape, dims_l=dims, dims_r=dims))
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     def __eq__(self, other):
         return super().__eq__(other) and np.allclose(
@@ -241,7 +241,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, DensityMatrix):
             other = DensityMatrix(other)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.kron(self._data, other._data)
         ret._op_shape = self._op_shape.tensor(other._op_shape)
         return ret
@@ -260,7 +260,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, DensityMatrix):
             other = DensityMatrix(other)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = np.kron(other._data, self._data)
         ret._op_shape = self._op_shape.expand(other._op_shape)
         return ret
@@ -281,7 +281,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         if not isinstance(other, DensityMatrix):
             other = DensityMatrix(other)
         self._op_shape._validate_add(other._op_shape)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = self.data + other.data
         return ret
 
@@ -299,7 +299,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = other * self.data
         return ret
 
@@ -356,7 +356,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         Returns:
             DensityMatrix: the state with reversed subsystem order.
         """
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         axes = tuple(range(self._op_shape._num_qargs_l - 1, -1, -1))
         axes = axes + tuple(len(axes) + i for i in axes)
         ret._data = np.reshape(
@@ -523,7 +523,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         """
         if qargs is None:
             # Resetting all qubits does not require sampling or RNG
-            ret = copy.copy(self)
+            ret = _copy.copy(self)
             state = np.zeros(self._op_shape.shape, dtype=complex)
             state[0, 0] = 1
             ret._data = state
@@ -715,7 +715,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         new_shape._dims_r = new_shape._dims_l
         new_shape._num_qargs_r = new_shape._num_qargs_l
 
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         if qargs is None:
             # Evolution on full matrix
             op_mat = other.data
@@ -792,7 +792,7 @@ class DensityMatrix(QuantumState, TolerancesMixin):
         """Return a new statevector by applying an instruction."""
         if isinstance(obj, QuantumCircuit):
             obj = obj.to_instruction()
-        vec = copy.copy(self)
+        vec = _copy.copy(self)
         vec._append_instruction(obj, qargs=qargs)
         return vec
 
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index b13ccde217..df39ba42f9 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -14,13 +14,14 @@
 Statevector quantum state class.
 """
 from __future__ import annotations
-import copy
+import copy as _copy
 import math
 import re
 from numbers import Number
 
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
@@ -104,10 +105,9 @@ class Statevector(QuantumState, TolerancesMixin):
                 raise QiskitError("Invalid input: not a vector or column-vector.")
         super().__init__(op_shape=OpShape.auto(shape=shape, dims_l=dims, num_qubits_r=0))
 
-    def __array__(self, dtype=None):
-        if dtype:
-            return np.asarray(self.data, dtype=dtype)
-        return self.data
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = self.data.dtype if dtype is None else dtype
+        return np.array(self.data, dtype=dtype, copy=copy)
 
     def __eq__(self, other):
         return super().__eq__(other) and np.allclose(
@@ -277,7 +277,7 @@ class Statevector(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, Statevector):
             other = Statevector(other)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._op_shape = self._op_shape.tensor(other._op_shape)
         ret._data = np.kron(self._data, other._data)
         return ret
@@ -318,7 +318,7 @@ class Statevector(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, Statevector):
             other = Statevector(other)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._op_shape = self._op_shape.expand(other._op_shape)
         ret._data = np.kron(other._data, self._data)
         return ret
@@ -339,7 +339,7 @@ class Statevector(QuantumState, TolerancesMixin):
         if not isinstance(other, Statevector):
             other = Statevector(other)
         self._op_shape._validate_add(other._op_shape)
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = self.data + other.data
         return ret
 
@@ -357,7 +357,7 @@ class Statevector(QuantumState, TolerancesMixin):
         """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         ret._data = other * self.data
         return ret
 
@@ -382,7 +382,7 @@ class Statevector(QuantumState, TolerancesMixin):
             qargs = getattr(other, "qargs", None)
 
         # Get return vector
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
 
         # Evolution by a circuit or instruction
         if isinstance(other, QuantumCircuit):
@@ -448,7 +448,7 @@ class Statevector(QuantumState, TolerancesMixin):
         Returns:
             Statevector: the Statevector with reversed subsystem order.
         """
-        ret = copy.copy(self)
+        ret = _copy.copy(self)
         axes = tuple(range(self._op_shape._num_qargs_l - 1, -1, -1))
         ret._data = np.reshape(
             np.transpose(np.reshape(self.data, self._op_shape.tensor_shape), axes),
@@ -619,7 +619,7 @@ class Statevector(QuantumState, TolerancesMixin):
         """
         if qargs is None:
             # Resetting all qubits does not require sampling or RNG
-            ret = copy.copy(self)
+            ret = _copy.copy(self)
             state = np.zeros(self._op_shape.shape, dtype=complex)
             state[0] = 1
             ret._data = state
diff --git a/qiskit/visualization/array.py b/qiskit/visualization/array.py
index b076e38b17..3a8ef29171 100644
--- a/qiskit/visualization/array.py
+++ b/qiskit/visualization/array.py
@@ -33,7 +33,7 @@ def _num_to_latex(raw_value, decimals=15, first_term=True, coefficient=False):
     """
     import sympy  # runtime import
 
-    raw_value = np.around(raw_value, decimals=decimals)
+    raw_value = np.around(raw_value, decimals=decimals).item()
     value = sympy.nsimplify(raw_value, rational=False)
 
     if isinstance(value, sympy.core.numbers.Rational) and value.denominator > 50:
diff --git a/releasenotes/notes/numpy-2.0-2f3e35bd42c48518.yaml b/releasenotes/notes/numpy-2.0-2f3e35bd42c48518.yaml
new file mode 100644
index 0000000000..3595f2f936
--- /dev/null
+++ b/releasenotes/notes/numpy-2.0-2f3e35bd42c48518.yaml
@@ -0,0 +1,5 @@
+---
+features:
+  - |
+    This release of Qiskit finalizes support for NumPy 2.0.  Qiskit will continue to support both
+    Numpy 1.x and 2.x for the foreseeable future.
diff --git a/requirements.txt b/requirements.txt
index 0a4f0f32c5..539f958799 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,9 +1,9 @@
 rustworkx>=0.14.0
-numpy>=1.17,<2
+numpy>=1.17,<3
 scipy>=1.5
 sympy>=1.3
 dill>=0.3
 python-dateutil>=2.8.0
 stevedore>=3.0.0
 typing-extensions
-symengine>=0.11
\ No newline at end of file
+symengine>=0.11
diff --git a/test/python/providers/test_backend_v2.py b/test/python/providers/test_backend_v2.py
index 6b974df902..70330085b1 100644
--- a/test/python/providers/test_backend_v2.py
+++ b/test/python/providers/test_backend_v2.py
@@ -19,6 +19,7 @@ from test import combine
 
 from ddt import ddt, data
 
+from numpy.testing import assert_array_max_ulp
 from qiskit.circuit import QuantumCircuit, ClassicalRegister, QuantumRegister
 from qiskit.circuit.library.standard_gates import (
     CXGate,
@@ -66,9 +67,11 @@ class TestBackendV2(QiskitTestCase):
     def test_qubit_properties(self):
         """Test that qubit properties are returned as expected."""
         props = self.backend.qubit_properties([1, 0])
-        self.assertEqual([0.0001697368029059364, 0.00017739560485559633], [x.t1 for x in props])
-        self.assertEqual([0.00010941773478876496, 0.00014388784397520525], [x.t2 for x in props])
-        self.assertEqual([5487811175.818378, 5429298959.955691], [x.frequency for x in props])
+        assert_array_max_ulp([0.0001697368029059364, 0.00017739560485559633], [x.t1 for x in props])
+        assert_array_max_ulp(
+            [0.00010941773478876496, 0.00014388784397520525], [x.t2 for x in props]
+        )
+        assert_array_max_ulp([5487811175.818378, 5429298959.955691], [x.frequency for x in props])
 
     def test_legacy_qubit_properties(self):
         """Test that qubit props work for backends not using properties in target."""
@@ -82,9 +85,11 @@ class TestBackendV2(QiskitTestCase):
                 return [self.target.qubit_properties[i] for i in qubit]
 
         props = FakeBackendV2LegacyQubitProps(num_qubits=2, seed=42).qubit_properties([1, 0])
-        self.assertEqual([0.0001697368029059364, 0.00017739560485559633], [x.t1 for x in props])
-        self.assertEqual([0.00010941773478876496, 0.00014388784397520525], [x.t2 for x in props])
-        self.assertEqual([5487811175.818378, 5429298959.955691], [x.frequency for x in props])
+        assert_array_max_ulp([0.0001697368029059364, 0.00017739560485559633], [x.t1 for x in props])
+        assert_array_max_ulp(
+            [0.00010941773478876496, 0.00014388784397520525], [x.t2 for x in props]
+        )
+        assert_array_max_ulp([5487811175.818378, 5429298959.955691], [x.frequency for x in props])
 
     def test_no_qubit_properties_raises(self):
         """Ensure that if no qubit properties are defined we raise correctly."""
diff --git a/test/python/transpiler/test_dynamical_decoupling.py b/test/python/transpiler/test_dynamical_decoupling.py
index 6460847b2e..5f11ede3a3 100644
--- a/test/python/transpiler/test_dynamical_decoupling.py
+++ b/test/python/transpiler/test_dynamical_decoupling.py
@@ -447,7 +447,9 @@ class TestPadDynamicalDecoupling(QiskitTestCase):
             representation to satisfy PadDynamicalDecoupling's check.
             """
 
-            def __array__(self, dtype=None):
+            def __array__(self, dtype=None, copy=None):
+                if copy is False:
+                    raise ValueError("cannot produce matrix without calculation")
                 return np.eye(2, dtype=dtype)
 
         # A gate with one unbound and one bound parameter to leave in the final