^†^†thanks: https://doi.org/10.1103/PhysRevE.107.015104

Data-Driven Constitutive Relation Reveals Scaling Law for Hydrodynamic Transport Coefficients

Candi Zheng Corresponding author: [email protected] Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Xueyuan Rd 1088, Shenzhen, China Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong SAR, China Yang Wang [email protected] Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong SAR, China Shiyi Chen [email protected] Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Xueyuan Rd 1088, Shenzhen, China

Abstract

Finding extended hydrodynamics equations valid from the dense gas region to the rarefied gas region remains a great challenge. The key to success is to obtain accurate constitutive relations for stress and heat flux. Data-driven models offer a new phenomenological approach to learning constitutive relations from data. Such models enable complex constitutive relations that extend Newton’s law of viscosity and Fourier’s law of heat conduction by regression on higher derivatives. However, the choices of derivatives in these models are ad-hoc without a clear physical explanation. We investigated data-driven models theoretically on a linear system. We argue that these models are equivalent to non-linear length scale scaling laws of transport coefficients. The equivalence to scaling laws justified the physical plausibility and revealed the limitation of data-driven models. Our argument also points out that modeling the scaling law could avoid practical difficulties in data-driven models like derivative estimation and variable selection on noisy data. We further proposed a constitutive relation model based on scaling law and tested it on the calculation of Rayleigh scattering spectra. The result shows our data-driven model has a clear advantage over the Chapman-Enskog expansion and moment methods.

^†^†preprint: APS/123-QED

I Introduction

Multiscale physics is widely encountered in fluid dynamics [1], soft matter systems [2], and quantum chemistry [3]. One of the typical multiscale physics problems is rarefied gas dynamics [4]. Rarefied gas flow simulation is known to be difficult due to the non-negligible dynamics at the mesoscopic scale. Simulation resolving these scales is computationally expensive for continuous and transitional flows, such as the Direct Simulation Monte Carlo (DSMC) [5] method. Instead, extended hydrodynamics equations at a coarse-grained macroscopic scale are efficient substitutes to reduce the computational cost. What lies within the heart of extended hydrodynamics is constitutive relations. Constitutive relations summarize mesoscopic scale dynamics as macroscopic phenomena, such as viscosity and heat conduction. Traditionally they are modeled by perturbation or polynomial expansion around the equilibrium of dense gas. The perturbation models, such as the Burnett-type equations [6, 7], utilize high-order spatial derivatives according to the Hilbert-Chapman-Enskog expansion [8, 9, 10]. However, difficulties exist in stability [11] and non-guaranteed convergence [12]. The polynomial expansion methods are represented by the Grad moment method [13] and its extensions [14] modifying the equilibrium distribution with orthogonal polynomials. Nevertheless, issues appear in unphysical solutions [15] and hyperbolicity [16]. The traditional methods’ perturbation or expansion nature limited their viability near dense equilibrium, constraining their applicable Knudsen numbers range [14].

Data-Driven models offer a new phenomenological approach to obtaining machine-learned constitutive relations from data. Compared to perturbation or expansion from existing theory, data provids an alternative source of information and is expected to expand the applicable range of extended hydrodynamics equations [17]. There have been attempts to learn constitutive relations from mesoscopic results [18] or to find proper moment equations [17]. Data-Driven models are also used in related areas such as learning the unknown governing of physical systems [19, 20, 21, 22], simulating physical dynamics [23, 24, 25] and solving the Boltzmann equation [26]. These attempts have proven the concept of data-driven modeling. However, the advantage over traditional models like Chapman-Enskog and Grad moment method hasn’t been established yet. Limitations for data-driven models include derivative estimation [27], determining input quantities (variable selection) [21], and modeling across a range of Knudsen numbers. Besides, the rather ad-hoc linear or neural network regression in data-driven models lack a clear physical explanation.

In this paper, we seek the physical explanation of data-driven models by investigating linear systems. We focus on the conservation laws and analyze data-driven constitutive relation models that extend Newton’s viscosity law and Fourier’s heat conduction law. We argue that these linear models are equivalent to non-linear length scale scaling laws of viscosity and heat conduction coefficients. These length scale scaling laws describe the change of viscosity and heat conduction coefficients, as we concern with dynamics at different length scales described by Knudsen numbers. The equivalence between data-driven constitutive relations and scaling laws justified the physical plausibility of data-driven models.

Based on our argument, we suggest modeling scaling laws explicitly in data-driven models. In doing so, we could involve high-order derivatives implicitly in constitutive relations without calculating them. It helps to avoid practical difficulties in data-driven models like derivative estimation and variable selection. We further modeled the constitutive relation based on our suggestion.

We apply our model to calculate the Rayleigh scattering spectra as the numerical benchmark. The Rayleigh scattering have been well studied [28] and used in Lidar wind measurement [29]. However, it remains difficult to correctly model the spectra shape in the transition region for today’s extended-hydrodynamic equations [30]. The numerical results show that our data-driven model can capture the spectra shape at the transition region. To our knowledge, it is the first time that the data-driven hydrodynamic model significantly outperforms the traditional Chapman-Enskog expansion and Grad moment methods.

II Methods

We consider the linearized extended hydrodynamics for one-dimensional homogeneous rarefied ideal gas. The hydrodynamics equations govern the dynamics of gas. The most important hydrodynamics equations are mass, momentum, and energy conservation laws. They form a one-dimensional (1D) linear system of density $\rho$ , velocity $v$ , and temperature $T$ , respectively. The non-dimensionalized linear system for conservation laws is

\begin{split}\frac{\partial{\rho}}{\partial{t}}+\frac{\partial{v}}{\partial{x}}&=0\\ \frac{\partial{v}}{\partial{t}}+\frac{\partial{T}}{\partial{x}}+\frac{\partial{\rho}}{\partial{x}}&=-\mbox{Kn}\frac{\partial{\sigma}}{\partial{x}}\\ \frac{3}{2}\frac{\partial{T}}{\partial{t}}+\frac{\partial{\rho}}{\partial{x}}&=-\frac{15}{4}\mbox{Kn}\frac{\partial{q}}{\partial{x}},\end{split}

(1)

in which Kn is the Knudsen number describing how rarefied the gas is. Detailed descriptions of non-dimensionalization and the definition of the Knudsen number are in Appendix A.

However, the hydrodynamics equations are not closed with two extra unknown terms: the stress ${\sigma}$ and the heat flux ${q}$ that encodes the mesoscopic dynamics. To close the equations, constitutive relations that model the stress and the heat flux with known quantities are necessary.

II.1 Data-Driven Constitutive Relations and Its Equivalence with Scaling Laws

We adopt a general form of data-driven constitutive relations consisting of derivatives of various orders similar to other data-driven models for physical systems [19, 20, 21, 22, 23, 24]. It is also motivated by the Hilbert-Chapman-Enskog expansion. The Hilbert-Chapman-Enskog expansion is a systematic way to generate constitutive relations for conservation laws at small Knudsen numbers. The leading order of expansion yields the well-known Newton’s law of viscosity and Fourier’s law of heat conduction and defines the viscosity coefficient $\mu_{0}$ and heat conduction coefficient $\kappa_{0}$ . However, they are not valid for rarefied gas effects at large Knudsen numbers [31]. For large Knudsen numbers, higher-order expansions extend the capability of constitutive relations by incorporating high-order spatial derivatives of density, velocity, and temperature. If we consider linear systems, these high-order spatial derivatives are combined linearly by coefficients determined by the Hilbert-Chapman-Enskog expansion. However, the Hilbert-Chapman-Enskog expansion guarantee neither convergence nor stability of the system [6]. Similar to the Hilbert-Chapman-Enskog expansion, we consider constitutive relations linear combinations of high-order spatial derivatives. However, we aim to determine combinations coefficients via a data-driven regression approach. Therefore we adopt the following general form of the data-driven constitutive relation

\begin{split}{\sigma}&=-\sum_{n=1}^{\infty}\left({a}_{n}\frac{\partial^{n}{v}}{\partial{x}^{n}}+{c}_{n}\frac{\partial^{n}{\rho}}{\partial{x}^{n}}+{e}_{n}\frac{\partial^{n}{T}}{\partial{x}^{n}}\right)\\ {q}&=-\sum_{n=1}^{\infty}\left({b}_{n}\frac{\partial^{n}{T}}{\partial{x}^{n}}+{d}_{n}\frac{\partial^{n}{\rho}}{\partial{x}^{n}}+{f}_{n}\frac{\partial^{n}{v}}{\partial{x}^{n}}\right),\end{split}

(2)

where $x$ is the non-dimensional spatial coordinate, $a_{n},b_{n},c_{n},d_{n},e_{n},f_{n}$ are unknown regression coefficients. The constitutive relation (2) has the same functional form obtained from Hilbert-Chapman-Enskog expansion since both are combinations of high-order spatial derivatives. But the regression coefficients in (2) are to be determined via a data-driven approach.

There are practical difficulties in directly applying constitutive relation (2). First is the problem in variable selection. This problem arises because we only have limited data in practice to determine the infinitely many regression coefficients in (2). Consequently, we could only determine a selected subset of regression coefficients. Choosing the best subset of regression coefficients is a challenging variable selection problem we wish to avoid. The second problem is density estimation. Constitutive relation (2) contains high-order spatial derivatives, which are difficult to estimate in practice. A naive attempt at estimating high-order spatial derivatives using the finite difference method requires a highly dense mesh and is very sensitive to noise. It completely fails on data generated by the DSMC method since they contain strong statistical noise. Finally, constitutive relation (2) does not guarantee the stability of hydrodynamic equations. The reason is there are no constraints on entropy production yet to respect the second law of thermodynamics. Fortunately, it turns out that reformulating the problem in the Fourier space with proper constraints on entropy production enables us to bypass the practical difficulties in variable selection and derivative estimation.

Now we reformulate the constitutive relations with the help of the Fourier transform and entropy production constraints. Fourier transform allows us to convert the derivatives in constitutive relations into algebraic expressions. Meanwhile, constraints on entropy production eliminate undesired terms and imaginary parts that appear in the Fourier transformation. The outline of the reformulation is as follows: Firstly, the entropy constraint reduces the constitution relations to the form that stress $\sigma$ consists of only velocity derivatives and heat flux $q$ consists of only temperature derivatives. This is because the stress and velocity, the same as heat flux and temperature, must be correlated to produce non-increasing entropy as follows.

\begin{split}\dot{s}&=-\frac{q}{T^{2}}\partial_{x}T-\frac{\sigma}{T}\partial_{x}v,\end{split}

(3)

where $\dot{s}$ is the entropy change rate per volume [32]. The only possibility is stress depends on velocity only, the same as heat flux depending on temperature, since density, velocity, and temperature are statistically independent [33]. Secondly, the non-increasing constraint on the entropy production eliminates undesired imaginary parts in the Fourier transform of constitutive relations. This constraint requires that each Fourier mode of the density, velocity, and temperature must produce non-negative entropy. It is necessary if we wish the linear system to be stable. As a result, constitutive relations are expressed as a summation of infinite polynomial series in the Fourier space. Finally, collecting and reformulating the summation in the constitutive relations leads to the following constitutive relations in the Fourier space

\begin{split}\tilde{\sigma}(k)&=-ik\frac{4}{3}\frac{\mu(k)}{\mu_{0}}\tilde{v}(k);\quad\mu(k)\geq 0\\ \tilde{q}(k)&=-ik\frac{\kappa(k)}{\kappa_{0}}\tilde{T}(k);\quad\kappa(k)\geq 0,\end{split}

(4)

where $\mu_{0}$ is the viscosity coefficient, $\kappa_{0}$ is the heat conduction coefficient, $k$ is the non-dimensional wavenumber for each Fourier mode, $\tilde{\sigma}(k),\tilde{q}(k),\tilde{v}(k),\tilde{T}(k)$ are corresponding spatial Fourier transforms of $\sigma,q,v,T$ . The detailed derivation from the constitutive relation (2) to (4) are shown in Appdendix B. The derivation also shows that the functions $\mu(k),\kappa(k)$ are even functions satisfy the natural constraints $\mu_{0}=\lim_{k\rightarrow 0}\mu(k)$ and $\kappa_{0}=\lim_{k\rightarrow 0}\kappa(k)$ . As we will discuss later, they describe the length scaling law for viscosity and heat conduction coefficient. Therefore we have shown that the data-driven constitutive relation (2) transformed into the form (4) containing scaling laws under the constraint of non-increasing entropy. This established the equivalence between data-driven constitutive relations and scaling laws.

The functions $\mu(k),\kappa(k)$ are the length scaling laws of viscosity and heat conduction coefficients. They describe the relative change of viscosity and heat conduction coefficients w.r.t length scale changes of the system. This is because $k$ is closely related to the Knudsen number $\mbox{Kn}=l/L$ that characterize the length scale of a rarefied gas system, in which $l$ is the mean free path of gas molecules and $L$ is the representative length scale of the system. Particularly, as defined in Appendix A, the Knudsen number of a Fourier mode is proportional to its non-dimensionalized wavenumber $k$

\mathrm{Kn}\propto|k|.

(5)

Therefore the even functions $\mu(k),\kappa(k)$ are also functions of the Knudsen number, hence are length scaling laws. These scaling laws could be measured experimentally [34]. However, we can not use such experimental results directly because the definition of the Knudsen number is not unified but varies according to the experiment setting. Alternatively, scaling laws could be learned through a data-driven approach from data like fluctuation spectra [35] containing information on viscosity and heat conduction.

Scaling laws $\mu(k),\kappa(k)$ in (4) are much easier to be determined than regression coefficients in (2). These coefficients may lead to divergence at large Knudsen numbers, making it ill-conditioned to determine regression coefficients valid for large Knudsen numbers. Instead, we could learn scaling laws $\mu(k),\kappa(k)$ uniformly from data at various Knudsen numbers without worrying about convergence. Learning scaling laws also eliminates the demand in variable selection, which refers to choosing a subset of regression coefficients. It is because all regression coefficients are now summarized in the function $\mu(k),\kappa(k)$ . Moreover, learning scaling laws is robust against noisy data since it avoids using estimated derivatives in constitutive relations (2). Therefore learning the scaling laws, compared to regression coefficients in (2), avoid practical difficulties in convergence issue, variable selection, and derivative estimation.

II.2 Modeling Scaling Laws using Neural Network

One difficulty that remains is learning scaling laws from data turns out to be a non-convex optimization problem that is difficult to solve. We overcome it by approximating scaling laws $\mu(k),\kappa(k)$ using neural networks, taking advantage of their stochastic optimization technique designed for non-convex optimizations [36].

Neural network modeling functions $\mu,\kappa$ must be constrained to obtain correct asymptotics and symmetry for hydrodynamics. Asymptotically, the function values of $\mu,\kappa$ must be specified to the equilibrium values $\mu_{0},\kappa_{0}$ at $\mathrm{Kn}=0$ to guarantee the constitutive relation’s consistency with the Navier-Stokes equation. In addition, we couple $\kappa$ and $\mu$ together

\kappa(k)=\frac{5k_{B}}{2m}\frac{\mu(k)}{\mbox{Pr}},\quad\mbox{Pr}=\frac{2}{3}

(6)

to constraint the Prandtl number Pr to the Chapman-Enskog result $\mbox{Pr}=\frac{2}{3}$ . While this coupling is not necessary, we find it accelerates the learning process without undermining the accuracy in practice. As for symmetry, homogeneity in space also requires the scaling laws to be even functions of $k$ . Homogeneity means there is no preferred direction in space. Hence the direction in space coordinate or the corresponding wavenumber $k$ should not make a difference in the scaling laws. For the one-dimensional case, the direction of k is its sign. Therefore, the scaling laws must be even functions independent of the sign of $k$ .

To satisfy all these constraints, we design the following non-dimensional constrained neural network for $\mu$ satisfying $M(\mathrm{Kn})=\frac{4}{3}\frac{\mu(\mathrm{Kn})}{\mu_{0}}$ with the architecture

\begin{split}M(\mathrm{Kn})&=\frac{4}{3}\left(1+\mathbf{W}_{2}\cdot\mbox{Tanh}(\mathbf{W}_{1}\cdot\mathcal{H}(\mbox{20Kn}))\right)\\ \mathcal{H}(x)&=\frac{1}{2}x^{2},x<1;\quad|x|-\frac{1}{2},x\geq 1,\end{split}

(7)

in which Kn are proportional to $|k|$ as in (5), $\mathbf{W}_{1},\mathbf{W}_{2}$ are the one-dimension weight vectors of the neural network, with the activation function Tanh acting element-wise on the vector input. The function $M(\mbox{Kn})$ is even and satisfies $M(0)=\frac{4}{3}$ and $M^{\prime}(0)=0$ . It guarantee the consistency with the NS equation. With the modeling of the scaling laws, we are prepared to investigate the capability of scaling laws in describing rarefied gas dynamics.

Refer to caption — Figure 1: Rayleigh scattering of electromagnetic(EM) waves with wave vector $\mathbb{k}_{i}$ and frequency $\omega_{i}$ through gas with density fluctuation $\delta\rho$ . $\hat{\mathbb{n}}$ is a unit vector. The intensity $I(\omega)$ of scattered EM wave with frequency shift $\omega$ is the Rayleigh scattering spectra. The Rayleigh scattering spectra are proportional to and determined by the density fluctuation spectra $\left<\delta\rho^{2}\right>$ . Therefore to compute the Rayleigh scattering spectra, we need only to compute the density fluctuation spectra of gas.

II.3 Rayleigh Scattering as Benchmark case

We will test the capability of scaling laws $\mu,\kappa$ in describing rarefied gas dynamics by calculating the Rayleigh scattering spectra. The Rayleigh scattering describes the refraction of electromagnetic waves (EM waves) passing through media with stochastic density fluctuation [37, 38, 39]. Such fluctuation usually appears as density fluctuation waves and happens spontaneously with the thermal motion of gas molecules. The Rayleigh scattering spectra are defined as the intensities $I(\omega)$ of scattered EM waves after the Rayleigh scattering with frequency shifts $\omega$ , as shown in Fig 1. They are proportional to the density fluctuation spectra of gas

I(\omega)\propto\left<\delta\rho^{2}\right>(k(\omega),\omega),

(8)

where $k$ is the wavenumber change of the scattered EM wave determined observation position and incident wave frequency, $\left<\delta\rho^{2}\right>$ is the density fluctuation spectra, which describes the intensity of density fluctuation waves at each wavenumber $k$ and frequency $\omega$ . A detailed description on the relation between the Rayleigh scattering spectra and density fluctuation spectra is shown in Appendix C. As a consequence of the proportionality between the Rayleigh scattering spectra and density fluctuation spectra, calculating the Rayleigh scattering spectra only needs to compute the density fluctuation spectra of the gas media.

The density fluctuation spectra $\left<{\delta\rho}^{2}\right>(k,\omega)$ describe the amplitude of density fluctuation waves caused by the collective motions of gas molecules. The wavenumber $k$ in the density fluctuation spectra specifies the wavelength of density fluctuation waves. It also sets the Knudsen number of density fluctuation waves since the wavenumber $k$ is proportional to the Knudsen number. Given a Knudsen number by specifying $k$ , the spectra $\left<{\delta\rho}^{2}\right>(k,\omega)$ could be calculated from macroscopic governing equations (1) using constitutive relation (4) (details in Appendix D). Hence the values of scaling laws $\mu,\kappa$ in the constitutive relation affect the shape of spectra $\left<{\delta\rho}^{2}\right>$ as a function of $\omega$ . It means that the density fluctuation spectra contain information on scaling laws which we aim to extract by training the neural network models. In practice, we train the neural network modeling scaling laws on density fluctuation spectra data $\left<\delta\rho^{2}\right>_{dsmc}$ computed by the DSMC method (Appendix E).

The density fluctuation spectra are not enough to confirm the capability of scaling laws in describing rarefied gas dynamics. It is because there is the risk of overfitting. Overfitting refers to the neural network learning the scaling law by rote from density fluctuation spectra. In other words, the neural network learns a scaling law that fails in predicting quantities other than density fluctuation spectra. To eliminate the risk of overfitting, we need to prepare test data to examine the neural network’s generalization ability: the ability to predict quantities that the neural network has not seen in the training process.

We examine the generalization ability of the neural network on test data consisting of velocity fluctuation spectra $\left<{v}^{2}\right>(k,\omega)$ . Similar to density fluctuation spectra, velocity fluctuation spectra describe the amplitude of velocity fluctuation waves caused by the collective motions of gas molecules. Velocity fluctuation spectra serve as ideal test data because of the following reasons: first, velocity fluctuation is consistent with the scaling law discussed in our paper since velocity fluctuation also obeys the hydrodynamic equations (1); second, velocity fluctuation corresponds to a different physical scenario compared to density fluctuation. In detail, velocity fluctuations are solved from the hydrodynamic equations with an initial condition ( (58) in Appendix D ) completely different from density fluctuation ( (50) in Appendix D ). The ‘consistent but different’ characteristic of velocity fluctuations makes them ideal for examing the generalization ability of our neural-network-modeled scaling laws.

II.4 Training the Neural Network on Density Fluctuation Spectra Data

We train the neural network models for scaling laws on the density fluctuation spectra data $\left<{\delta\rho}^{2}\right>$ . Specifically, it refers to learning the weight vectors $\mathbf{W}_{1},\mathbf{W}_{2}$ in the neural network (7) from data. This requires a loss function as the learning target. In our case, the loss function compares the difference between the observed spectra $\left<\delta\rho^{2}\right>_{dsmc}$ and predicted spectra $\left<\delta\rho^{2}\right>$ . The former are training data obtained from the DSMC computation (Appendix E), while the latter are the predictions of the governing equation. We define the loss function for any input weight vector $\mathbf{W}=\mathbf{W}_{1},\mathbf{W}_{2}$ as

\begin{split}L(\mathbf{W})&=\mathbb{E}_{\mathrm{Kn}\sim U}\mathbb{E}_{\omega\sim p(\omega|\mathrm{Kn})}\\ &\left|\left<\delta\rho^{2}\right>_{dsmc}(\mathrm{Kn},\omega)-\left<\delta\rho^{2}\right>(\mathrm{Kn},\omega;\mathbf{W}))\right|^{2},\end{split}

(9)

in which the predicted spectra $\left<\delta\rho^{2}\right>$ is a function on the weight vectors $\mathbf{W}$ because it depends on the neural network $M(\mathrm{Kn})$ . The symbol $\mathbb{E}_{\mathrm{Kn}\sim U}$ represents taking the expectation numerically by sampling $\mathrm{Kn}$ from a uniform distribution $U$ . Meanwhile, $\mathbb{E}_{\omega\sim p(\omega|\mathrm{Kn})}$ represents taking the expectation by sampling $\omega$ from a conditional distribution $p(\omega|k)$ , which is proportional to the amplitude of the DSMC spectra. Sampling $\omega$ in this way makes the sample point lies more in the peak region. After defining the loss function, we use the ADAM [40] optimizer to minimize the loss function and determine the weight vectors (Appendix F).

We take extra caution on the finite domain effect of the DSMC computed spectra data. The DSMC simulates gas confined in 1-D space of finite domain length $L_{d}$ . However, we aim to compute the density fluctuation spectra for Fourier modes with infinite spatial span. Therefore, finite domain length inevitably affects the spectra, especially for Fourier modes with a wavelength comparable to domain length. Such finite domain effect is proportional to the mean free path and domain length ratio $\frac{l}{L_{d}}$ which vanishes as $L_{d}$ tends to infinity. As a solution, we use a large domain length much greater ( $>200$ times) than the mean free path of the gas, eliminating the finite domain effect in the DSMC computed spectra.

In total, the numerical experimental setting could be divided into two processes: inference and training. The inference process calculates the density fluctuation spectra using the governing equations with the constitutive relations (4). The constitutive relations contain neural networks $M,K$ defined in (6) and (7) with weights to be determined. The training process determines the weights of neural networks by minimizing the loss function (9). The flow chart Fig. 2 summarizes the entire procedure.

III Results

We compare the density fluctuation spectra calculated by our model with the results of the NS equation and the Grad 13 method. For various Knudsen numbers, spectra $\left<\delta\rho^{2}\right>(\tilde{\omega})$ are shown in Fig. 3 as function of the non-dimensionalized frequency $\tilde{\omega}$ . At a small Knudsen number, all models give consistent spectra. However, at large Knudsen number, our model result matches accurately with the DSMC result, while the shape and amplitude of the NS equation and Grad 13 moments method deviate. Therefore, compared with the NS equation and the Grad 13 method, our model gives the most accurate spectra which are close to the DSMC result in both shape and amplitude.

We test the generalization ability of our model performance by predicting velocity spectra. The generalization ability ensures our model learns the rarefied dynamic physics rather than being forced to reproduce the DSMC density spectra data. As a linear benchmark, we predict the velocity fluctuation spectra of rarefied gas. Our model predicted these spectra in Fig 4(a), which matches with DSMC result much better than the NS equation. Moreover, to demonstrate the robustness of our model, we also plotted the $95\%$ confidence interval in Fig 4(b), estimated using multiple runs on randomly sampled training data. We claim our model has a robust generalization ability for rarefied gas fluctuations based on these benchmarks.

The potential risk of our model overfitting the density fluctuation spectra in training data is negligible. Overfitting refers to the phenomenon that the neural network is too powerful to remember the exact shape of spectra. However, given a Knudsen number, our neural network only models the viscosity scaling law, whose output is a number. Such a number is not enough to record the exact shape of spectra, which is a function of $\omega$ . Hence it is impossible for our neural network to learn the spectra by rote, making its potential risk of overfitting negligible. However, the negligible risk of overfitting does not mean our model generalizes well to all situations.

As we will discuss later, our model does not generalize to boundary regions. We demonstrate this by comparing our model with results on microchannel flow [41, 42] in Fig 4(b). The effective viscosity $\frac{\mu(\mathrm{Kn})}{\mu_{0}}$ of our model has a similar trend compared with microchannel flow results. However, it deviates from microchannel flow even at small Knudsen numbers. The reason is microchannel flow contain boundary region with large flow property gradients. Flow properties in such boundary regions are not governed by (1) and do not admit a Fourier decomposition. Therefore our model could not describe flows in boundary regions.

IV Discussion

Data-driven models modeling physics systems typically learn PDEs consisting of derivatives of various orders [19, 20, 21, 22, 23, 24]. The general form of data-driven models for linear constitutive relation is a linear regression of all derivatives as in (2). We have given it a clear physical explanation by pointing out the equivalence between constitutive relation and scaling law for transport coefficients. Our discussion also reveals that high-order derivatives enable constitutive relations to model more accurate scaling laws. The reason is additional terms of high-order derivatives in (2) contribute additional polynomials terms to (26) (Appendix B). These additional terms make scaling laws in (4) more flexible and hence more accurate. Therefore scaling laws helped us explain how high-order derivatives contribute to constitutive relations.

Scaling laws not only gives a physical explanation but also helps to avoid practical problems in learning the constitutive relation. Instead of regression on derivatives, we suggest directly modeling the scaling functions. It helps to avoid two major problems in regression on derivatives: derivative estimation and variable selection. Derivative estimation encounter stability and accuracy issue for high-order derivatives. Directly modeling the scaling functions avoids this without undermining its flexibility. Variable selection from infinite coefficients in (2) is difficult even if sparsity methods are involved. Modeling scaling functions replaced these coefficients with neural networks. In total, our argument suggests a better formulation for data-driven modeling.

Our model is implicitly connected with the Grad 13 method while avoiding its shortcomings by learning from data. Our model implicitly relies on the Grad 13 one-particle distribution $f(\mathbf{c})$ for gas molecule

f=f_{0}\left(1+\frac{3}{4}\mbox{Kn}({c}_{i}{c}_{j}-\frac{1}{3}\delta_{ij}\mathbf{{c}}^{2}){\sigma}_{ij}+\frac{\mbox{Kn}}{\mbox{Pr}}(\frac{{\mathbf{c}}^{2}}{2}-\frac{5}{2}){c}_{i}{q}_{i}\right)

(10)

where $\mathbf{{c}}$ is the non-dimensional peculiar velocity, and $f_{0}$ is the Maxwell distribution. It is because this distribution is the most probable form that admits arbitrary stress ${\sigma}$ and heat flux ${q}$ as its moments, which is the pre-requisite of modeling stress and heat flux as in (2). However, as we have shown in the result section, our model outperforms the NS equation and the Grad 13 method using only three conservation laws. The reason is our model uses the data-driven approach that learns the rarefied gas dynamics from data points equally, which converges uniformly for various Knudsen numbers. Contrarily, higher-order perturbation or moment expansion benefits little for more accurate spectra at large Knudsen numbers [30], limited by their slow convergence rate at large Knudsen numbers. In conclusion, by learning uniformly from data, the data-driven approach has demonstrated a clear advantage over the traditional Chapman-Enskog expansion and Grad’s moment method in handling rarefied gas dynamics.

The implicit connection with the Grad 13 method also reveals the limitation of our model. Distribution (10) is unsuitable for strong non-Maxwellian dynamics, such as shock waves [43]. Counter-intuitively, a data-driven model requires constraints rather than flexibility to resolve this, because strong non-Maxwellian dynamics tend to have correlated stress and heat flux that (10) can not handle, due to irregularly shaped distribution. Proper constraints on such correlation may be a future direction for data-driven strong non-Maxwellian models.

Similarly, our model relates to the Hilbert-Chapman-Enskog expansion hence bearing the same weakness at boundary regions. Theoretically, we could determine the coefficients in the constitutive relation model (2) by the Hilbert-Chapman-Enskog expansion. Therefore our model could be treated as a reformulation of it with data-driven enhanced convergence. However, the Hilbert-Chapman-Enskog expansion fails in boundary regions where the solutions have large gradients, such as boundary, shock, and initial layers, because the residual of expansion is proportional to gradients [44]. Hence our model also fails in such regions. Extending our model to boundary regions demands solving additional connection problems [45] concerning the gradients of gas flow, which changes the system equation. We expect those gradients to affect the scaling laws of transport coefficients by breaking the homogenous symmetry. Correspondingly, the neural network will no longer be an even function on Knudsen numbers. Flow gradients may further introduce non-local effects into scaling laws. Therefore extra equations describing the such non-local effects of transport coefficients may be required to extend our model to boundary regions.

Finally, we clarify our model’s valid scenario from machine learning’s point of view. Similar to other data-driven approaches, the training data scope limits our model’s viability. Specifically, the limitation is in two aspects: the range of Knudsen numbers and the governing equation. Our model could only capture physics within the Knudsen number range of the training data. In our case, it covers Knudsen numbers from 0 to 0.25. Our model is unreliable outside this Knudsen number range since it extrapolates the data. This limitation on the Knudsen number range does not undermine the utility of our model because we only need to train the model once on the desired Knudsen number range before applying it to other physical scenarios. As for the governing equation, our model works only for physical scenarios governed by the system equation (1). However, it admits different physical scenarios correspond to distinct initial conditions, such as velocity fluctuation in our test data.

In summary of the paper, we argued that the data-driven regression models for constitutive relations are equivalent to length scaling laws of transport coefficients. Our argument not only provides a theoretical justification for data-driven models but also helps to avoid practical problems. We further modeled constitutive relations based on our argumentation. On calculating the Rayleigh scattering spectra, our model significantly outperforms the Chapman-Enskog expansion and Grad moment methods. Our argumentation also reveals the implicit assumption and limitation of data-driven constitutive relations. Further constraints and modifications are necessary for it to accommodate strong non-equilibrium dynamics and boundary layers.

Acknowledgements.

We wish to acknowledge Professor Lei Wu (Southern University of Science and Technology) in offering suggestions on Rayleigh scattering modeling. We thank Ms Yuan Lan (Hong Kong University of Science and Techonology) for comments on the manuscript.

Appendix A The Governing Equation and Its Non-Dimensionalization

In this section, we hydrodynamics governing equations. Then we list the detailed non-dimensionalization for these conservation laws.

The hydrodynamics equations governing the macroscopic dynamics of gases are equations of gas’ statistical quantities, such as the number density $n$ , mass density $\rho=mn$ , the velocities $\mathbf{v}$ , temperature $T$ , stresses $\sigma$ , heat fluxes $q$ , etc. The most fundamental hydrodynamics equations are the conservation laws for mass, momentum, and energy, as described in [46]. In our paper, we consider linearized hydrodynamics equations for one-dimensional gas. These equations describe small fluctuations of statistical quantities around a specific equilibrium state of stationary gas with density $\rho_{0}$ and temperature $T_{0}$ . Specifically, one could obtain the linearized conservation laws via the first-order expansion of gases’ density, velocity, and temperature at the equilibrium. Here we omit the details of the expansion and give the linearized system directly as

\begin{split}\frac{\partial\bar{\rho}}{\partial\bar{t}}+\frac{\partial\bar{v}}{\partial\bar{x}}&=0\\ \frac{\partial\bar{v}}{\partial\bar{t}}+\frac{\partial\bar{T}}{\partial\bar{x}}+\frac{\partial\bar{\rho}}{\partial\bar{x}}&=-\mbox{Kn}\frac{\partial\bar{\sigma}}{\partial\bar{x}}\\ \frac{3}{2}\frac{\partial\bar{T}}{\partial\bar{t}}+\frac{\partial\bar{\rho}}{\partial\bar{x}}&=-\frac{15}{4}\mbox{Kn}\frac{\partial\bar{q}}{\partial\bar{x}},\end{split}

(11)

in which we use quantities with bars to represent the non-dimensionalized quantities.

In the main text, we omitted bars for the simplicity of notations. We also use non-dimensionalized quantities with bar omitted in appendix B and appendix D. However, we use dimensionalized quantities in SI unit in appendix C, E, F to simplify the computation and ensure the consistency with reference.

Now we describe the details of the non-dimensionalization. Suppose we aims to compute properties of the Fourier mode of fluctuations with wavenumber $k_{0}$ . In the non-dimensionalization process, it is natual to set the reference length scale as the wavelength $L=\frac{2\pi}{|k_{0}|}$ of the Fourier mode. Correspondingly, the reference time scale $T=\frac{L}{\sqrt{k_{B}T_{0}/m}}$ is time used by the sound wave to travel the distance of the reference length scale. We denote explicitly the non-dimensionalization of other quantities here. The non-dimensionalized time and spatial position in $x$ direction are $\bar{t}=\frac{t}{T}$ , and $\bar{x}=\frac{x}{L}$ . The non-dimensionalized Fourier wavenumber and angular frequency are $\bar{k}=kL=\frac{2\pi k}{|k_{0}|}$ , $\bar{\omega}=\omega T$ . The non-dimensionalized velocity ( $x$ direction component), density, and temperature are $\bar{v}_{x}=\frac{Tv_{x}}{L}$ , $\bar{\rho}=\frac{\rho-\rho_{0}}{\rho_{0}}$ , $\bar{T}=\frac{T-T_{0}}{T_{0}}$ . Moreover, the non-dimensionalized stress and heat flux in $x$ direction are $\bar{\sigma}_{xx}=\frac{L}{\mu_{0}}\sqrt{\frac{m}{k_{B}T_{0}}}\sigma_{xx}$ , $\bar{q}_{x}=\frac{L}{\kappa_{0}T_{0}}q_{x}$ , in which $\mu_{0}$ and $\kappa_{0}$ are the viscosity and heat conduction coefficients at equilibrium satisfying the relations $\kappa_{0}=\frac{5k_{B}}{2m}\frac{\mu_{0}}{\mbox{Pr}}$ , $\mbox{Pr}=\frac{2}{3}$ .

We define the Knudsen number $\mathrm{Kn}$ in (11) of the Fourier mode of fluctuations with wavenumber $k_{0}$ as

\begin{split}\mathrm{Kn}&=\frac{l}{2\pi/|k_{0}|}\\ l&=\sqrt{\frac{m}{k_{B}T_{0}}}\frac{\mu_{0}}{\rho_{0}},\\ \end{split}

(12)

in which $l$ is the mean free path and $\mu_{0},\kappa_{0}$ are the viscosity and heat conduction coefficient of the gas at equilibrium. For other Fourier modes non-dimensionalized wavenumber $\bar{k}$ , their Knudsen number $\mathrm{Kn}(\bar{k})$ are proportional to $\bar{k}$ as follows

\begin{split}\mathrm{Kn}(\bar{k})&=\frac{l}{2\pi/|k|}=\frac{|\bar{k}|}{2\pi}\mathrm{Kn}\end{split}

(13)

Finally, we give the formula to change the reference length scale. It is helpful to change the wavenumber $k_{0}$ of interest to another wavenumber, which corresponds to changeing the reference length scale $L$ to the wavelength of another Fourier mode. Suppose a non-dimensionalized physical quantity $\bar{f}$ are obtained from its dimensionalized version $f$ by $\bar{f}(\bar{x})=\frac{1}{D}f(\bar{x}L)$ . Define the spatial Fourier transformation of $\bar{f}(\bar{x})$ as $\bar{g}(\bar{k})$ , while the spatial Fourier transformation of $f(x)$ as $g(k)$ . Here we use the Fourier transformation in the symmetrical form [47, Eq 13.5]. Then $\bar{g}(\bar{k})$ and $g(k)$ satisfies

\bar{g}(\bar{k})=\frac{1}{DL}g(\frac{\bar{k}}{L})

(14)

With the help of (14), we could easily handle Fourier transformation from one non-dimensionalization of reference length scale to other reference length scales.

Appendix B The Equivalence Between Data-Driven Constitutive Relations and Scaling Laws

This section shows the equivalence between (2) and (4) in the main text. Before the derivation, we first discuss if the constitutive relation (2) is well-defined via dimension analysis.

As shown in the non-dimensionalization process in Appendix A, our system’s only degree of freedom is the Knudsen number. As a result, any non-dimensional numbers, such as the coefficients $a,b,c,d,e,f$ in the data-driven constitutive relation model (2) are functions depending on the Knudsen number which complex our analysis. For simplicity, in this section, we use the non-dimensionalization with $k_{0}=\frac{2\pi}{l}$ which makes $L=l$ . Under this non-dimensionalization, the Knudsen number of the system is fixed to $\mathrm{Kn}=1$ , while the non-dimensionalized wavenumber $k$ corresponds to the ‘relative’ Knudsen number of each Fourier mode. This non-dimensionalization makes the constitutive relation

\begin{split}{\sigma}&=-\sum_{n=1}^{\infty}\left({a}_{n}\frac{\partial^{n}{v}}{\partial{x}^{n}}+{c}_{n}\frac{\partial^{n}{\rho}}{\partial{x}^{n}}+{e}_{n}\frac{\partial^{n}{T}}{\partial{x}^{n}}\right)\\ {q}&=-\sum_{n=1}^{\infty}\left({b}_{n}\frac{\partial^{n}{T}}{\partial{x}^{n}}+{d}_{n}\frac{\partial^{n}{\rho}}{\partial{x}^{n}}+{f}_{n}\frac{\partial^{n}{v}}{\partial{x}^{n}}\right)\end{split}

(15)

well-defined with coefficients $a,b,c,d,e,f$ independent of the global $\mathrm{Kn}$ . Note that this constitutive relation is exactly (2) in the main text.

We now show that the above data-driven constitutive relation (15) is equivalent to the constitutive relation with scaling laws (4) in the form

\begin{split}\tilde{\sigma}(k)&=-ik\frac{4}{3}\frac{\mu(k)}{\mu_{0}}\tilde{v}(k);\quad\mu(k)\geq 0\\ \tilde{q}(k)&=-ik\frac{\kappa(k)}{\kappa_{0}}\tilde{T}(k);\quad\kappa(k)\geq 0,\end{split}

(16)

The constraint on entropy production is the key to achieving the equivalence. According to [32], The total entropy production rate of the conservation laws of the mass, moment, and energy is

\begin{split}\dot{S}&=\int-\frac{q_{j}}{T^{2}}\partial_{j}T-\frac{\sigma_{ij}}{2T}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)d^{3}x,\end{split}

(17)

, in which we use the Einstein’s summation convention with $i,j\in\{x,y,z\}$ as the dummy indices. As a fundamental constraint, the second law of thermodynamics requires the total entropy change rate $\dot{S}\geq 0$ . However, it is not enough to deduce the equivalence we desired.

Linearization of the system dramatically helps us by simplifying the constraints. First, we could replace the temperature $T$ with the equilibrium temperature $T_{0}$ as a first-order approximation since $\delta T=T-T_{0}$ is a small quantity.

\begin{split}\dot{S}&\approx\int-\frac{q_{j}}{T_{0}^{2}}\partial_{j}T-\frac{\sigma_{ij}}{2T_{0}}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)d^{3}x\geq 0,\end{split}

(18)

Second, the two terms for velocity and temperature field must be non-negative separately

\int q_{j}\partial_{j}Td\mathbb{x}\leq 0;\quad\int\sigma_{ij}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)d\mathbb{x}\leq 0,

(19)

because we know from statistical mechanics that the deviations of velocity and temperature from equilibrium are statistically independent [33]. If we reduce our problem to the 1D case and denote $v_{x}$ as $v$ , $\sigma_{xx}$ as $\sigma$ , and $q_{x}$ as $q$ . A Fourier transformation on (19) with the help of the Paseval’s theorem lead us to

\int\tilde{q}(k)(ik\tilde{T}(k))^{*}dk\leq 0;\quad\int\tilde{\sigma}(k)(ik\tilde{v}(k))^{*}dk\leq 0,

(20)

in which symbol with tilde represents the Fourier transform of corresponding function in $x$ and $*$ indicates complex conjugate. Finally, the fluctuation of statistical quantities generally behaves like white noise that spreads over the entire spectra. Therefore the linear system should be stable for each wavenumber $k$ , which means the entropy production from each wavenumber must be non-negative.

\tilde{q}(k)(ik\tilde{T}(k))^{*}\leq 0;\quad\tilde{\sigma}(k)(ik\tilde{v}(k))^{*}\leq 0,

(21)

Now we discuss what constraint (21) imposes on the constitutive relations. First, stress $\sigma$ depends on the velocity field $v$ only, while the heat flux $q$ depends solely on the temperature field $T$ . This is the only way to ensure (21) since density $\rho$ , velocity $v$ , and temperature $T$ fluctuations are statistically independent [48, 33] with no guarantee in their mutual products. Therefore the linear constitutive relations (15) reduce to the from

\begin{split}{\sigma}&=-\sum_{n=1}^{\infty}{a}_{n}\frac{\partial^{n}{v}_{x}}{\partial{x}^{n}}\\ {q}&=-\sum_{n=1}^{\infty}{b}_{n}\frac{\partial^{n}{T}}{\partial{x}^{n}},\end{split}

(22)

under the condition of non-decreasing entropy. A Fourier transformation of the above forms gives

\begin{split}\tilde{\sigma}(k)&=-\sum_{n=1}^{\infty}a_{n}(ik)^{n}v(k)\\ \tilde{q}(k)&=-\sum_{n=1}^{q}b_{n}(ik)^{n}T(k),\end{split}

(23)

Combine it with the constraint 21, leads us to

\begin{split}\sum_{n=1}^{\infty}(ik)^{n+1}a_{n}\left|\tilde{v}(k)\right|^{2}&\leq 0\\ \sum_{n=1}^{\infty}(ik)^{n+1}b_{n}\left|\tilde{T}(k)\right|^{2}&\leq 0\\ \end{split}

(24)

There should not be any imaginary part appearing in the LHS of (24). Hence all terms with odd powers on $ik$ vanishes. What remains is the following

\begin{split}(k^{2}a_{1}-k^{4}a_{3}+k^{6}a_{5}-k^{8}a_{7}\cdots)\left|\tilde{v}(k)\right|^{2}&\geq 0\\ (k^{2}b_{1}-k^{4}b_{3}+k^{6}b_{5}-k^{8}b_{7}\cdots)\left|\tilde{T}(k)\right|^{2}&\geq 0\\ a_{0}=a_{2}=\cdots=a_{{2n}}=\cdots=&0\\ b_{0}=b_{2}=\cdots=b_{{2n}}=\cdots=&0\\ \end{split}

(25)

Substitute (25) into (23) and take derivative w.r.t $x$ gives us the constitutive relation in Fourier space

\begin{split}\widetilde{\partial_{{x}}{\sigma}}({k})&=({k}^{2}{a}_{1}-{k}^{4}{a}_{3}+{k}^{6}{a}_{5}\cdots)\tilde{v}_{{k}}={k}^{2}M({k})\tilde{v}{(k)}\\ \widetilde{\partial_{{x}}{q}}({k})&=({k}^{2}{b}_{1}-{k}^{4}{b}_{3}+{k}^{6}{b}_{5}\cdots)\tilde{T}_{{k}}={k}^{2}K({k})\tilde{T}{(k)},\\ \end{split}

(26)

in which function $M,K$ are the infinite sum of series and should be non-negative even functions of $k$ . Note that the Knudsen number $\mathrm{Kn}(k)$ corresponds to the Fourier mode with wavenumber $k$ is $Kn(k)=\frac{|k|}{2\pi}$ according to (12). Therefore functions $M,K$ are even functions of Knudsen numbers, hence are well defined in the sense of dimensional analysis.

Functions $M,K$ are closely related to viscosity and heat conduction coefficients. They could be rewritten in the following form

\begin{split}M({k})&=\frac{4}{3}\frac{\mu({k})}{\mu_{0}}\\ K({k})&=\frac{\kappa({k})}{\kappa_{0}}\\ \end{split}

(27)

in which $\mu({k}),\kappa({k})$ are scaling laws of viscosity and heat conduction coefficients satisfying $\mu(0)=\mu_{0}$ and $\kappa(0)=\kappa_{0}$ . Under this notation, $\mu({k})=\mu_{0},\kappa({k})=\kappa_{0}$ exactly corresponds to the constitutive relation for the NS equation.

We could further deduce the stress and heat flux under the same notation with (27) by removing the spatial derivative in (26). The result is exactly (16), which completes our derivation from (15) to (16). Hence we have shown the equivalence between (2) and (4) in the main text.

In addition, we introduce the constitutive relation under the same non-dimensionalization with Appendix A, which is useful in the computation of spectra. Recall that we have used the reference length scale $L=l$ in this section instead of $L=\frac{2\pi}{|k_{0}|}$ in Appendix A. If we change the reference length scale to $L=\frac{2\pi}{|k_{0}|}$ . The constitutive relation (26) will becomes

\begin{split}\widetilde{\partial_{{x}}{\sigma}}({k})&={k}^{2}M({k}\mathrm{Kn})\tilde{v}{(k)}\\ \widetilde{\partial_{{x}}{q}}({k})&={k}^{2}K({k}\mathrm{Kn})\tilde{T}{(k)},\\ \end{split}

(28)

in which $\mathrm{Kn}=\frac{l}{2\pi/|k_{0}|}$ . This result could be derived with the help of (14). The merit of choosing $L=\frac{2\pi}{|k_{0}|}$ as the reference length scale is that the non-dimensionalized $k=2\pi$ exactly corresponds to the dimensionalized wavenumber $k_{0}$ . Therefore we could substitute $k$ with $2\pi$ everywhere if we are only interested in the dynamics at wavenumber $k_{0}$ .

Appendix C Rayleigh Scattering and Density Fluctuations

This section we introduce the Rayleigh scattering spectra and show that it is proportional to the density fluctuation spectra. The Rayleigh scattering was discovered by Lord Rayleigh back in the nineteenth century. It is the reason for the blue color of the sky in daytime and twilight. Specifically, the Rayleigh scattering is due to the refraction of electromagnetic waves (EM waves) passing through media with density fluctuations. Such fluctuation leads to changes in the dielectric constant hence generating refracted EM waves. This section only gives a rough introduction emphasizing the physical picture of Rayleigh scattering and its connection to density fluctuation spectra. One could refer to [37, 38, 39] for a detailed treatment of the Rayleigh scattering spectra. In addition, we do not use non-dimensionalization in this section for simplicity in discussing the related electrodynamics.

The Rayleigh scattering describes the refraction of incident electromagnetic waves passing through gas media with stochastic density fluctuation $\delta\rho$ . Considering the incident electromagnetic wave as plain EM wave with given wavevector $\mathbf{k}_{i}$ ,

\begin{split}\mathbf{E}_{inc}&=\bm{\xi}_{0}\exp(i\mathbf{k}_{i}\cdot\mathbf{r}+i\omega_{i}t)\\ \left|\mathbf{k}_{i}\right|&=\frac{\sqrt{\epsilon_{0}}\omega_{i}}{c}\end{split}

(29)

With $\bm{\xi}_{0}$ be the polarization vector, $\mathbf{k}_{i}$ the incident wave vector, $\omega_{i}$ the incident wave frequency, $c$ is the speed of light in vacuum, and $\epsilon_{0}$ is the dielectric constant of gas. The propagation of the incident wave is governed by the Maxwell’s equations in matter without source [49, Eq 10.21]. We approximate the permeability $\mu$ of gas with the permeability of the vacuum $\mu_{0}$ since they are very close for most materials. Under this approximation, the Maxwell’s equations reduces to

\begin{split}\nabla\cdot\mathbb{D}&=0\\ \nabla\times\nabla\times\mathbb{E}&=-\frac{1}{c^{2}}\frac{\partial^{2}\mathbb{D}}{\partial t^{2}},\\ \end{split}

(30)

in which $\mathbb{D}=\epsilon\mathbb{E}$ , $\epsilon$ is the dielectric constant of gases. This equation governs the propagation of the incident wave in gas media.

Now we analyze how the stochastic density fluctuation $\delta\rho$ affects the propagation of the incident wave $\mathbf{E}_{inc}$ . The dielectric constant of gases $\epsilon$ is a known function of the gas density. Therefore small fluctuations $\delta\rho$ in the gas density leads to the perturbation in $\epsilon$ and $\mathbb{E}$

\begin{split}\epsilon(\mathbb{r},t)&=\epsilon_{0}+\epsilon_{1}(\mathbb{r},t)+\cdots;\ \epsilon_{1}(\mathbb{r},t)=\frac{\partial\epsilon}{\partial\rho}(\rho_{0})\delta\rho(\mathbb{r},t)\\ \mathbb{E}(\mathbb{r},t)&=\mathbb{E}_{inc}(\mathbb{r},t)+\mathbb{E}_{1}(\mathbb{r},t)+\cdots,\\ \end{split}

(31)

in which $\rho_{0}$ is the equilibrium density of gas. Substitute the above expansion in to (30) yields perturbation equation of different orders obtained via perturbation theory in [37]. Specifically, the first order perturbation of the electric field $\mathbb{E}_{1}$ could be calculated by solving the following Helmholtz equation

\nabla^{2}\mathbb{D}_{1}-\frac{\epsilon_{0}}{c^{2}}\frac{\partial^{2}\mathbb{D}_{1}}{\partial t^{2}}=-\nabla\times\nabla\times(\epsilon_{1}\mathbb{E}_{inc})

(32)

in which $\mathbb{D}_{1}=\epsilon_{0}\mathbb{E}_{1}+\epsilon_{1}\mathbb{E}_{0}$ .

Under specific scenarios, we could find analytical solutions to (32) with the help of the Born approximation. Specifically, we consider observing the EM wave at position $\mathbb{r}$ scattered from gases of a certain volume $V$ centered at the origin. In addition, we assume $\epsilon_{1}=\delta\rho=0$ outside the volume $V$ . The distance between the observation position $\mathbb{r}$ and the volume $V$ is so large compared to the radius of the volume $V$ that it allows us to adopt the Born approximation [39, Eq 117.4] [50, Eq 10.53, 10.73] to compute the electric field

\begin{split}\mathbb{E}_{1}(\mathbb{r},\omega_{f})&=\frac{-\hat{\mathbb{r}}\times\hat{\mathbb{r}}\times\bm{\xi}_{0}}{\sqrt{2/\pi}c^{2}}\frac{\omega_{f}^{2}e^{ik_{f}r}}{r}\tilde{\epsilon}_{1}(\mathbb{k}_{f}-\mathbb{k}_{i},\omega_{f}-\omega_{i})\end{split}

(33)

in which $\mathbb{E}_{1}(\mathbb{r},\omega_{f})$ is the Fourier transform of $\mathbb{E}_{1}(\mathbb{r},t)$ w.r.t time, $k_{f}=\sqrt{\epsilon_{0}}\omega_{f}/c$ , $\mathbb{k}_{f}=k_{f}\hat{\mathbb{r}}$ , $\hat{\mathbb{r}}$ is the unit vector along the direction of $\mathbb{r}$ , $r=|\mathbb{r}|$ , $\tilde{\epsilon}_{1}(\mathbb{k},\omega)$ is the Fourier transformation of the function $\epsilon_{1}(\mathbb{r},t)$ w.r.t spatial and temporal coordinate. We call $\mathbb{E}_{1}$ as the electric field of scattered EM wave, whose amplitude is proportional to the perturbation of the dielectric constant.

Next, we define the Rayleigh scattering spectra as the intensity spectra of the scattered electric field $\mathbb{E}_{1}$ and connect it to the density fluctuation spectra. The Rayleigh scattering spectra refers to the intensity spectra $I$ of the scattered EM wave $\mathbb{E}_{1}$

I(\mathbb{r},\omega_{f})=\left<\mathbb{E}_{1}^{2}\right>(\mathbb{r},\omega_{f})=\frac{\sqrt{2\pi}}{T}|\mathbb{E}_{1}(\mathbb{r},\omega_{f})|^{2}\\

(34)

in which we assume $\mathbb{E}_{1}$ is a periodic function in time with period $T$ and $\left<\mathbb{E}_{1}^{2}\right>$ is the spectra of $\mathbb{E}_{1}$ w.r.t time.

We explain the definition of spectra in (34) in this paragraph. The spectra of a real function $f(t)$ with period $T$ is defined as

\begin{split}\left<f^{2}\right>(\omega)&=\left<\mathcal{F}_{t}\left[f(s)f(s+t)\right]\right>_{s}\\ &=\frac{1}{\sqrt{2\pi}T}\int_{-\frac{T}{2}}^{\frac{T}{2}}ds\int_{-\frac{T}{2}}^{\frac{T}{2}}dtf(s)f(s+t)e^{-i\omega t},\\ \end{split}

(35)

in which $\left<\right>_{s}$ represents the ensemble average w.r.t s defined as $\left<f(s)\right>_{s}=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(s)ds$ , $\mathcal{F}_{t}$ represents the Fourier transformation for periodic function w.r.t $t$ defined as $\mathcal{F}_{t}(f)(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-i\omega t}dt$ . The Wiener–Kinchin theorem [51] simplify the definition (35) to

\left<f^{2}\right>(\omega)=\frac{\sqrt{2\pi}}{T}|\tilde{f}(\omega)|^{2},

(36)

where $\tilde{f}(\omega)=\mathcal{F}_{t}(f)(\omega)$ is the Fourier transformation of $f$ . Similarly, the spectra for vector valued function $\mathbb{v}(t)$ with period $T$ are

\begin{split}\left<\mathbb{v}^{2}\right>(\omega)&=\left<\mathcal{F}_{t}\left[\mathbb{v}(s)\cdot\mathbb{v}(s+t)\right]\right>_{s}=\frac{\sqrt{2\pi}}{T}|\tilde{\mathbb{v}}(\omega)|^{2},\end{split}

(37)

which is exactly the definition we used in (34).

We are ready to show that the Rayleigh scattering intensity spectra (34) is proportional to the density fluctuation spectra. The term $|\mathbb{E}_{1}(\mathbb{r},\omega_{f})|^{2}$ in (34) could be calculated from (33) as

\begin{split}|\mathbb{E}_{1}(\mathbb{r},\omega_{f})|^{2}&=\frac{|\bm{\xi}_{0}|^{2}\sin(\psi)^{2}\omega_{f}^{4}}{2c^{4}r^{2}/\pi}|\tilde{\epsilon}_{1}(\mathbb{k}_{f}-\mathbb{k}_{i},\omega_{f}-\omega_{i})|^{2},\end{split}

(38)

where $\psi$ is the angle between $\mathbb{r}$ and $\bm{\xi}_{0}$ . The perturbation of the dielectric constant $\epsilon_{1}$ is proportional to the perturbation of density $\delta\rho$ by its definition in (31), hence we have

\begin{split}|\tilde{\epsilon}_{1}(\mathbb{k},\omega)|^{2}&=\left(\frac{\partial\epsilon}{\partial\rho}\right)^{2}|\delta\tilde{\rho}(\mathbb{k},\omega)|^{2}\end{split}

(39)

in which the $\delta\tilde{\rho}(\mathbb{k},\omega)$ is the Fourier transformation of $\delta\rho(\mathbb{r},t)$ in both the spatial and temporal coordinates.

The equation (39) establish the connection between intensity spectra (34) and the density fluctuation spectra. Assuming $\delta\rho(\mathbb{r},t)$ to be periodic in time with period T and noting that it vanishes outside volume $V$ , the explicit formula for $\delta\tilde{\rho}$ as the Fourier transformation of $\delta\rho$ is

\delta\tilde{\rho}(\mathbb{k},\omega)=\frac{1}{(2\pi)^{2}}\int_{-T/2}^{T/2}dt\int_{V}d\mathbb{r}^{3}\delta\rho(\mathbb{r},t)e^{-i\omega t}e^{-i\mathbb{k}\cdot\mathbb{r}},

(40)

Note that $\delta\tilde{\rho}$ is equal to the Fourier transform of the periodic extension of $\delta\rho$ with unit cell $V$ . Therefore we could treat $\delta\rho$ as if it is a periodic function in space, hence define the density fluctuation spectra by extending (37) to both the spatial and temporal coordinates

\begin{split}\left<\delta\rho^{2}\right>(\omega,\mathbb{k})&=\frac{(2\pi)^{\frac{d+1}{2}}}{TV}|\delta\tilde{\rho}(\mathbb{k},\omega)|^{2}\\ \end{split}

(41)

in which $d$ is the spatial dimension of volume $V$ . Finally, we get the intensity spectra $I$ of the scattered EM wave $\mathbb{E}_{1}$ in terms of the density fluctuation spectra by combining (34) (38) (39) (41)

I(\mathbb{r},\omega_{f})=\frac{V}{\sqrt{2\pi}^{3}}\frac{|\bm{\xi}_{0}|^{2}\sin(\psi)^{2}\omega_{f}^{4}}{2c^{4}r^{2}/\pi}\left(\frac{\partial\epsilon}{\partial\rho}\right)^{2}\left<\delta\rho^{2}\right>(\delta\omega,\delta\mathbb{k}),

(42)

in which $\delta\mathbb{k}=\mathbb{k}_{f}-\mathbb{k}_{i},\delta\omega=\omega_{f}-\omega_{i}$ .

The Rayleigh scattering intensity spectra (42) follow the famous $\omega_{f}^{4}$ frequency dependence, which means blue light (high $\omega_{f}$ ) is scattered more strongly than red light (low $\omega_{f}$ ). This frequency dependence is responsible for the blue color of the sky in daytime and twilight. Moreover, the most intense scattering happens at $\psi=\frac{\pi}{2}$ when one observes the scattering light in the direction perpendicular to the incident wave. Therefore the zenith is bluer than the sky at the horizon in the daytime. However, the detailed shape of the intensity spectra $I$ proportionally depends on the density fluctuation spectra $\left<\delta\rho^{2}\right>$ , which is determined by the hydrodynamics of gases in the scattering region.

Appendix D Calculating the density fluctuation spectra

This section computes the density fluctuation spectra $\left<\delta\rho^{2}\right>(\omega,k)$ introduced in the previous section. For simplicity, We omit the $\delta$ symbol throughout this section, making the notation consistent with Appendix A. Moreover, we non-dimensionalize the density fluctuation spectra as a function of the Knudsen number $\mathrm{Kn}$ and the non-dimensionalized $\omega$ in the form $\left<\rho^{2}\right>(\omega,\mathrm{Kn})$ . The computation largely follows [52].

First, we investigate the symmetries of spectra we shall exploit in computing the density fluctuation spectra. For a real function $f(t)$ with period $T$ , its correlation function

\left<f^{2}\right>(t)=\left<f(s)f(s+t)\right>_{s}

(43)

is an even function satisfying $\left<f^{2}\right>(-t)=\left<f^{2}\right>(t)$ . This could be deduced from the definition of ensemble average (35) in the previous section. It also holds when $T$ tends to infinity if $f$ is a stationary process.

The spectra $\left<f^{2}\right>(\omega)$ defined in (35) is exactly the Fourier transformation of the correlation function $\left<f^{2}\right>(t)$ . From now on we distinguish them by denoting the spectra $\left<f^{2}\right>(\omega)$ as $\left<f^{2}\right>_{\omega}$ while keeping the notation $\left<f^{2}\right>$ for $\left<f^{2}\right>(t)$ .

We could describe the spectra $\left<f^{2}\right>_{\omega}$ in terms of a one-sided Fourier transformation since the correlation function $\left<f^{2}\right>(t)$ is even. Specifically, we define the one-sided Fourier transformation of $\left<f^{2}\right>(t)$ as

\left<f^{2}\right>_{\omega}^{(+)}(\omega)=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}\left<f^{2}\right>(t)e^{-i\omega t}dt.

(44)

It is enough to obtain the full spectra $\left<f^{2}\right>_{\omega}(\omega)$ because the even correlation function $\left<f^{2}\right>(t)$ gives us the property

\left<f^{2}\right>_{\omega}(\omega)=2\mbox{Re}\left[\left<f^{2}\right>_{\omega}^{(+)}(\omega)\right],

(45)

in which Re represents the real part of complex numbers. Note that the one-sided Fourier transformation we use here is just the complex version of the Fourier cosine transformation. The one-sided Fourier transformation also have the property

\left(\partial_{t}g\right)_{\omega}^{(+)}(\omega)=i\omega g_{\omega}^{(+)}(\omega)-\frac{g(0)}{\sqrt{2\pi}}

(46)

if $g$ vanishes at infinity.

Another important property of the correlation is that the correlation $\left<fg\right>=\left<f(s)g(s+t)\right>_{s}$ between periodic function $f$ and $g$ is a linear functional acting on $g$ . It commute with the derivation $\partial$ w.r.t $g$

\left<f\partial g\right>=\partial\left<fg\right>

(47)

Consequently, the stress $\sigma(v)$ and heat flux $q(T)$ as linear functionals of $v$ and $T$ in the form of (22) satisfies

\begin{split}\left<f\sigma(v)\right>&=\sigma(\left<fv\right>)\\ \left<fq(T)\right>&=q(\left<fT\right>)\end{split}

(48)

because $\sigma,q$ are linear combinations of spatial derivatives $\partial_{x}$ that satisfies (47).

With these properties, we are ready to compute the density fluctuation spectra from the linear system (11) and constitutive relation (4). Taking the correlation between density $\rho$ and the governing equation (11), we obtain the linear system

\begin{split}\frac{\partial\left<{\rho}^{2}\right>}{\partial{t}}+\frac{\partial\left<{\rho v}\right>}{\partial{x}}&=0\\ \frac{\partial{\left<{\rho}v\right>}}{\partial{t}}+\frac{\partial\left<{\rho T}\right>}{\partial{x}}+\frac{\partial\left<{\rho^{2}}\right>}{\partial{x}}&=-\mbox{Kn}\frac{\partial\sigma(\left<{\rho v}\right>)}{\partial{x}}\\ \frac{3}{2}\frac{\partial\left<{\rho T}\right>}{\partial{t}}+\frac{\partial\left<{\rho^{2}}\right>}{\partial{x}}&=-\frac{15}{4}\mbox{Kn}\frac{\partial q(\left<{\rho T}\right>)}{\partial{x}},\end{split}

(49)

The above linear system governs the correlations of density with density, velocity, and temperature. However, the initial conditions for the linear system are required to determine the correlations completely.

The initial conditions of (49) describe the two-point correlations of densities, velocities, and temperatures between two simultaneous locations separated by distance $x$ at $t=0$ . Such correlation vanishes if the distance $x$ is non-zero since changes in one place require time to propagate to another. Therefore initial conditions should be delta function $\delta(x)$ multiplied by some amplitude constants. These amplitude constants could be determined by the fluctuation theory in statistical mechanics. The initial condition for $\left<{\rho}^{2}\right>(0,x)$ could be deduced from [48, Eq 88.2] with non-dimensionalization and DSMC’s Monte Carlo effects considered. As for $\left<{\rho v}\right>$ and $\left<{\rho T}\right>$ , they vanishes at $t=0$ since fluctuations of density $\rho$ , velocity $v$ , and temperature $T$ are statistically independent [48, 33]. Therefore we have

\begin{split}\left<{\rho}^{2}\right>(0,x)&=\frac{mN_{eff}}{\rho_{0}L}\delta(x)\\ \left<{\rho v}\right>(0,x)&=0\\ \left<{\rho T}\right>(0,x)&=0,\\ \end{split}

(50)

in which $m$ is the mass of gas molecule, $N_{eff}$ is the effective number of molecules per particle used in the DSMC simulation taking its Monte Carlo fluctuation into account, $\rho_{0}$ is the equilibrium gas density, $L$ is the reference length scale used in the non-dimensionalization.

To solve the linear system (49), we take the Fourier transform on the spacial coordinate and the one sided Fourier transformation on the temporal coordinate. With the help of the constitutive relation (28) we obtain

\begin{split}i\omega\left<{\rho}^{2}\right>_{\omega,k}^{+}+ik\left<{\rho}{v}\right>_{\omega,k}^{+}&=\frac{mN_{eff}}{2\pi\rho_{0}L}\\ i\omega\left<{\rho}{v}\right>_{\omega,k}^{+}+ik\left<{\rho}{T}\right>_{\omega,k}^{+}+ik\left<{\rho}^{2}\right>_{\omega,k}^{+}&=-k^{2}A_{k}(\mathrm{Kn})\left<{\rho}{v}\right>_{\omega,k}^{+}\\ \frac{3}{2}i\omega\left<{\rho}{T}\right>_{\omega,k}^{+}+ik\left<{\rho}{v}\right>_{\omega,k}^{+}&=-k^{2}B_{k}(\mathrm{Kn})\left<{\rho}{T}\right>_{\omega,k}^{+}\\ \end{split}

(51)

in which we define $A_{k}(\mathrm{Kn})=\mbox{Kn}M(k\mathrm{Kn})$ , $B_{k}(\mathrm{Kn})=\frac{15}{4}\mbox{Kn}K(k\mathrm{Kn})$ , and $f^{+}_{\omega,k}$ for arbitrary function $f(t,x)$ is obtained by taking the one sided Fourier transformation on time coordinate $t$ and the Fourier transform on space coordinate $x$ .

Finally, the solution of $\left<\delta{\rho}^{2}\right>_{\omega,k}^{+}$ is a function of the Knudsen number $\mathrm{Kn}$ and the frequency $\omega$ only

\left<{\rho}^{2}\right>_{\omega,k=2\pi}^{+}(\mathrm{Kn})=-\frac{ikN_{eff}m\left(-2k^{4}A_{k}(\text{Kn})B_{k}(\text{Kn})-3ik^{2}\omega A_{k}(\text{Kn})-2ik^{2}\omega B_{k}(\text{Kn})-2k^{2}+3\omega^{2}\right)}{(2\pi)^{2}\rho_{0}\left(-2k^{4}\omega A_{k}(\text{Kn})B_{k}(\text{Kn})-3ik^{2}\omega^{2}A_{k}(\text{Kn})+2ik^{4}B_{k}(\text{Kn})-2ik^{2}\omega^{2}B_{k}(\text{Kn})-5k^{2}\omega+3\omega^{3}\right)},

(52)

it does not depend on $k$ because that the non-dimensionalized wavenumber $k$ will be fixed to $2\pi$ by choosing the reference length scale $L=\frac{2\pi}{|k_{0}|}$ if we consider the Fourier mode at wavenumber $k_{0}$ . The density fluctuation spectra $\left<\rho^{2}\right>(\omega,\mathrm{Kn})$ is two times of the real part of (52) according to (45). This concludes the derivation of density fluctuation spectra for the constitutive relation (28).

Next we determine the density fluctuation spectra for the NS equation and the Grad 13 moment method. The density fluctuation spectra computed from the NS equation could be obtained by simply replace $A_{k},B_{k}$ in (52) by

A_{k}(\mathrm{Kn})=\frac{4}{3}\mathrm{Kn};\quad B_{k}(\mathrm{Kn})=\frac{15}{4}\mathrm{Kn}

(53)

As for the Grad 13 method, $\sigma$ and $q$ are unknown quantities determined by two addition equation for the stress and heat flux [30, Eq 35] as

\begin{split}\mathrm{Kn}\partial_{{t}}{\sigma}+\frac{4}{3}\partial_{{x}}{v}_{x}+\frac{8}{15}\partial_{{x}}{q}&=-\sigma\\ \mathrm{Kn}\partial_{{t}}{q}+\frac{4}{15}\partial_{{x}}{\sigma}+\frac{2}{3}\partial_{{x}}{T}&=-\frac{2}{3}q\end{split}

(54)

Note that the non-dimensinoalization used in [30] differ with us for stress and heat flux, specifically we have $\sigma_{[26]}=\sigma\mathrm{Kn}$ , $q_{[26]}=\frac{15}{4}q\mathrm{Kn}$ .

We apply the one sided Fourier transformation on time and Fourier transformation on space coordinates to the linear system again. The resulting governing equation for correlations of density between density, velocity, temperature, stress and heat fluxes for the Grad 13 are as follows

\begin{split}i\omega\left<{\rho}^{2}\right>_{\omega,k}^{+}+ik\left<{\rho}{v}_{x}\right>_{\omega,k}^{+}&=\frac{mN_{eff}}{2\pi\rho_{0}L}\\ i\omega\left<{\rho}{v}_{x}\right>_{\omega,k}^{+}+ik\left<{\rho}{T}\right>_{\omega,k}^{+}+ik\left<{\rho}^{2}\right>_{\omega,k}^{+}&=-ik\mathrm{Kn}\left<{\rho}{\sigma}\right>_{\omega,k}^{+}\\ \frac{3}{2}i\omega\left<{\rho}{T}\right>_{\omega,k}^{+}+ik\left<{\rho}{v}_{x}\right>_{\omega,k}^{+}+ik\left<{\rho}{q}\right>_{\omega,k}^{+}&=0\\ i\omega\mathrm{Kn}\left<{\rho}{\sigma}\right>_{\omega,k}^{+}+\frac{4}{3}ik\left<{\rho}{v}_{x}\right>_{\omega,k}^{+}+\frac{8}{15}ik\left<{\rho}{q}\right>_{\omega,k}^{+}&=-\left<{\rho}{\sigma}\right>_{\omega,k}^{+}\\ i\omega\mathrm{Kn}\left<{\rho}{q}\right>_{\omega,k}^{+}+\frac{4}{15}ik\left<{\rho}{\sigma}\right>_{\omega,k}^{+}+\frac{2}{3}ik\left<{\rho}{T}\right>_{\omega,k}^{+}&=-\frac{2}{3}\left<{\rho}{q}\right>_{\omega,k}^{+}\\ \end{split}

(55)

The spectra is hence calculated in the same way as NS equation and linear constitutional relation model. The result for the Grad 13 method is

\left<\delta{\rho}^{2}\right>_{\omega,k=2\pi}^{+}(\mathrm{Kn})=\frac{mN_{\text{eff}}k\left(-36ik^{4}\text{Kn}^{2}+189ik^{2}\text{Kn}^{2}\omega^{2}+165k^{2}\text{Kn}\omega-20ik^{2}-45i\text{Kn}^{2}\omega^{4}-75\text{Kn}\omega^{3}+30i\omega^{2}\right)}{(2\pi)^{2}\rho_{0}\left(135k^{4}\text{Kn}^{2}\omega-75ik^{4}\text{Kn}-234k^{2}\text{Kn}^{2}\omega^{3}+240ik^{2}\text{Kn}\omega^{2}+50k^{2}\omega+45\text{Kn}^{2}\omega^{5}-75i\text{Kn}\omega^{4}-30\omega^{3}\right)}

(56)

Finally, we compute the velocity fluctuation spectra as the test data. Taking the correlation between velocity $v$ and the governing equation (11) gives

\begin{split}\frac{\partial\left<{v\rho}\right>}{\partial{t}}+\frac{\partial\left<{v^{2}}\right>}{\partial{x}}&=0\\ \frac{\partial{\left<v^{2}\right>}}{\partial{t}}+\frac{\partial\left<{vT}\right>}{\partial{x}}+\frac{\partial\left<{v\rho}\right>}{\partial{x}}&=-\mbox{Kn}\frac{\partial\sigma(\left<{v^{2}}\right>)}{\partial{x}}\\ \frac{3}{2}\frac{\partial\left<{vT}\right>}{\partial{t}}+\frac{\partial\left<{v\rho}\right>}{\partial{x}}&=-\frac{15}{4}\mbox{Kn}\frac{\partial q(\left<{vT}\right>)}{\partial{x}},\end{split}

(57)

The initial condition for this linear system is obtained following a similar argument with the density spectra case. At time $t=0$ , the only non-vanishing correlation function is $\left<{v^{2}}\right>$ , which is proportional to $\delta(x)$ . The initial condition for $\left<{v}^{2}\right>(0,x)$ we used here is deduced from [48, Eq 88.5] with non-dimensionalization and DSMC’s Monte Carlo effects considered.

\begin{split}\left<{v\rho}\right>(0,x)&=0\\ \left<{v^{2}}\right>(0,x)&=\frac{mN_{eff}}{\rho_{0}L}\delta(x)\\ \left<{vT}\right>(0,x)&=0,\\ \end{split}

(58)

\left<{v}^{2}\right>_{\omega,k=2\pi}^{+}(\mathrm{Kn})=-\frac{iN_{eff}m\omega k\left(3\omega-2ik^{2}B_{k}(\text{Kn})\right)}{(2\pi)^{2}\rho_{0}(-2k^{4}\omega A_{k}(\text{Kn})B_{k}(\text{Kn})-3ik^{2}\omega^{2}A_{k}(\text{Kn})+2ik^{4}B_{k}(\text{Kn})-2ik^{2}\omega^{2}B_{k}(\text{Kn})-5k^{2}\omega+3\omega^{3})}

(59)

The velocity fluctuation spectra $\left<v^{2}\right>(\omega,\mathrm{Kn})$ is two times of the real part of (59) according to (45).

Appendix E DSMC Calculation Details

We use the DSMC0F program by A.Bird [53] to simulate the fluctuation of 1D homogeneous gas. Its geometry is a one-dimensional gas of unit cross section between two plane specularly reflecting walls that are normal to the x-axis. The computation domain between these two plane has spatial span $4.8m$ in the x direction and is divided uniformly into $1281$ cells. Each cell contains $8$ subcells utilized in determining collision pairs in the DSMC computation.

The initial condition of our DSMC computation uses particle velocity sampled as in [54] from the Maxwell distribution at $T=300K$ and zero mean velocity. The particle position is uniformly distributed in each cell. More details about the properties of the gas are shown in Table 1 using SI units.

The merit of the DSMC calculation is that no driven physical conditions are required for simulating fluctuations. It is because the DSMC method uses Monte Carlo samples to mimic the real gas molecules. Hence statistical quantities computed from such Monte Carlo samples naturally fluctuate in the same way as the real gas except for an enlarged fluctuation amplitude. Specifically, if one sample in the DSMC simulation represents $N_{eff}$ real gas molecules, the variances of fluctuations in statistical quantities computed from the DSMC simulation are $N_{eff}$ times larger than those of real gass. In our DSMC computation we have $128100$ simulation particle samples representing gases of number density $10^{20}m^{-3}$ , therefore in our computation each sample particles representing $N_{eff}=3.75\times 10^{15}$ real gas molecules.

No global Knudsen number is defined for our DSMC simulation of homogeneous gas since there is no mean flow across the simulation domain. However, fluctuations in density, velocity, and temperature exist and propagate according to hydrodynamics with well-defined Knudsen numbers. Specifically, the Knudsen numbers are defined for various Fourier modes (Phonons) of the fluctuations according to their wavelength.

The molecular model is crutial in DSMC calculations. It describes how two molecule collide with each other and determines the viscosity of the gas. The molecular model gives the relation between two characteristic quantities of classical binary collision problem [55, 56, 57]: the impact parameter $b$ and scattering angle $\theta$ . One of the typical molecular model used in DSMC is the variable hard/soft sphere model [58]

\theta=2\cos^{-1}((\frac{b}{d})^{\frac{1}{\alpha}})

(60)

in which $d$ is the effective diameter of the gas molecules and $\alpha$ is a parameter mainly effecting the diffusion coefficient. The diffusion parameter describes mass diffusion between components of gas mixtures and is irrelevant in our single species case. Therefore we use the default value $\alpha=1$ in the DSMC0F program corresponding to the variable hard sphere model. The effective diameter $d$ varies with the relative velocity between colliding molecules according to [59, Eq 4.63]

d=d_{ref}(\frac{(2k_{B}T_{ref}/(\frac{1}{2}mv_{r}^{2}))^{\omega-1/2}}{\Gamma(5/2-\omega)})^{1/2}

(61)

in which $m$ is the mass of gas molecule, $v_{r}$ is the relative velocity between the two colliding molecules, $\Gamma$ represents the Gamma function, $T_{ref}$ is the reference temperature, $d_{ref}$ is the reference molecule diameter, and $\omega$ is a parameter determines how viscosity coefficient changes w.r.t temperature. In our computation we use the default value $m=5\times 10^{-26}kg$ , $T_{ref}=273K$ , and $d_{ref}=3.5\times 10^{-10}m$ in the DSMC0F program. Note that (61) differ from the equation in [59] since they use the reduced mass $m_{r}=\frac{1}{2}m$ instead of molecule mass $m$ in our case.

The parameter $\omega$ in (61) determines the power law between viscosity coefficient $\mu$ and the temperature $T$ in the form $\mu\propto T^{\omega}$ [60, Eq 3.66]. However, viscosity coefficients appears only in the stress as a production with the velocity gradients. Such a change in viscosity is of second order $\delta v\delta T$ in the perturbation expansion hence is not important in our first-order linear case. As a result, we again use the default value $\omega=0.5$ in the DSMC0F program.

We compute the viscosity coefficient and heat conduction coefficient of our DSMC simulated gas using the Chapman-Enskog theory [60, Eq 3.73]

\begin{split}\mu_{0}&=\frac{5(\alpha+1)(\alpha+2)(\pi mk_{B})^{1/2}(4k_{B}/m)^{\omega-1/2}T^{\omega}}{16\alpha\Gamma(\frac{9}{2}-\omega)\sigma_{T,ref}v_{r,ref}^{2\omega-1}}\\ \kappa_{0}&=\frac{15k_{B}}{4m}\mu_{0}\end{split}

(62)

in which the reference total cross section $\sigma_{T,ref}=\pi d_{ref}^{2}$ and the reference velocity $v_{r,ref}=\sqrt{\frac{4k_{B}T_{ref}}{m\Gamma(5/2-\omega)^{\frac{1}{\omega-1/2}}}}$ .

To ensure the resolution at relative large Knudsen numbers, our DSMC computation, use the cell width to be 5 times smaller than the mean free path of the gas, while the time step to be 10 times smaller than the mean free time of the gas. We compute mean free path and mean free time from the collision rate per gas molecule according to [61, Eq 4.64]

f=4n\pi^{1/2}d_{ref}^{2}(\frac{T}{T_{ref}})^{1/2-\omega}(\frac{k_{B}T}{m})^{1/2}

(63)

in which $n$ is the number density. The mean free time of gas molecules is $t_{m}=\frac{1}{f}$ , while the mean free path of gas molecules is $l_{m}=\sqrt{\frac{8k_{B}T}{\pi m}}t_{m}\approx 1.28l$ , in which $l$ is the mean free path used in non-dimensionalization in Appendix A.

There is no need to worry about the convergence in the mean flow since our computation simulates homogeneous gas using homogeneous initial conditions. However, the finite simulation domain in our DSMC computation may introduce deviation in the spectra from theoretical results in Appendix D. To eliminate this finite domain effect, we use a domain length much (300 times) larger than the mean free path of the gas. Moreover, random perturbations that possibly appear in the initial condition are fully relaxed since we use the total simulation time of $30$ times more than the transverse time of sound speed over the computation domain.

A snap shot will be stored for every 5 time steps. Then the macroscopic quantities for each cell are calculated by averaging corresponding quantites of particles in each cell. The density fluctuation spectra used to train the neural network is ensemble average from 27 independent DSMC run, computed via discrete Fourier transformation using the equation (41).

Table 1: The coefficients and configuration used in DSMC0F program in SI unit.

Domain Length	4.8	Collision Model	VHS¹¹1Variable hard sphere model
Power law²²2The viscosity-temperature power law used in variable hard sphere model	0.5	Diameter³³3The reference molecule diameter	$3.5\times 10^{-10}$
Num of Cell	1281	Mean Free Path	$1.8\times 10^{-2}$
Simulation Particle	128100	Mean Free Time	$4\times 10^{-5}$
Density	$5\times 10^{-6}$	Temperature	300
Molecule Mass	$5\times 10^{-26}$	Sound Speed	371.56
Heat Conduction ⁵⁵5The heat conduction coefficient	0.0214	Viscosity ⁶⁶6The viscosity coefficient	$2.07\times 10^{-5}$
Time step size	$4\times 10^{-6}$	Cell width	$3.7\times 10^{-3}$
Subcell ⁴⁴4The number of subcell per cell used in particle collision process	8	Simulation time	$4\times 10^{-1}$

Appendix F Neural Network Training Details

In this section, we adopt the dimensionalized quantities instead of non-dimensionalized version in Appendix A to make this section consistent with the DSMC computation which is computed in SI unit. In the dimensionalized notation, the density fluctuation spectra is of the form $\left<\delta\rho^{2}\right>(\omega,k)$ , in which $\omega$ is the frequency and the wavenumber $k$ . The wavenumber $k$ directly determines the Knudsen number of phonons (Fourier modes of $\rho$ ). Given $k$ , the spectra $\left<\rho^{2}\right>_{k}(\omega)$ is a function of the angular frequency $\omega$ .

The training data set consists of 40000 $(k,\omega,\left<\rho^{2}\right>)$ tuples draw. While the validation set consists 400 such tuples. We generate these tuples by draw $k$ uniformly from the interval $\left[0,\frac{\pi\rho_{0}}{2\mu}\sqrt{\frac{k_{B}T_{0}}{m}}\right]$ (corresponds to $\mathrm{Kn}\in[0,0.25]$ ). Then for a given $k$ , we sample $\omega$ from the range $[-3ck,3ck]$ ( $c$ is the speed of sound) with probability proportional to the value of $\left<\rho^{2}\right>_{k}(\omega)$ calculated using DSMC. The merit of such sampling strategy is it emphasis the peak region of the spectra.

The common practice of choosing test set is to sample tuples of $(k,\omega,\left<\rho^{2}\right>)$ from the same distribution with the training set. However, such test set only test how good the neural network fit the density fluctuation spectra, not its ability to generalize to other physics scenarios. Instead, we use the model trained on density fluctuation spectra to predict the velocity fluctuation spectra $\left<v^{2}\right>(\omega,k)$ to test its generalization ability.

The function $M$ is modeled as a fully connected neural network without bias, as shown in the paper. The weights to be trained are $\mathbf{W}_{1,2}$ . These weights are initialized by Pytorch’s default uniform initialization.

The loss is defined as the mean square difference between spectra $\left<\rho^{2}\right>_{DSMC}$ computed using DSMC and the spectra $\left<\rho^{2}\right>$ predicted by linear constitutiont relation model. The optimizer we use is the Adam optimizer with learning rate $\alpha=0.005$ , Beta parameter $\beta_{1}=0.9$ and $\beta_{2}=0.999$ , and the parameter $\epsilon=10^{-8}$ . For each training epoch, the batch size for each step is 64. The training process stops if the loss obtained on validation set increases.

References

Chen and Doolen [1998] S. Chen and G. D. Doolen, Lattice boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30, 329 (1998).
Peter and Kremer [2009] C. Peter and K. Kremer, Multiscale simulation of soft matter systems – from the atomistic to the coarse-grained level and back, Soft Matter 5, 4357 (2009).
Lin and Truhlar [2006] H. Lin and D. G. Truhlar, Qm/mm: what have we learned, where are we, and where do we go from here?, Theoretical Chemistry Accounts 117, 185 (2006).
Kogan [1969] M. N. Kogan, Introduction, in Rarefied Gas Dynamics (Springer US, Boston, MA, 1969) pp. 1–27.
Bird [1994a] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994).
Agarwal et al. [2001] R. K. Agarwal, K.-Y. Yun, and R. Balakrishnan, Beyond navier–stokes: Burnett equations for flows in the continuum–transition regime, Physics of Fluids 13, 3061 (2001).
García-Colín et al. [2008] L. S. García-Colín, R. M. Velasco, and F. J. Uribe, Beyond the navier-stokes equations: Burnett hydrodynamics, Physics Reports 465, 149 (2008).
Hilbert [1912] D. Hilbert, Begründung der kinetischen Gastheorie, Mathematische Annalen 10.1007/BF01456676 (1912).
Cha [1916] VI. On the law of distribution of molecular velocities, and on the theory of viscosity and thermal conduction, in a non-uniform simple monatomic gas, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 10.1098/rsta.1916.0006 (1916).
Enskog [1917] D. Enskog, Kinetische theorie der vorgänge in mässig verdünnten gasen. i. allgemeiner teil, Uppsala: Almquist & Wiksells Boktryckeri (1917).
Bobylev [1982] A. Bobylev, The chapman-enskog and grad methods for solving the boltzmann equation, in Akademiia Nauk SSSR Doklady, Vol. 262 (1982) pp. 71–75.
McLennan [1965] J. A. McLennan, Convergence of the chapman‐enskog expansion for the linearized boltzmann equation, Physics of Fluids 8, 1580 (1965).
Grad [1949] H. Grad, On the kinetic theory of rarefied gases, Communications on Pure and Applied Mathematics 2, 331 (1949).
Struchtrup and Taheri [2011] H. Struchtrup and P. Taheri, Macroscopic transport models for rarefied gas flows: A brief review, IMA Journal of Applied Mathematics (Institute of Mathematics and Its Applications) 76, 672 (2011).
Weiss [1995] W. Weiss, Continuous shock structure in extended thermodynamics, Physical Review E 52, R5760 (1995).
Torrilhon [2009] M. Torrilhon, Hyperbolic moment equations in kinetic gas theory based on multi-variate pearson-iv-distributions, Communications in Computational Physics 7, 639 (2009).
Hana et al. [2019] J. Hana, C. Ma, Z. Ma, and E. Weinan, Uniformly accurate machine learning-based hydrodynamic models for kinetic equations, Proceedings of the National Academy of Sciences of the United States of America 10.1073/pnas.1909854116 (2019).
Zhang and Ma [2020] J. Zhang and W. Ma, Data-driven discovery of governing equations for fluid dynamics based on molecular simulation, Journal of Fluid Mechanics 10.1017/jfm.2020.184 (2020).
Raissi et al. [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Machine learning of linear differential equations using Gaussian processes, Journal of Computational Physics 10.1016/j.jcp.2017.07.050 (2017), arXiv:1701.02440 .
Rudy et al. [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 10.1126/sciadv.1602614 (2017), arXiv:1609.06401 .
Schaeffer [2017] H. Schaeffer, Learning partial differential equations via data discovery and sparse optimization, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 10.1098/rspa.2016.0446 (2017).
Raissi et al. [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 10.1016/j.jcp.2018.10.045 (2019).
Long et al. [2018] Z. Long, Y. Lu, X. Ma, and B. Dong, PDE-Net: Learning PDEs from data, in 6th International Conference on Learning Representations, ICLR 2018 - Workshop Track Proceedings (2018).
Long et al. [2019] Z. Long, Y. Lu, and B. Dong, PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network, Journal of Computational Physics 399, 10.1016/j.jcp.2019.108925 (2019), arXiv:1812.04426 .
LeCun et al. [1995] Y. LeCun, Y. Bengio, et al., Convolutional networks for images, speech, and time series, The handbook of brain theory and neural networks 3361, 1995 (1995).
Dal Santo et al. [2020] N. Dal Santo, S. Deparis, and L. Pegolotti, Data driven approximation of parametrized PDEs by reduced basis and neural networks, Journal of Computational Physics 416, 109550 (2020), arXiv:1904.01514 .
Baydin et al. [2018] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of machine learning research 18 (2018).
Yip and Nelkin [1964] S. Yip and M. Nelkin, Application of a kinetic model to time-dependent density correlations in fluids, Physical Review 135, 10.1103/PhysRev.135.A1241 (1964).
Fiocco and DeWolf [1968] G. Fiocco and J. DeWolf, Frequency spectrum of laser echoes from atmospheric constituents and determination of the aerosol content of air, Journal of Atmospheric Sciences 25, 488 (1968).
Wu and Gu [2020] L. Wu and X.-J. Gu, On the accuracy of macroscopic equations for linearized rarefied gas flows, Advances in Aerodynamics 2, 10.1186/s42774-019-0025-4 (2020).
Chapman and Cowling [1990] S. Chapman and T. G. Cowling, The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases (Cambridge university press, 1990).
Lifshitz and Pitaevskii [2013a] E. M. Lifshitz and L. P. Pitaevskii, Statistical physics: theory of the condensed state, Vol. 9 (Elsevier, 2013) Chap. 88, p. 372.
Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz, Course of theoretical physics, Vol. 5 (Elsevier, 2013) Chap. 112, p. 343.
Michalis et al. [2010] V. Michalis, A. Kalarakis, E. Skouras, and V. Burganos, Rarefaction effects on gas viscosity in the knudsen transition regime, Microfluidics and Nanofluidics 9, 847 (2010).
Ghaem-Maghami and May [1980] V. Ghaem-Maghami and A. D. May, Rayleigh-brillouin spectrum of compressed he, ne, and ar. ii. the hydrodynamic region, Physical Review A 22, 698 (1980).
Goodfellow et al. [2015] I. J. Goodfellow, Y. Bengio, and A. C. Courville, Deep learning, Nature 521, 436 (2015).
Pecora [1964] R. Pecora, Doppler shifts in light scattering from pure liquids and polymer solutions, The Journal of Chemical Physics 40, 1604 (1964).
Jackson [1998a] J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, 1998) Chap. 10.
Landau et al. [2013] L. D. Landau, J. Bell, M. Kearsley, L. Pitaevskii, E. Lifshitz, and J. Sykes, Electrodynamics of continuous media, Vol. 8 (elsevier, 2013) Chap. 117.
Kingma and Ba [2015] D. P. Kingma and J. L. Ba, Adam: A method for stochastic optimization, in 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings (2015) arXiv:1412.6980 .
Roohi and Darbandi [2009] E. Roohi and M. Darbandi, Extending the navier–stokes solutions to transition regime in two-dimensional micro- and nanochannel flows using information preservation scheme, Physics of Fluids 21, 082001 (2009).
Karniadakis et al. [2001] G. E. Karniadakis, A. Beskok, N. R. Aluru, and C.-M. Ho, Microflows and nanoflows: Fundamentals and simulation (2001).
Grad [1952] H. Grad, The profile of a steady plane shock wave, Communications on Pure and Applied Mathematics 5, 257 (1952).
Cercignani [1988a] C. Cercignani, Small and large mean free paths, in The Boltzmann Equation and Its Applications (Springer New York, New York, NY, 1988) p. 238.
Cercignani [1988b] C. Cercignani, Small and large mean free paths, in The Boltzmann Equation and Its Applications (Springer New York, New York, NY, 1988) p. 248.
Kardar [2007] M. Kardar, Statistical physics of particles (2007) Chap. 3, pp. 78–81.
Riley et al. [2006a] K. F. Riley, M. P. Hobson, and S. J. Bence, Integral transforms, in Mathematical Methods for Physics and Engineering: A Comprehensive Guide (Cambridge University Press, 2006) Chap. 13, p. 435, 3rd ed.
Lifshitz and Pitaevskii [2013b] E. M. Lifshitz and L. P. Pitaevskii, Statistical physics: theory of the condensed state, Vol. 9 (Elsevier, 2013) Chap. 88, p. 370.
Jackson [1998b] J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, 1998) Chap. 10.2, p. 463.
Griffiths and Schroeter [2018a] D. J. Griffiths and D. F. Schroeter, Introduction to Quantum Mechanics, 3rd ed. (Cambridge University Press, 2018) Chap. 10.4.
Riley et al. [2006b] K. F. Riley, M. P. Hobson, and S. J. Bence, Integral transforms, in Mathematical Methods for Physics and Engineering: A Comprehensive Guide (Cambridge University Press, 2006) Chap. 13, p. 450, 3rd ed.
Lifshitz and Pitaevskii [2013c] E. M. Lifshitz and L. P. Pitaevskii, Statistical physics: theory of the condensed state, Vol. 9 (Elsevier, 2013) Chap. 89, pp. 373–377.
Bird [1994b] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 11.6.
Bird [1994c] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. Appendix C, p. 426.
Bird [1994d] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 2.2, p. 33.
Griffiths and Schroeter [2018b] D. J. Griffiths and D. F. Schroeter, Introduction to Quantum Mechanics, 3rd ed. (Cambridge University Press, 2018) Chap. 10.1.1, p. 376.
Landau and Lifshitz [1982] L. Landau and E. Lifshitz, Mechanics: Volume 1, v. 1 (Elsevier Science, 1982) Chap. 18, p. 48.
Bird [1994e] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 2.6-2.7.
Bird [1994f] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 4.3, p. 93.
Bird [1994g] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 3.5, pp. 68–70.
Bird [1994h] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, 1994) Chap. 4.3, p. 93.