The Water Deuteration as a Function of Radius in Massive Star-forming Region W33A with Herschel HIFI and the
IRAM 30m Telescope
Hylke D. Hoogland June 29, 2018
Abstract
Context: Deuteration of water occurs most efficiently in cold (<30K) environments. By measuring the deuteration throughout the protostellar cloud, we can understand the conditions that are prevalent in the protostellar cloud. The deuteration as a function of radius for high- mass protostellar sources has been studied only a few times.
Aims: In my master thesis, I aim to determine the degree of water deuteration as a function of radius for the protostellar source W33A. This will give insight into how the water has formed in this source. Furthermore I aim to discover how the turbulent and in-fall velocities change throughout the cloud.
Method: I use Herschel HIFI and IRAM 30m spectral data in the sub-millimeter and far infrared from 203 GHz to 1910 GHz. I will see which isotopologues of water are present in the data, derive abundances for these and derive a D/H value by means of radiative transfer modeling with the RADEX and RATRAN codes.
Results: In total, 33 spectral features of water were found, though some uncertain: 7 H2O, 7 H218
O, 9 HDO, 1 HD18O and 1 D2O lines. This is one of the first times that D2O was found in a high-mass protostellar source. Some lines show a P Cygni profile indicative of an expanding envelope. Radiative transfer modeling with RATRAN shows that a step function abundance profile reproduces the observations best, with the abundances [HDOin]=10-6.5±1.0, [HDOout]=10-9±1.0, [H218
Oin]=10-6±1.0, [H218
Oout]=10-7.5±1.0and therefore D/Hin=10-3.7±1.0 and D/Hout=10-4.5±1.0. The observations show little evidence of in-falling material.
Conclusions: A high inner D/H ratio may be explained by a large amount of water sublimating from the dust grains as the protostar heats the surrounding material. The water may be formed in the pre-stellar phase. The cold conditions will have favoured a high D/H ratio in the water.
1 Introduction
Water is a very important molecule in the ISM for various reasons. A significant fraction of the oxygen in grain mantle ice is contained in water, from ∼40% to ∼65%, though in the gas phase this is much less (<1%) (Bergin et al.,2000). Oxygen occurs in numerous molecules, in ices (e.g. H2O, CO2), rocks (e.g. silicates) and organics (e.g. methanol, acetaldehyde) (see also Huss & Draine (2007)). Therefore water plays a role in the chemistry of these molecules. Due to its asymmetry, water has many rotational transitions which can cool molecular clouds. Cool gas is necessary to allow a cloud to collapse to form stars. Water is also a very important solvent in biological processes and thus very important for life. Hence it is scientifically important to understand where, when and under what conditions water forms (see for exampleCeccarelli et al. (2014)).
The heavy isotope fractionation observed in water is helpful in determining this. Due to a lower zero point vibrational energy, the reaction barrier of heavy isotopes is lower than of lighter isotopes at temperatures lower than ∼30K (Ceccarelli et al.,2014). In several lab experiments, it has been shown that water can form very efficiently on cold dust grains (e.g. Hiraoka et al.(1998),Mokrane et al. (2009), Dulieu et al. (2010), Oba et al. (2012)). At low temperatures (<30K), the water molecules can be more deuterated than would be the case in the gas phase. Therefore, finding a high deuterium fraction may be a sign of water formed on dust grains. A lower deuterium fractionation would indicate more formation of water in the gas phase. However, Furuya et al.
(2016) have shown that the (D2O/HDO)/(HDO/H2O) ratio is a better indicator for the formation
history of water than the D/H ratio. D2O is however very rarely detected in massive star forming regions (see for example,Neill et al.(2013)), making determining this ratio very difficult.
Many computer simulations of deuterium fractionation predict higher D/H ratios in the cold outskirts of the cloud (e.g. Cazaux et al.(2011),Aikawa et al.(2012),Furuya et al.(2016)). Only at very low densities, below 104 cm-3, does the inner D/H ratio become comparable to the outer D/H ratio (Cazaux et al., 2011). These models and observations from, for example, Roberts &
Millar(2007), also show that water has a low deuterium fractionation compared to other molecules.
Organic molecules can reach a D/H ratio as high as a few 10%.
Chemical abundances are derived from the spectrum using radiative transfer codes. These radiative transfer codes depend on molecular collision data determined from lab experiments (for a complete overview of lab experiments for astrochemistry, seeSmith(2011)). Although this allows for a calculation of the chemical abundances, the codes are all limited by the completeness of the observations and the accuracy of the molecular collision data. See van der Tak (2011) for a review on this topic. Before starting radiative transfer calculations, it is often useful to assume Local Thermodynamic Equilibrium (LTE). The energy levels are then distributed according to the Boltzmann distribution. This allows for analysis with rotational diagrams (Goldsmith & Langer, 1999). LTE is not always a good assumption. In this thesis I use the non-LTE radiative transfer code RADEX (van der Tak et al.,2007) to estimate line intensities. RATRAN (Hogerheijde & van der Tak,2000) is used here to determine the chemical abundances as a function of radius because it allows for a source model to be used.
Although there have been many studies performed on low-mass sources (e.g. van Dishoeck et al.
(1995), Jørgensen & van Dishoeck (2010), Taquet et al. (2013), Coutens et al. (2013a), Coutens et al.(2013b),Coutens et al.(2014), or the review ofLuhman(2012) on low-mass star formation), due to the scarcity of high-mass sources there have been significantly less studies on high-mass sources (e.g. van der Tak et al.(2006),Herpin et al.(2012),Neill et al.(2013); also see the recent review byMotte et al. (2017) on massive star formation). Therefore, it is very useful to perform this research on more high-mass sources, to study whether there is a difference between high- and low-mass protostars.
Deuterium fractionation has also been studied in comets. Altwegg et al. (2017) have ana- lyzed data from the Rosetta mission, allowing them to determine the HDO/H2O and D2O/HDO ratio of 67P/Churyumov-Gerasimenko: 1.04x10-4 and 1.8x10-2, respectively. Revealing that the (D2O/HDO)/(HDO/H2O) ratio mentioned in section 1 is 17 for this comet, much higher than for protostellar objects which range from 0.5 to 7. For many other comets, this was not possible.
But the isotopic ratios of, for example, N, C and O could be determined (e.g. Eberhardt et al.
(1995),Biver et al.(2016)). It appears that solar system comets generally have a lower D/H than protostellar objects by a factor ∼10, though the D/H values show some overlap. Some comets (e.g. 67P/Churyumov-Gerasimenko, C/2012 F6 Lemmon) have comparable D/H ratios (∼5x10-4) to some protostellar objects (e.g. IRAS 4B, NGC 1333-IRAS 2A at ∼3x10-4) (Altwegg et al., 2017). For a review on the chemical composition of comets and a comparison with hot corinos, circumstellar disks and dense cloud cores, seeMumma & Charnley(2011).
My master thesis is structured as follows. Section2gives a short overview of the source W33A, after which the observations used in my thesis are described in section3. Data reduction methods are described in section4. Results follow from section5onwards, showing the line fitting procedure I used. Section6will then describe an attempt to model the data assuming LTE and optically thin lines according the the method described inGoldsmith & Langer (1999). More realistic attempts are described in sections7and8, where non-LTE radiative transfer codes RADEX and RATRAN are used respectively (Hogerheijde & van der Tak, 2000). Discussion and conclusions follow in sections9 and 10. The appendix shows a list of spectral lines found (Table 3), the python code used to automate RATRAN and a list of observations (Table4).
2 Source description
W33A is a well-studied Massive Young Stellar Object (MYSO), containing at least 2 main cores, named MM1 and MM2 (van der Tak et al.,2000). 3 43 GHz sources were detected in W33A later on (van der Tak & Menten,2005). The dust structure was determined by de Wit et al.(2007) on scales down to 120 AU using VLTI/MIDI data. The source is located at a distance of 2.40 kpc in the larger W33 cloud complex, located in the Scutum spiral arm. There are several water masers in the W33 cloud complex, that have allowed its distance to be determined by parallax measurements
(Immer et al., 2013). The gas in the envelope of W33A shows complex kinematics. Maud et al.
(2017) have identified multiple cores, rotating filaments and a protostellar outflow. Until recently the stellar mass of the MYSO in W33A was estimated to be about 10 M (Galván-Madrid et al., 2010), but ALMA observations show that the main core of W33A, MM1, consists of 7 YSO’s of which the most massive has a stellar mass of 7 M and a luminosity of 104 L (Izquierdo et al., 2018).
The chemical composition of W33A has been a subject of study several times before. A great variety of ice-mantle species was identified byGibb et al.(2000). They have determined the column density for H2O and NH3. A summary of column densities for many ice species (e.g. H2O, HDO, CO, CO2, CH4and organic compounds) found in W33A by earlier studies can also be found in their paper. Solid HDO was detected in dust grain ice mantles around W33A (Teixeira et al., 1999), supporting the idea that highly deuterated water forms on dust grains. This detection is however disputed (Dartois et al.,2003). Also the gas phase chemistry has been well-studied, for example by Immer et al.(2014). They found that the chemical complexity increases as the protostar evolves.
They identified two line flux ratios which may be used as chemical clocks to trace the evolution of protostars: N2H+(3-2)/CS(6-5) and N2H+(3-2)/H2CO(42,2-32,1).
3 Observations
The data used in this thesis were obtained by the HIFI instrument (Heterodyne Instrument for the Far-Infrared) (de Graauw & Helmich,2001) aboard the Herschel Space Observatory (Pilbratt, 2003) and the IRAM 30m telescope. HIFI data was taken from the Herschel Science Archive (HSA1). In total 131 HIFI single point pipeline products of W33A were used. The integration times range from 3 minutes to 1 hour and 13 minutes. These observations were taken between March 2010 and October 2012 as part of the Water In Star-forming regions with the Herschel space observatory (WISH) program (van Dishoeck et al., 2011)). The observations are listed in Table4. The beam size of the Herschel HIFI data ranges from 44.2" at 480 GHz to 11.1" at 1910 GHz (Roelfsema et al., 2012). The spectral resolution of the HIFI data ranges from 0.8 km/s at low frequencies (around 500 GHz) to 0.3 km/s at high frequencies (around 1800 GHz). The system temperature of the instrument ranges from 100K at 480 GHz to 1000K at 1280 GHz and constant thereafter. The noise level of the observations depend on the integration time, number of overlapping observations and frequency, ranging from 0.01K to 0.1K. A calibration uncertainty of 10% was applied (Roelfsema et al.,2012).
The HIFI instrument is a heterodyne receiver, meaning that the data consists of two sets of observations at slightly different frequency, called the Lower Side Band (LSB) and Upper Side Band (USB). These two sets are convolved together to form the final dataset. This causes features from different frequencies to overlap. A sideband deconvolution was attempted by removing lines from observations where other observations show only continuum at the same frequency. In this way, the probability of false identifications is reduced (see section4).
4 Data reduction
A baseline correction was performed for all of the data products and the continuum set to zero.
The data products were all averaged to obtain a final spectrum used for the analysis in this thesis.
All observations were taken at different frequencies. Most of the misplaced spectral features (see section 3) could be removed by excluding spectral features in observations at frequencies where all other observations show only continuum. An example of this procedure is shown in Figure1.
The ’true spectrum’ is a collection of random Gaussian line profiles, not actual observations. The green and red lines show what the LSB and USB detectors of HIFI could observe of this spectrum.
Below, the overplotted green and red lines show how these data are convolved to create the final dataset shown in yellow. Notice how lines are filtered out in the overlapping regions if they appear at the wrong frequency. In most parts of the spectrum, I have 4-6 observations, which allow the true spectrum to be reconstructed.
1http://archives.esac.esa.int/hsa/whsa/
Figure 1: An example of the HIFI sideband deconvolution procedure for 2 observations.
The HIFI instrument has a main beam efficiency which changes throughout the frequency range which still needed to be corrected for. This efficiency is described by the following equation (Roelfsema et al.,2012):
ηmb= ηmb,0· e−(4πσλ )2 (1)
where ηmb,0 = 0.76 ± 0.02 and σ = 3.8 ± 0.9µm for λ in µm. For band 5 (1120-1410 GHz) these numbers change to ηmb,0= 0.66 ± 0.02 and σ = 3.8 ± 0.9µm.
The IRAM 30m data, which contains the H218O line at 203 GHz was added to the final spectrum after averaging all the HIFI data.
5 Line Fitting
5.1 Procedure
After data reduction, the whole spectrum is searched by a python programme for emission and absorption lines. The data for spectral lines was taken from the Splatalogue, which makes use of the CDMS (Müller et al.,2005), JPL (Pickett et al.,1998), Lovas/NIST (Lovas & Dragoset,2004) and ToyaMa (Kobayashi,2014) spectral line databases. It was assumed that any line within 10 km/s of the systemic velocity is due to (an isotopologue of) water. This is a reasonable first assumption because water is a very common molecule in the ISM. The systemic velocity was taken fromImmer et al.(2013) to be ∼37 km/s, based on the water maser G012.90-0.26 in the W33A cloud. Lines found by this method revealed that a systemic velocity of ∼39 km/s was more consistent with the peaks of the emission lines, though this is only slightly different from the value reported byImmer et al.(2013).
Individual components in the line profiles were identified by eye and a manual first guess was given for all components. A python code was used to find the best fitting Gaussian parameters for each component, starting from these first manual guesses. The identification of some rare isotopologues (H217O, D2O and HD18O) is uncertain and therefore the analysis described in the rest of my thesis will be needed to check whether the identification is correct (see sections6,7and 8).
5.2 Qualitative examination of line profiles
Spectral features from various water isotopologues have been identified: 7 H2O, 7 H218O, 0 or 1 H217O, 9 HDO, 1 to 3 HD18O and 1 or 2 D2O lines. A list of emission and absorption lines with fitting parameters can be found in Table3. In this table, the frequency is given in GHz, the Doppler shift and line width in km/s and the line height (peak intensity) in K. ID1 is the ID given to a single spectral feature, while ID2 is a sub-ID given to the individual components of a spectral feature.
A wide variety of line shapes can be seen, indicating the different conditions that prevail in the regions of the cloud that the photons are being emitted from. Most notable are the P Cygni profiles which can be seen in many lines. These are indicative of an expanding envelope. As the expanding material is moving towards Earth, it absorbs emission from the source at a lower velocity (with respect to Earth). An example is shown in Figure2. Notice that for all figures showing spectral lines in my thesis, the continuum has been removed. A few remarkable examples are highlighted in Figures2,3,4 and5, all other lines are shown in Figures6and7.
Figure 2: The IRAM 30m data of the H218O line at 203 GHz shows a clear P Cygni profile
Figure 3: The transition of H2O at 557 GHz shows very intricate absorption patterns and a clear distinction between the broad emission component assumed to originate from the outflow and a narrow emission component assumed to originate from the envelope
Very intricate absorption patterns can be seen in the transition of H2O at 557 GHz. This is because this transition is easily absorbed due to its lower level energy being zero (Figure3). Also noteworthy is the great range in line widths seen in the spectrum (see Figures4 and 5). This is likely due to turbulence.
Figure 4: The HDO line at 849 GHz shows a large broadening, likely due to turbulence.
Figure 5: This D2O line shows clearly in comparison with Figure4 that the range of line widths is large.
Figure 6: Other water lines found which were not already mentioned in section5.
Figure 7: Continuation of Figure6. Together these figures show all water lines found.
6 LTE modeling: rotational diagrams
As a first step the lines which were found in the spectrum were fitted with a Gaussian profile, consisting of a narrow, low velocity component and a broad, high velocity component (see section 5.1. The fitting parameters were used to make rotational diagrams according to the method from Goldsmith & Langer(1999):
ln(Nu
gu
) = ln N − ln Z − Eu
kBT (2)
where Nu is the upper level column density, N the total molecular column density (of the species), Euthe upper level energy, guthe degeneracy of the upper energy level, kBthe Boltzmann constant, T the temperature and Z the partition function:
Z =X
gue−EuT (3)
where the sum is over all energy levels of the species.
This method assumes Local Thermal Equilibrium (LTE) and lines which are optically thin. It was however quickly determined that these assumptions do not correspond to the data, as can be seen in Figure8. The data points derived from the data do not follow a straight line. Therefore, the column densities and temperatures derived by this method are unreliable for this data set.
Assuming temperatures in the range 50-100K showed that the H217O abundance would be higher than H2O, which made me conclude that the H217O identification is incorrect. There was no reason yet (see sections7 and8) to exclude any other lines.
Figure 8: Fit to the HDO data, assuming LTE and optically thin lines. Numbers above correspond to ID1, numbers below to ID2 in Table3.
7 Non-LTE modeling
The non-LTE radiative transfer and excitation solver RADEX (van der Tak et al.,2007) was used to obtain estimates for the temperature, molecular column density and collision density. RADEX does not assume that all lines are optically thin, as is the case for the rotational diagram analysis (see section 6). Instead, it computes the optical depth for each line. RADEX uses molecular collision data from the LAMDA database (Schöier et al.,2005).
First guesses for the temperature, molecular column density and collision density were taken from the rotational diagram analysis. RADEX uses a single value model to calculate the flux of the emission lines. Although RADEX modeling provided estimates for the density and temperature of the emitting region, the models did not fit the data very well, hinting at a far more complex cloud structure. The differences between models and observations are given in the following as the Root Mean Squared difference (RMS):
RM S = rD
(model − observations)2E
(4) The RMS difference (in units of log(Ngu
u)) for HDO was roughly equal for the rotational diagram analysis and the RADEX models: 2.0±0.8 for the rotational diagrams and 2.2±0.8 for the RADEX models. For H218O this was 3.5±1.3 and 2.5±1.3, respectively. The same general trend is observed for H2O: 1.0±0.2 for the rotational diagrams and 0.7±0.2 for the RADEX models. Despite not assuming LTE, the difference between the models and the observations is only slightly less than for the rotational diagram analysis. A combination of 2 RADEX models was attempted to obtain a better fit. These did however not result in significantly lower residuals and were therefore not considered further.
Table 1 shows the parameters found with RADEX. These were not significantly different for the various isotopologues, so a general view of the results is presented here. It is however expected that, especially for the column density, the difference between the isotopologues would be larger.
Therefore, RATRAN modeling was performed afterwards. This takes into account a more realistic model. An example of a rotational diagram made from the output of RADEX is shown in Figure 9.
Table 1: Overview of the parameters found with RADEX
Parameter Search range Best models
Temperature (K) 10-500 250-300
Molecular column density (log cm-2) 8.5-15 12-13.5 Collision density (of H2) (log cm-3) 3.5-12 8-10
With RADEX, I determined that one of the D2O lines was incorrectly identified. It was not possible to construct a model that would reproduce both D2O lines I had tentatively identified.
Only the 3(0,3)-2(1,2) line at 851 GHz could be reproduced.
Figure 9: RADEX output for HDO. There was one line for which a broad and narrow component could be observed. This did not influence the RADEX results.
8 Radiative transfer modeling with RATRAN
With RATRAN, the shapes of the emission lines were modeled using the physical structure of W33A determined by van der Tak et al. (2013), as described in section 4.1 of that paper. This physical structure consists of 50 concentric shells spaced logarithmically with radius.
Figure 10: The density of ortho-H2from the physical structure used for RATRAN modeling.
Figure 11: The density of para-H2from the physical structure used for RATRAN modeling.
Figure 12: The temperature from the physical structure used for RATRAN modeling. The black and blue lines show the upper level energy (in K) of the transitions observed for H218O and HDO, respectively.
Like RADEX, RATRAN uses molecular collision data from the LAMDA database (Schöier et al., 2005). For a complete description of the RATRAN code, seeHogerheijde & van der Tak (2000). The density, Doppler b parameter (a measure of turbulence) and in-fall velocity of this model were adjusted to fit the narrow Gaussian component of the line fits (where these could be separated from the broad component). RATRAN assumes spherical symmetry, while it is assumed that the broad component is due to outflows. Therefore the broad component is ignored and the outflows are not further discussed in my thesis. For a recent review on protostellar outflows, see Bally(2016).
The abundance of water and its isotopologues was determined by comparing the total line flux of the narrow component of the line fits with the total flux predicted by the RATRAN model. The kinematics were determined first by setting the means of both the line fit and the model to 1. This allows for direct comparison. In this case I took the root mean square of the difference between the model and the fit for both the abundance and kinematics.
8.1 Methodology
The methodology used for RATRAN is as follows. I have used three different shapes for the abundance profiles: a power-law, a step function and a constant profile. These were chosen because they have at most two parameters that define the entire profile. This makes it easier to check whether the best fitting model has been found. The step of the step function is where the physical structure fromvan der Tak et al.(2013) is ∼100K. This is assumed to be the temperature at which water sublimates from the dust grains into the gas phase (see alsoCeccarelli et al.(2000),van der Tak et al.(2006),Coutens et al.(2012)). There are some studies which favour higher sublimation temperatures (e.g. Coutens et al.(2014)), but due to limited time this was not further explored in my thesis.
The profiles for the kinematics (turbulent and in-fall velocity) are all power-law profiles. This is because this seems more physical to me than a step function. Having the same profile shape for the kinematics will also allow me to compare the kinematics profiles for different isotopologues and abundance profiles. If the kinematics profiles can be well determined, these should be comparable.
The best fitting profile finding procedure is fully automatic (see Appendix for Python code).
A brute force approach is adopted. All profiles start flat; the abundance profiles all start at log10(o−HX
2) = −10 (where X is the isotopologue being studied), the kinematics profiles at 0 km/s.
After calculation of this first model, the difference between the flux and line shape of the RATRAN model and the fits to the observed lines determine how much the model profiles will be changed for the next iteration. The change is random within this range. Only the shape of the profile is conserved (power-law, step function or constant). If the new model corresponds better to the observations, the new model is saved. Afterwards this new model is changed again randomly the same way as in the first step. The smaller the difference with the observations, the smaller the random change of the profiles. In this way, the code will search for the best fitting profiles of each shape for each isotopologue automatically.
Of all the isotopologues found, only H218O, HDO, HD18O and D2O were attempted to model.
H2O was excluded because the expected high optical depth of these lines would be very com- putationally expensive. Therefore, the D/H ratio will be determined from the H218O and HDO models. HD18O and D2O resulted in higher abundances of these isotopologues than HDO, which is very unlikely. Therefore, I decided that 2 of the HD18O and 1 of the D2O line identifications were unreliable. RADEX modeling had already shown that this D2O line was likely incorrectly identified (see section7). This leaves only 1 certain identification for both of these isotopologues and therefore no abundance profiles can be derived. For H217O, only 1 line was found. This is the only line which appears in just one observation, making it uncertain whether it has come from the LSB or the USB (see section4for details). However, even if the identification would be certain, 1 line is not enough to derive an abundance profile for the entire range of the physical structure.
8.2 Results
8.2.1 Abundance
The Figures 13 and 14 show the derived best fitting profiles for both the power-law and step function profiles for H218O and HDO, respectively. The constant profile is no longer considered because it was not as well in agreement with the observations as the power-law and step function profiles for both isotopologues (see Table2). Because the power-law and step function shapes are based on the abundance relative to ortho-H2, there is a slight curve in the models when comparing the abundance of the species to the total H2abundance. The difference between the flux from the RATRAN models and the observations is given in Table2, as the Root Mean Square:
RM S = rD
(RAT RAN − observations)2E
(5) For both these isotopologues, the RMS for the step function models is 30% lower than for the power-law models and 55% lower than for the constant models. This indicates that there is an abundance jump in both H218O and HDO around 100K in the protostellar cloud.
As expected, the RATRAN modeling predicts more of these water isotopologues in the centre of the cloud (>100K) than at the edge of the cloud (<100K). This is because the water is thought to sublimate from the dust grains on which it formed at these temperatures. Both the power- law and step function profiles indicate this. The uncertainty in the abundances for H218O and
HDO are estimated from the range of profiles which yield similar residuals with the observations.
The uncertainty is thus determined to be about a factor 10, or one order of magnitude, for both isotopologues.
Table 2: RMS difference RATRAN and observations, in units of log10(K km s-1) Species Step function Power-law Constant
H218O 0.468 0.672 1.033
HDO 0.390 0.558 0.703
Figure 13: The best fitting step function (red) and power-law (blue) profiles for the abundance of H218O, as determined by RATRAN.
Figure 14: The best fitting step function (red) and power-law (blue) profiles for the abundance of HDO, as determined by RATRAN.
As the step function agrees most with the observations (see Table2), the flux diagrams for these models is shown in Figures15and16for H218O and HDO, respectively. Clearly, the models do not represent the observations perfectly (also see Table2). This is likely due to the simplicity of the model, being restricted to a power-law or step function shape. Using a more realistic model, like the model derived inIzquierdo et al.(2018), would be an interesting next step for future research.
Their model includes a 3D gas structure with filaments feeding the main protostar, which is more realistic than the assumed spherical symmetry and simple abundance profiles I used here.
Figure 15: The flux of each line from the RATRAN step function model for H218O compared to the observations.
Figure 16: The flux of each line from the RATRAN step function model for HDO compared to the observations.
Figure 17: The ratio of the step function (red) and power-law (blue) profiles for the abundance of H218O and HDO as shown in Figures13and14, as determined by RATRAN.
From the H218O and HDO power-law and step function models, the H218O/HDO ratio can be derived, shown in Figure17. This will allow us to derive the D/H ratio. The step function and power-law profiles do not agree. However, the uncertainty in the abundances is estimated from Figures18 and19to be a factor of ∼10, or one order of magnitude. This is larger than the difference between the ratios derived from the power-law and step function abundance profiles.
Therefore, Figure17suggests that the H218O/HDO ratio is more or less constant throughout the cloud.
To assess the significance of the result, the D/H ratio is calculated for all RATRAN models of H218O and HDO, assuming16O/18O=500. This is shown in the Figures18and19for the power- law and step function profiles, respectively. The y-axis shows the D/H ratio as calculated from the innermost shell of the physical structure, while the x-axis shows this for the outermost shell. From these figures I estimate the uncertainty for both profiles to be about 1 order of magnitude, i.e. a factor ∼10.
Figure 18: A comparison of all power-law abundance profiles of H218O and HDO reveals that these predict that the D/H ratio is around 10-4for both the inner envelope as well as the outer envelope.
The color coding shows the RMS difference in the log of the fluxes between the RATRAN model and the observations, where red is a large RMS and blue a small RMS. White areas have no data.
Figure 19: A comparison of all step function abundance profiles of H218O and HDO reveals that these predict that the D/H ratio is around 10-3.7 for the inner envelope and 10-4.5 for the outer envelope. The color coding shows the RMS difference in the log of the fluxes between the RATRAN model and the observations, where red is a large RMS and blue a small RMS. White areas have no data.
8.2.2 Kinematics
Besides deriving the abundance profiles as described above, RATRAN also enables us to determine kinematics for the emitting gas in the cloud. As stated above (see section8.1), only power-law profiles were employed for determining the turbulent and in-fall velocities, as this seemed more physical to me than step function profiles. Despite the great variety of line shapes (see section 5.2), a clear result can not be obtained by the RATRAN modeling. All kinematics profiles derived in correspondence with the H218O and HDO power-law and step function abundance profiles are depicted in Figure20.
It seems that most models agree that the envelope is expanding (positive in-fall velocity), though one of them has an inverted slope. An expansion of the outskirts was expected from the qualitative analysis of the line shapes in section5.2. The disagreement between these in-fall velocity profiles may be explained by the lack of clearly asymmetrical emission lines found in the spectrum. Only the HDO line at 509 GHz shows clear asymmetry. For the other lines, only the lack of asymmetry can be used as an argument, namely that the in-fall velocity must be low. The general trend is that the in-fall velocity is slowly increasing with radius (larger expansion at larger radius) and the turbulent velocity is slowly decreasing with radius (see Figure21).
Figure 20: The profiles for the kinematics (turbulent velocity in blue, in-fall velocity in red) derived by RATRAN.
Figure 21: The best fitting kinematics profiles (turbulent velocity in blue, in-fall velocity in red) with uncertainty margins, to show the general trend.
Another explanation maybe an oversimplification of the models used in my thesis. FromMaud et al.(2017), it is clear that the structure of W33A is far more complex than a simple spherically symmetric power-law. A more realistic model would be the one derived byIzquierdo et al.(2018), which includes a 3D gas structure with gas flows in the feeding filaments. Furthermore it may be that if there is a part with a high in-fall velocity, it is not well represented in the observations. It is likely (and observed, e.g. Galván-Madrid et al. (2010), Maud et al. (2017)) that the centre of the cloud is collapsing. The density in the centre is very high, and therefore emission lines from this region may be obscured.
9 Discussion
9.1 On the abundance
The total abundance found for HDO is very similar to other sources studied, e.g. IRAS 16293-2422 (Ceccarelli et al.(2000),Parise et al.(2005),Coutens et al.(2012)), NGC 1333 IRAS2A (Liu et al.,
2011) and AFGL 2591 (van der Tak et al.,2006). The IRAS sources are low-mass protostars while AFGL 2591 is a high-mass protostar. The inner HDO abundances (with respect to H2) range from 8x10-8to 1.8x10-7, the outer HDO abundances range from 8x10-11to 7x10-10for the IRAS sources, for AFGL 2591 this is a bit higher at 4x10-9. The inner HDO abundance found for W33A falls within this range, at 10-6.5±1. The outer HDO abundance for W33A is a bit higher than for the IRAS sources, but very similar to AFGL 2591, at 10-9±1. The HDO abundance seems to be a bit higher in the outer envelopes of the high-mass sources, but this difference with the low-mass sources is very small (less than 1 order of magnitude).
There is a difference in the D/H ratio with radius with NGC 1333 IRAS2A and AFGL 2591, where it is increasing with radius. In the case of NGC 1333 IRAS2A, the HDO/H2O ratio increases from 0.01 to 0.07, while for AFGL 2591 this same ratio increases from 5x10-4to 4x10-2. The results derived in my master thesis are more similar to IRAS 16293-2422, which has been observed to have a higher D/H ratio in the inner part of the protostellar cloud (inner part: 0.03, outer part: 0.005) (Parise et al.,2005).
These results may be explained if large amounts of water formed in the pre-stellar phase. The cold conditions will have enhanced deuterium fractionation. As the protostar evolves and heats up the cloud, the water will be released from the dust grains, enriching the gas with highly deuterated water.
9.2 On the kinematics
Although the results from the kinematics in my master thesis are questionable due to the problems described earlier (see section8.2.2), a comparison is made with the results achieved inHerpin et al.
(2012) on W43-MM1 and Herpin et al. (2016) on NGC 6334I(N), DR21(OH), IRAS 16272-4837, and IRAS 05358+3543. They have found that a multi-step model for the kinematics works best to reproduce their data. Because W43-MM1 is a high mass source, I will focus the discussion mostly on that source. A comparison of the results from Herpin et al. (2012) for W43-MM1 with the results from this work for W33A are shown in Figure22.
Figure 22: The results fromHerpin et al.(2012) show decreasing in-fall velocity (pink dashed) and increasing turbulent velocity (black solid) for W43-MM1, contrary to what was found for W33A in this work (red for in-fall velocity, blue for turbulent velocity).
Similarly to W43-MM1, in all of the other sources mentioned the turbulent velocity is increasing with radius. W33A seems to show a different behavior, where the turbulent velocity is decreasing with radius, though not significantly. However, the turbulent velocities for all of these sources, including W33A, is around a few km/s, showing no large fluctuation between the sources. Except for IRAS 05358+3543, a constant model would predict negative in-fall velocity (collapse) (Herpin et al.,2016). Therefore, IRAS 05358+3543 seems more similar to W33A than the other sources in that respect.
10 Summary and Conclusions
In my master thesis, I have used spectral data from the Herschel HIFI instrument and the IRAM 30m telescope to determine the abundance of water isotopologues in the massive star-forming region W33A. I have identified 7 H2O, 7 H218O, 9 HDO, 1 HD18O and 1 D2O lines. Only H218O and HDO lines were suitable for the radiative transfer modeling I did.
I used 3 shapes for the abundance profiles: a power-law, a step function and a constant function.
I developed a python code to run RATRAN models automatically. The shape that reproduces the observations best was found to be the step function. This indicates that there is a sublimation point for water at ∼100K. The abundances (with respect to H2) found for this profile are:
[HDOin]=10-6.5±1.0 [HDOout]=10-9±1.0 [H218Oin]=10-6±1.0 [H218Oout]=10-7.5±1.0
which results in the following D/H with radius:
D/Hin=10-3.7±1.0 D/Hout=10-4.5±1.0
For the power-law profiles, which correspond slightly less to the observations, the D/H ratio is constant with radius:
D/Hin=10-4.0±1.0 D/Hout=10-4.0±1.0
The roughly constant D/H ratio is different from the expectation that the D/H ratio should be higher on the outside due to lower temperatures. This can be explained if large amounts of water formed in the pre-stellar phase on cold dust grains. This water is released into the gas phase as the star heats the surrounding gas and dust, enriching the inner cloud with highly deuterated water.
This master thesis was supervised by prof. Floris van der Tak.
References
Aikawa Y., Wakelam V., Hersant F., Garrod R. T., Herbst E., 2012,Astrophysical Journal, 760, 40
Altwegg K., et al., 2017,Monthly Notices of the Royal Astronomical Society,469, S130 Bally J., 2016,Annual Reviews Astronomy and Astrophysics,54, 491
Bergin E. A., et al., 2000,Astrophysical Journal, Letters,539, L129 Biver N., et al., 2016,Astronomy and Astrophysics,589, A78
Cazaux S., Caselli P., Spaans M., 2011,Astrophysical Journal, Letters,741, L34
Ceccarelli C., Castets A., Caux E., Hollenbach D., Loinard L., Molinari S., Tielens A. G. G. M., 2000, Astronomy and Astrophysics,355, 1129
Ceccarelli C., Caselli P., Bockelée-Morvan D., Mousis O., Pizzarello S., Robert F., Semenov D., 2014, Protostars and Planets VI,pp 859–882
Coutens A., et al., 2012,Astronomy and Astrophysics, 539, A132 Coutens A., et al., 2013a,Astronomy and Astrophysics,553, A75 Coutens A., et al., 2013b,Astronomy and Astrophysics,560, A39
Coutens A., et al., 2014,Monthly Notices of the Royal Astronomical Society,445, 1299
Dartois E., Thi W.-F., Geballe T. R., Deboffle D., d’Hendecourt L., van Dishoeck E., 2003,As- tronomy and Astrophysics,399, 1009
Dulieu F., Amiaud L., Congiu E., Fillion J.-H., Matar E., Momeni A., Pirronello V., Lemaire J. L., 2010,Astronomy and Astrophysics,512, A30
Eberhardt P., Reber M., Krankowsky D., Hodges R. R., 1995, Astronomy and Astrophysics,302, 301
Furuya K., van Dishoeck E. F., Aikawa Y., 2016,Astronomy and Astrophysics, 586, A127 Galván-Madrid R., Zhang Q., Keto E., Ho P. T. P., Zapata L. A., Rodríguez L. F., Pineda J. E.,
Vázquez-Semadeni E., 2010,Astrophysical Journal,725, 17 Gibb E. L., et al., 2000,Astrophysical Journal,536, 347
Goldsmith P. F., Langer W. D., 1999,Astrophysical Journal,517, 209 Herpin F., et al., 2012,Astronomy and Astrophysics,542, A76 Herpin F., et al., 2016,Astronomy and Astrophysics,587, A139
Hiraoka K., Miyagoshi T., Takayama T., Yamamoto K., Kihara Y., 1998,Astrophysical Journal, 498, 710
Hogerheijde M. R., van der Tak F. F. S., 2000, Astronomy and Astrophysics,362, 697 Huss G. R., Draine B. T., 2007,Highlights of Astronomy,14, 353
Immer K., Reid M. J., Menten K. M., Brunthaler A., Dame T. M., 2013,Astronomy and Astro- physics,553, A117
Immer K., Galván-Madrid R., König C., Liu H. B., Menten K. M., 2014,Astronomy and Astro- physics,572, A63
Izquierdo A. F., Galván-Madrid R., Maud L. T., Hoare M. G., Johnston K. G., Keto E. R., Zhang Q., de Wit W.-J., 2018,Monthly Notices of the Royal Astronomical Society,
Jørgensen J. K., van Dishoeck E. F., 2010,Astrophysical Journal, Letters,725, L172 Kobayashi K., 2014
Liu F.-C., Parise B., Kristensen L., Visser R., van Dishoeck E. F., Güsten R., 2011,Astronomy and Astrophysics, 527, A19
Lovas F. J., Dragoset R. A., 2004
Luhman K. L., 2012,Annual Reviews Astronomy and Astrophysics,50, 65
Maud L. T., Hoare M. G., Galván-Madrid R., Zhang Q., de Wit W. J., Keto E., Johnston K. G., Pineda J. E., 2017,Monthly Notices of the Royal Astronomical Society,467, L120
Mokrane H., Chaabouni H., Accolla M., Congiu E., Dulieu F., Chehrouri M., Lemaire J. L., 2009, Astrophysical Journal, Letters, 705, L195
Motte F., Bontemps S., Louvet F., 2017, preprint, (arXiv:1706.00118)
Müller H. S. P., Schlöder F., Stutzki J., Winnewisser G., 2005, Journal of Molecular Structure, 742, 215
Mumma M. J., Charnley S. B., 2011,Annual Reviews Astronomy and Astrophysics, 49, 471 Neill J. L., Crockett N. R., Bergin E. A., Pearson J. C., Xu L.-H., 2013, Astrophysical Journal,
777, 85
Oba Y., Watanabe N., Hama T., Kuwahata K., Hidaka H., Kouchi A., 2012,Astrophysical Journal, 749, 67
Parise B., et al., 2005,Astronomy and Astrophysics,431, 547
Pickett H. M., Poynter R. L., Cohen E. A., Delitsky M. L., Pearson J. C., Müller H. S. P., 1998, Journal of Quantitiative Spectroscopy and Radiative Transfer,60, 883
Pilbratt G. L., 2003, in Gry C., Peschke S., Matagne J., Garcia-Lario P., Lorente R., Salama A., eds, ESA Special Publication Vol. 511, Exploiting the ISO Data Archive. Infrared Astronomy in the Internet Age. p. 31
Roberts H., Millar T. J., 2007,Astronomy and Astrophysics,471, 849 Roelfsema P. R., et al., 2012,Astronomy and Astrophysics,537, A17
Schöier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J. H., 2005, Astronomy and Astrophysics, 432, 369
Smith I. W. M., 2011,Annual Reviews Astronomy and Astrophysics,49, 29
Taquet V., López-Sepulcre A., Ceccarelli C., Neri R., Kahane C., Coutens A., Vastel C., 2013, Astrophysical Journal, Letters, 768, L29
Teixeira T. C., Devlin J. P., Buch V., Emerson J. P., 1999, Astronomy and Astrophysics,347, L19 de Graauw T., Helmich F. P., 2001, in Pilbratt G. L., Cernicharo J., Heras A. M., Prusti T., Harris R., eds, ESA Special Publication Vol. 460, The Promise of the Herschel Space Observatory. p. 45 de Wit W. J., Hoare M. G., Oudmaijer R. D., Mottram J. C., 2007,Astrophysical Journal, Letters,
671, L169
van Dishoeck E. F., Blake G. A., Jansen D. J., Groesbeck T. D., 1995,Astrophysical Journal,447, 760
van Dishoeck E. F., et al., 2011,Publications of the Astronomical Society of the Pacific,123, 138 van der Tak F., 2011, in Cernicharo J., Bachiller R., eds, IAU Symposium Vol. 280, The Molecular
Universe. pp 449–460 (arXiv:1107.3368),doi:10.1017/S1743921311025191 van der Tak F. F. S., Menten K. M., 2005,Astronomy and Astrophysics,437, 947
van der Tak F. F. S., van Dishoeck E. F., Evans II N. J., Blake G. A., 2000,Astrophysical Journal, 537, 283
van der Tak F. F. S., Walmsley C. M., Herpin F., Ceccarelli C., 2006,Astronomy and Astrophysics, 447, 1011
van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007,Astronomy and Astrophysics, 468, 627
van der Tak F. F. S., et al., 2013,Astronomy and Astrophysics,554, A83
Appendix
A: Line list
Table 3: List of fitting parameters for the water lines found in the spectrum of W33A. Those with the same ID1 belong to the same spectral feature, for which ID2 identifies the various components.
The parameters are for Gaussian fits to each single component. See for details section 5.1.
ID1 ID2 Species Freq. (GHz) Shift (km/s) Err. (km/s) Width (km/s) Err. (km/s) Height (K) Err. (K) QN log(A (/s)) Eu(K) gu
0 0 H218O 203.408 39.306 0.640 3.045 1.028 0.306 0.410 3(1,3)-2(2,0) -5.318 203.684 7
0 1 H218O 203.408 35.640 5.602 4.314 1.467 -0.228 0.293 3(1,3)-2(2,0) -5.318 203.684 7
1 0 H218O 489.054 48.320 0.049 2.311 0.049 0.162 0.003 4(2,3)-3(3,0) -4.162 429.645 27
2 0 HDO 490.597 38.951 0.152 2.439 0.153 0.071 0.004 2(0,2)-1(1,1) -3.280 66.432 5
3 0 HD18O 501.567 24.733 0.066 1.740 0.066 0.109 0.004 1(1,0)-1(0,1) -2.660 46.255 3
4 0 HDO 509.292 31.929 3.515 4.703 2.586 0.015 0.004 1(1,0)-1(0,1) -2.635 46.755 3
4 1 HDO 509.292 38.393 0.303 1.957 0.478 0.040 0.012 1(1,0)-1(0,1) -2.635 46.755 3
5 0 H218O 547.676 39.047 0.605 2.960 0.495 0.062 0.004 1(1,0)-1(0,1) -2.483 60.462 9
5 1 H218O 547.676 33.664 0.806 1.708 0.679 -0.024 0.011 1(1,0)-1(0,1) -2.483 60.462 9
6 0 H2O 556.936 38.160 0.224 13.739 0.360 0.243 0.019 1(1,0)-1(0,1) -2.465 60.964 9
6 1 H2O 556.936 36.617 0.612 5.753 0.546 0.442 0.119 1(1,0)-1(0,1) -2.465 60.964 9
6 2 H2O 556.936 23.972 0.030 0.750 0.041 -0.237 0.011 1(1,0)-1(0,1) -2.465 60.964 9
6 3 H2O 556.936 29.871 0.269 2.272 0.195 -0.544 0.075 1(1,0)-1(0,1) -2.465 60.964 9
6 4 H2O 556.936 33.844 0.201 1.561 0.194 -0.758 0.124 1(1,0)-1(0,1) -2.465 60.964 9
6 5 H2O 556.936 37.245 0.145 1.431 0.301 -0.867 0.154 1(1,0)-1(0,1) -2.465 60.964 9
6 6 H2O 556.936 39.764 0.122 1.042 0.155 -0.698 0.168 1(1,0)-1(0,1) -2.465 60.964 9
6 7 H2O 556.936 41.770 0.103 0.749 0.068 -0.374 0.071 1(1,0)-1(0,1) -2.465 60.964 9
6 8 H2O 556.936 44.179 0.042 0.592 0.065 -0.159 0.021 1(1,0)-1(0,1) -2.465 60.964 9
8 0 HDO 599.927 38.329 0.137 2.131 0.134 0.104 0.006 2(1,1)-2(0,2) -2.462 95.224 5
10 0 H2O 752.033 36.978 0.158 7.786 0.222 0.713 0.030 2(1,1)-2(0,2) -2.156 136.938 5
10 1 H2O 752.033 37.479 0.046 2.049 0.065 1.265 0.034 2(1,1)-2(0,2) -2.156 136.938 5
11 0 HDO 753.411 38.367 0.090 1.810 0.089 0.159 0.007 3(1,2)-3(0,3) -2.229 167.561 7
13 0 HDO 848.962 38.642 0.065 3.041 0.065 0.167 0.003 2(1,2)-1(1,1) -3.033 83.631 5
14 0 D2O 850.758 43.196 0.079 1.348 0.081 0.080 0.004 3(0,3)-2(1,2) -1.774 101.358 7
15 0 HDO 893.639 33.749 0.113 1.411 0.113 -0.199 0.013 1(1,1)-0(0,0) -2.078 42.888 3
15 1 HDO 893.639 41.037 0.432 2.161 0.460 0.064 0.011 1(1,1)-0(0,0) -2.078 42.888 3
16 0 HDO 919.311 38.663 0.075 2.366 0.075 0.159 0.004 2(0,2)-1(0,1) -2.806 66.432 5
19 0 HDO 984.138 38.627 0.101 2.243 0.101 0.170 0.006 4(1,3)-4(0,4) -1.967 263.271 9
20 0 H2O 987.927 37.442 0.189 9.444 0.238 0.546 0.019 2(0,2)-1(1,1) -2.237 100.846 5
20 1 H2O 987.927 38.725 0.032 2.127 0.042 1.499 0.024 2(0,2)-1(1,1) -2.237 100.846 5
21 0 H218O 994.675 39.565 0.331 2.470 0.326 0.162 0.018 2(0,2)-1(1,1) -2.220 100.609 5
22 0 HDO 995.412 39.130 0.202 1.340 0.193 0.168 0.022 3(0,3)-2(1,2) -2.152 131.403 7
18
25 0 H2O 1097.365 36.243 0.171 8.578 0.214 0.465 0.018 3(1,2)-3(0,3) -1.788 249.436 21
25 1 H2O 1097.365 37.931 0.031 2.330 0.043 1.250 0.020 3(1,2)-3(0,3) -1.788 249.436 21
26 0 H218O 1101.698 34.404 1.358 2.658 0.785 -0.222 0.392 1(1,1)-0(0,0) -1.748 52.873 3
26 1 H218O 1101.698 38.144 6.645 3.454 2.139 0.161 0.255 1(1,1)-0(0,0) -1.748 52.873 3
27 0 H2O 1113.343 41.534 0.457 9.029 0.233 0.485 0.018 1(1,1)-0(0,0) -1.739 53.432 3
27 1 H2O 1113.343 24.157 0.026 0.639 0.029 -0.654 0.025 1(1,1)-0(0,0) -1.739 53.432 3
27 2 H2O 1113.343 29.132 0.179 1.698 0.101 -1.061 0.079 1(1,1)-0(0,0) -1.739 53.432 3
27 3 H2O 1113.343 33.773 0.106 1.981 0.252 -1.684 0.041 1(1,1)-0(0,0) -1.739 53.432 3
27 4 H2O 1113.343 37.337 0.163 1.161 0.193 -1.392 0.237 1(1,1)-0(0,0) -1.739 53.432 3
27 5 H2O 1113.343 39.509 0.105 0.914 0.099 -1.577 0.206 1(1,1)-0(0,0) -1.739 53.432 3
27 6 H2O 1113.343 41.608 0.070 0.668 0.049 -0.766 0.058 1(1,1)-0(0,0) -1.739 53.432 3
27 7 H2O 1113.343 44.185 0.027 0.320 0.029 -0.462 0.030 1(1,1)-0(0,0) -1.739 53.432 3
29 0 H218O 1655.868 34.139 0.081 1.560 0.082 -0.387 0.017 2(1,2)-1(0,1) -1.263 113.646 15
30 0 H2O 1661.008 33.062 0.152 2.134 0.163 -0.506 0.030 2(2,1)-2(1,2) -1.518 194.095 15
30 1 H2O 1661.008 40.083 0.061 1.502 0.063 1.065 0.036 2(2,1)-2(1,2) -1.518 194.095 15
32 0 H2O 1669.905 40.769 1.031 8.482 0.496 0.607 0.066 2(1,2)-1(0,1) -1.256 114.378 15
32 1 H2O 1669.905 24.329 0.043 0.672 0.050 -1.040 0.068 2(1,2)-1(0,1) -1.256 114.378 15
32 2 H2O 1669.905 28.747 0.168 1.457 0.200 -1.720 0.612 2(1,2)-1(0,1) -1.256 114.378 15
32 3 H2O 1669.905 33.777 0.777 2.838 1.313 -2.794 0.225 2(1,2)-1(0,1) -1.256 114.378 15
32 4 H2O 1669.905 38.456 0.282 1.711 0.715 -2.283 1.350 2(1,2)-1(0,1) -1.256 114.378 15
32 5 H2O 1669.905 39.940 0.093 0.597 0.165 -0.968 0.643 2(1,2)-1(0,1) -1.256 114.378 15
32 6 H2O 1669.905 41.751 0.125 0.769 0.106 -1.435 0.543 2(1,2)-1(0,1) -1.256 114.378 15
32 7 H2O 1669.905 44.336 0.046 0.415 0.054 -0.746 0.080 2(1,2)-1(0,1) -1.256 114.378 15
B: Automated code for RATRAN
The python code for running RATRAN automated consists of three parts. The program that is run in the terminal is runratran.py:
# −∗− c o d i n g : u t f −8 −∗−
"""
C r e a t e d on Mon Apr 23 2 2 : 2 3 : 1 8 2018
@author : H y l k e D. Hoogland
"""
from __future__ import d i v i s i o n import numpy a s np
import s c i p y . ndimage . f i l t e r s a s s f import o s
import i m p o r t l i b import t i m e
import m a t p l o t l i b . p y p l o t a s p l t
23
s f r = i m p o r t l i b . import_module ( ’ s k y f i t s r e a d ’ ) def f ( x , a , b , c , d ) :
return a ∗x∗∗3+b∗x∗∗2+ c ∗x+d
s p e c i e s = ’ m o l e c u l e ’ #e . g . hdo , h218o , e t c .
o r i g d a t a = np . l o a d t x t ( ’ path / t o / o r i g i n a l /RATRAN/ model . mdl ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ ) metadata = o r i g d a t a [ : 6 ]
rmax = f l o a t ( metadata [ 0 ] . s p l i t ( ’= ’ ) [ − 1 ] ) columns = metadata [ 3 ] . s p l i t ( ’ , ’ )
columns [ 0 ] = columns [ 0 ] . s p l i t ( ’= ’ ) [ − 1 ] d a t a = [ ]
f o r i in xrange ( len ( o r i g d a t a [ 6 : ] ) ) :
d a t a . append ( o r i g d a t a [ 6 : ] [ i ] . s p l i t ( ) ) d a t a = np . t r a n s p o s e ( d a t a ) . a s t y p e ( f l o a t )
#c o u n t t h e number o f s k y i n p u t f i l e s , i f c o m b i n i n g m o l e c u l a r d a t a f i l e s s k i e s = o s . l i s t d i r ( ’ . / ’ )
s q = 0
f o r i in s k i e s :
i f i [ : 3 ] == ’ sky ’ and i [ − 8 : ] == ’ n a m e o f i n p u t f i l e . i n p ’ : s q += 1
#t h i s p a r t i s t o e x t r a c t t r a n s i t i o n i n f o r m a t i o n from t h e LAMDA d a t a f i l e s i f s q == 2 :
Eup = [ ] l i n e s = [ ] t r a n s = [ ]
f o r z in [ ’ o ’ , ’ p ’ ] :
s k y d a t a = np . l o a d t x t ( ’ sky ’+z+ ’ n a m e o f i n p u t f i l e . i n p ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ )
f o r i in xrange ( len ( s k y d a t a ) ) :
i f s k y d a t a [ i ] . s p l i t ( ’= ’ ) [ 0 ] == ’ t r a n s ’ :
linum = ( s k y d a t a [ i ] . s p l i t ( ’= ’ ) [ 1 ] . s p l i t ( ’ , ’ ) ) l i n e s 0 = [ ]
t r a n s 0 = [ ]
l d a t = np . l o a d t x t ( ’ path / t o / molec / ’+z+ ’ h 2 o @ d a n i e l . d a t ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ )
k = 0
f o r j in xrange ( len ( l d a t ) ) : i f ’ ! ’ in l d a t [ j ] :
24
d l t r a n s 0 = j e l i f k == 5 :
d l t r a n s 1 = j e l i f k == 6 :
d l t 0 = j e l i f k == 7 :
d l t 1 = j f o r j in np . a r a n g e ( d l t 0 +1 , d l t 1 , 1 ) :
l i n e s 0 . append ( l d a t [ j ] . s p l i t ( ) ) f o r j in np . a r a n g e ( d l t r a n s 0 +1 , d l t r a n s 1 , 1 ) :
t r a n s 0 . append ( l d a t [ j ] . s p l i t ( ) )
l i n e s 1 = np . a r r a y ( [ l i s t ( x ) f o r x in zip ( ∗ l i n e s 0 ) ] ) t r a n s 1 = np . a r r a y ( [ l i s t ( x ) f o r x in zip ( ∗ t r a n s 0 ) ] ) l i n e s . append ( l i n e s 1 )
t r a n s . append ( t r a n s 1 ) Eup0 = [ ]
f o r i in linum :
y e s = l i n e s 1 [ 0 ] == i
Eup0 . append ( l i n e s 1 [ 5 ] [ y e s ] [ 0 ] ) Eup . append ( np . a r r a y ( Eup0 ) . a s t y p e ( f l o a t ) ) Eup = np . c o n c a t e n a t e ( Eup )
e l i f s q == 1 :
s k y d a t a = np . l o a d t x t ( ’ s k y n a m e o f i n p u t f i l e . i n p ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ ) f o r i in xrange ( len ( s k y d a t a ) ) :
i f s k y d a t a [ i ] . s p l i t ( ’= ’ ) [ 0 ] == ’ t r a n s ’ :
linum = s k y d a t a [ i ] . s p l i t ( ’= ’ ) [ 1 ] . s p l i t ( ’ , ’ ) l i n e s = [ ]
t r a n s = [ ]
i f s p e c i e s in [ ’ hdo ’ , ’ hd18o ’ ] :
l d a t = np . l o a d t x t ( ’ path / t o / molec / hdo . d a t ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ ) e l i f s p e c i e s in [ ’ d2o ’ ] : #add o t h e r s p e c i e s t h i s way
l d a t = np . l o a d t x t ( ’ path / t o / molec / o h 2 o @ d a n i e l . d a t ’ , dtype=str , d e l i m i t e r= ’ a s d r j k r 5 ’ )
k = 0
f o r j in xrange ( len ( l d a t ) ) : i f ’ ! ’ in l d a t [ j ] :
k += 1 i f k == 4 :
d l t r a n s 0 = j e l i f k == 5 :
d l t r a n s 1 = j
25
e l i f k == 7 : d l t 1 = j f o r j in np . a r a n g e ( d l t 0 +1 , d l t 1 , 1 ) :
l i n e s . append ( l d a t [ j ] . s p l i t ( ) ) f o r j in np . a r a n g e ( d l t r a n s 0 +1 , d l t r a n s 1 , 1 ) :
t r a n s . append ( l d a t [ j ] . s p l i t ( ) )
l i n e s = np . a r r a y ( [ l i s t ( x ) f o r x in zip ( ∗ l i n e s ) ] ) t r a n s = np . a r r a y ( [ l i s t ( x ) f o r x in zip ( ∗ t r a n s ) ] ) Eup = [ ]
f o r i in linum :
y e s = l i n e s [ 0 ] == i
Eup . append ( l i n e s [ 5 ] [ y e s ] [ 0 ] ) Eup = np . a r r a y ( Eup ) . a s t y p e ( f l o a t )
#l o o p RATRAN m o d e l s f o r e v e r while True :
constmod = np . l o a d t x t ( ’ . . / h2 ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ )
c o n v d a t = np . l o a d t x t ( ’ . . / c o n v d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ ) . T kinmin = c o n v d a t [ 8 ] == np . nanmin ( c o n v d a t [ 8 ] )
i f len ( c o n v d a t [ 8 ] [ kinmin ] ) > 0 and np . i s f i n i t e ( np . nanmin ( c o n v d a t [ 8 ] ) ) :
d o p p l e r = ( c o n v d a t [ 5 ] [ kinmin ]− c o n v d a t [ 4 ] [ kinmin ] ) / 5 0 ∗ np . a r a n g e (50)+\
c o n v d a t [ 4 ] [ kinmin ]
v = ( c o n v d a t [ 7 ] [ kinmin ]− c o n v d a t [ 6 ] [ kinmin ] ) / 5 0 ∗ np . a r a n g e (50)+\
c o n v d a t [ 4 ] [ kinmin ] e l s e :
d o p p l e r = np . z e r o s ( 5 0 ) v = np . z e r o s ( 5 0 )
#d o p p l e r = np . l o a d t x t ( ’ . . / db ’+ s p e c i e s +’ p r o f i l e s h a p e . d a t ’ ) #uncomment i f db and v s h o u l d
#v = np . l o a d t x t ( ’ . . / v ’+ s p e c i e s +’ p r o f i l e s h a p e . d a t ’ ) #b e f o r c e d t o h a v e t h e s e v a l u e s a l l s t d = np . l o a d t x t ( ’ . . / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . da t ’ )
a l l s t d 2 = np . l o a d t x t ( ’ . . / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e 2 . da t ’ ) i f a l l s t d [ − 1 ] > 3 :
a l l s t d [ − 1 ] = 3
#uncomment i f random g a u s s i a n c h a n g e s s h o u l d b e made t o p h y s i c a l s t r u c t u r e
"""
r = np . random . r a n d i n t ( l e n ( constmod ) ) #t o t a l l y f r e e p a r a m e t e r s R = np . random . normal ( )
p r i n t r , R
constmod [ r ] += R∗ a l l s t d [ −1]
constmod = s f . g a u s s i a n _ f i l t e r 1 d ( constmod , np . random . normal ( ) )
26
"""
rhon = np . random . normal ( ) ∗ a l l s t d [ −1] #power−l a w r h o x = np . random . normal ( ) ∗ a l l s t d [ −1]
dy = ( rhox−rhon ) / 5 0 ∗ np . a r a n g e (50)+ rhon constmod += dy
"""
#uncomment i f random s t e p f u n c t i o n c h a n g e s s h o u l d b e made t o p h y s i c a l s t r u c t u r e
"""
c o r e = np . random . normal ( ) ∗ a l l s t d [ −1] #s t e p e n v e = np . random . normal ( ) ∗ a l l s t d [ −1]
lncm = l e n ( constmod ) constmod [ : lncm / 2 ] += c o r e constmod [ lncm / 2 : ] += e n v e
"""
#uncomment i f random c o n s t a n t c h a n g e s s h o u l d b e made t o p h y s i c a l s t r u c t u r e
"""
dy = np . random . normal ( ) ∗ a l l s t d [ −1] #c o n s t a n t constmod += dy
"""
print constmod [ 0 ] , constmod [ − 1 ]
ki nr ms = np . mean ( a l l s t d 2 [ np . a r a n g e ( 0 , 2 ∗ len ( Eup ) , 2 ) ] ) i f np . i s n a n ( ki nrm s ) :
ki nr ms = a l l s t d [ − 1 ] / 4
#c h a n g e s t o k i n e m a t i c s a r e a l w a y s power−l a w s h a p e d dbn = np . random . normal ( ) ∗ ki nr ms ∗3
dbx = np . random . normal ( ) ∗ ki nr ms ∗3 dy = ( dbx−dbn ) / 5 0 ∗ np . a r a n g e (50)+ dbn
#n e g a t i v e t u r b u l e n t v e l o c i t y d o e s n o t make any s e n s e while any ( d o p p l e r+dy < 0 ) :
dbn = np . random . normal ( ) ∗ ki nr ms ∗3 dbx = np . random . normal ( ) ∗ ki nr ms ∗3 dy = ( dbx−dbn ) / 5 0 ∗ np . a r a n g e (50)+ dbn d o p p l e r += dy
print d o p p l e r [ 0 ] , d o p p l e r [ − 1 ] vn = np . random . normal ( ) ∗ ki nr ms ∗3 vx = np . random . normal ( ) ∗ ki nr ms ∗3 dy = ( vx−vn ) / 5 0 ∗ np . a r a n g e (50)+ vn v += dy
print v [ 0 ] , v [ − 1 ]
27
i f s q == 2 :
t 0 = t i m e . t i m e ( )
#run RATRAN f o r combined m o l e c u l a r d a t a f i l e s o s . system ( ’ amc␣amcow33a . i n p ’ )
o s . system ( ’ amc␣amcpw33a . i n p ’ ) o s . system ( ’ sky ␣ skyow33a . i n p ’ ) o s . system ( ’ sky ␣ skypw33a . i n p ’ ) t 1 = t i m e . t i m e () − t 0
#c a l l t h e c o d e t h a t r e a d s t h e s k y f i t s f i l e s
a l l s t d o , r o t r m s l o g o , r m s s t d l o g o , a l l E o , al l W f o , allWeo , allWmo , f i g s o = \ s f r . readshow ( l i n e s [ 0 ] , t r a n s [ 0 ] , s p e c i e s , ’ o ’ )
a l l s t d p , r o t r m s l o g p , r m s s t d l o g p , a l l E p , allWfp , allWep , allWmp , f i g s p = \ s f r . readshow ( l i n e s [ 1 ] , t r a n s [ 1 ] , s p e c i e s , ’ p ’ )
a l l s t d = np . h s t a c k ( ( a l l s t d o , a l l s t d p ) )
r o t r m s l o g = np . h s t a c k ( ( r o t r m s l o g o , r o t r m s l o g p ) ) r m s s t d l o g = np . h s t a c k ( ( r m s s t d l o g o , r m s s t d l o g p ) ) a l l E = np . h s t a c k ( ( a l l E o , a l l E p ) )
a l l W f = np . h s t a c k ( ( a l l W fo , a l l W f p ) ) allWe = np . h s t a c k ( ( allWeo , allWep ) ) allWm = np . h s t a c k ( ( allWmo , allWmp ) )
f i g s = np . h s t a c k ( ( f i g s o , f i g s p ) ) wl = 1/ r m s s t d l o g
u b e r s t d = np . s q r t ( np . nansum ( r o t r m s l o g ∗∗2∗ wl ) / np . nansum ( wl ) ) a l l s t d = np . h s t a c k ( ( a l l s t d , r o t r m s l o g , r m s s t d l o g , u b e r s t d ) ) e l s e :
t 0 = t i m e . t i m e ( )
#run RATRAN f o r a s i n g l e m o l e c u l a r d a t a f i l e o s . system ( ’ amc␣amcw33a . i n p ’ )
o s . system ( ’ sky ␣ skyw33a . i n p ’ ) t 1 = t i m e . t i m e () − t 0
#c a l l t h e c o d e t h a t r e a d s t h e s k y f i t s f i l e s
a l l s t d , r o t r m s l o g , r m s s t d l o g , a l l E , allWf , allWe , allWm , f i g s = \ s f r . readshow ( l i n e s , t r a n s , s p e c i e s , ’ ’ )
wl = 1/ allWe
a l l s t d = np . h s t a c k ( ( a l l s t d , r o t r m s l o g , r m s s t d l o g ,
np . s q r t ( np . nansum ( r o t r m s l o g ∗∗2∗ wl ) / np . nansum ( wl ) ) ) ) ki nr ms = np . nansum ( a l l s t d [ np . a r a n g e ( 0 , 2 ∗ len ( Eup ) , 2 ) ] ∗ wl ) / np . nansum ( wl )
#s a v e f o r a l l m o d e l s how l a r g e t h e d i f f e r e n c e w i t h t h e o b s e r v a t i o n s i s c o n v d a t = np . l o a d t x t ( ’ path / t o / c o n v d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ )
np . s a v e t x t ( ’ path / t o / c o n v d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ , np . v s t a c k ( ( convdat , np . h s t a c k ( ( constmod [ 0 ] , constmod [ − 1 ] , a l l s t d [ − 1 ] , t1 ,
d o p p l e r [ 0 ] , d o p p l e r [ − 1 ] , v [ 0 ] , v [ − 1 ] , ki nr ms ) ) ) ) ) print a l l s t d [ − 1 ]
28
p l t . f i g u r e ( )
p l t . e r r o r b a r ( a l l E , np . l o g 1 0 ( a l l W f ) , fmt= ’ ∗ ’ , ms=10 , y e r r=allWe ) p l t . p l o t ( a l l E , np . l o g 1 0 ( allWm ) , ’ o ’ )
p l t . x l a b e l ( ’ Eup␣ (K) ’ )
p l t . y l a b e l ( ’ l o g 1 0 ␣ Flux ␣ (K␣km␣ s ^−1) ’ ) p l t . g r i d ( )
p l t . s a v e f i g ( ’ f l u x . png ’ , bbox_inches= ’ t i g h t ’ , pad_inches =0) p l t . c l o s e ( )
d i f f d a t = np . l o a d t x t ( ’ path / t o / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ ) d i f f d a t 2 = np . l o a d t x t ( ’ path / t o / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e 2 . d a t ’ ) d i f f k i n = np . mean ( d i f f d a t 2 [ np . a r a n g e ( 0 , 2 ∗ len ( Eup ) , 2 ) ] )
#s a v e e x t r a d a t a f o r t h e b e s t f i t t i n g model s o f a r i f d i f f d a t [ − 1 ] >= a l l s t d [ − 1 ] :
np . s a v e t x t ( ’ path / t o / h2 ’+s p e c i e s+ ’ p r o f i l e s h a p e . d at ’ , constmod ) np . s a v e t x t ( ’ path / t o / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ , a l l s t d )
p l t . f i g u r e ( )
p l t . e r r o r b a r ( a l l E , np . l o g 1 0 ( a l l W f ) , fmt= ’ ∗ ’ , ms=10 , y e r r=allWe ) p l t . p l o t ( a l l E , np . l o g 1 0 ( allWm ) , ’ o ’ )
p l t . x l a b e l ( ’ Eup␣ (K) ’ )
p l t . y l a b e l ( ’ l o g 1 0 ␣ Flux ␣ (K␣km␣ s ^−1) ’ ) p l t . g r i d ( )
p l t . s a v e f i g ( ’ path / t o / f l u x ’+s p e c i e s+ ’ p r o f i l e s h a p e . png ’ , bbox_inches= ’ t i g h t ’ , pad_inches =0)
p l t . c l o s e ( )
f o r i in np . a r a n g e ( 0 , len ( f i g s ) , 2 ) :
f i g s [ i ] . s a v e f i g ( ’ . . / ’+s p e c i e s+s t r ( i )+ ’ p r o f i l e s h a p e . png ’ , bbox_inches= ’ t i g h t ’ , pad_inches =0)
i f d i f f k i n >= kin rm s or np . i s n a n ( d i f f k i n ) :
np . s a v e t x t ( ’ path / t o / d i f f d a t ’+s p e c i e s+ ’ p r o f i l e s h a p e 2 . d a t ’ , a l l s t d ) np . s a v e t x t ( ’ path / t o /db ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ , d o p p l e r ) np . s a v e t x t ( ’ path / t o /v ’+s p e c i e s+ ’ p r o f i l e s h a p e . d a t ’ , v )
f o r i in np . a r a n g e ( 1 , len ( f i g s ) , 2 ) :
f i g s [ i ] . s a v e f i g ( ’ . . / ’+s p e c i e s+s t r ( i )+ ’ p r o f i l e s h a p e . png ’ , bbox_inches= ’ t i g h t ’ , pad_inches =0)
It will call a script to modify the physical structure, denshower.py:
# −∗− c o d i n g : u t f −8 −∗−
"""
C r e a t e d on Mon Mar 19 1 1 : 1 8 : 1 4 2018
29
"""
from __future__ import d i v i s i o n import numpy a s np
def d e t e r d e n s ( Eup , constmod , d o p p l e r , v , data , metadata , columns ) : f o r i in xrange ( len ( columns ) ) :
i f columns [ i ] == ’ nh ’ : C = i
e l i f columns [ i ] == ’nm ’ : M = i
e l i f columns [ i ] == ’ db ’ : B = i
e l i f columns [ i ] == ’ v r ’ : V = i
d a t a [M] = d a t a [ C]∗10∗∗ − constmod d a t a [ B ] = d o p p l e r
d a t a [V] = v
f o r i in np . t r a n s p o s e ( d a t a ) . a s t y p e ( s t r ) :
metadata = np . h s t a c k ( ( metadata , ’ ␣ ’ . j o i n ( i ) ) )
np . s a v e t x t ( ’ path / t o / m o d i f i e d /RATRAN/ model . mdl ’ , np . t r a n s p o s e ( [ metadata ] ) , fmt= ’%s ’ )
and a script which reads the output and determines the fit with the observations, skyfitsread.py:
# −∗− c o d i n g : u t f −8 −∗−
"""
C r e a t e d on F r i Apr 20 1 3 : 5 4 : 1 0 2018
@author : H y l k e D. Hoogland
"""
from __future__ import d i v i s i o n import numpy a s np
import s c i p y . o p t i m i z e a s s o import a s t r o p y . i o . f i t s a s p y f i t s
import a s t r o p y . c o n v o l u t i o n a s a s t r o c o n import m a t p l o t l i b . p y p l o t a s p l t
import o s
30
