Effects of vertical nonlinearity on the superconducting gravimeter CT #036 at Ishigakijima, Japan

One of the characteristic features of the gravity recordings produced by the superconducting gravimeter CT #036 at Ishigakijima, Japan, is that it indicates gravity increase when a typhoon (hurricane) approaches the island. Since we are trying to detect small gravity signals associated with the long-term slow slip events in this region, it is very important in the interpretation of the observed data whether such gravity changes are of natural or instrumental origin. In this paper, we investigate whether or not nonlinearity in the sensor of the superconducting gravimeter is responsible for this phenomenon. Here we take the same theoretical approach as taken by our previous study which investigated the effect of coupling between horizontal and vertical components of the gravity sensor in order to understand the noise caused by the movements of a nearby VLBI antenna. From theoretical and experimental approaches, we prove that the gravity increase observed by CT #036 at the times of high background noise level cannot be explained by instrumental effects, such as the nonlinearity in the vertical component or the coupling between horizontal and vertical components of the gravity sensor. This implies that the observed gravity increases are real gravity signals of natural origin.


Introduction
Understanding the instrumental properties of a superconducting gravimeter (SG) is essential for fully making use of the high-quality gravity data it produces. The SG (CT #036) installed at Ishigakijima (Fig. 1), Japan, indicates several kinds of unusual gravity signals which are not seen at other SG sites, including Matsushiro, Kamioka, and Mizusawa. One of such signals are the systematic step changes caused by the motion of a nearby VLBI antenna (Honma et al. 2000). For any strange behavior of the instrument, there must be a physical cause that accounts for it. Imanishi et al. (2018) investigated this phenomenon from the viewpoint of the static and dynamic properties of the gravity sensor based on magnetic suspension. In the gravity sensor of the SG, there are two main superconducting coils, and the sphere is stably levitated near the plane of the upper coil (Fig. 2a). The levitation force is generated due to the interaction between the magnetic field and the induced currents in the superconducting sphere (Goodkind 1999). At the same time, the sphere is constrained to the central axis by the magnetic field, but there is room for movement in the horizontal plane. As a result of detailed analysis of data from a collocated seismometer, Imanishi et al. (2018) succeeded in giving a quantitative explanation to the observed gravity steps in terms of the coupling between the horizontal and vertical components of the gravity sensor. It was also shown by theoretical analysis that the mechanical eigenfrequency for horizontal translation of the levitating sphere is approximately 3 Hz in the case of CT #036.
Another characteristic feature of the CT #036 at Ishigakijima is that it indicates gravity increase when the background noise level is extremely high in stormy weather conditions. Due to its geographical location ( Fig. 1), the Ishigakijima island is often hit by typhoons (hurricanes). Figure  July through October 2012, as well as atmospheric pressure and the noise level recorded by the SG. The effect of atmospheric pressure on gravity has been corrected by the usual method with a single admittance. The noise level has been calculated by integrating the power of gravity signals in the frequency range 0.15-0.25 Hz. Note that this is significantly underestimated because of the attenuations by the analog anti-aliasing filter (GGP1 filter) in this frequency range. The events of low atmospheric pressure seen in this interval correspond to typhoons or similar meteorological events. We can see that when the atmospheric pressure becomes very low, the noise level becomes high, certainly because of stronger winds and higher oceanic waves. In such events, the gravity (Fig. 3b) also appears to be larger, with up to a few µGal excursion from the long-term trend (1 µGal = 10 -8 ms −2 ). For an ideal gravimeter, the DC output of the instrument should be independent of the background noise level at these frequencies. Because we are interested in detecting possible gravity signals associated with the long-term slow slip events occurring near Ishigakijima (Heki and Kataoka 2008), whether such gravity changes are real signals or not is extremely important in interpretation of the gravity recordings. At first, we suspected that this was another phenomenon to be explained with the coupling between the horizontal and vertical components of the gravity sensor as in the case of the VLBI antenna studied by Imanishi et al. (2018). However, it has turned out that this effect solely is not sufficient to provide the reason for the observed gravity changes, as discussed later in this paper. The purpose of this paper is to examine the effect of nonlinearity in the vertical component that may be existent in the gravity sensor of the SG in relation to the phenomenon mentioned above. Given the geometry Fig. 2 The gravity sensor of the superconducting gravimeter. a Superconducting sphere and coils. b Vertical force applied to the levitating sphere as a function of vertical position. The vertical force gradient is made weak for a limited interval of the vertical coordinate (z). Ideally, the force should be linearly dependent on z. In reality, there exists a small but finite quadratic term, which gives rise to nonlinear responses of the gravimeter Apparent correlation between noise level and gravity. a Atmospheric pressure. b Gravity residual. c Power of background noise in the frequency range 0.15-0.25 Hz calculated from gravimeter records. Notice that gravity changes are partly similar to the noise level, especially when the pressure is low inside the sensor and the currents in the superconducting coils, the total upward force depends on the position of the sphere as drawn in Fig. 3 of Goodkind (1999). By adjusting the ratio of the currents in the two coils, the vertical force gradient slightly above the plane of the upper coil can be made very weak so that the test mass levitated there is highly sensitive to small changes of gravity. Ideally, the levitation force around the position of balance should be linearly dependent on the vertical displacement, guaranteeing linearity of the gravity sensor as an acceleration transducer. Figure 2b shows schematically the result of our own computation of the levitation force by a finite element method based on the actual parameters for CT #036 (Imanishi and Takamori 2021). The overall curve of the levitation force is found to be very well approximated by a thirdorder polynomial function of the vertical coordinate z . For a particular value of z , the second derivative of the force equals zero, meaning that the force gradient is ideally constant. Below (above) that point, the graph of the function is convex downward (upward), and the coefficient of the quadratic term is negative (positive), as shown by a blue (red) curve. Existence of such quadratic dependence will give rise to a finite shift in the DC component in response to sinusoidal signal input, in other words, an apparent gravity change. Although this effect will be suppressed at lower frequencies thanks to the feedback control of the gravimeter, it might show up in the higher frequency range (above 0.1 Hz) where the feedback control is insufficient. Understanding the response of the SG to high-frequency disturbances will be necessary also for applying SG data to detection of earthquake-induced prompt gravity signals (Vallée et al. 2017;Kimura et al. 2019).
In the following sections, we will present the theory and measurement on this vertical nonlinearity. As a result of an experiment of artificial signal injection, a small but significantly positive estimate of the coefficient of the quadratic term was obtained. Based on this result and the data from a collocated seismometer, it will be shown that the effect of nonlinearity, if any, is negligibly small in the case of CT #036 at Ishigakijima. This implies that the observed gravity changes are real signals of natural origin.

Theory
To investigate the vertical nonlinearity of the SG, here we take the same theoretical approach as adopted by Imanishi et al. (2018). We shall begin with a brief review of the theory on the magnetic suspension in the SG. In a Cartesian coordinate system, the potential U sensed by the superconducting sphere (mass m ) can be expanded up to the third order as where α H , α V , β H , and β V are the coefficients of Taylor expansion series of the potential. The z is upward positive, and the origin of the coordinate system coincides with the mean position of the sphere. The coefficients of the second-order terms, α H and α V , denote the "spring constants" in the horizontal and vertical directions, respectively. The coefficients of the third-order terms, β H and β V , denote the deviation from a purely harmonic potential. Equations of motion for the sphere subject to external and frictional forces as well as the forces due to the magnetic field are written as where F x , F y , and F z are the components of the external force exerted on the sphere as viewed from the frame fixed to the gravity sensor. ω H and ω V are the angular eigenfrequencies of the sphere in the horizontal and vertical directions, respectively, given by and 2mh H ω H and 2mh V ω V are the coefficients of the friction for the horizontal and vertical components, respectively. Note that in the equations of motion (2)-(4) the terms containing α H and α V are first order of the coordinates, whereas the terms containing β H and β V are second order. By defining and assuming Equations (2) through (4) are rewritten as Here we have retained the terms up to the first order in the horizontal components and the terms up to the second order in the vertical component. Equations (12)-(14) are the basic equations to be studied in our analysis of the gravity sensor. The fifth term (containing β ′ V ) in the left-hand side of Eq. (14) comes from the third-order term of z in the potential. Imanishi et al. (2018) neglected this term in analyzing the effect of horizontal acceleration on the gravity sensor, because in the absence of input vertical acceleration z 2 can be regarded as negligibly small compared with x 2 , y 2 , and z . In the present study, we are interested in the effect of large vertical acceleration applied to the instrument, and therefore, possible nonlinear response of the gravimeter arising from the term containing β ′ V in Eq. (14) is the main subject to be studied. In the following, we derive the formulae that will be used as theoretical basis for experimentally measuring the coefficient β ′ V by applying artificial vertical acceleration to the superconducting sphere in the gravity sensor.
Let us consider the case where the applied vertical acceleration is a sinusoidal function of time with a single angular frequency ω . In a complex notation, the right-hand side of Eq. (14), the upward acceleration exerted on the sphere, can be written as where µ and ν are the real constants in the unit of acceleration denoting the cosine and sine parts of the complex magnitude of the input acceleration, respectively. Taking the real part of the right-hand side of Eq. (15), Eq. (14) can be written as where we have neglected the term containing β ′ H . Our task here is to find a solution of the differential Eq. (16).
Before we seek for a general solution to Eq. (16), let us consider the special case where β ′ V = 0 . In this case, Eq. (16) reduces to a linear differential equation which can be solved exactly. We assume a solution of the form: where a 1 and b 1 are real constants in the unit of displacement. Substituting Eq. (17) to Eq. (16), we obtain a set of equations: or in a matrix form as

Solving this, we have
This solution for z is the well-known linear response of a damped harmonic oscillator to an external sinusoidal force with a single frequency. The squared sum of a 1 and b 1 is given by We go on to the general case where β ′ V is finite. Because there is a z 2 term in the left-hand side of Eq. (16), it is no longer a linear differential equation, and z as the response to the external force would involve higher harmonics of ω . Let us assume an approximate solution of the form: where a 0 , a 1 , b 1 , a 2 , and b 2 are real constants in the unit of displacement. Substituting Eq. (22) into Eq. (16), and equating the coefficients for 1 (constant term), cos ωt , sin ωt , cos 2ωt , and sin 2ωt for the left-hand and righthand sides, we have the following five equations: (20) Regarding µ and ν as first-order quantities, Eqs. (24) and (25) indicate that a 1 and b 1 are first-order quantities, whereas Eqs. (23), (26), and (27) indicate that a 0 , a 2 , and b 2 are second-order quantities. Therefore, retaining up to first-order quantities, Eqs. (24) and (25) Equations (28), (29), and (31) give the approximate solutions for the coefficients a 0 , a 1 , b 1 , a 2 , and b 2 in equation (22). It is readily seen from Eqs. (29) and (31) that a 0 , a 2 , and b 2 vanish when β ′ V = 0 . For small but finite β ′ V , the response of the system slightly deviates from the purely harmonic solution, with the DC offset given by Eq. (29) and the second harmonics given by Eq. (31).
From Eqs. (21) and (33), the amplitude ratio between the ω and 2ω terms is given by This quantity is larger for lower frequencies. This means that the unknown parameters for the 2ω terms will be better determined using lower frequencies when we measure the response of the gravity sensor to injected artificial signals.

Measurements
The mechanical response of the gravimeter can be measured by applying artificial force to the superconducting sphere (Imanishi et al. 1996;GWR Instruments 1997;Van Camp et al. 2000). The method of measurement is very simple; one applies an external signal (in voltage) V in through the feedback circuit of the gravimeter with the normal feedback control disabled, and records a resultant signal V out in the error output of the gravimeter (the channel called Gravity Balance) as well as the input signal V in . Depending on the purpose of the measurement, the input signal V in can be an arbitrary function of time, including a step or a sinusoid.
Our experiment of measuring the instrumental response of CT #036 took place on June 7, 2018. We used an Agilent 33210A function generator to generate artificial signals and a Hakusan DATAMARK LS-8800 data logger to record both input and output signals. The sampling frequency was 200 Hz. Sinusoidal functions with 0.5 V amplitude at eight discrete frequencies from 0.001 to 0.2 Hz were applied (Table 1). Figure 4 shows the recorded time series of the input and output signals. For each frequency, more than six cycles of oscillations were recorded. As shown in Fig. 4b the output signal is affected by natural gravity changes, mostly the earth tides. We noticed that the power-linecycle (60 Hz) noise became large in the output signal from the Gravity Balance when it exceeded ± 3 V. The cause of this noise is unknown. The output signal from the Gravity Balance channel is filtered by a built-in analog lowpass filter with its corner frequency of 0.2 Hz. In order to extract the response of the gravimeter, amplitude and phase of the wave elements recorded on the Gravity Balance channel must be corrected for the response of this filter. We made an experiment to calibrate the response of the Gravity Balance lowpass filter for the CT #036. The calibrated response was found to be slightly different from that theoretically calculated with the nominal constants given in the circuit diagram (GWR Instruments, 1997). In the analysis below, we use the calibrated response to correct for the Gravity Balance filter. See Additional file 1 for more details on the calibration of the filter response.
We will divide the analysis of the measurement results into two stages. In the first stage, we adopt the linear model, corresponding to Eqs. (17,18,19,20,21), to interpret the response of the gravimeter. The purpose of this treatment is to obtain an internally consistent set of parameters for CT #036, which will be needed to combine measured voltage with physical quantities. Then, in the second stage, we adopt the model in which nonlinearity is taken into account, corresponding to Eqs. (22,23,24,25,26,27,28,29,30,31,32,33).

Case 1: linear response
In this subsection, we treat the linear case ( β ′ V = 0 ). In this case, the input and output time-domain signals, V in and V out , are linearly connected through an impulse response. In the frequency domain, this is described as where Ṽ in and Ṽ out are the Fourier transforms of V in and V out , respectively. φ 1 , φ 2 , and φ 3 are three stages of transfer functions as explained below. First, φ 1 , in the unit of ms −2 V −1 , translates voltage into acceleration. A signal (in voltage) of unit magnitude input to the feedback circuit of the gravimeter generates an acceleration applied to the superconducting sphere, whose magnitude is equal to φ 1 . Therefore, φ 1 stands for the DC sensitivity of the gravimeter, and is often called a scale factor. Calibration of the scale factor is usually made through parallel registration with an absolute gravimeter (e.g., Imanishi et al., 2002;Crossley et al., 2018). Absolute gravity measurements at Ishigakijima was performed in January 2015 (Miyakawa et al., 2020) for the purpose of calibrating instrumental drift as well as the sensitivity of CT #036. The resultant scale factor (for the normal 100 kΩ feedback resistor) was which falls within the typical range (i.e., (50 − 100) × 10 −8 ms −2 V −1 ) of sensitivity for an SG. In this study, we treat φ 1 as a known quantity.
Next, φ 2 is the mechanical response of the superconducting sphere against applied acceleration of unit magnitude. Neglecting higher-order terms in Eq. (14), φ 2 as the linear response of a damped harmonic oscillator is derived for angular frequency ω as This transfer function translates acceleration into a displacement of the sphere. The unit of φ 2 is s 2 . Note that among the three transfer functions φ 2 is the only frequency-dependent part. Note also that when ω = 0 , φ 2 is equal to 1/ω 2 V . Finally, φ 3 converts a displacement of the superconducting sphere of unit magnitude into voltage through the position detector. This has a unit of Vm −1 . In this paper, we do not take possible nonlinearity of the position detector into account, and assume that φ 3 is invariant within the range of displacement of the sphere in our experiment. As far as we are aware, this factor has gathered much less attention than the scale factor of the gravimeter among the users of SG. One of the purpose here is to calibrate φ 3 which will be used later to convert measured voltage into displacement of the sphere.
Putting together the three stages, the total transfer function is given by is a new variable in the unit of s −2 . ω V , η , and γ in the right-hand side of Eq. (38) are the parameters to be measured by our experiment.
For each frequency f , the input signal V in is fit to the model function: where ω = 2π f , and A in 0 , A in 1 , and B in 1 are the free parameters (in the unit of V) to be adjusted. If multiplied by φ 1 , the right-hand side of Eq. (40) corresponds to the righthand side of Eq. (16). In other words, The term A in 0 in the right-hand side of Eq. (40) is included to account for a small offset of the DC component that may exist in actual experiments.
Similarly, the output signal V out is fit to the function: where A out 0 , A out 1 , and B out 1 are the free parameters (in the unit of V) to be adjusted. If divided by φ 3 , the right-hand side of Eq. (42) corresponds to the right-hand side of Eq. (17), in other words, A out 1 and B out 1 must be corrected for the response of the analog lowpass filter for the Gravity Balance channel at the angular frequency ω.
In a real experiment, special care must be taken in estimating A out 0 , A out 1 , and B out 1 , because the gravimeter is subject to the gravity changes of natural origin as well as the injected artificial signal during the experiment. Then, instead of Eq. (42), the formula to be used is where g calc (t) is the known gravity signals in the unit of acceleration. Here we take into account the theoretical tides and the effect of atmospheric pressure as known natural signals. We assume that these are long periodic enough, so that there is no need to take the frequencydependent response of the gravity sensor into account. We also assume that there are no other agents of natural gravity signals affecting the measurements. The parameter C in the right-hand side of Eq. (44) is the DC sensitivity (in other words, the "scale factor") of the gravimeter operated without feedback. It is given by Imanishi et al. Earth, Planets and Space (2022) 74:73 where we have used Eqs. (37) and (39). Because C is dependent on ω V and γ , determination of C as well as ω V , η , and γ must be done in an iterative way as follows. First, a provisional value of C , say C 0 , is assumed to correct for the known natural signals in Eq. (44). Once the unknown parameters in the right-hand sides of Eqs. (42) and (44) are adjusted, ω V , η , and γ in Eq. (38) can be estimated. This is actually done in terms of the cosine and sine parts of the both sides of Eq. (38) for angular frequency ω as and Then, we can calculate C based on Eq. (45), which must reproduce the initial value C 0 so that the solution is internally consistent. After some trials, we found that the estimates of ω V , η , and γ , and therefore of C , are highly invariant for a plausible range of the choice of C 0 . Figure 5 shows the estimated value of C for different initial values of C 0 . As a result, C = 27.65 × 10 −8 m s −2 V −1 was found to be the most appropriate value of C . Comparing this value with the usual scale factor given in Eq. (36), we can say that the gravimeter without feedback control is approximately twice as sensitive as that under feedback control in the case of CT #036.  Figure 6 shows the magnitude and phase of the transfer function given by Eq. (38). The dots are the values for discrete frequencies measured in our experiment, and the curves are the values theoretically calculated using the results of Eqs. (48) through (50). The measured values are very well approximated by the model function. Note that the magnetic suspension of the SG as a mechanical pendulum system is overdamped as seen from Fig. 6(a) (Imanishi et al., 1996).
It will be useful now to verify the physical magnitude of the applied acceleration and the resultant displacement of the sphere. We take the case of frequency = 0.001 Hz. For the input signal, A in 1 = −0.0268V and B in 1 = −0.4984V (note that A in 1 2 + B in 1 2 ∼ = (0.5) 2 ). Therefore, from Eqs. (36) and (41) (52) φ 3 = γ φ 1 = 2.5 × 10 6 Vm −1 .  (45)) by an iterative estimation process Imanishi et al. Earth, Planets and Space (2022) 74:73 Case 2: nonlinear response Now we move on to the more general case where the third-order term of the potential is treated as finite. In this case, the linear transfer function φ 2 , and therefore the relation (38), does not apply. The purpose here is to directly estimate the magnitude of the second harmonics of the sphere excited by a sinusoidal input. Knowledge of φ 3 is necessary to convert observed amplitude in voltage to physical quantities in terms of displacement.
Following Eq. (26), the function to be fit to the output signal V out is where A out 0 , A out 1 , B out 1 , A out 2 , and B out 2 are the free parameters to be adjusted. The parameter C is treated as known from the analysis of the linear case. Once these parameters are estimated, they are divided by the known parameter φ 3 to obtain a 0 , a 1 , b 1 , a 2 , and b 2 in the right-hand (53) side of Eq. (22). In particular, the coefficients for the second harmonics are given by The coefficients A out 2 and B out 2 must be corrected for the response of the analog lowpass filter for the Gravity Balance channel at the angular frequency 2ω.
It can be seen from Eqs. (29) and (31) that the parameters a 0 , a 2 , and b 2 (in other words, A out 0 , A out 2 , and B out 2 ) contain information on the desired coefficient β ′ V . However, estimating β ′ V based on the measured value of a 0 is not practical, because it may be seriously affected by possible DC offsets in the measurement of both input and output signals as well as contamination of tidal and other long-periodic geophysical signals. Instead, one can determine a 2 and b 2 in a more robust way, from which β ′ V can be estimated. Let a 2 ± a 2 and b 2 ± b 2 be the estimates thus obtained ( a 2 and b 2 are the uncertainties). The theoretical predictions of a 2 and b 2 , given by Eq. (31), can be rewritten as and where and Both A and B can be calculated from the parameters already obtained in the linear case. Therefore, we could estimate β ′ V in two ways by for the cosine part and for the sine part. Considering that the relative estimation errors of a 2 and b 2 are much larger than those of the other parameters, the estimation error of β ′ V , �β ′ V , for Eqs. (59) and (60) may be obtained by and respectively. However, estimation by Eqs. (59) and (60) can be unstable because a 2 and b 2 depends on the initial phase of the applied oscillations. So, a more robust estimate of β ′ V can be obtained by taking a weighted mean of the two results as with its uncertainty given by Table 2 summarizes the results of fitting for frequencies f = 0.001, 0.002, 0.005, and 0.01 Hz. Note that for each frequency the magnitude of the 2f terms ( A out 2 and B out 2 ) is smaller than those of the f terms ( A out 1 and B out 1 ) by three to four orders of magnitude. Here, AIC 1 and AIC 2 are the Akaike Information Criterion (AIC; Akaike 1972) for the models based on Eqs. (44) and (53), respectively. For frequencies f = 0.001 and 0.002 Hz, AIC 2 is smaller than AIC 1 , meaning that existence of the second harmonics in the measured data is statistically significant from the viewpoint of AIC. We shall exclude the frequencies f = 0.005 and 0.01 Hz from the following discussion, because AIC 2 is larger than AIC 1 for these frequencies.
For f = 0.001 and 0.002 Hz, we have obtained the estimates of β ′ V for the cosine and sine parts based on Eqs. (59) and (60), as listed in Table 2. Note that all the four estimates of β ′ V are positive. Also, for each frequency, the estimates of β ′ V from the cosine and sine parts are marginally consistent with each other within the estimation error. Using Eqs. (63) and (64), our final estimate is for f = 0.001 Hz and for f = 0.002 Hz. These results from the two frequencies are also consistent with each other within the estimation error. Thus, we have obtained statistically significant results indicating that β ′ V is positive and its magnitude is around 2 × 10 3 in the case of CT #036.

Discussion
Based on the estimate of the coefficient β ′ V , we derive the expected offset in the DC component in gravity signal for general vertical acceleration input. Let Z be the vertical displacement of the ground with respect to the inertial frame. In the equation of motion (14) in the vertical component, the right-hand side is the applied acceleration and may be identified as −Z . Therefore, the equation of motion becomes where we have neglected the term containing β ′ H as before. This differential equation is nonlinear with respect to z . Based on the potential given by Eq. (1), here we may assume that β ′ V is sufficiently small, so that a condition holds for a possible range of displacement of the sphere. Under this condition, we can expand the solution for z into a power series of β ′ V as This is substituted into Eq. (67). For the zeroth order of The solution to this, ζ 0 , is the usual linear response in the vertical direction, expressed in the frequency domain as where ζ 0 and Z are the Fourier transforms of ζ 0 and Z , respectively. For the first order of β ′ V , we have Since we are interested in the DC component of ζ 1 , we adopt approximations �ζ 1 � = 0 and �ζ 1 � = 0 , where denotes temporal average. Then, taking temporal averages of Eq. (72), we obtain This means that, to the first order of β ′ V , the apparent gravity change due to this effect is approximated as Equation (74) indicates that a gravity increase is predicted if β ′ V > 0 . As was done in Imanishi et al. (2018), we can estimate ζ 2 0 using the data of a 1.0-Hz velocity transducer which is collocated with the gravimeter (Ohtaki and Nawa, 2013). First, the raw data in the vertical component of the seismometer are converted into ground displacements by deconvolving the transfer function of the sensor. Second, the spectrum of vertical ground displacements is converted to that of sphere's vertical motion with the transfer function given by Eq. (71). Finally, the power spectrum is integrated with respect to frequency to give the total power, i.e., ζ 2 0 . Here we use the data acquired when the typhoon "Neoguri" (https:// en. wikip edia. org/ wiki/ Typho on_ Neogu ri_ (2014)) approached Ishigakijima on July 8th, 2014. The track of this typhoon is shown in Fig. 1. Figure 7 shows the data of atmospheric pressure, sealevel, and gravity during the six-day period from July 5th till c Gravity July 10th. The atmospheric pressure data were recorded at the SG site, whereas the sealevel data were taken at the tide gauge station at the Ishigaki port, southern part of the island. Figure 7a and b clearly shows that the sealevel rose by about 40 cm due to the drop of atmospheric pressure by about 27 hPa. The gravity data have been corrected for the earth and oceanic tides as well as the atmospheric effect. For the latter correction, a single admittance (-0.338 µGal/hPa) has been used. After these corrections, an event of gravity increase of approximately 2 µGal around July 8th becomes evident in the residual (Fig. 7c). We can conclude that this gravity change is not explained solely by the loading effect of the atmosphere or the ocean, for the following reasons. First, a simple calculation shows that a uniform rise in the sealevel by 1 cm around Ishigakijima (100 km × 100 km area) causes a gravity increase of 0.0015 µGal at the SG station by the Newtonian effect. Therefore, the maximal rise in the sealevel of 40 cm in that area should have caused a gravity increase of about 0.06 µGal, and slightly larger when the deflection of the ground is taken into account. However, this is much smaller than the observed gravity change.
Next, it will be worth considering the possibility that the effective atmospheric admittance was different from usual when the typhoon approached the island, because it can be variable with temporal and spatial scales of atmospheric disturbances (Hinderer et al. 2014). Given the magnitude of the drop of the atmospheric pressure as well as the gravity increase, an atmospheric admittance as large (in the absolute sense) as -0.43 µGal/hPa would be required so that the gravity change associated with the typhoon would almost vanish after atmospheric correction. This large admittance is close to that derived for a simplistic Bouguer plate model (Zürn and Wielandt 2007), and is unlikely to be the case, even if we consider that the typhoon Neoguri was a huge one spanning about a few hundreds of kilometers. The fact that the gravity change is lagged behind the changes in the atmospheric pressure and the sealevel by a few hours, as shown in Fig. 7, provides another evidence that the residual gravity change is not explained solely by the loading effect. Figure 8a shows the spectrum of vertical displacements of the sphere obtained by converting the seismometer data. Numerical integration of this spectrum gives �ζ 2 0 � = 1.82 × 10 −13 m 2 . Then, taking the estimate of β ′ V obtained for f = 0.001 Hz, we obtain from Eq. (74) g ∼ 1.8 × 10 −10 ms −2 , in other words, about 0.02 µGal. This is roughly two orders of magnitude smaller than the recorded gravity change. Therefore, we can conclude that the effect of the second-order term in the equation of motion (in other words, the third-order term in the potential) cannot account for the gravity increase observed in stormy weather.
To confirm our conclusion further, we examine additionally whether or not the coupling between horizontal and vertical components of the SG (Imanishi et al. 2018) is responsible for the observed gravity increase. In contrast to the noise from the VLBI antenna as investigated in Imanishi et al. (2018), the background noise has continuous spectra with a peak around 0.2 Hz. In order to estimate the gravity effect according to the model of Imanishi et al. (2018), we need to calculate the mean squared horizontal displacements of the sphere in the right-hand side of Eq. (36) of Imanishi et al. (2018). To do so, we use the transfer function (34) of Imanishi et al. (2018) to convert the spectrum of ground displacements into the spectrum of sphere's displacements, and then integrate the resultant spectrum with respect to frequency to obtain the estimate of total power. Figure 8b and c shows the spectra of the horizontal displacements of the sphere converted from the seismometer records using the parameters estimated in the previous section.
Because the horizontal eigenfrequency of the sphere is about 3 Hz, the spectral power at 0.2 Hz is diminished in this spectrum, compared with the power at 4-5 Hz. By integrating these spectra in the frequency range 0-20 Hz, we obtain �x 2 � = 6.11 × 10 −15 m 2 and �y 2 � = 3.68 × 10 −14 m 2 , where x and y correspond to east-west and north-south components, respectively. Substituting these values into Eq. (36) of Imanishi et al. (2018) and using the revised estimate of β ′ H (Additional file 2), an apparent gravity change would be This is negligibly small compared with the observed gravity change. Note that this is of the same order of magnitude as, but larger than, the effect of the higher-order term estimated above.
As discussed in Introduction, the sign and magnitude of β ′ V depend on the exact position of the levitating sphere as well as the currents in the superconducting coils (Imanishi and Takamori 2021). It should be noted that our results obtained here apply only to CT #036 at Ishigakijima, and various installations of other SGs may indicate different characteristics regarding the response to a high level of background noise.

Conclusion
We have investigated the gravity increase which the SG CT #036 at Ishigakijima indicates at the times of high background noise level. From theoretical and experimental analyses, we have proved that this phenomenon cannot be explained by instrumental effects, such as the nonlinearity in the vertical component or the coupling between horizontal and vertical components of the gravity sensor. So we have concluded that those events of gravity increase are real gravity signals of natural origin. Now that we have reached this conclusion, the next challenge will be to investigate the geophysical cause of the observed gravity signal. Possible origins of the signal include the oceanic loading effects on the shelf of the island (Fratepietro et al. 2006;Mouyen et al. 2017) and the effects of underground water. These subjects will be investigated in future.