Analysis and Computation of Solutions for a Class of Nonlinear SBVPs Arising in Epitaxial Growth

: In this work, the existence and nonexistence of stationary radial solutions to the elliptic partial differential equation arising in the molecular beam epitaxy are studied. Since we are interested in radial solutions, we focus on the fourth-order singular ordinary differential equation. It is non-self adjoint, it does not have exact solutions, and it admits multiple solutions. Here, λ ∈ R measures the intensity of the ﬂux and G is stationary ﬂux. The solution depends on the size of the parameter λ . We use a monotone iterative technique and integral equations along with upper and lower solutions to prove that solutions exist. We establish the qualitative properties of the solutions and provide bounds for the values of the parameter λ , which help us to separate existence from nonexistence. These results complement some existing results in the literature. To verify the analytical results, we also propose a new computational iterative technique and use it to verify the bounds on λ and the dependence of solutions for these computed bounds on λ .


Introduction
Epitaxy means the growth of a single thin film on top of a crystalline substrate. It is crucial for semiconductor thin film technology, hard and soft coatings, protective coatings, optical coatings, etc. The epitaxial growth technique is used to produce the growth of semiconductor films and multilayer structures under high vacuum conditions [1]. The major advantages of epitaxial growth are reducing the growth time, better structural and superior electrical properties, eliminating waste caused during growth, wafering cost, cutting, polishing, etc. Several types of epitaxial growth techniques, such as hybrid vapor phase epitaxy [2], chemical beam epitaxy [3], and molecular beam epitaxy (MBE), have been used for the growth of compound semiconductors and other materials. In this work, we strictly focus on MBE, and we restrict our attention to the differential equation model, which was proposed by Escudero et al. [4][5][6][7]. The mathematical description of epitaxial growth is carried out by means of a function σ defined as σ : Ω ⊂ R 2 × R + → R, which describes the height of the growing interface in the spatial point x ∈ Ω ⊂ R 2 at time t ∈ R + . The authors in [4][5][6][7] show that the function σ obeys the fourth-order partial differential equation where η(x, t) models the incoming mass entering the system through epitaxial deposition, λ measures the intensity of this flux, and the determinant of Hessian matrix is (2) The stationary counterpart of the partial differential Equation (1) subject to the homogeneous Dirichlet boundary condition (4) and homogeneous Navier boundary condition (5) is defined as (see [6]) where η(x, t) ≡ G(x) is a stationary flux, and n is the unit out drawn normal to ∂Ω.
Using the transformation r = |x| and σ(x) = φ(|x|), as a result of symmetry, the above set of equations are transformed into the following set of equations: where = d dr .
In this paper, we also impose the following boundary conditions which complements the work in [6]: For simplicity, we take G(r) = 1, which physically means that the new material is being deposited uniformly on the unit disc. Now, using lim r→0 rφ (r) = 0, w = rφ and integrating parts from Equation (6), we have Using the transformation t = r 2 2 and u(t) = w(r), it is possible to reduce Equation (10) into the following equation: Corresponding to (11), we define the following three boundary value problems: Problem 1: Problem 2: Problem 3: The BVPs (12), (13) and (14) can be equivalently described as the following integral equations (IE): • IE corresponding to Problem 1: • IE corresponding to Problem 2: • IE corresponding to Problem 3: We assume that u ∈ C 2 loc 0, 1 2 ; R , where C 2 loc 0, 1 2 ; R is defined as In [6], Escudero et al. proved the existence and nonexistence of solutions of Problems 1 and 3 using upper and lower solution techniques. Corresponding to Problems 1 and 3, they have also provided the rigorous bounds of the values of the parameter λ, which helps us to separate existence from nonexistence. In [8], Verma et al. provide numerical illustrations via VIM to verify the results of Escudero et al. [6]. To verify their numerical results, they provided other iterative schemes based on homotopy [9] and the Adomian decomposition method [10].
Equation (13) has not been investigated theoretically in the existing literature to the best of our knowledge. Moreover, many investigations are still pending relating to BVPs (12), (13), and (14). Here, we focus on both theoretical and numerical work. We derive the sign of the solution and prove its existence in continuous space. We also compute the bounds of the parameter λ. The results of this paper complements existing theoretical results. We also provide an iterative scheme based on Green's function to compute the bounds and solutions to demonstrate the existence and nonexistence, which is dependent on λ.
To prove the existence of the solutions, we use the monotone iterative technique [11][12][13][14][15][16][17]. Recently, many researchers applied this technique on the initial value problem (IVP) for the nonlinear noninstantaneous impulsive differential equation (NIDE) [18], p-Laplacian boundary value problems with the right-handed Riemann-Liouville fractional derivative [19], etc. to prove the existence of the solution. Here, we also present numerical results to verify the theoretical results. To develop the iterative scheme based on Green's function, we consider Equations (12)- (14). Recently, many authors have used numerical approximate methods like the VIM [8], the Adomian decomposition method (ADM), the homotopy perturbation method (HPM) etc. to find approximate solutions for different models involving differential equations [20,21], integral equations [22][23][24], fractional differential equations [25,26], the Stefan problem [27][28][29][30], system of integral equations [31], etc. Thereafter, Waleed Al Hayani [32] and Singh et al. [33] applied ADM with Green's function to compute the approximate solution. Recently, Noeiaghdam et al. [34] proposed a technique based on ADM for solving Volterra integral equation with discontinuous kernels using the CESTAC method. To find out more about this method, please see [35,36]. They focused on the BVPs which have a unique solution. The major advantage of our proposed technique is its ability to capture multiple solutions together with a desired accuracy.
The remainder of the paper is focused on both theoretical and numerical results. We prove some of the basic properties of the BVPs in Section 2. The monotone iterative technique is presented in Section 3 to prove the existence of a solution. A wide range of λ for Equation (6), corresponding to different types of boundary conditions, is shown in Section 4. In section 5, we apply our proposed technique to the integral equations and show a wide range of numerical results. Finally, in Section 6, we draw our main conclusions.

Lemma 7.
Let u(t) ∈ C 2 loc (0, 1 2 ], R be the solution of Problem 2, then u(t) can be written in the following form: and also satisfies Proof. Using the boundary condition and properties of Green's function, we have From Equation (29) and Problem 2, we can easily derive Equation (27). Now, using the result of Lemma 1, we have Therefore, from Equation (30) and using a similar analysis as that in Lemma 6, we can prove the result (28).
, R be the solution of Problem 1, then u(t) can be written in the following form: and satisfies Proof. The Green's function of Problem 1 is given by Again, from Equation (33) and Problem 3, we derive integral Equation (31). Furthermore, using the result of Lemma 1, we have Again, using a similar analysis as that in Lemma 7, we get the inequality (32).

Existence of Solutions
In this section, we apply the monotone iterative technique coupled with lower and upper solutions to prove the existence of at least one solution for Problems 1-3. For this purpose, we need to prove some lemmas, which help us to prove the main results of this paper.

Construction of Green's Function
To investigate the Problems 1-3, we consider the corresponding nonlinear singular boundary value problems, which are given by Problem 1(a): Problem 2(a): Problem 3(a): where Throughout the paper, we assume the following conditions: Lemma 9. Let k satisfy H 0 and u(t) ∈ C 2 loc (0, 1 2 ], R be the solution of Problem 1(a), then where Green's function G(s, t) is given by and G(s, t) ≤ 0 for all 0 ≤ s ≤ 1 2 and 0 ≤ t ≤ 1 2 .
Proof. Using the boundary condition of Problem 1(a) and the properties of Green's function, we can easily prove Equation (39). Furthermore, we have G(s, t) ≤ 0 for all 0 ≤ s ≤ 1 2 and 0 ≤ t ≤ 1 2 .
Proof. In a similar manner to that in Lemma 9, we can easily obtain Equation (41) and prove G(s, t) ≤ 0 for all 0 ≤ s ≤ 1 2 and 0 ≤ t ≤ 1 2 .
Lemma 11. Let k satisfy H 1 and u(t) ∈ C 2 loc (0, 1 2 ], R be the solution of Problem 3(a), then where Green's function G(s, t) is given by where Proof. Again, using a similar analysis, we can easily derive the Green's function. Now, Hence, from (43), we have G(s, t) ≤ 0 for all 0 ≤ s ≤ 1 2 and 0 ≤ t ≤ 1 2 .
Proof. The proof is similar to that shown in Lemma 9.
Proof. The proof is similar to that shown in Lemma 10.

Lemma 14.
Let k satisfy H 4 and u(t) ∈ C 2 loc (0, 1 2 ], R be the solution of Problem 3(a), then where Green's function G(s, t) is given by where Proof. The proof is similar to that shown in Lemma 11. Proposition 1. Let k satisfy H 0 or H 2 (respectively, H 0 or H 3 and H 1 or H 4 ) and h(t) ∈ C 2 loc (0, 1 2 ], R is such that h(t) ≥ 0, then the solution of Problem 1(a) (respectively, Problem 2(a) and Problem 3(a)) is nonpositive.

Monotone Iterative Technique
Here, we define lower and upper solutions corresponding to Problems 1-3.

Theorem 1.
Assume H 0 (respectively, H 0 and H 1 ) is true, there exist α 0 , and β 0 ∈ C 2 loc (0, 1 2 ], R are upper and lower solutions of Problem 1 (respectively, Problem 2 and Problem 3), which satisfy the properties P 1 and P 2 such that β 0 ≤ α 0 = 0, then the Problem 1 (respectively, Problem 2 and Problem 3) has at least one solution in the region D 0 and the sequences {α n }, {β n } defined by (52)-(55) converge to solutions u, v uniformly and monotonically, respectively, such that Proof. We divide the proof into three parts. In the first part, we prove that β n is a lower solution of problem 1, β n ≤ β n+1 , and β n+1 ≤ α 0 ∀n ∈ N ∪ {0}.
In the second part of the proof, we have to show that α n is an upper solution of Problem 1 and α n+1 ≤ α n ∀n ∈ N.
Now, from (52) and (53), we have Therefore, by using (50), we have Hence, by Proposition 1, we have α 1 ≤ α 0 . Therefore, our assumptions are true for n = 0. Let our assumptions be true up to n = m. Then, we find that α n is an upper solution of Problem 1 and α n+1 ≤ α n for n = 1, 2, . . . , m.
In the last part of the proof, we want to show β n ≤ α n for all n ∈ N. Again, from (71) and (86), we have Since β n ≤ α n ≤ 0, we have and Hence, by Proposition 1, β n+1 ≤ α n+1 . Finally, we have Let t n ∈ 0, 1 2 for n ∈ N such that t n+1 < t n for n ∈ N, lim n→∞ t n = 0.

Estimations of λ
The objective of this section is to derive some qualitative bounds of the parameter λ, from which we can conclude about the nonexistence of solutions. Equation (11) can be written in the following form: Put v(t) = − u(t) t and integrating from 0 to t, Equation (100) becomes In view of the transformation, the boundary condition at r = 1 becomes BC of Problem 1: v BC of Problem 2: v Escudero et al. in [6] prove the following two lemmas: The set of numbers λ ≥ 0, for which there exists a solution u(t) ∈ C 2 loc (0, 1 2 ], R of Equation (11) satisfying lim t→0 u(t) √ t = 0 and u(t) ≤ 0, is nonempty and bounded from above.
We present the following results which complement the results proved by Escudero et al. [6].

Theorem 3.
Let λ 0 ∈ R + . If 0 ≤ λ < λ 0 , then Equation (10) corresponding to different types of boundary conditions are solvable. Moreover, there is no solution to these problems if λ > λ 0 . Furthermore, every solution w(r) of a governing equation corresponding to these three types of boundary condition satisfy w(r) ≤ 0, r ∈ (0, 1] and lim Proof. The proof of this can be deduced from Lemma 15, Lemma 16, Lemma 1, Lemma 2, Lemma 3, and Lemma 5.

Numerical Results and Discussion
Here, we present the numerical data to validate our derived theoretical results. In Section 5.1, we derive the numerical estimation of the bounds computed by ADM. In Section 5.3, we numerically show the existence of at least one solution.

ADM
To find the approximate solutions, we develop the iterative numerical schemes with the help of the Fredholm integral Equations (15), (16), and (17), respectively. Now, we decompose the solution u(t) of the form u(t) = ∑ ∞ i=0 u i (t), and approximate the nonlinear term in terms of Adomian's polynomials [38], which is given by where Therefore, from integral Equation (15), we define We compute the arbitary constant c using the Mathematica 10.9 program. For better understanding, we present the algorithm of our proposed technique corresponding to Equation (15) below.
Residue Error: Here, we define the residue error [39] corresponding to Equation (15) for error analysis, which is given by where λ is the parameter. Therefore, the maximum absolute residue error can be defined as
Step 2. Identify the constant term, and approximate the nonlinear term by Equation (143).
Step 4. Approximate the term Step 5. Compute the values of the constant and the approximate solutions u(t) = n+1 ∑ i=0 u i .
Step 6. Determine the residue error R(t) and set the stopping criteria L ∞ < , where is the tolerance.
Again, we apply the algorithm Section 5.2 to Equations (16) and (17), and we define the following iterative schemes: and Scheme of Problem 3 = A n 2s 2 ds, . . ., and c = − Approximate solutions for Equations (16) and (17) can be written as provided the series is convergent for n → ∞. Recently, the convergence of ADM was established by Verma et al. in [9]. Now, by using the transformation t = r 2 2 , u(t) = w(r), w(r) = rφ (r), and φ(1) = 0, we get the solutions of Equation (6). We arrive at two cases: Case (a): λ ≥ 0. For λ = 0, we get one trivial and one nontrivial solution. For 0 < λ ≤ λ critical , we always find two nontrivial solutions. We may refer to them as upper and lower solutions, respectively. Corresponding to Equations (9), (8), and (7), we find the critical values of λ, i.e., λ critical , are 31.94, 11.34, and 168.76, respectively. For λ > λ critical , we do not find any numerical solutions, as the value of c become imaginary. In Section 5.2.1, we tabulate residual errors of the approximate solutions corresponding to some λ. Case (b): λ < 0.
In this case, we always have two nontrivial numerical solutions corresponding to three types of boundary conditions. One solution is negative (namely, the negative solution) and the other solution is positive (namely, the positive solution). We do not find any negative critical λ. Please refer to Section 5.2.1 as regards residue errors.

Tables
Here, we have placed some numerical data of approximate solutions of φ(r) corresponding to different types of boundary conditions below. If we are increasing the value of λ, we see that the residue error of the lower solution is increasing and the residue error of the upper solution is decreasing (see : Table 1). Similarly, if we are decreasing the value of negative λ, we see that the residue error of both positive and negative solutions are decreasing (see: Table 2). The same can be seen in Tables 3-6.

Monotone Iterative Method
Here, we compute the monotone iterations using Equations (52)-(55) corresponding to three types of boundary condition.
Corresponding to problem (12): By using Lemma 19, we chose the lower and upper iterations Therefore, it is easy to show that α 0 and β 0 both satisfy the inequalities (50), (51), (56), and (57), such that β 0 ≤ α 0 . We consider that k = −1. Hence, by using Theorem 1, we have a monotonically and uniformly convergent sequence {β n } and {α n }, which are converging to the solution v and u of the problem (12). We denote u ADM as the approximation of the solution of (12) computed by ADM. By using the transformation t = r 2 2 , w(r) = u(t), and w(r) = rφ (r), we have the fourth-order iterations corresponding to the monotone iteration α i and β i . φ α i and φ β i are the fourth-order solutions corresponding to the monotone iteration α i and β i for i = 0, 1, . . ., respectively.
In Figure 1a, we plotted β 0 , β 1 , u ADM , α 1 , and α 0 corresponding to problem (12) for λ = 2. We have seen that the lower sequence {β n } and upper sequence {α n } always satisfies the inequality β n ≤ α n . In Figure 1b, we have placed the monotone iterations of fourth-order SBVP corresponding to problem (12) for λ = 2. Here, we also observed the existence of at least one solution for a fourth-order SBVP corresponding to the Dirichlet boundary condition. Corresponding to problem (13): From Lemma 18, we chose the initial monotone iterations as follows: The same remarks follow as discussed above.
In Figure 2a,b, we present the numerical results for k = −1 and λ = 2. We noticed that the approximate solution u ADM and φ ADM always lies between the lower sequence {β n } and upper sequence {α n }. Corresponding to problem (14): Here, we consider the initial monotone iterations We also included the same remarks as those stated above. Monotone lower and upper iterations corresponding to the second-order and the fourth-order differential equation are plotted in Figure 3.

Conclusions
In this work, we derived some qualitative properties of the singular boundary value problems that arise in the theory of epitaxial growth. Moreover, we proved the existence of a solution and discovered a range of parameter k, for which the nonlinear problem has multiple solutions in the region D 0 . We established the bounds of the parameter λ, from which we confirmed the nonexistence of solutions. Furthermore, the boundary value problems have multiple solutions, therefore it is challenging for researchers to obtain a suitable scheme to capture both solutions with the desired accuracy. However, we successfully developed iterative schemes and captured both solutions with a high accuracy. From Tables 1-4, we can see that the approximate solutions computed by our proposed method converge to the exact solutions very quickly. Corresponding to the boundary conditions (7), we notice that the positive approximate solution converges to the exact positive solution very slowly (See Table 6). We verified that our numerical results matched well with our theoretical results as well as the existing numerical results [9]. We conclude that our proposed technique is relatively powerful and efficient. Furthermore, this technique is an effective tool to solve BVPs, which have multiple solutions.