Interior solution of azimuthally symmetric case of Laplace equation in orthogonal similar oblate spheroidal coordinates

Curvilinear coordinate systems distinct from the rectangular Cartesian coordinate system are particularly valuable in the field calculations as they facilitate the expression of boundary conditions of differential equations in a reasonably simple way when the coordinate surfaces fit the physical boundaries of the problem. The recently finalized orthogonal similar oblate spheroidal (SOS) coordinate system can be particularly useful for a physical processes description inside or in the vicinity of the bodies with the geometry of an oblate spheroid. Such shape is aproximating well objects investigated within astrophysics. The solution of the azimuthally symmetric case of the Laplace equation was found for the interior space in the orthogonal SOS coordinates. In the frame of the derivation of the harmonic functions, the Laplace equation was separated by a special separation procedure. A generalized Legendre equation was introduced as the equation for the angular part of the separated Laplace equation. The harmonic functions were determined as relations involving generalized Legendre functions of the first and of the second kind. Several lower-degree functions are reported. Recursion formula facilitating determination of the higher-degree harmonic functions was found. The general solution of the azimuthally symmetric Laplace equation for the interior space in the SOS coordinates is reported.


Introduction
Curvilinear coordinate systems different from the rectangular Cartesian coordinate system are of value in simple geometric considerations as the cumbersome mathematical expressions for determination, for example, of areas and volumes can be overcome.Moreover, curvilinear coordinate systems are particularly valuable in the field theory and differential equations solution.In order to express boundary conditions of differential equations in a reasonably simple way, coordinate surfaces that fit the physical boundaries of the problem are essential [1].Therefore, a range of field problems that can be handled effectively depend upon the number of well-developed coordinate systems.Particularly, the orthogonal curvilinear coordinate systems [2] appear to be the most useful.
The recently finalized similar oblate spheroidal (SOS) orthogonal coordinate system [3,4] can be a powerful tool for a description of physical processes inside or in the vicinity of the bodies with geometry of an oblate spheroid.Such bodies range (but are not limited to) from planets with a small oblateness (like the Earth with the ratio of the semi-axes difference to the major semi-axis length ≈1:300), through elliptical galaxies up to significantly flattened objects like disk galaxies.
In the field of atmospheric physics, it was stated [3] that the SOS coordinates could be of help for better modeling of geopotential surfaces, allowing for a better description of the spatial variation of the apparent gravity.Still more advantageously, the SOS coordinates could be used for modeling of the potential and atmosphere of significantly more oblate celestial bodies, e.g., the gas giants.Further, similar oblate spheroids are frequently used for modeling of iso-density levels inside galaxies [5,6].Therefore, the SOS coordinates can be helpful also in this field.
In electrostatics and solid-state physics, the SOS coordinates could find application in description of electric field potential in ferroelectric materials (e.g., ferroelectric nanocomposites) containing dielectric inclusions (or vice versa-ferroelectric inclusions in dielectric matrix) [7], particularly the inclusions of spheroidal shape [8].
The similar oblate spheroidal coordinates are distinct from all the standard orthogonal coordinate systems [2,9] including the confocal oblate spheroidal system, as one coordinate surfaces family of the SOS coordinates is not of the second degree (and even not of the fourth degree) but it is formed by general power functions z ≈ x 1+μ (i.e., with a real-number exponent 1 + μ) rotated around the z-axis [3,4].They are orthogonal to the similar oblate spheroids representing the first set of the coordinate surfaces.(For terminological clarity, a spheroid means in this text an ellipsoid of revolution or rotational ellipsoid.An oblate spheroid is a quadric surface obtained by rotating an ellipse about the shorter principal axis).(2024) 139:409 Already reported [4] were the analytical coordinate transformation from the SOS coordinates to the Cartesian system as well as the metric scale factors and the Jacobian determinant.Further, the formulas necessary for the transformation of a vector field between the SOS system and the Cartesian coordinates were recently published [10].Generalized sine and cosine applicable in the transformation were introduced as well.The derived analytical relations cannot be expressed in a closed form; nevertheless they employ convergent infinite power series (with generalized binomial coefficients) and are thus still analytical.
Nevertheless, there are still missing significant parts of algebra connected with the SOS coordinate system.Predominantly, solutions of standard differential equations in the SOS coordinates are not derived.It is known that several important scalar differential equations can be reduced to the Helmholtz equation or to its special case, the Laplace equation [2].Therefore, the scalar fields in several areas of physics can be based on the solutions of the Helmholtz and of the Laplace equations.
This article deals with the task to solve the simpler of the two types of the abovementioned differential equations, the Laplace equation, in the SOS coordinate system.A solution of the Laplace equation would represent an important step in the development of the SOS coordinate system and its applicability in astrophysics (e.g., for modeling of gravitational potential) and in other fields of physics (cosmology, geodesy, climatology, electromagnetism).Determination of harmonic functions could also help to find solutions of the other, more complex, differential equations in the SOS coordinates.
The organization of the text is following.First, a short summary of the SOS coordinate system is provided.Then, the azimuthally symmetric Laplace equation in the SOS coordinates is formulated and a possible shape of the harmonic function is modeled.Special separation procedure is reported next, and the radial part of the equation is solved.Simple but informative enough solutions of the angular part of the separated equation are found.The found simple solutions as well as the Laplace equation itself are rewritten to the form enabling solution generalization.A new type of polynomials-generalizing the Legendre ones-are found to be a part of the harmonic functions in the SOS coordinates for the interior space.A recursion formula for the polynomials is derived.Finally, the complete interior solution of the Laplace equation in the azimuthally symmetric case in the SOS system is reported.
Occasionally during the derivation, the found interim results for the limiting case of oblate spheroids, i.e., for the spherical case, are compared with the well-known solution for the spherical coordinates.
A number of formulas needed within the article and derived on the basis of binomial identities are reported in Appendix A. Nevertheless, due to the fact that some derivations with an extensive use of combinatorial identities are very lengthy when written in a full detail, a large part of them are moved to Supplements to this article.

The SOS coordinates
The SOS coordinates were introduced by White et al. [3].The full description of the analytical solution of SOS coordinates can be found in [4].A shortened, but more overviewable summary is reported in [10].It is recommended to find details of the system in the abovementioned references.A limited description (to the extent needed for the tasks of this article) is given also in this section.
For the SOS coordinate system (R, ν, λ) [4], the basic coordinate surfaces of the R coordinate are similar oblate spheroids given in 3D Cartesian coordinates (x 3D , y 3D , z 3D ) by the formula The R coordinate value is equal to the equatorial radius of the particular spheroid from the family.The parameter μ characterizes the whole family of the similar oblate spheroids.The parameter μ > 0 for oblate spheroids, and the minor and major semi-axes of each member of the spheroid family have the ratio (1 + μ) −1/2 .As a limit (when μ 0), a sphere (and spherical coordinate surfaces) is determined by (1).A special reference spheroid is introduced with the equatorial radius R 0 , usually coinciding with the reference surface of the object for which the SOS coordinate system is to be applied.
The second set of the coordinate surfaces, orthogonal to the similar oblate spheroids defined above, are power functions in 3D of the shape [3,4] The labeling, i.e., the coordinate ν corresponding to these surfaces, is equivalent to the so-called parametric latitude [4].The coordinate ν is also equivalent to the parameter used for the standard parametric equation of the ellipse, representing the meridional section of the reference spheroid mentioned above, having a special equatorial radius (major semi-axis) equal to R 0 , i.e., Using this definition, the coordinate ν is equal to zero in the equatorial plane while it is equal to ± π/2 on the rotation axis (+ in the north,-in the south).From Eqs. ( 3) and ( 17) of [4], the relation for the power function coordinate surfaces in 3D (2) in dependence on the coordinate ν can be detailed as (the correctness of the relation can be easily tested when coszν and sinzν according to (3) are plugged into this relation).Finally, the third set of the coordinate surfaces, orthogonal to the previous two, are semi-infinite planes containing the rotation axis.The associated coordinate is the longitude angle λ, which is the same as its equivalent coordinate in the spherical coordinate system.
A key role in the derivation of the SOS coordinate algebra (e.g., the coordinates transformation between the SOS and the Cartesian system [4]) plays the dimensionless parameter W defined as which will be used frequently in this article.The constant-W surfaces are straight half-lines starting at the origin and rotated around the system axis, i.e., cones [4].In what follows, the calculation is restricted to the parametric latitude ν ∈ 0, π 2 .With this limitation, the parameter W , see Eq. ( 5), is always nonnegative which simplifies further derivations.Due to the symmetry (reflection with respect to the equator), expressions and solutions relevant for the SOS system in the complementary range ν ∈ − π 2 , 0) can be easily obtained.
The SOS coordinates can be transformed to the Cartesian coordinates using analytical expressions including infinite power series with generalized binomial coefficients [4,10].Relations enabling such approach were first reported by Pólya and Szegö [11].Due to the convergency limits of the involved power series, the expressions were derived separately in so-called "the small-ν region" and separately in "the large-ν region" [4].The border between the two regions is defined by the parameter W value (see Eq. ( 5)) fulfilling the following relation: This border surface is a straight half-line starting at the origin and rotated around the system axis (and forming thus a surface of a cone) [4].Although W border is a constant for a fixed μ, the coordinates R and ν vary along the half-line.

Important basic relations of the SOS coordinate system
In order to fulfill the aim of this article, some already derived relations for the SOS coordinates [4,10] are needed.Therefore, the relevant ones are listed in this section.The formulae are listed separately for the small-ν region and for the large-ν region.
First needed relations are the metric scale factors [4,10].Below are listed the metric scale factors for R coordinate, as well as for ν coordinate, , respectively.The left side (or the first one) of the above relations shows the formula valid for the small-ν region while the right side shows the formula valid for the large-ν region.The partial derivative in the h ν relations can be expressed as For completeness, we also report the partial derivative of W with respect to R, which is easily obtained from ( 5) Further, the already derived Jacobian determinant h R h ν h λ [4,10] is listed: Again, the first relation shows the analytic formula valid for the small-ν region, while the second relation shows the formula valid for the large-ν region.The Jacobian divided by the square of the h R metric scale factor is expressed as for the small-ν region and for the large-ν region, respectively [4,10].
The Jacobian divided by the square of the h ν metric scale factor is derived in Appendix A of this article both for the small-ν region (the first relation, (A8) in the appendix) and for the large-ν region (the second relation, (A9) in the appendix): Finally, the functions for the generalized cosine (f C ) and for the generalized sine (f S ) [10] in the frame of the SOS coordinate system are reported: The left sides (the first realtions) show the formulae valid for the small-ν region, while the right sides (the second relations) for the large-ν region.As f S , f C depend solely on W (not separately on R and ν), these functions are constant at the same surfaces in 3D where W is constant.As (see the text after ( 5)) the constant-W surfaces are cones with the tip at the origin and the axis coinciding with the rotation axis of the coordinate system, so are also f S , f C constant on the same surfaces.
For μ 0, the generalized cosine and the generalized sine become cosine function and sine function, respectively.Similarly as for the standard cosine and sine, the following relation is valid for the generalized cosine and the generalized sine: Further (see [10] or Appendix A, relation (A13), and using ( 16)) the relations of the generalized cosine and of the generalized sine to the metric scale factor h R are The functions listed in this section all contain infinite power series with generalized binomial coefficients.Operations with them are not as straightforward as a simple multiplication or summation of two expressions.Nevertheless, it is possible to deal with them with a help of known combinatorial identities [11][12][13][14][15].The important ones, used in the derivations also in this article, are listed in [10].With their help, a number of more complex relations needed for the Laplace equation solution are derived in Appendix A and subsequently used in the following sections.

Laplace equation in SOS coordinates in azimuthally symmetric case
The Laplace equation, V 0 (18) in the SOS coordinates can be derived using the general formula for the Laplace operator .In general curvilinear orthogonal coordinates, it is reported in [16], Eq. ( 15).Applied to the SOS case and the symbols for the SOS coordinates R, ν, λ, the operator can be written as where V (R,ν) denotes the potential function (or generally any function with physical meaning on which the Laplace operator is to be applied) and h R h ν h λ is the Jacobian determinant.The Laplace Eq. ( 18) in the SOS coordinates in azimuthally symmetric case (i.e., when the potential function is independent of λ) is then The basic strategy for its solution is to modify the equation algebraically in such a way that it enables variables separation, and thus transformation of the partial differential equation to a set of ordinary differential equations.
It follows from the shape of the involved functions (the metric scale factors and the Jacobian) in the SOS case that a standard separation procedure using Stäckel matrix (see e.g., [2]) is not applicable.Therefore, either separation is not possible, or it is a non-standard one, and another way for separation thus has to be searched for.
To test which of the abovementioned options is correct, a model potential was input into the Laplace Eq. (20).With the experience with the Laplace equation in the spherical coordinates, the model of the solution of the Laplace equation in the similar oblate spheroidal coordinates was first tested in the form of a product of three functions plus a constant V C : It was found in the first round of derivation that at least one (but certainly not the only one) solution in this form exists, and it is equal to where c 1 is an arbitrary real constant, f S is the generalized sine (15) and h R is the metric scale factor of the R coordinate in the SOS coordinate system [4].Both functions f S and h R depend on the composed variable W (see Eq. ( 5)), not individually on the coordinates R and/or ν.Although the derivation of the result ( 22) is non-trivial, it is not reported here in detail as the sequence of the steps leading to it is very similar to the one which will be reported in the second applied model further in the text.The found result (22) indicates that-at least for some solutions of the Laplace equation in the SOS coordinates, or for some region in the space-the solution can be written without an explicitly expressed dependence on ν coordinate, only with the explicit dependence on R and on W , moreover appearing in two separate functions, a product of which (plus a constant) represents the solution.Therefore, the second applied model function reflecting this finding and to be tested was of the form where r and F are functions of R and W , respectively.This new model excludes G(ν) function contained in the model (21); on the other side, it tests a more general function r(R) for the radial part of the model.
When we insert the model potential (23) into the Laplace equation in the azimuthally symmetric case (20), we obtain the equation in the form By an intensive use of the product and of the chain rules, and by the use of relations (A96), (A97) and (A38) derived in Appendix A, the equation can be rewritten into the form where f C and f S are the generalized cosine and sine reported in the previous chapter.This extensive derivation can be found in Supplement B, Sect.B1.
It can be easily checked that when setting μ 0 (i.e., in the spherical case) the equation simplifies (thanks to the characteristics of the generalized cosine and sine, as well as the metric scale factor h R ) to Further taking into account the form of derivatives of F(W ), see Eqs. (A44) and (A47) in Appendix A, we obtain-in this special case, i.e., for μ 0-the Laplace equation in the spherical coordinates: (it has to be taken into account that ν is latitude, not the polar angle, and that W , according to its definition by (5), does not depend on R in the case when μ 0, but solely on ν).

Separation of Laplace equation in the SOS coordinates
For a solution of Eq. (25), i.e., the Laplace equation in the SOS coordinates, we can apply a special variable separation approach.It can be seen that the first two terms in (25), containing the first and the second derivative of F(W ) and closed by curl brackets, depend only on the composed variable W , not explicitly on R, while the third and the fourth terms (the middle ones, containing the first and the second derivative of r(R)) depend only on R. The last term (also closed by curl brackets) is a product of a pure W -dependent function and a pure R-dependent function.It does worth to notice that there is no explicit dependence on the coordinate ν, as also the functions h R , f C and f S depend only on the composed variable W .The coordinate ν enters the equation only through W .
Equation (25) can be thus written in the following simplified form: with a clear correspondence of the individual terms in (28) with the terms in (25).In order this equation being fulfilled for any W and for any R, either b(R) and d(R) has to be equal to constants for any R, or a(W ) and c(W ) has to be equal to constants for any W . Otherwise, two different equations would arise for a(W ), c(W ) or for b(R), d(R) functions when keeping R equal to a particular value R 1 , and when keeping R equal to another particular value R 2 (or also similarly when keeping W W 1 , and when keeping W W 2 ).If b(R 1 ) is not equal to b(R 2 ), then the equations would be generally two different equations for a(W ) and c(W ), and they could not be fulfilled at once for all possible W .It cannot be compensated by different values d(R 1 ) and d(R 2 ).The same reasoning can be done for other possible cases.When the first possibility-i.e., a(W ) and c(W ) (according to (25)) are constants-is tested, it is concluded (after extensive algebraic operations not reported here) that the Laplace equation cannot be fulfilled in this case.
The second possibility-i.e., the case when b(R) and d(R) are constants-leads to the following relations (easily obtained by comparison of (28) and (25)): where K d , K b are separation constants.The first differential Eq. (31) has the solution for r(R) in the form This relation tells us that the radial part of the Laplace equation (at least in some solution of the Laplace equation, or in some region in space) is basically the same as for the spherical case solution, except the equatorial radius R is used instead of radius.
It follows from the second Eq.(32), see Supplement B, Sect.B2 for details, that the following equation connects the separation constants K b and K d : Then, the remaining part of the Laplace Eq. (25), i.e., (in the abbreviated notation) can be written-with a help of the separation constant K d -in the form (see (B36) in Supplement B, Sect.B2) Although it is not a fully rigorous terminology, this remaining part of the Laplace equation (and its algebraic modifications) in the SOS coordinates-depending solely on the composed parameter W -will be sometimes denoted as "angular part of the Laplace equation" in the following text, as it is the part of the Laplace equation depending also on the meridional coordinate ν hidden in W .

Relatively simple but informative solution of the angular part of the Laplace equation
The further strategy for the solution is to modify algebraically (36) to the form containing only the metric scale factor h R , not f C and f S , and to see how that form can be solved.With the help of relations (17), and also (A22) from Appendix A, and after rather extensive algebraic manipulations (see Supplement B, Sect.B2 for details), ( 36) is modified to the form containing-in the pre-factors of the derivatives and of the function F(W ) itself-only powers of h R function and the parameter W : Equation ( 37) is a rather complicated ordinary second order differential equation for F(W ), moreover depending also on the separation constant K d .Its solution is not trivial and-before finding a general solution-it would be of advantage to find at least one relatively simple solution.There is certainly a trivial solution (F(W ) equal to zero) and-for K d 0 which makes the last term zero-another clear solution is also when F(W ) is a constant.However, these solutions are not informative enough to enable steps forward in finding all the possible solutions.Therefore, still relatively simple, but a non-trivial non-constant solution has to be searched for.
Equation (37) can be possibly written also with a help of an unknown g(W ) function as follows [17]: where f (W ) is the factor by the first derivative of F(W ) in (37), i.e., By a comparison of the factor by the function F(W ) in ( 37) and (38), and by substituting f (W ) according to (39), we obtain the first order differential equation for g(W ) in the form: Then, a solution of the second order differential equation for F(W ) is transferred to a solution of the first order differential equation for g(W ).If such g(W ) could be found, then an exact solution of the partial Laplace Eq. (37) would be possible to find as well (see [17]).Nevertheless, the above first order nonlinear differential Eq. (40) for g(W ), which is the Ricatti general equation [17], has to be solved first for g(W ).
Note that writing the Laplace equation (the separated part of it depending only on W , to be precise) in the special form of (38) equation is possibly not the only option how to find a solution.Probably, other solutions could be found in cases when the equation cannot be written in the form of (38).Nevertheless, (38) could lead to one of the simplest solutions.The search for another/others would be facilitated when at least one solution is known.
Further, we still simplify the task by setting K d 1, which could lead to still more specific solution of (40).Nevertheless, even in such case the equation is not trivial to solve.Again, extensive trials and algebraic derivations (see Supplement B, Sect.B3 for details), which use (A75), (A111), (A29) derived in Appendix A, are employed for solving (40).The particular solution for K d 1 is then found with the shape where the generalized cosine function f C is given by (14).Equation ( 38), in which we know the function g(W ), has then the particular solution (see [17]): According to Appendix A, relation (A111), the integral containing the metric scale factor h R can be expressed in the following form: and the particular solution of the angular part of the Laplace equation is thus Thanks to the identities of [11], see Appendix A, (A1), (A2), (A3), (A4), (A5), this relation can be transformed to a still simpler expression (the derivation shown in Supplement B, Sect.B3), particularly where f S /h R is the generalized sine divided by the metric scale factor of the coordinate R. According to ( 7), ( 14), (15), and also (A7) from Appendix A, we can write for f S /h R the analytical expressions in the small-ν region and in the large-ν region.Notice that f S /h R depends solely on W (not separatelly on R and ν).Therefore, f S /h R function is constant on the same surfaces in 3D where W is constant, i.e., on the cones with the tip at the origin and the axis coinciding with the rotation axis of the coordinate system.For an axial 2D section through the coordinate system, these cones become straight half-lines radiating from the center.For a better imagination, x-z plane map of the constant levels of f S /h R function is shown in Fig. 1.
As f S /h R appears to be a very important ratio in the Laplace equation solution, it does worth to derive some identities valid for this ratio.From the abovementioned derivation in Supplement B, Sect.B3, the following identity can be deduced:  153) and ( 157) in [10]) that the ratio f S /h R in a point is proportional to the ratio of the Cartesian coordinate z 3D and the equatorial radius R of the corresponding SOS coordinate spheroid surface for that point, i.e., For the spherical case, where μ 0, h R 1, f S sin (ν) and R meaning the radius of a sphere, this relation becomes trivial: sin (ν) z 3D /R.
We can also derive Cartesian coordinates x 3D , y 3D in 3D as a function of f S /h R : The derivation is done in Appendix A (see Eqs. (A25), (A26)).Therefore, the squared distance from the axis of rotation expressed in the SOS coordinate system is Finally, the squared position vector magnitude can be determined using ( 49) and (52) as where (17) was used for the rightmost expression derivation.When we return to Eq. ( 45), we see that the first determined non-trivial and non-constant solution of the angular part of the Laplace equation is The f S /h R function is here multiplied by the arbitrary constant c A1 , as such multiplication evidently also leads to a solution of (37).The index 1 denotes that this is the solution valid for the separation constant K d equal to 1, while the index A denotes that this is the first-kind solution of the expected two different kinds of solution (as the equation to be solved is the second-order differential equation).( 2024) 139:409 Note that for μ 0, f S /h R results in sin (ν).In this limit, the solution function is This function is basically the solution of the angular part of the Laplace equation in the spherical coordinates, particularly the first-degree Legendre polynomial (considering that we do not have polar coordinates-the origin of ν coordinate is here at the equator).It does worth to note that powers of sin(ν) form (within the Legendre polynomials) the angular part of the solution of the Laplace equation in the spherical coordinates.

Simple solution of the second kind
It has to be noted that (54) is only a particular solution of Eq. (38).Nevertheless, the general solution for K d 1 can be obtained as well.The general solution (for K d 1) is given by (see [17]) The first part, i.e., by c A1 , is a solution of its own, as it is the already found solution (54).The second part, i.e., the function F B1 (W ), is thus a new, separate, second-kind solution: To find it in an analytic form, we have to determine the integrals (inner and outer) in (57).The integration and then the derivation of the second-kind solution for K d 1 is done in Supplement B, Sect.B4 (the relations ( 14), ( 15), (17), and also (A75), (A1), (A2), (A3), (A4), (A5), (A7) from Appendix A, as well as Eq. ( 82) from [4] are employed in the supplement).The second possible real solution of the angular part of the separated Laplace equation for the given model (K d 1) is then derived in the form or, equivalently, with inverse hyperbolic sine instead of logarithm, (taking into account the well-known relation between the inverse hyperbolic sine and the logarithm functions).A correctness of this more complex solution (denoted "the second-kind solution" in what follows) of the full Laplace Eq. ( 24), in which the found radial part solution (33) with K d 1 is used, is successfully tested in Supplement B Sect. B5 (( 17), (10), (16), and also (A91), (A93), (A18), (A98), (A78), (A92), (A38) listed in Appendix A are employed there).It is confirmed that (58) (as well as (59)) is one of the solutions.

The found solutions rewritten to the form suitable for further generalization
The found second-kind solution for K d 1, i.e., (58), contains relations with f S /f C .Nevertheless, it can be written in terms of f S /h R as well, as (see Eq. (A24) in Appendix A) The power series expressions of s are recorded in (46) and (47).As h 2 R is equal to 1 on the equator and to 1/(1 + μ) on the polar axis (see Eq. ( 91) in [4] and employ the fact that h R is constant on straight lines going through the origin), and as f S is 0 on the equator and 1 on the polar axis (see [10], Eqs. ( 114) and ( 115)), the ratio f S /h R is 0 on the equator, while it is equal to √ 1 + μ on the polar axis, i.e., Further (see Eq. ( 91) in [4] and Eq. ( 123) in [10]), on the reference spheroid, s can be expressed in a closed form as With a help of (61), the second-kind solution (58) in terms of s is (see the derivation in Supplement C, part C1) or equivalently Note that the expression in the square brackets simplifies, for μ 0 (i.e., for the spherical case), to which is the first degree Legendre function of the second kind.
Similarly, the first determined solution of the angular part of the Laplace equation, i.e., the first-kind solution, is written in terms of s as which is (for μ 0) the Legendre function of the first kind (or Legendre polynomial) of the first degree (see Eq. ( 55)).It has to be noted that the range of the variable s is, in the SOS case, larger than in the spherical case (where the range is < -1, 1 >).According to (62), the upper limit of s is √ 1 + μ in the SOS case.

The angular part of the Laplace equation rewritten to the form suitable for further generalization
We now know that solutions exist resembling Legendre functions of the first degree and both of the first kind and of the second kind.Therefore, it does worth to test if the angular part of the Laplace equation in the SOS coordinates can be written in the form similar to the Legendre equation.Further step is thus to rewrite the angular part of the Laplace equation in the SOS coordinates (37) to the form containing only the variable s.The derivation is reported in Supplement C, part C2.With the help of the relations (A19), (A20) and (A31) from Appendix A, we arrive at the equation termed "generalized Legendre equation" in what follows.Indeed, for a special case when μ 0 (i.e., for the spherical case) the equation reduces to the well-known Legendre equation 10 Solutions for the special case K d 0 We already know the solution of (68) in the special case K d 1 (see Eq. ( 67) and ( 65)).We can relatively easily find the solutions in the case when the separation constant K d 0. (68) then simplifies to It is clear that one solution (the first-kind solution) is when F(s) equals to a constant, here denoted c A0 : The second-kind solution is more complicated.Nevertheless, if we use-as a test function-the part of the already found secondkind first-degree solution (64), particularly the logarithm-based function, we get its first and second derivatives as follows.The first derivative is , and the second derivative is The derivatives inserted to (70) result in the equation After cancellation of a number of terms in the equation, what remains is which is obviously fulfilled for any s.Then, the given function (72) multiplied by a constant, i.e., is the searched solution of the second kind for K d 0.

Polynomials as a solution of the angular part of the Laplace equation in the SOS coordinates
It was shown in the previous text that s f S /h R function is a basis for the found simple solutions of the angular part of the Laplace equation in the SOS coordinates for K d 0 and K d 1.Moreover, f S /h R function becomes sin (ν) in the limit μ 0. At the same time, it is known that powers of sin (ν) form (within the Legendre polynomials) the solution of the Laplace equation in the spherical coordinates.Therefore, a possibility is further tested, if the angular part of the Laplace equation in the SOS coordinates for various K d could be written generally in the form of Legendre-like polynomials, which employ f S /h R (i.e., the generalized sine divided by the metric scale factor of R coordinate) instead of the standard sine as the argument, i.e., Here, c j are unknown constants, and β 1 , β 2 , β 2 ≥ β 1 , are some integer numbers determining the number and power of the solution terms.This model, hopefully, will lead to the first-kind solutions when proper c j constants are found.
The angular part of the Laplace Eq. ( 68) can be-with this model-written in the form The first and second derivatives of F Aβ 2 (s) (which are easily acquired for a polynomial function), as well as the polynomial F Aβ 2 (s) itself, are inserted to (79).The equation is subsequently simplified.This is done in Supplement D, part D1.Then, the following equation is obtained: This is the angular part of the Laplace equation in the SOS coordinates when (78) polynomial is used as a solution model.It can be seen that the term with the pre-factor j( j − 1)(1 + μ) 3 in (80) will be always the term with the lowest power when j is equal to β 1 , and-moreover-it cannot be compensated by any other term in the sum as no other has that low power.Therefore, it has to be zero for j equal to β 1 .Similarly, this term cannot be compensated by any other term (as no other has that power of s) also forj equal to β 1 + 1.This restricts the bottom limit of the sum, β 1 , to 0. In such case (β 1 0), both abovementioned terms are zero thanks to the pre-factor j( j − 1).Thus, we select β 1 0 in what follows.
Further, it can be seen that the term with the fourth power in s will always give in the sum the term with the largest power in s f S /h R when j is equal to β 2 .Moreover, it cannot be compensated by any other term in the sum as no other has that large power.Therefore, the pre-factor has to be zero itself (in order the particular c β 2 could be nonzero).It can be easily shown that in order to fulfill , the roots have to be Moreover, it can be seen that the term with the fourth power in s in (80) will always give in the sum the term with the second largest power in s when j is equal to β 2 -1.It cannot be compensated by any other term as no other has that large power.It can be easily shown that in order to fulfill , the roots have to be β 2 K d ± 1, or c β 2 −1 has to be zero.As the β 2 roots for the abovementioned cases j β 2 and j β 2 -1 are mutually incompatible, it is clear that c β 2 −1 has to be zero, and that the upper limit of the sum, β 2 , is either equal to K d or to K d -2.We conservatively select the larger one, i.e., K d .Moreover, as the upper limit of the sum, β 2 , has to be an integer number, then K d has to be an integer, too.Further, it has to be positive or zero (i.e., K d ≥ 0), according to our condition β 2 ≥ β 1 .Thus, K d cannot be negative.This is a significant result, as it is different from the spherical case (i.e., when μ 0, also discussed in Supplement D part D1) for which negative constants K d are also allowed.It means that (78) model polynomial cannot lead to a solution for any negative separation constant K d when μ > 0. Further search for solutions with negative exponents K d in the radial part (i.e., in fact for the exterior points up to the infinity) has still to continue in the future for the spheroidal (SOS) case.
For integer K d ≥ 0, we further employ the above derived limits in the sum of (80), and we then have the basic equation for the model (78), i.e., for the first-kind solutions, in the form 12 Solutions of the first kind for the special cases K d 2, 3, 4, 5 and 6 We already know the solution in the form of polynomial for K d 0 and 1 (see Eqs. ( 71) and ( 67)).In Supplement D, part D2, the polynomial solutions of (83) for K d 2, 3, 4, 5 and 6 are derived, and they are reported here.For K d 2, the polynomial solution is while it is for K d 3. The polynomial solution for K d 4 is and the solution for K d 5 is Finally, the polynomial solution for K d 6 is Notice that-for μ 0-the derived polynomials reduce to the Legendre polynomials of degree 2 to 6.For the found polynomial solutions with various K d , normalization factors which would enable employing a proper recursion formula have still to be determined.Such formula is the topic of the next section.

Generalized Legendre polynomials and Bonnet-like recursion formula
Derivation of the higher-degree polynomial solutions of the angular part of the separated Laplace equation in the SOS coordinates (similarly as in the previous section) is increasingly more and more time demanding for the extent of calculations involved.Therefore, it would be of advantage to have a recursion formula which can derive relatively easily the higher-degree polynomial from a limited number of the lower-degree polynomials.For the spherical coordinates case, such recursion is known as Bonnet formula.Its derivation is based on a formal expansion of the following generating function in powers of t, where P n (x) are the Legendre polynomials.The Bonnet recursion formula for the spherical case Legendre polynomials is then For SOS coordinates, derivation based on the already known seven Legendre-like polynomials (for K d 0, 1, 2, 3, 4, 5 and 6, see the previous sections) results in the similar formula which reproduces the first seven previously derived SOS Legendre-like polynomials n K d , P SI n (s) ∼ F An (s) , except different normalization factors.These polynomials are summarized in Table 1.We will call them the generalized Legendre polynomials (or generalized Legendre functions of the first kind), and use for them the notation P SI n (s) where the upper index SI discriminates them from the standard Legendre polynomials P n (s) applicable in the spherical case.The "I" part of the upper index tells that these polynomials are valid only for the interior solution of the angular part of the separated Laplace equation.
For μ 0, it can be easily seen that the recursion formula (91) reduces to the well-known Bonnet formula (90) for the spherical case Legendre polynomials.
In Supplement D, part D3, the recursion formula (91) is tested.It is used to determine higher-degree polynomials (K d 2, 3, 4, 5 and 6) from the first two (K d 0, 1).The results of the use of the Bonnet-like recursion formula (91) correspond to the ones determined in the previous section.The calculation by-products are the correct normalization coefficients of the polynomials, which are then included into the polynomials recorded in Table 1.In Table 1, therefore, the generalized Legendre polynomials P n SI for the SOS coordinates with the proper normalization factors for the interior first-kind solution of the generalized Legendre equation up to 6th degree are listed.
It can be seen from Table 1 that for μ 0, the SOS case polynomial solutions are the same as for the spherical case (as s f S /h R is then equal to sin ν sin χ).It can be also seen that the Legendre-like polynomials P n SI are finite for any nonnegative μ.When the generalized Legendre polynomials P n SI (s) are compared (see Table 1) with the Legendre polynomials P n (s) applicable for the spherical case (μ 0), we see that there are two main differences.First, the factors by the powers of s are no more simple integers k 0 , but they are more complex factors containing μ, which are (without a mathematical rigor) of the form (k 0 + k 1 μ + k 2 μ 2 + …)(1 + μ) l .For μ 0, these factors result in the same integer number as for the spherical case Legendre polynomials, as expected.Further, for μ 0, a normalization factor exists differing from the spherical case one by 1/(1 + μ) n for each generalized The substitution s f S /h R is used.All the polynomials of degree 0 to 6 were derived by a solution of the generalized Legendre equation.The polynomials of degree 2 to 6 were also checked by the Bonnet-like recursion formula (91).Comparison with the classical Legendre polynomials is provided in the last column.χ in the last column is latitude 123 Legendre polynomial.In the limit of the spherical case (μ → 0), the factor 1/(1 + μ) n becomes equal to one for all the polynomials.For the spheroidal case (μ > 0), nevertheless, it differs from 1 and, therefore, is essential when employing the recursion formula.
This above analysis indicates (although it is not a full rigorous proof for an arbitrary polynomial degree) that the generalized Legendre polynomials can be derived by the use of the recursion formula (91) also for the higher degrees than those listed in Table 1.It also indicates (although it is not proven rigorously) that the generalized Legendre polynomials determined in this way are solutions of the generalized Legendre equation (which is equivalent to the angular part of the Laplace equation in the SOS coordinates).
While the Legendre polynomials are orthogonal (with respect to the weighting function equal to 1), simple integration tests show that the generalized Legendre polynomials are not orthogonal with respect to the same weighting function.It is, therefore, more complicated to find proper coefficients for a solution given by particular boundary conditions.This is, nevertheless, out of the scope of this paper, and we leave finding of a proper weighting function on experienced mathematicians.(We only note that the known procedures used for Bessel functions weighting [18,19] do not work as the K d separation constant is present in (68) equation also in the factor by the first derivative, not only by the function itself.)

The second-kind solutions of the second degree and of the third degree
Let us title the second solution of the generalized Legendre equation (or equivalently of the angular part of the Laplace equation in the SOS coordinates) the generalized Legendre function of the second kind.We will denote it Q n S (s), where n denotes degree.Two second-kind solutions were previously found: for K d 0 (see Eq. ( 77)), to which Q 0 S (s) is proportional, i.e., and for K d 1 (see Eq. ( 64)), to which The above relations for the generalized Legendre functions of the second kind are with the normalization factor set in such a way that the initial part of the relations equals to P 0 SI (s)Q 0 SI (s) and P 1 SI (s)Q 0 SI (s).According to (62), the upper limit (i.e., on the pole) of s is √ 1 + μ.When assessing the denominator value in the logarithm function, it is seen that the generalized Legendre functions of the second kind diverge on the axis of rotation, similarly as in the spherical case.Therefore, the second-kind solution is probably of less applicability in physics than the first-kind solution.Nevertheless, as there could be some special applications also for the second-kind solution in case of singularity, e.g., when solving the Laplace equation in Schwarzschild spacetime [20], it does worth to derive the higher-degree second-kind solutions of the generalized Legendre equation in the SOS case as well.
The second degree solution of the second kind was found by a tedious way of testing a model function.The derivation is done in Supplement E, part E1.The result is as follows: It can be easily checked that-in the special case when μ 0-it would simplify to the classical second-kind Legendre function of the second degree.
For further test of a prospective recursion procedure, it does worth to have also the third degree solution.Extensive derivation and algebraic manipulations shown in detail in Supplement E, part E2, lead to the following formula for the solution: When μ 0, then the function simplifies to the classical second-kind Legendre function of the third degree, as can be easily found by the substitution μ 0: 15 The second-kind solutions The already found second-kind solutions of the degree 0 to 3 are listed in Table 2.In Table 2, there are-for comparison-also listed the classical Legendre functions of the second kind for the spherical case up to the degree six.It has to be noted that standardly, these functions are reported only up to the degree 3. Therefore, Q 4 (x) and Q 5 (x) were obtained as the corresponding coefficient of expansion of the generating function A generating function different from (97) which can help to determine the higher-degree Legendre functions of the second-kind Q n (x) is based on the generating function suggested in [21].U n (x) functions are here the polynomials forming the second part (i.e., the one not by the logarithm) of the Legendre functions of the second kind, and thus The classical Legendre functions of the second-kind Q 5 (x) and Q 6 (x) in Table 2 were obtained using (99).
It is known that for the spherical case solution (μ 0), the Bonnet recursion formula (90) is valid also for the second-kind solutions.With its help, higher-degree second-kind Legendre functions can be relatively easily derived when μ 0. It indicates that the Bonnet-like formula (91) could similarly work also for the second-kind generalized Legendre functions, i.e., for the SOS coordinates case.
It can be deduced by comparing P n SI (s) in Table 1 with Eqs. ( 92), ( 93), ( 94) and ( 95) that the following relation holds for the second-kind functions: where T n (s) is a polynomial (especially, T 0 (s) 0, T 1 (s) 1/(1 + μ)).The first term in (100) has identical form as for the standard Legendre functions (see e.g., [22,23]), and the second term contains special function (1 + μ) 2 − μs 2 which equals to 1 for μ 0. In Supplement E part E3, a test is performed if the recursion formula (91) can be used (similarly as for the standard Legendre functions) also for the generalized Legendre function of the second kind, i.e., and which normalization coefficients are to be used.The derivation confirmed that the recursion formula reproduces the found solutions for the second (Eq.( 94)) and third (Eq.( 95)) degree function.The functions with a proper normalization coefficients for the use in the recursion formula are reported in Table 2. Further, (101) is used in Supplement E, part E3, also for determination of the second-kind functions Q 4 SI (s), Q 5 SI (s) andQ 6 SI (s) which are then also included into Table 2.It can be observed that Q n SI (s) Q n (x) when μ 0 in all cases.
As the formula for Q n SI (s) determination contains in the pre-factor by logarithm (i.e., by Q 0 SI (s)) the term P n SI (s), which is the first-kind solution of the generalized Legendre equation, the range of the separation constant K d for which the second-kind solutions exist is the same as for the first-kind solutions.Therefore, neither the second-kind solution (similarly as the first-kind solutions) derived in the form shown above does exist for negative separation constants K d .

The complete solution in interior points
Finally, the complete interior solution of the azimuthally symmetric case of the Laplace equation in orthogonal similar oblate spheroidal coordinates (20) with applied variables separation introduced by ( 23) can be written in the form where the radial part solution (33) is used, and P n SI (s) and Q n SI (s) are the generalized Legendre functions of the first and of the second kind, respectively.Several low-index functions can be seen in Tables 1 and 2. Note that the zero degree solution a 0 R 0 P SI 0 (s) is equal to the constant V C (see Eq. ( 23)). a n and b n are arbitrary real coefficients, which can be determined for the particular boundary conditions.(102) thus represents the harmonic functions in the SOS coordinates.However, as shown in the previous text, this solution exists only for the nonnegative integer values of the separation constant K d n.

Conclusions
The interior solution of the azimuthally symmetric case of the Laplace equation in the orthogonal similar oblate spheroidal coordinates (20) was found (see Eq. ( 102)).During the derivation, the Laplace equation was separated by a special separation procedure to the radial (see Eqs. ( 31), ( 32)) and the so-called angular (see Eq. ( 36)) parts.While the solution (33) of the radial part is simple power of the equatorial radius of the oblate spheroidal coordinate surfaces, the separated angular part of the Laplace equation is significantly more complicated and leads to the generalized Legendre Eq. ( 68).Its solution then leads to the generalized Legendre functions of the first and of the second kind with argument s f S /h R , i.e., the generalized sine divided by the metric scale factor of the coordinate R. Analytical expressions with infinite power series involving generalized binomial coefficients were written for this argument (see Eqs. ( 46) and ( 47)).
Several lower-degree generalized Legendre functions of both kinds were derived and are listed in Tables 1 and 2. The recursion formula for the generalized Legendre polynomials (see Eq. ( 91)) was found.The general solution of the rotationally symmetric Laplace equation for the interior space is reported (see Eq. ( 102)).
The presented solution is a step forward in the establishment of the SOS coordinates as a useful tool in the field physics.Further, the determination of the harmonic functions could facilitate finding solutions of more complex differential equations in the SOS coordinates, for example the Poisson equation or the Helmholtz equation.
Within the Laplace solution determination process, relevant identities valid for the SOS coordinate system in general and for the generalized sine and cosine functions in particular were derived as well (see Appendix A).
The found solution exists only for the nonnegative integer values of the exponent of the equatorial radius in the solution.Therefore, the solution has physical meaning only for the interior space (i.e., it is restricted to a certain annulus) and cannot be used for a region extended up to the infinity, as the solution would diverge.In order to find solutions in the exterior points, another approach has to be used.The exterior points solution is in fact already partially known to the author.However, due to the extent of the present text, it will be published elsewhere together with other considerations on the Laplace equation solution in the SOS coordinates.A2.1.Jacobian divided by the square of the h ν scale factor for small-ν region.A2.2.Relations between metric scale factor h R and the generalized sine and cosine.A2.3.Relations between W , h R , f S and f C .A2.4.Relations involving h ν metric scale factor.A3.Derivatives of W and F(W ) with respect to ν and R. A3.1.Derivatives of W with respect to ν and R. A3.2.Derivatives of F(W ) with respect to W and ν.A4.Derivatives of the Pólya-Szegő-type power series.A4.1.Derivative of S C type of Pólya-Szegő power series.A4.2.Derivative of h 2 R .A4.3.Derivative of S A type of Pólya-Szegő power series.A4.4.Summary for the derivative of S C and S A type of Pólya-Szegő power series.A4.5.Derivatives of power series contained in metric scale factors.A4.6.Derivatives of the generalized cosine and sine.A4.7.Derivatives of functions containing Jacobian.A5.Integral identities involving the Pólya-Szegő-type power series.

A2. Basic operations with the expressions containing Pólya-Szegő-type power series
A2.1.Jacobian divided by the square of the h ν metric scale factor For the use in the Laplace equation solution, we derive here the Jacobian divided by the square of the h ν scale factor, first in the small-ν region.This quantity is needed in the case of Laplacian calculation in the SOS coordinates.We start with a division of (11) by the squared (8): For the determination of the ratio of the power series in (A6), the combination of the Hagen-Rothe identity with the Cauchy product (see Eq. ( 50) in [10]) is to be used.Here, ε 2, b -μ, c -(μ + 2) and a (μ + 1)/2, thus a + c -(μ + 3)/2.Then / h 2 ν ratio derivation for the large-ν region is carried out in a similar way by using (A7

. Relations between metric scale factor h R and the generalized sine and cosine
Several relations between the metric scale factors and the generalized sine and cosine (see Eqs. ( 7) and ( 14), ( 15)) are needed for the Laplace equation solution.An important relation is 1 − h 2 R /μ, which can be-for the small-ν region-written as where we used the substitution M k-1.When we further use (subsequently) two well-known binomial identities (see, e.g., [10], Eqs. ( 53) and ( 54)), we arrive at If we use the notation coming from (14), which introduced the generalized sine f S , in the last part of (A12), we obtain the relation Similarly, for the large-ν region, we can write where we used the substitution M k-1.When we further use (subsequently) the binomial identities listed in (A11), we get Using the definitions of the generalized cosine ( 14) as well as (16), we arrive at ) and (A16) are formally the same.Therefore, the same relation can be used in both the small-ν region and the large-ν region.
For completeness, we report also the expression (derived with the help of Eq. ( 91) of [4]) valid at the reference surface (i.e., for R R 0 ): By using ( 16) and ( 17) relations, the following relations valid in both regions can be easily derived between h R and the generalized sine and cosine: Further useful relations between the metric scale factor h R and and f S /h R or f C /h R are as follows: Thanks to (A19), we also obtain the following relation between h R and f S /h R : Further, Another useful relation comes from (17): A very important relation is the one between f S /f C and f S /h R .The derivation employs (A19).We obtain from ( 17) (1 + μ) Thus, We can also derive Cartesian coordinates in 3D as function of f S /h R .With the help of Eq. ( 47) of [4], of the relations (A2), (A7) and (A19), as well as of ( 14) and ( 16) from the main text of the present article, we can write squared x 3D in the small-ν region as Accordingly, squared y 3D can be derived as

A2.3. Relations between W, h R , f S and f C
When we look at the form of the generalized sine and cosine (14), (15) and use (A2), (A7), we see that-in the small-ν region-we have In the large-ν region, we obtain in a similar way It does worth to note that this expression works also for μ 0, as √ 1 + μ 1 μ has the limit for μ → 0 which is square root of e.In this limit, the expression is thus e sin ν.Therefore, we can write an important relation between the generalized trigonometric functions, the metric scale factor h R and the parameter W , valid in both regions, as The combination of (A24) with (A29) leads then to a possibility to express W in terms of f S /h R : Then (using the definition s≡f S /h R ), It can be also easily found by substituting (17) to (A31) that we can express W 2 using only h R , which, with a subsequent use of ( 17), leads to The W definition (5), and a subsequent use of the binomial theorem, results in the identity A2.4.Relations involving h ν metric scale factor We will also need relations involving h ν metric scale factor.For the further derivations of derivatives, it does worth to define the relations with power series contained in the h ν metric scale factor (8) by where the left side relation holds for the small-ν region, and the right side relation for the large-ν region.Such definition makes the following derivations easier to overview and having lower number of multiplication factors.It also does worth to notice that S ν depends only on W , not separately on R and ν.Then, h ν metric scale factor (8) can be written as in both the small-ν region and the large-ν region.That means the expression is formally the same for the small-and the large-ν regions.
With the help of ( 7), (A35), of Jensen's identity from [10] (particularly Eq. ( 52) with t 1) and of the relations for the generalized sine and cosine ( 14), ( 15)), the following identities can be derived in the small-ν region: With the help of the same relations but using t 1/(1 + μ)) instead, we arrive in the large-ν region to the same relation as the one valid in the small-ν region.Then, the following relations in this sub-section are valid for both regions.When applying (A36) onto the left side of (A37), a very useful relation between the metric scale factors and the generalized sine and cosine arises: Further, with a help of ( 16) and ( 17), we get from (A37) When we use the relation (A36) between S ν and h ν , (A39) leads to the following relation between the metric scale factor h ν and powers of h R : which is a very significant relation for the solution of the Laplace equation.

A3 Derivatives of W and F(W) with respect to ν and R
A3.1.Derivatives of W with respect to ν and R The parameter W (see Eq. ( 5)) and its derivatives (see Eqs. ( 9), (10)) are reported in the main text.On that basis, the second derivative with respect to ν is calculated as For completeness, the second derivative with respect to R is A3.2.Derivatives of F(W ) with respect to W and ν With a help of ( 9), the first derivative of a function F depending on W with respect to ν is and the derivative with respect to W is thus The second derivative with to W is Then (with a help of (A41)), Therefore, the second derivative of F(W ) with respect to W is finally

A4 Derivatives of the Pólya-Szegő-type power series
An essential element of differential calculations in this article is the knowledge of the derivatives of the Pólya-Szegő-type power series present in the relations (A1), (A2) of this Appendix, with respect to W .These series, are, in what follows, restricted (see the conditions (A4) and (A5)) to the cases with b -μ, ε 2, and with b μ/(1 + μ), ε -2/(1 + μ), respectively, as these correspond to the series used in the small-ν region and in the large-ν region, respectively, in the SOS coordinate system.Therefore, they are particularly interesting for this article.We define the following expressions and symbols for the two types of the series with generalized binomial coefficients: and where the left side relations hold for the small-ν region, and the right side relations for the large-ν region.These expressions cover all the expressions with the power series employing generalized binomial coefficients which are used in the SOS coordinate system.Just a parameter has to be set differently for the various elements of the system (relations for the coordinates transformation, the metric scale factors, the Jacobian, see [4], and the generalized sine and cosine, see [10]).Using (A48), (A49), (A7) and ( 7), it can be seen that in both regions.
A4.1.Derivative of S C type of Pólya-Szegő power series We can write for both cases (i.e., for the small-ν region as well as for the large-ν region) where b -μ, ε 2, β 0in the small-ν region, and b μ/(1 + μ), ε -2/(1 + μ), β 2a in the large-ν region.Then, the derivative can be relatively easily obtained by multiplication and division by a proper power of W followed by the product rule use: Further simple algebraic manipulation leads to When the parameters for the small-ν region are input, we get By applying (A50) onto the second term, we get in the small-ν region, where the rightmost expression comes from the relations between f S and h R (17).
When the parameters for the large-ν region are plugged into (A53), we get By applying (A50) on the second term, we obtain in the large-ν region, where the rightmost expression comes from the relations between f S and h R (17).

A4.2. Derivative of h R
More complicated is the derivative of the second-type power S A (see Eq. (A49)).For this purpose, we will need to find first the derivative of h R 2 with respect to W .We can compare Eqs.(F19) and (F20) of Supplement F of [4] to determine a helpful identity applicable in the small-ν region Similarly, by comparison of Eqs.(F39) and (F40) of Supplement F of [10], another helpful identity applicable in the large-ν region is determined Both these identities lead, in the respective regions, to the identity (see Eq. (F24) and (F44) of Supplement F of [4]) which is valid in both regions (small-ν, large-ν) when the proper series for the functions h R , h ν (see Eqs. ( 7), ( 8)) are used for the respective regions.Then (see Eqs. ( 8) and (A35)) in the small-ν region, while in the large-ν region.A correct relation for S ν according to (A35) has to be used for the particular region.Then, the formal expression for the derivative is the same as in the small-ν region.The final formulae in (A61), (A62) have an advantage with respect to the initial formulae that they do not contain explicitly R, which is a correct result (thanks to the h ν dependence on R).The given formulae for the derivative dh 2 R dW thus does not depend explicitly on R (it depends directly only on W ). (The use of h ν symbol and division by R 2 -as is the case in the initial formula-would be rather misleading).

A4.3. Derivative of S A type of Pólya-Szegő power series
Now, we have the tools for derivation of the S A power series with respect to W in terms of the already known functions.
Thanks to the identity (A50), we can write A4.5.Derivatives of power series contained in the metric scale factors For special values of the parameter a in the derivatives of S A expressions, we get derivatives of the metric scale factors.The first interim relation provides derivative of S ν (see Eq. (A35), connected with h ν metric scale factor by (A36)).It can be calculated from the S A derivative (see Eqs. (A65), (A69)) taking into account that S ν S A for a -(μ + 2) in the small-ν region, and for a -(μ + 2)/(1 + μ) in the large-ν region.Input of the parameter a for the small-ν region leads (using (A65)) to The same action in the large-ν region leads-using (A69)-to As can be seen, the relations are the same for the small and the large-ν regions when proper relations for S ν and h R are used in the respective region.
The second set of the relations should provide derivative of h R .Nevertheless, these were already derived above in (A61) and (A62) with the following result: Again, these relations are the same for the small-and the large-ν region when the proper relations for S ν and h R are used in the respective region.With a help of (A39), it can be also written in the following forms: valid also in both the small-and the large-ν region.The last expression is obtained with the help of the relation (A22).Further, dh R dW More generally, we have valid for any integer j.

A4.6. Derivatives of the generalized cosine and sine
We can also derive the derivatives of the squared generalized cosine f C and sine f S .With a help of ( 14) and (A66), S A according to the definition (A49) and with a -1, we get in the small-ν region where the last formula follows from algebraic manipulation using identities (A61) and (A37).In the large-ν region, by comparing ( 14) and (A49) with a -1/(1 + μ)), and by using (A70), The relations are formally the same in both regions.Thanks to the (16) relation, the derivative of the square of the second generalized trigonometric function, f S , with respect to W is equal but of the opposite sign than the derivative of the squared generalized cosine f C : Also for this generalized trigonometric function, the relation is valid for both the small-and the large-ν region.This simplifies all further derivations containing derivatives of the generalized sine and cosine, which can be done with the same notation in both regions.
It should be noted that the above derivatives are carried out with respect to the parameter W , not with respect to the coordinate ν.It means that when we want to compare limiting values of derivatives of generalized trigonometric functions with the derivatives of the normal trigonometric functions, we have to calculate the derivative with respect to ν first, and only afterward to set μ 0. The chain rule tells that where ( 9) was used for the partial derivative.For the limiting value μ 0 corresponding to the spherical case, we have (as h R 1, f C being equal to cosine, f S being equal to sine in such case) ∂ f 2 S ∂ν 2.1.sin 2 ν cos 2 ν 1 + 0. sin 2 ν sin ν cos ν 2 sin ν cos ν (A82) which is equal to the derivative of the squared sine with respect to ν.
The derivatives of the generalized cosine and sine (i.e., not squared) are and Still more generally, the derivatives of the powers of the generalized sine and cosine are Other important derivatives of the special composed functions containing generalized sine or cosine are as follows: Further useful derivative is the derivative of the ratio of generalized sine and cosine: and also derivative of its square: The reversed ratio of generalized goniometric functions has derivative Finally, a very important derivative is the one of the ratio of the generalized sine f S and the metric scale factor h R .First, let us find the derivative of the square of it.With the help of ( 17) and (A75), we easily obtain Thus and generally (also with the help of ( 17)) A remarkable property of (f S /h R ) m derivative is that the derivative (and the second derivative) gives the same function (f S /h R ) m multiplied by m and m-independent function.Remind also that the above derivatives are with respect to W , not ν.
The power m generalized cosine over metric scale factor derivative can be derived accordingly:

. Derivatives of functions containing Jacobian
For the Laplace equation solution, several derivatives of the expressions containing Jacobian are needed.First, we derive the derivative of the Jacobian divided by the square of the h R divided by ∂ W ∂ν for the small-ν region with respect to W .With the help of (12), and using (A55), i.e., the derivative of S C type of power series, we can derive the following: Further, we will need the derivative of the Jacobian divided by the square of the h ν , multiplied by ∂ W ∂ν for the small-ν region with respect to W .With the help of (13) and using (A55), we arrive at The above derivative with respect to the coordinate ν is A5 Integral identities involving the Pólya-Szegő-type power series The integral can be expressed in the following form: Now, we can employ the basic binomial identity (see [12]) and arrive at (we skip for the moment the integration constant and attach it again at the end of derivation only)

Fig. 1
Fig.1Constant levels of f S /h R function (the generalized sine divided by the metric scale factor h R ) in one quadrant of x-z plane for the SOS coordinate system with μ 2. The set of the orthogonal SOS coordinate system lines is used here as a background for a better overview

Table 1
Table of the first several lower-degree generalized Legendre polynomials P for SOS coordinates with the proper normalization factors for the first-kind solution of the generalized Legendre equation (i.e., the angular part of the separated Laplace equation in the SOS coordinates) up to the sixth degree n SI (s)

Table 2
Generalized Legendre functions of the second-kind Q n SI applicable for the solution of the Laplace equation in the SOS coordinates up to the degree six with proper normalization factors

Table 1
S /h R is used.Comparison with the Legendre functions of the second-kind Q n (the last column) is done as well.The polynomials P n SI (s) and P n (x) appearing in the table can be found in