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;
rho[r_] = If[r < a, 0,If[a < r<b,ρo Exp[ - α r],If[r > b, 0]]]
Plot[rho[r],{r,0,6},AxesLabel->{"r","ρ[r]"}]

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.)  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.

Normally we simply state that ρ is spherically symmetric, so that (OverVector[E]) is radial and the magnitude of (OverVector[E]), Emag, depends only on the distance from the origin O, i.e., on the scalar distance r.  But how do we actually know this?   (Hang on to to your hat!)

First, verify the radial direction of (OverVector[E]).  Although not stated, the above figure is a planar SLICE through the spherical charge distribution, passing though the origin.  With some added and subtracted features, this is the same slice:

We know that if we take a small (differential) volume element, say dτ1' within the shell of charge at the position (OverVector[r1']) [These ' indicate coordinates associated with source points - where there is some charge], the volume element will contain a differential charge dq1 = ρ dτ1'.  This charge produces a differential (OverVector[E]) field, d(OverVector[E1]), at the field point P (at position (OverVector[r]), in this case outside of the charged shell), where d(OverVector[E1]) = k(dq1/((ℜ1)^2)) (OverHat[ℜ])1, where (OverVector[ℜ1]) is the displacement vector between dq1 and P, i.e., (OverVector[ℜ1])  = (OverVector[r]) - (OverVector[r1']) and (OverHat[ℜ])1 is the unit vector in the direction of  (OverVector[ℜ1]).  

For a radially symmetric ρ, we can find a unique volume element dτ2' with the same differential volume and value of ρ, as shown in the diagram.  [Since the source vectors  (OverVector[r1']) and  (OverVector[r2']) have the same lengt, ρ(r1') = ρ(r2').]  Therefore dq1 = dq2.  Forming  d(OverVector[E2]) = k(dq1/((ℜ2)^2)) (OverHat[ℜ2]) and noting from the symmetry we realize that R1 = R2 (the distances), we conclude that  d(OverVector[E1]) and  d(OverVector[E2]) have equal z components and equal BUT OPPOSITE components perpendicular to z.  Therefore their sum, d(OverVector[E]) = d(OverVector[E1]) + d(OverVector[E2]) has only a z component.  As drawn, the z axis is in the (OverHat[r ])direction, therefore is normal to the spherical surface drawn through P (i.e., with center O and radius r).

All the points on the left side of the 3D charge distribution can be paired with corresponding points on the right side of the distribution, so that when the total (OverVector[E]) is found (by integration), the resulting (OverVector[E]) on the surface of a sphere of radius r will be normal to the spherical surface.

IF you are as tired as I am, you will nod enthusiastically at the following:  This process can be applied to P located anywhere (0 ≤ r ≤ ∞, OR saying the same thing: for r < a;  a ≤ r ≤ b; and r > b).    


What about the magnitude of (OverVector[E]) at the point P being a function of r the scalar r — the length of (OverVector[r])) only?

In spherical coordinates, the variables upon which |(OverVector[E])| could depend are {r, θ, ϕ}. If you pick a particular point P and then rotate the charge distribution to ANY other θ and ϕ, the symmetry says the distribution of charge is not changed.  Thus, |(OverVector[E])| cannot depend on θ or ϕ.  Obviously the r dependence remains.  (Think r = 0 vs r ≫ b in our problem.)  So the magnitude of (OverVector[E]) depends on r only.  (This means that  |(OverVector[E])| is constant on any spherical surface centered at O.

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
In order to efficiently and correctly apply Gauss’s Law.  On spherical Gaussian surfaces centered at the origin O, the integrals that arise using Gauss’s Law become greatly simplified.  In fact, the integrals are usually impossible without these simplifications.

(d) Find (OverVector[E])[(OverVector[r])] for r < a  (Probably will not need to use M)

Let’s call this region the “cavity”.  Choose a spherical Gaussian Surface (S1) centered on O, with r < a.

On this surface, (OverVector[E])[(OverVector[r])] = Ecavity[r] (OverHat[r])  and the differential area element is  (OverVector[dS])1 = dS1 (OverHat[r])   (Carefully keep track of what is a vector and what is a scalar in all of the expressions.)    
 
(We use the dummy variable rprime for the radius if it looks like we will be integrating over rprime.)

The charge density inside S1 is ρ[rprime] = 0; the differential volume element for this particular spherical symmetry (where ρ does not depend on θ or ϕ) is dτ1 = 4 π (rprime^2) drprime

Gauss's Law for this surface is:   (∫_S1)(OverVector[E])[(OverVector[r])]·d(OverVector[S1])  = (1/(ϵo)) (∫_(τ1))ρ[rprime]dτ1

OR:  (∫_S1)Ecavity[r]  (OverHat[r]).dS1 (OverHat[r])  =  (1/(ϵo))(∫_(τ1))ρ[rprime] 4 π(rprime^2) drprime

which can be rewritten:    (∫_S1)Ecavity[r]( (OverHat[r])·(OverHat[r]) )dS1  = (1/(ϵo)) (∫_(τ1))(0)4 π (rprime^2) drprime

First, (OverHat[r])· (OverHat[r]) = 1.  Second, on the spherical surface S1, r is a constant (equal to the radius of S1), so Ecavity[r] is constant and can come outside of the lhs integral.  Third, the RHS = 0.   

Therefore:   Ecavity[r] * (∫_S1)dS1  =0.       (∫_S1)dS1 = the area of S1, a sphere of radius r, which is non-zero. (We ignore the possibility that this is not true at r = 0.)   Therefore, for r < a, Ecavity[r] = 0.  Since Ecavity[r] is the magnitude of the vector (OverVector[E])[(OverVector[r])] inside the cavity, (OverVector[E])[(OverVector[r])] = 0 for r < a.

As seen, M is of not much use for this part of the problem.

(* Input Cell -- you really don't need to use it *)

(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:

We apply Gauss’s Law to a spherical Gaussian surface S2 with radius r such that a ≤ r ≤ b (r inside the thick spherical shell)

Since the electric field everywhere in this geometry is constant on and perpendicular to any sphere centered on the origin, the LHS of Gauss’s Law for S2, using the above symmetry arguments is simply:  LHS = Eshell[r] * (∫_S2)dS2  = Eshell[r] * 4 π (r^2) .

(The integral (∫_S2)dS2  is simply the area of S2, a sphere of radius r, i.e., 4 π (r^2).)

The RHS =  (1/(ϵo))(∫_(τ2))ρ[rprime] 4 π(rprime^2) drprime, where τ2 is the volume consisting of the empty spherical cavity (rprime < a) plus the spherical shell contained between a ≤ rprime ≤ r.  

As above, ρ[rprime] = 0 for rprime < a, so the empty cavity does not contribute to the RHS integral.  

Thus, the remaining contribution is:    RHS = (1/(ϵo))(∫_a^r)(Subscript[ρ, o]) Exp[-α rprime] 4 π (rprime^2)drprime, where again, dτ2 = 4 π (rprime^2) drprime

We define both sides in M and then use Solve[ ] to find the function Eshell[r] which represents the magnitude of the non-zero component of E (the radial component).  (Here, Solve[ ] is simply performing the trivial division by 4 π (r^2).)

(* Input Cell *)
ClearAll["`*"]
ρ[rprime_] = ρo Exp[-α rprime];
RHS = (1/(ϵo))(∫_a^r)ρ[rprime]4 π (rprime^2)\[DifferentialD]rprime;
LHS =Eshell[r]* 4 π (r^2) ;
sol = Solve[LHS ==RHS, Eshell[r]];
Eshell[r_]= Eshell[r]/.sol[[1]]//Simplify;
Print["Eshell[r] = ",Eshell[r]]

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

(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

We  apply Gauss’s Law using a spherical Gaussian surface S3 with a radius r > b (r outside of all the charge in  the thick spherical shell).

The LHS of Gauss’s Law for S3, using the above symmetry arguments is :  LHS = Eoutside_shell[r] * (∫_S3)dS3  = Eoutside_shell[r]* 4 π (r^2) 

	 (The integral (∫_S3)dS3  is simply the area of S3,which is a sphere of radius r, i.e., 4 π (r^2).)

The RHS =  (1/(ϵo))(∫_(τ3))ρ[rprime] 4 π(rprime^2) drprime, where τ3 is the volume consisting of the empty spherical cavity (rprime < a) plus the spherical shell contained between a ≤ rprime ≤ b plus the volume of the region for b ≤ rprime ≤ r 

As above, ρ[rprime] = 0 for rprime < a, so the empty cavity does not contribute to the RHS integral.  Likewise, ρ[rprime] = 0 for rprime > b, so the empty outer region does not contribute to the RHS integral.

Thus, the remaining contribution is:    RHS = (1/(ϵo))(∫_a^b)(Subscript[ρ, o]) Exp[-α rprime] 4 π (rprime^2)drprime.  Note the integral represents all the charge in the totality of the thick shell.  

We again define both sides in M and then use Solve[ ] to find the function Eoutsideshell[r] which represents the magnitude of the non-zero component of E (the radial component).  [Again, Solve[ ] is simply performing the trivial division by 4 π (r^2)]

(* Input Cell *)
RHS = (1/(ϵo))(∫_a^b)ρ[rprime]4 π (rprime^2)\[DifferentialD]rprime
LHS =Eoutsideshell[r]* 4 π (r^2) ;
sol = Solve[LHS ==RHS, Eoutsideshell[r]];
Eoutsideshell[r_]= Eoutsideshell[r]/.sol[[1]];
Print["Eoutsideshell[r] = \!\(\*StyleBox["\"",FontFamily->"Source Code Pro",FontWeight->"SemiBold",FontColor->RGBColor[0., 0., 0.]]\), Eoutsideshell[r]]

(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 *)
EE[r_] = If[r < a,0,If[a≤r≤b,Eshell[r],If[r > b,Eoutsideshell[r]]]];
a = 1; b = 2; α = 0.5; ρo = 1; ϵo=1;
plotEE = Plot[EE[r],{r,0,6},AxesLabel->{"r","Eradial[r]"}];
lines1 = Graphics[{Line[{{a,0},{a,1}}],Line[{{b,0},{b,1}}],Text["r = a",{a,0.2},Background->LightRed],Text["r = b",{b,0.2},Background->LightGreen]}];
Show[plotEE,lines1]

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]
EcavityCartesian[x_,y_] = {0,0}
EoutsideshellCartesian[x_,y_]= TransformedField["Polar"->"Cartesian",{Eoutsideshell[r],0},{r,θ}->{x,y}]
EshellCartesian[x_,y_]= TransformedField["Polar"->"Cartesian",{Eshell[r],0},{r,θ}->{x,y}]
EECartesian[x_,y_] = If[Abs[((((x^2) + (y^2)))^(1/2)) ]< a,EcavityCartesian[x,y],If[a<Abs[((((x^2) + (y^2)))^(1/2)) ] < b, EshellCartesian[x,y],If[Abs[((((x^2) + (y^2)))^(1/2)) ]> b,EoutsideshellCartesian[x,y]]]];
a = 1; b = 3; α = 0.5; ρo = 1; ϵo=1;
vecplotEE = VectorPlot[EECartesian[x,y],{x,-6,6},{y,-6,6}];
spheres = Graphics[{Green,Circle[{0,0},a],Red,Circle[{0,0},b],Text["r = a",{a,-0.3},Background->LightRed],Text["r = b",{b,0.3},Background->LightBlue]}];
vecplotEEfinal = Show[vecplotEE,spheres]

(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)
From our symmetry arguments and Gauss’s Law used to derive the E[r] functions for these three regions, we know that:

(1) In all regions, (OverVector[E])[r] = Er[r] (OverHat[r]).  Consistent with this, the vector plot shows all E vectors are radial; at any given r, the LENGTH of (OverVector[E])[r] is the same.

(2) No charge is contained within spherical Gaussian surfaces with r < a, so all E vectors inside this cavity = (OverVector[0]).  No vectors seen in the VectorPlot for this region because they are all of zero length.

(3) Inside the shell (a ≤ r ≤ b), the charge contained within spherical Gaussian surfaces starts at 0 from r = a and reaches a maximum at r = b.  Consequently the magnitude of Er[r] increases with r from 0 at r = a to a maximum at r = b.  The VectorPlot shows this increase.

(4) For r > b, the charge contained within spherical Gaussian surfaces of radius r remains fixed at 

	Qtotal  = (∫_a^b)(ρo Exp[-α rprime])(4 π (rprime^2))drprime = ((4 π ((e^(-a α)) (2+a α (2+a α))+(e^(-b α)) (-2-b α (2+b α))) ρo)/(α^3))

(5)  From past work with E-Fields due to spherically symmetric charge distributions, we know that the field outside the charge distribution will be that of a point charge = Qtotal located at the origin.  It should be a radial field whose magnitude falls as (1/(r^2)) .  The VectorPlot indicates such a field.

(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, 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

We generally take the zero of potential at r = ∞, so that:  

V[r] - V[∞] = – (∫_∞^r)E[rprime]drprime,  where rprime, again, is a dummy variable.  (This is a LINE INTEGRAL).

This expression gives the work per unit charge to move a test charge from infinity to the position r.

Just a reminder—this all holds for radial dependent functions.  CAREFULLY NOTE: the – sign.  If troubled by this, go here:  

http://www.phys.ufl.edu/~acosta/phy2061/lectures/ElectricPotential.pdf

For  (OverVector[E])[r] = Er[r] (OverHat[r]), we have proven that the only component of (OverVector[E])[r] that contributes to the line integral (and the potential) is the radial component; any motion at constant r does not require work.  

FIRST  r > b:  
 
TO FIND Voutsideshell[r] (corresponding to r > b) we form the following integral:  

Voutsideshell[r] =  - (∫_∞^r)Eoutside[rprime] drprime, 

where Eoutside[rprime] is the magnitude of the radial vector field.

(*  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;
Voutsideshell[r_] = - (∫_∞^r)Eoutsideshell[rprime] \[DifferentialD]rprime   (* (∫_∞^r)Eoutsideshell[r] is defined above *)

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;
Vinsideshellonly[r_] =-(∫_b^r)Eshell[rprime]\[DifferentialD]rprime
Vshell[r_] = Voutsideshell[b] + Vinsideshellonly[r];
Vshell[r_] = %//Simplify

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!).

(*  Input cell  *)
Vcavity[r_] = Vshell[a]

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

(*  Input cell *)
V[r_] = If[r < a, Vcavity[r], If [a ≤  r≤ b,Vshell[r],If[r > b, Voutsideshell[r]]]];
a = 1; b = 3; α = 0.5; ρo = 1; ϵo=1;
lines3 = Graphics[{Line[{{a,0},{a,1.4}}],Line[{{b,0},{b,1.4}}],Text["r = a",{a,0.7},Background->LightRed],Text["r = b",{b,0.7},Background->LightGreen]}];
plotv = Plot[ V[r], {r,0,10},PlotRange->{0,1.5},AxesLabel->{"r","V[r]"}];
Show[plotv,lines3]

(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 *)
pp = Plot[{ EE[r],0.2 V[r],0.5ρ[r],0.01 Qenclosed[r]}, {r,0,10},PlotRange->All,PlotLegends->{"Eradial[r]", "V[r]","ρ[r]","Qenclosed[r]"}];
lines4 = Graphics[{Line[{{a,0},{a,0.4}}],Line[{{b,0},{b,0.4}}],Text["r = a",{a,0.1},Background->LightRed],Text["r = b",{b,0.1},Background->LightGreen]}];
Show[pp,lines4]

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

(this is a text cell - type in your interpretation)
Although in some cases repetitive, here are a few more interpretive comments:

•  Outside the shell (r > b), E[r] falls as (1/(r^2)) and V falls as ((1   )/r) corresponding to the (OverVector[E]) and V of a point charge Qtotal located at the origin, 
  
 			where Qtotal = (∫_a^b)ρ[rprime]4 π (rprime^2)drprime
 			
Qtotal is the constant value of Qenclosed for r > b.

Inside the shell (a ≤ r ≤ b), Qenclosed[r] increases from zero (at r = a) to Qtotal (at r = b), due to the integration of the Exp[-α r] factor in ρ[r]; the E field increases in magnitude accordingly.   Because there is no charge contained in the empty cavity, Ecavity = 0.

As r decreases from r = b towards r = a, one must do work to move a test charge because there continues to be charge inside the corresponding Gaussian surface (and thus a non-zero radial electric field).  Thus, V[r] continues to increase with decreasing r.  Once r ≤ a, no more work is needed to move the test charge.

	Ecavity = 0; Vcavity is a constant.

Both E and V are continuous at r = a and r = b.  Such continuity is consistent with the Boundary Conditions of E and V at an uncharged interface.  (The spherical shell with nonzero ρ has no net free charge attached to it’s inner and outer surfaces.)  

The BCs relevant here are:

▫  (Subscript[E, r < a]^⊥) - (Subscript[E, r > a]^⊥) = ((σ[a])/(ϵo))  (where σ[a] = 0) 
▫  (Subscript[E, r < b]^⊥) - (Subscript[E, r > b]^⊥) = ((σ[b])/(ϵo))  (where σ[b] = 0)
▫  (Subscript[V, r < a]) = (Subscript[V, r > a]),  and
▫  (Subscript[V, r < b]) = (Subscript[V, r > b])

(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.

I am rewriting and re-executing vecplotEE below to spiff up the final plot.

(*  Input Cell *)
vecplotEE = VectorPlot[EECartesian[x,y],{x,-6,6},{y,-6,6},VectorColorFunction->Hue];
spheres = Graphics[{Green,Circle[{0,0},a],Red,Circle[{0,0},b],Text["r = a",{a,-0.3},Background->LightRed],Text["r = b",{b,0.3},Background->LightBlue]}];
vecplotEEfinal = Show[vecplotEE,spheres];
conplotV= ContourPlot[V[r]/.r-> ((((x^2) + (y^2)))^(1/2)),{x,-6,6}, {y,-6,6},Contours->30];
Show[conplotV,vecplotEEfinal]

(o) Interpret your combined plot.

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

Reminder,  the 2D plot contour is on a plane that cuts through the z axis of the original problem.

The charge density ( ρ[r] = ρo Exp[-α r] ) is non-zero only between r = a and r = b. 

•  Due to the symmetry of ρ[r],  (OverVector[E]) is radial:    (OverVector[E])[ (OverVector[r]) ] = E[ r ] (OverHat[r])  for all r.  Since no work is done moving along a surface of constant r (E remains constant), the surfaces of constant V are spheres (here, circles).    

Note :  IF you point your mouse at any contour in the plot, the value of V along that contour will pop up.

•  Brain-check:  V is constant (and non-zero) in all of the region where r < a.  Because the increment in V between contours is NOT infinitesimal, you cannot see any contours between r = a and the first plotted contour, which appears at larger r.  BUT V is changing (here decreasing) between r = a and the first plotted contour.

•  Since (OverVector[E])[ (OverVector[r]) ] = –  the gradient of V[r], which is in the direction of steepest descent in the potential, we see that all of the E vectors are perpendicular to the contours.