Explicit Equilibrium Solutions For the Aggregation Equation with Power-Law Potentials

Despite their wide presence in various models in the study of collective behaviors, explicit swarming patterns are difficult to obtain. In this paper, special stationary solutions of the aggregation equation with power-law kernels are constructed by inverting Fredholm integral operators or by employing certain integral identities. These solutions are expected to be the global energy stable equilibria and to characterize the generic behaviors of stationary solutions for more general interactions.


Introduction
The mathematical analysis of equilibria and traveling wave type patterns in many body descriptions of collective behavior models is one of the basic scientific problems in the explanation of coherent structures in mathematical biology and technology. They naturally appear in swarming of animal species, cell movement by chemotaxis, granular media interaction and selfassembly of particles, see for instance [33,16,31] and the references therein. Consensus in orientation patterns have been reported in the seminal papers [37,25] in which they numerically show the asymptotic stability of compactly supported groups of particles moving in a fixed direction. The amazing feature of these flock patterns is that the consensus in orientation can be established based only on attraction and repulsion effects. More precisely, the movement of each particle is assumed to follow the equations of motion (1) d dt where f is a function modeling the self-propulsion at low speed and friction at high speed, and K is the radial pairwise interaction potential. As the total number of individuals is large, the system of differential equations is difficult to analyse and usually a continuum description based on mean-field limits is adopted, either at the kinetic level for the particle distribution function [15,16] or at the hydrodynamic level for the macroscopic density and velocities [24,15]. Different stationary patterns were observed for the system at the equilibrium speed s 0 (with f (s 0 ) = 0), for instance coherent moving flocks and single or double rotating mills [37,25,15,19,20]. The stability of these patterns at the discrete particle level was analysed in [7,2,18]. At the continuum level, these patterns are characterized by searching for continuous probability densities or probability measures ρ of particle locations such that the total force acting on each individual balances out. This is equivalent to finding probability densities or measures ρ such that (2) ∇K * ρ = 0 on supp(ρ) .
Being the problem posed on the support of the unknown density ρ implies that the equation (2) is highly nonlinear. In fact, characterizing the interaction potentials K such that these profiles are continuous or regular in their support has recently triggered lots of attention. The richness of the qualitative properties of particular steady solutions depending on the potential K was reported in [34]. They numerically investigate the shape of the most stable solution for the first order dynamics associated to (1) given by   , i = 1, · · · , N , and they characterize several bifurcations starting from the particular solution of the Delta ring, a solution concentrated on a circle. However, getting conclusions about the qualitative properties of the continuum problem (2) turned out to be mathematically involved. The existence of explicit formulas for flock and mill patterns for particular potentials such as the Morse potentials, originally used in [25], was possible due to the particular properties of associated differential operators [36,6,21,17]. Let us point out that mill patterns are also characterized by a probability density or measure satisfying a similar equation to (2) giving the balance between attractive-repulsive and centrifugal forces of the form ∇K * ρ = −s 0 ∇ ln |x| in supp(ρ). The continuum evolution model associated to the particle system (3) via the mean-field limit [12] is given by where ρ(t, x) is the mass density function or measure of particles/individuals in space at time t ≥ 0. This equation is known as the aggregation equation and it has been proposed for swarming modeling in [38] as the inertia-less continuous system associated to (1). This model has been thoroughly studied with attractive and repulsive-attractive potentials, see [12,23] and the references therein. Depending on the singularity of the potential at the origin, the solutions enjoy different regularity properties or propagate certain regularity in time (even sometimes uniformly in time). All these early results imply that the regularity of the stationary states depends on the singularity or the strength of the repulsion of the potential K at the origin. In other words, the stationary states' regularity is deeply connected with the dissipation properties of the evolution (4). In fact, this equation has a gradient flow structure related to probability measures [22,41,14,4], whose implication for our problem here is that stable stationary states of the dynamics of (4) should be among the local minimizers in suitable topology of the interaction energy functional defined over all positive measures of mass M 0 . Actually, it is shown [3] that the local minimizers ρ of the interaction energy (5) in suitable transport distance topology have to satisfy the Euler-Lagrange equations: . The previous conditions have to be understood ρ-a.e for the first condition and a.e. Lebesgue for the second one in case ρ/M 0 is just a mere probability measure instead of a density function. Using these Euler-Lagrange conditions, it was shown that the dimension of the support of local minimizers depends strongly on the singularity of the potential at origin [3]; that local minimizers can be bounded, compactly supported or continuous on the support if the singularity is strong enough [13]. The last point was shown by using a hidden connection discovered in [13] between the Euler-Lagrange conditions (6) and the classical obstacle problem for elliptic operators. Finally, the existence of global minimizers of the interaction energy has been obtained just recently under quite sharp conditions [9,40]. Therefore, we know the existence of solutions to the Euler-Lagrange equations (6) for a quite large family of potentials including the Morse potentials under the non H-stability condition [25,9] and all power law potentials Let us remind the reader that the convention |x| 0 0 := log(x) is used. Because of the simpler topology, the one-dimensional case is in general better understood, see [27,28,11] and the references therein.
In this paper, we consider certain special exact stationary solutions in the sense that they satisfy the first condition of the Euler-Lagrange equations (6) K * ρ = E on the support of density ρ for certain ranges of the power-law potentials. In our examples, we will find compactly supported, radially symmetric densities ρ which are continuous inside its support. The power-law potentials might be considered a very restrictive case, however they are the typical examples obtained by asymptotics expansions and they should show the generic boundary behavior near the support of the stationary states. They also allow direct connections to more classical problems involving the Laplacian or its fractional counterparts. In fact, certain range of interactions are intimately connected to fractional diffusions [23,32]. Finally, having explicit solutions allows for excellent test cases for numerical schemes [8,10].
In this paper, we will always look for radial symmetric solutions associated with the power-law kernel (7), which are supported on a ball B R and are continuous inside. Explicit solutions are constructed in the range of parameters corresponding to −d < b < a, b ≤ 2 when either a or b is an even integer. These solutions are derived from the condition K * ρ = E on the support of density ρ with mass M 0 , and are reasonable candidate to be global minimizers of the interaction energy. The paper is organized as follows. Some basic integral equations and their solutions in one dimension are discussed in Section 2, but detailed derivations are deferred in in the Appendix. The one dimensional case with either a = 2 or b = 2 is studied thoroughly in Section 3, where some key observations allow a unified yet simpler approach in multidimensional case treated in Section 4. We finally conclude discussing possible extensions, and comment on certain open problems and conjectures of the generic behavior of solutions to (6).

Fredholm integral equations and integral identities
To derive the exact steady states governed by (6) with either a or b being an even integer, the key is to invert the Fredholm integral operator of the form that is, to find the solution of L[ρ] = f . In one dimension, singular integral equations with power-law kernel like this have been studied extensively [26], usually in connection with the theory of complex variables. The first step in solving L[ρ](x) = f (x) is to differentiate both sides of the equation such that the kernel |x − y| p becomes weakly singular, with exponent p between −1 and 0. Depending on the precise value of p, we have the following two cases.
If p ∈ (2k −1, 2k) for some non-negative integer k, upon taking derivatives 2k times, the integral equation becomes with ν = 2k − p ∈ (0, 1). The solution is given by (see [26] for the derivation) including a Cauchy principal value integral in the second term. Explicit solutions can be obtained for special right hand side f like polynomials. For example, if f (x) = 1, The detailed derivation, mainly involving the evaluation of the principal value integrals, is given in Appendix B.
If p ∈ (2k, 2k +1) for some non-negative integer k, upon taking derivatives 2k + 1 times, the integral equation becomes with ν = 2k + 1 − p ∈ (0, 1). The solution is given by (see [26] for the derivation) In addition to similar principal value integrals as in (9), the solution consists of a non-trivial null space spanned by the function (R 2 − x 2 ) ν 2 −1 . Explicit solutions can also be obtained when f is a polynomial. For example, if f (x) = x, then (14) ρ Besides the appearance of solutions to (8) or (12), certain integral identities have also to be established when the in order to eliminate free parameters in the final steady states. In many cases, the energy has to be calculated, to select the unique stable state. These special integrals turn out to be identities intimately connected to the solutions given above. For example, from the solution (10) for f (x) = 1, we obtain the identity valid for ν ∈ (−1, 1) and x ∈ (−R, R). Other identities with higher exponents can be obtained by successive integrations. For example, integrating twice on both sides of (15) gives where the integration constant (by evaluating at x = 0) is represented as special integrals from Appendix A. Since the exponent of the kernel |x − y| is usually fixed in the problem, it is more convenient to change ν to ν + 2 in the previous identity to have the same interaction kernel |x − y| −ν , that is, which is valid for ν ∈ (−3, −1). In fact, this identity holds for a wider range of parameter ν ∈ (−3, 1), precisely when the integral on the right hand side of (16) is well-defined. As we will see, these identities are more useful than the general solutions formula (9) and (13), because the former can provide a unified approach for the construction of the steady states, regardless of the range of the parameters. When p is an integer, the solution for the integral equation (8) has to be derived differently, but usually can be established equivalently from the above results by taking the limit. For this reason, the special cases of integer exponents are not discussed separately in the following sections.

Exact steady states in one dimension
In this section, we derive explicit expressions of the steady states in one dimension from the governing equation when either a or b is 2. These compactly supported symmetric steady states exist only in the parameter regime Figure 1: the potential integral is well-defined only for b > −1, but the steady states becomes two delta masses as b exceeds 2 (see [4,27,28]). Upon repeated differentiation on both sides of (GE), different integral equations of the form (8) or (12) appear, as summarized in Figure 1. In the special cases when either a or b is 2 of our interests, only the solutions (10), (11) and (14) are employed below. The solutions are assumed to be radially symmetric and supported on some ball B R with unknown radius R. The precise expressions are derived first, and then the radius R is determined from various selfconsistency conditions like the definitions of certain moments. Without loss of general, we take the total mass M 0 given. 3.1. The case when b = 2, a > 2. We focus on the case a ∈ (2, 3) first. Taking derivative three times to (GE), we get the homogeneous equation From (13), the solution is given by (17) ρ To determine the parameter c, we first use the definition of the total mass M 0 , that is, On the other hand, certain information is lost during the differentiation and has to be retrieved upon substituting the solution (17) back into (GE). By the special integral (42) in Appendix A, we infer that From the identity (16) Therefore, the left hand side of the governing equation (GE) is reduced to cos aπ/2 π .
For this expression to be a constant, the coefficient of x 2 should vanish, leading to the condition M 0 = c(a − 1). This, together with (18), can be simplified to determine the radius of the support and therefore the solution (17).
It is also instructive to calculate the associated energy in this case, because the parameter b = 2, a ∈ (2, 3) lies on the boundary separating compactly supported solutions that we are seeking and solutions consisting of two delta masses [3,7]. From (19) and (20), the energy can be simplified as On the other hand, the more singular solutions consisting of two delta masses can be written as Here the radius R 0 can be calculated in different ways, for example by minimizing the total energy This optimization procedure gives R 0 = 1/2 with the total energy E δ = (2 − a)M 0 /4a. Though difficult to compare analytically, numerical calculation shown in Figure 2(a) implies that E given by (21) is always less than E δ above, indicating the solution given by (25) is more stable and it is believed to be the global stable solution in the sense of the global minimizer of the interaction energy (5), see [9]. The convolution K * ρ, shown in Figure 2(b) for a = 5/2, is constant on [−R, R] (the solid line between two dots) and is larger outside the support. In contrast, K * ρ δ for a = 5/2 behaves differently: although R 0 is the optimal radius, the fact K * ρ δ has a local maximum at ±R 0 implies that the energy can be reduced by making the density less concentrated. We can proceed in the same way to find solutions for the other case when a ≥ 3 and b = 2, by differentiating the governing equation (GE). If a ∈ (2k − 1, 2k) for some integer k ≥ 2, the resulting homogeneous solution is of the form (8) and has only zero solution. If a ∈ (2k, 2k + 1) for k ≥ 2, the solution ρ(x) is proportional to (R 2 − x 2 ) (2k−1−a)/2 , whose coefficient must be zero to make sure no higher order terms appear in (GE). In fact, densities like (22) consisting two delta masses are the only stable steady states in this case.
3.2. The case when a = 2, −1 < b < 2. For a = 2, the governing equation becomes Depending on the value of b, derivatives of different orders have to be taken.  (17) compared with E δ from two Dirac masses; (b) The function K * ρ and K * ρ δ .
If b ∈ (1, 2), we take derivative twice on (23) to get From (10) with ν = 2 − b, the solution is given by Using the definition of the total mass the radius of the support is then given by If b ∈ (0, 1), we take derivative once on (23) to get From (14) with ν = 1 − b, the solution is can be written as (28) ρ where c is the free parameter. Because the component ( near the boundary |x| = R, we must have the constraint c ≥ 0 for the nonnegativity of the solution. The definition for the total mass is equivalent to from which c and R are related. Therefore, we have a one family of steady states parameterized by c ≥ 0. To select the most stable one of the family, we have to calculate the energy again as a function of c. The integration with the quadratic potential is straightforward, Using (15) and (16) Collecting all the terms above, the energy can be written as Using the constraint (29), the energy can be simplified as Now, we show that c = 0 is the global minimizer of E for c ≥ 0. First from (29), we can compute the derivative After some tedious algebra, we can show that for any c > 0 and b ∈ (0, 1). Therefore, the energy E is minimized precisely when c = 0 and there is no singular component (R 2 − t 2 ) − b+1 2 in the solution (28). The radius R is given by the same formula (26).
Finally, if b ∈ (−1, 0), no derivative has to be taken and the original governing equation (GE) can be written as with the rescaled second order moment From (10) and (11) with ν = −b, the solution is given by (30) ρ then by exactly the same procedure as for the case b ∈ (0, 1), one can show that the energy E is an increasing function of c for c > 0. Therefore, the stable steady state is given by (25) for c = 0 with the radius (20).
for the coefficient of the second term in (30), it is tempting to conclude that the energy is minimized if and only if this coefficient is zero. However, when E is taken as the only free parameter, both R and M 2 solved from the definition of the first two even moments depend on E too. As a result, it is not obvious that the coefficient (bM 2 (E) − M 0 )R(E) 2 /2 − bE is an increasing function of E. One illuminating counterexample is to exam the case b ∈ (0, 1), where −bE becomes an decreasing function of E but the whole expression (bM 2 − M 0 )R 2 /2 − bE has the opposite monotonicity.
Although the approaches are different for b in different ranges of the interval (−1, 2), the final stable solution takes the same form. This uniformity of the solution is a direct consequence of the fundamental governing equation (GE), since the integral equations (24) and (27) are obtained from (GE) by taking derivatives. This observation is the key to derive general solutions in higher dimensions later.
3.3. Extension with larger even integers of a and in higher dimensions. As we can see from the previous two subsections, the derivation of the exact steady states is in general quite involved, with different solutions to integral equations and the establishment of integral identities in different parameter regimes. The situation is further complicated by the appearance of free parameters that usually have to selected by minimizing the energy.
The above approach using explicit solutions of the integral equations (8) or (12) cannot be easily generalized to higher dimensions either, because the explicit dependence on angular integrals. One exception is in three dimensions, because of the following fact with r = |x| and s = |y|. This reduction is not useful yet, as there are two nonlocal terms associated with one single power. But these two terms can be combined into one integral by extending the radial density ρ(r) evenly on the interval [−R, R], i.e., Therefore, the governing equation becomes from which the steady states can be derived similarly as in one dimension, by taking derivatives on both sides first. Since higher powers appear in the kernels, the derivative of the steady states is more cumbersome than that in one dimension, hence is not pursued here. In summary, the above approach by solving integral equations like (8) or (12) is getting more involved with larger exponents in one dimension, and is unlikely to be generalized to higher dimensions. However, certain features like the universal form of the solutions across different parameter regimes indicate other unified ways, that we exploit next.

Exact solutions in higher dimensions
Now we focus on the construction of the steady states in general dimensions, based on some insights from the case of a = 2 above: the solutions take the same form in different parameter regimes and certain component never appears.
Assuming the radial solution ρ is supported on the ball B R = {x ∈ R d , |x| ≤ R} and a = 2k is an even integer, the left hand side of the equivalent governing equation is an even polynomial of degree 2k, whose coefficients can be written as moments of ρ. It turns out that all the solutions of (31) can be built successively from the fundamental identity which is well-known in potential theory [35,Appendix]. This identity (and the ones with higher exponents below) is expected from dimensional analysis, by counting the powers of R on both sides. Once the identity (32) is established, others with larger exponents can be derived by integration as those in the end of Section 2. Denote ∆ x the Laplace operator, then By the radial symmetry of the integrals, the Laplace operator ∆ x can be inverted easily, giving where the integration constant can be obtained from special integrals like (43) in the Appendix, by evaluating both sides at the origin. Replacing b with b − 2 to keep the same kernel |x − y| b , we arrive at The range of validity for b inherited from (32) is (2 − d, 4 − d), but can be increased to (−d, 4 − d) by continuation. Integral identities with higher powers on (R 2 − |y| 2 ) can be obtained similarly, which are special cases (for k being integers) of the relation in view of the connection to the (negative) fractional Laplacian [32,Eq (9)]. When k is a non-negative integer, the Gauss hypergeometric function 2 F 1 on the right hand side is a finite polynomial and (33) can be written as with the following expressions Here the Pochhammer symbol (r) n = Γ(r + n)/Γ(r) = r(r + 1) · · · (r + n − 1) is used. The range of b also increases with larger exponents, which becomes (−d, 2 + 2k − d) for (34), exactly when the components (R 2 − |y| 2 ) k− b+d 2 and |x − y| b of the integrand are integrable.
Equipped with these integral identities, we are now in a position to construct the steady state for general even integer a = 2k. First, the governing equation (31) can be written as for some coefficients F 0 , F 1 , · · · , F k depending on the rescaled moments M 2j = R −2j B R |x| 2j ρ(x)dx. More precisely, the coefficient F j only depends on the moments M 0 , M 2 , · · · , M 2(k−j) , but not on any higher order moments. From the integral identities above, the steady states take the form (36) for some constants A 0 , A 1 , · · · , A k to be determined. Here appropriate powers of R are used to make the coefficients F j and A j dimensionless and to simplify the calculations later. Motivated by the solutions from the special case derived in the last section, we now make some formal argument for A 0 vanishing identically in (36). First, recall that we expect the radial steady states of our interest here to exist only when as determined in [4]. In fact, for b > b max the preferred asymptotic radial stationary state seems to be the uniform measure or Dirac Delta on a particular (d − 1)-dimensional sphere [34], although it has not been rigorously proved. Under the constraint a > b, it is easy to see that b max is always less than two and is exactly the reason we only focus on a = 2k, an even integer.
is no longer integrable near the boundary |x| = R. Therefore, A 0 must vanish identically for the density to be well-defined. For b belonging to the second interval (−d, 2 − d), the density ρ is shown to be zero in [13] on the boundary |x| = R by formulating the governing equation (31) as an equivalent obstacle problem, which implies the coefficient A 0 of the unbounded yet integrable component (R 2 − |x| 2 ) − b+d 2 must be zero. In this range, the corresponding interaction energy is minimized by making a smooth transition to the zero density outside the support, instead of creating an unbounded jump when A 0 is positive.
With this new information A 0 ≡ 0, we can now construct the steady states in a straightforward way. We show the calculation for a = 2 and a = 4 first to demonstrate the procedure, and then go to general cases when a is an even integer. For a = 2, the governing equation is is determined uniquely by matching the coefficient of |x| 2 on both sides of (37), to obtain (38) ρ The radius R is fixed from the definition of the total mass, that is In one dimension, (38) reduces to the solution (25). When a = 4, the governing equation becomes Upon substituting ρ(x) = (R 2 − |x| 2 ) 1− b+d 2 (A 1 R 2 + A 2 (R 2 − |x| 2 )) followed by the general identity (34), the matching condition of the coefficients of |x| 2 and |x| 4 on both sides of (39) becomes Exacting C kj from (35), the solution of this linear system is given by independent of M 4 and R. To decide the radius of the support, we use the definition of the first two moments, From the explicit expressions (40) of A 1 and A 2 , these two definitions of the moments are homogeneous equations in M 0 and M 2 , and is equivalent to the eigenvalue problem The matrix on the right hand side has all positive entries. By Perron-Frobenius theorem [5], there is only one positive eigenvector M = [M 0 , M 2 ] T corresponding to the unique positive eigenvalue. In this case, the special eigenvalue can be obtained explicitly, and the resulting radius is given by The corresponding eigenvector can also be parameterized by the total mass M 0 , that is M 2 = dM 0 / (2 − b)(6 − b) and the coefficient A 1 can be simplified as Therefore, the solution ρ is completely determined. Now we can proceed with the general case a = 2k, and the steps are outlined below: Step a) The solution takes the following form Step by matching the coefficients of |x| 2 , · · · , |x| 2k on both sides of the governing equation.
Step c) The definitions of the (rescaled) moments can be formulated as an eigenvalue problem R b−2k M = D M . The radius R is the 1/(b − 2k)-th power of the eigenvalue of D and the eigenvector, or equivalently the moments, can be parameterized by the total mass M 0 .
Remark 2. Explicit forms of the coefficients can be worked out using symbolic packages for the general case a = 2k. But the calculation is more and more involved, and numerical softwares using floating point are preferred.  However, the solutions constructed in this way are not always physically relevant, because negative densities could appear. From the expressions (40) of A 1 and A 2 for a = 4, it is easy to see that the solution ρ(x) = (R 2 − |x| 2 ) 1− b+d 2 (A 1 R 2 + A 2 (R 2 − |x| 2 )) becomes negative at the origin when Sinceb is less than b max , ρ(0) is indeed negative for b betweenb and b max (see Figure 4). When a is an even integer larger than 4, numerical experiments show that the solutions constructed in the above approach is still negative near the origin, when b is close to its upper bound b max . In fact, extensive particle simulations indicate that a void (a region with zero density) starts to appear near the origin [3]. As b continues to increase, the density becomes concentrated more and more towards its outer boundary, and eventually collapses on a sphere as b passes b max . Because the solution is no longer supported on a ball, the above approach using integral identities like (33) does not work, and precise form of the solutions remains unknown.

Further discussions, conclusions and conjectures
Besides providing more exact solutions for theoretical analysis and numerical testing, these newly discovered exact solutions confirm many widely observed but yet unproved phenomena about boundary regularity even with the simplest power-law attractive-repulsive kernels. From the general form (41), it is easy to see that all constructed stationary solutions in this paper behave like (R − |x|) 1−(b+d)/2 near the boundary, depending only on the singular repulsive part of the kernel. More precisely, the Newtonian repulsive kernel dictated by b = 2 − d is critical: the solutions is zero on its boundary for b < 2 − d and become infinity (yet integrable) for 2 > b > 2 − d, while at the critical value b = 2 − d, the solution is finite on the boundary. We do expect this behavior to be generic for many potentials with power-law like repulsive component. As a result, if this boundary behavior conjecture is true, any converging solution of the evolution equation (4) can not be uniformed bounded in time for 2 − d < b < 2. Therefore, uniform (in time) L p bounds are expected only for 1 ≤ p < 2/(b + d − 2).
The same boundary behaviors for other general potentials seem to be supported by extensive analytical and numerical studies. For Newtonian repulsion and general attraction, bounded density on the boundary is established as an eigenvalue problem in [30,29]. Finite non-zero density on the boundary is proved for Newtonian-like repulsion like the Quasi-Morse potential [6]. The inverse square root singularity of the solutions for the Morse potential in two dimension (with b = 1) is also confirmed numerically in [39]. The boundary behaviors are rigorously proved for Newtonian or more singular than Newtonian repulsions and quite general attractions in [13] using the equivalence to classical obstacle problems in elliptic equations of the necessary conditions (6). This work may motivate further investigation towards the regularity of the solutions near boundary.
However, important questions like boundary behaviors are only inferred from special solutions and observed from limited numerical simulations. The extend to which these statements remain true are still widely open, and rigorous proofs may require sophisticated techniques from harmonic analysis and potential theory. Because the solutions stop to be valid due to the appearance of negative density as b is close to b max , the exact way in which the solution collapses to a co-dimension one sphere is not clear either. The construction of these nearly singular solutions is another interesting problem.
The Beta function can be represented using the Gamma function, that is B(α, β) = Γ(α)Γ(β)/Γ(α + β). Certain special integrals on the ball B R also appear repeatedly in this paper, with in one dimension and more generally Appendix B. Derivation of the solutions to the singular integral equations When the right hand sides of two singular integral equations (8) and (12) are polynomials, the solutions using the complicated formula (9) and (13) can be simplified significantly. The key step is to evaluate the principal integral of the form with n = 0 in (44). Singular integrals like these can be evaluated using Plemelj Formulae for sectionally analytic functions, which are analytic on the complex plane except a few contours [1,26]. is analytic on C \ L. For any real number x ∈ (−R, R), Ψ has two limiting values when approaching x from above and below of L, denoted as Ψ + (x) and Ψ − (x) (see Figure 5). The Plemelj Formulae relates ψ to the jump and average of Ψ + (x) and Ψ − (x) in the following way, Now we can calculate principal integrals like (44) using (46b), by finding appropriate complex function Ψ whose jump across L is exactly (R + y) n [(R − y)/(R + y)] (1−ν)/2 . from the relation (46a), the sectionally analytic function Ψ given by (45) has the same jump ψ across L as Φ. Therefore, Ψ − Φ is continuous across L and is analytic on the whole complex plane. On the other hand, since both Ψ and Φ are zero at infinity, so is their difference Ψ − Φ at infinity. Hence by Liouville's Theorem, Ψ is identical to Φ, and the second Plemelj formula (46b) becomes sin νπ 2 As a result, the principal value integral (44) with n = 0 and ν ∈ (0, 1) is obtained and the solution of (8) with f (x) = 1 can be shown to be (10). The same procedure can be generalized to larger powers n ≥ 1 and ν ∈ (0, 1), by choosing where P n (z) is a polynomial of degree n such that Φ(z) becomes zero at infinity. For example, we take to get (for n = 2) Once the principal value integral is obtained for any non-negative integer n and ν ∈ (0, 1) in (44), the solution (9) can be evaluated for any polynomial left hand side of (8).
The solutions (13) for any polynomial right hand side of (12) require the evaluation of a slightly different principal value integral P.V.