Sine and Cosine Integrals and the Delta Function

Trigger warning: This post is for mature audiences only; it may contain Fourier transforms and Dirac delta functions. 

In Chapter 11 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I examine some properties of Fourier transforms. In particular, we consider three integrals of sines and cosines. After some analysis, we conclude that these integrals are related to the Dirac delta function, δ(ω – ω’), equal to infinity at ω = ω’ and zero everywhere else (it’s a strange function consisting of one really tall, thin spike).

Are these equations correct? I now believe that they’re almost right, but not entirely. I propose that instead they should be 

You’re probably thinking “what a pity, the second three equations looks more complicated than the first three.” I agree. But let me explain why I think they’re better. Hang on, it’s a long story.

Let’s go back to our definition of the Fourier transform in Eq. 11.57 of IPMB

The first thing to note is that y(t) consists of two parts. The first depends on cos(ωt), which is an even function, meaning cos(–ωt) = cos(ωt). There’s an integral over ω, implying that many different frequencies contribute to y(t), weighted by the function C(ω). But one thing we know for sure is that when you add up the contributions from all these many frequencies, the result must be an even function (the sum of even functions is an even function). The second part depends on sin(ωt), which is an odd function, sin(–ωt) = – sin(ωt). Again, when you add up all the contributions from these many frequencies weighted by S(ω), you must get an odd function. So, we can say that we’re writing y(t) as the sum of an even part, y_even(t), and an odd part, y_odd(t). In that case, we can rewrite our Fourier transform expressions as

We should be able to take our expression for y_even(t), put our expression for C(ω) into it, and then—if all works as it should—get back y_even(t). Let’s try it and see if it works. To start I’ll just rewrite the first of the four equations listed above

Now for C(ω) I’ll use the third of the four equations listed above. In that expression, there is an integral over t, but t is a dummy variable (it’s an “internal” variable; after you do the integral, the result does not depend on t), so to keep things from getting confusing we’ll call the dummy variable by another name, t'

Next switch the order of the integrals, so the integral over t' is on the outside and the integral over ω is on the inside

Ha! There, inside the bracket, is one of those integrals were’re talking about. Okay, the variables ω and t are swapped, but otherwise it’s the same. So, let’s put in our new expression for the integral

The 2π’s cancel, and a factor of one half comes out. An integral containing a delta function just picks out the value where the argument of the delta function is zero. We get

But, we know that y_even(t) is an even function, meaning y_even(–t) equals y_even(t). So finally

It works! We go “around the loop” and get back our original function.

You could perform another calculation just like this one but for y_odd(t). Stop reading and do it, to convince yourself that again you get back to where you started from, y_odd(t) = y_odd(t).

Now, you folks who are really on the ball might realize that if you had used the old delta function relationships given in IPMB (the first three equations in this post), they would also work. (Again, try it and see.) So why use my fancy new formulas? Well, if you have an integral that adds up a bunch of cos(ωt), you know you’re gonna get an even function. There’s no way it can be equal to δ(ω – ω’), because that function is neither even nor odd. So, it just doesn’t make sense to say that summing up a bunch of even functions gives you something that isn’t even. In my new formula, that sum of two delta functions is an even function. The same argument holds for the integral with sin(ωt), which must be odd.

Finally (and this is what got me started down this rabbit hole), you often see the delta function written as

Jackson even gives this equation, so it MUST be correct. (For those of who aren’t physicists, John David Jackson wrote the highly regarded graduate textbook Classical Electrodynamics, known by physics graduate students simply as “Jackson.”)

In Jackson’s equation, i is the square root of minus one. So, this representation of the delta function uses complex numbers. You won’t see it in IPMB because Russ and I avoid complex numbers (I hate them).

Let’s use the Euler formula e^iθ = cosθ + i sinθ to change the integral in Jackson’s delta function expression to

Now use a couple trig identities, cos(A–B) = cosA cosB + sinA sinB and sin(A–B) = sinA cosB –cosA sinB, to get

This is really four integrals,

Then, using the relations between these integrals and the delta function given in IPMB (the first three equations at the top of this post), you get that the sum of these integrals is equal to

which is obviously wrong; we started with 2πδ and ended up with 4πδ. Even worse, do the same calculation for δ(ω + ω') with a plus instead of a minus in front of the ω'. I’ll leave it to you to work out the details, but you’ll get zero! Again, nonsense. However, if you use the integral relations I propose above (second set of three integrals at the top of this blog), everything works just fine (try it).

Gene Surdutovich, my new coauthor for the sixth edition of IPMB, and I are still deciding how to discuss this issue in the new edition (which we are hard at work on but I doubt will be out within the next year). I don’t want to get bogged down in mathematical minutia that isn’t essential to our book’s goals, but I want our discussion to be correct. Once the sixth edition is published, you can see how we handle it.

I haven’t seen my new delta function/Fourier integral relationships in any other textbook or math handbook. This makes me nervous. Are they correct? Moreover, Intermediate Physics for Medicine and Biology does not typically contain new mathematical results. Maybe I haven’t looked hard enough to see if someone else published these equations (if you’ve seen them before, let me know where…please!). Maybe I’ll find these results in Morse and Feshbach (another one of those textbooks known to all physics graduate students) or some other mathematical tome. I need to make a trip to the Oakland University library to look through their book collection, but right now its too cold and snowy (we got about four to five inches of the white stuff in the last 48 hours).

A Toy Model of Climate Change

Introduction

A screenshot of the online book
Math for the People.

In Introductory Physics for Medicine and Biology, Russ Hobbie and I make use of toy models. Such mathematical descriptions are not intended to be accurate or realistic. Rather, they‘re simple models that capture the main idea without getting bogged down in the details. Today, I present an example of a toy model. It’s not related to medicine or biology, but instead describes climate change. I didn’t originally derive this model. Much of the analysis below comes from other sources, such as the online book Math for the People published by Mark Branson and Whitney George.

Earth Without an Atmosphere

First, consider the earth with no atmosphere. We will balance the energy coming into the earth from the sun with the energy from the earth that is radiated out into space. Our goal will be to calculate the earth’s temperature, T.

The power density (energy per unit time per unit area, in watts per square meter) emitted by the sun is called the solar constant, S. It depends on how far you are from the sun, but at the earth’s orbit S = 1360 W/m². To get the total power impinging on our planet, we must multiply S by the area subtended by the earth, which is πR², where R is the earth’s radius (R = 6.4 × 10⁶ meters). This gives SπR² = 1.8 × 10¹⁷ W, or nearly 200,000 TW (T, or tera-, means one trillion). That’s a lot of power. The total average power consumption by humanity is only about 20 TW, so there’s plenty of energy from the sun.

We often prefer to talk about the energy loss or gain per unit area of the earth’s surface. The surface area of the earth is 4πR² (the factor of four comes from the total surface area of the spherical earth, in contrast to the area subtended by the earth when viewed from the sun). The power per unit area of the earth’s surface is therefore SπR²/4πR², or S/4.

Not all of this energy is absorbed by the earth; some is reflected back into space. The albedo, a, is a dimensionless number that indicates the fraction of the sun’s energy that is reflected. The power absorbed per unit area is then (1 – a)S/4. About 30% of the sun’s energy is reflected (a = 0.3), so the power of sunlight absorbed by the earth per unit of surface area is 238 W/m².

What happens to that energy? The sun heats the earth to a temperature T. Any hot object radiates energy. Such thermal radiation is analyzed in Section 14.8 of Intermediate Physics for Medicine and Biology. The radiated power per unit area is equal to eσT⁴. The symbol σ is the Stefan-Boltzmann constant, σ = 5.7 × 10^–8 W/(m² K⁴). As stated earlier, T is the earth’s temperature. When raising the temperature to the fourth power, T must be expressed as the absolute temperature measured in kelvin (K). Sometimes it’s convenient at the end of a calculation to convert kelvin to the more familiar degrees Celsius (°C), where 0°C = 273 K. But remember, all calculations of T⁴ must use kelvin. Finally, e is the emissivity of the earth, which is a measure of how well the earth absorbs and emits radiation. The emissivity is another dimensionless number ranging between zero and one. The earth is an excellent emitter and absorber, so e = 1. From now on, I’ll not even bother including e in our equations, in which case the power density emitted is just σT⁴.

Let’s assume the earth is in steady state, meaning the temperature is not increasing or decreasing. Then the power in must equal the power out, so 

(1 – a)S/4 = σT⁴. 

Solving for the temperature gives

T = ∜[(1 – a)S/4σ] .

Because we know a, S, and σ, we can calculate the temperature. It is T = 254 K = –19°C. That’s really cold (remember, in the Celsius scale water freezes at 0°C). Without an atmosphere, the earth would be a frozen wasteland.

Earth With an Atmosphere

Often we can learn much from a toy model by adding in complications, one by one. Now, we’ll include an atmosphere around earth. We must keep track of the power into and out of both the earth and the atmosphere. The earth has temperature T_E and the atmosphere has temperature T_A.

First, let’s analyze the atmosphere. Sunlight passes right through the air without being absorbed because it’s mainly visible light and our atmosphere is transparent in the visible part of the spectrum. The main source of thermal (or infrared) radiation (for which the atmosphere is NOT transparent) is from the earth. We already know how much that is, σT_E⁴. The atmosphere only absorbs a fraction of the earth’s radiation, e_A, so the power per unit area absorbed by the atmosphere is e_AσT_E⁴.

Just like the earth, the atmosphere will heat up to a temperature T_A and emit its own thermal radiation. The emitted power per unit area is e_AσT_A⁴. However, the atmosphere has upper and lower surfaces, and we’ll assume they both emit equally well. So the total power emitted by the atmosphere per unit area is 2e_AσT_A⁴.

If we balance the power in and out of the atmosphere, we get 

e_AσT_E⁴ = 2e_AσT_A⁴. 

Interestingly, the fraction of radiation absorbed by the atmosphere, e_A, cancels out of our equation (a good emitter is also a good absorber). The Stefan-Boltzmann constant σ also cancels, and we just get T_E⁴ = 2T_A⁴. If we take the forth root of each side of the equation, we find that T_A = 0.84 T_E. The atmosphere is somewhat cooler than the earth.

Next, let’s reanalyze the power into and out of the earth when surrounded by an atmosphere. The sunlight power per unit area impinging on earth is still (1 – a)S/4. The radiation emitted by the earth is still σT_E⁴. However, the thermal radiation produced by the atmosphere that is aimed inward toward the earth is all absorbed by the earth (since the emissivity of the earth is one, e_E = 1), so this provides another factor of e_AσT_A⁴. Balancing power in and out gives

(1 – a)S/4 + e_AσT_A⁴ = σT_E⁴ .

Notice that if e_A were zero, this would be the same relationship as we found when there was no atmosphere: (1 – a)S/4 = σT_E⁴. The atmosphere provides additional heating, warming the earth.

We found earlier that T_E⁴ = 2T_A⁴. If we rewrite this as T_A⁴ = T_E⁴/2 and plug that into our energy balance equation, we get

(1 – a)S/4 + e_AσT_E⁴/2 = σT_E⁴ .

With a bit of algebra, we find

(1 – a)S/4 = σT_E⁴ (1 – e_A/2) .

Solving for the earth’s temperature gives

T_E = ∜[(1 – a)S/4σ] ∜[1/(1 – e_A/2) ] .

If e_A were zero, this would be exactly the relationship we had for no atmosphere. The fraction of energy absorbed by the atmosphere is not zero, however, but is approximately e_A = 0.8. The atmosphere provides a dimensionless correction factor of ∜[1/(1 – e_A/2)]. The temperature we found previously, 254 K, is corrected by this factor, 1.136. We get T_E = 288.5 K = 15.5 °C. This is approximately the average temperature of the earth. Our atmosphere raised the earth’s temperature from –19°C to +15.5°C, a change of 34.5°C.

Climate Change

To understand climate change, we need to look more deeply into the meaning of the factor e_A, the fraction of energy absorbed by the atmosphere. The main constituents of the atmosphere—oxygen and nitrogen—are transparent to both visible and thermal radiation, so they don’t contribute to e_A. Thermal energy is primarily absorbed by greenhouse gases. Examples of such gases are water vapor, carbon dioxide, and methane. Methane is an excellent absorber of thermal radiation, but its concentration in the atmosphere is low. Water vapor is a good absorber, but water vapor is in equilibrium with liquid water, so it isn’t changing much. Carbon dioxide is a good absorber, has a relatively high concentration, and is being produced by burning fossil fuels, so a lot of our discussion about climate change focuses on carbon dioxide.

The key to understanding climate change is that greenhouse gasses like carbon dioxide affect the fraction of energy absorbed, e_A. Suppose an increase in the carbon dioxide concentration in the atmosphere increased e_A slightly, from 0.80 to 0.81. The correction factor  ∜(1/(1 – e_A/2) ) would increase from 1.136 to 1.139, changing the temperature from 288.5 K to 289.3 K, implying an increase in temperature of 0.8 K. Because changes in temperature are the same if expressed in kelvin or Celsius, this is a 0.8°C rise. A small change in e_A causes a significant change in the earth’s temperature. The more carbon dioxide in the atmosphere, the greater the temperature rise: Global warming.

Feedback

We have assumed the earth’s albedo, a, is a constant, but that is not strictly true. The albedo depends on how much snow and ice cover the earth. More snow and ice means more reflection, a larger albedo, a smaller amount of sunlight absorbed by the earth, and a lower temperature. But a lower temperature means more snow and ice. We have a viscous cycle: more snow and ice leads to a lower temperature which leads to more snow and ice, which leads to an even lower temperature, and so on. Intermediate Physics for Medicine and Biology dedicates an entire chapter to feedback, but it focuses mainly on negative feedback that tends to maintain a system in equilibrium. A viscous cycle is an example of positive feedback, which can lead to explosive change. An example from biology is the upstroke of a nerve action potential: an increase in the electrical voltage inside a nerve cell leads to an opening of sodium channels in the cell membrane, which lets positively charged sodium ions enter the cell, which causes the voltage inside the cell to increase even more. The earth’s climate has many such feedback loops. They are one of the reasons why climate modeling is so complicated.

Conclusion

Today I presented a simple description of the earth’s temperature and the impact of climate change. Many things were left out of this toy model. I ignored differences in temperature over the earth’s surface and within the atmosphere. I neglected ocean currents and the jet stream that move heat around the globe. I did not account for seasonal variations, or for other greenhouse gasses such as methane and water vapor, or how the amount of water vapor changes with temperature, or how clouds affect the albedo, and a myriad of other factors. Climate modeling is a complex subject. But toy models like I presented today provide insight into the underlying physical mechanisms. For that reason, they are crucial for understanding complex phenomena such as climate change.

Taylor Series

The Taylor series is particularly useful for analyzing how functions behave in limiting cases. This is essential when translating a mathematical expression into physical intuition, and I would argue that the ability to do such translations is one of the most important skills an aspiring physicist needs. Below I give a dozen examples from Intermediate Physics for Medicine and Biology, selected to give you practice with Taylor series. In each case, expand the function in the dimensionless variable that I specify. For every example—and this is crucial—interpret the result physically. Think of this blog post as providing a giant homework problem about Taylor series.

Find the Taylor series of:

1. Eq. 2.26 as a function of bt (this is Problem 26 in Chapter 2). The function is the solution for decay plus input at a constant rate. You will need to look up the Taylor series for an exponential, either in Appendix D or in your favorite math handbook. I suspect you’ll find this example easy.
2. Eq. 4.69 as a function of ξ (this is Problem 47 in Chapter 4). Again, the Taylor series for an exponential is required, but this function—which arises when analyzing drift and diffusion—is more difficult than the last one. You’ll need to use the first four terms of the Taylor expansion.
3. The argument of the inverse sine function in the equation for C(r,z) in Problem 34 of Chapter 4, as a function of z/a (assume r is less than a). This expression arises when calculating the concentration during diffusion from a circular disk. Use your Taylor expansion to show that the concentration is uniform on the disk surface (z = 0). This calculation may be difficult, as it involves two different Taylor series. 
4. Eq. 5.26 as a function of ax. Like the first problem, this one is not difficult and merely requires expanding the exponential. However, there are two equations to analyze, arising from the study of countercurrent transport. 
5. Eq. 6.10 as a function of z/c (assume c is less than b). You will need to look up or calculate the Taylor series for the inverse tangent function. This expression indicates the electric field near a rectangular sheet of charge. For z = 0 the electric field is constant, just as it is for an infinite sheet.
6. Eq. 6.75b as a function of b/a. This equation gives the length constant for a myelinated nerve axon with outer radius b and inner radius a. You will need the Taylor series for ln(1+x). The first term of your expansion should be the same as Eq. 6.75a: the length constant for an unmyelinated nerve with radius a and membrane thickness b.
7. The third displayed equation of Problem 46 in Chapter 7 as a function of t/t_C. This expression is for the strength-duration curve when exciting a neuron. Interestingly, the short-duration behavior is not the same as for the Lapicque strength-duration curve, which is the first displayed equation of Problem 46.
8. Eq. 9.5 as a function of [M']/[K]. Sometimes it is tricky to even see how to express the function in terms of the required dimensionless variable. In this case, divide both sides of Eq. 9.5, to get [K']/[K] in terms of [M']/[K]. This problem arises from analysis of Donnan equilibrium, when a membrane is permeable to potassium and chloride ions but not to large charged molecules represented by M’.
9. The expression inside the brackets in Eq. 12.42 as a function of ξ. The first thing to do is to find the Taylor expansion of sinc(ξ), which is equal to sin(ξ)/ξ. This function arises when solving tomography problems using filtered back projection.
10. Eq. 13.39 as a function of a/z. The problem is a little confusing, because you want the limit of large (not small) z, so that a/z goes to zero. The goal is to show that the intensity falls off as 1/z² for an ultrasonic wave in the Fraunhoffer zone.
11. Eq. 14.33 as a function of λk_BT/hc. This problem really is to determine how the blackbody radiation function behaves as a function of wavelength λ, for short wavelength (high energy) photons. You are showing that Planck's blackbody function does not suffer from the ultraviolet catastrophe.
12. Eq. 15.18 as a function of x. (This is Problem 15 in Chapter 15). This function describes how the Compton cross section depends on photon energy. Good luck! (You’ll need it).

Brook Taylor

Who was Taylor? Brook Taylor (1685-1731) was an English mathematician and a fellow of the Royal Society. He was a champion of Newton’s version of the calculus over Leibniz’s, and he disputed with Johann Bernoulli. He published a book on mathematics in 1715 that contained his series.

Numerical Solution of the Quadratic Formula

Homework Problem 7 in Chapter 11 of Intermediate Physics for Medicine and Biology examines fitting data to a straight line. In that problem, the four data points to be fit are (100, 4004), (101, 4017), (102, 4039), and (103, 4063). The goal is to fit the data to the line y = ax + b. For this data, one must perform the calculation of a and b to high precision or else you get large errors. The solution manual (available to instructors upon request) says that

This problem illustrates how students can run into numerical problems if they are not careful. With modern calculators that carry many significant figures, this may seem like a moot point. But the idea is still important and can creep subtly into computer computations and cause unexpected, difficult-to-debug errors.

Numerical Recipes

Are there other examples of numerical errors creeping into calculations? Yes. You can find one discussed in Numerical Recipes that involves the quadratic formula.

We all know the quadratic formula from high school. If you have a quadratic equation of the form

the solution is

For example,
 

has two solutions
 

so x = 1 or 2.

Now, suppose the coefficient b is larger,

The solution is

so x = 300 or 0.00667.

This calculation is susceptible to numerical error. For instance, suppose all numerical calculations are performed to only four significant figures. Then when you reach the step

you must subtract 8 from 90,000. You get 89992, which to four significant figures becomes 89990, which has a square root of (again to four significant figures) 300.0. The solutions are therefore x = 300 or 0. The large solution (300) is correct, but the small solution (0 instead of 0.00667) is completely wrong. The main reason is that when using the minus sign for ± you must subtract two numbers that are almost the same (in this case, 300 – 299.98667) to get a much smaller number.

You might say “so what! Who uses only four significant figures in their calculations?” Okay, try solving

where I increased b from 300 to 3000. You’ll find that using even six significant figures gives one nonsense solution (try it). As you make b larger and larger, the calculation becomes more and more difficult. The situation can cause unexpected, difficult-to-debug errors.

What’s the moral to this story? Is it simply that you must use high precision when doing calculations? No. We can do better. Notice that the solution is fine when using the plus sign in the quadratic equation. We need make no changes. It’s the negative sign that gives the problem,
 

Let’s try a trick; multiply the expression by a very special form of one:

Simplifying, we get

Voilà! The denominator has the plus sign in front of the square root, so it is not susceptible to numerical error. The numerator is simplicity itself. Try solving x² – 300x + 2 = 0 using math to four significant figures,

 

No error, even with just four sig figs. The problem is fixed!

I should note that the problem is fixed only for negative values of b. If b is positive, you can use an analogous approach to get a slightly different form of the solution (I’ll leave that as an exercise for the reader).

So, the moral of the story is: if you find that your numerical calculation is susceptible to numerical error, fix it! Look for a trick that eliminates the problem. Often you can find one.

Stirling's Approximation

I've always been fascinated by Stirling’s approximation,

ln(n!) = n ln(n) − n,

where n! is the factorial. Russ Hobbie and I mention Stirling’s approximation in Appendix I of Intermediate Physics for Medicine and Biology. In the homework problems for that appendix (yes, IPMB does has homework problems in its appendices), a more accurate version of Stirling’s approximation is given as

ln(n!) = n ln(n) − n + ½ ln(2π n) .

There is one thing that’s always bothered me about Stirling’s approximation: it’s for the logarithm of the factorial, not the factorial itself. So today, I’ll derive an approximation for the factorial. 

The first step is easy; just apply the exponential function to the entire expression. Because the exponential is the inverse of the natural logarithm, you get

n! = e^{n ln(n) − n + ½ ln(2π n)}

Now, we just use some properties of exponents

n! = e^{n ln(n)} e⁻ⁿ e^{½ln(2π n)}

n! = (e^ln(n))ⁿ e⁻ⁿ √(e^{ln(2π n)})

n! = nⁿ e⁻ⁿ √(2π n) 

And there we have it. It’s a strange formula, with a really huge factor (nⁿ) multiplied by a tiny factor (e⁻ⁿ) times a plain old modestly sized factor (√(2π n)). It contains both e = 2.7183 and π = 3.1416.

Let's see how it works.

n	n!	nn e−n √(2π n)	fractional error (%)
1	1	0.92214	7.8
2	2	1.9190	4.1
5	120	118.02	1.7
10	3.6288 × 106	3.5987 × 106	0.83
20	2.4329 × 1018	2.4228 × 1018	0.42
50	3.0414 × 1064	3.0363 × 1064	0.17
100	9.3326 × 10157	9.3249 × 10157	0.083

For the last entry (n = 100), my calculator couldn’t calculate 100¹⁰⁰ or 100!. To get the first one I wrote

100¹⁰⁰ = (10²)¹⁰⁰ = 10^{2 × 100} = 10²⁰⁰.

The calculator was able to compute e⁻¹⁰⁰ = 3.7201 × 10⁻⁴⁴, and of course the square root of 200π was not a problem. To obtain the actual value of 100!, I just asked Google.

Why in the world does anyone need a way to calculate such big factorials? Russ and I use them in Chapter 3 about statistical dynamics. There you have to count the number of states, which often requires using factorials. The beauty of statistical mechanics is that you usually apply it to macroscopic systems with a large number of particles. And by large, I mean something like Avogadro’s number of particles (6 × 10²³). The interesting thing is that in statistical mechanics you often need not the factorial, but the logarithm of the factorial, so Stirling's approximation is exactly what you want. But it’s good to know that you can also approximate the factorial itself. 

Finally, one last fact from Mr. Google. 1,000,000! = 8.2639 × 10^5,565,708. Wow!

Stirling’s Approximation

https://www.youtube.com/watch?v=IJ5N28-Ujno

The Helmholtz Coil and the Maxwell Coil

To do magnetic resonance imaging, you need a static magnetic field that is uniform and a switchable magnetic field that has a uniform gradient. How do you produce such fields? In this post, I explain one of the simplest ways: using a Helmholtz coil and a Maxwell coil.

Both of these are created using circular coils. The magnetic field B_z produced by a circular coil can be calculated using the law of Biot and Savart (see Chapter 8 of Intermediate Physics for Medicine and Biology)

where μ₀ is the permeability of free space (the basic constant of magnetostatics), I is the coil current, N is the number of turns, R is the coil radius, and z is the distance along the axis from the coil center.

The Helmholtz Coil

The Helmholtz coil consists of two circular coils in parallel planes, having the same axis and the same current in the same direction, that are separated by a distance d. Our goal will be to find the value of d that gives the most uniform magnetic field. By superposition, the magnetic field is 

 

To create a uniform magnetic field, we will perform a Taylor expansion of the magnetic field about the origin (z = 0). We will need derivatives of the magnetic field. The first derivative is

(The reader will have to fill in the missing steps when calculating these derivatives.) At z = 0, this derivative will go to zero. In fact, because the magnetic field is even about the z axis, all odd derivatives will be zero, regardless of the value of d.

The second derivative is

At z = 0, the two terms in the brackets are the same. Our goal is to have this term be zero, implying that the second order term in the Taylor series vanishes. This will happen if

or, in other words, d = R. This famous result says that for a Helmholtz pair the coil separation should equal the coil radius.

A Helmholtz coil produces a remarkably uniform field near the origin. However, it is not uniform enough for use in most magnetic resonance imaging machines, which typically have a more complex set of coils to create an even more homogeneous field. If you need a larger region that is homogeneous, you could always just use a larger Helmholtz coil, but then you would need more current to achieve the desired magnetic field at the center. A Helmholtz pair isn’t bad if you want to use only two reasonably sized coils.

The Maxwell Coil

The Helmholtz coil produces a uniform magnetic field, whereas the Maxwell coil produces a uniform magnetic field gradient. It consists of two circular coils, in parallel planes having the same axis, that are separated by a distance d, but which have current in the opposite directions. Again, our goal will be to find the value of d that gives the most uniform magnetic field gradient. The magnetic field is

The only difference between this case and that for the Helmholtz coil is the change in sign of the second term in the bracket. If z = 0, the magnetic field is zero. Moreover, the magnetic field is an odd function of z, so all even derivatives also vanish. The first derivative is

This expression gives us the magnitude of the gradient at the origin, but it doesn’t help us create a more uniform gradient. The second derivative is

This derivative is zero at the origin, regardless of the value of d. So, we have to look at the third derivative.

At z = 0, this will vanish if

or, in other words, d = √3 R = 1.73 R. Thus, the two coils have a greater separation for a Maxwell coil than for a Helmholtz coil. The Maxwell coil would be useful for producing the slice selection gradient during MRI (for more about the need for gradient fields in MRI, see Chapter 18 of IPMB).

Conclusion

Below is a plot of the normalized magnetic field as a function of z for the Helmholtz coil (blue) and the Maxwell coil (yellow). As you can see, the region with a uniform field or gradient is small. It depends on what level of accuracy you need, but if you are more than half a radius from the origin you will see significant deviations from homogeneity.
 

Russ Hobbie and I never discuss the Helmholtz coil in Intermediate Physics for Medicine and Biology. We don’t mention the Maxwell coil by name, but Problem 33 of Chapter 18 analyzed a Maxwell pair even if we don’t call it that.

The Maxwell coil is great for producing the magnetic field gradient dB_z/dz needed for slice selection in MRI, but how do you produce the gradients dB_z/dx and dB_z/dy needed during MRI readout and phase encoding? That, my friends, is a story for another post.

Decay Plus Input at a Constant Rate Revisited

In Chapter 2 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss the problem of decay plus input at a constant rate.

Suppose that in addition to the removal of y from the system at a rate –by, y enters the system at a constant rate a, independent of y and t. The net rate of change of y is given by

Then we go on to discuss how you can learn things about a differential equation without actually solving it.

It is often easier to write down a differential equation describing a problem than it is to solve it… However, a good deal can be learned about the solution by examining the equation itself. Suppose that y(0) = 0. Then the equation at t = 0 is dy/dt = a, and y initially grows at a constant rate a. As y builds up, the rate of growth decreases from this value because of the –by term. Finally when a – by = 0, dy/dt is zero and y stops growing. This is enough information to make the sketch in Fig. 2.13.

The equation is solved in Appendix F. The solution is

… The solution does have the properties sketched in Fig. 2.13, as you can see from Fig. 2.14.

Figure 2.13 looks similar to this figure

Sketch of the initial slope a and final value a/b of y when y(0) = 0. In this figure, a=b=1.

 And Fig. 2.14 looks like this

A plot of y(t) using Eq. 2.26, with a=b=1.

However, Eq. 2.26 is not the only solution that is consistent with the sketch in Fig. 2.13. Today I want to present another function that is consistent with Fig. 2.13, but does not obey the differential equation in Eq. 2.25.

Let’s examine how this function behaves. When bt is much less than one, the function becomes y = at, so it’s initial growth rate is a. When bt is much greater than one, the function approaches a/b. The sketch in Fig. 2.13 is consistent with this behavior.

Below I show both Eqs. 2.26 and 2.26’ in the same plot.

A plot of y(t) using Eq. 2.26 (blue) and Eq. 2.26' (yellow), with a=b=1.

The function in Eq. 2.26 (blue) approaches its asymptotic value at large t more quickly than the function in Eq. 2.26’ (yellow).

The moral of the story is that you can learn a lot about the behavior of a solution by just inspecting the differential equation, but you can’t learn everything (or, at least, I can’t). To learn everything, you need to solve the differential equation. 

By the way, if Eq. 2.26’ doesn’t solve the differential equation in Eq. 2.25, then what differential equation does it solve? The answer is

 How did I figure that out? Trial and error.