A fluid in a nonequilibrium state exhibits long-ranged correlations of its hydrodynamic fluctuations. In this article, we examine the effect of a transpiration interface on these correlations—specifically, we consider a dilute gas in a domain bisected by the interface. The system is held in a nonequilibrium steady state by using isothermal walls to impose a temperature gradient. The gas is simulated using both direct simulation Monte Carlo (DSMC) and fluctuating hydrodynamics (FHD). For the FHD simulations, two models are developed for the interface based on master equation and Langevin approaches. For appropriate simulation parameters, good agreement is observed between DSMC and FHD results with the latter showing a significant advantage in computational speed. For each approach, we quantify the effects of transpiration on long-ranged correlations in the hydrodynamic variables. The principal effect of transpiration is a suppression of the correlations, an outcome largely explained by a reduction in the temperature gradient due to the interface. We also observe a distortion of the temperature correlations, specifically the appearance of a new peak located near the interface.

## I. INTRODUCTION

Recent research in the field of fluid dynamics has seen an increased focus on nanoscale and microscale applications and, in particular, the role of thermal fluctuations. These fluctuations are important in areas ranging from microfluidic engineering^{1} to molecular biology.^{2} Fluctuations also yield information about the transport properties of the fluid itself^{3} and are important in dynamic processes such as phase transitions, combustion and ignition, nucleation, hydrodynamic instabilities, and a host of Brownian dynamics phenomena.^{4–11}

Traditional computational fluid dynamic approaches have focused on solutions to the deterministic Navier-Stokes equations, which neglect thermal fluctuations. To capture this important feature of nanoscale flows, particle based methods are typically employed. These techniques operate by replicating, with varying degrees of approximation, the dynamics of the fluid’s constituent molecules. Molecular dynamics^{12} takes the approach of deterministically simulating the interactions between molecules, often using a relatively complex interaction potential. Although this approach accurately simulates fluids at the nanoscale, its computational cost makes it inappropriate for mesoscale applications where thermal fluctuations are of importance. The direct simulation Monte Carlo^{13,14} (DSMC) method greatly simplifies the particle interactions using a stochastic approach, allowing the simulation of much larger systems while still reproducing the correct thermal fluctuations. DSMC solves the Boltzmann equation,^{15} so it can be used for flows of arbitrary Knudsen numbers,^{16} whereas solutions of the Navier-Stokes equations are valid only when the Knudsen number is sufficiently small. However, this technique applies only to dilute gasses and is still relatively computationally expensive.

An alternative method for simulating mesoscopic systems is given by fluctuating hydrodynamics (FHD). Following the approach of Landau and Lifshitz,^{17,18} FHD extends the deterministic Navier-Stokes equations by incorporating stochastic fluxes that lead to fluctuations that are consistent with statistical mechanics. The validity of fluctuating hydrodynamics for nonequilibrium systems has been confirmed experimentally by light scattering^{19} and shadowgraph experiments.^{20} Numerical methods for solving the compressible and incompressible FHD equations are discussed in Refs. 21–23, and the approach has been further extended to multispecies mixtures,^{7,24–29} multiphase flows,^{30} reacting flows,^{31} and electrokinetic flows.^{32,33} Hybrid methods combining FHD and DSMC are investigated in Refs. 34 and 35.

When a fluid is held in a nonequilibrium state, long range correlations in the thermal fluctuations^{18} can result in macroscopic effects.^{27} These nonequilibrium fluctuations have been studied in many scenarios in which the correlations extend through the bulk of the fluid; however, less is known regarding the effect of an interface, such as a porous wall or a thin membrane. This type of interface has previously been studied using molecular dynamics^{36} in the context of two separate reservoirs of gas, each held in equilibrium at different temperatures and densities. In this article, we consider a dilute gas in a volume bisected by a simple porous transpiration interface, with the system driven out of equilibrium by inducing a temperature gradient. The partitioning interface is taken to be a thin adiabatic wall, with holes that are small compared with the mean free path of the gas. This ensures that the pores themselves do not have a significant hydrodynamic effect, yielding a simple system for the study of correlations. Particles passing through the interface conserve their energy, but not their momentum, consistent with the particles experiencing randomizing elastic collisions as they pass through the membrane.

In this article, both DSMC and FHD are used to simulate this system, yielding information on long range correlations and validating the use of FHD for transpiration interfaces. In Sec. II, we begin with a brief review of DSMC and of the theory underlying fluctuating hydrodynamics including the numerical implementation of the FHD equations. In Sec. III, we present several models for transpiration interfaces, using both master equation and Langevin approaches in the FHD context, and the kinetic model applicable to DSMC simulations. In Sec. IV, we validate these approaches by comparing the long range correlations occurring in both FHD and DSMC simulations, in the presence of interfaces with varying porosity. The role of the interface in distorting and suppressing these correlations is discussed, and a brief comparison of the relative computational performance of the DSMC and FHD simulations is made. Some conclusions are given in Sec. V.

## II. FLUID SIMULATION MODELS

For the study of nonequilibrium fluctuations, we use two independent approaches for the bulk fluid: a particle model, Direct Simulation Monte Carlo (DSMC), and a continuum model, fluctuating hydrodynamics (FHD). This section summarizes these numerical models; the treatment of the transpiration interface is described in Sec. III.

### A. Direct simulation Monte Carlo (DSMC)

Direct simulation Monte Carlo is a well-known method developed by Graeme Bird for computing gas dynamics from the perspective of kinetic theory.^{13} As in molecular dynamics, the state of the system in DSMC is given by the positions and velocities of particles. At each time step, the particles are first independently moved ballistically, imposing boundary conditions if they strike a surface. Generally, isothermal boundary conditions are implemented using a diffuse surface interaction, while adiabatic boundaries are simulated using specular reflection; surfaces with partial accommodation can also be simulated. Collisions are evaluated by a stochastic process, selecting the post-collision angles from their kinetic theory distributions while conserving momentum and energy. We employ the hard sphere collision operator in the present work, but many other collision models are available in DSMC. Pedagogical explanations of DSMC are contained in Refs. 37–39 with a complete description given in Ref. 14.

We note that the statistical variation of the hydrodynamic variables arises from the inherently discrete, particle character of the fluid. Consequently, although DSMC is a Monte Carlo algorithm, randomness is not required to capture the statistical variation of the hydrodynamic; the same fluctuations also occur in deterministic methods, such as molecular dynamics. Under a relatively weak set of conditions (e.g., detailed balance), discrete event processes reproduce the correct physical fluctuations at hydrodynamic scales.^{40,41} When each DSMC particle represents one molecule in the gas, the magnitude of the hydrodynamic fluctuations is in agreement with equilibrium statistical mechanics.^{42} For both equilibrium and nonequilibrium systems, DSMC yields the physical spectra of spontaneous fluctuations as predicted by theory^{43,44} and observed in experiments.^{45,46} In this paper, the numerical measurements from DSMC simulations are used to validate the fluctuating hydrodynamic simulation results.

### B. Fluctuating hydrodynamics (FHD)

The equations of fluctuating hydrodynamics are obtained by adding white noise to each dissipative flux in the deterministic Navier-Stokes equations^{47} to represent thermal fluctuations. We consider the compressible fluctuating hydrodynamic equations for a pure monatomic ideal gas; generalizations are found in Refs. 7, 24–32. The continuity, momentum, and energy conservation equations for FHD are

Here, *ρ* is the mass density of the gas, **u** is the velocity, *P* is the pressure, *T* is the temperature, and $E=cVT+12|u|2$ is the total specific energy. For a monatomic ideal gas, the specific heat at constant volume is $cV=32kB/m$, where *k*_{B} is Boltzmann’s constant and *m* is the molecular mass; the equation of state is *P* = *nk*_{B}*T*, where *n* = *ρ*/*m* is the number density. From Eqs. (1)–(3), we define the fluxes in the mass, momentum, and total energy as *J*_{M}, $JP$, and $JE$, respectively.

In the momentum and energy equations [Eqs. (2) and (3)], the deterministic viscous stress tensor **Π** and deterministic heat flux vector **Q** are given by

where *δ*_{ij} is the Kronecker delta, *η*(*T*) is the shear viscosity, and *κ*(*T*) is the thermal conductivity; consistent with a monatomic gas, we take the bulk viscosity to be zero.

The terms that appear with a tilde in Eqs. (2) and (3) denote stochastic contributions from Gaussian random fields. The stochastic viscous tensor, $\Pi \u0303$, and the stochastic heat flux tensor, $Q\u0303$, have zero mean and covariances^{17,18}

and

where *δ*(..) is the Dirac delta function.

## III. TRANSPIRATION INTERFACE MODELS

As discussed in Sec. I, we consider an adiabatic, stationary interface that is porous to the gas molecules (e.g., a specular wall with small holes). Particles encountering the interface are allowed to pass with an effusion probability, *f*. This probability can be viewed as the fraction of the wall area comprised of holes. The conditions in the gas are chosen such that a particle that crosses the interface experiences a large number of intermolecular collisions before crossing again—this is called the effusion condition.^{13,48}

### A. DSMC interface model

The implementation of the transpiration interface is straightforward in DSMC. Particles reaching the interface are transmitted with probability *f* and specularly reflected with probability 1 − *f*. The transmitted particles are assumed to interact with the adiabatic porous interface such that the particle’s energy, *ϵ*, is conserved but not its momentum. Therefore, the velocity of a transmitted particle, **v**, is reset by randomly selecting it from a hemisphere of radius $|v|=2\u03f5/m$.

### B. FHD interface models

The specular condition for the surface of the interface, discussed above, is modeled hydrodynamically as a full slip adiabatic boundary condition. We present two methods for modeling the effusion of mass and energy through the interface, based on master equation and Langevin approaches; as discussed above, the effusive momentum flux is set to zero. We first discuss the kinetic theory underlying both these techniques.

#### 1. Kinetic theory

For a gas at equilibrium, the time between particle crossings through a transpiration interface of area *A* is exponentially distributed, with the mean time between crossings given by^{13}

By comparison, the mean time between collisions for a hard sphere particle^{13} is

where *d* is the collisional diameter of the particle. The effusion condition described above is *nAλ*_{c}*τ*_{e} ≫ *τ*_{c}, where *λ*_{c} is the mean free path, which is satisfied when *f* ≪ 1. The probability distribution for the kinetic energy, *ϵ*, of particles passing through the interface is

this is a gamma distribution, Γ(2, *k*_{B}*T*).^{49} Equations (7) and (9) form the basis of the stochastic simulation algorithm given in Sec. III B 2.

From the above, we find that the probability per unit time of observing a particle of energy *ϵ* crossing the interface is

where the arrows indicate the direction of the crossing and the subscripts *L* and *R* refer to the gas on the left and right hand sides of the interface. This leads to the master equation for $P(M,E)$, i.e., the probability density of the state on the right of the interface with mass *M* and total energy $E$,

Following from the above,^{36} the average one-way flux of mass passing from left to right through the interface is

with the arrow indicating the direction of flux. The total average mass flux is given by

The average energy of a particle crossing the interface is 2*k*_{B}*T*, so the average total energy flux is

As discussed above, we assume that there is no momentum flux across the interface, so $JP=0$. Since the nonequilibrium conditions considered in this article result in zero average net momentum transport, this difference in treatment between DSMC and FHD appears to have a negligible impact on the long-ranged correlations of interest (see Sec. IV C).

The variances and covariance of the mass and energy fluxes are

The resulting Langevin equations are given by

where $J\u0303M$ and $J\u0303E$ are Gaussian noise terms, with

These equations form the basis of the Langevin method discussed in Sec. III B 3.

#### 2. Master equation implementation

Following from the theory presented above, we define a stochastic simulation algorithm (SSA)^{50} implementation of the interface. This operates by repeatedly generating particle crossing times from an exponential distribution, until the cumulative sum of these crossing times exceeds a single FHD time step, Δ*t*. The energy of each crossing particle is selected from the appropriate gamma distribution, and a tally of the total mass and energy transfer is used to update the hydrodynamic cells on each side of the interface.^{51}

The implementation of the master equation model over a single time step is as follows:

Take

*ρ*_{L}and*ρ*_{R}to be the density at the cell centers to the left and to the right of the interface. Similarly, the temperature and specific energy in said cells is indicated by*T*_{L},*T*_{R},*E*_{L}, and*E*_{R}.From Eq. (7), the rate of particles crossing from each direction is

*r*_{L}= 1/*τ*_{e}(*ρ*_{L},*T*_{L}) and*r*_{R}= 1/*τ*_{e}(*ρ*_{R},*T*_{R}).Initialize the accumulation of mass and energy flux,

*F*_{M}and $FE$, across the interface to zero.Set the interface time counter

*t*_{ssa}= 0. Loop while*t*_{ssa}< Δ*t*:- Given an exponential distribution of time between particle crossings, the time until the next particle crossing is given bywhere $Ri\u2208(0,1]$ indicates an independent uniform random number. [We note that Eq. (21) is a standard approach for generating samples from an exponential distribution based on the inverse cumulative distribution function.(21)$\theta =\u2212log(R1)rL+rR,$
^{38,49}]where $T*$ =(22)$\u03f5=\u2212kBT*\u2061log(R3R4),$*T*_{L}if the particle is traveling left to right or $T*$ =*T*_{R}otherwise. Increment the interface time counter,

*t*_{ssa}+=*θ*. If*t*_{ssa}≥ Δ*t*, then break out of the loop.Select the direction for the particle crossing the interface. If $rL/(rL+rR)>R2$, then the particle crosses from left to right; otherwise, it crosses from right to left.

Following from Eq. (9), the energy of the effusing particle is gamma distributed and generated by

^{49}If the particle is traveling left to right,

*F*_{M}+=*m*and $FE\u2009+=\u2009\u03f5$. Otherwise,*F*_{M}−=*m*and $FE\u2009\u2212=\u2009\u03f5$.Finally, the left cell is updated as

*ρ*_{L}−=*F*_{M}/Δ*V*and $\rho LEL\u2009\u2212=\u2009FE/\Delta V$, and the right cell is updated as*ρ*_{R}+=*F*_{M}/Δ*V*and $\rho RER\u2009+=\u2009FE/\Delta V$ with Δ*V*being the cell volume.

Note that the interface is modeled as a slip wall since the flux of momentum crossing the interface is zero.

#### 3. Langevin equation implementation

In contrast to the SSA approach, which considers discrete particle crossing events, the Langevin algorithm generates the total mass and energy transfer based on the expected values. A Gaussian approximation is used for the fluctuations in these quantities. This is similar in principle to the stochastic fluxes generated between cells in the FHD algorithm, noting that the central limit theorem shows the equivalence of event-driven and Langevin approaches for large numbers of events. The theoretical foundation for this Langevin formulation comes from the Fokker-Planck equation that results from a second-order truncation of the Kramers-Moyal expansion of the master equation.^{52,53}

The implementation of the Langevin model over a single time step is as follows:

Take

*ρ*_{L}and*ρ*_{R}to be the density at the cell centers to the left and to the right of the interface. Similarly, the temperature and specific energy in said cells is indicated by*T*_{L},*T*_{R},*E*_{L}, and*E*_{R}.Calculate the mean mass and energy fluxes given by Eqs. (14) and (15) and the variances and covariance given by Eqs. (16)–(18).

Calculate the correlation coefficient, given by

Generate correlated random variables

where $G1$ and $G2$ are independent Gaussian random numbers of mean zero and variance one. [We note that Eq. (24) is a standard approach for generating a bivariate Gaussian with a desired covariance from independent Gaussian random variables.^{49}]

Following from Eq. (19), the mass and energy flux accumulations are given by

Finally, the left cell is updated as

*ρ*_{L}−=*F*_{M}/Δ*V*and $\rho LEL\u2009\u2212=\u2009FE/\Delta V$, and the right cell is updated as*ρ*_{R}+=*F*_{M}/Δ*V*and $\rho RER\u2009+=\u2009FE/\Delta V$ with Δ*V*being the cell volume.

Again, note that the interface is modeled as a slip wall since the flux of momentum crossing the interface is zero.

## IV. SIMULATION RESULTS

In this section, we begin by validating the use of FHD for transpiration by comparison of the FHD and DSMC simulations. For the majority of this section, the master equation method of Sec. III B 2 is used to implement the interface in the FHD simulations; several comparisons are shown in Sec. IV C to validate the Langevin method of Sec. III B 3. The effect of transpiration is demonstrated with several different effusion probabilities, *f*, by examining the long range spatial correlations in the fluid. This is followed by an analysis of the relative computational efficiency of the DSMC and FHD simulations.

We consider a quasi-one-dimensional domain with a gradient induced in the gas by setting different temperatures at *x* = 0 and *x* = *L* (see Fig. 1). These temperatures are set in the DSMC simulations using a fully diffuse boundary condition and in the FHD simulations with an isothermal no-slip boundary condition. Rather than using a temperature slip boundary condition, the Dirichlet boundary condition for temperature is set to match the mean extrapolated temperatures from the DSMC simulations—this is discussed further below. For the DSMC simulations, the *y* and *z* boundaries are taken to be periodic; as discussed above, the FHD simulations are performed in 1D. The transpiration interface is at *x* = *L*/2, the center of the domain. Simulation parameters are given in Table I.

Molecular diameter | 3.66 × 10^{−8} |

Molecular mass | 6.63 × 10^{−23} |

Reference temperature | 273 |

Sound speed | 30, 781 |

Specific heat c_{$v$} | 3.12 × 10^{6} |

Reference mean free path | 6.26 × 10^{−6} |

Reference mean free time | 1.64 × 10^{−10} |

System length (x dir) | 2.5 × 10^{−4} |

Interface x location | 1.25 × 10^{−4} |

Cell size Δx | 3.125 × 10^{−6} |

Time step Δt | 1.0 × 10^{−12} |

Number of cells (x dir) | 80 |

System height, width (y, z dir) | 1.25 × 10^{−6} |

Average particles per cell | 131 |

Molecular diameter | 3.66 × 10^{−8} |

Molecular mass | 6.63 × 10^{−23} |

Reference temperature | 273 |

Sound speed | 30, 781 |

Specific heat c_{$v$} | 3.12 × 10^{6} |

Reference mean free path | 6.26 × 10^{−6} |

Reference mean free time | 1.64 × 10^{−10} |

System length (x dir) | 2.5 × 10^{−4} |

Interface x location | 1.25 × 10^{−4} |

Cell size Δx | 3.125 × 10^{−6} |

Time step Δt | 1.0 × 10^{−12} |

Number of cells (x dir) | 80 |

System height, width (y, z dir) | 1.25 × 10^{−6} |

Average particles per cell | 131 |

### A. Temperature profiles and effusion condition

In Fig. 2, the temperature profiles are compared between the DSMC and FHD simulations (using the master equation algorithm of Sec. III B 2). Comparisons are shown for a range of *f* values and in the absence of a membrane. In all cases, the temperature boundary conditions for the DSMC simulations are set to 273 K and 519 K on the left and right hand boundaries, respectively, giving a large gradient of 7604 K/cm. Due to this large gradient, the gas at the boundary is unable to completely accommodate to the wall temperature, resulting in a significant Knudsen layer^{16} and “temperature slip.” This feature is correctly captured by the DSMC simulations but must be accounted for in the FHD simulations with a modified boundary condition. As mentioned above, the FHD boundary conditions have been extrapolated from these new temperatures; see Table II.

f
. | T(0)
. | T(L)
. |
---|---|---|

No membrane | 280.8 | 508.5 |

0.1 | 277.6 | 509.9 |

0.05 | 276.2 | 515.4 |

0.025 | 275.0 | 516.9 |

0.0125 | 274.1 | 517.9 |

0.006 25 | 273.6 | 518.4 |

f
. | T(0)
. | T(L)
. |
---|---|---|

No membrane | 280.8 | 508.5 |

0.1 | 277.6 | 509.9 |

0.05 | 276.2 | 515.4 |

0.025 | 275.0 | 516.9 |

0.0125 | 274.1 | 517.9 |

0.006 25 | 273.6 | 518.4 |

In all cases, the DSMC and FHD simulations show good agreement aside from a small discrepancy near the interface, which becomes prominent at *f* ≳ 0.1. As discussed in Secs. I and III, the transpiration models introduced for the FHD simulations are valid only when the effusion condition holds. Similar to the temperature boundary condition discussed above, this is highlighted by the appearance of a Knudsen layer at the surface of the interface, with the relevant Knudsen number being proportional to the effusion probability *f*. For *f* ≳ 0.1, this causes the solution near the interface to be poorly approximated by the FHD simulations. This occurs because in this regime there is a substantial chance of a particle crossing the interface in one direction and being reflected back in the other direction after only a small number of collisions, preventing the regions close to the interface from maintaining local equilibrium. This effect is accurately captured by DSMC, which does not rely on a local equilibrium assumption, but not by the FHD simulations.

### B. Measurement of hydrodynamic fluctuations

The fluctuations of the hydrodynamic variables are measured by accumulating statistical samples of the conserved variables. For example, in DSMC, the mass and momentum densities in the cell centered at *x* at time *t* are

where *N* is the number of DSMC particles in the cell and **v**_{i} is the velocity of particle *i*. In FHD, these conserved densities, *ρ* and **p** = *ρ***u**, are computed directly [see Eqs. (1)–(3)]. Means and second moments are computed from these statistical samples, for example,

and

where *S* is the number of time samples. Correlations of hydrodynamic fluctuations are obtained from these conserved variables, e.g.,

The calculation of temperature fluctuations is similar; see Ref. 54 for details. For all the results shown in this article, the simulation was run until a statistical steady state was reached, and recording of the quantities of interest was then begun. The value *S* is given in the caption for each result. In practice, each statistical sample was constructed from an ensemble of shorter simulations to give a total of *S* time samples.

### C. Hydrodynamic fluctuations in DSMC and FHD

In Fig. 3, long range spatial correlations, in the absence of an interface, are compared between DSMC and FHD simulations. Three correlations are shown, $\u27e8\delta ux*\delta \rho \u27e9$, ⟨*δT***δρ*⟩, and ⟨*δT***δT*⟩, in (a), (b), and (c), respectively. In each case, the variable with the asterisk is measured at the fixed location *x** = *L*/4 − Δ*x*/2, i.e., in the cell near the center of the left hand side of the domain; see Fig. 1. Note that in the absence of a temperature gradient (i.e., at thermodynamic equilibrium) $\u27e8\delta ux*\delta \rho \u27e9=\u27e8\delta T*\delta \rho \u27e9=0$ and ⟨*δT***δT*⟩ = 0 for *x* ≠ *x**.

We see good agreement between the FHD and DSMC results, observing that the FHD simulations have difficulty capturing the sharp peak in $\u27e8\delta ux*\delta \rho \u27e9$; this confirms the observations of Ref. 21. In (c), the expected Kronecker delta peak at *x* = *x** has been removed, from both the FHD and DSMC results, for ease of viewing. The FHD simulations again have difficulty in this case and exhibit significant undershoot in the two adjacent cells. This is a well known effect,^{55} and these data points have also been removed from the FHD results for ease of viewing—they are shown on the inset plot.

In Figs. 4 and 5, the same correlations are shown in the presence of interfaces with effusion probability *f* = 0.1 and *f* = 0.05, respectively. In the FHD simulations, the master equation method of Sec. III B 2 has been used to simulate the interface. For the *f* = 0.05 case, subject to the same caveats as above, good agreement is observed between the FHD and DSMC results in most of the domain, with some scattering in the FHD results in the vicinity of the interface. In the *f* = 0.1 comparison, a small systematic difference is observed between the two sets of results, particularly visible in Figs. 4(a) and 4(c)—this represents the breakdown of the effusion condition.

In Fig. 6, the master equation and Langevin approaches for modeling the interface, given in Secs. III B 2 and III B 3, are compared. The correlations ⟨*δT***δρ*⟩ and ⟨*δT***δT*⟩ are shown for the cases *f* = 0.1 and *f* = 0.006 25. Excellent agreement is observed, demonstrating the equivalence of the two approaches.

As can be seen from Figs. 3–5, the general impact of the interface is to reduce the magnitude of long range correlations in the gas. The value $\u27e8\delta ux*\delta \rho \u27e9$ approaches zero for *x*/*L* > 0.5 because the interface does not conserve momentum. Interestingly, a new peak in ⟨*δT***δT*⟩ develops immediately to the left of the membrane—this behavior is borne out for lower values of *f* not shown here. This distortion is also observed in a simple Fourier model, with the interface modeled as a region of low thermal conductivity.

Aside from these features, the reduction in magnitude appears to arise principally from the change in temperature gradient. In Fig. 2, a reduction in the gradient is visible with decreasing *f*, where the interface behaves increasingly like an adiabatic surface. In Fig. 7, the magnitude of correlations is shown as a function of the gradient. The top row shows measurements of the value $|\u27e8\delta T*\delta \rho \u27e9|\xaf$, where the overline indicates that the value has been smoothed by averaging over a spatial interval (see the figure caption). Figure 7(a) shows measurements from the left hand side of the interface, while Fig. 7(b) shows the values from the right. In the bottom row, Figs. 7(c) and 7(d) show the left and right side measurements of $|\u27e8\delta T*\delta T\u27e9|\xaf$. The results displayed are taken from DSMC and FHD simulations in the range 0.006 25 ≤ *f* ≤ 0.1, and the subsets shown in each plot are those where a statistically converged solution was obtained. Where both DSMC and FHD results are shown, the DSMC solution should be taken to be the more accurate; we have included the FHD results to illustrate the divergence of the solutions due to the violation of the effusion condition. For lower values of *f*, the magnitude of the correlations becomes smaller relative to equilibrium noise, requiring a larger sample size to obtain a result. FHD results were obtainable in this regime because this approach exhibits greater computational speed; this is discussed further in Sec. IV D.

In Figs. 7(a) and 7(b), ⟨*δT***δρ*⟩ exhibits a linear response with respect to the gradient; this agrees with the findings of Refs. 43 and 56, where a linear response to the temperature gradient is observed in the absence of a membrane. In the same references, a quadratic response is shown in ⟨*δT***δT*⟩. This appears to be supported by Fig. 7(d); in Fig. 7(c), measurement of this effect is confounded by the appearance of the additional peak near the surface of the interface, discussed above. Finally, it is clear from Figs. 7(b) and 7(d) that the fluctuations on opposite sides of the interface remain correlated even for small values of *f*.

### D. Comparative performance

The cell size used to make the comparisons shown in Figs. 3–5, given in Table I, has been selected based on the kinetic length scale that constrains the DSMC simulations. The time step has been selected to minimize the chance of fluctuations causing negative density or energy values in the FHD simulations—these may also be avoided by increasing the cell size, which is not bound by the kinetic scale as in the DSMC simulations.

In Fig. 8, the effect of changing cell size and time step is compared between the two methods, using the case *f* = 0.1. To maintain a constant diffusive stability value, *C* ∼ Δ*t*/Δ*x*^{2}, for every doubling of the cell size, the time step is increased by a factor of four. Due to the very small time steps used, the error appears to be dominated by the cell size, with both methods showing the expected quadratic convergence. The FHD algorithm is deterministically third order in time and second order in space and weakly second order for stochastic systems; DSMC is at best second order. We note also that, depending on the simulation parameters, one FHD time step executes approximately 10–30 times faster than the equivalent DSMC time step.

## V. CONCLUDING REMARKS

In this article, we have investigated the effect of a simple transpiration interface on long range hydrodynamic correlations in a dilute gas, where a nonequilibrium steady state has been induced by a temperature gradient. Several qualitative differences in the correlations were observed; however, the principal effect of the interface was to reduce their magnitude proportionally to the effusion probability, *f*. This outcome can largely be explained by changes in the temperature gradient induced by the interface. We also observe a distortion of the ⟨*δT***δT*⟩ correlations, specifically the appearance of a new peak located near the interface.

The investigation used two distinct numerical techniques, DSMC and FHD, which were found to be in agreement for appropriate values of *f*; for *f* ≳ 0.1, noncontinuum effects reduce the accuracy of the FHD method. The FHD simulations used both master equation and Langevin approaches for modeling the interface, with the two methods producing equivalent results. This validates the use of FHD for interfaces with low effusion probabilities, giving a computationally faster alternative to DSMC.

There are a number of potential extensions to the work presented here. Currently, both methods used to model the interface in the FHD simulations assume $JP=0$, i.e., there is zero momentum flux. Future extensions to the interface model could include momentum fluctuations around a zero mean, a more realistic model of thermal processes, and incorporation of multicomponent effects such as species dependent permeability. It would also be interesting to use DSMC measurements of the velocity distribution function near the interface to quantify the departure from local equilibrium. This information could aid in the reformulation of the master equation transition probabilities to incorporate nonequilibrium effects.

These developments would provide the starting point for numerical modeling of a range of membrane behaviors. More complex models could include models for transport of ions through membranes and the incorporation of chemical reactions that enable active transport. These capabilities would enable the methodology be used to model transport in biological membranes and transport process important to artificial photosynthesis.

## ACKNOWLEDGMENTS

The authors wish to acknowledge many years of fruitful discussions with Graeme Bird. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

## REFERENCES

_{2}subjected to temperature gradients with and without the influence of gravity