Exploring Cosmological Parameters with Functions#
Learning Objectives
Define Python functions and use their outputs within other functions, which is useful for calculations dependent on previous calculations.
Understand how to apply functions over array inputs both with loops and with vectorized operations.
Characterize and identify different models for the universe with the Hubble-Lemaitre Law, a common exercise in introductory astronomy courses.
Document functions properly with type decorators and docstrings.
Background: The Hubble-Lemaitre Law#
The Hubble-Lemaitre law (often referred to as just “Hubble’s Law”) describes the relationship between the velocity of a galaxy and its distance from us:
where \(v\) is is the recession velocity of a galaxy (in km/s), \(H_0\) is the Hubble constant (in km/s/Mpc; assume 70 km/s/Mpc), and \(d\) is the distance of the galaxy (in Mpc). \(H_0\) is the present-day expansion rate of the universe. In general, we cannot measure direct distances to very distant galaxies, and so we measure the recessional velocity \(v\) and can use that to find \(d\). The fact that this correlation exists is a consequence of the fact that our universe is expanding.
Exercise 1
Define a function
hubble_distanceto calculate the distance to a galaxy based on this formula. The function should take one required argument, the recessional velocity in km/s, and should return the distance to the galaxy in Megaparsecs (Mpc).Add the Hubble constant as an optional argument that takes a default value of 70 km/s/Mpc.
Add a print statement of the form: Given a recessional velocity of [] and a Hubble constant of [], a galaxy’s distance would be [] (fill in the units here, and use f-strings to handle the values).
Verify your code by checking that you recover 20 Mpc as the distance for a galaxy moving at 1400 km/s, assuming the default \(H_0\) value. (for reference, 20 Mpc is in the local universe, and space telescopes like JWST can resolve individual stars in galaxies at this distance).
Now, create an evenly-spaced array (hint: np.arange or np.linspace is your friend here!) of 20 recessional velocities called recess_vels in the range [1000,100000]. Using a loop, evaluate your function over each element. Store the results in a list.
Without using a loop, evaluate your function over recess_vels.
The output from 6 just prints everything in a rather sloppy way. Go back and add another optional argument print_res that is a boolean defaulting to False. Using an if statement, make it so that the results statement only prints if this variable is set to True.
Add type decorators to the inputs.
Lastly, add a docstring detailing the purpose of the function, its arguments, and the return type.
Questions 5-6 above get at a core concept in Python: vectorization. In Python, it will essentially always be faster to carry out an operation over a vector (array) simultaneously (part 5) as opposed to using a loop (part 4). This might not matter for small datasets but will become increasingly relevant as your dataset grows.
Now we can make our approach even more sophisticated! We will now build subsequent functions that expand upon the one above. This incremental function-building is a core component of research, since you rarely write a single function that does everything you want.
Background: The Cosmic Critical Density#
The cosmic critical density is the average density of matter and energy needed for the universe to be flat. If the actual density is higher, the universe is “closed” and may eventually collapse; if the density is lower, the universe is “open” and will expand forever. The density can be calculated from the Hubble Constant via \(\rho_c = \frac{3H_0^2}{8\pi G}\)
Background: The Age of the Universe#
Because the universe has been continuously expanding since the origin of the universe, the present-day expansion rate can be used to approximate the age of the universe. Put explicitly, the age of the universe can be estimated roughly as \( t = \frac{1}{H_0}\), but you need to be careful about units: you need to first convert \(H_0\) to units of inverse seconds for this to work.
Exercise 2
In this exercise, you will take galaxy distance and recession velocity data and use it to derive the density, age, and shape of the universe.
Write a function get_hubble_constant which takes in an array of galaxy distances in Mpc and an array of recessional velocities in km/s. Compute \(H_0\) for each element and then return a single average. Don’t use loops for anything in this exercise!
Write a function get_critical_density that computes the critical density according to the above formula, with only a single argument of \(H_0\). Determine and print the “shape” of the universe (open, closed, flat) by comparing against a model where the Hubble constant is fixed to \(H_0 = 73\) km/s/Mpc (which is realistic, and we think we live in a flat universe).
Write a function get_universe_age that takes in \(H_0\) in km/s/Mpc and returns the age of the universe in years. Note that you need to include two unit conversions: first, you will need to convert km/s/Mpc to just 1/s. After inverting this, you can then convert the time in seconds to a time in years.
Write a function cosmology_solver that wraps all of the three functions, i.e., it calculates the critical density, shape, and age of the universe. Specifically, it should take the same inputs as get_hubble_constant. Within cosmology_solver, you should run each of the functions you created above such that you can figure out the properties of the univese for a given dataset of galaxy distances/recession velocities. Print the critical density.
Test your complete cosmology solver with the following data: distances: [45,100,200,400] Mpc and velocities [3400, 7600, 15200, 31800] km/s. What kind of univese do you end up with? Is it older or younger than our \(H_0 = 70\)ish universe?