Multi-line Bayesian Modeling Pipeline

This directory includes the scripts and parameters used in [1] for analyzing and cross comparing the galaxy centers of NGC3351, NGC3627, and NGC4321. Most of the scripts are adaped from those in ../multiline-ngc3351, as the modeling approach is similar to that in [2].

This document lists the general workflow and explains the usage of each script. For detailed instructions on setting up the model and environment, please check the README file in ../multiline-ngc3351. The information (steps 1-4 in particular) here also applies to the scripts in ../multiline-ngc3351.

1. Construct RADEX model grids

• radex_pipeline.py: return model grids of line integrated intensity (and ratio) for all six lines (more details here)
• flux_model_6d.py: convert all intensity model grids to 5D and then add the 6th dimension (beam-filling factor)
• radex_par_extract.py: extract a specified parameter output (other than intensity; e.g. tau, Tex) from RADEX and return the corresponding 6D model grids

2. Test the modeling on selected pixels

• bestfit_onepix.py: calculate the best-fit (minimum chi^2) solution for a given pixel

  • generate and save a plot that shows how the solution is constrained by all the observed line intensities (e.g., Figure 6 in [1])
  • save a grid with chi^2 values that will be loaded into corner.py
  • print the minimum chi^2 value and best-fit solution

• corner_plot.py: calculate the marginalized probability distribution functions (PDFs) for all the modeled parameters for a given pixel (e.g., Figure 7 in [1])

  • generate and save a corner plot showing 1D and 2D likelihood distributions
  • optional settings: overplot vertical lines showing bestfit/1DMax solutions

• For indirectly modeled parameters such as CO-to-H2 factor (alpha_CO) and optical depth, the 1D PDFs at a designated pixel can be generated using Xco_marginalize.py, tau_marginalize.py, and sqrt_n_over_T_marginalize.py.

3. Handle observational data and masks (basic)

• flux_recover.py: calculate the 12m-only/12m+7m+TP flux recovery ratio and returns a mask that excluded regions with lower ratio than the desired threshold
• create all the required masks for S/N cutoff (e.g., > xxx sigma for xxx lines)
• mom0_subplots.py: show moment 0 maps of all the lines as in Figures 1 & 2 of [1]
• ratio_plot.py: generate line ratio maps that are properly masked
• radius_projection.py: define regions based on projected galactocentric radii (while also considering gas morphology and line ratio variations), and create a mask for each defined region
• region.py: generate the figure showing each defined region as in Figure 3 of [1]
• region_stat.py: return all kinds of statistics within each defined region

4. Run pixel-by-pixel modeling over the whole maps

• model_maps_nocov.py: return maps of the solutions for all directly modeled paramters as numpy arrays

  • Various options of output map type: 'bestfit', '1dmax', 'median_interp', 'median', and 'mean'.
  • 'median_interp' extracts the interpolated value of each parameter that corresponds to the PDF median, while 'median' returns the original value of the discretely sampled input parameter space.

• model_maps_nocov_testall.py: same as the above, but used for testing the results by inputting various combination of lines (i.e., Section 5.5 in [1])

• model_maps_plot.py: show properly masked maps of all directly modeled paramters returned from model_maps_nocov.py or model_maps_nocov_testall.py (e.g., Figures 8 & 9 in [1])

• For indirectly modeled parameters, the pixel-by-pixel maps are generated by e.g. model_maps_nocov_alpha.py, model_maps_nocov_tau.py, model_maps_nocov_sqrt_n_over_T.py.

• alpha2var_plot.py: show scatter plots of both median and 1DMax alpha_CO vs. projected galactocentric radius, color-coded by the defined regions (e.g., Figures 10b & 11b in [1])

5. Cross comparison of various quantities in all three galaxies

• cross_compare_plot.py: can generate various plots that incorporate all three galaxies. All figures in Section 5 of [1] are generated using this script. There are three options of figure types:

  • 'scatter': create scatter plots as in Figures 13a, 14b, 15, 16, 19 of [1]. 'colorcode' can be set to True to colorcode the data points by an additional variable, or set to False to simply colorcode the points by galaxies.
  • 'run_med': plot the x-axis binned medians as solid lines and the 25-75th percentile range as shaded areas, as in Figures 12, 14a, 17, 18 of [1]. 'compare_alpha' is set to True to generate Figure 18, which compares our results with past simulation studies.
  • bin_2d: create a scatter plot where the data points are colorcoded by another parameter that is 2D-binned (e.g., Figure 13b in [1]). The 2D binning is to bring out more clearly the overall trend of color gradient.

• alpha_fitting.py: perform power-law fit to a single or multi-variate function, evaluate the uncertainty of fitting via 1000 times refitting, and calculate the 25/75 percentiles of data scatter with respect to the fit. This script is used to determine Equations 4-6 in [1], which describes alpha_CO as functions of optical depth, temperature, and the observed line width or line ratio.

Citation

If you use or reference these script(s) in your work, please cite the following papers:

[1] Teng et al., "The Physical Drivers and Observational Tracers of CO-to-H2 Conversion Factor Variations in Nearby Barred Galaxy Centers", 2023, The Astrophysical Journal (ApJ), 950, 119. [paper]

[2] Teng et al., "Molecular Gas Properties and CO-to-H2 Conversion Factors in the Central Kiloparsec of NGC 3351", 2022, The Astrophysical Journal (ApJ), 925, 72. [paper]

More information on this project can be found in this page.