ReHILAE: Is Re-ionisation of Hydrogen due to Lyman-alpha emitters (LAEs)?

Week 5 Research Update:

Following on from last week’s session, we started today’s session with implementing the redshift dependent coefficient in the LAE UV model which we discussed last week. On page 21 of David’s paper here, you can see a plot of f_{esc} \times \xi_{ion} against redshift using this data and the equation shown on the plot we were able to implement the coefficient previously discussed.

Our main focus for today’s session was to focus on replacing all of our polyfit functions with SciPy curves. Not only would using SciPy return the best fit parameters for the fitted curves, it would also return the errors on the curves to us. The two main functions that use a fit derived from observational data in our new LyC framework are the f_{esc,LyC} and \rho_{L_{Ly\alpha}} functions. Using SciPi’s curve fit function we were able to to produce the following plots and derive the estimated functions for these parameters:

The errors on the on the initial data was estimated as 20% and 15% correspondingly. Using the equations shown on the plots we were able to replace the f_{esc,LyC} and \rho_{L_{Ly\alpha}} functions with their new equations. This reduced the run-time of our code considerably; instead of fitting a polyfit to the same data to generate a function for every single time-step, we now have predefined function to call. The code snippets below show how different our previous and new \rho_{L_{Ly\alpha}} are.

# Previous version
def P_L_Lya(z) # data from SC4K Sobral
    x = pylab.array([2.2,2.5,2.8,3.0,3.2, 3.3, 3.7, 4.1, 4.6, 4.8, 5.1, 5.3, 5.8 ])
    y = pylab.array([0.52, 0.74, 0.77, 0.88, 0.84, 0.85, 1.01, 0.87, 1.19, 1.12, 1.27, 1.08, 1.10])  # data from SC4K Sobral 
    p2 = pylab.polyfit(x, y, 2.0)
    p = pylab.poly1d(p2)
    if p(z) < 0:
        return 0
        return p(z) * 10**40
# New version
def P_L_Lya(x, P1, P2, P3):
    return 10**(P1*x**2 + P2*x +P3)

In the new version, P1 , P2 , P3 represent the coefficients shown on the log(\rho_{L_{Ly\alpha}}) plot above and x corresponds to the redshift value that would be parsed into the function.

Once we had replaced all of our polyfits we SciPy curves we moved onto the final phase of our project: running our re-ionisation model multiple times and varying our curve fit parameters for every iteration, resulting in a range of re-ionisation plots on our final graph. When selecting the parameters for our functions we use a normal distribution. On the plots above, it can be seen that our parameters are of the form a\pm b . We can then let a be our mean and b be our standard deviation and sample our parameters for each run from a normal distribution, within one standard deviation from our mean. After running this simulation with ten iterations, we produced the following plot:

When running our simulation for one thousand iterations:

Whilst both generate a line where the epoch of re-ionisation starts at z=10 and finish at z=6 , we cannot currently conclude that LAEs are all that was needed to re-ionise the universe. From the plot with one thousand iterations, it can be seen that we have lines starting across a range of redshifts from z=14 and z=5. As well as lines that re-ionise the universe almost immediately after starting and some that do not even finish before reaching z=0. This means we need to try and reduce our errors; reducing these errors will reduce the standard deviation in our normal distributions and therefore reduce the variation in our parameters for each iteration.