DRAFT September 29, 2000

The Imaging Capabilities of HESSI

Gordon J. Hurford, Edward J. Schmahl, Andrew Conway, Brian R. Dennis, Richard A. Schwartz, Lab for Astronomy and Solar Physics, NASA Goddard Space Flight Center, Greenbelt, MD 20771

Robert.P.Lin, Jim McTiernan, David Smith, Chris Johns-Krull, and Andre Csillaghy, Space Science Lab, University of California, Berkeley , CA 94720

Markus J. Aschwanden and Thomas R. Metcalf, Lockheed Martin Advanced Technology Center, Solar & Astrophysics Laboratory, Dept. L9-41 Palo Alto, CA 94304

Changes included: Completed references GH: Ref to HXIS ES: Subsection on Count rates and visibilities ES: Added more on statistics ES: Updated decimation scheme (thanks, DC) MA: Make definitions consistent

TBD: GH: Figure for various slit-slat ratios and energy dep’ce GH: Change (4) to ”harmonic motion” GH: More on statistics MA: Show simulation results for Bproj, CLEAN, MEM, FwdFit

Contents

1 Introduction 1 1.1 The HESSI Imager . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Expected Imaging Capabilities . . . . . . . . . . . . . . . . . 2

2 Modulation Principles 4 2.1 One-dimensional Modulation . . . . . . . . . . . . . . . . . . 4 2.2 Rotational Modulation . . . . . . . . . . . . . . . . . . . . . . 8 2.3 The HESSI Event List . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Visibilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Count Rates and Visibilities . . . . . . . . . . . . . . . . . . . 10

3 Image Reconstruction 13 3.1 Linear Methods: Backprojection and Fourier Transforms … 13 3.2 Nonlinear Reconstruction Methods . . . . . . . . . . . . . . . 14

3.2.1 CLEAN (Cartesian, Polar) . . . . . . . . . . . . . . . 16 3.2.2 Maximum Entropy Methods (MEM Sato, Polar MEM) 17 3.2.3 Forward Fitting (pixelized, non-pixelized) . . . . . . . 19 3.2.4 The PIXON Method . . . . . . . . . . . . . . . . . . . 20

3.3 Undeveloped Techniques . . . . . . . . . . . . . . . . . . . . . 22 3.3.1 Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3.2 Non-Negative Least Squares . . . . . . . . . . . . . . . 22 3.3.3 Two-Dimensional Autocorrelation . . . . . . . . . . . 22

4 Statistical Limits 22 4.1 Fourier Coefficient Distribution . . . . . . . . . . . . . . . . . 23 4.2 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . 24

5 Figure Captions 26

6 Appendix 29

Abstract

The High Energy Solar Spectroscopic Imager (HESSI) will use rotational modulation synthesis for imaging hard X-ray and gamma-ray flares with spatial resolution of 2.3′′ and spectral resolution of ∼ 1 keV. Like the Yohkoh/HXT hard X-ray telescope, HESSI relies on Fourier methods to produce images, but HESSI has many more sampled points in the Fourier plane, and is expected to produce higher-resolution maps with greater dynamic range. We summarize the methods to be used for HESSI imaging, touching on crucial details of modulation such as the Calibrated Event List, and describing the Back Projection method. Beyond these, we outline four basic tools for image reconstruction: CLEAN, MEM, PIXONS, and Forward Fitting, with simulated examples of each.

1 Introduction

The only practical method of obtaining ∼arc-second angular resolution in hard X-rays and gamma-rays within the cost, mass, and launch constraints of a small satellite is to use Fourier-transform imaging. (See Prince et al. 1988 for a review of imaging techniques.) One of the most powerful of the Fourier family techniques is rotational modulation synthesis, first proposed by Mertz (1967) and implemented by Schnopper et al. (1968). Rotating Modulation Collimators (RMC) were first constructed and used for solar flare X-ray imaging by the Japanese Hinotori telescope (Makashima et al. 1977, Ohki et al. 1982, Enome 1982), with angular resolution of 28′′ in the 20-40 keV energy range. The next hard x-ray solar imager was the Hard X-ray Imaging Spectrometer (HXIS), not a Fourier-transform imager, but important for developments in the ∼ 3 − 30 keV range (Van Beek et al. 1980). A later Japanese telescope, HXT, on the Yohkoh satellite, used nonrotating Fourier synthesis with angular resolution of ∼ 8′′ in the 20-100 keV energy range (Kosugi et al. 1991). A prototype of a rotation modulation telescope—the High Energy Imaging Device (HEIDI)— was flown in 1993 on a balloon, without successful flare imaging, but with a definitive test of a novel solar aspect system (SAS), proving 0.5-arcsecond performance at balloon altitudes. (Crannell et al. 1994)

The High Energy Solar Spectroscopic Imager (HESSI) also uses rotational modulation synthesis imaging. With a launch shortly after the peak of flare activity, HESSI has been designed to image hard X-ray and gamma-ray flares with spatial resolution of 2.3′′ and spectral resolution of ∼ 1 keV (Lin et al. 1993, 1994, 1998, Holman et al. 1997). We summarize briefly the HESSI instrumentation used to exploit Fourier synthesis imaging.

1.1 The HESSI Imager

HESSI uses nine bi-grid collimators, each consisting of a pair of widely separated grids in front of an X-ray/gamma-ray detector. Each grid consists of a planar array of equally-spaced, X-ray-opaque slats separated by transparent slits. The slits of each pair of grids are parallel to each other and their pitches (p) are identical, so the transmission through the grid pair depends on the direction of the incident X-rays. (See Figure 1.1) For slits and slats of equal width, the transmission is modulated from zero to 50% and back to zero for a change in source angle to collimator axis (orthogonal to the slits) of p/L where L is the separation between grids (1550 mm). The angular resolution is then defined as p/(2L). For HESSI, the transmission of the source photons through the grids is modulated by the rotation of the spacecraft at about 15 revolutions per minute. The detector records the arrival time and energy of individual photons from anywhere on the Sun, allowing the modulated counting rate to be determined as a function of rotation angle.

Note that the detectors have no spatial resolution and hence have been optimized for high sensitivity and energy resolution. The nine segmented Ge detectors (GeDs), one behind each RMC, detect photons from ∼3 keV to 15 MeV. The GeDs are cooled to ∼ 75 K by a space-qualified long-life mechanical cryocooler. As the spacecraft rotates, the RMCs convert the spatial information from the source into temporal modulation of the photon counting rates of the GeDs. The instrument electronics amplify, shape, and digitize the GeD signals, provide power, format the data, and interface to the spacecraft electronics.

HESSI exploits one of the fundamental efficiencies of Fourier imaging: precision spacecraft pointing can be traded for exact knowledge of pointing. With a Fourier imager, it is not necessary to stabilize the spin axis to any better than a few arc minutes, as long as the system gives sufficiently precise and high-bandwidth pointing information. For HESSI, this information is provided by the solar aspect system (SAS) — a heritage of the HEIDI SAS — and the roll angle system (RAS). The SAS consists of 3 linear diode arrays on which the full solar image is projected. The solar limb is determined at 6 positions every 10 ms, giving pitch and yaw to ∼ 1.5′′. The RAS consists of a star-imaging linear photodiode, providing roll angle to 2.7′ every minute.

1.2 Expected Imaging Capabilities

In 2-3 years of mission life, HESSI is expected to obtain observations of tens of thousands of microflares, thousands of hard x-ray flares, and of order a hundred gamma-ray line flares. Detected photons will be tagged with 1 μs time resolution up to a rate of ∼ 10^{5} photons/s. The intense 3-150 keV Xray fluxes that usually accompany large gamma-ray line flares are absorbed by the front segment of the detectors, so the rear segment will always count at moderate rates. This is essential for gamma-ray line measurements with optimal spectral resolution and high throughput.

Except in extreme circumstances, every photon is stored and then sent down. If an extended, extremely active period occurs, a decimation scheme is enabled which digitally decimates the incoming photon event stream, so that only a fraction of the events below an energy threshold get sent to telemetry.

3 Both the fraction and the energy threshold are programmable. There are decimators for both the front and the rear segments of the detectors. The front decimator thresholds are controlled by software to avoid overflowing the spacecraft recorders, and the rear decimators can be controlled from the ground to limit the background rate.

In addition to decimation, there are mechanical attenuators, which if required, reduce the flux below the saturation level of ∼ 50, 000 events/s. Together, decimation and attenuation provide minimum loss of imaging and no loss in gamma-ray spectroscopy.

It is expected that HESSI will provide detailed imaging spectroscopy in each of ten energy intervals with ∼ 2s time resolution for one event every ∼ 5 days with ∼ 10^{4} counts s^{−1} above 15 keV. With lower energy resolution (∼ 100 keV), crude imaging information could be obtained in tens of milliseconds. HESSI provides spatial resolution of 2.3 arcseconds at X-ray energies below ∼ 200 keV, 7 arcseconds to 400 keV, and 36 arcseconds for gamma-ray lines and continuum above 1 MeV.

2 Modulation Principles

HESSI, like all RMCs, relies on temporal modulation by its grids to provide spatial information about incoming photons. The imaging process may be thought of as a superposition of elemental processes of photon passage through a sub-collimator. Since the time of arrival of each photon is determined to 1 μs, and the aspect system determines the pointing of the telescope with arcsecond accuracy, the known geometry of the subcollimators provides well-defined information about the photon’s origin.

2.1 One-dimensional Modulation

We illustrate how this works by considering a subcollimator with two identical grids, each parallel to the x,y plane, whose slats are, at a given instant, parallel to the y axis, and whose spin lies parallel to the z axis. Let a photon travelling in the x-z plane strike the plane of the front grid. The incidence phase (measured in the x direction) of the photon relative to the slats of the front grid may be considered to be random since all photons arriving from a single point at infinity lie on parallel lines with a continuum of phases, and since the detectors have no spatial resolution, the modulation profile is found by integrating over all incidence phases.

The probability of passage of the photon through the front grid is proportional to the projected area of the slits of that grid, as viewed from the direction of the photon. For an ideal grid of finite thickness and infinite extent, this value is independent of the y direction, but dependent on the photon direction in a fashion characteristic of the grid thickness and the slit-slat ratio.

In the case of a single thin grid with equal slits and slats, the projected area, plotted as a function of incidence phase, is a periodic, rectangular waveform with equally long maxima and minima. For a thick grid, there is internal shadowing by the sides of the slats, and the projected open areas are reduced from the thin case, so the area waveform is still periodic, but with maxima and minima of unequal lengths. In the case of HESSI grids, the projected slit-slat ratios are optimized to have equal projected slit and slat widths at incidence angles of about 600′′. The aspect ratio of the slats (x-width/z-thickness) is 1:50. For photons with nearly normal incidence, the effective slit/pitch ratio is about 0.6, and for photons incident at angles of 415′′ in the x-z plane, the effective slit/pitch ratio after shadowing is close to 1:2.

After photons pass through the front grid, they have a finite probability of passing through the rear grid. If the incidence angle is less than about 900′′, the probability of incidence on the rear grid is almost unity. Only for photons entering at angles of greater than 1 degree does this probability become significantly less than 1 due to non-overlap of the front and rear grid areas. The passage of the photon through the second grid is, again, dependent on its direction in exactly in the same way as for the front grid. The fraction of projected open area for an incident photon is therefore the convolution of the transmission functions of each grid. It is easy to verify that in the ideal case, this is a triangle waveform. For thin grids, the triangle function has pointed peaks and valleys. But if the projected slit-slat ratio is larger or smaller than 1, the valleys of the waveforms will be flattened. The sharp corners of the profile’s peaks are also rounded by partial transmission effects. Figure 2 shows examples of the equivalent area function for various slit/pitch ratios, and also for (over-simplified) energy-dependent effects that smooth the corners of the profile.

It is useful and efficient to characterize the subcollimator transmission functions in terms of the first few harmonics. For slit/pitch ratios equal to 1/2, and for energies low enough that the slats are opaque (E < 100 keV), and for energies high enough that diffraction effects are negligible (E > 6 keV), the effective subcollimator area function is approximately a triangle waveform, 6 represented by the following function of θ, the photon’s angle of incidence relative to the imager axis:

S_{triangle}(θ) = \frac{A}{2} · (\frac{1}{2} + \frac{4}{π^{2}} cos(θ) + \frac{4}{9π^{2}} cos(3θ) + \frac{4}{25π^{2}} cos(5θ) + ⋯) (1)

This is precisely the form used by Schnopper et al. (1970) for their observations of galactic-center with a rocket-borne RMC. But for HESSI, the even harmonic terms are usually not negligible, and the Fourier coefficients are not those of a pure triangle waveform, so we must use the more general form:

S(θ) = A·(c_{0} +c_{1} cos[θ −θ_{1}]+c_{2} cos(2[θ −θ_{2}])+c_{3} cos(3[θ −θ_{3}])+⋯) (2)

where c_{0}A is the effective area averaged over θ. In practice, we keep terms only up to the 5th harmonic. In general, the collimator transmission function is the convolution of the transmission functions of the front and rear grids. If the front and rear grids of a subcollimator are geometrically identical (which is very nearly true for HESSI), this means that the Fourier coefficients of S(θ) are the squares of the transmission functions of a single grid. For the case of ideal grids of effective slit/pitch ratio q, the Fourier coefficients of the grid transmission are

g_{n} = \frac{1}{2π} ∫_{−qπ}^{qπ} cos(nθ)dθ = sin[nπq]/(nπq)  n = 1, 2, 3, ⋯  (3)

So the amplitudes c_{n} are:

c_{n} = 4sin^{2}[nπq]/(nπq)^{2}  n = 1, 2, 3, ⋯  (4)

When q = 1/2, the case of a simple triangle profile (which occurs for HESSI when θ = 412′′), all the even amplitudes vanish, but for q = 0.4, (photons at normal incidence on the HESSI grids), c_{2} = 0.1093, and for q = 0.6 (equivalent to angular incidence of 824′′), c_{2} = 0.0486. In general, for real conditions, the relative amplitudes c_{1}, c_{2}, ⋯ and the relative phases θ_{1}, θ_{2}, ⋯ are dependent on the grid characteristics, the roll angle, and the angle of incidence. The latter dependence is sufficiently small and slow that for sources on the Sun, that they can be considered uniform within a map in a given time bin. These parameters have been computed from the grid characterizations, and are included in the HESSI reconstruction algorithms.

2.2 Rotational Modulation

So far we have described the transmission only for one roll angle, namely the instant at which the pattern projected on the sky from an ensemble of sight lines from the detectors — the so-called “modulation pattern” — has contours of constant phase perpendicular to the vector on the plane of the sky from the source to the projected telescope axis. As the spacecraft spins, the projected modulation pattern spins with it. To the lowest order of approximation, if the grids are infinitely wide and infinitely thin, and the spin axis does not move, the modulation pattern would rotate unchanged. Since, however, the grids are thick, of finite extent, and the spin axis is not necessarily constant during a rotation, the projected modulation pattern changes with roll angle in a predictable way, which is reflected in a slow, smooth change of the phase and amplitude coefficients of the sub-collimator pattern.

It is useful to take the telescope as our reference frame, and let a point source move in a circle (roll angle α = 0−2π) on the celestial sphere around the spin axis. Then, because no modulation occurs in the y direction parallel to the slats, and modulation in the x-direction follows equation (2), the effective pitch changing in proportion to sec(α). During a rotation, the fundamental coefficients (c_{1}, θ_{1}) change least, and the function S is a simple function of roll angle α. For an ideal RMC (idealized grids, steady spin), the argument θ in equation (2) is proportional to the radial distance R (arcsec) of the source and inversely proportional to angular pitch p):

θ = 2πRcos(α)/p  (5)

The flux of photons on a HESSI detector from a point source a distance R off axis is proportional to the effective area, which leads to the modulation profile as a function of rotation angle α:

P_{0}(α) = A · (1 + c_{1}cos(2πRcos(α)/p − θ_{1}) + ⋯)  (6)

where p is the angular pitch (arc sec) of the grids in the subcollimator. For a given map center, the quantities (c_{n}, θ_{n}) are computed from the aspect system and the grid response matrices, and are returned in the HESSI calibrated event list. The above function is computable for any time interval and any point source on the Sun. Therefore, since the the signal expected for a superposition of sources is the superposition of point source signals, any model source can be used to predict a modulation profile. This is the basis for all HESSI reconstruction algorithms.

2.3 The HESSI Event List

The set of photon arrival times taken from the HESSI telemetry can be histogrammed to yield a raw count rate profile. The HESSI user can select an almost arbitrary set of time bin sizes {Δt_{i} } for histogramming. (Δt_{i} is necessarily ≥ 1 μs, and may be as large as ∼ 1 s, but is more typically ∼ 0.5 ms). After the time bin selection, the count/bin profile becomes a function of roll angle α, and one can proceed to calibrate the rates.

For phase and amplitude calibration, one must select a phase center for mapping. This point can be found in several ways: (a) determine the azimuth by the roll angle of fastest modulation, and the radius by the number of modulation cycles per rotation; (b) make a rough map of a large region (say a solar quadrant) and find the position of the maximum; or (c) use context data or (d) previous flare position data.

Using the aspect system data, the distance of the line of maximum transmission from map center is computed for each roll angle α_{i} and each subcollimator pitch p_{k}. This gives the phase at map center Φ_{ik}. (See figure (??.) For the case of a fixed spin axis, Φ_{ik} = 2πRcos(α_{i})/p_{k} − θ_{k}. But in practice, Φ_{ik} is computed without reference to R, using aspect data and the subcollimator geometries alone. One must incorporate the dead time, a function τ_{ij} of roll angle and energy E_{j}, the detector number k, and the subcollimator transmission T_{ijk}, as well as the modulation amplitudes c_{n}^{ijk} and modulation phases Φ_{n}^{ijk} to determine the model modulation profile:

S_{ijk} = F_{0} · {τ_{ijk}T_{ijk}[1 + ∑_{n} c_{n}^{ijk}cos[n(Φ_{ik} − Φ_{n}^{ijk})]}  (7)

When the HESSI analysis software computes the quantities τ_{ijk}, T_{ijk}, c_{n}^{ijk} and Φ_{n}^{ijk}, they are returned in a structure called the calibrated event list, the starting point for image reconstruction.

2.4 Visibilities

In radio astronomy, the quantity analogous to the modulation pattern is the fringe pattern, which is the response of an interferometer to point sources in the sky. This is the Fourier transform, and it is represented by a sum of complex exponentials in the u,v plane:

V_{i} = f e^{2πi(u_{i}x+v_{i}y)}  (8)

for each point source of flux f at coordinate (x, y). The variables u_{i}, v_{i} are the coordinates in the Fourier plane, and for HESSI, u_{i} = cosϕ_{i}/p and v_{i} = sinϕ_{i}/p, so the points (u_{i}, v_{i}) lie on circles in the Fourier plane. The importance of visibilities to HESSI is that they represent an intermediate between count-rate profiles and images that is independent of the pointing. The visibilities have all of the spin-axis wobble, slit-shadowing effects, transmission and deadtime accounted for. Visibilities, unlike countrate profiles, can be added from 1/2 rotation to another.

A significant simplification for imaging occurs when one changes from the usual Cartesian coordinates to polar coordinates in both the image plane and the Fourier plane:

x = r cosθ  u = k cosϕ
y = r sinθ  v = k sinϕ

Since the visibilities are found on circles in the Fourier plane, the radial integration drops out (V is proportional to a delta function of radius k), and after a little algebra, the Fourier transform for one subcollimator of wavenumber k becomes:

F (r, θ) = k ∫_{0}^{2π} V(ϕ)e^{ikr cos(θ−ϕ)}dϕ = kV ⊗ e^{ikr cos(ϕ)}  (9)

The nice thing about this, is that for each radius, F (r, ϕ) is a convolution of the visibility and a complex exponential. Since convolutions are orders of magnitude faster than matrix multiplications, polar coordinate representations of HESSI data provide a great increase in speed of “backprojection” (section 3.1) or other reconstructions.

2.5 Count Rates and Visibilities

Given an incident flux of photons on a HESSI subcollimator, after the selection of a map center, the predicted count rate C is given by equation (7) with subscripts i, j, k dropped:

C = F_{0} T τ · {1 + c_{1} cos[Φ] + c_{2} cos[2(Φ − Φ_{2})] + ⋯}  (10)

where F_{0} is the incident photon flux on front grid, T is the subcollimator transmission, τ is the livetime, Φ is the phase at map center, c_{n} is the modulation amplitude, and Φ_{n} is the modulation phase. (All these and subsequent variables are defined in the Appendix.)

In general, the phase at map center depends on the results of the aspect solution, the time binning, the subcollimator, and the position of the map center. The phase at map center Φ depends on the the projection of the vector between the subcollimator axis (x_{coll}, y_{coll}) and the map center (x_{map}, y_{map}):

Φ = 2π{[x_{map} − x_{coll}]cosα + [y_{map} − y_{coll}]sinα)/p}  (11)

where

(x_{coll}, y_{coll}) = Subcollimator direction in inertial coordinates
(x_{map}, y_{map}) = Map center in inertial coordinates
α = Azimuth angle of slits
p = Pitch of subcollimator grids

In the simplest case, where, the spin axis is at the origin and is stationary, the phase at map center is a function only of the azimuth ϕ, the pitch p, and the location of the map center (R_{0}, ϕ_{0}):

Φ = 2πR_{0}cos(ϕ − ϕ_{0})/p  (12)

But in general, regardless of the position of the spin axis, the predicted count rate is given by equation (10). The visibility of a point source at the map center is:

V = F_{0} e^{i2π[x_{map} cosα+y_{map} sinα)]/p}  (13)

Note that the visibility does not depend on the subcollimator coordinates (x_{coll}, y_{coll}), but the count rate is subject to phase shifts caused by telescope motion. The phase shifts due to this motion are:

ΔΦ = 2π{x_{coll} cosα + y_{coll} sinα)/p  (14)

And these are reflected in additional variations of the count rate profile. (The quantity ΔΦ is computed from the aspect solution.) We can separate out this variation explicitly in the visibility:

V = F_{0} e^{i[Φ+ΔΦ]}  (15)

If we subtract out the mean of the count rate (equation 2), and divide out the T τ factor, keeping only the fundamental term,

(C − F_{0} T τ)/(T τc_{1}) = F_{0} cosΦ  (16)

So, for this particular case of a point source at map center the relation between the countrate profile and the visibilities are found by comparing equations 15 and 16:

Re\{V e^{−iΔΦ}\} = \frac{C − F_{0}T τ}{c_{1}T τ} = \frac{C− < C >}{c_{1}T τ}  (17)

It is important to note that, by equation 17, the peak-to-peak variation of the visibility is several times larger than the peak-to-peak variation of C. This is because the visibility reflects the transform of the incident photons, while the count rate is the modulated profile of only the detected photons. In fact, when the livetime is unity (T = 1), and when the subcollimator pattern is perfectly triangular, (c_{1}τ = 2/π^{2}), the peak-to-peak variation of Re{V } is π^{2}/2 = 4.9 times that of the count profile.

For a more general case, where the source is an arbitrary superposition of N point sources of flux F_{j} and location (R_{j}, ϕ_{j}), the visibility is:

V = e^{iΔΦ} ∑_{j=1}^{N} F_{j} e^{iΦ_{j}}  (18)

where the Φ_{j} are the phases of the point sources at (R_{j}, ϕ_{j}). Then, the corresponding predicted count rate is:

C = T τ · ∑_{j=1}^{N} F_{j} {1 + c_{1} cos[Φ_{j})] + ⋯}  (19)

Comparing equations 18 and 19, we find that equation 17 remains valid, providing that F_{0} is defined as the total flux of the N sources.

The generalization to higher harmonics is straightforward. The visibility for the n^{th} harmonic is given by:

Re\{V^{(n)}e^{−iΔΦ}\} = \frac{C^{(n)}}{c_{n}T τ}  (20)

where C^{(n)} is the n^{th} term in the Fourier expansion of the count rate (equation 10). Note that equation (20) determines only a projection of the complex visibility at each instant of time. Since, however, each modulation cycle contains phase as well as amplitude information, equation (20) can be used to determine both components for each modulation cycle. Each visibility profile so derived from equation 20 then defines a portion of the entire Fourier transform obtained on a circle of radius n/p in the Fourier plane.

3 Image Reconstruction

3.1 Linear Methods: Backprojection and Fourier Transforms

Back projection (Mertz, Nakano, and Kilner 1986) is the most straightforward and basic method of image reconstruction. It is equivalent to a Fourier transform (Kilner and Nakano, 1989.) A map constructed by this method is called the ”dirty map”, being the analogue of the radio astronomer’s initial Fourier transform of the observed visibilities. Back projection is a perfectly linear process: maps for arbitrary time intervals may be added together, and maps for different pitches and harmonics may be summed, generally leading to improvement of the image. Further improvements to the image by CLEAN or MEM (for example) do not share this property of linearity.

The principle of back projection is to consider the quantities S_{ijk} in equation (7) as circular slices of the two-dimensional modulation patterns P_{im}. The modulation profile (7) may be thought of as the amplitude at an inertial position of a series of two-dimensional modulation patterns that rotate around at the same rate as the telescope. The modulation patterns must shift in phase and amplitude, depending on the relative position of the telescope axis as a function of roll angle.

The brightness D_{m} of each pixel (m) in the back-projection map is defined by the following linear combination of the count rates O_{i}:

D_{m} = ∑ O_{i}W_{im}  (21)

where the weights W_{i} are linearly related to the modulation patterns P_{im}:

W_{im} = (P_{im}− < P_{m} >)/(< P_{m}^{2} > − < P_{m} >^{2}) · T_{i}τ_{i},  (22)

and the averages < P_{m} > and < P_{m}^{2} > are taken over the range of roll angles α_{i}. This particular normalization of the backprojection weights W_{i} was first derived by Durouchoux et al (1983). The explanation for the normalization of the weights W_{i} derives from the fact that the predicted count rate O_{i} for a point source of strength S at sky position m is:

O_{i} = S T_{i}τ_{i} P_{im}  (23)

Inserting the last equation (23) into the equation for the dirty map (21) shows that the peak of the dirty map equals the strength of the point source,

D_{m} = S,  (24)

since ∑ W_{im}T_{i}τ_{i} P_{im} = 1. Although (24) is strictly true only for a single source, simulations show that it produces nearly optimum back-projection maps for more complicated sources.

Figure 3.1 shows an example of a back-projection image for a simulated double source (Gaussian widths = 1 × 1′′) and the corresponding countrate profiles.

Although it is possible to use only one modulation pattern per roll angle to produce back-projection images, it turns out to be more efficient to construct ”universal” modulation patterns in the form of cosine and sine pairs. Sinusoidal modulation patterns of arbitrary phase can then be computed by appropriate linear combinations. This formulation has the advantage that neighboring half rotations may use the same ”universal” set of patterns, even if there are different aspect solution profiles.

It is worth mentioning that in early versions of RMC design, there were both sine and cosine subcollimators (e.g., Mertz et al. 1986), Murphy, 1990), but this turns out not to be necessary in hardware, since with sufficient sampling of the modulation pattern, relative amplitudes of the sine and cosine components can be determined from the modulation profile itself.

Each modulation pattern has contours of constant phase perpendicular to the roll angle azimuth vector, and has the periods and amplitudes used in equation (7). The modulation patterns so computed for each angle bin, are added together with weights given by the count rate divided by the product of the transmission factors and deadtimes. The sum is the backprojection map. For HESSI, simulations suggest that summation of the 9-subcollimator back-projection maps gives a dynamic range of order 10:1 for rates ∼ 10^{5}counts^{−1}s^{−1}.

3.2 Nonlinear Reconstruction Methods

All of the practical image reconstruction algorithms for improving a backprojection (”dirty”) map are nonlinear, since linear deconvolution algorithms such as Wiener filtering and inverse filtering are inapplicable to applications with incomplete sampling of the Fourier plane, as is the case for HESSI modulation synthesis. A number of reconstruction methods already exist, or are being improved, or are being considered for addition to the HESSI reconstruction suite. We summarize several of these here.

15 3.2.1 CLEAN (Cartesian, Polar)

CLEAN is an iterative algorithm which deconvolves a Point Spread Function— the imager’s response to a delta function source—from a “dirty map”. This algorithm is of fundamental importance in radio astronomy, where it is used to create images of astronomical sources obtained with interferometers. The basic CLEAN method was developed by Högbom (1974). It was originally designed for point sources, but it has been found to work well for extended sources as well when given a reasonable starting model. The Högbom CLEAN constructs discrete approximations to the CLEAN Map in the plane from the convolution equation

P ⊗ I_{source} = D  (25)

where P is the HESSI PSF for one or many subcollimators and/or harmonics, D is the dirty map, and ⊗ denotes a convolution. The CLEAN algorithm starts with an initial approximation I_{0} = 0 to the residual map. At the n^{th} iteration, it then searches for the largest value in the residual map. A delta function is then centered at the location of the largest residual flux and given an amplitude μ (the so-called “loop gain”) times this value. P is then subtracted from I_{n−1} to yield I_{n}.

I_{n} = I_{n−1} − μP  (26)

Iteration continues until a specified iteration limit is reached, or until the peak residual or RMS residual decreases to some level. The resulting final map is denoted I_{final}, and the position of each delta function is saved in a “CLEAN component” table. At the point where component subtraction is stopped, it is assumed that the residual brightness distribution consists mainly of noise.

To damp out high spatial frequency features which may be spuriously created in the iteration, each CLEAN component is convolved with the so-called CLEAN PSF, P_{clean} (the “CLEAN beam” in radioastronomy), which is simply a suitably smoothed version or inner portion of the PSF, P. A CLEAN map is produced when the final residual map is added to the approximate solution,

I_{clean} = I_{final} ⊗ P_{clean} + {D − P ⊗ I_{final}}  (27)

in order to include the noise.

The main disadvantage of CLEAN is that it does not, at least in the Högbom version, compare the observed modulation profile with a model modulation profile, to assess the “goodness of fit” during the iteration. Two variants of CLEAN exist for HESSI. One is for rectangular coordinates, and the other is for polar coordinates. In most circumstances, the latter has a significant advantage of speed over the former.

Figure ?? shows an example simulated source, the “dirty map”, the “CLEAN map”, and the residuals.

3.2.2 Maximum Entropy Methods (MEM Sato, Polar MEM)

The principle of Maximum Entropy is to find the map which is maximally noncommittal with regard to the unavailable information (i.e. the regions of the Fourier plane for which HESSI has no data). The image found is that which fits the data to within the noise level and also maximizes the entropy. In one view, the entropy is something which, when maximized, produces a positive image with a compressed range in pixel values. (Cornwell 1984) Image entropy thus defined is therefore not to be confused with a “physical entropy”. Several kinds of ”entropy” have been proposed to achieve this, but according to Skilling (1984) and Sivia (1996) the only function which guarantees that no untoward correlation is imposed is:

H = − ∑_{j=1}^{M} f_{j} ln(f_{j})  (28)

The idea is to maximize the objective function Q:

Q = H − λ χ^{2} − β ∑_{1}^{M} f_{j}  (29)

for some values of λ, where λ is a parameter giving a measure of the tug-of-war between maximizing entropy and minimizing the χ^{2} statistic and β is a Lagrange parameter multiplying the flux constraint. At maximum,

f_{j} = F − j e^{−λ∂χ^{2}/∂f_{j}−α}  (30)

One may solve the equations (11) iteratively or by directed search of Q.

One of the HESSI MEM reconstruction algorithms is MEM Sato, derived from the Maximum Entropy Method used for HXT imaging (Sato, 1998, Sato, Kosugi, and Makishima, 1999).

18 Figure 3.2.2 shows an example simulated source, its count profile, the “MEM map”, and the residuals.

3.2.3 Forward Fitting (pixelized, non-pixelized)

The forward-fitting method is based on models that represent a spatial map by a superposition of multiple source structures, which are quantified by circular Gaussians (4 parameters per source), elliptical Gaussians (6 parameters per source), or curved ellipticals (7 parameters), designed to characterize real solar flare hard X-ray maps with a minimum number of geometric elements.

In one realization of this scheme, maps of superimposed Gaussians are created iteratively, at each step constructing a model modulation profile from the map. In this “pixelized” form, the map evolves into a best fitting image. A second form of forward-fitting is the unpixelized version, where the parameters of the Gaussians are used to find the “equivalent point source”, which is equivalent to an amplitude and phase at each roll angle. The “equivalent point source” is then used to create the model modulation profile for comparison with the observed count profile and continued iteration. The unpixelized form requires at least an order of magnitude less memory and computation than the pixelized form, but the programs are more complex and their robustness is as yet unknown.

3.2.4 The PIXON Method

The PIXON method is another technique removes the sidelobe pattern of a telescope while mitigating the problems of correlated residuals and spurious sources which are commonly seen in Fourier deconvolution, chi-square fitting, and Maximum Entropy approaches.

The goal of the Pixon method is to construct the simplest, i.e. smoothest, model for the image that would be consistent with the data, i.e. have an acceptable chi-square fit. Being the simplest model, the derived image would be artifact free, i.e. there would be no spurious sources. In addition, the model is necessarily be most tightly constrained by the data, and consequently have the most accurately determined parameters.

The PIXON method changes the global smoothness idea of MEM into a local condition in which local maximum smoothness of the image is imposed. From an information science point of view, one selects a model with the minimum information content from the family of multiresolution basis functions (pixons) and which statistically fit the data.

Since the model has minimum complexity, spurious sources are unlikely to arise. Each parameter is determined using a larger fraction of the data, hence it is determined more accurately. This usually results in superior photometric and positional accuracy. And since the minimum number of parameters are used, the data cannot be over fitted.

Figure 3.2.1 shows an example simulated source, its count profile, the “PIXON map”, and the residuals.

21

3.3 Undeveloped Techniques

The available HESSI reconstruction algorithms in no way span the universe of techniques. A few of the likely methods to be added to the existing suite of programs are outlined here.

3.3.1 Wavelets

HESSI has pitches at multiplicative intervals of \sqrt{3}, and so, in a way it is a wavelet telescope, where the wavelet scales are logarithmically spaced by \sqrt{3}. To date, no one has exploited wavelet mathematics for HESSI imaging.

3.3.2 Non-Negative Least Squares

A promising technique used for medical imaging (ref) is Non-Negative Least Squares (NNLS). In general, linear Least Squares is not expected to be useful for HESSI, but the enforcement of positivity in NNLS provides an essential nonlinear constraint that may lead to good images in the HESSI parameter domain.

3.3.3 Two-Dimensional Autocorrelation

Mertz (1996) has suggested an autocorrelation technique for imaging with RMCs. It has the desirable properties of positivity and smaller sidelobes and shows promise for fast mapping of single sources. However, as Mertz points out, difficulties in processing complex sources seem to limit its applicability. It may, however, be useful for centroid finding before final reconstruction.

4 Statistical Limits

Since HESSI is a photon-counting instrument, imaging involves stochastic processes. During the time during which the data for an image are acquired, the incident photon rate λ may be considered a constant. The number of photons per time bin is given by Poisson statistics, where the probability density function p_{poi}(λ) is the expected number of photons per time bin incident on the front grid. The process of modulation is also stochastic, since the incident phase relative to the collimator fringes is random, and the probability of a photon detection is proportional to the function (7) normalized to unity at its maximum. The counts per bin is therefore given by the product of these two probabilities:

C_{i} = p_{poi}(λ) S_{norm}  (31)

This is the ansatz used for simulations of HESSI count rate profiles (Schmahl, 1998). Examination of equation (7) for S_{ijk} shows that if its Fourier expansion is truncated at the fundamental (as is often very useful), there are situations when S can go negative. This is not really a problem for backprojection or for Fourier mapping, neither of which depend on positivity. But it is of importance for CLEAN, MEM, and Pixons. Of singular importance, if it is negative, the normalized version (S_{norm}) of equation (7) fails to have any meaning as a probability.

4.1 Fourier Coefficient Distribution

The problem arises due to an expectation that S_{ijk} should represent a count rate. If, however, the Fourier coefficients are considered as primary elements, then the statistics can be done with them, regardless of their number and whether they produce a negative sum.

Let us consider the fundamental only. After choosing a map center, at any given rotation angle, one can use the calibrated phase at map center to determine the angular length of each modulation cycle. For the fundamental, the Fourier coefficients can be computed from 4 angular bins, assumed to be equally spaced within one cycle. Let these bins have counts {s_{0}, s_{1}, s_{2}, s_{3}}. Then the DC, cosine and sine coeffcients, derived either by the method of least squares or Fourier analysis, may be written as:

c_{0} = [s_{0} + s_{1} + s_{2} + s_{3}]/4  (32)
a_{1} = [s_{0} − s_{2}]/2  (33)
b_{1} = [s_{1} − s_{3}]/2  (34)

Then, since σ^{2} = λ for Poisson statistics, the standard deviation for the DC coefficient is, as expected, σ_{0} = \sqrt{λ}/2. Given a value for c_{0} (the mean cycle flux), the allowed range of a_{1} and b_{1}, due to the nonnegativity of the countrates, is −k/2 to k/2, where k is the number of photons in the cycle. The standard deviations for the sine and cosine coefficients a_{1} and b_{1} are σ_{0}/\sqrt{2}.

⋯ ⋯ ⋯ ⋯

A similar analysis can be done for the higher harmonics (n > 1), where of course it is necessary to have n bins per fundamental cycle for adequate sampling.

4.2 Regularization

When the count rates are low, the time bins will contain very small numbers of photons, perhaps all ones and zeroes, particularly for the finest subcollimators, where the time bins must be small enough to resolve the modulation cycles. Even though the modulation cycles may not be apparent to the eye, use can be made of the process called “regularization” to average cycles in phase. The binning that regularizes the modulation profile is found by (1) sorting the map-center phase so that it is increasing in time, (2) assigning cycle boundaries whereever the phase is a multiple of 2π, (3) dividing each cycle into 4 bins, and (4) taking the mean count rate in each bin. This process is illustrated in Figure 4.2

After regularization, one may make in-phase averages of the modulation profile with successively smaller “super bins”. Starting with ∼ 1/4 the rotation period, one tests c_{0} − c_{2} and c_{1} − c_{3} for significance, then steps down to smaller superbins, testing at each level until one reaches the lowest timescale limited by counting statistics. Each of the differences so produced gives a Fourier representation of the amplitude and phase of the modulation profile for use in modeling (e.g. Forward Fitting) the source.

It follows logically from the above scheme that data from the finer grids contain as much information about point sources as the coarser grids, regardless of the count rate (provided only that the grid transparancies are the same). This can also be seen by the following argument. Given one subcollimator of angular pitch p_{n}, (n = 1 − 9), the modulation profile of two point sources separated by a distance Δx = p_{n} = 4.3′′ × \sqrt{(}3)^{n−1} shows slow (Δt ∼ 1 s) amplitude modulation caused by “beating” of the sources. During one half rotation, the amplitude has maxima at angles where the modulation is fastest (i.e. where the phase at map center is near zero), and minima at angles where the modulation is slowest (i.e. where the phase at map center is near maximum). By simple geometrical scaling of both the instrument and the source geometry, this statement is true for all p_{n}.

25 That is, for the finest subcollimator, two point sources 4.3′′ apart produce the same amplitude modulation as two point sources 4.3 × 3^{4} = 348′′ apart would produce for the coarsest subcollimator. Extraction of the imaging information for the finest grids is done automatically by the backprojection, FT, CLEAN and MEM algorithms, so that even if counts/timebin are of the order of unity and modulation is not apparent to the eye, the algorithms make use of the information contained in the data. For Forward Fitting, the regularization scheme described above seems to be necessary. · · · ⋯ ⋯ ⋯ ⋯ ⋯

5 Figure Captions

Figure 1.1

Schematic geometry of the HESSI subcollimators, showing representative incidence of photons with respect to the collimator axis.

Figure 2.1

Profiles of the equivalent area function for various slit/pitch ratios, and also showing energy-dependent effects that smooth the corners of the profile.

Figure ?? Diagram showing the relative positions of the spacecraft coordinates, heliocentric coordinates and the mapped region. Note that the spacecraft rotates clockwise as viewed against the Sun (i.e. α increases with time). The projected distance Θ from the line of subcollimator peak response to the map center is the phase at map center used in the calibrated event list.

Figure 3.1

An example of a back-projection image for a simulated source and the corresponding countrate profile.

Figure 3.2.1

An example simulated source, its count profile, the “dirty map”, the “CLEAN map”, and the corresponding count profile derived from the CLEAN map.

Figure 3.2.2

An example simulated source, its count profile, the “MEM map”, and the corresponding count profile derived from the MEM map.

Figure 3.2.4

An example simulated source, its count profile, the “PIXON map”, and the corresponding count profile derived from the PIXON map.

Figure 4.2 The binning that regularizes the modulation profile.

References

Cornwell, T. J., in Indirect Imaging, Proc. IAU/URSI Symp., ed. J.A. Roberts, Cambridge Univ. Press, 291, 1984.

Crannell , C.J., Hurford,G.J., Orwig,L.E., and Prince,T.A. 1986, SPIE 571, 142. Crannell,C.J., American Institute of Aeronautics and Astronautics: Washington DC, AIAA-94-0299, 1994.

Durouchoux, P., Hudson, H., Hurford, G., Hurley, K., Matteson, J., and Orsal, E., Ann. Astron. 120, (1983).

Enome, S., Adv. Space Res. 2/11, 201 1982.

Gull, S.F., and Skilling, J., 131, 646, 1984.

Högbom, J.A., Astron Astrophys 15, 417, 1974.

Holman, G. D.; Lin, R. P.; Dennis, B. R.; Crannell, C. J.; Ramaty, R.; Rosenvinge, T. T.; Canfield, R. C.; Emslie, A. G.; Hudson, H. S.; Hurford, G. J.; Madden, N. W.; van Beek, H. F.; Benz, A.; Bornmann, P. L.; Brown, J. C.; Enome, S.; Kosugi, T.; Vilmer, N.; Zehnder, A. Bull. AAS, 191, 7416, 1997.

Kosugi, T., Masuda, S., Makishima, K., Inda, M., Murakami, T., Dotani, T., Ogawara, Y., Sakao, T., Kai, K., Nakajima, H. Solar Phys. 136, 17 , 1991.

Lin, R.P., Dennis, B.R., Emslie, A.G., Ramaty, R.,Canfield, R., and Doschek, G., Adv.Space Res. 13/9,401. 1993.

Lin, R.P., Dennis, B.R., Ramaty, R., Emslie, A.G., Canfield, R., and Doschek,G., Geophysical Monograph 84, Washington DC:AGU, 283, 1994.

Lin, R. P., Hurford, G. J., Madden, N. W., Dennis, B. R., Crannell, C. J., Holman, G.D., Ramaty, R., von Rosenvinge, T. T., Zehnder, A., van Beek, H. F., Bornmann, P. L., Canfield, R. C., Emslie, A. G., Hudson, H. S., Benz, A., Brown, J.C., Enome, S., Kosugi, T., Vilmer, Nicole, Smith, David M., McTiernan, J., Hawkins, I., Slassi-Sennou, S., Csillaghy, A., Fisher, G., Johns-Krull, C., Schwartz, R., Orwig, L. E., Zarro, D., Schmahl, E., Aschwanden, M., Harvey, P., Curtis, D., Pankow, D., Clark, D., Boyle, R. F., Henneck, R., SPIE 3442, p.2, 1998.

Makishima, K., Miyamoto,S., Murakami,T., Nishimura,J., Oda,M., Ogawara,Y., and Tawara,Y., in New Instrumentation for Space Astronomy, (K.A. van der Hucht and G.Vaiana eds.), New York: Pergamon Press. 1977.

Mertz, L.N., Proc. Symp. on Modern Optics, v. 17 of the Microwave research Institute Symposia Series (Polytechnic Institute of Brooklyn), New York, 1967.

Mertz, L.N., Nakano, G.H., and Kilner, J.R., J. Opt. Soc. Am. 3, 2167, 1977.

Murphy, M.J., Nucl. Instr. and Methods in Physics Research, A290, 551-558, 1990.

Ohki, K., Tsuneta,S., Takakura,T., Nitta,N., Makishima,K., Murakami,T., Ogawara,Y., Oda,M., Miyamoto,S., Proc. Hinotori Symp. on Solar Flares, ISAS, Tokyo, Japan, p. ..., 1982

Palmer, D., and Prince, T.A., IEEE Trans. Nucl. Sci., NS-34, No. 1, 71, 1987.

Sato, Jun, PhD Thesis, NAO, 1998.

Sato, J., Kosugi, T., Makishima, K., PASJ, 51, 127, 1999.

Schmahl, E.J., AGU meeting, May, 1998.

Sivia, D.S., Data Analysis, a Bayesian Tutorial, Clarendon Press, Oxford, 1996.

Schnopper, H.W., Thompson, R.I., and Watt, S., Space Sci. Rev. 8, 534, 1968.

Schnopper, H.W., H. W.; Bradt, H. V.; Rappaport, S.; Boughan, E.; Burnett, B.; Doxsey, R.; Mayer, W.; Watt, S., ApJ 161L, 161, 1970.

Van Beek, H.F., Hoyng, P., Lafleur, B., and Simnett, G.M., Solar Phys., 65, 39, 1980.

6 Appendix

F_{0} = Incident photon flux on front grid
T = Subcollimator transmission
τ = Livetime
Φ = Phase at map center
c_{n} = Modulation Amplitude, n = 1, 2, 3, 4, 5
Φ_{n} = Modulation Phase, n = 1, 2, 3, 4, 5
p = angular pitch
L = separation between grids
p/(2L) = angular resolution
S_{triangle} = effective subcollimator area for ideal thin grids
θ = the photon′s angle of incidence relative to the imageraxis
c_{n} = modulation amplitude for harmonic n
θ_{n} = modulation phases for harmonic n
g_{n} = Fourier coefficients of the grid transmission
q = effective slit/pitch ratio
θ_{n} = relative phases for harmonic n
α = roll angle
P_{0} = modulation profile
Δt_{i} = time binsizes
Φ = phase at map center
τ = dead time
T = subcollimator transmission
E = energy
u, v = coordinates in the Fourier plane
ϕ = azimuthal coordinate in the Fourier plane
k = radial coordinate in the Fourier plane

V = visibility
C = predicted count rate
x_{coll} = subcollimator axis
x_{map} = map center
R_{0}, ϕ_{0} = location of the map center
ΔΦ = phase shift from telescope axis to sun center
P_{im} = two − dimensional modulation patterns
D_{m} = brightness of pixel m in the back − projection map
O_{i} = predicted count rate in bin i
P = HESSI Point Spread Function
μ = the loop gain in CLEAN
I_{n} = residual map intensity at nth iteration of CLEAN
H = Entropy
λ = parameter in MEM
p_{poi} = Poisson variate
S_{norm} = probability of a photon detection
a_{1}, b_{1} = cosine and sine Fourier coeffcients