Low-frequency molecular fluctuations in the translational nonequilibrium zone of one-dimensional strong shock waves are characterized for the first time in a kinetic collisional framework in the Mach number range $2\u2264M\u226410$. Our analysis draws upon the well-known bimodal nature of the probability density function (PDF) of gas particles in the shock, as opposed to their Maxwellian distribution in the freestream, the latter exhibiting two orders of magnitude higher dominant frequencies than the former. Inside the (finite-thickness) shock region, the strong correlation between perturbations in the bimodal PDF and fluctuations in the normal stress suggests introducing a novel two-bin model to describe the reduced-order dynamics of a large number of collision interactions of gas particles. Our model correctly predicts two orders of magnitude differences in fluctuation frequencies in the shock vs those in the freestream and is consistent with the small-amplitude fluctuations obtained from the highly resolved direct simulation Monte Carlo computations of the same configuration. The variation of low-frequency fluctuations with changes in the conditions upstream of the shock revealed that these fluctuations can be described by a Strouhal number, based on the bulk velocity upstream of the shock and the shock-thickness based on the maximum density gradient inside the shock, that remains practically independent of Mach number in the range examined.

## I. INTRODUCTION

It is well known that numerical solutions of the deterministic Euler and Navier–Stokes equations (NSE) require special treatment if one is to obtain reliable flow features containing shocks. Solutions of differential equations describing fluid flow motion, in either a shock-capturing^{1,2} or a shock-fitting context,^{3,4} can predict both location and unsteadiness of shocks accurately, but the description of the internal shock structure and predictions of the gas properties inside a shock must be modeled in the context of such solutions. Information about the interior of shocks may be obtained from the numerical solution of systems of equations that are more general than the Navier–Stokes equations. The first documentation of differences in shock structure predicted by the Navier–Stokes equations and solution of the more general Bhatnagar–Gross–Krook model for the Boltzmann equation of transport is due to Liepmann *et al.*,^{5} who showed that the Navier–Stokes equations are inadequate in describing the internal structure of a shock wave at high Mach numbers ($M\u22653$) and that the differences between the two solutions increase in the low-pressure region of the shock layer with an increase in Mach number. Soon afterward, Bird^{6} solved the exact Boltzmann equation using the stochastic DSMC method^{7} and quantified the level of strong translational nonequilibrium in the interior of shocks.

The significance of the correct description of not only the location and motion of shocks but also of their internal dynamical structure cannot be overstated in fluid mechanics. Among other classes of flows, shocks are responsible for the production of sound in supersonic jets as shown by the extensive studies of their interaction with turbulence and isolated vortices.^{8–20} The shock-wave/boundary-layer interactions (SBLIs) on compression ramps, cones, and flat plates are widely investigated for their role in generating separation bubbles, unsteadiness, and large surface heat and pressure fluxes.^{21–23} It is also recognized that the study of receptivity of shocks to freestream or induced disturbances is of importance in the investigation of transition in hypersonic boundary layers.^{24–28}

However, the vast majority of research in laminar to turbulent transition and SBLIs focuses on the physics of transition in the boundary layer and not the shock.^{29–33} In fact, research on the effect of kinetic fluctuations on triggering transition from laminar to turbulent flows has been explored using Landau–Lifshitz's theory of fluctuating hydrodynamics^{34} in several works on the receptivity of boundary layers in the incompressible^{35,36} and high-speed compressible^{37,38} regimes. In the phenomenological approach of fluctuating hydrodynamics, stochastic fluxes are added to the stress tensor and heat flux vector in the NSE formulation to account for the effect of molecular fluctuations on the flow field. Using this approach, Luchini^{36} showed that the presence of kinetic fluctuations is a non-negligible source of disturbances in a boundary layer and sets an upper bound on the transition Reynolds number. Fedorov and Tumin^{37} showed that in a compressible flat-plate boundary layer, the kinetic fluctuations could trigger random wave packets of Tollmien–Schlichting waves or a Mack second-mode instability. However, this approach cannot be easily extended to account for the kinetic fluctuations in shocks, where the fluctuation amplitude is known to be larger than the predictions of the equilibrium theory.^{39}

Therefore, to understand the role of molecular fluctuations in the shock layer, kinetic methods that provide insight into the molecular nature of the shock in terms of particle velocities and energy distribution function are important. However, understanding detailed molecular collisions in two- and three-dimensional (2D and 3D) flows is challenging because of the tremendous length and time scales related to boundary layers, instabilities, and unsteadiness of SBLIs. Our study aims at closing this theoretical gap by examining the origin of molecular fluctuations in the well-known one-dimensional (1D) shock of a monatomic gas, argon, in a kinetic framework. We show the surprising result that even simple one-dimensional shocks exhibit *low-frequency fluctuations* that have been, up until now, ignored, yet their presence may well contribute to laminar to turbulent transition in supersonic and hypersonic flows.

Historically, solutions of the internal structure of strong shocks were sought using kinetic models^{5} as they allow for anisotropy of stresses and heat fluxes, which are not accounted for by the traditional Navier–Stokes–Fourier constitutive relations. Bird^{6} successfully modeled 1D argon shock structures using the DSMC method and obtained shock-thicknesses that were in good agreement with the experiments of Alsmeyer.^{40} Since then, the method has been widely used for modeling the steady-state behavior of strong shocks^{41–45} as well as unsteady flows without^{46–50} and with shocks.^{51} Stefanov *et al.*^{39} simulated Mach 26 bow shock over a 2D cylinder using DSMC and attributed an increase in velocity fluctuations inside the shock layer to “thermal nonequilibrium effects”; however, they did not delve deeper into the reason for these effects, the changes in timescales of fluctuations in comparison with an equilibrium state, and their dependence on the strength of the shock wave. Our work attempts to answer these questions by investigating the fluctuations in molecular velocity and energy distribution functions obtained from the solution of the Boltzmann equation. This work will show that the major role of bimodality in the distribution function inside the shock is to change the dominant frequencies of molecular fluctuations compared to an equilibrium freestream.

Our past work has exploited the fidelity of DSMC to understand the critical role of shocks in hypersonic SBLIs. Tumuklu *et al.*^{52,53} simulated laminar SBLIs in a Mach 16 axisymmetric flow over a double-cone and observed strong coupling between the complex shock structure and the separation bubble. For the freestream unit Reynolds number of $Re1=3.74\xd7105$ m^{−1}, they found oscillations of the detached (bow) and separation shocks characterized by a Strouhal number of 0.078. In this work, we will show that the Strouhal number associated with the low-frequency fluctuations in an isolated 1D shock falls within a similar range of $St=$ 0.001–0.02. Sawant *et al.*^{54} extended the capability of DSMC in modeling the computationally challenging 3D laminar SBLIs of a Mach 7 flow over a double-wedge at near-continuum input conditions through the use of adaptively refined octree grids. In addition, we have investigated the stability of the spanwise-periodic laminar separation bubble in the double-wedge flow to self-excited, small-amplitude, spanwise homogeneous perturbations.^{55,56} We demonstrate the existence of the strong coupling of the shock structure with the separation bubble and show that linear amplification of the most unstable 3D flow perturbations leads to synchronization of low-frequency unsteadiness of the separation bubble with the unsteadiness of the manifold of shock structures corresponding to a Strouhal number of 0.028.

The remainder of the paper is organized as follows: Sec. II describes the DSMC numerical implementation and the nonequilibrium aspects of a 1D shock structure of argon. It then describes the fluctuations in the overall stress inside the shock with that in the freestream and contrasts the differences in their frequencies. From additional analyses of cross-correlations, it is hypothesized that the differences in fluctuations are caused by the long-time collision interaction of particles in two modes of the bimodal energy distribution of particles. To prove this hypothesis, we construct a simplified two-energy-bin ordinary differential equation (ODE) model in Sec. III similar to the predator–prey model of Lotka–Volterra^{57–59} but with modified terms accounting for intermolecular collisions between the two bins. The evaluation of rate coefficients used in the model is also described in this section, and the simplification of the collision processes is discussed in detail in Appendix B. In Sec. IV, we discuss the dynamics of the two-energy-bin model and compare with the broadband frequencies obtained from the power spectral density (PSD) analyses of fluctuations in the DSMC residuals. Section V establishes a range of Strouhal numbers for Mach numbers 2–10 as well as for variations in upstream temperature at a given Mach number. Using this range, an estimate of low frequencies is made and qualitatively compared with experiments. Finally, Sec. VI summarizes the present findings.

## II. A PARTICLE REPRESENTATION OF NONEQUILIBRIUM FLUCTUATIONS INSIDE A ONE-DIMENSIONAL SHOCK

### A. DSMC simulation methodology and properties of a 1D shock

The DSMC solution of the shock is obtained using a 1D version of the Scalable Unstructured Gas-dynamic Adaptive mesh-Refinement (SUGAR) DSMC solver,^{54} which makes use of a binary adaptive mesh refinement structure, where each computational Cartesian “root” cell is recursively refined into smaller cells until its size is smaller than the local mean-free-path. The smallest cells, referred to as “leaf” or “collision” cells, are used to select neighboring collision partners to perform binary elastic collisions between particles of a monatomic gas using the majorant frequency scheme.^{60} The macroscopic flow and transport parameters are computed based on the statistical equations of kinetic theory^{7} and shown on larger root cells. In terms of traditional DSMC solvers, each reference to a computational cell in this paper means a Cartesian root cell. The gas is assumed to follow a variable hard sphere (VHS) molecular model^{7} with viscosity index $\omega =0.81$, mass $m=6.637\xd710\u221226$ kg, and reference diameter $dr=4.17\xd710\u221210$ m at reference temperature *T _{r}* = 273 K. We compared our DSMC profiles of normalized viscous stress and heat flux along the direction of shock propagation with the results of Ohwada

^{61}obtained by using a finite-difference Boltzmann solver for a 1D, Mach 3 flow of argon and obtained excellent agreement.

The DSMC simulation was initialized using the Rankine–Hugoniot jump conditions at $X/\lambda 1=0$, where *X* is the streamwise direction normal to the shock and *λ*_{1} is the upstream mean-free-path. At the upstream boundary, $X/\lambda 1<$ 0, computational particles are introduced using a directed local Maxwellian flow at a number density, bulk velocity, and temperature of $n1=1\xd71022$ m^{−3}, $ux,1=3572.24$ m s^{−1}, and $Ttr,1=710$ K, respectively, to generate hypersonic flow. At the subsonic downstream boundary, $X/\lambda 1>$ 0, a specularly reflecting surface moving at a bulk velocity of $ux,2$, equal to the downstream Rankine–Hugoniot velocity is used (see Ref. 7, Chap.12). The domain length is 0.02 m, that is, $X/\lambda 1=\u221257.45$ to 57.45, which is sufficiently large so that the boundaries are placed sufficiently far from the shock layer but small enough so that the shock does not move from its initialized central location during the early transient period. In this domain, the flow takes only ∼6 *μ*s to transition from the jump condition to form a typical shock structure, a time much shorter than the low-frequency molecular fluctuations that we will discuss. Bird's STABIL boundary condition^{7} is used to prevent random walk effects from causing the shock to move; however, because of the high number density and the use of a large number of computational particles ($\u2248$475 000), the drift was found to be very small. Nevertheless, the STABIL algorithm has not been previously used in a 1D shock simulation where we wish to characterize molecular fluctuations that generate low-frequency unsteadiness. For this reason, additional clarification is provided in Appendix A to demonstrate that the frequencies that we will observe from DSMC in Secs. II B and V are not influenced by the STABIL algorithm.

Using the aforementioned computational approach, Fig. 1 shows time-averaged profiles of macroscopic flow and transport parameters obtained from the DSMC simulation of a 1D, Mach 7.2 shock in argon. Unless otherwise indicated, subsequent figures and results are based on these 1D shock conditions. The time-averaged or mean quantities discussed in this work were calculated after the 6-*μ*s transient period for additional time steps until ∼0.2 ms. Figure 1(a) shows deviation of the $X\u2212$ and average lateral directional temperatures, $Ttr,x$ and based on axial symmetry, $Ttr,n=(Ttr,y+Ttr,z)/2$, respectively, from the overall translational temperature, *T _{tr}*, within a region of $\u221216<X/\lambda 1<6$, indicating the presence of translational nonequilibrium in the shock. In the equilibrium regions ($X/\lambda 1<\u221216$ and $X/\lambda 1>6$), it can be seen that all directional temperatures are equal to the overall temperature. Within the region $2<X/\lambda 1<6$, the gradient of overall translational temperature asymptotes to zero, yet there is a strong deviation of directional temperatures from each other, similar to that in the DSMC simulation of a Mach 8 argon shock by Bird.

^{6}Figure 1 also shows the locations of the numerical probes

*F*($X/\lambda 1=\u221229$),

*P*($X/\lambda 1=\u22121$), and

*D*($X/\lambda 1=29$), based on which later discussions will follow.

In addition, Fig. 1(b) shows the deviation of time-averaged pressure, *p*, from the *X*-directional overall stress, *σ _{xx}*, inside the shock as well as the normal viscous stress components,

*τ*and $\tau nn=(\tau yy+\tau zz)/2$, calculated based on the relationship

_{xx}where *δ _{ij}* is the Kronecker delta function. The pressure and stress components are normalized by the upstream parameter $\rho 1\beta 1\u22122$, where $\beta =m/2\kappa bTtr$ is the inverse of the most probable speed of molecules,

*m*is mass,

*κ*is the Boltzmann constant, $\rho =nm$ is the mass density, and

_{b}*n*is the number density. The existence of these non-zero stresses leads to a finite-thickness of the shock wave, which, when modeled by accounting for molecular thermal fluctuations using a stochastic method, such as DSMC, reveals nonequilibrium bimodal velocity and energy distributions of particles. This paper will show that these bimodal distributions exhibit a dominant frequency that is two orders of magnitude lower than those found in the freestream.

With respect to the equilibrium regions of the shock, Figs. 1(c) and 1(d) show excellent agreement between the DSMC-computed velocity and number density fluctuations and theory, as was observed by Stefanov *et al.*^{39} For a dilute gas, theory predicts that at equilibrium, statistical fluctuations in macroscopic quantities of bulk velocity, *u _{x}*, and number of particles,

*N*, have a standard deviation of $\kappa b\u27e8Ttr\u27e9/m\u27e8N\u27e9$ and $\u27e8N\u27e9$, respectively, where the brackets denote time or ensemble average of macroscopic quantities.

^{34,62}The DSMC standard deviations are computed for each computational cell of width $\Delta x=1\xd710\u22124$ m for instantaneous data collected at every time step of $\Delta t=3$ ns from 6

*μ*s to 0.2 ms. Note that since the velocity fluctuations shown in Fig. 1(c) are expressed in terms of DSMC computational particles, the amplitude of actual thermal fluctuations is obtained by division by a factor of $Fn$, where $Fn=107$ is the number of dilute gas molecules represented by each DSMC computational particle. Therefore, the 0.344% standard deviation of velocity fluctuations in the freestream in Fig. 1(c) corresponds to an actual value of $1.09\xd710\u22124$ %, a small yet significantly large number on the scale of small amplitude perturbations considered in shock-dominated flows.

### B. Examination of nonequilibrium fluctuations

We are particularly interested in the nonequilibrium zone of the shock layer, where a significant deviation from equilibrium is seen in DSMC-computed standard deviations, as shown in Figs. 1(c) and 1(d). Note that the standard deviations peak at the location of absolute maximum gradients of respective flow parameters, that is, at location $X/\lambda 1=0$. In contrast to the findings of Stefanov *et al.*,^{39} the density fluctuations inside a shock also deviate from the equilibrium Poisson law, although their magnitude is much smaller than the velocity or energy fluctuations.

To understand the frequencies of fluctuations in the shock layer, we focus on probe *P* located in the vicinity of maximum absolute bulk velocity gradient in the shock to obtain the time history of fluctuations of the normalized *X*-directional overall stress, $\sigma xx\beta 12\rho 1\u22121$, about the time-averaged mean. These data are important because fluctuations in *σ _{xx}* correspond to fluctuations in $Ttr,x$ given that $\sigma xx=\rho RTtr,x$ and density fluctuations are negligible. The original instantaneous signal is window-averaged with a moving time-window of 0.3

*μ*s (100 timesteps) such that frequencies greater than 3.33 MHz are removed. The power spectral density (PSD) of this mean-subtracted, window-averaged signal reveals frequency content in fluctuations of

*σ*, as shown in Fig. 2. For spectral estimation here and elsewhere in the paper, Welch's method

_{xx}^{63,64}is used in SciPy

^{65}software with two Hann-window weighted segments of data sampled with a frequency of 333 MHz prior to the Fast Fourier Transform (FFT) such that the frequency resolution is 0.9 kHz. In Appendix A, the sensitivity of the PSD obtained using DSMC particle data to the value of

*F*is discussed and it is shown that when

_{n}*F*is set to unity (i.e

_{n}*.,*no possible numerical amplification of noise), the main features of the non-dimensionalized low-frequency spectrum do not change.

At probe *P*, the PSD shows a broadband of low frequencies that goes up to 90 kHz frequency, as shown in Fig. 2(a). This upper bound of the broadband is defined as the frequency at which the normalized cumulative energy (NCE), obtained from normalizing and cumulatively summing the PSD spectrum, exhibits an inflection point. This low-frequency broadband is defined by its weighted average of 37.5 kHz and a standard deviation of 21.4 kHz. Figure 2(b) shows the contours of PSD plotted on the axes $X/\lambda 1$ vs frequency created by interpolating the PSD data at each $X/\lambda 1$ location spatially separated by $X/\lambda 1=1$ in the Tecplot-360^{66} software using the inverse-distance algorithm with default parameters (exponent = 3.5, point selection=Octant, Number of points = 8). The figure reveals that such low-frequency broadband is observed within $\u22126<X/\lambda 1<3$, the region of strong translational nonequilibrium seen earlier in Fig. 1.

In contrast to the nonequilibrium region of the shock, Fig. 2(c) shows that the PSD at probe *F* in the freestream contains widely distributed energy across the whole spectral limit of 166 MHz and does not exhibit a noticeable inflection point within this limit. It also shows that most of the spectral energy is concentrated near the frequency of 3.5 MHz, which is the inverse of the mean-collision-time of 0.284 *μ*s. This is consistent with the analysis of Kogan (see Ref. 67, pp. 82 and 124) for a gas in absolute equilibrium (zero bulk velocity, constant temperature, and number of particles), which showed that small perturbations to the velocity distribution function and its moments decay exponentially with a characteristic timescale of mean collision time. In addition, from the comparison of PSDs in Figs. 2(a) and 2(c), it can be noted that the low-frequency fluctuations in the shock have about two orders of magnitude larger amplitude than the freestream, and therefore, low-frequency fluctuations of shock may be even more important in laminar to turbulent transition of high-speed boundary layers than, as commonly believed, freestream kinetic fluctuations. The order of magnitude differences in these timescales of fluctuations will be explained by the new dynamics model discussed in Sec. IV.

The key to understanding the fluctuations in the overall stress, *σ _{xx}*, is to correlate the fluctuations in the PDF $f\xi x$ of the normalized

*X*-directional energy of particles,

*ξ*, and its mean,

_{x}*μ*. Note that $\xi x=(vx2\u2212ux2)\beta 2$, where

*v*is the

_{x}*X*-directional instantaneous molecular velocity of particles, which is the sum of their bulk and thermal components,

*u*and

_{x}*c*, respectively. The mean,

_{x}*μ*, is related to overall stress,

*σ*, as

_{xx}Figure 3 shows that the behavior of the PDF changes at different locations in the flow, that is, from the upstream to downstream regions relative to the shock, $f\xi x$ changes from a nearly symmetric equilibrium distribution [Fig. 3(a)] to a one-sided, asymmetric equilibrium distribution [Fig. 3(b)]. At the location of maximum density gradient in the shock [Fig. 3(c)], it can be seen that the density function is bimodal with an inflection point at $\xi x=1.33$, where the portions $\u22120.86<\xi x<1.33$ and $1.33<\xi x<3.52$ have shapes similar to PDFs in the subsonic and hypersonic regions, shown in Figs. 3(b) and 3(a), respectively. Note that the bimodal PDF can be expressed as a linear combination of upstream and downstream contributions of equilibrium PDFs,^{68} similar to the Mott–Smith model of the bimodal velocity distribution.^{69}

Toward that end, the cross-correlation coefficient of *μ* and the number of particles as a function of *ξ _{x}* are defined as

where

Note that *ξ _{x}* is discretized into 200 energy bins from its minimum to maximum value. $N(\xi x,w)$ is the total number of DSMC particles within the normalized

*X*-directional energy space

*ξ*and $\xi x+\Delta \xi x$ from time window

_{x}*w*to

*w*+

*1, and $\u27e8N(\xi x)\u27e9w$ denotes the number of particles in the same energy space but averaged over time-windows,*

*W*. Similarly, $\mu (w)$ is the instantaneous mean of the PDF $f\xi x$, from time-window

*w*to

*w*+

*1, whereas $\u27e8\mu \u27e9w$ is the mean computed by averaging over all time-windows. $\Sigma N(\xi x)$ and $\Sigma \mu $ are the standard deviations in the fluctuations of $N(\xi x,w)$ and $\mu (w)$ about their respective means. Figures 4(a) and 4(b) show the fluctuations in number of particles $[N(\xi x,w)\u2212\u27e8N(\xi x)\u27e9w]$ in the energy space*

*ξ*at probes

_{x}*P*and

*F*, respectively.

^{70}

The analysis of the correlation functions combined with particle distributions suggests an approach for grouping of particles into a finite number of energy bins whose dynamics can then be further analyzed. Starting with Fig. 5, the distribution of $c(\xi x)$ at probes *P*, inside a shock, and *F*, in the freestream, is compared. $c(\xi x)$ is calculated using a total number of time windows, *W*, of 646, each representing a time interval of 0.3 *μ*s starting from *t *=* *6 *μ*s. At probe *P*, Fig. 5(a) shows that $c(\xi x)$ has a strong negative correlation between $\u22120.86<\xi x<1.33$ and positive correlation for $1.33<\xi x<3.52$. Opposite signs of $c(\xi x)$ means that particles in these parts lose and gain particles, respectively, when *μ* and *σ _{xx}* increase and vice versa, and understanding the underlying mechanisms in these regions, such as intermolecular collisions and flux exchange with adjacent cells, is necessary to understand fluctuations in

*σ*. These two zones in the energy space are denoted as energy bins “

_{xx}*A*” and “

*B*,” which are coarser than the size of the histogram bins, $\Delta \xi x=0.0115$, used to obtain smooth variations of $f\xi x$. We ignore the energy space beyond $\xi x>3.52$, because $c(\xi x)$ remains negative but small in comparison with other two bins and the total average number of particles in this region is only 3.5% of total particles at probe

*P*, which can be calculated from the cumulative sum of the distribution of average number of particles shown in Fig. 6(a). The demarcation between the coarse energy bins

*A*and

*B*is at $c(\xi x)=0$, which is also the location of inflection point. Examination of $c(\xi x)$ for probe

*F*in the freestream [see Fig. 5(b)] shows that, in contrast, the maximum magnitude of $c(\xi x)$ is not more than 0.2, indicating that the fluctuations are random in nature and even if a strong correlation exists, it must be on timescales smaller than the window size of 0.3

*μ*s used to obtain $c(\xi x)$. Comparison of the correlation functions at probes

*P*and

*F*also shows that the distribution of $c(\xi x)$ at the latter location is almost symmetric about the inflection point at $\xi x=0$, which is consistent with the

*X*-directional energy distribution having a nearly symmetric shape, as shown in Fig. 3(a).

If we now group particles into two *X*-directional energy bins, *A* and *B*, which are denoted in Figs. 5(a) and 5(b) for probes *P* and *F*, respectively, we can study the ratios of the instantaneous to time averaged particles, $NA/\u27e8NA\u27e9t$ and $NB/\u27e8NB\u27e9t$ as a function of time at these two probes, as shown in Figs. 5(c) and 5(d). The average number of DSMC computational particles per computational cell in energy bins *A* and *B* is $\u27e8NA\u27e9t=1088$ and $\u27e8NB\u27e9t=611$ at probe *P*, and $\u27e8NA\u27e9t=510$ and $\u27e8NB\u27e9t=493$ at probe *F*, which can be converted to number density $(m\u22123)$ by multiplication of an $Fn=1\xd7107$ and division by the volume of computational cell, $\Delta v=1\xd710\u221212\u2009m3$. From the comparison of these two figures, it can be seen that fluctuations in these ratios are larger in magnitude and exhibit longer timescales (0.01–0.05 ms) at probe P inside the shock in comparison to probe F in the freestream, where fluctuations appear to be random as their timescales are on the order of mean collision time of 0.284 *μ*s, smaller than the time interval used to time-average the signal (0.3 *μ*s). The standard deviations in the fluctuations of the ratios $NA/\u27e8NA\u27e9t$ and $NB/\u27e8NB\u27e9t$ are 4.5 and 2.1%, respectively, at probe *P*, vs 1.8 and 1.6% at probe *F*. In addition, a noticeable negative correlation between the fluctuations (out-of-phase fluctuations) of two energy bins at probe *P* can also be seen from Fig. 5(c), while such dynamics are absent at probe *F*.

Finally, to demonstrate the generality of the difference in properties between probes *P* and *F*, Fig. 6 shows a comparison of the time-averaged number of particles as a function of *ξ _{x}* at probes

*P*and

*F*, and at their respective left-adjacent computational cells. For probe

*P*[Fig. 6(a)], two inflection points can be identified at $\xi x=1.33$ and 3.52, where $c(\xi x)$ changes sign. We can also see that the average number of particles in energy bin

*A*is larger at probe

*P*than in the cell to its left, whereas in bin

*B*it is lower. This is consistent with the fact that from the upstream to downstream region inside a shock, the contribution of the upstream symmetric distribution decreases, while the downstream asymmetric distribution increases. At probe

*F*[Fig. 6(b)], however, the PDF of energies does not change between adjacent cells and there is no inflection point.

## III. COLLISION PROCESSES THAT DEFINE THE TWO-BIN MODEL

Although the distribution of *X*-directional particle energies is continuous, we propose a two-bin model to understand the role of particle collisions between the two bins, *A* and *B*, and how they induce low-frequency time dynamics in fluctuations of macroscopic flow parameters such as normalized overall stress, *σ _{xx}*, inside a shock. We first develop a simple two energy-bin model similar to the Lotka–Volterra's predator–prey model.

^{57–59}We then evaluate the collisions rate coefficients for different energy transfer processes using the DSMC particle distribution data and show that instead of sustained particle oscillations, the solution to our ODE contains dampening terms that cause the periodic fluctuations to die out on timescales an order of magnitude longer than the period of oscillation.

### A. The two-bin energy model

The dynamics of the number of particles in energy bins *A* and *B* may be written as

where the first and second terms on the right are the change in number of particles due to collisions and fluxes, respectively. In addition, since the average number of particles in energy bins *A* and *B* constitutes 96.5% of the total particles, we can assume them to be mutually exclusive and write

The flux coefficient *F _{j}* for bin “

*j”*is defined as the difference between the net average influx from the left boundary and outflux from the right boundary of the “

*j”*-type particles per second divided by the average number of “

*j”*-type particles,

where $\Delta t$ is the DSMC time step, superscripts “*l”* and “*r”* refer to the left and right boundaries of the computational cell, respectively, and “in” and “out” refer to the incoming and outgoing particles, respectively, as shown in Fig. 7 for location *P* and energy bins *A* and *B*.

Flux coefficients are required in the evaluation of Eq. (3). Using the values of *A-* and *B*-type particles given in Table I, flux values of *F _{A}* and

*F*= −1 277 925.8 and 2 500 681.6 s

_{B}^{–1}, respectively, at probe

*P*and zero in the freestream (as expected) are obtained. To make sure that the predictions of the model are not affected by instantaneous fluctuations or their numerical amplification, the DSMC-derived fluxes provided in Table I are temporally averaged from 6

*μs*to 0.21 ms, a time period ∼8 times longer than the timescale corresponding to the average low-frequency observed in the PSD. Note that at probe

*P*, $FA<0$ because on average more

*A*-type particles travel to a cell downstream than came in from the upstream, as also seen from Fig. 6(a). The opposite is true for type-

*B*particles, resulting in $FB>0$. This is synonymous with the fact that if there were no collisions, that is, $[dNAdt]coll=0$, the number of

*A*particles, those that mainly represent the subsonic part of the bimodal distribution, would decrease, while the number of

*B*particles, those that mainly represent the upstream hypersonic flow, would increase.

Probes . | $NAl,in$ . | $NAl,out$ . | $NAr,out$ . | $NAr,in$ . | $NBl,in$ . | $NBl,out$ . | $NBr,out$ . | $NBr,in$ . |
---|---|---|---|---|---|---|---|---|

$P$ | 41.46 | 7.19 | 47.40 | 8.96 | 69.84 | 0.68 | 65.37 | 0.80 |

$F$ | 49.28 | 0 | 49.28 | 0.0 | 58.63 | 0.0 | 58.63 | 0.0 |

Probes . | $NAl,in$ . | $NAl,out$ . | $NAr,out$ . | $NAr,in$ . | $NBl,in$ . | $NBl,out$ . | $NBr,out$ . | $NBr,in$ . |
---|---|---|---|---|---|---|---|---|

$P$ | 41.46 | 7.19 | 47.40 | 8.96 | 69.84 | 0.68 | 65.37 | 0.80 |

$F$ | 49.28 | 0 | 49.28 | 0.0 | 58.63 | 0.0 | 58.63 | 0.0 |

To evaluate the collision terms in Eq. (3), we need to identify all possible collision processes that can cause a loss or gain of type-*A* and type-*B* particles. Note that different collision processes can occur at different rates, which could change the dynamics of particle exchange between the two bins and consequently the timescale of fluctuations in *σ _{xx}*. Such collision processes are listed in column two of Table II, where the

*type*of

*A*particle is labeled based on the collision process that it undergoes, denoted by a subscript “

*i*,” the relevant pre-and post-collisional states are specified where the latter is denoted by a primed superscript, and the rate coefficient for each

*fundamental*collision process,

*k*, unless it is a “compound” rate, which is denoted by a tilde. The table shows that there are six fundamental processes (second column) and eight types of

_{i}*A*particles, $Aa,Ab\u2032,Ac,Ac\u2032,Ad,Ad\u2032,Ae\u2032,$ and

*A*, that can cause a net change in the number of particles in energy bins

_{f}*A*and

*B*. For example, $Ab\u2032$ particles are the postcollisional particles located in energy bin

*A*that are formed due to the collision process

*b*, while

*A*are the precollisional particles in energy bin

_{c}*A*that take part in collision process

*c*. Using the DSMC-derived PDFs for all sub-types of type-

*A*particles of the normalized

*X*-directional energy, $\xi x=(vx2\u2212ux2)\beta 2$, and total transverse energy, $\xi y+\xi z=(vy2+vz2)\beta 2$, denoted by $f\xi xA$ and $f\xi y+\xi zA$, respectively (Fig. 8), we can reduce the number of types of

*A*particles further, which allows us to write processes

*p*(column two) in the form of processes

_{i}*P*(column three). In what follows, we only describe the grouping process for processes

_{i}*a*and

*b*, while details for the other processes are given in Appendix B. Then, the

*P*processes are further simplified to processes

_{i}*Q*(column 4), as described in Appendix B.

_{i}Id . | Detailed collision process . | Collision process after grouping alike A-particles . | Simplified collision process . |
---|---|---|---|

$i$ . | $pi$ . | $Pi$ . | $Qi$ . |

a | A_{a} + B $ka\u2192$ B + B | A_{a} + B $ka\u2192$ B + B | A + B $k\u0303a\u2192$ B + B |

b | B + B $kb\u2192$ $Ab\u2032$ + B | B + B $kb\u2192$ A_{a} + B | B + B $kb\u2192$ A + B |

c | A_{c} + B $kc\u2192$ $Ac\u2032$ + $Ac\u2032$ | A_{c} + B $kc\u2192$ $Ac\u2032$ + $Ac\u2032$ | A + B $k\u0303c\u2194$ $Ac\u2032$ + $Ac\u2032$ |

d | A_{d} + A_{d} $kd\u2192$ $Ad\u2032$ + B | $Ac\u2032$ + $Ac\u2032$ $kd\u2192$ A_{c} + B | ⋯ |

e | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ |

f | A_{f} + A_{f} $kf\u2192$ B + B | $Ae\u2032$ + $Ae\u2032$ $kf\u2192$ B + B | ⋯ |

g^{a} | A_{g} + B $kg\u2192$ B + $Ag\u2032$ | A_{g} + B $kg\u2192$ B + A_{g} | ⋯ |

h^{a} | A_{h} + A_{h} $kh\u2192$ $Ah\u2032$ + $Ah\u2032$ | A_{h} + A_{h} $kh\u2192$ A_{h} + A_{h} | ⋯ |

i^{a} | B + B $ki\u2192$ B + B | B + B $ki\u2192$ B + B | ⋯ |

Id . | Detailed collision process . | Collision process after grouping alike A-particles . | Simplified collision process . |
---|---|---|---|

$i$ . | $pi$ . | $Pi$ . | $Qi$ . |

a | A_{a} + B $ka\u2192$ B + B | A_{a} + B $ka\u2192$ B + B | A + B $k\u0303a\u2192$ B + B |

b | B + B $kb\u2192$ $Ab\u2032$ + B | B + B $kb\u2192$ A_{a} + B | B + B $kb\u2192$ A + B |

c | A_{c} + B $kc\u2192$ $Ac\u2032$ + $Ac\u2032$ | A_{c} + B $kc\u2192$ $Ac\u2032$ + $Ac\u2032$ | A + B $k\u0303c\u2194$ $Ac\u2032$ + $Ac\u2032$ |

d | A_{d} + A_{d} $kd\u2192$ $Ad\u2032$ + B | $Ac\u2032$ + $Ac\u2032$ $kd\u2192$ A_{c} + B | ⋯ |

e | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ | B + B $ke\u2192$ $Ae\u2032$ + $Ae\u2032$ |

f | A_{f} + A_{f} $kf\u2192$ B + B | $Ae\u2032$ + $Ae\u2032$ $kf\u2192$ B + B | ⋯ |

g^{a} | A_{g} + B $kg\u2192$ B + $Ag\u2032$ | A_{g} + B $kg\u2192$ B + A_{g} | ⋯ |

h^{a} | A_{h} + A_{h} $kh\u2192$ $Ah\u2032$ + $Ah\u2032$ | A_{h} + A_{h} $kh\u2192$ A_{h} + A_{h} | ⋯ |

i^{a} | B + B $ki\u2192$ B + B | B + B $ki\u2192$ B + B | ⋯ |

^{a}

These processes cause no net change in the number of particles of energy bins $A$ and *B*.

Collision process *p _{a}* describes a mechanism for type-

*A*particles in energy bin

_{a}*A*that collide with particles in energy bin

*B*and now belong to bin

*B*resulting in the loss of type-

*A*particles. Simultaneously, process

*p*describes the collision process of two bin

_{b}*B*particles that causes one of them to move to energy bin

*A*leading to a gain of type-

*A*particles, denoted by $Ab\u2032$. However, Figs. 8(a) and 8(b) show that at location

*P*the energy distributions $f\xi xA$ and $f\xi y+\xi zA$ for

*A*and

_{a}*A*particles are the same. Therefore, we can group

_{b}*A*- and

_{b}*A*-type particles and write the second reaction

_{a}*P*as shown in column three of Table II. The same holds true for probe

_{b}*F*as shown in Figs. 8(c) and 8(d); however, the distributions are different from those at probe

*P*.

We can also group all collisions between type-*A* particles that do not lead to a net change in the energy distribution of particles within bin *A* as process *p _{h}* and see from Fig. 8 that the energy distributions of all such type-

*A*particles before and after collisions,

*A*and $Ah\u2032$, is the same at probe

_{h}*P*as well as

*F*. Note that at probe

*F*, the distribution $f\xi y+\xi zA$ for

*A*-type particles matches with the

_{a}*A*-type, that is, in the freestream the transverse energy of particles taking part in processes

_{h}*P*and

_{a}*P*is the same as those in

_{b}*P*. The differences at probe

_{h}*P*indicate the role of transverse translational modes in the collision processes

*P*and

_{a}*P*. Similar to process

_{b}*p*, process

_{h}*p*, is a bin-exchange process, which causes no net change in the number of particles in bins

_{g}*A*and

*B*. Particles simply swap energy bins by exchanging

*X*-directional energies. Using similar logic, additional simplifications can be made to the detailed collision processes,

*p*, to construct the

_{i}*P*column of Table II based on the energy distribution functions at locations

_{i}*P*and

*F*, as described in Appendix B.

To complete the dynamical model for particles moving between energy bins *A* and *B*, an expression is required for the net change in the number of particles in energy bin *A* due to collisions from all possible processes related to *A*-type particles, that is,

where $NAa,\u2009NAc,\u2009NAc\u2032$, and $NAe\u2032$ are the number of particles of *A _{a}*,

*A*, $Ac\u2032$, and $Ae\u2032$ types, respectively, which are defined based on

_{c}*P*processes. Using standard rate equation formulations for elementary reactions (see Ref. 71, Sec. 16.15), the terms on the right hand side can be expressed as

_{i}The conditions in the 1D shock allow us to make additional simplifications to Eq. (8), as described in Appendix B, so that using Eqs. (3), (4), and (B9) we obtain the final system of ordinary differential equations used to study perturbation dynamics,

### B. Evaluation of rate coefficients

Finally, to solve the dynamical system of equations we need to evaluate the rate coefficients in Eq. (9). A summary of the expressions for the rate coefficients used in Eq. (9) and their values are presented in Table III for probes *P* and *F* based on the DSMC collision data collected during the simulation, summarized in Table IV. Note that $\u27e8Ci\u27e9t$ is the average number of DSMC collisions per time step per cell volume that take part in the *i ^{th}* detailed collision process of particle types

*R*

_{1}and

*R*

_{2}[see. Eq. (B1)]. Similar to the flux coefficients, $\u27e8Ci\u27e9t$ were temporally averaged from 6

*μs*to 0.21 ms. The rate coefficients

*k*/($s\u22121$) shown in Table III can be converted to the traditional units of $K/(m3\u2009s\u22121)$ by multiplying them with the volume of computational cell, $\Delta v=1\xd710\u221212\u2009m3$, and dividing by the parameter, $Fn=1\xd7107$. By this conversion, it can be easily shown that the rates $K\u0303a$,

*K*range from $1.8\u22125.0\xd710\u221217$ m

_{b}^{3}s

^{−1}at probes

*P*and

*F*.

Rate coefficient . | Formula . | Probe P
. | Probe F
. |
---|---|---|---|

$k\u0303a$ | $\u27e8Ca\u27e9t/\Delta t\u27e8NA\u27e9t\u27e8NB\u27e9t$ | 188.9 | 497.0 |

k _{b} | $\u27e8Cb\u27e9t/\Delta t\u27e8NB\u27e9t2$ | 417.5 | 513.5 |

$k\u0303c$ | $(\u27e8Cc\u27e9t\u2212\u27e8Cd\u27e9t)/\Delta t\u27e8NA\u27e9t\u27e8NB\u27e9t$ | 1843.5 | 0 |

k _{e} | $\u27e8Ce\u27e9t/\Delta t\u27e8NB\u27e9t2$ | 75.0 | 0 |

Rate coefficient . | Formula . | Probe P
. | Probe F
. |
---|---|---|---|

$k\u0303a$ | $\u27e8Ca\u27e9t/\Delta t\u27e8NA\u27e9t\u27e8NB\u27e9t$ | 188.9 | 497.0 |

k _{b} | $\u27e8Cb\u27e9t/\Delta t\u27e8NB\u27e9t2$ | 417.5 | 513.5 |

$k\u0303c$ | $(\u27e8Cc\u27e9t\u2212\u27e8Cd\u27e9t)/\Delta t\u27e8NA\u27e9t\u27e8NB\u27e9t$ | 1843.5 | 0 |

k _{e} | $\u27e8Ce\u27e9t/\Delta t\u27e8NB\u27e9t2$ | 75.0 | 0 |

Probes . | $\u27e8Ca\u27e9t$ . | $\u27e8Cb\u27e9t$ . | $\u27e8Cc\u27e9t$ . | $\u27e8Cd\u27e9t$ . | $\u27e8Ce\u27e9t$ . | $\u27e8Cf\u27e9t$ . | $\u27e8Cg\u27e9t$ . | $\u27e8Ch\u27e9t$ . | $\u27e8Ci\u27e9t$ . |
---|---|---|---|---|---|---|---|---|---|

$P$ | 0.377 | 0.468 | 5.074 | 1.396 | 0.084 | 0.0079 | 2.980 | 8.849 | 1.481 |

$F$ | 0.375 | 0.374 | 0.382 | 0.383 | 0.0 | 0.0 | 1.028 | 0.893 | 0.887 |

Probes . | $\u27e8Ca\u27e9t$ . | $\u27e8Cb\u27e9t$ . | $\u27e8Cc\u27e9t$ . | $\u27e8Cd\u27e9t$ . | $\u27e8Ce\u27e9t$ . | $\u27e8Cf\u27e9t$ . | $\u27e8Cg\u27e9t$ . | $\u27e8Ch\u27e9t$ . | $\u27e8Ci\u27e9t$ . |
---|---|---|---|---|---|---|---|---|---|

$P$ | 0.377 | 0.468 | 5.074 | 1.396 | 0.084 | 0.0079 | 2.980 | 8.849 | 1.481 |

$F$ | 0.375 | 0.374 | 0.382 | 0.383 | 0.0 | 0.0 | 1.028 | 0.893 | 0.887 |

Note that in the freestream, $\u27e8Ca\u27e9t\u2248\u27e8Cb\u27e9t$, due to detailed balance in a local equilibrium condition, which dictates that each detailed collision must be balanced by its inverse collision process (see Ref. 72, pp. 38, 39, and 42). Furthermore, the net rate, $k\u0303c$, is nearly zero in the freestream because again by principle of detailed balance, $\u27e8Cc\u27e9t\u2248\u27e8Cd\u27e9t$, as shown in Table IV. On the other hand, in the shock, a factor of 3.63 difference between the average values of these two types of collisions results in $K\u0303c=1.83\xd710\u221216$ m^{3} s^{−1}. Also, from Table IV it can be seen that both $\u27e8Ce\u27e9t$ and $\u27e8Cf\u27e9t$ are zero in the freestream, which means that these processes take place only in conditions of translational nonequilibrium. As a result, the rate *K _{e}* is zero in the freestream, whereas in the shock it is two orders of magnitude slower, $7.5\xd710\u221218$ m

^{3}s

^{−1}, than $K\u0303c$.

The rate coefficients given in Table III can be compared with theoretical estimate of a rate coefficient, *K _{th}*, for a forward collision process of a quasi-equilibrium gas (see Ref. 72, pp. 213–216), which is defined as

The rate of change of *R*_{1}-type particles per unit volume due to collisions is written as

where $nR1$ and $nR2$ are number densities of *R*_{1} and *R*_{2} type particles, respectively. For a gas following the VHS molecular model, the rate coefficient *K _{th}* can be expressed as

where

where $\u27e8\sigma \u27e9t$ is the average equilibrium collision cross section for the VHS gas, $\u27e8gr\u27e9t$ is the time-averaged relative velocity, *γ* is the fraction of collisions for which the relative translational energy along the line of centers of colliding molecules exceeds the activation energy of *E _{a}* [see Ref. 7, Sec. 6.2 and Eq. (4.72))], and Γ is the gamma function. Also,

*m*is the reduced mass [see Ref. 7, Eq. (2.7)],

_{r}*d*is the VHS molecular diameter (see Ref. 7, Eq. (4.63)], and

*ε*is a symmetry factor equal to two because R

_{1}= R

_{2}. Table V shows that the estimated parameters in Eq. (10) and average rate coefficients at probes

*P*and

*F*located inside the shock and the freestream, respectively, are consistent with the values presented in Table III.

Parameters . | At probe P
. | At probe F
. |
---|---|---|

$\u27e8Ttr\u27e9t$/(K) | 10 440 | 710 |

$\u27e8gr\u27e9t$/(m s^{−1}) | 3326.2 | 867.6 |

$\u27e8\sigma \u27e9t$/(m^{2}) | $1.81\xd710\u221219$ | 4.16 × 10^{–19} |

E/(J) _{a} | 1.83 × 10^{–19} | 1.25 × 10^{–20} |

γ | 0.351 | 0.351 |

$Kth1$/(m^{3} s^{−1}) | 1.05 × 10^{–16} | 6.33 × 10^{–17} |

Parameters . | At probe P
. | At probe F
. |
---|---|---|

$\u27e8Ttr\u27e9t$/(K) | 10 440 | 710 |

$\u27e8gr\u27e9t$/(m s^{−1}) | 3326.2 | 867.6 |

$\u27e8\sigma \u27e9t$/(m^{2}) | $1.81\xd710\u221219$ | 4.16 × 10^{–19} |

E/(J) _{a} | 1.83 × 10^{–19} | 1.25 × 10^{–20} |

γ | 0.351 | 0.351 |

$Kth1$/(m^{3} s^{−1}) | 1.05 × 10^{–16} | 6.33 × 10^{–17} |

## IV. DYNAMICS OF ENERGY FLUCTUATIONS INSIDE A SHOCK USING THE TWO-ENERGY-BIN MODEL

The numerical solution of the ODE system in Eq. (9) at probes *P* and *F* is shown in Fig. 9 as a function of time and in Fig. 10 as path lines in the phase space of *N _{A}* vs

*N*. Starting with location

_{B}*P*, Fig. 9(a) shows that the solution converges to the non-zero critical values of $NA,crit$ and $NB,crit$, obtained by setting $dNA/dt=dNB/dt=0$, which are listed in Table VI. There are 18.2 and 7.5% differences between the critical values and the average number of DSMC particles per computational cell volume in bins

*A*and

*B*at probe

*P*, which may be attributed to the simplifications of the two-energy-bin model discussed above. Yet, as will be shown, this simple model reveals two orders of magnitude disparity in frequencies at probes

*P*and

*F*. Figure 10(a) shows the streamlines of a vector field defined by the normalized growth rate vector $G=[GNA,GNB]T$ with components,

Solution . | Probe P
. | Probe F
. |
---|---|---|

$NA,crit$ | 1286 | 511.5 |

$NB,crit$ | 657.2 | 495 |

$L1$ | −281 762 + 1 765 303.3 $j$ | 0.0 |

$L2$ | −281 762 − 1 765 303.3 $j$ | −499 976 |

Solution . | Probe P
. | Probe F
. |
---|---|---|

$NA,crit$ | 1286 | 511.5 |

$NB,crit$ | 657.2 | 495 |

$L1$ | −281 762 + 1 765 303.3 $j$ | 0.0 |

$L2$ | −281 762 − 1 765 303.3 $j$ | −499 976 |

We can see from the formation of a stable spiral that the critical point does not depend on the initial values. The figure shows that the streamlines are tangent to the direction of the maximum growth vector at a given location of (*N _{A}*,

*N*) and show eventual decay to the stable critical point of 1245 and 681. The overlaid path line in Fig. 10(a) describes these dynamics for the instantaneous solution shown in Fig. 9(a).

_{B}The Jacobian matrix of Eq. (9) is given by

The eigenvalues of matrix $J(NA,crit,NB,crit)$ are listed in Table VI. At probe *P,* they are complex conjugates with a negative real part in the shock that gives the rate of decay from the initial to critical values. This is consistent with Fig. 9(a), which shows that at 0.01, 0.015, and 0.02 ms, the solution converges to 1%, 0.3%, and 0.07% of the critical values, respectively. These percentage values correspond to an approximate difference between the critical and instantaneous values of the number of *A*-type particles of 13, four, and one, respectively.

More importantly, these results suggest that inside the shock, the number of particles are subjected to recurrent microscopic thermal fluctuations, which always decay with timescales on the order of 0.01–0.02 ms with corresponding frequencies on the order of 100–50 kHz. These frequencies are close to the weighted average of the low-frequency broadband of 37.5 kHz observed in Fig. 2(a) from the PSD analysis. The instantaneous particles in the energy bins *A* and *B*, shown in Fig. 5(c), are a result of superposition of responses to recurrent low-frequency microscopic fluctuations. Therefore, even with the underlying simplifications of the two-energy-bin model, it is able to predict the existence of the low frequencies that originate due to the interactions in the average number of particles between the two modes of the bimodal PDF of particle energies in the presence of translational nonequilibrium in a 1D shock. In addition, the imaginary part of the pair of eigenvalues listed in Table VI at probe *P* also indicates oscillations with a time period of $2\pi /Im(L1)$ (frequency of 281 kHz), which may correspond to fluctuations observed in Fig. 2(a) in the mid-frequency range of ∼100–300 kHz outside of the low-frequency broadband.

Finally, we can relate the two timescales predicted by Eq. (9) in the shock to collision processes *Q _{e}* and

*Q*that describe the relaxation of high energy particles in

_{b}*B*to low energy bin

*A*(see Table II). The frequency associated with process

*Q*is on the order $f\u223cO(keNB)$, or equivalently $f\u223cO(KenB=51\u2009\u2009kHz)$, where $nB=m\u27e8NB\u27e9t=6.81\xd71021$ m

_{e}^{−3}is the average number density of particles in energy bin

*B*. Similarly, the frequency associated with process

*Q*is $f\u223cO(KbnB=284\u2009\u2009kHz)$. Therefore, it is seen that the two timescales present in the shock are a direct consequence of the bimodal energy distribution function. Also, note that although the analysis in this work is for a steady shock, these results hold for unsteady shock as well, because the PDF inside the shock remains bimodal.

_{b}Turning to probe *F* in the freestream, we want to demonstrate that the two-bin energy model is able to predict different dynamics than observed in the shock. Since the flux coefficients *F _{A}* and

*F*and the rate coefficients $k\u0303c$ and

_{B}*k*are zero, we obtain

_{e}which is satisfied by an infinite number of critical points corresponding to different freestream conditions. That is, the number densities are different but the temperature is the same since all solutions have the same shape of the PDF shown in Fig. 3(c), and therefore, the same cross-correlation coefficient is shown in Fig. 5(b). In this case, the solution depends on the initial values of *N _{A}* and

*N*, and converges to a critical point that satisfies Eq. (13). Figures 9(b) and 10(b) show the solution of Eq. (9) as a function of time for a set of initial values of (

_{B}*N*,

_{A}*N*) and a streamplot along with an overlaid path line of the solution, respectively. In Fig. 10(b), the dependence on the initial values can be clearly seen. For the given conditions, the solution converges to critical values of $NA,crit=511.5$ and $NB,crit=495$, which are within 0.5% of the average number of DSMC particles per computational cell volume, of $\u27e8NA\u27e9t=509.9$ and $\u27e8NB\u27e9t=492.8$, and satisfy the ratio of 1.033 given by Eq. (13).

_{B}For the freestream conditions, there are two real eigenvalues, as shown in Table VI with the zero eigenvalue indicating again a non-unique solution of Eq. (13). The negative eigenvalue of -499 976 gives the linear rate of decay, the inverse of which gives the characteristic timescale of 2 *μ*s. This fast decay rate is consistent with the PSD at probe *F*, discussed in Sec. II, which showed that ∼94% of the total spectral energy is distributed in frequencies greater than 1 MHz. More importantly, in the freestream the model predicts the absence of two orders of magnitude lower dominant frequencies having significantly large energy, which appear to be unique to the region of strong nonequilibrium inside a shock.

## V. STROUHAL NUMBER SCALING

With the observation of dominant low-frequency perturbations of macroscopic flow parameters inside the shock, a natural question arises as to whether one can define a nondimensional Strouhal number,

that would be constant for different shock strengths. We propose to define the timescale as that required for the flow to traverse a distance equal to the shock-thickness with the upstream bulk velocity. To justify this timescale, cases ranging from Mach two to 10 were run by varying the upstream bulk velocity but keeping the upstream temperature constant ($Ttr,1$=710 K). Starting with the selection of frequency, *f*, the PSD of the instantaneous pressure spectrum at the location of the maximum gradient of bulk velocity in the shock for different Mach numbers is shown in Fig. 11(a). A broadband of low frequencies is seen, the boundary of which is defined by the inflection point in the NCE, shown in Fig. 11(b).

The characteristic length, *L _{s}*, defined by the density-gradient shock-thickness is calculated as the overall density change divided by the maximum density gradient (e.g., see Ref. 72, Chap. X, Sec. 9) and is shown in Fig. 12(a). We note that

*λ*

_{1}, the upstream mean-free-path, is obtained through the VHS model with a viscosity index of $\omega =0.81$. The SUGAR-1D DSMC shock-thickness values match within 2% with the DSMC calculations of Macrossan and Lilley

^{42}and Bird (see Ref. 7, p. 294) for the same viscosity index. We note that their results are scaled by 76% to account for the differences in

*λ*

_{1}, which they defined based on the hard-sphere model.

^{7}The noticeable discrepancy between the SUGAR-1D DSMC results for $\omega =0.81$ with the experiments of Alsmeyer

^{40}is due to the choice of viscosity index, while qualitatively the variation of shock-thickness with Mach number is consistent with the experimental results. Good agreement is observed between SUGAR-1D DSMC and experiment for a Mach 8 shock simulated with $\omega =0.75$. Based on these values, the calculated ratio of shock-thickness to upstream bulk velocity, $Lsux,1\u22121$, decreases with Mach number, as shown in Fig. 12(b).

The Strouhal number defined based on the above quantities and the weighted average frequency, *f*, shown in Fig. 11(a), is relatively constant within a range of *St *=* *0.007–0.011 and a standard deviation from 0.001 to 0.02. For Mach 8 case simulated with $\omega =0.75$, *St* was found to be in the same range and equal to 0.0092. This range of Strouhal numbers is similar to those corresponding to low-frequency unsteadiness of shock and separation bubble, as observed in the literature for highly compressible flows over embedded bodies, indicating that the two phenomena could be related. For example, in the case of shock-dominated separated flows the Strouhal number associated with low-frequency shock motion ranges from 0.02 to 0.05 (e.g., Refs. 2, 23, 53, and 73–75), assuming that the characteristic length and velocity scales are given by the length of the separation bubble and the upstream bulk velocity, respectively. In the study of oblique SBLIs, Nichols *et al.*^{76} defined the Strouhal number based on the upstream boundary-layer thickness and freestream velocity and found it to be within a range of 0.0003–0.05. Therefore, we hypothesize that these fluctuations may play a major role in shock-dominated flows having shock-thickness comparable to other important length scales in the flow, such as the size of the boundary layer in SBLIs and turbulent length scale in shock-turbulence interaction.

To test whether the established range of Strouhal numbers holds for variations in upstream temperature, another set of DSMC cases were simulated in which the upstream temperature was changed from 710 K by fractions of 1/8, 1/4, 1/2, and 2, while the upstream bulk velocity was varied accordingly to maintain a constant Mach number of 7.2. Figure 13 shows a decrease in the ratio of $Lsux,1\u22121$ with increase in temperature, due primarily to the increase in $ux,1$ rather than the shock-thickness. The Strouhal numbers obtained using the aforementioned approach are within the range of St = 0.005–0.011, with a small reduction with decrease in temperature. Nonetheless, the ±1 standard deviation of broadband frequencies is contained within the limits of *St *=* *0.001 to 0.02. Appendix A also shows that the Strouhal number remained in the same range for a case run by increasing the freestream number density, *n*_{1}, by three orders of magnitude, where the mean-collision-time of molecules decreased, frequency increased, and shock-thickness decreased by the same order. In Ref. 68, we show that a semi-empirical relationship to predict low-frequency oscillations in shocks can be derived from bimodal probability distribution functions based on noncentral Chi-squared functions for a broad range of Mach numbers and freestream temperatures.

Finally, our Strouhal number formulation was used to qualitatively explain the ∼300 kHz frequencies observed inside a shock in the experiments of Chynoweth *et al.* (see Ref. 77, Fig. 7) of a Mach 6 air flow and $Re1=6.73\xd7106$ m^{−1} over a 7° half-angle cone attached with a 10° fin. At freestream conditions experienced by the fin downstream of the leading-edge shock, that is, $\rho 1=4.636\xd710\u22122$ kg m^{−3}, Mach 5, and freestream temperature, $Ttr,1$=71.67 K, the normal shock-thickness, $Ls\u223c$ 6.6 *μ*m, was estimated from Fig. 12(a) and *λ*_{1} calculated using the VHS molecular model. Using Eq. (14) and estimating *St* from Fig. 12(b), the low-frequency broadband was predicted to have a weighted average frequency of 1.2 MHz with one standard deviation limit from 0.4 to 2.1 MHz, which was also verified from a separate DSMC simulation. The broadband from the simulation also revealed high spectral energy near 50 kHz frequency, which is also present in the experimental PSD spectrum. This estimate is qualitatively reasonable, even though it ignores factors that could reduce the range of dominant frequencies, such as differences in gas properties, reduced strength of the oblique shock generated by the fin compared to the normal shock, and surface effects.

## VI. CONCLUSION

The investigation of macroscopic fluctuations in a DSMC-computed Mach 7.2 shock layer revealed two orders of magnitude lower dominant frequencies of fluctuations than found in the freestream. In the peak overshoot region in the shock, a low-frequency broadband with a weighted-average of ∼37 kHz was observed as opposed to the freestream region where peak frequencies are concentrated at ∼3.5 MHz corresponding to the inverse mean-collision-time. These disparities were attributed to the differences in particle distribution functions in the nonequilibrium zone of shock vs the equilibrium regime upstream. The fluctuations in the *x*-directional normalized overall stress component, which is the mean of the PDF $f\xi x$, were found to be correlated with perturbations in the entire *X*-directional energy space of PDF. Two distinct energy bins were identified depending on the sign of the cross-correlation coefficient, and a Lotka–Volterra type two-energy bin ODE model was constructed.

The model accounted for the interaction of two energy bins through intermolecular collisions and particle fluxes from neighboring computational cells. At a location inside the shock, the model predicted two disparate timescales: a longer timescale (50–100 kHz) associated with the time of decay of fluctuations in the number of particles in energy bins to critical values and an order of magnitude smaller timescale (281 kHz) associated with oscillations in the number of particles. In the freestream, the model predicted faster decay of fluctuations, $O(MHz)$, and more importantly, an absence of low-frequency fluctuations consistent with the DSMC spectral analysis.

Finally, a Strouhal number was defined based on the density gradient shock-thickness and upstream bulk velocity to nondimensionalize the low-frequency broadband characterized by a weighted average frequency and ±1 standard deviation across a wide range of Mach numbers from 2 to 10. The Strouhal number was found to range from 0.007 to 0.011 with ±1 standard deviation between 0.001 and 0.02, consistent with values found in the literature corresponding to low-frequency unsteadiness of shock and separation bubbles in shock-dominated flows. The range of Strouhal numbers was also found to hold for a Mach 7.2 shock simulated with different upstream temperatures and number density, and explain qualitatively the frequencies inside a shock observed in the work of Chynoweth *et al.*^{77}

The presence of low-frequency fluctuations of shock layers suggests that these disturbances may play a key role in the receptivity process of transition in hypersonic flows, SBLIs, and shock–turbulence interactions, especially when the shock-thickness is comparable to other important length scales in the flow. Our work implies that the leading-edge shock interactions with freestream disturbances, such as tunnel-induced turbulence, can result in selective amplification of disturbances in the low-frequency range, and consequent excitation of the downstream boundary layer. Future numerical investigations focused on understanding the effect of freestream disturbances on transition should not ignore shock-induced perturbations in the boundary layer within the range of Strouhal numbers established in this work.

## ACKNOWLEDGMENTS

The research conducted in this paper is supported by the Office of Naval Research under Grant No. N000141202195 titled “Multi-scale modelling of unsteady shock-boundary layer hypersonic flow instabilities,” with Dr. Eric Marineau as the program officer.

This work used the STAMPEDE2 supercomputing resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE) at the Texas Advanced Computing Center (TACC) through Allocation No. TG-PHY160006. In addition, S.S.S. would like to thank colleague Nakul Nuwal for helpful discussions regarding the two-energy bin model.

## DATA AVAILABILITY

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

### APPENDIX A: ADDITIONAL NUMERICAL TESTS

First, we discuss the frequencies associated with potential disturbances introduced due to the STABIL algorithm for the Mach 7.2 shock case described in Sec. II B. The STABIL algorithm works as follows: when the number of particles in the domain boundary increases by a small, pre-set number, then all particles are translated by a small distance, *δ _{x}*, in the direction opposite to the drift of the shock. After translation, molecules that move outside of the boundary are deleted and those that are within a

*δ*distance from the opposite boundary are replicated. One-way spurious numerical frequencies could be introduced at the boundaries is due to the recurring frequency at which additional particles are added within a

_{x}*δ*distance from the boundaries to stabilize the shock. Starting with the upstream side, it was found that particles were not added near the upstream boundary during stabilization. On the downstream side, Fig. 14 shows the number of particles added during stabilization as a function of time and time passed since the last instance of particle addition. The mean values of these quantities are 236 (or ∼6.2% of the downstream computational particles in a root cell) and 0.0018 ms, respectively. We note that mean time between adding particles corresponds to a frequency of ∼553 kHz, which is much higher than the frequency of low-frequency fluctuations that are discussed in Sec. II. It is possible that the disturbances generated due to these particle additions could propagate to the shock at the downstream speed of sound, but the frequency associated with such disturbances would be 205 kHz, which is again outside of the low-frequency broadband. Therefore, this analysis supports the conclusion that molecular frequencies at the peak over shoot of the shock are not numerical artifacts of the STABIL algorithm.

_{x}Next, we analyze the possible effect of *numerical* amplification of fluctuations on our low-frequency predictions caused by the stochastic nature of DSMC for $Fn>1$. The best way to assess this is to study its effect on the PSD, as shown in Fig. 15. The DSMC simulation, discussed in Sec. II B, was re-run first with *F _{n}* reduced by four orders of magnitude, which now requires 33 million particles, 1.2 × 10

^{6}time steps with $\Delta t$ = 1 ns, for the same number density of 10

^{22}m

^{−3}. An additional DSMC simulation was also performed for the case of

*F*= 1 (i.e

_{n}*.,*no numerical amplification of fluctuation amplitudes) with the freestream number density now raised by three orders of magnitude, $n1=1025$ m

^{−3}. We note that we have changed the abscissa from frequency to Strouhal number using Eq. (14) because when we change number density we do not expect the “line-by-line” PSD spectrum to remain the same (due to the molecular nature of the flow), but we do expect the mean Strouhal numbers to be similar if we are looking at a phenomena related to fluctuations induced by the motion of the shock and not computational noise. A comparison of the new PSD spectra in the shock with the original case for the 10

^{22}m

^{−3}upstream condition is shown in Fig. 15(a) where indeed we find only a change of 2.7% in the weighted average frequency of the low-frequency broadband in the shock such that the mean Strouhal number is $St\u223c0.012$. Finally, Fig. 15(b) shows the PSD obtained from the simulation run by increasing the freestream number density by three orders of magnitude, $n1=1025$ m

^{−3}and

*F*= 1. Comparison of this result with Fig. 15(a) shows that both conditions generate PSDs consistent with the Strouhal number scaling discussed in Sec. V and produce similar quantities such as mean Strouhal number and broadband limits.

_{n}### APPENDIX B: SIMPLIFICATION OF COLLISION PROCESSES

We first describe the grouping of detailed collision processes *p _{i}* into like-

*P*processes of Table II based on the energy density functions at locations

_{i}*P*and

*F*. Starting with the reduction of collision processes

*p*to like-

_{i}*P*processes, Fig. 8 shows that the energy distributions of

_{i}*A*and $Ad\u2032$ are very similar at probes

_{c}*P*and

*F*, so that $Ad\u2032$-type particles can be considered to be the same as those of

*A*. In addition, although the distributions $f\xi xA$ and $f\xi y+\xi zA$ of $Ac\u2032$ and

_{c}*A*particles are the same at probe

_{d}*F*, they are slightly different at probe

*P*at high

*ξ*and low $\xi y+\xi z$ energies. This occurs because when

_{x}*A*- and

_{c}*B*-type particles collide, the

*B*particle loses its

*X*-directional energy, which results in the increase in not only the

*ξ*of the

_{x}*A*-type particle but also a noticeable increase in its $\xi y+\xi z$ energy, as can be seen in the transverse energy distributions of $Ac\u2032$-type particles in Fig. 8(b) at probe

_{c}*P*. Also, when two

*A*particles collide, one of them is switched from energy bin

_{d}*A*to bin

*B*at the expense of higher transverse energy of the other

*A*particle that remains in bin

_{d}*A*. In the freestream, since these two

*A*-type particles are the same $Ac\u2032$-type particles that were generated from process

_{d}*p*, their transverse energy distributions exactly match, as seen in Fig. 8(d). However, at probe

_{c}*P*, it is also possible that only one of the two

*A*particles is an $Ac\u2032$-type, and the other one has lower transverse and higher

_{d}*X*-directional energy than the $Ac\u2032$-type particle, as can be seen in the

*slightly*larger

*X*-directional energy distribution of

*A*-type particles than the $Ac\u2032$-type in Fig. 8(a) between $0.5<\xi x<1.33$. We demonstrate that small distinctions such as this do not affect the dynamics of the system and assume that $Ac\u2032\u223cAd$ to construct processes

_{d}*P*and

_{c}*P*given in Table II.

_{d}Similarly, we can describe the process *p _{e}*, where the particles of type-$Ae\u2032$ have a very different transverse energy distribution shown in Fig. 8(b) than the other type-

*A*particles at probe

*P*. This process does not occur at probe

*F*. Ignoring the differences at low transverse energies, the distribution of $Ae\u2032$ and

*A*particles is the same, as shown in Fig. 8(b), and therefore, they can be denoted by the same identifier $Ae\u2032$.

_{f}We move on to present the derivations for the simplifications to Eq. (8) that allow us to reduce the six collision processes (*P _{i}*) given in Table II to four processes (

*Q*). To construct compound collision processes,

_{i}*Q*we start with the rate coefficient

_{i}*k*for process

_{i}*P*, which can be evaluated from the DSMC simulation as

_{i}where the collision takes place with between particle types *R*_{1} and *R*_{2}. $\u27e8Ci\u27e9t$ for collision processes defined in Table II are listed in Table IV. Using the definition of rates *k _{a}* and

*k*, the first kinetic equation in Eq. (7) can be written as

_{b}where we have used the fact that

as shown in Fig. 16(a). The ratio $NAa/\u27e8NAa\u27e9t$ is estimated from Fig. 8(a) by obtaining the number of particles between $0.0<\xi x<1.33$ at probe *P*. At probe *F*, the ratio of $NAa/\u27e8NA\u27e9t$ cannot be directly estimated from Fig. 8(c) because the distribution of *A _{a}*-type particles overlaps with other

*A*-type particles. Yet, Eq. (B3) is assumed to hold true under the assumption that the entire

*X*-directional energy zone $\u221228<\xi x<0$ is coarsened into a single bin

*A*.

Using the definitions of rates *k _{c}* and

*k*, the second kinetic equation in Eq. (7) can be written as

_{d}Since the average number of $Ac\u2032$ particles is not known, it is difficult to estimate the second term on the right-hand side. However, this term can be simplified by observing from Fig. 16(b) that the ratio of rates $Cc/Cd$ does not vary significantly from the average rates at probe *P*, that is,

and the standard deviation in the fluctuation of the ratio of $Cc/Cd$ is only 3.7% of the average ratio of 3.63 for probe *P* and 4.2% of an average ratio of one for probe *F*. Using Eq. (B1) for *C _{c}* and

*C*in B5, and the definition of

_{d}*k*based on the average rate, we obtain

_{c}In addition, we can see from Figs. 8(a) and 8(b) that the *A _{c}* particles have the same energy distributions as the

*A*at probe

_{h}*P*, which implies that any particle that takes part in process

*P*can take part in process

_{c}*P*, and together, they constitute 72.37% of total collisions in which the particles in energy bin

_{h}*A*take part. Therefore, we can expect the changes in the ratio $NAc/\u27e8NAc\u27e9t$ are correlated with the changes in ratio $NA/\u27e8NA\u27e9t$ and write

The same assumption holds at probe *F* because the *A _{c}*-type particles constitute the majority of particles from $\u221220<\xi x<\u22125$, which is a significant portion of energy bin

*A*($\u221228<\xi x<0$), as seen from Fig. 8(c).

By using the aforementioned assumptions, the third kinetic equation in Eq. (7) becomes

We also drop the second term on the right-hand side of the fourth kinetic equation in Eq. (7) because $\u27e8Ce\u27e9t\u226b\u27e8Cf\u27e9t$. Note that both $\u27e8Ce\u27e9t=\u27e8Cf\u27e9t=0$ at probe *F*.