An effective method in investigating structures of polytropic protoplanets formed via gravitational instability.

In this article, we have reinvestigated the initial distribution of thermodynamic variables inside the protoplanets formed via gravitational instability having mass range 0.3 - 10 M J ( 1 M J = 1.8986 × 10 30  gm ) by an embedded RKACeM(4,4) method assuming that the polytropic gas law holds in the protoplanets. The findings attained by our numerical experiments are recognized to be consistent with the results acquired through other notable investigations in this regard. Furthermore, the model is easily computable. The used method is found to be efficient in investing the structures of polytropic protoplanets in their initial stages in terms of accuracy, stability, computational cost, and solving endpoint constraints.

Introduction

Certainly, numerical computing plays an essential and important role in solving real-time problems arising in physics, mathematics, engineering, and other branches of science to provide an optimal and efficient solution. In such computing, three stages are of interest, namely the formation of a suitable numerical method, the application of the scheme to achieve an efficient solution, and the confirmation of the acquired findings [1]. But before selecting and/or constructing novel methods, it is a prerequisite to figuring out different aspects, viz. type of the equation of interest, the availability of equipment, programming, information about the speed of execution, the validity of the attained results with reliability and accuracy [1], [2]. It is worth noting that RKACeM(4,4) is an embedded technique, dubbed with the 4th order Runge-Kutta (RK) arithmetic mean (RKAM(4,4)) method and 4th order RK based centroidal mean (RKCeM(4,4)) method, which provides the facility of selecting step size to control local truncation error. The method of interest can be used to address various types of problems involving ordinary differential equations (ODEs) with initial or boundary values. The method also provides an advantage to use a larger step size in integrating an initial value problem (IVP) or boundary value problem (BVP) in ODEs, for the RKCeM(4,4) method has a slightly larger stability region than that of the RKAM(4,4) method as will see later. Further, if the local truncation error of an RK-based method approaches zero, it is said to be consistent. It is well known that in RKAM(4,4) method, there is no facility for estimating local truncation errors. This leads to the formulation of a variety of embedded techniques, namely the RKARMS(4,4) method, RKACeM(4,4) method, RKAHeM(4,4) method, RKF(4,5) method, etc. But not all of the methods are suitable for solving all the types of the IVPs/BVPs arising in mathematical physics [1], [2].

It is worth mentioning that the interest of researchers in planetary science has been rekindled with the detection of a planet (51 Pegasi B) outside our solar system. Thereby, a significant volume of works has been conducted inside our solar system and elsewhere by using the physical conditions regulating the interior of such planets [3] and the researches are still being carried out towards the same. However, the details about the evolution process of the planets are still debated [4]. But it is believed widely that the solar system planets or extrasolar ones are formed from materials having high orbital angular momentum left over from the star's birth [5]. The two competing paradigms, found in the literature, to explain the evolution process of gas giants are (1) core accretion (CA) [6], [7] and (2) gravitational (disk) instability (GI) [8], [9]. In accordance with the first paradigm, giant planets are formed by the accumulation of solid bodies followed by the accretion of a gaseous envelope [5]. This paradigm has, so far, been accepted as the standard one in explaining the process of evolution of the solar planets as well as that of extrasolar ones. But there are some problems with the CA scenario. One major problem is that the formation of gas giants via the CA takes a very long time. It is also believed that the gas neighbouring a young star disappears before the solid core formation [10]. As in [11], the CA paradigm's main difficulty lies in its very beginning of the growth. It is, so far, unclear how metre-sized rocks stick together while colliding at high speeds, subject to high radial drifts into the parent star. Further, the CA scenario cannot efficiently describe the process of formation of gas giants at radial distances out to ∼100 AU from the parent star [12]. However, many authors (see e.g., Ormel and Klahr [13]) used the pebble accretion mechanism to solve the limitations of the CA scenario. But the subsequent pebble growth into planetesimals is not clearly understood [14]. As in [15], the streaming instability mechanism can make available an acceptable solution in this regard. Upon satisfaction of the threshold of the dust-to-gas ratio, the pebble clumps can collapse straight into planetesimals. It is seen that the growth of the core via the accretion of pebbles is much quicker; however, the mechanism is not so effective [16], [17] and may necessitate more pebbles than the observations indicate [18]. Furthermore, the formation of terrestrial planets is relatively fast [19]. With the problems in the CA mechanism, the GI paradigm, an alternative to the CA mechanism of planetary formation, has been reformulated with fragmentation from massive protoplanetary discs (PPDs) [11], [18], [20], [21], [22], [23]. But, as in [20], this paradigm has also received a lot of criticism. Some noteworthy problems of the scenario are (i) PPDs are unable to form planetary embryos at the present location of Jupiter [24], conflicting with the earlier results by Boss [9]; (ii) sedimentation of dust is so a slow process that it is failed to yield observed solid cores in giant embryos of mass [25]; (iii) OB stars (hot, massive stars of spectral types O or early-type B) are too rare to explain the abundance of Jupiter like planets observed in exosolar planetary systems [20] and (iv) it is believed that the formation of Earth-sized terrestrial planets is rarely possible by the GI [4], [20]. However, many authors have tried to solve the problems found in the GI scenario [11], [18], [20], [21], [22], [26], [27], [28], [29]. The tidal downsizing mechanism, a revised form of the GI paradigm with planet migration inward and tidal disruptions of GI fragments in the inner regions of the disc [27], could account for all observational facts relevant to the process of planetary evolution [22]. With the theoretical work as well as simulations, Humphries and Nayakshin [18], [19] exhibited that giant planets assembled by the GI can be very metal-rich as required by the solar and exoplanetary data and that these planets can migrate inward and explain the closer-in data as well [29]. Recently, Atacama Large Millimeter/submillimeter Array (ALMA) has been used widely to investigate the PPDs, which makes available information on the planetary systems orbiting stars other than the Sun during all stages of evolution [23], [30]. As in [18], the age of the ALMA planets is not a challenge for the GI mechanism. Because in this mechanism, planets form within the first ∼0.1 Myr [12].

Thus, the GI mechanism with some amendments can be a promising route to the rapid formation of gas giants in the exterior of our solar system, including our own. But one of the problems lies within the estimation of the initial configuration [4]. Through literature survey, it is seen that diverse numerical models present different structures for initially formed protoplanets [4], [31], [32], and no author has, so far, shown that these protoplanets exit, in reality, with their definite structures [2]. This leads researchers to conduct more research on the estimation of initial structures of the protoplanets formed by the GI envisioned by Boss [9], and planetary formation through this mechanism is still a subject for ongoing research. However, Boss [33], in his investigation, presumes the protoplanets in their initial stage to be in radiative equilibrium, whereas in [25], [31], [34], the gas blob of the protoplanets was assumed to be fully convective with a thin outer radiative zone, which is consistent with the assumption made in [35]. In their investigations, Paul et al. [36], [37] considered such a protoplanet to be convective fully, which is consistent with that found in [38]. On the other hand, Paul and Bhattacharjee [39] and Paul et al. [4], [40], [41] performed their numerical experiments considering the conductive-radiative heat transport. Paul et al. [42], [43] also used a polytropic equation of state, where they concluded that the polytropic protoplanets (PPs) with polytropic index as well as , the distributions of thermodynamic and other variables are nearer to reality. Here, the initial structure of the planet formation approach is on the basis of the renowned polytrope conjecture, which is used to characterize gaseous planets, main-sequence and fully convective stars, and even compact objects like neutron stars and white dwarfs [44]. A polytrope is a simple structural assumption between pressure and density, assumed to be valid throughout a gas sphere and can provide significant insight into the structure and evolution of stars. Our intention is that whether the law can provide significant insight into the structure of a protoplanet. However, this law leads to the formulation of the Lane–Emden (LE) equation, which is a dimensionless form of the Poisson equation for the gravitational potential of a Newtonian self-gravitating, spherically symmetric, polytropic fluid [45]. Such an equation, in most cases, may be solved only numerically. It is known that the LE equation has exact solutions for , 1, and 5; for the rest, one may rely on numerical or power series methods. The key interest when solving the LE equation using power series is how to make the series converge to the outer surface of a gas sphere. However, various endeavours have been made to enumerate what the index of polytrope signifies in the restricted context of the ideal gas law [46]. There exist other endeavours to import physical meaning to the index of polytrope [47], [48], but limited in scopes in their results. This work will, to the best of our knowledge, be the first to reveal the configuration of polytropic protoplanets with the physical meaning of the index of polytrope in the usual case for the protoplanets formed via GI with an efficient approach.

Therefore, this communication aims to investigate initial configurations of PPs, formed via GI, with mass range 0.3-10 by an embedded RKACeM(4,4) method for having optimal and efficient results.

Mathematical foundation

Protoplanetary structure

We assume isolated spherically symmetric gaseous objects of solar composition, formed via GI, having a mass limit of 0.3-10. The reason behind the choice is that the mass limit covers most of the detected extrasolar giant planets [31], [37]. As in [4], we presume that during the pre-collapse stage, the object contract quasi-statically, and the energy source is only the gravitational contraction. Each of the objects is assumed to behave as a polytrope meaning that the local pressure and density are related through a power law, known as the polytropic gas equation of state [49]: In Eq. (1), P designates the pressure inside the PP at a distance interior to a radius r; ρ is its density at the distance r from the centre; is a constant of proportionality, which is related to the total mass of the configuration; the ratio of specific heat at constant pressure to that at constant volume, , n is the polytropic index (not necessarily an integer). It is of interest to note that for a gas sphere with adiabatic convective heat transport, , which leads to , from which an estimation of the index n in the case of PP can be made. It is worthy to mention at this juncture that the polytropic index n appropriately signifies the behaviour of polytropic gas, where for initial PPs, n should significantly be small () [43]. Because at the initial stages, the PPs remain less centrally condensed [42], [43], as is to be expected.

The polytropic gas equation of state given by Eq. (1) by combining with the equations of hydrostatic support and mass conservation gives rise to the LE [45], named after astrophysicists Jonathan Homer Lane and Robert Emden. It describes the density profile of a gaseous self-gravitating object [45]. The equation is of importance in astrophysics, for the values of the polytropic index n between 0 and 5, the equation approximates to a reasonable accuracy the structures of a variety of realistic stellar models. Following Chandrasekhar [50], the LE equation governing the equilibrium structure of a self-gravitating isothermal gas sphere can be given by Eq. (2) is a second-order non-linear ODE in which the independent variable ξ is a dimensionless radius, the dimensionless independent variable θ is the polytropic temperature. The polytropic index n of the gas sphere is a free parameter ( and for cases of practical interest in problems of planetary structure, [43]. The centre () and surface (, is the value of ξ where the first zero of the dependent variable θ occurs) are two singularities of this ODE. Physical conditions at the centre and surface require Eq. (2) to satisfy the standard boundary conditions (BCs) , and the slope at . Physical conditions in the interior of a polytropic model require θ and remain finite for every value of ξ lies between 0 and . The polytropic temperature, θ, mentioned above is linked to T, ρ, and P through In Eq. (3), , , and specify the central pressure, density, and temperature, respectively, of a PP. The central values for , 1.5 are, respectively, given by , , and , where K can be determined by the mass radius relation, as we will see later. But for , the pressure-density relation cannot be used explicitly, as diverges when . However, for , , and can be estimated by means of the above corresponding mentioned equations, whereas can be attained by [43]. For , the calculation for central values is presented in subsection 2.2.

In the equations presented above, G symbolizes Newton's universal gravitational constant, M represents the mass inside radius R, i.e., the entire mass of a PP; , , and represent constants partaking different values for different n. For such n, the values of and are obtainable in any standard book [51]. The value of for a specific n can be attained via , where designates the mean molecular weight of the gas of the protoplanets containing hydrogen, helium and heavy elements, H ( gm) represents the mass of a hydrogen atom, and k () represents the Boltzmann constant.

The equation specifying the mass distribution in a PP is [43] where denotes the mass of the PP interior to r.

Mass-radius relationship for polytropes

Using Eqs. (2) and (4), we see where . Now, Eq. (5) on eliminating with yields the mass-radius relation presented below: From Eq. (6), we see that if we specify the mass and radius of a PP, and supply the assumed polytropic index n, except and 1, K can be determined, which is essential because K summarizes the LE equation's thermodynamical property. Our calculated K values for and are presented in Table 1.

Table 1

Polytropic constant K for n = 0.5 and n = 1.5.

Mass	R	K
		n = 0.5	n = 1.5
0.3 MJ	3.5 × 1012	1.3971 × 1026	8.2153 × 1014
1 MJ	5.3 × 1012	3.3371 × 1026	1.8583 × 1015
3 MJ	7.8 × 1012	7.6797 × 1026	3.9444 × 1015
5 MJ	8.4 × 1012	6.6745 × 1026	5.0364 × 1015
7 MJ	9.1 × 1012	7.1138 × 1026	6.1037 × 1015
10 MJ	11.0 × 1012	1.2852 × 1027	8.3095 × 1015

It is to be noted here that for , the mass and radius of a protoplanet are independent of each other (see Eq. (6)). In this case, K can be generated for a given R from Eq. (6) as . Then with that constant value of K, the central density of each of the PP having assumed masses can be obtained from , which is obtained from Eq. (5) for , and thereby with the relation , the central value of pressure, , for each PP can be obtained. The central temperature for the constant μ can be obtained from with as K and are now known.

Numerical method

Nondimensionalization

For numerical treatment, we use the transformation , which brings Eq. (2) to the following form: where is the first zero of the LE function.

The conditions for solving Eq. (7) then become and .

The use of transformations and in Eq. (4) brings it to the following form:

A description of the RAKCeM(4,4) method

The RAKCeM(4,4) method

Let us consider the following ODE: with the initial condition , , where is sufficiently differentiable function in a neighbourhood of the exact solution (). With a view to solving Eq. (9) effectively, an embedded RK technique can be employed [1], [2], [52]. A general s–stage RK pair can be characterized by means of the extended Butcher tableau of parameters, tabulated in Table 2 [32].

Table 2

The Butcher tableau.

C	A
	bT
	bˆ^T
	ET

In Table 2, , , and and . Then the respective two approximations and at using the two methods can be expressed as [52] where represents the step size, and the slopes are given by , ,  .

It is to be noted that from the embedded form, one can estimate the local truncation error (τ) by , which can be used to control h [53].

The butcher array form for the 4-stage technique can be given in a form presented in Table 3.

Table 3

The Butcher array form corresponding to the 4-stage technique.

0
c_2	a21
c_3	a31	a32
c_3	a41	a42	a43
	b1	b2	b3	b4

The Butcher tableau form for the RKAM(4,4) method can be presented in Table 4 as [52]

Table 4

The Butcher array form of the RKAM(4,4) technique.

0
12	12
12	0	12
1	0	0	1	0
1	16	13	13	16

and the equivalent equations defining the RKAM(4,4) method as [54] with the slopes The RKCeM(4,4) method is [55] with the slopes In Butcher array form, the RKCeM(4,4) method can be put in Table 5 appended below [55]:

Table 5

The Butcher tableau of the RKCeM(4,4) technique.

0
12	12
12	124	1124
1	112	−25132	2566
	13	13	13

Combination of the RKAM(4,4) method given by Eq. (12) with that of RKCeM(4,4) given by Eq. (13) gives rise to an embedded method, referred to as RKACeM(4,4) method, and can be formulated as under: The Butcher tableau form for the RKACeM(4,4) method is given in Table 6 as [55]

Table 6

The Butcher tableau form for the RKACeM(4,4) technique.

0
12	12
12	0	12
1	0	0	1
	…	…	…
12	124	1124
1	124	−25132	7366
	13	13	13
	13	13	13
	ET

Hence by Butcher array presented in Table 2, and the local truncation error estimation, . To approximate the solution through the RKACeM(4,4) technique, one needs four stages, which share the same set of vectors and using and , but and use while and use .

Error control in the RKACeM(4,4) method

Error estimate (ERREST) for the embedded RKACeM(4,4) technique is obtained through the local truncation errors provided by the RKAM(4,4) and RKCeM(4,4) approaches. The following subsection discusses the local truncation errors.

Local truncation error in the RKACeM(4,4) method

Local truncation error at the point is the difference between the computed value and the value at the point on the solution curve that goes through the point .

The local truncation error, τ, for the RKACeM(4,4) method can be obtained from Eqs. (14) and (15) as that can be used to control h. As in [56], [57], an ERREST for the RK(4,4) method can be given by , where and are constants.

For the RKAM(4,4) technique, τ is given by and that for the RKCeM(4,4) method is conferred by , where and are the respective approximated outcomes given by Eqs. (16) and (17) at , and and are the respective local truncation errors obtained via the methods, RKAM(4,4) and RKCeM(4,4). An ERREST for the outcomes at the point by the methods is given by . The τ includes the y derivatives of the RKAM(4,4) method and can be set as [55] while the same for the RKCeM(4,4) technique is given by

Then one can see using Eqs. (18) and (19) that As in [56], if it is assumed that the following bounds for f and its partial derivatives hold for and , one can obtain where p represents the order of the method, and and are constants.

In this case, we have . Hence, with the use of Eq. (21), one can see [55] Thus, from Eqs. (20) and (21), we have .

Hence, .

Now, if the tolerance is assumed to be , then by setting , the EC and the selection of step size can be established to yield

Algorithm for solving 1st order IVP by RKACeM (4,4) technique

The algorithm for finding an approximate solution to the IVP given by Eq. (9) is presented below:

INPUT: endpoints ; initial condition α; tolerance ToL; initial step size h.

OUTPUT: x, w, h, where w approximates , and the step size h was used.

Step 1 set , , .

OUTPUT ().

Step 2 while (FLAG = 1) do steps 3-11.

Step 3 set

Step 4 set

Step 5 If , then do steps 6 and 7.

Step 6 set

Step 7 OUTPUT ()

Step 8 set

Step 9 if then set

else if then set .

Step 10 (The procedure is complete.)

STOP.

Results and discussion

In the case of polytropic protoplanets, their initial structures are directly dependent on the solution of Eq. (7), which necessitates the specification of the parameter n. It is found that the higher is the value of n, the greater is the value of central temperature [43]. But an initial protoplanet should have a low central temperature and hence n is supposed to be small [42]. Following Paul et al. [43], the four values of n, namely 0, 0.5, 1.0, and 1.5 are considered. It is to be noted down here that for all the values of n, LE equation does not possess analytic solutions. So, one may rely on numerical techniques. It is also to be noted down here that for and , the LE equation yields the respective analytic solution, and . It is perceived that for , the solution is monotonically decreasing towards the surface which is physically reasonable. This is also true for . However, for solving Eq. (7), for all the assumed values of n, numerical calculations can be carried out. In our calculation, we have used the embedded RKACeM(4,4) technique to integrate Eq. (7) effectively for all the assumed values of n in estimating the distribution of θ as well as for . Then with the help of Eq. (3), the initial configurations of the thermodynamic variables are attained. The remaining mass distribution was obtained by Eq. (8) using the distribution of for the mentioned x. Our calculated initial profiles of the thermodynamic and other variables for different n are shown in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9.

Figure 1

The distribution of θ for different n.

Figure 2

Comparison of our computed results with exact solutions for n = 0 and n = 1; (a) for n = 0 and (b) for n = 1.

Figure 3

Temperature distribution inside the PPs with masses 0.3 M and 10 M for different n; (a) with 0.3 M and (b) with 10 M.

Figure 4

Temperature distribution inside some PPs with n = 1.5.

Figure 5

The configurations for pressures inside the PPs with masses 1 M and 7 M for assumed n; (a) with 1 M and (b) with 7 M.

Figure 6

Pressure distribution inside some PPs with n = 1.

Figure 7

Density distribution inside the PPs with masses 0.3 M and 10 M for the assumed values of n. (a) with 0.3 M and (b) with 10 M.

Figure 8

Density distribution inside the assumed PPs with n = 1.5.

Figure 9

Mass distribution inside a PP with mass 5 M for the assumed values of n.

Fig. 1 shows the distribution of θ that immediately regulates the distribution of pressure, temperature, and density for different n. The results can be realized to be acceptable to that achieved in [42], [43]. The computed distribution of θ for and with the analytic solutions is illustrated in Fig. 2. We have also estimated the respective RMSE (root mean square error) values to show the accuracy of our obtained results and are presented in Table 7, from which it can be observed that our computed outcomes agree fairly well with the exact solutions.

Table 7

RMSE values of our computed results with that of the exact solution.

Polytropic index	RMSE value
n = 0	2.9860 × 10−9
n = 1	2.4168 × 10−7

Fig. 3 delineates temperature profiles inside the PPs having masses 0.3 and for the assumed values of n. It is seen from Fig. 3 that a higher value of n leads to a higher central temperature of a PP. Our computed temperature distribution that came through our numerical experiment compares well with those obtained in [4], [32], [43], [54]. Fig. 4 displays initial profiles of the temperature inside the assumed PPs for . Fig. 4 shows that a massive protoplanet has a hotter interior. Our estimated profiles for temperature for all the considered values of n inside the assumed PPs compare well with those attained in [31], [32], [39], [43]. It is also interesting to note here that our estimated temperature profile for one Jupiter mass PP with the polytropic equation of state for agrees reasonably well with the configuration of temperature attained in the study due to Helled et al. [58] that was conducted with SCvH (Saumon, Chabrier, Van Horn) equation of state.

Fig. 5 signifies the profiles for pressure for the PPs having masses and for all the assumed values of n and Fig. 6 depicts the same for the PPs with assumed masses for . It is seen from Fig. 5 that the pressure profile of a protoplanet depends on n, for a large n, central pressure of a PP increases, while its surface pressure decreases. The distribution of pressure for with attained in this study is found to be in a reasonable agreement with that attained in [58]. It is aforementioned that the study due to Helled et al. [58] was conducted with SCvH equation of state. However, on the other hand, Fig. 6 shows that as massive as the PP is taken within the assumed mass so higher the central pressure is attained, except for the PP of mass , which is an excellent agreement with those of Paul et al. [4], [32], [40], [43].

Fig. 7 represents density distribution inside the initially formed PPs of masses as well as for all n, and Fig. 8 depicts the distribution of density inside the assumed PPs for .5. Fig. 7 shows that a higher value of n yields a higher central condensation. It is noteworthy that relates to homogeneous models, it yields a sphere of constant density. Also, it can be observed from Fig. 7 that for , the distributions are flatter and almost similar to a constant density model. On the other hand, one can observe from Fig. 8 that the PP with mass (Saturn) is denser over the and PPs, and the PP is less denser over the and PPs at a specific distance from the centre, which is found to be an excellent agreement with the corresponding findings attained in Paul et al. [4], [32], [37], [43]. It is also apparent from the figures (Figure 7, Figure 8) that in our model, the distribution of matter is not uniform, which is as to be expected. However, for and , the density profiles attained in this study for the considered PPs compare fairly well with those presented in [4], [43], [54]. But our attained density distribution considerably differs from those obtained in [31]. In reality, initial profiles of the protoplanets formed through the GI are still unidentified, and diverse initial configurations for them can be found to be predicted by different numerical models [32].

Fig. 9 illustrates our calculated mass distribution in a PP for the assumed values of n. The mass distribution for is seemed to be closer to reality, which is also true for . This nature is found to be valid for all the PPs with and 1.5 and as a result of the same reason mentioned above, the corresponding outcomes are excluded. As in [43], if the shock wave is the trigger for fragmentation of the nebula, the protoplanets in their initial stages are expected to be convective. It is noteworthy that for convection, [43]. Though, Figure 2, Figure 5, Figure 7, Figure 9 show that for , the PPs have a little envelope, but the distributions of pressure, temperature, density, and mass are realistic, which are also true for . However, the system possesses a unique solution for all n signifying that the GI scenario is a reasonable paradigm. The outcomes emanated in the study can be significant in studying the evolution of giant planets from the protoplanets (precursior of planets) formed by GI that recognizes gravity as the only force capable of creating structures by accumulating material in space.

For drawing a comparison, we used RKAM(4,4) and RKCeM(4,4) methods in producing the results with that of the RKACeM(4,4) method. To have a better understanding and to present them meaningfully, the initial configurations of the thermodynamic variables inside only PP for produced by the mentioned methods are tabulated in Table 8. The presentation of other such results for other PPs is overlooked, as we think they would only be space-consuming. It can be contingent from Table 8 that the outcomes obtained by dint of the method of our interest are as good as those attained via the RKAM(4,4) and RKCeM(4,4) methods.

Table 8

Comparison of our computed results for varying x (0.99 − 0.001) by the RKAM(4,4), RKCeM(4,4), and the novel RKACeM(4,4) techniques for the configuration of thermodynamic variables inside a 10 M PP with n = 1.

r/R	RKAM(4,4) method			RKCeM(4,4) method			RKACeM(4,4) method
	P (dynes cm−2)	T (°K)	ρ (gm cm−3)	P (dynes cm−2)	T (°K)	ρ (gm cm−3)	P (dynes cm−2)	T (°K)	ρ (gm cm−3)
0.99	4562.3286	1191.4634	1.6132e-08	4562.3286	1191.4634	1.6132e-08	4561.8384	1191.3994	1.6131e-08
0.9	5161.3694	1267.2722	1.7158e-08	5161.3694	1267.2722	1.7158e-08	5160.8885	1267.2131	1.7157e-08
0.8	5813.1070	1344.9049	1.8209e-08	5813.1070	1344.9049	1.8209e-08	5812.6480	1344.8518	1.8208e-08
0.7	6435.2524	1415.0449	1.9159e-08	6435.2524	1415.0449	1.9159e-08	6434.2164	1414.9310	1.9157e-08
0.6	7010.8850	1476.9774	1.9997-08	7010.8849	1476.9774	1.9997e-08	7009.9375	1476.8776	2.0000e-08
0.5	7525.4508	1530.2193	2.0718e-08	7525.4508	1530.2194	2.0718e-08	7525.1043	1530.1841	2.0718e-08
0.4	7966.6785	1574.4399	2.1317e-08	7966.6784	1574.4398	2.1317e-08	7966.3915	1574.4115	2.1317e-08
0.3	8322.8260	1609.2475	2.1788e-08	8322.8260	1609.2475	2.1788e-08	8322.3033	1609.1969	2.1787e-08
0.2	8583.6253	1634.2662	2.2127e-08	8583.6253	1634.2662	2.2127e-08	8583.4766	1634.2521	2.2125e-08
0.1	8743.1226	1649.3800	2.2332e-08	8743.1226	1649.3800	2.2332e-08	8743.0325	1649.3715	2.2331e-08
0.01	8796.2695	1654.3854	2.2399e-08	8796.2695	1654.3854	2.2399e-08	8796.2682	1654.3853	2.2399e-08
0.001	8796.8249	1654.4376	2.2400e-08	8796.8249	1654.4376	2.2400e-08	8796.8249	1654.4376	2.2400e-08

With the intention of comparing the computational efficiency in producing the results by the methods, the codes with the mentioned methods were run on the same computer (4th generation Intel(R) Core(TM) i5-4570 processor) setting various initial time steps with different starting points, namely 0.05, 0.01, and 0.001. The cost of computation for the RKACeM(4,4) framework was noticed to be slightly higher than that of the RKAM(4,4) and RKCeM(4,4) methods for each scenario except when x is taken as 0.05 and 0.01. Regarding computing factors, the outcomes of our calculations strike the findings of Paul and Senthilkumar [2]. The selection of tolerance may be the reasoning for the fact. Because, for the embedded approach such as RKACeM(4,4), the adapted step size can be resized, if required, based on the accuracy of the outcomes at each step of the solution process. In such a case, it may necessitate repeated calculations before the intended result is obtained therein. However, there is no way to change the step size in either the RKAM(4,4) or RKCeM(4,4) approaches, and therefore the current step size is used in a specific integration domain. Table 9 elucidates the above-mentioned point that the RKACeM(4,4) requires fewer time steps, but it takes a fraction of a second longer to compute the results than the RKAM(4,4) and RKCeM(4,4) methods do. But the time taken by the suggested model for the computing purpose might not be major bargaining. Nevertheless, it is worth noting that the unknowns of a system of ODEs are interdependent. Hence, if an error is sustained in the case of an unknown, it will affect the estimates of the unknowns in the next subsequent steps, and errors can thereby be piled up. Depending on how the tolerance is set up in the RKACeM(4,4) method, an outcome with a specified precision may be obtained. But there is not such an advantage in both the RKAM(4,4) and RKCeM(4,4) methods. The findings produced by the RKAM(4,4) and RKCeM(4,4) methods rely completely on the step size selection, where there is no way to advance the size of the step. This contributes to an escalation in the overall number of time steps, which can tend to increase computational error. Therefore, in terms of accuracy, the method chosen in this analysis is considered to be superior to any of the two methods.

Table 9

Comparison of our computed results for central values calculated by the RKAM(4,4), RKCeM(4,4), and RKACeM(4,4) methods. The calculations are made for n = 1 and a PP with 1 M with different initial time steps. Starting values are different but the calculations are made upwards to the point 0.99 in each case.

Method	Initial time step and starting values	Total step needed	Computational time (second)	P (dynes cm−2) at the endpoint	T (°K) at the endpoint
RKAM(4,4)	0.05	9400	0.0061	600.5622	227.2759
	0.01	9800	0.0102	612.0836	229.4456
	0.001	9890	0.0133	612.9525	229.6084
RKCeM(4,4)	0.05	9400	0.0071	600.5622	227.2759
	0.01	9800	0.0119	612.0836	229.4456
	0.001	9890	0.0158	612.9525	229.6084
RKACeM(4,4)	0.05	25	0.0012	600.4087	227.2468
	0.01	31	0.0078	612.0151	229.4327
	0.001	33	0.0160	612.8868	229.5961

In order to test the efficiency of the RKACeM(4,4) method, ERRESTs for θ are computed while estimating them from the point 0.001 upwards to 0.99 for a PP with and are displayed in Fig. 10. The figure depicts the effectiveness of the method adopted in this analysis. It is worth mentioning here that in the case of both the RKAM(4,4) and RKCeM(4,4) methods, there is no way of quantifying errors.

Figure 10

ERRESTs in θ. Here a polytropic protoplanet having mass 1 M with n = 1 is considered.

We have carried out our experiments for varying endpoints using all of the methods we considered to see how well they worked in solving endpoint constraints. Major discrepancies in outcomes are noticed at the presumed endpoints while using the RKAM(4,4) and RKCeM(4,4) methods, but the RKACeM(4,4) method yields consistent results in the case of the variable endpoints. The results while estimating the distributions of temperature and pressure of a PP with mass for are displayed in Fig. 11 for a deeper interpretation. Other similar figures are omitted for brevity. But the same analyses may be carried out with other protoplanets of presumed masses. Thus, concerning the accuracy, efficacy, and solving endpoint constraints, the RKACeM(4,4) method is found to be more appropriate for solving LE than that of the established RK(4,4) and RKCeM(4,4) methods. This embedded method can raise the computational expense somewhat regarding the choice of the starting step size and tolerance, which is its one limitation. But in today's high-performance machines, computational cost as estimated via the experiment cannot be a serious concern. However, the computational cost can be reduced by setting up a suitable time step, starting point, and tolerance.

Figure 11

Profiles for P and T with varying endpoints obtained via the RKAM(4,4), RKCeM(4,4), and RKACeM(4,4) methods for a polytropic protoplanet with mass 10 M for the polytropic index n = 1.

We have analyzed the stability of the RKAM(4,4) and RKCeM(4,4) techniques. It is to be noted here that while solving an IVP in ODEs numerically, at each integration step, a local truncation error is induced due to the incorrectness of the formula. The GTE (global truncation error) of a numerical method may become massive for the intensification and growth of the local truncation error, even when the error is little at each integration step of the numerical method. The growth phenomenon of the GTE in the numerical method is known as the numerical instability of that method. The stability regions of the RKAM(4,4) and RKCeM(4,4) methods are examined wherein, as a test equation, the first order IVP, , is employed. Hence, from Eqs. (10) and (11), the stability polynomial for each of the RKAM(4,4) and RKCeM(4,4) techniques is obtained, respectively, as and where .

To regulate the stability regions of the RKAM(4,4) and RKCeM(4,4) methods using Eqs. (22) and (23), respectively, on the complex plane, the following conditions are used:

The stability regions of the RKAM(4,4) and RKCeM(4,4) methods are shown in Fig. 12. From the figure, it can be viewed that the stability region of the RKCeM(4,4) method is slightly bigger than that of the RKAM(4,4) method, which means that larger step sizes can be considered in integrating the problem by the RKCeM(4,4) method over the RKAM(4,4) method, which, in turn, may reduce computational cost.

Figure 12

Stability regions of RKAM(4,4) and RKCeM(4,4) methods.

Conclusion

In this study, the RKACeM(4,4) method, due to its advanced characteristics, is employed to test its efficiency in investigating the initial distribution of thermodynamic and other variables in gas giant protoplanets formed by GI, assuming that the polytropic law holds well in them. The model computed polytropic temperature, on which the distribution of thermodynamic and other variables are dependent, is found to agree well with the analytic solutions for and on the basis of the attained RMSE values. The simulated outcomes by the present investigation are also found to be in a reasonable agreement with some corresponding published results. The advantage of the method used in the study is that there is a scope of using tolerance that helps to control the error occurred in the output data, and the stability region of this method is larger than that of the RKAM(4,4) method. That means the deliberated method is more stable than the RKAM(4,4) method. Therefore, it can be concluded that the RKACeM(4,4) method can be an advantageous as well as efficient to study and analyze the initial structures of the PP formed via the GI in terms of ERREST, stability, solving endpoint constraints as well as computational cost with setting up appropriate initial time step, starting point, and tolerance over the RKAM(4,4) method. The method used in this study, therefore, can be helpful for other research subjects, such as the N-body problem or gravitational instability. The problems require the strict conservation of energy and momentum, which means that a reasonable accuracy of the results attained by the numerical method is essential.

Declarations

Author contribution statement

Gour Chandra Paul: Conceived and designed the experiments; Performed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Mrinal Chandra Barman, Hafijur Rahman: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

No data was used for the research described in the article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.