COUPLED ESCAPE PROBABILITY FOR AN ASYMMETRIC SPHERICAL CASE: MODELING OPTICALLY THICK COMETS

We begin our modeling of IR rovibrational spectra of a coma by initially following the method used by Chin & Weaver (1984), Crovisier (1987), and others with some minor improvements. As in those models, we assume a constant expansion velocity, thus linearly relating any radial distance to a specific time since a "parcel" of gas was released from the surface of the nucleus. Therefore, we numerically integrate over time the linear differential equations defined by the Einstein coefficients and collisional rate coefficients to get fractional molecular energy level populations for each distance from the nucleus, from which we calculate emission spectra:

$$\begin{eqnarray} \frac{dn_{k}}{dt} &=& \sum _{i}(A_{ik} n_{i} + J_{\nu _{ik}}(B_{ik} n_{i} - B_{ki} n_{k})) + \nonumber \\ && - \sum _{l}(A_{kl} n_{k} + J_{\nu _{kl}}(B_{kl} n_{k} - B_{lk} n_{l})) \nonumber \\ && +\sum _{j} C_{jk}(n_{k} - n_{j} e^{-E_{jk}/k_{B} T})\protect. \protect \end{eqnarray} \tag{ 1 }$$

Here i, j, k, and l indicate energy levels of a molecule, with the n's with those indices being the corresponding level populations. The Einstein coefficients between levels x and y are A_xy and B_xy, and C is a similar collisional coefficient. We use $J_{\nu _{xy}}$ for the mean intensity of radiation at the frequency corresponding to the transition between x and y, E_xy is the energy difference between levels x and y, and k_B and T have their usual meanings. The summations are over all levels i, l, or j that have a transition into or out of level k. For collisional coefficients $C = n_{_{{\rm H_{2}O}}} \sigma \bar{v}$, where $n_{_{{\rm H_{2}O}}}$ is the number density of H₂O (assumed to be the dominant collisional partner), σ is the collisional cross section for a given transition, and $\bar{v}$ is the mean (thermal) molecular speed. This allows us to include a time-variable production rate. We use such a coma integration as the initial basis for our coma model before including potential optical depth effects.

Although this simple approach to expansion velocity is not strictly accurate near the nucleus (where the gas is undergoing acceleration before reaching a constant velocity), it is a commonly used approximation in conjunction with a Haser model and is often relied on by cometary scientists (e.g., to derive production rates). For the majority of the observations modeled here as hypothetical examples, it should be a sufficiently reasonable approximation, since most of the modeled observations are farther from the nucleus than the extent of major acceleration (see, e.g., Tenishev et al. 2008). In future modeling, this approximation can be replaced with a more accurate treatment for modeling Deep Impact observations of the very innermost comae of comets Tempel 1 and Hartley 2.

Our primary improvement is the inclusion of radiative transfer calculations using our spherical adaptation of the Coupled Escape Probability method (hereafter "CEP"; see Elitzur & Asensio-Ramos 2006, hereafter, "CEP06") to more correctly model optically thick (or potentially thick) regions of cometary comae. This is described in detail below and is the main part of this paper. We use the coma integration results to provide the "initial guess" values for populations used in the subsequent radiative transfer calculations using CEP.

For the purposes of the initial coma model, we treat the comet as spherically symmetric and as having a uniform and constant gas production rate over its entire surface. The outward speed of the gas is also assumed to be constant, as in Chin & Weaver (1984). While this is not strictly physically accurate (see, e.g., Combi 1996), the variation over the majority of the coma is relatively small. We use a temperature profile that varies with radial distance from the nucleus, having closely followed Combi's (1989) model (see Figure 1), and we can scale the initial gas temperature (at the surface) to alter our profile. These approximations make integration over time equivalent to calculating these values over increasing distances from the comet nucleus for a "shell" of gas expanding outward from the nucleus. (See, e.g., Bockelée-Morvan & Crovisier 1987, Combi 1996, Tenishev et al. 2008, etc., for more detailed and accurate treatments of coma temperature and velocity, parameterized over production rates, etc. For the present purpose of presenting our adaptation of CEP, the specifics of the coma integration model do not represent a limitation on our approach to radiative transfer.)

Zoom In Zoom Out Reset image size

We ignore the photodestruction of CO in our coma model, due to its long lifetime (see Crovisier 1994). The lifetime is >10⁷ s, and we are integrating out to 10⁵ km with an expansion velocity of 0.8 km s⁻¹. (Note that others, e.g., Huebner et al. 1992 and Morgenthaler et al. 2011, find a shorter lifetime, but still ≳ 4 × 10⁵ s, or ≳ 1.4 × 10⁵ s, which is still large enough that it can be neglected in our modeling out to 10⁵ km. Even for the smaller of these values, only in the outer coma, ≳ 2.5 × 10⁴ km, would the densities be reduced by ≳10%. As we are primarily interested in optical depth effects in this study, the outer coma, where optical depths will be lower, is of less interest.)

The model can include coma morphology features as well, each modeled with its own coma integration using separate conditions. Such features, as seen in the Deep Impact and EPOXI encounters, are a main motivation for creating this model to better understand possible optical depth effects in the near-nucleus regions of the coma (see, e.g., Feaga et al. 2007; A'Hearn et al. 2011). After describing our method in Section 2, we present (in Section 3) some results demonstrating its use in better understanding possible optically thick spectra for the carbon monoxide 1-0 (X¹Σ⁺) band in spherically symmetric comae. These may be useful for ground-based (or space-based) observations of comets. In Section 4, we apply our CO model to actual spectra of comet Garradd, observed at 2 AU by the Deep Impact spacecraft. Forthcoming model results for CO₂, H₂O, and near-nucleus morphological features will follow (in other papers).

In CEP06 Equation (7), Elitzur & Asensio-Ramos derive a purely analytical expression for the net radiative bracket in a plane-parallel slab, based on the formal solution of the radiative transfer equation and the definition of the net radiative bracket:

$$\begin{eqnarray} p(\tau) &=& 1 - \frac{1}{2 S(\tau)} \int _{0}^{\tau _{t}} S(t) dt \int _{-\infty }^{\infty } \Phi ^2(x) dx \nonumber\\ && \times \int _{0}^{1} e^{-|\tau -t|\Phi (x)/\mu } \frac{d\mu }{\mu }\protect, \protect \end{eqnarray} \tag{ 2 }$$

where τ and t refer to optical depths for a specific wavenumber, with τ_t being the overall optical thickness of the slab, S(τ) or S(t) to the source function for a given line ($S_{\nu _{21}} = ({A_{21} n_{2}}/{B_{12} n_{1} -B_{21} n_{2}})$), Φ(x) is a dimensionless line profile, and μ = cos θ, where θ is the angle of a given ray measured from the normal to the plane.

We use a spherical analog to this theoretical expression (i.e., as opposed to the discrete expression introduced later in CEP06 involving a number of "zones") for the net radiative bracket:

$$\begin{eqnarray} &&p(\tau (r,\theta,\phi)) = 1 - \int _{4\pi } \left [ \frac{1}{2 S(\tau (r,\theta,\phi, \Omega))}\right. \nonumber\\ &&\left.\quad\quad \times \int _{0}^{\tau (r,\theta,\phi, \Omega)} S (t(r,\theta,\phi, \Omega)\right. \nonumber \\ &&\left.\quad\quad \times \int _{-\infty }^{\infty } \Phi ^2(x) (e^{-|\tau (r,\theta,\phi, \Omega)-t(r,\theta,\phi, \Omega)|\Phi (x)} ) dx dt \right] d\Omega, \quad\quad \end{eqnarray} \tag{ 3 }$$

where p(τ(r, θ, ϕ)) refers to the net radiative bracket at any point labeled by the coordinates (r, θ, ϕ) in spherical coordinates with the origin at the center of some sphere of interest. (In our study, the sphere is centered on the comet nucleus; it could be the center of any spherical astronomical object for any given case.)

As in CEP06, $\Phi (x) = (\pi \, \Delta \nu _{D})^{-1/2}\, e^{-x^{2}}$ is a dimensionless line profile, normalized so that ∫Φ(x)dx = 1, where x = (ν − ν₀)/Δν_D is the dimensionless line shift from line center ν₀ and Δν_D is the Doppler line width, ν₀/c(2kT/m)^1/2. In the present work, we extend the CEP06 treatment of line profiles (which for the sake of presentation of a simple example used a Φ(x) that was the same throughout). We have included variations in line width (due to temperature) and Doppler shifts between different regions of the coma. The latter have been calculated for all three relevant aspects of the model with differing line-of-sight velocities: incident solar radiation, emergent flux along the observers' line of sight, and the calculation of the net radiative bracket between regions in the coma.

Both τ(r, θ, ϕ, dΩ) and t(r, θ, ϕ, dΩ) refer to optical depths and S(τ(r, θ, ϕ, dΩ)) to the source function as viewed from the coordinates (r, θ, ϕ) along a "pencil" of solid angle dΩ (see Figure 2).

Zoom In Zoom Out Reset image size

The optical depth along a "pencil" (or cone) of solid angle dΩ is, of course, highly dependent on the particular direction of a given solid-angle element dΩ, i.e., dependent on θ', ϕ', the direction of a vector centered on (r, θ, ϕ) pointing along the centerline of dΩ (see Figure 2).

The above CEP06 expression for p(τ) (Equation (2)) is turned into a discrete expression by dividing the slab into multiple "zones" (labeled with indices i and j). This yields CEP06 Equation (14):

$$\begin{eqnarray} p^{i} &=& 1 - \frac{1}{2 \tau ^{i,i-1} S^{i}} \sum _{j=1}^{z} S^{j} \times \int _{\tau ^{i-1}}^{\tau ^{i}} d\tau \int _{\tau ^{j-1}}^{\tau ^{j}} dt \nonumber \\ && \times \int _{-\infty }^{\infty } \Phi ^2(x) dx \int _{0}^{1} e^{-|\tau -t|\Phi (x)/\mu } \frac{d\mu }{\mu }\protect. \protect \end{eqnarray} \tag{ 4 }$$

This expression for p can be calculated for each zone (and each wavenumber/frequency) and the resultant p values included in the equations of statistical equilibrium, which are then solved. (See also CEP06 Equations (32)–(36), and accompanying CEP06 text, for a more complete discussion.)

Similar to CEP06 Equation (14), for the purposes of our integration of a discretized pⁱ the sphere is divided into spherical shells (analogous to the plane-parallel zones of the original CEP) where i (or j) is a shell index (see Figure 3). The integration is broken down into a sum of integrals along different "cones" of solid angle (the aforementioned dΩ's). Note that although dΩ may seem somewhat conceptually analogous to dμ in the plane-parallel case, we cannot use the convenient exponential integrals as in the plane-parallel situation, but instead must integrate along each dΩ separately:

$$\begin{eqnarray} p^{i} &=& 1 - \int _{4\pi } \left[ \frac{1}{2 \tau ^{i,i-1}(\Omega) S^{i}} \sum _{j=1}^{z} S^{j} \int _{\tau ^{i-1}(\Omega)}^{\tau ^{i}(\Omega)} d\tau \int _{\tau ^{j-1}(\Omega)}^{\tau ^{j}(\Omega)}\right. dt \nonumber \\ && \left.\times \int _{-\infty }^{\infty } \Phi ^2(x) dx (e^{-|\tau (\Omega)-t(\Omega)|\Phi (x)} ) \right] d\Omega. \qquad \end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Here we have dropped the coordinate subscripts r, θ, ϕ for clarity/simplicity and use shell indices instead (where some shell i contains the point defined by r, θ, ϕ). Note that the "z" in the summation, the maximum number of shells along a given dΩ, will be different for each dΩ. For each "cone" of solid angle we will sum S^j(Ω)e^{−|τ − t|} from "z" (the outermost shell) to i + 1, where i + 1 is the shell adjacent to shell i.

In actual practice (i.e., computer implementation), the integration of the source function of dΩ over 4π steradians will also be done by a discrete summation.

Along any particular element of dΩ originating in shell i, this sum will be

$$\begin{eqnarray} &&\sum _{j=1}^{z} S^{j}(d\Omega)\left(1 - e^{-\tau ^{j,j-1}(d\Omega)}\right) e^{-\tau ^{j-1,i}(d\Omega)} \qquad \qquad \nonumber \\ &&= \sum _{j=1}^{z} S^{j}(d\Omega)\left(1 - e^{-\tau ^{j,j-1}(d\Omega)}\right) \prod _{j^{\prime }=j-1}^{i} e^{-\tau ^{j^{\prime }}(d\Omega)}\protect. \protect \end{eqnarray} \tag{ 6 }$$

Note that unlike the plane-parallel case, each dτ must be calculated explicitly from the geometry and local molecular energy level populations for each shell along the "pencil":

$$\begin{equation} d\tau ^{j}(d\Omega) = \alpha ^{j}(d\Omega) \times ds^{j}(d\Omega), \end{equation} \tag{ 7 }$$

where α^j(dΩ) is the absorption coefficient in region j and ds^j(dΩ) is the distance through that region. Note that this distance may vary over the width of a given dΩ, but we can either approximate using the centerline of dΩ or derive a θ'-dependent expression and calculate a proper integration to get a mean value over the region. (In our implementation we use the former option, for the sake of simplicity.)

This integration over solid angle is more tedious than in the plane-parallel case (which was able to use exponential integrals over a zone's dτ) and more computationally costly, but straightforward enough to be feasible.

Turning this all into a fully discrete expression for pⁱ (for shell "i"), we get

$$\begin{eqnarray} p^{i} &=& 1 - \sum _{\omega =1}^{Z} \left[ \frac{1}{2 \tau ^{i,i-1}_{\omega } S^{i}_{\omega }} \sum _{j=1}^{z} S^{j}_{\omega } \left(1 - e^{-(\alpha ^{j}_{\omega } \times \Delta s^{j,j-1}_{\omega })}\right) \right.\nonumber \\ && \left. \times \prod _{j^{\prime }=j-1}^{i} e^{-(\alpha ^{j^{\prime }}_{\omega } \times \Delta s^{j^{\prime }, j^{\prime }-1}_{\omega })} \times \Delta \Omega _{\omega }\right]\protect, \protect \end{eqnarray} \tag{ 8 }$$

where Z is the maximum numbered spherical shell and all the quantities indicated by the subscript ω are dependent on a particular direction viewed from a given point (or shell i, in a spherically symmetric case).

Our motivating interest in this study is comets' comae, which are not spherically symmetric cases.

To adapt CEP to this asymmetry, we divide each shell further into "regions." In the case of comets, one source of asymmetry is incident solar radiation coming from one direction (outside the outermost shell). Therefore, the natural way to further divide shells into regions is along lines parallel to the direction of solar radiation (i.e., the center-of-comet-to-center-of-Sun line, which we will arbitrarily label as the z-axis). Thus, we superimpose a set of co-axial cylinders (with radii equal to corresponding shell radii, for the sake of simplicity) on the shells to divide the coma into regions bounded by two cylinders and two spherical shells. (Note that some regions, specifically those along the z = 0 plane perpendicular to the solar radiation, are only bounded by an inner cylinder and outer sphere. See Figures 4 and 5.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These regions form annuli or rings of unusual, but easily envisioned, cross sections (see Figure 5).

Incident solar radiation is parallel to the z-axis (due to our choice of the z direction). Hence, each ray of sunlight travels along the axial direction of a specific cylinder. For each region, the solar radiation absorbed is calculated based on the relevant optical depths along that direction (the dτ's) of those regions in the same cylinder that are closer to the Sun than the given region. This is similar to the case of external radiation described in CEP06 Appendix A, Equations (A3) and (A4), in which we simply set μ₀ = 1, due to the above constraint of cylinders being co-axial with the incident solar radiation:

$$\begin{eqnarray} & \bar{J_{e}^{i}} = J_{e} \frac{1}{\tau ^{i,i-1}} \left[ \gamma (\tau ^{i}) - \gamma (\tau ^{i-1}) \right], \\ &&\nonumber \\ & {\rm where} \nonumber \\ &&\nonumber \\ & \gamma (\tau) = \int _{-\infty }^{\infty } [ 1 - e^{-\tau \Phi (x)} ] dx \nonumber \protect, \protect \end{eqnarray} \tag{ 9 }$$

where $\bar{J_{e}^{i}}$ is the average over a region i of $J_{e}^{i}$, the contribution of external radiation, J_e, to the mean intensity of the region, which is to be included in the equations of statistical equilibrium (Equation (1)) by addition to the "B" term.

From a purely radiative standpoint, assuming that within each region there exists uniform density, temperature, and other physical conditions, the radiative excitation of molecules (hence the emission and absorption) in each region/annulus should be equal throughout the region. In our expanded CEP implementation, these regions are analogous to zones in the plane-parallel CEP. Each region's radiative effect or contribution to each other region, i.e., the net radiative bracket, must be calculated. Note that self-irradiation from around an annulus must also be taken into account, as well as irradiation from other regions. Once this calculation is done, the entire region will be equal with respect to radiative processes (i.e., there is symmetry around the z-axis).

If models of distantly observed comets were all we needed, this might be sufficient. However, we are motivated by the desire to better understand Deep Impact and EPOXI spectral observations that have very high spatial resolution around the comets' nuclei (see, e.g., Feaga et al. 2007, 2011). Therefore, the above radiative treatment alone is insufficiently asymmetric to fully model a cometary coma, when coma morphology is included in the model. The inclusion of morphology undoes the aforementioned symmetry around the z-axis within each annulus/region. These observations are one of the primary motives for this study, and therefore these morphological asymmetries must also be dealt with appropriately in this model.

To model morphological features, we use a cone shape superimposed over the aforementioned divisions into regions. (Other geometric shapes could also have been used. We chose to implement a cone due to its similarity in shape to many observed coma features.) A cone of arbitrary orientation and size with its vertex at the center of the sphere creates intersections with the above-described regions. Each of these is then added as a sub-region, which may have different properties from the surrounding (or subsumed/replaced) region.

Each sub-region can possess different initial conditions from the surrounding region. Thus, morphological features, which by their nature tend not to be axisymmetric around our z-axis, can be included in the model. It should be noted that these sub-regions are only included as necessary. Thus, for those annuli that do have constant axisymmetric conditions (i.e., no interesting morphological features impinging on them), we can save time and memory computationally by leaving them undivided, as they would have been originally.

Given the above geometric divisions, we have implemented Asymmetric Spherical CEP as follows. For each region (or sub-region, as applicable) we take representative population values from the coma integration and use these as the basis of an "initial guess." We then make an immediate improvement to the initial guess values by recalculating each region's populations (individually) taking into account the attenuation of incident solar radiation by intervening regions in the solar direction (as per Equation (9)). These recalculated populations are the values we then use as the initial guess (required by the implementation of Newton's method in Press et al. 1992) for CEP calculations. Based on these populations, we calculate the necessary source functions, delta-taus, and net radiative brackets, "p," as above, for each wavenumber (or line, transition, etc.). Following the above discretized Equation (8) for each region, which in this context we will call the "recipient" region, we iterate over all other regions to calculate their contribution to the recipient's p. Each other region's contribution is essentially its own source function attenuated over the optical depth of all intervening regions along the line of sight between itself and the recipient region, integrated over (or, to simplify, multiplied by) the solid angle subtended by one region from the other. This is then divided by the recipient region's source function and optical depth (along the given line of sight).

To implement this in a practical algorithm of manageable complexity, we make several simplifying approximations.

Due to the z-axis symmetry that exists (before adding sub-regions), we can partially simplify to a two-dimensional (2D) diagram in which a region is represented by the cross section of the annulus in the (arbitrarily chosen) y = 0 plane. We calculate a region's "centroid," i.e., the centroid of its 2D projected area in this plane that (in our approximation) corresponds to a point (r, θ, ϕ) in spherical coordinates.

For (cone-shaped) sub-regions, which in general do not have their centerline on the X–Z plane, we must use a different centroid. We use the midpoint along the cone's centerline (within region boundaries).

We also choose a series of points distributed evenly along circles parallel to the X–Y plane around each region, which will be the "starting points" for calculation of lines of sight (which will terminate at the centroids). These are chosen by rotating a region's centroid around the z-axis by multiples of some angle that depends on the size (radius) of the region. The choice of angle is such that the larger a region's size, the more starting points it will have, and thus the region will be divided into more elements of solid angle.

We use the line of sight between the "centroids" of regions and this series of starting points to calculate the contributions of every other region (or the region to itself) to a given recipient region's p. We calculate the optical depth of each intervening region, along the line of sight, based on the molecular population levels of the intervening regions (see Figure 6). These "integration lines" encapsulate the main part (within the square brackets in Equation (8)) of the calculation of p.

Zoom In Zoom Out Reset image size

We approximate the solid angle subtended by the integration lines from another region's centroid by integrating dΩ = dϕ sin θdθ from zero up to the mean value of the angles between the starting point of that region, the centroid of the recipient region, and the multiple "corners" of that region around the starting point (see Figure 7). Effectively, this gives a solid angle between regions of 2π(1 − cos θ_mean).

Zoom In Zoom Out Reset image size

Note that due to these approximations, the sum of solid angles over all integration lines between a region and all regions in a given shell exterior to the region is not necessarily constrained to exactly equal 4π, as it should be in reality. Therefore, in calculating a "p" value, we sum the solid angles involved and average over that sum instead of 4π steradians (as Equations (3) and (5) dictate we should do).

In the limit of arbitrarily small (and numerous) regions, these approximations would approach a physical situation of arbitrarily precise accuracy. Thus, we maintain the "exactness" of the CEP method.

Unlike the plane-parallel situation, the flux exiting the surface of the coma (or other sphere of interest) is not simply a single value (per wavenumber) that has been integrated over angle. In the spherical situation, the resultant intensities form a 2D mapping (in a plane perpendicular to the observer's line of sight; see Figure 8).

Zoom In Zoom Out Reset image size

In our implementation, this plane is specified by rotation angles, θ, ϕ, and ψ with the comet's center at the origin, and is assumed to be at a distance ⩾R_coma, the maximum radius of the coma. We can also specify the density of and interval between points on this plane for which the output intensities will be calculated. Each point in this planar mapping shows the intensity (or surface brightness) integrated along a specific line of sight, perpendicular to the plane, through the coma from one side to the other (again, for each wavenumber):

$$\begin{equation} I_{{\rm surf}} = \sum _{i=1}^{z} \left(S_{i} \Delta \nu _{i}\left(1 - e^{-\tau ^{i}}\right) \prod _{j=i-1}^{j=1} e^{-\tau ^{j}} \right)\protect, \protect \end{equation} \tag{ 10 }$$

where S_i is the source function of a region i, Δν_i is the line width of wavenumber/frequency ν in region i, and τⁱ, or τ^j, represents the optical depth of wavenumber ν in region i or j along the relevant line of sight. Indices i and j run from 1 to z, where z equals the number of regions along a given line of sight.

Thus, the spherical CEP algorithm produces results that could be described as a four-dimensional data "hypercube": for each point in the above 2D mapping, there is a complete (2D flux versus wavenumber) spectrum. These data can then be presented in multiple formats. Several forms of data presentation for simulating observations are described in the following section.

This is also precisely the output needed to compare with the Deep Impact and EPOXI observations that have been displayed as 2D brightness maps for specific wavelengths or bands.

We present here one example of a brightness map of a modeled coma (see Figure 9). This format of output presentation is most similar to the radiance maps of Feaga et al. (2007) and A'Hearn et al. (2011). For a spherical coma with no morphological features, it is rather uninteresting. It is nevertheless included here as a demonstration and used to illustrate the azimuthal angles of the radial profiles presented below with the addition of overlaid lines. Note that for the Q_CO = 10²⁸ s⁻¹ case, there is some difference in brightness noticeable to the eye between the sunward side (azimuthal angles nearer to zero) and the anti-sunward side, especially in the near-nucleus portion of the image. This is due to the optical depth along the solar direction reducing the excitation and emission from one side of the coma to the other.

Zoom In Zoom Out Reset image size

Abundances of cometary species are frequently calculated from observed fluxes using fluorescence efficiencies, or g-factors. In an optically thin case, the brightness of a given line or band is proportional to the column density of the relevant molecule. In such cases, F_band = g_band × N and g_band = ∑_bandA_u × n_u, where F_band is the band total flux (or brightness), g_band the band g-factor, N the total column density, A_u the Einstein A coefficient for the relevant transition originating in upper level "u," and n_u the column density of the population of a specific upper level "u" (which in our model is numerically approximated as the sum over all regions along a line of sight of the fractional population of level "u" times each region's column density).

However, large optical depths will spoil this simple linear relation between column density and brightness. With radiative transfer modeling, it is possible to get a calculated g-factor (g_band = ∑_bandA_u × n_u) from the model and the "effective g-factor," g_eff = F_band/N, which is the actual ratio of brightness to column density.

In Figures 10, 11, and 12 we present all these values together as radial profiles, for theoretical comets of three different production rates, Q_CO = 10²⁶, 10²⁷, and 10²⁸ s⁻¹. (All are observed at 1 AU, at a phase angle of 90° and multiple azimuthal angles. In all cases $Q_{{\rm H_2O}} = 10 \times Q_{{\rm CO}}$.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our results, we typically see that g_eff does tend toward the calculated g-factor values at larger impact parameters, where optical depth effects are negligible, as would be expected. (The actual numerical values of the "asymptotic" band g-factors produced by our model, 2.4 × 10⁻⁴ s⁻¹ per molecule for CO at 1 AU, also agree well with other published values, such as those calculated by Chin & Weaver 1984; Crovisier & Le Bourlot 1983; and Weaver & Mumma 1984.) The actual radii at which this convergence occurs depend primarily on the production rate of a comet. We can use the distance at which g_eff = 0.9 g_thin as a very rough measure of the point where a coma can be considered to transition from optically thick to thin. For the "thin" and "intermediate" coma models (Q_CO = 10²⁶ and Q_CO = 10²⁷ s⁻¹) the convergence is fairly close to the nucleus, within ∼100–200 km. However, for the "thick" model (Q_CO = 10²⁸ s⁻¹) with its high production rate, the "optically thick regime" can extend as far as O(10³) km, which can be spatially resolved even in some remote observations. Note that at radial distances very near the nucleus, worrying about optical depth effects may be relevant even for lower production rates.

Another optical depth effect is variation of brightness (and corresponding g-factor) with varying angles, both phase angle (of observer) and azimuthal angle within any given observation.

The radial profiles in Figures 10–12 demonstrate the azimuthal variation for a single phase (observing) angle. It may seem somewhat surprising and counterintuitive that the radial profile lines in the sunward direction are not the brightest (nor those closer to sunward on either side). However, when velocities and Doppler shifts are accounted for, this is readily explained. Along the sunward (or anti-sunward) direction there is no change in the sunward component of the expansion velocity, and thus no relative Doppler shifting of the line-center frequencies, all along that direction. This causes greater effective optical depths and less solar excitation in that direction. Conversely, the profiles of directions perpendicular to sunward show greater excitation. Even for the only moderately thick case of Q_CO = 10²⁷ s⁻¹ for radii ⩽100 km, there is a difference of as much as ∼5%–10% between sunward and anti-sunward directions. The effect is even more pronounced (as much as ∼20%–30%) for the thicker case of Q_CO = 10²⁸ s⁻¹.

In Figures 12–16, we present profiles of brightness (and column density) for model results observed from different phase angles, ranging from 0° to 180° , for the optically thicker case of Q_CO = 10²⁸ s⁻¹. Each plot includes multiple azimuthal angles, which would all be visible simultaneously in a wide-field observation (i.e., including the entire coma) from each given phase angle. (The above plot for 90° phase angle for Q_CO = 10²⁸ s⁻¹, Figure 12, should be considered part of this series as well.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Observations with a slit spectrometer with sufficient spatial resolution (see, e.g., DiSanti et al. 1999) might observe along one specific azimuthal angle and thus see possible variations along the slit with sufficiently high spatial resolution.

The most obvious effect seen at a glance in these figures is the spread among azimuthal angles for a given phase angle. As would be expected, the 0° and 180° phase angles (sunward and anti-sunward) have no real azimuthal variation. From phase angle 45° to 90° to 135° there is a progression: the azimuthal lines get spread out farther from each other, as well as noticeably dropping in brightness for those lines farther from the sunward side.

With respect to total brightness, the phase angles 0° , 45° , 90° , and 135° produce roughly equal peak brightness for their strongest azimuthal profiles (the more sunward directions, at the nucleus grazing radius) of about ∼3 × 10¹² photons s⁻¹ cm² sr⁻¹. Most significantly, for the 180° phase angle, the peak values are about ∼2 × 10¹² – significantly less bright than at other phase angles.

With respect to "effective" (or "observed") g-factors, the minimum values in the most optically thick regions (again, at the nucleus-grazing radii) also show a trend from sunward to anti-sunward. For the 0° through 135° phase angles, the minimum values of the ratio of the band total flux over the column density are roughly 3 × 10⁻⁵ photons s⁻¹ molecule⁻¹,while for 180° the value decreases to a minimum of about 2 × 10⁻⁵. This trend is due to a combination of two optical depth effects. The first is attenuation of incident solar light from the sunward to anti-sunward sides of the coma, leading to less fluorescent pumping, and thus less emission, by the anti-sunward direction. Second, whatever emission there is more likely to "escape" the coma along lines of sight along the azimuthal angles closer to where it is emitted, experiencing less total optical depth on its path out of the coma. Thus, the already greater emission of the sunward regions is also more likely to be observed along azimuthal directions closer to sunward.

However, the similarly located values for "actual" calculated g-factors do not follow a similar simple monotonic trend. The greatest values for a given phase angle rise from 0° through 45° and peak for phase angle 90° . From 90° through 135° down to 180° they fall through the same values, creating a symmetric peak around 90° . This symmetrical and non-monotonic pattern of the calculated values is less intuitive than the trend of g_eff over angles in the same cases. However, it is clearly understood in light of the fact that these values are based only on the actual population distributions in different regions and do not include optical depth effects on the emergent radiation. Thus, observing from phase angles 0° and 180° are sampling exactly the same lines of sight and regions' populations, including the least illuminated parts (i.e., least excited populations) of the anti-sunward side of the coma. The same is essentially true for 45° and 135° (due to symmetry around the z-axis), but they do not sample the darkest parts of the anti-sunward directions (and the differences in populations are less between outermost regions on the sunward side among azimuthal angles between ±45—they are all experiencing direct solar illumination). For 90° the higher azimuthal profiles are sampling from more excited and higher emission populations than the profiles with lower values, and consistently so all the way along their lines of sight (which is not true for 45° and 135° ). Thus, the more sunward lines for 90° are the brightest seen, and the anti-sunward values fall between values of the 0° or 180° and the 45° or 135° lines.

Finally, the profiles for Q_CO = 10²⁸ s⁻¹ have a deviation from linear slopes in their brightness in the vicinity of ∼100 ∼1000 km, most easily visible for the 180° plot, but also present for the other phase angles (and growing in size from 0° up to 180° ). This is mostly due to the temperature profile reaching its minimal values at these radii in conjunction with the higher density of the Q_CO = 10²⁸ s⁻¹ case. The higher density leads to this still being a collisionally dominated regime, and the low temperatures lead to the lowest population levels being most highly populated. These levels also have the highest Einstein A values, thus leading to a higher overall number of photons emitted for the same number of molecules. The effect is greatest for the 180° view due to a cumulative effect—the lines of sight all sample the most dense and cold regions at these radii. The same extreme effect is not seen for the 0° phase angle due to the overall greater fluorescent excitation of the sunward regions dominating it. (Note that the difference in total brightness between the two azimuthal angles at these radii is about a factor of two.)

If one is observing with high spectral resolution but low spatial resolution, the spectra observed will be the sum of as much of the coma as fills the field of view. To model this, we have simulated "aperture-averaged" spectra, where the "aperture" controls the area of the coma sampled. Our apertures are square boxes and are all centered exactly on the center of the comet, and sample a nucleus-centered area equal to the square of the "aperture size" over which we average the brightness. We present a series of apertures from 2 × 10¹ km (very near the nucleus) through 2 × 10⁵ km (the whole coma) for each of the three production rates. See Figures 17–19. (All these example spectra are modeled at a heliocentric distance of 1 AU and phase angle 90° . $Q_{{\rm H_{2}O}} = 10 \times Q_{{\rm CO}}$, as in Chin & Weaver 1984.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Band shape for apertures including the outer coma (approximately 10⁴ < R_ap < 10⁵ km) does not change significantly for different production rates. The total brightness for this regime increases approximately in linear proportion to production rate. This is due to spectra with such large aperture sizes being dominated by the fluorescence-dominated optically thin outer coma, with optical depth effects playing a minimal role (but not entirely non-existent: note the small, ≲ 3%, reduction in g-factor for the highest production rate model for R_ap = 2 × 10⁴ km).

In the inner coma, however, optical depth effects can be very striking. The "thickest" spectra (for higher production rates and smaller R_ap) have remarkably altered band shapes from the optically thin spectra.

First, however, a word about changes that are not specifically caused by optical depth effects. It is clear that there is much variation, even within a given production rate, from large to small aperture sizes. Not all of this is due to optical depth effects. Even in our optically thin case, the band shape changes noticeably in breadth. This occurs primarily due to the temperature profile. We have used a fairly simplified profile, which can be scaled to a surface temperature parameter but does not vary much otherwise between different model cases. This provides a straightforward "control" for this aspect of spectral change with aperture size.

In the innermost coma near the nucleus, the temperature is quite warm (∼100–200 K), which leads to a broader band in the 10 km spectrum. The coma gas cools to a minimum (∼20 K) around 100–200 km out from the nucleus, which produces a much narrowed band. Since the Einstein A coefficient for the lowest J lines is higher than for the lines in the "wings" of the band, the cold temperature also increases the g-factor, even in optically thick regions. At larger radii, the temperature rises again, but becomes less significant since the coma gets less dense and tends toward fluorescent equilibrium. Between these regimes, in a "transition region," there are still optical depth effects, which can be more easily isolated as g-factors are less temperature controlled.

Temperature is also a factor in determining Doppler broadening and line width, which is proportional to T^1/2, so the ratios between line widths for the coldest and the warmest regions of the coma are about 2–3, for a given wavenumber. This may lead to temperature playing a significant role in the optical thickness of the coma to incident solar radiation.

Temperature effects notwithstanding, the spectra from the denser near-nucleus regions of a coma show optical depth effects in several aspects. In addition to the total brightness no longer increasing linearly with production rate (and a corresponding reduction of g-factors), energy is also dramatically shifted between lines within the band.

The notable shifting of flux from R branch to P branch (evident in many of the thicker spectra) and shifting to lower wavenumbers in both branches (as is most evident in the R_ap = 20–200 km spectra for Q = 10²⁸ s⁻¹) are very noticeable optical depth effects. (See Sahai & Wannier 1985 for an analytical discussion of similar effects.)

This effect appears due to the branching ratio of a given pair of P- and R-branch lines originating in the same upper level ((J' + 1)/J' for the pair of lines from upper level J'), which generally (slightly) favors emission in the P-branch line. In optically thick cases, repeated absorption and emission of photons lead to a cumulative effect that favors the P branch over the R branch much more than in optically thin conditions (where it is probable that any emitted photon will not be re-absorbed before escaping the coma). See Figure 20.

Zoom In Zoom Out Reset image size

Similarly, flux is "pushed" outward in the branches, and more so in the P branch due to combination with the above effect. This is due to the lines closer to the center of the band becoming optically thick before those in the wings (both due to their higher Einstein coefficients and generally being more populated). Flux initially emitted in lines that are optically thick will through repeated absorption and emission be forced out into lines that are less thick.

As seen in Figures 17–19, the P-branch total brightness and the R-branch total brightness vary with respect to each other over different optical depths (as well as other factors, such as temperature distribution). The ratio of the sums of P/R branches' brightnesses can be useful to alert an observer (or anyone analyzing observations) that they must in a given case beware of, and if possible account for, optical depth effects.

This alone would not be sufficient, as temperatures along a given line of sight are also a significant factor in controlling the P/R ratio, in the collisionally dominated inner coma. Colder population distributions will emit more in the lower lines (in both branches) for which the ratio of P/R for each pair of lines originating in the same upper state is greater. Also, the values of τ and dτ in Equations (2)–(10) will effectively vary inversely with line width, other factors being equal.

Use of a model like ours can show where the P/R ratio is large due to temperature and where (its excess beyond that value is) due to optical depth. In our optically thin, Q_CO = 10²⁶ s⁻¹, model the P/R ratio does not exceed ∼1.4, even for aperture sizes dominated by the coldest portion of the coma. (Note, however, that this is an aperture-averaged value. In Figure 21 below, the peak value is slightly higher, ∼1.5.) However, the ratio for corresponding aperture sizes in the Q_CO = 10²⁷ s⁻¹ and Q_CO = 10²⁸ s⁻¹ cases is ∼1.6 and ∼1.9, respectively. Furthermore, in the thickest case modeled, even the spectrum with aperture size of 2000 km has a ratio of ∼1.4. Note that in all cases the 2 × 10⁵ km aperture, which is dominated by the outer coma in fluorescent equilibrium, has a ratio of ∼1.12. All of this indicates that a P/R ratio in excess of ∼1.4–1.5 is a warning sign of optical depth effects involved.

Zoom In Zoom Out Reset image size

While it would be ideal to be able to derive a simple correction factor from the P/R ratio in such cases, alas, it is not exactly possible. However, a rough estimate of the degree of optical depth effects can be derived.

To do so, we create radial profiles of the P/R ratio, for both the observed emergent flux/brightness and the calculated value based on underlying populations without attenuation of emergent light, as shown in Figures 21, 22, and 23. By cross-referencing the observed P/R ratio for a given radial distance with the corresponding g-factor in Figures 10, 11, and 12, one can ascertain the "real" g-factor to use to calculate a correct column density from the observed flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This heuristic is, of course, limited in use to the carbon monoxide X¹Σ⁺ band. Other spectra with P and R branches will have their own ratios, which can be derived by similar modeling. More complicated spectra may also, but such ratios would be more complicated to find than for cases with a simple two-branch structure.

The DI-HRI observations of Garradd had sufficient signal-to-noise ratio (S/N) to detect the coma signal above the background and extract spectra from the data in a 5 × 9 pixel aperture, corresponding to 21,050× 18,945 km² aperture, centered on the unresolved nucleus. (Note that the pixels are rectangular, not square.) Smaller apertures of 1 pixel and 3 × 5 pixels, corresponding to apertures of 4210 × 2105 km² and 12,630 × 10,525 km², respectively, were also extracted from the data. The S/N for the CO band is ∼3 in the largest aperture and as low as ∼2 in the smallest. Thus, noise in the spectra (see Figures 27–29) is significant.

The P and R branches in the initial HRI Garradd CO spectra seemed to be somewhat but not very clearly distinguishable in some of the spectra. (See Figures 27–29.) For some aperture sizes the most obvious features could possibly be noise. (e.g., The single pixel wide "peaks" on the edges of the CO band's wavelengths.) However, the peak in the P branch at around 4.65 ∼ 4.675 μm is visible although it is a very narrow one and very close to the band center.

Comparison with our theoretical models implies that the observations are dominated by cold and collisionally dominated regions of the coma.

In addition to considerations of temperature, the large P/R branch ratio of Garradd CO fluxes also indicated that some optical thickness was involved (as per Figure 19 above). The peak of the P branch is much greater than that of the R branch (in those Garradd spectra where the R-branch peak can be distinguished).

The smallest aperture Garradd spectrum appears qualitatively similar in band shape and branch ratio to our hypothetical 200 km aperture CO spectrum for Q = 10²⁸ s⁻¹ (see Figure 19(c)) and the largest aperture spectrum somewhere between the theoretical 200 km and 2000 km (see also Figure 19(d)) spectra. The 200 km aperture spectrum of the theoretical model samples both the coldest (∼20 K) part of the coma and a fairly optically thick part, with a P/R ratio ≳1.8 (some of which is due to low temperature, not only optical thickness). Thus, in conjunction with our low temperature, a factor of a few times greater density than that theoretical model (and corresponding increase in optical thickness) seemed a reasonable starting point to match the observations. This would be somewhat greater than the estimate based on observations and the use of the "optically thin" g-factor (see Table 3).

Based on preliminary production rates derived from those observations (see Table 3), we modeled all three species for Garradd, at a distance of 2 AU. The production rates and other model parameters used are listed in Tables 4 and 5.

Table 4. Model Input Parameters for Modeling of Actual Observations of Comet Garradd for CO

Model Parameters Used for C/2009 P1 Garradd
Mean nucleus radius	3 km
Heliocentric distance	2 AU
Heliocentric speed	14.5 km s−1
Tsurface	40 K (tested 20–100 K)
Expansion speed Vexp	0.5 km s−1
QCO	3.2 × 1028 s−1
$Q_{{\rm H_{2}O}}$	4.6 × 1028 s−1
Band center wavenumber	2149 cm−1
Band center Einstein A	33 s−1
Highest rot. level, Jmax	20
Solar flux (at 2 AU)	6.25 × 1012 photons cm−2 s−1 cm−1
σrot	3.2 × 10−15 cm2

Download table as:  ASCIITypeset image

We used a simple spherical model without including any coma morphology. We used a convolution and pixel binning simulating that of the Deep Impact HRI instrument, with λ/δλ ∼ 200–330 (see Hampton et al. 2005).

We tested a range of CO production rates beginning with the large-aperture optically thin value of Q_CO ∼ 2.7 × 10²⁸ s⁻¹ and up to 4 × 10²⁸ s⁻¹. In conjunction with the best-fit temperature, Q_CO = 3.2 × 10²⁸ yielded the closest total brightness values to the data.

Fitting the temperature proved crucial to fitting the data with the model. Villanueva et al. (2012) measured a rotational temperature around 40 ± 7 K (1σ) for Garradd at 2 AU for several species (within ∼3000 km of the center of the comet, judging from the spatial profiles in their Figure 2; note, however, that their measurement was pre-perihelion and using other species). We tested temperatures within the range of the 40 ± 7 K measurements of Villanueva et al. (2012), and even some below the lower end of their 1σ range of 33 K. We found that these changes in temperature made some relatively small differences in the resulting spectra.

The best fit was achieved using a constant temperature of 40 K throughout the coma, in place of our initial profile (Figure 1). This is a gross but not unreasonable approximation. It is in rough agreement with Villanueva et al. (2012). Furthermore, our initial temperature profile has a fairly broad trough around its minimum of ∼20 K that extends out to several hundred kilometers from the nucleus. It is reasonable that the warming of the coma after the minimum should be considerably slower/farther out in the coma at 2 AU than at 1 AU (since this warming is primarily due to photodissociation of H₂O) and probably stretches out the low-temperature zone discussed to a few thousand kilometers, or perhaps even farther. It should also be noted that the outer coma, where fluorescence dominates, is mostly unaffected by this approximation. Thus, our use of a constant temperature of 40 K makes the inner coma cold so as to fit the band profile in the inner coma, but has no effect on spectra of the outer coma where fluorescence dominates.

Due to noise considerations, the aperture-averaged band total brightness was deemed the most important measure of quality of model fitting, followed by peak values and band shape. See Table 3 for these brightness values and calculated column densities and production rates derived from them (see Feaga et al. 2014 for further details). Corresponding model results are shown in Table 5.

We present our model results in Figures 24, 25, and 26, as well as the model overplotted on the actual data in Figures 27, 28, and 29. Also see the values shown in Table 5 (and compare with observational data in Table 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our "best-fit model" for CO used a production rate of Q_CO = 3.2× 10²⁸ s⁻¹. This value is greater than the value that would be calculated (from the largest aperture) using the optically thin g-factor of g_CO = 2.4 × 10⁻⁴ s⁻¹ (1 AU value), which would yield Q_CO ∼ 2.7 × 10²⁸ s⁻¹ (for which a model was tried, but it produced a less close match to the data). Thus, some optical depth effects are still indicated for CO, even in the largest aperture size.

Figure 30 shows radial profiles for our "best-fit" model. The locations of the vertical lines indicating where the effective g-factor is 90% and 98% of the optically thin g-factor also demonstrate that none of the spectra observed out to a radius of approximately 10,000 km are entirely free from optical depth effects. Furthermore, although these effects are relatively small at 10,000 km, they are considerably greater for the smaller aperture spectra.

Zoom In Zoom Out Reset image size

4.3.1. Garradd Model Conclusions

The DI HRI Garradd spectra show optical depth effects, even for the largest aperture, which includes approximately a 10,000 km radius from the center of the comet. Our model results show that even a box aperture integrated out to that radius would underestimate the production rate by about 15% if an optically thin g-factor were assumed.

Numerically, these effects are also shown in Table 5, in the section showing percentage of optically thin g-factors; the values for the smallest aperture are about 50%–60% of the optically thin values. Thus, were one to naively use the optically thin g-factor to calculate production rates based on brightness values from the coma up to about 1000–2000 km radius, the resulting production rates would be incorrect by almost a factor of two!

In Figures 31, 32, and 33, we plot the residuals of the model and data, for the CO band. As can be readily seen, the level of noise on either side to the left or right of the band is of a similar magnitude to the residuals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our models have fit the data quite well in band total brightness, and fairly closely in peak values and band shape. The differences between the models and the data are at levels comparable to the noise and uncertainty of the observations themselves, particularly with respect to band shapes. The production rates used in these models are somewhat higher, but not dramatically so, than those that would be derived from the data based on the largest aperture assuming that the observations were optically thin. However, the production rate is significantly higher than what would be derived from the smallest aperture data assuming optical thinness.

Temperature also proved important in modeling band shape, to the noise-limited extent that shape fitting was possible. Our best-fit models use a fairly low temperature of 40 K throughout the inner coma, which at first may seem quite low at 2 AU, but it is in agreement with other observations (Villanueva et al. 2012.) This low temperature, in the collisionally dominated inner coma, works to fit the relatively narrow band shapes observed.