task: metadata: install: pip install -e .
patch: diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 89a35c3b0..d11095c95 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -167,6 +167,10 @@
 * The `PauliRot` operation now gracefully handles single-qubit Paulis, and all-identity Paulis
   [(#860)](https://github.com/PennyLaneAI/pennylane/pull/860).
 
+* Fixes a bug whereby binary Python operators were not properly propagating the `requires_grad`
+  attribute to the output tensor.
+  [(#889)](https://github.com/PennyLaneAI/pennylane/pull/889)
+
 <h3>Contributors</h3>
 
 This release contains contributions from (in alphabetical order):
diff --git a/pennylane/numpy/tensor.py b/pennylane/numpy/tensor.py
index 57e4bef3e..c8ab8803c 100644
--- a/pennylane/numpy/tensor.py
+++ b/pennylane/numpy/tensor.py
@@ -113,6 +113,55 @@ class tensor(_np.ndarray):
         out_arr = tensor(obj, requires_grad=self.requires_grad)
         return super().__array_wrap__(out_arr)
 
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        # pylint: disable=no-member,attribute-defined-outside-init
+
+        # unwrap any outputs the ufunc might have
+        outputs = [i.view(onp.ndarray) for i in kwargs.get("out", ())]
+
+        if outputs:
+            # Insert the unwrapped outputs into the keyword
+            # args dictionary, to be passed to ndarray.__array_ufunc__
+            outputs = tuple(outputs)
+            kwargs["out"] = outputs
+        else:
+            # If the ufunc has no ouputs, we simply
+            # create a tuple containing None for all potential outputs.
+            outputs = (None,) * ufunc.nout
+
+        # unwrap the input arguments to the ufunc
+        args = [i.unwrap() if hasattr(i, "unwrap") else i for i in inputs]
+
+        # call the ndarray.__array_ufunc__ method to compute the result
+        # of the vectorized ufunc
+        res = super().__array_ufunc__(ufunc, method, *args, **kwargs)
+
+        if ufunc.nout == 1:
+            res = (res,)
+
+        # construct a list of ufunc outputs to return
+        ufunc_output = [
+            (onp.asarray(result) if output is None else output)
+            for result, output in zip(res, outputs)
+        ]
+
+        # if any of the inputs were trainable, the output is also trainable
+        requires_grad = any(
+            isinstance(x, onp.ndarray) and getattr(x, "requires_grad", True) for x in inputs
+        )
+
+        # Iterate through the ufunc outputs and convert each to a PennyLane tensor.
+        # We also correctly set the requires_grad attribute.
+        for i in range(len(ufunc_output)):  # pylint: disable=consider-using-enumerate
+            ufunc_output[i] = tensor(ufunc_output[i], requires_grad=requires_grad)
+
+        if len(ufunc_output) == 1:
+            # the ufunc has a single output so return a single tensor
+            return ufunc_output[0]
+
+        # otherwise we must return a tuple of tensors
+        return tuple(ufunc_output)
+
     def __getitem__(self, *args, **kwargs):
         item = super().__getitem__(*args, **kwargs)
 

problem_statement: Arithmetic operations between pennylane.numpy arrays have inconsistent `requires_grad` output 
**Issue**

The PennyLane wrapper if the numpy library, which is used to optimise circuits with the standard autograd interface, has an inconsistent behaviour when using arithmetic functions such as +, -, * (here using the sum to demonstrate). 

Summing two `pennylane.numpy.array` objects, one having gradients enabled and the other not, returns a differentiable array or not, depending on the order.

``` python
import pennylane as qml
from pennylane import numpy as np

a = np.array(1., requires_grad=True)
b = np.array(0., requires_grad=False)

c = a+b
print(c.requires_grad) # True

d = b+a
print(d.requires_grad) # False
```
**Expected behaviour**
The `requires_grad = True` flag should always be overwriting, and the output of the two print statements should be true in both cases. This is important for example in loss function where target outputs (which do not require gradients) and predictions (which depend on the parameters and require gradients) are compared, i.e. via `np.sum(np.abs((target - prediction))**2)`.

**System information**
Reproducible with PennyLane v0.12 stable release, and with the current master, in tape mode or not in tape mode.

repo: PennyLaneAI/pennylane
test_cmd: pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
test_patch: diff --git a/tests/test_numpy_wrapper.py b/tests/test_numpy_wrapper.py
index 9d4f824f6..e77b57c03 100644
--- a/tests/test_numpy_wrapper.py
+++ b/tests/test_numpy_wrapper.py
@@ -350,6 +350,58 @@ class TestNumpyIntegration:
         assert isinstance(res, np.tensor)
         assert not res.requires_grad
 
+    def test_binary_operator_mixed_trainable_left(self):
+        """Test that binary operators on one trainable and one non-trainable
+        arrays result in trainable output."""
+        x = np.array([[1, 2], [3, 4]], requires_grad=True)
+        y = np.array([[5, 6], [7, 8]], requires_grad=False)
+
+        res = x + y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x - y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x / y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x * y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x @ y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+    def test_binary_operator_mixed_trainable_right(self):
+        """Test that binary operators on one non-trainable and one trainable
+        arrays result in trainable output."""
+        x = np.array([[1, 2], [3, 4]], requires_grad=False)
+        y = np.array([[5, 6], [7, 8]], requires_grad=True)
+
+        res = x + y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x - y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x / y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x * y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
+        res = x @ y
+        assert isinstance(res, np.tensor)
+        assert res.requires_grad
+
 
 class TestAutogradIntegration:
     """Test autograd works with the new tensor subclass"""