Automated Lyman-Alpha Emission-Line Fitting for
IUE Short-Wavelength Prime Spectra of 3C 273

Introduction

The quasar 3C 273 (RA 12h 29m 06.7s, Dec +02° 03′ 09″; redshift \(z = 0.158\)) is among the brightest and most studied quasi-stellar objects in the sky. Its UV continuum and broad emission lines have been monitored extensively with the International Ultraviolet Explorer (IUE) satellite, which operated from 1978 to 1996 and accumulated a uniquely dense, long-baseline UV spectral archive for this source.^[1,2]

The Short-Wavelength Prime (SWP) camera covers 1150–1975 Å in low-dispersion mode at a resolving power of \(R \approx 300\) (\(\Delta\lambda \approx 4\)–6 Å), capturing Lyman-Alpha at its redshifted observed wavelength \(\lambda_\text{obs} \approx 1408\) Å. The campaign files have been divided by \(1+z\), so this pipeline works in the rest frame and the line appears near 1213–1218 Å.

Measuring Ly\(\alpha\) flux accurately is critical for: (a) reverberation mapping of the broad-line region (BLR), (b) variability analysis of the accretion disc, and (c) quantifying the relationship between continuum and emission-line fluxes on year-to-decade timescales.^[4,5] Automated, reproducible fitting across all 255 SWP spectra requires a pipeline robust to varying S/N, continuum slope, and line asymmetry.

End-to-End Pipeline Overview

The pipeline is implemented in Python and exposed through a Flask web application. Each spectrum file (two-column ASCII: wavelength in Å, flux in erg s⁻¹ cm⁻² Å⁻¹) is processed independently. The eight principal stages are:

Data Ingestion & Pre-processing

Each IUE SWP low-dispersion spectrum is a two-column text file: wavelength \(\lambda\) (Å) and flux density \(F_\lambda\) (erg s\(^{-1}\) cm\(^{-2}\) Å\(^{-1}\)). The 255 campaign files are rest-frame spectra: they were divided by \(1+z\) with \(z = 0.158\) for 3C 273. Two facts confirm this: the pixel spacing is \(1.448\) Å \(= 1.676/(1+z)\) Å, where 1.676 Å is the SWP MXLO sampling in the observed frame, and geocoronal Ly\(\alpha\) airglow (1215.67 Å observed) appears at \(1215.67/1.158 = 1049.8\) Å.

Usable data run from 993.6 Å to ≈1832 Å (pixels with zero flux at the blue edge are ignored). No rebinning is applied. Observation dates (JD) are taken from the SWP log jd.xlsx; five spectra (251–255) have no recorded date.

Continuum Modeling: Power-Law Fit

The UV continuum of a quasar/AGN is well described by a power law \(F_c(\lambda) = A\,\lambda^\alpha\), where \(A\) is the continuum amplitude and \(\alpha\) is the spectral index (typically \(\alpha \approx -1.5\) to \(-2.5\) for UV continua).^[6] Taking the natural logarithm linearises the problem:

Ordinary least squares is applied to \((\ln\lambda_i, \ln F_i)\) for all pixels in five line-free rest-frame windows: 1100–1150, 1325–1370, 1425–1475, 1590–1620 and 1660–1700 Å. These are the windows of the campaign tables; with them the spectral indices reproduce the tabulated values for all 255 spectra (maximum difference \(5\times10^{-5}\)). The 1325–1370 Å window contains the weak C II λ1335 blend, which is undetected in most epochs.^[7]

Least-Squares Solution

Let \(x_i = \ln\lambda_i\) and \(y_i = \ln F_i\) for \(n\) continuum pixels. The OLS estimators are:

The best-fit continuum \(\hat{F}_c(\lambda) = \hat{A}\,\lambda^{\hat{\alpha}}\) is then evaluated on the full wavelength grid, and the continuum-subtracted (emission-only) spectrum is:

Uncertainties

With residual variance \(s^2 = \sum (y_i - \hat\alpha x_i - \ln\hat A)^2/(n-2)\) and \(S_{xx} = \sum (x_i - \bar x)^2\), the standard errors of the OLS estimators are^[10]

The campaign tables list \(\sigma_A\) values that are \(\sqrt{n}\) (≈12) times larger, consistent with the \(1/n\) term being omitted; the pipeline reports Eq. 3a. For Spectrum 70: \(\hat\alpha = -0.933 \pm 0.055\), \(\hat A = (1.57 \pm 0.62)\times10^{-10}\).

Stage A: Discrete Candidate Grid Search

5.1 Design Philosophy

Rather than feeding a continuous solver a single starting point and trusting it to find the global minimum, the candidate engine takes a decision-making approach: enumerate a physically constrained set of Gaussian models and rank them by a multi-criteria score that mirrors how a human expert would judge a fit. This avoids local-minimum traps common in UV emission-line fitting where low S/N or asymmetric profiles create false optima.^[8]

5.2 Gaussian Emission-Line Model

Each candidate emission-line model is a single Gaussian:

where \(A\) is the peak amplitude (erg s\(^{-1}\) cm\(^{-2}\) Å\(^{-1}\)), \(\mu\) is the line center (Å), and \(\sigma\) is the Gaussian width (Å).

5.3 Parameter Grids

The search space is defined over three parameters:

Why these ranges. A Gaussian has FWHM \(= 2\sqrt{2\ln 2}\,\sigma \approx 2.3548\,\sigma\), so the Ly\(\alpha\) FWHM range of 10–24 Å corresponds to \(\sigma = 4.25\)–10.19 Å (2500–5900 km s\(^{-1}\) at 1216 Å). This brackets the approved fits (FWHM 10.6–20.0 Å, 2600–4950 km s\(^{-1}\)) and published broad-line widths of 3C 273 (Paltani & Türler 2005). No line can appear narrower than the instrument: the IUE SWP low-dispersion resolution is ≈6 Å FWHM observed, i.e. 6/(1+z) = 5.2 Å in these de-redshifted spectra (z = 0.158), giving the floor \(\sigma \ge 2.2\) Å used for the weaker lines. The grid only seeds the fit: \(\sigma\) may take any value between its bounds in the least-squares refinement.

5.4 Noise Estimation (Robust MAD)

A robust noise estimate is required both for scoring (\(\chi^2\)) and for SNR computation. The pipeline uses the Median Absolute Deviation (MAD), scaled to match the standard deviation of a normal distribution:^[9]

where \(r_i = F_\text{sub}(\lambda_i) / \bar{F}_c\) are the normalised continuum-subtracted residuals over a ±20 Å window around the line. The factor 1.4826 converts the MAD to a Gaussian-equivalent \(\sigma\). The absolute noise level is then \(\sigma_n = \hat{\sigma}_\text{noise}\,|\bar{F}_c|\).

Composite Scoring Function (seed selection)

Each grid candidate \((\mu, \sigma, W)\), with its amplitude from linear least squares, receives a composite score \(S \in [0, 1]\). The score only chooses the starting point of the least-squares refinement (§7); the reported parameters and errors come from that fit. Weights as implemented in src/candidate_engine.py:

There is no penalty for a low \(\chi^2_\text{red}\): with the coarse IUE sampling, values well below 1 occur for legitimate fits (approved fits reach 0.05).

Stage B: Weighted Least-Squares Refinement

The grid candidate seeds a continuous fit of the observed flux over a per-line window (e.g. 1140–1290 Å for Ly\(\alpha\) + N V, 998–1085 Å for Ly\(\beta\) + O VI, 1495–1595 Å for C IV; the Ly\(\alpha\) and C IV windows are those of the reference hand fits). The model contains the line, any overlapping companion line and a local linear continuum, all fitted simultaneously:

\(\chi^2\) is minimised with the trust-region reflective algorithm (scipy.optimize.least_squares)^[8] under bounds \(A_k \ge 0\), \(|\mu_k - \lambda_{0,k}| \le 8\) Å and \(\sigma_\text{min} \le \sigma_k \le \sigma_\text{max}\) (§5). Fluxes are fitted in units of \(10^{-13}\) so all parameters have comparable scale. Masked ranges (geocoronal airglow at 1043–1056 Å) are excluded. Fitting the continuum together with the line avoids the bias of subtracting a baseline that itself contains line flux.

7.1 Pixel noise

\(\sigma_\text{pix}\) is the DER_SNR estimate^[16] \(\sigma = 1.482602/\sqrt{6}\;\text{median}\,|2F_i - F_{i-2} - F_{i+2}|\), evaluated on the fit residuals \(F_i - M(\lambda_i)\) so that the curvature of the line itself does not inflate it. Because the weights are uniform, the best-fit parameters do not depend on \(\sigma_\text{pix}\); it sets \(\chi^2\) and the uncertainties.

7.2 Uncertainties

The covariance matrix is \(C = (J^\mathsf{T}J)^{-1}\,\max(1, \chi^2_\text{red})\), where \(J\) is the Jacobian of the normalised residuals; scaling by \(\chi^2_\text{red}\) when it exceeds 1 inflates errors for imperfect models.^[10] Parameters held fixed (§7.3) are left out of the fit and carry no variance; \(\nu\) counts only the free parameters.

7.3 Narrow + broad decomposition (Ly\(\alpha\), C IV)

The Ly\(\alpha\) and C IV profiles of 3C 273 have a sharp core on a wide base, which a single Gaussian cannot follow. Each of these lines is therefore modelled as two Gaussians of the same line, a narrow component (N) and a broad component (B), fitted together with the local continuum (and N V for Ly\(\alpha\)) in Eq. 8:

When all six line parameters are free, the two components trade flux and the broad flux is poorly constrained. In a free fit of 51 spectra the broad component reached \(F/\sigma_F \ge 3\) in only 11 (Ly\(\alpha\)) and 21 (C IV). The broad centre and width are therefore held fixed at the median of the reference hand fits (references/CT-3C273-LUCK-EMMY.xls, sheet Check, 15 spectra). Per spectrum, only \(A_B\), \(A_N\), \(\mu_N\) and \(\sigma_N\) are fitted. \(\mu_N\) may move ±3 Å around the observed peak of that spectrum.

The two σ ranges do not overlap, so the components cannot swap roles. The values are set in config/Lya.yaml and config/CIV.yaml (profile: narrow_broad). derive_broad_shape.py gives the alternative measured from the data alone: the median of the free fits with broad \(F/\sigma_F \ge 5\), no parameter at a bound and \(\Delta\text{BIC} > 10\), which gives \(\sigma_B \approx 10.4\) Å (Ly\(\alpha\)) and 12.9 Å (C IV). Which of the two to use is an open question (questions.txt).

Is the second Gaussian justified? Each spectrum is also fitted with one Gaussian over the same pixels, and \(\Delta\text{BIC} = \text{BIC}_1 - \text{BIC}_2\) is reported, with \(\text{BIC} = n\ln(\text{RSS}/n) + k\ln n\). \(\Delta\text{BIC} > 10\) strongly favours narrow + broad (Kass & Raftery 1995).

Naming. Every line carries the short label of the reference sheets (label in its YAML): Ly-B, Lya-N, Lya-B, NV, OI, CII, SIV, CIV-N, CIV-B, HeII, OIII]. Each component is reported, plotted and tabulated as its own emission line, with its own light curve and \(F_\text{var}\).

Physical Statistics Computation

8.1 Integrated flux and its error

For a sum of Gaussians (narrow + broad, §7.3) the same propagation is written as \(F = \sum_k \sqrt{2\pi}A_k\sigma_k\), \(\sigma_F^2 = \mathbf{g}^\mathsf{T} C\, \mathbf{g}\) with \(\mathbf{g} = \partial F/\partial \mathbf{p}\), which includes the covariance between the components. Each component's own flux uses only its own terms.

This is the full Gaussian flux, independent of any window. For comparison with the campaign tables, flux_window also reports the Gaussian integrated only over \([\mu - W/2, \mu + W/2]\) (the convention of the reference tool), which captures \(\operatorname{erf}\!\left(W/(2\sqrt2\,\sigma)\right)\) of the line, i.e. 72–92 % for the approved Ly\(\alpha\) windows.

8.2 Significance

A line is detected when \(F/\sigma_F \ge 3\).

8.3 Equivalent width

where \(F_c(\mu) = c_0 + c_1(\mu - \bar\lambda)\) is the fitted local continuum at the line centre.^[11]

8.4 FWHM

with \(\text{FWHM}_\text{inst} = 6/(1+z) = 5.18\) Å (IUE SWP low dispersion, rest frame).^[14,21] Both observed and intrinsic widths are reported.

8.5 Reduced chi-squared

\(n_\text{par}\) counts every fitted parameter (3 per Gaussian + 2 continuum terms). \(\chi^2_\text{red} \approx 1\) means the residuals match the pixel noise; a low value is not penalised. Its formal significance is tested in §9.

Quality Evaluation & Acceptance Thresholds

A fit is classified as ACCEPTED, MARGINAL, or REJECTED based on three criteria:

A low \(\chi^2_\text{red}\) is not penalised. Detection uses the flux significance \(F/\sigma_F\) of the least-squares fit (Eq. 9).

9.1 Verification step

Every refined fit is checked with standard statistical tests. A fit is VERIFIED when all required checks pass; otherwise it is flagged CHECK together with the failed tests. A failed check does not discard the measurement; it states that a single Gaussian does not fully describe that profile.

\(\rho\) is the lag-1 autocorrelation of the pixel noise, measured per spectrum in the line-free continuum windows (typically \(\rho \approx 0.5\), because the MXLO spectra are resampled). For AR(1) noise the variance of \(\chi^2\) is inflated by \(c\); comparing \(\chi^2/c\) with \(\nu/c\) degrees of freedom is Satterthwaite's moment-matching approximation.^[20] \(F_\text{direct}\) is the trapezoidal integral of data minus fitted continuum and companion lines.

9.2 Diagnostic graphs

1. Rest-frame spectrum: observed flux, power-law continuum, shaded continuum windows, fitted line centres.
2. Line fit: data with ±1σ pixel-noise error bars, local continuum, each Gaussian component and the total model; masked pixels shaded.
3. Normalised residuals \((F - M)/\sigma\) with ±1σ and ±2σ reference lines.
4. Residual distribution: histogram of normalised residuals with the standard normal N(0,1).

The same panels are written for every line to a publication figure (PNG) per spectrum.

9.3 Profile Diagnostic Viewer (/ and /v2)

This is the home page. For Ly\(\alpha\) and C IV the line-fit graph shows the broad component (solid, shaded) and the narrow component (dashed) separately, and the parameter table adds the broad and narrow fluxes and FWHMs. The old single-spectrum fitter is kept at /v1 (footer link “Old version (v1)”).

After a spectrum is loaded the page shows, from top to bottom:

1. Header: file name, number of detected and verified lines, and the continuum \(\alpha \pm \sigma_\alpha\), \(A \pm \sigma_A\) (one per spectrum); buttons for the Excel export and the all-lines PNG.
2. Graph 1 and the emission-line summary (window, min/max λ, \(\chi^2_\text{red}\), status per line).
3. For the selected line: Graphs 2–4 and the parameter table (all values ±1σ).
4. Collapsible sections, each showing its key result while closed: Verification (§9.1); Graph 5 — variability context, which dates the spectrum from the SWP log by its SWP number and plots the selected line's detected fluxes of the same calendar year with this spectrum highlighted, plus a table of this spectrum's flux vs. that year's \(\bar F\) and \(F_\text{var}\) for every line; and Graph 6 — full-spectrum decomposition, the observed spectrum with every selected line's Gaussian in its own colour (Okabe & Ito colour-blind-safe palette), followed by a fit and light-curve panel per selected line. A line selector chooses which lines are shown and downloaded.
5. Campaign reference (collapsed): \(F_\text{var}\) per line per year for all dated spectra.

Every graph downloads as PNG (2× resolution) or SVG.

9.3a Campaign Results (/v5)

/v5 shows every line of every spectrum, fitted with the method of §7 by run_variability.py in parallel worker processes. A line selector drives all views: Graph A (all selected lines over the campaign as \(F/\bar F\)), a spectrum browser with the full-spectrum decomposition and per-line fit and residual panels, and a line-by-line results table. The campaign workbook adds a Flux matrix sheet (one row per spectrum, all lines side by side). Uploaded batches of .txt spectra are fitted the same way and dated by their SWP image number.

9.3b Emission-line results (/results)

/results shows the precomputed campaign without any upload:

1. Spectrum picker (1–255, with ◀ ▶). It draws the spectrum in the style of the reference figures: the continuum-subtracted flux \(F - F_c\) with every fitted component drawn alone and named (Lya-N, Lya-B, NV, CIV-N, CIV-B, SIV, Ly-B, …). A table below gives that spectrum's values for every line.
2. Emission-line workbook, displayed like the Excel file: sheet tabs along the bottom, one tab per emission line, one row per spectrum. Clicking a row opens that spectrum above. The buttons download the same workbook by line or by spectrum (§9.5).

The table is served by /api/results/sheets, which reads outputs/line_fluxes.csv with the standard csv module so it also runs on the deployed server, where pandas is not installed.

9.4 Manuscript method (/v4)

A separate viewer, /v4, and batch run (run_variability.py --method manuscript, output in outputs/manuscript/) measure the lines exactly as described in the 3C 390.3 manuscript. Everything else (continuum, peak search, seed grid, line windows and ranges, verification, dating, \(F_\text{var}\)) is shared with the default method. The differences are:

Without a local continuum term, structure the power law does not follow (e.g. around Ly\(\alpha\) + N V) remains in the residuals; the verification panel reports this.

9.5 Excel outputs

Campaign Variability Statistics

Over 250 dated IUE SWP epochs spanning 1978–1996, the emission-line fluxes of 3C 273 varies on timescales of weeks to years, driven by changes in the accretion rate and BLR response. Two standard variability metrics are computed:^[3,12,13]

10.1 Fractional Variability Amplitude \(F_\text{var}\)

The fractional variability amplitude quantifies intrinsic variability relative to the measurement uncertainties:

where \(S^2 = \frac{1}{N-1}\sum(F_i - \bar{F})^2\) is the sample variance, \(\langle\sigma_\text{err}^2\rangle = \frac{1}{N}\sum\sigma_{F,i}^2\) is the mean squared measurement error, and \(\bar{F}\) is the mean flux. When \(S^2 \le \langle\sigma_\text{err}^2\rangle\) the scatter is consistent with the measurement errors alone and \(F_\text{var}\) is undefined ("not variable").

\(F_\text{var}\) is computed separately for every emission line and every calendar year, using the light curve \((\text{JD} - \text{JD}_\text{min},\ F,\ F_\text{err})\) of the detected (\(F/\sigma_F \ge 3\)) fluxes; observation dates come from the SWP log (jd.xlsx). Its uncertainty follows Vaughan et al. (2003) and Sukanya et al. (2018):

10.2 Flux Ratio \(R_\text{max}\)

The simplest variability measure is the ratio of maximum to minimum observed flux:

For the 3C 273 Ly\(\alpha\) campaign, typical values are \(F_\text{var} \approx 0.1\)–0.3 and \(R_\text{max} \approx 2\)–4, consistent with the moderate variability amplitude characteristic of radio-loud quasars in this luminosity range.

Worked Example: Spectrum 70 (SWP 35476)

Ly\(\alpha\) in 70dcdr2dswp35476mxlo.txt, fitted jointly with N V over 1180–1265 Å (51 degrees of freedom):

IUE SWP Instrument Characteristics

The Short-Wavelength Prime camera of the IUE satellite employed a Vidicon detector with a caesium telluride photocathode, covering 1150–1975 Å. In low-dispersion mode (the standard for monitoring campaigns), the spectral resolving power is \(R = \lambda/\Delta\lambda \approx 300\), corresponding to a resolution element of approximately 5–6 Å FWHM. Individual spectra are extracted and calibrated by the IUE reduction pipeline (NEWSIPS for post-1992 data, or the original IUESIPS for earlier data) to produce flux-calibrated 1-D spectra in absolute units.^[14]

Two numbers matter for this pipeline. (i) Resolution: ≈6 Å FWHM in the observed frame, 5.18 Å in these rest-frame spectra; it is the lower limit for line widths and is removed in quadrature from the reported intrinsic FWHM. (ii) Sampling: 1.676 Å per MXLO pixel observed, 1.448 Å rest frame. Neighbouring pixels are correlated (lag-1 \(\rho \approx 0.5\)), which is accounted for in the goodness-of-fit test (§9.1).^[14]

Sky background is negligible for bright point sources in the UV; scattered-light corrections are included in the NEWSIPS flux calibration. No additional background subtraction is applied by this pipeline.

Software Architecture

The web application exposes the pipeline through these pages: Profile Diagnostic Viewer (/ and /v2; one spectrum, all graphs and checks), Results (Excel) (/results; any of the 255 precomputed spectra and the emission-line workbook, §9.3b), Campaign Results (/v5; light curves, \(F_\text{var}\) and downloads), Manuscript Method (/v4), Multi-Epoch Ensemble (/v3; client-side viewer of the reference Excel) and this documentation. The original dark-theme fitter is kept at /v1, linked from the footer as “Old version (v1)”.

References