## 3 Numerical Procedure
### Finding the Junction Point Based on the Type of Meniscus  
There are two general approaches to finding the junction point \(c\) according to the direction of integration of the meniscus. The first approach is to integrate the meniscus from \(X=X_{\infty}\) to \(X_{c}\) in the negative \(X\)-direction [29, 31]. Here, \(X_{\infty}\) is sought in an iterative manner so that the derivatives of the meniscus and cavity profiles match at the junction point \(c\). However, it is difficult to compute the shape of the meniscus using this method because the initial slope becomes zero when the exact initial condition, \(Z_{\infty}=\alpha_{\infty}=0\) at \(X=X_{\infty}\), is applied. To circumvent this, [PERSON] [29] and [PERSON] and [PERSON] [40] employed \(Z_{a}=\alpha_{a}=0.00001\) at \(X=X_{a}<X_{\infty}\) as an initial condition. [PERSON] and [PERSON] [31] further improved this method by using an approximate analytical solution of \(Z_{a}\) with the Bessel function. Nevertheless, as these methods do not give an exact initial condition at \(X=X_{\infty}\), the accuracy of the solution is still limited. The second approach is to integrate the meniscus from \(X_{c}\) to \(X_{\infty}\) in the positive \(X\)-direction [23, 30, 32]. Although this approach has the advantage of obtaining the \(X_{\infty}\) position more precisely, there is another difficulty in that the meniscus shape changes dramatically around the junction point \(c\), as shown in Figure 4.  
Figure 4: Three types of the meniscus shape according to the location of the guessed junction point. Here, \(\varphi_{c}\) is the angle at the exact junction point. Figure 3: Comparison of the bubble shapes obtained according to the tolerance of the ODE solver for a large bubble (\(B=10^{-16}\)) using (**a**) the conventional method, and (**b**) the improved method involving the root-finding technique. The insets plot the cavity profiles up to 45 degrees and are not scaled. Hence, a method for effectively locating the junction point \(c\) is proposed herein by using a specific function that distinguishes the type of meniscus shape. A conventional bisection method for finding \(c\) narrows the interval by comparing the sign of \(H-\widetilde{H}\) at both ends of the solution interval. At each iteration step, it is necessary to find the height \(\widetilde{H}\), which is calculated differently depending on the meniscus type: if the meniscus profile is convex, its minimum value is chosen; if the meniscus profile is monotonically decreasing, an inflection point is selected. Interestingly, without calculating \(\widetilde{H}\), the bisection method can become simpler by considering the meniscus type only. Hence, the following function is defined to determine whether a meniscus is convex at a guessed junction point:  
```
function is_meniscus_convex(\(c\))  Integrate Equations (15) and (16) from \(X=X_{c}\) to \(X_{max}\) to find \(\sin\alpha_{max}\),  where \(X_{max}>X_{\infty}\) is any large value return the logical value of \((\sin\alpha_{c}\cdot\sin\alpha_{max}<0)\) endfunction
```  
**Algorithm 1** Modified bisection method  
This function is then used to modify the bisection method as Algorithm 1:  
```
Input:\(P_{L},P_{R}\) Output:\(P_{c}\) for\(n=1\) to \(n_{max}\) \(P_{c}=(P_{L}+P_{R})/2\) ifis_meniscus_convex(\(P_{c}\)) \(P_{R}=P_{c}\) else \(P_{L}=P_{c}\) endif endfor
```  
**Algorithm 2** Modified bisection method  
where \(P\) represents the \(X\) or \(Y\) value on the cavity. The maximum number of iterations is given by \(n_{max}=\log_{2}\left(\left|P_{L}^{0}-P_{R}^{0}\right|/\epsilon\right)\). Here, \(\epsilon\) and \(\left[P_{L}^{0},\ P_{R}^{0}\right]\) are the error tolerance and initial interval of the bisection method, respectively. If \(\varphi_{c}\) is between \(90^{\circ}\) and \(135^{\circ}\), the initial interval is \([Y_{90},\ Y_{135}]\); if \(\varphi_{c}\) is between \(135^{\circ}\) and \(180^{\circ}\), the initial interval is \([X_{135},\ X_{180}]\), where \(X_{180}\
eq 0\).