=============== Nebulous

What is Nebulous

Nebulous is a reversible jump Markov Chain Monte Carlo that takes Cloudy (Ferland et al, 2013) outputs and fits line ratios to find the optimal :math:\Phi and :math:n_H parameter space the broad line region occupies under a similar assumption of the Locally Optimally Emitting Cloud (LOC) model (Baldwin et al., 1995). Cloudy is a well-established and extensive spectral synthesis code designed to simulate realistic physical condition inside gas clouds. It predicts the thermal, photoionisation, and chemical structure of a non-equilibrium cloud illuminated by an external source of radiation and predicts its resulting spectrum.

The programs first calculates line EW values (with respect to continuum flux at 1215 Å) given a line intensity grid from Cloudy and :math:\Phi and :math:n_H parameter space, described by a shape (that is fit using the reversible jump MCMC). Then uses these EW values to create line ratios. The shape can have :math:n>3 vertices and occupy any region in the allowed parameter space (example shown in Figure 1).The RJMCMC fit boths the number of vertices, the position of the vertices in :math:\Phi and :math:n_H space, and the orientation of the BLR.

.. figure:: docs/polygon_image_for_docs.jpg :scale: 60% :align: center

Fig. 1 : Example of a shape and the grid points that fall within shape and that will be used in LOC-like summation to calculate line ratios.

In the LOC approach, the observed emission line spectrum is the sum of emission line contributions from a weighted distribution of 'clouds' with a range of gas densities at a range of radii. The resulting line strength is given as,

.. math:: L_{line} \propto \int_{r_{min}}^{r_{max}} \int_{n_{H,min}}^{n_{H,max}} W_{1215}(r,n_H) f(r) g(n_H) dn_H dr,

where :math:f(r) = r^{\Gamma}, :math:g(n_H) = n_H^{\beta} are the weighting functions or 'covering fraction' of the various clouds (can be thought of as the an number density of clouds at radius :math:r and density :math:n_H), :math:r_{ min} and :math:r_{ max} are the minimum and maximum radii of the BLR, and :math:n_{ H,min} and :math:n_{ H,max} are the minimum and maximum cloud densities considered. :math:W_{1215} is the equivalent width of the line referenced to the incident continuum at 1215Å. This expression can be rewritten in term of :math:\log \Phi_H and :math:\log n_H.

.. math:: L_{line} \propto \int_{\log \Phi_{H,min}}^{\log \Phi_{vmax}} \int_{\log n_{H,min}}^{\log n_{H,max}} W_{1215}(\log \Phi,\log n_H) 10^{(\beta+1)\log n_H-0.5(\Gamma+1)\log \Phi} d \log n_H d \log \Phi.

and the :math:EW_{line} is given as

.. math:: EW_{line} \propto \frac{ \int_{\log \Phi_{H,min}}^{\log \Phi_{vmax}} \int_{\log n_{H,min}}^{\log n_{H,max}} W_{1215}(\log \Phi,\log n_H) 10^{(\beta+1)\log n_H-0.5(\Gamma+1)\log \Phi} d \log n_H d \log \Phi }{ \int_{\log \Phi_{H,min}}^{\log \Phi_{vmax}} \int_{\log n_{ H,min}}^{\log n_{ H,max}} 10^{(\beta+1)\log n_H-0.5(\Gamma+1)\log \Phi} d \log n_H d \log \Phi }.

In this analysis, instead of integrating between :math:\log \Phi_{H,min} & :math:\log \Phi_{H,max}, and :math:\log n_{H,min} & :math:\log n_{H,max}, we sum over the contributions within the prescribed shape, that is,

.. math:: EW_{line} \propto \frac{\sum_i^N W_{1215}(\log \Phi_i,\log n_{H_i}) 10^{(\beta+1)\log n_{H_i}-0.5(\Gamma+1)\log \Phi_i} \Delta \log n_H \Delta \log \Phi}{\sum_i^N 10^{(\beta+1)\log n_{H_i}-0.5(\Gamma+1)\log \Phi_i} \Delta \log n_H \Delta \log \Phi}.

where i is the ith grid point enclosed within the shape and N is the number of grid points enclosed by shape (See Figure 1 for illustration).

It is common practice to adopt the a gas density distribution weighting function :math:\beta = -1, described in \citet{Baldwin1995}, which roughly corresponds with the gas density distribution of fragmenting BLR clouds resulting from magnetohydrodynamic instabilities. For simplicity I also set :math:\Gamma=-1. Meaning that each grid point in log space is given equal weighting. These parameters can be altered in nubulous code.

Line strength is highly dependent on the value of the integration limits. Therefore we can use the combination of observed ratios to constrain the integration parameter space.

There are two cases that we can consider when summing over the :math:\Phi_H and :math:n_H parameter space defined by the shape: first we can sum over the total emissivity, or we can assume an inclination and define an observed emissivity given the contributions from the inward and outward emissivity of the cloud, such that,

.. math:: \epsilon_{\rm obs}(r,\theta) = \frac{\epsilon_{\rm in}(r)}{2}\left(1-\cos \theta\right)+ \frac{\epsilon_{\rm out}(r)}{2}\left(1+\cos \theta\right),

where :math:\epsilon_{\rm total}(r)=\epsilon_{\rm in}(r)+\epsilon_{\rm out}(r), :math:\epsilon_{\rm in}(r) is the emissivity of the inward cloud face (i.e. towards the ionising continuum source), and :math:\epsilon_{\rm out}(r) is the outward emissivity (i.e. away from the ionising continuum source). This equation can be rewritten in terms of the anisotropy distribution, :math:F(r)=\epsilon_{\rm in}(r)/\epsilon_{\rm total}(r) (shown for all lines in Figure \ref{fig:aniso}),

.. math:: \epsilon_{\rm obs}(r,\theta) = \frac{\epsilon_{\rm total}(r)}{2}\left(1-[2F(r)-1]\cos \theta\right).

An isotropically emitting cloud will have an anisotropy distribution of :math:F(r)=0.5, while :math:F(r)=1 corresponds to a cloud that emits only towards the ionising continuum source, i.e. completely reflected. The anisotropic nature of each emitting cloud is highly dependent on the ionising parameter :math:U_Hc=\Phi_H/n_H. As shown in Figure \ref{fig:aniso}, the Helium lines are emitted fairly isotropically across the whole parameter space except for a small section that is dominated by inward emission. The Hydrogen lines are generally isotropic for ionisation conditions, :math:U_Hc>10^{11} and dominated by emission towards the ionising source below this limit. My investigations have suggested that inclination is required to achieve reasonable fits so it is fit along with the shape.

Reversible Jump MCMC (RJMCMC)

As we do not know the shape or size of the best fit :math:\Phi and :math:n_H parameter space we chose to model it by a polygon, where the position and number of vertices is allowed to vary. This makes this a problem where the number of unknowns are unknown. We adopt a Bayesian approach to solve this problem, specifically we adopt the reversible jump Markov Chain Monte Carlo (RJMCMC) algorithm proposed by Green et al (1995). We follow a very similar approach to Luo et al. (2010) who used RJMCMC to constrain the shape of a gravity anomalous body.

This inverse problem is typically non-linear and non-unique and it is possible that significantly different :math:n_H-\Phi parameter space regimes can fit the data equally well. We wish to obtain an :math:n_H-\Phi parameter space model that fits the data, and but we also want to measure the flexibility in the :math:n_H-\Phi parameter space that give reasonable results.

We fit the shape of the polygon, the inclination of the system and the covering fraction (in cases where the continuum flux near 1215Å is included in the fit).

Given a prior distribution of parameter values :math:\pi(\theta) and a likelihood of the observed data given the parameters :math:\pi(y|\theta) given by,

.. math:: \ln \pi(y|\theta) = -\frac{ 1 }{ 2 }\chi^2= -\frac{1}{2} \sum_{i=1}^{N}\frac{({\rm ratios}{\rm real[i]}-{\rm ratios }{\rm model[i]})^2}{{\sigma}_{\rm ratios[i]}^2}

the distribution of the parameters :math:\theta conditional on the observations :math:y, the posterior distribution, is determined by Bayes theorem

.. math:: \pi(\theta|y) = \frac{\pi(y|\theta)\pi(\theta)}{\int \pi(y|\theta}\pi(\theta) d\theta)

Markov chain Monte Carlo is a method that samples the posterior ... A Markov process is...balanced if the likelihood of the the transition from state :math:s to :math:s' is as likely as from :math:s' to :math:s, that is the transition is reversible. Or in more technical terms the transition rates between each pair of states :math:s and :math:s' in the state space obey

.. math:: q(s,s')\pi(s) = q(s',s)\pi(s')

where :math:q(s,s') is the Markov chain transition kernel, which is effectively the likelihood of moving to state :math:s' given the present state. In our case, we do not know how many parameters should be used to specify our model. Let :math:\left\{ \mathcal{M}_k\right\} denote a collection of candidate models. Model :math:\mathcal{M}_k has a vector :math:\theta_k of unknown parameters :math:\theta_k\in \mathds{R}^{n_k}, where the dimension :math:n_k depends on the model. In our case, these unknown parameters are the coordinates of the polygon vertices(, the inclination and the covering fraction), therefore :math:n_k = 2k(+3). Under a Bayesian framework, inference on the model and model parameters is carried out using the point posterior :math:\pi(\theta,\mathcal{M}_K|Y).

A detailed balanced is satisfied for a Markov chain if the proposed move from (:math:\theta_i,\mathcal{M}_i) to (:math:\theta_i,\mathcal{M}_i) is accepted with probability :math:\alpha = min\left\{1,\alpha_{i\rightarrow j}(\theta_i,theta_j)\right\} \citep{Green1995}, where,

.. math:: \alpha_{i\rightarrow j}(\theta_i,\theta_j) = \frac{\pi(\theta_j,\mathcal{M}j)r{j\rightarrow i}(\theta_j)q_{j\rightarrow i}(\theta_j,\theta_i)}{\pi(\theta_i,\mathcal{M}i)r{i\rightarrow j}(\theta_i)q_{i\rightarrow j}(\theta_i,\theta_j)}

where :math:r_{i\rightarrow j} is the probability that a proposed jump from model :math:\mathcal{M}_i to model :math:\mathcal{M}_j is attempted, and :math:q_{i\rightarrow j}(\theta_i,\theta_j) is the density from which the proposed parameter :