Solar fluorescence model of CH3D as applied to comet emission

Journal of Molecular Spectroscopy 291 (2013) 118–124

Contents lists available at SciVerse ScienceDirect

Journal of Molecular Spectroscopy journal homepage: www.elsevier.com/locate/jms

Solar fluorescence model of CH3D as applied to comet emission Erika L. Gibb a,⇑, Boncho P. Bonev b,c, Geronimo Villanueva b,c, Michael A. DiSanti c, Michael J. Mumma c a

Department of Physics & Astronomy, University of Missouri – St. Louis, 503 Benton Hall, One University Blvd., St. Louis, MO 63121, USA Department of Physics, The Catholic University of America, Washington, DC 20064, USA c Goddard Center for Astrobiology, Solar System Exploration Division, NASA Goddard Space Flight Center, Code 690, Greenbelt, MD 20771, USA b

a r t i c l e

i n f o

Article history: Available online 25 May 2013 Keywords: Methane Fluorescence Comets Model

a b s t r a c t We developed a solar fluorescence emission model for the m4 band of CH3D for application to low-density, optically thin environments like cometary atmospheres. Our model utilizes transition frequencies, energy levels, and line strengths from the HITRAN 2008 database. We calculated the statistical weights of states with A-type and E-type spin symmetry, the partition functions for each spin manifold as a function of rotational temperature, and the Einstein A- and B-coefficients for individual rotational–vibrational transitions. We used these parameters to construct a database of fluorescence efficiencies (g-factors) at rotational temperatures relevant to cometary atmospheres. The effects of nuclear spin temperature and rotational temperature on the resulting g-factors are disentangled and modeled separately. The primary application of this model is to constrain the abundance ratio CH3D/CH4 in cometary ice, a possible cosmogonic parameter. We present the upper limit for CH3D/CH4 in comet C/2007 N3 (Lulin) using our model and illustrate the ability of current spectrometers to stringently constrain astrochemical model predictions for CH3D/CH4 in the early solar system. Ó 2013 Elsevier Inc. All rights reserved.

1. Motivation Methane and its isotopologues (CH3D, 13CH4) play a fundamental role in the organic chemistry of planetary atmospheres, interstellar material, and the origins of the solar system. The methane abundance in cometary bodies can serve as a revealing cosmogonic indicator. For interstellar clouds, current models for formation of molecules at very low temperatures show that equilibrium abundance ratios for some species, particularly for deuterated species, vary with temperature and evolutionary stage [1,2]. Consequently, abundances of deuterated species stored in cometary ices may provide insight into the physical conditions present at the time of chemical synthesis. If their volatiles formed in cold (<20 K) environments, such as on interstellar dust grains via H (or D) atom addition reactions, a strong D/H enrichment (relative to molecular hydrogen) would be expected. However, if volatiles formed primarily in warmer regions such as a protoplanetary disk atmosphere, then the D/H ratio in molecules would be much lower. In the solar system today, methane is stored as ice in cometary nuclei and is released via sublimation as comets approach the Sun. CH4 infrared emission has been detected through rotational–vibration transitions (m3 band) in the expanding atmospheres (comae) of many comets. The m4 band of CH3D and the m3 band of CH4 fall in the same spectral region (3.3 lm). Thus, both molecules can be ⇑ Corresponding author. Fax: +1 (314) 516 6152. E-mail address: [email protected] (E.L. Gibb). 0022-2852/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jms.2013.05.007

sampled simultaneously with available high-resolution astronomical spectrometers, in particular with NIRSPEC at the Keck II observatory and CRIRES at ESO’s VLT. Simultaneous sampling (i.e., covering the frequencies) of the two isotopic forms of methane eliminates many sources of systematic uncertainty, resulting in well-constrained CH3D/CH4 ratios (or upper limits). However, searches for monodeuterated methane (CH3D) have been conducted in only a few comets, thus far providing only upper limits for the abundance ratio CH3D/CH4. Nonetheless, these can enable significant tests of alternative scenarios for methane synthesis in the early solar system. The quality of such tests depends critically on details of the fluorescence models, solar spectrum, atmospheric transmittance model, and other such parameters used in the analysis. Kawakita and Watanabe [3] developed the first fluorescence model for the m4 band of CH3D and later applied it to several comets. We recently developed improved models for non-LTE fluorescence emission from multiple molecular band systems (primarily for application to cometary comae) utilizing the standard HITRAN molecular database with our customized updates. The details of the technique and its application to several molecules of importance to cometary science are discussed in Villanueva et al. [4]. Here, we present our fluorescence model based on improvements in the molecular parameters for the m4 band of CH3D and also in the pumping solar spectrum. The new model permits improved interpretation of archival cometary spectra and more accurate predictions for sensitivity limits (in terms of minimum detectable line

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124

flux) in future spectroscopic searches for monodeuterated methane. The model produces a set of robust fluorescent efficiencies (g-factors, expressed in photons s1 molecule1) at select rotational temperatures that are used to predict and interpret the intensities of CH3D lines in the infrared spectral region. In this paper, we discuss how the HITRAN 2008 database (Section 2) is utilized to develop a fluorescence emission model for CH3D (Section 3), including how we decouple extractions of spin temperature and rotational temperature at low temperatures (Section 3.3). We present the implementation of the model to synthesize emissions originating from cometary comae spanning a broad range of rotational temperatures (5–150 K) (Section 3), and we lastly discuss its application to interpret the infrared spectra of comet C/2007 N3 (Lulin) (Section 4, [5]). 2. The database For our CH3D fluorescence model, we utilized the following parameters from the HITRAN database: frequency (m [cm1]), line strength at room temperature (S(296 K) [cm1/(molecule cm2)]), lower state energy (E1 [cm1], the subscript ‘1’ indicates lower state, while ‘2’ denotes upper state), upper and lower ro-vibrational quantum numbers. We then calculated the statistical weights (Section 2.2), the partition functions (see Section 2.3 and Table 1), and the Einstein A- and B-coefficients for each rotational–vibrational transition (Sections 2.4 and 3.1). Using these parameters, we constructed a database of fluorescence efficiencies (g-factors, see Section 3.2and Table 2) for each transition (gi) for a range of rotational temperatures appropriate for cometary studies (5–200 K). The details of the calculations are provided below. In the most recent (2008) update of the HITRAN database, the parameters of 12CH4 were updated and new bands of CH3D (near 8, 2.9, and 1.56 lm) were added [6] using data from Nikitin et al. [7]. Spectroscopic information (frequencies and line intensities) for the CH3D band used in this work (m4) is from Nikitin et al. [8].

CH3D is a symmetric top molecule with C3v symmetry and has three non-degenerate vibrational modes [m1(A1), m2(A1), m4(A1)] and three degenerate modes [m4(E), m5(E), m6(E)]. The vibrational structure of CH3D can be seen in Fig. 2 of Nikitin et al. [7]. The rotational structure of the ground vibrational state is presented in Fig. 1. As mentioned above, we utilized the ro-vibrational energies and transition frequencies of individual lines in the HITRAN database when generating our g-factors. The rotational structure for the ground state of the m4 band of CH3D is shown in Fig. 1. For symmetric top molecules, J and K are the local quantum numbers that determine the rotational structure, where J is the total

Table 1 Partition functions for selected Trot. Trot (K)

20 50 75 100 125 150

angular momentum quantum number and K is the projection of the total angular momentum on the molecular symmetry axis – therefore, J can only have values greater than or equal to K. Fig. 1 shows that for a given sequence of levels with the same values of K, J P K. The other rotational quantum number is the total symmetry (C of the molecular symmetry group. For molecules like CH3D that have a threefold axis of symmetry, levels with K = 0, 3, 6, 9, . . . have a larger statistical weight than those with K = 1, 2, 4, 5, 7, 8, . . . (see Section 2.2). The stronger levels are denoted as A1 and A2 and the weaker levels as E. In Fig. 1, the A states are indicated in red and the E states in blue. The selection rules for the m4 band of CH3D are: DJ = J2  J1 = 1, 0, +1, DK = 0, K 6 J. Further, the symmetry selection rules require that only transitions between A , A and E , E are permitted while exchange between different spin species (A , E) are not. 2.2. Statistical weights With three symmetrically placed H nuclei, CH3D has two nuclear spin ladders (A and E) of which A is doubly degenerate. We calculated the statistical weights explicitly for our modeling due to an error found in the HITRAN database where the statistical weights for A states are missing a factor of two (HITRAN personnel are aware of this issue and it will be corrected in a forthcoming release). The total rotational statistical weight (generally denoted as ‘g’, however, we use ‘w’ in this work to distinguish it from the fluorescence efficiency g-factor) is given by wtot = wJwKwDwCwI, where wJ = 2J + 1 and wK = 1 (for K = 0) and wK = 2 (for K > 0). The deuterium nuclear spin weight (2ID + 1) is equal to 3 (ID = 1), and the carbon nuclear spin weight (2IC + 1) is equal to 1 (IC = 0). The hydrogen spin weighting factors are taken from Herzberg [9, p. 28] for the case of a symmetric top molecule with three identical atoms. In a different functional form, these weights are also given in McDowell [10]:

For A-states :¼ wI ¼

2.1. Energy levels and selection rules for CH3D

Total partition functionsa Direct sum

Analyticalb

HITRANc

94.8431 336.8470 612.5233 938.5211 1307.9950 1716.3854

91.9025 337.0910 612.9719 939.2116 1308.9600 1717.6538

– – 613.359596 939.169799 1308.897438 1717.547178

a The partition functions are calculated assuming Trot = Tspin (see Section 2.3 and [6]). b From [10]. c From HITRAN 2008 [6].

119

For E-states : wI ¼

1 ð2IH þ 1Þð4I2H þ 4IH þ 3Þ ¼ 4 3

1 ð2IH þ 1Þð4I2H þ 4IH Þ ¼ 2 3

ð1Þ

ð2Þ

where IH = 1/2. 2.3. Calculating the partition function HITRAN includes a file of total partition functions for molecules in the database from 70 to 3000 K. However, since cometary emission often must be described by temperatures lower than 70 K, the partition functions need to be calculated for lower temperatures. We calculated the partition functions explicitly for all temperatures in our modeling using the following procedure. Since cometary comae are cold (usually <150 K), generally only the lowest energy vibrational and electronic levels are populated and only the rotational partition function has a non-negligible contribution to the total partition function (Qtot(T) = QrotQvibQe, with Qvib  1.0 and Qe  1.0). The total rotational partition function is given by Herzberg (p. 505) [9] and McDowell [10] as

Q rot ¼

X

wi expðhcEi =kTÞ:

In general, if the HITRAN database is sufficiently complete, direct summation can be used to calculate the partition function using the statistical weights as calculated in Section 2.2. Over the temperature ranges applicable to comet observations (up to 150 K), this is often the case. To verify this, we also used the analytical solution (Eq. (24)) in McDowell [10] for a prolate symmetric rotator to calculate the partition function:

120

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124

Table 2 Fluorescence efficiencies of selected CH3D m4 lines at 1 AU for different models. Assignmenta

R

R3(4) R3(3) R2(2) R R1(2) R R0(2) R R1(1) R R0(1) R R0(0) P P2(2) P P3(3), PP1(4) P P3(4) R R

Wavenumber (cm1)

50 K g-factors

3061.41 3053.84 3044.29 3042.34 3040.41 3034.69 3032.74 3025.04 2997.82 2988.04 2980.21

75 K g-factors

3.63 9.31 4.67 4.80 5.56 7.81 7.06 6.63 4.93 7.78 2.89

5.00 9.33 3.91 3.66 4.16 4.94 4.51 3.83 3.57 7.33 3.77

g-Factors for 100 Kb Kawakita et al. [3]

This work (blackbody approximation of solar flux)

This work (solar spectrum)

4.72 7.73 4.14 2.65 2.87 3.20 2.89 2.40 2.53 5.94 4.06

5.19 8.23 3.22 2.83 3.21 3.45 3.20 2.57 2.67 6.26 3.78

5.21 8.25 3.23 2.84 3.22 3.46 3.21 2.57 2.68 6.27 3.79

a The superscript in the designation of individual lines indicates DK = +1 or 1 for R and P sub-bands, respectively. The subscript indicates the K value in the lower state. For example, the RR3(4) indicates the J1 = 4 line of the R branch of the sub-band K2 = 2 ? K1 = 3. See Chapter IV.2 of [9] for more details. b In units of 106 photons s1 molecule1.

800

A

E

E

A

E

Einstein A coefficients may be calculated explicitly using the formalism presented in Simecková et al. [12]:

E

A21 ¼ 8pcm2 Q tot ðT rot ÞS12 ðT rot Þ=f½1  expðhcm=kTÞ 600

 ½expðE1 hc=kTÞIa w2 g;

Energy [cm-1]

10 10

5

5

10

10

10

10 10

5

5

5

5

10

10

10

400

200

0 0

2

4

6

8

10

where E1 is the lower state energy, m is the line frequency, and Qtot(Trot) is the partition function calculated in Section 2.3. We also calculated the Einstein coefficients for induced absorption using B12 = S12(296 K)/[hm f(296 K)], where S12(296 K) is the line strength tabulated at 296 K [in cm1/(molecule cm2)] and scaled by the isotopic abundance for CH3D (6.1575  104), and f is the fractional population given by f = w1 exp(hcE1/kTrot)/Qtot(Trot), while w1 and w2 are the statistical weights of the lower and upper state, respectively. At a temperature of 296 K, the total band intensity for m4 is Sm4(296 K) = 3.97  1021 cm1 (molecule cm2)1. In Section 3, we use B12 coefficients and fractional population to compute pumping rates and A21 coefficients to compute branching ratios in building the CH3D fluorescence model.

K-Ladder Fig. 1. Rotational structure of the ground state of CH3D m4. The red lines corresponding to K = 0, 3, 6, 9 ladders indicate A states. The blue dashed lines, corresponding to the K = 1, 2, 4, 5, 7, 8 ladders, indicate E states. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Q rot  r ðpmÞ1=2 exp½bð4  mÞ=12b3=2 ½1 þ b2 ð1  mÞ2 =90 þ   ð1 þ dÞð1 þ q0 þ q1 b1 þ q2 b2 þ q3 b3 Þ; where b = hcB/kT, m = B/A, and r = (2IH + 1)3/3, d is the nuclear spin correction and q are the centrifugal distortion constants given in Table 1 of [10]. B and A are the rotational constants (B = 3.8801950 cm1 and A = 5.250821 cm1 for CH3D, [10]). Sample total partition functions for a direct summation over the HITRAN database, the analytical solution of [10], and the partition function file provided by HITRAN are compared in Table 1. Numerical results from these three methods are in close agreement.

3. CH3D m4 fluorescence emission model for cometary comae Emission processes in cometary atmospheres are non-LTE (where LTE denotes ‘‘local thermodynamic equilibrium’’). When comets enter the inner solar system and are warmed by the Sun, the ices sublimate, thereby releasing native (primary) volatiles that drag dust grains with them to form the coma. Solar radiation then induces infrared fluorescent emission. Interpreting these emission lines requires development of full quantum mechanical models for each molecular band used. 3.1. Fluorescence pumping In cometary atmospheres, molecules are pumped into excited vibrational states by solar radiation, which then decay via resonant fluorescence to levels in the ground vibrational state. Fluorescence emission rates (g-factors) are determined by first calculating the pumping rate from solar radiation (gpump [s1]) for each line in the database following the standard formalism:

2.4. Line intensities and Einstein coefficients

g pump ðJ 0 ; K 0 Þ ¼ rB12 F  ðms Þf :

CH3D has nine bands in the 3–5 lm regions, including three fundamentals (m1, m2, m4), three overtones (2m3, 2m5, 2m6), and three combination bands (m3 + m6, m3 + m5, m5 + m6) [11]. The line intensities at 296 K are provided in the HITRAN database. From these the

The source of pumping radiation is the solar flux (F  ) at the Doppler shifted frequency of each line as seen by each comet. Historically (c.f., Gibb et al. [13]), a 5760 K blackbody was assumed. More recently, a synthetic solar spectrum combining an empirical

ð3Þ

121

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124 6

g-factors (x10 -5 photons molecule -1 s -1)

model [14] and a stellar flux model [15] has been adapted (for details, see Appendix B in Villanueva et al. [4]). The synthetic solar spectrum in the region of interest is shown in Fig. 2. For comparison, we also overplot a 5770 K blackbody. We find, in general, that while an accurate solar model is important for constraining relative line strengths of some molecules, the solar spectrum has relatively few lines in our region of interest and so makes only a minor difference to the g-factors for CH3D. This is well illustrated for 100 K in Table 2 where g-factors for the blackbody and full spectrum model agree to within 0.4%. The total fluorescence pump for the m4 band of CH3D at 100 K was calculated to be 2.21  104 s1.

100 K 5

4 3

2

1

3100

3050

3.2. Fluorescence emission line intensities, and Einstein coefficients

g i ¼ g pump ðJ 0 ; K 0 Þ  A21 =Atot : Fig. 3 shows a plot of g-factors for the m4 CH3D band at 100 K (top) and 40 K (bottom). Also shown, for comparison, are g-factors in that spectral region for the m3 band of CH4. Selected g-factors for strong CH3D lines at three rotational temperatures (50 K, 75 K, and 100 K) are given in Table 2. The strongest line at rotational temperatures applicable to most comets is the RR3(3) line, where the superscript in the designation of individual lines indicates DK = +1 or 1 for R and P sub-bands, respectively. The subscript indicates the K value in the lower state and the number in parentheses is J1. For example, the RR3(3) indicates the J1 = 3 line of the R branch of the sub-band K2 = 2 ? K1 = 3. See Chapter IV.2 of [9] for more details. 3.3. Disentangling the g-factor dependency on spin temperature and rotational temperature

Solar Flux at 1 AU (10-27 J s cm-3)

Molecules like CH4 and CH3D contain multiple hydrogen nuclei that can have their nuclear spins oriented either parallel or antiparallel. This gives rise to what are called spin states, denoted A, E, or F for CH4 and A or E for CH3D [9]. Radiative transitions between different spin states are forbidden. For CH3D, the lowest energy level in the A-ladder lies 9.1 cm1 (13 K) below the lowest E level, so the A to E ratio is temperature dependent if we assume

2.4 2.2 2.0 1.8 1.6 1.4 3100

3050

3000

2950

2900

Frequency [cm-1] Fig. 2. Solar spectrum over the region of the m4 band of CH3D (black trace), used to generate the g-factors in Table 2. For comparison, a 5770 K Blackbody is also shown (solid red line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2950

6

g-factors (x10 -5 photons molecule-1 s-1)

Once the line-by-line pumping rates were calculated, the branching ratio for each downward transition (2,1) was calculated from the ratio of the probability of spontaneous emission for that transition relative to the total from that level, A21/Atot. Atot is the total of Einstein A21 coefficients contributing to the upper ro-vibrational level. The fluorescence emission rates (g-factors) for each line, i, were then calculated from

3000

Frequency (cm -1)

40 K 5

4 3

2 1 0 3100

3050

3000

2950

Frequency (cm -1) Fig. 3. Fluorescence g-factors (at heliocentric distance 1 AU) for lines in the m4 band of CH3D (black) and the m3 band of CH4 (red) shown at 100 K (upper panel) and 40 K (lower panel).

that this ratio was realized under LTE- conditions (or if it is intentionally made so in the laboratory). Thus a nuclear spin temperature can be formally derived from measured A/E ratios. The relationship between Tspin and A/E is strongly non-linear at low temperatures (below 40 K, Fig. 4). The interpretation of spin temperatures for cometary molecules has been a subject of long-term debate. Measurements of Tspin in comets are relevant to the fields of chemical physics (where understanding nuclear spin conversion is of interest) and to cosmogony (since the spin temperature is thought to indicate the chemical formation temperature of a molecule [16,17]). Importantly, the spin temperature, which reflects nuclear spin symmetries, is distinct from the gas rotational temperature (Trot). The latter is defined as the Boltzmann temperature that best approximates the rotational population distribution in the ground vibrational state of a molecule and reflects local conditions in the inner coma of a comet at the time of the observations. The g-factors of molecules with identical nuclei have separate dependencies of spin temperature and rotational temperatures. Without deconvolving spin and rotational temperature in the resulting g-factors, the relative line strengths cannot be accurately reproduced by our models, and this can influence the calculated abundances. The g-factor dependencies from Trot and Tspin have been quantified for CH4 [5,18] and H2O [18,19], C2H6 [4], and CH3OH [20,21]. For CH3D, we followed a similar procedure (described in detail by [18]) and derived how the g-factor would vary with E/ A, respectively Tspin, while Trot is treated as a parameter:

g A ðT rot ; E=AÞ ¼ g A ðT rot ; E=A ¼ 1Þ  2=ðE=A þ 1Þ;

ð4Þ

g E ðT rot ; E=AÞ ¼ g E ðT rot ; E=A ¼ 1Þ  2  ðE=AÞ=½ðE=AÞ þ 1;

ð5Þ

122

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124

the dominant source of the difference is likely the solar radiation density used. We used the synthetic solar spectrum discussed in Section 3.3 and in Villanueva et al. [4]. Kawakita and Watanabe [3] used the radiation density from Thekaekara [22], which covers the near-UV to visible range. In more recent work, their solar radiation field has been updated and in general, g-factors for molecules are now closer agreement (see for example Kawakita and Mumma [23]). Both models for m4 CH3D emission are still to be rigorously tested, pending a secure detection of multiple lines of CH3D in a comet.

1.0

E/A Ratio

0.8

0.6

0.4

0.2

4. Application to comets 0

10

20

30

50

40

60

Tspin (K) Fig. 4. The relative populations of E and A spin states (A/E) are shown as a function of temperature (Tspin = Trot). For comparison to other molecules with spin symmetries, see Fig. 4 of Villanueva et al. [4].

where gA (Trot; E/A) and gE (Trot; E/A) are the g-factors for A and E transitions for a given rotational temperature and an arbitrary E/A ratio, gA (Trot; E/A = 1) and gE (Trot; E/A = 1) are the g-factors for the same rotational temperature, but for statistical equilibrium (i.e., for E/A = 1, Tspin above 40 K). Note that the g-factors in Table 2 are calculated for Tspin = Trot. Eqs. (4) and (5) can be used to calculate g-factors for a different spin temperature.

3.4. Comparison with previous model for deuterated methane

1.0 0.8 0.6 0.4 0.2 0 CH4 at 68 K OH*

1000 800

Flux Density (10

-20

2

-1

W/ m / cm )

Telluric Transmittance

Kawakita and Watanabe [3] developed an earlier fluorescence emission model for the m4 band of CH3D. A comparison of g-factors from that work and this work for 100 K is shown in Table 2. With the exception of the RR2(2) and PP3(4) transitions, the Kawakita and Watanabe [3] g-factors are systematically lower than ours by about 7–10%. A small part of this difference could be attributed to updates in the CH3D molecular line database since 2000. However,

Comets are remnants from the early formation of the solar system. They are composed of a mix of volatiles (ices) and refractory grains that have been stored at very cold temperatures in either the Edgeworth–Kuiper Belt or the Oort cloud reservoirs. When they enter the inner solar system, solar radiation induces infrared fluorescence emission from coma gases. The near-infrared is ideal for studying comets since many simple, abundant molecules, including CH4 (which is a symmetric molecule having no dipole moment and so cannot be studied in the radio) and CH3D, emit in the 1–5 lm region and can be observed with sensitive, high-resolution spectrometers, such as NIRSPEC on Keck 2, CRIRES on the VLT, and IRCS on Subaru. Each of these instruments has sufficient spectral coverage to simultaneously sample CH4 m3 and CH3D m4 (see Figs. 3 and 5). CH4 has been commonly observed via high-resolution infrared spectroscopy in comets since its 1996 discovery in Hyakutake (C/1996 B2) [24]. CH3D has not been detected in comets to date, but we show below that significant upper limits are challenging models of the chemical evolution during solar system formation. We applied the model described in this paper to the moderately bright comet C/2007 N3 (Lulin). Since CH3D has not been detected, we generally assume statistical equilibrium A/E and determine the rotational temperature in the coma from other molecules. When possible, the uncertainty in gas rotational temperature is

600 400 200

100

CH3 D at 68 K, 10% model prediction

50 0 3080

3070

3060

3050

3040

3030

Observed Frequency [cm-1 ] Fig. 5. An illustration of fluorescence emission modeling in comets. The top panel shows the telluric transmittance over the spectral region of interest, convolved to the appropriate instrumental resolving power (k/Dk  24,000). The middle panel shows the residual spectrum (black) for comet Lulin on 31 January 2009 after subtracting the transmittance model (see Gibb et al. [5] for details). The red and blue lines are fluorescence emission models for CH4 and OH, respectively, at the best fit temperature of 68 K. The green is the 1-sigma error envelope. The bottom panel shows the residual spectrum after subtracting the CH4 and OH synthetic emission models with the 1-sigma error envelope overplotted in green, illustrating the quality of the fit. A CH3D/CH4 = 10% model is shown to illustrate the expected line positions and how our results challenge current models. See Section 4 for more details.

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124

propagated to the uncertainty in the derived upper limit for CH3D. We determined stringent upper limits of CH3D/CH4 <2.9% and <3.0% (3-sigma) on 31 and 30 January 2009, respectively [5]. A sample spectrum from our NIRPSEC observations of comet C/ 2007 N3 (Lulin) is shown in Fig. 5. One challenge to near-infrared ground-based studies is the synthesis of an accurate telluric model (top panel, Fig. 5). We used an implementation of the LBLRTM (Line-By-Line Radiative Transfer Model, [25]) to synthesize telluric spectra (see details in Villanueva et al. [26]). The best-fit synthetic fluorescence emission models for CH4 (red) and OH (blue) (Doppler shifted to the 54 km s1 velocity of Lulin at the time of the observations) are compared with the residual (i.e., continuum-subtracted) spectrum (middle panel). The bottom panel shows the residual for comet Lulin after subtracting the (modeled) CH4 and OH emission. If all fluorescence emission features are accounted for, this plot should represent the stochastic noise of the data (the 1-sigma error envelope is shown in green). For comparison, above the residual we show a synthetic model for a CH3D/CH4 ratio of 10%, consistent with predictions from astrochemical models by Aikawa and Herbst [27], Willacy and Woods [28], and Albertsson et al. [2] for dark molecular clouds and protoplanetary disks. At the 10% CH3D/CH4 abundance level, the intensities of the strongest CH3D lines are well above the error envelope for the observations (green line in the bottom panel), especially that of the line near 3054 cm1. In other words, if the astrochemical models correctly predicted CH3D abundance (relative to either methane or water), we should have had a clear detection of CH3D in comet Lulin. This illustrates that neither the restrictive upper limits we find for CH3D in C/2007 N3 (Lulin) nor the comparable upper limits that have been reported for C/2004 Q4 (Machholz) (<2%, [18,29]) and C/2001 Q2 (NEAT) [30] can be explained by current astrochemical modeling. Significant constraints on the CH3D/ CH4 ratio can be obtained even in the case of such moderately productive comets. Our results therefore provide a stringent test of existing chemical models of solar system formation.

5. Conclusion We present a fluorescence emission model for the m4 band of CH3D for application to cometary comae. Earlier fluorescence emission models were generated using this methodology for HCN [31], C2H6 m7 [4], C2H6 m5 [32], C2H2 (m3 and m2 + m4 + m5) [33], NH3 [23], NH2 [23], CH4 m3 [13], CO m1 [33], H2CO (m1, m5, m3 + m6) [20], CH3OH m3 [21] and m2 [34], H2O and HDO [26], and OH [19]. Application of our models to high-resolution, near-infrared spectra of comets permits retrieval of rotational and (where applicable) spin temperatures required for establishing accurate molecular production rates. We note that this technique could easily be applied to the other ro-vibrational bands of CH3D and used for any environment in which fluorescent emission may be important. For example, CH3D and the other molecules discussed above have been detected in some planetary atmospheres. The m2 + m3 band of CH3D was detected in Saturn’s atmosphere by Kim et al. [35] and in Titan’s atmosphere by Kim et al. [36] using high-resolution, near-infrared spectroscopy. The fluorescence excitation modeling described in this paper can be part of a comprehensive treatment of the complicated atmospheric chemical modeling used to explain features in these outer solar system bodies. While CH4 and CH3D have not yet been detected in protoplanetary disks, several other molecules have been, including H2O, HCN, C2H2, and CO2 (see for example Carr and Najita [37,38]; Mandell et al. [39]). As instrument sensitivities improve, additional molecules will be discovered. Disk atmospheres are subject to stellar radiation, and fluorescence modeling may be needed to understand

123

their emission spectra. The methodology described in this paper can be applied to such systems by replacing the solar flux in Eq. (3) with that for a T Tauri star of the appropriate spectral type. Fluorescence emission modeling of molecules has important applications for atmospheric science. The technique presented in this paper can be applied to comets, planetary atmospheres, exoplanet atmospheres, and protoplanetary disk studies. CH3D is an important molecule as it can be diagnostic of deuterium fractionation processes and low temperature chemistry. Acknowledgments We greatly appreciate the efforts of the HITRAN molecular database project and the numerous laboratories that have contributed experimental data and theoretical calculations that made this work possible. BPB and ELG were supported by NSF Planetary Astronomy grants AST-0807939 and AST-1211362. GLV, MAD and MJM acknowledge support from NASA’s Planetary Atmospheres and Astronomy Programs (08-PATM08-0031 (PI: GLV), 09-PATM09-0080 (PI: MAD), 08-PAST08-0033/34 (PI: MJM), 09-PAST09-0034 (PI: MAD)), and NASA’s Astrobiology Institute (NAI5/NNH08ZDA002C, PI: MJM). The spectra presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. References [1] Y. Aikawa, V. Wakelam, F. Hersant, R.T. Garrod, E. Herbst, Astrophys. J. 760 (2012) 40–58. [2] T. Albertsson, D.A. Semenov, Th. Henning, 2011, arXiv1110.2644A. [3] H. Kawakita, J.-I. Watanabe, Astrophys. J. 582 (2003) 534–539. [4] G.L. Villanueva, M.J. Mumma, K. Magee-Sauer, J. Geophys. Res. 116 (2011) 08012. [5] E.L. Gibb, B.P. Bonev, G. Villanueva, M.A. DiSanti, M.J. Mumma, E. Sudholt, Y. Radeva, Astrophys. J. 750 (2012) 102–115. [6] L.S. Rothman, I.E. Gordon, A. Barbe, D.C. Benner, P.F. Bernath, et al., J. Quant. Spectrosc. Radiat. Transfer 100 (2009) 533–572. [7] A.V. Nikitin, J.P. Champion, L.R. Brown, J. Mol. Spectrosc. 240 (2006) 14–25. [8] A.V. Nikitin, L.R. Brown, L. Féjard, J.P. Champion, V.G. Tyuterev, J. Mol. Spectrosc. 216 (2002) 225–251. [9] G. Herzberg, Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules, Krieger, Malabar, FL, 1945. [10] R.S. McDowell, J. Chem. Phys. 93 (1990) 2801–2811. [11] A.V. Nikitin, J.P. Champion, V.G. Tyuterev, L.R. Brown, J. Mol. Spectrosc. 184 (1997) 120–128. [12] M. Simecková, D. Jacquemart, L.S. Rothman, R.R. Gamache, A. Goldman, J. Quant. Spectrosc. Radiat. Transfer 98 (2006) 130–155. [13] E.L. Gibb, M.J. Mumma, N. Dello Russo, M.A. DiSanti, K. Magee-Sauer, Icarus 165 (2003) 391–406. [14] F. Hase, P. Demoulin, A.J. Sauval, G.C. Goon, P.F. Bernath, A. Goldman, J.W. Hannigan, C.P. Rinsland, J. Quant. Spectrosc. Radiat. Transfer 102 (2006) 450– 463. [15] Kurucz, 1997, . [16] J. Crovisier, Astron. Astrophys. 130 (1984) 361–372. [17] M.J. Mumma, H.A. Weaver, H.P. Larson, Astron. Astrophys. 187 (1987) 419– 424. [18] B.P. Bonev, M.J. Mumma, E.L. Gibb, M.A. DiSanti, G.L. Villanueva, K. MageeSauer, R.S. Ellis, Astrophys. J. 699 (2009) 1563–1572. [19] B.P. Bonev, Towards a Chemical Taxonomy of Comets: Infrared Spectroscopic Methods for Quantitative Measurements of Cometary Water (With an Independent Chapter on Mars Polar Science), Ph.D. Thesis, , 2005. [20] M.A. DiSanti, B.P. Bonev, K. Magee-Sauer, N. Dello Russo, M.J. Mumma, D.C. Reuter, G.L. Villanueva, Astrophys. J. 650 (2006) 470–483. [21] G.L. Villanueva, M.A. DiSanti, M.J. Mumma, L.-H. Xu, Astrophys. J. 747 (2012) 37–47. [22] M.P. Thekaekara, Appl. Opt. 13 (1974) 518–522. [23] H. Kawakita, M.J. Mumma, Astrophys. J. 727 (2011) 91–101.

124

E.L. Gibb et al. / Journal of Molecular Spectroscopy 291 (2013) 118–124

[24] M.J. Mumma, M.A. DiSanti, N. Dello Russo, M. Fomenkova, K. Magee-Sauer, C.D. Kaminski, D.X. Xie, Science 272 (4266) (1996) 1310–1314. [25] S.A. Clough, M.W. Shephard, E.J. Mlawer, J.S. Delamere, M.J. Iacono, K. CadyPereira, S. Boukabara, P.D. Brown, J. Quant. Spectrosc. Radiat. Transfer 91 (2005) 233–244. [26] G.L. Villanueva, M.J. Mumma, B.P. Bonev, R.E. Novak, R.J. Barber, M.A. DiSanti, J. Quant. Spectrosc. Radiat. Transfer 113 (2012) 202–220. [27] Y. Aikawa, E. Herbst, Astron. Astrophys. 351 (1999) 233–246. [28] K. Willacy, P.M. Woods, Astrophys. J. 703 (2009) 479–498. [29] H. Kawakita, H. Kobayashi, Astrophys. J. 693 (2009) 388–396. [30] H. Kawakita, J.-I. Watanabe, R. Furusho, T. Fuse, D. Boice, Astrophys. J. 623 (2005) L49–L52. [31] H. Lippi, G.L. Villanueva, M.A. DiSanti, H. Bohnhardt, M.J. Mumma, B.P. Bone, D. Prialnik, A&A 551 (2013) 51–59.

[32] Y.L. Radeva, M.J. Mumma, G.L. Villanueva, M.F. A’Hearn, ApJ 729 (2011) 135– 143. [33] G.L. Villanueva, M.J. Mumma, M.A. DiSanti, B.P. Bonev, E.L. Gibb, K. MageeSauer, G.A. Blake, C. Salyk, Icarus 216 (2011) 227–240. [34] M.A. DiSanti, B.P. Bonev, G.L. Villanueva, M.J. Mumma, Astrophys. J. 763 (2013) 1–18. [35] J.H. Kim, S.J. Kim, T.R. Geballe, S.S. Kim, L.R. Brown, Icarus 185 (2006) 476–486. [36] S.J. Kim, T.R. Geballe, K.S. Noll, R. Courtin, Icarus 173 (2005) 522–532. [37] J.S. Carr, J.R. Najita, Science 319 (2008) 1504–1506. [38] J.S. Carr, J.R. Najita, Astrophys. J. 733 (2011) 102–119. [39] A.M. Mandell, J. Bast, E.V. van Dishoeck, G.A. Blake, C. Salyk, M.J. Mumma, G. Villanueva, Astrophys. J. 747 (2012) 92–105.