๐กโ๏ธ๐ Photon Noiseยถ
Light travels as photons. When we talk about the luminosity of a star or the brightness of the Sun in the sky, we're talking about so many photons that it makes sense to treat light as a continuous stream. However, when light spreads out of interstellar distances or when we're imaging something that's intrinsically faint, we need to start paying attention to individual photons. When we're measuring light, even if we neglect all other sources of noise, the randomness of the times with which these photons arrive at our telescope create an inescapable source of uncertainty called "photon noise" or "Poisson noise."
This page presents a pedagogical simulation of a cartoon telescope, and uses it to explore the concept of photon noise.
from astr3510 import catch_photons_in_bucket
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
Telescopes as Light Bucketsยถ
Let's imagine we focus our telescope on a star that has a particular brightness, meaning that there is a particular rate at which photons from this star rain down on our telescope, in units like $\mathrm{photons/s/m^2}$. We can imagine our telescope as a light bucket, that we hold out in this rain of photons. The larger our bucket (telescope aperture) or the longer we let it collect rain (exposure time), the more raindrops will fall into the bucket (photons).
If we know the collecting area of the telescope, the exposure time, and the brightness of the source we're observing, we can calculate an exact average number of photons we should expect to see in an exposure. In reality, we won't alway detect exactly that expected number, and not just because this expectation value might not be a non-integer number of photons. The discrete nature of light means that the actual photons we detect in a given observation will be drawn from a Poisson probability distribution centered at the expected number of photons.
Let's experiment a few times with some default settings:
N = catch_photons_in_bucket()
N = catch_photons_in_bucket()
N = catch_photons_in_bucket()
The left panel visualizes positions where individual photons might have entered the telescope; the upper right shows the numerical details of how many photons we expect and why; and the lower right shows how the actual number of photons detected compares to the Poisson probability distribution for the given expectation value.
In three different exposures, even with the same aperture, exposure time, and object brightness, we detected a different number of photons due the randomness of the times at which individual photons arrive at the telescope.
$\sqrt{N}$ and $1/\sqrt{N}$ยถ
Let's play with the input parameters of our observation a little by setting some keywords in the catch_photons_in_bucket function. We'll change the incoming photon flux (rate), the telescope mirror diameter (diameter), and the exposure time (time). We'll tweak the values to get us to an expectation value of 100 photons.
N = catch_photons_in_bucket(rate=1*u.photon/u.s/u.m**2,
diameter=2/np.sqrt(np.pi)*u.m,
time=100*u.s)
Now, let's compare what we see how things changes for different expectation values for the number of photons. We'll do this by changing the exposure time, collecting photons for a shorter or longer duration.
for t in [1, 10, 100, 1000, 10000]*u.s:
N = catch_photons_in_bucket(rate=10*u.photon/u.s/u.m**2,
diameter=2/np.sqrt(np.pi)*u.m,
time=t)
If we expect $N$ photons, the standard deviation of the Poisson distribution will be about $\sqrt{N}$. Looking closely at the numbers, we'll see that the absolute values of the numbers tend to be farther away from the expected values for larger $N$. However, it's important to consider what these means for the fractional uncertainty:
- At small values, $N$ might be $10\pm\sqrt{10}$, which means there's a very large fractional uncertainty on how many photons we'll actually detect, as $\sqrt{10}/10 = 31.6\%$.
- At large values, with $N$ being like $10^4 \pm \sqrt{10^4}$, the fractional uncertainty will be much lower, as $\sqrt{10^4}/10^4 = 1\%$ The uncertainty on the number of photons $N$ will increase as $\sqrt{N}$, but the fractional uncertainty on the number of photons we're detecting (which is what we'll ultimately translate into a statement of "how bright is this object") will go down as $\sqrt{N}/N = 1/\sqrt{N}$.
We can test this behavior by calling our photon-catching functions for many exposure times and seeing how the numbers change with exposure time. We'll do this with visualizations turned off, to avoid clutter.
# try a bunch of different exposure times
expectations = []
actuals = []
times = np.logspace(0, 2, 200)
for t in times:
# store the expected number
brightness = 1*u.photon/u.s/u.m**2
D = 2/np.sqrt(np.pi)*u.m
expected = np.pi*(D/2)**2*brightness*t*u.s
expectations.append(expected)
# store to observe number
N = catch_photons_in_bucket(rate=brightness,
diameter=D,
time=t*u.s,
visualize=False)
actuals.append(N)
N_expected = u.Quantity(expectations).value
N_actual = u.Quantity(actuals).value
# make a grid to store the results plot
fi, ax = plt.subplots(2, 2, figsize=(6,4), dpi=300, sharex='col', constrained_layout=True)
# plot the actual numbers of photons linearly
plt.sca(ax[0, 0])
plt.scatter(times, N_actual, s=10)
plt.fill_between(times,
N_expected-np.sqrt(N_expected),
N_expected+np.sqrt(N_expected),
zorder=-1)
plt.ylabel('$N_{actual}$')
# plot the actual numbers of photons logarithmically
plt.sca(ax[0, 1])
plt.scatter(times, N_actual, s=10)
plt.fill_between(times,
N_expected-np.sqrt(N_expected),
N_expected+np.sqrt(N_expected),
zorder=-1)
plt.yscale('log')
plt.xscale('log')
plt.ylabel('$N_{actual}$')
# plot the actual divided by the expected, linearly
plt.sca(ax[1,0])
plt.scatter(times, N_actual/N_expected, s=10)
plt.fill_between(times,
1-1/np.sqrt(N_expected),
1+1/np.sqrt(N_expected),
zorder=-1)
plt.ylabel('$N_{actual}/N_{expected}$')
plt.ylim(0, 2)
# plot the actual divided by the expected, logarithmically
plt.sca(ax[1,1])
plt.scatter(times, N_actual/N_expected, s=10)
plt.fill_between(times,
1-1/np.sqrt(N_expected),
1+1/np.sqrt(N_expected),
zorder=-1)
plt.ylabel('$N_{actual}/N_{expected}$')
plt.ylim(0, 2)
plt.xscale('log')
fi.supxlabel('Time (s)');
With our little numerical experiments, we can see that the Poisson distribution gets wider as $\sqrt{N}$ for increasing expectation value $N$, but that the fraction difference from the true underlying expected value goes down as $1\sqrt{N}$.
Learn more! ๐งโ๐ซ๐โ๏ธยถ
Play around with the catch_photons_in_bucket function to develop your understanding of how photon-counting statistics behave in different limits. Have fun!