The optical depth distribution of an absorber that corresponds to the \(i\) cell with neutral hydrogen density \(n_{\mathrm{HI},i}\), Doppler parameter \(b_i\) and velocity \(v_i\) can be expressed as a Gaussian given by:

\[\Phi_i(v) = n_{\mathrm{HI},i} \frac{1}{\sqrt{\pi} b_i} \mathrm{exp} \left[ -\left(\frac{v-v_i}{b_i} \right)^2 \right]\]

Note that the velocity of the absorber \(v_i\) includes the Hubble flow velocity at the position of cell \(i\) and the peculiar velocity of the gas at that position:

\[v_i = v_{\mathrm{H},i} + v_{\mathrm{LOS},i}\]

This distribution is shown in the figure below

To compute the contribution of the absorber \(i\) to the optical depth at cell \(j\) located at the velocity coordinate \(v_j\) we must integrate the Gaussian profile over the cell \(j\), this is:

\[\tau_{j,i} = \int_{v_{j-1/2}}^{v_{j+1/2}} \Phi_i(v) dv\]

That integral corresponds to the area shown in the figure below:

The definition of the error function is given by:

\[\operatorname{erf} x=\frac{1}{\sqrt{\pi}} \int_{-x}^{x} e^{-t^{2}} d t\]

Defining \(y_{j+1/2} = ( v_{j+1/2} - v_i )/b_i\), then the error function evaluated at \(y_{j+1/2}\) corresponds to the area in the figure below:

Similarly for \(y_{j-1/2} = ( v_{j-1/2} - v_i )/b_i\), then the error function evaluated at \(y_{j-1/2}\) corresponds to the area in the figure below:

Then subtracting the integrals results in the area below:

For this reason there is a factor of 2 that has to be included, this results in the contribution to the optical depth at cell \(j\) from cell \(i\) given by:

\[\tau_{j,i} = \frac{\pi e^{2} \lambda_0 f_{12}}{m_{e} c H} n_{\mathrm{HI},i} \frac{ erf(y_{j+1/2,i}) - erf(y_{j-1/2,i}) }{2}\] \[y_{j+1/2,i} = ( v_{j+1/2} - v_{\mathrm{H},i} - v_{\mathrm{LOS},i} )/b_i\] \[y_{j-1/2,i} = ( v_{j-1/2} - v_{\mathrm{H},i} - v_{\mathrm{LOS},i} )/b_i\]

Finally the optical depth at cell \(j\) is given by the summation over all \(i\) cells:

\[\tau_{j} = \frac{\pi e^{2} \lambda_0 f_{12}}{m_{e} c H} \sum_i n_{\mathrm{HI},i} \frac{ erf(y_{j+1/2,i}) - erf(y_{j-1/2,i}) }{2}\]