Metric deformation and boundary value problems in 3D

A novel perturbative method, proposed by Panda {\it et al.} [1] to solve the Helmholtz equation in two dimensions, is extended to three dimensions for general boundary surfaces. Although a few numerical works are available in the literature for specific domains in three dimensions such a general analytical prescription is presented for the first time. An appropriate transformation is used to get rid of the asymmetries in the domain boundary by mapping the boundary into an equivalent sphere with a deformed interior metric. The deformed metric produces new source terms in the original homogeneous equation. A deformation parameter measuring the deviation of the boundary from a spherical one is introduced as a perturbative parameter. With the help of standard Rayleigh-Schr{\"o}dinger perturbative technique the transformed equation is solved and the general solution is written down in a closed form at each order of perturbation. The solutions are boundary condition free and which make them widely applicable for various situations. Once the boundary conditions are applied to these general solutions the eigenvalues and the wavefunctions are obtained order by order. The efficacy of the method has been tested by comparing the analytic values against the numerical ones for three dimensional enclosures of various shapes. The method seems to work quite well for these shapes for both, Dirichlet as well as Neumann boundary conditions. The usage of spherical harmonics to express the asymmetries in the boundary surfaces helps us to consider a wide class of domains in three dimensions and also their fast convergence guarantees the convergence of the perturbative series for the energy. Direct applications of this method can be found in the field of quantum dots, nuclear physics, acoustical and electromagnetic cavities.


I. INTRODUCTION
The three dimensional Helmholtz equation is encountered frequently by physicist and engineers in different areas − like the eigenanalysis in acoustic and electromagnetic cavities, transmission of acoustic waves through ducts and in quantum mechanics. The main concern to pursue these problems is to calculate the eigenspectrum of the linear homogeneous  [2,3]. Constructing on the boundary a co-ordinate system suitable to it often helps in finding the solution. For the Helmholtz equation in three dimensions, there are 11 such 'special' co-ordinate systems [4]. However, even in these sacred classes, the multi-parameter dependence of the solution makes them difficult to use. So, for a general domain, finding solutions to the Helmholtz equation becomes a formidable job. However, many physical problems demand us to solve the equation for an arbitrary domain which is significantly deviated from the above mentioned idealised scenario. In particular, the geometry of the nanoscale second-phase particles are described analytically as superellipsoids [5] and experimentally verified from just one micrograph [6]. To explore the shapes of these particles more accurately, we need a good estimation of the eigenspectra of these geometries which urge us to solve the Helmholtz equation for these arbitrary boundaries. Recently in the field of medical science for automatic prostate segmentation using deformable superellipses have been efficiently applied [7]. The analytic approach towards this problem came recently by Panda and Hazra [8] where they have generalised Rayleigh's theorem in three dimensions and used the standard time independent perturbation method of quantum mechanics to calculate the eigenfunction and eigenvalue corrections in closed forms up to the second order of perturbation for general shapes in three dimension in terms of the spherical harmonics expansion coefficients. Most of the earlier attempts made towards this problem were via numerical means in the area of chaos for general three dimensional billiards [9][10][11][12][13][14] and also in some recent experiments to study the wave chaos in microwave [15] and resonant optical cavities [16]. Among the numerical schemes, the finite difference method (FDM), the finite element method (FEM) [17] and the boundary element method (BEM) [18] are popular ones but they consume huge amount of time to generate the mesh for a complicated geometry.
Recently as alternatives to the above, some meshless methods have been developed, like the boundary collocation method [19], the radial basis function [20], the method of particular solutions [21] and the method of fundamental solutions [22]. These methods also have the drawback of getting spurious eigensolutions.
In this paper, we have prescribed an alternative method to solve the boundary value problem for the Helmholtz equation for a general simply connected convex domain in three dimensions. This analytic formulation is an extension of our earlier work for two dimensions [1] where a suitable diffeomorphism in terms of a Fourier series was chosen to map the general problem into an equivalent one where the boundary was circular but the equation gets complicated due to the deformation of the metric in the interior and was solved by perturbation technique of quantum mechanics. In this method, an arbitrary domain in three dimensions is mapped to a regular closed region (in which the Helmholtz equation is exactly solvable) by a suitable co-ordinate transformation resulting the deformation of the interior metric. As a result, the homogeneous Helmholtz equation gets modified and can be written as the original one with additional source terms present on the right hand side of it. These extra terms can now be treated as a perturbation to the original equation and is solved using the standard perturbation technique. The corrections to the eigenfunction are obtained in closed form at each order of perturbation irrespective of the boundary conditions. The eigenvalue correction at each order is calculated by applying the proper boundary conditions on the respective eigenfunction corrections. To be able to express the eigenfunction corrections independent of the boundary conditions is a unique advantage of this method over the existing ones. Also, in this method the boundary conditions maintain a simple form at each order of perturbation and is easy to apply on some regular closed surface. In this analysis, we have used spherical harmonics to represent the asymmetries of the boundary from its equivalent sphere which allows us to implement the method for a large variety of geometries in three dimensions. To test the efficacy of the perturbative method we have tabulated the analytical values with the corresponding numerical ones for spheroids, supereggs, stadium of revolution, rounded cylinder and pear shaped enclosures.
The method seems to work extremely well for these cases for both Dirichlet and Neumann conditions. It has a potential of applicability for various general shapes in three dimensions, where the deviations are large from a spherical shape.
The paper is organised as follows : the section II describes the general formalism in the abstract sense. In section III, we have applied the method to various enclosures. The last section IV is reserved for results, conclusions and comments.

II. FORMULATION
The homogeneous Helmholtz equation reads, where g ij is the background metric component on three dimensional space (V) and ∇ represents a covariant derivative. △ is the Laplacian operator in three dimensions. Our interest lies in finding the eigenspectrum and the solutions in the interior of a closed simply connected convex region satisfying either the Dirichlet boundary condition (DBC), ψ = 0, or the Neumann boundary condition (NBC), ∂ψ ∂n = 0, on the boundary ∂V, where ∂ψ ∂n is the derivative along the outward normal direction to ∂V. Depending on the boundary condition the parameter k 2 can be matched to a physical quantity (for the DBC it is proportional to the energy of a quantum particle confined within ∂V, where as for the NBC it will determine the square of the resonating frequencies for acoustical cavities).
For convenience, we choose to work in spherical polar coordinate system (R, Θ, Φ). Now, We consider a general arbitrary surface of the form R = R(Θ, Φ). The periodicity condition implies that R(Θ, Φ + 2π) = R(Θ, Φ). We assume that any general shape (which is not very elongated in one direction) can be expressed as a deformation around an effective spherical boundary. We choose spherical boundary because the Helmholtz equation in 3D is exactly solvable for this boundary (the analysis in principle can work for perturbation around any simple surface for which the Helmholtz equation is exactly solvable). We make a coordinate transformation from the old one (R, Θ, Φ) to a new one (r, θ, φ) of the form where ǫ is a deformation parameter. Now, for an intelligent choice of the well behaved function f (θ, φ), our arbitrary surface will be transformed into a sphere of average radius, in the (r, θ, φ) system. As a result of the deformation, there will be changes in the underlying metric component g ij (R, Θ, Φ). Henceforth, we use the notation The dependence on the arguments (θ, φ) is not shown explicitly for brevity. The flat background metric in the old coordinate system (R, Θ, Φ) is given by g = diag (1, R 2 , R 2 sin 2 Θ).
Under the coordinate transformations (2) it takes the form The non-vanishing components of connection Γ for the flat background metric (g) given by, As a result of the coordinate transformation (2) no spurious curvature is induced in the manifold (i.e. Riemann tensor, R i jkl = 0 ∀ i, j, k, l). Under the map (R, Θ, Φ) → (r, θ, φ), Eq. (1) takes the following form where ψ (i,j,k) ≡ ∂ i+j+k ψ ∂r i ∂θ j ∂φ k . After some simplification, Eq. (4) reduces to ∞ n=0 ǫ n H n ψ + Eψ = 0, where the operator H i 's are given by . . . and We will now implement the standard Rayleigh-Schrödinger perturbation theory to solve for ψ and E. Perturbed eigenfunction ψ and the corresponding eigenvalue E are expanded in a power series of the perturbation parameter, ǫ, as where superscripts denote the orders of perturbation.
Substituting Eq. (10) in Eq. (5), and setting the coefficients of different orders of ǫ to zero, yields We can infer from the above equations, (11), that new source terms have been generated at each order in ǫ to the unperturbed homogeneous Helmholtz equation as a by product of the mapping given by (2). So, we started with an arbitrary boundary with the absence of sources inside the domain and with a mapping effectively generated a spherical boundary with non-vanishing sources inside it. In order to maintain the simple boundary condition, we have incorporated new source terms in the equations. Now, the eigenvalue corrections for different orders can be calculated by The corresponding boundary conditions for the DBC and the NBC are respectively and for all i ∈ N, where the radius of the sphere R 0 is defined by (3).
The general solution of the Eq. (11a) is given by n,l is calculated using the n th zero [23] of j l , denoted by β n,l , and for the NBC, it is dictated by the n th zero of j ′ l (i.e. derivative of j l with respect to its argument), denoted by α n,l . The boundary conditions Eq. (13) and Eq. (14) where all the levels with non-zero l are (2l + 1)-fold degenerate.
In this prescription, the energy corrections can be obtained in two ways. Primarily, E (i) is extracted out by imposing the respective boundary conditions on ψ (i) given by Eqs. (13) and Eqs. (14), which as a bonus give the coefficients of Bessel functions (in ψ (i) ). Alternatively, it can be verified using Eqs. (12) from the information of ψ (m) (∀ m < i). In principle, this formulation can be applied to obtain correction at all order of perturbation. In the following we calculate the eigenvalue as well as eigenfunction corrections for both the cases of boundary conditions. Till now the formalism was proceeding in an abstract sense as it did not require a particular form of f . From now on without a loss of generality we choose the following specific form for f in terms of spherical harmonics where C b a are the expansion coefficients. The constant part C 0 0 can always be absorbed by redefining the R in (2).
A. Non-degenerate states (l = 0) The first order correction to the eigenfunction is obtained by solving the Eq. (11b). Thus, where the unknown coefficients E (1) ) and A q p will be calculated by imposing the respective boundary conditions. The terms containing N n,0 ρj 1 (ρ) make the particular integral of the Eq. (11b). Now, applying the boundary conditions, given by Eq.
(13) and Eq. (14), for i = 1, we extract out the first order eigenvalue corrections as well as the unknown expansion coefficients as E (1) n,0 = 0 (for both the cases); It is evident that the first order eigenvalue correction is zero for both the cases of boundary conditions which is in confirmation with Eq. (12b). The orthogonality relation between n,0,0 and ψ (1) n,0,0 dictates that the remaining constant A 0 of Eq. (19) is zero for both the boundary conditions. The eigenfunction correction for the second order is given by solving Eq. (11c) as The first non-zero eigenvalue correction E (2) n,0 is obtained with the knowledge of ψ (0) n,0,0 and ψ (1) n,0,0 and imposing the boundary conditions, Eq. (13) and Eq. (14), for i = 2, one extracts out and The coefficients B q p are determined as where n 1 n 2 m 1 m 2 |nm gives the Clebsch-Gordan coefficient for the decomposition of |n 1 , n 2 , n, m in terms of |n 1 , m 1 |n 2 , m 2 and |a−p| |q−b| ≡ Maximum (|a − p|, |q − b|). The unknown coefficient B 0 can be fixed by the normalisation condition of the eigenfunction corrected up to second order. These results are consistent with the earlier work done by the other method [8]. The similarities of these results with the two dimensional metric deformation results [1] can also be noticed. With this simplification the expansion for f given in (18) will reduce to where we denote C a ≡ C 0 a . The above expression implies that the boundary has azimuthal symmetry. Due to this simplification the boundary conditions for NBC in (14) will reduce further and all the terms containing φ-derivative of f will go to zero for the axisymmetric geometries. For l = 0 case, the first order eigenfunction correction, given by is a solution of Eq. (11b) with ψ n,l,m given by (15b). The first order eigenvalue correction is evaluated by setting the respective boundary conditions given in Eq. (13) and Eq. (14) for i = 1 and it is also confirmed with the Eq. (12b). Thus, we have where the corresponding unperturbed eigenvalues E n,l are given by Eq. (16) and Eq. (17) respectively. There is a non-zero correction to the eigenvalue even at the first order unlike the non-degenerate case. Further, the coefficients A q p s (for p = l) are extracted as, where rest of the coefficients A q p for q = m are zero. The remaining coefficient A m l is calculated from the normalisation condition and found to be By solving Eq. (11c) we get the second order correction to the eigenfunction as n,l . Imposing the boundary conditions (13) and (14), for i = 2, we obtain the second order correction to the eigenvalue for the DBC and the NBC respectively as follows n,l = α 2 n,l − 3l(l + 1) α 2 n,l − l(l + 1) 2l k=|a−s| (2a + 1)(2s + 1) 2π C a C s as00|k0 2 kl00|l0 kl0m|lm × sp00|l0 sp0m|lm 1 + k(k + 1) + l(l + 1) − p(p + 1) 2{α 2 n,l − l(l + 1)} × × 1 + 2α 2 n,l + s(s + 1) − p(p + 1) − l(l + 1) 4 j p (α n,l ) α n,l j ′ p (α n,l ) .
The expansion coefficients B q p s are algebraically complicated to calculate and are not needed for our present purpose. These results are matching with the results obtained in our earlier paper [8] by a different method.

III. EXAMPLES
In the previous section we have described the formalism in an abstract sense and now we apply it to estimate the energyspectra for various axisymmetric boundary surfaces like spheroidal, superegg, stadium of revolution, rounded cylinder and pear shaped enclosures.
These geometries are naturally encountered in nuclear physics [24] and in the experiments on nanoscale structures [6]. The analytic results have been compared against the numerical ones obtained by using finite element method (with the help of Matlab and Mathematica) and are tabulated below for the above mentioned domains satisfying different boundary conditions.

A. Supereggs
We consider superegg [25] shaped enclosures which are surface of revolution of supercircles [26] about either of its in-plane axes. The representation of it in the spherical polar coordinate is given by r(θ, φ) = 1 (| cos θ| n + | sin θ| n ) 1/n ,

B. Rounded cylinder and Stadium of revolution
Next, we consider a rounded cylindrical enclosure in three dimensions. If we consider a rounded rectangular region with circular arc of radius R and the separation in between them of d in a plane, then the surface of revolution around the major axis will be a rounded cylinder of height and base diameter equal to 2(d+R). The equation of the rounded rectangle in polar coordinate is and the form of the rounded cylinder in spherical polar coordinate is for θ ∈ [0, π] and φ ∈ [0, 2π].
The shape of the rounded cylinder is displayed in fig (2a). Further, we consider the stadium of revolution in three dimensions. It is a surface of revolution attained by rotating a half stadium shape in two dimensions about its major axis. The equation of the half stadium in polar coordinate is where R is the radius of the circular arc and d is the separation between them. The angular co-ordinate θ is in the usual sense. The form of stadium of revolution in three dimensions in spherical polar coordinate is simply given by for θ ∈ [0, π] and φ ∈ [0, 2π].
The shape of the stadium of revolution in three dimensions is depicted below in fig (2b).

C. Spheroids
A spheroid in the spherical polar coordinates is given by where r a and r c (> 0) are called equatorial and polar radii respectively. For r c < r a the spheroid is known as oblate while for r c > r a it is called prolate. We have considered

D. Pear shaped enclosures
Finally, we consider the following pear shaped domains motivated by the recent work [24] in the experimental nuclear physics for the 220 Rn and 224 Ra nuclei. The equation of the pear shape in spherical polar coordinate is given by where R 0 is the radius of equal volume spherical shape, which we set to one. The values for the expansion coefficients C a s are identified with the parameter β a s and chosen as the same given in [24,27]. The shapes with these values are shown in fig (4).

IV. RESULTS AND DISCUSSIONS
These particular forms of r(θ, φ), given in (25) imply that this perturbative formalism works considerably well even for higher excited states and also for domains which are highly deformed. For this extreme cases addition of higher order corrections will surely improve the matching but they are algebraically complicated.
Due to the deformation from the sphere the degenerate states, of multiplicities (2l + 1) for a given l value, will be splitted. It is observed that for degenerate states the energy  to the accidental closeness of the zeros of derivative of spherical Bessel functions j 0 and j 3 , which are α 1,0 (= 4.49341) and α 1,3 (= 4.51409) respectively. In the calculation of the second order energy correction for l = 0 due to the non-zero finite value of C 3 , the term of the 2π , picks up a very high negative value (as the denominator is very close to zero) and causes the problem in the levels marked by stars (⋆) for both the shapes.
Similarly, for l = 3 and m = 0 the second order energy correction term, containing , produces a large positive correction as the denominator is again very close to zero and which creates the error for the levels marked by daggers ( †) for both the cases. Other than the above mentioned exception the matching between the perturbative and numerical values are outstanding. The agreement between the numerical and perturbative values are better for supereggs than that for spheroids for both the cases of boundary conditions and is analogous to the behaviour of their two dimensional counterparts i.e. between supercircle and ellipse as reported in earlier works [28,29]. We would like to point out that the eigenvalue correction for the state l = 2, m = 0 behaves in a peculiar fashion and disrupts the matching between different levels and produces the instances of maximum deviations for most of the cases in the above mentioned examples for both the boundary conditions. This abnormality for l = 2 is an issue which calls for further investigations. Another unique feature of this method is that it can capture the degeneracy patterns for apparently similar looking shapes having different degeneracy patterns and also distinct shapes having similar type of degeneracies.
As an example, the prolate, the stadium of revolution and the superegg (n = 2.5) are looking similar apparently but the first two shapes have similar degeneracy pattern only for the DBC.
The distinction is made in their NBC degeneracy patterns. In contrast, rounded cylinder and superegg (n = 2.5) are having different structures but they produce similar degeneracy pattern for the DBC. So, our perturbative method traces the degeneracy patterns (provided by the numerical results) in a correct fashion for low-lying levels for different shapes and both the cases of boundary conditions. On the other hand, there are a few instances where three very close energylevels show degeneracy mismatch (viz. superegg (n = 1.7) for DBC levels 5, 6 and 7; oblate for NBC levels 13, 14 and 15; rounded cylinder for NBC levels 11, 12 and 13). We expect that inclusion of higher order correction term may provide a positive correction to the non-degenerate state whereas the two degenerate states may get a negative correction and which might bring them in order.
In conclusion, we point out that expansion of the boundary asymmetry in terms of spherical harmonics makes the method completely general to encompass a large variety of domains in three dimensions. The closed form corrections for wavefunctions and energies at each order of perturbation in terms of expansion coefficients help us to apply the prescription for a general boundary conditions also. Since the solutions are boundary condition free this method can tackle variety of problems. Also the boundary conditions maintain a simple form at each order of perturbation at the cost of complicating the equations. This method can treat both degenerate and non-degenerate states in the same way for axisymmetric boundaries. This method has an edge over the other method discussed by Panda and Hazra recently [8]. In the present method applying boundary condition is easy in contrast to the other method where it is extremely cumbersome. In this case all order equations as well as energy corrections are also written.