Green Function for Dirichlet Problem on Sphere

Φ(a,θ) = V,  for 0<θ< π/2 and  Φ(a,θ) = -V,  for π/2<θ< π.

Clear["`*"]

The Green function is

G(r,θ,ϕ; a,ψ,χ) = ((a ((r^2)-(a^2)))/(((((a^2) + (r^2) - 2 a r x))^(3/2)))) 

where  x = Cos[(Subscript[θ, 12])] = Cos[θ]Cos[ψ]+Sin[θ]Sin[ψ] Cos[χ-ϕ].
 
 In our case, there is azimuthal symmetry so Φ must be independent of the angle ϕ. We may, therefore, simply evaluate the potential at ϕ = 0.

greenf  =  
      a (r^2 - a^2)/(a^2 + r^2 - 2 a r x)^(3/2)

We can express green in terms of the ratio ρ=(r/a)

g =  
      Simplify[greenf /.  r ->  a ρ, 
        Assumptions  ->  a > 0]

Expand g in an asymptotic series in ρ and retain the "n" leading terms

n = 10

g = Normal[Series[g, {ρ, Infinity, n}]]

Substitute for the angle variable x in terms of the integration and observation angles

x =  
      Cos[θ] Cos[ψ]  +  
        Sin[θ] Sin[ψ] Cos[χ]

Expand g so that it is simpler for mathematica to integrate

eg =  Expand[g]

Do the integral over χ

part1 =  Integrate[eg, {χ, 0, 2  π}]

Now, do the integral over ψ, divide out the factor of 4π, and expand the result

part2  =  
      Expand[(Integrate[Sin[ψ] part1, {ψ, 0, π/2}] -  
              Integrate[
                Sin[ψ] 
                  part1, {ψ, π/
                    2, π}])/(4  π)]

Substitute to express your answer in terms of μ = Cos[θ]

Φ0  =  
      TrigExpand[
        FullSimplify[part2 /.  θ ->  ArcCos[μ]]]

Pick out leading term

len =  Length[Φ0]

coef  =  Φ0[[len]]/μ

Factor out the leading power

res  =  Expand[Φ0/coef]

Put the result back in the form of an asymptotic expansion

res  =  Normal[Series[res, {ρ, Infinity, n}]]

lnr  =  Length[res]

Express the coefficients of (ρ^(-2 n)) in terms of Legendre polynomials

Do[t[i] =  
        Simplify[ρ^(2 i) res[[lnr + 1 - i]]/
              LegendreP[2 i + 1, μ]], {i, 1, lnr - 1}]

Table[t[i], {i, 1, lnr - 1}]

Put it back together in a neat form

For[{k = 1, fac =  (Subscript[P, 1])[μ]}, k < lnr, k++, 
      fac = fac +  
          t[k] (Subscript[P, 2 k + 1])[μ]/ρ^(2 k)]

fac  /.  ρ ->  r/a

Here is the resulting asymptotic expansion of Φ in terms of Legendre Polynomials

Φ  = 
      V  Expand [coef fac  /.  ρ ->  r/a]

V (Normal[
          Series[(1 - (z^2 - a^2)/(z Sqrt[z^2 + a^2])), {z, 
              Infinity, 10}]])