Spherical Shell of Polarized Dielectric Material

A thick spherical shell has a dielectric filling the region for R1 < r < R2.  The dielectric has a frozen-in Polarization (OverVector[P]) [(OverVector[r])] = Pzero f[r] (OverHat[r]), where Pzero is a constant and r is the (scalar) distance from the center of this spherical system.  

We will examine bound charge distributions and E fields for a few f[r].

(1) Let f[r] = Pzero r

We'll be interested in doing the following:

(a) Finding the E-field everywhere (i.e., all (OverVector[r])).

(b) Making plots in 2D of P and the E field (separately). 2D implies a CUT through the sphere going through the z (polar) axis.

(c)  Commenting on how they relate.

(d) making plots of the radial component of P, the radial component of the E-field, and the magnitude of any relevant charge distributions.

NOTE:  To aid identifying lines of output, I’ve inserted in several places Print Statements of the form:   Print[“u[x] = “,u[x]] which perform like this:

ClearAll["`*"]
u[x_] := uzero Exp[- a x];
Print["u[x] = ", u[x]]

Close enough.

(a) Obviously this problem is Loaded with spherical symmetry.  What can we say about the expected functional forms of any resulting bound charge distributions and electric fields?

<Text Cell.   Type in Answer Here:>

We expect the surface bound charge densities (on the inner and outer surfaces) on the dielectric to be independent of θ and ϕ (since r will either be R1 or R2) they are likely constants.  We expect the bulk bound charge density inside the dielectric to be a function of r only.  

We expect the E field to be of the form {(Subscript[E, r])[r], 0, 0} (i.e., the (OverHat[r]) component is the only non-zero component;  (Subscript[E, r])   depends only on r).

(b) First, let us find the bound charge densities on the surfaces and inside the polarized dielectric.

ClearAll["`*"]
f[r_] := Pzero r;
P = {f[r],0,0};   (* 3D vector in Spherical coordinates; only an r component that depends on the magnitude r only *)
Print["P = ",P,"  REMEMBER:  this is only for R1 < r < R2; otherwise P = 0 "]
ρρρ[r_] =- (1/(r^2))(∂_r)((r^2)P[[1]])//Simplify;    (*  This is the only non-zero component of -Div P in Spherical coordinates; here, we take the derivatives by "hand"*)
Print["ρbound(r)('by hand') = ",ρρρ[r]]
ρbound[r_] := -Div[P,{r,θ, ϕ},"Spherical"]//Simplify    (*Using M's Div function in spherical coordinates  *)
Print["ρbound[r] ('using M') = ",ρbound[r]]
σboundatR1 = P.{-1,0,0}/.r-> R1;
σboundatR2 = P.{1,0,0}/.r-> R2;
Print["σbound on inner surface (r=R1) = ",σboundatR1]
Print["σbound on outer surface (r=R2) = ",σboundatR2]

So we see that from P we can get the bound charge densities; ρbound is non-zero and constant; we get the same result (–3 Pzero) both ways.

The two σbounds have opposite signs and are not equal in magnitude.  We will examine the total bound charge later.

(c)  Find the radial component (only non-zero one) of the E field inside the dielectric.  Use Gauss’s Law.

(A little repetitive):  From the Spherical Symmetry of the problem we can conclude that for R1 < r < R2, the E field will be radial and will depend on r only.

I.e.,   (OverVector[EE])radialinside[r] will be of the form:    (OverVector[EE])radialinside[(OverVector[r])] = { EEradialinside[r] , 0, 0} =  EEradialinside[r] (OverHat[r])    (Note d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r]))

Thus when we write Gauss’s Law:   (∫_S)(OverVector[E])·(OverVector[dS]) = ((q^enclosed)/(ϵo))  for S a spherical surface of radius r (R1 < r < R2), the lhs becomes:

lhs = EEradialinside  (∫_S)dS  = EEradialinside  (4 π (r^2) )      
 (OverVector[EE])radialinside || d(OverVector[S]),  so (OverVector[EE])radialinside · (OverVector[dS]) = EEradialinside dS.)

and EEradialinside is constant for a given r so it comes outside of integral.

(* Apply Gauss's Law *)
lhs = EEradialinside  4 π (r^2);  (* for programming convenience, I've dropped the explicit r dependence; EEradialinside STILL depends on r*)
Print["Bound charge inside dielectric (R1 < rprime < r) = ",(∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r ]  (* NOTE: actual calculation embedded into the Print Statement *)
Print["Bound charge on inner surface (r = R1) = ",σboundatR1 * (4π(R1^2))]
rhs = (1/(ϵo))((∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r + σboundatR1 * (4π(R1^2)))//Simplify;
Print["Right Hand Side of Gauss's Law = ", rhs]
Print["Reminder:  Left hand side of Gauss's Law = ",lhs]
Print["We find EEradial by simply solving for it from lhs = rhs"]
(* Solve for the E field in region where P(and therefore ρ ) are ≠ 0 *)
sol = Solve[lhs == rhs, EEradialinside];
EEradialinside = EEradialinside/.sol[[1]];
Print["EEradialinside =  ",EEradialinside]
Print["Therefore, the vector (OverVector[EE])radial inside the dielectric (in spherical coordinates) = ", {EEradialinside,0,0}]

Note that (OverVector[EE])radialinside[r] = –(((OverVector[P])[r])/(ϵo))        REASON: By definition,  (OverVector[D]) = ϵo (OverVector[E]) + (OverVector[P]).  

For any spherical Gaussian Surface of radius r,

			 (∫_S)(OverVector[D])·(OverVector[dS]) = (Subscript[q, free]^enclosed)  = 0;  Therefore (OverVector[D]) = 0
Then, from definition of (OverVector[D]) we see that (Overscript[E,  ⇀]) =  –(((OverVector[P])[r])/(ϵo)) ( which is (OverVector[EE])radialinside[r] ).

(d)  Find the total bound charge for the entire region 0 ≤ r ≤ R2.

From above, recall that:

"ρbound[r] = "-3 Pzero
"σbound on inner surface (r=R1) = "-Pzero R1           
"σbound on outer surface (r=R2) = "Pzero R2

[ALL THREE OF THESE QUANTITIES ARE CONSTANTS FOR THE GIVEN constants:  Pzero, R1, and R2)

We can therefore find TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound 4 π (r^2)dr

Print["Bound (surface) Charge on Inner Surface (r = R1) = ",σboundatR1*4π(R1^2)]
Print["Bound (surface) Charge on Outer Surface (r = R2) = ",σboundatR2*4π(R2^2)]
Print["Total Bound (bulk charge) inside dielectric  (R1 < r < R2) = ",(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify]
TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify;
Print["Total Bound charge = ",TotalQ]

(e)  Comment on your result for Total Bound Charge:

The result that the total bound charge = 0 is not too surprising.  If we imagine applying an E field to a neutral, initially unpolarized dielectric, the E field would cause molecules (or atoms) to undergo internal charge separation and/or permanent dipoles (certain molecules) to partially align in the field.  NO charge was created or released in the process.  The surface bound charge can be thought of as being due to the tips and tails of the dipoles at the surfaces (tips +; tails –).  The bulk bound charge arises because the dipole density is not constant but has a non-zero divergence.

Mathematically, the Divergence Theorem comes to our rescue:  

If we apply it to P:   (∫_S)P·dS (OverHat[n])  = (∫_S)(P·(OverHat[n]) )dS = (∫_τ)Div[P]dτ 

where S is the surface and τ is the volume of the polarized dielectric.  But P·(OverHat[n]) is the surface bound charge density and Div[P] is equal to – (volume charge density)

So we can write (∫_S)(surface bound charge density )dS + (∫_τ)(volume charge density)dτ = 0

The lhs is the sum of the total surface bound charge and the total volume bound charge which = THE TOTAL BOUND CHARGE

Therefore  The total bound charge = 0      QED

(f) Find the electric field outside of the dielectric, i.e., r > R2.

Spherical symmetry implies that for r > R2, (OverVector[E])outside will be of the form {Eoutside[r], 0, 0} = Eoutside[r] (OverHat[r]) .
      
d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r])

Applying Gauss’s Law on a concentric sphere of radius r > R2 (utilizing  (q^(enclosed )) = 0)

(∫_S)(OverVector[E])·d(OverVector[S]) = ((q^enclosed)/(ϵo))        ⇒  EEoutside = 0.

(g) Find the electric field inside the spherical cavity;  i.e., r  <  R1.

Applying Gauss’s Law on a concentric sphere of radius r < R1 and utilizing  (q^(enclosed )) = 0, we find that Einsidecavity = 0.

Now let’s do some plotting.  First for vector plots, we can cut the spheres with a plane (say at ϕ = 0).  Since there is no ϕ dependence, all of the resulting 2 D plots will look the same.  You will want to transform to Cartesian Coordinates to keep M a happy plotter.

(h) Create a 2D VectorPlot of the polarization (OverVector[P]).   [2D is sufficient because of the spherical symmetry — any plane cut through the polar (z) axis exhibits the same vector plot.]

P is non zero only for R1 < r < R2.  We can use If statements (or the Piecewise function) to generate a function that covers all space:

Clear[Pzero]
Pcart = TransformedField["Polar"-> "Cartesian",{f[r],0},{r, θ} -> {x,z}]
Pzero = 1; a = .3; R1 = .6; R2 = 2.5; ϵo = 1;
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
vecP = VectorPlot[If[((((x^2) + (z^2)))^(1/2))<R1,{0,0},If[R1 < ((((x^2) + (z^2)))^(1/2))<R2,Pcart,If[((((x^2) + (z^2)))^(1/2)) > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["2D Vector Plot of the Polarization Vector, P"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecP,sphereboundaries]

P is non-zero only inside the dielectric (where there are molecules that can be polarized).   

∴  Inside the hollow region (r < R1),  P = 0.

Inside the dielectric (R1 < r < R2):  With P = {Pzero r, 0, 0}, we see radial P vectors increasing in length as r increases.  

For r > R2 (outside the sphere), P = 0.

(i)  Create a 2D VectorPlot of the Electric Field vector.

From above we can say that EE = 0 for r < R1,  EE = EEradialinside =  -((Pzero r)/(ϵo)) for R1 < r < R2,  and EE = 0 for r > R2.

Clear[Pzero,a,R1,R2,ϵo]
EEradialinside
EEcartinside = TransformedField["Polar"-> "Cartesian",{EEradialinside,0},{r, θ} -> {x,z}]
Pzero = 1; a = .8; R1 = .6; R2 = 2.5; ϵo = 1;
r = ((((x^2) + (z^2)))^(1/2))
vecE = VectorPlot[ If[r < R1, {0,0}, If[R1 < r < R2, EEcartinside, If[r > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["2D Vector Plot of the E-Field Vector, E"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecE,sphereboundaries]

From our determination of Inside the hollow region, no E field; with P = {Pzero r, 0, 0}, we got ρbound = –Div[P] = –3;  This means that |E| should increase with r and point towards the origin.  E outside the sphere = 0.  

Check –Div P:

Clear[r, Pzero]
-Div[P,{r,θ, ϕ},"Spherical"];
Print["ρbound using -Div P = ", %]

(j)  For clarity, plot the magnitudes of P, the E-Field, and ρbound vs. r.   The vectors are one D (radial component only) so simply plot those components vs. r; ρbound is a scalar function of r only, i.e.,  ρbound[r].

Clear[r]
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
Pzero = 1; a = 0.4; R1 = .6; R2 = 2.5; ϵo = 1;
PPradialplot = Plot[PPradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of P, the Polarization Vector(radial component), vs. r"],12,Red,Background->Lighter[Green]]]

P is given as radial; P is non zero only inside the shell (R1 < r < R2) and increases linearly from Pzero R1 to Pzero R2 inside the shell.

rho = ρbound[r]
ρρ = If[r<R1,0,If[R1 < r<R2, rho,If[r> R2,0]]];
ρplot = Plot[ρρ,{r,0,2R2},PlotRange-> All,PlotLabel->Style[Framed["Bound Volume Charge Density, ρbound[r]"],12,Red,Background->Lighter[Green]]]

ρbound is zero outside of the shell (P is zero outside, so outside:   – Div P is zero).  P ~ r means ρbound is constant inside the shell (= – 3 Pzero)

EEEEradial = If[r<R1,0,If[R1 < r<R2,
EEradialinside,If[r> R2,0]]];
EEradialplot = Plot[EEEEradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of the E-Field Vector(radial component)"],12,Red,Background->Lighter[Green]]]

The E field is zero inside the hollow cavity of the shell (r < R1) and outside the shell (r > R2) for reasons given above.  Inside the shell (R1 < r < R2), we showed that EEradial = -((Pzero r)/(ϵo)).  So, inside the shell, (OverVector[E]) points toward the origin.

Put them all together

Plot[{ PPradial+0.06,.8ρρ-0.06,2EEEEradial},{r,0, 2 R2},PlotRange-> All, PlotLegends->{"Polarization (radial)","Bound Charge Density ((Subscript[ρ, bound]))", "Electric Field (radial)"},PlotLabel->Style[Framed["Magnitudes of P, ρbound, and the E-Field Vector(radial component)"],11,Red,Background->Lighter[Green]] ](*  I've displaced P and E slightly so you can see the curves  *)
Print["Bound Surface Charge Density on inner surface = ",σboundatR1]
Print["Bound Surface Charge Density on outer surface = ",σboundatR2]

Note the signs of the bound surface charge densities, negative at r = R1 and positive at r = R2.

You should be able to reason why these signs arise from this VectorPlot of (OverVector[P]).  (Remember:  σbound = (OverVector[P]) · (OverVector[n]))

2.  f[r_] = Pzero (r^2)

Basically all we have to do is take the code above and simply replace f[r]:   f[r_] = Pzero (r^2)   

Most of what you see below is copied from above.

(a) Again, we note that this problem is Loaded with spherical symmetry.  What can we say about the expected functional forms of any resulting bound charge distributions and electric fields?

<Text Cell.   Type in Answer Here:>

Same answer as before:  We expect the surface bound charge densities (on the inner and outer surfaces) on the dielectric to be independent of θ and ϕ (since r will either be R1 or R2) they are likely constants.  We expect the bulk bound charge density inside the dielectric to be a function of r only.  

We expect the E field to be of the form {(Subscript[E, r])[r], 0, 0} (i.e., the (OverHat[r]) component is the only non-zero component;  (Subscript[E, r])   depends only on r).

(b) First, let us find the bound charge densities on the surfaces and inside the polarized dielectric.

ClearAll["`*"]
f[r_] := Pzero (r^2);
P = {f[r],0,0};   (* 3D vector in Spherical coordinates; only an r component that depends on the magnitude r only *)
Print["P = ",P,"  REMEMBER:  this is only for R1 < r < R2; otherwise P = 0 "]
ρρρ[r_] =- (1/(r^2))(∂_r)((r^2)P[[1]])//Simplify;    (*  This is the only non-zero component of -Div P in Spherical coordinates; here, we take the derivatives by "hand"*)
Print["ρbound(r)('by hand') = ",ρρρ[r]]
ρbound[r_] := -Div[P,{r,θ, ϕ},"Spherical"]//Simplify    (*Using M's Div function in spherical coordinates  *)
Print["ρbound[r] ('using M') = ",ρbound[r]]
σboundatR1 = P.{-1,0,0}/.r-> R1;
σboundatR2 = P.{1,0,0}/.r-> R2;
Print["σbound on inner surface (r=R1) = ",σboundatR1]
Print["σbound on outer surface (r=R2) = ",σboundatR2]

So we see that from P we can get the bound charge densities; ρbound is non-zero = –4 Pzero r.  We get the same result (–4 Pzero r) both ways.

The two σbounds have opposite signs and are not equal in magnitude.  We will examine the total bound charge later.

(c)  Find the radial component (only non-zero one) of the E field inside the dielectric.  Use Gauss’s Law.

(A little repetitive):  From the Spherical Symmetry of the problem we can conclude that for R1 < r < R2, the E field will be radial and will depend on r only.

I.e.,   (OverVector[EE])radialinside[r] will be of the form:    (OverVector[EE])radialinside[(OverVector[r])] = { EEradialinside[r] , 0, 0} =  EEradialinside[r] (OverHat[r])    (Note d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r]))

Thus when we write Gauss’s Law:   (∫_S)(OverVector[E])·d(OverVector[S]) = ((q^enclosed)/(ϵo))  for a spherical surface S of radius r (R1 < r < R2), the lhs becomes:

lhs = EEradialinside  (∫_S)dS  = EEradialinside  (4 π (r^2) ) .     
[ (OverVector[EE])radialinside || d(OverVector[S])  (so (OverVector[EE])radialinside · d(OverVector[S]) = EEradialinside dS), where EEradialinside is constant for a given r so it comes outside of integral.]

(* Apply Gauss's Law *)
lhs = EEradialinside  4 π (r^2);  (* for programming convenience, I've dropped the explicit r dependence; EEradialinside STILL depends on r*)
Print["Bound charge inside dielectric (R1 < rprime < r) = ",(∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r ]  
(* NOTE: actual calculation embedded into the Print Statement *)
Print["Bound charge on inner surface (r = R1) = ",σboundatR1 * (4π(R1^2))]
rhs = (1/(ϵo))((∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r + σboundatR1 * (4π(R1^2)))//Simplify;
Print["Right Hand Side of Gauss's Law = ", rhs]
Print["Reminder:  Left hand side of Gauss's Law = ",lhs]
Print["We find EEradial by simply solving for it from lhs = rhs"]
(* Solve for the E field in region where P(and therefore ρ ) are ≠ 0 *)
sol = Solve[lhs == rhs, EEradialinside];
EEradialinside = EEradialinside/.sol[[1]];
Print["EEradialinside =  ",EEradialinside]
Print["Therefore, the vector (OverVector[EE])radial inside the dielectric (in spherical coordinates) = ", {EEradialinside,0,0}]

Note that (OverVector[EE])radialinside[r] = –(((OverVector[P])[r])/(ϵo))        REASON:  (OverVector[D]) = ϵo (OverVector[E]) + (OverVector[P]);  for any spherical Gaussian Surface of radius r,
			 (∫_S)(OverVector[D])·(OverVector[dS]) = (Subscript[q, free]^enclosed)  = 0;  Therefore (OverVector[D]) = 0
Then, from definition of (OverVector[D]) we see that (Overscript[E,  ⇀]) =  –(((OverVector[P])[r])/(ϵo)) (which is (OverVector[EE])radialinside[r] ).

(d)  Find the total bound charge for the entire region 0 ≤ r ≤ R2.

From above, recall that:

ρbound = –4 Pzero r
σbound on inner surface (r = R1) =  –Pzero (R1^2)
σbound on outer surface (r = R2) =  –Pzero (R2^2)

We can therefore find TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound 4 π (r^2)dr

Print["Bound (surface) Charge on Inner Surface (r = R1) = ",σboundatR1*4π(R1^2)]
Print["Bound (surface) Charge on Outer Surface (r = R2) = ",σboundatR2*4π(R2^2)]
Print["Total Bound (bulk charge) inside dielectric  (R1 < r < R2) = ",(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify]
TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify;
Print["Total Bound charge = ",TotalQ]

(e)  Comment on your result for Total Bound Charge:

Again, the result that the total bound charge = 0 is not too surprising.  If we imagine applying an E field to a neutral, initially unpolarized dielectric, the E field would cause molecules (or atoms) to undergo internal charge separation and/or permanent dipoles (certain molecules) to partially align in the field.  NO charge was created or released in the process.  The surface bound charge can be thought of as being due to the tips and tails of the dipoles at the surfaces (tips +; tails –).  The bulk bound charge arises because the dipole density is not constant but has a non-zero divergence.

Mathematically, the Divergence Theorem comes to our rescue:  

If we apply it to P:   (∫_S)P·dS (OverHat[n])  = (∫_S)(P· (OverHat[n]) )dS = (∫_τ)Div[P]dτ 

where S is the surface and τ is the volume of the polarized dielectric.  But P. (OverHat[n]) is the surface bound charge density and Div[P] is equal to –(volume charge density)

So we can write (∫_S)(surface bound charge density )dS + (∫_τ)(volume charge density)dτ = 0

The lhs is the sum of the total surface bound charge and the total volume bound charge which = THE TOTAL BOUND CHARGE

Therefore  The total bound charge = 0      QED

(f) Find the electric field outside of the dielectric, i.e., r > R2.

Spherical symmetry implies that for r > R2, (OverVector[E])outside will be of the form {Eoutside[r], 0, 0} = Eoutside[r] (OverHat[r]) .
d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r]).

Applying Gauss’s Law on a concentric sphere of radius r > R2 (utilizing  (q^(enclosed )) = 0)

(∫_S)(OverVector[E])·(OverVector[dS]) = ((q^enclosed)/(ϵo))        ⇒  EEoutside = 0.

(g) Find the electric field inside the spherical cavity;  i.e., r < R1.

Applying Gauss’s Law on a concentric sphere of radius r < R1 and utilizing  (q^(enclosed )) = 0, we find that Einsidecavity = 0.

Now let’s do some plotting.  First for vector plots, we can cut the spheres with a plane (say at ϕ = 0).  Since there is no ϕ dependence, all of the resulting 2 D plots will look the same.  You will want to transform to Cartesian Coordinates to keep M a happy plotter.

(h) Create a 2D VectorPlot of the polarization (OverVector[P]) in 2D.

P is non zero only for R1 < r < R2.  We can use If statements (or the Piecewise function) to generate a function that covers all space:

Clear[Pzero]
Pcart = TransformedField["Polar"-> "Cartesian",{f[r],0},{r, θ} -> {x,z}]
Pzero = 1; a = .3; R1 = .6; R2 = 2.5; ϵo = 1;
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
vecP = VectorPlot[If[((((x^2) + (z^2)))^(1/2))<R1,{0,0},If[R1 < ((((x^2) + (z^2)))^(1/2))<R2,Pcart,If[((((x^2) + (z^2)))^(1/2)) > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the Polarization Vector, P"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecP,sphereboundaries]

P is non-zero only inside the dielectric (where there are molecules that can be polarized).   

∴  Inside the hollow region (r < R1), P = 0.

Inside the dielectric (R1 < r < R2):  With P = {Pzero r, 0, 0}, we see radial P vectors increasing in length as r increases.  

For r > R2 (outside the sphere), P = 0.

(i)  Create a 2D VectorPlot of the Electric Field vector.

From above we can say that EE = 0 for r < R1,  EE = EEradialinside =  -((Pzero r)/(ϵo)) for R1 < r < R2,  and EE = 0 for r > R2.

Clear[Pzero,a,R1,R2,ϵo]
EEradialinside
EEcartinside = TransformedField["Polar"-> "Cartesian",{EEradialinside,0},{r, θ} -> {x,z}]
Pzero = 1; a = .8; R1 = .6; R2 = 2.5; ϵo = 1;
r = ((((x^2) + (z^2)))^(1/2))
vecE = VectorPlot[ If[r < R1, {0,0}, If[R1 < r < R2, EEcartinside, If[r > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the E-Field Vector, E"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecE,sphereboundaries]

From our determination of Inside the hollow region, no E field; with P = {Pzero r, 0, 0}, we got ρbound = –Div[P] = –4 Pzero r;  This means that |E| should increase with r and point towards the origin.  E outside the sphere = 0.   

Check –Div P:

Clear[r, Pzero]
-Div[P,{r,θ, ϕ},"Spherical"];
Print["ρbound using -Div P = ", %]

Plots:   magnitudes vs. r of PPradial, EEradial, and ρbound[r]

Clear[r]
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
Pzero = 1; a = .4; R1 = .6; R2 = 2.5; ϵo = 1;
PPradialplot = Plot[PPradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of P, the Polarization Vector(radial component), vs. r"],12,Red,Background->Lighter[Green]]]

rho = ρbound[r]
ρρ = If[r<R1,0,If[R1 < r<R2, rho,If[r> R2,0]]];
ρplot = Plot[ρρ,{r,0,2R2},PlotRange-> All,PlotLabel->Style[Framed["Bound Volume Charge Density, ρbound[r]"],12,Red,Background->Lighter[Green]]]

EEEEradial = If[r<R1,0,If[R1 < r<R2,
EEradialinside,If[r> R2,0]]];
EEradialplot = Plot[EEEEradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of the E-Field Vector (radial component) vs. r "],12,Red,Background->Lighter[Green]]]

Put them all together

Plot[{ PPradial+0.03,.8ρρ-0.03,2EEEEradial},{r,0, 2 R2},PlotRange-> All, PlotLegends->{"Polarization (radial)","Bound Charge Density ((Subscript[ρ, bound]))", "Electric Field (radial)"},PlotLabel->Style[Framed["Magnitudes of P, ρbound, and the E-Field Vector(radial component)"],11,Red,Background->Lighter[Green]]]
Print["Surface Charge Density on inner surface = ",σboundatR1]
Print["Surface Charge Density on outer surface = ",σboundatR2]

3.  f[r_] = Pzero Exp[–a r]

(b) First, let us find the bound charge densities on the surfaces and inside the polarized dielectric.

ClearAll["`*"]
f[r_] := Pzero Exp[- a r];
P = {f[r],0,0};   (* 3D vector in Spherical coordinates; only an r component that depends on the magnitude r only *)
Print["P = ",P,"  REMEMBER:  this is only for R1 < r < R2; otherwise P = 0 "]
ρρρ[r_] =- (1/(r^2))(∂_r)((r^2)P[[1]])//Simplify;    (*  This is the only non-zero component of -Div P in Spherical coordinates; here, we take the derivatives by "hand"*)
Print["ρbound(r)('by hand') = ",ρρρ[r]]
ρbound[r_] := -Div[P,{r,θ, ϕ},"Spherical"]//Simplify    (*Using M's Div function in spherical coordinates  *)
Print["ρbound[r] ('using M') = ",ρbound[r]]
σboundatR1 = P.{-1,0,0}/.r-> R1;
σboundatR2 = P.{1,0,0}/.r-> R2;
Print["σbound on inner surface (r=R1) = ",σboundatR1]
Print["σbound on outer surface (r=R2) = ",σboundatR2]

So we see that from P we can get the bound charge densities; ρbound is non-zero and constant; we get the same result ( (((e^(-a r)) Pzero (-2+a r))/r)) both ways.

The two σbounds have opposite signs and are not equal in magnitude.  We will examine the total bound charge later.

(c)  Find the radial component (only non-zero one) of the E field inside the dielectric.  Use Gauss’s Law.

(A little repetitive):  From the Spherical Symmetry of the problem we can conclude that for R1 < r < R2, the E field will be radial and will depend on r only.

I.e.,   (OverVector[EE])radialinside[r] will be of the form:    (OverVector[EE])radialinside[(OverVector[r])] = { EEradialinside[r] , 0, 0} =  EEradialinside[r] (OverHat[r]).    
(Note d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r]).)

Thus when we write Gauss’s Law:   (∫_S)(OverVector[E])·(OverVector[dS]) = ((q^enclosed)/(ϵo))  for S a spherical surface of radius r (R1 < r < R2), the lhs becomes:

lhs = EEradialinside  (∫_S)dS  = EEradialinside  (4 π (r^2) )      
[ (OverVector[EE])radialinside || d(OverVector[S])  (so (OverVector[EE])radialinside · (OverVector[dS]) = EEradialinside dS.]

and EEradialinside is constant for a given r so it comes outside of integral].

(* Apply Gauss's Law *)
lhs = EEradialinside  4 π (r^2);  (* for programming convenience, I've dropped the explicit r dependence; EEradialinside STILL depends on r*)
Print["Bound charge inside dielectric (R1 < rprime < r) = ",(∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r ]  (* NOTE: actual calculation embedded into the Print Statement *)
Print["Bound charge on inner surface (r = R1) = ",σboundatR1 * (4π(R1^2))]
rhs = (1/(ϵo))((∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r + σboundatR1 * (4π(R1^2)))//Simplify;
Print["Right Hand Side of Gauss's Law = ", rhs]
Print["Reminder:  Left hand side of Gauss's Law = ",lhs]
Print["We find EEradial by simply solving for it from lhs = rhs"]
(* Solve for the E field in region where P(and therefore ρ ) are ≠ 0 *)
sol = Solve[lhs == rhs, EEradialinside];
EEradialinside = EEradialinside/.sol[[1]];
Print["EEradialinside =  ",EEradialinside]
Print["Therefore, the vector (OverVector[EE])radial inside the dielectric (in spherical coordinates) = ", {EEradialinside,0,0}]

Note that (OverVector[EE])radialinside[r] = – (((OverVector[P])[r])/(ϵo))        REASON:  (OverVector[D]) = ϵo (OverVector[E]) + (OverVector[P]);  for any spherical Gaussian Surface of radius r,
			 (∫_S)(OverVector[D])·d(OverVector[S]) = (Subscript[q, free]^enclosed)  = 0;  Therefore (OverVector[D]) = 0.
Then, from definition of (OverVector[D]) we see that (Overscript[E,  ⇀]) =  –(((OverVector[P])[r])/(ϵo)) ( which is (OverVector[EE])radialinside[r] ).

(d)  Find the total bound charge for the entire region 0 ≤ r ≤ R2.

From above, recall that:

ρbound[r] = ((Exp[-a r] Pzero (-2 + a r))/r)
σbound on inner surface (r = R1)  = – Exp[– a R1] Pzero
σbound on outer surface (r = R2)  = + Exp[– a R2] Pzero

We can therefore find TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound 4 π (r^2)dr

Print["Bound (surface) Charge on Inner Surface (r = R1) = ",σboundatR1*4π(R1^2)]
Print["Bound (surface) Charge on Outer Surface (r = R2) = ",σboundatR2*4π(R2^2)]
Print["Total Bound (bulk charge) inside dielectric  (R1 < r < R2) = ",(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify]
TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify;
Print["Total Bound charge = ",TotalQ]

(e)  Comment on your result for Total Bound Charge:

One more time:  The result that the total bound charge = 0 is not too surprising.  If we imagine applying an E field to a neutral, initially unpolarized dielectric, the E field would cause molecules (or atoms) to undergo internal charge separation and/or permanent dipoles (certain molecules) to partially align in the field.  NO charge was created or released in the process.  The surface bound charge can be thought of as being due to the tips and tails of the dipoles at the surfaces (tips +; tails –).  The bulk bound charge arises because the dipole density is not constant but has a non-zero divergence.

Mathematically, the Divergence Theorem comes to our rescue:  

If we apply it to P:   (∫_S)P·dS (OverHat[n])  = (∫_S)(P· (OverHat[n]) )dS = (∫_τ)Div[P]dτ 

where S is the surface and τ is the volume of the polarized dielectric.  But P·(OverHat[n]) is the surface bound charge density and Div[P] is equal to – (volume charge density).

So we can write (∫_S)(surface bound charge density )dS + (∫_τ)(volume charge density)dτ = 0

The lhs is the sum of the total surface bound charge and the total volume bound charge which = THE TOTAL BOUND CHARGE.

Therefore  The total bound charge = 0      QED

(f) Find the electric field outside of the dielectric, i.e., r > R2.

Spherical symmetry implies that for r > R2, (OverVector[E])outside will be of the form {Eoutside[r], 0, 0} = Eoutside[r] (OverHat[r]) ;       d(OverVector[S]) = {dS, 0, 0} = dS (OverHat[r]).

Applying Gauss’s Law on a concentric sphere of radius r > R2 (utilizing  (q^(enclosed )) = 0)

(∫_S)(OverVector[E])·d(OverVector[S]) = ((q^enclosed)/(ϵo))        ⇒  EEoutside = 0.

(g) Find the electric field inside the spherical cavity;  i.e., r < R1.

Applying Gauss’s Law on a concentric sphere of radius r < R1 and utilizing  (q^(enclosed )) = 0, we find that Einsidecavity = 0.

Now let’s do some plotting.  First for vector plots, we can cut the spheres with a plane (say at ϕ = 0).  Since there is no ϕ dependence, all of the resulting 2 D plots will look the same.  You will want to transform to Cartesian Coordinates to keep M a happy plotter.

(h) Create a 2D VectorPlot of the polarization (OverVector[P]) in 2D .

P is non zero only for R1 < r < R2.  We can use If statements (or the Piecewise function) to generate a function that covers all space:

Clear[Pzero]
Pcart = TransformedField["Polar"-> "Cartesian",{f[r],0},{r, θ} -> {x,z}]
Pzero = 1; a = .3; R1 = .6; R2 = 2.5; ϵo = 1;
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
vecP = VectorPlot[If[((((x^2) + (z^2)))^(1/2))<R1,{0,0},If[R1 < ((((x^2) + (z^2)))^(1/2))<R2,Pcart,If[((((x^2) + (z^2)))^(1/2)) > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the Polarization Vector, P"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecP,sphereboundaries]

P is non-zero only inside the dielectric (where there are molecules that can be polarized).   

∴  Inside the hollow region (r < R1), P = 0.

Inside the dielectric (R1 < r < R2),  where P = {Pzero Exp[– a r], 0, 0}, we see radial P vectors decrease in length as r increases.  

For r > R2 (outside the sphere), P = 0.

(i)  Create a 2D VectorPlot of the Electric Field vector.

From above we can say that EE = 0 for r < R1,  EE = EEradialinside =  -(((e^(-a r)) Pzero)/(ϵo))for R1 < r < R2,  and EE = 0 for r > R2.

Clear[Pzero,a,R1,R2,ϵo]
EEradialinside
EEcartinside = TransformedField["Polar"-> "Cartesian",{EEradialinside,0},{r, θ} -> {x,z}]
Pzero = 1; a = .8; R1 = .6; R2 = 2.5; ϵo = 1;
r = ((((x^2) + (z^2)))^(1/2))
vecE = VectorPlot[ If[r < R1, {0,0}, If[R1 < r < R2, EEcartinside, If[r > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the E-Field Vector, E"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecE,sphereboundaries]

Again, inside the hollow region, no E field; Within the polarized material, P = {Pzero Exp[– a r], 0, 0}, we got ρbound = –Div[P] = (((e^(-a r)) Pzero (-2+a r))/r);  This means that |E| should increase with r and point towards the origin.  E outside the sphere = 0.   

Check –Div P:

Clear[r, Pzero,a]
-Div[P,{r,θ, ϕ},"Spherical"];
Print["ρbound using -Div P = ", %]

Plots:   magnitudes vs. r of PPradial, EEradial, and ρbound[r]

Clear[r]
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
Pzero = 1; a = .4; R1 = .6; R2 = 2.5; ϵo = 1;
PPradialplot = Plot[PPradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of P, the Polarization Vector(radial component), vs. r"],12,Red,Background->Lighter[Green]]]

rho = ρbound[r];
ρρ = If[r<R1,0,If[R1 < r<R2, rho,If[r> R2,0]]];
ρplot = Plot[ρρ,{r,0,2R2},PlotRange-> All,PlotLabel->Style[Framed["Bound Volume Charge Density, ρbound[r]"],12,Red,Background->Lighter[Green]]]

EEEEradial = If[r<R1,0,If[R1 < r<R2,
EEradialinside,If[r> R2,0]]];
EEradialplot = Plot[EEEEradial,{r,0,2R2},PlotLabel->Style[Framed["Magnitude of the E-Field Vector(radial component)"],12,Red,Background->Lighter[Green]]]

Put them all together

Plot[{ PPradial+0.03,.8ρρ-0.03,2EEEEradial},{r,0, 2 R2},PlotRange-> All, PlotLegends->{"Polarization (radial)","Bound Charge Density ((Subscript[ρ, bound]))", "Electric Field (radial)"},PlotLabel->Style[Framed["Magnitudes of P, ρbound, and the E-Field Vector(radial component)"],11,Red,Background->Lighter[Green]]]
Print["Surface Charge Density on inner surface = ",σboundatR1]
Print["Surface Charge Density on outer surface = ",σboundatR2]

Here I Copy and merge the blocks of code — for laughs, I put in one more variation of f[r]; if you are highly motivated, check out the resulting equations and plots, SIGNS in particular.

ClearAll["`*"]
f[r_] := Pzero Sin[π r];
P = {f[r],0,0};   (* 3D vector in Spherical coordinates; only an r component that depends on the magnitude r only *)
Print["P = ",P,"  REMEMBER:  this is only for R1 < r < R2; otherwise P = 0 "]
ρρρ[r_] =- (1/(r^2))(∂_r)((r^2)P[[1]])//Simplify;    (*  This is the only non-zero component of -Div P in Spherical coordinates; here, we take the derivatives by "hand"*)
Print["ρbound(r)('by hand') = ",ρρρ[r]]
ρbound[r_] := -Div[P,{r,θ, ϕ},"Spherical"]//Simplify    (*Using M's Div function in spherical coordinates  *)
Print["ρbound[r] = ",ρbound[r]]
σboundatR1 = P.{-1,0,0}/.r-> R1;
σboundatR2 = P.{1,0,0}/.r-> R2;
Print["σbound on inner surface (r=R1) = ",σboundatR1]
Print["σbound on outer surface (r=R2) = ",σboundatR2]
(* Apply Gauss's Law *)
lhs = EEradialinside  4 π (r^2);  (* for programming convenience, I've dropped the explicit r dependence; EEradialinside STILL depends on r*)
Print["Bound charge inside dielectric (R1 < rprime < r) = ",(∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r ]  (* NOTE: actual calculation embedded into the Print Statement *)
Print["Bound charge on inner surface (r = R1) = ",σboundatR1 * (4π(R1^2))]
rhs = (1/(ϵo))((∫_R1^r)ρbound[r] 4 π (r^2) \[DifferentialD]r + σboundatR1 * (4π(R1^2)))//Simplify;
Print["Right Hand Side of Gauss's Law = ", rhs]
Print["Reminder:  Left hand side of Gauss's Law = ",lhs]
Print["We find EEradial by simply solving for it from lhs = rhs"]
(* Solve for the E field in region where P(and therefore ρ ) are ≠ 0 *)
sol = Solve[lhs == rhs, EEradialinside];
EEradialinside = EEradialinside/.sol[[1]];
Print["EEradialinside =  ",EEradialinside]
Print["Therefore, the vector (OverVector[EE])radial inside the dielectric (in spherical coordinates) = ", {EEradialinside,0,0}]
Print["Bound (surface) Charge on Inner Surface (r = R1) = ",σboundatR1*4π(R1^2)]
Print["Bound (surface) Charge on Outer Surface (r = R2) = ",σboundatR2*4π(R2^2)]
Print["Total Bound (bulk charge) inside dielectric  (R1 < r < R2) = ",(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify]
TotalQ = σboundatR1*4π(R1^2)+σboundatR2*4π(R2^2)+(∫_R1^R2)ρbound[r] 4 π (r^2)\[DifferentialD]r//Simplify;
Print["Total Bound charge = ",TotalQ]
Clear[Pzero]
Pcart = TransformedField["Polar"-> "Cartesian",{f[r],0},{r, θ} -> {x,z}]
Pzero = 1; a = .3; R1 = .6; R2 = 2.5; ϵo = 1;
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
vecP = VectorPlot[If[((((x^2) + (z^2)))^(1/2))<R1,{0,0},If[R1 < ((((x^2) + (z^2)))^(1/2))<R2,Pcart,If[((((x^2) + (z^2)))^(1/2)) > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the Polarization Vector, P"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecP,sphereboundaries]
Clear[Pzero,a,R1,R2,ϵo]
EEradialinside
EEcartinside = TransformedField["Polar"-> "Cartesian",{EEradialinside,0},{r, θ} -> {x,z}]
Pzero = 1; a = .8; R1 = .6; R2 = 2.5; ϵo = 1;
r = ((((x^2) + (z^2)))^(1/2));
vecE = VectorPlot[ If[r < R1, {0,0}, If[R1 < r < R2, EEcartinside, If[r > R2,{0,0}]]],{x,-4,4},{z,-4,4},PlotLabel->Style[Framed["Plot of the E-Field Vector, E"],16,Blue,Background->Lighter[Yellow]]];
sphereboundaries = Graphics[{Red,Thick,Circle[{0,0},R1],Purple,Thick,Circle[{0,0},R2]}];
Show[vecE,sphereboundaries]

Clear[r]
PPradial = If[r<R1,0,If[R1 < r<R2,f[r],If[r> R2,0]]];
Pzero = 1; a = .4; R1 = .6; R2 = 2.5; ϵo = 1;
PPradialplot = Plot[PPradial,{r,0,2R2}]
rho = ρbound[r]
ρρ = If[r<R1,0,If[R1 < r<R2, rho,If[r> R2,0]]];
ρplot = Plot[ρρ,{r,0,2R2},PlotRange-> All]
EEEEradial = If[r<R1,0,If[R1 < r<R2,
EEradialinside,If[r> R2,0]]];
EEradialplot = Plot[EEEEradial,{r,0,2R2}]
Plot[{ PPradial+0.03,.8ρρ-0.03,2EEEEradial},{r,0, 2 R2},PlotRange-> All, PlotLegends->{"Polarization (radial)","Bound Charge Density ((Subscript[ρ, bound]))", "Electric Field (radial)"}]
Print["Surface Charge Density on inner surface = ",σboundatR1]
Print["Surface Charge Density on outer surface = ",σboundatR2]