Return the Jacobian from `ews_to_stellar_parameters`

By the time `ews_to_stellar_parameters` returns, it has the full Jacobian of the
four residual equations at the solution. It uses that once to make four error
bars and then drops it. I'd like it back.

The four classical conditions are strongly correlated, so four independent
error bars are a poor summary of what was actually fit. Everything needed for
the joint,

Σ_θ = J⁻¹ Σ_r J⁻ᵀ

is in scope at `fit_via_EWs.jl:424` and only `abs.(J \ σ_r)` comes out:

```julia
J = DiffResults.jacobian(J_result)[.!fix_params, .!fix_params]
param_uncertainties = zeros(4)
param_uncertainties[.!fix_params] .= abs.(J \ residual_uncertainties[.!fix_params])

params, param_uncertainties
```

`abs.()` on a linear solve throws away the sign structure that carries the
correlation, so it cannot be reconstructed from the four numbers afterwards.

We rebuild `J` ourselves by re-running ForwardDiff over
`_stellar_param_equation_residuals`, and the correlations are large:
ρ(Teff, [m/H]) = 0.87 on a solar HARPS EW table (169 Fe lines), 0.83 ± 0.09
across 56 FGK stars from the Gaia FGK Benchmark Stars library. Isochrone and SED
fitting routinely take the four parameters as independent, and at ρ ≈ 0.85 that
is a materially wrong likelihood.

What would help, in order of preference:

1. Return `J` and `residual_uncertainties` alongside the parameters — optional
named tuple field, or behind a keyword so the two-value return stays default.
2. Return `Σ_θ`. `fit_spectrum` already has a `covariance` field, so this would
make the two entry points consistent.
3. Document that the returned uncertainties are marginals of a correlated
posterior and should not be used as independent errors.

Recomputing `J` externally costs a full ForwardDiff pass over the exact
`ews_to_abundances` — 31% of our per-star wall time, duplicating a matrix Korg
has already built at that point.

I can PR (1). Smallest change, non-breaking.