OBSERVATIONAL STUDIES ON THE NEAR-INFRARED UNIDENTIFIED EMISSION BANDS IN GALACTIC H ii REGIONS

The present studies are based mainly on the slit spectroscopic data obtained by AKARI/IRC in the framework of a director's time program during the AKARI posthelium mission phase (Phase 3), during which only the NIR channel (2.0–5.5 μm) was in operation. All of the observations were performed by the use of the "Ns" or "Nh" slit (08 length × 5'' width or 1' length × 3'' width, respectively) for diffuse sources (see Onaka et al. 2007 for more detailed information on the instrument). The director's time program was planned for the purpose of calibration. In the present observations, for the relative calibration of the NIR, disperser grism and prism spectroscopy were done separately for the same slit area in the first and the latter half of a single pointing observation. This procedure provides two kinds of spectra for the same slit area with different wavelength coverages and spectral resolutions (2.5–5.0 μm, R ∼ 100 for the grism and 1.7–5.4 μm, R ∼ 20–40 for the prism; see Figure 1). The target sources were selected from Galactic ultracompact (UCH ii) and giant H ii (GH ii) regions cataloged by Crowther & Conti (2003) and Conti & Crowther (2004) according to the visibility of the satellite, except for two unclassified infrared sources, G331.386-0.359 and G345.528-0.051. These sources are unclassified, but their spectra show H ii region-like properties, and we classify them as H ii region-like objects. Generally, two pointing observations were performed for each target, whereas some targets were observed only once. We selected 61 pointing observations of 36 different targets in total that are not severely affected by instrumental artifacts or stray light. The observation logs and target parameters are summarized in Table 1. We narrowed the selection further at the spectral analysis stage by removing those without a sufficient signal-to-noise ratio to quantify the PAH emission features and extinction and, in some cases, those without available complementary MIR imaging data (see details below) from the following analysis. The accurate slit position is determined from the 3.2 μm image (N3) taken during each pointing observation by referring to the positions of point sources in the 2MASS catalog.

Zoom In Zoom Out Reset image size

The data reduction basically follows the official IDL pipeline for AKARI Phase 3 data (ver. 20111121; Ohyama et al. 2007; Onaka et al. 2009). Since the wavelength is fixed for a given detector column in the IRC slit spectroscopy, a systematic pattern from the detector anomaly or a dark pattern in the column direction sometimes persists even after the standard pipeline process. The grism spectroscopic data are especially sensitive to the fluctuation of those patterns. Therefore we perform additional dark-subtraction for the grism spectroscopic data. In the grism mode, the brink of the image can be used to correct for some systematic error. It is blocked by the aperture mask and allows us to estimate the residual dark currents. We obtain median values of 10 pixels corresponding to a given wavelength in this region and subtract them from the spectrum.

The spectral data generally show a spatial variation within the entire slit area. The FWHM of the point-spread function of the IRC is approximately 3.2 pixels (∼46) at the NIR wavelengths during Phase 3 period (Onaka et al. 2010). We extract spectra taken with the grism and the prism mode, respectively, from the same slit area with the length of 6 pixels (∼9'') one by one along the slit direction, avoiding those affected by point-like sources. As a result, we obtain several sets of grism and prism spectra from each pointing observation (see Appendix B). Even though the telescope is pointed toward the same position on the sky due to the different optical alignments, spectral images are shifted systematically by about four pixels along the slit direction between the grism and the prism modes. When we create a pair of grism and prism spectra of the same position, this position shift is taken into consideration. Figure 1 gives an example of the obtained spectra. In the short wavelength region, the spectral resolution of the prism is too low to discuss features quantitatively. Furthermore, saturation frequently occurs in the prism spectra because of its extremely high sensitivity, particularly in this wavelength range. For these reasons, we truncate the prism spectra and only use a 4.35–5.4 μm wavelength region in the following analysis. Appendix B exhibits all of the resultant spectra together with the individual S9W band images taken by the AKARI MIR all-sky survey (see details below), in which the location of the slit and each spectra-extracted region are depicted with colored boxes.

We secondarily utilize the diffuse data obtained in the AKARI MIR all-sky survey (Ishihara et al. 2010). The AKARI MIR all-sky survey was carried out by the IRC at two MIR wide band filters, S9W (λ_eff = 9 μm, 6.7–11.6 μm) and L18W (λ_eff = 18 μm, 13.9–25.6 μm). In the NIR regime, the effect of the zodiacal light is negligible, but its intensity increases dramatically by a factor of approximately 100 in units of Jy at MIR wavelengths. We therefore subtract the medium value of a relatively dark area within the same frame from the processed images as the background, most of which is expected to arise from the zodiacal emission, and perform aperture-photometry within a 5'' radius circle centered on each spectrum-extracted region. From the nature of the present targets, the background signals that we subtract are considerably small relative to the source ones in most cases. The pixels where saturation occurs due to the presence of bright astronomical sources are automatically masked through the data reduction process. Some of the present target objects (e.g., M17a) contain such masked regions, which are excluded from the analysis involving the MIR imaging data.

First, we discuss the analysis of spectra taken with the grism. The ice absorption features can be regarded as a "foreground-screen." Except for the broad absorption feature around 3.05 μm associated with H₂O ice, the CO and CO₂ ice absorption features cannot be completely resolved with the spectral resolution of the AKARI/IRC slit spectroscopy. We fit the laboratory spectrum of the pure H₂O ice at 10 K to the observed spectra, while the CO and CO₂ ice absorption features are modeled with Gaussian functions with negative signs. Then, the function we fit to the observed spectra is given by

$${\begin{eqnarray} \nonumber F_{{ \lambda }}{ (\lambda)}&=&\left [ \sum _{ k=0}^{ 5} a_{ k} \lambda ^{ k} + \sum _{ k_l=1}^{ 2} b_{ k_l} f_{ {\rm l}}{(\lambda ;\,\lambda _{ k_l};\,\gamma _{ k_l})} + \sum _{ k_g=1}^{ 12}c_{ k_g} f_{ {\rm g}}{ (\lambda ;\,\lambda _{ k_g};\,\gamma _{ k_g})} \right ] \\ &&\cdot \exp \left(-d\,{\rm ln10}\,A_{ \rm H_2O}{ ({\rm \lambda })}\right) - \sum _{ k^{\prime }_g=1}^{ 2}e_{ k^{\prime }_g} f_{ {\rm g}}{(\lambda ;\,\lambda _{ k^{\prime }_g};\,\gamma _{ k^{\prime }_g})},\qquad \end{eqnarray}} \tag{ 1 }$$

where f_l(λ; λ₀; γ) and f_g(λ; λ₀; γ) denote Lorentzian and Gaussian functions with the central wavelength λ₀ and the FWHM γ, respectively. $A_{\rm H_2O}({\rm \lambda })$ denotes the absorbance of the H₂O ice. Here we employ the laboratory data, which are retrieved from the Leiden atomic and molecular database Ehrenfreund et al. (1996). The second term represents the CO and CO₂ absorption features. The FWHM $\gamma _{k^{\prime }_g}$ and the center of the wavelength $\lambda _{k^{\prime }_g}$ are fixed to the best-fit values, which are estimated from the spectra with good signal-to-noise ratios. Due to the different slit width (5'' for the "Ns" slit and 3'' for the "Nh" slit), the spectral resolution is slightly different between the spectra taken with "Ns" slit and those with "Nh" slit, and thus the best-fit values are derived for the "Ns" slit and "Nh" slit, respectively (see Table 2). The terms in the first brackets of Equation (1) represent emission features, which are made up of three components, a continuum, PAH, and line emission features. The continuum is modeled with a polynomial function of the 5th order and constrained to be nonnegative. We express both the PAH 3.3 μm and the adjacent 3.4 μm subfeature with a Lorentzian function, whereas the subfeatures at around the 3.5 μm are altogether modeled by one Gaussian function centered at 3.48 μm, even though they are supposed to consist of multiple components. It is because the spectral resolution achieved by AKARI/IRC is not sufficient to resolve each of faint features around the 3.5 μm. We attempt other combinations of Gaussians and Lorentzians, but the adopted combination provides the best fit. The center of the wavelength and the FWHM of these three components are fixed to the best-fit values estimated for the spectra obtained with AKARI/IRC (see Table 3). The emission lines are modeled with Gaussian profiles with the FWHM $\gamma _{k_g}$ fixed to match with the spectral resolution, 0.031 μm for "Ns" and 0.025 μm for "Nh" (see Table 4). In the fitting, only a_k, $b_{k_l}$, $c_{k_g}$, d, and $e_{k^{\prime }_g}$ are free parameters. The fitting is nonlinear and thus we use the Levenberg–Marquardt method (Press et al. 2002). We note that the wavelength region of 2.9–3.1 μm is not used in the fitting, since this region is at the bottom of the H₂O absorption feature and is easily contaminated by weak foreground emissions. The left panel of Figure 3 shows an example of the fitting for the grism spectra.

Zoom In Zoom Out Reset image size

Table 2. Gaussian Profile Parameters for Ice Features

Band	$\lambda _{k^{\prime }_g}$ (μm)	$\gamma _{k^{\prime }_g}$ (μm) for "Ns" Slit	$\gamma _{k^{\prime }_g}$ (μm) for "Nh" Slit
CO2 Ice	4.26	0.060	0.048
CO Ice	4.67	0.031	0.025

Download table as:  ASCIITypeset image

Table 3. Lorentzian and Gaussian Profile Parameters for UIR Bands

$\lambda _{k_l}$ (μm)	$\gamma _{k_l}$ (μm) for "Ns" Slit	$\gamma _{k_l}$ (μm) for "Nh" Slit
3.29	0.048	0.045
3.41	0.044	0.043
$\lambda _{k_g}$ (μm)	$\gamma _{k_g}$ (μm) for "Ns" Slit	$\gamma _{k_g}$ (μm) for "Nh" Slit
3.48	0.113	0.100

Download table as:  ASCIITypeset image

Table 4. Gaussian Profile Parameters for Emission Lines

Line	$\lambda _{k_g}$ (μm)	$\gamma _{k_g}$ (μm) for "Ns" Slit	$\gamma _{k_g}$ (μm) for "Nh" Slit
H i Hu15	3.91	0.031	0.025
H i Hu13	4.18	0.031	0.025
H i Hu12	4.38	0.031	0.025
H i Pf12	2.76	0.031	0.025
H i Pfη	2.88	0.031	0.025
H i Pfγ	3.75	0.031	0.025
H i Pfβ	4.65	0.031	0.025
H i Brβ	2.63	0.031	0.025
H i Brα	4.05	0.031	0.025
He i (3D1–3F0)	4.30	0.031	0.025
H2 0–0 S(13)	3.85	0.031	0.025

Download table as:  ASCIITypeset image

In the prism spectroscopy, the spectral dispersion varies with wavelength. Taking into account this non-uniform spectral dispersion in the prism mode, we succeed in reproducing the prism spectrum from the grism spectrum by a simulation (see details in Appendix A). On the basis of this analytical result, we simulate the spectral profile in a spectrum taken with the prism mode, $f_{k_f}$(λ), for each fitting component. As the fitting components, three emission features are selected: Pfβ at 4.65 μm, Huδ at 5.13 μm, and the 5.25 μm band, which are clearly seen in most of the resultant spectra taken with the prism. Then, we fit the spectra with the function

$$\begin{eqnarray} &&F_{\scriptscriptstyle \lambda }{\scriptstyle (\lambda)} = \sum _{\scriptscriptstyle k=0}^{\scriptscriptstyle 5}a_{\scriptscriptstyle k} \lambda ^{\scriptscriptstyle k} + \sum _{\scriptscriptstyle k_f=1}^{\scriptscriptstyle 3}b_{\scriptscriptstyle k_f} \cdot f_{\scriptscriptstyle k_f}{\scriptstyle (\lambda)}, \end{eqnarray} \tag{ 2 }$$

where the first term represents the continuum, and the $f_{k_f}$(λ) in the second term is the modeled function of each fitting component. In the formulation of the $f_{k_f}$(λ), we assume that Pfβ and Huδ lines have a Gaussian profile and the 5.25 μm band has a Lorentzian profile in grism spectra. The FWHM of the 5.25 μm band is set to the literature value reported by Boersma et al. (2009), 0.12 μm, and the central wavelength is slightly modified to match with the best-fit value for the AKARI/IRC prism spectroscopy, 5.22 μm. An example of the fits is given by the right panel of Figure 3. As described above, the wavelength calibration in the prism mode is no better than ∼0.03 μm in this spectral range, but the assumed spectral profile fits the observed spectrum fairly well.

We adopt the dust model of "Milky Way, R_V = 3.1" of Weingartner & Draine (2001) as the extinction curve and estimate the value of the visual extinction A_V from the observed intensity ratio of Brβ to Brα with the assumption of the Case B conditions of T_e = 10⁴ K and n_e = 10⁴ cm⁻³ (Storey & Hummer 1995) and the foreground-screen extinction (see Figure 2). It is remarked that the line and band intensities derived from the fitting are corrected for the extinction. The A_V value ranges widely from 0 to 40 mag. In several regions, the effect of the extinction cannot be ignored. In addition to the extinction correction, we also correct for a contribution from the unresolved emission line Pfδ at 3.30 μm to the UIR 3.3 μm band. According to the Case B condition, we assume that the intensity of Pfδ is equal to 9.3% of the extinction-corrected intensity of Brα and subtract it from the extinction-corrected intensity of the 3.3 μm band.

Based on the results of their laboratory experiments and quantum-chemical calculations, previous studies indicate that the astronomical 5.25 μm band could be a blend of overtone, difference, and combination bands of the fundamental C–H stretching and bending vibration modes (Allamandola et al. 1989a; Boersma et al. 2009). Boersma et al. (2009) also suggest that C–C modes of large ionized PAHs contribute to those two bands. From the point of view of observations, Boersma et al. (2009) report that the 5.25 μm band strength is strongly correlated with the 11.3 μm band strength, which supports their close connection with C–H vibration modes. However, their observational study is based on only four spectra obtained by the ISO/SWS for HD 44179 (post-AGB star), NGC 7027 (planetary nebula), and two positions in the Orion bar (H ii regions). A larger collection of astronomical spectra is important to make a firm conclusion.

Here, in order to consider the nature of the 5.25 μm band, we investigate their relative variation to the 3.3 μm band and the 3.4–3.6 μm subfeatures (the summation of the 3.41 and 3.48 μm components). Figures 4(a) and (b) show the relative strength of the 5.25 μm band with respect to the 9 μm surface brightness against those of the 3.3 μm band and the 3.4–3.6 μm subfeatures, respectively. There is a tight correlation of the 5.25 μm band both with the 3.3 μm band and the 3.4–3.6 μm subfeatures (the weighted correlation coefficient r is 0.92 and 0.88, respectively). As mentioned above the 3.3 μm band is assigned to an aromatic C-H stretching mode, and the 3.4–3.6 μm subfeatures are assigned to aliphatic C–H vibrational modes (Duley & Williams 1981), or overtone of the fundamental aromatic C–H stretching mode. The present result confirms that the 5.25 μm band has a close relation to C–H vibration modes. Given that the emission feature at a short wavelength like the 3.3 μm band requires high excitation, it is compatible with the hypothesis that the overtone bands of fundamental C-H vibration modes contribute to the 5.25 μm band. A slightly higher correlation is found with the 3.3 μm band than with the 3.4–3.6 μm subfeatures, suggesting that the 5.25 μm band has stronger connection to aromatic ones. Since the 3.4–3.6 μm subfeatures correlate with the 3.3 μm band to some extent, the correlation between the 3.4–3.6 μm subfeatures and the 5.25 μm band may be a secondary relation. However, the difference between those two correlations is not significant, and thus we cannot draw a clear conclusion on this point from the present data.

Zoom In Zoom Out Reset image size

We fit a linear function, y = A_intercept + B_slope · x, where y is the relative strength of the 5.25 μm band to the 9 μm surface brightness and x is that of the 3.3 μm band for all the data points and for those classified into UCH ii and GH ii regions separately, using the IDL routine "FITEXY" (Press et al. 2002). Th fitting results are summarized in Table 5. In all cases, the intercept A_intercept is very small, and the regression lines pass through the origin within uncertainty. The slope B_slope does not differ significantly between the UCH ii and H ii regions and is in agreement with the intensity ratio of the 5.25 μm band to the 3.3 μm band of HD 44179, NGC 7027, and two positions in the Orion bar reported by Boersma et al. (2009), which ranges from 0.16 to 0.38 within a factor of two. These results suggest that the relative intensity of the 5.25 μm band to the 3.3 μm band is rather constant over a wide range of Galactic objects regardless of their evolutionary stage.

It is widely known that the 3.3 μm PAH band is commonly accompanied by adjacent satellite features at 3.41, 3.51, and 3.56 μm, as well as a broad plateau extending to 3.6 μm for various kinds of astronomical objects (e.g., Geballe et al. 1985; Joblin et al. 1996). However, the origin of these satellite features has not been clearly identified. Aliphatic hydrocarbon materials have a fundamental C–H stretching mode at 3.4 μm (Pendleton & Allamandola 2002). Hydrogenated PAHs also cause the 3.4 μm band. Additional H atoms to peripheral C atoms of PAHs convert aromatic rings to aliphatic rings, which newly create aliphatic C–H stretching bands, especially near 3.4 μm and 3.5 μm (Bernstein et al. 1996). Aliphatic side-groups at the periphery of PAHs like methyl (–CH₃), methylene (–CH₂–), and ethyl (–CH₂CH₃) also give rise to the satellite features around 3.4–3.5 μm as aliphatic C–H stretching modes (e.g., Duley & Williams 1981). On the other hand, overtone bands of the aromatic C–H stretching mode also appear in that spectral range due to the anharmonicity (at 3.41, 3.47, 3.5, and 3.56 μm, Barker et al. 1987). Joblin et al. (1996) report the 3 μm spectrum variation within the reflection nebulae NGC 1333 SVS3 and NGC 2023, suggesting that the observed 3.4 μm band is too intense to originate from overtone bands of the aromatic C–H stretching mode alone, and that the relative variation of the 3.4 μm band with respect to the 3.3 μm band can be explained by photochemical erosion of aliphatic side-groups attached to PAHs due to the strong FUV radiation field.

In spectra of reflection nebulae, Sellgren (1984) indicate the presence of the continuum emission with a high color temperature of ∼1000 K in the NIR range. From observations of the Orion Nebulae with the Wide Field Cryogenic Telescope-II, Haraguchi et al. (2012) suggest that after subtraction of the free–free emission from the 3.7 μm continuum there remains a significant residual component, and that its relative strength to the 3.3 μm band emission reflects the ionization degree of PAHs. In the present study, the relative intensity of the 3.4–3.6 μm subfeatures to the 3.3 μm band, I_{3.4–3.6 μm}/I_3.3 μm, lies in a range of relatively small values, 0.0–0.5, which implies that the aliphatic structure is a small part of the band carriers, at least for the present target sources (e.g., Li & Draine 2012), but exhibits a small, but clear, systematic variation against the intensity ratio of the 3.7 μm continuum intensity to the 3.3 μm band, I_{cont, 3.7 μm}/I_3.3 μm (see Figure 5). The 3.7 μm continuum intensity, I_{cont, 3.7 μm}, is the integrated flux over a 3.65–3.71 μm range. Hydrogen free-free emission and free-bound and helium free-free emission are estimated from the intensity of Brα and are subtracted from the observed continuum. As suggested by Haraguchi et al. (2012), we assume that the ratio of I_{cont, 3.7 μm}/I_3.3 μm mirrors the ionization degree of PAHs. As the ionization degree of PAHs increases, the ratio of I_{3.4–3.6 μm}/I_3.3 μm slightly decreases with a large scatter.

Zoom In Zoom Out Reset image size

The aliphatic fraction of the band emitters is thought to be a major factor that influences the observed intensity ratio of I_{3.4–3.6 μm}/I_3.3 μm. Aliphatic C-H bonds are less resilient than aromatic ones, and they need less energy to break. The laboratory experiments reveal that irradiation of UV light and heating drive hydrocarbon solid materials to evolve to aromatic-rich materials (e.g., Iida et al. 1985; Smith 1984; Sakata et al. 1990). Recently, Jones (2012a, 2012b, 2012c) presented a theoretical model for the evolution of amorphous hydrocarbon materials and suggested that UV-photolysis leads to the ultimate transition toward aromatic-dominated materials, especially for small grains (a ≲ 20 nm). Given that the grains in this size range dominate the 3.3–3.6 μm band emission (e.g., Schutte et al. 1993; Draine & Li 2007), his prediction can be applicable to our study. The PAH ionization degree is controlled by a balance between photoionization and recombination with ambient electrons. Among the present target regions, the variation of the UV radiation field is expected to be the most contributing factor in the difference of the PAH ionization degree. The observed ratio of I_{3.4–3.6 μm}/I_3.3 μm decreases with the increase of the PAH ionization degree. This trend is consistent with the scenario that aliphatic C–H bonds in the band emitters are efficiently destroyed prior to aromatic C–H bonds in ionized gas and/or its boundary with the evolution of the UV radiation field. If a photothermal process is the main cause of the decrease of the 3.4–3.6 μm subfeatures relative to the 3.3 μm band, we can investigate the evolutionary history of hydrocarbon materials in interstellar environments from the intensity ratio of I_{3.4–3.6 μm}/I_3.3 μm. Kaneda et al. (2012) and Yamagishi et al. (2012) report the enhancement of the 3.4–3.6 μm subfeatures relative to the 3.3 μm band in the CO molecular loop near the Galactic center and in the halo of M82. They conclude that this enhancement arises from shattering of carbonaceous grains driven by shock winds (Jones et al. 1996). In contrast, the present study indicates no such enhancement in the ionized gas-dominated regions. Therefore, it can be inferred that shattering is not an effective process in the interior and the periphery of H ii regions, which is consistent with the model prediction of Micelotta et al. (2010). Some previous studies suggest that the variation of the 3.4–3.6 μm subfeatures relative to the 3.3 μm band reflects the processing of interstellar hydrocarbon materials due to thermal annealing, but they are discussed only for a limited number of astronomical objects including H ii regions, planetary nebulae, and reflection nebulae (e.g., Joblin et al. 1996; Goto et al. 2003; Boulanger et al. 2011). The present study is the first investigation based on a large quantity of high-quality spectra of various Galactic H ii regions, and it gives the first clear evidence for dust processing in an ionized medium.

The AKARI S9W band filter covers a wavelength range of 6.7–11.6 μm, which includes the prominent PAH emission band features at 6.2, 7.7, 8.6, and 11.3 μm. Therefore, it can be conjectured that the 9 μm band intensity mostly originates in these PAH emission features (e.g., Ishihara et al. 2007; Kaneda et al. 2012). On the other hand, the origin of the 18 μm emission is still ambiguous. In the report of AKARI imaging observations of the reflection nebulae IC 4954 and IC 4955 region, Ishihara et al. (2007) suggest that the continuum emission originating from stochastic heating of very small grains (VSGs) dominates spectra in the wavelength range of the L18W band filter. However, this assumption does not always hold in the case of H ii regions, the target objects of this study. In H ii regions, the intensity of radiation field reaches 10⁴–10⁵ times that in the solar vicinity (e.g., Tielens et al. 1993; Berné et al. 2009; Salgado et al. 2012). In such a harsh environment, the equilibrium temperature of big grains (BGs) becomes higher, and thermal emission can be dominant even at short wavelengths around 18 μm.

As shown in Figure 6, the MIR color of AKARI 9 μm to 18 μm band steeply declines against the ratio of I_Brα/I_3.3 μm. The ratio of I_Brα/I_3.3 μm is considered as a good indicator of the fraction of the ionized gas along the line of sight. The present result suggests that the MIR excess around 18 μm becomes stronger relative to the PAH 9 μm emission with the transition from PDRs to ionized-gas dominated regions. The observed trend can be explained by PAH destruction and either VSG or BG replenishment inside the ionized medium, which are proposed by several authors as a possible interpretation of Herschel and Spitzer multiwavelenght observations (e.g., Paradis et al. 2011; Paladini et al. 2012). On the other hand, it is expected that the temperature of BGs is raised with the increase of the incident radiation field from PDRs to ionized-gas, which can also contribute to the observed trend. The present result may reflect the variation in the incident radiation field between PDRs and ionized gas rather than that of dust abundance. These two effects cannot be distinguished from the present dataset, and we cannot draw a clear conclusion at the moment. In comparison with the GH ii regions, the UCH ii regions are distributed over a different region in Figure 6 in the relatively left side of the diagram. This might result from the different geometric structures between them, that is, larger contamination from their surroundings in the UCH ii sources. The AKARI MIR survey covers 99% of the whole sky. The wide scope of the AKARI MIR survey will be very helpful to investigate the dust processing in the ISM.

Zoom In Zoom Out Reset image size

We summarize the NIR spectra and the S9W image of each target object in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (36 panels)

Tarball of figure set images