Fitting · Payne Zero

Fit the full spectrum because abundance responses are coupled

A stellar abundance is not measured only by the isolated lines carrying that element's name. Blends allow several species to affect one pixel. Magnesium and other electron donors can change the H⁻ continuum. Carbon, nitrogen, and oxygen share molecular equilibria, and a changed mixture can alter line blanketing and the atmospheric structure. A full-spectrum fit lets all retained pixels constrain these connected responses together.

Direct synthesis also makes a discrepancy easier to locate. A mismatch can enter through the stellar state, the atmosphere, line data, continuum treatment, velocity, broadening, the instrument line-spread function, or the observations themselves. Payne Zero keeps those operations separate, so a fit can test the relevant physical or instrumental cause rather than placing an empirical correction directly on the flux.

Begin with a model and observation on clearly defined grids

The synthesis grid is chosen to resolve the intrinsic opacity and line profiles. The observed grid is set by detector pixels and the reduction. A fit should not compare those arrays until the model has passed through the same velocity, broadening, line-spread, and sampling assumptions used for the observation.

NormalizedSpectrum stores the observed wavelength, normalized flux, inverse variance, and Boolean mask. Invalid or excluded pixels remain explicit in the mask rather than being removed silently, so the model, data, and weights retain a stable one-to-one indexing.

Project total and continuum flux through the instrument

ObservedSpectrumOperator starts from the intrinsic logarithmic wavelength grid and produces flux on any strictly increasing observed grid inside that interval. It can apply a residual Doppler shift, Gaussian broadening, a constant resolving power or supplied shift-invariant line-spread kernel, and final resampling.

from fitter import ObservedSpectrumOperator

instrument = ObservedSpectrumOperator(
    native_wavelength_nm,
    observed_wavelength_nm,
    resolving_power=22_500,
    device=device,
    dtype=torch_dtype,
)
instrument.set_parameters(
    radial_velocity_km_s=residual_velocity,
    broadening_sigma_km_s=broadening_sigma,
)

This instrument model is passed to synthesis, which transforms total and continuum flux consistently before returning normalized flux. Wavelength-dependent or detector-specific line-spread functions use the same interface through an output wavelength array and a convolve_fluxes method. The APOGEE adapter uses that form.

Project a generated spectrum to a simple observed grid

This example starts from the converged spectrum generated in the tutorial, builds 180 output pixels, and applies a constant resolving power together with residual Gaussian broadening:

observed_wavelength_nm = np.linspace(
    physical_spectrum.wavelength_nm[3],
    physical_spectrum.wavelength_nm[-4],
    180,
)

instrument = ObservedSpectrumOperator(
    physical_spectrum.wavelength_nm,
    observed_wavelength_nm,
    resolving_power=25_000,
    device=device,
    dtype=torch_dtype,
)
instrument.set_parameters(
    radial_velocity_km_s=0.0,
    broadening_sigma_km_s=2.0,
)

observed_model = synthesize(
    atmosphere_path,
    wavelength_start_nm=500.0,
    wavelength_end_nm=502.0,
    resolution=300_000,
    spectral_operator=instrument,
    device=device,
    dtype=torch_dtype,
)

A converged solar spectrum from 500 to 502 nanometers on its intrinsic grid and after line-spread convolution and sampling to 180 observed pixels — A converged solar spectrum over 500–502 nm on its intrinsic R = 300,000 grid, followed by line-spread convolution and sampling onto 180 pixels at R = 25,000. Broadening reduces the contrast of features that are not resolved by the delivered instrument model.

Use a sampled line-spread kernel

If one shift-invariant instrument kernel is already known, provide its odd-length nonnegative weights instead of a resolving power. The samples are spaced in native log-wavelength pixels and are normalized internally.

lsf_kernel = np.array([0.02, 0.12, 0.32, 0.44, 0.32, 0.12, 0.02])

instrument = ObservedSpectrumOperator(
    native_wavelength_nm,
    observed_wavelength_nm,
    lsf_kernel=lsf_kernel,
    device=device,
    dtype=torch_dtype,
)

Rotational broadening when it belongs in the model

Rotation can be applied on the intrinsic grid before the instrument line-spread function. The Gray-profile implementation remains in the Torch graph for flux derivatives:

from fitter import RotationalBroadening

rotation = RotationalBroadening(
    native_wavelength_nm,
    maximum_vsini_km_s=100.0,
    limb_darkening=0.6,
    device=native_flux.device,
    dtype=native_flux.dtype,
)

rotated_flux = rotation(native_flux, vsini_km_s=12.0)

The scalar vsini_km_s value is intended for bounded or finite-difference fitting. The current interface does not calculate its derivative automatically.

Weights, masks, and continuum terms define the comparison

The residual vector is evaluated only where the mask is true and is weighted by inverse variance. The fitter can also profile a smooth multiplicative continuum. A user-supplied basis with shape(N, K) is multiplied by the physical model, and its coefficients are solved by weighted linear least squares for every trial set of stellar parameters.

Profiling keeps smooth detector or normalization structure out of the fitted stellar parameters. It also makes the modeling assumption visible: the chosen basis determines which broad differences between model and data the continuum is allowed to absorb.

Fit the stellar parameters within declared physical bounds

FitConfiguration fixes parameter order, starting values, lower and upper bounds, the step used to measure each spectral response, the largest permitted update, and the stopping conditions. A user-supplied model function receives one trial set of parameters and returns normalized model flux on the observed wavelength grid.

import numpy as np
from fitter import (
    FitConfiguration,
    NormalizedSpectrum,
    fit_normalized_spectrum,
)

observation = NormalizedSpectrum(
    wavelength=observed_wavelength_nm,
    flux=measured_normalized_flux,
    inverse_variance=measured_inverse_variance,
    mask=np.isfinite(measured_normalized_flux),
)

configuration = FitConfiguration(
    names=("M_H",),
    initial=np.array([-0.05]),
    lower=np.array([-1.0]),
    upper=np.array([0.3]),
    derivative_steps=np.array([0.05]),
    trust_half_width=np.array([0.30]),
    maximum_iterations=4,
)

def model(parameters):
    labels = {
        **solar_labels,
        "metallicity": float(parameters[0]),
    }
    return synthesize_from_labels(
        **labels,
        wavelength_start_nm=500.0,
        wavelength_end_nm=502.0,
        r_grid=300_000,
        spectral_operator=instrument,
        device=device,
        dtype=torch_dtype,
    ).normalized_flux

result = fit_normalized_spectrum(
    observation,
    configuration,
    model,
)
result.save("runs/metallicity_fit")

At each step, the fitter measures how the spectrum responds to the active parameters, either from supplied derivatives or by evaluating small one-sided changes. It solves a weighted local least-squares problem, limits the proposed update to the declared bounds and maximum step, and accepts a trial only when it improves the weighted fit to the data.

Adding an abundance generally adds another response spectrum to a finite-difference iteration. The cost therefore grows with the active parameters and their synthesis evaluations, not with a precomputed dense label grid.

In the controlled example below, [M/H] is the logarithmic metal-to-hydrogen ratio relative to the Sun. The symbol σ denotes one standard deviation of the injected pixel noise.

A controlled noisy mock spectrum and the fitted Payne Zero model, with the data-minus-model residual compared with the injected one-sigma noise — This controlled calculation uses the same one-label fitting interface demonstrated by the public notebook. It synthesizes a spectrum at [M/H] = −0.300, adds Gaussian noise with 1σ = 0.0025, and starts at [M/H] = −0.050. Four bounded iterations recover [M/H] = −0.301. This is an internal-consistency example, not a measurement of accuracy on observed stars.

Fit with a fixed atmosphere, then check the physical atmosphere

Converging a new atmosphere at every trial would spend most of the fit updating structures that will be discarded. The intended strategy separates exploration from acceptance. Repeated fixed-atmosphere synthesis evaluations first locate a candidate; the atmosphere solver is then called at that candidate, and the converged atmosphere is synthesized through the same instrument model and continuum treatment.

A stellar-parameter fit followed by a converged-atmosphere evaluation and acceptance check — The search proposes a candidate. A separately converged atmosphere then tests whether convergence materially changes the flux or the fit before a result is retained.

refine_with_physical_atmosphere records two separate checks: whether the fixed-atmosphere and converged spectra agree within the declared tolerance after continuum profiling, and whether the converged model preserves the weighted fit statistic. If a material discrepancy remains, the fixed-atmosphere model proposes a small parameter correction within the declared bounds. Another atmosphere is solved only when that correction is expected to improve the fit, and it is retained only if the converged model actually improves.

from tempfile import TemporaryDirectory

from fitter import (
    PhysicalAtmosphereConfiguration,
    refine_with_physical_atmosphere,
)
from payne_zero_atmosphere import solve_structured_atmosphere
from payne_zero_synthesis import synthesize

def converged_model(parameters):
    labels = {
        **solar_labels,
        "metallicity": float(parameters[0]),
    }
    with TemporaryDirectory() as directory:
        atmosphere_path = solve_structured_atmosphere(
            **labels,
            out_dir=directory,
        )
        return synthesize(
            atmosphere_path,
            wavelength_start_nm=500.0,
            wavelength_end_nm=502.0,
            resolution=300_000,
            spectral_operator=instrument,
            device=device,
            dtype=torch_dtype,
        ).normalized_flux

physical_result = refine_with_physical_atmosphere(
    observation,
    configuration,
    result,
    model,
    converged_model,
    PhysicalAtmosphereConfiguration(
        maximum_discrepancy_rms=2.0e-3,
        maximum_objective_degradation=0.1,
        minimum_predicted_objective_improvement=1.0e-3,
        maximum_physical_evaluations=4,
    ),
)
physical_result.save("runs/metallicity_fit_physical")

The final scientific product is therefore the model evaluated from a converged atmosphere. The earlier search result remains useful as a candidate and as a record of how the fitting routine reached that state, but it does not show that its atmosphere has converged.

APOGEE pixels require a wavelength-dependent instrument model

The packaged APOGEE DR14 example applies the same sequence to 7,514 retained pixels. In place of one constant resolving power, it uses a wavelength-dependent line-spread function, detector sampling, and the survey mask. It also fits a small residual velocity and an effective Gaussian broadening width, while solving a smooth continuum separately on each detector.

The public array interface fits effective temperature, surface gravity, bulk metallicity, alpha enhancement, and microturbulence; an option adds independent C, N, and O abundances. Use the generic fitter for a different element-by-element mixture. The packaged line-spread function is a representative DR14 all-slit mean, not an official visit-, fiber-, or star-specific ASPCAP kernel.

The default notebook loads a retained fit so the complete data preparation and model comparison can be inspected without rerunning the longer optimization. Setting PAYNE_ZERO_RUN_APOGEE=1 executes the fit on the packaged arrays.

The retained Payne Zero model and an APOGEE DR14 giant spectrum across all three detector wavelength ranges, with residuals in units of the pixel uncertainty — The public retained example for APOGEE DR14 object 2M08002084+4044415 across all three detectors. Black is the observed spectrum and orange is the seven-parameter model: effective temperature, surface gravity, bulk composition, alpha enhancement, microturbulence, residual velocity, and Gaussian broadening. The lower panel shows data minus model in units of each pixel's quoted 1σ uncertainty; the gray band marks ±1σ. The remaining structure is visible rather than removed by selecting a favorable line window.

The installed command accepts the same retained-pixel arrays:

payne-zero-fit-apogee \
  examples/data/apogee_dr14_example.npz \
  runs/apogee_example \
  --object-id 2M08002084+4044415 \
  --reference-labels 4858.537 2.426797 -0.3255714 0.2527028 1.278733 \
  --reference-vmacro 3.617606 \
  --synthesis-r-grid 300000 \
  --fresh-jacobian-rounds 0 \
  --force

The command uses the default intrinsic grid density R = 300,000 before the APOGEE line-spread function and detector sampling. This R describes numerical wavelength sampling, not the delivered resolving power of the spectrograph.

The paper tests an element-by-element survey fit

The manuscript applies the same instrument and fitting machinery to a separate 1,600-giant experiment. [Fe/H] is the logarithmic iron-to-hydrogen ratio relative to the Sun. [X/Fe] states how one element X differs from the solar element-to-iron ratio. The calculation fits [Fe/H] and eleven independent [X/Fe] ratios—C, N, O, Mg, Al, Si, S, Ca, Ti, Mn, and Ni—while retaining the published DR14 effective temperature, surface gravity, and microturbulence in its main configuration. A sensitivity calculation also frees those three stellar quantities. Each accepted result is recalculated with a physically converged atmosphere.

Calibrate line dataMove from stellar-parameter fitting to catalog-specific atomic-parameter corrections.