Homotopy Analysis Method for the Rayleigh Equation Governing the Radial Dynamics of a Multielectron Bubble

The homotopy analysis method is used to obtain analytical solutions of the Rayleigh equation for the radial oscillations of a multielectron bubble in liquid helium. The small order approximations for amplitude and frequency fit well with those computed numerically. The results confirm that the homotopy analysis method is a powerful and manageable tool for finding analytical solutions of strongly nonlinear dynamical systems.


Introduction
Nonlinear equations are widely used for modeling complex phenomena in various fields of sciences and engineering.Nonlinear problems are in general difficult to solve analytically.In recent years, the tool for finding analytical solutions of nonlinear equations known as the homotopy analysis method HAM has been developed by  This nonperturbation technique, independent of small/large physical parameters, allows to effectively control the convergence and accuracy of the series solution to the model under consideration.HAM has been applied successfully to many nonlinear problems such as free oscillations of self-excited systems 4 , heat radiation 5 , finding travelling-wave solutions of the Kawahara equation 6 , finding solitary wave solutions for the fifth-order KdV equation 7 , exponentially decaying boundary layers 8 , free vibrations of tapered beams 9 , and finding solutions for Duffing equations with cubic and quintic nonlinearities 10 .
HAM method is applied in this paper for finding a periodic solution of the Rayleigh equation, which describes the radial free oscillations of a multielectron bubble MEB in liquid Helium 11, 12 : where R is the bubble radius, ρ and σ are the mass density and surface tension of the liquid Helium, respectively, e is the electron charge, Z is the number of electrons in the bubble, ε is the dielectric constant of helium, ε 0 is thevacuum permittivity, R c is the equilibrium Coulomb radius of the MEB, and ξ is a dimensionless parameter.Another aim of this work is to obtain analytical solutions of the Rayleigh equation 1.1 with an additional external force driving term, a pressure step function.With this kind of forcing, Tempere et al. 13, 14 have shown that the angularly undeformed MEB undergoes cyclic collapses and reexpansions with high frequency in the order of MHz .During the last stage of collapse, the radial acceleration is sharp and large greater than 10 6 m/s 2 resulting in pulses of electromagnetic radiation or sonoluminescence.

Foundations of the HAM
The text sequence used in 9, 10 to show the basic ideas of the HAM is closely followed here.Consider a nonlinear differential equation expressed by the following: where N is a nonlinear differential operator, t denotes the independent variable; r t is an unknown function.
Applying the transformation τ ωt, 2.1 can be expressed by the following: Now a homotopy in general form can be constructed as follows: where q ∈ 0, 1 is an embedding parameter, φ is a function of τ and q, h is a nonzero auxiliary parameter, H τ is a nonzero auxiliary function, and L denotes an auxiliary linear operator.When the parameter q increases from 0 to 1, ω q varies from ω 0 to ω, and the solution φ τ, q varies from the initial approximation r 0 τ to the exact solution r τ .Put differently, φ τ, 0 r 0 τ is the solution of H φ, q, h, H τ | q 0 0, and φ τ, 1 r τ is the solution of H φ, q, h, H τ | q 1 0. Setting H φ, q, h, H τ | 0, the zero-order deformation equation is obtained as under the initial conditions: where k is a constant.The functions φ τ, q and ω q can be expanded as power series of q using Taylor's theorem as follows where r m τ and ω m are called the mth-order deformation derivatives.Differentiating the zero-order deformation equation with respect to q and then setting q 0, yields the first-order deformation equation which gives the first-order approximation of r τ as follows: subject to the initial conditions: Differentiating Equations 2.4 and 2.5 m times with respect to q, then setting q 0 and finally divided them by m!, yields the so called mth-order m > 1 deformation equation: subject to the initial conditions where D m r m−1 , ω m−1 , r m−1 and ω m−1 are defined as follows: 12 2.14 where α, θ, and A are dimensionless parameters given by the following:

Radial Free Oscillations
Introducing the new independent variable τ ωt, where ω is the dimensionless natural frequency, 3.2 and the initial conditions 3.3 can be rewritten as follows: In order to reduce the number of recursions and save computational efforts, this will be clarified later, an additional transformation is introduced R τ 3 r τ , 3.7 then, 3.5 becomes under the initial conditions r 0 ξ 3 , dr 0 dτ 0. 3.9 The free oscillation of the conservative system represented by 3.8 and 3.9 is a periodic motion, which can be expressed by the following base functions 1 : With the purpose of satisfying the initial conditions 3.9 , the initial guess of r τ for zeroorder deformation equation is chosen as follows: Under the Rule of Solution Expression denoted by 3.11 , the linear operator is selected as: r τ; q , 3.13 with the property where C 1 and C 2 are the integral constants.From 3.8 , the nonlinear operator is written as N r τ; q , ω q 1 3 ω q 2 r τ; q ∂ 2 r τ; q ∂τ 2 − 1 18 ω q 2 ∂r τ; q ∂τ For complying with the general form of the base functions 3.10 , the auxiliary function H τ is chosen as follows: H τ 1.

3.16
Substituting 3.15 into 2.11 and performing the homotopy-derivatives 3 , yields: 3.17 where a prime indicates first derivative with respect to τ.At this point, it is worth comparing 3.17 with the mth-order homotopy derivative corresponding with 3.5 , which reads as follows:

3.18
Clearly, the number of recursions for computing 3.18 is bigger than that using 3.17 .Thus, this last equation is used to save long calculations.
According to the property 3.14 of the auxiliary linear operator L, if the term cos τ exists in D m τ , the so-called secular term τ sin τ will appear in the final solution disobeying the rule solution expression 3.11 .To avoid the secular terms, the coefficient of cos τ in D m τ must be zero.
From 3.17 , the first-order approximation m 1 of HAM for D 1 τ is as follows:

3.20
Solving 2.7 subject to conditions 2.8 r 1 τ is obtained as:

3.21
Similarly, from the coefficient of cos τ in R 2 τ , ω 1 is obtained as follows: Solving 2.9 under conditions 2.10 for m 2, r 2 τ results in:

3.24
The higher-order approximations of r m τ and ω m−1 τ , for m > 2, are calculated similarly.

Driven Oscillations
The Rayleigh equation 1.1 with a driving pressure step p is written as follows: A nondimensional radius x and time t are defined as follows: x R R , t t t .

3.27
In terms of the dimensionless variables 3.27 , the equation of motion 3.25 becomes and the initial conditions 3.26 become where α, β, θ, and a are dimensionless parameters given by the following:

3.30
Under the transformation τ ωt, where ω is the dimensionless natural frequency, 3.28 and the initial conditions 3.29 can be rewritten as follows: x 0 a, dx 0 dτ 0.

3.32
Equation 3.31 can be integrated once 13 , yielding: Journal of Applied Mathematics 9 The form of 3.33 is simpler than that of 3.31 implying less computations.
In other words, the homotopy-derivative of the first term of 3.31 i.e., r 0 ω q−r ω r , which is at the order of magnitude n 8 is no longer present in 3.33 .Only remains the second term of 3.31 without the numeric factor 3/2 whose homotopy-derivative i.e., r 0 ω q−r ω r is at the order of magnitude n 8 .Even though the term α 2/3 β 2θ x appears in 3.33 , its homotopyderivative i.e., α 2/3 β 2θ x m−1 is only at the order of magnitude n 1 .
An additional simplification is found when 3.33 is substituted into 3.31 , resulting in Now, the first term of 3.31 is recovered instead of the first term of 3.33 .This is not advantageous because both terms involve the same number of recursions, but the second term in 3.33 and its corresponding recursions i.e., 2/3 From now on, the system represented by 3.34 and the initial conditions 3.32 , is adopted to be solved by HAM.
The oscillation of the conservative system 3.34 and 3.32 after the application of the step pressure is a periodic motion which can be expressed by the following base functions:

3.36
It should be noted that the system 3.34 and 3.32 oscillates between the initial radius x 0 a and the minimum radius x min , that is, when bubble wall velocity becomes zero dx/dτ 0; then from 3.33 : where Φ is given by the following:

10 Journal of Applied Mathematics
To satisfy the initial conditions 3.32 , the initial guess of x τ for zero-order deformation equation is chosen as follows: x 0 τ a − x min 2 cos τ a x min 2 .

3.39
Under the Rule of Solution Expression denoted by 3.36 , the linear operator for this case is identical to that chosen for free oscillations, that is, 3.13 and its corresponding property denoted by 3.14 .From 3.34 , the nonlinear operator is defined by the following: N x τ; q , ω q ω 2 q x 5 τ; q ∂ 2 x τ; q ∂τ 2 − 1 2 αx 3 τ; q 3 2 α β 3θ x τ; q − 4θ.

3.40
The solution must comply with the general form of the base functions 3.35 .Hence, the auxiliary function H τ must be chosen as follows:

3.42
As before, to avoid the secular terms, the coefficient of cos τ in D m τ must be zero.Thus, the frequency ω 0 is obtained by setting to zero the coefficient of cos τ in D 1 τ , yielding:

3.43
Then solving 2.7 under the conditions 2.8 , x 1 τ is obtained as follows:

3.45
The higher order approximations for ω m−1 and x m τ when m > 1 can be similarly calculated.Two sets of results are presented.The first one is for free oscillations Figures 1-7 and the second one is for driven oscillations Figures 8-11 .
The following characteristic scales were chosen to study the effect of varying both the initial radius of the MEB by changing ξ and the quantity of electrons in its surface by changing θ :

3.47
Figure 1 shows the curves of ω as function of h for different order approximations.The curves indicate that the valid regions of h i.e., the interval where ω versus h is a horizontal line are −3 < h < −2 and −3 < h < −1 for the 2nd and 5th-order approximations, respectively.Clearly, the convergence region of the series of the frequency ω increases with increasing the order approximation.The 2nd-order approximation for a proper value of h chosen in the valid interval −3 < h < −2 is compared with the numerical solution by Runge-Kutta scheme as illustrated in Figure 2.
Solutions in Figure 2 correspond to a MEB containing 10 4 electrons, as that numerically simulated by Salomaa and Williams 11 , accordingly its equilibrium Coulumb radius R c is about 1.064 microns.The initial radius of the MEB is 1.063 precisely the value of ξ times bigger than the equilibrium radius R c ; thus, after time zero, the bubble begins to oscillate with frequency around MHz. Clearly the 2nd-order HAM approximation is in excellent agreement with numerical solutions.This is quantified by computing the relative error for the Mth-order approximation as follows Relative Error r τ N − M m 0 r m τ r τ N × 100, 3.48 where r τ N is the solution by Runge-Kutta method.The maximum relative error, indicated with a filled triangle in Figure 3, between 2ndorder approximation and numerical solution is about 0.129%.When increasing the order of the approximation, the maximum relative error becomes smaller; accordingly, the maximum relative error between 5th-order HAM approximation and numerical solution is reduced to 0.063% indicated in Figure 3 with a hollow triangle , which is acceptably low.The effect of incrementing the initial radius of the MEB is depicted in Figure 4.As before, a bubble containing 10 4 electrons is considered but now the initial dimensionless radius is fixed to ξ 1.15.Seemingly, the 2nd-order approximation in Figure 4 a describes well the behavior of both the radius and the bubble wall velocity.Nevertheless, some deviation between HAM approximation and numerical solutions for the bubble wall acceleration is observed.This is overcome by increasing the order of approximation; consequently, the Journal of Applied Mathematics 1, θ 0.999939, h −2.3 .Dash-dotted line with symbols: relative error between 2nd-order HAM approximation and numerical results; solid line with symbols: relative error between 5th-order HAM approximation and numerical results.
5th-order HAM solution exhibits good agreement with the numerical solution as shown in Figure 4 b .
Figure 5 illustrates the first, second and fifth-order approximations of ω versus numerical results from Runge-Kutta method.As before, higher-order HAM solutions agree better with numerical results.The comparisons indicate that the 5th-order approximation agrees well with the numerical results even for a large initial dimensionless radius ξ 1.2.In fact, the maximum relative error between 5th-order approximation and numerical solution is about 0.066%; instead, the maximum error between the 2nd-order HAM solution and the numerical one is around 0.948%, as shown in Figure 6.It is interesting to note that the exponential growth tendency of the error curve for the 2nd-order approximation is dramatically reduced by increasing the order of HAM series.Clearly, the error curve for the 5th-order approximation shows slight upward trend with regard to that for the 2nd-order.Figure 7 shows the validity of the 5th-order HAM solutions for r τ when compared with numerical integrations.The effect of varying the quantity of electrons in the MEB while fixing its initial radius ξ is analyzed via changing the value of θ. Figure 7 a illustrates that smaller bubbles collapse more violently implying both smaller minimum radiuses and higher oscillation frequencies, for example, a small bubble containing 10 electrons attains at the first main collapse a minimum dimensionless radius ≈0.23 and develops an oscillation dimensionless frequency ω ≈ 0.33, whereas a big bubble containing 10 8 electrons attains at the first main collapse a minimum dimensionless radius ≈0.72 and develops a dimensionless frequency ω ≈ 0.27.Clearly, the analytical approximations for different values of θ closely approach the numerical results.Accordingly, in quantitative terms, the maximum relative errors 0.418, 0.154, 0.152, and 3.068% for the respective values of θ 0.617, 0.833, 0.999, and 1.05 are acceptably low; even for the biggest bubble which contains 10 8 electrons.
In Figure 7 b , the robustness of the 5th-order HAM solution is analyzed for high values of the initial radius ξ while fixing the quantity of electrons in the MEB.For ξ ≈ 1.07, the bubble implodes reaching both high acceleration ≥1 × 10 6 m/s 2 and a minimum radius more and more close to zero.This strong collapse conditions might involve sonoluminescence phenomenon as suggested by Tempere et al. 13 .Even for ξ as high as 1.4 the analytical solution agrees to some extent with the numerical results, the maximum relative error in this case is around 5%.It is worthy of note that the parameter h was varied to ensure the convergence of the analytical solution, the freedom of changing conveniently h is precisely one of the virtues of the HAM 1, 6, 9 .
In order to study the effect of varying the driving pressure p in 3.25 , the following scales where chosen:  Figure 8 depicts comparisons between the 4th-order HAM solutions and numerical integrations of: x τ , x τ , and x τ for two values of the pressure step.Excellent agreement between analytical and numerical solutions is clearly seen for p 300 Pa β 0.443495 ; accordingly, the maximum relative error between numerical and analytical results for x τ is 0.304%, indicated by a hollow triangle in Figure 9 a .Nevertheless, when the pressure step is increased to p 1000 Pa β 1.47832 , some disagreement between HAM approximations and numerical solutions of: x τ , x τ , and x τ is observed.In this case the maximum relative error for x τ is 4.664%, which might be acceptable if rough estimations are required.Even higher disagreement is clearly observed between HAM approximation and numerical integration for x τ .This discrepancy might be reduced by incrementing the order of the HAM approximation; but this would imply larger and larger analytical solutions, and therefore more and more difficult to handle with.For the current example, an explicit expression for the 4th-order x τ solution would require several pages to write upon.As pointed out by Liao 3 , in present time of computer with huge data storage capacities and high-speed CPU, an analytical expression with many terms might be accepted by most researchers.However, an important factor to be taken into account is the CPU time for evaluating an analytical solution with many terms of course this time tends to increase together with the number of terms .The CPU time consumed in evaluating the 4th-order approximation for x τ is 12.933 s computed by Timing command of Mathematica 5.1 in a Windows PC with Intel Core i7 processor at 3.4 GHz , which is around 275 times longer than time for computing a numerical solution of 3.27 and 3.28 for x τ using the fourth order Runge-Kutta algorithm this time is computed by Timing and NDSolve commands of Mathematica 5.1 in a Windows PC with Intel Core i7 processor at 3.4 GHz .Of course, this result does not invalidate the 4th-order analytical approximations but it makes them unsuitable for certain situations.For instance, an analysis of shape instabilities of bubbles using the 4th-order approximations for x τ , x τ and x τ would imply more computing time and less accuracy than those obtained by a numerical scheme.Figure 10 illustrates the variation of dimensionless bubble frequency ω versus β.Clearly, the analytical approximation for ω solid line is very close to the numerical calculations symbols .Even though the relative error has a tendency to grow exponentially as shown in Figure 11, its values are acceptably low for large magnitude of the step forcing, for example, for β 2 p 1360 Pa the relative error is just about 1.15%.

Conclusions
In this study the HAM has been used to obtain analytical solutions of the Rayleigh equation for the radial oscillations of a MEB in liquid helium.The small-order HAM approximations for freely oscillating bubbles agree very well with numerical solutions even for bubbles with initial radial amplitudes as high as 1.6 times the equilibrium Coulomb radius.The analytical solutions for radius, velocity and acceleration of the freely oscillating bubble wall are accurate enough to accomplish surface stability studies both parametric and Rayleigh-Taylor instabilities could be computed with the possibility of both saving calculations and giving a bigger understanding of bubble shape instabilities when compared to solutions from a numerical scheme.
In the case of forced oscillations, the fourth-order HAM solutions for displacement and velocity of the bubble wall agree well with those computed numerically.Nevertheless, when the magnitude of pressure step is large enough β 1.47832, p 1000 Pa , noticeable differences are observed between analytical and numerical curves for acceleration.It was shown that higher order approximations can be adapted to increase the convergence of the solution but at the expense of a huge number of terms in it, which might be impractical.Even this inconvenience, current analytical approximations show the great potential of the HAM for complex problems with strong nonlinearities.HAM offers the possibility of controlling in a convenient way the convergence of approximation series; this fundamental characteristic is what makes the HAM more powerful than other nonperturbation techniques such as Adomian decomposition method 15 and the homotopy perturbation method HPM 16 .It was confirmed here that convergence of solutions is ensure by choosing h parameter in an appropriate way.

Figure 9 :Figure 10 :
Figure 9: Relative error between 4th-order HAM solution and the corresponding numerical integration for displacement x τ ; a and b parts of this figure match with those in Figure 8.

Figure 11 :
Figure 11: Relative error between 4th-order HAM approximation and numerical solutions for the frequency ω.