## 3 Numerical Procedure
### Improving the Accuracy of the Cavity Shape  
The cavity shape can be obtained by integrating Equations (8) and (9) or Equations (11) and (12) numerically. However, previous studies have shown that it is necessary to divide the cavity into two [29; 32; 39] or three [31] parts in order to avoid divergence of Equations (8) and (9) at \(\varphi=\pi/2\). In the present study, when the cavity was divided into two parts, the numerical error was found to increase rapidly as the bubble size decreased (see Figure 6). Hence, the cavity shape was divided into the following three parts or intervals: (I) \(0\leq\varphi<\pi/4\), (II) \(\pi/4\leq\varphi<3\pi/4\), and (III) \(3\pi/4\leq\varphi<\pi\), and Equations (8) and (9) were applied to intervals (I) and (III), while Equations (11) and (12) were applied to interval (II). However, because \(\varphi\) is not an independent variable in these equations, the numerical integrations are not completely finished at \(\varphi=\pi/4\) and \(\varphi=3\pi/4\). Hence, a conventional method is to set an arbitrary far point \(X_{f}\) (or \(Y_{f}\)) as a boundary condition of the ODEs and to finish the calculation when \(\varphi\) is just greater than \(\pi/4\) (or \(3\pi/4\)). Unfortunately, this leads to potential errors at the endpoint \(X_{45}\) (or \(Y_{135}\)) of the divided cavity depending on ODE solvers. The conventional method works well on medium-sized bubbles, but produces a severe error in large bubbles, especially at the position of \(X_{45}\), because the bottom of the cavity becomes substantially flat. This error is not only affected by the boundary points \(X_{f}\) (or \(Y_{f}\)), but also by the tolerance of the ODE solver. To remedy this problem, an improved method was proposed for accurately finding the location of the connection points \(X_{45}\) and \(Y_{135}\) regardless of the ODE solvers and their boundary points. In this approach, a monotonically increasing function \(f(x;\varphi)\) is defined as follows:  
```
function\(f(x;\varphi)\)  Integrate Equations (8) and (9) from \(X=0\) to \(x\) to find \(\sin\varphi\) at the endpoint return\((\sin\phi-\sin\varphi)\) endfunction
```  
An accurate position \(x=X_{45}\) is then found by solving \(f(x;\pi/4)=0\) in the interval \(0<x<X_{f}\) using a root-finding technique. Here, the boundary point \(X_{f}\) is any value greater than \(X_{45}\); usually, \(X_{f}\approx 100\) is appropriate for most of the bubble sizes considered in this study. As shown in Figure 2b, \(X_{f}\) does not affect the accuracy of the solution obtained using this method at all. Once \(X_{45}\) has been determined, Equations (8) and (9) are integrated again from \(X=0\) to \(X_{45}\) to obtain the first part of the cavity profile along with the remaining \(Y_{45}\) and \(\varphi_{45}\). Similarly, the second part of the cavity and its endpoint \((X_{135},\ Y_{135},\ \varphi_{135})\) can be computed using the following function:  
```
function\(g(y;\phi)\)  Integrate Equations (11) and (12) from \(Y=Y_{45}\) to \(y\) to find \(\cos\varphi\) at the endpoint return\((\cos\varphi-\cos\varphi)\) endfunction
```  
**Algorithm 2** Function \(g(y;\phi)\)  
The computed bubble shapes according to the position of the boundary point \(X_{f}\) for a large bubble (\(B=10^{-16}\)) are presented in Figure 2. Here, the bubble shapes are seen to vary significantly depending on \(X_{f}\) when the conventional technique is used (Figure 2a). When the root-finding method is used, however, the bubble shapes not affected by \(X_{f}\) (Figure 2b). Similarly, the effect of ODE solver tolerances for the same size of bubble are presented in Figure 3. Here, the bubble shape changes slightly according to the tolerance using the  
Figure 2: Comparison of the bubble shapes obtained according to the boundary point \(X_{f}\) of Equations (8) and (9) for a large bubble (\(B=10^{-16}\) ) using (**a**) the conventional method, and (**b**) the improved method involving the root-finding technique. The tolerance of the ODE solver is set to \(10^{-9}\).