Uniform Asymptotic Approximation Method with P¨oschl-Teller Potential

,


I. INTRODUCTION
A century after the first claim by Einstein that general relativity (GR) needs to be quantized, the unification of Quantum Mechanics and GR still remains an open question, despite enormous efforts [1].Such a theory is necessary not only for conceptual reasons but also for the understanding of fundamental issues, such as the big bang and black hole singularities.Various theories have been proposed and among them, string/M-Theory and Loop Quantum Gravity (LQG) have been extensively investigated [2,3].Differences between the two approaches are described in [4,5].
LQG was initially based on a canonical approach to quantum gravity (QG) introduced earlier by Dirac, Bergmann, Wheeler, and DeWitt [6].However, instead of using metrics as the quantized objects [6], LQG is formulated in terms of densitized triads and connections, and is a non-perturbative and background-independent quantization of GR [7].The gravitational sector is described by the SU(2)-valued Ashtekar connection and its associated conjugate momentum, the densitized triad, from which one defines the holonomy of Ashtekar's connection and the flux of the densitized triad.Then, one can construct the full kinematical Hilbert space in a rigorous and well-defined way [3].An open question of LQG is its semiclassical limit, that is, are there solutions of LQG that closely approximate those of GR in the semiclassical limit?
Although the above question still remains open, concrete examples can be found in the context of loop quantum cosmology (LQC) (For recent reviews of LQC, see [8][9][10][11][12][13][14][15][16][17] and references therein).Physical implications of LQC have also been studied using the effective descriptions of the quantum spacetimes derived from coherent states [18], whose validity has been verified numerically for various spacetimes [19,20], especially for states sharply peaked on classical trajectories at late times [21].The effective dynamics provide a definitive answer on the resolution of the big bang singularity [22][23][24][25][26][27], replaced by a quantum bounce when the energy density of matter reaches a maximum value determined purely by the underlying quantum geometry.
One of the major challenges in the studies of cosmological perturbations in LQC is how to solve for the mode functions µ k from the modified Mukhanov-Sasaki equation.So far, it has mainly been done numerically [8][9][10][11][12][13][14][15][16][17].However, this is often required to be conducted with highperformance computational resources [41], which are not accessible to general audience.
In the past decade, we have systematically developed the uniform asymptotic approximation (UAA) method initially proposed by Olver [42][43][44], and applied it successfully to various circumstances [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]  1 .In this paper, we shall continuously work on it by considering the case in which the effective potential has only zero points but without singularities.To be more concrete, we shall consider the Pöschl-Teller (PT) potential, for which analytical solutions are known [66].The consideration of this potential is also motivated from the studies of cosmological perturbations in dressed metric and hybrid approaches [67,68], in which it was shown explicitly that the potentials for the mode functions can be well approximated by the PT potential with different choices of the PT parameters.In particular, in the dressed metric approach, the mode function satisfies the following equation [67] in which V (η) serves as an effective potential.During the bouncing phase it is given by where γ B is a constant introduced in [67], and m Pl and t Pl are respectively, the Planck mass and time.This potential can be well approximated by a PT potential with Here η is the conformal time related to the cosmic time t by dη = dt/a(t).On the other hand, in the hybrid approach, the effective potential during the bouncing phase is given by which can be also modeled by the PT potential (1.3) but now with [68] V For more details, we refer readers to [67,68].The rest of the paper is organized as follows: In Sec.II we provide a brief review of the UAA method with two turning points, and show that the first-order approximate solution will be described by the parabolic cylindrical functions.In Sec.III we construct the explicit approximate analytical solutions with the PT potential, and find that the parameter space can be divided into three different cases: A)

and C)
to cosmological perturbations in GR was carried by Habib et al, [64,65].
where k and β are real constants.After working out the error control function T [cf.Appendix C] in each case, we are able to determine the parameter q 0 , introduced in the process of the UAA method in order to minimize the errors.Then, we show the upper bounds of errors of our approximate solutions with respect to the exact one, given in Appendix B. In particular, in Case A), the upper bounds are ≲ 0.15%, while in Case B) they are no larger than 10%.In Case C), the errors are also very small, except the minimal points [cf.Fig. 10], at which the approximate solutions deviate significantly from the analytical one.The causes of such large errors are not known, and still under our investigations.In each of these three cases, we also develop our numerical codes, and find that the numerical solutions trace the exact one very well, and the upper bounds of errors are always less than 10 −4 % in each of the three cases.The paper is ended in Sec.IV, in which our main conclusions are summarized.There are also three appendices, A, B, and C, in which some mathematical formulas are presented.

II. THE UNIFORM ASYMPTOTIC APPROXIMATION METHOD
Let us start with the following second-order differential equation It should be noted that all second-order linear homogeneous ordinary differential equations (ODEs) can be written in the above form by properly choosing the variable y and µ k (y).Instead of working with the above form, we introduce two functions g(y) and q(y), so that the function f (y) takes the form [42] f (y) = λ 2 g(y) + q(y), where λ is a large positive dimensionless constant and serves as a bookmark, so we can expand µ k (y) as After all the calculations are done, one can always set λ = 1 by simply absorbing the factor λ −n into µ (n) k (y).It should be noted that there exist cases in which the above expansion does not converge, and in these cases we shall expand µ k (y) only to finite terms, say, N , so that µ k (y) is well approximated by the sum of these N terms.On the other hand, the main reason to introduce two functions g(y) and q(y), instead of only f (y), is to minimize errors by properly choosing g(y) and q(y).
In general, the function g(y) has singularities and/or zeros in the interval of our interest.We call the zeros and singularities of g(y) as turning points and poles, respectively.The uniform asymptotic approximate (UAA) solutions of µ k (y) depend on the properties of g(y) around their poles and turning points [42][43][44].The cases in which g(y) has both poles and turning points were studied in detail in [47,49,54], so in this paper we shall focus ourselves on the cases where singularities are absent and only turning points exist.As to be shown below, the treatments of these cases will be different from the ones considered in [47,49,54].In particular, in our previous studies the function q(y) was uniquely determined by requiring that the error control function be finite and minimized at the poles, while in the current cases no such poles exist.So, to fix q(y), other analyses of the error control function must be carried out.

A. The UAA Method
The UAA method includes three major steps: (i) the Liouville transformations; (ii) the minimization of the error control function; and (iii) the choice of the function y(ζ), where ζ is a new variable.In the following, we shall consider each of them separately.

The Liouville Transformations
The Liouville transformations consist of introducing a new variable ζ(y), which is assumed that the inverse y = y(ζ) always exists and is thrice-differentiable.Without loss of the generality, we also assume that y(ζ) is a monotonically increasing function [cf.Fig. 1].Then, in terms of U (ζ), which is defined by Eq.(2.1) takes the form, where and It should be noted that Eqs.(2.1) and (2.5) are completely equivalent, and so far no approximations are taken.However, the advantage of the form of Eq.(2.5) is that, by properly choosing q(y), the term |ψ(ζ)| can be much smaller than λ 2 ẏ2 g , that is, so that the exact solution of Eq.(2.1) can be well approximated by the first-order solution of Eq.(2.5) with ψ(ζ) = 0.This immediately raises the question: how to choose q(y) so that the condition (2.8) holds.To explain this in detail, let us move onto the next subsection.

Minimization of Errors
To minimize the errors, let us first introduce the error control function [42-44, 47, 49, 54] Then, introducing the free parameters a n and b n into the functions g(y) and q(y), so we have where n = 1, 2, ..., N , with N being an integer.It is clear that for such chosen g(y) and q(y), the error control function T (ζ) will also depend on a n and b n .To minimize the errors, one way is to minimize the error control function by properly choosing a n and b n , so that On the other hand, the errors also depend on the choice of y(ζ), which in turn sensitively depends on the properties of the functions g(y) and q(y) near their poles and turning points.In addition, it must be chosen so that the resulting equation of the first-order approximation (obtained by seting ψ(ζ) = 0) can be solved explicitly (in terms of known functions).Considering all the above, it has been found that y(ζ) can be chosen as [42-44, 47, 49, 54] in the cases with zero, one and two turning points, respectively.Here sgn(g) = 1 for g > 0 and sgn(g) = −1 for g < 0.
In the rest of this paper, we shall consider only the cases with two turning points.

B. UAA Method for Two Turning Points
For the cases with two turning points, we can always write g(y) as where y 1 and y 2 are the two turning points, and p(y) is a function of y with p(y i ) ̸ = 0, (i = 1, 2).In general, according to the properties of y 1 and y 2 , we can divide all the cases into three different subclasses: 1. y 1 and y 2 are two distinct real roots of g(y) = 0; 2. y 1 = y 2 , a double real root of g(y) = 0; and 3. y 1 and y 2 are two complex roots of g(y) = 0. Since g(y) is real, in this case these two roots must be complex conjugate, y 1 = y * 2 .To apply the UAA method to Eq.(2.6), we assume that the following conditions are satisfied [47,49,54]: • When far away from any of the two turning points, we require q(y) g(y) ≪ 1. (2.14) • When near any of these two points, we require provided that the two turning points are far away from each other, that is, when • If the two turning points are close to each other, |y 1 − y 2 | ≃ 0, then near these points we require It should be noted that, when |y 2 − y 1 | ≫ 1, the two turning points are far away, and each of them can be treated as an isolated single turning point [42,43].In addition, without loss of generality, we assume that g(y) < 0 for y > y 2 or y < y 1 , when y 1 and y 2 are real.When y 2 and y 1 are complex conjugate, we assume that g(y) < 0 [cf.Fig. 2].Then, in this case we adopt a method to treat all these three classes listed above together [44,47,49,54].In particular, we choose ẏ2 g as (2.17) FIG. 2. The function g(y) defined by Eq. (3.3) for different choices of k and β.In particular, the dotted black line denotes the case k 2 < β 2 , and the solid blue line denotes the case k 2 = β 2 , while the dash-dotted red line denotes the case so that ζ is an increasing function of y [cf.Fig. 1] and When we integrate the above equation, without loss of the generality, we shall choose the integration constants so that Then, we find that > 0, y 1,2 real, and y 1 ̸ = y 2 , = 0, y 1,2 real, and with where " + " corresponds to the cases that the two turning points y 1 and y 2 are both real, and " − " to the cases that the two turning points y 1 and y 2 are complex conjugate.When y 1 and y 2 are complex conjugate, the integration of Eq.(2.21) is along the imaginary axis [44].When the two real roots are equal, we have ζ 0 = 0. To proceed further, let us first derive the relation between ζ(y) and y by first integrating the right-hand side of Eq.(2.18).To this goal, it is found easier to distinguish the case in which y 1 and y 2 are real from the one in which they are complex conjugate.
(2.18) we find Now let us turn to consider the case when y 1 and y 2 are complex.For this case ζ 2 0 is always negative, ζ 2 0 < 0, thus from Eq. (2.4) we find [44]

The First-order Approximate Solutions
With the choice of Eq.(2.17), we find that Eq. (2.6) reduces to where we assume that ζ ∈ (−ζ 2 , ζ 2 ), with ζ 2 being a real and positive constant, which can be arbitrarily large Neglecting the ψ(ζ) term, we find that the approximate solutions can be expressed in terms of the parabolic cylinder functions W ( 1 2 λζ 2 0 , ± √ 2λζ) [44], and are given by from which we have where α k and β k are two integration constants, ϵ 1 and ϵ 2 are the errors of the corresponding approximate solutions, whose upper bounds are given by Eqs.(A.1) and (A.2) in Appendix A.
For the choice of Eq.(2.17), we find that the associated error control function defined by Eq.(2.9) now takes the form for g < 0, and for g > 0.

III. UAA SOLUTIONS WITH THE P ÖSCHL-TELLER POTENTIAL
To study the case in which only turning points exist, in this paper we consider the second-order differential equation (2.1) with a Pöschl-Teller (PT) potential [67,68] as in this case exact solutions exist, where k is the comoving wavenumber, and β 0 is a real and positive constant.The two parameters β 0 and α determine the height and the spread of the PT potential, respectively.Under the rescaling αy → y, the α parameter can be absorbed into the wavenumber k and β 0 by redefining (k/α → k, β 0 /α → β 0 ).As a result, there is no loss of generality to set α = 1 from now on.Then, the exact solutions in this case exist, and are presented in Appendix B.
On the other hand, to apply the UAA method to this case, and to minimize the errors of the analytic approximate solutions, we tentatively choose q as where q 0 is a free parameter, to be determined below by minimizing the error control function (2.9) with the choice of ẏ2 g given by Eq.(2.17).Then, we have where β ≡ β 2 0 − q 2 0 .In this paper, without loss of generality, we shall choose q 0 so that β is always real, that is Thus, from g(y) = 0 we find that the two roots are given by It is clear that, depending on the relative magnitudes of β 0 and k, as well as the choices of q 0 , two turning points can be either complex or real.In Fig. 2, we plot out the three different cases, k 2 < β 2 , k 2 = β 2 , and k 2 > β 2 , from which it can be seen clearly that the two turning points are real and different for k 2 < β 2 , real and equal for k 2 = β 2 , and complex conjugate for k 2 > β 2 , respectively.Then, from Eqs.(3.2) and (3.3), we find that for |y| ≫ 0, and for |y| ≃ |y i | and |y 1 − y 2 | ≫ 1, and for |y| ≃ |y 1 | and |y 1 − y 2 | ≃ 0. In the following, let us consider the three cases: (a) k 2 ≫ β 2 ; (b) k 2 ≃ β 2 ; and (c) In this case, we have g(y) is always negative, g(y) < 0, so that the two turning points of g(y) = 0 are complex conjugate and are given by As discussed in the last section, now ζ 2 0 < 0, for which Eq.(2.26) can be cast in the form where ζ2 0 ≡ −ζ 2 0 > 0. Note that in writing down the above equation, we had replaced U by W .In addition, the new variable ζ is related to y via from which we find that ζ0 is given explicitly by Moreover, in the case of the PT potential, the integration of Eq. (3.11) can be carried out explicitly, giving where ϵ y denotes the sign of y with x ≡ 1/ cosh(y), and AppellF 1 is the Appell hypergeometric function.Ignoring the ψ term in Eq. (3.10), we find the general solution where W denotes the Weber parabolic cylinder function [69], and α k and β k are two integration parameters which generally depend on the comoving wavenumber k.
The validity of the analytic solution (3.14) depends on the criteria given by Eqs.(2.14) -(2.16), while its accu-racy can be predicted by the error control function T .In the current case, we find that T of Eq.(2.29) can be written as a combination of three terms as that given by Eqs.(C.3) 2 , where where A is given by Eq.(C.5).It should be noted that T 1 , T 2 and T 3 given in Eq. (3.15) all vanish when y = 0 (for which we have x = 1 and ζ = 0), that is, Besides, as the PT potential is an even function, the er-ror control function is antisymmetric about the origin, namely, T (−y) = −T (y).As a result, we will study its behavior only on the positive y axis, y ≥ 0. With the help of Eq. (3.11), the numeric value of the error control function at any point y > 0 can be found from Eq.(3.15).In particular, for β/k ≪ 1, we find that  In Fig. 3 we plot the functions, |q/g|, |q(y − y 1 )/g|, and |q(y − y 1 )(y − y 2 )/g|, together with the error control function defined by Eqs.(C.3) -(C.5) for (k, β) = (5.0,1.0), with q 0 being given by Eq.(3.18).(Recall . From these figures it is clear that the conditions (2.14) -(2.16) are well satisfied, and the error control function remains small all the time.In particular, it decreases as β/k decreases.
In Fig. 4, we plot the mode functions µ k (y), µ N k (y), µ E k (y), and the relative difference δ A (y) defined by where A = (N, E), µ k (y) denotes the mode function ob-tained by the UAA method given by Eq.(3.14), µ N k (y) is the numerical solution obtained by integrating Eq.(2.1) directly with the same initial conditions, while µ E k (y) is the exact solution given by Eq.(B.5).From these figures we can see that the maximal errors occur in the region near y = 0, but the upper bound is no larger than 0.15% at any given y, including the region near y ≃ 0.
It is interesting to note that this analytical approximate solution is only up to the first-order approximation of the UAA method.With higher order approximations, the relative errors are even smaller.
To check our numerical solutions, in Fig. 4, we also plot the relative differences ϵ(y) between µ N k (y) and µ E k (y), defined by From these figures it can be seen that ϵ(y) is no larger than 10 −7 , and our numerical code is well tested and justified.
It is also interesting to note that the mode functions are oscillating for y ≲ −10, and these fine features are captured in all three mode functions, although there are some differences in the details.Again, as shown by their relative variations, these differences are very small.In addition, we also consider other choices of β and k, and find that they all have similar properties, as long as the condition k 2 ≫ β 2 .In this case, depending on k ≳ β or k ≲ β, the function g(y) has different properties, as shown in Fig. 2. Therefore, in the following subsections let us consider them separately.

k ≳ β
When k ≳ β the function g(y) is always non-positive for y ∈ (−∞, ∞).Then, from Eqs.(C.3) and (3.15) we find that as y → ∞, but now with ϵ ≡ (k − β)/k.Thus, to have the error control function be finite at y = ∞, now we must set instead of the value given by Eq.(3.18) for the case k ≫ β 2 .In Fig. 5, we plot the quantities |q/g|, |q(y − y 1 )/g|, |q(y − y 1 )(y − y 2 )/g|, and the error control function T for k = 5.0, β = 4.9 and q 0 = 1/2, for which we have y 1 = y * 2 = 0.200335i.From these figures we can see clearly that the conditions (2.14) -(2.16) are well satisfied, and the error control function remains small all the time.Then, the corresponding quantities µ k (y), µ N k (y), µ E k (y), δ A (y) and ϵ(y) are plotted in Fig. 6.From the curves of δ N (y) and δ E (y) we can see that now the errors of the first-order UAA solution are ≤ 4%, which are larger than those of the last subcase.This is mainly because of the fast oscillations of the solution in the region y < 0. Therefore, in order to obtain solutions with high precision, high-order approximations for this case are needed.However, we do like to note that our numerical solution still matches to the exact one very well, as shown by the curve of ϵ(y), which is no larger than 6.0 × 10 −6 .ϵ FIG. 6. Plots of the mode functions µ k (y), µ N k (y), µ E k (y) and their relative differences δ N (y), δ E (y) and ϵ(y) for k = 5, β = 4.9 and q0 = 1/2.

k ≲ β
In this case, we find that On the other hand, from Eqs.(2.30), (C.3) and (C.6) we find that where ζ(0) ≡ ζ(y)| y=0 < ζ 0 .Note that in calculating the error control function near the turning point y ≃ y 2 , we have used the relation so that the divergence of the second term of T 2 cancels exactly with that of T 3 .Eq.(3.25) can be obtained directly from the relation Similarly, it can be shown that It is clear that to minimize the errors, in the present case q 2 0 must be also chosen to be as that given by Eq. (3.22).In Fig. 7 Then, the corresponding quantities µ k (y), µ N k (y), µ E k (y), δ A (y) and ϵ(y) are plotted in Fig. 8. From the curves of δ N (y) and δ E (y) we can see that now the errors of the first-order UAA solution are ≲ 10%.Similar to the last subcase, this is mainly because of the fast oscillations of the solution in the region y < 0. Therefore, in order to obtain high precision, high-order approximations for this case are needed, too.In addition, our numerical solution still matches well to the exact one, as shown by the curve of ϵ(y), which is no larger than 2.0 × 10 −6 .
In this case, two real turning points appear, given, respectively by Then, we find that Eqs.(3.23) and (3.24) still hold in the current case, while Eq.(3.25) is replaced by as y → ∞, but now with ϵ ≡ k/β.Combining Eqs.(3.23), (3.24) and (3.29), we find that currently the proper choice of q 0 is still that given by q 0 = 1/2, as those given in the last two subcases.
In Fig. 9, we plot the quantities |q/g|, |q(y − y 1 )/g|, |q(y − y 1 )(y − y 2 )/g|, and the error control function T for k = 0.6, β = 4.0 and q 0 = 1/2, for which we have y 1 = −y 2 ≃ −2.58459.From this figure we can see that the preconditions (2.14)-(2.16)are well satisfied.Then, to the first-order approximation of the UAA method, the solution can be approximated by Eq.(2.28),where ζ 2 0 is given by Eq.(3.23), α k and β k are two integration constants, and ϵ 1 and ϵ 2 the errors of the corresponding approximate solutions, whose upper bounds are given by Eqs.(A.1) and (A.2) in Appendix A. In Fig. 10 (a), we plot out our first-order approximate solution, while Fig. 10 (b) to compare the approximate solution with the exact one, we plot both of them.In particular, the solid line represents the exact solution, while the red dotted line the approximate solution.From this figure it can be seen that except the minimal points, the two solutions match well.However, at these extreme minimal points, they deviate significantly from each.The causes of such errors are not clear, and we hope to come back to this issue in another occasion.
Finally, similar to all other cases, our numerical solution still matches well to the exact one, as shown by the curve of ϵ(y), which is no larger than 8.0 × 10 −6 .

IV. CONCLUSIONS
In this paper, we have applied the UAA method to the mode function µ k with a PT potential, for which it satisfies the second-order differential equation where k and β 0 are real constants.In this case, the exact solution is known and given by Eq.(B.5).The implementation of the UAA method includes the introduction of an auxiliary function q(y), which is taken as where q 0 is a free parameter.Then, we carry out the integration of the error control function, defined by where Clearly, the error control function T (ζ) will depend on q 0 .After working out the details, we find that it is convenient to distinguish the three cases: A) , where β 2 ≡ β 2 0 − q 2 0 > 0. In particular, in Case A), a proper choice of q 0 is q 0 = 1/ √ 24, while in Cases B) and C), it is q 0 = 1/2.Once q 0 is fixed, the analytical approximate solutions are uniquely determined by the linear combination of the two parabolic cylindrical functions W (ζ 2 0 /2, ± √ 2ζ), as shown by Eq.(2.28).In particular, in Case A) the upper bounds of errors are ≲ 0.15%, as shown in Fig. 4. In Case B), two subcases are considered, one with k ≳ β, and the other with k ≲ β.In the first case, the upper bounds of errors are ≲ 4%, while in the second case they are ≲ 10%, as shown, respectively, by Figs. 6 and 8.In Case C), the approximate solutions trace also very well to the exact one, except the minimal points, as shown in Fig. 10.This might be caused by the fact that at these points the mode function µ k is almost zero, and very small non-zero values will cause significantly deviations.We are still working on this case, and hope to come back to this point in another occasion.
As mentioned in the Introduction, the potentials of the mode functions in both dressed metric and hybrid approaches can be well modeled by PT potentials.Therefore, the current analysis on the choice of the function q(y) and the minimization of the error control function shall shed great light on how to carry out similar analyses in order to obtain more accurate approximate solutions in these models.We have been working on it recently, and wish to report our results soon in another occasion.
In addition, the differential equations for the quasinormal modes of black holes usually also take the form of Eq.(2.1) with potentials that have no singularities 3 , but normally do have turning points [70,71].For example, the effective potential for the axial perturbations of the Schwarzschild black hole is given by where ω denotes the qusinormal mode.Clearly, for r ≥ 2m, this potential also has no poles, but in general f (r) ≡ V (r) − ω 2 have two turning points.From [70,71], it can be seen that the properties of this potential are shared by many other cases, including those from modified theories of gravity.Thus, one can equally apply the analysis presented here to the studies of quasi-normal modes of black holes.
we find that Eq. (B.5) Here 2 F 1 (a 1 , a 2 , a 3 , x) denotes the hypergeometric function, and a k and b k are two independent integration con-stants, and are uniquely determined by the initial conditions.

APPENDIX C. COMPUTING THE ERROR CONTROL FUNCTION
In this appendix, we collect some useful formulae for working out the error control function explicitly.Considering the particular form of the PT potential, it is easier to compute the error control function by using the new variable x = sech(y), thus where ϵ y denotes the sign of y.In terms of the new variable, To calculate the error control function explicitly, let us consider the cases g < 0 and g > 0 separately.

g < 0
In this case, the error control function is defined by Eq.(2.29), which can be written as where In this case, the error control function is defined by Eq.(2.30), which can be also written as Eq.(C.3), but now with

FIG. 1 .
FIG.1.The function ζ(y) vs y, which is assumed to be always an increasing function of y.

1 .
When y1,2 Are Real Let us first consider the case when y 1 and y 2 are real.Then, when y > y 2 , we have ζ(y) > ζ 0 [cf.Fig.1].Hence, from Eq. (2.18) we find y ) as x → 0 (or y → ∞).Note that ζ → ∞ as y → ∞, which can be seen clearly from Eq.(3.11).Thus, to minimize the2 In this case, the associated error control function is V ζ 1 ,ζ (T ) for any given ζ 1 , where ζ 1 ∈ (−∞, ∞) [44].In this paper, we choose ζ 1 = 0, so the integrations will be carried out in the interval ζ ∈ [0, ∞), corresponding to y ∈ [0, ∞).Due to the symmetry of error control function for very large values of y, we must the equation, one can easily obtain the solutions for the region y ∈ (−∞, 0] by simply replacing y by −y (or ζ by −ζ).