A Perturbative Approach to the Tunneling Phenomena

The double-well potential is a good example, where we can compute the splitting in the bound state energy of the system due to the tunneling effect with various methods, namely WKB or instanton calculations. All these methods are non-perturbative and there is a common belief that it is difficult to find the splitting in the energy due to the barrier penetration from a perturbative analysis. However, we will illustrate by explicit examples containing singular potentials (e.g., Dirac delta potentials supported by points and curves and their relativistic extensions)that it is possible to find the splitting in the bound state energies by developing some kind of perturbation method.


Introduction
Most real quantum mechanical systems can not be solved exactly and we usually apply some approximation methods, the most common one being perturbation theory, to get information about the energy levels and scattering amplitudes. However, not all quantum systems can be analyzed by perturbative methods. There are various problems for which one can not deduce any information by simply using perturbation theory since these problems are inherently non-perturbative phenomena like the formation of bound states or penetration through a potential barrier. For such non-perturbative phenomenon, other tools, such as WKB [1,2] or instanton calculations [3], are particularly useful. The particle moving in a one-dimensional anharmonic potential V (x) = λ 2 8 (x 2 − a 2 ) 2 is a classic example, where we can study the barrier penetration through the WKB analysis. When Figure 1: Anharmonic Potential the energy scale determined by the length scale a is extremely small compared with the binding energy of the system, i.e., 2 /2ma 2 << E B , or λa 2 >> 1, the potential separates into two symmetrical wells with a very high barrier (see Fig. 1). In this extreme regime, as a first approximation, each well has separately quantized energy levels and these energy levels are degenerate due to the symmetry. However, once the large but finite value of the coupling constant λ is taken into account, the particle initially confined to one well can tunnel to the other well so the degeneracy in the energy levels disappear. The splitting in the resulting energy levels (between the true ground state and the first excited level due to the tunneling) is given by [1,2] where S 0 = 2mωa 2 3 and ω 2 = λ 2 a 2 m . The above exponentially decaying factor with respect to the separation between the wells illustrates the tunneling effect. The true ground state corresponds to a symmetric combination and the excited level corresponds to the anti-symmetric combination of the WKB corrected wave functions.
Among the exactly solvable potentials in quantum mechanics, Dirac delta well potentials are the most wellknown text book example [4]. Moreover, it has been studied extensively in mathematical physics literature from different point of views, in particular in the context of self-adjoint extension of symmetric operators [5]. Although it is easier to define it rigorously in one dimension through the quadratic forms, one possible way to define it in higher dimensions is to consider the free symmetric Hamiltonian defined on a dense domain excluding the point, where the support of the Dirac delta function is located, and then apply the self-adjoint extension techniques developed by J. Von Neumann (see the monograph [6] for the details and also for the historical development with extensive literature in the subject). Then, the formal (or heuristic) definition of one-dimensional Dirac delta potentials in the physics literature is understood as the one particular choice among the four parameter family of the self-adjoint operators, where the matching conditions of the wave function are just obtained from the boundary conditions (which define the domain of our self-adjoint operator) constructed through the extension theory. Another way to introduce these point interactions uses the resolvent method, developed by M. Krein, and it is based on the observation that for such type of potentials the resolvent can be found explicitly and expressed via the so-called Krein's formula [7]. Within this approach, the Hamiltonian for point interaction (in two and three dimensions) is first approximated (regularized) by a properly chosen sequence of self-adjoint operators H and then the coupling constant (or strength) of the potential is assumed to be a function of the parameter in such a way that one obtains a non-trivial limit. This convergence is actually in the strong resolvent sense, so the limit operator is self-adjoint [8]. Since the Dirac delta potentials in two and three dimensions require renormalization, it is usually considered as a toy model for the renormalization originally developed in quantum field theories and it helps us to better understand various ideas in field theory such as renormalization group and asymptotic freedom [9,10,11,12]. Furthermore, point like Dirac delta interactions have also been extended to various more general cases. For our approach, to illustrate the main ideas, we are mainly concerned with the delta potentials supported by points on flat and hyperbolic manifolds [13,14,15], and delta potentials supported by curves in flat spaces, and its various relativistic extensions in flat spaces [16,17,18,19].
In this paper, we explicitly demonstrate for a class of singular potential problems that the splitting in the energy levels due to the tunneling can be realized by simply developing some kind of perturbation theory. We have two basic assumptions here: 1) Binding energies of individual Dirac delta potentials are all different. Otherwise we need to employ degenerate perturbation theory. Actually, we briefly discuss a particular degenerate case, namely the two center case to compare with the double well potential. 2) The support of singular interactions are sufficiently separated from one another, as a result the bound state wave functions decay rapidly over the distances between them.
All the findings about the splitting in the bound state energies for singular potentials on hyperbolic manifolds treated here could be applied to the two dimensional systems such as graphite sheets. We can model impurities in these systems as attractive centers in some approximation and these sheets can be put in various shapes. This is especially true for surfaces with variable sectional curvature which is not completely negative. The negatively curved surfaces, of course, cannot be realized as embedded surfaces in three dimensions due to Hilbert's wellknown theorem. Nevertheless, we may envisage these models as an effective description of unusual quasi-particle states of some two dimensional systems. Due to the interactions, the system may develop a gap in the spectrum and the effective description may well be best understood through a negative sectional curvature space. Models related to point interactions on Lobachevsky plane have been studied from variety of different perspectives in [20,21]. The point interactions can be extended on more general class of manifolds as well [22]. In particular, they have been studied on some particular surfaces in R 3 , namely on the infinite planar strip as a natural model for quantum wires containing impurities [23] and on the torus [24]. A more heuristic approach for point interactions on Riemannian manifolds has been constructed through the heat kernel in [13,14]. The physical motivation behind studying the Dirac delta potentials supported by curves is based on the need for modelling semiconductor wires [25]. They could be considered as toy models for electrons confined to narrow tube-like regions.
The paper is organized as follows: In Section 2, we formally summarize the resolvent formulae, called Krein's formulae, for Hamiltonians perturbed by singular potentials including Dirac delta potentials supported by points and curves. The principal matrix for each case is given explicitly. The relativistic and the field theoretical extension of it has been also reviewed in the subsections of this section. In Section 3, we briefly discuss the analytic structure of the principal matrix and the bound state spectrum for such singular interactions. In Section 4, we discuss how the off-diagonal terms of the principal matrix change in the tunneling regime. Section 5 and Section 6 contain the formulation of the perturbative analysis and explicit calculations of the splitting in the bound state energy when these singular potentials are placed far away from each other, which is the main result of the paper. We finally discuss the degenerate case and wave functions, and compare these results results with the exact result in Section 7 to get a feeling for the accuracy of our approximation.

Krein's Formulae for Free Hamiltonians Perturbed by Singular Interactions
Before we are going to discuss the perturbative analysis of singular interactions for large separations of the supports of these potentials, we first present the basic results about our formulation of the singular Hamiltonians.
In this paper, we are mainly concerned with the Dirac delta potentials supported by finitely many points and finitely many curves in flat spaces, and their extension to the hyperbolic manifolds. Moreover, we also consider some relativistic extensions of these singular interactions.
Since we study the spectral properties of different kinds of Dirac delta potentials, we first introduce the notation for Dirac delta functions of interest. The Dirac delta distribution δ a supported by a point a in R n is defined as a continuous linear functional whose action on the test functions ψ is given by Similarly, Dirac delta distribution δ γ supported by a curve Γ in R n is defined as a continuous linear functional whose action on the test functions ψ is given by [26] δ γ , ψ = Γ ds ψ(γ(s)) .
The left hand sides in the definitions (2) and (3) can be expressed in the Dirac's bra-ket notation, most common in physics literature, as a|ψ and γ|ψ , respectively. As we have already emphasized in the introduction, there are several ways to define rigorously the Hamiltonian for Dirac delta potentials. Here, we start with a finite rank perturbations of self-adjoint free Hamiltonian H 0 (e.g., H 0 = P 2 /2m in the non-relativistic case and H 0 = √ P 2 + m 2 in the semi-relativistic case): where ϕ i ∈ H and ., . denotes the sesqui-linear inner product in the Hilbert space H. Then, it is well-known that the resolvent of H can be explicitly found in terms of the resolvent of the free part by simply solving the inhomogenous equation [7] (H − z)ψ = ρ , for a given ρ ∈ H and ψ ∈ D(H 0 ) = D(H). Here D stands for the domain of the operator and we assume that (z) > 0. It is well-known that H is self-adjoint on D(H 0 ) due to the Kato-Rellich theorem [5]. The resolvent of H could be found in two steps: First, we apply the resolvent of the free part to the equation (5) and project this on the vector ϕ j , we can then find the solution ϕ i , ψ so that the resolvent R z (H) = (H − z) −1 of the Hamiltonian H at z is: where Actually, the resolvent formula (7) is valid even in the case where the vectors ϕ i 's do not belong to the Hilbert space. Such perturbations represent the singular type of interactions, e.g. Dirac delta potentials supported by points or curves [6,16]. In Dirac's bra-ket notation, one can also express the above resolvent formula as: The explicit expression of the resolvent (7) or (9) is known as Krein's resolvent formula. Alternatively, these singular interactions can be defined directly through von Neumann's self-adjoint extension theory (or quadratic forms in some cases). Since our aim is the spectral behaviour and especially the bound state problem of such singular interactions, Krein's explicit formula is much more useful. Throughout the paper, following the terminology introduced by S. G. Rajeev in [27] we call the matrix Φ as the principal matrix (this is equivalent to the matrix Γ used in [6]). Actually, one can also develop the above resolvent formula (9) to relativistic and field theoretical extensions of the singular models, as we will discuss in the next subsections. Let us now summarize explicitly the resolvent formulae and principal matrices in all classes of singular interactions that we are going to discuss in this paper:

Point-like Dirac delta interactions in R
The Hamiltonian for the non-relativistic particle moving in fixed N point like Dirac delta potentials in one dimension can be expressed in terms of the formal projection operators given by the Dirac kets |a i where H 0 is the non-relativistic free Hamiltonian, and λ j 's are positive constants, called coupling constants or strengths of the potential. Throughout this paper, we will use the units such that = 2m = 1 for non-relativistic cases and = c = 1 only for the relativistic case. Since we have fairly complicated expressions, this simplifies our writing, hoping that this does not lead to any further complications. It is well-known in the literature that there are different ways to make sense of this formal Hamiltonian in a mathematically rigorous way (strictly speaking, the above expression (10) has no meaning as an operator in L 2 (R)). Let us define R z (H) := R(z) and R z (H 0 ) := R 0 (z) for simplicity. Even though it is hard to make sense of the Hamiltonian, one can find the resolvent of this formal operator algebraically and the result is consistent with the one given by a more rigorous formulation. Choosing ϕ i as the Dirac kets |a i formally in the previous section, the resolvent is explicitly given by Here R 0 (a i , a j ; z) = a i |(H 0 − z) −1 |a j is the free resolvent kernel. It is useful to express the principal matrix in terms of the heat kernel K t (a i , a j ) -the fundamental solution to the Cauchy problem associated with the heat equation -using Then, we obtain These expressions should be considered as analytical continuations of the formulae beyond their regions of convergence in the variable z. From the resolvent (11), one can also write down the resolvent kernel Using the explicit expression of the integral kernel of the free resolvent we have Here √ z is defined as the unambiguous square root of z with √ z is positive. Since we study the bound state spectrum, it is sometimes convenient to express the above matrix Φ(z) in terms of a real positive variable

Point-like Dirac delta Interactions in R 2 and R 3
We assume that the centers of the Dirac delta potentials do not coincide, that is, a i = a j whenever i = j. If we follow the same steps outlined above, we find exactly the same formal expression for the resolvent for point interactions in two and three dimensions except for the fact that the explicit expression of the integral kernel of the free resolvent in R 2 and R 3 [6] are given by respectively. Here H (1) 0 is the Hankel function of the first kind of order zero and √ z > 0. Unfortunately, the diagonal part of the free resolvent kernels are divergent so the diagonal part of the principal matrices are infinite. This is clear for the three dimensional case from the asymptotic behavior of the Hankel function [28] H (1) as x → 0. This difficulty can be resolved by the so-called regularization and renormalization method. Instead of starting with the higher dimensional version of the formal Hamiltonian (10), we first consider the regularized Hamiltonian through the heat kernel where r|a i = K /2 (r, a i ). The heat kernel associated with the heat equation ∇ 2 ψ − ∂ψ ∂t = 0 in R n is given by It is important to note that as → 0 + in the distributional sense. Then, we can easily find the resolvent kernel associated with the regularized Hamiltonian (22) where If we choose 1 where E i B < 0 (the spectrum of the free Hamiltonian only includes the continuous spectrum: [0, ∞]) is the bound state energy of the particle to the i th center in the absence of all the other centers and take the formal limit → 0 + we find where From the explicit form of the heat kernel formula (23), we obtain in two dimensions and in three dimensions.
Since we deal with the bound states in this paper, it is convenient to express the principal matrices in terms of the real positive variable ν = −i √ z: in two dimensions and in three dimensions. Here we have used K 0 (z) = iπ 2 H 1 0 (iz) with −π < arg(z) < π/2 and K 0 (z) is the modified Bessel function of the third kind [28].

Point-like Dirac delta Interactions in H 2 and H 3
Here we assume that the particle is intrinsically moving in the manifold. Our heuristic approach to study such type of interactions on Riemannian manifolds is based on the idea of using the heat kernel as a regulator for point interactions on manifolds [13,14]. Thanks to the fact (24), the regularized interaction is chosen as the heat kernel on Riemannian manifolds. Once we have regularized the Hamiltonian, one can follow essentially the same steps outlined in the previous section, and obtain exactly the same form of the resolvent and principal matrix as in (28) and (29), respectively. In this paper, we only consider the particular class of Riemannian manifolds, namely two and three dimensional hyperbolic manifolds for simplicity. The heat kernel on hyperbolic manifolds of constant sectional curvature −κ 2 can be analytically calculated and given by [29] where d(x, y) is the geodesic distance between the points x and y on the manifold. The explicit form of the principal matrix in H 3 can then be easily evaluated [15]: Similarly, the principal matrix in H 2 can simply be evaluated by interchanging the order of integration with respect to t and s where ψ is the digamma function with its integral representation [28] for (w) > 0, and Q is the Legendre function of second type [28] with its integral representation for real and positive a and (α) > −1.
Since the spectrum of the free Hamiltonian in H n includes only the continuous part starting from (n−1) 2 κ 2 /4, it is natural to assume E i B < (n − 1) 2 κ 2 /4.

Two Types of Relativistic Extensions of Point-like Dirac delta Interactions
We first consider the so-called semi-relativistic Salpeter type free Hamiltonian (also known as relativistic spin zero Hamiltonian) perturbed by point like Dirac delta potentials in one dimension. This problem for the single center case has been first studied in [30] from the self-adjoint extension point of view. The formal Hamiltonian is exactly in the same form as in (10), except for the free part in the units where = c = 1. This non-local operator is a particular pseudo-differential operator and defined in momentum space as multiplication by p 2 + m 2 [31], which is known as the symbol of the operator. After following the renormalization procedure outlined above for the point interactions in two and three dimensions, the resolvent and the principal matrix is exactly the same form as in (28) and (29), respectively. However, the explicit expression of the heat kernel in this case is given by [31] K t (x, y) = mt where K 1 is the modified Bessel function of the first kind. Due to the short-time asymptotic expansion the diagonal term in the principal matrix (29) is divergent. In contrast to the one-dimensional case for point Dirac delta potentials, this problem therefore requires renormalization, as noticed by [18,32]. Choosing the coupling constants as in (27) by substituting the heat kernel (40) and taking the limit → 0 + , we obtain the resolvent in the form of the Krein's formula (11). The explicit form of the diagonal principal matrix is given by where Its off-diagonal part is given by where E i B is the bound state energy to the i th center in the absence of all the other centers. Since the spectrum of the free Hamiltonian includes only the continuous spectrum starting from m, it is natural to expect that E i B < m. An alternative relativistic model can be introduced from a field theory perspective in two dimensions. If we take very heavy particles interacting with a light particle, in the extreme limit of static heavy particles one recovers the following model: where a i refer to the locations of static heavy particles. Here where † denotes the adjoint. Since this model was worked out in [17], we will be content with the resulting formulae only referring to the original paper for the details. We can compute the diagonal principal matrix as and the off-diagonal part as for −m < (z) < m. Moreover, the binding energy of the single center should satisfy −m < E i B < m, and the lower bound is due to the stability requirement, to prevent pair creation to reduce the energy further thus rendering the model unrealistic in single particle sector.

Dirac delta Interactions supported by curves in R 2 and in R 3
We consider N Dirac delta potentials supported by non-intersecting smooth curves γ j : [0, L j ] → R n of finite length L j (n = 2, 3). Each curve is assumed to be simple, i.e., γ j (s 1 ) = γ j (s 2 ) whenever s 1 = s 2 , where s 1 , s 2 ∈ (0, L j ). Our formulation also allows the simple closed curves.
The Hamiltonian of the system is given by where γ i |r = Γi ds δ(r − γ i (s)). Then, the Schrödinger equation (H|ψ = E|ψ ) associated with this Hamiltonian is In contrast to the point-like Dirac delta interactions, this equation is a generalized Schrödinger equation in the sense that it is non-local. The resolvent kernel of the above Hamiltonian is explicitly given in the same form associated with point like Dirac delta potentials, namely where or if we express it in terms of the heat kernel Using the explicit form of the heat kernel in two dimensions, the above principal matrix becomes The spectrum of the free Hamiltonian includes only continuous spectrum starting from zero, so we expect that the bound state energies must be below z = 0. For this reason, we restrict the principal matrix to the negative real values, i.e., z = −ν 2 , ν > 0. Then, we have For non self-intersecting curve γ i , we can expand it around the neighbourhood of s i = s i in the Serret-Frenet frame at s i [33]: is the torsion of the curve). In the first approximation, keeping only the linear terms in s i − s i , and translating and rotating the Serret-Frenet frame attached to the coordinate system Oxy in such a way that t i (s i ) = (1, 0) and n i (s i ) = (0, 1), we have Then, the integral in the diagonal part of the principal matrix (55) around s i = s i in the first approximation is By making change of coordinates ξ i = (s i +si) 2 , the above integral becomes Using and the integrals of modified Bessel functions [34] ∞ 0 dx x n K 0 (ax) = 2 n−1 a −n−1 Γ 2 1 + n 2 , where n = 0, 1 and Γ is the gamma function, it is easy to see that the integral that we consider is finite around η i = 0 (s i = s i ). For non self-intersecting curves, the integrals in the diagonal and off-diagonal terms in (55) are finite whenever s i = s i due to the upper bounds of the Bessel functions [14] K 0 (x) In three dimensions, the Dirac delta potentials supported by curves requires the renormalization. Using the explicit formula of the heat kernel (23) for three dimensions, we find One can show that the the diagonal part of the above principal matrix (53) includes a term which is divergent around s i = s i . This can be immediately seen using the similar method outlined above, that is, the above integral includes the following integral in the new variable η i : which is divergent around η i = 0. Similar to the non-relativistic and relativistic point interactions, we first regularize the resolvent and then by choosing the coupling constant as a function of the cut-off parameter : and taking the formal limit → 0 + , we obtain the resolvent which is exactly the same form as in (51) except the matrix Φ is given by Here, E i B is the bound state energy of the particle to the delta interaction supported by ith curve in the absence of all the other delta interactions. Since the spectrum of the free Hamiltonian only includes the continuous part starting from zero, we have E i B < 0. Using the explicit form of the heat kernel, the principal matrix turns out to be a finite expression: A semi-relativistic generalization of particles interacting with curves is presented in [19]. The formal Hamiltonian can be written as We refer to this work for the details and we are content with writing down the resulting Φ matrix, since for tunneling corrections to the bound spectra this is all we need: As usual, these formulae must be analytically continued in z outside of their region of convergence. In our approach we are interested in the bound states for which these formulae are valid.

Analytic Structure of the Principal Matrices and the Bound State Spectrum
It is well-known that the bound state spectrum is determined by the poles of the resolvent, so the bound state spectrum should only come from the points z below the spectrum of the free Hamiltonian, where the matrix Φ is not invertible, i.e., the bound state energies are the real solutions of the equation where E < σ(H 0 ). Looking at the explicit forms of the principal matrices Φ ij (z), we see that they are all matrix-valued holomorphic functions (defined on their largest possible set of the complex plane). The analytical structure of the principal matrices can be determined by using the generalized Loewner's theorem [35], which simply states that if f 0 is a real valued continuously differentiable function on an open subset ∆ of (−∞, ∞), then the following are equivalent: • There exists a holomorphic function f with f ≥ 0 on the upper half-plane of the complex plane such that f has an analytic continuation across ∆ that coincides with f 0 on ∆.
• For each continuous complex valued function F on ∆ that vanishes off a compact subset of ∆, where for ζ, η ∈ ∆, For simplicity, let us explicitly show the analytical structure of the principal matrix associated with the Dirac delta potential supported by a single curve in two dimensions. In this case, the principal matrix (52) is just the diagonal part, say Φ(E), and continuously differentiable function of E, where E is on the negative real axis. Then, we have where ζ, η is on the negative real axis and L is the length of the curve Γ. Using the resolvent identity for the free resolvent, i.e., where R † 0 (η) = R 0 (η) = R 0 (η). The positivity is preserved in the limiting case ζ → η as well. This shows that the analytically continued function, sayΦ is a Nevallina function. We denote the analytically continued function by the same letter Φ for simplicity. The aforementioned theorem can be generalized to the matrix valued function Φ ij (E), as a result to ensure the holomorphicity we verify that: and the principal matrix in all the other cases, including the relativistic extension of these problems, can be similarly analyzed. Hence, for a large region of the complex plane, which contains the negative real axis, the principal matrix is a matrix-valued holomorphic function so that its eigenvalues and eigenprojections are holomorphic near the real axis [36]. In fact, we get poles on the real axis for the eigenvalues and the residue calculus can be used to calculate the associated projections. Let us consider the eigenvalue problem for the principal matrix depending on the real parameter E: where k = 1, 2, . . . , N and we assume there is no degeneracy for simplicity (we consider the generic case). In order to simplify the notation, we sometimes suppress the variable E in the equations, e.g., A k (E) = A k and so on. Then, the bound state energies can be found from the zeroes of the eigenvalues ω, that is, for each k. Thanks to Feynman-Hellmann theorem [37,38], we have the following useful result where ., . denotes the inner product on C N . Using the expression of the principal matrices in all class of singular interactions described above and using the positivity of the heat kernel, it is possible to show that This is an important result, since it implies that every eigenvalue cuts the real axis only once, that particular value gives us a bound state if it is below the spectrum of the free part. Moreover we deduce that the ground state energy corresponds to the smallest eigenvalue of Φ.

Off-Diagonal Terms of the Principal Matrices in the Tunneling Regime
For simplicity, we assume that all binding energies E i B 's or/and λ i 's are different. We consider the situation where the Dirac delta potentials (supported by points and curves) are separated far away from each other in the sense that the de Broglie wavelength of the particle is much smaller than the minimum distance d between the point Dirac delta potentials or than the minimum distance between the delta potentials supported by non intersecting regular curves with finite length, namely or in the semi-relativistic case, this can be stated as d λ Compton . This regime can also be defined in terms of the energy scales, namely where E B is the minimum of the binding energies to the single delta potentials in the absence of all the others (recall that = 2m = 1). In the non-relativistic problem for point interactions in one and three dimensions, it is clear from the explicit form of the principal matrices (18), (33) all the off-diagonal terms are getting exponentially small as d increases, i.e., and as d → ∞. For point interactions in two dimensions, thanks to the upper bound of the Bessel function [14], for all x, the off-diagonal terms of the principal matrix (32) is going to zero exponentially as d → ∞. In the above expressions for principal matrices, we have expressed them in terms of a real positive variable ν for simplicity. Not all the bound state spectra of the potentials we consider in this paper are negative, so it is not always useful to express the principal matrix in terms of a real positive variable ν. For that purpose, we will consider the principal matrices restricted to the real values, namely z = E, where E is the real variable (not necessarily negative). For point interactions in three dimensional hyperbolic manifolds, the off-diagonal principal matrix restricted to the real values E < κ 2 is exponentially small as d → ∞. Here d is the minimum geodesic distance between the centers. As for the point interactions in two dimensional hyperbolic manifolds, the off-diagonal principal matrix restricted to the real values E < κ 2 /4 becomes Using the series representation of the Legendre function of second kind [28] Q where v = 1 2 + − E κ 2 + 1 4 > 1 and α = κd(a i , a j ), and splitting the sum, we obtain Since Gamma function is increasing on [2, ∞), < 1 for all k ≥ 2, and v > 1, we can find an upper bound for the above the infinite sum as which is simply a geometric series. All these show that the off-diagonal principal matrix in two dimensional hyperbolic manifolds is exponentially small as d → ∞ and the leading term is given by the first term of the series expansion. As for the delta interactions supported by curves, the minimum of the pairwise distances between the supports of Dirac delta potentials always exists since d ij (s, s ) = |(γ i (s) − γ j (s ))| is a continuous function on compact interval s ∈ [0, L], so we have for i = j. Then, Due to the upper bound of the Bessel function (85), the off-diagonal principal matrix is going to zero as d → ∞.
Similarly, the explicit forms of the off-diagonal parts of the principal matrices (44) and (48) in the relativistic cases are exponentially going to zero as d → ∞ (by assuming the order of the limit and the integral can be interchanged). For the other relativistic cases (including the relativistic delta potentials supported by curves), the off diagonal terms of the principal matrices can also be shown to be exponentially small. Therefore, we see that the principal matrices for all the above models are diagonally dominant in the "large" separation regime. However, the exponentially small off-diagonal terms are not analytic in the small parameter ( 1 E B d 2 ). Nevertheless, we can keep track of small values of the off-diagonal terms by introducing an artificial parameter in order to control the orders of terms in the perturbative expansion, that we are going to develop in the next section.

Splitting in Bound State Energies through Perturbation Theory
Let us consider the family of principal matrices restricted to the real axis E: where Φ 0 is the diagonal part of the principal matrix, and δΦ is off-diagonal part of it and this is the "small" correction (perturbation) to the diagonal part. Since Φ(E) is symmetric (Hermitian), we can apply standard perturbation techniques to the principal matrix [36,39,40]. For this purpose, let us assume we can expand the eigenvalues and eigenvectors as follows: for each k.
The solution to the related unperturbed eigenvalue problem is given by Once we have found the eigenvalues and eigenvectors of the diagonal part of the principal matrix or unperturbed eigenvalue problem, we can perturbatively solve the full problem. The standard perturbation theory gives us the eigenvalues ω k up to second order: and the first order correction to the eigenvectors A k is given by Since the bound state energies are determined from the solution of equation (78), the bound state energies in the zeroth order approximation can easily be found from ω k 0 (E) = 0. The solution is given by and the corresponding eigenvector is where 1 is located in the kth position of the column and other elements of it are zero or we can write Here e k i s form a complete orthonormal set of basis.
The bound state energies to the full problem up to the second order is then determined by solving the following equation where we have used the first order result ω k 1 = 0 (106) from the Equation (98).
Let us now expand ω k 0 (E) and Φ kl (E) for k = l around E = E k B : where ω k 0 (E k B ) = 0. If we substitute (107) into (105) and (99), and use Feynman-Hellman theorem given in previous section, the condition (105) up to the second order turns out be If we also expand the last factor in the powers of (δE k ) and ignore the second order terms and combine the terms using the symmetry property of principal matrix, we find Ignoring the second and third terms on the left hand side of the equality (this is guaranteed by the assumption Φ kk (E k B ) |Φ kl (E k B )|) and setting = 1, we get the change in E k (to first order) as, This is our main formula for all types of singular interactions we consider. It is striking that it contains the information about the tunneling regime.

Explicit Examples for the Splitting in the Energy
Let us now compute explicitly how the bound state energies change in the tunneling regime for the above class of singular potentials. For point Dirac delta potentials in one dimension, the bound state energies are negative so E k B = −|E k B | and in the tunneling regime d |E B | 1. For point Dirac delta potentials in two dimensions, the bound state energies are negative and again in the tunneling regime. Here we have used the asymptotic expansion of the modified Bessel function of the third kind K 0 (x) ≈ π 2x exp(−x) for x 1 [28]. In three dimensions, we have For point interactions in three dimensional hyperbolic manifolds, the bound state energies are below κ 2 (see [15] for details) and in the tunneling regime. Here we have used sinh 2 x ≈ e 2x 4 as x 1. For point interactions in two dimensional hyperbolic manifolds, the bound state energies are below κ 2 /4 (see [15]) and where ψ (1) is the polygamma function and we have used the infinite series representation of the Legendre function of second kind (89). For semi-relativistic point interactions in one dimensions, the bound state energies are below m. Let us first find explicitly integrals in the off-diagonal part of the principal matrix asymptotically in the tunneling regime md 1. For this purpose, let us rescale the integration variable s = µ/m so that the above integral becomes m 2 π ∞ 1 e −sm|a k −a l | √ s 2 −1 m 2 (s 2 −1)+(E k B ) 2 . Note that −s in the exponent has its maximum at s = 1 on the interval (1, ∞). Then, only the vicinity of s = 1 contributes to the full asymptotic expansion of the integral for large m|a k − a l |. Thus, we may approximate the above integral by m 2 π 1 e −sm|a k −a l | √ s 2 −1 m 2 (s 2 −1)+(E k B ) 2 , where > 1 and replace the function √ s 2 −1 m 2 (s 2 −1)+(E k B ) 2 in the integrand by its Taylor expansion [41]. It is important to emphasize that the full asymptotic expansion of this integral as m|a k − a l | → ∞ does not depend on since all other integrations are subdominant compared to the original integral. Hence, we find where we have used the fact that the contribution to the integral outside of the interval (1, ) is exponentially small. Substituting this result into Eq. (110), we find when E k B < 0 and For the field theory motivated relativistic version we can use a saddle point approximation, assuming that tunneling condition, given by m 2 − (E k B ) 2 d ij >> 1 is satisfied. Here it is enough to consider the function m(1 + s 2 ) 1/2 − E k B s and expand it around the maximum E k B / m 2 − (E k B ) 2 . The denominator can be replaced by its value at the maximum, we find that the leading behaviour goes as (assuming that E k B d ij 's remain large) evaluating the integral we end up with, Once we obtain the off-diagonal terms responsible for the tunneling contributions, calculating the derivatives of the diagonal parts are simple, Substituting these expressions into the general formulae we have derived, gives the tunneling contribution to energy levels that leads to small shifts in the binding energies. For Dirac delta potentials supported by curves in two dimensions: we define a kind of center of mass by and write |γ( in the argument of the functions in the principal matrix. When we evaluate the expressions we expand these terms by keeping only first order terms in the small quantities. The resulting Bessel functions can be expanded again to find the leading corrections for the curve to curve interaction terms. We use the expression above for the off diagonal terms and define d ij = |x i − x j | for simplicity and introduce a unit vector asd ij in a similar way. As a result we have the leading order expansion, When we insert this into Φ ij expression and integrate over the curve, we find and similarly for the other part. Thus we see that the only contribution comes from the second order which we neglect for our purposes. However a systematic expansion in powers of 1 dij can be developed for higher order correction as described. Using the asymptotic expansion of K 0 (z) for large values of z [28], for all ν ≥ 0 we get from (110) a more elegant expression, where Φ ll and its derivative at E k B can be computed from the explicit expression of the principal matrix (54). For Dirac delta potentials supported by curves in three dimensions, there is really no change, since renormalization is required only for the diagonal parts, we have the off-diagonal expressions already in a simpler form, as a result of the above analysis, the leading order expression is found to be, where Φ ll and its derivative at E k B can be computed from the explicit expression of the principal matrix (67). In a similar way, we look at the tunneling correction to bound state energies for relativistic particle coupled to Dirac potentials supported over curves. Again we use the approximation that the separation of the curves are large and the extend of the curves compared to these distances are small. This is not the only possible approximation, one can envisage a situation in which the separations are large but the extend of the curves are also large. The essential ideas are captured by our example so to achieve technical simplicity we keep this approximation. Essential point is to expand the off-diagonal terms in the leading order. By scaling t variable in the integral we can write Φ ij (E k B ) term as, where in the second line we used the asymptotics of K 1 for large argument (127). We may now use the same argument by means of the center of mass of the curves to define center to center distances and expand around the center of mass, not surprisingly we again find that the first order corrections become zero, only the center to center distance matters. Therefore, to leading order we have a simpler expression, This is of the type we have worked out for the semi-relativistic particle, and in the same manner, a saddle point approximation can be applied in a simple way, resulting We may now employ our general expressions to find the tunneling corrections. The derivative of the diagonal term can be simplified by means of ∂K0(z) ∂z = −K 1 (z).

Degenerate Case and Wave Functions for Point Interactions
Let us now compute the energy splitting of two equal strength delta functions supported by the points −a and a in two dimensions. This is very similar to the double well problem we discuss in the introduction, yet this version can be solved exactly. The approximation we use corresponds to the standard WKB approach. Let us recall that when we have two degenerate eigenvalues the degeneracy is lifted by the diagonal perturbation and as is well known, diagonalizing the perturbation matrix in the degeneracy subspace gives us the first order correction: If we call the common bound state as E B , for k = 1, 2 to get the first order correction we truncate the eigenvalue equations as, where we have used the asymptotic expansion of K 0 given by (127). Thus the splitting is given by which should be compared with the usual one-dimensional double well potential splitting given in the introduction. Note that in the former case, the strength of each harmonic well is proportional to the square of the separation therefore the initial energy level is not independent as in the delta function case and is proportional to the square of the separation. the exponent thus gets the square of the distance as the suppression factor, if we assume that E B ∼ |a| 2 one can see that the exponents behave in a similar way. Actually, one can also compare the first order perturbation result for the splitting δE 1 with the numerical result by solving det Φ(ν) = ln(ν/µ) − ±K 0 (2aν) = 0 numerically for each a by Mathematica (see Figure 2). We assume that a > e γ in order to guarantee the existence of the second bound states, where γ is the Euler's constant.
The same method can also be applied to the one-dimensional case. In the symmetrically placed Dirac delta potentials with equal strengths λ, the exact bound state energies when they are sufficiently far away from each other (when a > 1/λ, there are two bound state energies) can analytically be computed [42] where W is the Lambert W function [43], which is defined as the solution y(x) of the transcendental equation ye y = x. From (17), the principal matrix in this case reads Then, the first order perturbation result following the above procedure gives where we have used well-known result E B = − λ 2 4 . Then, one can easily find the error between the exact result δE exact = E + − E − and the first order perturbation result δE 1 in the splitting of the energy, see the Figure 3.
The three dimensional case can also be studied in this way and we can similarly solve det Φ(ν) = (ν − µ) − ± 1 2a e −2aν in terms of the Lambert W function and compare with the first order order perturbation result for the splitting in the energy (Figure 4): Here we assume that a > 1/2µ in order to guarantee the existence of second bound states. Let us emphasize that in the usual WKB approach one constructs the wave functions in classically allowed and forbidden regions respectively and use a subtle argument to connect the different regions. In this case, there is really no forbidden region, except the supports of the attractive regions. Indeed right there classically there is no sensible way to define the motion of a particle. Nevertheless, it is possible to find the effect of tunneling for the wave functions from our formalism. It relies on the first order corrections to the eigenstates of the principal operator, notice that an expansion of the eigenstates of the principal operator can be found in the non-degenerate case as Note that to this order normalization is not important, moreover we do not need to use a subtle argument about the shift of the eigenvalues since the change of eigenvalue is already second order in the exponentially small quantities, any such correction will be of lower order as we have seen in the shift of energy calculations. It is well-known that the wave function of the system associated with the bound states can be found from the explicit expression of the resolvent formula. Since the eigenvalues are isolated we can find the projections onto the subspace corresponding to this eigenvalue by the following contour integral (Riesz Integral representation) [40]: where C k is a small contour enclosing the isolated eigenvalue, say E k . We note that the free resolvent does not contain any poles on the negative real axis for the Dirac delta potentials supported by points, so all the poles on the negative real axis will come from the poles of inverse principal matrix Φ −1 (z). Since the principal matrix is self-adjoint on the real axis, we can apply the spectral theorem. Moreover, its eigenvalues and eigenprojections are holomorphic near the real axis, as emphasized in section 3. Then, we can write the spectral resolution of the inverse principal matrix, where P k (z) ij = A ik (z)A jk (z), A ki (z) is the normalized eigenvector corresponding to the eigenvalue ω k (z). Then, from the residue theorem, we find the square integrable wave function associated with the bound state energy E k as where α = (− ∂ω k ∂E E k ) −1/2 is the normalization constant. This is actually a general formula for the bound state wave function for the Dirac delta potentials supported by points in R n . For n = 2, we have Let us recall that the eigenstates for the unperturbed levels are given by unit vectors (103), we write this into the formula for the wave function (145). As a result, using the first order correction (100) to the eigenstate A k where we use 1 This form of the wave function clearly shows the tunneling nature of the solutions. It is now quite straightforward to compute the wave functions in this approximation for all the other cases we consider.

Conclusion
In this paper, we have first reviewed the basic results about some singular interactions, such as the Dirac delta potentials supported by points on flat spaces and hyperbolic manifolds, and delta potentials supported by curves in flat spaces. Moreover, the results in the relativistic extensions of the above-mentioned potentials have been also reviewed which was essentially given in [13,14,15,17,18,19]. The main result of this paper is to develop some kind of perturbation theory applied to a class of singular potentials in order to find the splitting in the energy due to the tunneling. This was only developed extensively for Dirac delta potentials supported by points in [14], here we extend the method to various kinds of Dirac delta potentials as well as to their relativistic versions. It is possible to give some bounds over the error terms if we assume that the errors in perturbation theory can be estimated. Typical perturbative expansions are asymptotic therefore a truncation is needed to get more accurate results, one knows that it gets worse beyond a few terms. The more accurate thing to do is to obtain a Borel summed version but that is beyond the content of the present paper, it will depend very much of the specifics of the model whereas we prefer to give a broader perspective.
The comparison with conventional methods certainly would be very useful, nevertheless at present we do not know how a more conventional approach, such as WKB or instanton calculus can be performed in these singular problems. Since the potentials are localized at points or along the curves, the variation of the potential relative to any wavelength is always much more important. Indeed this unusual behavior changes the problem completely. We need to give a meaning to these potentials first and redevelop the WKB analysis. Our main point here is that in this description of the singular potentials via resolvents, the WKB's reincarnation is given by a perturbative analysis of the eigenvalues of the principal operator for large separations of the supports.