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.

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 λ / 4,6–10 comb-shaped,14–18 tree-branching,8,20 diffusive,19,22 and complex filamentary structures.19,22 In the λ / 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 λ / 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.

FIG. 1.

Schematics of the structure of ionization-front observed in past experiments.16 (a) Front view, namely, from the viewpoint of the beam source. (b) Side view.

FIG. 1.

Schematics of the structure of ionization-front observed in past experiments.16 (a) Front view, namely, from the viewpoint of the beam source. (b) Side view.

Close modal

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 λ / 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 model23–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.

FIG. 2.

Hypotheses of ionizing-front propagation of the comb-shaped structure. These are cross-sectional views across the central axis. (a) Previous studies have considered that granular plasma spots propagate continuously even in the central region. (b) To explain the discharge structure (Fig. 1) and propagation speed depending on the place, we consider that the plasma spots propagate discretely in the central region, but they propagate continuously in the peripheral region.

FIG. 2.

Hypotheses of ionizing-front propagation of the comb-shaped structure. These are cross-sectional views across the central axis. (a) Previous studies have considered that granular plasma spots propagate continuously even in the central region. (b) To explain the discharge structure (Fig. 1) and propagation speed depending on the place, we consider that the plasma spots propagate discretely in the central region, but they propagate continuously in the peripheral region.

Close modal
Microwave reflectometry is a suitable method that does not rely on an optical camera to observe dense plasma in high-resolution time.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 λ / 2,
(1)
where λ is the wavelength of the millimeter wave, T p is a period in which interference is enhanced when phases of the incident and reflected millimeter waves are aligned, f p = 1 / T p is the peak frequency that is obtained by conducting a Fourier transform of interference wave. However, it is unclear whether Eq. (1) can be used when forming discrete structures. Moreover, it has not been experimentally confirmed whether the location method can be applied for discharge structure identification.

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 λ / 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.

FIG. 3.

Schematics of the proposed method to identify discharge structures.

FIG. 3.

Schematics of the proposed method to identify discharge structures.

Close modal

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.

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 model24,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 studies23–25 if the compressible fluid module is off.

The electromagnetic-wave of the transverse magnetic (TM) mode propagation is described using Maxwell’s equations
(2)
(3)
where E is the electric-field vector, μ 0 is the vacuum magnetic permeability, H is the magnetic field vector, ε 0 is the vacuum electric permittivity, and J is the current density vector which is estimated by Ohm’s law as J = e n e v e. Here, e is the elementary charge, n e is the electron number density, and v e is the mean velocity vector of the electrons, which is estimated using the following equation of motion of electrons:
(4)
where m e is the electron mass and ν m is the elastic collision frequency between the electrons and neutral particles. Here, ν m is calculated using ν m = 5.3 × 10 5 p, where p is the pressure of neutral gas. Equations (2) and (3) were discretized using the finite-difference time-domain (FDTD) method with total field/scattered field formulation.48 Equation (4) is numerically integrated using an explicit Euler method.
The time scale of electromagnetic-wave propagation is much shorter than that of fluid dynamics; thus, the electric field E is treated as a root mean square (RMS) value during one millimeter-wave period in the fluid models to reduce the computational time.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:
(5)
where n i is the number densities ( i = e , p , n), Γ i is the density flux vector, and S i is the source term. Here, the subscripts e , p , n are electrons, positive ions, and negative ions, respectively. The positive and negative ion density fluxes were assumed to be zero, as in a previous study.23 The electron density flux Γ e is obtained by the following equation:
(6)
(7)
where D eff is the effective diffusion coefficient,23, ζ = ν i τ m = λ D 2 / L 2 is the transition parameter. Here, ν i is the ionization frequency, λ D is the Debye length, and L is the scale length of the electron distribution. τ m = ε 0 / [ e n e ( μ i + μ e ) ] is the Maxwell relaxation time, μ i and μ e are the ion and electron mobility, respectively. D e is the free diffusion coefficient of the electron, which is calculated by D e = μ e k 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 μ e = e / m e ν m, and the ion mobility is assumed to be μ e / 200 in this study.23 The source terms S = S e , S p, and S n are obtained using the following equations:
(8)
(9)
(10)
where S ion is the electron-impact ionization rate, and S att is the electron-attachment rate. Other chemical reactions were ignored in this study because they were not important at higher beam intensities, even under subcritical conditions.40 Reaction rates were calculated using the following approximation:25 
(11)
(12)
where E eff is the effective electric field as follows:
(13)
where E rms is the RMS value of local electric-field intensity and ω is the angular frequency of the millimeter wave. E rms was calculated as follows:
(14)
where T is the period of the millimeter wave and E z is the z component of the electric-field vector E. The x and y components of E were zero because the irradiated electromagnetic-wave mode was the TM mode in this simulation. A calculation of E rms is conducted along with the FDTD solver. Equation (5) was numerically integrated using the explicit Euler method. The diffusion term in Eq. (5) is discretized using the central difference scheme.
Neutral fluid dynamics were introduced using the following 2D Euler equation to evaluate the effect of gas expansion by gas heating via plasma:
(15)
Here,
(16)
(17)
(18)
(19)
where ρ is the mass density of the neutral gas, ε is the total energy per unit volume, u is the x component of the flow velocity, v is the y component of the flow velocity, and W g is the heat source term which is calculated by the following equation:
(20)
where J E is the Joule heating. It is assumed that 30% of the RMS Joule heating energy is transferred to the translational energy.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
The frequency of the millimeter-wave beam f was set to 170 GHz following past experiments.14,16 The wavelength λ 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 λ × 5 λ, and the grid size Δ x and Δ y are Δ x = Δ y = λ / 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 simulations23,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 / λ , y / λ ) = ( 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 × 10 15 m 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 3, to avoid the initial electron effects on the discharge structure. The time step of the plasma and the neutral fluid solver Δ T fluid was calculated as below:
(21)
The time step of the electromagnetic-wave solver Δ T EM must satisfy the Courant–Friedrichs–Lewy (CFL) condition as below:
(22)
where c is the speed of light. In this study, Δ T EM = Δ T fluid 160 was used to satisfy Eq. (22). While the calculation of the plasma and neutral fluid solvers is conducted one time, the calculation of electromagnetic-wave solver is conducted 160 times, and E rms and J E are passed to the plasma and the neutral fluid solvers. The boundary conditions of the electromagnetic waves were set to Mur’s first order absorption condition. Those of reaction–diffusion and Euler equations were set as the outflow boundary conditions. The computational code was parallelized using a message passing interface (MPI) library. The computations were conducted using ten nodes (480 cores) of the JSS3 TOKI-SORA supercomputer system of Japan Aerospace Exploration Agency (JAXA).
FIG. 4.

2D calculation domain in this study.

FIG. 4.

2D calculation domain in this study.

Close modal

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.

FIG. 5.

Distributions of the electron number density n e. (a) Comb-shaped structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

FIG. 5.

Distributions of the electron number density n e. (a) Comb-shaped structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

Close modal

Figures 6(a) and 6(b) are x t diagram of the electron number density n e along y / λ = 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 × 10 9 for the comb-shaped structure and 1 × 10 13 cm 3 for the discrete structure. Figure 6(a) shows that the ionization-front propagates continuously toward the beam-source direction ( x direction). However, the ionization-front of the discrete structure propagates as a leapfrog, as shown in Fig. 6(b).

FIG. 6.

The x t diagram of the electron number density along y / λ = 2.5 axis. (a) Comb-shaped structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

FIG. 6.

The x t diagram of the electron number density along y / λ = 2.5 axis. (a) Comb-shaped structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

Close modal

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 / λ = 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.

FIG. 7.

The x t diagram of the RMS value of the local electric-field intensity E rms along y / λ = 2.5 axis. (a) Comb-shape structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

FIG. 7.

The x t diagram of the RMS value of the local electric-field intensity E rms along y / λ = 2.5 axis. (a) Comb-shape structure at E 0 , rms = 1.4 MV / m. (b) Discrete structure at E 0 , rms = 2.2 MV / m.

Close modal

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 / λ , y / λ ) = ( 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.

FIG. 8.

(a) The time-variation waveform of the RMS value of local electric-field intensity E rms at ( x / λ , y / λ ) = ( 0.01 , 2.5 ). (b) The Fourier spectrum of the waveform in the comb-shape structure at E 0 , rms = 1.4 MV / m.

FIG. 8.

(a) The time-variation waveform of the RMS value of local electric-field intensity E rms at ( x / λ , y / λ ) = ( 0.01 , 2.5 ). (b) The Fourier spectrum of the waveform in the comb-shape structure at E 0 , rms = 1.4 MV / m.

Close modal

Figure 9(a) shows a time-variation of the standing-wave intensity E rms / E 0 , rms at the position of ( x / λ , y / λ ) = ( 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.

FIG. 9.

(a) The time-variation waveform of the RMS value of local electric-field intensity E rms at ( x / λ , y / λ ) = ( 0.01 , 2.5 ). (b) The Fourier spectrum of the waveform in the discrete structure at E 0 , rms = 2.2 MV / m.

FIG. 9.

(a) The time-variation waveform of the RMS value of local electric-field intensity E rms at ( x / λ , y / λ ) = ( 0.01 , 2.5 ). (b) The Fourier spectrum of the waveform in the discrete structure at E 0 , rms = 2.2 MV / m.

Close modal

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.

In the numerical simulation in this study, the continuous propagation occurred at E 0 , rms = 1.4 , 1.6 , 1.8 MV / m. In particular, Fig. 8(a) shows the time-varying waveform at ( x / λ , y / λ ) = ( 0.01 , 2.5 ) at E 0 , rms = 1.4 MV / m. The movement of the ionization-front induces a periodic waveform because the standing waves also move with the ionization-front, which behaves as a fixed end for wave reflection. The local maximum and the local minimum in the time-varying waveform correspond to the anti-node and node of the standing waves, respectively. Thus, the ionization-front propagation speed u ion in the numerical simulation is expected to have a relationship of Eq. (1) between the wavelength of the millimeter-wave λ and the period of the waveform of the standing wave T p because the distance between the node and the anti-node is λ / 2. Table I shows the ionization-front propagation speed obtained from the x t diagram ( u ion , m), the peak frequency measured from the Fourier spectrum ( f p , m), and the theoretical peak frequency obtained by assuming Eq. (1) ( f p , th). Therefore, f p , th was evaluated using the below equation:
(23)
Here, u ion , m for each E 0 , rms is evaluated from the slope of the black line in x- t diagram of n e, which is calculated using the least square method. In the range of 1.4 E 0 , rms 1.8 MV / m, Table I shows that f p , m have a good agreement with f p , th. Therefore, it was verified that Eq. (1) worked in our numerical simulation, similar to the experiment for a continuous structure. It is expected that the location method can be applied to discrete structures because a standing wave is also formed, which moves along with the ionization-front. As the next step, we will verify the location method using numerical simulations when a discrete structure is induced.
TABLE I.

Measured propagation speed uion,m, the measured peak frequency of the waveform fp,m, and theoretical peak frequency of the waveform fp,th calculated by Eq. (1) in each E0,rms.

E0,rms (MV/m)uion,m (m/s)fp,m (MHz)fp,th (MHz)
1.4 374 0.43 0.424 
1.6 416 0.47 0.472 
1.8 426 0.53 0.483 
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 
E0,rms (MV/m)uion,m (m/s)fp,m (MHz)fp,th (MHz)
1.4 374 0.43 0.424 
1.6 416 0.47 0.472 
1.8 426 0.53 0.483 
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.

FIG. 10.

Peak frequency measured from the Fourier spectrum, f p , m, against the propagation speed of the ionization-front measured from the x t diagram, u ion , m. The black line is the peak frequency f p , th theoretically calculated by Eq. (1).

FIG. 10.

Peak frequency measured from the Fourier spectrum, f p , m, against the propagation speed of the ionization-front measured from the x t diagram, u ion , m. The black line is the peak frequency f p , th theoretically calculated by Eq. (1).

Close modal

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.

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 μ 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.

FIG. 11.

Waveform around t = 2.2 μ s in the case of the discrete structure is formed at E 0 , rms = 2.2 MV / m.

FIG. 11.

Waveform around t = 2.2 μ s in the case of the discrete structure is formed at E 0 , rms = 2.2 MV / m.

Close modal

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.

FIG. 12.

Distribution of electron number density n e and normalized local electric-field intensity E rms / E 0 , rms. (a) n e at t = 2.00 μ s. (b) E rms / E 0 , rms at t = 2.00 μ s. (c) n e at t = 2.14 μ s. (d) E rms / E 0 , rms at t = 2.14 μ s. (e) n e at t = 2.16 μ s. (f) E rms / E 0 , rms at t = 2.16 μ s. (g) n e at t = 2.21 μ s. (h) E rms / E 0 , rms at t = 2.21 μ s.

FIG. 12.

Distribution of electron number density n e and normalized local electric-field intensity E rms / E 0 , rms. (a) n e at t = 2.00 μ s. (b) E rms / E 0 , rms at t = 2.00 μ s. (c) n e at t = 2.14 μ s. (d) E rms / E 0 , rms at t = 2.14 μ s. (e) n e at t = 2.16 μ s. (f) E rms / E 0 , rms at t = 2.16 μ s. (g) n e at t = 2.21 μ s. (h) E rms / E 0 , rms at t = 2.21 μ s.

Close modal

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 / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 1.75 , 2.06 ) and ( 1.75 , 2.93 ) in front of the central plasma spot at ( x / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 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 / λ = 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 μ s in Fig. 11. At this time, Fig. 12(c) shows that two diagonal plasma spots are formed at ( x / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 0.01 , 2.5 ) because the front of the cut-off region moved λ / 4 toward beam source ( x / λ = 0) along the x axis.

At t = 2.16 μ 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 / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 0.01 , 2.5 ) is decreased at 2.16 μ s (A2) compared with 2.14 μ s (A1) because of the absorption by the central plasma spot during the formation.

FIG. 13.

1D distributions of normalized electric-field intensity E rms / E 0 , rms on the central axis at 2.14, 2.16, and 2.21  μ s.

FIG. 13.

1D distributions of normalized electric-field intensity E rms / E 0 , rms on the central axis at 2.14, 2.16, and 2.21  μ s.

Close modal

In contrast, the standing-wave intensity takes the local maximum value again at t = 2.21 μ s as point A3 in Fig. 11. A new central plasma spot has been completely formed at ( x / λ , y / λ ) = ( 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  μ s, the standing intensity is decreased to x / λ = 0 from x / λ = 1.6, which indicates that the reflected millimeter wave is diverging. However, the standing-wave intensity is constant between x / λ = 0 to 1 at 2.21 μ 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 μ 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 μ 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 / λ , y / λ ) = ( 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 λ / 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.

FIG. 14.

Distribution of electron number density n e and normalized local electric-field intensity E rms / E 0 , rms. (a) n e at t = 2.32 μ s. (b) E rms / E 0 , rms at t = 2.32 μ s. (c) n e at t = 2.35 μ s. (d) E rms / E 0 , rms at t = 2.35 μ s. (e) n e at t = 2.40 μ s. (f) E rms / E 0 , rms at t = 2.40 μ s.

FIG. 14.

Distribution of electron number density n e and normalized local electric-field intensity E rms / E 0 , rms. (a) n e at t = 2.32 μ s. (b) E rms / E 0 , rms at t = 2.32 μ s. (c) n e at t = 2.35 μ s. (d) E rms / E 0 , rms at t = 2.35 μ s. (e) n e at t = 2.40 μ s. (f) E rms / E 0 , rms at t = 2.40 μ s.

Close modal

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.

In Secs. IV AIV C, the observation point was set on the axis at ( x / λ , y / λ ) = ( 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 / λ , y / λ ) = ( 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 / λ = 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 / λ , y / λ ) = ( 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 μ s, whereas the waveform observed at (0.01,4.5) has a local maximum near t = 11.58 μ 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.

FIG. 15.

Observation points selected in this section and the normalized electric-field intensity distribution E rms / E 0 , rms at t = 11.58 μ s for the beam intensity of E 0 , rms = 1.4 MV / m.

FIG. 15.

Observation points selected in this section and the normalized electric-field intensity distribution E rms / E 0 , rms at t = 11.58 μ s for the beam intensity of E 0 , rms = 1.4 MV / m.

Close modal
FIG. 16.

Time-varying waveform of the standing-wave intensity E rms / E 0 , rms at each observation points. (a) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 2.5 ), ( 0.01 , 4.5 ), and ( 0.01 , 4.99 ) for E 0 , rms = 1.4 MV / m. (b) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 4.99 ), ( 0.5 , 4.99 ), and ( 1 , 4.99 ) for E 0 , rms = 1.4 MV / m. (c) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 2.5 ), ( 0.01 , 4.5 ), and ( 0.01 , 4.99 ) for E 0 , rms = 2.2 MV / m. (d) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 4.99 ), ( 0.5 , 4.99 ), and ( 1 , 4.99 ) for E 0 , rms = 2.2 MV / m.

FIG. 16.

Time-varying waveform of the standing-wave intensity E rms / E 0 , rms at each observation points. (a) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 2.5 ), ( 0.01 , 4.5 ), and ( 0.01 , 4.99 ) for E 0 , rms = 1.4 MV / m. (b) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 4.99 ), ( 0.5 , 4.99 ), and ( 1 , 4.99 ) for E 0 , rms = 1.4 MV / m. (c) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 2.5 ), ( 0.01 , 4.5 ), and ( 0.01 , 4.99 ) for E 0 , rms = 2.2 MV / m. (d) E rms / E 0 , rms observed at ( x / λ , y / λ ) = ( 0.01 , 4.99 ), ( 0.5 , 4.99 ), and ( 1 , 4.99 ) for E 0 , rms = 2.2 MV / m.

Close modal

Figures 17(a)17(d) show the Fourier spectrum of the waveform observed at ( x / λ , y / λ ) = ( 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 λ / 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 λ / 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.

FIG. 17.

Fourier spectrum of the time-varying waveform E rms / E 0 , rms observed at each observation point when the comb-shaped structure is formed at E 0 , rms = 1.4 MV / m. (a) ( x / λ , y / λ ) = ( 0.01 , 4.5 ). (b) ( x / λ , y / λ ) = ( 0.01 , 4.99 ) . (c) ( x / λ , y / λ ) = ( 0.5 , 4.99 ). (d) ( x / λ , y / λ ) = ( 1 , 4.99 ).

FIG. 17.

Fourier spectrum of the time-varying waveform E rms / E 0 , rms observed at each observation point when the comb-shaped structure is formed at E 0 , rms = 1.4 MV / m. (a) ( x / λ , y / λ ) = ( 0.01 , 4.5 ). (b) ( x / λ , y / λ ) = ( 0.01 , 4.99 ) . (c) ( x / λ , y / λ ) = ( 0.5 , 4.99 ). (d) ( x / λ , y / λ ) = ( 1 , 4.99 ).

Close modal

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
Yu. Ya.
Brodskii
,
S. V.
Golubev
,
V. G.
Zorin
,
A. G.
Luchinin
, and
V. E.
Sernenov
, “
New mechanism of gasdynamics propagation of a discharge
,”
Sov. Phys. JETP.
57
,
989
993
(
1983
).
2.
N. A.
Bogatov
,
S. V.
Golubev
, and
V. G.
Zorin
, “
Ionizing radiation from a microwave discharge
,”
Sov. Tech. Phys. Lett.
9
,
382
383
(
1983
).
3.
Yu. Ya.
Brodskii
,
I. P.
Venediktov
,
S. V.
Golubev
,
V. G.
Zorin
, and
I. A.
Kossy
, “
Nonequilibrium microwave discharge in air at atmospheric pressure
,”
Sov. Tech. Phys. Lett.
10
,
77
79
(
1984
).
4.
N. A.
Bogatov
,
Yu. V.
Bykov
,
N. P.
Venediktov
,
S. V.
Golubev
,
V. G.
Zorin
,
A. G.
Eremeev
, and
V. E.
Semenov
, “
Gasdynamic propagation of a nonequilibrium microwave discharge
,”
Sov. J. Plasma Phys.
12
,
416
420
(
1986
).
5.
G. M.
Batanov
,
S. I.
Gritsinin
,
I. A.
Kossy
,
A. N.
Magunov
,
V. P.
Silakov
, and
N. M.
Tarasova
, “High-pressure microwave discharges,” in Plasma Physics and Plasma Electronics, edited by L. Kovrizhnykh (Nova Science Publishers, 1989), p. 241.
6.
Y.
Hidaka
,
E. M.
Choi
,
I.
Mastovsky
,
M. A.
Shapiro
,
J. R.
Sirigiri
, and
R. J.
Temkin
, “
Observation of large arrays of plasma filaments in air breakdown by 1.5 MW 110 GHz gyrotron pulses
,”
Phys. Rev. Lett.
100
,
035003
(
2008
).
7.
Y.
Hidaka
,
E. M.
Choi
,
I.
Mastovsky
,
M. A.
Shapiro
,
J. R.
Sirigiri
,
R. J.
Temkin
,
G. F.
Edmiston
,
A. A.
Neuber
, and
Y.
Oda
, “
Plasma structures observed in gas breakdown using a 1.5 MW, 110 GHz pulsed gyrotron
,”
Phys. Plasmas
16
,
055702
(
2009
).
8.
A.
Cook
,
M.
Shapiro
, and
R.
Temkin
, “
Pressure dependence of plasma structure in microwave gas breakdown at 110 GHz
,”
Appl. Phys. Lett.
97
,
011504
(
2010
).
9.
A. M.
Cook
,
J. S.
Hummelt
,
M. A.
Shapiro
, and
R. J.
Temkin
, “
Observation of plasma array dynamics in 110 GHz millimeter-wave air breakdown
,”
Phys. Plasmas
18
,
100704
(
2011
).
10.
S. C.
Schaub
,
J. S.
Hummelt
,
W. C.
Guss
,
M. A.
Shapiro
, and
R. J.
Temkin
, “
Electron density and gas density measurements in a millimeter-wave discharge
,”
Phys. Plasmas
23
,
083512
(
2016
).
11.
K. V.
Khodataev
, “Propagation of microwave subcritical streamer discharge against radiation by brunching and looping,” AIAA Paper 2008-1405, 2008.
12.
K. V.
Artem’ev
,
G. M.
Batanov
,
N. K.
Berezhetskaya
,
V. D.
Borzosekov
,
A. M.
Davydov
,
L. V.
Kolik
,
E. M.
Konchekov
,
I. A.
Kossyi
,
A. E.
Petrov
,
K. A.
Sarksyan
,
V. D.
Stepakhin
, and
N. K.
Kharchev
, “
Discharge in a subthreshold microwave beam as an unusual type of ionization wave
,”
Plasma Phys. Rep.
44
,
1146
1153
(
2018
).
13.
K. V.
Artem’ev
,
G. M.
Batanov
,
N. K.
Berezhetskaya
,
V. D.
Borzosekov
,
A. M.
Davydov
,
L. V.
Kolik
,
E. M.
Konchekov
,
I. A.
Kossyi
,
A. E.
Petrov
,
K. A.
Sarksyan
,
V. D.
Stepakhin
, and
N. K.
Kharchev
, “
Location of the front of a subthreshold microwave discharge and some specificities of its propagation
,”
Plasma Phys. Rep.
45
,
965
972
(
2019
).
14.
Y.
Oda
,
K.
Komurasaki
,
K.
Takahashi
,
A.
Kasugai
, and
K.
Sakamoto
, “
Plasma generation using high-power millimeter-wave beam and its application for thrust generation
,”
J. Appl. Phys.
100
,
113307
(
2006
).
15.
Y.
Oda
,
K.
Kajiwara
,
K.
Takahashi
,
A.
Kasugai
,
K.
Sakamoto
, and
K.
Komurasaki
, “
In-tube shock wave driven by atmospheric millimeter-wave plasma
,”
Jpn. J. Appl. Phys.
48
,
116001
(
2009
).
16.
Y.
Oda
,
M.
Takahashi
,
N.
Ohnishi
,
K.
Komurasaki
,
K.
Sakamoto
, and
T.
Imai
, “
A study on the macroscopic self-organized structure of high-power millimeter-wave breakdown plasma
,”
Plasma Sources Sci. Technol.
29
,
075010
(
2020
).
17.
Y.
Oda
,
M.
Takahashi
,
K.
Tabata
,
N.
Ohnishi
,
K.
Komurasaki
, and
K.
Sakamoto
, “Frequency dependence of atmospheric millimeter wave breakdown plasma,” in Proceedings of 43rd International Conference on Infrared, Millimeter, and Terahertz Waves (IEEE, 2018), pp. 1–2.
18.
M.
Fukunari
,
S.
Tanaka
,
R.
Shinbayashi
,
Y.
Yamaguchi
,
Y.
Tatematsu
, and
T.
Saito
, “
Observation of a comb-shaped filamentary plasma array under subcritical condition in 303 GHz millimetre-wave air discharge
,”
Sci. Rep.
9
,
17972
(
2019
).
19.
K.
Tabata
,
Y.
Harada
,
Y.
Nakamura
,
K.
Komurasaki
,
H.
Koizumi
,
T.
Kariya
, and
R.
Minami
, “
Experimental investigation of ionization front propagating in a 28 GHz gyrotron beam: Observation of plasma structure and spectroscopic measurement of gas temperature
,”
J. Appl. Phys.
127
,
063301
(
2020
).
20.
K.
Shimamura
,
J.
Yamasaki
,
K.
Miyawaki
,
R.
Minami
,
T.
Kariya
,
J.
Yang
, and
S.
Yokota
, “
Propagation of microwave breakdown in argon induced by a 28 GHz gyrotron beam
,”
Phys. Plasmas
28
,
033505
(
2021
).
21.
K.
Tabata
,
Y.
Oda
,
K.
Komurasaki
,
A.
Manabe
, and
R.
Kawashima
, “
Non-equilibrium aerodynamics between ionization-wave and shock-wave fronts in millimetre-wave supported detonation
,”
Jpn. J. Appl. Phys.
62
,
116001
(
2023
).
22.
K. V.
Artem’ev
,
G. M.
Batanov
,
N. K.
Berezhetskaya
,
V. D.
Borzosekov
,
A. M.
Davydov
,
L. V.
Kolik
,
E. M.
Konchekov
,
I. A.
Kossyi
,
D. V.
Malakhov
,
I. V.
Moryakov
,
A. E.
Petrov
,
K. A.
Sarksyan
,
V. D.
Stepakhin
, and
N. K.
Kharchev
, “
Changes in structure of subthreshold discharge in air occurring with decreasing microwave radiation intensity
,”
Plasma Phys. Rep.
48
,
170
177
(
2022
).
23.
J. P.
Boeuf
,
B.
Chaudhury
, and
G. Q.
Zhu
, “
Theory and modeling of self-organization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure
,”
Phys. Rev. Lett.
104
,
015002
(
2010
).
24.
B.
Chaudhury
,
J. P.
Boeuf
, and
G. Q.
Zhu
, “
Pattern formation and propagation during microwave breakdown
,”
Phys. Plasmas
17
,
123505
(
2010
).
25.
G. Q.
Zhu
,
J. P.
Boeuf
, and
B.
Chaudhury
, “
Ionization-diffusion plasma front propagation in a microwave field
,”
Plasma Sources Sci. Technol.
20
,
035007
(
2011
).
26.
K.
Kourtzanidis
,
J. P.
Boeuf
, and
F.
Rogier
, “
Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air
,”
Phys. Plasmas
21
,
123513
(
2014
).
27.
K.
Kourtzanidis
,
F.
Rogier
, and
J. P.
Boeuf
, “
ADI-FDTD modeling of microwave plasma discharges in air towards fully three-dimensional simulations
,”
Comput. Phys. Commun.
195
,
123513
(
2015
).
28.
K.
Kourtzanidis
and
L.
Raja
, “
Limitations of the effective field approximation for fluid modeling of high frequency discharges in atmospheric pressure air: Application in resonant structures
,”
Phys. Plasmas
24
,
112105
(
2017
).
29.
P.
Shu
and
P.
Zhao
, “
Effect of ambient gas species on microwave breakdown pattern
,”
Jpn. J. Appl. Phys.
60
,
126001
(
2021
).
30.
P.
Ghosh
and
B.
Chaudhury
, “
Efficient dynamic mesh refinement technique for simulation of HPM breakdown-induced plasma pattern formation
,”
IEEE Trans. Plasma Sci.
51
,
66
76
(
2023
).
31.
M.
Takahashi
and
N.
Ohnishi
, “
Plasma filamentation and shock wave enhancement in microwave rockets by combining low-frequency microwaves with external magnetic field
,”
J. Appl. Phys.
120
,
063303
(
2016
).
32.
M.
Takahashi
, “
Development of plasma fluid model for a microwave rocket supported by a magnetic field
,”
J. Phys.: Conf. Ser.
905
,
012024
(
2017
).
33.
M.
Takahashi
,
Y.
Kageyama
, and
N.
Ohnishi
, “
Joule-heating-supported plasma filamentation and branching during subcritical microwave irradiation
,”
AIP Adv.
7
,
055206
(
2017
).
34.
M.
Takahashi
and
K.
Komurasaki
, “
Discharge from a high-intensity millimeter wave beam and its application to propulsion
,”
Adv. Phys.: X
3
,
1417744
(
2018
).
35.
M.
Takahashi
and
N.
Ohnishi
, “
Gas-species-dependence of microwave plasma propagation under external magnetic field
,”
J. Appl. Phys.
124
,
173301
(
2018
).
36.
Y.
Nakamura
,
K.
Komurasaki
,
M.
Fukunari
, and
H.
Koizumi
, “
Numerical analysis of plasma structure observed in atmospheric millimeter-wave discharge at under-critical intensity
,”
J. Appl. Phys.
124
,
033303
(
2018
).
37.
M.
Takahashi
and
N.
Ohnishi
, “
Gas propellant dependency of plasma structure and thrust performance of microwave rocket
,”
J. Appl. Phys.
125
,
163303
(
2019
).
38.
M.
Takahashi
and
N.
Ohnishi
, “Numerical study of discharge and thrust generation in a microwave rocket,” AIAA Paper 2019-1242, 2019.
39.
Y.
Nakamura
and
K.
Komurasaki
, “
Theory and modeling of under-critical millimeter-wave discharge in atmospheric air induced by high-energy excited neutral-particles carried via photons
,”
Plasma Sources Sci. Technol.
29
,
105017
(
2020
).
40.
S.
Suzuki
,
K.
Hamasaki
,
M.
Takahashi
,
C.
Kato
, and
N.
Ohnishi
, “
Numerical analysis of structural change process in millimeter-wave discharge at subcritical intensity
,”
Phys. Plasmas
29
,
093507
(
2022
).
41.
S.
Suzuki
,
C.
Kato
,
M.
Takahashi
, and
N.
Ohnishi
, “
Plasma propagation via radiation transfer in millimeter-wave discharge under subcritical condition
,”
J. Phys.: Conf. Ser.
2207
,
012046
(
2022
).
42.
A.
Manabe
,
H.
Britz
,
F.
Gutierrez
,
A.
van de Wetering
,
T.
Kinoshita
,
K.
Komurasaki
,
R.
Kawashima
,
H.
Sekine
, and
H.
Koizumi
, “
Numerical analysis on multi-cycle operation of microwave rocket with reed valve air-breathing system (in Japanese)
,”
J. Jpn. Soc. Aeronaut. Space Sci.
72
,
29
40
(
2024
).
43.
E.
Mazzucato
, “
Microwave reflectometry for magnetically confined plasmas
,”
Rev. Sci. Instrum.
69
,
2201
2217
(
1998
).
44.
R.
Grosso Ferreira
,
J.
Vicente
,
F.
Silva
,
B.
Gonçalves
,
L.
Cupido
, and
M.
Lino da Silva
, “
Reflectometry diagnostics for atmospheric entry applications: State-of-the-art and new developments
,”
CEAS Space J.
16
,
1
18
(
2024
).
45.
G. M.
Batanov
,
V. D.
Borzosekov
,
L. V.
Kolik
,
E. M.
Konchekov
,
D. V.
Malakhov
,
A. E.
Petrov
,
K. A.
Sarksyan
,
V. D.
Stepakhin
, and
N. K.
Kharchev
, “
Self-action of a gaussian beam of microwaves in the subthreshold field generated by the waves in air
,”
Plasma Phys. Rep.
47
,
598
602
(
2021
).
46.
A. M.
Cook
,
J. S.
Hummelt
,
M. A.
Shapiro
, and
R. J.
Temkin
, “
Millimeter wave scattering and diffraction in 110 GHz air breakdown plasma
,”
Phys. Plasmas
20
,
043507
(
2013
).
47.
G. J. M.
Hagelaar
and
L. C.
Pitchford
, “
Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models
,”
Plasma Sources Sci. Technol.
14
,
722
733
(
2005
).
48.
K. S.
Kunz
and
R. J.
Luebbers
,
The Finite Difference Time Domain Method for Electromagnetics
(
CRC Press
,
1983
).
49.
Y.
Wada
and
M. S.
Liou
, “A flux splitting scheme with high-resolution and robustness for discontinuities,” AIAA Paper 1994-83, 1994.
50.
B.
van Leer
, “
Toward the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method
,”
J. Comput. Phys.
32
,
101
136
(
1979
).