Taken from here Why Parallelism? A parallel computer is a collection of compute modules which work together to perform a task efficiently In principle, more processors means more power However, a major bottleneck is communication speed Think of a room of people working together on a calculation. Information needs to travel between them for them to cooperate efficiently. Better communication == faster results If communication takes longer than the task at hand, then the upsides of parallelism go away Some processors might have easier tasks than others, or are more efficient than other nodes. So even if you have a bunch of people, the task moves as slow as the hardest subtask We want to evenly distribute the work amongst all processors Given the above, we want to write parallel programs that scale well Each piece can safely be done in parallel Keep all the processors busy Keep synchronization/ comms to minimum Superscalar Processors A program is really just a series of instructions encoded in binary that a particular CPU can understand. The rate at which these instructions are executed is defined by the clock For simplicity, assume we have a CPU that executes exactly one instruction per cycle At it’s core then, CPU need to: Fetch/Decode: Figure out what instruction to execute Execute: Load the required data and do the thing Write Back: Update he internal state of the system Suppose that we have the program: $a = x^{2}+y^{2}+z^{2}$ The above a the sequential version of the program What happens if you could execute two instructions at the same time? You could run $x^{2}$ and $y^{2}$ in parallel, then add them in parallel with doing the z multiplication, then add up the final results. So now it takes 3 clock cycles This is the idea behind instruction level parallelism (ILP)/ superscalar processors The key idea: some instructions can be run at the same time without interfering with each other. Hence, you can run multiple constructions at the same time! A real world example: the Pentium 4 Single Thread Wall Unfortunately, the ILP gains seems to empirically saturate after 4 instructions In addition, there are hard power limits based on the frequency and voltages (EMI and all that) If chips get too hot, the don’t function anymore The way forward is multiple processors working in tandem Or alternatively, specialized processors for specific workloads Memory There’s the old not-very-funny C joke: “The world is just a giant array of undifferentiated bits” In the programmer’s mental model, memory is exactly this Since most types of memory are not stored on the processor, you need to go fetch it from somewhere This takes a long time if you are getting this from the disk (Memory access times can be 100’s of cycles) Caches Caches act as a way to help with these latencies: By providing an intermediate storage of highly accessed values, you reduce latency Caches should not alter the correctness of a program, only the performance Cache has two methods of “data locality”: Spatial locality: nearby addresses are more likely to be accessed (think array iteration) Temporal locality: high frequency addresses will more likely be accessed in the future Cache can have some hierarchy (L1,L2, L3) it could be implemented as DRAM (variable cycle count) but the important part is that whatever the caching system should be opaque (most) programmers The smallest unit of a cache is a “line” which has some width. If you get 1 address from memory, the adjacent memory values are put into a cache line. Subsequent memory addresses check the cache first before going to memory There are a bunch of policies for how the cache decides what data is high frequency. For now, just assume that a cache of size N bytes stores the last N addresses accessed New values follow LRU policy (least recently used) ie. the oldest access gets replaced first As a concrete example, say you iterate through addresses 0x0 and 0xF with a cache width of 8 bytes and a line width 4 bytes. The above shows the access pattern You get a cache miss every 4 bytes, and the least recently used cache address gets evicted if the cache is full Modern CPUs make predictions about what data will be fetched (so called “pre-fetching”). If the CPU guesses wrong, you can get reduced performance Stalling Suppose that you miss the cache, and are now spending a very long time getting memory from DRAM/disk While you are waiting for the access, you can swap to a different program and do work there There might be a period where the first thread finished stalling and is in a runnable state, but is not being ran by the processor This is OK if aiming for a high throughput system. However, there is a “no free lunch theorem” trade off between throughput (ie. pipelining a bunch of work) and latency (ie. reducing the time each task takes ) A CPU’s utilization rate is defined by the time spend doing actual work instead of stalling. This heavily depends on the number if instructions you can run before stalling. Adding more threads once you are at 100% utilization rate doesn’t increase the utilization. It just delays when you get back to the original task Moving Past Single Core Since you hit a wall in single core performance, one way to increase performance is to use the transistors to make multiple cores instead of making the ILP logic properly Another way is to use SIMD (Single Instruction Multiple Data) You have one register with a bunch of bits in them (for concreteness, say 256 bits). You then subdivide each subsection of this large register (say 64 bit chunks). Each chunk then gets passed into their own ALU Hence the name: a single instruction can work on the multiple chunks of data in a register For SIMD, you need to ensure that each subsequent loop executes the same set of instructions across the same elements (called “coherent execution” or “instruction stream coherence”) Coherent execution is NOT a requirement for parallelism across multiple cores If a program lacks instruction stream coherence, it is called “divergent execution” On a GPU, you have N instances of a program running in parallel on a large number of different cores. Divergent execution in this context can cause massive performance hits SIMD can be expressed at the compiler level (via AVX2 instructions on intel for instance) or can be done implicitly by the hardware Branches Remember that you are assuming that all iterations of the loop you are optimizing are independent of each other. How do you deal with conditional branches (re: if/else kind of statements)? Idea: Have all ALU cores try to run the if branch. If the condition is satisfied for that particular iteration, enable the computation, otherwise disable the ALU core. After that, repeat for the else branch, but disable the other set of registers Latency Issues Think of driving down the highway. If you wanted to increase the throughput (ie. how many cars per hour can move along the highway), you need to either: Increase the speed of each car Increase the number of lanes Memory is the same. There is some maximum throughput/bandwidth which you can transfer data on a channel Memory Latency is defined as the time it takes to transfer any one item You can increase the throughput of a system while maintaining a fixed latency by pipelining your system (ie. when one subprocess finishes, immediately begin a new subprocess) This approach is latency limited: you can only go as fast as the slowest subprocess Incidentally, this is how multiplication takes “a single cycle”. This refers to instruction throughput, not latency Ie. it takes like 4 cycles to run a multiply instruction end to end, but only 1 cycle assuming proper pipelining To go to the extreme, this of a GPU. The max bandwidth is in the TB/sec. Transferring enough memory to saturate this bandwidth is physically impossible However, even a 1 percent utilization rate is still much faster than an 8 core CPU Modern computing therefore is constrained by the memory bandwidth! Performant programs fetch data from memory less often Doing more math instead of more memory requests can be beneficial to performance (up to a point) Parallel Programming Idioms/ ISPC ispc is the Intel SPMD Program Compiler which provides programming primitives Calling a ISPC function spawns a “gang” of ISPC program instances which run concurrently Each instance has it’s own copy of local variables Some keywords unique to ISPC: uniform: type modifier which asserts that all instances use the same value for this variable (non necesssary for correctness) programCount: the number of simultaneously executing instances in a gang (uniform value) programIndex: the id of the current instance in the gang foreach: a parallel loop keyword. This demands that work required to execute all of these loops be done by the gang as a whole, and not each instance It’s the compiler’s job to assign iterations to program instances in the gang There are a bunch of caveats where this will either be undefined or give a compiler error If multiple iterations of a loop body can write to the same memory location (race conditions are bad) Cannot return many copies of a variable which expects one return value cannot have a uniform value which gets assigned a different value across instances The above program interleaves the assignment of program instances (Ex: instance 0 works on indices 0, 8, 16 …) Temporally, at time t=0, each instance is working on a contiguous set of indices, so you can run a packed vector load You could also do “blocked” assignment, where each instance does a contiguous set of loop indices With this setup, you need to access disparate memory locations at t=0 Requires “gather” instructions (vgatherdps) to implement (which can be costly SIMD) If you need to talk between instances, ISPC gives you some intrinsics: uniform int64 reduce_add(int32 x);: compute sum of a variable across all program instances in a gang uniform int32 reduce_min(int32 a);: compute the min of all values in a gang int32 broadcast(int32 value, uniform int index);: Broadcast a value from one instance to all instances in a gang int32 rotate(int32 value, uniform int offset);: For all i, pass value from instance i to instance i+offset% programCount There is a difference between the semantics of an abstraction (ie. what should something do?) and the implementation (ie. how is something implemented) A great deal of caution should be made to disentangle an implementation from the semantics For ISPC, the abstraction is that the compiler spawns programCount logical instruction streams The implementation is the compiler emitting SIMD instructions with the appropriate masking There is some leakage with the uniform keyword (in that you need to think about if a value is shared across all instances) This abstraction/implementation is not unique! What if you didn’t want to provide fine-grain control of what happens in each instance? You would remove programIndex and programCount and be forced to use uniform outside the foreach loop You could also disallow array indexing and force the programmer too think about collections of data (think numpy!)

Cell Overview Plasma Membrane Outer shell of the cell responsible for seperating the cell from the enviornment Regulates transport of nutrients and waste of the cell Consists of a phospholipic bilayer (seperates hydrophobic materials on the interior of the membrane from the aqeous enviornment The heads are hydrophilic, while the tails are hydrophobic Lots of things are embedded in and pierce the bilayer (peripheral proteins, protein channels, Alpha-Helix proteins etc.) The proteins which pierce the membrane can be asymetric (re: inside has one structure, outside has a different structure) Cytoplasm refers to the intracellular content 80% of interior is water consists of cytosol (intracellular fluid) and organelles (internal structures) cytosol is a non-newtonian fluid and is colorless Cytoskeleton Protein filaments in the cytoplasm Holds cellular structure, aids in migration, signalling, and chromosome segregation during cytokinesis There are more specialized cytoskeleton structures: flagella, cilia, filopodia, lammelipodia etc. There are 3 board classes of cytoskeleton: actin, microtubules, and intermediate filaments. Actin are physically the smallest, microtubules are physically the largest Cell shape dictated by microtubules main branches, with smaller actin branches doing fine shape Organelles Nucleus Nucleus is present in eukaryotic cells It’s enclosed by double membrane (re: two bilayers) supported by lamina (re: cytoskeleton support) The outermost membrane folds onto itself into the endoplasmic reticulum (ER) nuclear pores pieces into the interior of the nucleus the interior is filled with chromatin (linear DNA molecules in complex with other proteins) has other suborganelles (nucleolus/speckles) Purpose: stores genetic material site of DNA replication and translation ribosome synthesis Endoplasmatic Reticulum Purpose: lipid synthesis membrane protein synthesis Ca++ ion storage detoxification Network of interconnected closed membrane tubules and vesicles composed of smooth ER and rough ER (called “rough” because it has ribosomes) The rough ER helps transpose hydrophobic materials within their bilayers as well as ribosomes The smooth ER generates new lipids and membranes Golgi Apparatus Looks like the rough ER Packages ER produces into small vesicles (exocytotic, secretory, lysosomal) which get transported to their final destination Lysosomes single membrance vesicle which decomposes things pH of 5 Peroxisomes Like lysosomes, but specifically fatty acids and toxic compounds single membrance They have a crystalline core Mitochondrium drives ATP production for aerobic metabolism has an double membrane, cristae, and a matrix Chloroplast Photosynthesis site for plants Plants also have mitocondria for some reason Tissue Cultures The term “tissue culture” is used to denote a population of cells to be examined “Culturing” a tissue means to grow them in a lab Can be “In Vitro” (in a test tube /glass) or “In vivo” (in a living organism) Depends on the subject of study Why use cultures? Controlled enviornment Ease of access Problems with tissue culture: They are outside the body, so they aren’t a perfect replica of the system Cells grow in 2D in vitro, while they grow in 3D in humans Human cells divide roughly once a day in vitro, but rarely do in vivo First proof of concept of tissue culture in 1885 by Wilhelm Roux was able to maintain chicken embryos in vitro for a couple of days HeLa cells was the first cell culture line (1951) From Henrietta Lack’s cervical cancer cells Tissue Cultures are grown in dishes with a medium that has the correct nutrients, pH buffer ,indicator etc. These dishes are stored in incubators To work with these cells, you need to shove them into a sterile fume hood Microscopy in Cell Biology Light Microscopy resolution of 200 nm Oldest form of microscopy. Passes light through a thin section of cell tissue More modern versions involve more lenses, but the same idea holds For liquid sames, we need to invert the microscope (light coming from above, objective lenses below) for many lense systems, you can insert “tube lenses” to circumvent the finite focal lengths of lenses We can modulate the incoming light of a microscope to generate constrast Humans are good at seeing changes in wavelength and amplitude, but bad at seeing phase and polarization changes Most biological specimens respond to phase, so the trick is to convert the phase shifts to something that we can see The idea is that: you let light pass through your cells. Most will pass through just fine, but some will get phase shifted as they pass through things Using a phase ring, you shift the diffracted light by a factor of $\frac{\pi}{4}$, leaving the phase of the unpeturbed light unchanged Damp the amplitude of the unpeturbed light to the amplitude of the phase-modulated version Look at the interference patterns between the two Another method (Differential Interference Constrast Microscopy (DIC)) Send linearly polarized light through your sample You have two equal polarization components Use a Wollaston prism to seperate the two components away from each other Pass the components through your specimen, which shifts the two polarizations a different amount Recombine the components with a prism, then pass through an analyzer lens Fluorescense Microscopy Idea: Shine one wavelength in, the dye you put into the cells reacts and produces a difference wavelength which you can see The shift between the maximum of the emission and the absorption spectra is called the Stokes shift You are limited to 4 colors: Constrained to visible spectrum for obvious reasons UV regime can be cyto-toxic to your cells These is also some spread associated with each source How do you make biological molecules fluorescent? You get lucky and there is some compound which can naturally bind to the molecule you want to study Phalloidin: a mushroom toxin that binds to actin Mitotracker: synthetic binder which binds mitochondria Design florescent anti-bodies which bind to the target structure You permeate the cell (poke a bunch of holes into it), flush a bunch of these antibodies in and let them bind, and then flush a clean buffer solution to clear up the free floaters obviously, this only works on dead cells Direct chemical labelling for your Molecule of Interest GFP (Green Florescent Protein) Jellyfish proteins which can attached to living cell structures Accomplished by modifying the genes which produce a protein of interest by attaching the genome of GFP in an appropriate place Can attach to nearly any protein No need for the microinjection with anti-bodies Electron Microscopy Resolution of 1 nm Uses DeBroglie wavelength of electron to increase the resolution compared to optics Uses tunnelling to produce a modulating current which gets reconstructed to an image Disadvantages: Very high vacuum required Specimens need to be fixed, embedded, sectioned an stained with an electron-dense material Central Dogma of Molecular Biology DNA (deoxyribonucleic acid) to RNA (ribonucleic acid) via transcription RNA to Protein via translation DNA A polymer of highly charged polyelectrolytes the monomers of DNA (nucleotides) are made of a phosphate group, a sugar, and an organic base The sugar for DNA is deoxyribose, while for RNA it’s ribose The organic bases are thymine (T), cytosine(C), adenine(A), guanine(G), and uracil(U) T for DNA and U for RNA T and A are purines (double ring structures) while C, G and U are pyrimadine (single ring) A and T pair up for DNA (A and U for RNA), while C and G pair up in both You match a purine and pyrimadine together to maintain the width of DNA The pairing occurs via hydrogen bonding (A and T have 2 hydrogen bonds, while C and G have 3 hydrogen bonds) When describing the double band structure, there is a 3’ and a 5’ end The 3’ side is attached to an oxygen of a phosphate The 5’ side is connected to the $CH_{2}$ side The 3’ and 5’ sides have opposite polarity on each strand There are 3 variants of DNA (A,B and Z) DNA is not in isolation in the nucleus. They are complexed with histones Since DNA is negatively charged from the phosphates, the histones are positively charged The charge on the histones is modulated by the acetylation of its’ tails RNA RNA is shares a lot of similarities to DNA One big difference is that RNA has a flexible backbone The three main types of RNA are messenger RNA, transfer RNA and ribosomal RNA RNA Polymerase There are 3 polymerases which synthesis the RNA (RNA polymerase I synthesises rRNA. II for mRNA and II for tRNA) These polymerases need a couple of things to do transcription: A start sequence The TATA box is a common one, but it’s not unique Transcription factors are used to unwind and prepare the DNA A stop sequence like, 200 A base pairs a 5’ cap The polymerase starts translating from the 3’ to the 5’ end (the strand that it works on is called the leading strand) For mRNA, there are exon (protein encoding regions) and introns (non-protein encoding regions) Proteins Consists of polymers A string of 20 monomers Each monomer has 3 components: The amino group ($NH_{3}$) The acid group ($CO_{2}H$) A central carbon atom, which bridges the acid and amino groups (via single bond to N and C) The R group: Gives all of the variance between the monomers The monomers get chained together via peptide bonds The process forms a single covalent bond between the carbon and nitrogen of the acid and amino group, producing water as a biproduct The menagerie of Amino Acids side groups: Nonpolar, aliphatic Aromatic Positively charged Polar, uncharged Negatively charged There is a hierarchy to Proteins Primary: the sequence of monomers Secondary structure: $\alpha$ helices an $\beta$ pleated sheets Tertiary structure: the 3D shape of your protein Quaternary structure: multiple proteins complexing together There is a well known mapping from RNA base pairs to monomers (re: codons) Lipids Lipids are partitioned into polar heads and nonpolar tails This partitioning allows lipid to self assemble (align the heads towards water and tails away from water) There is also a whole zoo of lipids: phospholipids gangliosides free fatty acids Cholesterol Steroids Carbohydrates Vary in size: monosaccarides, oligosacharides, and polysacchardies the number of rings determines which category they fall into Consists of hydrogens and carbons (and some oxygens) They can be energy sources, or serve a structural function (cellulose and chitin) Mesoscopic Forces Van der Waals These are dipole-dipole interactions There are 3 types: between 2 permanent dipoles (Keesom orientation energy) between 2 induced dipoles (London’s dispersion force) between 1 permanent and 1 induced (Debye energy) The interaction potential of a Van der Waals force between two point particles scales like $r^{-6}$ This does not hold for other geometries. You can derive other geometries via taking infinitesimal masses which do scale like $R^{-6}$ The proportionality constants come from quantum mechanic (Hamaker constant) They are long range (> 10 nm) They are weak (~ 1 $kJmol^{-1}$) Hydrogen Bonds Occurs between a donor (strong polar group) and a proton acceptor (slightly electronegative) They are weaker than ionic and covalent bonds (~ 10-40 $kJmol^{-1}$) Used in molecular self-assembly Hydration forces The polar nature of water leads to an electrostatic interaction The large dipole moment of water + water’s capacity for hydrogen bonding gives rise to a short range-order of water molecules Described by an exponential decay ($f = f_{0} exp(-\frac{t}{\lambda})$) These are very short range (~1 nm) Electrostatics point to point (Coulomb): scales like $r^{-2}$ point to dipole: scales like $r^{-3}$ dipole to dipole: scales like $r^{-3}$ the dominant long-range interaction Debye-Hueckel Theory There is some screening of counter ions around the macro ion The distribution of the counter ions as a function of the potential is $\rho_{i}(\phi) = \rho_{0} exp(-\frac{\phi}{kT})$ You can plug this charge distribution into the Poisson equation to yield: $\Delta \phi(\vec{r}) = \frac{1}{\epsilon_{0}\epsilon} \Sigma_{i=1}^{N} z_{i} e \rho_{i}(\infty) exp(-\frac{ - z_{i}e \phi(\vec{r})}{kT})$ DLVO Theory A Theory which describes forces between charged colloidal particles in a solution $\frac{V(r)}{B} = -\frac{A}{12\pi r^{2}} + \frac{64 k T \epsilon_{0} \Gamma_{0}^{2}}{\kappa} exp(\kappa T)$ There is once again, competition between Van der Walls attraction and electrostatic repulsion Depletion Forces Mediated by polymers or small particles Imagine that you have an isotropic distribution of smaller particles (polymer chains for instance) If you have two particles, then they feel an isotropic pressure from these particles If these two particles are close enough (re: smaller than the average size of the polymer chains), then this will break the isotropy and push the particles together generally attractive driven by entropy gain of polymers leads to colloidal aggregation The joining pressure is $\Pi = \frac{N}{V} kT$ Undulation Forces Imagine that you have a very think membrane which you stretch out to sizes that are several orders of magnitude larger than the thickness An energy of size kT is sufficient to substantially deform the system Steric Forces mediated by polymers The polymers act as a spring which repells particles away from a surface The repulsive energy per unit area is $W(r) \approx 36 k T e^{-\frac{r}{R_{g}}}$ Bridging Flocculation The polymers get attached to colloids, and between two colloids, the polymers intertwine Aggregating Self-Assembly Systems try to self-assemble so as to minimize $F_{min}$ errors in assembly can lead to disease: sicle cell anemia amyloidosis (prions) Alzheimer’s disease (self-assembled beta sheets amyloidplaques) Self-Assemblies: aggregating (micellization of lipids) non-aggregating (folding of globular proteins) morphologies produced by phase separation Self-Organization: actively bringing/moving parts together (re mitosis) Entropy loss of the assembling molecule is counteracted by entropy gain in solvent molecules Self assemblies aim to minimize the surface energy, which runs contrary to Entropy maximization Ostwald Ripening You have two particles which exchange particles via the solvent This is a purely thermodynamic process Surfactants they are ambiphiles (they have both hydrophilic and hydrophobic effects) Phospholipids are an example of this type of molecules These lipids tend to get attached to surfaces first prior to forming self assemblies These phospholipids can get joined at their tails to self-assemble Liquid Crystals Mesostate (between solid and liquid) Anisotropic (preferred direction in behavior) Ability to flow (liquid-like) Lots of biological molecules for these crystals biologically important (membranes, spider slik, cellulose etc.) Roughly 4 types: Rod-like, disk-like, banana-like, and conical-like Nematic: molecules tend to align in the same direction Semetic: molecules tend to arrange themselves in layers The order parameter $S = <P_{2}(\cos\theta)>$ determines which category the LC falls into $\theta$ refers to the angle between the z-axis and the primary axis of the LC $P_{2}$ is the 2nd order Legendre polynomial Diffusion Brownian Motion Discovered by Brown in 1826 via pollen wiggling around Einstein in 1905 gave kinematic formulation for Brownian motion Random Walk You have a process described by: $x_{i}(n) = x_{i}(n-1)_\epsilon$ On average, the displacement is 0 The mean squared displacement is non-zero ($n\epsilon^{2}$) This can be related to the diffusion constant: $D = \frac{\epsilon^{2}}{2\tau}$ $tau$ is the time step of your simulation On a macroscopic level, Fick’s Law states: $J_{x} = -D \frac{\partial c}{\partial x}$ J is the particle flux and c is the concentration These concentrations are very dependent on the boundary conditions of your system

Following Misner, Wheeler and Thorne. Foundations of General Relativity All laws of physics can be expressed geometrically Local reference frames are flat gravitational mass is equivalent to inertial mass All reference frames are equivalent (re: there is no preferred reference frame) Minimal coupling: work in SR to find physical laws , then go to general frame to include gravity Basic Definitions and Math Things Take naturalized units (c=G=1) when working with GR Events: Things that happen at a particular point in spacetime. Does not depend on the underlying spacetime structure (re: could be flat or curved) Curves: A curve is some path through spacetime $P(\lambda)$. Once again, this is independent of the notion of any underlying spacetime (although a metric is needed to ascribe a length to the curve) Coordinates: How you label events in spacetime. There is no unique way of assigning coordinates Written as $x_{\mu}$, where $\mu$ ranges from 0 to 3 A coordinate transformation is a set of 4 equations. Each one determines how the new coordinate depends on the old 4 coordinates Coordinate singularities can occur (think of the north pole of the earth on latitude and longitude) These can be dealt with via multiple patches of coordinates Vectors: The seperation between two events in spacetime. In flat spacetime, this is the difference between the two coordinates assigned to each event. This falls apart in curved spacetime, but is an good approximation for arbitrarily close events These transform under coordinate transformations as: $\epsilon^{\beta} = \frac{\partial x^{\alpha}}{\partial x^{\beta}} \epsilon^{\beta}$ Follows from Taylor expansion Can reduce the notion of a vector from one that requires 2 events to one that requires 1 Imagine a parameterized line between the two events: $P(\lambda) = A + \lambda(B-A)$, where $0\leq \lambda \leq 1$ Taking the derivative w.r.t. $\lambda$ evaluated at $\lambda=0$ gives P(1)-P(0). This construction $\frac{dP}{d\lambda}_{\lambda=0}$ is a 1 point object called the tangent vector Summation Notation: $\epsilon_{\alpha}\gamma^{\alpha} = \Sigma_{\alpha=0}^{3} \epsilon_{\alpha}\gamma^{\alpha}$ The metric tensor: a machine for computing scalar products of vectors. It takes in two vectors and spits out a scalar It’s symmetric in it’s arguments: $g(u,v) = g(v,u)$ It’s linear w.r.t. it’s arguments: $g(a\vec{u}+b\vec{v},\vec{w}) = a g(\vec{u},\vec{w}) + b g(\vec{v},\vec{w})$ If you know the basis vectors of the frame ($\vec{e_{\alpha}}$), then you can calculate it’s output for any input (follows from linearity) Define the metric coefficients to be $g(e_{\alpha}, e_{\beta}) = \vec{e_{\alpha}} \cdot \vec{e}_{\beta} $ The scalar product then becomes: $u^{\alpha}v^{\beta}g_{\alpha\beta}$ In special relativity, this is called $\eta_{\alpha\beta}$, which is diagonal with time being -1 and the space coordinate being 1 1-forms: The 3rd class of geometric objects Think of the 4-momentum vector $p_{alpha} = m\mu_{\alpha}$ The de Broglie wavelength gives an alternative interpretation for momentum. If you diffract the wave on a lattice and observe the diffraction pattern, you can then map surfaces of equal phase. This series of surfaces are a one-form If you run a vector through this one form from ove event to another, then you can calculate the phase difference between the point via $<\tilde{k},\vec{v}>$ You can think of the 1-form as the local form of these equal phase surfaces (re: the best linear approximation of the phase $\phi$ near an event) More mathematically, you can think of a 1 form as a linear, real-valued funtion of vectors (re: something that maps vectors to scalars) the set of all 1-forms at a given event is a “vector space” For each vector $\vec{p}$, the is a unique 1-form defined by $\vec{p}\cdot \vec{v} <\tilde{p},\vec{v}>$ which holds for all $\vec{v}$ In words: the projection of v onto the 4 momentum equals the number of surfaces the vector v pieces of the 1-form The gradient of a scalar $\mathbf{d}f$ is the simplest example of a 1-form The more familiar vector notion of the gradient is the unique vector associated with the 1-form notion of the gradient Define $\partial_{\vec{v}} = (\frac{d}{d\lambda}) $ at $\lambda=0$ along the curve $P(\lambda)-P(0) = \lambda v$ as the directional derivative The gradient describes the first order changes in f in the neighborhood of $P_{0}$: $f(P) = f(P_{0}) <\mathbf{d}f, P-P_{0}>$ The gradient is related to the directional derivative. Apply the directional derivative to a scalar function f: $\partial_{\vec{v}} f = <\mathbf{d}f, \frac{dP}{d \lambda} = <\mathbf{d}f, \vec{v}>$ Tensors: The generalization of these lower dimensional objects. They can have multiple indices which transform under the Poincare group Think of tensors as multilinear maps which take in an appropriate number of one-forms and vectors and then outputs as scalar The order in which you insert vectors/1-forms matters If one knows the value of a tensor in 1 frame for some basis vectors/1-forms, then you can transform to any other rest frame (just transform each of vectors/1-forms with their appropriate matrices, then tack those onto the tensor to get the transformation law) We define the “rank” of the tensor as the number of indices used to describe it Tensor Operations (for example’s sake, we will deal with a 4-tensor): Gradient: $\nabla S(u,v,w,\epsilon)$ is $\partial_{\epsilon} S(\vec{u},\vec{v},\vec{w})$ with u,v and w fixed This is roughly the difference of S(u,v,w) from the tip of $\epsilon$ to the tail of $\epsilon$ In component notation: $\nabla S(u,v,w,\epsilon) = (\frac{\partial S_{\alpha\beta\gamma}}{\partial x^{\alpha}} \epsilon^{\delta}) u^{\alpha}v^{\beta}w^{\gamma} = S_{\alpha\beta\gamma,\delta} u^{\alpha}v^{\beta}w^{\gamma}\epsilon^{\delta}$ This increases the rank of the tensor by 1 Contraction: You can reduce the rank of a tensor by 2 via contraction $M(u,v) = \Sigma_{a=0}^{3} R(e_{a}, u, w^{a}, v) $ You need to have one index up and one index down prior to summing them Divergence: You take the gradient of a tensor, and then contract it with one of the original slots Hence $\nabla \cdot S = S^{\alpha}_{\beta\gamma,\alpha}$ If you take the divergence of S on the first slot Tranpose: interchange two slots for each other Symmetrization and Anti-symmetrization A tensor is symmetric if any transposition of it’s indices yields the same tensor A tensor is antisymmetric if swapping indicies causes the tensor to reverse sign Wedge product: a way to construct a completely antisymmetric tensor from a set of vectors and 1-forms Define the “bivector” $u \wedge v = u\otimes v - v\otimes u$ Define the “2-form” $\alpha \wedge \beta = \alpha \otimes \beta - \beta \otimes \alpha$ If a vector is a linear combination of other vectors in the set, the wedge product between all the vectors is 0 You can define a generalized p-form as a completely anti-symmetric tensor of rank p, which you need to normalize by $p!$. Wheather there is a - or + is determined by the Levi-Civiti symbol Wedge products obey distributive and addition, and commute with scalar multiplication and addition If you are commuting wedge products (say that you have a p-form a, and a q-form b): $\alpha \wedge \beta = (-1)^{pq} \beta \wedge \alpha$ Duals: For vectors, antisymmetric rank 2 tensors, and antisymmetric rank 3 tensors: $*J_{\alpha\beta\gamma} = J^{\mu}\epsilon_{\mu\alpha\beta\gamma}$ $*F_{\alpha\beta} = \frac{1}{2} F^{\mu\nu}\epsilon_{\mu\nu\alpha\beta}$ $*B_{\alpha} = \frac{1}{3!} B^{\lambda\mu\nu}\epsilon_{\lambda\mu\nu\alpha}$ More generally, the dual is denoted as Hodge star (*) and describes the one to one correspondence between a p form and an n-p form Requires specifying the metric for this to make sense $*(\epsilon^{i_1}\wedge … \wedge \epsilon^{i_p}) = \frac{1}{(n-p)!}{\epsilon^{i_1…i_{p}}}_{j{p+1}… j{n}}\epsilon^{j{p+1}}\wedge … \epsilon^{jn}$ Uses the metric to raise the first p indices of $\epsilon_{1…n} = \sqrt{-g} \tilde{\epsilon}_{1…n}$ where $\tilde{e}$ is the Levi Civiti symbol The $\sqrt{g}$ is needed to make it transform like a tensor Reference frames: a system to relate two coordinate systems to each other Intertial reference frames are the preferred class since they require the minimal set of forces needed to describe motion. Non-interital frames tack on additional forces (re: forces that arise from accelerating w.r.t. the inertial one) Special Relativity Postulates: All inertial frames are equivalent There is a maximum speed that information can travel at (happens to be c) Lorentz Transformations Define some set of matrices $\Lambda_{\alpha}^{\beta}$ which transforms your coordinates from one frame to another $x^{a} = \Lambda_{b}^{a} x^{b}$ and $x^{b} = \Lambda_{a}^{b} x^{a}$ Going from a primed frame to an unprimed frame and back is the identity. Hence the Lorentz transform must be inverses of each other: $\Lambda_{\beta}^{\alpha} \Lambda_{\gamma}^{\beta} = \delta^{\alpha}_{\gamma}$ The Lorentz transforms define how the basis vectors and 1 forms transforms via the Einstein summation notation In matrix notation, suppose that your frame is moving along the x axis. Then your Lorentz matrix takes the form $\begin{pmatrix}\gamma & -\beta \gamma & 0 & 0\\ -\beta \gamma & \gamma &0 &0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}$ Apply this matrix to some vector (t,x,y,z) to get the transformations to a primed system moving at some velocity v w.r.t. the unprimed Can derive by noting that: Time dialates: ($\Delta t’ = \gamma \Delta t$) where ’ is the moving frame Lengths along the motion contract ($\Delta x’ = \frac{1}{\gamma}\Delta x$ Lengths perpendicular to the motion are unchanged This can be framed in hyperbolic geometry by making the associations $\beta = \tanh (\alpha)$, $\sinh (\alpha) = \beta \gamma$ and $\cosh (\alpha) = \gamma$ This lets you easily add boosts together by simply adding the $\alpha$ Spacetime Interval The interval is defined as $ds^{2} = -dt^{2} +dx^{2}$. This is invariant under Lorentz transformation This of this as a frame-invariant way of stating the distance between two spacetime events This can also be written as $ds^{2} = \eta_{ab} dx^{a} dx^{b}$ The interval is also invariant under spatial rotations (can think of boosts as rotations in t-x plane) The interval has 3 cases: timelike (<0), spacelike (>0) or null (=0) Timelike implies a causal lightcone structure Spacelike involves a non-causal structure Nulllike implies moving at light speed SR Mechanics Define $\tau$ to be the proper time to be the time experienced by an observer when at rest (ie. $d\tau^{2} = -dt^{2} +dx^{2}$ This proper time will be less than or equal to every other frame (b/c time dialation) $\tau$ is Lorentz invariant Suppose that you parameterize your interval via the proper time $\tau$ $\tau$ is the 4-magnitude of your vector (Think of it as the frame where no velocity is happening) With a set of 4 orthonormal basis vectors $\vec{e_{\alpha}}$, you can construct the 4-velocity u as $\vec{u} = u^{0}\vec{e_{0}} + u^{1}\vec{e_{1}} + u^{2}\vec{e_{2}} + u^{3}\vec{e_{3}}$ $u^{0} = \frac{dt}{d\tau} = \frac{1}{\sqrt{1-v^{2}}}$ $u^{j} = \frac{dx^{j}}{d\tau} = \frac{v^{j}}{\sqrt{1-v^{2}}}$ $v^{2}$ is the ordinary 3d velocity, and $v^{j}$ is that paricular coordinate of the velocity (doesn’t need to be Cartesian ) The magnitude of $v^{2}$ dictates the angle that the 4-velocity makes with the temporal axis This transforms like the coordinate $x^{a}$ since $\tau$ is invariant The norm of u is -1 U is tangent to the particle’s path From the above, we can construct a 4-momentum for massive objects: $p^{a} = mu^{a} = (\gamma m, \gamma m \vec{v})^{a} = (E, \vec{p})^{a}$ The norm of the 4-momentum is $-m^{2}$ For a massless particle, we demand that $p^{a} = (\frac{h}{\lambda}, \frac{h}{\lambda} \vec{e})^{a}$ $\lambda$ is the deBroglie wavelength $\vec{e}$ points along the direction of propagation The four-force is defined as $\vec{F} \frac{d\vec{p}}{dt}$ where t is the coordinate time and p is the relativistic 3 momentum $\vec{p} = \gamma m \vec{v}$ The power is $\vec{v}\cdot\vec{F}$, which can be seen by differentiating $E^{2} = \vec{p}^{2}+m^{2}$ w.r.t. time Flux Densities The relativistic flux of some substance q is given by $j^{a} = (\rho, \vec{J})$ $\rho = frac{dq}{d^{3}\vec{r}}$ $j^{i} = \frac{dx^{i} dq}{d^{4}x}$ $d^{4}x$ is Lorentz invariant ($\det(\Lambda)=1$ is the volume element) If q is conserved, then we have $\partial_{a} j^{a} = 0$ E&M Fields We want to write the Lorentz Force law in some frame-invariant way The standard form: $\frac{\vec{p}}{dt} = e(\vec{E} + \vec{v}\times\vec{B})$ This depends on a lot of frame dependent quantities (3-momentum, the time of a specific frame, E and B fields etc.) The frame-invariant form should only depend on frame invariant quantities. Swap the 3-momentum with the 4 momentum, and swap t for $\tau$. We want to find some machine $\mathbf{F}(u^{\alpha})$(re: tensor) which takes in the 4-velocity and outputs the appropriate derivative of the 4 momentum divided by the charge In index notation: $\frac{dp^{\alpha}}{d\tau} = e F^{\alpha}_{\beta}u^{\beta}$ You can calculate the spatial parts of F via comparing against the standard form of the Lorentz force law You can calculate the time components as $\frac{dp^{0}}{d\tau} = \frac{1}{\sqrt{1-\vec{v}^{2}}} \frac{dE}{dt} = e \vec{E}\cdot\vec{u}$ The end result is that the force law generalizes to $\frac{dp^{\mu}}{d\tau} = e F^{\mu}_{\nu} u^{\nu}$ The transformation of the fields under Lorentz boost can easily be seen by applying a Lorentz transformation on the Faraday tensor People typically take about the “covariant components” of F, which can be found by lowering an index: $F_{\alpha\beta} = \eta_{\alpha\gamma}F^{\gamma}_{\beta}$ Maxwell’s equations can also be written in terms of the Faraday tensor: Define the 4-current J where the time component is $\rho$ and the spatial components are the current density We have that $\mathbf{d}F = F_{\alpha \beta, \gamma} + F_{\beta \gamma, \alpha} + F_{\gamma \alpha, \beta} = 0$ (re: The exterior derivative of the tensor vanishes) where $F_{a,b} = \frac{\partial}{\partial x^{b}} A_{a}$ This encode $\nabla \cdot \vec{B} = 0$ and $\frac{\partial B}{\partial t}+\nabla \times E = 0$ We have that $\mathbf{d*} F = 4\pi * J \rightarrow F^{\alpha\beta}_{,\beta} = 4\pi J^{\alpha}$ (or the exterior derivative of the dual of F is proportional to the dual of J) This encodes $\nabla \cdot E = 4\pi \rho$ and $\frac{\partial E}{\partial t}-\nabla\times b = -4\pi J$ Exterior Calculus Exterior Calculus Formulation of E&M Recall that the Faraday tensor is $F = F_{\alpha\beta} \mathbf{d}x^{\alpha} \otimes \mathbf{d}x^{\beta}$ The Faraday tensor can instead be written in terms of wedge products: $F = \frac{1}{2} F_{\alpha\beta} dx^{\alpha} \wedge dx^{\beta}$ $dx^{\alpha} \wedge dx^{\beta} = dx^{\alpha} \otimes dx^{\beta} -dx^{\beta} dx^{\alpha}$ Any antisymmetric, second-rank tensor can be expanded like this Observe that $B_{2}-E^{2} = \frac{1}{2} F_{\alpha \beta}F^{\alpha\beta}$ and that $\vec{E}\cdot \vec{B} = \frac{1}{4} F_{\alpha\beta}* F^{\alpha\beta}$ The 4-density current is defined as $J^{\mu} = e \int \delta^{4}(x^{nu}-a^{\nu}(a)) \frac{da^{\mu}}{d\alpha} d\alpha$ Can think of a simple 2-form as a “honeycomb” structure with circulation in each of the cells. You can generate this simple 2-form by intersecting two one forms together In general, you can’t think of general 2-forms like this. You need to think of a general 2-form as a sum of simple 2-forms general 2-forms may or may not simplify to a simple 2-form Combinatorics dictates the number of 2-forms you can construct from a set of 1 forms (for 4 basis 1-forms, you have 6 unique basis 2-forms) Just like how 1-forms can convert vectors to numbers via the number of surfaces pierced by the vector, 2 forms convert 2 vectors to a number by the number of cells spliced by the parallelogram generated by the wedge product Hence, a general f-form can be written as $\phi = \frac{1}{f!} \phi_{\alpha_{1}\alpha_{2},…\alpha_{f}} dx^{\alpha_{1}} \wedge … dx^{\alpha_{f}}$ The exterior derivative of a general f-form $d\phi = \frac{1}{f!} \frac{\partial \phi_{\alpha_{1}…\alpha_{f}}}{\partial x^{\alpha_{0}}} dx^{\alpha_{0}}\wedge … dx^{\alpha_{f}}$ Notationally, this is a boldface d We can write $F = \mathbf{d} A$ to automatically satisfy $\mathbf{d} F = 0$ In alternative notation: $F_{\alpha\beta} = \partial_{\alpha}A_{\beta}-\partial_{\beta}A_{\alpha}$ You can use the Lorentz gauge ($\frac{\partial A^{\nu}}{\partial x^{\nu}} = 0$) and the following definition of the 4-density current $J^{\mu} = e\int \delta^{4}(x^{\nu}-a^{\nu}(\alpha))\dot{a^{\mu}}(\alpha) d\alpha$ to get: $\square A_{\mu} = - 4\pi J_{\mu}$ $\square A = (-\frac{\partial^{2} \phi}{\partial t^{2}} + \nabla^{2} )A $ Stress-Energy Tensor and Conservation Laws Spacetime has a river of 4-momentum, which consists of particles moving along their worldlines We can quantify the flow of this river via the stress energy tensor $T_{\mu\nu}$ Define the volume 1-form to be $\Sigma_{\mu} = \epsilon_{\mu\alpha\beta\gamma}A^{\alpha}B^{\beta}C^{\gamma}$ A,B, and C are the sides of the parallelpiped $T(…,\Sigma)$ is the momentum flowing through the box $T^{00}$ is the energy density $T^{j0}$ is the momentum flux is the 4-momentum per unit volume $T^{0k}$ is the energy flux $T^{jk}$ iis the stress tensor (re jth component of force produced by fields at x^{k}$ Define the number-flux vector $S_{A} = N_{A} \vec{u}_{A}$ This is a 4 vector whose time component is the Lorentz contracted number density, and whose spatial components is the flux of particles Perfect Fluids We have a box of particles whose velocities are distributed isotropically The off-diagonal elements are zero, the on diagonal elements are just the rest-energy density and a pressure along each axis In index notation: $T_{\alpha\beta} = \rho u_{\alpha}u_{\beta} + P (\eta_{\alpha\beta}+ u_{\alpha} u_{\beta})$ In geometric language: $T = P \mathbf{g} + (\rho+P) \mathbf{u} \otimes \mathbf{u}$ Maxwell Stress Energy Tensor $T^{\mu\nu} = F^{\mu\alpha}F_{\alpha}^{\nu} - \frac{1}{4}\eta^{\mu\nu} F_{\alpha\beta}F^{\alpha\beta}$ Symmetry of Stress Energy Tensor We know that $T^{0j} = T^{j0}$ (re: momentum density equals energy flux) due to the equivalence of mass and energy For the stress tensor portion, make the following argument: WLOG, examine the torque around the z-axis for a cube of side L $\tau^{z} = -T^{yx} L^{2} \frac{L}{2} + T^{yx} L^{2} (-\frac{L}{2}) - (-T^{xy}L^{2}) \frac{L}{2} - T^{xy}L^{2}(-\frac{L}{2})$ In words, $\tau^{z}$ is the y component of the force on the +x face times the lever arm on the +x face. plus the y-component of force on the -x face times the lever arm to the -x face minus the x component of force on the +y face times the lever arm to +y face times the lever arm to the +y face minus the x-component of force on the -y face times the lever arm to the -y face $\tau^{z} = (T^{xy}-T^{yx})L^{3}$ The moment of inertia on the cube is $T^{00} L^{3} * L^{2}$ Torque decreases as $L^{3}$ and moment of inertia decreases as $L^{5}$, hence the torque could create an arbitrarily high angular acceleration, which is absurd The paradox is avoided if $T^{xy} = T^{yx}$ The same arguments hold for the other axes Conservation of 4-Momentum We know that $\oint_{\partial V} T^{\mu\alpha} d^{3}\Sigma_{a} = \int_{V} T^{\mu\alpha}_{,\alpha} dV$ In words, the flux of 4-momentum outwards across a closed 3-surface must vanish Since the volume we choose is arbitrary, we know that $T^{\mu\alpha}_{,\alpha} = 0$ Conservation of Angular Mommentum We can define a conserved angular momentum $J^{\alpha\beta}$ from the symmetry of the stress energy tensor Define some event A to be your origin which your angular momentum is defined around $x^{\alpha}(A) = a^{\alpha}$ $J^{\alpha\beta\gamma} = (x^{\alpha}-a^{\alpha})T^{\beta\gamma}-(x^{\beta}-a^{\beta})T^{\alpha\gamma}$ $x^{\alpha}-a^{\alpha}$ is the vector separation of the field point x from the origin a $J^{\alpha\beta\gamma}_{,\gamma} = 0$ from the stress-energy tensor symmetry From the vanishing divergence, we know that $\oint_{\partial V}J_{,\gamma}^{\alpha\beta\gamma} d^{3} \Sigma_{\gamma}= 0$ We can see that the spacelike surface of a constant time t is the standard angular momentum $J^{\alpha \beta} = \int J^{\alpha\beta0}dx dy dz$ Accelerated Observers SR works for accelerating observers. Just imagine an interpolation of constant velocity frames which you transition between Define $a^{\alpha} = \frac{du^{\alpha}}{d\tau}$ From $u^{2} = -1$, we can see that $0 = a_{\alpha}u^{\alpha}$ In the rest frame of the passenger (re: instantaneous inertial frame), we have that $a^{\alpha} = (0, \vec{a})$ This frame is the co-moving frame Fermi-Walker Tetrad This comoving frame by construction cannot be global, because you run into paradoxes If you confine this frame locally (Fermi-Walker Tetrad), then you are fine This tetrad is defined as follows: Define a set of 4 basis vectors $e_{\alpha}$ (the subscript denotes vectors,not components of one vector!) Let the $e_{0}$ vector be aligned with the 4-velocity. Define all the other vectors such that they are orthogonal and nonrotating The orthogonal condition can be encoded as $e_{\mu}\cdot e_{\nu} = \eta_{\mu}{\nu}$ For the non-rotating condition, we know that the basis vectors at two successive instants must be related by a Lorentz transformation Since the unit 4-velocity has a constant magnitude, then accelerations therefore imply some “rotation” of the 4-velocity Lorentz boosts can be thought of as “rotations” in the x-t plane Hence, “non-rotating” means that there is no additional boosts or rotations other than the bare minimum required to move the 4-velocity Tensor Algebra With flat spacetime, a global Lorentz frame is impossible The basis vectors change as you move around spacetime Define your new metric to be $g_{\alpha\beta} = e_{\alpha}\cdot e_{\beta}$ $g_{\alpha\beta} g^{\beta\gamma} = \delta_{\alpha}^{\gamma}$ Our Lorentz transformation matrices now become any arbitrary, nonsingular transformation matrices $e_{\beta} = e_{\alpha}L^{\alpha}_{\beta}$ $p^{\beta} = L^{\beta}_{\alpha} p^{\beta}$ $L_{\alpha}^{\beta} L_{\gamma}^{\alpha} = \delta^{\beta}_{\gamma}$ We can define coordinate bases: $e_{\alpha} = \frac{\partial P}{\partial x^{\alpha}}$ $L^{\alpha}_{\beta} = \frac{\partial x^{a}}{\partial x^{b}}$ The Levi-Civita tensor has components which now depend on the bases: $\epsilon_{\alpha \beta \gamma \delta} = \sqrt{-g} [\alpha\beta\gamma\delta]$ g is the determinant of $g_{\alpha\beta}$ $[…]$ means the antisymmetric symbol, where $0123 = 1$ Each spacetime point resides in it’s own tangent plane to the manifold. You need to “parallel transport” a vector at another spacetime point you can associate some tangent vector u with the corresponding directional derivative $\partial_{u}$ in the tangent plane Since a global Lorentz frame is impossible, the best we can do is demand that $g_{\alpha\beta}$ acts like $\eta_{\alpha\beta}$ in the neighborhood of some event $P_{0}$ Hence, we demand that $g_{\mu\nu,\alpha}(P_{0}) = 0$ Define the covariant derivative $\nabla_{u} T$ along a curve $P(\lambda)$ whose tangent vector is $u = \frac{dP}{d\lambda}$ $lim_{\epsilon\rightarrow 0} \frac{T(P(\epsilon)-T(P(0))}{\epsilon}$ where you need to parallel transport $P(\epsilon)$ to $P(0)$ Define the Christofel symbols as $\Gamma_{\beta\gamma}^{\alpha} = <\omega^{\alpha}, \nabla_{\gamma} e_{\beta}>$ A more convenient form of the symbols for a given basis in terms of the metric is: $\Gamma_{\mu\beta\gamma} = \frac{1}{2}(g_{\mu\beta,\gamma}+g_{\mu\gamma,\beta}-g_{\beta\gamma,\mu})$ The geodesic equation is defined as parallel transporting the tangent vector onto itself: $\nabla_{u} u = 0$ In coordinates: $\frac{d^{2}x^{\alpha}}{d\lambda^{2}}+\Gamma_{\mu\gamma}^{\alpha} \frac{dx^{\mu}}{d\lambda}\frac{dx^{\gamma}}{d\lambda}=0$ In component form, we can define the components of $\nabla_{u} T$ as $\frac{D T_{\alpha}^{\beta}}{d \lambda} = T_{\alpha;\gamma}^{\beta}u^{\gamma} = \frac{dT_{\alpha}^{\beta}}{d\lambda}+(\Gamma_{\mu\gamma}^{\beta} T_{\alpha}^{\mu}-\Gamma_{\alpha\gamma}^{\mu}T^{\beta}_{\mu}) \frac{dx^{\gamma}}{d\lambda}$ Differential Topology You don’t need a metric to describe curvature; it just makes calculation useful Many of our prior definitions still hold without a metric Events are points in spacetime, irrespective of the underlying geometry Curves are still paths through spacetime; without a metric, you can’t ascribe a proper length to the path though The parameter $\lambda$ is arbitrary up to some origin and unit change (ie. the transformation $\lambda’ = a\lambda +b$ traces out the same path in spacetime Vectors in flat spacetime could be though of as the difference between two events In curved coordinates, this doesn’t make sense Vectors instead become the directional derivative along some curve: $u = \partial_{u} = \frac{d}{d\lambda}$ along the curve The directional derivatives for a tangent space, which is isomorphic to the tangent vectors you normally think of Another useful construction is to embed your manifold in some higher dimensional flat space, take some curve that is tangent to your point of interest, and take the limit as $\lambda=0$ to arrive at the tangent vector on the manifold This is a bit of a hack, it that you shouldn’t need to refer to some higher dimensional space to extract properties of the manifold, but it works for computation In the tangent space at event $P_{0}$, we can define a coordinate basis as $e_{\alpha} = \frac{\partial}{\partial x^{\alpha}}$ This let’s us write any vector as $v^{\nu}e_{(\nu)}$ If we are transforming between coordinate bases, then we know that $\frac{\partial}{\partial x^{\alpha}} = \frac{\partial x^{\beta}}{\partial x^{\alpha’}}\frac{\partial}{\partial x^{\beta}}$ In general, we have some non-singular transformation matrix $L^{\beta}_{\alpha}$ One forms also reside in the tangent space A set of one-forms $w^{\beta}$ obeys the identity $<w^{\beta},e_{\alpha}> \delta^{\beta}_{\alpha}$ This doesn’t require the metric to be interpretable! Just think of this as demanding that each $e_{\alpha}$ is parallel to the associated $w^{\alpha}$ and is perpendicular to every other one-form In flat spacetime, we know that $<u,v> = u\cdot v$ Without the metric, tensors can’t change their rank. They are the same otherwise Since vectors can be mapped to directional derivatives, one can ask of the order you act these operators on a function matter. Define $[u,v] f = u(v(f))-v(u(f))$ as the commutator In a coordinate basis (re: $u = \frac{\partial}{\partial x^{a}}$ and $v = \frac{\partial}{\partial x^{b}}$), then the commutator is 0 due to partial derivatives commuting In general, this does not hold in arbitrary coordinate systems Affine Geometry Remember that events and trajectories are independent of the underlying spacetime structure How you trace the curve (ie. the $\lambda$ parameter) is not unique You can “start your clock” at some arbitrary point You can change the units of your timescale (years to seconds for example) Combining the above lets you define a more general affine parameter $\alpha = a\lambda +b$ Given some initial event location, with some initial velocity, the worldline of the particle will be unique Up to the affine reparameterization of course In order to define the length of a curve, you need a metric Parallel Transport GR dictates that two test bodies, initially falling on parallel neighboring geodesics, get pushed towards each other by spacetime curvature. What does it mean for two geodesics to be parallel? Imagine that you have two points along some curve $r(\lambda)$. Two points along that curve A and B reside in their own tangent spaces. There is no way a priori to relate a vector in the tangent space of A to a vector in the tangent space of B What you can do is define some local coordinate system at each point along the curve which obeys Minkowski spacetime Take the tangent vector in A that you want to compare with the tangent vector in B (both of them lying along the same curve) As you move along the curve and transition between each coordinate system, keep the tangent vector unchanged as seen in each coordinate system Once you arrive at B, compare the two vectors. If they are the same, then you are in a flat spacetime. Otherwise, you have curvature! This follows from the equivalence principle (Minkowski spacetime is valid locally!) Schild’s ladder let’s you visualize this process: Imagine that you have some curve through A and B (need not be a geodesic) Take some tangent vector and propagate along that vector to some nearby point X Define some nearby point R which lies along the curve from A to B Take the geodesic from X to R with some affine parameterization $\lambda$. Define the point Y along XR at the point $\lambda_{Y} = \frac{1}{2}(\lambda_{x}+\lambda_{R})$ Take the geodesic from A to Y with some affine parameter $\lambda$. You need increase $\lambda$ by $\lambda_{Y}$. To arrive at Z, let $\lambda_{Z} = 2*\lambda_{X}$ The curve RP gives vector AX “parallel transported” along the curve AB Rinse and repeat with R being the new A and P being the new X In flat spacetime, AX and RP are parallel, and all vectors here are local, so this holds in curved spacetime (equivalence principle) Ask how rapidly a vector field $\vec{v}$ is changing along a curve with tangent vector $\vec{u} = \frac{d}{d\lambda}$ Call this $\nabla_{u} v$ (ie. the covariant derivative of v along u) Construct this as follows: Take the geodesic along the tangent vector $\vec{u}$ originating at some point defined at $\lambda_{0}$ Take the vector from the vector field defined at the point $\lambda = \lambda_{0}+\epsilon$ along the geodesic Parallel transport v back to $\lambda_{0}$. Define this parallel transported vector as $\vec{w}(\lambda_{0}+\epsilon)$ Define $\nabla_{\vec{u}} \vec{v} = lim_{\epsilon\rightarrow 0} \frac{\vec{w}(\lambda_{0}+\epsilon)-\vec{v}(\lambda_{0})}{\epsilon}$ Schild’s ladder let’s you pictorally show: Symmetry: $\nabla_{u} v \nabla_{v} u = [u,v]$ chain rule: $\nabla_{u}(f \vec{v}) = f \nabla_{u} v+ v \partial_{u} f$ addition: $\nabla_{u}(v+w) = \nabla_{u} v + \nabla_{u} w$ From linearity: $\nabla_{a\vec{u}+b\vec{n}} \vec{v} = a \nabla_{u} v + b \nabla_{n} v$ a and b are some scalar functions (or numbers) Given the above, we can describe parallel transport as: $\frac{d\vec{v}}{d\lambda} = \nabla_{u} \vec{v} = 0$ implies that the vector field V is parallel transported along the vector $u = \frac{d}{d\lambda}$ You can test if a curve is geodesic if it’s tangent vectors covariant derivative along itself vanishes: $\nabla_{u} u = 0$ $\nabla$ is not a tensor! Parallel Transport: Component Notation Since each point has it’s own tangent space, the basis vectors and dual basis vectors change from point to point To quantify how these basis change, use the covariant derivative. Define $\nabla_{e_{\beta}} = \nabla_{beta}$ to denote the covariant derivative along some basis vector $e_{\beta}$ This rate of change is a vector, so it can be expanded in the vector basis: $\nabla_{beta} e_{\alpha} = e_{\mu} \Gamma^{\mu}_{\alpha,beta}$ $\Lambda_{\alpha\beta}^{\mu} = (\omega^{\mu}, \nabla_{\beta} e_{\alpha})$ Similarly, you can show that $\nabla_{\beta} \omega^{\nu} = -\Gamma^{\nu}_{\alpha\beta} = \omega^{\alpha}$ This implies that $(\nabla_{\beta} \omega^{\nu}, e_{\alpha}) = -\Gamma^{\nu}_{\alpha\beta}$ This can be easily generalized to covariant derivatives of tensors: You have the standard partial derivative part You add terms of the form $S*\Lambda$ For an upper index, contract that upper index of S with the non-differentiating lower index in $\Lambda$ For a lower index, contract that index with the upper index of $\Lambda$, and preprend a minus sign As an example, let $S_{\beta\gamma}^{\alpha}$. Denote $\nabla_{\delta} S_{\beta\gamma}^{\alpha} = S^{\alpha}_{\beta\gamma;\delta}$ $ S_{\beta\gamma;\delta}^{\alpha} = S_{\beta\gamma,\delta}^{\alpha} + S_{\beta\gamma}^{\mu} \Lambda_{\mu \delta}^{\alpha} - S_{\mu\gamma}^{\alpha} \Lambda_{\beta \delta}^{\mu} - S_{\beta\mu}^{\alpha} \Lambda_{\gamma \delta}^{\mu}$ Make sure that you contract on the non-differentiating index! Spacetime Curvature Relative Acceleration of Neighboring Geodesics Imagine a family of geodesics. You can distinguish between two geodesics by a continuous selector parameter n Hence, you can describe a point via two parameters: $\lambda$ is your affine parameter for each geodesic $n$ determines which geodesic you are on From this, we can define: The tangent vector $\frac{\partial}{\partial \lambda}$ The separation vector $\frac{\partial}{\partial n}$ which denotes the separation between some fiducial geodesic at n and some typical nearby geodesic $n+\Delta n$ geodesic deviation means that: You parallel transport some seperation vector n along the fiducial geodesic This tip of this vector traces out some cannonical trajectory The actual trajectory of a test geodesic might deviate from this cannonical path Imagine that you displace the separation vector along the fiducial geodesic by $+\frac{\lambda}{2}$ and $-\frac{\lambda}{2}$ Calculate the covariant derivative at these two displacements (ie. $\nabla_{u} n$) Move these derivative vectors back to $\lambda$ via parallel transport and subtract them after dividing through by $\delta \lambda$. The resulting 2nd covariant derivative is thus $\nabla_{u}\nabla_{u} n$ This is the “relative acceleration” between the geodesics Now imagine that you go around in a “square” loop: You parallel transport the seperation vector by $\Delta \lambda$, the you parallel transport to an adjacent geodesic at $n+\Delta n$, parallel transport back along this new geodesic by $-\delta \lambda$, then move back to the start by parallel transporting by $-\Delta n$ This gives you that $\nabla_{u} \nabla_{u} n + [\nabla_{n},\nabla_{u}] u = 0$ This follows since when you “parallel transport”, you use the geodesic equation $\nabla_{u} a = 0$ where u is the tangent vector of the curve. We are also parallel transporting the initial tangent vector u. $[\nabla_{n},\nabla_{u}]$ is almost a valid curvature operator The problem is that the above depends on the derivatives at the point of evaluation. The actual definition of the curvature tensor is $R(A,B) = [\nabla_{A}, \nabla_{B}]-\nabla_{[A,B]}$ $\nabla_{[A,B]}$ is the derivative along the vector $[A,B]$ This gives the Riemann curvature tensor as $R(\sigma, C,A,B) = <\sigma, R(A,B)C>$ This is a type 1,3 tensor (takes in a one-form $\sigma$ and 3 vectors and spits out a number) In a coordinate basis, we have $R_{\beta\gamma\delta}^{\alpha} = \frac{\partial \Gamma_{\beta\delta}^{\alpha}}{\partial x^{\gamma}} - \frac{\partial \Gamma_{\beta\gamma}^{\alpha}}{\partial x^{\delta}}+\Gamma_{\mu\gamma}^{\alpha}\Gamma_{\beta\delta}^{\mu}-\Gamma_{\mu\delta}^{\alpha}\Gamma_{\beta\gamma}^{\mu}$ A zero Reimann tensor implies a flat space (re: geodesics are straight lines) Reimann Normal Coordinates In a curved spacetime, you can’t define a global coordinate system such that $\Lambda_{\beta\gamma}^{\alpha} = 0$ everywhere To define a local inertial system, use Riemann normal coordinates Pick a point P in your manifold, and define an orthonormal basis at this point Radiate all the geodesics out from point P. Each geodesic has a tangent vector v at point P associated to it This overcounts by a lot Suppose that you step $\lambda$ along a geodesic with tangent vector v you reach the same point if you step $\frac{1}{2}\lambda$ along a geodesic with tangent vector v This means that if we fix $\lambda$ and vary v, we can reach the same points as all possible geodesics Find a tangent vector v at $P_{0}$ for which $P = \mathbf{P}(1,v)$ Expand v in terms of the basis $e_{a}$: $P = \mathbf{P}(1; x^{a}e_{a})$ From the above construction, we can see that: $e_{a}(\mathbf{P_{0}}) = \frac{\partial}{\partial x^{a}}_{P}$ $\Gamma_{\beta\gamma}^{\alpha}(P_{0}) = 0$ $\Gamma_{\beta\gamma,\mu}^{\alpha}(P_{0}) = -\frac{1}{3}(R_{\beta\gamma\mu}^{\alpha}+R_{\gamma\beta\mu}^{\alpha})$ Reimannian Geometry In order to define a notion of distance between points on a manifold, we define a metric This is a rank 0-2 tensor that is symmetric in it’s indices Explicitly, this is $g = g_{ij} dx^{i} \otimes dx^{j}$ Define $g^{ab}$ to be the inverse ie. $g^{\alpha\mu}g_{\mu\beta} = g_{\beta}^{\alpha} = \delta_{\beta}^{\alpha}$ With the metric, you can define a one-to-one mapping between the cotangent and the tangent space Use $g_{ab}$ to convert (1,0) tensors to (0,1) tensors Use $g^{ab}$ to convert (0,1) tensors to (1,0) tensors To specialize to general relativity, we demand that at every point on the manifold, you need to be able to define an orthonormal basis such that $g_{ab} = e_{a}\cdot e_{b} = \eta_{ab}$ We also want the metric to be as “flat” as possible: $g_{ab,y}(P_{0}) = 0$ We won’t be able to set 2nd derivatives to 0 generically though In light of the above, you can also derive the geodesic equations from extermizing the action $S = \int (g_{\mu\nu}u^{\mu}u^{\nu})^{\frac{1}{2}} d\lambda$ You do a trick by reparameterizing the curve in terms of the differential of the action itself: $dS = (g_{\mu\nu}u^{\mu}u^{\nu})^{\frac{1}{2}} d\lambda$ You get back the geodesic equation via Eular-Lagrange It’s useful to calculate Christoffel coefficients Metric Induced Symmetries of Reimann Tensor $R_{\alpha\beta\gamma\delta}$ is antisymmetric between the last two indices It is symmetric under swapping adjacent pairs of indices: $R_{\alpha\beta\gamma\delta} = R_{\gamma\delta\alpha\beta}$ $R_{\beta\gamma\delta}^{\alpha}$ satisfies a Bianchi identity: $R_{\beta\gamma\delta}^{\alpha}+R_{\gamma\delta\beta}^{\alpha}+R_{\delta\beta\gamma}^{\alpha} = 0$ You can notate this as $r_{[\beta\gamma\delta]}^{\alpha} = 0$ Another Bianchi identity is $R_{\beta[\gamma\delta;\epsilon]}^{\alpha} = 0$ Reminder that ; denotes a covariant derivative w.r.t. the index to the right of it All the above implies that $R_{[\alpha\beta\gamma\delta]} = 0$ Bianchi Identities From a physics POV, we demand the divergence of the stress energy tensor to be zero since this encodes a conservation law Alternatively, we can say that the exterior derivative of T vanishes What tensorlike feature of the geometry of spacetime can provide this conservation? (ie. What tensor derived from curvature objects has an exterior derivative of 0?) If this exists, we can make it proportional to the stress-energy tensor The metric is a set of 10 potentials which can be combined into a single tensor From the metric, you can define the curvature operator $R =\frac{1}{4} e_{\mu} \wedge e_{nu} R_{\alpha\beta}^{\mu\nu} dx^{\alpha} \wedge dx^{\beta}$ We know that R obeys two Bianchi Identities as listed above The differential identity ( the 2nd one), is a manifestation of the general principle “the boundary of the boundary is zero” Boundary of Boundary is 0 Look at a 3D case. Define some tensor field at each point in space (for clarity, this can be a scalar field or a vector field) Imagine that you sum this quantity over all of the edges of the cube, where each edge has some orientation to it There are no “exposed” edges . Namely, every edge has a corresponding edge on another face with the opposite orientation Summing across edges gives you zero due to this pairing of opposite orientation edges The idea generalizes to higher dimensions. Each volume Has a set of oriented hypersurfaces The “edges” of the hypersurfaces cancel out in pairs Moment of Rotation Imagine a 3-volume. For simplicity, take it to be a cube with cartesian coordinates The rotation associated with the front face $\delta y \delta z$ is $e_{\lambda} \wedge e_{\mu} {R^{|\lambda\mu|}}_{yz} \delta y \delta z$ We are in Reimann normal coordinates, and we demand $\lambda < \mu$ as denoted by $|\lambda \mu|$ Follows from parallel transporting a test vector A around the surface edge and seeing how much it differs from the start The moment of rotation is taking $P_{cff}-P$ and $\wedge$ prior to the rotation $P_{cff}$ is the center of the front face P is some arbitrary point in space that you are comparing the moment against The exact choice does not matter since the Bianchi identity will make it irrelevant only the back face is subtracted The difference between the two is a well defined vector, but the actual locations of the two points are a bit nebulous $P_{cff}$ and P need to be infinitesimally close though Do a similar calculation on the back face. We have that the difference between the two faces is $P_{cff} - P_{cbf} = \Delta x e_{x}$ Sum across all 6 faces to get $e_{\nu} \wedge e_{\lambda} \wedge e_{\mu} {R^{|\lambda \mu|}}_{|\alpha\beta|} dx^{\nu}\wedge dx^{\alpha} \wedge dx^{\beta}$ This is evaluated on the cube $e_{x} \wedge e_{y} \wedge e_{z} \Delta x \Delta y \Delta z$ Alternatively, we can write this as $dP \wedge R$ $dP = e_{\sigma} dx^{\sigma}$ $R = \frac{1}{4} e_{\mu} \wedge e_{\nu} {R^{\mu\nu}}_{\alpha\beta} dx^{\alpha}\wedge dx^{\beta}$ $dR =0$ is the 2nd Bianchi identity $dP \wedge R$ is a trivector. Taking the dual of this would just give a 1-form, which is equivalent and easier to work with In general, we can see that $*(e_{\nu} \wedge e_{\lambda} \wedge e_{\mu}) = {\epsilon_{\nu\lambda\mu}}^{\sigma} e_{\sigma}$ We also se that $dx^{\nu}\wedge dx^{\alpha} \wedge dx^{\beta} = \epsilon^{\nu\alpha\beta\tau} d^{3} \Sigma_{\tau}$ $d^{3}\Sigma_{\tau}$ is the 3 volume element Use the above to rewrite the moments across the entire cube as : $*(dP \wedge R) = e_{\sigma} {X_{\nu}}^{\sigma\nu\tau} d^{3} \Sigma_{\tau}$ $X_{\nu} = (*R*)_{\nu}$ ${X_{\nu}}^{\sigma\nu\tau}$ is the Einstein tensor it relates the initial starting volume with the final moment Observe that $d*G = 0$ (follows from dd = 0 and that * and d commute) Hence, *G is a candidate object whose exterior derivative vanishes The Field equations use G as the source of curvature because of this Einstein Field Equations The main idea is that the stress energy tensor should give rise to the geometry of the system and vice-versa. Written abstractly, we have that $G = \kappa T$ where T is the stress energy tensor, $\kappa$ is some proportionality constant, and G is a purely geometric tensor G must satisfy the following properties: G should vanish when the spacetime is flat G needs to purely be a function of the metric, and other tensors that can be contructed from the metric (like the Reimann tensor) G needs to be linear in Reimann (simplest non-trivial assumption. Could include higher order terms) G needs to be symmetric G needs to have a vanishing divergence $\nabla \cdot G = 0$ The only tensor which satisfies this is $G_{\mu\nu} = R_{\mu\nu}-\frac{R}{2}g_{\mu\nu}$ This can be though of as $G_{\mu\nu} = G_{\mu\alpha\nu}^{\alpha}$ where the rank (1,3) tensor is the double dual of the Reimann tensor This can be proven by making the general ansatz $G = aR_{\alpha\beta}+bRg_{\alpha\beta}+\Lambda g_{\alpha\beta}$, then applying each condition you want to satisfy To calculate kappa, the above equation should reduce to the Newtonian limit Assume that you are in a region of constant density This implies that the $T_{00}$ component dominates all other components of the tensor. Also, the tensor is diagonal In the weak field approximation, we have $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$ where the magnitude of $h_{\mu\nu}$ is small Examine the 00 component and grind through the algebra gives you a Poisson equation. Compare this against the actual Poisson equation ($\nabla\cdot\frac{F_{g}}{m} = -4\pi \rho$ )to see that $\kappa = 8\pi$ Hence we have $G_{\mu\nu} = 8\pi T_{\mu\nu}$ EFE Applications Thermodynamics Assume that we are dealing with equillibrium thermo of a perfect fluid with fixed chemical compositions We can charaterize the system via a set of potentials (n, $\rho$, P, T, s ,$\mu$) Potentials are frame independent functions (since they are scalars) However, we can evaluate the potentials in the fluid’s rest frame The potentials, evaulated in the rest frame, are: n: baryon number density: number of baryons per unit 3d volume of the rest frame. Anti-baryons are negative $\rho$: density of the total mass-energy rest frame P : isotropic pressure in rest frame T: temperature s: entropy per baryon in rest frame $\mu$: chemical potential of baryons We are assuming: the chemical composition is fixed by the baryon density and entropy per baryon We are assuming that chemical reaction rates are two slow to matter, or that they are so fast that equillibrium is reached We assume baryon conservation: $\nabla \cdot S = 0$ where $S = nu^{\mu}$ $\frac{d}{d\tau}(nV) = 0$ The 2nd law states entropy can be generated, but not destroyed Assume that entropy very slowly flows in and out of a fluid element $ $\frac{ds}{d\tau} \geq 0$ The 1st law states: In words: The differential of the energy in a volume element containing a fixed number A of baryons equals the negative pressure times the volume differential plus the temperature times the entropy differential $d(\frac{\rho A}{n}) = -P d(\frac{A}{n}) + T d(As) \rightarrow d\rho = \frac{\rho+P}{n} dn +nT ds$ “d” is a exterior derivative Can think of $\rho = \rho(n,s)$. We can then think of the first law as a total differential to state: $\frac{\rho+P}{n} = \frac{\partial \rho}{\partial n}_{s}$ and similar for entropy The two partials must obey the Maxwell relationship $(\frac{\partial P}{\partial s}) = n^{2} \frac{\partial T}{\partial n}_{s}$) for consistency The chemical potential is defined as follows: Take a sample of the simple fluid in a fixed thermodynamic state (fix n and s) Take a much smaller sample containing $\delta A$ baryons in the same state as the larger sample (fix n and s) Inject the smaller sample into the larger sample, holding the volume fixed during the injection process The total mass energy, plus the work required to perform the injection is $\mu \delta A$ $\mu \delta A = \rho(\frac{\delta A}{n})+ P \frac{\delta A}{n} = \frac{\partial \rho}{\partial n}_{s}$ Hydrodynamics is just enforcing that $\nabla \cdot T=0$ E&M The are unmodified from SR! Just change your metric. Recall that $F^{0j} = E^{j}$ and $F^{jk} = \epsilon^{jkl} B^{l}$ We have that: $F^{\alpha\beta}_{;\beta} = 4\pi J^{\alpha}$ $F_{\alpha\beta;\gamma} + F_{\beta\gamma;\alpha} + F_{\gamma \alpha; \beta} = 0$ $ma^{\alpha} = F^{\alpha\beta}q u_{\beta}$ Geometric Optics Define the following length scales : $\lambda$ is the wavelength of the wave involved L is the length over which the wave parameters vary R is the radiius of curvature through which the waves propagate Geometric optics is the regime where $\lambda «L$ and $\lambda « R$ Hence, we can assume that waves are locally plane waves The main results from these assumptions are that: Light rays are null geodesics The polarization vector is perpendicular to the rays and the parallel-propagated along the rays The amplitude is governed by an adiabatic invariant (ie. Number of photons is conserved)

Scales Black holes are described by 3 parameters: Mass, spin, and charge In practice, astrophysical black holes have no charge, since they are embedded in a plasma which “shorts out” the hole We can define some characteristic scales of black holes $R_{s} = \frac{2GM}{c^{2}}$ is the Schwartzchild radius $\Delta t_{LC} = \frac{R_{s}}{c}$ is the light crossing time (roughly the time needed to cross the black hole, up to a scale factor) How bright can be a black hole be? We would want to convert all of the black hole’s mass into radiation in the relavent time scale $L_{max} = \frac{Mc^{2}}{\frac{2GM}{c^{3}}} \approx \frac{c^{5}}{G} \approx 10^{59} \frac{erg}{s}$ For reference, supernovae are around $10^{52}$ What is the luminosity of the observable universe (excluding the blackholes)? This is Fermi estimation problem: Find the mass of the sum, calculate the mass to energy conversion, and calculate the time needed for a photon to escape the center The escape time is roughly a billion years, which gives a luminosity of around $10^{37}$ There are roughly $10^{10}$ stars in the universe. Assume all stars are like the sun There are also roughly $10^{10}$ galaxies in the universe Combining these gives that the luminosity of the universe is $10^{53}$ The point of this is that gravity is very efficient at converting mass to energy If the black hole is moving close to the speed of light towards you (pileup of photons), then the observed luminosity can exceed this threshold Gamma ray bursts (GRBs) are the prime example of this boosted luminosity Eddington Limit Imagine that you have a star which is gravitationally bound and emits radiation. These two forces oppose each other: gravity attracts to the center, while radiation pressure pushes the star boundary outwards Larger stars have more mass, but they also tend to produce more radiation pressure The Eddington limit is when these two forces balance each other Define $\dot{m}$ as the accretion rate (gas falling onto a compact object or star) The luminosity is $L = \frac{GM\dot{m}}{r}$, where $\frac{M}{r}$ is the compactness The associated Eddington luminosity is an upper bound on the brightness of these objects (there are exceptions though) Let’s look at the simplest case: We have a fully ionized Hydrogen gas (plasma) accreting isotropically onto a star/ compact object Electrons getting accelerated around protons produces radiation, which creates the outward pressure Thomson Scattering Imagine a EM wave propagating along the $\hat{z}$ direction. This will oscillate an electron in a direction transverse to z (fix this direction to be x) This oscillating generates the radiation Newton’s 2nd Law, and letting $E(z,t) = exp(i\omega t) \hat{x}$, the resulting motion is $x(t) = A exp(i\omega t)$ where $A = \frac{qE(z)}{m\omega^{2}}$ The cross section scales like $A^{2}$, which means that the cross section of the proton is roughly 1 million times smaller than the electron $\frac{\sigma_{Tp}}{\sigma_{Te}} = \frac{m_{e}^{2}}{m_{p}^{2}}$ The radiation force is defined as the energy flux times the cross section divided by the speed of light $F = \frac{L}{4\pi r^{2}} \rightarrow f_{rad} = \frac{\sigma_{T} L }{4\pi r^{2} c}$ Setting the gravitational force (ignoring the electron mass) equal to the radiative force (ignoring the proton cross section) gives the target luminosity (re: Eddington Luminoosity) as $L_{edd} = \frac{4\pi c G M m_{p}}{\sigma_{T}}$ Setting the mass scale to that of the sun, we get that $L_{Edd} = 1.3E38 (\frac{M}{M_{*}}) \frac{ergs}{s}$ For a black hole, $L_{edd} = 1.5E45 M_{7} \frac{erg}{s}$ where $M_{7} = \frac{M}{1E7 M_{*}}$ If we assume that the object is a perfect blackbody, we can define an effective temperature $L = 4\pi R^{2} \sigma_{B} T_{eff}^{4}$ Variability The acceleration is $\frac{GM}{r^2}$. Suppose that the radiation pressure suddenly vanished. How long would it take for the star to collapse? Naively, one would say that $\frac{1}{2} a t_{dyn}^{2} = R \rightarrow t_{dyn} = \sqrt{\frac{2R^{3}}{GM}}$ Assuming that force at star radius is the same for all smaller radius This is roughly $t_{dyn} = \frac{1}{\sqrt{G\rho}}$ $t_{dyn}$ sets the timescale for how quickly the luminosity of the star can change Sense of Scale The above outlines the rough energies, wavelengths, frequencies and temperatures for each band of radiation Atmospheric Attenuation Our atomsphere absorbs different amounts of radiation depending on the frequency This change is absorption is described by the opacity of the sky The above sets constraints on where expeiments/telescopes can be (radio and visible can be ground based, but other bands might demand high altitudes or be in space) Intensities Define $I_{\nu} = \frac{dE}{dA dt d\nu d\Omega}$ as the frequency dependent brightness of a source This is constant along rays in free space Imagine that you have 2 detectors which are coaxial with seperation R and area $dA_{1}$ and $dA_{2}$ respectively The energy change across each is $dE_{1} = I_{\nu_{1}} dA_{1} dt d\Omega_{1} d\nu_{1}$. Same for the 2nd detector Solid angle is related to area via $d\Omega_{1} = \frac{dA_{2}}{R^{2}}$ The frequency and time crossing are the same Use conservation of energy to see that the brightness is conserved Define the flux as $F_{\nu} = \int I_{\nu} \cos \theta d\Omega$ cosine modulation to account for when you don’t look at a source head on You can recover the inverse square law by integrating over a small solid angle Define the pressure as $P_{\nu} = \frac{1}{c} \int I_{\nu} \cos^{2}(\theta) d\Omega$ Additional cosine arises from projecting onto the area of interest Define the energy density $u_{\nu}$ as $dE = u_{\nu}(\Omega) dV d\Omega d\nu$, or in terms of the brightness: $u_{\nu} = \frac{I_{\nu}}{c}$ Define the mean intensity as $J_{\nu} = \frac{1}{4\pi} \int I_{\nu} d\Omega$ Particle acceleration Cosmic Rays There exists cosmic rays which are protons of energies 1E21 eV The LHC boosts protons to 1E13 eV When they hit the atmosphere, they create showers of high energy particles, which produce Cherenkov radiation cones in the air There is a hierachy of densities in our Universe: On earth, the density of air is roughly 1E20 per cc In the ISM, there is a density of 1E0 particles per cc In the IGM, there is a density of 1E-6 per cc In the ISM and IGM, the mean free path $\lambda = \frac{1}{n \sigma}$ is very long, which lets particles get accelerated for a long time before getting scattered Fermi Acceleration There are these objects called molecular clouds These are objects in the ISM which form in slightly over-dense regions They are called “star nurseries” Point being, despite being made up of a bunch of particles, collectively the clouds have some average velocity in some direction There are a lot of these molecular clouds scattered around. Fermi realized that a particle could bounce off of these clouds and potentially gain these massive energies This hinges on the fact that the velocity distribution of the clouds is assymetric A non-relativistic calculation of this phenomena: Let $\Delta$ be the distance from the particle to the nearest cloud Let $v$ and $V$ be the velocities of the particle and the cloud respectively Let’s calculate the probability of a collision The time until a head-on collision is $\tau_{u} = \frac{\Delta}{v+V}$ The velocities point in opposite directions The time until a catch-up collision is $\tau_{c} = \frac{\Delta}{v-V}$ The velocities point in same direction The probability of a head-on collision is $\frac{\frac{1}{\tau_{u}}}{\frac{1}{\tau_{u}}+\frac{1}{\tau_{c}}} = \frac{v+V}{2v}$ The probability of a catch-up collision is $P_{c} = \frac{v-V}{2v}$ What is the change in energy for each case? Assume that the mass of the cloud is much larger than the particle. This is just an elastic collision For a head-on collision: velocity after the collision is $v+2V$, hence $\Delta E_{h} = \frac{m}{2} (v+2V)^{2} -\frac{m}{2} v^{2}$ For a catch-up collision: velocity after the collision is $-v+2V$ What is the expected change is energy? Its $<\Delta E > = \Sigma_{i} P_{i} E_{i}$ Doing that calculation out gives $\frac{<\Delta E >}{E} = 4 (\frac{V}{v})^{2}$ This is a 2nd order process (b/c of the -2), which means it’s pretty inefficient This Fermi acceleration is a general concept. There are “converging magnetic flows” where the “clouds” are always moving towards you. Hence you get a more efficient energy gain Let $\frac{dE}{dt} \propto \nu < \Delta E>$, where $\nu$ is the frequency of collisions For the molecular clouds, we see that $\frac{dE}{dt} = \alpha E$ for some characteristic $\alpha$ of the ISM Let’s think of this problem in phase space, where our axes are E and x $\frac{d N}{d t} = -\frac{\partial \phi_{s}}{\partial x} - \frac{\partial \phi_{e}}{\partial E} + \frac{\partial N}{\partial t}$ Think Stokes Theorem, where the flow on the boundary and the production of particles determines the change in the number of particles in some volume The spatial flux is defined by $\phi_{s} = -D \frac{\partial N}{\partial x}$, where D is some diffusion constant The energy flux is $\phi_{e} = \frac{N \Delta E}{\Delta t}$ We define some timescale $\tau$, which dictates how quickly particles leave the system (for instance, the Larmor radius grows so large so as to leave the system). This then gives that $\frac{\partial N}{\partial t} \approx -\frac{N}{\tau}$ For molecular clouds, D=0, and the system is in equillibrium, so $\frac{dN}{dt} = 0$ Solving for N yields $N(E) \approx N_{0} exp(-(1+\frac{1}{\alpha}{\tau}))$ where $\alpha = 4 \nu (\frac{V}{v})^{2}$ In log-log space, this is a straight line (re: a power law) MHD (Magneto-hydrodynamics) We have Maxwell’s equations $\nabla \cdot E = 4\pi \rho_{e}$ $\nabla \cdot B = 0$ $\nabla \times E = -\frac{1}{c} \frac{\partial B}{\partial t}$ $\nabla \times B = \frac{4\pi}{c} J \frac{1}{c} \frac{\partial E}{\partial t}$ From Maxwell, we can derive charge conservation: $\frac{\partial \rho_{e}}{\partial t} = \nabla \cdot J$ We also make the ideal assumption (re: no resistivity) This gives that there is no Lorentz force in the plasma (re: if there was a force, then due to no resistivity, the charges would rearrange themselves and short out the force) Hence we have $\vec{E} \approx - \frac{\vec{v}}{c} \times \vec{B}$ We add on the additional constraint of conservation of mass: $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$ Writing down momentum conservation gives: $\rho \frac{d\vec{v}}{dt} = \rho \frac{\partial \vec{v}}{\partial t} + \rho (\vec{v} \cdot \nabla) \vec{v} = -\frac{\vec{B}}{4\pi} (\nabla \times \vec{B} - \frac{1}{c} \frac{\partial \vec{E}}{\partial t})$ In non-relativistic MHD, we can drop the $\frac{\partial E}{\partial t}$ term Suppose that you have some B field along the z axis, which experiences some sinusoidal peturbation in the x axis: $\vec{B} = B_{0} \hat{z} + B_{A} exp(ikz- i\omega t) \hat{x}$ This could arise via neutron star quakes, or acretion disk turbulances These quakes have some length scale $\lambda$, which modules the perturbation by $\sin(k_{y} y)$ The E field which gets created is $\vec{E} = \frac{ic}{\omega}(ikB_{a} \sin(k_{y}y) \hat{y} - k_{y} B_{A} \cos (k_{y} y) \hat{z}) exp(ikz-i\omega t)$ So this peturbation creates an E field along the magnetic field line, which can accelerate particles to relativistic speeds This speed up is reduced due to synchrotron radiation Acretion Disks We have gas around a object with a strong gravitational field Gravity can accelrate the charges, causing energy loses via radiation; angular momentum is preserved though Radially, The balance of centrifugal acceleration and gravity causes the gas to gather into a disk Vertically, the gravity gets balanced by the radiation pressure Thin Disks A thin disk is defined as being thin in the z-axis (due to being cold and having low pressure) The radial pressure gradient is negligible A thicker disk is hotter (hence larger pressure gradients) How does all the gas fall into the black hole when the portions of the disk have angular momentum? There is a net flux of angular momentum out of the disk Imagine a ring of gas of radial thickness dR at radius R. The angular momentum on the inside is $R\Omega(R)$, where $\Omega(R) = \frac{v_{\phi}}{r}$, and on the outside is $(R+dR)\Omega(R+dR)$ We imagine that there is some particle exchange from inside the disk to the outside, and vice versa The torque due to the outflow is $\tau = \dot{M_{1}} (R+dR) R \Omega(R)$ The torque due to the inflow is $\tau = \dot{M_{2}} R(R+dR) R \Omega(R+dR)$ In short, the particles keep their old angular momentum as they cross the boundary The inflow and outflow mass rate are the same ($\dot{M_{1}} = 2\pi R \rho(R) H \tilde{v}$). H is the height of the disk. $\rho(z)$ is the density. $\tilde{v}$ is the average radial velocity of the disk Hence, the net torque is $\tau = \dot{M_{1}} R(R+dR)(\Omega(R)-\Omega(R+dR)) = 2\pi R H \tilde{v} \rho(z) R(R+dR)(dR \Omega’(R))$ The kinematic viscosity is defined as $\nu = dR \tilde{v}$ (units of length squared over time) Hence: $\tau = -2\pi \nu \epsilon R^{3} \Omega'$ There is a characteristic length and velocity scales in the problem (height of disk and speed of sound). You can encode this in the kinematic viscosity via $\nu = \alpha c_{s} H$, where $0 \leq \alpha \leq 1$ Due to the viscosity of the gas, the friction generated causes the gas to radiate like a blackbody The mass flow of a ring is given by $R\frac{\partial \Epsilon}{\partial t} + \frac{\partial}{\partial R}(R\Epsilon v_{R}) = 0$ Imagine the mass flow at the inner and outer radii, then take the limit as dR goes to 0 $\Epsilon$ is the mass density of the ring, $v_{r}$ is the radial velocity Similarly, there is a flow of angular momentum: $R \frac{\partial}{\partial t}(\Epsilon R^{2} \Omega) + \frac{\partial}{\partial R}(R\Epsilon v_{r} R^{2}\Omega) = -\frac{1}{2\pi} \frac{\partial \tau_{out}}{\partial R}$ Same story: Take the difference between the inner and outer disks, and also add on the torques at the inner and outer rings Since the rings are Keplerian, we know that $\Omega = \frac{v}{R} \propto R^{-\frac{3}{2}}$ We can combine mass and angular momentum conservation, along with the Keplerian ring assumption to calculate the radial velocity: $v_{r} = -\frac{3}{\Epsilon R^{\frac{1}{2}}} \frac{\partial}{\partial R}(v \Epsilon R^{\frac{1}{2}})$ We can define $\dot{m} = 2\pi R \Epsilon (-v_{R})$ Thick Disks If the cooling is inefficient, or the mass accretion rate is too large, then disk has a tendency to become thick There are a lot of different regimes with funny names: ADAF: advection-dominated accretion flow RIAF: raidatively inefficient accretion flow CDAF: convection dominated accretion flow ADIOS: advection dominated inflow/outflow solution You get like a quadropole movement in the flow, where one axis is an inflow while the other axis is an outflow Magnetic fields become more important, and the ideal MHD approximation breaks down Analytic Assume that the angular velocity is now a function of z and R ($\Omega(z,R)$) For simplicity, assume that the radial velocity is much smaller than the angular velocity. the vertical velocity is also small Force balance gives $\frac{1}{\rho}\nabla P = -\nabla \Phi +\Omega^{2}\vec{R} = \vec{g_{eff}}$ To simplify the analysis, let’s examine isobaric surfaces $\nabla P = 0$ To further simplify the analysis, pretend that the z profile goes like $\pm R \tan \alpha $ The potential takes the form $\Phi = -\frac{GM}{(R^{2}+z^{2})^{\frac{3}{2}}}$ Assume that z « R Doing the algebra and calculus gives $\Omega^{2}(R) = \frac{GM \cos \alpha}{R^{3}}$ This reduces to a Keplerian disk in the small $\alpha$ limit Otherwise, the angular velocity is sub-Keplerian The radiation force becomes $\vec{F_{rad}} = -\frac{c}{\kappa \rho} \nabla P$ Integrating the force over the surface area of the disk gives you the luminosity: $L = \frac{4\pi c GM}{\kappa} -\frac{c}{\kappa} \int \nabla \cdot \Omega^{2} dV$ Doing the integral, you get $L = L_{edd}(1+\sin \alpha \frac{r_{2}}{r_{1}})$ This can be super-Eddington! This can explain how these super-massive black holes can be created in a Hubble time Acretion Columns Imagine that the central accretor has some non-zero magnetic dipole moment. The B field for a non-zero dipole is $\vec{B} = \frac{3\hat{n}(\hat{n}\cdot\hat{m})-\vec{m}}{r^{3}}$ In our case, $\hat{m}$ is along the rotation axis of the acretor, and $\hat{n}$ denotes the direction we are measuring the field at As a function of radius, assuming we are in the plane, we see that $B(R) = B_{0}(\frac{R_*}{R})^{3}$ The pressure of the gas is given by $\frac{\rho GM}{R}$ The magnetic pressure is given by $P = \frac{B^{2}}{2\mu_{0}}$ The magnetic pressure has dynamical consequences once these pressures are comparable to each other The radius where this interchange happens defines this transition on behavior The gas, once below this radius, will flow along the magnetic field lines towards the poles of the acretor. This forms an accretion column MHD We can describe our system as a function of pressure, velocity, and density at all x,y,z,t. These are not all conserved. We can construct a set of conserved variables from these: the total energy density (kinetic plus thermal), the momentum density and the density itself There are two pictures of flows that are useful The Eulerian description imagines that you fix a box in space. You then look at the flow of stuff through the faces of the box and any stuff which gets generated in the box The Lagrangian description follows a particular section of stuff around as a function of time and space With conservation of mass, we see that $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$ Can derive by looking at the mass flow in a little box is size $\delta V = \delta x \delta y \delta z$ With momentum conservation, we find that: $\frac{\partial}{\partial t}(\rho \vec{v}) + \nabla \cdot (\rho \vec{v}\vec{v}) = -\nabla \cdot \sigma + \vec{F_{b}}$ $\vec{v}\vec{v}$ and $\sigma$ are 3x3 matrices (re:tensors) You can split the stress tensor into diagonal pressure components and off-diagonal shear components You can dot the above with $\vec{v}$ to get an power conservation law The 1st law of thermodynamics states: $\frac{du}{dt} = T\frac{ds}{dt} + \frac{P}{\rho^{2}}\frac{d\rho}{dt}$ S is the entropy per unit mass Bondi Accretion This describes the accretion of cold gas onto a central mass where the flow is primary in the radial direction The interstellar medium (ISM) is such a cold gas Far away from the center accretor, the thermal energy density is given by $\rho c_{s}^{2}$. Define some radius $R_{a}$ where the thermal energy density equals the gravitational energy density Another relavent radii is $r_{s}$ (the sonic radius), where in infalling gas speed equals the sound speed Imagine that you have some steady radial flow: From mass conservation, we know that $\frac{1}{r^{2}} \frac{d}{dr}(r^{2}\rho v) = 0$ The mass accretion rate is then defined as $\dot{M} = 4\pi r^{2} \rho (-v)$ (minus sign since the gas is always moving inwards) From the momentum conservation, we have that $+\frac{\nabla P}{\rho}+ \frac{GM}{r^{3}} + v \frac{dv}{dr} = 0$ We assume an equation of state $P = D \rho^{\Gamma}$, which is representative of some ideal gas equation of state (or if $\Gamma=1$, it’s an isotherm) If we define $c_{s}^{2} = \frac{dP}{d\rho}$, do some chain rule to change $\frac{dP}{d r} = \frac{\partial P}{\partial \rho}\frac{d \rho}{dr}$, and incorporate the mass continuity equation, we see that: $\frac{1}{2}(1-\frac{c_{s}^{2}}{v^{2}}) \frac{d}{dr}(v^{2}) = -\frac{GM}{r^{2}}(1-\frac{2c_{s}^{2}r}{GM})$ This describes a transition from subsonic to supersonic speed Radii which this occurs at is when $v(r_{s}) = c_{s}(r_{s})$, which translates to $r_{s} =frac{GM}{2c_{s}^{2}(r_{s})$ We know the sound speed and the density at infinity. We need to relate these parameters to the values at $r_{s}$ Look at energy conservation: $\frac{1}{2}v^{2} + \int \frac{dP}{\rho} - \frac{GM}{r} = constant$ Assuming that $\Gamma \neq 1$, we know that v=0 at $r=\infty$. This let’s you write the sound speed as a function of the boundary conditions. Remembering that $v = c_{s}$ at $r_{s}$, we see that: $c_{s}(r_{s}) = c_{s}(\infty) (\frac{2}{5-3r})^{\frac{1}{2}}$ $\rho(r_{s}) = \rho(\infty) (\frac{c_{s}(r_{s})}{c_{s}(\infty)})^{\frac{2}{\Gamma-1}}$ Roche Lobes Most massive stars are in binary systems (~70%) We define the ratio between the two masses as $q = \frac{M_{2}}{M_{1}}$, where $M_{2} \leq M_{1}$ In an inertial frame, the potential is just the sum due to each mass In a co-rotating frame, you can peg the masses in place, and add in some additional potential to model this rotation For simplicity, we will only assume that the centrifugal force is relavent. Euler and Coreolis forces are small The centrifugal acceleration is $-\omega\ \times \omega \times r$, which means the potential is $\Phi(r) = -\frac{1}{2}(\vec{\omega}\times \vec{r})^{2}$ (just take the derivative w.r.t. r! to get back the force!) The relavent lengths of the system are: The binary separation a: related to the period via Kepler’s Law: $4\pi^{2}a^{3} = G(M_{1}+M_{2})P_{orb}^{2}$ $R_{i}$ the average radii of each mass $b_{i}$ the distance from the mass to the L1 Lagrange point The equipotentials of this system in the co-rotating frame looks like two lobes which meet at the Lagrange point L1 Some heutristics for how these parameters evolve are: $\frac{R_{2}}{a} 0.38+0.2 \log q$ $\frac{b_{1}}{a} = 0.5 -0.277 \log q$ Jets Imagine you have an electron moving through a constant magnetic field. Align your coordinate system along the electrons velocity, with the z axis orthogonal in the direction of the B field There is some pitch angle $\alpha$ which the electron makes with the B field the gyration radius is given by $r_{g} = \frac{v\sin\alpha}{\omega_{g}}$ where $\omega_{g} = \frac{eB}{m_{e}c}$ ( follows from Lorentz force law) The radiated power in the electron’s rest frame is then $P = \frac{2e^{2}}{3c^{3}}\omega_{g}^{2} v \sin^{2} \alpha$ (Larmor formula) The power emitted by the electron in some other frame is $P’ = \frac{4}{3} \sigma_{T} c \beta^{2} \gamma^{2} u_{B}$ $\sigma_{T}$ is the Thompson cross section (need Klein-Nishina if quantum is important) $U_{B} = \frac{B^{2}}{8\pi}$ You also need to average over the pitch angle: $\int 2\pi \sin^{3} \alpha d\alpha = \frac{8\pi}{3}$ Since energy and time both have the same transformation, then the above powers are the same This means that the energy gets relativistically beamed forwards As the electron whips around, the beam very quickly enters and leaves the line of sight of an observer Hence, the time we observe the beam is short, which means the frequency content is very large (hence why radio pulsars are a thing…) The width of the cone is roughly $\frac{1}{\gamma}$ The actual measured power depends on the velocity distribution of electrons

Following John Preskill’s 1998 Quantum Information and Computation notes. Introduction Landauer’s principle: It takes energy to erase information (since erasure always compresses phase space, such processes are irreversible) You can store a bit of information as one molecule in a box. If it’s on the left, it’s on, else it’s off. If you slowly compress the volume in half, you’re gaurenteed to be in the LHS. The change in entropy is k ln 2, which has some associated work that needs to be performed Logic gates used to perform computation are typically irreversible For instance, NAND is irreversible, since one bit of information is lost for each gate NOT is an example of a reversible gate Charles Bennett observed that any computation can in principle be done reversibly You can construct a Toffoli gate: Input is (a,b,c) Output is $(a,b, c \oplus (ab))$ So a and b get mirrored, and the third bit gets flipped if the first two bits are both 1 (otherwise, it mirrors the input) a Toffoli gate is universal (provided you throw away the ancilla bits when interpreting You can in principle do any computation up to the end, print out a copy of the answer (logically reversible process), then step back the computation back to the beginning Maxwell’s Demon The aforementioned ideas allows a resolution of the Maxwell’s demon paradox The original formulation is as follows: You have a partitioned box (split into A and B parts) and some demon which observes the molecules in the box If a fast particle is moving from A to B and it will cross the partition, then the demon allows it through If a fast particle is moving from B to A across the partition, the demon blocks it Over time, you will get faster particles on the right and slower ones on the left with minimal work. This has the effect of having heat flow from a cold place to a hot place at no cost, in violation of the 2nd law of thermodynamics The resolution to this is that the demon needs to have some memory/ keep information on the molecules in the box. If the demon has a finite memory capacity, eventually, information will need to be erased, which results in work being done Quantum Influence The quantum nature of reality changes the definition of information Quantum mechanics is a truly random process, which has no place in deterministic classical dynamics The uncertainty principle implies that the act of acquiring information from a system invariably disturbs the system The no-cloning theorem of quantum mechanics states that quantum information cannot be copied with perfect fidelity In classical computation, you can clearly copy a state with impunity The main fundamental difference of quantum mechanics is Bell’s theorem. This states that quantum mechanics is not a local hidden variable theory. All of the information in a quantum system is encoded in nonlocal correlations that have no classical analog Quantum Complexity Classically, the indivisible unit of information is the bit The quantum analog is the qubit, which is a vector in a 2D complex vector space with an inner product ie. $|\phi > = a|0> + b | 1>$ Performing measurements on $|\phi>$ gives a probabilistic output Generalizing to N qubits, you need $2^{N}$ basis vectors to completely specify the state This can be compactly written as $\Sigma_{x=0}^{2^{N}-1} a_{x} | x>$ where $a_{x}$ are complex numbers Thus, any quantum computation consists of applying unitary transformations onto this N qubit representation. You can then measure the state by projecting onto one of the basis vectors A quantum computation is probabilistic by nature, so multiple runs are not guaranteed to yield the same result All of these operations can be simulated on a classical computer (re: vector representations, matrix multiplications, inner products) The trouble arises if you want to simulate a large number of qubits. Even 100 qubits requires manipulating way more than $10^{30}$ complex numbers as a vector You can’t divide and conquer the complexity of a quantum system due to the nonlocal correlations imparting the vast majority of the information content Standard computers are Turing complete: If given an infinite amount of time and infinite memory, any computation can be done Problems can be classified as “hard” or “easy” depending on how much time and memory they consume Describe how “hard” a problem is should be universal: it should not depend on the hardware you’re running on The standard distinction is between polynomial time algorithms and exponential time algorithms Simulating a quantum computer on a classical computer is not a polynomial time algorithm In light of this physical reality, a classical Turing machine is not an appropriate model for quantum computers Quantum Parallelism (Deutsch Problem) Imagine a function f(x) which takes a long time to compute. We want to know if f(0)==f(1) (constant) or f(0)!=f(1) (balanced). A classical computer would need to calculate both values, then compare Define a 2 qubit system x,y. Define a unitary operator U such that $|x>y> \rightarrow |x>|y \oplus f(x)>$ ie. flip the bit if f(x) acting on x is 1. Otherwise, don’t do anything Imagine that you prepare $x = \frac{1}{\sqrt{2}} (|0> + |1>)$ and $y = \frac{1}{\sqrt{2}} (|0> -|1>)$ We then project the first qubit onto the basis $|\pm > = \frac{1}{\sqrt{2}}(|0>\pm |1>)$ We see that we get $|+>$ if the function is constant and $|->$ if the function is balanced. Hence, we only need to do 1 pass of the machine to get the answer! By feeding in a superposition of states, we can create a final state from which we can extract information (hallmark of quantum parallelism) Can generalize this readily to N bits: Set the input to a known superposition, calculate the function f(x), then choose some funny basis or do additional transformations to teach out global properties of f Errors Information is stored in nonlocal correlations of the system. A large quantum system can couple to its’ environment, which “spreads out” the correlation to the environment. You can’t measure the environment perfectly, so you inevitably lose information Can think of this as decoherence: the enviornment is continually “measuring” the state, which causes it to drift over time A separate problem is that we don’t have perfect quantum gates. Trying to implement perfect unitary transformation ain’t happening. There will always be some error of order $\epsilon$ from the ideal You need some sort of error correcting codes to deal with this reality Classically, there are a ton of error correcting codes (think Hamming codes) The simplest is just repetition: make N copies of your data, then do majority polling to determine the correct output Quantum mechanically, more things can go wrong You can have the standard bit errors like in classical computing You can have phase errors, where, for instance, $|0> \rightarrow - |0>$ This phase shift is continuous, unlike the jump discontinuities of classical bit flips You can’t clone a quantum system with perfect fidelity You can’t measure the system without disturbing it Unsurprisingly, Peter Shor developed the first quantum error-correcting code, which can be thought of as an extension of the N-bit repetition code. For simplicity, suppose that we want to encode one qubit as 3 qubits So the state $a|0> + b |1> = a|000> + b|111>$ We want to be able to detect errors without destroying this superposition If I measure the first qubit and get 0, then all the states collapse to 0 and I lose information about the coefficients a and b What if you instead measure pairs of qubits? Let the 3 qubit state be represented as $|x,y,z>$. Define $x \oplus y$ denote XOR-ing the bits Define the two-bit observable $(y \oplus z, x \oplus z)$ For our 3-qubit system, this observable will dictate which index an error occurred at 0 for no error, and then 1,2,3 (in binary) to denote that qubit 1,2,3 has been flipped (going left to right) What if there is a small deviation in state? Say that the following perturbations occur: $|000> \rightarrow |000> + \epsilon |100>$ and $|111> \rightarrow |111> + \epsilon |011>$ When you project onto the eigenstate of $(y \oplus z, x \oplus z)$, most of the time (with probability $1-|\epsilon|^{2}$), the state will get reprojected back to the original state. The $|\epsilon|^{2}$ outcome just sets the observable to (0,1), indicating a bit flip in position 1 This scheme fails with a probability of $|\epsilon|^{4}$ if multiple bit-flips happen The above scheme allows us to: Make a measurement of the system without damaging information ( you gain information about the error location, but not the exact configuration of your system) Small continuous errors either get corrected immediately, or produce a large discrete error which can be corrected Avoid the no cloning theorem ($a|000>+b|111>$ is not the same as $(a|0>+b|1>)^{3}$) The only source of error left are phase errors Do a similar trick used to address the other problems: Define a set of 9-qubit states $|0> = \frac{1}{2^{\frac{3}{2}}}(|000> + |111>)(|000> + |111>)(|000> + |111>)$ and $|1> = \frac{1}{2^{\frac{3}{2}}}(|000> - |111>)(|000> - |111>)(|000> - |111>)$ Define a “cluster” as state within each parentheses If a bit flip occurs within a cluster, you can correct it with the aforementioned scheme Observe that the phases between each cluster is aligned (ie. they are all + or all -). If a phase flip occurs, then the phases between each cluster won’t be aligned (ie. one is + and one is 1) You can generate a similar index observable, but now it’s a 6-bit observable. If this observable is non-zero, than you can detect which cluster has a different sign compared to the others and correct for it The most general single-qubit unitary transformation can be expanded to order $\epsilon$ in terms of the Pauli matrices: $U = 1+ i\epsilon_{x}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}+ i \epsilon_{y}\begin{pmatrix} 0 & -i \\ 0 & i \end{pmatrix}+ i\epsilon_{z} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$ Each term can be though of as a bit flip, a phase flip, and a combination of the two respectively The key takeaways of the quantum repetition code are: That the errors became digitized: Either you are in a state of no error, or you are in a discrete set of error states that you know how to recover from You can measure the error without measuring the data (re: sampling all of the qubits) The errors are local (re: uncorrelated), and the encoded information is nonlocal, so as long as you keep your measurements to single qubits, then you can’t get information out of the system States and Ensembles Quantum Mechanics Axioms A state is a ray in Hilbert space A Hilbert space is a vector space over the complex numbers You can define an inner product which maps an ordered pair of vectors to C that is strictly positive, linear, and skew symmetric (ie. swapping order of vectors is complex conjugation) An additional constraint is that there is some notion of a norm: $||\phi|| = <\phi| \phi>^{\frac{1}{2}}$ A ray is an equivalence class of vectors which differ by a multiplication by a nonzero complex scalar. We choose the representative of this class to have a unit norm An observable is a property of a physical system which can be measured Must be a linear, self-adjoint operator(ie. $<\psi| A \psi> = < A^{\dag} \psi | \psi >$ ) As a result of the self-adjoint property, we can write $ A= \Sigma_{n} a_{n} P_{n}$ where the projection operator $P_{n}$ satisfies $P_{n}P_{m} = \delta_{n,m} P_{n}$ and $P^{\dag} = P$ We can make a measurement of a state by using the projection operator The dynamics of the system is unitary (re: probability preserving) and is governed by $\frac{d}{dt}|\phi> = -i H |\phi>$ We can describe the time evolution of a state as $|\phi(t)> = U(t) | \phi(0)>$ If H is time independent, then $U = exp(-i t H)$ Schrodinger equation is linear in the Hamiltonian If we want to compose systems A and B, we take the tensor product of their Hilbert spaces. So states $|\phi_{A}>$ and $|\phi_{B}>$ get combined to the state $|\phi_{A}> \otimes |\phi_{B}>$ You can compose a basis of the combined system as a tensor product of the bases of each original state: $|i_{A}> \otimes \mu_{B}>$ You can define the operator $X = A \otimes I$ to only act on the first subsystem Qubits A qubit is a state in a 2D Hilbert space that can take the form $a|0> + b|1>$ Symmetries Any symmetry of a quantum system must leave the probabilities untouched This implies that a symmetry is a automorphism of the Hilbert space which preserves the absolute values of inner products for all members of the space Each symmetry maps onto either a unitary or anti-unitary operator Anti-unitarity doesn’t matter for continuous symmetries Compositions of symmetries should also be a symmetry (up to an overall phase factor). This follows from group theory Symmetries should commute with the dynamics of a system Let R be the symmetry in question. The following must hold: $U(R) exp(-itH) = exp(-itH) U(R)$ Namely, that symmetries should commute with the time evolution operator For a continuous symmetry, we can let R get arbitrarily close to the identity: $R = I + \epsilon T$. This implies that, to first order, $U = 1-i\epsilon Q$ where Q is unitary. This in turn implies that Q commutes with H We call Q the generator Since any finite transformation can be written as a product of infinitesimal ones: $R = (1+\frac{\theta}{N} T)^{N} \rightarrow U = exp(i\theta Q)$, knowing how the infinitesimal symmetry transformations are represented allows you do finite transformations Rotations A finite rotation is given by $R(\hat{n}, \theta) = exp(-i\theta\hat{n}\cdot J)$ $\hat{n}$ is the axis of rotation, $\theta$ is the rotation angle, and $\vec{J}$ is the angular momentum The associated commutation relationship for angular momentum is $[J_{k}, J_{l}] = i \epsilon_{klm} J_{m}$ The simplest non-trivial irrreducible representation of angular momentum is 2D, as given by the Pauli matrices: $J_{k} = \frac{1}{2}\sigma_{k}$ The Pauli matrices satisfy the anti-commutation relationship: $\sigma_{k}\sigma_{l}+\sigma_{l}\sigma_{k} = 2\delta_{lk} I$ We can use the Pauli matrices to write any finite rotation as $U(\hat{n},\theta) = exp(-i\frac{\theta}{2} \hat{n}\cdot \vec{\sigma})$ There is a $2\pi$ ambiguity, which gives rise to spinor representations For rotation, you have the action $U(\hat{n} ,\theta = 2\pi) = -1$ The components of angular momentum transform under rotations like a vector $U(R)J_{k} U(R)^{\dag} = R_{kl} J_{l}$ This implies that if state $|m>$ is an eigenstate of $J_{3}$, then $U(R)|m>$ is an eigenstate of $RJ_{3}$ with the same eigenvalue The above implies that we can construct eigenstates of angular momentum along the axis $\hat{n} = < \sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta >$ by applying a rotation through $\theta$ to the z axis More explicitly: all direction measurements can be performed by first rotating the $\hat{n}$ axis to the z axis, and then measuring along z The axis $\hat{n’} = < -\sin \phi, \cos \phi, 0>$ is the axis around which you perform the counter clockwise rotation by $\theta$ Density Operator In real quantum systems, we don’t have a closed system: there is always some interaction with the environment This can cause some axioms of quantum mechanics to appear to be violated States are not rays Measurements are not orthogonal projections Evolution is not unitary Say that you have a two qubit system A and B, but you can only observe A. How can you characterize the observations made on A alone? We can define the density matrix of the system as $\rho = \Sigma_{i} p_{i} |i><i|$ This of this as a representation of the ensemble of possible quantum states, each with their own probability We can then define expectation values of observable Q as $tr(Q \rho)$ This idea extends to any bipartite system: You can calculate the expectation value of one system by partial tracing over the subsystem In general, we have that $\rho_{A} = tr_{B}(|\phi><\phi|)$ A general density matrix (in diagonal form) can also be written as $\rho_{A} = \Sigma_{a} p_{a} |a>< a|$ The density matrix has the following properties: self-adjoint ($\rho_{A} = \rho_{A}^{\dag}$) $\rho_{A}$ is positive definite $tr(\rho_{A}) = 1$ For so called “pure states”, we have that $\rho^{2} = \rho$ You can think of a pure state as when the subsystem is also a ray (re: no mixing of components) Bloch Sphere We can write any 2x2 self-adjoint matrix in the basis of the Pauli matrices and the identity $\rho(\vec{P}) = \frac{1}{2}(I+\vec{P}\cdot \vec{\sigma})$ where $\vec{P}$ is some 3D vector The $\frac{1}{2}$ arises because $tr(\rho)=1$ and the Pauli matrices are traceless, hence we need to scale down the identity so that $\rho$ is a density operator The eigenvalues for $\rho$ must be non-negative (since they are interpreted as probabilities) Looking at the determinant of $\rho$ we get $det(\rho) = \frac{1}{4}(1-\vec{P}^{2})$ The above constrains $\vec{P}^{2} \leq 1$, which corresponds the the unit ball So Bloch “Sphere” is a bit of a misnomer The actual sphere (ie. $|\vec{P}| = 1$) corresponds to density matrices with a vanishing determinant. Since the trace is 1, the eigenvalues are either 0 or 1. Hence, the boundary of the ball are pure states Accordingly, we can write a pure state as: $\rho(\hat{n}) = \frac{1}{2}(I+\hat{n}\cdot \vec{\sigma})$ We can write this in spherical coordinates as well: $\rho(\theta,\phi) = \frac{1}{2} I + \frac{1}{2} \begin{pmatrix} \cos \theta & \sin \theta exp(-i\phi) \\ \sin \theta exp(i\phi) & -\cos(\theta) \end{pmatrix}$ This is easily derived from the vector $\Phi(\theta,\phi)> = \begin{pmatrix} exp(-i\phi/2) \cos(\theta/2) \\ exp(i\phi/2) \sin(\theta/2)\end{pmatrix}$ We remember that $\hat{n} = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$ You can construct the density matrix of a system by measuring $\hat{p}\cdot\vec{\sigma}$ across the 3 linearly independent axes Schmidt Decomposition The standard orthonormal basis for a bipartite system is $|\phi_{AB}> = \Sigma_{\alpha, \mu} \phi_{\alpha\mu} |a_{A}> \otimes |\mu_{B}>$ An alternative representation is the so called Schmidt decomposition: $|\phi_{AB}> = \Sigma_{i} \sqrt{p_{i}} |i_{A}> \otimes | i’_{B} >$ This is just the SVD of the original operator $\tilde{i_{B}}> = \Sigma_{\mu} \phi_{i\mu} |\mu_{B}>$ These turn out to be orthogonal to each other the transformation from unprimed to primed can be encoded in a unitary transformation $U_{B}$ The $p_{i}$ are the singular values This decomposition holds for any matrix Can derive as follows: Assume some general composite state: $|\phi> = \Sigma_{i,\mu} |i_{A}>\otimes \mu_{B}> = \Sigma_{i} |i> \tilde{i}>$ Choose a basis for which $\rho_{A}$ is diagonal: $\rho_{A} = \Sigma_{i} p_{i} |i><i|$ We know that $\rho_{A} = tr_{B}(|\phi><\phi|)$ Comparing the two results shows you that $<\tilde{i}|\tilde{j}> = p_{i} \delta_{ij}$ So $\tilde{i}>$ are actually orthogonal. We normalize them by rescaling each basis by $\sqrt{p_{i}}$ Calculating the density matrices $\rho_{A}$ and $\rho_{B}$ by tracing on the other subsystem. Applying the Schmidt decomposition to each, we find that $\rho_{A}$ and $\rho_{B}$ share the same non-zero singular values If A and B are different dimensions, the difference is made up with zeros Entanglement Define the Schmidt number as the number of non-zero singular values in a bipartite pure state If the schmidt number is greater than one, then the state is entangled. Otherwise, it’s seperable So if the state can be written as a direct product of the subspaces (re: $|\phi_{AB}> = |\phi_{A}>\otimes |\phi_{B}>$) then it’s seperable A seperable state is not necessarily uncorrelated: If you have $\uparrow_{A}>\uparrow_{B}>$, then the system is seperable but correlated An entangled state is fundamentally different in that there are non-local quantum correlations (re: there must be interactions between the subsystems) Ensemble Interpretation Ambiguity The density operator is self-adjoint, nonnegative and has $tr(\rho) = 1$ From the above, you can construct a density matrix as a convex linear combination between two other density matrices which still satisfies all the properties: $\rho(\lambda) = \lambda \rho_{1} + (1-\lambda) \rho_{2}$ where $0 \leq \lambda \leq 1$ The above implies that the density operators are a convex subset of the real vector space of dxd Hermetian operators Pure states can’t be defined as a complex sum of other states Define the pure state $\rho = |\phi><\phi|$ and define $\phi_{\perp}>$ as some orthogonal vector to $\rho$ Suppose that $\rho$ can be written as a convex linear combination $<\phi_{\perp}|\rho|\phi_{\perp}> = 0 = \lambda <\phi_{\perp}|\rho_{1}|\phi_{\perp} + (1-\lambda) <\phi_{\perp}|\rho_{2}|\phi_{\perp}>$ For the above to hold, both terms must vanish. If $\lambda$ is 0 or 1, then $\rho_{1}= \rho_{2} = \rho$ Otherwise, $\rho_{1}$ and $\rho_{2}$ are orthogonal to $|\phi_{\perp}>$. Since $\phi_{\perp}>$ can be any orthogonal vector, you arrive at the same conclusion These pure states are extremal points of the set You can see this structure on the Bloch sphere (extremal points are on the boundary) d=2 is a special case where the extremal points are all pure states. This does not hold for d>2 The convexity of the set of density matrices has a physical interpretation: Calculating the expectation value of some observable M, we have $<M> = tr(\rho(\lambda) M)$ Suppose that we have two states $\rho_{1}$ and $\rho_{2}$, where the former has a probability of $\lambda$ and the latter has probability of $1-\lambda$ Taking the expectation value of M on this state also gives $tr(M \rho(\lambda))$ So this preperation of $\rho$ gives the same observable, regardless of $\rho_{1}$ and $\rho_{2}$ used in preparation Since pure states can’t be a convex linear combination of other density matrices, there is an unambiguous way of preparing a pure state This ambiguity in mixed states stands in stark constrast to classical systems, where there is a unique preparation method to generate a probability distribution Faster Than Light Communication?! (No…) Suppose that qubit A has the density matrix $\rho_{A} = \frac{1}{2}(|\uparrow_{z}><\uparrow_{z}|+|\downarrow_{z}><\downarrow_{z}|$ This density matrix could arise from an entangled bipartite pure state $|\phi_{AB}>$ with the Schmidt decomposition $\rho_{A} = \frac{1}{2}(|\uparrow_{zA}><\uparrow_{zB}|+|\downarrow_{zA}><\downarrow_{zB}|)$ We can realize the ensemble interpretation of A via measuring qubit B Since $\rho_{A}$ has degenerate eigenvalues, this Schmidt basis is not unique, so any direction works (just apply a unitary tranformation V to A and $V^{*}$ to B to rotate your basis) Possible faster than light information propagation: Prepare many copies of this entangled state. Alice takes all the A qubits to one location, while Bob takes all the other B qubits to another Bob measures along the z axis for all of the prepared states. Alice, concurrently, measures all of her qubits to see which axis was read. For now, ignore the fact that the measurements need to be simultaneous (re: Bob calls Alice and tells her to measure something) The main problem with this is that density matrix $\rho_{A}$ is the same between the two setups, so it’s fundamentally impossible to distinguish between the two states Follows from the convexity of the density matrices set Quantum Erasers The density matrix $\rho_{A} = \frac{1}{2} I$ is an incoherent mixture between $|\uparrow>$ and $\downarrow>$ A coherent superposition would be $|\phi> = \frac{1}{\sqrt{2}}( |\uparrow > \pm \downarrow>)$ The distinction is that the relative phase of a coherent superposition is observable Entanglement causes decoherence Imagine an entangled state like that in previous section Bob makes a measurement along the x axis and sends his measurement result to Alice This forces Alice’s spin to be a pure state along the x axis, which in turn can be interpreted as a coherent superposition of z axis spins So even though Alice initially had an incoherent state, Bob’s measurement caused Alice’s state to become coherent Another thought experiment: Bob uses a Stern-Gerlach experiment to measure his z axis spin. This precludes Alice from having a coherent superposition along the z axis Bob refocuses the two beams of the Stern-Gerlach which then passes through a Stern-Gerlach along the x axis. Alice’s coherence along the z axis is now restored! (re: Alice is in a known x axis configuration, but the z axis position is lost) This situation is called a quantum eraser. Coherence of a state can be restored if an orthogonal measuremnt occurs Information is physical: measuring an system fundamentally changes the physical description of the system HJW Theorem How do you extend the quantum eraser to multiple qubits with a more general density matrix? Consider the general density matrix: $\rho_{A} = \Sigma_{i} p_{i} |\phi_{i}><\phi_{i}|$ where $\Sigma p_{i} = 1$ Don’t assume that $|\phi_{i}>$ are orthogonal, but do assume they are normalized Imagine some bipartite system such that performing a partial trace over the subsystem B yields $\rho_{A}$ $|\phi_{1AB}> = \Sigma_{i} \sqrt{p_{i}} |\phi_{iA}>\otimes |\alpha_{iB}>$ $<\alpha_{i}|\alpha_{j}> = \delta_{ij}$ $tr_{B}(|\phi_{1AB}><\phi_{1AB}|) = \rho_{A}$ The construction of this bipartite state is called purification. It can be though of as representing a mixed state as a pure state in a higher dimensional Hilbert space Performing a measurement on B is projecting onto a $|\alpha_{iB}>$ basis, which in turn forces system A to be in the pure state $|\phi_{i}><\phi_{i}|$ Is there a different ensemble interpretation of $\rho_{A}$ that we can construct by making a different measurement of B? Let $\rho_{A} = \Sigma_{\mu} q_{\mu} |\phi_{\mu}>< \phi_{\mu}|$ so a different ensemble of pure states There is a similar purification of A: $|\phi_{2AB}> = \Sigma_{\mu} \sqrt{q_{\mu}} |\phi_{\mu A}>\otimes |\beta_{\mu B}>$ How are $|\phi_{1}>$ and $\phi_{2}>$ related? Partial tracing over B yields the same density matrix in both cases. They both have Schmidt decompositioons: $\phi_{1}> = \Sigma_{k} \sqrt{\lambda_{k}} |k_{A}> \otimes | k_{1B}>$ $\phi_{2}> = \Sigma_{k} \sqrt{\lambda_{k}} |k_{A}> \otimes | k_{2B}>$ Alternatively: $|\phi_{1AB}> = (I_{A} \otimes U_{B}) | \phi_{2AB}>$ Since $k_{1}>$ and $k_{2}>$ are both orthonormal bases for B, there is a unitary transformation $U_{B}$ between the two Hence, $|\phi_{1}>$ and $|\phi_{2}>$ are the same purification. You just need to change which direction in B you measure along Suppose that we have many ensembles that realize $\rho_{A}$, where the max number of ensembles is d We can then choose a Hilbert space $H_{B}$ of dimension d and a pure state $|\phi_{AB}> \in H_{A} \otimes H_{B}$ such that any of the ensembles can be realized by measuring a suitable observable of B This is the HJW theorem The general density matrix mixes the pure states incoherently (re: can’t detect the relative phases of these states) So you can erase information by making a measurement in B, and restore the coherence by making a different measurement Fidelity Suppose that you have two density operators $\rho$ and $\sigma$. The fidelity is defined as $(tr(\sqrt{\rho^{\frac{1}{2}} \sigma \rho^{\frac{1}{2}}}))^{2}$ Think of this as the distance metric between two mixed states This is well defined since $\rho$ and $\sigma$ are positive definite matrices; hence, you can take the square root via the spectral theorem The fidelity is bound between 0 and 1. 1 occurs if $\rho$ and $\sigma$ are identical If $\rho = |\phi>< \phi|$ (ie. a pure state), the fidelity becomes $F(\rho,\sigma) = <\phi|\sigma|\phi>$ If both density matrices are pure states, the fidelity reduces to the Born interpretation An alterrnative definition of the fidelity is $F(\rho,\sigma) = ||\sigma^{\frac{1}{2}} \rho^{\frac{1}{2}}||_{1}$ $||A||_{1} = tr \sqrt{A^{\dag}A}$ For Hermetian matrices, take the sum of the absolute values of the eigenvalues The symmetry between the arguments of the fidelity is more manifest in this notation, since for any Hermetian matrices A and B, we have that $||AB|| $ $= ||BA||$ This follows from the fact that BAAB and ABBA have the same eigenvalues (hence the same traces) Uhlmann’s Theorem How does the fidelity of two density operators relate to the overlap of their purifications? Define $|\Phi_{AB}>$ as the purifcation of the density operator $\rho_{A}$ where $\rho_{A} = tr_{B}(|\Phi><\Phi|)$ Suppose that $\rho = \Sigma_{i} p_{i} |i><i|$ where $|i_{A}>$ is a orthonormal basis for system A The associated purifcation is then $\Phi_{\rho}> = \Sigma_{i} \sqrt{p_{i}} |i_{A}> \otimes |i_{B}>$ By the HJW theorem, the general purification is $\Phi_{\rho}(V)> = I \otimes V |\Phi_{\rho}>$ where V is unitary This can also be written as $(\rho^{\frac{1}{2}} \otimes V) |\tilde{\Phi}>$ where $\tilde{\Phi}> = \Sigma_{i} |i_{A}> \otimes | i_{B}>$, which is a normalized maximally entangled state Suppose that we have two density operators $\rho$ and $\sigma$ acting on A. What is the inner product of their purifications? $<\Phi_{\sigma}(W)|\Phi_{\rho}(V)> = <\tilde{\Phi}| \sigma^{\frac{1}{2}} \rho^{\frac{1}{2}} \otimes W^{\dag} V |\tilde{\Phi}>$ Using the fact that $U\otimes I |\tilde{\Phi}> = I \otimes U^{T} |\tilde{\Phi}>$, we can write the inner product of the purifications as $<\tilde{\Phi}|\sigma^{\frac{1}{2}}\rho^{\frac{1}{2}}U \otimes I | \tilde{\Phi}> = tr( \sigma^{\frac{1}{2}}\rho^{\frac{1}{2}}U )$ where $U = (W^{\dag}V)^{T}$ We use the polar decomposition: $A = U’ \sqrt{A^{\dag}A}$ This is the “obvious” fact that any square complex matrix can be factorized as $A=UP$, where U is a unitary matrix and P is a positive semi-definite Hermetian matrix so $ tr( \sigma^{\frac{1}{2}}\rho^{\frac{1}{2}}U ) = tr(UU’\sqrt{\rho^{\frac{1}{2}}\sigma \rho^{\frac{1}{2}}})$ This is very close to the square root of the fidelity. If we select $U’ = U^{-1}$, then they are the same This selection can be thought of as maximizing the inner product of the purifications (this can be seen by Schmit decomposing $ \sqrt{\rho^{\frac{1}{2}}\sigma \rho^{\frac{1}{2}}}$ into its’ non-negative eigenvalues $\lambda_{a}$ and eigenvectors $|a>$, and realizing that $<a|UU’|a>$ is maximized when $UU’ = I$ The final result of Uhlmann’s theorem is that $F(\rho,\sigma) = (tr( \sqrt{\rho^{\frac{1}{2}}\sigma \rho^{\frac{1}{2}}}))^{2} = max_{V,W} | <\phi_{\sigma}(W)|\phi_{\rho}(V)>|^{2}$ An immediate corrollary to this is the monotonicity of fidelity: $F(\rho_{AB},\rho_{AB}) \leq F(\rho_{A},\rho_{A})$. In words: tracing out a subsystem cannot decrease the fidelity of two density operators (follows from Uhlmann’s because purification of $\rho_{AB}$ also purfies $\rho_{A}$ Measurement and Evolution Suppose that we have some system S embedded in a larger system. We want to make a measurement on S. This process can be described as a orthonormal projection operator on the whole system. However, an orthonormal projection on the whole system might not be an orthonormal one on S Let’s generalize the notion of measurement: to measure an observable M, you need to modify the Hamiltonian of the world by coupling the observable M to to some other variable which represents our apparatus. This auxilliary system used to facililate the measurement is called “the pointer”, “the meter”, or “the ancilla” We can then prepare an eigenstate of M by observing the pointer More concretely, think of the pointer as a particle of mass m that propagates freely apart from the coupling We want to measure the position of the pointer, so we need the wavefunction to be spatially confined, but not too confined since otherwise the packet widens too quickly via uncertainty The width of the wavepacket will evolve like $\Delta x (t) \approx \Delta x + \frac{\hbar t}{m\Delta x}$ which minimizes when $\Delta x \approx \sqrt{\frac{\hbar t}{m}}$. This is a upper bound on the resolution of the pointer’s position $\Delta x = \sqrt{\frac{\hbar t}{m}}$ is called the standard quantum limit (SQL) For simplicity, assume the pointer mass is heavy enough for the SQL to not matter The Hamiltonian takes the form $H = \mathbf{H_{0}} + \frac{\mathbf{P^{2}}}{2m} + \lambda(t) \mathbf{M} \otimes \mathbf{P}$ where $\mathbf{H_{0}}$ is the unperturbed Hamiltonian, $\frac{\mathbf{P^{2}}}{2m}$ is the kinetic energy, $\lambda$ is the coupling Assume the kinetic energy term is negliglible due to the large mass of the pointer For simplicity, assume that $[H_{0}, M] = 0$, or that the measurement is quickly done so as to not perturb the evolution. The Hamiltonian then becomes $H \approx \lambda(t) M \otimes P$ Assume that $\lambda$ suddenly turns on a t=0 and turns off at t=T. The time evolution operator is then $U(T) \approx exp(i\lambda T \mathbf{M}\otimes \mathbf{P}) = \Sigma_{a} |a> exp(i\lambda T M_{a} \mathbf{P}) <a|$ where we diagonalize $\mathbf{M}$ Recall that $exp(-i x_{0} \mathbf{P}) \phi(x) = \phi(x-x_{0})$ (Follows from $P = -i\frac{d}{dx}$ and Taylor expansion) Suppose that our initial state is in a superposition of M eigenstates, initally unentangled with the position-space wavepacket $\phi(x)>$ (ie. $\Sigma_{a} \alpha_{a} | a> \otimes \phi(x)>$. Acting U(T) on this state gives $\Sigma_{a} \alpha_{a} |a> \otimes |\phi(x-\lambda T M_{a}) > $ The time evolution correlates the position of the pointer with the value of observable M Hence, observing the shift from x, we can figure out the associated eigenstate of M with probability $|\alpha_{a}|^{2}$ A classical example is the Stern-Gerlach experiment: you have some observable $\sigma_{3}$. You pass through a B field of the form $B_{3} = \lambda z$. The perturbed Hamiltonian is then $-\lambda \mu z \sigma_{3}$. This coupling imparts an impulse on the pointer (z is a translator of $P_{z}$) which can be observed as spin up or spin down Imagine that we can couple any observable to a pointer as described above. Observing the pointer will be some orthogonal projection in the Hibert space Let $[E_{a}]$ be some set of orthongonal projectors such that: $E_{a} = E_{a}^{\dag}$ $E_{a}E_{b} = \delta_{ab} E_{a}$ $\Sigma_{a=0}^{N-1} E_{a} = 1$ Define the unitary transformation $U = \Sigma_{a,b} E_{a} \otimes |b+a> <b|$ Addition of b + a is modulo N Act this unitary transformation on the initial state of the system: $|\Phi> = |\phi> \otimes |0>$ to get $U|\Psi> = \Sigma_{a} E_{a} |\phi> \otimes |a>$ The effect of this transformation on the density matrix due to measurement is to transform it like $\rho \rightarrow \Sigma_{a} E_{a} \rho E_{a}$ The probability of observing a then becomes $<\phi| E_{a}|\phi>$ So we can see that by coupling the system to a macroscopic apparatus, and doing appropriate unitary transformations, we can observe the pointer in its’ fiducial basis to perform any orthogonal measurement on the system

Referencing “Modern Compiler Implementation in ML” by Andrew W. Appel. Lexical Analysis This refers to breaking up the input source code into individual words/tokens A “token” is a sequence of characters that can be treated as a unit in the grammer of a programming language Some common tokens which are typical in programming languages: ID: think names of any kind (variable, function classes etc.) NUM: Think integers REAL: Think floating point RESERVED: Words which are reserved in the programming language. “If” is a typical one These reserved words cannot be used as identifiers (ID) macros and preprocessor directives are done prior to lexing Only “weak” languages require a macro preprocessor (what does that mean?) Some of these token hae semantic values attached to them (think identifiers and literals), which give additional information Regular Expressions You could describe the lexical rules of a programming language in a normal language That’s very verbose and not precise A language can be though of as a set of strings A string is a finite sequence of symbols the symbols are drawn from a finite alphabet A language can be a finite or infinite set The Pascal language is an infinite set, since it includes all possible Pascal programs The language of C reserved words is the set of all alphabetic strings that cannot be used as identifiiers in the C programming language In this viewpoint, we just need to figure out if a given string is within the language set or not Regular expressions provide a compact notation to represent a set of strings Symbol: each symbol the alphabet of the language, the regular expression $\mathbf{a}$ denotes the language containing just the string “a” Alternation: given two regular expressions M and N, the alternation operator ("|") creates a new regular expression $M | N$. A string is in this new language if it was in language M or language N Concatenation: Given regular expressions M and N, $M\dot N$ creates a new regular expression. A string is in $M \dot N$ if said string the concatenation of any two strings in the sub languages (so $(a|b) \dot a$ is the language of strings “aa” and “ba”) $\epsilon$ represents the languages whose only elements is the empty string Repetition: Given a regular expression M, the Kleene closure M* defines a new language. A string is in M* if said string is a concatenation of zero or more strings, all of which are in M (so $((a|b)\dot a)*$ is an infinite set of strings $("", “aa”, “ba”,“aaaa”,“baaa”, “aaba”, …)$) The above are the base operations. You can make some notational short hand: $[abcd] = (a| b|c|d )$ $[b-g] = [bcdefg]$ $M? = (M| \epsilon)$ $M+ = (M\cdot M*)$ There are also some convenience things: a “.” represents any single character except newline A string wrapped in quotes (re: “a.+*” stands for itself literally) Using the above, we can write some common tokens as: Regex Tokens if $250 [a-z][a-z0-9]* (ID) [0-9]+ (NUM) ([0-9]+"."[0-9]*) ([0-9]*"."[0-9]+) A given lexical specification is complete if it always matches some initial substring of the input Can achieve this by having a rule that matches any single character There is some ambiguity in these rules, like what if there are multiple matches? Lex uses two disambiguation rules Accept the longest match which satisfies the regex rules If multiple regular expressions match a longest substring, then the first regular expression determines it’s type This implies that the order you list the regular expressions in matters As an example, “if8” will get matched to an identifier, while “if” gets matched to the IF token Deterministic Finite Automata (DFA) Need to be able to implement regex in a computer program. One representation which lends itself well to this is a finite automaton Consists of a finite set of states, along which edges (labelled by a symbol) which transitions between each state One state is denoted as the start state, while there can be multiple states which are final states (re: a regex expression needs to end on one of these) In a deterministic finite automaton (DFA), no two edges leaving from the same state are labelled with the same symbol a DFA can take in a string, and either accept or reject said string. From the start state, each character in the string induces a transition to some other state. If there is no valid transition available, or if after consuming the string, the DFA is not in a designated final state, then the string is rejected You can convert this graphical representation of the DFA into a transition matrix This is a 2D array which is subscripted by the state number (re: all possible states) and the input character (re: all of the alphabet of the language) Each entry in the matrix tells you which state you need to jump to There is a “dead” state that just loops to itself on all characters (used to encode the absence of an edge) You also need a “finality” array, which maps state number to tokens If you want the longest match for a given input string, you also need to keep track of the last final state the DFA was in, as well as the associated character position Nondeterministic Finite Automata (NFA) Unlike DFAs, NFAs can have edges which share the same symbol. In addition, special $\epsilon$ edges can be traversed without consuming any input The non-determinism comes from the fact that you need to pick which branch you go down It’s easy to convert a regex into a NFA. Just follow the above rules. At the end, you should end up with some NFA with a tail and a head Each symbol gets converted into a tail (re: transition) and a head (re: a state) Same for an $\epsilon$ edge Concatenations are just chains Alternations are just split paths with an epsilon edge buffer Repetition is just a feedback loop This is not necessarily the most efficient realization of said regular expression, but it works Each regular expression is either a primitive (a single symbol or $\epsilon$) or it is made from smaller expressions. Ditto for NFAs Going from a NFA to DFA is harder In order to get the longest match, you need to be able to traverse all of the parallel $\epsilon$ edges at a given node (This applies recursively) This describes the notion of an $\epsilon$-closure Let edge(s,c) be the set of all NFA states reachable by following a single edge with label c from state s For a set of states S, closure(S) is the set of states that can be reached from a state in S without consuming any of the input (re: you cross all the reachable $\epsilon$ edges) This can be represented mathematically. closure(S) is the smallest set T such that $T = S \cup (\cup_{s\in T} edge(s, \epsilon)) $ T can be calculated in a straight forward manner: while $T_{new} \neq T$, let $T_{new} = T \cup (\cup_{s\in T} edge(s, \epsilon))$ This converges since T can only grow every iteration. There is a finite number of states, so the algorithm must terminate. If the cardinality doesn’t grow ($T_{new}= T$), this means that T includes $(\cup_{s\in T} edge(s, \epsilon))$ Suppose that you have a set d of NFA states you need to visit. Define a function DFAedge(d,c) where closure( $\cup_{s\in d}$ edge(s,c)) Can think of this as traversing the edges associated with a character, and then traversing all of the epsilons You can simulate an NFA via the following algorithm let d be the closure of your start state $s_{1}$ for each character in you input string, calculate DFAedge(d,$c_{i}$), then update you set d The above outlines how to construct a DFA using closures and DFAedge $\Sigma$ denotes the alphabet of the language In words, you start with the set generated by the closure of the start state. Iterate through this set. If a character matches one of the edges, check if the connecting state is already in the set. If it is, set that the associated element in the transition matrix to the next state. Otherwise, append the state to the queue, and set the transition matrix element. Do this until you have gone through all of the states This algorithm doesn’t reach unreachable states, which prevents exponential blowups This is not an optimal automaton

Following Baumann’s Cosmology. Conventions We are assuming a (- + + + ) signature for the metric. Greek letters for spacetime indices, and Latin letters for spatial indices overdots for physical time derivatives and primes for conformal time derivatives Qualitative Overview Meters are a cumbersome unit to work with at universe scales. We use lightyears and parsecs instead 1 ly is roughly 9.5E15 m 1 parsec (pc) is defined as the distance at which the radius of the Earth’s orbit around the Sun subtends an angle of one arcsecond (re: $\frac{1}{3600}$ of a degree) In light years, 1 parsec is around 3.26 lightyears For a sense of scale: The nearest star (Proxima Centauri) is 4.2 lightyears away Our Galaxy, the Milky Way, is around 30 kpc across The nearest galaxy (Andromeda) is around 2 million light years away (2 Mpc) The local group (~ 50 nearby galaxies) is around 10 Mpc in diameter The local supercluster (Laniakea) is 500 Mpc The observable universe is around 14 Gpc or 46.5 billion lightyears This is ostensibly larger than the age of the universe (13.8 lightyears) due to expansion The consituents of the universe are: Ordinary matter: This is us. It’s 5% Dark matter: This is the vast majority of the matter contents Major experimental evidence for this is increased rotational speeds of hydrogen gas in the outer reaches of galaxies (Vera Rubin) Only can be explained if galaxies embedded in halos of dark matter Another major piece is lensing of the CMB is in tension with measurements of element abundances without non-baryonic dark matter Also, the anisotropies of the CMB wouldn’t manifest as they do without dark matter Dark Energy: This is the biggest component, and is roughly 70% of the universe’s energy content Some negative pressure associated with the universe’s expansion Big Bang In the beginning, it was hot and dense, with all particle species in roughly equal abundances At some point, it expanded rapidly, and cooled off in the process We can map out the various phases by converting these temperatures to energies: EW phase transition (100 GeV): the electroweak symmetry is broken, allowing EM and the weak force to become distinct entities and particles gain mass via the Higgs As the temperature drops below the mass of particle species, particle antiparticle annihilation begins while pair production becomes inefficient QCD phase transition (150 MeV): the remaining quarks condense to hadrons At the grand old age of 1, the neutrinos decouple from equillibrium due to the expansion of the universe creating the cosmic neutrino background ($c\nu B$) Big Bang Nucleosynthesis (BBN) occurs around 1 minute in, creating elements lighter than lithium-7 370,000 years later, the universe cooled enough to allow recombination to occur. This allowed the photons to stream through the universe (CMB) There are two periods which we don’t know exactly when they occured dark matter production: lots of theories here (WIMPs, axions etc.) baryogenesis: The symmetry breaking between matter and antimatter (1 part in $10^{10}$) to give the matter to photon ratio observed Structure formation The first stars (Population III stars) are boorn 100 millions years after the Big Bang They were massive, and died off rapidly. The UV light emitted heated up the surrounding gas, leading to reionization They could have created the supermassive black holes They created the heavier elements in the universe The first galaxies formed around 1 billion years after the Big Bang The distribution of galaxies is not uniform (on small scales): there are spatial correlations Inflation We observe some objects in the universe to be older than the universe’s age Inflation helps solve this: if the distances between objects suddenly increases in a short time frame, you can explain these apparent causality violations This period would need to happen before the Big Bang In 1 billionth of a trillionth of a trillionth of a second, the universe doubled in size about 80 times The Expanding Universe GR Concept Refresher Spacetime is defined by some metric: $ds^{2} = \Sigma_{\mu,\nu}^{3} g_{\mu\nu}dx^{\mu}dx^{\nu}$ $ds^{2}$ is invariant to all observers In Minkowski space, this reduces to $ds^{2} = -c^{2}dt^{2} + \delta_{ij} dx^{i} dx^{j}$ the metric is used to convert between contravariant vectors ($A_{\mu}$) and covariant vectors ($A_{\mu}$) You can contract a convariant with a contravariant vector via the metric to produce an invariant scalar FRW Metric On large enough scales, the distribution of mass is isotropic (same in all directions) and homogeneous (same at every point in space) The metric then takes the form $ds^{2} = -c^{2} dt^{2} + a^{2}(t) d\vec{l}^{2}$ homogeneous and isotropic spacetimes must have constant intrinsic curvature This means that we can write $d\vec{l}^{2} = dx^{2}+ k \frac{(x \cdot dx)^{2}}{R_{0}^{2}-kx^{2}}$ $R_{0}$ denotes the intrinsic curvature of the space In spherical coordinates, we have $dl^{2} = \frac{dr^{2}}{1-k\frac{r^{2}}{R_{0}^{2}}}+ r^{2} d\Omega^{2}$ k can be 0 for flat space, 1 for spherical, and -1 for hyperbolic The RW metric becomes $ds^{2} = -c^{2} dt^{2} + a^{2}(t)(\frac{dr^{2}}{1-k\frac{r^{2}}{R_{0}^{2}}}+ r^{2} d\Omega^{2})$ The line element has a rescaling symmetry $a \rightarrow \lambda a$, $ r\rightarrow r/ \lambda$ and $R_{0} \rightarrow R_{0}/\lambda$ This allows $a(t_{0}) = 1$ (re: the scale factor at the present time) as well as $R_{0} = R(t_{0})$ to be interpreted as the physical curvature today we call r a comoving coordinate, while $r_{phys} = a(t) r$ is the physical coordinate The physical velocity is just the time derivative of $r_{phys}$ $v_{phys} = H r_{phys} + v_{pec}$ We define $H = \frac{\dot{a}}{a}$ to be the Hubble parameter $v_{pec}$ to be the velocity of a comoving observer (re: someone who follows the Hubble flow) We can make the change of coordinates $d\chi = \frac{dr}{\sqrt{1-\frac{kr^{2}}{R_{0}^{2}}}}$ The FRW metric becomes: $ds^{2} = -c^{2} dt^{2} + a^{2}(t) (d\chi^{2}+$ Kinematics Define a free particle lagrangian in an arbitrary coordinate system as $\mathbf{L} = \frac{m}{2} g_{ij}(x^{k}) \dot{x}^{i}\dot{x}^{j}$ Substitute into the Euler-Lagrange equution to get $\frac{d^{2} x^{i}}{dt^{2}} = - \Gamma_{ab}^{i} \frac{dx^{a}}{dt}\frac{dx^{b}}{dt}$ We define $\Gamma_{ab}^{i} = \frac{1}{2} g^{ij} (\partial_{a} g_{jb} + \partial_{b} g_{ja} + \partial_{j} g_{ab})$ These are the Christoffel symbols This is easily extended to 4 dimensions. Simply replace the latin indices for greek ones (re: add time equations in) The derivaion for massless particles require some care though If you define the 4-momentum as $P^{\mu} = m \frac{dx^{\mu}}{d\tau}$, the geodesic equation becomes $P^{\alpha}(\partial_{\alpha}P^{\mu} + \Gamma_{\alpha \beta}^{\mu} P^{\beta}) = 0$ The bracketed term is the covariant derivative, denoted as $\nabla_{\alpha} P^{\mu}$ Substituting in the FRW metric to the geodesic equation tells you how the 4-momentum depends on the scale factor Looking at $\mu=0$, we see that: For massless particles, we have $\frac{1}{E} \frac{dE}{dt} = - \frac{\dot{a}}{a}$ This implies an energy scalinig of $E = \frac{1}{a}$ For massive particles, we have that the momentum scales like $\frac{1}{a}$ Redshift Since spacetime is getting stretched out, the wavelength of light gets stretched as well. This increase in observed wavelength is called redshifting Since energy of photons scales like $\frac{1}{a}$, the wavelength is proportional to a. Hence $\lambda_{0} = \frac{a(t_{0})}{a(t_{1})} \lambda_{1}$ redshift is defined as the fractional shift in wavelength: $z = \frac{\lambda_{0}-\lambda_{1}}{\lambda_{1}}$ Using $a(t_{0})=1$ for today’s scale factor, we have that $1+z = \frac{1}{a(t_{1})}$ Conventionally, we labe events in the universe’s history by their redshift: the surface of last scattering happens at z=1100 and the first galaxies formed around z= 10 For nearby sources (z< 1), we can expand the scale factor around $t_{0}$ and define $H_{0} = \frac{\dot{a_{0}}}{a(t_{0})}$ as Hubble’s constant to get that $a(t_{1}) = 1+ (t_{1}-t_{0})H_{0}$ For close objects, $\Delta t$ is simply $frac{d}{c}$, which means that redshift increases linearly with distance Alternatively, the recession speed $v = cz$ becomes $v \approx H_{0}d$ The above is called the Hubble-Lemaitre law It’s convention to describe Hubble constant measurements as $H_{0} = 100 h \frac{km}{s Mpc}$ since the measurements originally have very large uncertainties h is $0.73 \pm 0.01$ from supernovae measurements while h is $0.674 \pm 0.005$ from CMB measurements There is a statistically significant difference between the two, giving rise to the “Hubble tension” Distances These are hard to measure. The metric distance is unobservable, and the physical distance assumes a fixed time To make measurements, we need definitions which take into account the expansion and the finite travel time of light Luminosity Distance We can use “standard candles” (re: objects of known intrinsic brightness) and compare against their observed brightness to calculate distances We have the Cepheids, which have a periodic brightness. This period is correlated with the intrinsic brightness of a star More explicitly, the brightness of the star versus the logarithm of the period shows a linear relationship So if you can measure the distance to a Cepheid via parallax, you can extrapolate the distances to the rest of the Cepheids The next rung up on the cosmic distance ladder are the Type 1a supernovae These happen when a white drawf accretes too much matter from a companion star and explode They are so bright that they outshine the stars in their local group These explosions occur at a precise moment (re: when the white drawf exceeds the Chandrasekhar limit) and hence have a fixed brightness To generalize, suppose that we have a source of known luminosity L (re: energy per unit time) There is some observed flux (re: energy per unit time per unit area) from with L can be inferred Suppose that the source is at redshift z. The comoving distance is $\chi(z) = c\int_{0}^{z} \frac{dz}{H(z)}$ The flux in a static space is $F = \frac{L}{4 \pi \chi^{2}}$ since the source is isotropic the expansion complicates this for 3 reasons: the radius of the sphere expands so that the area becomes $4\pi a^{2}(t_0) d_{M}^{2}$ The arrival rate of the photons is smaller than the rate at which they are emitted by a factor of $\frac{1}{1+z}$ The energy of the photons get red shifted, which contributes another factor of $\frac{1}{1+z}$ The end result is that $F = \frac{L}{4\pi d_{M}^{2} (1+z)^{2}} = \frac{L}{4\pi d_{L}^{2}}$ where $d_{L}$ is the luminosity distance $d_{L} = (1+z) d_{M}(z)$ where $d_{M}$ is the metric distance Angular Diameter Distance Standard rulers are objects of know physical size An example is the typical size of the hot and cold spots in the CMB from theory Once again, assume a static space. The angular size of an object that is of known length D at some comoving distance $\chi$ away is $\delta \theta = \frac{D}{\chi}$ With expansion, the formula becomes $\delta \theta = \frac{D}{a(t_{1}) d_{M}} = \frac{D}{d_{A}}$ $d_{A}(z) = \frac{d_{M}}{1+z}$ relates the angular diameter distance to the metric distance Dynamics To calculate the scale factor, we need to solve the Einstein equation: $G_{\mu\nu} = \frac{8\pi G}{c^{4}} T_{\mu\nu}$ G is the Einstein tensor (curvature measurement) which T is the energy momentum tensor First, let’s look at number density $N^{\mu}$ $\mu=0$ is the number density of particles $\mu=i$ component is the flux of the particles in direction $x^{i}$ What constraints must $N^{\mu}$ obey to satisfy homogeneity and isotropy? Isotropy implies that $N^{i} = 0$ while homogeneity implies $N^{0} = c n(t)$ In a comoving frame, we have $N^{\mu} = n U^{\mu}$ Particle number should be conserved, which implies $\partial_{\mu} N^{\mu} = 0$ In curved spacetime, we have $\nabla_{\mu} N^{\mu} = 0$ $\nabla_{\mu} A^{\nu} = \partial_{\mu}A^{\nu} + \Gamma_{\mu\lambda}^{\nu} A^{\lambda}$ $\nabla_{\mu} B_{\nu} = \partial_{\mu}B_{\nu} - \Gamma_{\mu\nu}^{\lambda} B_{\lambda}$ So we have $\partial_{\mu} N^{\mu} = -\Lambda_{\mu \lambda}^{\mu} N^{\lambda)$ In the rest frame, we have $ \frac{\dot{n}}{n} = -3 \frac{\dot{a}}{a}$ This implies that $n(t) \propto a^{-3}$ In an analagous way to the number density, we can constrain the form of $T_{\mu\nu}$ We can break the tensor into 3 components: an energy density (scalar), a vector (momentum density/energy flux. Same values because of symmetry), and a smaller 3x3 stress tensor Isotropy forces the vector components to be 0 Homogeneity requires the energy density to be indpendent of position, but dependent on time Isotropy around the origin constrains the mean value of the stress tensor to be proportional to $\delta_{ij}$. Homogeneity requires that the proportionality constant be only a function of time All of that allows us to write down the perfect fluid stress-energy tensor: $T_{\mu\nu} = (\rho+\frac{P}{c^{2}})U_{\mu}U_{\nu}+P g_{\mu\nu}$ Local energy and momentum conservation imply $\nabla_{\mu} T^{\mu}{\nu} = 0$ Expand out in terms of Christoffel symbols and examine the $\nu=0$ term The $T^{i}_{0}$ term vanishes by isotropy In what remains, $\Gamma_{\mu 0}^{\lambda}$ is 0 unless $\lambda$ and $\mu$ are the same spatial indices, which implies $\Gamma_{i0}^{i} = \frac{3}{c} \frac{\dot{a}}{a}$ From the continuity equation, we arrive at: $\dot{\rho} + 3 \frac{\dot{a}}{a} (\rho + \frac{P}{c^{2}}) = 0$ Since the time translation is broken in an expanding space, global energy conservation doesn’t hold Most cosmological fluids are parametered as $P = w(\rho c^{2})$ This implies $\rho \propto a^{-3(1+w)}$ Matter is a fluid whose pressure is much smaller than it’s energy density. w=0, which implies $\rho \propto a^{-3}$ This could be baryons or dark matter Radiation is anything which obeys $P = \frac{1}{3} \rho c^{2}$. So $w=\frac{1}{3}$ and $\rho \propto a^{-4}$ This can be very light particles, photons, neutrinos, gravitons Dark energy has negative pressure, which implies $\rho \propto a^{0}$ Einstein Tensor for FRW metric The Einstein tensor is defined as $G_{\mu\nu} = R_{\mu\nu}-\frac{1}{2}R g_{\mu\nu}$ $R_{\mu\nu}$ is the Ricci tensor and $R = g^{\mu\nu}R_{\mu\nu}$ is the Ricci scalar You can write the Ricci tensor in terms of the Christoffel symbols: $\partial_{\lambda} \Gamma_{\mu\nu}^{\lambda} - \partial_{\nu} \Gamma_{\mu\lambda}^{\lambda} + \Gamma_{\lambda\rho}^{\lambda} \Gamma_{\mu\nu}^{\rho} - \Gamma_{\mu\lambda}^{\rho}\Gamma_{\nu\rho}^{\lambda}$ $R_{0i} = R_{i0} = 0$ because of isotropy (re: vector component) $R_{00} = -\frac{3}{c^{2}} \frac{\ddot{a}}{a}$ $R_{ij} = \frac{1}{c^{2}} (\frac{\ddot{a}}{a} + 2 (\frac{\dot{a}}{a})^{2} + 2 \frac{kc^{2}}{a^{2} R_{0}^{2}}) g_{ij}$ The author doesn’t directly calculate this. Instead, he computes $R_{ij}$ around x=0, then argues that since this is a tensorial reslationship, the equation holds everywhere The spatial metric is $\gamma_{ij} = \delta_{ij} + \frac{k x_{i} x_{j}}{R_{0}^{2}-k(x_{k}x^{k})}$ They expand to 2nd order since you need derivatives of $\gamma_{ij}$, and $\gamma_{ij}(x=0) = \delta_{ij}$ doesn’t have this information $\gamma_{jk}^{i} = \frac{k}{R_{0}^{2}} x^{i} \delta_{jk}$ The Ricci scalar is then $R = \frac{6}{c^{2}}(\frac{\ddot{a}}{a} + (\frac{\dot{a}}{a})^{2} +\frac{kc^{2}}{a^{2} R_{0}^{2}})$ Friedmann Equations Plugging in the FRW metric into the Einstein equations gives the Friedmann Equations The time-time component gives: $(\frac{\dot{a}}{a})^{2} = \frac{8\pi G}{3} \rho - \frac{kc^{2}}{a^{2}R_{0}^{2}}$ This is often written in terms of the hubble parameter $H = \frac{\dot{a}}{a}$ The first Friedmann equation can also be written as $\frac{H^{2}}{H_{0}^{2}} = \Omega_{r} a^{-4} + \Omega_{m} a^{-3} \Omega_{k} a^{-2} + \Omega_{\Lambda}$ Define $\Omega_{i} = \frac{\rho_{i}}{\rho_{crit}}$ where $\rho_{crit} = \frac{3H_{0}^{2}}{8 \pi G}$ os the density of a flat universe evaluated at today $\Omega_{\Lambda} = -\frac{kc^{2}}{(R_{0}H_{0})^{2}}$ If you evaluate at $t_{0}$, you get that $1 = \Omega_{r} + \Omega_{m} + \Omega_{\Lambda} + \Omega_{k}$ The spatical component(s) yield: $\frac{\ddot{a}}{a} = -\frac{4\pi G}{3}(\rho+\frac{3P}{c^{2}})$ The 2nd Friedmann equation is also called the Raychaudhuri equation It can also be derived by taking the time derivative of the 1st equation and applying continuity for $\dot{\rho}$ Exact Solutions of the Friedmann Equations In general, you need numerics to solve the Friedmann equation Single Component Universes Assume a flat universe with just a single component You can parameterize the specific component by its associated w factor We have $\frac{d \ln a}{dt} \approx H_{0} \sqrt{\Omega_{i}} a^{-\frac{3}{2}(1+w_{i})}$ Tranforming to conformal time, we find $a(\eta} \propto \eta^{\frac{2}{1+3 w_{i}}$ The special case of $w_{i} = 1/3$ is a universe dominated by spatial curvature A pure matter universe is a called a Einstein-de Sitter universe This is a good approximation for long matter-dominated periods in the universe Using a Matter dominated universe implies that the Hubble constant is $\frac{2}{3} \frac{1}{t_{0}}$. This gives a universe age of 9 billion years, which is the famous age problem A pure curvature universe is only permitted if k=-1. This is called a Milne universe* A pure curvature universe is only permitted if k=-1. This is called a Milne universe If $w=-1$, then we have a de Sitter space, which is a good approximation in the far past and far future All of these solutions (baring de Sitter space) have a singularity at t=0. This singularity (which goes like $\rho \propto t^{-2}$ is a generic feature of all Big Bang cosmologies (Hawking and Penrose) This assumes a strong energy condition ($\rho c^{2} + 3P \geq 0 $ which implies that $\ddot{a} \leq 0$ at all times Gravity is attractive, so it slows down the expansion Our Universe There are a handful of parameters that parameterize our universe First off is the temperature of the cosmic microwave background (CMB), which is given as $T_{0} = 2.7255 \pm 0.0006 K$ The energy density of the photons is given by $\Omega_{\gamma} = 5.4E-5$ The energy density of the neutrinos is 68% that of the photons This is currently being constrained by $0.0012 < \Omega_{\nu} < 0.003$, where the lower limit comes from neutrino oscillation experiments and the upper bound comes from cosmology All of these measurements came from COBE The anisotropy of the CMB places constrains on the curvature energy density ($|\Omega_{k}| < 0.005$ This implies that curvature is a very negligable part of the universe’s energy contents (and more so in the past) We have the baryon density, which we can infer from the abundances of light chemical elements produced in the Big Bang. This coems to around $\Omega_{b} \approox 0.05$ The dark matter density is $\Omega_{c} \approx 0.27$ Dark energy, responsible for the expansion of the universe, makes up $\Omega_{\Lambda} \approx 0.68$ of the energy content of the universe Integrating the Friedmann equation yields: $H_{0} t = =\int_{0}^{a} \frac{da}{\sqrt{\Omega_{r}a^{-2}+\Omega_{m}a^{-1}+\Omega_{\Lambda}a^{2}+\Omega_{k}}}$ We can infer the age of the universe (a=1) to be $t_{0} \approx 13.8$ GYrs We can evaluate the time scales at which matter-radiation and matter-dark energy equality occur at: $t_{mr} \approx 50000$ yrs $t_{m\Lambda} \approx 10.2$ Gyrs The coincidence problem: Why did dark energy come to dominate the expansion so close to the present time? The Hot Big Bang Natural Units We define natural units, which set the speed of light and $\hbar$ to unity We also define the reduced Planck mass as $M_{PI} = \sqrt{\frac{\hbar c}{8\pi G}}$ The Boltzmann constant is also typicall set to unity, equating temperature and mass Statistical Mechanics Refresher The probability of a particle having some momentum and position at time t is described by a distribution function From homogeneity, f can’t depend on x From isotropy, f can only depend on the magnitude of p In equillibrium, we have that the gas is in a state of maximum entropy and has the following distribution function $f(p,T} = \frac{1}{e^{\frac{E(p)-\mu}{T}}\pm 1}$ Where the + is for fermions and the - is for bosons The density of state is the number of particles assigned to each probability bin For a box normalized system, we have $\frac{g}{h^{3}} = \frac{g}{(2\pi)^{3}}$ This gives rise to: The number density: $n(T) = \frac{g}{(2\pi)^{3}} \int d^{3} p f(p,T)$ The energy density: $n(T) = \frac{g}{(2\pi)^{3}} \int d^{3} p E(p) f(p,T)$ $E(p) = \sqrt{m^{2}+p^{2}}$ The pressure: $n(T) = \frac{g}{(2\pi)^{3}} \int d^{3} p f(p,T)\frac{p^{2}}{3E(p)}$ Pressure is momentum change per unit time per unit area. The change in momentum is 2|p| along one axis. The volume swept out in unit time is $v_{x} dA = p_{x} \frac{dA}{E}$. We then calculate $\frac{p_{x}^{2}}{E}$ over the particles moving in the correct direction (divide by 2) and along the correct axis (divide by 3) For now, we assume that the chemical potential is much smaller than the temperature The Primordial Plasma Using a relativistic energy relationship: $E = \sqrt{p^{2}+m^{2}}$ and setting $\mu=0$, we have: $n = \frac{g}{2\pi^{2}} \int_{0}^{\infty} dp \frac{p^2}{exp(\frac{\sqrt{p^{2}+m^{2}}}{T})\pm 1 }$ $\rho = \frac{g}{2\pi^{2}} \int_{0}^{\infty} dp \frac{p^2 \sqrt{p^{2}+m^{2}}{exp(\frac{\sqrt{p^{2}+m^{2}}}{T})\pm 1 }$ We can define the dimensionless variables $x = \frac{m}{T}$ and $\epsilon = \frac{p}{T}$ The integrals become $I_{\pm}(x) = \int_{0}^{\infty} d\epsilon \frac{\epsilon^{2}}{exp(\sqrt{\epsilon^{2} + x^{2}}\pm 1}$ and $J_{\pm}(x) = \int_{0}^{\infty} d\epsilon \frac{\epsilon^{2}\sqrt{\epsilon^{2}+x^{2}}}{exp(\sqrt{\epsilon^{2} + x^{2}}\pm 1}$ Relativistic limit The momentum dominates the rest mass, so x approaches 0 $I_{\pm}(x=0) \int_{0}^{\infty} d\epsilon \frac{\epsilon^{2}}{e^{\epsilon}\pm 1}$ Expand the denominator in a geometric series For bosons (-) we get $I_{-}(0) = 2 \zeta(3)$ with $\zeta$ being the Reimann zeta function For fermions, we find that $I_{+}(0) = \frac{3}{4} I_{-}(0)$ A cute trick to see this is to note that $\frac{1}{e^{\epsilon}+1} = \frac{1}{e^{\epsilon}-1}-\frac{2}{e^{2\epsilon}-1}$ Hence we have $I_{+}(0) = I_{-}(0)- 2 * (\frac{1}{2})^{3}* I_{-}(0) = \frac{3}{4} I_{-}(0)$ Hence we get that $n = \frac{\zeta(3)}{\pi^{2}} g T^{3}$ and $\rho = \frac{\pi^{2}}{30} gT^{4}$ with a multiplicative factor of 1, $\frac{3}{4}$ and 1, $\frac{7}{8}$ for bosons and fermions respectively This tells use the number density and energy density of relic photons Setting $p = E$ yields $P = \frac{1}{3} \rho$, as expected of a relativistic gas The early universe consisted of a collection of different species, which makes the total energy density $\rho = \Sigma_{i} \frac{g_{i}}{2\pi^{2}} T_{i}^{4} J_{\pm}(x_{i})$ It’s standard to write the density in terms of the photon temperature: $\rho = \frac{\pi^{2}}{30} g_{*}(T) T^{4}$ $g_{*}(T) = \Sigma_{i} g_{i} (\frac{T_{i}}{T})^{4} \frac{J_{\pm}(x_{i})}{J_{0}}$ Calculating $g_{*}$ If we assume only relativistic species, we have that $g_{*} = \Sigma_{i=b} g_{i} (\frac{T_{i}}{T})^{4} + \frac{7}{8} \Sigma_{f} g_{i} (\frac{T_{i}}{T})^{4}$ where b and f represent summing over bosons and fermions respectively How many internal degrees of freedom are there? A massive particle has g = 2s+1, while a massless particle oof any spin has g=2 There are 4 force carriers (re: gauge bosons) The photon is massless (g=2) The weak force carriers are all spin 1 gauge bosons, and there are 3 of them (g=9) The Gluons are massless, and there are 8 of them (because 8 generators of SU(3)) (g=16) The leptons are massive spin $\frac{1}{2}$. (g=12, including the antiparticles) The quarks have 6 flavors and each come in 3 colors. (g=72, including antiparticles) The nutrinos, despite being massive spin $\frac{1}{2}$ particles only contribute 1 degree of freedom (only left-helicity neutrinos have ever been observed) Non-relativistic limit we let x» 1 which makes the integral the same for bosons and fermions $I(x) = \int_{0}^{\infty} d \epsilon \frac{\epsilon^{2}}{e^{\sqrt{\epsilon^{2}+x^{2}}}}$ Taylor expand the square root to 1st order and perform the Gaussian integral Rational is that most of the integral contributions come from $\epsilon « x$ $I_{\pm}(x) = \sqrt{\frac{\pi}{2}} x^{\frac{3}{2}} e^{-x}$ $n = g (\frac{mT}{2\pi})^{\frac{3}{2}} e^{-\frac{m}{T}}$ The exponential suppresses massive particles at these low temperatures (called Boltzmann suppression) Entropy Conservation Recall the first law of thermodynamics: $TdS = dU+PdV$ (dropped chemical potential since it’s small) Recalling that $U = \rho V$, define the entropy density $s = \frac{S}{V}$ Since s and $\rho$ are intrinsic quantities (re: volume independent), we we re write the first law as $(Ts-\rho-P) dV + V(T\frac{ds}{dT}-\frac{d\rho}{dT}) = 0$ Both bracketed terms must vanish seperately since dV and dT can be arbitrary so we have $s = \frac{\rho+P}{T}$ and $\frac{ds}{dT} = \frac{1}{T} \frac{d\rho}{dT}$ We can use the continuity equation $\frac{d\rho}{dt} = -3H (\rho+P) = -3H Ts$ to rewrite the above as $\frac{d(sa^{3})}{dt} = 0$ This implies that the total entropy is conserved in equillibrium and that the entropy density scales as $s \propto a^{-3}$ Cosmic Neutrino Background Since the neutrinos are the most weakly interacting particles in the Standard Model, we expect them to decouple from the thermal plasma first Neutrinos couple to the bath via processes like: $\nu_{e} + \bar{\nu_{e}} \leftrightarrow e^{+} + e^{-}$ $e^{-} + \bar{\nu_{e}} \leftrightarrow e^{-} + \bar{\nu_{e}}$ The interaction rate (per particle) is $\Gamma = n \sigma |v|$ where n is the number density of the target particle, $\sigma$ is the cross section and v is the relative velocity v is c is the relativistic limit We can get a rough scale of the decoupling temperature: We know that $\sigma \approx G_{F}^{2} T^{2}$ where $G_{F} \approx 1.2E-5 GeV^{-2}$ is the Fermi constant n scales as $T^{3}$ The interaction rate becomes $\Gamma \approx G_{F}^{2} T^{5}$ The hubble rate scales like $\frac{T^{2}}{M_{Pl}}$ These scales are roughly equal to each other when T is around 1 MeV After decoupling, the neutrinos preserve the relativistic Fermi-Dirac distribution along all the geodesics The number density scales like $n \propto a^{-3} \int d^{3}q \frac{1}{exp(\frac{q}{aT_{\nu}})+1}$ q = ap is the a time-independent momentum (needed since $p \propto a^{-1}$) particle number conservation requires $n_{\nu} \propto a^{-3}$. This is only consistent if $T_{\nu} \propto a^{-1}$ As long as photon temperature still scales like the neutrino temperature, we will be in equllibrium. Particle annihilation breaks this scaling Electron-Positron Annihilation Shortly after neutrino decoupling, the temperature drops below the electron mass, which means that the electrons and positrons annihilate with the photons This causes the photons to be “heated” (re: the photon bath temperature decreases more slowly compared to the neutrinos) We can characterize this by comparing the effective number of degrees of freedom. For $T<m_{e}$, we have g=2, while for $T \geq m_{e}$ we have $2+\frac{7}{8}\otimes 4$ electron-positron annihilation causes the $e_{\pm} \gamma$ plasma to evolve quasistatically into a $\gamma$ only plasma due to the timescale of annihilation This implies that the $T_{nu} = (\frac{4}{11})^{\frac{1}{3}} T_{\gamma}$ (Follows from conservation of entropy ($S = g(aT)^{3}$) Cosmic Microwave Background At temperatures of 1 eV, the universe consisted of a plasma of free electrons and nuclei. Photons tightly coupled to electrons via Thomson scattering, while Coulomb scattering coupled electrons to protons At temperatures below 0.3 eV, electrons and nuclei formed neutral atoms, which caused the density of free electrons to decrease sharply This caused the mean free path to quickly exceed the Hubble length. Around 0.25 eV, photons decoupled from matter enough to allow the photons to escape (re: the universe became transparent) Chemical Potential Photons have 0 chemical potential, since their number is not conserved Ex: double Compton scattering ($e^{-} + \gamma \rightarrow e^{-}+\gamma + \gamma$) The chemical potential of an antiparticle is $\mu_{X} = -\mu_{\bar{X}}$ Follows from particle-antiparticle annihilation in equillibrium with a photon bath Recombination Saha Equation The formation of hydrogen occurs via the reaction: $e^{-}+p^{+} \leftrightarrow H + \gamma$ In equillibrium, when $T< m_{i}$ where i could be e, p, or H, the have equillibrium abundances: $n_{i}^{eq} = g_{i}(\frac{m_{i}T}{2\pi})^{\frac{3}{2}} exp(\frac{\mu_{i}-m_{i}}{T})$ Using $\mu_{e}+\mu_{p} = \mu_{H}$, we can eliminate the chemical potentials: $\frac{n_{H}}{n_{e}{n_{p}}} = \frac{g_{H}}{g_{e}g_{p}} (\frac{m_{H}}{m_{e}{m_{p}}}\frac{2\pi}{T})^{\frac{3}{2}} exp(\frac{m_{p}+m_{e}-m_{H}}{T})$ Make the approximation that $m_{H} \approx m_{p}$ in the prefactor In the exponent, define $E_{I} = m_{p}+m_{e}-m_{H}$ as the ionization energy The number of degrees of freedom are 2, 2 and 4 for e,p and H respectively e and p come from being spin 1/2 H comes from the singlet and triplet states We believe the universe to overall be electrically neutral (ie. $n_{e} = n_{p}) All the above gives: $\frac{n_{H}}{n_{e}^{2}} = (\frac{2\pi}{m_{e}T})^{\frac{3}{2}} e^{\frac{E_{I}}{T}}$ Define $X_{e} = \frac{n_{e}}{n_{p}+n_{H}} = \frac{n_{e}}{n_{p}+n_{H}}$ to be the free-electron fraction (or the ionization fraction) Neglecting the small amount of helium, $n_{p} + n_{H}$ is approximately the baryon density, which we define as $n_{b} = \eta n_{\gamma} = \eta \frac{\zeta(3)}{\pi^{2}} T^{3}$ We can combine the above 2 into 1 equation: $\frac{1-X_{e}}{X_{e}^{2}} = \frac{n_{H}}{n_{e}^{2}} n_{b}$ Subtituting in the number density ratio yields the Saha equation The solution to the Saha equation is $X_{e} = \frac{-1+\sqrt{1+4f}}{2f}$ where $f(T,\eta) = \frac{2\zeta(3)}{\pi^{2}} \eta (\frac{2\pi T}{m_{e}})^{\frac{3}{2}} e^{\frac{E_{I}}{T}}$ define $T_{rec}$ as the temperature at which $X_{e} = 0.5$ The shape of the Saha equation looks like a sigmoid function w.r.t. T (or z via the conversion $T = T_{0}(1+z)$) You find out that for $\eta \approx 6E-10$ that at $z \approx 1300$ or so is when $X_{e} = 0.5$ Since matter-radiation equality happens are $z\approx 3400$, we know that recombination occurred in the matter-dominated era Photon Decoupling Early on, photons are strongly coupled to free electrons via Thompson scattering: $e^{-} + \gamma \leftrightarrow e^{-} + \gamma$ At some point , the expansion rate eclipsed the interaction rate, causing the photons to decouple Denote the temperature where these two rates are equal to be $T_{dec}$ $\Gamma_{\gamma}(T_{dec}) = n_{b} X_{e}(T_{dec}) \sigma_{T} = \frac{2\zeta(3)}{\pi^{2}} \eta \sigma_{T} X_{e}(T_{dec}) T_{dec}^{3}$ $H(T_{dec}) = H_{0} \sqrt{\Omega_{m}} (\frac{T_{dec}}{T_{0}})^{\frac{3}{2}}$ Set the rates equal to each other, plug in the Saha equation and cosmological parameters to find that $T_{dec} \approx 0.27 eV$ This places decoupling shortly after recombination, but this small change creates a drastic drop in the ionization fraction Last Scattering The precise moment of photon decoupling is called last-scattering We define the probability that a photon will scatter within some time interval to be $\Gamma_{\gamma}(t)dt$ The optical depth is the integrated probability: $\tau(t) = \int_{t}^{t_{0}} \Gamma_{\gamma}(t)dt$ Taking $t_{0}$ to be the present time, the moment of last scattering is defined to be $\tau(t_{*}) = 1$ To a good approximation, $t_{*} \approx t_{dec}$ However, the free electron density at the end of recombination greatly influences when $t_{*}$ happens Boltzmann Equation We want to describe evolution away from equillibrium. This requires the Boltzmann Equation We recall from FRW equations that $\frac{dn_{i}}{dt}+3\frac{\dot{a}}{a} n_{i} = \frac{1}{a^{3}}\frac{d(n_{i}a^{3})}{dt} = 0$ This implies that the number of particles in a fixed physical volume is conserved. Hence $n \propto a^{-3}$ due to expansion To include collisions, we set the RHS to $C_{i}([n_{j}])$. This is the Boltzmann equation More technically, it’s the integrated Boltzmann equation over the distribution functions We will only consider 2 particle scattering. More concretely, we will look at reactions of the form $1+2 \leftrightarrow 3+4$ For these reactions, the total overall reaction rate is given by $-\alpha n_{1} n_{2}+ \beta n_{3} n_{4}$ You lose reactants 1 and 2, and this rate is proportional to the number density of them both Similarly, increase the reaction rate is from 3 reacting with 4 $\alpha = <\sigma v>$ is the thermally-averaged cross section $\beta = (\frac{n_{1}n_{2}}{n_{3}n_{4}})_{eq} \alpha$ due to detailed balance (ie. at equillibrium, the forward and backward rates should equal each other) Given the above, it’s nice to write the Boltzmann Equation in terms of the number of particle in a comoving volume ($N_{i} \propto n_{i} a^{3}$) $\frac{d \ln N_{1}}{d \ln a} = -\frac{\Gamma_{1}}{H}(1-(\frac{N_{1,eq}N_{2,eq}}{N_{3,eq}N_{4,eq}}) \frac{N_{3}N_{4}}{N_{1}N_{2}})$ $\Gamma_{1} = n_{2} <\sigma v>$ is the interaction rate of species 1 Define $\frac{\Gamma_{1}}{H}$ as the interaction efficiency Dark Matter Freeze-Out Define the following interaction $X+ \bar{X} \leftrightarrow l + \bar{l}$ X is some heavy fermion and l is some light (basically massless) particle X could be dark matter We will assume that the light particle remains coupled to the plasma, meaning that they maintain their equillibrium densities ($n_{l} = n_{eq}$) We also assume that there is no initial asymmetry between X and $\bar{X}$ densities ($N_{X} == N_{\bar{X}}$) The Boltzmann equation becomes: $\frac{1}{a^{3}} \frac{d(n_{X}a^{3})}{dt} = -<\sigma v > (n_{X}^{2} - (n_{X}^{eq})^{2})$ Define $Y_{X} = \frac{n_{X}}{T^{3}}$ Define a new time $x = \frac{M_{X}}{T}$. $M_{X}$ denotes the temperature scale associated with the mass of particle X, which is when most of the dynamics happen at The Hubble parameter can then be written as $H = H \frac{M_{X}}{x^{2}}$ Given the above, you get the Riccati equation: $\frac{dY_{X}}{dx} = - \frac{\lambda}{x^{2}} (Y_{X}^{2} - (Y_{X}^{eq})^{2})$ $\lambda = \frac{\Gamma(M_{X})}{H(M_{X})}$ is dimensionless and is taken to be constant This is not analytically tractable, but numerically, you see that there is a plateauing of Y for large values of $x_{f}$ This is “freeze-out”: massive particles initially are in eq. with the plasma, but are eventually become too far apart to interact fast enough to stay in EQ There is a rough freeze out time $x_{f}$ at which this detailed balance gets broken For large enough times, $Y_{X}» Y_{X}^{eq}$. Hence the Boltzmann equation reduces to $\frac{dY_{X}}{dx} \approx - \lambda \frac{Y_{X}^{2}}{x^{2}}$. Integrate from $x_{f}$ to $x = \infty$, and realize that $Y_{X}^{f} » Y_{X}^{eq}$ gives $Y_{X}^{\infty} = \frac{x_{f}}{\lambda}$ Hence, the freeze out density is inversely proportional to the interaction rate WIMP Miracle The number density of particles after freeze out scales like $N_{X} \propto a^{-3}$ $n_{X,0} = n_{X,1} (frac{a_{1}}{a_{0}})^{3} = Y_{X}^{\infty} T_{0}^{3} (\frac{a_{1}T_{1}}{a_{0}T_{0}})^{3}$ $a_{1}$ is an arbitrary time that is late enough for the species to have reached their relic abundance, but before particle annihilation has broken the $T \propto a^{-1}$ scaling From conservation of entropy $G_{*S}(aT)^{3} = const$ we have $n_{X,0} = Y_{X}^{\infty}T_{0}^{3} \frac{G_{*S}(T_{0})}{G_{*S}(M_{X})}$ We can calculate the energy density of the particles to be $\Omega_{X} = \frac{\rho_{X,0}}{\rho_{crit,0}}$ Plugging in the observed values yields $\Omega_{X} \approx 0.1 \frac{x_{f}}{\sqrt{g_{*}(M_{X})}} \frac{1E-8 GeV^{-2}}{<\sigma v>}$ This is insensitive to the mass of the new particles The required cross section to see the observed dark matter density is about 0.1 $\sqrt{G_{F}}$ This is the WIMP miracle, since the weak interaction scale gives the right order of magnitude for the dark matter abundance observed today Baryogenesis You can apply the above Boltzmann Equation analysis to baryon annihilation The cross section scales like the inverse mass squared of the pion (nucleon interations effectively mediated by pions) This gives a density of baryons to be $\Omega_{b} 1E-11$, which is much smaller than the observed value We need a theory of baryogenesis to produce the abundance of baryons in the early universe Sakharov Conditions In 1967, Sakharov identified three necessary conditions that any theory of baryogenesis must satisfy Violation of baryon number: Since we start with a equal number of baryons (B=0), we expect that at present time $B\neq 0$) Standard model lagrangian conserves baryon number via guage symmetry. Baryon violation occurs in SM via “triangle anomaly” C and CP violation C is charge conjugation. If C was an exact symmetry, then $i\rightarrow f$ would have equal probability to $\bar{i} \rightarrow \bar{f}$, and thus there would be no net baryon number generated CP (charge conjugation with parity) is more subtle CPT (Charge, parity, time reversal) is a symmetry of any relativistic quantum field theory So CP invariance is the same as time invariance. If you have a process $i(p_{a},s_{a}) \rightarrow f(p_{a},s_{a})$, where p and s are the momenta and spins respectively, then under time reversal, you have $f(-p_{a}, -s_{a}) \rightarrow i(-p_{a}, -s_{a})$. Integrating over all momenta and spin yields a vanishing net baryon number. Hence, CP must be violated Both C and CP are violated in the weak interaction Departure from equillibrium: processes that generate the baryon asymmetry must occur out of thermal equillibrium Big Bang Nucleosynthesis (BBN) Qualitatively: In the beginning, baryonic matter was mostly protons and neutrons, which were coupled to each other via $\beta$ and inverse $\beta$ decay $n+ \nu_{e} \leftrightarrow p^{+} + e^{-}$ $n+e^{+} \leftrightarrow p^{+} + \bar{\nu}_{e}$ At 1 MeV, the neutrons decoupled from the bath. These free neutrons decay, further reducing their abundances free neutrons also combine with protons to form deuterium (~0.1 MeV): $n + p^{+} \leftrightarrow D + \gamma $ Deuterium nuclei combine with free protons to produce helium: $D+p^{+} \leftrightarrow ^{3}He+\gamma$ $D+^{3}He \leftrightarrow ^{4}He + p^{+}$ All the promordial neutrons get converted to helium heavy element formation difficult, which is why the protons survive and make hydrogen In detail: In equillibrium, the abundances go like $\frac{n_{n}}{n_{p}} = exp(-\frac{m_{n}-m_{p}}{kT})$ Assuming chemical potentials are negligable. The prefactor of $(\frac{m_{n}}{m_{p}})^{\frac{3}{2}}$ gets dropped b/c of similar masses, but not in the exponential Let $Q = m_{n}-m_{p}$ Neutron Freeze out occurs during decoupling at 1 MeV. The neutron fraction is $ X_{n} = \frac{n_{n}}{n_{n}+n_{p}}$ Can use equillibrium ratio to get that $X(T)_{eq} = \frac{e^{-\frac{Q}{T}}}{1+e^{-\frac{Q}{T}}}$ Use the Boltzmann equation, with 1 being neutrons, 3 being protons, 2 and 4 being leptons (which remain in equillibrium) You find that the neutron fraction obeys: $\frac{X_{n}}{dt} = -\Gamma_{n}(X_{n} -(1-X_{n})exp(-\frac{Q}{T}))$ $\Gamma_{n}$ comes from a cross section analysis Make the change of variables to $x = \frac{Q}{T}$ Integrate this with the initial condition of equillibrium of neutrons and protons Numerically, this traces out a flipped sigmoid curve. We find that $X_{n}$ approaches 0.15 for large x One wrinkle to this is the neutron life time The lifetime causes the final “equillibrium” to decay If fusion doesn’t happen fast enough, then in principle, the neutrons would vanish. Lucky for us that this isn’t the case Deuterium is the only efficient intermediate to large nuclei (ie. it’s the first stable one to form). This is the so called deuterium bottleneck neutron to neutron is very unstable proton to proton needs to overcome the Coulomb barrier 3 way interactions don’t really happen on these small time scales (100 seconds) After this point, all of the free neutrons get trapped into tritium and helium via the following reactions: $D + p^{+} \leftrightarrow He^{3} +\gamma$ $ D + D \leftrightarrow H^{3} + p^{+}$ $D+D \leftrightarrow He^{3} + n$ $H^{3} +p^{+} \leftrightarrow He^{3} + n$ $H^{3} +D \leftrightarrow He^{4} +n$ $He^{3} +D \leftrightarrow He^{3} + p^{+}$ Evantually, the final equillibrium helium abundance is half that of the neutron abundance Elements with A larger than 7 don’t exist in large quantities from BBN since their lifetimes aren’t long enough and there aren’t enough of them to construct a chain reaction Inflation The universe is largely homogeneous and flat, which would require a very fine tuning of the initial conditions The horizon problem: Why is the universe homogeneous everywhere you look? Why are causally seperated regions of the CMB roughly the same temperature? Flatness problem is: why is the universe so close to being flat currently when it was not initially? A priori, the universe can’t relax to flatness fast enough Why are fluctuations in the CMB correlated over distances larger than the light travel time since the Big Bang (superhorizon correlations)? Inflation, namely the rapid stretching of initial quantum fluctuations, gives a generic way to solve all of these problems that is agnostic to the initial conditions Horizon Problem Define the comoving particle horizon as $d_{h} = \int_{t_{i}}^{t} \frac{dt}{a(t)}$ This is the maximum possible distance which light (re: nill geodesic) can influence a future event The event horizon is the opposite: how far in the past can light influence an event First, rewrite the particle horizon as $d_{h} = \int_{\ln a_{i}}^{\ln a} (aH)^{-1} d(\ln a)$ We call $(aH)^{-1}$ the comoving Hubble radius. For large times, this is large, hence $d_{h} \approx (aH)^{-1}$ The recombination time is much smaller than the age of the universe today Imagine that you have two photons which start streaming at recombination time which are sufficiently far apart In standard Big Bang cosmology, this implies that there is not enough time for the two events to be in causal contact with each other (their light cones don’t interact) Then why is the CMB homogeneous? This implies that all of the sky was in causal contact with each other How large of a patch do you need for this to be a problem? Depends on your cosmology model. Find $\theta = \frac{2d_{h}(t_{rec})}{d_{A}(t_{rec})}$ where $d_{A}$ is the angular diameter distance between the events The Friedmann equations tell you what H is One caveat to the above: We observe that the conformal time blows up at the origin if a=0 For a small $\delta t$ after the Big Bang, we are assuming that the conformal time here doesn’t contribute significantly to the overall conformal time Namely, we are assuming that GR doesn’t break down at these small time scales. We know that it does at scales comparable to the Planck time A quantum gravity theory would elucidate the geometry in this area and potentially avoid this blow up Inflation provides an answer in lieu of a new quantum gravity theory We have assumed that the Hubble radius monotonically increases as a function of time. Doing this causes points to spread out too fast What if we demand that, in the early universe, that the Hubble radius decreases over time? Namely: $\frac{d}{dt}(aH)^{-1} <0$. Using $H = \frac{\dot{a}}{a}$, this is equivalent to saying that $\ddot{a} >0$ (ie. accelerated expansion) If we do this, we can allow large $d_{H}$ while maintaining a Hubble radius that allows causal contact Making the Hubble radius decrease as a function of time causes the $d_{H}$ integral to be dominated by the early Hubble radius Unlike Big Bang cosmology, inflation allows allows the singularity to occur at a negative conformal time ($\eta = -\infty$), hence providing more spacetime for events to be causally connected Flatness Problem $\Omega_{k}(t) = \frac{\rho_{c}-\rho}{\rho_{c}} = (\frac{a_{0}H}{aH})^{2} \Omega_{k,0}$ This is the time dependent curvature parameter, where $\rho_{c} = 3M_{p}^{2}H^{2}$ is the time dependent critical density This denotes how curved the universe is as a function of time Since the Hubble radius $\frac{1}{aH}$ increases as a function of time, this implies that the curvature of the universe was much smaller in the past However, we observe that the universe is extremely flat. This implies that the early universe was extremely flat. Why is this? A shrinking Hubble radius in the early universe would explain the apparent extreme flatness: any initial curvature would be quickly flattened out, after which the curvature would very slowly increase Superhorizon Correlations The universe contains density fluctuations which are correlated over acausal distances (superhorizon) " These detailed correlations make it more difficult to appeal to a creator without sounding like a young Earth creationist, arguing that the fossil record was planted to deceive us” - David Tong Inflation solves this by letting apparently causally connected regions interact, then quickly seperate to their acausal locations Inflation Duration Of course, inflation needs to last long enough so that the slow ramp up matches the observed values We need the Hubble radius of the observable universe $(a_{0}H_{0})^{-1}$ to be smaller than the Hubble radius during inflation $(a_i H_i)^{-1}$ Define $ N = \ln (\frac{a_{e}}{a_i})$ be the number of e-foldings (think how many doublings of space happen, but instead of base 2 it’s base e) Assume that the Hubble rate during inflation stayed roughly constant ($H_{i} \approx H_{e}$) Inflation goes up until the Hot Big Bang, after which the Hubble radius grows to the current value Assume that the Hubble radius doesn’t change significantly between the end of inflation and the start of the Big Bang (ie. $(a_{e}H_{e})^{-1} \approx (a_{R}H_{R})^{-1}$) Assume a radiation dominated universe for simplicity. We see that $\frac{a_{0}H_{0}}{a_{R}H_{R}} = \frac{a_{R}}{a_{0}} \approx \frac{T_{0}}{T_{R}}$ $T_{R}$ is the reheating temperature Combining all of the above gives bounds on how many e-foldings you need Structure Formation All of the prior stuff was assumed a homogeneous universe. We now sprinkle in inhomogeneities and see how they evolve over time We first look at the Newtonian limit to build some intuition as to what happens non-relativistic velocities low gravity (curvature much smaller than the Hubble radius) P « $\rho$ Fluid Dynamics We have three equations which govern fluid flow: $\frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \vec{u}) = 0$ Follows from mass conservation $\rho \frac{d\vec{u}}{dt} = \rho(\frac{\partial}{\partial t} + \vec{u} \cdot \nabla ) \vec{u} = - \nabla P$ Follows from momentum conservation $P = P(\rho,T)$ is an equation of state which closes the system Imagine that we have a static fluid (u=0). This implies that $\rho$ and P and constant. What would happen if we peturb from this steady state? let $\rho = \bar{\rho}+\delta p$ and similarly for P Make this substitution and only keep terms up to 1st order to get: $\partial_{t} \delta \rho = - \nabla \cdot (\bar{\rho} \vec{u})$ $\bar{\rho}\partial_{t} \vec{u} = -\nabla \delta P$ Combining the above yields: $\partial_{t}^{2} \delta \rho - \nabla^{2} \delta P = 0$ For now, assume a barytropic fluid: $P = P(\rho)$. We can then write $\delta P = \frac{\partial P}{\partial \rho} \delta \rho = c_{s}^{2} \delta \rho$ $c_{s}^{2} = \frac{\partial P}{\partial \rho}$ is the sound speed With the above substitution, the above becomes a wave equation We can write $\delta \rho(x,t)$ as the sum of all the Fourier modes: $\delta \rho(x,t) = \int \frac{d^{3} k } {(2\pi)^{3}} exp(i k\cdot x ) \delta \rho (k,t)$ This then gets converted to a ODE of k , with solutions : $\delta(k,t) = C_{k} exp(i\omega_{k}t)+D_{k}exp(-i\omega_{k}t)$ where $\omega_{k} = c_{s}|k|$ Newtonian Fluid Dynamics with Gravity Keep the mass continuity the same, but now add a gravitational potential: $\rho \frac{d\vec{u}}{dt} = \rho(\frac{\partial}{\partial t} + \vec{u} \cdot \nabla ) \vec{u} = - \nabla P - \rho \nabla \Phi$ Linearize in an analagous fashion to the previous section to see that $(\frac{\partial^{2} }{\partial t^{2}} + c_{s}^{2}k^{2} -4 \pi G \bar{\rho})\delta \rho (k,t) = 0$ $\omega_{k} = c_{s}^{2}k^{2} - 4\pi G \bar{\rho}$ is the new frequency of the system. Setting this to 0 defines a critical wavenumber called the Jeans scale: $k_{J} = \frac{\sqrt{4\pi G\bar{\rho}}}{c_{s}}$ for k»$k_{J}$, we recover the previous behavior for k « $k_{J}$, we have that gravity dominates. The solutions look like $\delta \rho \propto exp(\pm \frac{t}{\tau})$ where $\tau = \frac{1}{\sqrt{4\pi G \bar{\rho}}}$. These exponential solutions are called Jean instabilities Adding Expansion Now, we imagine that our coordinates change over time: $r(t) = a(t) x$ $u(t) = \dot{r} = Hr+v$ where H is the Hubble flow and v is the physical peculiar velocity “peculiar velocity” is just relative velocity w.r.t. some baseline (In this case, the Hubble flow) Originally, the time and space derivatives were independent. Due to expansion, we now get: $\frac{\partial f}{\partial t_{r}} = ((\frac{\partial}{\partial t_{x}}) - H x \cdot \nabla_{x}) f$ $t_{i}$ refers to what variable is being held constant Chain rule $f(t,a(t))$ and transform from r to x Previously, we implicitly held r constant in the time derivatives for the continuity equation, momentum equation Now, we need to hold x (the comoving coordinate) constant instead The continuity and momentum equations are derived by subbing in the derivative relationship b/w x and r above, and then using the relationship for the velocity field u The Poisson equation picks up a factor of $a^{2}$ on the RHS Define the background solution to be when $\vec{v} = 0$. The equations become: $\frac{\partial \bar{\rho}}{\partial t} = -3H \bar{\rho}$ $\frac{\partial \vec{\bar{u}}}{\partial t} = -\frac{1}{a}\frac{\nabla \bar{P}}{\bar{\rho}} - \frac{1}{a} \nabla \bar{\Phi}$ The solution follows that of a homogeneous non-relativistic fluid $\bar{\rho} \propto a^{-3}$ Solutions to Euler equation take the form $\vec{u} = Ha \vec{x}$, $\bar{P} = const$, $\nabla \bar{\Phi} = -\ddot{a} a \vec{x}$ This solution holds as long as the scale factor a satisfies $\frac{\ddot{a}}{a} = -\frac{4\pi G}{3} \bar{\rho}$ This this the acceleration equation derived from the Friedmann equations Define the linear peturbations: $\rho = \bar{\rho}(1+\delta)$ $\delta = \frac{\delta \rho} {\bar{\rho}}$ $\vec{u} = Ha\vec{x} + \vec{v}$ The velocity $\vec{v}$ is the expansion $P = \bar{P} + \delta P$ $\Phi = \bar{\Phi} + \delta \Phi$ Doing the expansion, and keeping terms up to first order, we find: $\dot{\delta} = -\frac{1}{a} \nabla \cdot \vec{v}$ $\dot{v} + Hv = -\frac{1}{a \bar{\rho}} \nabla \delta P - \frac{1}{a} \nabla \delta \Phi$ Recall that $(\vec{v} \cdot \nabla) \vec{x} = \vec{v}$ Combine the above two into a single 2nd order equation Take time derivative of first, plug $\dot{v}$ from first and $\nabla \cdot v$ from second, then the perturbed Possion equation ($\nabla^{2} \delta \Phi = 4\pi G a^{2} \bar{\rho} \delta$, then the peturbed pressure $\delta P = c_{s}^{2} \bar{\rho} \delta$ (remember we are assuming a barotropic fluid!) You get $\ddot{\delta} + 2H\dot{\delta} - (\frac{c_{s}^{2}}{a^{2}}\nabla^{2} + 4\pi G \bar{\rho}(t)) \delta = 0$ This describes the growth of density flunctuations in an expanding universe Growth of Matter Perturbations Make the standard Fourier expansion ansatz: $\delta \rho(x,t) = \int \frac{d^{3}k}{(2\pi)^{3}} exp(i\vec{k}\cdot \vec{x}) \delta \rho(\vec{k},t)$ In Fourier space, we have: $\ddot{\delta}(\vec{k},t) + 2H \dot{\delta}(\vec{k},t) + c_{s}^{2} (\frac{k^{2}}{a^{2}} - k_{J}^{2}) \delta(\vec{k},t) = 0$ The physical Jeans wavenumber is $k_{J}(t) = \frac{\sqrt{4\pi G \bar{\rho}(t)}}{c_{s}(t)}$ It’s time dependent now! On small scales ($\frac{k}{a} » k_{J}$) you get a damped harmonic oscillator, with the hubble expansion acting as friction On large scales ($\frac{k}{a} « k_{J}$), you recover the static space case, but the extra friction term prevents the fluctuations from growing exponentially Explicitly: $\ddot{\delta} + 2H \dot{\delta} - 4\pi G \bar{\rho}(t) \delta = 0$ where $\delta$ is a function of k and t Hence, Jeans wavenumber sets the scale where you go from suppression at small k to damped growth at large k The time dependent nature of the Jean scale is important For baryons, prior to recombination, the baryons are strongly coupled to the photons, hence have a sound speed like the photons ($c_{s} \approx \frac{1}{\sqrt{3}}$ which gives a Jeans scale on the order of a Hubble radius. Hence, there is no growth in baryon fluctuations on subhorizon scales Post recombination, the sound speed drops and flunctuations start to grow For cold dark matter, the sound speed is negligable, hence subhorizon modes can grow earlier Linear Growth Function Since for the large scale case, the ODE involves only time derivatives, we can write an ansatz of the form: $\delta(k,t) = \delta_{+}(k) D_{+}(t) + \delta_{-}(k) D_{-}(t)$ Intuition: Need 2 independent functions to span the solution space, each of which only depends on time. Each of these functions can be modulated by some function of k (ie. The fourier modes) $D_{+}(t)$ is called the linear growth function and is normalized such that $D_{+}(t_{0}) = 1$ plugging in this ansatz to the density flunctuation equation we can determine how the flunctuations grow as a function of the scale factor If we have multiple components, then the density becomes $\Sigma_{\alpha} \bar{\rho_{\alpha}} \delta_{\alpha}$ Statistical Properties Correlation Functions By definition, the mean value of the density perturbations is zero ($<\delta> = 0$) Hence, we need to use the two-point correlation function to get any statistics on the density field $\Epsilon(|x-x’|, t) = <\delta(x,t)\delta(x’,t)> = \int D\delta P(\delta) \delta(x,t) \delta (x’,t)$ This is a functional integral over field configurations. $P(\delta)$ is the probability of realizing some config $\delta(x,t)$ $<…>$ creates an emsemble average of the stochastic process Homogeneity implies that the two-point function should be invariant under translations (ie. only depends on separation between two points) Isotropy implies that only the distance between separation matters We want to compute this two-point function, but it is not what we observe The observed universe is a single realization of a random process We assume ergoicity (ie. ensemble averages are equal to spatial averages as the volume approaches infinity) This allows different parts of the universe to be seen as different realization of the underlying random process Somce we have a finite volume, this introduces statistical fluctuations called sample variance We can also write the two point function in Fourier space: $<\delta(\vec{k}), \delta(\vec{k’})> = \int d^{3}x d^{3} x’ exp(-i\vec{k}\cdot \vec{x}) exp(-i\vec{k’}\cdot \vec{x’}) <\delta(x)\delta(x’)>$ This can be simplified: $\int d^{3}r d^{3}x’ exp(-i\vec{k}\cdot \vec{r}) exp(-i(\vec{k}-\vec{k’})\cdot x’) \Epsilon(r) = (2\pi)^{3} \delta_{D}(\vec{k}-\vec{k’}) P(k)$ $P(k) = \int d^{3}r exp(-i\vec{k}\cdot\vec{r}) \Epsilon(r)$ The delta function in k implies that each wavevector is independent of each other $P(k)$ is the power spectrum (ie. the 3D Fourier transform of the correlation function $\Epsilon(r)$ Due to rotational inveriance, we know the power spectrum only depends on the magnitude of the wavevector Hence, working in spherical coordinates: $P(k) = \frac{4\pi}{k} \int_{0}^{\infty} dr r\sin(kr) \Epsilon(r)$ Gaussian Random Fields For a Gaussian random field, the probability distribution is a Gaussian function of $\delta(x)$ the PDF for N points in space $x_1$ to $x_{N}$ iS $P(\delta_{1},…\delta_{N}) \propto \frac{1}{\sqrt{det(\epsilon_{ij})}} exp(-\delta_{i} \epsilon_{ij}^{-1} \delta_{j})$ where $\delta_{i} = \delta(x_{i})$ and $\epsilon_{ij} = \epsilon(|i-j|)$ Harrison Zeldovich Spectrum $P(k,t) = T^{2}(k) \frac{D_{+}^{2}(t)}{D_{+}(t_{i})} P(k,t_{i})$ (will be derived later on) The transfer function is defined as $T(k) = \frac{D_{+}(t_{i})}{D_{+}(t)} \frac{\delta(\vec{k},t)}{\delta(\vec{k},t)}$ $t_{i}$ is some initial time which is taken to be after inflation This is necessary since the dependency inside and outside the Hubble radius is different Harrison Zeldovich and Peebles argued that $P(k,t_{i}) = Ak^{n}$ We expect the gravitational potential to be scale invariant during inflation If there was some characteristic scale, then you would observe clumps of energy during inflation In addition, as inflation happens, larger and larger scales enter the Hubble radius. We need the magnitude of each new scale to be the same as the previous one, hence scale invariance From Fourier transforming the Poisson equation, we see that $P_{\Phi} \propto k^{-4} P_{\delta} \propto k^{n-4}$ To have the gravitational potential to be scale free, we need the correlation function/variance to not depend on the wavenumber We know that $\sigma_{\Phi}^{2} = \Epsilon_{\Phi}(0) = \int \frac{dk}{k} \frac{k^{3}}{2\pi^{2}} P_{\Phi} $ Hence, to keep scale invariance, we need n=1

Following the book “Understanding Digital Signal Processing” by Richard G. Lyons. Discrete Sequences and Systems Definitions and Notation Signal Processing: the science of analyzing time-varying physical processes Analog refers to waveforms that are continuous in time and can take on a continuous range o amplitude values Comes from analog computers, where physical electronic differentiators and integrators were used to model differential equations Discrete: the time variable is quantized This means that the signal at intermediate values is unknown. You can make pretty good guesses at intermediate values, but ultimately, you don’t know A discrete system is one which takes in some input stream(s) $x_{i}(n)$ and returns some output streams $y_{j}(n)$ where n denotes which discrete time step you’re dealing with A fundamental difference of discrete systems from analog systems is the sampling rate $f_{s}$, which is typically independent of the signal you’re sampling An alternative representation of any discrete signal is in the frequency basis via Fourier transforms. Hence, $x(n) \rightarrow X(m)$ where m indexes the frequency (m=0 is the DC component) Basic signal processing operation symbols Decibel Fun history lesson: The original power difference was in units of “bels”, where the power difference is $10 * log_{10}(\frac{P_{1}}{P_{2}})$ bels This was inconvenient, since you routinely measured power differences smaller than that. so the decibel was introduced to account for this (just multiply by 10!) The log scale lets us compare signals of wildly different orders of magnitude on the same plot Linear and Time-Invariant Systems A linear System obeys the homogeneity property: $c_{1} x_{1}(n)+c_{2} x_{2}(n) \rightarrow c_{1} y_{1}(n)+c_{2} y_{2}(n) $ This system is also commutative A time-invariant system is one where a time delay/shift in the input sequence causes an equivalent time delay in the output sequence. To wit: $x(n+k) \rightarrow y(n+k)$ Given a time-invariant linear system (LTI), we can determine the output response via a convolution of the input sequence with system response to a unit impulse signal (ie. a signal where $x(0) = 1$ and all other times are 0) The frequency response of a LTI system via a discrete Fourier transform of the impulse response Periodic Sampling Frequency Ambiguity When sampling at a rate of $f_{s}$ samples/second, if k is an integer, you can’t distinguish sampled values of a sinewave of $f_{0}$ Hz and a sinewave of $f_{0}+kf_{s}$ Hz Suppose that you have a Shark’s tooth frequency response The region of interest is based around 0, with a spread of $f_{s}$ Outside this region in both directions, the band of interest gets repeated Because of this copying, any signal you measure in the band of interest cannot be differentiated by an equivalent signal in some other band This arises due to convolution Low Pass Aliasing This frequency ambiguity shows up when dealing with low pass signals (ie. the spectral energy is all below some cutoff frequency B) The above demonstrates the problem graph A represents the continuous incoming signal graph B represents the sampled spectra when the sampling rate is $f_{s} \geq 2B$ (so called Nyquist criterion). All of the frequency bands are nicely separated from each other graph C represents the sampled spectra when the sampling rate is less than the Nyquist criterion. As you can see, you get crossing of bands, which corrupts the initial signal You won’t know the amplitude in the cross over regions a priori All of the frequency content outside of the band gets folded into it when sampling (regardless of Nyquist criterion) (see above) This necessitates analog low pass filtering prior to sampling BandPass Filtering Sometimes, we care more about the bandwidth of a signal (ie. the range of frequencies that we can measure) compared to the absolute frequency that we can measure As such, we want to be able to sample some range of frequencies around some non-zero $f_{c}$ such that the result is centered around 0 Band Pass Aliasing The above shows how aliasing can be a problem with band pass signals We want to downshift the band by $f_{c}-\frac{B}{2}$ at some sampling frequency $f_{s}$. This means that we need to close a distance of $2f_{c}-B$ between the upper edge of the negative frequencies and the lower edge of the positive frequencies Because of the spectral replication, we can fit $m*f_{s}$ in this gap, where m is some positive integer Hence: $m*f_{s} = 2f_{c}-B$ or $f_{s} = \frac{2f_{c}-B}{m}$ Increasing the frequency past this point for fixed m will shift the P band right and the Q band left, causing aliasing. Hence, $\frac{2f_{c}-B}{m}$ gives an upper bound on $f_{s}$ What happens if you decrease $f_{s}$? In Fig. b, the P lobe gets shifted down, and the Q lobe gets shifted up. Hence, there is some lower bound where the bands start to overlap again. This occurs when the lower edge of the negative frequencies and the upper edge of the positive frequencies match. This distance is $2f_{s}+B$ Because you’re shifting the other way, you need to add one more band in the mix. Hence, the lower bound is $f_{s} = \frac{2f_{c}+B}{m+1}$ Hence we have: $\frac{2f_{c}-B}{m} \geq f_{s} \geq \frac{2f_{c}+B}{m+1}$ where m being an arbitrary positive integer constrains $f_{s} \geq 2B$ This is called 1st-order sampling Spectral Inversion If you choose m to be an odd integer, then you can get spectral inversion, where the positive and negative frequency slopes get flip-flopped If this is a problem, simply multiply the the discrete time samples by $(-1)^{n}$ This has the effect of flipping the the bands around $\pm \frac{f_{s}}{4}$ (with spectral replication of course) Positioning Sampled Spectra at $\frac{f_{s}}{4}$ For use cases like spectral inversion, we can choose $f_{s} = \frac{4f_{c}}{2k-1}$ so that the spectra is centered at $\frac{f_{s}}{4}$ Noise in Bandpass Sampling Suppose that the signal has some noise outside the region of interest. The power of the noise of the sampled signal increases by a factor of (m+1) while the signal power remains the same. Hence, the Signal to Noise Ratio (SNR) decreases by this same factor DFTs The continuous Fourier transform has the form $X(f) = \int_{-\infty}^{\infty} x(t) exp(-i 2\pi f t) dt$ The discretized version looks like $X(m) = \Sigma_{n=0}^{N-1} x(n) exp(\frac{-i 2\pi nm}{N})$ Remember that each m can be mapped to some frequency via $f(m) = \frac{m f_{s}}{N}$ where m ranges from 0 to N-1 You can think of each X(m) as a point-for-point multiplication between the input values and some complex sinusoid The X(0) term simplifies to $X(0) = \Sigma_{n=0}^{N-1} x(n)$, which is just proportional to the average. This is the DC bias of the signal Occasionally, you will get some normalization term in front of the DFT The convention used here is to have a factor of $\frac{1}{N}$ on the inverse Fourier Transform Some code uses a factor of $\frac{1}{\sqrt{N}}$ for both the forward and inverse Fourier transform. The reasoning behind this is so that there is no scale change when transforming in either direction DFT Properties Most signals that we care about are real. Under this assumption: $X(m) = X^{*}(N-m) If x(n) is even, then X(m) is real and even if x(n) is odd, then X(m) is imaginary and odd DFTs are linear transformations If you shift x(n) by k samples ,then the resulting DFT is shifted by an associated phase More concretely: $x(n) \rightarrow x(n+k)$ implies that $X_{shifted}(m) = exp(\frac{i 2\pi k m }{N}) X(m)$ DFT Leakage What happens if your signal doesn’t exactly equal a harmonic of $\frac{f_{s}}{N}$? The power of that frequency gets spread across all of the spectra (or all of the bins) This leakage means that bins which contain low-level signals get corrupted by high-amplitude signals in neighboring bins Another thing to note is that is leakage is periodic (ie. the m=0 bin can leak onto the N-1 bin) Windows Leakage occurs due to the sidebands of the Fourier transform of the signal One way to mitigate this is via windowing For each datapoint in your Fourier transform, you add some weight to each sample Hence, the weighted Fourier transform becomes: $X_{w}(m) = \Sigma_{n=0}^{N-1} w(n) x(n) exp(\frac{-i 2\pi nm}{N})$ You pick your weights such that the beginning and end of the sample interval both go smoothly to some common amplitude value There are a bunch of different weighting schemes: Boxcar/uniform/rectangular: $w(n) = 1$ for all n Triangular: $w(n) = \frac{n}{N/2}$ for n from 0 to N/2 and $w(n) = 2-\frac{n}{N/2}$ for N between N/2+1 to N-1 Hanning: $w(n) = 0.5-0.5\cos(\frac{2\pi n}{N})$ for all n Hamming: $w(n) The standard Fourier transform uses boxcar. This results in a Fourier transform of sinc(x). The side lobes of sinc are due to the sharp edges of the sampling window. The alternative windowing functions try to smoothly taper off to 0 instead The trade-offs for having a reduced sideband response is that the mainband gets stretched out in frequency space Scalloping For simplicity, assume no windowing (re: boxcar). Each Fourier bin gets transformed via sinc(x) The overlap between adjacent bins creates this ripple in the overall magnitude of the DFT away from the bin centers (see below for visual) In practice, this apparently isn’t a big problem, since most applications have so many frequency bins that this ripple is not noticeable Zero-Padding The idea behind this is to inflate the number of sampled points by adding zeros to the end of a input This is typically done up to the nearest power of 2 for FFT reasons The net result of this is to get a better approximation of the CFT One side effect is that the index associated with each point gets changed as you zero pad This is not a bad thing, since the associated frequency of that bin is still the same as the original One gotcha that you might encounter is combining zero-padding and windowing. Make sure that your window doesn’t touch the zero-padding! DFT Processing Gain Can think of each DFT bin as a bandpass filter centered around some frequency This bandpass effect can allow one to extract single frequencies from some background noise As the number of points increases, the signal frequency rises above the noise Can define a signal-to-noise ratio (SNR) which is the amplitude of the signal over the average background noise Increasing the number of points tends to increase the SNR Can’t increase ad infinitum because the number of multiplications scale as $N^{2}$ An alternative approach is to average a smaller number of DFTs DFTs in Practice Use some FFT implementation instead of the naive DFT Cooley-Tukey’s radix-2 is the classic version Make sure that you same fast enough and long enough fast enough means meeting Nyquist criterion at minimum (typical applications aim for 2.5x to 4x the bandwidth) If you don’t know the input frequency, watch out for FFT results which have large spectral components near half the sampling rate (re: aliasing) The number of samples you need to collect is given by $N = \frac{f_s}{f_{r}}$ where $f_{r}$ is your frequency resolution Make sure to round up to the nearest power of 2 for FFT reasons Make sure to zero pad your data as needed If you’re zero padding, make sure that you do this after windowing your signal Leakage is still a problem even if you zero pad the DC bias can affect your upper spectrum due to warping To mitigate this, can do an averaging of the time sequence and subtract of this average, removing the DC bias Make sure to do this before windowing If you need to improve SNR, try taking multiple small FFTs and averaging them Apparently, this lets you detect signals which are buried in the noise? Finite Impulse Response (FIR) Filters Given a finite duration of nonzero input values, an FIR filter wil eventually output a finite duration of of nonzero output values Conversely, if the input to an FIR filter suddenly became all zeros, the output would eventually become all zeros Can think of averaging as a simple example of an FIR filter that demonstrates the above: If the input signal suddenly drops to zero, the average will gradually taper off to zero Can be thought of as a kind of low pass filter The “delay” modules can be though of as shift registers which temporarily store the values for the computation Due to these delay modules, any change in the input takes some time to be realized in the output. The intermediate output of the filter is called the “transient” period FIRs as Convolution Picture a general M-tap FIR filter. The nth output is then $y(n) = \Sigma_{k=0}^{M-1} h(k) x(n-k)$ where $h(k)$ is the weight which gets applies to the kth input. This is the definition of a convolution Using the convolution theorem, we can equal convolution in the time domain to multiplication in the frequency domain (ie $h(n) * x(n) = H(m) \times X(m)$, where * denotes the convolution) One way to determine h(k) is to send an impulse response as input, then see what the filter outputs Can take the Fourier transform of the above to get H(m) Lowpass FIR Filter Design What if we want to go the other way around? Given some cutoff frequency, what coefficients do we need to implement a lowpass FIR filter? There are two methods to figure out the coefficients: window design, and optimum methods For simplicity, assume a perfect low pass filter: unity gain below the threshold frequency, and zero gain above the threshold Window Design Take the inverse Fourier transform of your desired frequency response in order to get the coefficients Some nuances: you want to define the frequency response from 0 to $f_{s}$ instead of the symmetric $\mp \frac{f_{s}}{2}$ When you do this process, you get some ripple power in your reconstructed filter due to Gibb’s phenomena This arises because no finite set of sinuoids will be able to change fast enough to be equal to an instantaneous discontinuity If you’re doing windowing, there is a trade off between ripple power suppression and main lobe widening that you need to figure out Two windows which allow a fine-tuning of this trade off are the Chebyshev window and the Kaiser window Bandpass and High Pass Design Due to the negative frequencies, a “low pass” filter is roughly symmetric around the origin If you could shift this band up so that it’s centered at $\frac{f_{s}}{4}$, then you would get a bandpass filter around $\frac{f_{s}}{4}$! This can be achieved by multiplying you low pass filter coefficients by a sinusoid of frequency $\frac{f_{s}}{4}$ This takes the form of a sequence ${-1,0,1,0,-1,…}$ A high pass filter is analogous, but you instead multiply by by a sinusoid of frequency $\frac{f_{s}}{2}$ This takes the form of a sequence ${-1,1,-1,…}$ Parks-McClellan Exchange FIR Filter Design Method (Optimum Method) In one sentence, Parks-McClellan tries to minimize the error in the pass and stop bands by utilizing the Chebyshev approximation It’s the defacto FIR filter generation algorithm and can be used for most realistic response curves We define $f_{p}$ and $f_{s}$ as the upper frequency bound on the low pass side and the lower frequency bound on the stop band side We also define $\delta_{p}$ and $\delta_{s}$ as the ripple deviation in the low pass and stop bands Half-Band Filters The basic idea here is that we want the frequency response to be symmetric around $\frac{f_{s}}{4}$ This implies that the average of $f_{p}$ and $f{s}$ is $\frac{f_{s}}{4}$ If the number of coefficients is odd, then every other coefficient is zero (excluding the central tap) This halves the number of multiplications needed to do the filter The number of taps N must be such that N+1 is a multiple of 4. Otherwise, you reduce down to a smaller effective number of taps You can use the Optimum Method to generate Half-Band filters of this type FIR Filter Phase Considerations We need to introduce the idea of group delay. Think of it as $G = -\frac{d\phi}{df}$ where $\phi$ is the phase We want G constant for amplitude modulated (AM) signals Why? Because if G is constant, all of the frequency components get delayed by the same amount. This ensure that the shape of the incident wave packet is conserved There are two ways in which discontinuities can show up in a phase response plot the range of arctan is $-\frac{\pi}{2}$ to $\frac{\pi}{2}$. Lots of software will force results to lie within this region. This is not physical Remember that the phase response function is complex. For each frequency component of the response, if the polarity of one component changes while the other remains the same, then a jump discontinuity will occur. This is physical Infinite Impulse Response (IIR) Filters IIR Filters allow feedback from their output Think of FIR filters as Moore FSMs, and IIR filters as Mealy FSMs The “Infinite” comes from the fact that you can end up in a state where the filter continuously outputs stuff even if there is no input Why use IIRs? Because they can be more efficient that FIRs (re: fewer multiplications), which allows higher sampling rates compared to FIRs IIRs consist of feedforward coefficients $b(k)$ (like an FIR filter) and feedback coefficients $a(k)$, which couple past outputs to the current on You can’t calculate the IIR coefficients from the impulse response (not even $b(k)$!) Laplace Transforms Useful tool for solving differential equations Take Laplace transform, do algebra to get output function’s form in Laplace domain, and then inverse Laplace transform Defined as $F(s) = \int_{0}^{\infty} f(t) e^{-st} dt$ where s is a complex number: $s = \sigma+\omega i$, which can be though of as a complex frequency The Laplace transform allows you to convert the time series of the inputs and outputs into the complex frequency space ($X(s)$ and $Y(s)$ respectively) The transformation which enables you to go from $X(s)$ to $Y(s)$ is called the transfer function $H(s)$ ie. $Y(s) = H(s) X(s)$ In principle, given a transfer function, one can Laplace transform x(t), apply H(s), then inverse Laplace to get y(t) In practice, H(s) have enough information to give you the frequency response as well as if the filter is stable Stability TL;DR look at the poles of H(s) What are poles? They are complex numbers $s_{0}$ where the transfer functions blows up They are called poles because if you look at the magnitude of H(s) around $s_{0}$, it looks like a flag pole H(s) can have multiple poles Let $s_{0} = \sigma_{0}+\omega_{0} i$ be a pole of H(s). If $\sigma_{0}$ is positive, the system is unstable. If $\sigma_{0}=0$, then the system is conditionally unstable (re. it will oscillate back and forth) There are also zeros of H(s), which correspond to the numerator beinng zero People typically denote “x” for poles and “o” for zeros on the complex plane Z-Transform This is a discretized version of the Laplace transform that can be used to construct IIR filters Continuous Laplace equations get applied to differential equations, while the z-transform gets applied to difference equations It’s defined as $H(z) = \Sigma_{n=-\infty}^{\infty} h(n) z^{-n}$ If we write z in polar form: $H(z) = \Sigma_{n=-\infty}^{\infty} h(n) r^{-n} exp(-i\omega n)$ From the r term, we immediately see that the stability condition is that $|z| \leq 1$ The ranges of frequencies are also restricted by the sampling frequency (ie. $\frac{\omega_{s}}{2} \geq \omega \geq \frac{-\omega_{s}}{2}$) Hence, if all of the poles of a transfer function are located inside the unit circle, things are stable Adding a delay in your signal processing pipeline looks like $y(n) = x(n-1)$ Plugging that into the z-transform, changing variables to $k=n-1$ and rearranging, you get that $Y(z) = z^{-1} X(z)$ or more verbosely, having k delay units multiplies the Z transform by $z^{-k}$ Alternatively, a transfer funtion of $z^{-k}$ can be thought of as a delay of $kt_{s}$ seconds where $t_{s} = \frac{1}{f_{s}}$ Applications to IIR We can write any M-th order IIR filter with N feedforward states and M feedback stages in Direct Form I We can write the output Z-transform as: $Y(z) = X(z) \Sigma_{k=0}^{N} b(k) z^{-k} + Y(z) \Sigma_{k=1}^{M} a(k) z^{-k}$ We can then define the transfer function as $H(z) = \frac{Y(z)}{X(z)} = \frac{\Sigma_{k=0}^{N} b(k) z^{-k}}{1-\Sigma_{k=1}^{M} a(k) z^{-k}}$ The order of the transfer function is the larger of M or N You can make the substitution $z = e^{i\omega}$ (remember, stability means that the radius is 1!) to get the frequency response This can be though of as taking the intersection of the unit radius cylinder with $H(\omega)$, cutting the curve at the $\pi$ branch cut and then unwrapping everything Various Notations for IIR Filters For simplicity, assume that M and N are equal (same number of feedforward and feedbackwards stages) Since everything is linear and time invariant, we have some freedom to shuffle around the multiplications and additions Transposition Theorem To transpose a filter: Reverse the direction of all signal-flow arrows Convert all adders to signal nodes Convert all signal nodes to adders Swap the x(n) input and y(n) output labels Why would you do this? Say you’re transposing Direct Form I. The transpose first calculates the zeros, then the poles. The poles can lead to a blowup in intermediate calculations, so postponing this calculation until the end helps prevent truncation Finite Precision Effects in IIRs Since we can’t use infinite precision arithmetic, we need to used a finite precision representation of all of our computation The three major problems from this discretization are: Coefficient quantization: Lacking the necessary precision on coefficients causes the zeros and poles to wobble around, which can result in stability issues and a distorted frequency response Overflow errors: Intermediate calculations can go outside the dynamic range of the precision specified. This can lead to overflow oscillations Roundoff errors: To mitigate the overflow issues, you can round off intermediate calculations to some specified precision You can normally treat roundoff errors as some random process The roundoff noise becomes highly correlated with the signal (called limit cycles). This can cause infinite oscillations Chaining Filters Together If you cascade one filter after another, then overall transfer function of the system is the product of the two The associated coefficients are also the convolution of the two original filters If you run two filters in parallel with each other (ie. same inputs) then sum their outputs, the overall transfer function is the sum of the two individual transfer functions The coefficients are simply the sum of the two original filters This idea allows you to break down a higher order filter into smaller order ones, which are easier to work with One big problem is a combinatorial explosion of equivalent circuits (assuming infinite precision) To demonstrate this, say you want to break down a 2M order filter into M distinct 2nd-order filters there are M! ways of arranging these filters, and there are M! ways of pairing coefficients of the filters If you’re dealing with 2nd order filters, there is a simple (but not-optimal) algorithm to make these choices: First, factor a high-order IIR filter into the form $H(z) = \frac{\Sigma_{i}^{N} (z-z_{i})}{\Sigma_{j}^{N} (z-p_{i})^}$ where N denotes the number of poles and zeros After than, find a pole or pole pair closest to the unit circle, then find the zero which is closest to those poles. These poles and zeros form one of your 2nd order filters. Remove them from the pool, repeat this procedure until nothing is left Scaling Gain of IIR Filters For some applications, we want to be able to limit the gain of the data values in the IIR filter This is accomplished by modifying the DC gain of the filter The DC gain is defined as $\frac{\Sigma b(k)}{1-\Sigma a(k)}$ or in words: the sume of the feedforward coefficients divided by 1 minus the sum of the feedback coefficients If we want to have a DC gain of 1, we simply divide all of the coefficients by this initial DC gain Neither the frequency response nor the phase response does not change Differentiators, Integrators, and Matched Filters Differentiators Imagine a continuous analog sine wave. It’s derivative is a cosine times the angular frequency The transfer function of the ideal differentiator is $H(\omega) = i \omega$ Can easily derive by taking the Fourier transform of derivative of signal The naive way of implementing this is via a tapped-delay line structure A First difference differentiator simply subtracts adjacent sames from each other The First difference method is very sensitive to noise and high frequency signals This follows from the frequency response of the filter: $|H(\omega)| = 2 |\sin (\frac{\omega}{2})|$ A Central-difference differentiator takes the average of two consecutive first difference differentiators This works out to $\frac{x(n)-x(n-2)}{2}$ The central difference method is more robust against noise, but there is the linear regime is smaller than that of the finite difference method Follows from $|H(\omega)| = |\sin (\omega)|$ Both of these methods have a very narrow operating regime. Which means that you can only use these if the input has a low spectral content compared to the sampling rate The group delay of both of these filters goes like $\frac{D}{2}$ where D is the number of taps involved This follows from the anti-symmetry in the coefficients Weather the delay is an integer or not is important for specific applications In practice, you shove the frequency response through Park-McClellan to get the best performing filter Oddly enough, you get a more accurate differentiation when N is even compared to when N is odd You can also window the input sequence to reduce the ripples in frequency space Narrowband Differentiators Richard Hamming proposed a family of low-noise Lanczos differentiating filters, which have 2M+1 coefficients You generate the coefficients via $h(k) = \frac{-3k}{M(M+1)(2M+1)}$ where k ranges from -M to M This has a smaller linear operating frequency range compared to the naive way You can implement this filter with half the multiplications as the number of taps You can also multiply through by the common denominator in order to do integer multiplication instead of floating point Wideband Differentiators Suppose that we want to be able to differentiate a large number of frequencies up to some cutoff $\omega_{c}$, after which we have perfect attenuation This is impossible to do without an infinite number of taps All approximations will have some ripples at frequencies above the cutoff Taking the inverse transform of the ideal frequency response yields $h_{gen}(k) = \frac{\omega_{c}\cos(\omega_{c}k)}{\pi k} - \frac{\omega_{c}\sin(\omega_{c}k)}{\pi k^{2}}$ This is sufficient, but it has negatives, which are annoying to work with We can rewrite the coefficients as $h_{gen}(k) = \frac{\omega_{c} \cos(\omega_{c} \lceil k-M\rceil)}{\pi \lceil k-M\rceil}-\frac{\omega_{c} \sin(\omega_{c} \lceil k-M\rceil)}{\pi \lceil k-M\rceil^{2}}$ $M = \frac{N-1}{2}$ and $0 \leq k \leq N-1$ Integrators An ideal analog integrator has a frequency response proportional to $\frac{1}{\omega}$ Discretization prevents us from doing this, but there are several schemes which approximate this Rectangular Rule: Simply keep a running sum of all of the inputs $y_{Re}(n) = x(n) + y_{Re}(n-1)$ $H_{Re}(\omega) = \frac{1}{1-exp(-i\omega)}$ Trapezoidal Rule: Average the two most recent incoming signals $y_{Tr}(n) = 0.5 x(n) + 0.5 x(n-1) + y_{Tr}(n-1)$ $H_{Tr}(\omega) = \frac{0.5+0.5exp(-i\omega)}{1-exp(-i \omega)}$ Simpson’s Rule: Weight the three most recent incoming signals as though you fit a quadratic to them $y_{Si}(n) = \frac{x(n)+4x(n-1)+x(n-2)}{3} + y_{Si}(n-2)$ $H_{Si}(\omega) = \frac{1+4 exp(-i\omega)+exp(-2\omega i)}{3(1-exp(-2i \omega))}$ All of the above have z-domain transfer function poles at z=1, which corresponds to a cyclic frequency of 0 This means that the register widths need to be large enough to prevent overflow from corrupting an integrator The above shows the absolute error of each of the integrators Tick’s rule is some weird scheme which trades off drastically larger errors at high frequencies for more precise low frequency approximation Matched Filters Suppose that our input consists of a signal component and a noise component We want the filter to be able to output the signal to noise ratio (SNR), which can then be passed to some thresholding function to determine if a signal of interest was detected The design goal of the filter is to maximize the SNR For an analog signal, the optimum frequency response is $H(f) = S(f)^{\dag}e^{-2\pi f T i}$ where $\dag$ refers to complex conjugation and T is the time duration of the signal Physically, this is the complex conjugate of the frequency response times a linear phase shift which is proportional to T Can think of this as convolution: you have some reference template which you drag across the input signal. The maximum SNR occurs when the dot product is minimized Physically, the phase of the signal constructively interferes when the input matches the template, but the noise destructively interferes We can discretize under the following assumptions: The signal of interest is some N sample sequence $S(m)$ is the N-point DFT of the signal of interest m is some discrete frequency index $0 \leq m \leq N-1$ $x_{in}$ is sampled at a rate of $f_{s}$ The time of the N-length sample sequence is N-1 sample periods Hence we get $H(m) = S^{\dag}(m) exp(-\frac{2\pi m (N-1)}{N}i)$ the associated coefficients are $h(k) = x_{s}(N-k-1)$ where $0 \leq k \leq N-1$ A matched filter can be thought of as a digital correlator If the samplings are large, it may be more efficient to do multiplication in frequency domain than convolution in time domain Quadrature Signals What if, instead of a real signal, you have a complex signal? This complex signal is called a quadrature signal, and takes the form $e^{i 2\pi f t} = cos(2\pi f t) + i sin (2\pi f t)$ You can generate an “imaginary” signal via multiplication by i: delay the phase of the input by $\frac{\pi}{2}$ In order to generate a quadrature, you could have two function generators which generate cosine and sine respectively. As long as you keep the phase delay between the two pinned to $\frac{\pi}{2}$, then the two wires together are a quadrature Hence, you always need two wires to represent a quadrature Can imagine two classes of quadrature signals in the complex plane. One rotates clockwise, while the other goes counter-clockwise This gives rise to the notion of “negative frequencies”. They are simply signals which have the opposite polarization as the standard signals If you are representing a real signal in quadrature, the negative and positive frequencies are evenly symmetric around the DC point Similarly, purely imaginary signals have odd symmetry around the DC point Down Conversion Imagine that you have some real-valued time sequence x(n). This has a corresponding frequency representation which is symmetric around f=0 and is joined at $\frac{f_{s}}{2}$ If we multiply x(n) by $e^{-i 2\pi f_{c} n t_{s}}$, where $t_{s} = \frac{1}{f_{s}}$ and $f_{c}$ is some frequency shift less than $f_{s}$, then we have effectively shifted the spectra over by $f_{c}$ This multiplication involves splitting up the signal and multiplying each branch by cosine and negative sine respectively An analogous construction holds for quadrature signals One big sticking point of this type of down conversion is that there needs to be a very strict 90 phase delay between the wires of the quadrature signal If you digitize before you do the mixing, this helps lock the phase difference Complex Resonators You can think of the above IIR filter as a way of generating a frequency between 0 and $f_{s}$ Hz If you let x(n) be a unity impulse (re 1, then 0 forever after), this circuit will oscillate at a rate specified by $\omega_{r}$ $\omega_{r}$ is a normalized angle between 0 and $2\pi$ So if you want $\frac{f_{s}}{4}$ as the output frequency, set $\omega_{r} = \frac{\pi}{2}$ and take the real component of the output Hilbert Transform The Hilbert transformation let’s you convert real signals to complex signals Can think of this as multiplication by i. In physical terms, this is a 90 degree phase delay. This can be combined with the original signal to produce a quadrature signal So we need to design a circuit which, given a purely real input signal, outputs the same signal, but delays by 90 degrees What is the frequency response of the shifted output? $X_{ht}(\omega) = H(\omega) X_{r}(\omega)$ $H(\omega) = \pm i $ where -i is over positive frequencies and +i is over negative frequencies Suppose that you can perform the Hilbert transform on $x_{r}$. The Fourier transform of the quadrature signal is $X_{c}(\omega) = X_{r}(\omega)+ i X_{i}(\omega)$ This is zero in the negative frequencies (ie. a one sided signal) Impulse Response for Hilbert Transform We know what the response function in frequency space looks like: $H(\omega) = \pm i $ Taking the inverse Fourier transform yields: $h(t) = \frac{1-cos(\frac{\omega_{s} t}{2})}{\pi t}$ where $\omega_{s} = 2\pi f_{s}$ This can be digitized to $h(n) = \frac{f_{s}}{\pi n}(1-cos(\pi n))$ For $n=0$, let $h(n)=0$ FIR Implementation Naively, we could simply just use the h(n) coefficients and convolve them with the incoming time sequence. Little problem: there are an infinite number of coefficients, so we need to truncate the series If the truncation as an odd number of taps and the coefficients are antisymmetric (as in the case of the Hilbert transform), $|H(0)| = 0$ and $|H(\frac{\omega_{s}}{2})| = 0$. If there are an even number of taps, only the DC component of H needs to be 0 Look at the above 15 tap Hilbert transform above The magnitude response drops off at 0. This effectively turns the Hilbert transform into a bandpass filter There is some ripple in the passband. This comes from the finite number of taps in the coefficients. You can reduce the ripple via windowing, but this will reduce the bandwidth of the filter Due to the aforementioned drop off, it’s hard to compute the Hilbert Transform of low-frequency signals The phase response is linear (with a slight discontinuity at 0) To offset this phase response, we delay the real component of the quadrature signal by the group delay (which remember, is $G = \frac{D}{2}$ where D is the number of delay elements) Given a set of coefficients, remember to reverse the order in order to do the convolution. For most FIR filters, the coefficients are symmetric, so this doesn’t matter, but Hilbert filters aren’t Sample Rate Conversion Say that you have two hardware processors operating at two different sample rates. You need to exchange data between the two. How do you do this? One way this could be done is to do use a DAC and ADC setup to convert to an analog signal, then resample at the new rate. This limits the dynamic range of the signal and is not normally utilized What you do depends if you are increase the sample rate or decrease the sample rate As an aside, people in the DSP community use a lot of similar sounding words to describe various concepts of sampling rate conversion Use the above to make the transition between terms Decimation Decimation is how you go about downsampling. It consists of low pass filtering a signal, and then downsampling Downsampling by a factor of M is done by only keeping every Mth sample, then discarding the rest. This reduces the sampling rate by a factor of M Care needs to be taken to ensure that the new sampling rate still obeys the Nyquist criterion (ie. the sampling rate is twice the bandwidth of the signal), otherwise aliasing can occur The lowpass portion is essential when you drop below the Nyquist limit (ie. it prevents the upper frequencies of your signal from aliasing with the lower frequencies) The low pass portion of decimation can most easily be done with a FIR filter (due to the linear phase response) Two-Stage Decimation It may be more computationally efficient and stable to break a decimation block into multiple sub-blocks Ie. The number of taps that you is can be significantly lower than a single stage If you have two stage decimation (so the downsampling factor is $M = M_{1}M_{2}$), the optimal value for $M_{1}$ becomes: $M_{opt} = 2M \frac{1-\sqrt{\frac{MF}{2-F}}}{2-F(M+1)}$ $F = \frac{f_{stop}-B}{f_{stop}}$ This is the ratio of the single-stage low pass filter’s transition region width to that of the filter’s stopband frequency Downsampling Gotchas Downsampling is not a time invariant process. So you aren’t allowed to reorder time invariant operations with it While downsampling doesn’t cause time-domain signal amplitude loss, there is a drop in magnitude by a factor of M in the frequency domain Interpolation Interpolation can be thought of as upsampling, followed by a low pass filter Upsampling by a factor of L can be thought of as inserting L-1 samples between each of the original signal’s samples One easy way of upsampling is via zero stuffing: make all of the new intermediate values 0! You then need to lowpass filter this interpolated signal to smooth the signal Another potential way of interpolating is just to repeat the values of the original sequence until the next one The problem with this is that the frequency domain gets multiplied by a sinc function. This would require an egregious number of taps in the lowpass filter to smooth out Interpolation Gotchas An artifact of upsampling is that, in frequency space, you get phantom copies of the original spectra. An ideal lowpass would eliminate these copies completely, but due to the attenuation of real filters, there is always some residue Due to the zeros, the interpolation reduces the time series magnitude by a factor of L. You need a corresponding DC offset of L to maintain a unity transformation Conversely, the frequency magnitude will increase by a factor of L Polyphase Filters Naively, if you want a sampling rate that this not an integer, then you would need to chain decimations and interpolations of different integers Doing this by zero padding and discarding samples is excessively wasteful Polyphase filters helps eliminate this unnecessary computation Polyphase Interpolation and Decimation Interpolation is easier to understand. Suppose that you have a 12-tap low pass filter and you want to upsample by a factor of 4 For x(0), there are only three meaningful multiplications; the rest get zeroed out For x(1) (re: the next sample over), there are still only three meaningful multiplications; the difference is that the FIR coefficients you are using are all different This repeats for x(2) and x(3) If you perform all of these branches in parallel, and then multiplex the output based on the current sample, then you have effectively eliminated all of the zeroing The implementation shown in the image above is not really any more efficient than the zeroing scheme (3 out of the 4 branches are just getting wasted) Instead of multiplexing the output, you can multiplex the coefficients instead. This way, only one branch is active at a given time Some gotchas: Try and make sure that the initial FIR filter has an integer multiple of your interpolation factor for ease of implementation Polyphase interpolation still has drops in magnitude by a factor of L Can fix by scaling the coefficients by L, or multiplying the output by L Decimation is very similar In fact, you use the same FIR coefficients as interpolation! Now, instead of multiplexing the output, you multiplex your input down each branch sequentially, then sum the outputs of all the branches The z-transformation of a up/downsampled signal is simply $W(z) = X(z^{L})$ for upsampling and $W(z) = X(z^{\frac{1}{L}})$ for downsampling Because upsampling and downsampling are not time invariant, you can’t swap them around with other filters arbitrarily However, if the filters have some special form, then it’s alright Specifically, you can swap a transfer function of the form $H(z^{L})$ witth an upsampler and downsampler as per the noble identities Cascaded Integrator Comb (CIC) Filters CIC filters are computationally efficient implementations of narrowband lowpass filters Very useful for anti-aliasing filtering in decimation and interpolation A key advantage of these filters is that they require only additions and subtractions The above image shows D point moving averagers A running sum is an averager that doesn’t do the one multiplication Fig a shows a standard non-recursive averager Fig b shows an equivalent recursive setup: instead of summing all of the delayed elements, you delay the signal by D cycles, then subtract that from the input, multiply by 1/D, and then add the last sample Fig c is equivalent to Fig b, but without the multiplication Fig d is equivalent to Fig c, but uses more compact notation for the delay elements and annotates each section The governing difference equation of Fig. d is $y(n) = x(n)-x(n-D)+y(n-1)$ The feedforward section (re $x(n)-x(n-D)$) is the comb section, the delay unit is called differential delay, and the feedback section is called the integrator The associated transfer function is then $H(z) = \frac{1-z^{D}}{1-z^{-1}}$ The corresponding frequency response is $H(f) = exp(-i\pi f(D-1)) \frac{\sin(\pi f D)}{\sin(\pi f)}$ The above shows the characteristics of a D=5 CIC filter The magnitude of the frequency response is like a squared sinc function The phase response is piecewise linear, with jump discontinuities at the crossing points of the frequency magnitudes There is a pole on the edge of the unit circle, which ostensibly should lead to conditional stability. However, since there is no coefficient quantization error, this gives guaranteed stability to CIC filters The DC gain of the filter is D (use L’Hopital’s rule on frequency response at f=0) Signal Averaging Can think of the mean, standard deviation, and variance of a set of N samples as the DC signal, the AC signal, and the power of the AC signal All of these involve calculating the mean of a signal. How do you do that? Coherent Averaging An alternative name for this is time-synchronous averaging. This means that the start of each sample set is taken at the same phase of the target signal Take for instance, a sine wave with some noise on it. If the sine wave samples all start at the same phase, they will converge to the true sample, while the noise, which is decoherent, will go to zero This sort of averaging across a set of samples yields a SNR increase of $\sqrt{N}$, where $SNR = \frac{R}{\sigma}$ Incoherent Averaging This means that the start of each set of samples is not taken at phase (think of a oscilloscope without a proper trigger getting set) This is useless in the time-domain, but pretty useful in the frequency domain One particular application is a spectrum analyzer You want to figure out the spectra of a signal, but you don’t want a large tap DFT. Since addition is cheaper than multiplication, you can just take the average of multiple smaller DFTS You get a similar increase in SNR of the power measurements. This gain of SNR is called integration gain Averaging Multiple DFTs Suppose that you have k DFTs, each with N taps Simply average across all of the k bins You get a reduction of the noise by a factor of k In practice, there is an overlap between adjacent DFTs (so points in the first DFT also show up in the second DFT etc.) How much overlap is too much? Keep it less than 50% and you don’t see any noticable reduction in SNR Phase Averaging This is harder than time-domain signal averaging and frequency domain magnitude averaging due to the branch cut Treat the phase angles as the arguments of two complex numbers, add the numbers, then extract out the phase Averaging Implementations You can use the standard moving averager discussed in previous sections, in both the recursive and non-recursive setup The problem with these are that N (the number of sample sets) must be an integer The second is that abrupt amplitude changes in the input don’t result in abrupt changes in the output (lowpass filter behavior and all that…) One popular averaging scheme is exponential averaging Define the recursive lowpass filter: $y(n) = \alpha x(n) + (1-\alpha) y(n-1)$ where $0 < \alpha < 1$ There are three big benefits of this type of averaging: You can tune the frequency response of the filter by changing $\alpha$ It requires less computation than the standard moving averager You only need one delay element The above shows the response of the exponential averager when the input is a step function If $\alpha$ is 1, then the input is not attenuated at all and no averaging takes place This makes it so that the averager responds immediately to changes in the input As $\alpha$ decreases, the response time of the filter also increases. The expoential weighting of previous inputs averages the signal and reduces noise A trick people use is to have a large $\alpha$ value initially so that there is an immediate output, then taper $\alpha$ off to increase SNR Given a certain noise reduction factor R, we can determine what $\alpha$ we need via: $\alpha = \frac{2}{R+1}$ The frequency response has a pole at $z=1-\alpha$, so this is unconditionally stable The phase of this filter is nonlinear. We normally don’t care about this since averagers are more concerned with the DC response Digital Signal Processing Tricks Frequency Translation by $\frac{f_{s}}{2}$ Multiply x(n) by the sequence $(-1)^{n}$, where both sequences are run at $f_{s}$ sampling rate This is efficient since you don’t actually need to do multiplication; just flip the sign bit This can be shown by doing the convolution theorem use the form $(-1)^{n} = exp(i\pi n)$ End effect is to shift the top and bottom halves of the spectra down and up (respectively) by a factor of $\frac{f_{s}}{2}$ Frequency Translation by $-\frac{f_{s}}{4}$ Define the sequences $\cos(\frac{\pi n }{2}) = 1, 0, -1, 0, …$ $\sin(\frac{\pi n}{2} ) = 0, -1, 0, 1, …$ You can multiply these sequences with any x(n) to shift the frequency spectrum by $\frac{f_{s}{4}$ (sign depends on initial phase of cosine and sine sequences) Don’t acutally need multiplication, just sign flips and zeroing out The only difference between the outputs is the 90 degree phase shift Can use this to generate quadrature signal High Speed Vector Magnitude Approximation $(\alpha Max + \beta Min$) Suppose that $I^{2}$ and $Q^{2}$ (representing the in-phase and quadrature signal powers) are readily available. How do you quickly calculate the square root of the sum of the two? A decent approximation is $|V| \approx \alpha Max + \beta Min$, where max and min are the larger and smaller of I and Q Can quickly compare magnitudes by dropping sign bit and comparing values If $\alpha$ and $\beta$ are both 1, this is just another way of writing the triangle inequality $\alpha$ and $\beta$ are typically chosen to be powers of two for easy multiplication, but not a hard restriction Need to experimentally determine optimal hyperparameters The error is dependent on the relative magnitude (re: phase) between I and Q Conditional Faster Complex Multiplication In the very specific case where the hardware allows 3 addition operations to be faster than 1 multilplication, you can do complex multiplication faster by computing some intermediate values Define $k_{1} = a(c+d)$, $k_{2} = d(a+b)$, and $k_{3} = c(b-a)$ where a+ib and c+id are your two complex numbers then $R = k_{1} - k_{2}$ and $I = k_{1} + k_{3}$ This takes 3 multiplications and 5 additions instead of 4 multiplications and 2 additions Working with Real-Valued Data for FFTs The general FFT algorithm accepts complex numbers. Most data is real. You can trim the fat of the FFT to take advantage of this One way is to do the FFT of two seperate real sequences at the same time Suppose that you have two N point input sequences. Treat one sequence (a) as the real component, and the other (b) as the imaginary component Do the N point complex FFT. You can extract the FFTs of the original sequences via: $X_{a}(m) = \frac{X^{*}(N-m)+X(m)}{2}$ $X_{b}(m) = \frac{i(X^{*}(N-m)-X(m))}{2}$ There is a slight gotcha at m=0, where X(N) doesn’t exist. You just need to remember that $X(N) = X(0)$ due to the periodicity A similar trick works to calculate a 2N real FFT with a N sample complex FFT Break the sequence into even and odd samples, shoving the even into the real and odd into the imaginary Do the complex FFT, then define the follow 4 quantities $X_{r}^{+}(m) = \frac{X_{r}(m)+X_{r}(N-m)}{2}$ $X_{r}^{-}(m) = \frac{X_{r}(m)-X_{r}(N-m)}{2}$ $X_{i}^{+}(m) = \frac{X_{i}(m)+X_{i}(N-m)}{2}$ $X_{i}^{-}(m) = \frac{X_{i}(m)-X_{i}(N-m)}{2}$ The full FFT is then $X_{real}(m) = X_{r}^{+}(m)+\cos(\frac{\pi m}{N}) X_{i}^{+}(m) - \sin(\frac{\pi m}{N}) X_{r}^{-}(m)$ $X_{imag}(m) = X_{i}^{+}(m)+\sin(\frac{\pi m}{N}) X_{i}^{+}(m) - \cos(\frac{\pi m}{N}) X_{r}^{-}(m)$ This is a sequence of length N-1, but the other have can be gotten via complex cojugation due to the realness of the input Computing Inverse FFT using Forward FFT The first method conjugates signal, does the FFT, then conjugates again The second swaps the real and imaginary components, does the FFT, then swaps the output You can easily prove the above two are equivalent to the inverse FFT by chugging through the algebra Folded FIRs If an FIR filter has symmetric coefficients, then you can reduce the number of multiplications by roughly half Simply add up the signals which share the same coefficient (so if you have 5 taps, signal 0 adds to signal 4, signal 1 adds to signal 3 etc.), then multiply by the coefficient Reducing ADC Quantization Noise Option 1: simply increase the sampling rate The noise scales like the inverse sampling frequency Option 2: dither. This means add some noise to the analog signal prior to discretizing This seems like a stupid idea, but it’s apparently not Take for instance the case of clipping. If an analog signal goes outside the range of the ADC, then values will get pegged to the maximum/minimum outputs of the ADC This clipping creates artificial spectral harmonics Another way of saying this is that clipping causes the quantization noise to become coherent with the signal If you dither an analog signal, then the quantization noise becomes less coherent, which means that these stray spectra harmonics are dampened Your noise floor increases as as result, but the SNR also increases due to the elimination of the stray harmonics You can use dithering for low amplitude analog signals highly periodic analog signals slowly varying analog signals ADC Tests Coherence Sampling Imagine that you feed a single sine wave of unity amplitude into an ADC. This will span the full range ADC and let you characterize the response of the ADC via an FFT there is spectral leakage that needs to be mitigated for this to work. the main trick is to have the ADC sampling rate be coherent with the signal frequency. This means that $f_{sample} = \frac{mf_{sig}}{n}$ The choice of m and n matter. People have found that if you make m a odd prime number help minimize the redundant output codes of the ADC and make the quantization more random (desirable in this case) The FFT amplitude gets boosted by a factor of $10 * log_{10}(\frac{64}{2})$, so the actual SNR of the signal is smaller than the FFT would suggest SINAD (Signal to Noise and Distortion) Also called SNDR To get a better indicator of the SNR, calcualate the SINAD Do the N Point FFT of the sequence and discard the negative frequencies Compute the total signal spectral power (re: summ the square of the bins with signal in them) sum the total noise-only spectral power Calculate $SINAD = 10* log_{10}(\frac{P_{sig}}{P_{noise}})$ Dynamic Range Testing To test the dynamic range, feed two equal amplitude analog sinewaves into an ADC Increase the amplitudes until you get some spurious spectral component which rises above the noise floor Define the SFDR as the dB difference between the high level signal spectral magnitude sample and the spurious spectral signal Missing Codes Occasionally, ADCs will be incapable of producing a specific binary word (re: a code). so for an 8-bit ADC, it could happen that the number 33 is never output Just histogram all of the output values. Any missing codes will show up as 0 bins Time Domain Convolution as FFT In some cases, doing a time domain convolution in frequency space is faster than the FIR implementation Simply break up the signal into chunks (with a power of 2 size per chunk, zero padding as needed), FFT each of them, multiply in frequency space, then inverse FFT Sharpening FIR Filters Suppose that you want to improve the stopband attentuation of a FIR filter, but you can’t change its coefficients You could double up the filter (re: make a copy, then run it after the first), but this also increases the passband ripple Alternatively, you can do the above to keep the stopband attenuation, but reduce the ripple power G denotes the desired gain of the signal The Delay block is a consists ofo $\frac{N-1}{2}$ delays where N is the number of samples The x2 and x3 can be implemented via left shift by 1, and left shift by 1 added to original signal I don’t know why this works since the book doesn’t give a derivation Goertzel Algorithm Say you want to detect if a single frequency is in a signal. You could do an FFT, then check if the bin you want is above some threshold. But that requires computing the full FFT Goertzel’s algorithm lets you only compute a subset of the bins of a DFT As pictured above, Goertzel has two feedback coefficients, and one feedforward one The time domain difference equations are: $w(n) = 2\cos(\frac{2\pi m}{N}) w(n-1) -w(n-2) + x(n)$ $y(n) = w(n)-exp(-i \frac{2 \pi m}{N}) w(n-1)$ The transfer function of the IIR filter is as follows: $\frac{Y(z)}{X(z)} = \frac{1-exp(-i \frac{2\pi m }{N})z^{-1}}{1- 2 \cos (\frac{2\pi m }{N}) z^{-1}+z^{-2}}$ There is both a pole and zero at $exp(-i\frac{2\pi m}{N})$ (plus a conjugate pole which doesn’t cancel out) The frequency response is to become sharply peaked around the m frequency bin The one surviving pole is on the unit circle, which might cause stability problems Luckily, the algorithm resets w(n-1) and w(n-2) for after each block, preventing the filter from running for too long and hitting instability The advantages of Goertzel’s algorithm over the FFT: N doesn’t need to be a power of 2 The resonant frequency can by anything between 0 and $f_{s}$ Hz The number of coefficients is much smaller You don’t need to wait for the entire block of data to come in before starting the FFT You can implement M copies of Goertzel for multi-frequency detection and it will still be more efficient than an FFT if $M < log_{2}(N) $ When doing Goertzel, you only retain every Nth sample. This causes the frequency response to become a sinc function Sliding DFT Goertzel’s algorithm compute a single complex DFT bin for every N input samples The Sliding window algorithm lets you compute a DFT bin on a sample-by-sample basis You can evaluate the DFT for a single m bin: $X^{m}(q+1) = \Sigma_{n=0}^{N-1} x(n+q+1) exp(-i \frac{2\pi mn}{N})$ The moniker “sliding” comes from the fact that the x(n) values used to compute $X^{q}$ gets shifted over by 1 every time step You can do some index juggling to tease out a recurence relationship for $X^{m}$ $X^{m}(q+1) = exp(i \frac{2\pi m}{N})(X^{m}(q)+x(q+N)-x(q)$ This version appears to look ahead at the future. That’s impossible. You can index on n instead of q to get: $X^{m}(q+1) = exp(i \frac{2\pi m}{N})(X^{m}(n-1)+x(n)-x(n-N)$ This have a pole on the unit circle, which can result in stability issues This can be mitigated by multiplying by $-r^{N}$ instead of -1 where r<1 Zoom DFT Instead of taking the FFT over a large number of samples, you get get a better frequency resolution around some region centered at $f_{c}$ with the above scheme There is an inflection point where using a Zoom FFT is more advantageous compared to just using a higher resolution FFT Fast Arctangent Sometimes, you need to compute the phase of a quadrature signal let x = I+iQ. Can use above to do a dirty arctan caclutation depends on the Octant you’re in (start at real axis, increase index in a counter-clockwise manner) Frequency Demodulation Tying into the above, we can calculate the instantaneous frequency of a signal by: converting it to a quadrature signal ( $ x=X+iY = A \exp(i\theta) = A \exp(2\pi f t i)$) calculating the phase ($ \theta = tan^{-1}(\frac{Y}{X})$) differentiating the phase w.r.t. time ($\Delta \theta$) use $f = f_{s}\frac{\Delta\theta}{2 \pi}$ to map the change in phase to a frequency If you want to avoid doing the costly arctan calculation, just take the derivative of $\arctan(\frac{Y}{X})$ with respect to time Doing that yields the top schematic If you know the magnitude of your signal is constant (re purely frequency modulation), then the denominator can be dropped and you can use a basic central difference scheme for differentiation This is the bottom schematic DC Removal Simplest way is compute the average of a block of length N, then subtract this from the each original sample This is not real time The IIR Filters above do real time DC removal They are all equivalent Their response functions are $H(z) = \frac{1-z^{-1}}{1-\alpha z^{-1}}$ There is a zero at z=1, and a pole at $\alpha$ The zero means that there is infinite attenuation at DC As $\alpha$ approaches 1, then the attenuation around the DC point gets sharper The phase experiences a jump discontinuity at DC and tappers off away from DC Efficient Polynomial Evaluation Horner’s Rule A standard kth order polynomial looks like: $f_{k}(x) = \Sigma_{0}^{k} c_{k}x^{k} You can massage the above to be written in terms of FMA (fused multiplication and addition), which is dedicated hardware circuit on some devices For specificity, look at the cubic. We have $c_{3}x^{3}+c_{2}x^{2}+c_{1}x^{1}+c_{0} = x(x(c_{3}x + c_{2}) + c_{1}) + c_{0}$ The above is easily generalizable, and each group can be performed with FMA With fixed point, large left and right shifts can cause overflow A similar trick applies to this case, but swap multiplication for shifting Estrin’s Method Horner’s Rule is inherently serial. You can parallelize it via writing the polynomials as nested expressions with a base rule of $ax^{q}+b$ This maps well onto parallel MAC operations found in most modern CPUs High Order FIR Filters How do you create a high-order FIR filter when your DSP software using Parks-McClellan refuses to converge In brief, you make a lower order FIR filter whose response is a crude approximation of your desired response, then you zero pad the frequency responese and recalculate the new coefficients In detail: Define a target number of taps to be M*N Set the number of taps in the lower order FIR filter (hence forth the prototype filter) to be N After calculating the coefficients for the prototype, perform an N-point DFT to get frequency response of the prototype Insert M-1 zero-valued samples just before the $\frac{N}{2}$ sample of the frequency response to get the new frequency response Calculate the inverse DSP to get the new interpolated FIR coefficients Take the real part to deal with finite precision issues Multiply the coefficients by M to compensate for the 1/M amplitude loss Test the interpolated coefficients (re: compare the FFT of the interpolated coefficients to the FFT of a massively zero padded version of the interpolated coefficients) Automatic Gain Control (AGC) The simple linear AGC circuit compares subtracts off the signals output power from some desired reference power, then feeds back this difference with some scaling factor $\alpha$ There is some lag time before the signal returns to the constant power output. This lag/attack time is input dependent This lack of control over the attack time leads nicely to a logarithmic AGC design The logarithmic AGC gives you more control over the attack time of the filter The low pass filter helps suppress rapid changes in the input The logarithmic design also has more dynamic range in comparision to the linear one Quadrature Osscillator The difference equation is given as $y_{r}(n) = y_{r}(n-1)\cos\theta - y_{i}(n-1) \sin \theta$ $y_{i}(n) = y_{r}(n-1)\sin\theta + q_{i}(n-1) \cos\theta$ This is just multiplying the previous input by $e^{i\theta}$ One big problem is that the input can grow/decay out of control. introduce AGC into the mix Faster Variance The unbiased variance can be re-written as $\sigma_{ub} = \frac{\Sigma_{n=1}^{N} (x(n)^{2}) - Nx_{avg}^{2}}{N-1}$ This introduces two additional multiplications at the expense of N fewer additions In the standard formula, we need to retain the all of the N samples in order to computer the average. With this set up, we don’t need to do that Signal Transition Detection Suppose that you need to detect a transition in an input which spans many samples. A standard digital differentiator would have the memory to accomudate that. We can use the so called “time-domain slope filtering”, which is a FIR filter with the following coefficients: $C_{k} = -\frac{12k-6(N-1)}{N^{3}-N}$ where N is any integer The affect of the above is to create a linear ramp in coefficient value This approach is more robust against high noise levels in comparision to a standard digital filter

Notes following this MIT course and is based on the xv6 toy operating system. Operating System Interfaces Processes and Memory Forking Exec Operating System Organization Abstracting Physical Resources CPU Modes Processes xv6 startup Security Page Tables Kernel Address space Xv6 Implementation Physical Memory Allocation Xv6 Physical Memory Implementation Process Address Space Sbrk (Xv6) Exec (Xv6) Real World Paging Traps and System Calls RISC-V traps User Space Traps Syscalls Kernel Space Traps Page Fault Exceptions Interrupts ad Device Drivers Console Input Console Output Timer Interrupts Locking Races Spinlocks Using Locks Deadlocks Spinlocks Scheduling Context Switching Scheduler Sleep and Wakeup Semaphores Lost Wake-Up Problem Pipes Process Locks Real World Scheduling File System Buffer Cache Layer Logging Layer Block Allocator Inode Layer Directory Layer Path Layer File Descriptor Layer Operating System Interfaces Operating systems serve a number of purposes Abstracts away lower level hardware, so that programs don’t have to care about what disk you are writing to Allows multiple programs to run “at the same time” Provide well defined interface for programs to interact with each other For the interface, we want a simple interface that is easy to implement, but allows high sophisticated operations to happen UNIX philosophy of composing many simple programs together helps keep this problem trackable Can divide operating system into kernel space and user space kernel space is the program which provides services to other programs. It’s given hardware level privileges which user programs don’t have access to Each running program under the kernel’s manager (called a process) gets it’s own memory which houses the instructions (ie. the CPU assembly instructions), the data (variables created when running the program) and the stack (organizes program procedure calls) if process needs a kernel service, it invokes a system call, which briefly transfers ownership to the kernel, which executes the syscall, and then returns ownership to the user program Called user space and kernel space Typically, people use a shell to allow user input. The shell resides in user space, so you have flexibility in what shell you run Processes and Memory For each running process, there is some user-space memory (instructions, data, stack) and kernel-space memory (which the user’s can’t access) One of the kernel-space things is the PID (process identifier) Xv6 has it’s own custom syscalls, which are somewhat stripped down versions of the UNIX ones Unless otherwise stated, all of these calls return 0 for no error and -1 if there’s an error Call Sig Desc int fork() Create a process, return child’s PID. Makes an exact copy of user-space. returns in both parent and child process. Returns child’s PID in parent process, and 0 in child process int exit(int status) Terminate the current process; Releases all held resources. Exit status reported to wait(). No return. int wait(int *status) Wait for a child to exit; exit status in *status (out parameter); returns child PID. If caller has no children, immediately returns -1. If parent don’t care about child exit status, pass 0 address (ie. (int *) 0) to wait int kill(int pid) Terminate process PID. Returns 0, or -1 for error. int getpid() Return the current process’s PID. int sleep(int n) Pause for n clock ticks. int exec(char *file, char *argv[]) Load a file and execute it with arguments; only returns if error. This replaces the current process’s memory with that of the file. For Xv6, the ELF format is used. char *sbrk(int n) Grow process’s memory by n bytes. Returns start of new memory. int open(char *file, int flags) Open a file; flags indicate read/write; returns an fd (file descriptor). int write(int fd, char *buf, int n) Write n bytes from buf to file descriptor fd; returns n. int read(int fd, char *buf, int n) Read n bytes into buf; returns number read; or 0 if end of file. int close(int fd) Release open file fd. int dup(int fd) Return a new file descriptor referring to the same file as fd. int pipe(int p[]) Create a pipe, put read/write file descriptors in p[0] and p[1]. int chdir(char *dir) Change the current directory. int mkdir(char *dir) Create a new directory. int mknod(char *file, int, int) Create a device file. int fstat(int fd, struct stat *st) Place info about an open file into *st. int stat(char *file, struct stat *st) Place info about a named file into *st. int link(char *file1, char *file2) Create another name (file2) for the file file1. int unlink(char *file) Remove a file. Forking As a small case study, look at the following C snippet int pid = fork(); if(pid > 0){ printf("parent: child=%d\n", pid); pid = wait((int *) 0); printf("child %d is done\n", pid); } else if(pid == 0){ printf("child: exiting\n"); exit(0); } else { printf("fork error\n"); } After the fork(), the pid let’s you determine what the parent and the child do NOTE: While after forking, the memory and register content of both processes is the same, they located in different locations which can evolve independently from each other Don’t expect the output to be the same. Whichever process gets to printf first will output first (parent and child output could also be interleaved) Exec Another case study of exec char *argv[3]; argv[0] = "echo"; argv[1] = "hello"; argv[2] = 0; exec("/bin/echo", argv); printf("exec error\n"); If echo is functioning properly, then the last printf will never get executed Argv[0] is typically reserved for the program name Operating System Organization An operating system has three requirements multiplexing: multiple processes should be able to access the same hardware resources (time-sharing) isolation: each process should not be subject to the bugs of another process (ie. if one process goes under, the others keep running) interaction: processes should be able to interact with each other (eg. like UNIX pipes) xv86 runs on a multi-core RISC-V processor, and is written in “LP64 C” ,which means longs and pointers are 64 bits, but integers are 32 bits. Abstracting Physical Resources In theory, you don’t need an operating system. Just provide a library to interact with the hardware resources, then allow programs to link against this library real-time systems and embedded devices do this for speed reason The problem with this approach is that it requires each program to trust each other (ie. yield control of the CPU time to achieve multiplexing, called coperative time-sharing) and not have any bugs. Both of these are ludicrously naive UNIX uses a file system to abstract the hardware so that programs don’t have to fully trust each other This hardware abstraction also allows transparent switching between hardware CPUs (save registers when you switch, load them back in when you come back) CPU Modes Three modes of CPU machine mode: instructions have full privilege privilege mode: only the kernel can run in this mode. If a user program tries to execute one of these instructions, it’s gets terminated by the kernel user mode: most restricted mode of cpu Need to use a special function to transition to privilege mode Operating systems need to make a choice about what mode it will run in UNIX makes the choice of putting everything in privilege mode (monolithic kernel) Another choice is to have most of the code in user space though (called a microkernel), since if you mess up in kernel mode, everything crashes Processes A process is stored in a struct at kernel/proc.h inn xv86. This info lets the process act like it has it’s own CPU, and act like it has it’s own memory Each process has a thread of control, which holds the state needed to execute the process (p->thread) Each process has two stacks: a user space stack(p->stack), and a kernel space stack (p->kstack) in RISC-V, ecall is used to enter kernel space, and sret is used to return to user space p->state indicates wheather the process is allocaated, ready to run, currently running, waiting for I/O, or exiting p->pagetable holds the process table info, formatted for RISC-V The pagetable is responsible for translating virtual addresses (ie. the addreses a CPU instruction manipulates) to a physical address (an address that the CPU chip sends to main memory) xv6 startup CPU first reads and execute the bootloader from ROM, which loads the rest of the kernel The RISC-V starts with paging hardware disabled The bootloader puts the kernel at 0x80000000 (place here b/c I/O devices are in the lower range), and the starts executing at _entry at _entry, there is a start function which does things only machine mode: acts like it returns from supervisor mode (even though you are running in machine mode) sets the return address to that of main disables the virtual address space delegates all interrupts and exceptions to supervisor mode and then programs the clock chip to generate interrupts in the main function it creates the first user process the first syscall called is exec, replaces the current process with the shell the args to exec get put in a0 and a1, and puts the syscall number in a7 syscall number matches the entries in the syscalls array, a table of function pointers syscall retrieves the syscall number from a7 The first syscall is SYS_exec once sys_exec returns, it’s return value is put in p->trapframe->a0 This causes the original userspace call to return that value Security The kernel assumes that user code is malicious (ie. it will try it’s best to crash the system) The kernel space code is assumed to be safe\ Page Tables Xv6 uses page tables to implement isolation. Maps virtual addresses (address that RISC-V instruction manipulates) to physical addresses (address that CPU sends to main memory) the top page of the address houses the trampoline page and the trapframe page, each of 4096 bytes the trampoline page contains the code to transition in and out of the kernel the trapframe page houses the process’s user registers To elaborate on the mapping: only the bottom 39 bits out of the 64 bits are used for virtual address in RISC-V, there are $2^{27}$ page table entries (PTE). each PTE contains a 44-bit physical page number (PPN) These flags are used to state now that page can be used PTE_V indicates wheather a PTE exists PTE_R indicates if you can read the page PTE_W controls if you can write to a page PTE_U controls if user mode instructions can access the page The top 27 bits of the 39 bits are used to index into the page table to find a PTE You can then construct a 56-bit physical address whose top 44 bits come from the PPN in the PTE, and whose bottom 12 bits are from the virtual address RISC-V uses a hierarchy of pagetables, since that uses less data (for example, three sections of 9 bits are used to index 3 levels of page tables) Kernel Address space Each process gets one page table. The kernel gets it’s own page table to describe it’s address space The kernel uses direct mapping (ie. virtual addresses are equal to physical addresses). The cases where direct mapping isn’t used by the kernel are The trampoline page: it’s mapped to the top of the address space. The user page tables also have this same mapping. This allows you to easily transfers registers between layers The kernel stack pages: each process gets it’s own kernel stack, along with an unmapped guard page whose PTE is invalid. This helps guard against stack overflows Xv6 Implementation The central datastructure is pagetable_t, which is a pointer to a RISC-V root page-table page This can either be a kernel page table, or a preprocess page table Early on in booting, the kernel’s page table gets constructed prior to paging being enabled The root page-table pge gets allocated, then gets populated with the translations the kernel needs, then the kernel stack for each process gets allocated (leaving space for the guard pages) and the appropriate flags are set Physical Memory Allocation The kernel is responsible for allocating and freeing physical memory at runtime for page tables, user memory, kernel stacks, and pipe buffers xv6 uses physical memory from the end of ther kernel to PHYSTOP for run-time allocation Allocation consists of removing a page from the linked list, while freeing is adding a page to the list Xv6 Physical Memory Implementation the main data structure of the allocator is a free list of physical memory pages each entry in the list is a struct run. The struct stores each run structure in the free page itself, since nothing else is stored there The struct has a spin lock on it, which means you need to acquire and release the lock (don’t worry about this for now…) The free list is initialized to hold every page between the end of the kernel and and PHYSTOP Normal operating systems can figure out how much to allocate via some hardware configurations. Xv6 just assumes 128 Mb of RAM a PTE can only refer to a physical address that is aligned to the 4096-byte boundary, so PGROUNDUP macro does this alignment The allocator treats addresses in two different ways As integers in order to perform arithmetic As pointers to read and write memory The aforementioned reasons is why there are so many C type casts Xv6 sets all freed memory to 1 via kfree (to try and catch dangling references) kfree repends the page to the free list (so the list grows head first) Process Address Space Each process gets it’s own page table. When processes get switched, it also changes page tables When a process asks for more memory, xv6 uses kalloc to allocate physical pages, and then populates the page table with the correct PTEs. The OS also sets the correct permissions on the page This schema has some nice advantages Different process page tables translate user addresses to different pages of physical memory. This allows each process to have a private user memory Each process sees its memroy as having contiguous virtual addresses starting at 0, which physical memory is non-contiguous The kernel maps a page with trampoline code at the top of user address space, which means that a single page of physical memory shows up in add address spaces (hence, allows you to share between processes) Referring to above figure, we note that the user stack is a single page, followed by a guard page, which acts as a way to detect a stack overflow. Xv6 just panics, but a real OS would try to allocate more memory for the user stack Sbrk (Xv6) Sbrk is used by a process to grow or shrink memory. Depending on the sign of allocation, sbrk either allocates or deallocates a page at a time as needed (so if you allocate 1 byte, 1 whole page will be allocated) Exec (Xv6) Exec creates the user part of an address space If does this by reading a file stored in the file system in the form of an ELF file It first reads the first 4 bytes to check of the magic number matches ‘0x7F’, ‘E’, ‘L’, ‘F’. If it matches, it assumes that the file is formatted correctly A new empty page table is allocated, then each ELF segment is allocated with uvmalloc then loaded with loadseg walkaddr is used to find the physical addresses of the allocated memory exec then allocates and initializes one page for the user stack, along with a guard page to protect against overflow and too many arguments If any program segment is invalid, the program image is freed, and exec returns -1. Exec is risky, since addresses in the ELF file can refer to the kernel. So user programs can in theory gain kernel privileges Real World Paging real world operating systems deviate from Xv6 drastically Real OS’s don’t crash if a user stack overflows. They actually try to save a page=fault exception Xv6 assumes that physical RAM is always at 0x8000000. The real world can’t make this assumption. Real OS’s need to be able to map arbitrary hardware physical memory layouts into predictable kernel virtual address layouts On machines with large amounts of memory, you might want to use “super pages” (ie. make pages 4 Mb instead of 4kB) to avoid excess page-table manipulation the lack of malloc prevents the kernel from using data structures with require dynamic allocation Real kernels would allow allocation of variable size blocks instead of a hard-coded 4096 bytes Traps and System Calls “Traps” refer to any event which causes the CPU to set aside ordinary execution of instructions Could be a system call Could be an exception (ie. an instruction has done something illegal) Could be a device interrupt: when a device signals that it needs attention Traps often should be transparent to other processes (ie. other code shouldn’t know that they happened) The overall structure of trap handling is a trap forces control to transfer to the kernel the kernel saves registers and other state so that execution can be resumed later the kernel executes the appropriate handler code the original code resumes where it left off “Handler” refers to code which that processes a trap The first handler instructions written in assembly are called a “vector” RISC-V traps RISC-V CPUs have a set of control registers that the kernel writes to to tell the CPU how to handle traps the kernel can read from to find out about how a trap occurred The most important registers are: stvec: the address of it’s trap handler. RISC-V jumps to this address to handle the trap sepc: When a trap occurs, RISC-V save the program counter here the sret instruction copies sepc to pc. The kernel can write to sepc to change where sret goes scause: RISC-V puts a number here that describes the reason for the trap sscratch: the trap handler code uses this to help overwrite user registers before saving them (ie. a scratch pad) sstatus: the SIE bit in sstatus controls weather device interrupts are enabled. If the kernel clears SIE, RISC-V will defer device interrupts until the kernel sets SIE. The SPP bit indicates whether a trap came from user mode or supervisor mode and controls to what mode rset returns You can’t read or write these registers in user mode. A similar set of registers exist for traps handled in machine mode (xV86 only uses these for timer interrupts) RISC-V does the following sequence: if the trap is a device interrupt, ssstatus SIE is cleared. Don’t do any of the following Disable interrupts by clearing the SIE bit in sstatus Copy the pc to sepc Save the current mode in the SPP bit in status set scause to reflect the trap’s cause set the mode to supervisor copy stvec to the pc Start executing in supervisor mode In theory, you could skip some of these steps. In practice, you want to maintain process isolation (ie. prevent user code from running kernel level instructions) Note that the CPU doesn’t do the following operations. It’s the kernel’s job to do these tasks. The CPU does this to allow greater flexibility for the OS (for instance, some OS’s omit the page table switch in some situations for performance reasons) Doesn’t switch to the kernel page table Doesn’t switch to a stack in the kernel Doesn’t save any registers other than the pc User Space Traps In RISC-V, the hardware doesn’t switch page tables on a forced trap. This means that the trap handler address must have a valid mapping in the user page table and in the kernel page table each process automatically gets assigned a trampoline page, which is located at the top of the virtual address space above any user code Cruicially, it is located at the same address in user and kernel space, which makes the transition between the two easy The csrw instruction copies a register over to sscratch the TRAPFRAME is saved just below the TRAMPOLINE virtual address Once all of the registers are copied, you switch to the kernel page table, then figures out what trap occurred, then dispatches the correct subroutine Once the interrupt routine finishes, the RISC-V reloads the old program counter, then switches to the original user page table Syscalls the arguments for exec are placed in a0 and a1, and places the system call number in a7. This system call number indexes into a syscall array, which is an array of function pointers After the associated syscall returns, the return value is recorded in p->trapframe->a0 C calling convention on RISC-V places return values in a0 Conventionally, negative numbers are errors, positive number are success Kernel Space Traps Upon entering the kernel, all 32 registers are pushed to the stack for later restoration Importantly, they are saved on the interrupted kernel thread. So if the trap causes a switch to a different thread, then these interrupted threads will be preserved If an interrupt was due to a timer interrupt, and a process’s kernel thread is running, then the current thread yields to allow other threads to run Once a trap is done, it returns to the original code There in theory can be a race condition when the ‘stvec’ register is still set ot uservec instead of kernelvec. If an interrupt occurred in the interim, that would be bad RISC-V disables interrupts when it starts to take a trap, and they don’t get enabled until after it sets stvec Page Fault Exceptions For xv86, if an exception occurs in user space, the kernel kills the faulting process. If an exception happens in the kernel, the kernel panics. Real kernels can do more interesting things For instance, copy-on-write (COW) fork. The basic idea is that when you create a fork, you need to copy over the memory space of the parent to the child. Xv6 allocate a separate chunk of physical memory for the child. it would be nice if you could use the same physical memory, but a naive implementation would cause the two processes to step over each other the CPU can raise a page-fault exception if a virtual address has no mapping in the page table, or PTE_V is cleared, or a process doesn’t have the correct permissions Initially, parent and child shares all physical pages, but both are read-only. If a page-fault from an attempted write occurred, then a new page of physical memory is allocated for that process, but now with write enabled. Subsequent writes are unimpeded Many use cases of fork() doesn’t allocate new memory, so this save a bunch of page allocations Lazy allocation is another use case. If an application demands memory, the kernel keeps track of the increase in size, but doesn’t allocate any physical memory until a page fault occurs If an application asks for more memory than it needs, you save on allocating all those pages Even if the application massively grows the address space, you can spread out the cost of allocating that space over time This does incur the overhead of more paging requests, but you can allocate a small number of pages on page fault There is demand paging. Naively (like xv6 does), one loads in all of the application text and data into memory before starting the application. You can instead create a user page table with all PTEs marked invalid. Once the program starts and encounters a page for the first time, the kernel allocates that page Programs might need more memory than there is in RAM. the kernel splits the data between RAM and disk, and marks all the disk pages as invalid. If a disk page gets hit, then a new page is allocated in RAM, the page from the disk is copied over, and the PTE is pointed to RAM (this process is called ‘paging in’. ‘paging out’ goes the other way) If there is no free physical RAM, then the kernel must free a physical page (or evict it) from RAM to the disk, and mark the page as invalid again. This is expensive eviction is infrequent if an application only uses a subset of their memory pages which can all fit in RAM Most computers have little to no free memory (think multiplexing user on a cloud computer). When free physical memory is scarce, allocation is expensive Interrupts ad Device Drivers A “driver” refers to code in an operating system that manages a particular device. This means setting the device up, handling interrupts, and mediate processes waiting for I/O from the device Importantly this need to be done in an thread safe concurrent fashion Typically, you can split a driver into two halves The top half runs in a process’s kernel thread. This is called via system calls such as read and write that want to perform device I/O The bottom half executes at interrupt time. On interrupt, the bottom figures out what was done, wakes up a waiting process if needed, and then tells the hardware to start work on any queued operation Console Input The console driver takes inputs from the UART serial-port hardware and stores it in a line buffer. User processes then use the read syscall to fetch these lines from the console The UART appears to software via a set of memory-mapped control registers. That is, certain addresses map to the hardware instead of RAM Console Output The device driver maintains some output buffer so that writing processes don’t need to wait for the UART to finish sending You need to wait if the buffer is full though This is a general pattern: You decouple the device activity from the process activity via buffering and interrupts The console driver can process input even when no process is waiting to read it Similarly, processes can send output without having to wait for the device This decoupling technique is called I/O concurrency. It’s useful for slow devices (like UART) or if you need to do immediate action Timer Interrupts Xv6 needs timer interrupts to maintain the current time, and switch amoung compute-bound processes These are just roll-over counters with a fixed increment frequency and well defined thresholds for rollover (stimecmp in Xv6) Timer interrupts only increment time for one CPU (prevents time from passing faster if therer are multiple CPUs) The kernel code may be moved from one CPU to another without warning Locking Most kernels allow the execution of multiple activities Multiple CPUs can execute independently Multiple CPUs share physical RAM, which allows data structures to be shared between them A CPU, even on a single processor system, can switch between multiple threads of execution A device interrupt handler can interrupt at any time You need to manage this parallelism in order to ensure correctness of program execution (concurrency control techniques) Locks provide mutual exclusion, ensuring that only one CPU at a time can hold the lock Upside: Easy to understand Downside: Serializes concurrent operations Races Imagine that you have your bog standard linked list struct element { int data; struct element *next; }; struct element *list = 0; void push(int data) { struct element *l; l = malloc(sizeof *l); l->data = data; l->next = list; list = l; } What happens if two CPUs execute push on the same linked list? Only the one that occurs latter in time will be successfully realized. The first one will get lost. This is an example of a race Races occur when simultaneous memory access occurs with at least one being a write The outcome of races can depend on The machine code generated by the compiler The timing of the two CPUs involved How the two CPU memory operations are ordered by the memory system Locks are one way of avoiding races. Locks ensure mutual exclusion, so that only one CPU can execute the offending lines of code at a time These sections are often called critical sections Locks protect data by maintaining invariants on the data In the critical sections, these invariants might be violated, but at the end, the invariants should hold again Look at linked lists: The invariant is that the head of the list points to the first element, and the next field points to the next element push temporarily violates this invariant when it updates the pointers. Locks ensure that is section is run by only one CPU Locks decrease performance by serializing concurrent code If multiple processes need a lock at the same time, then they are in contention with each other A big part of kernel design is to avoid locks when possible via clever data structures and algorithms Locks also need to be placed correctly. If they are two wide, they can be serializing code which doesn’t need to be Linked lists again: if you locked the entirety of push, then even malloc would be serialized when it doesn’t need to be Spinlocks One implementation of locks is a spinlock. You have one bit of data which represents if a lock has been acquired or not. Naively, one might implement this like void acquire(struct spinlock *lk) // does not work! { for(;;) { if(lk->locked == 0) { lk->locked = 1; break; } } } The problem occurs at the if statement. If two CPUs both reach line 25, they both might see that lk->locked is 0, which means both of them acquire the lock. This violates mutual exclusion To get around this, most multi-core CPUs have some sort of atomic read and swap instruction which does the following in an uninterruptible way (using ampswap r, a in RISC-V as an example): reads the value at memory address a writes the contents of register r to address a puts the read value into r in Xv6, acquiring and releasing a lock basically boils down to using this instruction To acquire a lock, a CPU continually spams the ampswap instruction by writing a 1 to lk->locked If you get a 1 in the return value, then you know that some other CPU has the lock If you get a 0, then this means that no other CPU was holding the lock previously. Hence, the current CPU now has exclusive ownership of the lock and can move on with life Release works in a similar way, but spams 0 instead Using Locks It can be difficult to figure out what data and invariants each lock should protect Some principle to help guide you: Any time a variable can be written by one CPU at the same time that another CPU can read or write it, a lock should be used to keep the two operations from overlapping If an invariant involves multiple memory locations, typically all of these locations need to be protected by a single lock to ensure the invariance is maintained The above are rules for when locks are necessary, but they say nothing about when they are unnecessary You can divide locking into “coarse-locking” and “fine-locking” to handle different speed requirements For instance, in Xv6, there is a single free list with a single global lock (coarse grain) For the file system, each file has it’s own lock, so that processes manipulating different files don’t need to wait for each Deadlocks Locks need to be acquired and released in the same order, otherwise, you get deadlocks For instance, if one CPU acquires B then A, and the other acquires A then B, then both CPUs will be waiting for the other to finish and will stall For Xv6, there are well defined lock ordering which must be maintained daf In general, it’s difficult to maintain global deadlocking avoiding order For instance, module M1 might call module M2, but the lock order requires the lock in M2 to be acquired before M1 This sets constraints on how fine grained that you can set your locks One way that you might be able to avoid deadlocks is with re-entrant locks (or recursive locks) This means that if a process attempts to acquire a lock it already has, then the kernel just allows it instead of panicking struct spinlock lock; int data = 0; // protected by lock f() { acquire(&lock); if(data == 0){ call_once(); h(); data = 1; } release(&lock); } g() { aquire(&lock); if(data == 0){ call_once(); data = 1; } release(&lock); } h() { ... } Take a look at the above code, you would expect that call_once() would be called once, either by f or by g, but not by both If re-entrant locks are allowed and h happens to call g, call_once will be called twice Without recursive locks, you would get a deadlock. Either bug is bad, but a deadlock is easier to track down since the kernel panics on deadlock, but silently fails with reentrant locks If you wanted to use reentrant locks, then you would need to keep track of the depth of the lock You can get deadlocks with interrupts as well For instance, suppose that a CPU acquires a lock for sys_sleep, and the CPU gets a timer interrupt. The interrupt will try to get the lock, won’t be able to get it, and then will wait for it to be released. sys_sleep won’t every release until the timer interrupt is handled To avoid this, if a spinlock is used by an interrupt handler, a CPU must naver hold the lock while interrupts are enabled Xv6 is more conservative: if a CPU acquires any lock, then interrupts are disabled on that CPU You can also uncover deadlocks when the compiler does out of order execution Mentally, you can think of assembly code are running in the order that the code is written is This is not the case: compiler might reorder your code, omit extraneous memory stores and loads, run independent sections of code in parallel etc. That last one isn’t perfect. You can get code that works sequentially, but fails miserably when done out of order You can tell a compiler not to re-order via memory barriers (__sync_syncrhonize() in C library) Spinlocks Sometimes, a lock needs to be held for a long time Say in the file system, a file is held for tens of milliseconds while read and write operations happen Spinlocks are not good for this situation since a lot of CPU time is wasted spinning the lock Another drawback is that a process cannot yield the CPU while retaining the spick lock ie. Another process wants to use the CPU while your waiting for the slow disk if you yield with a spinlock, then you run into the possibility that another thread of the CPU tries to acquire the spinlock. You then have both threads waiting for each other to finish (deadlock!) The solution to this is sleep-locks This consists of a locked field which is protected by a spinlock, and acquiresleep’s call which does an atomic sleep on the CPU, yielding as needed Since sleep locks keep interrupts enabled, the cannot be used in interrupt handlers Scheduling Operating systems typically run more processes than the computer has CPUs, so you need some way to time-share the CPUs amongst the processes Ideally, this sharing should be transparent to the user A common way to do this is to give the illusion that each process has it’s own CPU by multiplexing processes (ie. the work of a process gets shared between CPUs) There are a bunch of challenges that you need to handle to get this working First, you need to be able to save and restore CPU registers (this needs to be done in assembly) How do you force switches in a way that is transparent to user processes? Xv6 uses the standard technique in which hardware timer interrupts drive context switches all CPUs switch amoung the same set of processes, so you need a locking plan to avoid races A process’s memory and other resources must be freed when the process exits, but it can’t do all of this itself for various reasons One of these is that it can’t free ints own kernel stack while still using it Each CPU must remember which process it is executing so that syscalls affect the correct process kernel Finally, sleep and wakeup allow a process to give up the CPU and wait to be woken up by another process or interrupt Context Switching The above shows how you switch from one user process to another a trap from user space to the old process’s kernel thread a context switch to the current CPU’s scheduler thread a context switch to a new process’s kernel thread a trap return to the new user-level process Xv6 has dedicated threads which execute the scheduler because it’s not safe for the scheduler to execute on any process’s kernel stack Another CPU might wake the process up and run it, and it would be bad to have the same stack on two different CPUs Hence, each CPU has it’s own scheduler thread to cope with situations in which more than one CPU is running a process Upon a context switch, the CPU scheduler thread will save the old process’s registers, stack and program counter, and then load the new process’s registers, stack and program counter The function swtch does this save and restore. It only focuses on saving and restoring as set of 32 RISC-V registers, collectively called “contexts” swtch only saves the callee registers. It does not save the program counter; it save the ra register (ie. the return address from which swtch was called) After that, swtch restores registers from the new context, which holds the register value of a previous swtch call. Upon return, the instrutions point to the restored RA (ie. the instruction which called swtch). It also returns on the new thread’s stack Scheduler Basic ideas is that if a process wants to give up the CPU, it must acquire its own process lock release any other locks it is holding update its own process state Disables interrupts then call the scheduler Each CPU runs a scheduler in it’s own dedicated thread The CPU checks that the above sequence was run, then context switches to the scheduler thread The scheduler then dispatches to any waiting processes NOTE: the process lock is held across context switches This must be done since context switching does not preserve invariants. The lock can only be released if the scheduler clears the cpu’s process the new process’s kernel is completely up and running The scheduler in Xv6 is very simple It loops through the process table and looks for a process that is RUNNABLE, sets that process to RUNNING, then context switches and waits for it to yields. It then continues on Sleep and Wakeup What if you need two threads to communicate with each other? Xv6 uses a mechanism called sleep and wakeup Sleep allows a kernel thread to wait for a specific event A thread can call wakeup to indicate that threads waiting on a specific event should resume Semaphores As an example of the utility of Sleep and Wakeup, let’s try and implement a higher-level syncrhonization mechanism called a semaphore. Semaphores maintain a spinlock and a reference count “Producers” can increment this count when ready “Consumers” busy wait for this count to be non-zeroo, then decrement the counter and return A naive implementation might look like 100 struct semaphore { 101 struct spinlock lock; 102 int count; 103 }; 104 105 void 106 V(struct semaphore *s) 107 { 108 acquire(&s->lock); 109 s->count += 1; 110 release(&s->lock); 111 } 112 113 void 114 P(struct semaphore *s) 115 { 116 while(s->count == 0) 117 ; 118 acquire(&s->lock); 119 s->count -= 1; 120 release(&s->lock); 121 } This implementation is expensive, since the consumer is busy waiting for potentially long periods of time To avoid the busy wait would require the consumer to yield the CPU, and only resume after function V increments the counter Suppose that there exists some functions sleep(chan) and wakeup(chan) sleep(chan) puts the calling process to sleep and releases the CPU for other work wakeup(chan) wakes all processes that are in sleep calls with the same channel number. If no processes are waiting on chan, wakeup does nothing The above allows us to change the semaphore implementation as follows 200 void 201 V(struct semaphore *s) 202 { 203 acquire(&s->lock); 204 s->count += 1; 205 wakeup(s); 206 release(&s->lock); 207 } 208 209 void 210 P(struct semaphore *s) 211 { 212 while(s->count == 0) 213 sleep(s); 214 acquire(&s->lock); 215 s->count -= 1; 216 release(&s->lock); 217 } The consumer is no longer busy waiting (it gives up the CPU instead) Lost Wake-Up Problem The above sleep and wakeup implementation runs into the “Lost Wake-Up” Problem Suppose that P finds that s->count ==0 on line 212. While P is between lines 212 and 213, V runs on another CPU, which changes s->count to be nonzero and calls wakeup, which finds no processes and does nothing. P executes lines 213, but missed the previous wake-up, so it continuously sleeps. You need to wait for V to send another wakeup, which might never happen This happens because the invariant that P sleeps only when s->count==0 is violated by V running at just the wrong moment Naively fixing this by moving the acquire call creates a deadlock THe actual fix is to change the interfact of sleep to be sleep(chan, cond_lock), where cond_lock is the semaphore’s spinlock. This allows the lock to be released only after the process has be properly put to sleep This forces a concurrent V to wait until P as finished sleeping Hence, the proper implementation is 400 void 401 V(struct semaphore *s) 402 { 403 acquire(&s->lock); 404 s->count += 1; 405 wakeup(s); 406 release(&s->lock); 407 } 408 409 void 410 P(struct semaphore *s) 411 { 412 acquire(&s->lock); 413 while(s->count == 0) 414 sleep(s, &s->lock); 415 s->count -= 1; 416 release(&s->lock); 417 } Pipes A pipe struct consists of a spinlock a data buffer nread and nwrite keep track of the number of bytes read/written The buffer is a circular buffer (ie. the next byte after $buf[PIPESIZE-1]$ is $buf[0]$). nread and nwrite do not wrap around This convention allows ups to distinguish between a full buffer (nwrite=nread+PIPESIZE) from an empty buffer(nwrite==nread) To index the buffer, you need to use an index of nread % PIPESIZE instead of nread The pipe works by busy-waiting for the lock to free If writing, the pipe writes to the buffer until it’s out of data, or the buffer is full. If the buffer gets full, then wakeup is called to alert any sleeping readers that there is data to be read The write pipe then sleeps and waits for a wakeup call indicating that the pipe has space sleep releases the pipe’s lock during this process If reading, then after acquiring the lock, the pipe reads from the buffer until the buffer is empty or the desired bytes have been recovered If the buffer is empty, then you read all the available bytes, then send a wakeup to all sleeping write pipes Process Locks p->lock must be held while reading or writing any of the following proc fields p-> state, p->chan, p->killed p->xstate, x->pid This is because these fields can be used by other processes On a higher level, p->lock also prevents race conditions (along with p->state) when allocating new processes Conceals a process from view while it is being created or destroyed Prevent a parent’s wait from collecting a process that has set its state to ZOMBIE but has not yet yielded the CPU It prevents another CPU’s scheduler from deciding to run a yielding process after it sets its state to RUNNABLE but before it finishes swtch Ensures only one CPU’s scheduler decides to run a RUNNABLE process prevents a timer interrupt from causing a process to yield while it is in swtch Prevents wakeup from overlooking a process that is calling sleep but hasn’t yielded the CPU yet Prevent the victim process of kill from exiting makes kill check and write p->state in an atomic fashion Real World Scheduling Xv6 uses round-robin shedueling, where each process is run in turn. Real schedulers might implement some priority scheme so that some processes are given preferential treatment. This is a hard thing to do (things like priority inversion and convoys can happen) sleep and wakeup is a simple synchronization method, but it’s not the only one: for example, the Linux kernel uses an explicit process queue Linear scans of the processes during wakeup is inefficient. There are a variety of data structures which could be more efficient. In the context of threading libraries, this structure is referred to as a condition variable sleep is mapped onto wait wakeup is mapped onto signal wakeup currently wakes up all sleeping processes on a particular channel. This causes all the processes to wake up and rush to grab the lock This sort of behavior is called “a thundering herd” A more sane way might be to have a way to wakeup only one process on a channel, or broadcast a wake up to all sleeping processes on a channel File System A file system’s purpose is to organize and store data; it allows one to share data amoung users and applications, as well as persistence so that data is still available after a reboot There are lots of challenges a file system must handle You need on-disk data structures to represent the tree of named directories, the identities of the blocks that hold a file’s content, and to record the free areas of the disk Support crash recovery: If a crash occurs, the file system should still work/ be self-consistent Different processes must be able to operate on the file system at the same time, so the file system must manage invariants The file system needs an in-memory cache of hot blocks, since disk access is orders of magnitude sower than accessing memory Xv6 uses a 7 layer file system. From bottom up the first layer is the disk the buffer cache layer caches disk blocks and synchronizes access to them, ensuring that only one kernel process at a time can modify the data stored in any particular block Logging layer allows higher layer to wrap updates as atomic transactions. This ensures that the blocks are updated atomically in the face of crashes The inode layer provides individual files Each inode as a unique i-number and some blocks which hold the file data The directory layer implements each directory as a special kind of inode whose content is a sequence of directory entries, which each contains a file’s name and i-number The pathname layer provides hierarchical path names and resolves then with recursive lookup The file description layer creates a top-level file system interface for user-level programs Disk hardware traditionally presents data on the disk in 512-bit wide blocks called sectors. These sectors are sequentially layed out For instance, sector 0 is the first 512 bits, sector 1 is the next 512 bits etc. Xv6 holds recently read sectors in memory This data sometimes is out of sync with the disk (perhaps it hasn’t read from the disk, or maybe the cached data has been updated by the software, but not the hardware) Xv6 splits the file system up into several blocks: Block 0 is the boot sector. The filesystem doesn’t touch this Block 1 is the superblock. It contains meta data bout the file system (file system size in blocks, the number of data blocks, the number of inodes, and the number of blocks in the log) Block 2 holds the rest of the file system Buffer Cache Layer Has two jobs Synchronize access to disk blocks to ensure that only one copy of a block is in memory and that only one kernel thread at a time uses that copy Cache popular blocks so that they don’t need to be re-read from the slow disk The interface of Xv6 uses bread obtains a buf struct which contains a copy of the block which can be read or modified in memory bwrite pushes this modified buffer to the appropriate block on the disk brelse frees the buffer once it’s done The buffer cache uses a doubly linked list of buffers This gets initialized with NBUF empty buf structs Each buffer has a valid field, which states wheather the buffer contains a copy of the block, and a disk field, which indicates that the buffer content has been handed to the disk If you reach the end of the list and you ask for a block that is not in the cache, the buffer cache boots out the least recently used buffer for the new block bget scans the buffer list for a buffer with the given device and sector numbers. If there are none in the cache, then the least recently used block in the cache is forcibly cleared, made invalid so that bread is forced to go to disk, then presented as the return from bread brelse releases the sleep lock on a buffer, then moves that buffer to the front of the linked list. This orders the list from front to back by how recently the buffer was released Logging Layer Suppose that you only perform a subset of writes before the computer crashes The crash may leave an allocated but unreference content block The crash might also leave an inode with a reference to a content block If xv6 supports multiple users, then the latter would allow the old file’s owner to read and write block in the new file, owned by a different user Xv6 performs logging by keeping a seperate log which describes all the write a syscall wants to perform. Once all of them are queued up, the syscall writes a commit record to the disk indicating that the log contains a complete operation. The syscall then copies the writes over to the disk. After the writes complete, the log on the disk gets erased Imagine a crash happens If a log is marked as containing a complete operation, then the recovery code copies the writes to where they belong If a log is not marked as complete, the recovery code ignores the log In either case, the log gets erased by the recovery code Let’s think about why this works: If the crash occurs before the operation commits, then the log in the disk will not be marked as complete. The recovery code will ignore it and the state of the disk will be as though the operation didn’t happen If the crash occurs after the operation commits, then recovery will replay all of the operation’s writes It might even repeat some that the operation has started prior to the crash The above makes disk writing atomic: either all the writes appear on the disk, or none of them appear On Xv6, the log resides in the superblock It contains a header block, followed by a sequence of updated block copies The header contains an array of sector numbers (one for each logged block) and the count of log blocks The count is either zero (indicating no blocks) or non-zero (indicatin the log contains a completed transaction with the indicated number of logged blocks) Xv6 write the header block when a transaction commits. The counter gets set to 0 after copying the logged blocks to the file system So a crash midway through a transation will result in a count of zero A crash after a commit will have a non-zero count The log can accumulate logs across many syscalls. To avoid splitting a syscall across transactions, the logging system only commits when no file-system syscalls are happening This batching of transactions is known as group commit The log has a fixed amount of disk space. If a syscall asks for too many blocks, then you need to split these blocks up so that everything fits (for instance, write() could asks for many blocks) Block Allocator We have a free pool of disk blocks which need to be allocated Xv6 block allocator uses a bitmap on disk. A zero bit in the mask indicates the block is free. Otherwise, it is in use Mkfs sets the bits corresponding to the boot sector, superblock, log blocks, inode blocks and bitmask blocks The block allocate provides two functions balloc allocates a new disk block. This loops over all possible blocks (from 0 to sb.size), looking for a free block. Once found, the bitmap is updated and the block is returned There is an optimization where the loop is split into two parts. The outer loop reads seach block of bitmap bits. The inner loop checks all Bits-Per-Block (BPB) in a single bitmap block A race might occur if two processes try to allocate a block at the same time, but the buffer cache prevents this by only letting one process use any one bitmap block at a time bfree finds the right bitmap block and clears the right bit Inode Layer Inode has two meanings The on-disk data structure containing a file’s size and list of data block numbers An in-memory inode, which contains a copy of the on-disk inode as well as extra information needed within the kernel The on-disk inodes are in a contiguous area called the inode blocks. Each inode is the same size. Each inode has a i-number which states where in this contiguous array that it lives at struct dinode contains a type field which distinguishes between files, directories, and special files (ie. devices). A type of zero indicates that the inode is free nlink field counts the number of directory entries that refer to this inode. Needed to know when the on-disk inode and its data blocks should be freed size field counting the number of bytes of content in the file addrs array records the block number of the disk blocks holding the file’s content The active inodes in memory are kept in an itable a struct inode is only stored if there is a C pointer referring to that inode the refs field in inode keeps track of these C pointers. Once the number of pointers hits 0, the kernel discards the inode from memory There are 4 locking mechanisms in the inode code itable.lock protects the invariant that an inode is present in the inode table at most once, as well as the invariant that an in-memory inode’s ref field counts the number of in-memory pointers to the inode Each in-memory inode has a lock field containing a sleep-lock, which ensure exclusive access to the inode’s fields a non-zero ref field forces the system to maintain that inode in the table nlink counts the number of directory entries that refer to a file. Xv6 won’t free an inode if the link count is greater than 0 iget() returns an inode pointer that is gaurenteed to be valid until the matching iput() is called In words, iget() provides non-exclusive access to the the inode, which some processes depend on (think knowing if a file is open or not) If you want the inode structure to be populated with the most recent update of the inode, then you need to call ilock() and iunlock() to get exclusive ownership, and then read from the disk This seperation helps avoid deadlock in some scenarios (say directory lookup) Directory Layer Directories are implemented internally like a file (ie. inodes) Directories have a linked list of names and inodes numbers. Directories with inode number zero are free the function dirlookup() searches a directory for an entry with a given name. If it finds one, it returns a pointer to the correesponding inode and sets *poff to the byte offset of the entry within the directory The function dirlink() write a new directory entry into the specified directory by looping through the directory entries; If an unallocated entry is found, the inode is created at that offset Path Layer Path name lookup involves sucessive calls to dirlookup. The starting point is either the root directory or the current directory After locking the current directory pointer, you scan through the entries, set the current directory to the found directory, and then interate on that new directory There are some concurrency issues: One kernel thread is looking up a pathname while another kernel thread may be changing the directory tree by unlinking a directory. A potential risk is that a lookup may be searching a directory that has been deleted by another thread Locking a directory when doing path lookup and reference counting inodes are some deadlock prevention mechanisms File Descriptor Layer Everything is a file in UNIX. Everything! This layer helps acheive this uniformity Each process gets its own table of open files (file descriptors) struct file is a wrapper around inode or pipe, plus an I/O offset each open() cll create a new open file Multiple processes can have the same open file, but different I/O offsets A single file can appear multiple times in one process’s file table; it can also appear across multiple processes This could happen with dup or fork All open files are kept in a global file table (ftable) Filealloc scans the file table for an unreferenced file and returns a new reference Filedup increments the reference count fileclose decrements this counter. When the counter hits 0, the underlying pipe or inode is released