Simulation of Patch Field Effect in Space-Borne Gravitational Wave Detection Missions

TL;DR

This paper introduces a new simulation method to study patch field effects in space-based gravitational wave detectors, improving sensor accuracy.

Contribution

A novel boundary element partitioning and octree-based simulation algorithm is proposed to model electrostatic and geometric impacts of patch fields.

Findings

ΔFx and ΔKxx show linear dependence on patch potential variation (Δu).

Patch effects can be modeled using a quartic polynomial dependent on patch radius (rₘ).

The new method efficiently models both electrostatic and geometric impacts with low computational complexity.

Abstract

Space-borne gravitational wave detection missions demand ultra-precise inertial sensors with acceleration noise below 3×10−15 m/s2/Hz. Patch field effects, arising from surface contaminants and nonuniform distribution of potential on the test mass (TM) and housing surfaces, pose critical challenges to sensor performance. Existing studies predominantly focus on nonuniform potential distributions while neglecting bulge effects (surface deformation caused by the adhesion of pollutants or oxides, production and processing defects, and other factors) and rely on commercial software with limited flexibility for customized simulations. This paper presents a novel boundary element partitioning and octree-based simulation algorithm to address these limitations, enabling efficient simulation of both electrostatic and geometric impacts of patch fields with low spatiotemporal complexity (O(n)). Leveraging this framework, we systematically investigate the influence of single patches on the TM electrostatic force (ΔFx) and stiffness (ΔKxx) through parametric studies. Key findings reveal that ΔFx and ΔKxx exhibit linear dependence on patch potential variation (Δu) and can be fitted by a quartic polynomial (which can be simplified in some cases, such as only a cubic term) about patch radius (rₘ). The proposed method’s capability to concurrently model geometric bulges and potential nonuniformity offers significant advantages over conventional approaches, providing critical insights for gravitational wave data analysis. These results establish a foundation for optimizing mitigation strategies against patch-induced noise in future space missions.

Figures

Full text

1. Introduction

The detection of gravitational waves has revolutionized astrophysics since the first direct observation of a binary black hole merger by LIGO in 2015 [1]. This breakthrough opened a new window for exploring the universe, particularly in the mid- and low-frequency bands, which are best observed from space using laser interferometry [2,3,4,5,6]. To explore these frequency bands, several space-based missions have been proposed, including LISA [7], Taiji [8], and Tianqin [9]. These missions aim to achieve unprecedented sensitivity, enabling the detection of gravitational waves from sources such as supermassive black hole mergers and compact binary systems. Inertial sensors are critical components of space gravitational wave detectors, featuring a sensitive structure composed of a test mass (TM) and a housing equipped with electrodes. These electrodes include injection electrodes and sensing/actuation electrodes for position measurement and feedback control [10]. Space gravitational wave detection demands extremely low acceleration noise levels (e.g., $[eqn]$ or better) [11,12]. To achieve this, it is essential to analyze and mitigate all potential noise sources, including those arising from patch fields of the inertial sensor.

Patch fields can arise on the surface of the TM or the inner surface of the housing due to factors such as surface contaminants, nonuniform material segregation, and variations in crystal orientation [13,14,15,16,17,18]. Patches lead to nonuniform potential distributions and localized bulges (surface deformation caused by the adhesion of pollutants or oxides, production and processing defects, and other factors), which can significantly affect sensor performance by introducing unwanted forces and stiffness changes. Numerous studies have been conducted on the patch field effect of inertial sensors in space gravitational wave detection missions [13,19,20,21,22,23,24]. These studies can be separated into two approaches: theoretical derivation and experimental investigations. The experimental approaches are further divided into physical experiments and computer simulations. Existing studies mainly focus on the effect of nonuniform potential distribution caused by patches, ignoring the effect of bulges (surface deformation). Additionally, many simulation-based studies rely on commercial software, which limits the flexibility (such as using GPU, convenient simulation of patches, and so on) to design customized algorithms for unique requirements.

In this paper, we propose a novel numerical simulation algorithm based on boundary element partitioning and octree structure. This method enables the efficient simulation of patch field effects under various conditions with reduced computational complexity. Using this approach, we investigate the impact of a single patch on the force and stiffness of the TM. The remainder of this paper is organized as follows: Section 2 presents the mathematical model and simulation method for patch field analysis, including the subdivision of surfaces and the calculation of electric field forces. Section 3 introduces a low-complexity algorithm based on charge equivalence and octree structure to improve computational efficiency. Section 4 details the experimental setup, results, and discussion, focusing on the effects of patch size, location, and potential distribution. Finally, Section 5 concludes the paper with a summary of key findings and potential future research directions.

2. Mathematical Modeling of Patch Field Simulation

2.1. Surface Subdivision and Electric Field Modeling

Figure 1 illustrates the schematic diagram of the sensitive structure. We define a rectangular coordinate system in space. With the housing as the reference frame, the geometric center of the housing is located at the origin, each face is perpendicular to a coordinate axis, and the $[eqn]$ -axis direction coincides with the sensitive axis direction. When the TM is in equilibrium, its geometric center is also at the origin, and each face of the TM is also perpendicular to a coordinate axis. The surfaces of both the TM and housing are discretized into multiple rectangular regions. Each rectangular region is further subdivided into smaller rectangles, which are then divided into four triangular elements, as illustrated in Figure 2. This fine-grained subdivision enables precise modeling of surface geometry and facilitates the simulation of localized features, such as bulges caused by contaminants. By leveraging the geometric flexibility of triangles, we can simulate bulges by displacing specific nodes, as demonstrated in Figure 3.

2.1.1. Potential and Charge Distribution

According to [25,26], the potential of the TM is

$$[eqn]$$

where $[eqn]$ represents the charge accumulated on the TM due to cosmic rays, $[eqn]$ is the capacitance between the TM and a portion of the housing, and $[eqn]$ is the potential of the corresponding electrode. If the portion is not an electrode, $[eqn]$ is set to zero, as the area can be approximated as ground.

The potential at the centroid of each triangular element is assigned based on its location. For instance, if one triangular element whose centroid lies on an electrode with a potential of 3 mV, the potential of the triangular element is set as 3 mV. The potential at the centroid of triangular element i due to triangular element j is given by:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ is the charge surface density, and $[eqn]$ is the vacuum dielectric constant. $[eqn]$ is the area of the triangular element j, and $[eqn]$ stands for the distance between point Q and the centroid of the triangular element i. The potential at the centroid of triangular element i is the sum of contributions from all triangular elements:

$$[eqn]$$

where n stands for the number of triangular elements. For each triangular element, its $[eqn]$ is unknown, while the potential at its centroid can be obtained directly. We can construct a system of linear equations e.g., Equation (4) to solve the charge distribution:

$$[eqn]$$

The matrix element $[eqn]$ is defined as:

$$[eqn]$$

where $[eqn]$ is the area of triangular element j, and $[eqn]$ is the distance between point $[eqn]$ on triangular element j and the centroid of triangular element i. Since the value of n can be very large, we cannot solve Equation (4) directly. Instead, we employ the Generalized Minimal Residual (GMRES) method, an iterative solver well-suited for large linear systems.

2.1.2. Calculation of Electric Field Force and Stiffness

In a moment, once the charge distribution is worked out, the electrostatic force and its stiffness acting on the test mass can be calculated. Specifically, we focus on $[eqn]$ , which stands for the electrostatic force acting on the TM along the X-axis direction and the electrostatic stiffness $[eqn]$ , which is defined as the rate of change of $[eqn]$ with respect to the displacement of the TM along the X-axis. These quantities can be expressed as:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ stands for the distance between points P and Q, and $[eqn]$ and $[eqn]$ represent the x coordinates of points P and Q, respectively. The number of triangular elements on the TM is $[eqn]$ , and of all the triangular elements (on all faces of the TM and the housing, a total of 12), the first $[eqn]$ are on the TM.

2.1.3. Algorithm Complexity

Storing the information of triangular elements requires a space complexity of $[eqn]$ . Due to the large value of n, the matrix in Equation (4) cannot be stored entirely. For example, if $[eqn]$ (the calculation accuracy is very low under such subdivision scale) and matrix elements are stored using 64-bit variables, the storage required to store the entire matrix will be about 7450.6 GB. Because of the extremely large storage requirement, matrix elements are computed on-the-fly when needed. The space complexity for calculating potential, force, and stiffness is $[eqn]$ .

When solving the equations iteratively, the time complexity required for one iteration of the solver is $[eqn]$ , and the same complexity applies to the calculation of electrostatic force and stiffness. The high time complexity limits the algorithm’s applicability, as it can only handle problems with relatively small values of n. Therefore, a low-complexity calculation method is urgently needed.

2.2. Simulation of a Patch

We first define the outer normal direction. At any point on a surface, the outer normal direction is perpendicular to the surface and points toward the gap between the housing and the TM. For the outer surface of the TM, the outer normal direction points to the outside of the TM, while for the inner surface of the housing, the outer normal direction points inward.

To simulate the single bulge, we introduce two parameters: $[eqn]$ (the base radius) and $[eqn]$ (the maximum displacement). A center point is selected, and all the nodes with a distance $[eqn]$ from the center are displaced along the outer normal direction by a distance d. The relationship between r and d is given by:

$$[eqn]$$

Equation (8) and Figure 4 illustrate the typical shape of the bulge, which resembles a cone with a base radius of $[eqn]$ and a height of $[eqn]$ . Specifically, if $[eqn]$ , the patch affects only the potential distribution without altering the surface geometry.

All the triangular elements that have at least one vertex belonging to the nodes of the patch are considered part of the patch. Each triangular element of the patch involves the same change in the amount of the potential of the centroid of the triangular element. The influence of the parameters, such as $[eqn]$ and $[eqn]$ , on the force and stiffness is investigated in detail in Section 4.

3. Low-Complexity Calculation Method

As mentioned in Section 2.1.3, we need to find an efficient algorithm to reduce computational complexity for large-scale simulations. Ying, L. has proposed a fast multipole method for particle simulations [27]. Many other improved algorithms [28,29,30] based on this method have also been developed. These algorithms inspired us. Aiming at the problems studied in this paper, we have made modifications and optimizations on the basis of these algorithms. This section will elaborate on the modified method and its application in our research.

3.1. The Theoretical Basis of Charge Equivalence

For the Laplace equation $[eqn]$ , which governs potential distributions in singularity-free regions, the first kind of boundary conditions can be used as a complete definite solution condition [31]. For an area without any singularity inside (or outside) it, if two groups of charges outside (or inside) the region produce the same potential at any point on the border of the area, the effects of the two groups on the inside (or outside) of the area are equivalent. See Figure 5.

To discretize the problem, we employ m sampling points distributed along the region’s boundary. Let the charge distribution outside (or inside) the region be represented by a charge vector $[eqn]$ . The resultant electric potential at these sampling points can be formulated as $[eqn]$ , where $[eqn]$ denotes the potential coefficient matrix mapping charges to potentials. To reconstruct the electric field inside (or outside) the region, we introduce an equivalent charge vector $[eqn]$ comprising k charges and a corresponding potential mapping matrix $[eqn]$ , where $[eqn]$ provides the potential vector at the sampling points. Then we can obtain

$$[eqn]$$

If E is invertible, then

$$[eqn]$$

otherwise

$$[eqn]$$

where $[eqn]$ is the generalized inverse of E. The process of multiplying T and S is called the first stage of the equivalency, and the process of multiplying $[eqn]$ (or $[eqn]$ ) and TS is called the second stage of the equivalency.

3.2. Building an Octree

For a three-dimensional region $[eqn]$ , we can divide it into eight small regions listed in Table 1, where $[eqn]$ , $[eqn]$ , and $[eqn]$ . After executing this process recursively, we can obtain an octree (see Figure 6). The recursive exit is when the volume of the current region is less than the threshold.

A pair of nodes are neighbors if and only if they are in the same layer of the octree and have common vertices, including the case where two nodes are the same node.

A pair of nodes are interactive nodes of the same layer if and only if their parent nodes are neighbors but they are not neighbors themselves.

We chose a cube with edges of length L whose centroid is at (0, 0, 0), which can surround the whole sensitive structure. We constructed an octree recursively with this cube as the root node. To ensure each octree node contains at least one triangular element, the octree should be pruned. If one region is inside the TM, or outside the housing, or in the gap between the TM and the housing, we will not continue to divide it, and the node will not be considered, such as in Figure 7.

For each node of the octree, we built two auxiliary surfaces, scaling up with the centroid of the node as the center. According to [27], the proportional coefficients can be determined as $[eqn]$ and $[eqn]$ , respectively. We refer to the inner (outer) one of the two surfaces as the inner (outer) auxiliary surface.

When analyzing the influence on the area outside the outer auxiliary surface of a node made by the charges inside the node, we should reconstruct the electric field outside the outer auxiliary surface by placing charges (equivalent point charges, defined in Section 3.1) on the inner auxiliary surface. The sampling points (defined in Section 3.1) are on the outer auxiliary surface.

When analyzing the influence on the area inside the inner auxiliary surface of a node made by part of the charges outside the outer auxiliary surface, we should reconstruct the electric field inside the inner auxiliary surface by placing charges on the outer auxiliary surface. The sampling points are on the inner auxiliary surface.

The placement of equivalent point charges and sampling points on the two kinds of auxiliary surfaces should be regular to facilitate the design of the algorithm. Each auxiliary surface is made up of six faces. We use the parameter p to describe the number of sampling points or charges. For one of the six faces, the number of sampling points or charges is $[eqn]$ . Each face is divided into $[eqn]$ regions by an imaginary grid, and the sampling points or charges are distributed on the grid nodes. For example, Figure 8 shows the distribution on a single face when $[eqn]$ .

3.3. The Operation Performed on the Octree

Here, we define four kinds of operations, including the equivalenting source operation, up operation, down operation, and horizontal operation. Figure 9 shows the schematic diagram of these operations. Among these operations, the equivalenting source operation and up operation are essentially the same as S2M and M2M in [27], while the down operation and horizontal operation are improved upon L2L and M2L in [27].

Equivalenting source operation: For a leaf, using the charges on its inner auxiliary surface to be equivalent to those triangular elements inside it and calling this process the equivalenting source operation for the leaf. The sampling points lie on the outer auxiliary surface of the leaf.

Up operation: For a node, excluding the root node, using the charges on its parent’s inner auxiliary surface to be equivalent to the charges on its inner auxiliary surface and calling this process the up operation for the node. The sampling points lie on the outer auxiliary surface of the node’s parent.

Down operation: For a node, excluding the root node, using the charges on its outer auxiliary surface to be equivalent to the charges on its parent’s outer auxiliary surface and calling this process the down operation for the node. The sampling points lie on the inner auxiliary surface of the node.

Horizontal operation: For a node A and a node B, which is one of the interactive nodes of the same layer of A, using the charges on the outer auxiliary surface of A to be equivalent to the charges on the inner auxiliary surface of B and calling this process the horizontal operation from B to A.

3.4. Main Body of the Algorithm

3.4.1. Calculating Equivalent Charges

The general flow of calculating equivalent charges is to iterate from bottom to top first, calculate the charges on the inner auxiliary surfaces, and then iterate from top to bottom, and calculate the charges on the outer auxiliary surfaces. See Algorithm 1 for details. Algorithm 1. Calculating equivalent chargesInput: Constructed octree, needed matrixes, triangular elements informationOutput: The equivalent charges on two auxiliary surfaces for each node1parfor each leaf, i2  Executing the first and second steps of equivalenting source operation for i in sequence3end parfor4for i = tree_hight : 25  parfor each node, j, in layer i of the octree6    Executing up operation for j7  end parfor 11end for12for i = 2 : tree_height13  parfor each node, j, in layer i of the octree14    for each interactive node of the same layer, k, of j15      Executing horizontal operation from k to j16    end for 17  end parfor 18   $[eqn]$ (if it exists) of the octree19    Executing down operation for j20  end parfor 24end for25return the equivalent charges on two auxiliary surfaces for each node

The reason why we split the equivalenting source operation (line 2) into two steps instead of executing the whole directly will be explained in Section 3.5.1.

Different from calculating potential distribution, when obtaining force and stiffness, we only need the electric field provided by those triangular elements on the housing. So, we should rebuild the octree after we set $[eqn]$ for each $[eqn]$ .

To accelerate the calculation, we can adopt multi-dimensional GPU parallel computing. The starting node and ending node of each operation, as well as the matrix operation corresponding to the operation, can all correspond to different dimensions of multi-dimensional threads.

3.4.2. Updating Expressions

Due to the equivalence of charges on triangular elements, those relevant expressions in Section 2 are no longer applicable and need to be modified accordingly. $[eqn]$ stands for the set composed of subscripts of the triangular elements in all neighbor nodes of the leaf where the i-th triangular element is located, and $[eqn]$ stands for the list composed of the ratio of charge quantity to vacuum dielectric constant of each equivalent charge on the outer auxiliary surface of the leaf where the i-th triangular element is located. In addition, $[eqn]$ stands for the set composed of the positions of each equivalent charge on the outer auxiliary surface of the leaf where the i-th triangular element is located. The expression from Equation (3) that describes the potential at the centroid of the i-th triangular element should be modified to

$$[eqn]$$

where $[eqn]$ is the value of parameter p for the equivalent charges of the four kinds of operations.

If $[eqn]$ is the set formed by the subscripts of the triangular element located on the housing in all the neighboring nodes of the leaf where the i-th triangular element is located, and $[eqn]$ is the list formed by the relative charge amounts of each equivalent charge on the outer auxiliary surface when the charge on the surface of the TM is set to 0, then, the expression (6), which describes the force on the TM in the X-axis direction, should be replaced by (13), and the expression (7), which describes the rate of change of the force on the TM along the x axis with respect to the displacement of the TM along the x axis, should be modified to (14).

$$[eqn]$$ $$[eqn]$$

With the help of the octree, we can improve the efficiency of calculating $[eqn]$ , $[eqn]$ , and $[eqn]$ . Just like in Section 3.4.1, calculating $[eqn]$ , $[eqn]$ , and $[eqn]$ can be regarded as performing matrix operations on each leaf node, and multi-dimensional GPU parallel computing can still be used for acceleration.

3.5. Algorithm Complexity

The volume of leaves is determined, so it can be approximated that the average amount of triangular elements, $[eqn]$ , inside each leaf is determined. The number of leaves is $[eqn]$ . According to the properties of the octree, the number of nodes is $[eqn]$ .

Let $[eqn]$ stand for the value of parameter p for the sampling points of all the equivalenting source operations. The following analysis is based on the assumption that $[eqn]$ , $[eqn]$ , and $[eqn]$ are constant.

3.5.1. Space Complexity

Each operation has a corresponding matrix which can be translated into two matrixes corresponding to the first and the second step of the operation. Both the up operation and down operation have at most eight different relative positions of one node and its parent. The horizontal operation contains at most $[eqn]$ different relative positions of two nodes in the same layer of an octree [29]. The function $[eqn]$ has the property that for each $[eqn]$ , $[eqn]$ . Making use of this property, we can execute the up operation, down operation, and horizontal operation with $[eqn]$ space complexity [29,30]. Because of the uncertain distribution of triangular elements in each leaf, if we want to store the matrixes of the equivalenting source operation, we need the space complexity of $[eqn]$ . However, we can store the matrix of the second step of the equivalenting source operation with $[eqn]$ , and if we do not precompute the matrixes of the first step, we can also finish the equivalenting source operation with $[eqn]$ space complexity. To optimize computational efficiency, we precompute the matrix of the second step, thereby circumventing the computationally intensive matrix inversion operation inherent in the solution procedure.

To store the information of triangular elements, the space complexity is $[eqn]$ . To store the information of the octree, the space complexity is $[eqn]$ .

In summary, the total space complexity is $[eqn]$ .

3.5.2. Time Complexity

For the convenience of analysis, the acceleration brought by the use of GPU parallel computing is not considered. First, we analyze the time complexity of calculating equivalent charges in one iteration. For each leaf, the corresponding equivalenting source operation should be executed exactly once in the whole procedure. So, the total time complexity of the equivalenting source operation is $[eqn]$ . For each node, the computational procedure requires a maximum of one up operation, one down operation, and $[eqn]$ times horizontal operation [30]. The total time complexity of the up operation, down operation, and horizontal operation is $[eqn]$ , $[eqn]$ and $[eqn]$ , respectively. In summary, the total time complexity of calculating equivalent charges in one iteration is $[eqn]$ .

Second, we analyze the time complexity of calculating potential, force, and stiffness. Each node has at most 27 neighbors and $[eqn]$ charges on the outer auxiliary surface, so the total time complexity is $[eqn]$ .

4. Experiments, Results, and Discussion

4.1. Parameter Setting

For computational simplicity and analytical convenience, we model the TM as an ideal cube (with uniform side lengths of 46 mm), while the inner surface of the housing is approximated as an ideal cuboid surface. Based on the coordinate system defined in Section 2.1, Figure 10 shows the distribution of electrodes.

The red electrodes are injection electrodes, and the green electrodes are sensing and actuation electrodes. The other area of the surface of the housing can be approximated as grounding. According to Ref. [26], alternating current was applied to each electrode. We select a moment and assume that at that moment, the potential on all sensing electrodes is 1000 mV, and the potential on all injection electrodes is 3 mV. In addition, we assume $[eqn]$ .

To achieve surface subdivision, we first divide the domain into grids of squares with side length 0.1 mm; each square is subsequently divided into four triangular elements (as shown in Figure 11), forming the fundamental mesh unit. This process yields a total of 11,802,080 triangular elements.

4.2. Single Contaminant Patch Bulge Impact Analysis

In this section, we will analyze the influence of a single contamination patch bulge on the system, focusing on three key parameters: contaminant size, location, and potential variation. To streamline the analysis process, we adopt $[eqn]$ in Equation (8), where $[eqn]$ represents the characteristic bulge height.

Before our experiments, we verified the correctness of our simulation algorithm. See Appendix A for detailed procedures.

To begin with, we calculated the electrostatic force $[eqn]$ and the electrostatic stiffness $[eqn]$ under the circumstance that no patch exists. In theory, $[eqn]$ . The value of $[eqn]$ can be thought of as the error due to the subdivision, algorithm, and rounding error of floating-point numbers.

4.2.1. Impact of Potential Δu

To characterize the electromechanical coupling effects under surface contamination, we conducted systematic variations of the electrostatic potential difference $[eqn]$ from $[eqn]$ to 150 $[eqn]$ with a step $[eqn]$ . Throughout the experiments, the patch base radius $[eqn]$ was fixed at $[eqn]$ mm, and the center of the bulge’s base was positioned at the coordinate $[eqn]$ As shown in Figure 12, a clear linear relationship is observed between $[eqn]$ and $[eqn]$ and the potential difference Δu. Notably, the electrostatic force change $[eqn]$ is defined as $[eqn]$ , while the stiffness change $[eqn]$ is characterized as $[eqn]$ .

We tried to fit the data using linear functions and obtained the following expressions

$$[eqn]$$ $$[eqn]$$

The actual meanings of their intercepts are the changes in the force or stiffness when the bulge has no influence on potential. Therefore, to reduce the error, we calculated the force and stiffness when the change in electric potential caused by the bulge was zero and assigned these values as the intercepts. During fitting, only the slope was adjusted.

4.2.2. Impact of Base Radius rm

To investigate the dependence of $[eqn]$ and $[eqn]$ on the geometric scale of the bulge, we selected $[eqn]$ mm, set $[eqn]$ , and chose point $[eqn]$ as the center of the bottom of the bulge. Figure 13 illustrates the relationship between $[eqn]$ and the deviations $[eqn]$ and $[eqn]$ .

Theoretically, when $[eqn]$ , both $[eqn]$ and $[eqn]$ are expected to vanish, so we attempted to fit $[eqn]$ and $[eqn]$ using cubic polynomials passing through the origin (see “Fitted curve (A)” in Figure 13). The resulting fitting expressions are:

$$[eqn]$$ $$[eqn]$$

In Equations (17) and (18), the coefficients of the cubic terms are several orders of magnitude larger than those of the remaining terms. Based on this observation, we inferred that $[eqn]$ and $[eqn]$ . To verify this conclusion, we re-fitted the data. The new fitting results are shown in “Fitted curve (B)” in Figure 13, and the formulas are:

$$[eqn]$$ $$[eqn]$$

The fitting results confirm that $[eqn]$ and $[eqn]$ exhibit a cubic dependence on $[eqn]$ .

4.2.3. Impact of Location

To further investigate the spatial dependence of contaminant path effects, we analyzed the positional sensitivity of patch effects by selecting four distinct bulge locations: Locations 1 and 2 on the TM surface at (−23 $[eqn]$ , 0 $[eqn]$ , 0 $[eqn]$ ) and (−23 $[eqn]$ , −11 $[eqn]$ , 0 $[eqn]$ ), with Location 2 positioned opposite an electrode, and Locations 3 and 4 on the housing surface at (−27 $[eqn]$ , 0 $[eqn]$ , 0 $[eqn]$ ) and (−27 $[eqn]$ , −11 $[eqn]$ , 0 $[eqn]$ ), with Location 4 on an electrode. For all cases, we set $[eqn]$ and $[eqn]$ .

Both Location 1 and Location 2 are at the surface of the TM. Location 2 is opposite one electrode. Both Location 3 and Location 4 are at the surface of the housing. Location 4 is at one electrode.

From Figure 14, it can be observed that patches opposite to an electrode or located at one electrode exhibit a bigger impact on $[eqn]$ compared to those situated in other regions. Conversely, patches on the housing surface demonstrate a significantly greater effect on $[eqn]$ than those implemented on the TM surface.

4.3. Single Planar Patch Impact Analysis

In this section, we will investigate the influence of the variation of potential, size, and position of a single planar patch. The condition $[eqn]$ is adopted in Equation (8).

4.3.1. Impact of Potential Δu

We vary $[eqn]$ from $[eqn]$ to $[eqn]$ in steps of $[eqn]$ , with $[eqn]$ and the patch centered at $[eqn]$ . Figure 15 shows $[eqn]$ and $[eqn]$ under different $[eqn]$ . We fitted the data using a proportional function, as theoretically, when $[eqn]$ , both $[eqn]$ and $[eqn]$ should equal zero. The fitted equations are:

$$[eqn]$$ $$[eqn]$$

According to the result, we can conclude $[eqn]$ and $[eqn]$ .

To further validate this conclusion, we conducted additional experiments in which the patch was located at $[eqn]$ and $[eqn]$ . Figure 16 shows the results, and the fitted proportional functions are given in Equations (23) and (24).

$$[eqn]$$ $$[eqn]$$

The conclusion that $[eqn]$ and $[eqn]$ has been further supported.

An additional observation is the reversal of monotonicity in $[eqn]$ and $[eqn]$ between Figure 15 and Figure 16. This phenomenon will be explained in Section 4.3.3.

4.3.2. Impact of Patch Radius rm and Position

To investigate the impact of patch radius $[eqn]$ and position on the electrostatic force and stiffness, we selected three distinct points, $[eqn]$   $[eqn]$ and $[eqn]$ as the center of the patch and conducted three sets of experiments. In each experiment, for each $[eqn]$ , we set $[eqn]$ and calculated $[eqn]$ and $[eqn]$ .

In the first set of experiments, the patch was centered at coordinates $[eqn]$ and the results are shown in Figure 17.

We tried to fit the data. Considering the accuracy of the fit and the simplicity of the expression, we obtained the following fitting results:

$$[eqn]$$ $$[eqn]$$

In the second set of experiments, the patch was centered at coordinates $[eqn]$ and the results are shown in Figure 18. We fitted the data and obtained the expressions:

$$[eqn]$$ $$[eqn]$$

In the third set of experiments, the patch was centered at coordinates $[eqn]$ and the results are shown in Figure 19. We fitted the data and obtained the following expressions:

$$[eqn]$$ $$[eqn]$$

As can be readily observed from the results above, the changes in the electrostatic force and stiffness depend on both the position and size of the patch. When the patch is positioned directly opposite the electrode, we find that changes in electric forces and stiffness satisfy $[eqn]$ , $[eqn]$ , given $[eqn]$ , and the magnitudes of such changes increase with the patch radius $[eqn]$ . When the patch is located opposite the non-electrode region of the housing, we find that $[eqn]$ and $[eqn]$ for a given value of $[eqn]$ and the magnitudes of such changes increase with the patch radius $[eqn]$ When the patch faces the edge of the electrode, the patch effects can be considered as the superposition of those from the above two cases. We will explain it in Section 4.3.3.

4.3.3. Explanation of the Curve of ΔFx and ΔKxx

To explain the curve shape of $[eqn]$ and $[eqn]$ mentioned above, we analyzed the charge redistribution induced by the patch and its interaction with the housing electrodes. Figure 20 illustrates the charge density ( $[eqn]$ ) distribution on the TM and housing surfaces.

To verify the charge distribution in Figure 20, we plotted the equipotential distribution on three planes ( $[eqn]$ , $[eqn]$ , and $[eqn]$ ) derived from the charge distribution, as shown in Figure 21.

Figure 21 shows that: (i) the potential in the ground region of the housing is zero, (ii) the potential on the electrode matches the predefined potential value, and (iii) the TM is an equal potential body. These results thereby verify the correctness of the charge distribution.

Baseline Charge Distribution (Figure 20a,b): Under a 1000 mV bias applied to the sensing and actuation electrodes, the housing electrode regions exhibit strong positive charge accumulation, while the ground regions display negative charge accumulation (Figure 20b). Electrostatic coupling induces complementary charge polarization on the TM: negative charges dominate regions opposite housing electrodes, whereas positive charges accumulate elsewhere (Figure 20a).

Patch-Induced Charge Redistribution (Figure 22a–f): When a patch is present with $[eqn]$ , additional positive charges accumulate within the patch region (Figure 22a,c,e), with negligible boundary effects due to the small boundary area. Additional negative charges accumulate on the surface of the housing opposite to the patch (Figure 22b,d,f).

From a physical point of view, $[eqn]$ , $[eqn]$ are determined by the interactions among the surface charges from the TM and housing. The dominant contributions are from the interactions between the baseline background charge distributions $[eqn]$ and $[eqn]$ . The leading corrections arise primarily from interactions between the patch-induced charge and the original background, i.e., two pairs of interactions between $[eqn]$ and $[eqn]$ and between $[eqn]$ and $[eqn]$ . As the magnitude of the patch-induced charge density is relatively small, the second-order corrections to the forces and stiffness from interactions between the patch-induced charge distribution can be ignored.

According to Equations (13) and (14), for the interaction between $[eqn]$ and $[eqn]$ ( $[eqn]$ and $[eqn]$ ), if $[eqn]$ ( $[eqn]$ ), then $[eqn]$ and $[eqn]$ ; if $[eqn]$ ( $[eqn]$ ), then $[eqn]$ and $[eqn]$ .

In Figure 15, the patch is opposite the non-electrode region and $[eqn]$ ; we observe that $[eqn]$ and $[eqn]$ , leading to $[eqn]$ and $[eqn]$ . On the contrary, in Figure 16, the patch is opposite an electrode, $[eqn]$ ; we observe that $[eqn]$ , $[eqn]$ > 0, resulting in $[eqn]$ and $[eqn]$ . The magnitude of these deviations scales linearly with $[eqn]$ .

In Figure 17, no matter how $[eqn]$ changes, the patch is always only opposite an electrode, and $[eqn]$ , so $[eqn]$ , so $[eqn]$ , $[eqn]$ . The curves become steeper and steeper because the area of the patch is a quadratic function of $[eqn]$ .

In Figure 18, no matter what the value of $[eqn]$ is, half of the patch is always opposite an electrode, and the other half is not opposite. The shape of the curve is determined by the superposition of the effects of the two parts.

In Figure 19, $[eqn]$ , which determines $[eqn]$ and $[eqn]$ when the patch is not big enough, it is only opposite the ground area on the housing, so $[eqn]$ and $[eqn]$ , $[eqn]$ . The curve gets steeper and steeper. As $[eqn]$ expands further, a part of the patch will be opposite an electrode. In this part, $[eqn]$ ; the superposition of this part with the former causes the concavity, monotonicity, and steepness of the curves to be affected.

4.4. Analysis of the Metric Requirements in Space-Borne Gravitational Wave Detection Missions

The acceleration noise of inertial sensors can be classified into two types: direct noise and coupled noise. The direct noise of the patch field effect needs to be calculated based on a large amount of time series data related to the patch parameters, which has extremely high requirements for the efficiency of the simulation algorithm and computing power resources. This section does not conduct research on direct noise. Coupled noise is formed by the coupling of the position fluctuation of the TM and stiffness.

This section will respectively analyze the parameter limitations of bulges or large planner patches. In extreme cases, the stiffness of the patch field effect is approximately $[eqn]$ [10]. In this section, we adopt a stricter value, with the stiffness of the patch not exceeding $[eqn]$ (that is, the stiffness of the force is $[eqn]$ when the mass of the TM is 1.93 kg).

4.4.1. Patch with Bulge

Figure 14 indicates that the change in stiffness at different positions is not significant. Therefore, we ignore the differences in the effects of bulges at different locations.

According to Equation (20), when $[eqn]$ , $[eqn]$ . If we let $[eqn]$ then we can obtain $[eqn]$ . This is a very broad range. It can be considered that the impact of the bulge of a single patch can be ignored. If there are N bulges that do not affect each other, and the value of $[eqn]$ of each bulge is 0.1 mm, we can obtain $[eqn]$ according to $[eqn]$ . Such a requirement is very lenient.

According to Equation (16), when $[eqn]$ , $[eqn]$ . If we assume there are N bulges that do not affect each other and let $[eqn]$ then we can obtain $[eqn]$ . Such restrictions can still be ignored.

4.4.2. Single Planar Patch

In this section, we look for the values of $[eqn]$ and $[eqn]$ that meet the stiffness limitation when the center of the patch is located at $[eqn]$ and $[eqn]$ , respectively.

Figure 23 shows the distribution of acceptable combinations of $[eqn]$ and $[eqn]$ . If the tuple $[eqn]$ is acceptable, the point will be marked as blue. With the increase in $[eqn]$ , the acceptable values of $[eqn]$ in Figure 23a,b are gradually decreasing. When $[eqn]$ increases to 150 mV, the acceptable values of $[eqn]$ in Figure 23a,b decrease to 3.16 mm and 1.62 mm, respectively. In a word, in order to ensure that the stiffness meets the requirements when $[eqn]$ , $[eqn]$ needs to be less than 1.62 mm.

5. Conclusions

In this paper, we provided a method that can be used to simulate the influences of a given patch field on the force and the stiffness on the TM and used this method to conduct some experiments. We can draw conclusions as follows:

We established a mathematical model based on a partition of boundary elements and the GMRES method to simulate the patch field. This model can solve the charge distribution according to the potential distribution and use the charge distribution to calculate the force and stiffness on the TM. The bulges and the nonuniform distribution of the potential led by patches can also be simulated with the model.To overcome the difficulty of excessive computational complexity, based on existing algorithms, we made improvements according to needs, and designed an algorithm with $[eqn]$ spatiotemporal complexity for calculating potential distribution, force, and stiffness.With the method mentioned above, we researched what impacts single bulge forms from contaminant attachment have on force, $[eqn]$ , and stiffness, $[eqn]$ . The control variable method is adopted, and $[eqn]$ , $[eqn]$ , and location are taken as variables, respectively. The results show that both $[eqn]$ and $[eqn]$ are linear functions of $[eqn]$ , approximately proportional to $[eqn]$ to the third power. The patch which is opposite one electrode or is at one electrode has a bigger impact on $[eqn]$ than the patch in the other area. The patch at the surface of the housing has a bigger influence on $[eqn]$ than the patch at the surface of the TM.In addition, we also studied a single patch without a bulge and found that the relation of $[eqn]$ and $[eqn]$ to $[eqn]$ can be approximated to a quartic curve passing through the origin, which can be simplified in some cases, and both $[eqn]$ and $[eqn]$ are approximately proportional to $[eqn]$ .With the help of Conclusion 3~4, we respectively studied the bulges and single planar patch for space-borne gravitational wave detection missions. If we limit the stiffness caused by the patch not to exceed $[eqn]$ , we found that, under normal circumstances, the impact of a bulge can be ignored. When $[eqn]$ , the value of $[eqn]$ of a single patch should be less than 1.62 mm.

This study focuses on patch effects in inertial sensors and establishes a simulation framework that is scalable to the case where alternating current (AC) voltage is applied to actuation electrodes, as well as to the multi-patch case. Future research will focus on more complex scenarios and exploring mitigation strategies to reduce the impact of patch fields on gravitational wave detection. We will also try to study the general nonuniform distribution of electric potential caused by patch field effects based on the Monte Carlo method, rather than single patches, but this requires a large number of experiments to reach a statistically significant conclusion. Due to the large scale of the partition of our simulation experiment, we will further optimize the algorithm and improve the computing power resources to complete this part of the work. Our partner team is conducting patch effect experiments simultaneously, and there will be an actual measurement of the patch distribution in the future.

In addition to the patch field effects described in this article, other interference sources affecting the TM include temperature variations, residual gas interactions, cosmic ray exposure, self-gravitation force, magnetic fields, and so on. The influence of the patch field effects may exhibit dependencies on these factors. Moreover, temporal variations must be systematically taken into consideration. In future work, we will analyze the coupled interactions and incorporate time-varying analyses during satellite orbital operations.