## 4 Computational Approach
### Numerical Method  
Implicit Large Eddy Simulation (iLES) is the method used in this study. A high-resolution sixth-order compact difference scheme with tenth-order compact filtering is used [31]. As the sixth-order compact differencing scheme has been shown to have 50-100 times higher spatial resolution than conventional second-order schemes, this approach saves on the numbers of computational grids. As a result, this saves computer time and memory [32]. Such a scheme was necessary even with the world's leading-edge supercomputer because even one specific case requires considerable time for computation. In this iLES method, tenth-order filtering plays a role as an alternative to the sub-grid model. For time integration, lower-upper symmetric and alternating-directional implicit symmetric Gauss-Seidel (ADI-SGS) [33] methods are used. To ensure time accuracy, a backward second-order difference formula is used for time integration, and 3-5 sub-iterations are adopted. Details of the numerical algorithms have been previously reported [26]. Freestream Mach numbers in the numerical simulations are set to be 0.1 or 0.2 for computational efficiency, but it only models low-speed flows assuming that the phenomena are essentially the same. The computer program used in this study is our original in-house computer code named LANS3D (originally developed by [PERSON] and [PERSON] in 1984 and modified by [PERSON] et al. at Institute of Space and Astronautical Science, Sagamihara, Japan) that has been used for many applications and is well validated. Recent applications for the flows over airfoils at similar Reynolds numbers have been reported [34; 35]. The code is modified in consideration of the source (body force) terms in Equations (2) and (3).  
Computational grids are carefully prepared to be appropriate for iLES simulations at this range of the Reynolds numbers. Namely, the grid sizes based on the wall unit, which were calculated for the attached flow over the NACA0015 airfoil, are \(\Delta\xi\approx 8\) (in the streamwise direction), \(\Delta\eta\approx 9\) (in the spanwise direction) and \(\Delta\xi\approx 1\) (in the wall-normal direction). Our previous study using the same compact differencing and filtering schemes showed that the grid resolutions to resolve near-wall turbulence are \(\Delta\xi\approx 36\) in the streamwise direction, \(\Delta\eta\approx 15\) in the spanwise direction, and \(\Delta\xi\approx 1\) in the wall-normal direction. Therefore, the grid spacing near the airfoil surface has a sufficient spatial resolution to resolve near-wall turbulence. A grid-independence study was also conducted [36].  
We employed an overset zonal method [37] to enhance the spatial grid resolution in the region around the plasma actuator. There are two zonal grids for the flow simulations: one for the overall airfoil region (zone 1) and the other for the actuator region (zone 2). The latter grid with much fine spatial resolution is overlapped onto the former grid. The body force terms are computed in advance using the grid with ultra-fine spatial resolution. Then, the result is mapped onto the fine grid region (zone 2) in the flow simulations. The span length of the simulation is 20%. The total number of grid points is approximately 20 million.