quantile fails on an empty dimension where other reductions return NaN

What happened?

quantile on a dimension of length zero fails three different ways, while the other reductions return nan.

empty.mean() / .median() / .std()   -> nan
empty.quantile(0.5)                 -> AxisError: source: axis 0 is out of bounds for array of dimension 0
empty.quantile(0.5, skipna=False)   -> IndexError: index -1 is out of bounds for axis 0 with size 0
DataArray(zeros((0,3))).quantile(0.5, dim="y") -> ValueError: applied function returned data with an unexpected number of dimensions

That last one reduces over y, which has length 3. Only x is empty.

median and quantile(0.5) are the same operation and disagree, and pandas returns nan:

pd.Series([], dtype=float).quantile(0.5)   # nan

Where it comes from

np.nanquantile drops the leading q axis when the reduced axis is empty:

np.nanquantile(np.zeros((3,)),  [0.5], axis=-1).shape   # (1,)
np.nanquantile(np.zeros((0,)),  [0.5], axis=-1).shape   # ()      <- q axis gone
np.nanquantile(np.zeros((2,3)), [0.5], axis=-1).shape   # (1, 2)
np.nanquantile(np.zeros((0,3)), [0.5], axis=-1).shape   # (0,)    <- q axis gone

Variable.quantile._wrapper does moveaxis(_quantile_func(npa, **kwargs), 0, -1), which assumes that axis is there, so apply_ufunc then rejects the result. The skipna=False path goes to np.quantile, which raises on empty input by itself.

What did you expect to happen?

nan, matching mean, median, std and pandas, with the same result for both skipna values.

Minimal Complete Verifiable Example

import numpy as np
import xarray as xr

empty = xr.DataArray(np.array([], dtype=float), dims="x", coords={"x": []})
print(empty.median().item())     # nan
print(empty.quantile(0.5))       # AxisError

# reached without building an empty array by hand
d = xr.DataArray([1.0, 2.0], dims="x", coords={"x": [0, 1]})
d.where(d > 99, drop=True).quantile(0.5)

Environment

xarray 2026.7.1.dev43+ga48f3152c (main at a48f315), numpy 2.5.2, pandas 3.0.5, python 3.12.11

Start at Variable.quantile._wrapper and reproduce the empty-dimension examples with both skipna values, including reduction over y. Check the quantile and reduction tests around these cases, then add coverage showing that empty quantiles return nan with consistent dimensions and match median, mean, std, and pandas.