Two Dimensional Treatment of the Magnetic Field due to a circular coil (ring of current)

If we consider the (OverVector[B]) Field created by a ring of current oriented normal to the z axis, we expect that planes passing through the z axis at various θ would exhibit the same (OverVector[B]) Field plots ((OverVector[B]) is symmetric with respect to rotations in the θ direction).  Thus, a 2D representation of (OverVector[B]) in any of these planes is sufficient to describe it's behavior.

Here we show the current loop normal to the x-z plane, flowing around the z axis and a possible representation of the (OverVector[B]) field streamlines in the x-z plane.   If you can imagine other planes passing through the z axis, similar (OverVector[B]) field streamlines would be expected to occur on these planes due to the θ symmetry.   Realize that the ring is actually in the x-y plane (the blue “ring” is misleading -- it’s intersection with the x-z plane would just be two points on the x axis (corresponding to the points {R, 0, 0} and {-R,0,0}).

After some interpretation, we will extend the problem to 3 Dimensions.  (you can relax; we do the work).

Using the Biot-Savart Law, we wish to derive the appropriate (OverVector[B])[x,z].   We will use an approach presented by Phil Duxbury (Michigan State U.) in this notebook:  http://computation.pa.msu.edu/phy201/worksheet8nb/.

(a)  We know that Biot-Savart looks like this:    (OverVector[B]) = ∫(((Subscript[μ, o])I)/(4π))(((OverVector[dl])×(OverVector[ℝ]))/(ℝ^3)) where (OverVector[ℝ]) = (OverVector[r])ring - (OverVector[r])fieldpoint.    

(OverVector[r])ring is the radius vector of the differential element of the ring (OverVector[dl]) over which we are integrating.

(OverVector[r])fieldpoint is a point in the plane of interest passing through the z axis and cutting the ring similar to the plane shown in the figure above.  (This is the x-z plane and therefore (OverVector[r])fieldpoint = {x, 0, z}).

Enter below the code in M to define  (OverVector[r])ring , (OverVector[r])fieldpoint,    (OverVector[ℝ]), and the argument of the Biot-Savart integral:     ((((Subscript[μ, o])I)/(4π))(((OverVector[dl])×(OverVector[ℝ]))/(ℝ^3))).  

HUGE HELPERS:    Letting Ro represent the (constant) radius of the ring, the vector rring = Ro {Cos[theta], Sin[theta], 0}   and the vector dl = {0,0,1} × rring  (I am leaving off the arrows above (OverVector[dl]) and (OverVector[r])ring.)  

Be sure you verify these ‘helpers’.  They are written in Cartesian Coordinates and theta is an angular parameter that takes you around the ring.  It will be our variable of integration.

Now enter the code:

ClearAll["Global`*"]
rring=Ro {Cos[theta],Sin[theta],0} ;
rfieldpoint = {x, 0, z};   (* we are examining B in the x-z plane; ∴ we set y = 0 *)
 dl={0,0,1} \[Cross] rring   (*  NOTE:  the smaller x:  \[Cross]  is for the vector cross-product *)
ℝ =  rfieldpoint - rring
ℝmag = (Sqrt[ℝ.ℝ])    (* ℝmag is the length of the VECTOR ℝ  *)
dB = ((μo II)/(4π))((dl\[Cross]ℝ)/((ℝmag)^3))  //Simplify    (*  You could also write for dl\[Cross]ℝ:  Cross[dl,ℝ];   dB is a VECTOR  *)

(b)  Hanging firmly onto your hat, integrate these three terms using M to obtain expressions for B in the x-z plane.  (M fearlessly attacks these integrals and will generate solutions involving Elliptical Integrals).

Important note:  It took M too long to do the definite integrals,  so I did the indefinites first, then evaluated the limits (0 to 2π);  worked like a champ.

I will first do each component separately, then combine into (hopefully) a 2D vector with x and z components only.  The integrals are over theta from 0 to 2 π.  M does the indefinite integral extremely fast whereas the definite integral is very slow to emerge.  Therefore I do the indefinite and then evaluate at the endpoints  (0 and 2 π)

(*  II = 1; Ro =1; μo = 1;  *)
dB[[1]]  (*  note that it is a function of theta *)
Bx[theta_] =Integrate [dB[[1]], theta]   (* dB[[1]] is the x component of dB;  this is the indefinite integral over theta *)
Bxx = Bx[2  π] - Bx[0]//Simplify  (*  Here is the final definite integral for the x component  *)

dB[[2]]  (*  note that it is a function of theta *)
By[theta_] =Integrate [dB[[2]], theta]   (* dB[[2]] is the y component of dB;  this is the indefinite integral *)
Byy = By[2  π] - By[0]//Simplify(*  Here is the final definite integral for the y component  *)

dB[[3]]  (*  note that it is a function of theta *)
Bz[theta_] =Integrate [dB[[3]], theta]   (* dB[[3]] is the z component of dB;  this is the indefinite integral *)
Bzz = Bz[2  π] - Bz[0]//Simplify  (*  Here is the final definite integral for the z component  *)

Messy! but M can evaluate the Elliptical Functions to very high precision, very quickly.  NOTE that Byy = 0 so the B vectors lie in the plane chosen ( the x-z plane).

(c)  From earlier work we know that (OverVector[B]) along the symmetry axis (the z axis) (OverVector[B])zaxis = ((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2))))(OverHat[z]).  Check to see if the above mess reduces to this result.

Bzaxis = Limit[ {Bxx, Bzz}, x-> 0]

(d)  From earlier work we know that (OverVector[B]) at the center of the current loop (the origin) = ((II μo)/(2  Ro))(OverHat[z]).  Check to see if the above reduces to this result.

Bcenter = Limit[{Bxx, Bzz}, {x-> 0,z-> 0},Assumptions-> Ro >0]

(e)  Generate a VectorPlot and StreamPlot for  B in the x-z plan using the results from the integrations above.  Be sure you understand where the current loop is sitting in space.

II = 1; Ro = 1; μo = 1;   
VectorPlot[{Bxx,Bzz},{x,-2,2},{z,-2,2},VectorPoints->22,VectorScale->{Medium,Scaled[0.7]},VectorStyle->Yellow,FrameStyle->White,Background->Black,FrameLabel->{x,z}]
StreamPlot[{Bxx,Bzz},{x,-2,2},{z,-2,2},StreamPoints->26,StreamColorFunction->Hue,FrameLabel->{x,z}]

Where is the current loop?:  It intersects the x-z plane inside the ‘circles’ (at {0, ± 1}) and lies normal to this plane.  It’s symmetry axis is the z axis.

(f)  Explore the Divergence of B and comment using M’s Div function.      THINGS BREAK APART!!!!   Div not appear to be 0!

Clear[II, Ro, μo]
divvy = Div[{Bxx,Byy, Bzz},{x,y, z}]//Simplify

Doesn’t look like zero;  if we do a plot of Div[(OverVector[B])] we see that it is not zero!   What is wrong???

II = 1; Ro = 1; μo = 1;   
Plot3D[divvy,{x,-2,2},{z,-2,2},PlotRange-> {-1.2,1.2}]

This type of behavior arises all the time when you reduce the problem to 2 (or 1) dimension(s).  

IN ONE D:  Recall that (OverVector[B]) along the symmetry axis of the ring (examined above) is (OverVector[B])[z] = ((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2))))(OverHat[z]).   Since this is a function of z only (and (OverVector[B]) points in the z direction only) we would assume that Div[(OverVector[B])[z]] is given by:   D[B[z], z]  ((OverVector[B])[z] is a vector; B[z] is the functional form of the length of (OverVector[B])[z]).  

DO IT:

D[((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2)))),z]

Here (in the 1D case) and in the 2D case we are ignoring the fact that there are components of B in the x and/or y direction adjacent to the points where we are evaluating the divergence.  E.g., for (OverVector[B]) along the symmetry axis (z), even though (OverVector[B]) is only in the z direction, if we are not careful we might conclude that Div B ≠ 0 (which means there are magnetic monopoles -- experimentally, so far, their ain’t any!).

Reminder:  Form of Div in Cartesian Coordinates.  

Take the “dummy” vector function BB = {BBx[x, y, z], BBy[x, y, z], BBz[x, y, z]} and take apply Div:

Div[{BBx[x,y,z], BBy[x,y,z], BBz[x,y,z]},{x,y,z}]

Converting back to (OverVector[B]) , The term(s) we are missing in the divergence is (By^((0,1,0)))[x, y, z] which = D[By[x,y,z], y] and (Bz^((0,0,1)))[x, y, z] which = D[Bz[x,y,z], z] 

We expect Div[(OverVector[B]) ] = 0 when we have all three components of (OverVector[B])  correct.  SO,  let’s try to solve the 3D problem for the ring of current.

We start just like above but we let the field point be 3D:    rfieldpoint = {x, y, z};
Again, theta is being treated like a parameter which takes us around the current loop (it actually is the polar angle in the x-y plane).

ClearAll["Global`*"]
rring=Ro {Cos[theta],Sin[theta],0} ;  (*parametric version of the source points on the ring *)
rfieldpoint = {x, y, z};   (* we are examining B NOW in 3D so rfieldpoint = {x,y,z} *)
 dl={0,0,1} \[Cross] rring   (*  NOTE:  the smaller x:  \[Cross]  is for the vector cross-product *)
ℝ =  rfieldpoint - rring
ℝmag = (Sqrt[ℝ.ℝ])    (* ℝmag is the length of the VECTOR ℝ  *)
dB = ((μo II)/(4π))((dl\[Cross]ℝ)/((ℝmag)^3))  //Simplify

M does better on the calculus if we convert the dB we have in Cartesian coordinates to dB = dBcyl in Cylindrical Coordinates and then integrate.  We use M’s TransformedField to do this coordinate transformation:

dBcyl = TransformedField["Cartesian"->"Cylindrical", dB,{x,y,z}->{s,ϕ,ζ}]//Simplify  (* a 3D Vector *)

(*  NOW  integrate to find the three components of B in cylindrical coordinates  Bss, Bϕϕ, and Bζζ   where ζ is the Greek z  *) 
$Assumptions = {II,Ro,μo,theta,s,ϕ, ζ} ∈ Reals && 0≤ {theta,ϕ}≤ 2π;
dBcyl[[1]]  (*  note that it is a function of theta *)
Bss[theta_] =Integrate[dBcyl[[1]], theta]//Simplify ; (*  this is the indefinite integral *)
Bss[s_,ϕ_,ζ_] = (Bss[2 π ] - Bss[0])  //Simplify   (*  the definite integral over 0 ≤ theta ≤ 2  π *)

dBcyl[[2]]  (*  note that it is a function of theta *)
Bϕϕ[theta_] =Integrate[dBcyl[[2]], theta]//Simplify ; (*  this is the indefinite integral *)
Bϕϕ[s_,ϕ_,ζ_] = (Bϕϕ[2 π ] - Bϕϕ[0])  //Simplify   (*  the definite integral over 0 ≤ theta ≤ 2  π *)

We should get Bϕϕ = 0 due to the cylindrical symmetry of the current distribution.

dBcyl[[3]]  (*  note that it is a function of theta *)
Bζζ[theta_] =Integrate[dBcyl[[3]], theta]//Simplify ; (*  this is the indefinite integral *)
Bζζ[s_,ϕ_,ζ_] = (Bζζ[2 π ] - Bζζ[0])  //Simplify   (*  the definite integral over 0 ≤ theta ≤ 2  π *)

Div[B] = 0 in (3D) cylindrical coordinates.  [Hurray!]   No Monopoles.

BB[s_,ϕ_,ζ_] = {Bss[s,ϕ,ζ],Bϕϕ[s,ϕ,ζ],Bζζ[s,ϕ,ζ]};
Div[BB[s,ϕ,ζ] ,{s,ϕ,ζ},"Cylindrical"]//FullSimplify

Expect (OverVector[B]) along the symmetry axis (the z axis) (OverVector[B])zaxis  to equal:   ((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2))))(OverHat[z]).   I.e., Bss and Bϕϕ = 0 
( we showed above that Bϕϕ = 0 for all points in space;  performing the following shows that Bss = 0 for s = 0.  

I.e., we get:  (OverVector[B])zaxis = {0,0,((II (Ro^2) μo)/(2 ((((Ro^2)+(ζ^2)))^(3/2))))}. )

Limit[BB[s, ϕ, ζ],s->0]   (* if you get a Warning about assumptions, ignore *)

We can find B at the center of the ring:

Limit[%, ζ -> 0,Assumptions->Ro>0]//Simplify

It is redundant but we can quickly show that Bss = 0 for s = 0

Limit[Bss[s, ϕ, ζ],s->0]  (* if you get a Warning about assumptions, ignore *)

Plot B in cylindrical coordinates for case of ϕ = 0.   This will be like a Cartesian plot where s -> x  and ζ -> z.  Note that s is always positive so we can only represent  “half” the x - z plane:

Clear[II,Ro,μo,s,ϕ,ζ]
II = 1; Ro = 1; μo = 1;  
BB2D[s_, ζ_]= {BB[s, ϕ, ζ][[1]],BB[s, ϕ, ζ][[3]]}/. ϕ -> 0
StreamPlot[BB2D[s,ζ] ,{s,0.01,2},{ζ,-2,2},AspectRatio->Automatic]

The current loop is coming out/going into the paper at the point {1,0}, inside the ‘circle’.  For the Field direction shown by the streamlines, the current is going into the paper.

Transform to Cartesian Coordinates - This is the full 3D solution for BB in {x,y,z} coordinates.  The functions can be simplified:   This article derives fairly simple expressions for Cartesian, Spherical, and Cylindrical coordinates in 3D

Clear[II,Ro,μo]
Bcart3D[x_,y_, z_] = TransformedField["Cylindrical"->"Cartesian",BB[s, ϕ, ζ],{s,ϕ,ζ}->{x,y,z}]//Simplify

Div[B] works in (3D) Cartesian coordinates just fine:

Div[Bcart3D[x,y, z] ,{x,y,z},"Cartesian"]//Simplify

Bcard3D evaluated along the z axis  is, as before  [we should get:    {0,0,((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2))))} =  ((II (Ro^2) μo)/(2 ((((Ro^2)+(z^2)))^(3/2)))) (OverHat[z])]
Letting M actually Do It (in Cartesian Coordinates):

Limit[Bcart3D[x,y,z],{x-> 0, y -> 0}]     (*  THIS TAKES A LONG TIME -on an iMac, ~2 minutes  *)

We can plot Bcard3D either in 3D:  (Blind use of VectorPlot3D -- it’s not pretty)

II = 1; Ro = 1; μo = 1;  
VectorPlot3D[Bcart3D[x,y,z], {x,-2,2},{y,-2,2},{z,-2,2}]  (* 3D:  this graph is a mess -- much playing around might help  *)

We can plot Bcard3D in 2D as we did above.  we set y -> 0.

(*  2D Plots:  Just set y -> 0;  Results will look just like the 2D results obtained above *)
Bxx = Bcart3D[x,y,z][[1]]/.y->0;   (* remove ; if you want to see functional forms  *)
Bzz = Bcart3D[x,y,z][[3]]/.y->0;
VectorPlot[{Bxx,Bzz},{x,-2,2},{z,-2,2},VectorPoints->22,VectorScale->{Medium,Scaled[0.7]},VectorStyle->Yellow,FrameStyle->White,Background->Black,FrameLabel->{x,z}]
StreamPlot[{Bxx,Bzz},{x,-2,2},{z,-2,2},StreamPoints->26,StreamColorFunction->Hue]

Exactly as the 2D results obtained above.