Partnership for Integration of Computation into Undergraduate Physics

**Exercise 1: Force Fields** A force field, like any field, associates a value with each point in space. Since force is a vector, the force field due to the gravitational wave can be represented with an arrow (to indicate both magnitude and direction) at every point in space. In the programs that follow, we will ignore the z direction and analyze force fields in the x-y plane. We will start with a program that places an arrow, representing force, at every point on a 15 x 15 grid. The point (0,0) will be in the middle of the grid. We will begin with a constant and uniform force field before modeling a gravitational wave. A physical example is gravity near the surface of the Earth. It is approximately constant in time[^footnote] and uniform on small scales. 1. Run the template provided for Exercise 1. The force used is $\vec{F} = <0,-mg,0>$, where $g = 9.8 ~m/s^2$. 2. The divergence of a vector field $\vec{A}$ is defined as $div \vec{A} = \frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y}+\frac{\partial A_z}{\partial z}$. The value of the divergence may depend on (x,y,z), in which case it varies from point to point in the field. Physically, it represents the net flux (flow) from an infinitesimal volume surrounding a point. One can get an intuitive feel for whether a field has non-zero divergence in a certain region by visually examining the flux. We will imagine a small volume centered on a certain point in the field. Roughly speaking, a field that appears to be (overall) "exiting" the volume has a positive divergence; a field with a net "entrance" into the volume has a negative divergence. If there is equal entrance and exit, there is no divergence. Before continuing, it may help to get some practice identifying positive and negative divergences. The following links give examples where the divergence is either positive at all points, or negative at all points: https://mathinsight.org/divergence_idea https://mathinsight.org/divergence_subtleties When you feel comfortable recognizing divergence, return to examining your template. * Looking at the output of this program, predict whether the divergence of the field is positive, negative, or zero. * Check your prediction with a direction calculation (on paper). * Modify the code to display a force field with positive divergence. (The force field does not have to have physical significance.) **Exercise 2: Gravitational Wave Polarization** 1. Change your code so that it visualizes the following gravitational wave: $$F_{x}= \frac{-m}{2} h_{+} x\omega^{2}\cos(\omega t),$$ $$F_{y}= \frac{m}{2} h_{+} y\omega^{2}\cos(\omega t).$$ Each grid point has a different $(x,y)$ value, and thus a different value of force. The following parameters are recommended: | Parameter | Symbol | Value | | ------------- |:-------------:| -----:| | Amplitude of gravitational wave | $h_{+}, h_{\times}$ | $ 0.2 $ | | Angular frequency of oscillation | $\omega$ | $50 ~rad/s $ | | Mass of object subject to force field | $m$ | $0.1~ kg$ | | Time | $t$ | $0$ | This will show the output at a time $\omega t =0$. 2. Let's check to make sure the force field is plausible. Suppose we examine the arrow at the upper right-hand corner. $F_x$ at this point depends on $x$, $\cos(\omega t)$, and a factor $\frac{-m}{2} h_{+} \omega^{2}$. Let's figure out whether $F_x$ is overall negative or positive by examining these in turn. First, what is (x,y) for this grid point? The middle of the plot is (0,0), and the points range from -7 to +7 in each direction, incrementing in integer steps. Therefore, the value of (x,y) at this point is (+7,+7). The factor $x$ is positive. * Is $\cos(\omega t)$ positive, negative, or zero at $\omega t =0$? * What about $\frac{-m}{2} h_{+} \omega^{2}$? * Use your answers above to determine if $F_{x}$ negative, zero, or positive. * Do the same for $F_{y}$. * Based on your answer, determine the direction of the force arrow at this point. Compare to the output of the code. If they disagree, resolve the discrepancy. Another way to check the plausibility of the output is to calculate the divergence. 3. Looking at the output of your program, predict whether the divergence is positive, negative, or zero. 4. Calculate the divergence of this force (assume the z-component is zero). Hint: for $F_x$, as far as the variable "x" is concerned, the force looks like: $ F_x = const*x$. The time-dependent cosine term has no explicit x dependence, so in calculating $\frac{\partial F_x}{\partial x}$, it is treated like a constant. 5. Finally, check that your plot of the force field resembles a plus sign. This wave is called plus-polarized. The term "polarization" refers to the direction of the force field (not the direction the wave is travelling). Here, the wave travels in the z-direction, but the polarization is in the x-y plane. Light waves work similarly - if an electromagnetic plane-wave is travelling in the z direction, the polarization (direction of the associated electric field) is in the x-y plane. These types of waves are known as "transverse." *Gravitational Waves: Cross Polarization* 6. Repeat the exercise above for the cross polarization, $$ F_{x}= \frac{-m}{2} h_{\times}y \omega^{2}\cos(\omega t),$$ $$ F_{y}= \frac{-m}{2} h_{\times}x \omega^{2}\cos(\omega t).$$ Note that the placements of $x$ and $y$ are reversed. As before, perform at least 3 plausibility checks on the output: 7. Pick a point and determine the direction of the force at that point using the above equations. Compare with the arrow direction at that point. 8. All gravitational waves in this Newtonian approximation are divergence-free. Check that the divergence is zero (visually and mathematically). 9. Check that the output resembles a cross. [^footnote]: The small changes of Earth's surface gravity with time (due to mass shifting around) are a major factor that limits the precision of ground-based gravitational wave detectors. For this reason, space-based gravitational wave detectors (eg, LISA) have been designed. **Exercise 3: Effect of Gravitational Waves on a Ring** In this exercise, we will allow the force due to the gravitational wave to evolve in time, and show the effect on a ring of test masses in the x-y plane. Each test mass will have a different value of $(x_0,y_0)$ that it will oscillate around as a result of being in the force field. In this Newtonian approximation, force determines the acceleration of the test masses through Newton's 2nd law. In the x direction, this gives (for the plus polarization): $$ F_x = m \ddot x.$$ which gives $$\ddot x= \frac{-1}{2} h_{+} x\omega^{2}\cos(\omega t). \label{ac1}$$ The motion of the test mass at location $x_0$ can be written $x(t) = x_0 + \delta x$, where the displacements $\delta x$ are small compared to $x_0$. Since $x_0$ has no time dependence, $\ddot{x_0} = 0$, and thus eq. \ref{acc} actually represents $\ddot{ \delta x}$: $$\ddot{\delta x}= \frac{-1}{2} h_{+} x\omega^{2}\cos(\omega t). \label{acc}$$ To simulate this system numerically, starting with acceleration and allowing the program to determine the subsequent motion, we need to set the initial conditions for each test mass. The point $(x_0,y_0)$ is not (necessarily) the initial position of each mass; it is the equilibrium position that the mass oscillates around. We will call the initial positions and velocities $(x_i,y_i)$ and $(v_{xi},v_{yi})$. A convenient choice for the initial conditions is to set them so that the resulting motion is purely oscillatory around the point $(x_0,y_0)$. In the calculations that follow, you will use calculus to determine what these initial values will need to be. So that we can solve this equation analytically, we will approximate $x(t) \approx x_0$ in eq. \ref{acc}, and take $\omega$ and $h_{+}$ to be constant. Thus, the equation to solve is: $$\ddot{ \delta x}= \frac{-1}{2} h_{+} x_0\omega^{2}\cos(\omega t).$$ * Integrate this equation twice to derive $\delta x$ as a function of time. Show that the solution has the form: $\delta x= \frac{h_{+}}{2}x_0\cos(\omega t)+v_{xi} t + c_2$, or equivalently: $$x(t) = x_0 + \delta x=x0+ \frac{h_{+}}{2}x_0\cos(\omega t)+v_{xi} t + c_2 \label{sol} $$ To get motion that oscillates around $x_0$, we need to set $v_{xi} = 0$ and $c_2 = 0$. But our simulation doesn't allow us to input $c_2$ directly; our inputs are the initial conditions. * Find the value of $x_i$ that corresponds to $c_2 = 0$ by setting $t = 0$ and $x(0) = x_i$ in eq.\ref{sol}. Hint: you should find that the test mass has to be initially displaced from its equilibrium value $x_0$ in order to get $c_2 = 0$. * Now, find the initial conditions for the y component of the motion similarly. Integrate the acceleration in the y direction twice, and find the value of $y_i$ that corresponds to $c_2 = 0$. As before, determine the value of $v_{yi}$ that eliminates the linear term in $y(t)$. * Simulate the effect of this gravitational wave on a ring of test masses using the template for Exercise 3. You'll need to put in the force and the initial conditions. Use $x_0$ instead of $x$, and $y_0$ instead of $y$, when you enter the force. (Otherwise, the force you are simulating will not exactly match the force you used to determine the desired initial conditions.) Plausibility checks: * Do you observe behavior consistent with a "plus" polarization? * Do the masses oscillate around $(x_0,y_0)$? If the motion is shifted, the initial position is not quite right. If they fly off to infinity, the initial velocity is not quite right. (Physically, there is nothing wrong with making the initial conditions anything you like, but for the purposes of understanding the effect of the force, it is helpful to remove the "flying off" behavior.) *Analysis* For each of the following parameters in the simulation, describe in words what effect the parameter has on the motion of the test masses. If it has no effect, explain why not. For example, try changing $R$, the radius of the ring of test masses. Does the amount of deformity of the ring change as a result? Why or why not? | Parameter | Effect on motion (if none, explain why not) | | ------------- |:----------------------------------------------------------:| | $h_{+}$ | $ $ | | $\omega$ | $ $ | | $R$ | $ $ | * Copy your code into a new file so that you don't overwrite your existing code. Modify your simulation so that it shows the cross polarization. Make sure to re-derive the initial positions, since they will be different than before. * Perform plausibility checks, as before, to confirm that your program is working correctly. **Exercise 4: Orbiting Systems and Circular Polarization** Objects that orbit each other create gravitational waves. In this section, our source of gravitational waves are two masses, $m_1$ and $m_2$, that orbit in the x-y plane, with the origin at the center of mass of the system. This could represent a system where the two masses are of comparable size, such as a binary star system, or one in which a planet orbits a star. The gravitational waves produced by this system will be the same as would arise from a single mass $\mu$ orbiting at a radius $A$. (We will assume that this orbit is circular). This effective radius $A$ is equal to the difference in the orbital radii of the physical masses ($\vec{r_i}$). Specifically, $$\mu = \frac{m_1 m_2}{m_1+m_2}\label{mu}$$ and $$A \equiv |\vec{r_1} - \vec{r_2}|$$ Let's imagine that our gravitational wave detector (which we will again model as a ring of test masses) is placed on the z axis, a far distance $r$ from the source of the gravitational waves. In this case, the resulting gravitational wave will have both plus and cross polarizations at the same time. The angular frequency of the gravitational wave turns out to have twice the angular frequency of the orbiting sources. The frequency of the source will be denoted $\omega_s$. $$h_{+}(t) = \frac{4G\mu \omega_{s}^{2} A^2}{rc^4}\cos(2\omega_s t), $$ $$ h_{\times}(t) = \frac{4G\mu \omega_{s}^{2} A^2}{rc^4}\sin(2\omega_s t). $$ Note that Newton's gravitational constant (G) appears, because these equations come from the Einstein equation. The Einstein equation relates the stress-energy tensor (information about mass/energy) to the curvature of spacetime. The constant $G$ tells spacetime how strongly to curve in response to the mass/energy present. * Explain what happens to the amplitudes $h_{+}$ and $h_{\times}$ as the detector gets farther from the source. (The distance between them is $r$). Why does this make sense? * Determine the initial conditions that will produce purely oscillatory motion for $\delta x$ and $\delta y$ (no constant or linear terms). You'll have to obtain acceleration from the force equations, and then integrate twice to find $\delta x$ and $\delta y$. As before, approximate $x$ and $y$ as $x_0$ and $y_0$. Hint: You should find that the initial velocity in the x direction must be $\dot {\delta x}(0)= \frac{G\mu \omega_{s}^{3} A^2}{rc^4} y_0$, in order to make the first constant of integration equal zero. Once you have determined the four initial conditions, you are ready to simulate this gravitational wave. * Without overwriting your previous program, create a program to model this gravitational wave. The net force will be the sum of the contributions from the two polarizations: $$ F_{x}= \frac{-m}{2}[ h_{\times}(t) y +h_{+}(t) x],$$ $$F_{y}= \frac{m}{2}[ h_{+}(t) y- h_{\times}(t)x ].$$ The following parameters are recommended (taken from measurements of the Alpha Centauri system, a group of stars about 4 light-years from Earth). | Parameter | Symbol | Value | | ------------- |:-------------:| -----:| | Reduced mass of system | $\mu$ | $ 0.50~ M_{\odot} $| | Angular frequency of source | $\omega_s$ | $2.5\cdot 10^{-9} ~rad/s $ | | Radius of source orbit | $A$ | $ 18~AU $ | | Distance from test masses to source| $r$ | $267,000~ AU $ | |Radius of ring of test masses |$R$ | $6.0~m$| Masses are conveniently expressed in "solar masses," $m_\odot = 1.989*10^{30}~kg$. An astronomical unit is $AU = 1.496\cdot10^{11}~m$. When you write your code, make sure to convert all units to SI to avoid inconsistencies. If you run your code with these parameters, you won't see any displacement of the test masses. They will just sit there. This makes sense because if ripples in spacetime were large enough to see the distortion of objects, we would see them in real life, and wouldn't need a kilometer-long detector to find them. In reality, gravitational waves, while omnipresent, are so small that they require heroic experimental efforts to detect even with 21st century instruments. To be able to get visible results in our program, we'll need to exaggerate the amplitude of the wave. * Scale up your parameters so that the apparent amplitude of the wave is $~0.2$ or higher. Plausibility checks * Do you observe that each mass is moving in a small circle? * Do the masses oscillate around $(x_0,y_0)$? If the motion is shifted, the initial position is not quite right. * Do the masses stay near $(x_0,y_0)$, or do they fly off to infinity? If they fly off, the initial velocity should be checked. Also, make sure that you have approximated $(x,y)$ in the force as $(x_0,y_0)$. This circling behavior is a signature of a circularly-polarized wave. As written, the masses circulate counter-clockwise; this is called "right-hand circularly polarized," since it matches the orientation of your right hand when your thumb points back at you. * Convert your wave to a left-hand circularly polarized wave by multiplying $F_y$, but not $F_x$, by a minus sign. (Make sure to modify $v_{iy}$ and $y_i$ accordingly). Check your output to make sure you get the masses circulating in small clockwise circles. We have been assuming that the detector is on the z-axis, so that it looks at the orbital plane "head on." If the orbit of this source had been oriented differently relative to the detector, the polarization would no longer be purely circularly-polarized. Thus, experimentalists are interested in determining the polarization of a gravitation wave, since it reveals the orientation of the orbit of the source. **Exercise 5: Detecting Gravitational Waves** In the previous section, we found that the gravitational waves produced by Alpha Centauri have tiny amplitudes - on the order of $10^{-23}$ when they reach Earth. Surprisingly, it is not this tiny amplitude that is the biggest limitation to detecting the stars' gravitational waves. Gravitational wave detectors, such as LIGO, can only sense gravitational waves with frequencies above $1~Hz$. The optimum frequency for Advanced LIGO is $10^2~ - ~10^3 ~Hz$. For our purposes, we'll assume that our detector can pick up a gravitational wave signal if its amplitude is at least $10^{-22}$ and its frequency falls within this range. Frequency, when measured in Hz, indicates how many cycles per second the system is making. How can a binary star system possibly make 10 revolutions *per second*? Would it have to go faster than the speed of light? Let's find out. The speed of an object moving in a circle is determined by both its frequency and the radius of the orbit ($A$). * Use the formula for tangential velocity ($v = \omega A$) and a frequency of $10~Hz$ to determine how large the radius of the orbit can be, while keeping $v$ no more than one-tenth the speed of light. (Recall that $\omega = 2\pi f$). Your answer will turn out to be smaller than the radius of the Earth (6,000 km). And yet, this is a typical scale for the orbit of "compact binaries" - binary systems composed of neutron stars and/or black holes. These are the sources of gravitational waves that LIGO is designed to find. Our final goal is to simulate the effect of a compact binary system. First, we need to make sure that the value used for the frequency of the orbit is consistent with the orbital radius - as stars get closer to each other, it takes less and less time for them to complete an orbit. In this simulation, we ignore the inward spiralling of the orbit, and approximate the radius as constant. Undergoing uniform circular motion at radius $A$, the centripetal acceleration of the system is $$a_c = \frac{v^2}{A} = \frac{G(m_1+m_2)}{A^2}. $$ This equation, combined with $v = \omega A$, requires $A =( 2 G m_s/\omega_s^{2})^{1/3}$. * With the values in the table below (and eq. \ref{mu}), simulate a binary neutron star system. Copy your code from Exercise 3 into a new file so that you don't overwrite it. | Parameter | Symbol | Value | | ------------- |:-------------:| -----:| | Mass of each star | $m_i$ | $1.4~ M_{\odot} $| | Angular frequency of source | $\omega_s$ | $2\pi(10~Hz) ~rad/s $ | | Radius of source orbit | $A$ | $ (2 G m_s/\omega_s^{2})^{1/3}$ | | Distance from test masses to source| $r$ |(variable) | * Experiment with the variable $r$ (keeping it much larger than the radius of the orbit, $A$). Can the detector get close enough to the star system to make the effect of gravitational wave visible to the naked eye? - How far away can a compact binary be and still be detectable by our detector - assumed to be sensitive to amplitudes as small as $10^{-22}$? Express the answer in Mega-parsecs ($1 Mpc = 3.086\cdot 10^{22}$ m).