Thanks @eendebakpt the change LGTM (unrelated return changes are also fine). I doubt it will make issues, but maybe add a very short release note anyway?

Arrg, I noticed one thing that maybe we should defer. But our scalars do the same dance and still return the Python boolean.
And while it is easy to change here, I am not actually sure that it can be changed for the scalars because they may have to adhere to the NotImplemented protocol.

Without any comments, I may just go ahead soon, but I want to let it sink in briefly for myself.

Which scalar operations return the Python bool? I checked with

import numpy as np
print(np, np.__version__)

for pair in [( np.array(1), np.array('s')), ( np.array(1), np.int16(2)), ( np.array(1), 2),  (np.float16(1), np.int16(1)),  (np.float16(1), 3)]:
    r= pair[0] == pair[1]
    print(f'{pair}: {r=},\t{type(r)=}')
    r= pair[1] == pair[0]
    print(f'{pair[::-1]}: {r=},\t{type(r)=}')

and they return a numpy bool.

For the NotImplemented protocol (I checked https://docs.python.org/3/reference/datamodel.html#object.__eq__) the method should either return NotImplemented , or an object that can be cast to True/False using bool. Returning a numpy.bool would be fine then.

Ah, I guess only the ones for which the comparison is not defined np.float64(3) != np.str_("a"). I.e. this path, but for scalars.

Thanks. So we have:

np.float64(3) != np.int16(4) # numpy.bool
np.float64(3) != np.str_("a") # python bool

This is not causing any issues (as far as i know), but it is a bit inconsistent.

Note to self: python bool might be returned here

OK, thanks @eendebakpt let's give this a try. FWIW, we should of course return a 0-D array in some future (which this does not make any harder, I suspect it just means deleting the code branch completely).
But for now the scalar is right. I dislike the difference to the normal scalars, but c'est la vie.