1 I. INTRODUCTION

The Rayleigh–Taylor instability (RTI)1,2 is a fundamental hydrodynamic instability found in many natural phenomena, such as supernova explosions3–6 and jet-driven lobes in galaxy clusters.7 RTI is also a critical concern in inertial confinement fusion (ICF).8 With the recently reported milestone achievements in ignition,9,10 deliberate ignition designs are required for better ICF implosion performance by increasing the coupling efficiency between the laser energy and the capsule and controlling the hydrodynamic instabilities. In ICF implosions,11 the laser irradiates either directly on the fusion pellet or indirectly on a high-Z hohlraum to transfer the laser energy into x-ray energy. When irradiated by laser or x-ray energy, the outer shell of the fusion pellet absorbs the radiation energy and transforms into high-temperature low-density plasma that is then ejected from the pellet at high speed. This ablation leads to the shell’s inward acceleration and compression due to momentum conservation, and during this acceleration phase the outer shell surface is unstable to RTI because a lighter fluid (ablated plasma) is pushing on a denser fluid (unablated plasma). The development of RTI seeded by short-wavelength target surface roughness or laser imprint12 will cause the interchange of fluids: the bubbles of lighter fluid rise through the denser fluid, while the denser fluid penetrates down through the lighter fluid in the form of narrow spikes. The integrity of the target shell is severely compromised by this material mixing, which can degrade the implosion performance and prevent thermonuclear ignition.

Mass ablation is a key feature in the RTI at the outer surface of an ICF pellet, and ablative RTI (ARTI) evolves very differently from classical RTI (CRTI) because of this ablation. Previous studies on ARTI or CRTI were focused mostly on RTI developing from initial sinusoidal perturbations in the density or velocities near the interface where the low-density (ρ_l) fluid is supporting the high-density (ρ_h) fluid against a gravitational field g. The linear theory2 for single-mode CRTI shows that initial infinitesimal sinusoidal modulations with wavelength λ grow exponentially in time with growth rate γCRTI=ATkg, where k = 2π/λ is the perturbation wave number and A_T = (ρ_h − ρ_l)/(ρ_h + ρ_l) is known as the Atwood number. It is generally considered that when the mode amplitude reaches ∼0.1λ, the linear growth stage ceases and the CRTI growth becomes nonlinear. The bubble velocity U_b—which is the speed at which the bubble vertex penetrates the dense fluid and so is an important quantity for assessing RTI risk—then saturates and approaches a constant terminal value U_cl predicted by inviscid potential-flow theory13,14 in the highly nonlinear phase. The 2D terminal velocity is Ucl2D=g(1−rd)/3k, where r_d = ρ_l/ρ_h is the density ratio of the light and heavy fluids.

As is well known, ablation lowers the linear RTI growth rate by supplying a restoring force when the interface is distorted, and this stabilizing effect is more substantial for short-wavelength modes. The linear growth rate of ARTI can be approximated by the Takabe–Bodner15,16 formulaγARTI=αkg−βkVa,(1)where Va=ṁ/ρa is the ablation velocity representing the penetration velocity of the ablation front into the heavy shell, ṁ is the mass ablation rate per unit area, and ρ_a is the density of the unablated shell at the ablation interface. The terms α and β are parameters that depend on the ablation material15 and are chosen to fit the numerical simulation data; for a deuterium–tritium (DT) ablator, a proper choice is α = 0.94 and β = 2.7.17 Equation (1) works well with a sharp ablation front (kL_m ≪ 1, where L_m is the minimum density gradient scale length near the ablation front)18 but is less precise when there is significant radiation energy transport, which leads to smooth density profiles.17 An improved formula that includes the finite thickness of the ablation front can be found in Refs. 19 and 20 as γ=kg/(1+kLm)−βkVa, which agrees well with experiments.21 More-complicated but more-precise self-consistent stability theories17,22–27 were developed for the limits of large or small kL₀ and small or large Froude number Fr=Va2/(gL0). For large Fr, the linear growth rate for ARTI can be formulated as17γARTI=ATkg−AT2k2Va2/rd−(1+AT)kVa,(2)where L₀ is the characteristic thickness of the ablation front, which is related to L_m and the power index for the thermal conduction (κ ∼ T) as L₀ = L_mν/(ν + 1)⁽⁺¹⁾. γ_ARTI is lower than γ_CRTI under the same density profile in the presence of V_a, which imposes an upper cutoff wave number k_c (i.e., a lower cutoff wavelength λ_c) on the instability spectrum, and all modes with k > k_c (i.e., λ < λ_c) are linearly stable.17,22,23,25 The ablative stabilization effects on RTI have been nicely demonstrated experimentally,6,28 with decreased RT growth observed if the interface is ablated by laser-driven energy fluxes.

Like CRTI, ARTI also saturates after its linear phase, and U_b approaches U_cl when the mode amplitude is sufficiently large for the ablation effect on the bubble vertex to be negligible. However, it has been found that ARTI bubble velocities can exceed U_cl because of the vortex acceleration mechanism in both 2D29 and 3D30 geometry, especially for the short-wavelength modes near the linear cutoff, and vortex acceleration has been verified in experiments.31 This destabilizing effect caused by ablation is a potential risk to target integrity because of the short-wavelength modes in ICF implosions. In its highly nonlinear phase, ARTI seeded by multiple-mode perturbations eventually reaches a self-similar32 turbulent-like regime, which is of intense research interest in both the ICF and astrophysical communities.33–38

Under ICF-relevant experimental conditions, ARTI can grow because of different types of imperfections. The nonuniformity of the irradiation (e.g., laser imprint) and the nonuniformity of the target (e.g., target roughness) can both cause perturbations at the ablation front. While the nonuniformity of the target seeds ARTI as an initial condition, the nonuniformity of the irradiation can serve as a continuous perturbation source coming from the ablation side. Mitigating laser imprints in ICF has been pursued in different ways. It was recently found that laser imprints can be mitigated using hybrid high-Z target designs,39,40 and another mitigation concept is to move the laser speckles around with time. In direct-drive ICF schemes, the required high irradiation uniformity is achieved via beam-smoothing techniques, including distributed phase plates,41,42 polarization smoothing,43,44 and 2D smoothing by spectral dispersion (SSD).45–47 In particular, SSD is known to shift the laser imprint dynamically; the strategy is to vary the interference pattern in the laser speckle on a time scale that is shorter than the characteristic hydrodynamic response time of the target so that the nonuniformity of laser intensity can be smoothed in time,45 which also leads to time-varying perturbations of the heat flux on the laser absorption zone in the ablated plasma.

Herein, we use numerical simulations to study ARTI growth induced by heat-flux perturbations. Both stationary and time-varying heat-flux perturbations are investigated, and the time period of the variation is set in a broad range that includes the SSD-relevant modulation periods used in the leading ICF facilities.46,48,49 The rest of this paper is organized as follows. In Sec. II, the simulation methods and configurations are introduced. In Secs. III and IV, the simulation results for ARTI evolution induced by stationary (Sec. III) and time-varying (Sec. IV) heat-flux perturbations are presented and discussed. Finally, the paper concludes with a summary in Sec. V.

2 II. SIMULATION SETUP

In this study, we used the hydrodynamic code ART,29,30 which was designed for simulating ARTI in ICF-relevant regimes. ART solves the single-fluid hydrodynamic equations together with the Spitzer–Harm (SH) model of thermal conduction50 on a Cartesian grid. The hydrodynamic equations are as follows:∂ρ∂t+∇⋅(ρu)=0,∂(ρu)∂t+∇⋅(ρuu)+∇P=ρg,∂E∂t+∇⋅[(E+P)u]=ρu⋅g,(3)where ρ, P, u, and g are the fluid density, pressure, velocity, and acceleration, respectively, E = P/(Γ − 1) + ρ|u|²/2 is the total energy, and Γ = 5/3 is the ratio of specific heats. The thermal conduction part in the energy equation is solved separately via Strang splitting51 to avoid the strict time step (Δt) required by the explicit heat diffusion equation solvers; a similar treatment is commonly used in high-energy-density physics simulations.52,53 The thermal conduction equationcvρ∂T∂t=∇⋅(κ∇T)(4)is solved after the hydrodynamic equations are advanced in each time step with the same Δt, where c_v is the specific heat at constant volume and κ is the thermal conductivity coefficient provided by the classical SH model with no flux limiter as κ ∼ T^2.5. The SH heat fluxes in the whole simulation domain are found to be smaller than the free-stream heat fluxes times a commonly used flux-limiting factor (0.065). The corrections to P and E are then applied after solving Eq. (4) by using the equation of state for a perfect gas. A fifth-order WENO-JS54 scheme in spatial finite difference is used with a third-order Runge–Kutta method advancing in time. ART has been benchmarked and used in a number of ICF-relevant ARTI studies.29,30,37,38,55–59

A typical simulation setup is shown in Fig. 1. In the upper half is the cold and dense unablated DT shell, and in the lower half is the ablated plasma with higher temperature but lower density [see Figs. 1(a) and 1(c)]. The initial profile is set to a quasi-equilibrium state that is typical for direct-drive experiments on the National Ignition Facility,60,61 which is calculated by integrating the hydrodynamic equations below and above the ablation front with different approximations based on the fluid parameters on the ablation front and the entropy profile in the shell. Figure 1(c) shows the profiles of initial quantities along the z direction, these being the density, temperature, and velocity in the negative z direction. The ablation velocity is set to be V_a = 2.0, 3.5, and 5.0 μm/ns for different cases, giving a Froude number in the range of 0.1 < F_r < 2.

Fig. 1. Simulation setup for λ = 70 μm and V_a = 3.5 μm/ns: (a) initial density profile; (b) temperature fluctuations at t = 0.08 ns; (c) initial profiles of ρ (solid), T (dot-dashed), and v_z (dashed) along z axis.

下载图片 查看所有图片

The density profile has no initial perturbations, so the ablation front is initially flat. To keep the ablation front approximately fixed in the box, the acceleration is set initially to g₀ = 100 μm/ns² and then adjusted dynamically to balance the ablation pressure during the simulation as g(t)=[(P+ρu2)bot−(P+ρu2)top]/Mtotal, where the subscripts “top” and “bot” indicate the integral values on the top and bottom boundaries, respectively, and M_total is the total mass in the whole computational region. This is equivalent to solving the fluid equations in the frame of reference of the accelerated shell. Because the mass of the total simulation volume decreases with time because of the ablation, g must be increased to keep the ablation front fixed.

The energy flux transported from the laser-absorption region toward the ablation front is simulated by a heat flux applied on the bottom boundary of the simulation box (z = 0). The initial hydrodynamic profiles including temperature are obtained by integrating the 1D hydrodynamic equilibrium equations in the frame of reference of the shell from the ablation front toward both sides, then the average heat flux on the bottom boundary (Q₀) is calculated initially as Q₀ = −κ∇T evaluated on the boundary and kept constant throughout a simulation. Equilibrium hydrodynamic boundary conditions are used for the vertical direction, and periodic boundary conditions are used for the horizontal direction. The top boundary is set to be adiabatic because the majority of the heat flow is stopped near the ablation front. The dense fluid in the upper half can be approximated well by a piston, so ∂_zv_z = 0 is used for the inflow boundary. A typical simulation box for the case λ = 70 μm is 70 × 120 μm² with a 700 × 1200 grid, and a grid convergence study was performed to ensure numerical convergence.

The perturbed heat flux on the bottom boundary is set as Q(x) = Q₀[1 + A · cos(kx)] for a stationary perturbation and Q(x, t) = Q₀[1 + A · cos(kx + ωt)] for a time-varying perturbation in the form of a “wave” traveling along the −x direction. Here, A is the ratio of perturbation amplitude to Q₀ and was set as A = 10% in all the simulations reported herein. Furthermore, we have ω = 2π/τ, where τ is the period of the wave and was set in the range of 0.05–1 ns in our simulations, covering the frequencies of the laser interference patterns in ICF-relevant SSD regimes.48,49

Figure 1(b) shows the temperature fluctuations at t = 0.08 ns after the stationary heat-flux perturbation of λ = 70 μm is assigned initially in the case of V_a = 3.5 μm/ns. The temperature fluctuations are calculated by subtracting (i) the quasi-steady temperature field in a reference simulation with a uniform heat flux Q₀ applied on the bottom boundary from (ii) the temperature field with the heat-flux perturbation. As can be seen, the assigned heat-flux perturbation on the bottom boundary induces considerable fluctuations in the temperature field inside the box, which then perturb the ablation front and induce ARTI growth.

3 III. ARTI INDUCED BY STATIONARY HEAT-FLUX PERTURBATION

We performed a series of stationary heat-flux (SHF) perturbation simulations with λ ranging from 7 μm (which is close to the linear cutoff) up to 80 μm, and the simulation results are shown in Fig. 2. A characteristic Peclet number P_e representing the rate ratio of advection to diffusion is calculated as P_e = LU/α ≈ 0.45, where L is chosen as the distance from the bottom boundary to the ablation front, U is a characteristic fluid advection velocity chosen as |v_z| at the bottom boundary, and α = κ/(ρc_p) is the thermal diffusivity at the bottom boundary. It is usually considered that advection dominates diffusion when P_e ≫ 1, which can prevent heat-flux perturbations from transporting upward to the ablation front. In our simulations, we found that the SHF perturbation applied on the bottom boundary could penetrate into the ablated plasma, perturb the interface (i.e., the ablation front), and induce ARTI growth. However, unlike the usual velocity or density perturbations applied initially on the interface, as the perturbed heat flux penetrates into the plasma toward the ablation front, the amplitude of the heat-flux perturbation decreases with the penetration distance because of the transverse heat conduction [see Fig. 2(c)]. Figure 2(c) also shows that this traverse thermal smoothing effect is enhanced for shorter-wavelength modes, which have finer transverse structures and thus mitigated actual perturbations at the interface.

Fig. 2. (a) Ratio of bubble velocity to corresponding classical value U_cl2D and (b) linear growth rate of ablative Rayleigh–Taylor instability (ARTI) induced by stationary heat-flux (SHF) perturbation for different λ with V_a = 3.5 μm/ns. The linear growth rates are shown for the simulation results (red triangles) and those obtained theoretically using the improved Takabe-like formula [Eq. (1)] (pink dot-dashed line) and the formula based on stability theory [Eq. (2)] (black solid line). (c) x-averaged temperature perturbation amplitude vs z penetration distance induced by SHF perturbation for three cases with different λ at t ≈ 0.5 ns.

下载图片 查看所有图片

The ratio of the ARTI bubble velocity U_b to the corresponding classical value U_cl2D for different λ is plotted in Fig. 2(a). U_b is measured as the velocity of the bubble vertex relative to the dense plasma in the simulations, and U_cl2D is calculated using the quantities obtained from the simulation data. As can be seen, all modes experience a linear growth phase and then saturate as U_b approaches U_cl2D, then the vortex acceleration mechanism29 drives U_b slightly beyond U_cl2D. The shorter-wavelength modes grow substantially later than the longer-wavelength modes, which is consistent with the enhanced transverse thermal smoothing shown in Fig. 2(c). In the simulations with λ ≤ 10 μm, ARTI did not grow until very late in time, which is why the corresponding results are absent from Fig. 2.

Although the heat-flux perturbation is supplied continuously at the bottom boundary, it is not expected to alter the linear ARTI growth rate significantly because the perturbation actually applied at the ablation front is rather small whereas ARTI is an intrinsically unstable process that grows exponentially in time. In Fig. 2(b), the linear growth rates computed in the simulations are compared with those obtained using Eqs. (1) and (2). The values of g and r_d used in the calculations are the average values during the linear phases in the simulations. Good agreement is reached between simulation and theory, verifying that ARTI induced by SHF perturbation has a similar growth rate as that induced by the velocity or density perturbations near the interface.

4 IV. ARTI INDUCED BY TIME-VARYING HEAT-FLUX PERTURBATION

When the heat-flux perturbation on the bottom boundary varies with time in the form of a wave traveling horizontally, the density and temperature fluctuations inside the ablated plasma react also in a wave-like manner, which leads to new features of the ARTI evolution. We performed a series of ART simulations with time-varying heat-flux (TVHF) perturbation for different values of λ, V_a, and τ to investigate how TVHF perturbation affects ARTI growth; the other simulation conditions were basically the same as those for SHF perturbation.

Figure 3 shows the density contours of ARTI induced by TVHF [(a)–(e)] and SHF [(f)] perturbations with λ = 70 μm and V_a = 3.5 μm/ns at t ≈ 5.0 ns, where τ increases from Figs. 3(a)–3(e). As can be seen, the ARTI amplitude (i.e., the vertical height between the bubble vertex and the spike tip) in the SHF case grows larger than in all the TVHF cases, which indicates that the TVHF perturbation mitigates the ARTI growth for a given wavelength. The ARTI amplitude is smaller for shorter τ, except for the special case of τ = 0.2 ns [see Fig. 3(c)]. Figure 3(c) also shows a shell that is extraordinarily denser than in the other cases, indicating a higher shell compression in this case. Similar phenomena were also seen in the simulations for λ = 50 and 30 μm with V_a = 3.5 μm/ns, in which the special values of τ for TVHF perturbation were ∼0.15 and 0.10 ns, respectively.

Fig. 3. Density contours of ARTI induced by time-varying heat-flux (TVHF) perturbation for different τ with λ = 70 μm and V_a = 3.5 μm/ns at t ≈ 5.0 ns. The perturbation periods are (a) τ = 0.05 ns, (b) τ = 0.1 ns, (c) τ = 0.2 ns, (d) τ = 0.5 ns, and (e) τ = 1.0 ns, and (f) is for SHF perturbation.

下载图片 查看所有图片

The temporal evolution of ARTI bubble velocity induced by TVHF and SHF perturbations with λ = 70, 50, and 30 μm and V_a = 3.5 μm/ns is plotted in Figs. 4(a)–4(c), respectively. As can be seen, in each case U_b grows from the same initial value around V_a and saturates at the classical terminal velocity U_cl2D after the linear growth stage. U_b grows the earliest in the SHF cases, while the ARTI growth in all the TVHF cases is postponed for all λ. For given λ, ARTI grows later for shorter τ, and this is believed to be because shorter τ leads to more-substantial thermal smoothing and therefore smaller ARTI initial perturbation amplitude on the ablation front, which is consistent with the simulation results showing that SSD with higher modulated frequency leads to a lower rms level of illumination nonuniformity45 and a decreased amplitude of the initial perturbation of ARTI.62 The TVHF cases with the special values of τ exhibit extraordinarily faster linear growth, which is also shown in Fig. 5(a), and exhibit larger saturation value of U_b for all λ.

Fig. 4. Temporal evolution of bubble velocity for ARTI induced by TVHF and SHF perturbations for (a) λ = 70 μm, (b) λ = 50 μm, and (c) λ = 30 μm with V_a = 3.5 μm/ns.

下载图片 查看所有图片

Fig. 5. (a) Simulation (red triangles) and theoretical (black solid line) results for linear growth rate and (b) average effective acceleration (blue dots) and ablation pressure (red squares) in linear stage of ARTI induced by TVHF perturbation for different τ with λ = 70 μm and V_a = 3.5 μm/ns. (c) Profiles of average heat flux along z direction for three different configurations at t ≈ 2.0 ns.

下载图片 查看所有图片

To address the abnormal behavior near the special period, we performed more simulations with finer τ resolution and measured the linear growth rate for each case. Figure 5(a) compares the measured linear growth rates with those calculated theoretically by Eq. (2) for the TVHF cases with λ = 70 μm and V_a = 3.5 μm/ns. The simulation results agree well with the theoretical ones, and both peak at τ ≈ 0.2 ns. Figure 5(a) also shows that for τ larger or smaller than 0.2 ns, γ approaches that in the SHF case [see the λ = 70 μm case in Fig. 2(b)]. In the calculation using Eq. (2), the acceleration is taken as the average value in the linear stage for each case, which is plotted in Fig. 5(b). A substantial acceleration enhancement of more than 40% is observed near τ = 0.2 ns despite the fact that the same average heat flux on the bottom boundary is assigned for all cases, which may facilitate implosion efficiency with the same energy flux imported.

The enhanced acceleration near the special period is found to be due to the heat-flux enhancement in the ablated plasma, which increases the ablation pressure on the ablation front. Figure 5(b) shows that the ablation pressure reaches its maximum near the resonant period, which is closely correlated with the behavior of g. Figure 5(c) compares the heat-flux profiles of three simulations with different τ at the same time t ≈ 2 ns when ARTI has barely developed. In the case with the special value of τ = 0.2 ns, the heat flux is larger than that in the 1D steady ablation case (with no ARTI perturbation), while the case with τ far from the special value shows a heat flux that is very close to that in the 1D steady ablation case. Although the heat-flux differences in Fig. 5(c) look mild, the temperature and pressure at the ablation front can be influenced significantly by integrating Q = −κ∇T.

The behavior near the special period demonstrates a resonance mechanism in the system. We performed more simulations with finer τ resolution for different values of λ and V_a, and the linear growth rates are plotted in Fig. 6. To identify the key parameter that leads to resonance, we changed the plasma conditions (temperature, pressure, etc.) in the ablation region by altering V_a. As can be seen, the maximum linear growth rates are reached near the points where the phase velocity of the traveling-wave TVHF perturbation (i.e., V_ph = ω/k = λ/τ) matches the characteristic sound speed c̄s in the ablation region, where c̄s=∫ρcs2dV/∫ρdV is a density-weighted average of the local sound speed c_s in the ablation region. Similar phenomena are observed for all V_a and λ despite the fact that the average temperature and characteristic sound speed in the ablation region differ substantially for different values of V_a, as shown in Figs. 6(a)–6(c).

Fig. 6. Linear growth rate of ARTI vs phase velocity of TVHF perturbation for different λ and V_a: (a) V_a = 2.0 μm/ns; (b) V_a = 3.5 μm/ns; (c) V_a = 5.0 μm/ns. Pink circles: λ = 30 μm; blue triangles: λ = 50 μm; red squares: λ = 70 μm. The red dashed line represents the corresponding characteristic sound speed in the ablation region for each V_a.

下载图片 查看所有图片

The term c̄s characterizes the traveling speed of the density fluctuations, while V_ph determines that of the temperature fluctuations along the x direction in the ablation region; c̄s≈Vph is then expected to cause a resonance between the density and temperature fluctuations. In the TVHF regime, we plot the relative phase (δϕ) between the two types of fluctuations vs V_ph with λ = 70 μm and V_a = 3.5 μm/ns in Fig. 7. Here, δϕ̄ is defined as the time average ofδϕ(t)=∫ρ[2π|xρmax(z,t)−xTmax(z,t)|/λ]dz∫ρdz(5)in the first 2 ns of the simulations, where x_ρmax and x_Tmax are the x locations of the maxima of the density and temperature fluctuations, respectively. As can be seen, the minimum δϕ̄<0.2π is reached when c̄s≈Vph. In the SHF regime, the x location of the maximum of the density fluctuation is expected to overlap with the minimum of the temperature fluctuation under the isobaric approximation along the x direction, which yields δϕ̄=π. The resonance changes the phase considerably by bringing the maxima of the density and temperature fluctuations closer and thus leads to more-efficient heating of more mass in the ablation region, which causes the heat-flux enhancement shown in Fig. 5(c).

Fig. 7. Average phase difference of density and temperature fluctuations traveling along x direction in ablation region with λ = 70 μm and V_a = 3.5 μm/ns during t = 0–2 ns.

下载图片 查看所有图片

5 V. CONCLUSION

In this study, the evolution of 2D ARTI induced by single-mode stationary and time-varying heat-flux perturbations was investigated via simulations. In the SHF cases, the shorter-wavelength modes grew substantially later than the longer-wavelength modes, which was reasoned to be due to the enhanced transverse thermal smoothing caused by their finer structures and the mitigated actual perturbation levels at the interface. The linear growth rate was barely affected by the continuously applied heat-flux perturbation on the bottom boundary and could still be predicted well by linear theory.17

The TVHF perturbation mitigated the ARTI growth because of the enhanced thermal smoothing induced by the wave-like heat-flux source. It was found that ARTI growth was postponed for shorter τ with given λ. The TVHF perturbation near the resonant period for which Vph≈c̄s in the ablation region minimized the relative phase between the density and temperature fluctuations by bringing their maxima closer, thus leading to more-efficient heating of more mass. Consequently, this resonance effect offers significantly enhanced electron thermal conduction and thus larger ablation pressure and effective acceleration, which consequently yield higher linear growth rate and saturated bubble velocity in ARTI growth.

The enhancement of ablation pressure and effective acceleration caused by near-resonance TVHF perturbation can facilitate increased implosion velocity but can also compromise the integrity of ICF shells by causing faster ARTI growth in both the linear and nonlinear phases. The positive and negative influences of this resonance effect deserve attention and further study in ICF-relevant research, and important topics should be explored in the future, such as ARTI evolution in the presence of both interfacial and heat-flux perturbations.

6 ACKNOWLEDGMENTS

Acknowledgment. This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant Nos. XDA25050400 and XDA25010200), by the National Natural Science Foundation of China (Grant No. 12175229 and 12388101), by the Frontier Scientific Research Program of the Deep Space Exploration Laboratory (Grant No. 2022-QYKYJH-HXYF-019), and by the Fundamental Research Funds for the Central Universities. The numerical calculations reported herein were performed on the supercomputing system in the Supercomputing Center of the University of Science and Technology of China.