The prescription of ​Bourrier et al is to calculate (at each pixel) the average optical depth \tau_i for range of velocity bins with spacing \Delta v (or frequency bins with size \Delta \nu, so that

\sum \tau_i \Delta \nu = \tau

So we have

\tau_i \Delta \nu = \tau_i \frac{\Delta v}{\lambda_0} = \displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \tau (v) dv

and

\tau (v) = \sigma_{v0} \int N(v') \phi ( v-v' ) dv'

Now to get the line profile in terms of velocities, we can start with the line profile in terms of frequencies:

\phi ( \nu - \nu' ) = \frac{\Gamma}{4 \pi^2 (\nu - \nu')^2 + \left ( \Gamma / 2 \right )^2}

and realizing that this is a distribution, we can use

\phi ( v - v') dv' = \phi ( \nu - \nu') d \nu

where

\nu = \nu_0 + \nu_0 \frac{v}{c}

so

\nu - \nu' = \frac{v - v'}{\lambda_0}

and

d \nu = \frac{dv}{\lambda_0}

to get

\phi ( v-v' ) =\frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2}

Putting this all together we have…

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int N(v') \frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv' dv

Now if we bin the column densities into velocity bins, we have

N_j = \displaystyle \int_{v_j-\frac{1}{2}\Delta v}^{v_j + \frac{1}{2}\Delta v} N (v') dv'

and

\int N(v') \frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv' = \sum N_j \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2}

and we have

\tau_i = \frac{\lambda_0}{\Delta v} \sigma_{v0} \sum N_j \displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv

or

\tau_i = \frac{\lambda_0}{\Delta v} \sigma_{v0} \sum N_j C_{ij}

where

C_{ij}=\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv

Now \Delta v is of order 20 km/s and \Gamma \lambda_0 is 76 m/s… so if i == j, then the integral may as well go to infinity and we have

C_{ii} = 1

and for i \ne j (wave arms in air for a bit)

C_{ij} = \frac{\Gamma \lambda_0}{4 \pi^2 (v_i-v_j)^2} \Delta v

so we have

\tau_i = \frac{\sigma_{v0} N_i \lambda_0}{\Delta v} + \displaystyle \sum_{i \ne j} N_j \sigma_{v0} \frac{\Gamma \lambda_0^2}{4 \pi^2 (v_i-v_j)^2}

Or at least that's how I get agreement with their equation 11…

So were I to implement this approach, I would create projections (integrations along LOS) of \kappa_i where

\kappa_i(n,v) = \delta_{ij}\frac{\sigma_{v0} n X \lambda_0}{\Delta v} + \left ( 1 - \delta_{ij} \right ) n \sigma_{v0} \frac{\Gamma \lambda_0^2}{4 \pi^2 (v_i-v_j)^2}

where j corresponds to the bin that the LOS cell velocity v lies within and X is the ionization fraction computing using either the ​Saha equation - or balancing ​radiative recombination with ​direct ionization… They give very different answers for the ionization fraction at the planet wind temp.

Remaining questions

1. How to chose the ionization fraction? - Need to balance photo ionization with radiative recombination
2. We should use thermal broadening (which dominates over natural broadening)..
3. Binning particle velocities before convolution seems numerically expedient - but binning errors are orders larger than convolution corrections…?
4. Why not just numerically approximate integral of line profile - and calculate contributions to each absorption bin for each cell.

Thermal broadening

It should be more accurate to calculate the LOS velocity and the thermally broadened line profile, and then calculate the contribution of that profile to each bin in frequency space. The thermal broadening profile is just a Gaussian, so using the error function we can calculate the contribution to each bin…

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int N(v') \phi(v-v') dv' dv

Now without thermal broadening, we have

N(v') = \sum n_j \delta(v'-v_j) \Delta z

however, with thermal broadening, we have

N(v') = \sum n_j \frac{1}{\sqrt{\pi} v_{th}} e^{-(v'-v_j)^2/v_{th}^2} \Delta z

where

v_{th}=\sqrt{\frac{2k_B T}{m}}

which gives us

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \sum n_j \Delta z \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int \frac{1}{\sqrt{\pi} v_{th}} e^{-(v'-v_j)^2/v_{th}^2} \phi(v-v') dv' dv

Now the integral over v' is a convolution with a lorentz and gaussian profile - which gives a Voight profile… however, the thermal broadening is much larger than the natural broadening, so we can approximate the natural broadening as a delta function…

\phi(v-v')=\delta(v-v')

which gives

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \sum n_j \Delta z \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{1}{\sqrt{\pi} v_{th}} e^{-(v-v_j)^2/v_{th}^2} dv

or upon integrating the exponential, we get

\tau_i = \frac{\sigma_{v0} \lambda_0}{2\Delta v}\displaystyle \sum n_j \Delta z \left [ \mbox{erf}\left ( \frac{v_j-v_i+\frac{1}{2}\Delta v)}{v_{th}} \right ) - \mbox{erf} \left ( \frac{v_j - v_i - \frac{1}{2}\Delta v)}{v_{th}} \right ) \right ]

So for each frequency bin, we can create a LOS integration of a corresponding integrated opacity

\kappa_i=\frac{n \sigma_{v0} \lambda_0}{2\Delta v} \left [ \mbox{erf}\left ( \frac{v_j-v_i+\frac{1}{2}\Delta v)}{v_{th}} \right ) - \mbox{erf} \left ( \frac{v_j - v_i - \frac{1}{2}\Delta v)}{v_{th}} \right ) \right ]

where n is the number density of neutral Hydrogen

Photoionization

Also need to update code to include photoionization rates as well as direct ionization and radiative recombination…

Photoionization rate can be calculated from

A=n_{HI} \displaystyle \int_{\nu_L}^\infty \frac{L_\nu \sigma_H(\nu)}{4 \pi r^2 h \nu} d\nu

We can approximate

\sigma_H(v) = \sigma_L \left ( \nu/ \nu_L \right )^{-3}

where

\sigma_L = 6.3 \times 10^{-18} \mbox{cm}^2

and

\nu_L=\frac{13.6 \mbox{eV} }{h}

and we can calculate the luminosity from Planck's law integrated over the area and solid angle…

L_\nu=\displaystyle \int_S \int_0^{2\pi} \int_0^{\pi/2} I(T, \nu) \cos \theta \sin \theta d\theta d\phi dA

L_\nu=4 \pi R^2 \left (\pi \right ) \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k T}}-1}

L_\nu=4 \pi R^2 \left (\pi \right ) \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{\nu}{\nu_S}}-1}

L_\nu=\frac{8 \pi^2 h R^2}{c^2} \frac{\nu^3}{e^{\frac{\nu}{\nu_S}}-1}

where

\nu_S = \frac{k T}{h}

and

\lambda_S = \frac{c}{\nu_S}

Putting this all together, we get

A=2 \pi n_{HI} \nu_S \frac{\nu_L^3}{\nu_S^3} \frac{R^2}{r^2} \frac{\sigma_L}{\lambda_S^2 } \displaystyle \int_{\frac{v_L}{v_S}}^\infty \frac{1}{x} \frac{1}{e^x-1} dx

Now assuming a temperature of 5778 K, this gives

A = \frac{\alpha n_{HI}}{r^2}

where

\alpha = 3.6774 \times 10^{17} \mbox{cm}^2/\mbox{s}

and a photoionization time scale (at the .047 AU) of .0435 yr

Here is the same plot as before, but included are the equilibrium ionization fraction as a function of temperature for material at the orbital separation and with densities of 1e7 and 1e9 particles/cm³

Now in Bourrier et al, they leave the photoionization rate a free parameter and adjust it between .5, 1, and 5 times the solar value.

They quote the ionizing flux at the solar minimum from

​Bzowski et al 2008

which cite

​http://www.aanda.org/articles/aa/pdf/2008/43/aa8810-07.pdf

which have a photoionization rate at 1 AU of 1e-7/s or 115 days. At .047 AU this would be 6 hours which is close to the 8 hours shown in Bourrier et al.

This would imply a photoionizing rate of 1.7165405 \times 10^{19} cm^2/s / r^2

currently I am using 3.5808d19