Potential Inside and Outside Sphere with Surface Potential V[θ] = Vo ((Sin[θ])^4)

PROBLEM:
Consider a hollow, nonconducting sphere of radius R with a surface potential given by:

	V[θ] = Vo ((Sin[θ])^4) .

Find the electric potential and the electric field V[r,θ] both inside and outside the sphere.  

Execute the following code to view a surface-density plot of this potential.

ClearAll["Global`*"] (*Leave code for user to execute*)
ParametricPlot3D[{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]},{θ,0,π},{ϕ,0,2π},ColorFunction->Function[{x,y,z,θ},ColorData["TemperatureMap"][((Sin[θ])^4)]],ColorFunctionScaling->False,Mesh->True, AxesLabel->{x,y,z}]
DensityPlot[x,{x,0,1},{y,0,0.05},ColorFunction->ColorData["TemperatureMap"],ColorFunctionScaling->False,AspectRatio->Automatic,PlotRangePadding->None,FrameTicks->{{None,None},{Automatic,None}}, PlotLabel->"Potential vs Color in Units of Vo"]

Because the boundary is a sphere, it is convenient to express the solution in spherical coordinates.  Further, the potential at the boundary does not depend on ϕ, so the resulting potential will not depend on ϕ.  In this sense, the solution is “two-dimensional” (depending on r and θ only), even though the potential and electric field are defined at every point in 3D space.  

Laplace’s equation in spherical coordinates, with no ϕ dependence, is given by:

(1/(r^2))(∂/(∂r))((r^2)((∂V[r,θ])/(∂r)))+(1/((r^2) sinθ))  (∂/(∂θ))(sinθ ((∂V[r,θ])/(∂θ)))=0

One can generate solutions for this partial differential equation by separating the variables.  One starts by assuming a solution of the form V[r,θ] = R[r] Θ[θ].  With some algebra, one can often reduce the result to two ordinary differential equations (easier!), one in R[r] and one in Θ[θ].

(a) Use M’s Laplacian command to calculate (∇^2)[R[r] Th[θ].  The syntax of the command is:

	Laplacian[V[r,θ], {r,θ,ϕ}, “Spherical”]
	
The Expand option produces a simpler result in this case.  You can also name the solution for future reference, for instance:

	Lap1 = Laplacian[V[r,θ], {r,θ,ϕ}, “Spherical”]//Expand

(* Input code here *)

ClearAll["Global`*"]
Lap1 = Laplacian[R[r] Th[θ], {r,θ,ϕ}, "Spherical"]//Expand

(b)  Multiplying each term by (r^2) and dividing by R[r] Th[θ] will eliminate the r-dependence in  the terms in Th’[θ] and Th’’[θ] and eliminate the θ-dependence from the terms in R’[r] and R[r]’’.  Use M to perform this multiplication.  The simplify option helps.  If you don’t want to copy the entire solution, you can refer to it by the name we gave it above.

(* Input code here *)
Lap2 = Lap1*(((r^2)/(R[r] Th[θ]))) //Simplify

If all went well, one set of terms will depend on r only, and the other set of terms will depend on θ only.  Since r and θ are independent variables, it is not likely that any function of r will sum to zero with a function of θ unless both functions are constant.  As you will see, it is useful to choose a constant of the form l(l+1), where l is a non-negative integer:  l = 0,1,2,3....  (In text files, I use an italics l from the Consolas font to avoid misreading it as a 1.)   

We can define an ordinary differential equation for R[r] by selecting its terms and setting them equal to  l(l+1).  To ensure a zero sum, we can define a second differential equation with Th[θ] and its terms and set it equal to –l(l+1).  This choice of “separation constants” will ensure that Laplace’s equation is satisfied by any solution of the from R[r] Θ[θ] for a given value of l.

(* Leave this code for students to execute *)
Clear[r];
rOnly = ((r (2 R'[r]+r R''[r]))/(R[r]))==l(l+1)
thOnly = ((Cot[θ] Th'[θ]+Th''[θ])/(Th[θ]))== -l(l+1)

(c)  Use DSolve to solve these differential equations.  Extract and simplify the solutions as necessary.

(* Input code below *)
Rsol = DSolve[rOnly, R[r],r]//Simplify
Thsol = DSolve[thOnly, Th[θ],θ]//Simplify

As usual, DSolve lists the solutions in brackets.  Solution functions can be extracted from the list using the replacement command.  For instance, R[r_] = R[r]/.Rsol[[1]].  The double brackets enclosing the “1” in Rsol[[1]] tell M to look for the first solution to the equation Rsol that is enclosed in two brackets.  The resulting R[r] can be used to define a function R[r_], although it is often safer to rename the function.  (M overwrites R[r].)

R[r_] = R[r]/.Rsol[[1]]
Th[θ_] = Th[θ]/.Thsol[[1]]

If your solution for R[r] looks like mine, it is complex.  M assumes l is complex, and that makes R[r] complex.  With pen and paper, you can show that the solution for R[r] can be written as: A[l](r^l) +B[l](r^(-1-2l)), provided that l is real and not negative.  The Simply command can show this, but not as a “PostFix” operation—that is, not using the //Simplify notation.  As a PostFix, Simplify does not recognize assumptions.  Assumptions can be added as options to the full-blown Simplify command.  For instance:

Simplify[R[r], Assumptions->{l ∈ Integers && l ≥ 0}]

Try it!  (O.K., I’ll do it for you):

(* Input code for Simplify operation here *)

Simplify[R[r],Assumptions->{l ∈Integers && l ≥0}]

(d)  The above solution can be rewritten as A[l](r^l) +B[l](r^(-1-2l)).  Which solutions (∝ (r^l) or ∝ (r^(-1-2l))) are appropriate for potentials inside the sphere?  Which are appropriate for solutions outside the sphere?  Do these equations make any assumptions about the location of the zero potential?

<Enter answers below>

Solutions of the form (r^l) are appropriate for solutions inside a spherical boundary, because they converge as r → 0 (the center of the sphere).  Conversely, solutions of the form (r^(-(l+1))) are appropriate for solutions outside a spherical boundary, because they converge as r → ∞.  [Since l is a non-negative integer, the exponent (l+1) in (r^(-(l+1))) never equals zero.  Therefore all of these terms go to zero as as r → ∞.]

The potential for Vout[r,θ] goes to zero as r→∞, no matter how the constant C[2] is chosen.  Like the Coulomb potential, the zero of potential in the series expansion is found as r→∞.

(e)  You should have obtained two solutions for the angle dependence, TTh[θ]:  Legendre polynomials of the first kind, LegendreP[l,Cos[θ]], and Legendre polynomials of the second kind, LegendreQ[l,Cos[θ]].  

Plot the first five LegendreP and the first five LegendreQ polynomials (for l = 0, 1, 2, 3, and 4) over the range {θ, 0, π} and answer the following questions:

Which of these functions are suitable solutions for the θ-dependence of the potential in spherical coordinates?  Why?

(* Input code here *)

ClearAll["Global`*"](* Leave ClearAll statement *)
Plot[{LegendreP[0,Cos[θ]],LegendreP[1,Cos[θ]], LegendreP[2,Cos[θ]], LegendreP[3,Cos[θ]], LegendreP[4,Cos[θ]]}, {θ,0,π}, PlotLabel-> "First Five Legendre Polynomials of the First Kind", PlotLegends->"Expressions"]

ClearAll["Global`*"]
Plot[{LegendreQ[0,Cos[θ]],LegendreQ[1,Cos[θ]], LegendreQ[2,Cos[θ]], LegendreQ[3,Cos[θ]], LegendreQ[4,Cos[θ]]}, {θ, 0, π}, PlotLabel-> "First Five Legendre Polynomials of the Second Kind", PlotLegends->"Expressions"]

<Answer questions in this text cell>

The potential should be continuous and finite over the range of possible θ-values: {θ, 0, π}.  

The Legendre polynomials of the first kind, LegendreP[l,θ], appear to be well-behaved (finite-valued and continuous] over the entire range.  They should work.  

The Legendre polynomials of the second kind, LegendreQ[l,θ], diverge as θ → 0 and as θ → π.  They are therefore unsuitable as solutions for the angular dependence of the potential over the range {0,π}.  They might be of use in a problem where the θ-dependence is somehow restricted.  Not today.

(f)  Write general, series solutions of Laplace’s equation in spherical coordinates for the potentials inside and outside the sphere: Vin[r,θ] and Vout[r,θ].  Each series will be a linear combination of terms including Legendre polynomials, LegendreP[l,Cos[θ]] with undetermined coefficients, say A[l] for Vin[r,θ] and B[l] for Vout[r,θ].

(* Input code here *)

ClearAll["`*"] (* Leave ClearAll statement *)
Assumptions-> l ∈ Integers && l ≥ 0 && {r ,θ} ∈ Reals && r > 0 && 0 ≤ θ ≤ π ;(* Leave Assumptions *)
Vin[r_,θ_] =(∑_(l=0)^∞)A[l] (r^l) LegendreP[l,Cos[θ]]
Vout[r_,θ_] =(∑_(l=0)^∞)B[l] (r^(-(l+1))) LegendreP[l,Cos[θ]]

(g) Choose a few terms from each solution and show that they satisfy Laplace’s equation in spherical coordinates.  Show your results for Vin[r,θ] and Vout[r,θ] for one of these l-values.  The Laplacian command with the //Simplify option is probably the easiest way to do this:

Laplacian[your function here, {r,θ,ϕ},”Spherical”]//Simplify

(* Input code here *)

ClearAll["`*"] (* Leave ClearAll statement *)
l=5;
Laplacian[A[l] (r^l) LegendreP[l,Cos[θ]], {r,θ,ϕ}, "Spherical"]//Simplify
Laplacian[B[l] (r^(-(l+1))) LegendreP[l,Cos[θ]], {r,θ,ϕ}, "Spherical"]//Simplify

<Input any comments in this text cell>

M says that the Laplacian for the l = 5 terms in each series is zero.  I haven’t tested all the possibe values of l, but if each term (by itself) satisfies Laplace’s equation, any sum of such terms will also satisfy Laplace’s equation.  We leave the general proof for another day.

(h) Write out the boundary conditions in terms of these expansions.

(* Input code here *)

Vin[R,θ] ==Vo ((Sin[θ])^4) 
Vout[R,θ] ==Vo ((Sin[θ])^4)

The constants {A[l]} and {B[l]} may be determined from the boundary conditions by using the Fourier trick.  You multiply both sides of the boundary condition equations by the mth order Legendre polynomial, LegendreP[m,Cos[θ]]×Sin[θ] and integrate from 0 to π.  (The factor of Sin[θ] appears when you transform the argument of LegendreP[m,x] from x to Cos[θ].)  You must not assume that the order m equals l in LegendreP[l,Cos[θ] above. 

The integrals for the Fourier trick have the following form:  

(∫_0^π)Vin[R,θ] LegendreP[m,Cos[θ]] Sin[θ]dθ   =(∫_0^π)Vo ((Sin[θ])^4)LegendreP[m,Cos[θ]] Sin[θ]dθ 
(∫_0^π)Vout[R,θ] LegendreP[m,Cos[θ]] Sin[θ]dθ =(∫_0^π)Vo ((Sin[θ])^4)LegendreP[m,Cos[θ]] Sin[θ]dθ 

You will also need the orthogonality relation for Legendre polynomials.

(∫_0^π)P[m,Cos[θ]]P[n,Cos[θ]]Sin[θ]dθ =({{{, ({{0, m ≠ n}, {(2/(2l+1)), m=n}})}})    = ((2 (Subscript[δ, mn]))/(2l+1))

(i) Use your series expansions for Vin[] and Vout[] and the orthogonality relation to simplify the LHS of the “Fourier Trick” expression.  It is safest to do this by hand.

<Input derivation in this text cell>

The left hand side of the Fourier trick expression is

(∫_0^π)Vin[R,θ] LegendreP[m,Cos[θ]] Sin[θ]dθ

where Vin[r,θ] =(∑_(l=0)^∞)A[l] (R^l) LegendreP[l,Cos[θ]] from Part (a) above.  Replacing Vin[R,θ] with the corresponding infinite series gives:

(∫_0^π)(∑_(l=0)^∞)A[l] (R^l) LegendreP[l,Cos[θ]] LegendreP[m,Cos[θ]] Sin[θ]dθ

The factors that do not depend on θ can be taken outside the integral.  This yields:

(∑_(l=0)^∞)A[l] (R^l)(∫_0^π)LegendreP[l,Cos[θ]] LegendreP[m,Cos[θ]] Sin[θ]dθ

The form of the expression inside the integral now matches the orthogonality relation, and therefore equals  ((2 (Subscript[δ, lm]))/(2l+1)).  The LHS then becomes:  

(∑_(l=0)^∞)((2 A[l] (R^l)(Subscript[δ, lm]))/(2l+1)) = ((2 A[m] (R^m))/(2m+1))

The only nonzero term in the sum is that for which l = m, leaving

LHS for Vin = ((2 A[m] (R^m))/(2m+1))

The boundary condition for Vout[r,θ] differs from the boundary condition for Vin[r,θ] in two respects:  A[l] is replaced by B[l], and (R^l) is replaced by by (R^(-(l+1))).  Repeating the same series of operations on the boundary conditions for Vout yields.  

LHS for Vout = ((2 B[m] (R^(-(m+1))))/(2m+1))

(j)  Now consider the RHS, which is the same for Vin[r,θ] and Vout[r,θ].  In general, the integral will depend on m, the order of the Legendre polynomial.  

Use M to evaluate the RHS for m = 0 through 10.  Put the results in a table.  Comment on your results.  For instance, some terms are zero.  If you can, explain why the zero terms might be expected to equal zero.

(* Input code here *)

ClearAll["`*"]  (* Leave ClearAll Statement *)
RHS[m_]=(∫_0^π)Vo ((Sin[θ])^4)LegendreP[m,Cos[θ]] Sin[θ]\[DifferentialD]θ
RHSTable = Column[Table[{m, RHS[m]},{m,0, 10}]]

<Input comments in this text cell>

The RHS is zero except for m = 0, 2, and 4.  These are the first three “even” numbered Legendre polynomials.  The even numbered Legendre polynomials include only even powers of Cos[θ].  

((Sin[θ])^4) can be expressed as the sum of three terms with the form ((Cos[θ])^n), where n is even. 

((Sin[θ])^4)=(((1-((Cos[θ])^2)))^2) = 1 – 2 ((Cos[θ])^2) + ((Cos[θ])^4) 

The same three powers of ((Cos[θ])^n) (n = 0, 2, 4) appear in LegendreP[4,Cos[θ]].  Note that the LegendreP[2,Cos[θ]] polynomial includes the first two powers of ((Cos[θ])^n) (n = 0, 2).  Likewise the LegendreP[0,Cos[θ]] polynomial includes only the n = 0 term.  So there ought to be a linear combination of LegendreP[0,Cos[θ]], LegendreP[2,Cos[θ]], and LegendreP[4,Cos[θ]] that sums to ((Sin[θ])^4).  The Fourier trick is one way of determining how much weight to give each Legendre polynomial in the sum.  

Since there are only three nonzero terms in the series, we may conclude that the sum of these three terms yields the “exact” solution to the potential.  The result is not an approximation.  Nice.

(k)  Derive expressions for A[l] and B[l].  Put the results in a table.  I left the integral from the RHS side in my expressions, as it is simple enough to recalculate the value of the integral.  Another approach would be to copy your results for the integration from Part (e).

<Input derivation in this text cell>

Equating LHS to RHS we have

LHSin =((2 (R^l))/(2l+1)) A[l]=  RHS = (∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]dθ

So that: 

A[l]=(( 2l+1)/(2 (R^l)))(∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]dθ

Similarly 

LHSout =(2/((2 l+1)(R^((l+1))))) B[l] = (∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]dθ  

Solving for B[l] gives:

B[l] = (((2 l+1)(R^((l+1))))/2)(∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]dθ 

We know that only the l = 0, 2, and 4 terms of the series are nonzero.  I calculate a few more that yield zeros 

The M-code I used follows.

(* Input code here *)

ClearAll["`*"]
Assumptions →{l ,m}∈ Integers && {l,m}≥0; (* Leave Assumptions *)
lmax = 10;
A[l_]=(( 2l+1)/(2 (R^l)))(∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]\[DifferentialD]θ
B[l_] = (((2 l+1)(R^((l+1))))/2)(∫_0^π)Vo ((Sin[θ])^4)LegendreP[l,Cos[θ]] Sin[θ]\[DifferentialD]θ
Column[Table[{l,A[l], B[l]},{l,0,lmax}]]

(l) Use your values of A[l] and B[l] to compute Vin[r,θ] and Vout[r,θ] in terms of Cos[θ] and r.  If you do not clear your constants, M will remember and use the above values.

(* Input code here *)
Vin[r_,θ_] =(∑_(l=0)^lmax)A[l] (r^l) LegendreP[l,Cos[θ]]
Vout[r_,θ_] =(∑_(l=0)^lmax)B[l] (r^(-(l+1))) LegendreP[l,Cos[θ]]

(m) Use Plot3D[] to plot the potential inside and outside the sphere for reasonable choices of R and Vo.  I set R = 1 m and Vo = 1 V.  Remember that Plot3D[] requires the function to be expressed in terms of Cartesian coordinates, {x,y,z}.  You will have to choose two of the three Cartesian coordinates, as you need one coordinate axis to show the potential values.

Will you plot V[x,y], V[y,z], or V[z,x]?  Why?  What substitutions are needed to express the above potential in terms of your Cartesian coordinates?  Write and run your M code.

<Describe your plan here>

I plotted V[y,z], assuming that x = 0.  One expects the potential to depend on x and y in similar ways because of the symmetry.  By symmetry, rotations about the z-axis will not change the potential anywhere on the sphere.  For instance, one can “rotate” the x-axis around the z-axis so that it lies on top of the old y-axis.  

The z-axis is special.  Cos[θ] = z/r.  With z, either x or y would be a good choice.  I plotted V[y,z].  The following “replacements” were required:

r → ((((x^2)+(y^2)))^(1/2))
Cos[θ] → z/r.

(* Input code here *)

Vo = 1; R = 1;
r = ((((y^2)+(z^2)))^(1/2));
Vin[y_,z_] =((8/15)-(8/21) (r^2) (-1+3 ((Cos[θ])^2))+(1/35) (r^4) (3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4)))/.Cos[θ]->z/r
Vout[y_,z_] =((8/(15 r))-((8 (-1+3 ((Cos[θ])^2)))/(21 (r^3)))+((3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4))/(35 (r^5))))/.Cos[θ]->z/r
Plot3D[If[r<R, Vin[y,z], Vout[y,z]],{y,-2,2},{z,-2,2}, AxesLabel->{y,z, V[y,z]}, PlotLabel->"Potential Inside and Outside Sphere with V[θ] = Sin[θ]^4"]

(n) Describe your plot below.  Where is the z-axis of the sphere?  Where is the equator of the sphere?  Do the peaks and troughs of the potential function match the maximum and minimum values of the ((Sin[θ])^4) function?

<Input your description in this text cell>

I plotted V[x=0,y,z]:  the potential as a function y and z for the special case of x = 0.  The x-values are not shown.  The z-axis is the conventional z-axis for spherical coordinates.  The N and S poles of the sphere are located at points {y,z} = {0,±1}.  Of the points lying on the equator, the plot shows only two, {y,z} = {±1,0}.  

The maximum potential of Vo = 1 is reached at points {y,z} = {±1,0}, which corresponds to θ =π/2, on the equator of the sphere.  This coincides the the maximum value of ((Sin[θ])^4) when θ = π/2.  Likewise, the minimum value of ((Sin[θ])^4) = 0 is reached at θ = 0 and π, which corresponds to points {y,z} = {0,±1} at the poles of the sphere.

It is remarkable how quickly the θ-dependence of the potential drops as you move away from the sphere.  In the limit of large r, we expect the potential to vary as (r^(-1)), like the potential of a point charge.  The first term in the series for Vout[r,θ] varies as (r^(-1)) does not depend on θ.  The higher order terms (which do depend on θ) both decay quickly with increasing r (as (r^(-3)) and (r^(-5))).  These higher-order terms are responsible for the fast decay observed above.

(o) Before we plot the electric field, it is useful to show that the field component that you do not plot is zero.  For instance, if you calculated V[y,z] above, you will want to plot the electric field in the y-z plane, where x = 0.  If the x-component of the field is not zero when x = 0, your plot will not show it. That would be bad.  

First, write Vin[r,θ] and Vout[r,θ] in terms of all three Cartesian components.  Then use the D[] command (partial derivative) to calculate the one component of the electric field that you will not show in your field plot.  For instance, if I hope to plot the field in the y-z plane, I hope that the x-component of the field in the y-z plane is zero.  Use the replace all command (for instance, /.x→0) to see if the x-component of the field zero when x = 0.

(* Input code here *)

Vo = 1; R = 1;
r = ((((x^2)+ (y^2)+(z^2)))^(1/2));
Vin[x_, y_,z_] =((8/15)-(8/21) (r^2) (-1+3 ((Cos[θ])^2))+(1/35) (r^4) (3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4)))/.Cos[θ]->z/r;
Vout[x_, y_,z_] =((8/(15 r))-((8 (-1+3 ((Cos[θ])^2)))/(21 (r^3)))+((3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4))/(35 (r^5))))/.Cos[θ]->z/r;
-D[Vin[x,y,z],x]/.x->0          (* x-component of Vin for x = 0 *)
-D[Vout[x,y,z],x]/.x->0        (* x-component of Vout for x = 0 *)

<Describe your results in this text cell>

I plan to plot the field in the y-z plane, where x is zero.  The x-component of the field is the partial derivative of the potential with respect to x.  When I set x to zero, both partial derivatives are zero.  When I plot the fields as functions of y and z (with x set to zero), I can be confident that my plot shows the true direction of the electric field.  When x = 0, the all the field vectors will lie in the y-z plane.

(p) Construct a ContourPlot of the potential and a StreamPlot of the electric field.  Superimpose them.  Describe the result.  Do you see  evidence for negative surface charge on the sphere?

(* Input code here *)

Vo = 1; R = 1;
r = ((((y^2)+(z^2)))^(1/2));
Vin[y_,z_] =((8/15)-(8/21) (r^2) (-1+3 ((Cos[θ])^2))+(1/35) (r^4) (3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4)))/.Cos[θ]->z/r;
Vout[y_,z_] =((8/(15 r))-((8 (-1+3 ((Cos[θ])^2)))/(21 (r^3)))+((3-30 ((Cos[θ])^2)+35 ((Cos[θ])^4))/(35 (r^5))))/.Cos[θ]->z/r;
EEin[y_,z_]=-Grad[Vin[y,z],{y,z}]//Simplify
EEout[y_,z_]=-Grad[Vout[y,z],{y,z}]//Simplify
contoury = ContourPlot[If[r<R, Vin[y,z], Vout[y,z]],{y,-2,2},{z,-2,2}, PlotLabel->"V[y,z] for Sin[(θ^4)] Potential Distribution"]
streamy = StreamPlot[If[r<R, EEin[y,z], EEout[y,z]],{y,-2,2},{z,-2,2}]
Show[contoury, streamy]

<Input your description in this text file>

As expected, the electric field is directed outward from parts of the sphere with high positive potentials, near the equator.  Interestingly, the electric field is directed inward toward parts of the sphere with low (but still positive) potentials, near the poles of the sphere.  

The presence of field lines directed inward, toward the surface, at the poles of the sphere indicates the presence of negative surface charge in these regions.  The absence of negative potential does not imply the absence of negative charge.  

The D-shaped contours at the poles of the sphere are interesting.  I don’t know what to make of them, but they are interesting.

(q) Griffiths claims that the potential at the center of a sphere is equal to the average potential on the surface of the sphere (Section 3.1.4).  

Find the average potential of the surface of the sphere and compare it to the potential at the center of the sphere.

The average potential on the surface of the sphere is given by

∫V[R,θ,ϕ]dA  ÷ ∫dA

where dA is the surface area element, (R^2) Sin[θ] dθ dϕ = 2π (R^2) Sin[θ] dθ.  The surface area of the sphere is given by ∫dA = 4π (R^2).  For Vo = 1, the average is given by:

(* Input code for average potential on sphere surface here *)

Average = Integrate[((Sin[θ])^4) 2π Sin[θ], {θ,0,π}]/(4π)

The potential at the center of the sphere is given by the first term in the Legendre expansion for Vin[r,θ], that is A[0].  (All the other terms in the Legendre series for Vin[r,θ] equal zero when r = 0.)  If the variable A[0] has not been cleared or overwritten, you can still query M for its value.  Otherwise, you can manually check your results above.

(* Input code to find potential at center of sphere here *)

A[0]

For R = 1 and Vo = 1, the average potential on the surface of the sphere is 8/15.  Under the same set of assumptions, the potential at the center of the sphere (=A[0]) also equals 8/15.  

These results are consistent with Griffiths’ theorem stating that the potential at the center of a sphere is equal to the average potential on the surface of the sphere (Section 3.1.4).   

If you had never done the calculation, this theorem would tell you that the potential at the center of the sphere must not be zero, because the potential on the surface of the sphere is zero or positive everywhere (so it’s average over that surface is not zero).  

Note that the electric field is a different beast:  it is zero at the center of the sphere (easily seen by symmetry).