A standing wave induced in front of the ionization-front of a millimeter-wave discharge was numerically investigated to develop an interferometric discharge structure identification method. The time-varying waveform of the standing-wave intensity obtained at a distant observation point was smooth when a continuous comb-shaped structure was formed, whereas it was noisy with high-frequency components when a discrete structure was formed. The peak frequency of the Fourier spectrum of the time-varying waveform was proportional to the ionization-front propagation speed. The rapid time-variation of the waveform was caused by an increase in millimeter-wave absorption in a new plasma spot formation in the discrete structure. The results suggest that discharge structure identification, measurement of ionization-front propagation, and timing of plasma spot formation can be conducted experimentally without using a high-speed camera.

## I. INTRODUCTION

Atmospheric discharge induced by high-power millimeter-wave beam irradiation has been actively investigated since the 1980s, along with the development of gyrotrons for nuclear fusion research.^{1–21} Discharge experiments have shown that an ionization-front, which is the region where the background gas ionizes, propagates toward the beam source. The propagation speed and shape of the trajectory of the ionization-front, that is, the discharge structure, depend on the beam power density of the millimeter-wave, millimeter-wave frequency, pressure, and background gas species.^{11–20} Past experiments have reported some discharge structures, such as $\lambda /4$,^{6–10} comb-shaped,^{14–18} tree-branching,^{8,20} diffusive,^{19,22} and complex filamentary structures.^{19,22} In the $\lambda /4$-structure, the ionization-front propagates discretely. Millimeter-wave discharge is classified into overcritical and subcritical conditions, depending on the local beam power density. The $\lambda /4$-structure is observed in the overcritical condition, whereas the comb-shaped, diffusive, and complex filamentary structures are observed in the subcritical condition. The overcritical condition is that the electron production by electron-impact ionization exceeds the loss of electrons by the attachment process owing to the high electric-field intensity. The propagation speed was generally greater than 10 km/s.^{6,7} The balancing condition was defined as the critical intensity $ E cr$, which was $2.45 MV / m$ in the case of 170 GHz at atmospheric pressure. Under subcritical conditions, the beam-power density is lower than $ E cr$, and long-time-scale phenomena appear, such as neutral-gas heating and chemical reactions.^{4} The propagation speed is in the range of tens to thousands of meters per second.^{16,18,19}

The millimeter-wave discharge of the 170 GHz band under subcritical conditions has been investigated using the gyrotron developed for ITER.^{14–17} This study also focuses on the 170 GHz band because abundant experimental data and past numerical simulation models are available. A discharge experiment using a 170 GHz Gaussian beam has reported that a comb-shaped structure is formed by granular plasma spots spreading in a diagonal direction.^{16}^{,} Figures 1(a) and 1(b) show schematic sketches of the granular plasma spots drawn based on snapshots captured in previous experiments.^{16} In fact, the plasma spots in the experiment are more abundant than those shown in the figure. Figure 1(a) shows that the plasma spots arranged in concentric circles are observed from the viewpoint of the beam source. In the side view of Fig. 1(b), the plasma spots are arranged in a bow-shaped arrangement with even spacing. Thus, the central axis region is hidden by the peripheral plasma spots from a side viewpoint. Moreover, the chronophotography of the high-speed camera shows that the propagation speed in the diagonal direction was less than 400 m/s, whereas that toward the beam source was approximately 1000 m/s.^{16} The order of the difference in propagation speed can be caused by a different mechanism of ionization-front formation between the central axis region and the peripheral region. In fact, it has not been observed whether granular plasma spots are formed in the central region because the luminescence of plasma spots arranged in concentric circles shields the emission from the central region of the comb-shaped structure. Moreover, the Gaussian beam provides a higher electric-field region in the central axis region and may cause changes in the ionization mechanism. In addition, an experiment using a flat-top beam has been conducted, and a comb-shaped structure in which the plasma spots are arranged in a row was reported.^{21} This difference in the structures indicates that a beam profile with a strong intensity at the central axis is important for forming a comb-shaped structure spreading diagonally.

As mentioned above, it is difficult to elucidate millimeter-wave discharge phenomena only through discharge experiments, owing to its supersonic propagation, microscopic millimeter-scale structure, and large gradients of plasma density. Therefore, numerical studies have been conducted actively. In the overcritical conditions, the numerical model coupling electromagnetic-wave propagation and reaction–diffusion of the plasma have simulated the propagation speed and $\lambda /4$-structure of the overcritical experiments.^{23–30} Past numerical simulations have been conducted in the 170 GHz band as a target under subcritical conditions because of abundant experimental data. However, the propagation speed of the comb-shaped structure observed at 170 GHz was not reproduced in the conventional numerical models. Therefore, it is necessary to modify the numerical model to design engineering applications, such as the beaming propulsion rocket.^{21,31–42}

Our research group conducted numerical simulations to investigate the formation mechanism of the 170 GHz millimeter-wave discharge under subcritical conditions.^{33,40} We proposed the combination of an overcritical model^{23–30} with compressible neutral fluid, detailed chemical reaction, and radiation transport modules to simulate subcritical discharge.^{40} The 1D simulation aimed to reproduce the one-dimensional discharge front propagation along the central axis; however, it cannot simulate the electric-field enhancement by multi-dimensional millimeter-wave interference. The simulation results using this model reproduced three different discharge structures: diffusive, sharp, and discrete, in the order of small electric-field intensity. The discrete structure is formed by electron-impact ionization in the anti-node of the standing wave created by the interference of the incident millimeter-wave and the millimeter-wave reflected by the ionization-front. The sharp structure is formed by electron-impact ionization assisted by a high-temperature neutral gas, which is heated through a fast gas-heating process. The diffusive structure was formed by thermal and associative ionization in a high-temperature neutral gas heated through a vibrational–translational relaxation process. The sharp structure exhibited large electron density gradients, indicating that the continuous plasma spots propagation in the experiment corresponds to the sharp structure in the 1D simulation. The sharp structure was obtained in the numerical simulation in the range of experimental electric-field intensity in which the comb-shaped structure was obtained in the past experiments (0.4–1.2 MV/m). The propagation speed was 1000 m/s along with the central axis in the experiment but was 100–200 m/s in the simulation, which was the same order as the diagonal plasma spots propagation in the experiment at the intensity range of 0.4–1.2 MV/m. In the higher electric-field intensity greater than 1.4 MV/m, a discrete structure was formed, and its propagation speed was of the order of 1000 m/s, which agreed with the propagation speed at the central region of the comb-shaped structure. The propagation speed of sharp and diffusive structures was subsonic because the propagation speed could not exceed the propagation speed of the shock waves; however, the propagation speed of the discrete structure was supersonic because the jumping of the ionization-front overtook the shock wave. These findings have suggested that the central region of the comb-shaped structure may propagate discretely; however, this cannot be concluded only using 1D simulations.

Based on these experimental and numerical findings, we propose a hypothesis for the composite discharge structure. In general, the comb-shaped structure in the 170 GHz experiment has been considered as granular plasma spots that spread continuously in the radial direction, as shown in Fig. 2(a). However, the propagation speed along the central axis in the experiment is of the same order as that of the discrete propagation reproduced in the numerical simulation. Therefore, it is consistent to consider that the comb-shaped structure propagating at supersonic speed can be a composite structure of discrete propagation in the central region and continuous propagation in the periphery region, as shown in Fig. 2(b). However, observations using high-speed cameras have not confirmed whether a composite structure exists for now because the emission of the peripheral plasma spot hides the central region.^{16} Thus, it is necessary to develop a fast observation method that does not rely on an optical camera to identify the discharge structures.

^{43}At the Prokhorov General Physics Institute, Russia Academy of Science (GPI-RAS), a 75 GHz discharge experiment has been conducted, and a comb-shaped structure, similar to the 170 GHz experiment, was also formed. They used an interferometer to measure the ionization-front propagation speed; this is known as the location method.

^{13}However, they did not consider the transition of discharge structures and have assumed that the ionization-front propagates continuously, and the below equation can represent the ionization-front propagation speed $ u ion$ because interference of incident and reflected millimeter waves is enhanced when the ionization-front propagates over a distance of $\lambda /2$,

Therefore, we propose a new interferometric discharge structure identification (IDSI) method using a time-varying waveform of the standing-wave intensity, which is formed by the interference of the incident millimeter-wave and the millimeter-wave reflected by the ionization-front. A one-dimensional schematic diagram of the proposed IDSI method is shown in Fig. 3. A high-power millimeter-wave is irradiated from the left and propagates toward the right in Fig. 3. A standing wave is formed in front of the ionization-front because part of the incident millimeter-wave is reflected in the ionization-front and overlaps with the incident millimeter-wave. As the ionization-front propagates, the standing wave also moves with the ionization-front. Therefore, the discharge structure can be identified by measuring the temporal variation in the standing-wave intensity at an observation point. A smooth time-variation waveform can be obtained when the ionization-front propagates continuously. The previous research considers this condition to measure the propagation speed and assumes that the ionization-front propagates over a distance of $\lambda /2$ in a period of the enhancement of interference. By contrast, a noisy time-variation waveform containing high-frequency components can be obtained when the ionization-front propagates intermittently or discretely. Moreover, it is unclear whether the propagation speed of the ionization-front can be measured as the same as the continuous propagation. Therefore, this study confirms whether the proposed IDSI method can identify discharge structures and measure the propagation speed using two-dimensional numerical simulations before conducting discharge experiments. The numerical simulation results can provide the relationships between the microscopic structure of the plasma spots and the time-varying waveform of the standing-wave intensity at the observation points. The IDSI method can be applied to reconstruct a three-dimensional discharge structure using multi-point measuring without high-speed cameras.

Although the optical system of the IDSI method is one among the special forms of homodyne reflectometry,^{43,44} the biggest difference between them is the usage of a high-power millimeter wave generated by a gyrotron. The IDSI method can integrate the mixer and the receiving antenna of ordinary reflectometry by using a high-power millimeter-wave, simplifying the millimeter-wave discharge experiment system. However, if this method is applied to other plasma discharges such as ExB devices, the advantage of non-intrusiveness to its dynamics can be lost because the high-power millimeter wave for diagnosis promotes ionization, thereby changing the physical properties. For such purposes, typical reflectometry using a low-power microwave is suitable. Moreover, although microwave reflectometry generally aims to measure the electron number density distribution or its fluctuation,^{43,44} the IDSI method is specialized in the structural identification of the ionization-front and the measurement of its propagation speed. This is because it is applicable when the electron number density is greater than the cut-off density, the electron number density gradients are very steep, and the spatial scale is smaller than the wavelength of the incident beam. The millimeter-wave discharge can meet this condition and the plasma front can be treated as a fixed end of reflection, making the IDSI method applicable.

Although this work focuses on the reflected millimeter wave and plasma structure in the subcritical conditions, previous studies revealed that wave refraction, scattering, and diffraction occur during the discharge process, in addition to reflections. Concretely, Batanov *et al.* reported the refracted millimeter-wave by the filamentary plasma and EUV haro region in the subcritical discharge experiment,^{45} inducing reionization plasma spots behind the ionization-front. However, in this study, the refracted millimeter-wave does not affect the measurement by the IDSI method because the observation points are set on the front side of the ionization-front. Additionally, Cook *et al.* conducted experimental research to investigate the scattering and diffraction by the millimeter-wave discharge plasma in the overcritical condition.^{46} They reported the reflected millimeter-wave signal toward the beam source to be a high-frequency time-varying waveform. The reason was stated in their study, but this waveform may be generated by the intermittent propagation of the ionization-front, according to the same principle of the IDSI method.

The purpose of this study is to investigate how ionization-front propagation affects the time-variation of standing-wave intensity at a distant point. A numerical simulation is used to obtain detailed information on the ionization-front propagation and the time-varying waveform of the standing-wave intensity. Using this information, we verify (1) whether the discharge structure can be identified, (2) whether the propagation speed of the ionization-front can be measured, and (3) whether the discrete propagation of the plasma spots can be observed using the standing-wave intensity.

## II. NUMERICAL MODEL

We performed 2D simulations using a numerical model coupling electromagnetic-wave propagation, reaction–diffusion of plasma, and compressible neutral fluid modules. This model was proposed by Takahashi *et al.* in which the transport coefficients were calculated using the Boltzmann equation solver, Bolsig+.^{33} Bolsig+ is a free software developed by Laboratoire Plasma et Conversion d’Energie (LAPLACE).^{47} However, our model uses approximate formulas to calculate the transport coefficients that have been used in the previous overcritical model^{24,25} to reduce computation time. The comb-shaped and discrete structures were reproduced in our model, as in Takahashi *et al.* under subcritical conditions.^{33} In the overcritical condition, the propagation speed and the discharge structure of our model agree with those of previous studies^{23–25} if the compressible fluid module is off.

^{48}Equation (4) is numerically integrated using an explicit Euler method.

^{23}The number densities of the electrons $ n e$, the positive ions $ n p$, and the negative ions $ n n$ were calculated by solving the following reaction–diffusion equation:

^{23}The electron density flux $ \Gamma e$ is obtained by the following equation:

^{23}

^{,}$\zeta = \nu i \tau m= \lambda D 2/ L 2$ is the transition parameter. Here, $ \nu i$ is the ionization frequency, $ \lambda D$ is the Debye length, and $L$ is the scale length of the electron distribution. $ \tau m= \epsilon 0/[e n e( \mu i+ \mu e)]$ is the Maxwell relaxation time, $ \mu i$ and $ \mu e$ are the ion and electron mobility, respectively. $ D e$ is the free diffusion coefficient of the electron, which is calculated by $ D e= \mu ek T e/e$, $k$ is the Boltzmann constant, and $ T e$ is the electron temperature, which is assumed to be 1.5 eV in this study.

^{24}

^{,}$ D a$ is the ambipolar diffusion coefficient. The electron mobility is given by $ \mu e=e/ m e \nu m$, and the ion mobility is assumed to be $ \mu e/200$ in this study.

^{23}The source terms $S= S e, S p$, and $ S n$ are obtained using the following equations:

^{40}Reaction rates were calculated using the following approximation:

^{25}

^{33}The time integration of Eq. (15) was conducted using the explicit Euler method and spatially discretized using the cell-centered finite volume method. The AUSM-DV scheme and the third-order MUSCL interpolation were used to evaluate the numerical flux.

^{49,50}

## III. CONDITION

^{14,16}The wavelength $\lambda $ is 1.76 mm. The incident RMS electric-field intensity $ E 0 , rms$, is set as 1.4, 1.6, 1.8, 2.0, 2.1, 2.12, 2.15, and 2.2 MV/m, which corresponds to the millimeter-wave power densities of 5.20, 6.79, 8.59 10.6, 11.7, 11.9, 12.3, and 12.8 $ GW / m 2$. Figure 4 shows the 2D calculation domain used in this study. The computational domain has $5\lambda \xd75\lambda $, and the grid size $\Delta x$ and $\Delta y$ are $\Delta x=\Delta y=\lambda /N$. Here, $N=100$ was used as in a previous study.

^{33}The millimeter-wave beam was irradiated from the left boundary; therefore, the ionization-front propagates toward the left direction. The incident millimeter-wave was set as a plane wave of the TM mode. The TM mode was chosen to reduce the computational time because the discharge structure in the TE mode requires a 1/10 smaller grid size compared with the TM mode.

^{24}However, discrete and continuous structures are formed either in the TE or TM mode, and the standing wave is expected to be formed in both cases. Therefore, the relationship between the time variation of the standing-wave intensity and the discharge structure can be investigated in both modes. Moreover, the experimental discharge structures were reproduced using the plane wave in the past numerical simulations

^{23,33,36}even though the Gaussian beam was used in the experiments.

^{6,16}This is because the discharge occurs at the central beam axis where the beam intensity is high and its spatial variation is sufficiently small to be treated as a plane wave for the plasma. As similar to the previous simulations, in this study, for performing a simple evaluation of our method, the plane wave has been used by assuming that the plasma spot is sufficiently smaller than the beam radius. However, it is expected that the IDSI method will be applicable when a more practical beam profile such as the Gaussian shape is used in the experiment because the phase of the reflected millimeter wave does not change regardless of the beam profile. The initial electron and expansion region of the neutral fluid is set as Gaussian distribution at $(x/\lambda ,y/\lambda )=(3.75,2.5)$ to ignite discharge at the subcritical condition. The full width at half maximum of the Gaussian distribution was 0.176 mm. The peak electron number density is $1.0\xd7 10 15 m \u2212 3$ and the peak expansion ratio of neutral gas is 33.3, which is the same as in previous research.

^{33}The initial electron number density is considerably smaller than the cut-off density, of the order of $ 10 21$– $ 10 22 m \u2212 3$, to avoid the initial electron effects on the discharge structure. The time step of the plasma and the neutral fluid solver $\Delta T fluid$ was calculated as below:

## IV. RESULTS

### A. Change of standing wave depending on discharge structures

Figures 5(a) and 5(b) show the distributions of the electron number density $ n e$ at $ E 0 , rms=1.4$ and $2.2 MV / m$, respectively. A comb-shaped structure was formed at 1.4 MV/m, as shown in Fig. 5(a), because of the neutral gas expansion through Joule heating, as reported in a previous work.^{33} However, a discrete discharge structure was formed at 2.2 MV/m, as shown in Fig. 5(b). A high electric-field intensity at the anti-node of the standing wave formed by millimeter-wave interference induces stepwise propagation of the ionization-front. This structure was the same as that in a previous work.^{33} Comb-shaped structures are induced below 1.8 $ MV / m$ and the discrete structures were formed at least 2.0 $ MV / m$.

Figures 6(a) and 6(b) are $x$– $t$ diagram of the electron number density $ n e$ along $y/\lambda =2.5$ at $ E 0 , rms=1.4$ and 2.2 $ MV / m$, respectively. The black lines in Figs. 6(a) and 6(b) indicate the location of the ionization-front, which is determined as the position where the electron density exceeds $1\xd7 10 9$ for the comb-shaped structure and $1\xd7 10 13 cm \u2212 3$ for the discrete structure. Figure 6(a) shows that the ionization-front propagates continuously toward the beam-source direction ( $\u2212x$ direction). However, the ionization-front of the discrete structure propagates as a leapfrog, as shown in Fig. 6(b).

Figures 7(a) and 7(b) are the $x$– $t$ diagrams of the normalized RMS value of electric-field intensity $ E rms/ E 0 , rms$ along $y/\lambda =2.5$ at $ E 0 , rms=1.4$ and $2.2 MV / m$, respectively. Figures 7(a) and 7(b) show that standing waves were formed in front of the ionization-front. They are induced by the interference of incident millimeter waves and reflected millimeter waves from the ionization-front. Here, the reflection point moves with the ionization-front propagation. The nodes and anti-nodes of standing wave move smoothly along with the continuous ionization-front propagation in the case of forming the comb-shaped structure as shown in Fig. 7(a). However, the nodes and anti-nodes of standing wave move intermittently in case of forming the discrete structure, as shown in Fig. 7(b), because of the stepwise ionization-front propagation. This finding suggests that the time-variation of the standing-wave intensity includes information on the ionization-front propagation and discharge structure.

The relationship between the time-varying waveform of the standing-wave and the discharge structure was then examined to distinguish between the comb-shaped and discrete structures. Figure 8(a) shows a time-variation of the standing-wave intensity $ E rms/ E 0 , rms$ at the observation position of $(x/\lambda ,y/\lambda )=(0.01,2.5)$ on the central axis at $ E 0 , rms=1.4 MV / m$. This figure shows that the waveform was smooth and did not contain high-frequency components. Figure 8(b) is a Fourier spectrum of the waveform in the case of $ E 0 , rms=1.4 MV / m$, which has one strong peak of 0.43 MHz without high-frequency components at more than the peak frequency.

Figure 9(a) shows a time-variation of the standing-wave intensity $ E rms/ E 0 , rms$ at the position of $(x/\lambda ,y/\lambda )=(0.01,2.5)$ in the case of $ E 0 , rms=2.2 MV / m$. Figure 9(b) is a Fourier spectrum of the time-varying waveform at $ E 0 , rms=2.2 MV / m$, which has one strong peak of 2.62 MHz with high-frequency components greater than the peak frequency. The high-frequency components originated from the stepwise propagation of the ionization-front because they were not formed in the case of the comb-shaped structure, as shown in Fig. 8. The detailed mechanism of the high-frequency components induced in the discrete structure is discussed in Sec. IV C. By contrast, we can conclude that the discrete propagation occurs if a high-frequency waveform is obtained. This IDSI method can be applied to millimeter-wave discharge experiments to distinguish discharge structures without a high-speed camera.

### B. Waveform and ionization-front propagation speed

In the previous experiments,^{13} the propagation speed was measured using the peak frequency of the interference wave when a comb-shaped structure was formed. However, it is unclear if this location method can be applied when a discrete structure is induced. In this section, we discuss whether the method can measure the propagation speed of the ionization-front in this case by investigating the relationship between the Fourier spectrum with a strong peak and the ionization-front propagation speed.

First, we focused on whether the location method could be applied to the numerical simulation when a continuous structure was induced. Previous experiments have reported that the peak frequency of the waveform is proportional to the propagation speed of the ionization-front.^{13} Therefore, if the location method is applicable to a numerical simulation, it is expected that similar relationships will be obtained. Although the time-varying waveform was a mixed signal of incident and reflected millimeter waves in the previous experiment, this waveform corresponds to the temporal evolution of the rms electric-field intensity $ E rms$ for the standing wave, which was captured at a specific observation point in the numerical simulation. Here, the standing-wave intensity was normalized by the incident electric-field intensity of the millimeter-wave $ E 0 , rms$ to compare the case of the discrete structure.

E_{0,rms} (MV/m)
. | u_{ion,m} (m/s)
. | f_{p,m} (MHz)
. | f_{p,th} (MHz)
. |
---|---|---|---|

1.4 | 374 | 0.43 | 0.424 |

1.6 | 416 | 0.47 | 0.472 |

1.8 | 426 | 0.53 | 0.483 |

2 | 1353 | 1.51 | 1.53 |

2.1 | 1704 | 1.97 | 1.93 |

2.12 | 1745 | 2.02 | 1.98 |

2.15 | 1836 | 2.22 | 2.08 |

2.2 | 2213 | 2.62 | 2.51 |

E_{0,rms} (MV/m)
. | u_{ion,m} (m/s)
. | f_{p,m} (MHz)
. | f_{p,th} (MHz)
. |
---|---|---|---|

1.4 | 374 | 0.43 | 0.424 |

1.6 | 416 | 0.47 | 0.472 |

1.8 | 426 | 0.53 | 0.483 |

2 | 1353 | 1.51 | 1.53 |

2.1 | 1704 | 1.97 | 1.93 |

2.12 | 1745 | 2.02 | 1.98 |

2.15 | 1836 | 2.22 | 2.08 |

2.2 | 2213 | 2.62 | 2.51 |

The relationship between the $ f p , th$ obtained by Eq. (23) and $ f p , m$ obtained by the numerical simulation in the discrete structure was investigated in the same manner as in the continuous structure. The discrete structure was formed at $ E 0 , rms=2.0,2.1,2.12,2.15,2,2 MV / m$. Table I shows that the structural transition from 1.8 to 2.0 $ MV / m$ leads to a sudden increase in the ionization-front propagation speed $ u ion , m$. The $ f p , m$ are plotted against the $ u ion , m$ in Fig. 10, which shows a proportional relationship, and the $ f p , m$ is in agreement with the theoretical value $ f p , th$ as the same as in the case of continuous propagation at less than 1.8 MV/m. These findings suggest that the peak frequency in the Fourier spectrum of the standing wave is induced by the ionization-front propagation in both comb-shaped and discrete structures.

Therefore, our numerical simulation suggests that the ionization-front propagation speed can be measured using the peak frequency in the frequency spectrum of the standing-wave intensity observed at a specific point even if a structural change of the millimeter-wave discharge occurs.

### C. Waveform in discrete structure

In Sec. IV B, the ionization-front propagation induced a low-frequency component of the time-varying waveform of the standing-wave intensity when comb-shaped and discrete structures were induced. However, the high-frequency component is formed only when a discrete structure is formed. The physical significance of the high-frequency component is discussed in this section.

A periodic waveform was formed when a discrete structure was induced as shown in Fig. 9(a). Figure 11 is a close-up view between $t=1.95$ and $t=2.45 \mu s$ of Fig. 9(a). The waveform has six local maximum and local minimum points: local maximum (A1), local minimum (A2), local maximum (A3), local minimum (N1), local maximum (N2), and local minimum (N3) points in sequence. In this section, we numbered the points from A1 to A3 for the anti-node and, additionally, N1 to N3 for the node, as shown in Fig. 11, to investigate the relationship between the waveform and plasma spots of the discretely ionization-front propagation. Point 0 in Fig. 11 is the last local minimum point in the previous cycle.

Figure 12 shows the distribution of the electron number density $ n e$ and the distribution of the normalized local $ E rms$ for the anti-node at each time when the local maximum and minimum points are formed, respectively.

First, the initial state (point 0) is examined before the local maximum of point A1 is formed. A central plasma spot is formed at $(x/\lambda ,y/\lambda )=(1.81,2.5)$ as shown in Fig. 12(a). At this time, Fig. 12(b) shows that a high electric-field intensity region is formed at diagonal positions at $(x/\lambda ,y/\lambda )=(1.75,2.06)$ and $(1.75,2.93)$ in front of the central plasma spot at $(x/\lambda ,y/\lambda )=(1.81,2.5)$ due to a focusing of the reflected millimeter wave. Figure 12(b) also shows that the node of the standing wave is formed at the observation point of the waveform at $(x/\lambda ,y/\lambda )=(0.01,2.5)$. Thus, the waveform exhibits a local minimum (initial state point 0), as shown in Fig. 11. Here, the incident millimeter-wave is completely cut off at $x/\lambda =2$ as shown in Fig. 12(b) because of a high number density of electron.

Subsequently, the local maximum of point A1 is formed at $t=2.14 \mu s$ in Fig. 11. At this time, Fig. 12(c) shows that two diagonal plasma spots are formed at $(x/\lambda ,y/\lambda )=(1.75,2.07)$ and $(1.75,2.93)$ where the high electric-field intensity regions have been formed in Fig. 12(b). This is because of a time lag between the high electric-field region formation and the plasma spot formation. Comparing Figs. 12(b) with 12(d), the node of the standing wave changes to the anti-node at the observation point of $(x/\lambda ,y/\lambda )=(0.01,2.5)$ because the front of the cut-off region moved $\lambda /4$ toward beam source ( $x/\lambda =0$) along the $x$ axis.

At $t=2.16 \mu s$ of the local minimum of point A2, the standing-wave intensity decreases temporarily. The reason is that as the new central plasma spot is being formed at $(x/\lambda ,y/\lambda )=(1.55,2.5)$, as shown in Fig. 12(e), the reflected millimeter-wave intensity is reduced because the millimeter-wave absorption is increased in the new plasma spot where the electron number density is less than the cut-off density $ n c$. Moreover, Fig. 12(f) indicates the electric-field intensity is still enhanced at $(x/\lambda ,y/\lambda )=(1.55,2.5)$ where the new central plasma spot is being formed, which indicates that the millimeter-energy is absorbed effectively. Comparing Figs. 12(d) with 12(f), the anti-node is maintained at the observation point in Fig. 12(f) because the distribution of the nodes and anti-nodes does not change. However, the 1D distribution of $ E rms/ E 0 , rms$ along the central axis (Fig. 13) shows that the electric-field intensity of the standing wave at $(x/\lambda ,y/\lambda )=(0.01,2.5)$ is decreased at $2.16 \mu s$ (A2) compared with $2.14 \mu s$ (A1) because of the absorption by the central plasma spot during the formation.

In contrast, the standing-wave intensity takes the local maximum value again at $t=2.21 \mu s$ as point A3 in Fig. 11. A new central plasma spot has been completely formed at $(x/\lambda ,y/\lambda )=(1.55,2.5)$ and reaches the cut-off density at this time as shown in Figs. 12(g) and 12(h). Thus, the absorption of the millimeter wave by the central plasma spot is smaller than the time when the local minimum of A2 is formed, which leads to an increase in the electric-field intensity at the observation point compared to that at point A2. Additionally, a local maximum of the waveform is formed because the diagonal plasma spots, which are larger than the central plasma spot, contribute to the formation of an anti-node at the observation points, as shown in Fig. 12(h). Figure 13 shows a 1D distribution of the electric-field intensity $ E rms$ along the central axis. At 2.14 and 2.18 $ \mu s$, the standing intensity is decreased to $x/\lambda =0$ from $x/\lambda =1.6$, which indicates that the reflected millimeter wave is diverging. However, the standing-wave intensity is constant between $x/\lambda =0$ to $1$ at $2.21 \mu s$, which suggests that the diagonal plasma spots arrangement shown in Fig. 12(g) contributes to the focus of reflected millimeter-wave toward the observation point. These processes increase the standing-wave intensity at $2.21 \mu s$ and induce the local maximum A3. Overall, the small and rapid variations in the standing-wave intensity were caused by the time lag between the formation of the central and diagonal plasma spots.

Subsequently, the standing-wave intensity drops rapidly at $t=2.32 \mu s$, as shown in Fig. 11 (N1). At this time, a node is formed at the observation points as shown in Fig. 14(b). Figure 14(a) shows that two diagonal plasma spots are formed at $(x/\lambda ,y/\lambda )=(1.52,2.07)$ and $(1.52,2.93)$. The arrangement of the plasma spots shown in Fig. 14(a) is the same as that in Fig. 12(c), which has the diagonal plasma spots in the head region when the local maximum of point N1 is formed. The reason why the local maximum and the local minimum are replaced is that the node and the anti-node are replaced due to the head of the ionization-front propagating $\lambda /4$ even though the same plasma spots arrangement is formed. Therefore, the local minimum and maximum of points N1–N3 were formed by the same mechanism as those of points A1–A3. The only difference is whether an anti-node or a node is observed at the observation point.

In this numerical simulation, an ideal plasma spot formation, such as staggered propagation, was considered. However, the high-frequency time-variation of the waveform synchronized with a new plasma spot formation can occur in any plasma spots arrangement such as the complex filamentary structure observed in the 28 GHz discharge,^{19} because millimeter-wave absorption is increased when a new plasma spot is formed. Therefore, the timing of plasma spot formation can be experimentally measured by investigating the high-frequency signal superimposed on the lower-frequency waveform formed by the ionization-front propagation. The time resolution for measuring the standing-wave intensity, which is the same as the electric-field intensity, is of the order of GHz using an oscilloscope, which is much higher than the time resolution of a high-speed camera (MHz). Thus, the plasma spot formation timing obtained using a high-speed camera can be complemented by measuring the standing-wave intensity using an antenna.

### D. Off-axis observation point

In Secs. IV A–IV C, the observation point was set on the axis at $(x/\lambda ,y/\lambda )=(0.01,2.5)$. However, the antenna and rectifier cannot be installed on the axis in the actual experiment because the millimeter-wave beam propagating on the axis has a sufficiently high intensity to destroy the rectifier. Since the incident beam power density decreases exponentially as it departs from the beam axis because a Gaussian beam is used for the experiment, not a plane wave, the antenna and rectifier are expected to operate without damage due to the weakened incident wave if they are set off-axis. Therefore, for conducting the experiment, it is necessary to investigate the waveform change when the observation points are set off-axis. Here, although our simulation model assumes the plane wave injection for simplification, the waveform change depending on the observation locations can be examined because the phase of the reflected millimeter-wave is not changed depending on the beam profile.

The observation points of $(x/\lambda ,y/\lambda )=(0.01,4.5)$, (0.01,4.99), (0.5,4.99), and (1,4.99) are selected as the off-axis location as shown in Fig. 15. Figures 16(a)–16(d) show the time-varying waveform of the standing-wave intensity observed at each observation point including the on-axis point ( $y/\lambda =2.5$) for comparison. Figures 16(a) and 16(b) show the case of the comb-shaped structure induced at $ E 0 , rms=1.4 MV / m$. Figures 16(c) and 16(d) show the discrete structure induced at $ E 0 , rms=2.2 MV / m$. Smooth waveforms were obtained when the comb-shaped structure was induced, while noisy waveforms were obtained when the discrete structure was induced in at each observation point, which was the same as when the observation point was set at $(x/\lambda ,y/\lambda )=(0.01,2.5)$. Therefore, the IDSI method is useful even if the observation point is not set on the axis. However, the timing when the node and anti-node are observed is changed depending on the observation points as shown in Figs. 16(a)–16(d). This is because the standing-wave distribution is arcuate-shaped due to an effect of the plasma spot arrangement as shown in Fig. 15. As a specific example, the waveform observed at (0.01,4.99) has a local minimum at $t=11.58 \mu s$, whereas the waveform observed at (0.01,4.5) has a local maximum near $t=11.58 \mu s$ as shown in Fig. 16(a). The reason can be clearly explained by investigating the standing-wave intensity distribution because Fig. 15 shows that the node is located at (0.01,4.99), but the anti-node is located at (0.01,4.5) due to the arc-shaped standing-wave distribution. This result suggests that the arrangement of the plasma spot can be observed by conducting multi-spot measurements of the standing-wave intensity.

Figures 17(a)–17(d) show the Fourier spectrum of the waveform observed at $(x/\lambda ,y/\lambda )=(0.01,4.5)$, (0.01,4.99), (0.5,4.99), and (1,4.99), respectively, at a beam intensity of $ E 0 , rms=1.4 MV / m$. The general shapes of the spectrum are almost the same, but the peak frequencies in Figs. 17(a)–17(d) are 0.4, 0.4, 0.34, and 0.34 MHz, respectively, which decreases as the observation point deviates from the on-axis or approaches the ionization-front. Here, the peak frequency is 0.43 MHz when the observation point is set on-axis, as shown in Fig. 8. Such a decrease in the peak frequency at the off-axis observation points can be explained using Fig. 15. The time interval of the local maximums of the time-varying waveform in Fig. 16 corresponds to the distance between anti-nodes, which is equal to $\lambda /2$ when the observation point is on the axis as shown in Fig. 15. However, the distance between anti-nodes along the $x$ axis observed at the off-axis point is greater than $\lambda /2$ as shown using red arrows in Fig. 15. This is due to the arc-shaped standing-wave distribution generated by the plasma spot arrangement. Therefore, in actual experiments, depending on the observation point, the peak frequency is expected to decrease compared with the frequency calculated by Eq. (23) using the actual ionization-front propagation speed. It is suitable to use a beam splitter on the beam axis to measure the propagation speed, similar to the previous study,^{13} but the arrangement of the plasma spot can be obtained using the off-axis observation points. Moreover, the off-axis observation point may be applicable for measuring the propagation speed if the apparent wavelength is corrected. This point will be discussed in the forthcoming paper.

## V. CONCLUSION

2D numerical simulation results of the subcritical millimeter-wave discharge are presented to confirm the IDSI method. The IDSI method can identify continuous comb-shaped and discrete discharge structures by examining the time-variation of the standing-wave intensity in the numerical simulation because a noisy waveform is obtained for the discrete structure, whereas a smooth waveform is obtained for the continuous comb-shaped discharge structure. For both the continuous and discrete discharge structures, the Fourier spectrum of the waveform has a strong peak, and its frequency is proportional to the propagation speed of the ionization-front, which agrees with the theoretical formula proposed in the previous experiment. Thus, the location method can also measure the ionization-front propagation speed by conducting a Fourier transform of the time-variation waveform of the standing-wave intensity even if the discharge structure is changed. Moreover, the high-frequency components observed in the discrete structure are induced because the millimeter-wave absorption during a new plasma spot formation causes a small time-variation of the waveform. These findings can be applied to measure the timing of the plasma spot formation in discharge experiments without using a high-speed camera. We plan to conduct a discharge experiment to measure standing-wave intensity, which will be presented in a forthcoming paper.

## ACKNOWLEDGMENTS

This research was supported by JSPS KAKENHI (Grant No. 23KJ0103). Numerical simulations were conducted using JSS3 TOKI-SORA supercomputer system of JAXA. We would like to acknowledge Editage (www.editage.jp) for English language editing.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**S. Suzuki:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). **M. Takahashi:** Resources (equal); Supervision (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Plasma Physics and Plasma Electronics*, edited by L. Kovrizhnykh (Nova Science Publishers, 1989), p. 241.

*Proceedings of 43rd International Conference on Infrared, Millimeter, and Terahertz Waves*(IEEE, 2018), pp. 1–2.

*The Finite Difference Time Domain Method for Electromagnetics*