diff --git a/qiskit_ionq/__init__.py b/qiskit_ionq/__init__.py
index 08af3927..e0fd426f 100644
--- a/qiskit_ionq/__init__.py
+++ b/qiskit_ionq/__init__.py
@@ -28,5 +28,5 @@
 
 from .ionq_provider import IonQProvider
 from .version import __version__
-from .ionq_gates import GPIGate, GPI2Gate, MSGate, ZZGate
+from .ionq_gates import GPIGate, GPI2Gate, MSGate, VirtualZGate, ZZGate
 from .constants import ErrorMitigation
diff --git a/qiskit_ionq/helpers.py b/qiskit_ionq/helpers.py
index cbb6f87c..8426b235 100644
--- a/qiskit_ionq/helpers.py
+++ b/qiskit_ionq/helpers.py
@@ -118,6 +118,7 @@
     "gpi",  # TODO All single qubit gates can transpile into GPI/GPI2
     "gpi2",
     "ms",  # Pairwise MS gate
+    "gz",  # VirtualZ gate (RZ)
     "zz",  # ZZ gate
 ]
 
diff --git a/qiskit_ionq/ionq_gates.py b/qiskit_ionq/ionq_gates.py
index 3f7335f9..0819ba5a 100644
--- a/qiskit_ionq/ionq_gates.py
+++ b/qiskit_ionq/ionq_gates.py
@@ -27,8 +27,6 @@
 """Native gateset for IonQ hardware."""
 
 from typing import Optional
-import cmath
-import math
 import numpy
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.parameterexpression import ParameterValueType
@@ -58,8 +56,8 @@ def __init__(self, phi: ParameterValueType, label: Optional[str] = None):
 
     def __array__(self, dtype=None):
         """Return a numpy.array for the GPI gate."""
-        top = cmath.exp(-self.params[0] * 2 * math.pi * 1j)
-        bot = cmath.exp(self.params[0] * 2 * math.pi * 1j)
+        top = numpy.exp(-self.params[0] * 2 * numpy.pi * 1j)
+        bot = numpy.exp(self.params[0] * 2 * numpy.pi * 1j)
         return numpy.array([[0, top], [bot, 0]], dtype=dtype)
 
 
@@ -87,9 +85,42 @@ def __init__(self, phi: ParameterValueType, label: Optional[str] = None):
 
     def __array__(self, dtype=None):
         """Return a numpy.array for the GPI2 gate."""
-        top = -1j * cmath.exp(self.params[0] * 2 * math.pi * -1j)
-        bot = -1j * cmath.exp(self.params[0] * 2 * math.pi * 1j)
-        return numpy.array([[1, top], [bot, 1]], dtype=dtype) / math.sqrt(2)
+        top = -1j * numpy.exp(self.params[0] * 2 * numpy.pi * -1j)
+        bot = -1j * numpy.exp(self.params[0] * 2 * numpy.pi * 1j)
+        return numpy.array([[1, top], [bot, 1]], dtype=dtype) / numpy.sqrt(2)
+
+
+class VirtualZGate(Gate):
+    r"""Single-qubit rotation about the Z axis.
+    **Circuit symbol:**
+    .. parsed-literal::
+             ┌───────┐
+        q_0: ┤ GZ(ϴ) ├
+             └───────┘
+    **Matrix Representation:**
+    .. math::
+
+        GZ(\theta) = \exp\left(-i\frac{\theta}{2}Z\right) =
+            \begin{pmatrix}
+                e^{-i\frac{\theta}{2}} & 0 \\
+                0 & e^{i\frac{\theta}{2}}
+            \end{pmatrix}
+    """
+
+    def __init__(
+        self,
+        theta: ParameterValueType,
+        label: Optional[str] = None,
+    ):
+        """Create new GZ gate."""
+        super().__init__("gz", 1, [theta], label=label)
+
+    def __array__(self, dtype=None):
+        """Return a numpy.array for the GZ gate."""
+        theta = self.params[0]
+        top = numpy.exp(-1j * theta / 2)
+        bot = numpy.exp(1j * theta / 2)
+        return numpy.array([[top, 0], [0, bot]], dtype=dtype)
 
 
 class MSGate(Gate):
@@ -134,15 +165,15 @@ def __array__(self, dtype=None):
         phi0 = self.params[0]
         phi1 = self.params[1]
         theta = self.params[2]
-        diag = numpy.cos(math.pi * theta)
-        sin = numpy.sin(math.pi * theta)
+        diag = numpy.cos(numpy.pi * theta)
+        sin = numpy.sin(numpy.pi * theta)
 
         return numpy.array(
             [
-                [diag, 0, 0, sin * -1j * cmath.exp(-1j * 2 * math.pi * (phi0 + phi1))],
-                [0, diag, sin * -1j * cmath.exp(-1j * 2 * math.pi * (phi0 - phi1)), 0],
-                [0, sin * -1j * cmath.exp(1j * 2 * math.pi * (phi0 - phi1)), diag, 0],
-                [sin * -1j * cmath.exp(1j * 2 * math.pi * (phi0 + phi1)), 0, 0, diag],
+                [diag, 0, 0, sin * -1j * numpy.exp(-1j * 2 * numpy.pi * (phi0 + phi1))],
+                [0, diag, sin * -1j * numpy.exp(-1j * 2 * numpy.pi * (phi0 - phi1)), 0],
+                [0, sin * -1j * numpy.exp(1j * 2 * numpy.pi * (phi0 - phi1)), diag, 0],
+                [sin * -1j * numpy.exp(1j * 2 * numpy.pi * (phi0 + phi1)), 0, 0, diag],
             ],
             dtype=dtype,
         )
@@ -174,13 +205,13 @@ def __init__(self, theta: ParameterValueType, label: Optional[str] = None):
 
     def __array__(self, dtype=None):
         """Return a numpy.array for the ZZ gate."""
-        itheta2 = 1j * float(self.params[0]) * math.pi
+        itheta2 = 1j * float(self.params[0]) * numpy.pi
         return numpy.array(
             [
-                [cmath.exp(-itheta2), 0, 0, 0],
-                [0, cmath.exp(itheta2), 0, 0],
-                [0, 0, cmath.exp(itheta2), 0],
-                [0, 0, 0, cmath.exp(-itheta2)],
+                [numpy.exp(-itheta2), 0, 0, 0],
+                [0, numpy.exp(itheta2), 0, 0],
+                [0, 0, numpy.exp(itheta2), 0],
+                [0, 0, 0, numpy.exp(-itheta2)],
             ],
             dtype=dtype,
         )
diff --git a/test/helpers/test_qiskit_to_ionq.py b/test/helpers/test_qiskit_to_ionq.py
index 38c4c558..dc6a5e6a 100644
--- a/test/helpers/test_qiskit_to_ionq.py
+++ b/test/helpers/test_qiskit_to_ionq.py
@@ -35,7 +35,7 @@
 
 from qiskit_ionq.exceptions import IonQGateError
 from qiskit_ionq.helpers import qiskit_to_ionq, decompress_metadata_string_to_dict
-from qiskit_ionq.ionq_gates import GPIGate, GPI2Gate, MSGate, ZZGate
+from qiskit_ionq.ionq_gates import GPIGate, GPI2Gate, MSGate, VirtualZGate, ZZGate
 from qiskit_ionq.constants import ErrorMitigation
 
 
@@ -273,6 +273,7 @@ def test_native_circuit_transpile(simulator_backend):
     circ = QuantumCircuit(3, name="blame_test")
     circ.append(GPIGate(0.1), [0])
     circ.append(GPI2Gate(0.2), [1])
+    circ.append(VirtualZGate(0.3), [1])
     circ.append(MSGate(0.2, 0.3, 0.25), [1, 2])
     circ.append(ZZGate(0.4), [0, 2])
 
@@ -291,6 +292,7 @@ def test_full_native_circuit(simulator_backend):
     qc = QuantumCircuit(3, name="blame_test")
     qc.append(GPIGate(0.1), [0])
     qc.append(GPI2Gate(0.2), [1])
+    qc.append(VirtualZGate(0.3), [1])
     qc.append(MSGate(0.2, 0.3, 0.25), [1, 2])
     qc.append(ZZGate(0.4), [0, 2])
     ionq_json = qiskit_to_ionq(
@@ -325,6 +327,7 @@ def test_full_native_circuit(simulator_backend):
             "circuit": [
                 {"gate": "gpi", "target": 0, "phase": 0.1},
                 {"gate": "gpi2", "target": 1, "phase": 0.2},
+                {"gate": "gz", "targets": [1], "phase": 0.3},
                 {"gate": "ms", "targets": [1, 2], "phases": [0.2, 0.3], "angle": 0.25},
                 {"gate": "zz", "angle": 0.4, "targets": [0, 2]},
             ],
diff --git a/test/ionq_gates/test_gates.py b/test/ionq_gates/test_gates.py
index f3465a9b..9292e507 100644
--- a/test/ionq_gates/test_gates.py
+++ b/test/ionq_gates/test_gates.py
@@ -24,16 +24,15 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-"""Tests for the IonQ's GPIGate, GPI2Gate, MSGate, ZZGate."""
+"""Tests for the IonQ's GPIGate, GPI2Gate, MSGate, VirtualZGate, ZZGate."""
 # pylint: disable=redefined-outer-name
 
-import math
 import numpy
 
 import pytest
 
-from qiskit.circuit.library import XGate, YGate, RXGate, RYGate
-from qiskit_ionq import GPIGate, GPI2Gate, MSGate, ZZGate
+from qiskit.circuit.library import XGate, YGate, RXGate, RYGate, RZGate
+from qiskit_ionq import GPIGate, GPI2Gate, MSGate, VirtualZGate, ZZGate
 
 
 @pytest.mark.parametrize("gate,phase", [(XGate(), 0), (YGate(), 0.25)])
@@ -44,7 +43,7 @@ def test_gpi_equivalences(gate, phase):
 
 
 @pytest.mark.parametrize(
-    "gate,phase", [(RXGate(math.pi / 2), 1), (RYGate(math.pi / 2), 0.25)]
+    "gate,phase", [(RXGate(numpy.pi / 2), 1), (RYGate(numpy.pi / 2), 0.25)]
 )
 def test_gpi2_equivalences(gate, phase):
     """Tests equivalence of the GPI2 gate at specific phases."""
@@ -52,7 +51,16 @@ def test_gpi2_equivalences(gate, phase):
     numpy.testing.assert_array_almost_equal(numpy.array(gate), numpy.array(gpi))
 
 
-@pytest.mark.parametrize("phase", [0, 0.1, 0.4, math.pi / 2, math.pi, 2 * math.pi])
+@pytest.mark.parametrize("gate,angle", [(RZGate, 0), (RZGate, numpy.pi)])
+def test_virtualz_equivalence(gate, angle):
+    """Tests equivalence of the VirtualZGate gate."""
+    virtualz = VirtualZGate(angle)
+    numpy.testing.assert_array_almost_equal(
+        gate(angle).to_matrix(), virtualz.to_matrix()
+    )
+
+
+@pytest.mark.parametrize("phase", [0, 0.1, 0.4, numpy.pi / 2, numpy.pi, 2 * numpy.pi])
 def test_gpi_inverse(phase):
     """Tests that the GPI gate is unitary."""
     gate = GPIGate(phase)
@@ -60,7 +68,7 @@ def test_gpi_inverse(phase):
     numpy.testing.assert_array_almost_equal(mat.dot(mat.conj().T), numpy.identity(2))
 
 
-@pytest.mark.parametrize("phase", [0, 0.1, 0.4, math.pi / 2, math.pi, 2 * math.pi])
+@pytest.mark.parametrize("phase", [0, 0.1, 0.4, numpy.pi / 2, numpy.pi, 2 * numpy.pi])
 def test_gpi2_inverse(phase):
     """Tests that the GPI2 gate is unitary."""
     gate = GPI2Gate(phase)
@@ -75,9 +83,9 @@ def test_gpi2_inverse(phase):
         (0, 1, 0.25),
         (0.1, 1, 0.25),
         (0.4, 1, 0.25),
-        (math.pi / 2, 0, 0.25),
-        (0, math.pi, 0.25),
-        (0.1, 2 * math.pi, 0.25),
+        (numpy.pi / 2, 0, 0.25),
+        (0, numpy.pi, 0.25),
+        (0.1, 2 * numpy.pi, 0.25),
     ],
 )
 def test_ms_inverse(params):
@@ -90,7 +98,7 @@ def test_ms_inverse(params):
 
 @pytest.mark.parametrize(
     "angle",
-    [0, 0.1, 0.4, math.pi / 2, math.pi, 2 * math.pi],
+    [0, 0.1, 0.4, numpy.pi / 2, numpy.pi, 2 * numpy.pi],
 )
 def test_zz_inverse(angle):
     """Tests that the ZZ gate is unitary."""