E-Field and Electric Potential for a thick shell of charge with a  ρ[r] = (Subscript[ρ, o]) Exp[-α r]

A thick, non-conducting spherical shell of inner radius a and outer radius b has a volume charge density

 ρ[r] = (Subscript[ρ, o]) Exp[-α r] 	for a ≤ r ≤ b; assume (Subscript[ρ, o]) and α are positive.
 ρ[r] = 0 		elsewhere (for r < a and r > b).
Find  (OverVector[E]) and V at all r.

(a)Just for reference, make a plot of ρ[r] (we need to assume some values for the constants;  I’ll use a = 1; b = 2; α = 0.5; ρo = 1)

ClearAll["`*"]
a = 1; b = 2; α = 0.5; ρo = 1;  (*  just to get you started *)

We want to Find  (OverVector[E]).  Let’s not beat around the bushes - we will use Gauss’s Law to find (OverVector[E]).

(b) But first, in the text cell below, show that the vector (OverVector[E]) is radial in direction and that the magnitude of (OverVector[E]) depends only on the distance r from the origin O, where r = |(OverVector[r])|, the magnitude of (OverVector[r]).  Remember r is a distance, therefore a scalar.  This allows us to write:   (OverVector[E])[ (OverVector[r]) ] = E[ r ] (OverHat[r]).  (Note carefully where the symbols have vector vs. scalar status.)  Motivation:  This greatly simplifies applying Gauss’s Law to find (OverVector[E]).

this is a text cell — type in your response  You might want to include a figure.
First, verify the radial direction of (OverVector[E]).

Show that the magnitude of (OverVector[E]) at the point P being a function of r the scalar r — the length of (OverVector[r])) only.

For obvious reasons, ONCE we go through this lengthy exercise (or similar arguments for other symmetries) we generally and happily say "because of symmetry, . . . " and quickly move on.

(c) State why we make these arguments.

this is a text cell — type in your response

(d) Find (OverVector[E])[(OverVector[r])] for r < a  (Probably will not need to use M).  A diagram would be useful.

<this is a text cell>

(e)  Find (OverVector[E])inside_shell[(OverVector[r])], i.e., for a ≤ r ≤ b  (Suggest you set up both sides of Gauss’s Law, then use M.)

this is a text cell — type in your analysis of the LHS and RHS of Gauss’s Law here:

HINTS:  I wrote out both sides of Gauss’s Law and then used Solve to find the function Eshell[r] (the magnitude of E inside the shell of charge)

The Solve statement looked like this:  sol = Solve[LHS ==RHS, Eshell[r]].  

The vector field (OverVector[E])inside_shell[(OverVector[r])] is simply Eshell[r] (OverHat[r])  (the magnitude times the unit vector).

(* Input Cell *)
ClearAll["`*"]
ρ[rprime_] = ρo Exp[-α rprime];  (* just a starter  *)

(f)  Now find (OverVector[E]) outside the shell, (OverVector[E])outside_shell[(OverVector[r])], that is, (OverVector[E]) for r > b.

this is a text cell — type in your reasoning here

More Hints:  Again, HINTS:  I wrote out both sides of Gauss’s Law and then used Solve to find the function Eoutsideshell[r] (the magnitude of E outside the shell of charge)

The Solve statement looked like this:  sol = Solve[LHS ==RHS, Eoutsideshell[r]].  

The vector field (OverVector[E])outside_shell[(OverVector[r])] is simply Eoutsideshell[r] (OverHat[r])  (the magnitude times the unit vector).

(* Input Cell *)

(g)  Use If functions to define the magnitude of (OverVector[E])[r],(which I write as EE[r]), for all r   (An alternative approach is to use PieceWise[ ]. )

( ok, as a gift — this will work:  EE[r_] = If[r < a,0,If[a≤r≤b,Eshell[r],If[r > b,Eoutsideshell[r]]]] )

Chose values for the constants and Plot EE[r] from r = 0 to some r > b. 

(I used the following values:    a = 1; b = 2; α = 0.5; ρo = 1; ϵo=1;)

(* Input Cell *)

Just for comparison, below is a plot of the Charge Enclosed (Qenclosed[r]) by a spherical Gaussian Surface of radius r for our particular ρ[r]

To beat a very dead horse, for a spherical symmetric ρ[r]:    EE[r] *(4 π (r^2)) = ((Qenclosed[r])/(ϵo))

Clear[a ,b ,α , ρo , ϵo]
Qenclosed[r_] = If[r < a,0,If[a≤r≤b, (∫_a^r)ρ[rprime]4 π (rprime^2)\[DifferentialD]rprime,If[ r > b,(∫_a^b)ρ[rprime]4 π (rprime^2)\[DifferentialD]rprime]]];
a = 1; b = 3; α = 0.5; ρo = 1; ϵo=1;
plotQ = Plot[Qenclosed[r],{r,0,6},AxesLabel->{"r","Qenclosed[r]"}];
lines2 = Graphics[{Line[{{a,0},{a,35}}],Line[{{b,0},{b,35}}],Text["r = a",{a,7},Background->LightRed],Text["r = b",{b,7},Background->LightGreen]}];
Show[plotQ,lines2]

The relation between this curve and the above plot of E[r] is:    E[r]  = ((Qenclosed[r])/(4 π ϵo  (r^2)))= k((Qenclosed[r])/(r^2))

Note that the constant value of Qenclosed for r ≥ b equals the total charge contained in the whole shell which is:

		Qtotal = (∫_a^b)ρ[rprime]4 π (rprime^2)drprime

(h)  Produce a VectorPlot of EE[r]·(Overscript[r, ⋀]) in two dimensions (i.e., in a plane passing through the polar axis of the sphere).  Remember, you will need to explicitly define 2D vector functions (for the three regions:  r < a; a ≤ r ≤ b, and r > b) and convert them to Cartesian coordinates to use M’s VectorPlot function.

(* Input Cell *)
Clear[a ,b ,α , ρo , ϵo]   (* this simply clears the values assigned to these constants  *)

(i) Interpret the resulting 2D Vector Plot for the given spherically symmetric charge distribution:

	     0 				for r < a
ρ[r] =       ρo Exp[-α r]         for a ≤ r ≤ b
                0                         	for r > b

(this is a text cell - type in your interpretation)

(j)  Now derive the electric potential V[r] for all regions.  It is very easy to show that if (OverVector[E])[r] = Er[r] (OverHat[r]), V will be a function of r only.  We will go with that.

(this is a text cell - type in your interpretation)
      First, it might be helpful if we copy and paste clean spherical coordinate versions of Er[r] (the magnitude of (OverVector[E])[r] = Er[r] (OverHat[r]) ) for the three regions:
      
      	       Eoutsideshell[r] =  ((((e^(-a α)) (2+a α (2+a α))+(e^(-b α)) (-2-b α (2+b α))) ρo)/((r^2) (α^3) ϵo))     	for r > b
      		 
Er[r] =         Eshell[r] =  ((((e^(-a α)) (2+a α (2+a α))+(e^(-r α)) (-2-r α (2+r α))) ρo)/((r^2) (α^3) ϵo))          		for a ≤ r ≤ b

                    Ecavity[r]  = 		    0                              				for r < a

FIRST do  r > b:

(*  Input cell *)
 Clear[a, α, b, ϵo, ρo]
Eoutsideshell[r_] =  ((((E^(-a α)) (2+a α (2+a α))+(E^(-b α)) (-2-b α (2+b α))) ρo)/((r^2) (α^3) ϵo));   (* a reminder *)
$Assumptions =(b≤  rprime ≤ ∞) && rprime ∈ Reals && r∈ Reals && r ≥ b;  (* useful for helping M integrate *)

SECOND, inside the shell (a ≤ r ≤ b):

To find Vshell[r], we integrate Eoutside[r] from ∞ to b (which equals Voutsideshell[b]) and then integrate Eshell[r] from b to r (to get what I call Vinsideonly[r]).  (Of course we insert the – sign.)  

Vshell[r] is the sum of the two contributions.  The limits on these integrals are consistent with moving the test charge from ∞ to an r inside the shell, a ≤ r ≤ b.  Note:  I have changed the assumptions to aid the Vinsideonly[r] integration.  I recommend using Simplify near or at end.

Third, inside the empty cavity (r < a):   Since there is no charge enclosed, the electric field in the cavity is zero and the extra work to move the test charge from r = a to anywhere in the empty cavity will also be zero.  Thus, the potential inside the cavity Vcavity[r] = Vshell[a] (a constant!).

Here are the calculations to obtain Vshell[r]

(*  Input cell  *)
  Eshell[r_] =  ((((E^(-a α)) (2+a α (2+a α))+(E^(-r α)) (-2-r α (2+r α))) ρo)/((r^2) (α^3) ϵo));
$Assumptions =( a ≤  rprime ≤ b) && rprime ∈ Reals && ( a ≤  r ≤ b) && r∈ Reals && r > 0&& α ∈Reals &&α > 0;

THIRD, inside the empty cavity (r < a):

(*  Input cell --enter expression for Vcavity[r_] *)

(k) Now plot V[r] for all three regions combined.

(*  Input cell *)

(l) Plot the magnitude of (OverVector[E])[r], V[r], the charge density ρ[r], and the Qenclosed[r] together in a single graph. You may want to adjust some of the amplitudes so they stay usefully scaled.

Plotting the magnitude of E[r], V[r], the charge density ρ[r], and the Qenclosed[r] over all the regions to compare their behavior (amplitudes have been adjusted to simplify comparison of curve behaviors):

(*  Input Cell *)
ρ[r_] = If[r < a, 0, If [a ≤  r≤ b,ρo Exp[-α r],If[r>b,0]]]; (* just some help *)

(m)  Make some interpretive remarks about these plots (worry not about repeating comments from above).

(this is a text cell - type in your interpretation)

(n)  Generate a 2D Contour Map of the potential over all three regions.  {Hint: ContourPlot uses Cartesian Coordinates—you can use your general V[r] above and simply replace it with V[r]/.r-> ((((x^2) + (y^2)))^(1/2))}

Using Show, superimpose this ContourPlot with the VectorPlot of EE[r]·(Overscript[r, ⋀]) above.

(*  Input Cell *)

(o) Interpret your combined plot.

(this is a text cell - type in your interpretation)