On the finite-size effects in two segregated Bose-Einstein condensates restricted by a hard wall

TL;DR

This paper investigates finite-size effects in two segregated Bose-Einstein condensates confined by a hard wall, analyzing boundary conditions, stability, and phase transitions using the Gross-Pitaevskii equations within the double-parabola approximation.

Contribution

It identifies all possible boundary conditions and their effects on condensate stability and phase transitions, highlighting the stability of Neumann BC and the wetting transitions from Robin and Dirichlet BCs.

Findings

Neumann boundary condition yields a stable ground state.

Robin and Dirichlet boundary conditions lead to unstable states.

Wetting phase transitions occur from the unstable states, with Robin BC being more favorable.

Abstract

The finite-size effects in two segregated Bose-Einstein condensates (BECs) restricted by a hard wall is studied by means of the Gross-Pitaevskii equations in the double-parabola approximation (DPA). Starting from the consistency between the boundary conditions (BCs) imposed on condensates in confined geometry and in the full space, we find all possible BCs together with the corresponding condensate profiles and interface tensions. We discover two finite-size effects: a) The ground state derived from the Neumann BC is stable whereas the ground states derived from the Robin and Dirichlet BCs are unstable. b) Thereby, there equally manifest two possible wetting phase transitions originating from two unstable states. However, the one associated with the Robin BC is more favourable because it corresponds to a smaller interface tension.

Equations

$$L{A}_j={\xi }_j/\sqrt{K-1}$$

$${\xi }_j=\hslash /\sqrt{2m_j{\mu }_j}$$

$$H=\underset{V}{\int }\mathrm{d}V{H}_{\text{GP}}+\underset{W}{\int }\mathrm{d}S{H}_{\text{W}},$$

$$H_{\text{GP}}=\sum _{j=1}^2{\psi }_j^{\ast }\left(-\frac{{\hslash }^2}{2m_j}{\nabla }^2\right){\psi }_j+V\left({\psi }_1,{\psi }_2\right),$$

$$V\left({\psi }_1,{\psi }_2\right)=\sum _{j=1}^2\left(-{\mu }_j{\psi }_j^{\ast }{\psi }_j+\frac{g_{jj}}{2}{\left|{\psi }_j\right|}^4\right)+g_{12}{\left|{\psi }_1\right|}^2{\left|{\psi }_2\right|}^2,$$

$$H_{\text{W}}=\sum _{j=1}^2\frac{{\hslash }^2}{2m_j{\lambda }_j}{\psi }_{\text{W}j}^{\ast }{\psi }_{\text{W}j},$$

$${\psi }_1(\vec{a},z=-h^{\prime })=0$$

$${\psi }_2\left(\vec{a},z=-h^{\prime }\right)=0,$$

$${\psi }_2\left(\vec{a},z=-h^{\prime }\right)\ne 0.$$

$$\underset{h^{\prime }\to \infty }{\lim }{\psi }_2\left(\vec{a},z=-h^{\prime }\right)={\psi }_2\left(\vec{a},-\infty \right)=0,$$

$$H_{\text{W}}=\frac{{\hslash }^2}{2m_2{\lambda }_2}{\psi }_2^{\ast }\left(\vec{a},z=-h^{\prime }\right){\psi }_2\left(\vec{a},z=-h^{\prime }\right)$$

$$\displaystyle\bigg{(}-\frac{\hbar^{2}}{2m_{1}}\Delta-\mu_{1}+g_{11}|\psi_{1}|^{2}+g_{12}|\psi_{2}|^{2}\bigg{)}\psi_{1}=0,$$

$$\displaystyle\bigg{(}-\frac{\hbar^{2}}{2m_{2}}\Delta-\mu_{2}+g_{22}|\psi_{2}|^{2}+g_{12}|\psi_{1}|^{2}\bigg{)}\psi_{2}=0,$$

$$\vec{n}\nabla {\psi }_{\text{W}j}=\frac{1}{{\lambda }_j}{\psi }_{\text{W}j},$$

$${\left(\frac{\partial {\psi }_2\left(\vec{a},z\right)}{\partial z}\right)}_{z=-h^{\prime }}=\frac{1}{{\lambda }_2}{\psi }_2\left(\vec{a},z=-h^{\prime }\right).$$

$${\left(\frac{\partial {\psi }_2\left(\vec{a},z\right)}{\partial z}\right)}_{z=-h^{\prime }}=0,$$

$${\psi }_2\left(\vec{a},z=-h^{\prime }\right)=0,$$

$${\psi }_j={\psi }_j(z),$$

$$\displaystyle H_{\text{GP}}=\sum\limits_{j=1}^{2}\psi_{j}\bigg{(}-\dfrac{\hbar^{2}}{2m_{j}}\frac{\mathrm{d}^{2}}{\mathrm{d}z^{2}}\bigg{)}\psi_{j}+V(\psi_{1},\psi_{2})$$

$$\displaystyle V(\psi_{1},\psi_{2})=\sum\limits_{j=1}^{2}\bigg{(}-\mu_{j}\psi_{j}^{2}+\frac{g_{jj}}{2}\psi_{j}^{4}\bigg{)}+g_{12}\psi_{1}^{2}\psi_{2}^{2}.$$

$$\tilde{H}\left({\varphi }_1,{\varphi }_2\right)=\sum _{j=1}^2{\varphi }_j\left(-\frac{{\xi }_j^2}{{\xi }_1}\frac{{\mathrm{d}}^2}{\mathrm{d}{\varrho }^2}\right){\varphi }_j+\tilde{V}\left({\varphi }_1,{\varphi }_2\right),$$

$$\tilde{V}\left({\varphi }_1,{\varphi }_2\right)=\sum _{j=1}^2\left(-{\varphi }_j^2+\frac{{\varphi }_j^4}{2}\right)+K{\varphi }_1^2{\varphi }_2^2.$$

$$-\frac{{\mathrm{d}}^2{\varphi }_1}{\mathrm{d}{\varrho }^2}-{\varphi }_1+{\varphi }_1^3+K{\varphi }_2^2{\varphi }_1=0,$$

$$-{\xi }^2\frac{{\mathrm{d}}^2{\varphi }_2}{\mathrm{d}{\varrho }^2}-{\varphi }_2+{\varphi }_2^3+K{\varphi }_1^2{\varphi }_2=0,$$

$$H_{\text{W}}$$

$${\tilde{V}}_{\text{DPA}}({\varphi }_j,{\varphi }_2)=2({\varphi }_j-1{)}^2+{\beta }^2{\varphi }_{j^{\prime }}^2-\frac{1}{2},\{\begin{array}{r}\varrho >\ell ,(j,j^{\prime })=(1,2),\\ \varrho <\ell ,(j,j^{\prime })=(2,1),\end{array}$$

$${\left(\frac{\mathrm{d}{\varphi }_j\left(\varrho \right)}{\mathrm{d}\varrho }\right)}_{\varrho =\ell -0}={\left(\frac{\mathrm{d}{\varphi }_j\left(\varrho \right)}{\mathrm{d}\varrho }\right)}_{\rho =\ell +0},$$

$${\varphi }_j\left(\varrho =\ell -0\right)={\varphi }_j\left(\varrho =\ell +0\right).$$

$${\varphi }_1\left(\varrho =\ell ;\ell ,K,\xi \right)={\varphi }_2\left(\varrho =\ell ;\ell ,K,\xi \right),$$

$$\ell =f\left(K,\xi \right).$$

$$-\frac{{\mathrm{d}}^2{\varphi }_1}{\mathrm{d}{\varrho }^2}+2({\varphi }_1-1)=0,$$

$$-{\xi }^2\frac{{\mathrm{d}}^2{\varphi }_2}{\mathrm{d}{\varrho }^2}+{\beta }^2{\varphi }_2=0.$$

$$-\frac{{\mathrm{d}}^2{\varphi }_1}{\mathrm{d}{\varrho }^2}+{\beta }^2{\varphi }_1=0,$$

$$-{\xi }^2\frac{{\mathrm{d}}^2{\varphi }_2}{\mathrm{d}{\varrho }^2}+2({\varphi }_2-1)=0.$$

$${\varphi }_1\left(-h\right)=0.$$

Full text

Key words: finite-size effect, Bose-Einstein condensates, boundary condition, double-parabola approximation

PACS: 03.75.Hh, 05.30.Jp, 68.03.Cd, 68.08.Bc

Ключовi слова: åôåêòè ñêiíчåííîãî ðîçìiðó, êîíäåíñàòè Áîçå-Åéíøòåéíà, ãðàíèчíà óìîâà, íàáëèæåííÿ ïîäâiéíî¿ ïàðàáîëè

1 Introduction

The physics of binary Bose-Einstein condensate mixture in infinite space has explored many interesting instabilities. It was proved theoretically that the well-known instabilities in classical fluid also take place in quantum fluid, including the Rayleigh-Taylor instability [1, 2], the Kelvin-Helmholtz instability [3, 4], the Richtmyer-Meshkov instability [5], and the Rayleigh-Plateau instability [6] for immiscible case. For the miscible case, the counter-superflow instability was discovered in [7, 8, 9] and the dynamical instability was recently predicted to occur in both single-component [10] and two-component Bose-Einstein condensates (BECs) [11]. This instability is associated with the superfluidity breaking of BECs and, therefore, it is purely a quantum effect. In the semi-infinite space where the system of two segregated BECs is restricted by a hard wall, its ground state is fully determined by the boundary conditions (BCs) at hard wall. In [12, 13], with the Dirichlet BCs imposed for both condensates, the authors indicated impeccably that the system of interest can undergo a wetting phase transition. In the present paper we explore the unstable states associated with the BCs at a hard wall, and then we proceed to the study of wetting phase transitions. Here, we prove that the wetting phase transition associated with the Robin BC at hard wall is more favourable.

To this end, we make use of the double-parabola approximation (DPA) method which was proposed in [14, 15] and was developed in [16, 17, 18]. The boundary problem we deal with is schematically given in figure 1 where the hard wall (optical wall) is placed at $z=-h^{\prime}$, the interface at $z=L$ and the curve $j$ $(j=1,2)$ represent the condensate profile ${\phi_{j}}(z)$ which will be introduced in the next section. The interval $L{A_{j}}$ is the so-called penetration depth of condensate $j$ [14] inside the condensate $j^{\prime}\neq j.$ Its expression reads

[TABLE]

with

[TABLE]

being the healing length of condensate $j$. The parameters $K$, ${m_{j}}$ and ${\mu_{j}}$ will be cleared up in what follows. Therefore, the interval ${A_{1}}{A_{2}}$ is the coexistence region of two condensates.

Following [19, 20], $H$ can be written in the form

[TABLE]

here, the Gross-Pitaevskii (GP) Hamiltonian ${H_{\text{GP}}}$ in the bulk and the wall Hamiltonian ${H_{\text{W}}}$ are given as follows:

– The GP Hamiltonian

[TABLE]

here, ${\psi_{j}}={\psi_{j}}(\vec{x})$, ${m_{j}}$ and ${\mu_{j}}$ $(j=1,2)$ are respectively the wave function, the atomic mass and the chemical potential of each species $j$, with $\vec{x}=(\vec{a},z)$, $\vec{a}=(x,y)$ to be denoted in what follows. For inhomogeneous state, in which two condensates are separated by interface, we have $K=\frac{g_{12}}{\sqrt{{g_{11}}{g_{22}}}}>1$ and hereafter we will consider the case of all coupling constants being positive.

– The wall Hamiltonian. Usually, the wall Hamiltonian is chosen in the phenomenological form [19, 20]

[TABLE]

where ${\psi_{\text{W}j}}$ $(j=1,2)$ is the wall field induced by ${\psi_{j}}$ on the wall and ${\lambda_{j}}$ is its extrapolation length, introduced by de Gennes [21]. Next, let us pay special attention to the wall fields when the wall is flat. For the first condensate, figure 1 indicates that the Dirichlet BC

[TABLE]

is justified when $h^{\prime}+L\geqslant\frac{\xi_{1}}{\sqrt{K-1}}$ or $K\geqslant\frac{{\xi_{1}^{2}}}{{{{\left({h^{\prime}+L}\right)}^{2}}}}+1$.

Moreover, the Dirichlet BC is also consistent with the condition in the infinite space ${\psi_{1}}({\vec{a},-\infty})=0$ [16, 18]. As to the second condensate, we have two options:

either

[TABLE]

or

[TABLE]

It is easily seen that the consistency between the BCs in the semi-infinite space and in the infinite space eliminates the first option. Indeed, if $h^{\prime}\to\infty$ one obtains

[TABLE]

which contradicts the BC for this condensate in the infinite space [14, 22], ${\psi_{2}}\left({\vec{a},-\infty}\right)\neq 0.$

Hence, in reality the wall Hamiltonian (1.3) is simplified to

[TABLE]

for the flat wall.

The fact that the equilibrium values of the fields minimizes the total Hamiltonian given in (1.1), (1.2) and (1.4) leads to the time-independent GP equations in the bulk

[TABLE]

together with the BCs at the wall

[TABLE]

where $\vec{n}$ is the unit vector normal to the wall and pointing inside the system. For a flat wall, the foregoing condition turns out to be the Robin BC

[TABLE]

It is easily seen that (1.6) is reduced to the Neumann BC

[TABLE]

when ${\lambda_{2}}$ tends to infinity and to the Dirichlet BC

[TABLE]

when ${\lambda_{2}}\to 0$.

The present paper is organized as follows. In section 2, employing the DPA, we calculate the analytical solutions of the GP equations obeying the Robin BC at a hard wall. Section 3 is devoted to the determination of the interface tension in the grand canonical ensemble (GCE) and the wetting phase transition. The conclusion is given in section 4.

2 Ground state in DPA

For two segregated BECs, the translation along the $Oz$ axis is spontaneously broken by the presence of the domain wall-the interface, and one usually focuses on the condensate profiles that only depend on $z$ and for simplicity we do not change their symbols

[TABLE]

which are assumed to be real without loss of generality. In terms of ${\psi_{j}}\left(z\right)$, the GP Hamiltonian (1.2) becomes

[TABLE]

It is convenient to write all quantities in the dimensionless forms by introducing new quantities $\varrho=z/{\xi_{1}}$, $h=h^{\prime}/{\xi_{1}}$, $\ell=L/{\xi_{1}}$ and ${\phi_{j}}={\psi_{j}}/\sqrt{{n_{j}}},$ where the densities ${n_{j}}$ $(j=1,2)$ were determined in [16]. In the dimensionless forms, the GP Hamiltonian (2.1) is rewritten as follows:

[TABLE]

In the dimensionless form, equations (1.5) read

[TABLE]

here, $\xi={\xi_{2}}/{\xi_{1}}$.

Correspondingly, the wall Hamiltonian (1.4) is rewritten in the form

[TABLE]

Now, let us present in detail the formulation of the DPA [16]. The DPA method consists of three constituents:

a) We temporarily introduce $\varrho=\ell$ the interface location which separates the system into two regions as plotted in figure 1.

b) The foregoing potential (2.2b) is linearized by the DPA potential ${V_{\text{DPA}}}({{\phi_{1}},{\phi_{2}}})$ in appropriate half spaces separated by the interface

[TABLE]

where $\beta=\sqrt{K-1}$.

c) The smooth variation and the continuity of both condensate profiles when passing the interface are assumed. Mathematically, this assumption is expressed by

[TABLE]

It is important to note that the introduction of the interface position $\ell$ is merely a mathematical tactics based on which one sets up the DPA potential. Indeed, the interface is defined as the intersection of two condensates

[TABLE]

which indicates that $\ell$ is a function of $K$ and $\xi$

[TABLE]

Therefore, in reality, both condensates depend only upon $\varrho$ and the model parameters.

Taking into account (2.2a), (2.4) and (2.5), we arrive at the total Hamiltonian in the DPA, ${H_{\text{DPA}}}$. Minimizing it, we are led to the GP equations in DPA together with the BCs, namely:

– GP equations.

Equations in the right-hand side of the interface

[TABLE]

– The BCs that both condensate profiles fulfil are:

The Dirichlet BC at hard wall for the first condensate

[TABLE]

The Robin BC at hard wall for the second condensate

[TABLE]

with $c=\frac{{{\xi_{1}}}}{{{\lambda_{2}}}}.$

Besides these, we still have the conditions for both condensates at infinity

[TABLE]

The analytical expressions of the solutions of equations (2.10) which obey the conditions (2.11), (2.12) and (2.13) are in order:

– In the RHS of the interface

[TABLE]

– In the LHS of the interface

[TABLE]

Substituting (2.14), (2.15) into (2.6), (2.8), we find the analytical expressions of ${A_{j}}$, ${B_{j}}$ $(j=1,2)$ which are given in appendix A.

The most important requirement for any approximation method is its reliability. To check the reliability of DPA, let us plot in figure 2 the condensate profiles corresponding to several values of $\xi$ and at $c=1/\sqrt{2}$, $h=0$: a) $\xi=1.0$, $K=3$; b) $\xi=0.5$, $K=2.5$; c) $\xi=1.5$, $K=3.5$. The solid lines denote the solutions (2.14), (2.15) of the DPA equations (2.10) and the dashed lines denote the numerical solutions of the GP equations (2.3), with the same BCs (2.11), (2.12) and (2.13).

Figure 2 tells us that:

– The graphs obtained in DPA are close to those found in the GP theory in the whole variation region of $\varrho$. This implies that the DPA is a reliable method.

– The condensates calculated respectively in DPA and in GP equations provide the same interface position $\ell$.

Note that the foregoing results remain valid for other values of $c$, including $c=0$ (Neumann BC) and $c=\infty$ (Dirichlet BC).

– Finally, let us calculate the value $c$ appearing in the Robin BC (2.12). At first, let us look for ${\phi_{2}}\left({-h}\right)$ from (2.15)

[TABLE]

The expressions of ${M_{j}}$, $1\leqslant j\leqslant 4$ are given in appendix B.

When $h$ tends to infinity, (2.16) becomes

[TABLE]

which must coincide with the BC for ${\phi_{2}}\left(\varrho\right)$ in the infinite space [14, 16], ${\phi_{2}}\left({-\infty}\right)=1$. Hence, we get $c=0$ and then the Robin BC (2.12) turns out to be the Neumann BC.

3 Interface tension in grand canonical ensemble

The interface tension was determined in grand canonical ensemble (GCE) and canonical ensemble (CE) by [23, 24] and [14], accordingly. The system can be viewed as having a direct contact with the bulk reservoirs of condensates and we have ${\mu_{j}}={g_{jj}}{n_{j}}$. The interface tension defined in GCE are equal to $4$ times their counterparts in CE [13] in the infinite space and different from $4$ times in the semi-infinite space [17]. To begin with, let us follow closely [23, 24] to start from the grand potential of the interface in the dimensionless form calculated at the bulk coexistence of two condensates

[TABLE]

where $P_{1}=P_{2}=P_{0}=g_{jj}n_{j}^{2}/2$, $A$ is the interface area and $\tilde{V}$ is given by (2.2b).

Replacing $\tilde{V}$ by the DPA potential (2.5) in equation (3.1), we arrive at the grand potential in DPA

[TABLE]

Then, taking into account (2.5), (2.10) and (3.2), the interface tension is determined as the excess grand potential per unit area,

[TABLE]

with ${\Omega_{\text{b}}}=-PV$ being the bulk grand potential.

Substituting (2.14), (2.15) into (3.3), leads to

[TABLE]

in which

[TABLE]

(3.4) is the interface tension with Robin BC.

When $c$ tends to infinity, equation (3.4) turns out to be

[TABLE]

Expression (3.5) totally coincides with the interface tension calculated imposing the Dirichlet condition at a hard wall [17]. When $c$ tends to zero, equation (3.4) reduces to the interface tension with Neumann BC

[TABLE]

In order to better understand the physical significance of the Neumann BC it is convenient to let $h=0$ and rewrite the interface tension (3.3) in the following form

[TABLE]

where

[TABLE]

The expression $\int\nolimits_{0}^{\ell}{{\phi_{2}}}\mathrm{d}\varrho$ is exactly the area of the shaded domain limited by the condensate $2$ and the interval $0\leqslant\varrho\leqslant\ell$ in the $\varrho$ axis, as is shown in figure 3. It is easily realized that the area monotonously decreases for $c$ changing from Neumann BC (blue line), $c=0$, to Robin BC (green line), $c=1$, for example, and then to Dirichlet BC (red line), $c\rightarrow\infty$. This implies that the interface tension (3.7) obeys the inequality

[TABLE]

The above inequality is illustrated in figure 4 which displays the $K$ dependence of interface tensions (3.4), (3.5) and (3.6), corresponding, respectively, to the Robin, Dirichlet and Neumann BCs at $\xi=1$. For other values of $\xi$, the above inequality (3.8) remains unchanged (illustrated in figure 5).

Combining the inequality (3.8) and the figure 4, we arrive at the important statement:

a) The state derived from the Neumann BC is stable while those which are derived from the other BCs are unstable.

b)

[TABLE]

Next, let us consider whether or not the system would undergo a transition from immiscible to miscible states when $K$ tends to $1$. To this end, we numerically compute the evolution of the function $\ell=f(K,\xi)$ given in (2.8) versus $K$ and the result is visualized in figure 6, which gives

[TABLE]

For other values of $\xi$, the preceding results remain unchanged. Based on (3.10), we immediately have

[TABLE]

which is valid for all interface tensions (3.4), (3.5) and (3.6).

Combining (3.9) and (3.11) proves that as $K$ tends to $1$, there occurs a first-order phase transition from immiscible to miscible states.

To proceed to the wetting phase transition, let us remark that this phase transition could be realized in the system which was initially in the unstable state. In our case, the wetting phase transition could equally depart from one of the two unstable states as follows:

a) Ground state derived from the Dirichlet BC. We consider the behaviour of the corresponding interface tension given in (3.5) in strong segregation

[TABLE]

On the other hand, for strong segregation, the behaviour ${\gamma_{12}}$ is expressed through the sum of two wall tensions ${\gamma_{\text{W}j}}$ $(j=1,2)$

[TABLE]

Comparing (3.12) and (3.13), immediately gives

[TABLE]

The wetting phase transition takes place when the Antonov rule is realized [13, 17]

[TABLE]

The phase diagram of (3.15) is plotted in figure 7 with blue line.

b) Ground state derived from the Robin BC. Analogously, the strong segregation behaviour of the interface tension derived from the Robin BC (3.4) reads

[TABLE]

with

[TABLE]

Combining (3.13) and (3), yields

[TABLE]

The phase diagram of the wetting phase transition determined by the Antonov rule (3.15) is depicted in figure 7 by red line. There is easily seen a distinction between two diagrams corresponding to two different phase transitions. Although both phase transitions are possible, but the one associated with the Robin BC is more favourable because it corresponds to the smaller interface tension (3.8).

4 Conclusion and outlook

Based on the GP equations in the DPA, we investigated the system of two segregated Bose-Einstein condensates restricted by a hard wall. Our main results are in order:

a) Making use of the consistency between the the BCs imposed on the condensates in confined geometry and in infinite space, we determined the BCs for both condensates at a hard wall, namely, the Neumann BC for the second condensate and Dirichlet BC for the first condensate.

b) The analytical expressions for both condensate profiles, calculated in different BCs, are very close to those derived from the GP equations by numerical computations. This proves the reliability of the DPA.

c) Based on the analytical expressions of interface tension, calculated in different BCs, we showed that

[TABLE]

which proves that the ground state corresponds to the Neumann BC and is stable while other two BCs lead to unstable states.

d) The above mentioned result leads to an important consequence: In our system there exist two scenarios of wetting phase transitions which correspond to two different unstable states, but the one associated with the Robin unstable state is more favourable because it corresponds to a smaller interface tension.

e) As $K$ tends to $1$, the system undergoes a first-order phase transition from immiscible to miscible states.

The results mentioned in items a), c) and d) are our major successes, they are very meaningful because the system of binary BECs mixture is also useful for many other systems of great interest, such as multi-component superconductors, Mott-insulators and so on [25]. Especially, it is very interesting to extend the method presented in this work, in a straightforward manner, to two-gap superconductors [26] whose Hamiltonians have the same structure of our system. It is expected that the application of the DPA to these systems could provide new interesting effects.

Last but not least, it is worth to remark that our previous results were derived within the framework of the Gross-Pitaevskii theory based on the mean field approximation, whereby the thermal and quantum fluctuations were completely neglected [27, 28]. Although a great number of the effects related to the Bose-Einstein condensations can be successfully explained by the Gross-Pitaevskii theory, in recent years a more complete theory is most welcome in order to comprehend various experiments associated with quantum fluctuations, solitons, vortices and so on [29, 30, 31].

Acknowledgements

This paper is supported by the Vietnam Ministry of Education and Training in the framework of the Scientific Research Project under Grant B 2015-25-33. The discussion with Pham The Song is acknowledged with thanks.

Appendix A The analytical expressions of $A_{j}$, $B_{j}$ $(j=1,2)$

[TABLE]

[TABLE]

Appendix B The analytical expressions of ${M_{j}}$, $1\leqslant j\leqslant 4$

[TABLE]

[TABLE]