Plausible Detection of Feasible Cave Networks Beneath Impact Melt Pits on the Moon Using the Grail Mission

: In the future, when humans build their bases on terrestrial planets and their moons, caves will be the safest place for inhabitation. Large holes, believed to be cave entrances, have been discovered on the Moon, along with small features called “impact melt pits.” In the Gravity Recovery and Interior Laboratory (GRAIL) gravity model, which is expressed in spherical harmonics (SH), it is difﬁcult to express the gravity anomaly created by a small empty space below the surface. Nevertheless, we propose that a cave network, akin to an anthill, exists under the impact melt pits discovered on the Moon. This is because we think it is natural to apply a network created by Earth’s small caves to the Moon. We obtained accurate Bouguer gravity measurements by calculating regional crustal density using localized admittance of the study area and detected weak gravity (mass deﬁcit) information. By increasing the degrees and order of SH at regular intervals, we estimated the change in gravity at a speciﬁc position at high degrees and order, thereby extracting shallow depth information. To validate our method, we compared our results with those of existing studies that analyzed the previously known Marius Hills Hole (MHH) area. The analysis of seven regions in our study area revealed a mass deﬁcit in some impact melt pits in four lunar regions (Copernicus, King, Stevinus, and Tycho). We propose that there is a cave network in this region, indicated by the gravitation reduction in the impact melt pits region. Our results can be useful for the selection of landing sites for future in situ explorations of lunar caves.


Introduction
When humans go into space, the Moon, in particular, a lunar cave, will be suitable as an outpost. There are typically three types of caves on Earth: solutional caves, sea caves, and lava tubes. Solutional and sea caves are created with water as the main medium, and lava tubes are created as a medium for lava. Recent studies have revealed the existence of water on Mars, and spectroscopic observations of the Moon have also revealed trace amounts of water on the lunar surface [1,2]. However, the presence of water on a large enough scale to be considered as a "water resource", such as those observed on Earth, have not been reported on both the Moon and Mars. As suggested by Lauro et al. [1], if there is liquid water on Mars, the existence of solutional and sea caves can also be considered. However, because there is no observational evidence of liquid water on the Moon, we can only consider the presence of lava tubes for investigating the existence of caves. The existence of lava tubes on the Moon has already been predicted from traces of lava activity on the lunar surface [3][4][5]. Lava tube regions have been proposed as suitable sites for human inhabitation on the lunar surface [6,7].
The role of a subsurface base is to defend against the physical damage that humans may experience on the lunar surface. Typical damages that we may experience on the lunar surface include diurnal temperature fluctuations, meteorite impacts, cosmic radiation, high energy particles, and dust. The daily temperature difference of the Moon varies according to latitude, with a change of about 100-400 K at the equator and a constant temperature of the surface with a void height of 75 m roughly. However, Kobayashi et al. [28] denied the results of Kaku et al. [27], with respect to presenting the conditions and data processing methods required when searching for lava tubes using SELENE LRS data; to understand the use of this technology and its use in the detection of cave openings and its internal structure, the results of lunar cave exploration using GPR need to be investigated further.
In this study, we explored a lunar cave via a geophysical approach using the GRAIL gravity exploration data. However, we only targeted impact melt pits on the Moon, and not mares and highland pits, which have been the most common topics of several previous studies, because we think that it is natural to apply the network created by small caves on Earth rather than large single individual caves to the moon.

Data
The GRAIL's gravity observation data are distributed as a gravity model consisting of spherical harmonics (SH) coefficients. From the early model, GRAIL420C1A, to the current model, GL1500E, gravity models of various degrees and orders have been developed over the past years [23,[29][30][31][32][33][34]. In the gravity model, as a characteristic of SH, the higher the degree and order, the shorter the wavelength; therefore, shallow depth information can be effectively expressed using this model [35]. Additionally, when a gravity model is expressed as a map, it has a high spatial resolution, with high degree and order; therefore, such models can be used to deduce the gravity of a small area. The GL1500E model offers the highest degree and order among the currently published gravity models and is suitable for exploring lunar caves. However, Goossens et al. [36] improved the GRGM1200B gravity model (with a degree and order of 1200) and developed the GRGM1200B RM1 model, which provided a high degree and order; notably, GRGM1200B RM1 showed a significantly higher correlation between the observed gravity and the gravity derived from topography, compared to the GL1500E model. Thus, GRGM1200B RM1 expresses shallow depth gravity information more effectively, even though GL1500E has a high spatial resolution. Therefore, in this study, we proceeded with the GRGM1200B RM1 for λ = 1 model, which provides more detailed shallow depth information, even if it compromised the spatial resolution to some extent.

Hypothesis
Before analyzing the gravity model, it is necessary to know how the existence of a small lunar cave changes the gravity field. In this study, we determined whether the information of a small subsurface void can be expressed in a gravity model. Lowrie [37] describes the change in gravity in a simple form, with respect to the differences in the density, radius, and depth of a mass, when an object having a different density from the perceptual density exists below the surface. The equation of Lowrie [37] was used, modified for our study: where R is the radius of the void, x, y is the distance from the center (surface above the void position), z is the depth, ∆ρ is the density variation, and G is the gravitational constant. Based on this, the changes in gravity, with respect to the density of the lunar crust (2550 kg/m 3 ), density of empty space (0 kg/m 3 ), radius (1-13 m), and depth, are shown in Figure 1a. When the radius of the empty space was more than 5 m, the gravity change appeared at a depth of several tens of meters, and this amount of change was expressed in the gravity model effectively. Next, we checked the spatial change in gravity that appeared under specific conditions. Figure 1b portrays a gravity anomaly at void center g 0 when there was a spherical empty space having a radius of 5 m at a depth of 20 m, with the lunar crust as a background. As the distance from the empty space increased, the reduced gravity recovered, and the shape of the wavelength showed a sharp pin, rather than a gentle curve. This wavelength corresponded to SH; notably, the degree and order must be approximately 100,000 or more to be expressed effectively. A small lunar cave corresponding to this degree and order cannot be expressed in the currently used gravity model. Therefore, we hypothesized that if a small lunar cave exists, it will form a network (akin to an anthill), with multiple entrances, rather than resembling a single entity, such as the MHH, MIH, and MTH. If subsurface voids are entangled similar to a net in a certain area, such as the lava tubes observed on Earth, the change in gravity in the area having empty spaces can be sufficiently detected in current the gravity field model. In addition, because the lunar crustal structure is very complex, the small signal of the lunar cave will most likely be expressed in a combination of various high degrees and orders. in the gravity model effectively. Next, we checked the spatial change in gravity that appeared under specific conditions. Figure 1b portrays a gravity anomaly at void center g0 when there was a spherical empty space having a radius of 5 m at a depth of 20 m, with the lunar crust as a background. As the distance from the empty space increased, the reduced gravity recovered, and the shape of the wavelength showed a sharp pin, rather than a gentle curve. This wavelength corresponded to SH; notably, the degree and order must be approximately 100,000 or more to be expressed effectively. A small lunar cave corresponding to this degree and order cannot be expressed in the currently used gravity model. Therefore, we hypothesized that if a small lunar cave exists, it will form a network (akin to an anthill), with multiple entrances, rather than resembling a single entity, such as the MHH, MIH, and MTH. If subsurface voids are entangled similar to a net in a certain area, such as the lava tubes observed on Earth, the change in gravity in the area having empty spaces can be sufficiently detected in current the gravity field model. In addition, because the lunar crustal structure is very complex, the small signal of the lunar cave will most likely be expressed in a combination of various high degrees and orders. ∆ is the surface above the void position where the change in gravity is strongest, and the radius is indicated in the legend by specifying the color from 1 m to 13 m (at intervals of 2 m). Lunar crustal density was set as the background density ( = 2550 kg/m 3 ), and the density of the empty space was set as = 0 kg/m 3 . (b) Gravity anomaly when an empty sphere having radius of 5 m exists at depth of 20 m. Gravity anomaly variation decreases sharply as the distance between the empty space and the surface increases and decreases by 95 % when compared to the point at 100-m depth. Density conditions of (b) are same as those of (a).

Research Target
Our study area is an impact melt pit consisting of a terrain that has a high possibility of a subsurface lunar cave network. Unlike the pits in the mare and highland areas, this terrain is small, has an irregular shape, consists of several clusters that form near fractures, and is distributed in impact melt deposits [14]. In addition, because it has a similar shape to lava tube ceilings observed on Earth, it may connect to subsurface voids [15,38,39]. If there is a cave network beneath the impact melt pits, they can affect the gravity of the impact melt pits area, so we set the analysis target as an impact melt pit. Wagner and Robinson [14] found the presence of impact melt pits in a total of 29 craters on the Moon. We performed the analysis on seven craters, containing more than 10 impact melt pits (Table 1). In the latest lunar impact pit catalog released by Wagner and Robinson [21], there was a change in the number of pits by region. In the case of the Lalande crater, there were eight pits, not ten; however, our targets were based on the results of [14]. ∆g 0 is the surface above the void position where the change in gravity is strongest, and the radius is indicated in the legend by specifying the color from 1 m to 13 m (at intervals of 2 m). Lunar crustal density was set as the background density (ρ 0 = 2550 kg/m 3 ), and the density of the empty space was set as ρ 0 = 0 kg/m 3 . (b) Gravity anomaly when an empty sphere having radius of 5 m exists at depth of 20 m. Gravity anomaly variation decreases sharply as the distance between the empty space and the surface increases and decreases by 95 % when compared to the g 0 point at 100-m depth. Density conditions of (b) are same as those of (a).

Research Target
Our study area is an impact melt pit consisting of a terrain that has a high possibility of a subsurface lunar cave network. Unlike the pits in the mare and highland areas, this terrain is small, has an irregular shape, consists of several clusters that form near fractures, and is distributed in impact melt deposits [14]. In addition, because it has a similar shape to lava tube ceilings observed on Earth, it may connect to subsurface voids [15,38,39]. If there is a cave network beneath the impact melt pits, they can affect the gravity of the impact melt pits area, so we set the analysis target as an impact melt pit. Wagner and Robinson [14] found the presence of impact melt pits in a total of 29 craters on the Moon. We performed the analysis on seven craters, containing more than 10 impact melt pits (Table 1). In the latest lunar impact pit catalog released by Wagner and Robinson [21], there was a change in the number of pits by region. In the case of the Lalande crater, there were eight pits, not ten; however, our targets were based on the results of [14]. Table 1. Seven craters containing more than ten pits [14].

Estimation of Localized Crustal Density
In this study, we used SHTOOLS (https://shtools.oca.eu/shtools/public/index.html (accessed on 15 October 2021)), which is a software developed for a potential field operation (represented by SH) for gravity model analysis [40]. Previous representative studies of the lunar gravity field have mainly focused on subsurface structure and evolution, such as crustal thickness determination and dike exploration for the entire Moon [41,42]. In the case of gravity field analysis for the entire Moon, we set one value of crustal density when calculating the Bouguer gravity, which expressed only the gravity beneath the surface (by removing the influence of topography). However, due to the difference in the composition of the mare and highland crust of the Moon, the actual density varied (depending on the location). This difference affects the Bouguer gravity of a specific area, resulting in inaccurate analysis. Therefore, to analyze the localized Bouguer gravity, it is necessary to estimate the localized crustal density of the corresponding area.
In this study, we calculated the admittance to estimate the perceptual density at the seven craters. Admittance refers to the gravity anomaly generated by a change in topography [43]. Watts [43] explained the relationship between crustal density and admittance using four isostasy models (uncompensated, Airy, Pratt, and flexure). In the form of admittance, as the wavenumber (replaced by the degree and order of spherical harmonics, k = (l + 1/2)/R; k is wavenumber, l is spherical harmonics degree, and R is radius) increased, the three models converged to the uncompensated model. Our target is a small terrain, and we needed to consider high degrees and orders. Therefore, in our study, we applied the uncompensated model, which is the simplest form and does not require assumptions about parameters to fit the Moon. The equation of Watts [43] was used: where ρ c is oceanic crustal density, ρ w is density of water, d is depth, k is wavenumber, and G is gravitational constant. However, this form of equation applies only to Earth. We need to modify this equation for Moon as follows: where ρ c is the lunar crustal density. When admittance is expressed as an equation, it is the value obtained by dividing the cross-power spectrum of gravity and topography by the power spectrum of topography. The equation of admittance spectrum is as follows: where l is the spherical harmonics degree, S hg (l) is the cross-power spectrum of gravity and topography, and S hh (l) is the power spectrum of topography.
To determine the localized admittance, we computed a localized multitaper power spectrum for gravity and topography [44][45][46]. The MoonTopo2600p SH model was used to express the topographic data [47]. We set a spherical cap of a radius of 10 • for the seven regions to calculate a localized multitaper power spectrum. The localizing window bandwidth was set at 100, and we used 44 windows, with a concentration factor > 0.99. Figure 2 portrays the localized admittance for the seven regions. By combining Equations (3) and (4), we estimate the local crustal density and depth through least square fitting. Table 2 portrays the crustal densities calculated for the seven regions for the range of degree and order being 250-600, with only a small difference between the isostasy models in each adjustment.
bandwidth was set at 100, and we used 44 windows, with a concentration factor > 0.99. Figure 2 portrays the localized admittance for the seven regions. By combining Equations (3) and (4), we estimate the local crustal density and depth through least square fitting. Table 2 portrays the crustal densities calculated for the seven regions for the range of degree and order being 250-600, with only a small difference between the isostasy models in each adjustment. , and all other areas are located at mare or near mare. The admittances of Crookes and King increased due to the difference in topography (according to location). Marius Hills (yellow) was not a target area but an area for validating our method. The red color in the background was the range used to calculate the localized crustal density. Note. In depth, negative value means higher elevation rather than lunar radius.

New Processing Method of the Gravity Model
In this study, we used a gravity gradient tensor analysis method suitable for shortwavelength identification, using the Bouguer gravity derived from the crustal density of each region [41,48]. In addition, it is necessary to apply degree and order filtering to the gravity model to detect gravity changes at shallow depths. We used a high-pass filter, which remove the degree and order values of 1-60 to exclude the regional trend [41]. The gravity was calculated by increasing the range by a regular interval, from the minimum degree and order. Because the filtering range increased, the calculated results incorporated the shallow depth information. If each result was combined, such as a structure of spectral image (e.g., visible or infrared), the changes in gravity can be observed according to the depth of any region. We calculate the horizontal eigenvalues for effective analysis of residual trends [41]. The eigenvalue of gravity gradient tensor indicated negative , and all other areas are located at mare or near mare. The admittances of Crookes and King increased due to the difference in topography (according to location). Marius Hills (yellow) was not a target area but an area for validating our method. The red color in the background was the range used to calculate the localized crustal density. Note. In depth, negative value means higher elevation rather than lunar radius.

New Processing Method of the Gravity Model
In this study, we used a gravity gradient tensor analysis method suitable for shortwavelength identification, using the Bouguer gravity derived from the crustal density of each region [41,48]. In addition, it is necessary to apply degree and order filtering to the gravity model to detect gravity changes at shallow depths. We used a high-pass filter, which remove the degree and order values of 1-60 to exclude the regional trend [41]. The gravity was calculated by increasing the range by a regular interval, from the minimum degree and order. Because the filtering range increased, the calculated results incorporated the shallow depth information. If each result was combined, such as a structure of spectral image (e.g., visible or infrared), the changes in gravity can be observed according to the depth of any region. We calculate the horizontal eigenvalues for effective analysis of residual trends [41]. The eigenvalue of gravity gradient tensor indicated negative features for high density and positive features for low density. Calculating eigenvalues in this way indicated which feature change trend appeared in which degree and order, and the gravity changes at shallow depths were also identified efficiently. Figure 3 shows a simple conceptual diagram of this. We calculated the maximum degree and order starting from 65, increasing in intervals of 5. features for high density and positive features for low density. Calculating eigenvalu this way indicated which feature change trend appeared in which degree and order the gravity changes at shallow depths were also identified efficiently. Figure 3 sho simple conceptual diagram of this. We calculated the maximum degree and order sta from 65, increasing in intervals of 5. Concept of the method devised in our study. In the gravity model of spherical harm (SH), as the degree and order increase, a shallow short wavelength is expressed. Calculatio specifying a degree and order range can indicate the gravity at a specific depth level. It is ea find out what information is added at a shallower depth by increasing the degree and order at regular intervals. In arbitrary SH degrees -, if a mass exists in a specific degree and (correspond to depth), with a structure similar to that shown in (a), it can be expressed as (b) thr eigenvalues of gravity gradient tensor. In the range of degree and order corresponding to high sity mass, eigenvalues indicated a negative trend, and in the low-density mass section, eigenv indicated a positive trend.
In this study, we used gravity model GRGM1200B RM1 for λ = 1 (degree and o 1200). However, it is not possible to calculate the gravity field up to degree and order blindly to obtain shallow depth information. Because the gravity model express shorter wavelength with increasing degree and order, even though shallow depth i mation can be known, the effect of errors may increase. In this study, we determine maximum degree and order. In previous studies, the power spectrum of Bouguer gr revealed that as the degree and order increase, the power first decreases and the creases, which is a phenomenon that occurs when the error is larger than the main s [23,[29][30][31][32][33][34]. If there is a cave network beneath the impact melt pits on the Moon, inte tation at a high degree and order is required. Therefore, to exclude this error, we c lated the power spectrum of the localized Bouguer gravity for seven regions, and th gree and order of the point where the power increased was set to the maximum val be used for filtering; the values are shown in Table 3.  Concept of the method devised in our study. In the gravity model of spherical harmonics (SH), as the degree and order increase, a shallow short wavelength is expressed. Calculations by specifying a degree and order range can indicate the gravity at a specific depth level. It is easy to find out what information is added at a shallower depth by increasing the degree and order range at regular intervals. In arbitrary SH degrees n 1-4 , if a mass exists in a specific degree and order (correspond to depth), with a structure similar to that shown in (a), it can be expressed as (b) through eigenvalues of gravity gradient tensor. In the range of degree and order corresponding to high density mass, eigenvalues indicated a negative trend, and in the low-density mass section, eigenvalues indicated a positive trend.
In this study, we used gravity model GRGM1200B RM1 for λ = 1 (degree and order 1200). However, it is not possible to calculate the gravity field up to degree and order 1200 blindly to obtain shallow depth information. Because the gravity model expresses a shorter wavelength with increasing degree and order, even though shallow depth information can be known, the effect of errors may increase. In this study, we determined the maximum degree and order. In previous studies, the power spectrum of Bouguer gravity revealed that as the degree and order increase, the power first decreases and then increases, which is a phenomenon that occurs when the error is larger than the main signal [23,[29][30][31][32][33][34]. If there is a cave network beneath the impact melt pits on the Moon, interpretation at a high degree and order is required. Therefore, to exclude this error, we calculated the power spectrum of the localized Bouguer gravity for seven regions, and the degree and order of the point where the power increased was set to the maximum value to be used for filtering; the values are shown in Table 3.

Validation
We verified the proposed method by analyzing the MHH region, which has been significantly studied by Chappaz et al. [26]. The crustal density of the spherical cap centered on the MHH determined from the admittance at an angular radius of 10 • is 2569.72 kg/m 3 . To ignore the error that appears as the degree and order increase in SH, the gradient of the power of the localized Bouguer gravity was calculated, and the maximum degree and order applied to the analysis are 507, with a gradient exceeding 0. The eigenvalue map for the 60 to 510 section including the maximum order 507 is shown, and the positive feature can be seen along with the rille form including the pit (Figure 4a,b) is the eigenvalue variation for the two pit regions, and the range of 50 is set based on 510 including the maximum order 507, and the red dot is MHH. The reason for placing the maximum degree and order over the range of 50 was to identify the positive trend that indicated the mass deficit at the maximum degree and order, which portrayed the shallow depth information, excluding errors, within a narrow range.

Validation
We verified the proposed method by analyzing the MHH region, which has been significantly studied by Chappaz et al. [26]. The crustal density of the spherical cap centered on the MHH determined from the admittance at an angular radius of 10° is 2569.72 kg/m 3 . To ignore the error that appears as the degree and order increase in SH, the gradient of the power of the localized Bouguer gravity was calculated, and the maximum degree and order applied to the analysis are 507, with a gradient exceeding 0. The eigenvalue map for the 60 to 510 section including the maximum order 507 is shown, and the positive feature can be seen along with the rille form including the pit (Figure 4a,b) is the eigenvalue variation for the two pit regions, and the range of 50 is set based on 510 including the maximum order 507, and the red dot is MHH. The reason for placing the maximum degree and order over the range of 50 was to identify the positive trend that indicated the mass deficit at the maximum degree and order, which portrayed the shallow depth information, excluding errors, within a narrow range. The eigenvalue variation for the two pit points shows positive trends as the degree and order increase. Still, in the case of MHH (red dot), within the specified range, the trend does not increase but increases after 483 SH degrees. On the contrary, the increasing trend of the green point is more pronounced. From this result, it seems that the underground mass deficit of the green point is more apparent than that of MHH. Since it is not a cluster region of pits corresponding to our hypothesis, the cave network, there is a limit to the eigenvalue trend of MHH, strongly supporting our hypothesis.
Considering common sense, a mass deficit can naturally appear even in an area where there is no pit, so the verification result may be weak. However, we note that the cave network may contribute to several causes of mass deficit in the pit area. It is not possible to accurately distinguish caves from current data. Thus, in the best way, we verified the existence of the cave network in the pit area, which is one clue that can be accessed in the current situation. The eigenvalue variation for the two pit points shows positive trends as the degree and order increase. Still, in the case of MHH (red dot), within the specified range, the trend does not increase but increases after 483 SH degrees. On the contrary, the increasing trend of the green point is more pronounced. From this result, it seems that the underground mass deficit of the green point is more apparent than that of MHH. Since it is not a cluster region of pits corresponding to our hypothesis, the cave network, there is a limit to the eigenvalue trend of MHH, strongly supporting our hypothesis.
Considering common sense, a mass deficit can naturally appear even in an area where there is no pit, so the verification result may be weak. However, we note that the cave network may contribute to several causes of mass deficit in the pit area. It is not possible to accurately distinguish caves from current data. Thus, in the best way, we verified the existence of the cave network in the pit area, which is one clue that can be accessed in the current situation.

Results and Discussion
In the existing literature, no study has explored small-scale voids, such as cave networks in terrestrial subsurfaces, by analyzing gravity field data. Therefore, we derived the results by presenting a new method for analysis. Fortunately, this method has been verified to produce the same results, compared to the results of Chappaz et al. [26]. The magnitude of change in eigenvalues in the analysis results of the study area was small, and any errors may interfere with our interpretation. Therefore, we found the maximum degree and order, wherein the error would not be stronger than the gravity, through power spectrum analysis, and analyzed the change in gravity caused by the mass in the lunar subsurface, and thus, excluding the error as much as possible.
Then, we analyzed the variations in the eigenvalues over a specific range from the maximum degree and order for the seven regions. When interpreting the results, our main focus was to identify the variation trend of the eigenvalues, according to the degree and order, and not the eigenvalues distribution in the impact melt pit region. The color of the eigenvalues map in each region applied different stretches for improving the detection of eigenvalue variations and displayed the range of color bars.
As a result of the analysis, there are only four areas (Copernicus, King, Stevinus, and Tycho) where the increasing trend of the eigenvalue appears. The number of pits with increasing eigenvalue trends in each region is five in Copernicus, seventeen in King, four in Stevinus, and only one in Tycho. In the eigenvalue map, due to spatial resolution (2.2 km/pixel), a pit cluster is included in one pixel, so Copernicus and King are represented in three places, and Stevinus and Tycho are equal to the number of pits. In addition, the eigenvalues trend in the other pits area is not consistent with our hypothesis or it is not easy to interpret clearly in both the primary model and the clone model. In that example, the eigenvalues trend of all the pits in the King area ( Figure A1) where the pits are most concentrated is shown in Appendix A.
For the GRGM1200B RM1 for λ = 1 model, there are a total of 100 clone models. Although we applied to filter to exclude high-order errors from the spherical harmonic function, we randomly selected five replicated models for more detailed error analysis and applied the same research method to check whether the same results are obtained. As a result of the clone model analysis, an increasing pattern of eigenvalues appears, as in the result of the main model, and it is confirmed that the error of the gravity model does not contribute to the SH range of our study. This shows that our results are not affected by errors and are reliable.
The clone model results ( Figure A2) are shown in the Appendix A. Additionally, the clone analysis results of the Marius Hills region ( Figure A3) are also attached to the Appendix A.
Copernicus is a complex crater on the near side of the Moon and is 96 km in diameter. The crater is located on the wall of the buried crater that has a different topography compared to the surface [49][50][51]. In Copernicus, most of the pits are clustered in the northern part of the crater plain, so a mass deficit was expected. Contrary to expectations, however, the analysis showed mass deficits in three regions that were not related to the main cluster (Figure 5a,b). All of these pits are in the red region corresponding to the buried crater wall, and in particular, four pits are gathered in the area with green and blue points, suggesting the possibility of a cave network.
Considering the presence of olivine in the central peak of Copernicus [52], we were able to deduce that the mantle material was ejected through both the mare and highland crusts. If so, the mantle should also be uplifted, and the eigenvalues of the region should show a negative trend. However, in Figure 5a, Copernicus has a positive feature (red color) for most of it. It is reasonable to view this as the influence of buried craters. The buried crater center has a crustal thickness of 15 to 20 km and the Copernicus center has a crustal thickness of 30 to 40 km [42]. Regarding crustal thickness, we can say there is definitely mantle uplift in buried craters, but not in Copernicus. However, since the olivine of the Copernicus central peak exists, the mantle uplift can be expected, but the trace seems to be hidden due to the buried crater wall. Considering the presence of olivine in the central peak of Copernicus [52], we were able to deduce that the mantle material was ejected through both the mare and highland crusts. If so, the mantle should also be uplifted, and the eigenvalues of the region should show a negative trend. However, in Figure 5a, Copernicus has a positive feature (red color) for most of it. It is reasonable to view this as the influence of buried craters. The buried crater center has a crustal thickness of 15 to 20 km and the Copernicus center has a crustal thickness of 30 to 40 km [42]. Regarding crustal thickness, we can say there is definitely mantle uplift in buried craters, but not in Copernicus. However, since the olivine of the Copernicus central peak exists, the mantle uplift can be expected, but the trace seems to be hidden due to the buried crater wall.
The King crater is located on the far side and is an impact melt deposit region; this includes a satellite crater called King Y located north of the main crater, with a diameter of 20 km. The impact melt pits in this crater are mainly distributed in the central and eastern parts of the impact melt deposit. There is a total of 62 pits in King, and a positive trend appears only in three regions (Figure 6a,b). Three (red), ten (green), and four (blue) pits are clustered in each area. Among them, it can be said that there is a high possibility of the existence of a cave network in the green point area. The King crater is located on the far side and is an impact melt deposit region; this includes a satellite crater called King Y located north of the main crater, with a diameter of 20 km. The impact melt pits in this crater are mainly distributed in the central and eastern parts of the impact melt deposit. There is a total of 62 pits in King, and a positive trend appears only in three regions (Figure 6a,b). Three (red), ten (green), and four (blue) pits are clustered in each area. Among them, it can be said that there is a high possibility of the existence of a cave network in the green point area. Unlike the other target areas, King Y has a small area of impact melt deposit and the largest number of impact melt pits per area [14]. Additionally, the distribution of fractures and the impact melt pits overlap [53]. Considering the results of gravity gradient analysis in this region and the distribution of a large number of impact melt pits in a small area, the possibility of a potential lunar cave network seems to be most likely. Unlike the other target areas, King Y has a small area of impact melt deposit and the largest number of impact melt pits per area [14]. Additionally, the distribution of fractures and the impact melt pits overlap [53]. Considering the results of gravity gradient analysis in this region and the distribution of a large number of impact melt pits in a small area, the possibility of a potential lunar cave network seems to be most likely.
The Stevinus crater has a diameter of 75 km. There was a mass deficit in only 4 out of 26 pits. Similar to Copernicus, most of the pits are clustered in the eastern part of the crater, but a positive trend was observed in the location unrelated to the cluster (Figure 7a,b). Additionally, unlike Copernicus and King, the three pits (red, green, and blue) are located hundreds of meters apart, but it is difficult to see these pits clustering. Unlike the other target areas, King Y has a small area of impact melt deposit and the largest number of impact melt pits per area [14]. Additionally, the distribution of fractures and the impact melt pits overlap [53]. Considering the results of gravity gradient analysis in this region and the distribution of a large number of impact melt pits in a small area, the possibility of a potential lunar cave network seems to be most likely.
The Stevinus crater has a diameter of 75 km. There was a mass deficit in only 4 out of 26 pits. Similar to Copernicus, most of the pits are clustered in the eastern part of the crater, but a positive trend was observed in the location unrelated to the cluster ( Figure  7a,b). Additionally, unlike Copernicus and King, the three pits (red, green, and blue) are located hundreds of meters apart, but it is difficult to see these pits clustering.

Conclusions
In this study, we analyzed the eigenvalues of gravity gradient tensor changes for seven lunar craters to determine the possibility that an empty space, such as a lunar cave, may exist at shallow depths below the surface. We hypothesize the possibility of a network of small-scale caves, rather than large-scale ones. We selected an impact melt pit as a candidate and analyzed these pits to understand and determine gravity changes related

Conclusions
In this study, we analyzed the eigenvalues of gravity gradient tensor changes for seven lunar craters to determine the possibility that an empty space, such as a lunar cave, may exist at shallow depths below the surface. We hypothesize the possibility of a network of small-scale caves, rather than large-scale ones. We selected an impact melt pit as a candidate and analyzed these pits to understand and determine gravity changes related to depth. The existing literature does not explore small-scale subsurface voids using GRAIL's gravity model; therefore, we devised a new method that could support such explorations. For the validation of the method, we compared our results to those of existing studies that analyzed known areas, such as the MHH.
From the analyses of the seven regions using our newly devised method, a mass deficit was detected beneath the pits in the four regions of the Copernicus, King, Stevinus, and Tycho craters. This study revealed the possibility of the existence of a cave network in shallow depth of the crust. Moreover, the possibility of the existence of a cave network is one of several reasons that can contribute to the mass deficit detected in the region.
We obtained meaningful results using the latest data on human space exploration conducted till date; however, our argument may seem weak due to the limitation of the absolute resolution of the images used in our study. However, because in situ exploration is essential for the reliable exploration of lunar caves, remote sensing can play a crucial role for similar future lunar explorations. In addition, when observing gravity and seismic waves in situ, more detailed internal information can be obtained than remote sensing, so it is necessary to actively consider it as a science payload for a landing mission.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
As mentioned in Section 3, It shows an example in which it is not easy to interpret the eigenvalues' increasing trend. There are cases where the oscillation is large, both an increasing and a decreasing pattern appear, and a case where there is little change in the trend. Additionally, in the case of the clone model, it is difficult to interpret because it shows a similar pattern to Figure A1.
We describe the analysis results of 5 clone models randomly selected out of 100 for each study area ( Figure A2). In Clones 3 and 5 of Stevinus, green and blue are not clearly increasing. However, it can be argued that a positive trend appears for almost all clones. the eigenvalues' increasing trend. There are cases where the oscillation is large, both an increasing and a decreasing pattern appear, and a case where there is little change in the trend. Additionally, in the case of the clone model, it is difficult to interpret because it shows a similar pattern to Figure A1.
We describe the analysis results of 5 clone models randomly selected out of 100 for each study area ( Figure A2). In Clones 3 and 5 of Stevinus, green and blue are not clearly increasing. However, it can be argued that a positive trend appears for almost all clones.   Figure A2. Clone model analysis results for the four regions in Section 3. A total of 5 were randomly selected from among 100 clone models. There are some points where the increasing trend is weak in Stevinus' clone 3 and clone 5, but it can be seen that a positive trend appears in almost all of the clone models analyzed. Each plot is normalized as shown in Figures 5-7. Figure A2. Clone model analysis results for the four regions in Section 3. A total of 5 were randomly selected from among 100 clone models. There are some points where the increasing trend is weak in Stevinus' clone 3 and clone 5, but it can be seen that a positive trend appears in almost all of the clone models analyzed. Each plot is normalized as shown in Figures 5-7. Figure A2. Clone model analysis results for the four regions in Section 3. A total of 5 were randomly selected from among 100 clone models. There are some points where the increasing trend is weak in Stevinus' clone 3 and clone 5, but it can be seen that a positive trend appears in almost all of the clone models analyzed. Each plot is normalized as shown in Figures 5-7.