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 the script(s) in your study, 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.