### Week 3 Research Update:

Today we finally started implementing implementing the Lyman alpha methodology after successfully creating our own re-ionisation model and reproducing the results shown in the original masters paper. This methodology will be introduced by implementing the following two equations:

Where:

Note, is the number of LyC photons produced per second predicted by observables, not the number of photons produced per second.

Our first step in introducing the Lyman alpha methodology was to write a function that could convert the escape fraction to a Lyman continuum escape fraction.

Originally we tried to use the equation shown here (p. 32 right-hand side figure):

The plot from last week (repeated below) showed that for the LAEs we have data for, they have an average of 140 angstroms. Using the equation above, we ended up with a Lyman continuum escape fraction of 96%. This seems extremely high as LAEs only have a Lyman alpha escape fraction between 30% – 50%. After speaking to David, we were pointed towards A. Verhamme et al.‘s paper. Table one in this paper provided us with data for several galaxies, this data contained their equivalent widths and LyC escape fractions. This plot is also shown below.

The blue data point in the right-hand plot is the average equivalent width shown by the near-constant line on the left-hand plot. Using this data, it can be seen that for our average equivalent width LAEs would have a LyC escape fraction of approximately 10%.

The other term in the equation that we need to determine via data is (data will be added to our GitHub repository soon). We hardcoded the data from the fits file provided and the used pylab’s polyfit module to generate as a function of redshift, as shown below.

```
def L_Lya(z): # Data from Table_C2_Sobral18_SSC4K_compilation.fits
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])
p2 = pylab.polyfit(x, y, 2.0)
p = pylab.poly1d(p2)
return p(z)*10**40
```

Once we had implemented the methodology, we ran our model and created the following plot:

This is very different to the original plots we produced in the previous weeks. Firstly we are now (incorrectly) ionising the universe as soon as it’s born. Secondly, the epoch of re-ionisation is finishing around . As there are many dependencies throughout our model, it is quite hard to determine what is causing the error by looking at this plot alone. These dependencies have been mapped out below.

Identifying all these relationships tells us which plots to produce when we run our simulation so we can visualise these dependencies and identify our errors. We are working on a plotting module which will automatically plot all these relationships when we solve the final ODE, this will hopefully be finished by the end of next week. In the mean time, these plots have been created in a stand-alone file. The plots showing the relationships between , and redshift are shown at the top of this entry.

The plot that stands out the most here is the luminosity – redshift plot. Our current function is returning negative luminosities for redshifts greater than approximately seven. The plots for and also show that in our model, the availability of ionising photons rapidly drops after redshift five, which is why our re-ionisation process is taking so long.

In the next session we will look into improving our function and determining why our number of ionising photons drops to zero.