forward mode differentiation issue with scipy.special.xlogy
### Description

A few days ago, a [commit](https://github.com/google/jax/commit/dd023e266e6616494cdd590959103ddc646109c4) to jax.scipy.special seems to be causing my previously-working code to bomb out with NaN errors.  I believe the root problem is that you've added a custom jvp to xlogy and xlog1py which need the same safety features as the original functions.

```@custom_derivatives.custom_jvp
@_wraps(osp_special.xlogy, module='scipy.special')
def xlogy(x: ArrayLike, y: ArrayLike) -> Array:
  # Note: xlogy(0, 0) should return 0 according to the function documentation.
  x, y = promote_args_inexact("xlogy", x, y)
  x_ok = x != 0.
  safe_x = jnp.where(x_ok, x, 1.)
  safe_y = jnp.where(x_ok, y, 1.)
  return jnp.where(x_ok, lax.mul(safe_x, lax.log(safe_y)), jnp.zeros_like(x))

def _xlogy_jvp(primals, tangents):
  (x, y) = primals
  (x_dot, y_dot) = tangents
  result = xlogy(x, y)
  return result, (x_dot * lax.log(y) + y_dot * x / y).astype(result.dtype)
xlogy.defjvp(_xlogy_jvp)
```
My code is obtaining NaNs in _xlogy_jvp when I call scipy.special.entr(p) which then calls xlogy(p,p) and one of the elements of p is zero.

### What jax/jaxlib version are you using?

jax v0.4.9

### Which accelerator(s) are you using?

TPU

### Additional system info

Colab

### NVIDIA GPU info

_No response_
diff --git a/jax/_src/scipy/special.py b/jax/_src/scipy/special.py
index fb33d9e92..d8bedc771 100644
--- a/jax/_src/scipy/special.py
+++ b/jax/_src/scipy/special.py
@@ -155,12 +155,24 @@ def _xlog1py_jvp(primals, tangents):
   return result, (x_dot * lax.log1p(y) + y_dot * x / (1 + y)).astype(result.dtype)
 xlog1py.defjvp(_xlog1py_jvp)
 
+@custom_derivatives.custom_jvp
+def _xlogx(x):
+  """Compute x log(x) with well-defined derivatives."""
+  return xlogy(x, x)
+
+def _xlogx_jvp(primals, tangents):
+  x, = primals
+  x_dot, = tangents
+  return  _xlogx(x), x_dot * (lax.log(x) + 1)
+_xlogx.defjvp(_xlogx_jvp)
+
+
 @_wraps(osp_special.entr, module='scipy.special')
 def entr(x: ArrayLike) -> Array:
   x, = promote_args_inexact("entr", x)
   return lax.select(lax.lt(x, _lax_const(x, 0)),
                     lax.full_like(x, -np.inf),
-                    lax.neg(xlogy(x, x)))
+                    lax.neg(_xlogx(x)))
 
 
 @_wraps(osp_special.multigammaln, update_doc=False)
diff --git a/tests/lax_scipy_test.py b/tests/lax_scipy_test.py
index 3ecf24e67..3e8c84398 100644
--- a/tests/lax_scipy_test.py
+++ b/tests/lax_scipy_test.py
@@ -243,6 +243,11 @@ class LaxBackedScipyTests(jtu.JaxTestCase):
 
     jtu.check_grads(lsp_special.xlog1py, (x0, y0), order=2)
 
+  def testGradOfEntrAtZero(self):
+    # https://github.com/google/jax/issues/15709
+    self.assertEqual(jax.jacfwd(lsp_special.entr)(0.0), jnp.inf)
+    self.assertEqual(jax.jacrev(lsp_special.entr)(0.0), jnp.inf)
+
   @jtu.sample_product(
     [dict(order=order, z=z, n_iter=n_iter)
      for order, z, n_iter in zip(
['tests/lax_scipy_test.py', 'jax/_src/scipy/special.py']