Tutorial 3: Surface Brightness

:::info Migrated from Old Wiki This tutorial has been enhanced with content from the old GLEE wiki. Content should be double-checked for accuracy. :::

In this tutorial, we will learn how to model the full surface brightness distribution of a lens system. We minimize the following $\chi^2$ expression:

$$\chi^2_{\mathrm{light}} = \sum_{x=1}^{N_x} \sum_{y=1}^{N_y} \frac{[I_{\mathrm{obs}}(x,y) - \mathrm{PSF} \otimes I_{\mathrm{mod}}(x,y)]^2}{\sigma(x,y)^2}$$

Where:

• $N_x$ and $N_y$ are the dimensions of the science image cutout.
• $I_{\mathrm{obs}}(x,y)$ is the observed intensity at pixel $(x,y)$.
• $I_{\mathrm{mod}}(x,y)$ is the modelled intensity.
• $\mathrm{PSF}$ is the point spread function.
• $\otimes$ denotes the convolution operator.
• $\sigma(x,y)$ is the uncertainty at pixel $(x,y)$, provided by the uncertainty map.

note

Only unmasked pixels (those with a value of 1 in the lens mask, as discussed in Tutorial 0: Preliminary Data) contribute to the $\chi^2$ calculation.

Step 1: Setting Up the Configuration File

To begin the modelling, configure the esources section, placed directly after the sources section. We start by modelling the brightest components—the quasars—using a PSF light profile. To accurately capture their light without masking them out, we supply a 1-pixel arc mask.

esources 1
    Dds/Ds    1.000000  exact:
    esource_ngy          35
    esource_ngx          240
    esource_dx           0.08000
    esource_data         path/to/sci.fits
    esource_err          path/to/errormap.fits
    esource_arcmask      path/to/1pixmask.fits
    esource_lensmask     path/to/lensmask.fits
    esource_mod_light    LensOnly
    esource_psf          path/to/psf.fits
    esource_sub_agn_psf  path/to/psf_agn.fits
    esource_sub_agn_psf_factor     1
    esource_sub_esr_psf  path/to/psf_esr.fits
    esource_sub_esr_psf_factor     1
    esource_regopt       SpecRegPrecSigFigOnce
    esource_reglampre    1
    esource_reglamnup    1000
    esource_regtype      curv
    esource_reglam       247.388
    esource_reglamlo     1e-05
    esource_reglamhi     100000
    esource_light  4
    psf
        4.160000  #x-coord   flat:3.36,4.96  label:x_qso_1   step:0.001
        3.930000  #y-coord   flat:3.13,4.73  label:y_qso_1   step:0.001
       80.000000  #amp       noprior:  min:0  step:0.01
    psf
        4.970000  #x-coord   flat:4.17,5.77  label:x_qso_2   step:0.001
        5.590000  #y-coord   flat:4.79,6.39  label:y_qso_2   step:0.001
       70.000000  #amp       noprior:  min:0  step:0.01
    psf
        4.090000  #x-coord   flat:3.29,4.89  label:x_qso_3   step:0.001
        5.310000  #y-coord   flat:4.51,6.11  label:y_qso_3   step:0.001
       30.000000  #amp       noprior:  min:0  step:0.01
    psf
        5.430000  #x-coord   flat:4.63,6.23  label:x_qso_4   step:0.001
        4.660000  #y-coord   flat:3.86,5.46  label:y_qso_4   step:0.001
       80.560000  #amp       noprior:  min:0  step:0.01

Run a siman/MCMC cycle until convergence is achieved. Once the quasar positions and amplitudes are well sampled, fix these parameters to facilitate the lens light modelling.

Step 2: Modelling the Lens Light

Next, we incorporate the central lens galaxy's light using two Sersic profiles. At this stage, we mask the extended arcs and quasars by pointing esource_arcmask to the real arc mask.

esources 1
    ...
    esource_arcmask      path/to/arcmask.fits
    ...
    esource_light  6
    psf
        4.164410  #x-coord   exact:  label:x_qso_1   step:0.01
        3.931569  #y-coord   exact:  label:y_qso_1   step:0.01
       84.381816  #amp       exact:  min:0  step:0.1
    ...[Other PSFs set to exact] ...
    sersic
        4.675963  #x-coord   link:lens_x  a:0,1,1
        4.798325  #y-coord   link:lens_y  a:0,1,1
        0.595251  #q         exact:  label:sersic0q  step:0.04
        2.820239  #PA        exact:  label:sersic0pa  step:0.1
      7.0691e-04  #amp       exact:  min:0  label:814amp    step:0.001
        2.741933  #r_eff     exact:  min:0  label:sersic0reff  step:0.1
        0.480960  #n_sersic  exact:  label:sersic0n  step:0.1
    sersic
        4.675963  #x-coord   link:lens_x  a:0,1,1
        4.798325  #y-coord   link:lens_y  a:0,1,1
        0.904249  #q         exact:  label:sersic1q  step:0.5
        2.549768  #PA        exact:  label:sersic1pa  step:0.1
        0.012060  #amp       link:814amp  a:0,17.06,1
        0.602792  #r_eff     exact:  min:0  label:sersic1reff  step:0.1
        4.289665  #n_sersic  exact:  label:sersic1n  step:0.1

Finally, run another siman/MCMC cycle until convergence is reached.

You can iterate several times over these steps to ensure parameters are well sampled. Once both the quasar and lens galaxy light components are properly modelled, you will have a robust surface brightness model suitable for further extended arc source reconstruction.

:::tip What's next? With a converged surface brightness model you can:

• Reconstruct the lensed source by expanding the esource_arcmask to include the full arc and adding a pixellated source grid.
• Run a full Bayesian analysis with MCMC to obtain posterior distributions on all lens model parameters.
• Combine with time delays (chi2type 128) for Hubble constant constraints.

See the GLEE Configfile reference for all available keywords. :::

---

Understanding Masks: Arcmask vs. Lensmask

GLEE uses two distinct mask types to separate lensed arc light from foreground lens light:

esource_arcmask

Marks pixels (value = 1) containing lensed arc emission to be reconstructed on the source plane.

Purpose: Identify which image-plane pixels map to the source grid.

Processing:

1. Subtract all esource_light profiles from observed data inside arcmask: $$I_{\text{arc},j} = I_{\text{data},j} - I_{\text{light profiles},j}$$
2. Use residual $I_{\text{arc},j}$ to reconstruct source surface brightness

:::tip Creating Arcmasks from DS9 Regions Use GLEE's getMaskInFitsFromDS9reg.csh script to convert DS9 polygon regions to FITS masks:

getMaskInFitsFromDS9reg.csh arcmask.reg 100 120 arcmask.fits

• arcmask.reg: DS9 region file (image coordinates, in pixels)
• 100, 120: Output FITS dimensions (nx, ny)
• Ensure region file has 4 header lines before the polygon (add #dummyline if needed) :::

esource_lensmask

Marks pixels (value = 1) containing foreground lens light to be modeled via esource_light profiles.

Purpose: Define where to fit lens galaxy light (outside the arcmask).

Processing:

• GLEE automatically removes overlap: lensmask_effective = lensmask - arcmask
• Fit lens light (Sérsic, Gaussian, etc.) to lensmask_effective pixels
• Compute $\chi^2_{\text{lens light}}$ only on these pixels

:::note Overlap is Allowed You can provide overlapping masks. GLEE will prioritize arcmask pixels for source reconstruction and remove them from lensmask. :::

---

Noise Modeling

Instead of a static uncertainty map, GLEE can model pixel noise dynamically based on intensity:

Power-Law Noise Model

esource_noise_model      powerlaw
    0.20      #scale      exact:
    1.00      #power      exact:
    0.03      #const      exact:

Formula:

$$\sigma_i^2 = \text{scale} \cdot (d_i)^{\text{power}} + \text{const}$$

where:

• $d_i$: observed intensity at pixel $i$
• First term: Astrophysical source noise (Poisson when power=1, scale=1 for counts)
• Second term: Background noise (read noise, sky background)

Calibration

1. Power: Set to 1.0 for Poisson noise
2. Const: Estimate background variance from an empty region: $\text{const} = \sigma_{\text{bkgd}}^2$
3. Scale: Adjust for image units
   • If image is in counts: scale = 1.0
   • If image is in counts/sec: scale = 1.0 / t_exp

Example (HE0435, $t_{\text{exp}} = 9940$ sec, units = counts/sec):

esource_noise_model      powerlaw
    1.006e-04  #scale   (= 1/9940)
    1.00       #power
    0.03       #const   (measured from background)

Noise Mask

Optionally exclude problematic regions (e.g., PSF core mismatch):

esource_noise_model_mask    mask_qso4.fits

Pixels in this mask get artificially high $\sigma$ → effectively ignored in $\chi^2$.

:::warning Priority If esource_noise_model is specified, it overrides any esource_err or esource_wht file. :::

---

PSF Correction (Advanced)

GLEE can simultaneously optimize PSF corrections during source reconstruction, correcting for systematic PSF errors (e.g., pixelation, undersampling, wings).

Syntax

esource_psfc_ngx     9
esource_psfc_dx      0.080000
esource_psfc_regtype      grad
esource_psfc_reglam       2.02951e+06

Keyword	Description	Notes
esource_psfc_ngx	Size of PSF correction grid (pixels)	Must be odd (e.g., 7, 9, 11)
esource_psfc_dx	Pixel scale of correction grid	Defaults to image pixel scale if omitted
esource_psfc_regtype	Regularization type	grad (gradient) or curv (curvature)
esource_psfc_reglam	Regularization strength	Tune to avoid overfitting

How It Works

1. PSF correction grid $\mathbf{P}_c$ is optimized on pixels inside arcmask
2. Updated PSF = Original PSF + PSF correction (interpolated and renormalized)
3. AGN amplitudes automatically rescaled to match updated PSF

:::caution Initial Guess Required PSF correction requires a good starting regularization constant. If your initial guess is poor:

1. Run glee -f 2 iteratively to optimize esource_reglam
2. Watch for warnings about residuals being too large
3. Once stable, enable PSF correction :::

:::note Correction Scope PSF corrections only affect PSF light blobs (AGN, quasars). Extended source PSF (esource_sub_esr_psf) uses the central part of the corrected PSF. :::

---

Regularization Options

Regularization prevents source reconstruction from overfitting noise. GLEE supports several schemes:

Regularization Types (esource_regtype)

Type	Formula	Use Case
curv	Curvature: $R = \sum (\nabla^2 s)^2$	Recommended for smooth sources (galaxies)
grad	Gradient: $R = \sum (\nabla s)^2$	Compact sources (quasars)
quad	Quadratic: $R = \sum s^2$	Suppress spurious emission

Regularization Optimization (esource_regopt)

Option	Description	When to Use
SpecRegVal	Fix $\lambda$ to specified value	When optimal $\lambda$ is known
SpecRegPrecSigFig	Optimize $\lambda$ to match precision	Explore $\lambda$ every iteration
SpecRegPrecSigFigOnce	Optimize $\lambda$ once every N steps	Recommended (balance speed & accuracy)

Configuration:

esource_regopt       SpecRegPrecSigFigOnce
esource_reglampre    1          # Precision: 1 significant figure
esource_reglamnup    1000       # Optimize every 1000 MCMC steps
esource_reglam       2.5e-05    # Initial guess
esource_reglamlo     1e-07      # Search lower bound
esource_reglamhi     0.001      # Search upper bound
esource_regfact      10         # Over-regularization factor

Tuning Regularization

• Too low $\lambda$ → Noise amplification (checkerboard source)
• Too high $\lambda$ → Over-smoothing (loss of structure)

Best practice: Use SpecRegPrecSigFigOnce with broad search bounds (reglamlo, reglamhi) and let GLEE find the optimal value.