diff --git a/README.md b/README.md
index bb7e4f15..407aa960 100644
--- a/README.md
+++ b/README.md
@@ -498,6 +498,8 @@ The (current) list of properties and methods is:
    * `tau` (`M`, `mass`): see note above
    * `tau2` (`M2`, `mass2`)
    * `beta`: scalar(s) between $0$ (inclusive) and $1$ (exclusive, unless the vector components are infinite)
+   * `deltaRapidityPhi`: $Delta R_{\mbox{rapidity}} = \Delta\phi^2 + \Delta \mbox{rapidity}^2$
+   * `deltaRapidityPhi2`: the above, squared
    * `gamma`: scalar(s) between $1$ (inclusive) and $\infty$
    * `rapidity`: scalar(s) between $0$ (inclusive) and $\infty$
    * `boost_p4(four_vector)`: change coordinate system using another 4D vector as the difference
diff --git a/docs/api/vector._compute.lorentz.rst b/docs/api/vector._compute.lorentz.rst
index 701674d2..49238826 100644
--- a/docs/api/vector._compute.lorentz.rst
+++ b/docs/api/vector._compute.lorentz.rst
@@ -130,6 +130,24 @@ vector.\_compute.lorentz.boost\_p4 module
    :show-inheritance:
    :private-members:
 
+vector.\_compute.lorentz.deltaRapidityPhi module
+--------------------------------------
+
+.. automodule:: vector._compute.lorentz.deltaRapidityPhi
+   :members:
+   :undoc-members:
+   :show-inheritance:
+   :private-members:
+
+vector.\_compute.lorentz.deltaRapidityPhi2 module
+---------------------------------------
+
+.. automodule:: vector._compute.lorentz.deltaRapidityPhi2
+   :members:
+   :undoc-members:
+   :show-inheritance:
+   :private-members:
+
 vector.\_compute.lorentz.dot module
 -----------------------------------
 
diff --git a/docs/usage/intro.ipynb b/docs/usage/intro.ipynb
index cd8304c0..9a852a8c 100644
--- a/docs/usage/intro.ipynb
+++ b/docs/usage/intro.ipynb
@@ -2310,6 +2310,8 @@
     "   * `tau` (`M`, `mass`): see note above\n",
     "   * `tau2` (`M2`, `mass2`)\n",
     "   * `beta`: scalar(s) between $0$ (inclusive) and $1$ (exclusive, unless the vector components are infinite)\n",
+    "   * `deltaRapidityPhi`: $Delta R_{\mbox{rapidity}} = \Delta\phi^2 + \Delta \mbox{rapidity}^2$\n",
+    "   * `deltaRapidityPhi2`: the above, squared\n",
     "   * `gamma`: scalar(s) between $1$ (inclusive) and $\\infty$\n",
     "   * `rapidity`: scalar(s) between $0$ (inclusive) and $\\infty$\n",
     "   * `boost_p4(four_vector)`: change coordinate system using another 4D vector as the difference\n",
diff --git a/tests/compute/lorentz/test_deltaRapidityPhi.py b/tests/compute/lorentz/test_deltaRapidityPhi.py
new file mode 100644
index 00000000..fb4c1633
--- /dev/null
+++ b/tests/compute/lorentz/test_deltaRapidityPhi.py
@@ -0,0 +1,120 @@
+# Copyright (c) 2019-2021, Jonas Eschle, Jim Pivarski, Eduardo Rodrigues, and Henry Schreiner.
+#
+# Distributed under the 3-clause BSD license, see accompanying file LICENSE
+# or https://github.com/scikit-hep/vector for details.
+
+import numpy
+import pytest
+
+import vector._backends.numpy_
+import vector._backends.object_
+
+
+def test_lorentz_object():
+    v1 = vector._backends.object_.MomentumObject4D(
+        vector._backends.object_.AzimuthalObjectXY(1.0, 1.0),
+        vector._backends.object_.LongitudinalObjectZ(1.0),
+        vector._backends.object_.TemporalObjectTau(1.0),
+    )
+    v2 = vector._backends.object_.MomentumObject4D(
+        vector._backends.object_.AzimuthalObjectXY(-1.0, -1.0),
+        vector._backends.object_.LongitudinalObjectZ(-1.0),
+        vector._backends.object_.TemporalObjectTau(1.0),
+    )
+    expected_result = numpy.sqrt(
+        # phi
+        numpy.pi**2
+        # rapidity
+        + ((0.5 * numpy.log(3 / 1) - 0.5 * numpy.log(1 / 3)) ** 2)
+    )
+    assert v1.deltaRapidityPhi(v2) == pytest.approx(expected_result)
+
+    for t1 in (
+        "xyzt",
+        "xythetat",
+        "xyetat",
+        "rhophizt",
+        "rhophithetat",
+        "rhophietat",
+        "xyztau",
+        "xythetatau",
+        "xyetatau",
+        "rhophiztau",
+        "rhophithetatau",
+        "rhophietatau",
+    ):
+        for t2 in (
+            "xyzt",
+            "xythetat",
+            "xyetat",
+            "rhophizt",
+            "rhophithetat",
+            "rhophietat",
+            "xyztau",
+            "xythetatau",
+            "xyetatau",
+            "rhophiztau",
+            "rhophithetatau",
+            "rhophietatau",
+        ):
+            tr1, tr2 = getattr(v1, "to_" + t1)(), getattr(v2, "to_" + t2)()
+            assert tr1.deltaRapidityPhi(tr2) == pytest.approx(expected_result)
+
+
+def test_lorentz_numpy():
+    v1 = vector._backends.numpy_.VectorNumpy4D(
+        [(1.0, 1.0, 1.0, 1.0)],
+        dtype=[
+            ("x", numpy.float64),
+            ("y", numpy.float64),
+            ("z", numpy.float64),
+            ("tau", numpy.float64),
+        ],
+    )
+    v2 = vector._backends.numpy_.VectorNumpy4D(
+        [(-1.0, -1.0, -1.0, 1.0)],
+        dtype=[
+            ("x", numpy.float64),
+            ("y", numpy.float64),
+            ("z", numpy.float64),
+            ("tau", numpy.float64),
+        ],
+    )
+    expected_result = numpy.sqrt(
+        # phi
+        numpy.pi**2
+        # rapidity
+        + ((0.5 * numpy.log(3 / 1) - 0.5 * numpy.log(1 / 3)) ** 2)
+    )
+    assert v1.deltaRapidityPhi(v2) == pytest.approx(expected_result)
+
+    for t1 in (
+        "xyzt",
+        "xythetat",
+        "xyetat",
+        "rhophizt",
+        "rhophithetat",
+        "rhophietat",
+        "xyztau",
+        "xythetatau",
+        "xyetatau",
+        "rhophiztau",
+        "rhophithetatau",
+        "rhophietatau",
+    ):
+        for t2 in (
+            "xyzt",
+            "xythetat",
+            "xyetat",
+            "rhophizt",
+            "rhophithetat",
+            "rhophietat",
+            "xyztau",
+            "xythetatau",
+            "xyetatau",
+            "rhophiztau",
+            "rhophithetatau",
+            "rhophietatau",
+        ):
+            tr1, tr2 = getattr(v1, "to_" + t1)(), getattr(v2, "to_" + t2)()
+            assert numpy.allclose(tr1.deltaRapidityPhi(tr2), expected_result)
diff --git a/tests/compute/lorentz/test_deltaRapidityPhi2.py b/tests/compute/lorentz/test_deltaRapidityPhi2.py
new file mode 100644
index 00000000..151174be
--- /dev/null
+++ b/tests/compute/lorentz/test_deltaRapidityPhi2.py
@@ -0,0 +1,120 @@
+# Copyright (c) 2019-2021, Jonas Eschle, Jim Pivarski, Eduardo Rodrigues, and Henry Schreiner.
+#
+# Distributed under the 3-clause BSD license, see accompanying file LICENSE
+# or https://github.com/scikit-hep/vector for details.
+
+import numpy
+import pytest
+
+import vector._backends.numpy_
+import vector._backends.object_
+
+
+def test_lorentz_object():
+    v1 = vector._backends.object_.MomentumObject4D(
+        vector._backends.object_.AzimuthalObjectXY(1.0, 1.0),
+        vector._backends.object_.LongitudinalObjectZ(1.0),
+        vector._backends.object_.TemporalObjectTau(1.0),
+    )
+    v2 = vector._backends.object_.MomentumObject4D(
+        vector._backends.object_.AzimuthalObjectXY(-1.0, -1.0),
+        vector._backends.object_.LongitudinalObjectZ(-1.0),
+        vector._backends.object_.TemporalObjectTau(1.0),
+    )
+    expected_result = (
+        # phi
+        numpy.pi**2
+        # rapidity
+        + ((0.5 * numpy.log(3 / 1) - 0.5 * numpy.log(1 / 3)) ** 2)
+    )
+    assert v1.deltaRapidityPhi2(v2) == pytest.approx(expected_result)
+
+    for t1 in (
+        "xyzt",
+        "xythetat",
+        "xyetat",
+        "rhophizt",
+        "rhophithetat",
+        "rhophietat",
+        "xyztau",
+        "xythetatau",
+        "xyetatau",
+        "rhophiztau",
+        "rhophithetatau",
+        "rhophietatau",
+    ):
+        for t2 in (
+            "xyzt",
+            "xythetat",
+            "xyetat",
+            "rhophizt",
+            "rhophithetat",
+            "rhophietat",
+            "xyztau",
+            "xythetatau",
+            "xyetatau",
+            "rhophiztau",
+            "rhophithetatau",
+            "rhophietatau",
+        ):
+            tr1, tr2 = getattr(v1, "to_" + t1)(), getattr(v2, "to_" + t2)()
+            assert tr1.deltaRapidityPhi2(tr2) == pytest.approx(expected_result)
+
+
+def test_lorentz_numpy():
+    v1 = vector._backends.numpy_.VectorNumpy4D(
+        [(1.0, 1.0, 1.0, 1.0)],
+        dtype=[
+            ("x", numpy.float64),
+            ("y", numpy.float64),
+            ("z", numpy.float64),
+            ("tau", numpy.float64),
+        ],
+    )
+    v2 = vector._backends.numpy_.VectorNumpy4D(
+        [(-1.0, -1.0, -1.0, 1.0)],
+        dtype=[
+            ("x", numpy.float64),
+            ("y", numpy.float64),
+            ("z", numpy.float64),
+            ("tau", numpy.float64),
+        ],
+    )
+    expected_result = (
+        # phi
+        numpy.pi**2
+        # rapidity
+        + ((0.5 * numpy.log(3 / 1) - 0.5 * numpy.log(1 / 3)) ** 2)
+    )
+    assert v1.deltaRapidityPhi2(v2) == pytest.approx(expected_result)
+
+    for t1 in (
+        "xyzt",
+        "xythetat",
+        "xyetat",
+        "rhophizt",
+        "rhophithetat",
+        "rhophietat",
+        "xyztau",
+        "xythetatau",
+        "xyetatau",
+        "rhophiztau",
+        "rhophithetatau",
+        "rhophietatau",
+    ):
+        for t2 in (
+            "xyzt",
+            "xythetat",
+            "xyetat",
+            "rhophizt",
+            "rhophithetat",
+            "rhophietat",
+            "xyztau",
+            "xythetatau",
+            "xyetatau",
+            "rhophiztau",
+            "rhophithetatau",
+            "rhophietatau",
+        ):
+            tr1, tr2 = getattr(v1, "to_" + t1)(), getattr(v2, "to_" + t2)()
+            assert numpy.allclose(tr1.deltaRapidityPhi2(tr2), expected_result)