Array API: subtract_overscan model fitting goes through numpy-only astropy.modeling

### Problem

`subtract_overscan` with `model=...` fits via `astropy.modeling`:

https://github.com/astropy/ccdproc/blob/9d25eeefda8a38fc442c1a18e79ca2fa8ca40559/ccdproc/core.py#L616

astropy's fitters operate on numpy internally (`LinearLSQFitter` ends up in `np.linalg.lstsq`), so CuPy input hits the implicit-conversion error. This matches the `test_subtract_overscan_model` failures noted in #910.

### Options

- Treat model fitting as a CPU-only operation under the policy in #935 (explicit host transfer + convert the evaluated model back to `xp`).
- For the common low-order polynomial case, an xp-native least-squares fit (`xp.linalg.lstsq` is in the array API linalg extension) could keep everything on-device.
- Longer term: array-API support in `astropy.modeling` (tracked in #940).

---
Found during a review of the array API implementation from #885; follow-up to #909 / #910.

Start in ccdproc/core.py around line 616 and review the test_subtract_overscan_model failures noted in #910. Determine whether the fix should use the CPU-transfer policy in #935 or an xp-native fit, then verify that CuPy input avoids implicit conversion and the relevant model-fitting tests pass.

Written by the indexing model from the issue text.