{'fraction_non_alphanumeric': 0.03283322900470489, 'fraction_numerical': 0.012297316398441128, 'mean_word_length': 4.708818443804034, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Delineation of the flow and mixing induced by Rayleigh-Taylor instability through tracers

Rayleigh-Taylor-instability(RTI) induced flow and mixing are of great importance in both nature and engineering scenarios. To capture the underpinning physics, tracers are introduced to make a supplement to discrete Boltzmann simulation of RTI in compressible flows. Via marking two types of tracers with different colors, the tracer distribution provides a clear boundary of two fluids during the RTI evolution. Fine structures of the flow and thermodynamic nonequilibrium behavior around the interface in a miscible two-fluid system are delineated. Distribution of tracers in its velocity phase space makes a charming pattern showing quite dense information on the flow behavior, which opens a new perspective for analyzing and accessing significantly deep insights into the flow system. RTI mixing is further investigated via tracer defined local mixedness. The appearance of Kelvin-Helmholtz instability is quantitatively captured by mixedness averaged align the direction of the pressure gradient. The role of compressibility and viscosity on mixing are investigated separately, both of which show two-stage effect. The underlying mechanism of the two-stage effect is interpreted as the development of large structures at the initial stage and the generation of small structures at the late stage. At the late stage, for a fixed time, a saturation phenomenon of viscosity is found that further increase of viscosity cannot see an evident decline in mixedness. The mixing statues of heavy and light fluids are not synchronous and the mixing of a RTI system is heterogenous. The results are helpful for understanding the mechanism of flow and mixing induced by RTI.


I. Introduction
Rayleigh-Taylor instability (RTI) occurs at the perturbed interface of two fluids when the pressure gradient points from the heavy to the light fluid [1][2][3] . RTI widely exists in both scientific and engineering scenarios, including supernova evolution in astrophysics 4 , stratigraphic dynamics in geophysics 5 , and material mixing in chemical engineering 6 . The inertial confinement fusion (ICF), which is expected to be a viable approach to develop fusion energy, control of RTI is important for the successful ignition. If it allows free development of RTI, the applied compression is limited and the integrity of the fuel shell will be badly destroied, leading to the failure of ignition [7][8][9][10] . Although with the importance, the control of RTI is still a challenge due to unclearness of this strong nonlinear process, and better understanding of RTI will benefit both fundamental science and engineering practice 11,12 . Flow-induced by RTI is time-dependent. It undergoes different development stages, with both linear and nonlinear features, showing typical patterns of fingering and mixing of two fluids. In some cases, the flow can transit to a fully turbulent state. The RTI phenomenon was observed and upward to a scientific level by Rayleigh 1 , and further analyzed by Taylor 2 , a) Corresponding author: xu_aiguo@iapcm.ac.cn Lewis 3 , and et al. The concept of RTI was invented after then and a physical quantity of RTI, the growth rate was derived from physical modeling 13 . However, it is hard to generalize the whole process of RTI with a simple theory. Most theoretical works only focused on either the early linear (weakly non-linear) growing stage or the fully turbulent stage [14][15][16][17] . Intervals between stages have yet not well understood 18 . Especially, the accompanied Kelvin-Helmholtz instability (KHI) induces large-scale vortices and mixing, which will enhance the non-linearity and make the instability even harder for the theoretical analysis. Besides theoretical analysis, experimental and numerical efforts was also made for better understanding RTI. In recent experiments, advanced apparatuses were employed, including planer laser-induced fluorescence to capture the RTI-induced mixing layer Waddell, Niederhaus, and Jacobs 19 , Wilkinson and Jacobs 20 and high-speed video camera to visualize the abating effect of RTI on Richtmyer-Meshkov instability (RMI) by Luo et al. 21 . Although experiments provide good insight into the physics of RTI, the disadvantage is relatively expensive and time cost. The obtained information is limited by the ability of measuring techniques. Therefore, numerical simulation are widely adopted in RTI researches. Different methods in computational physics, including level-set method 22 , front tracking method 23 , volumeof-fluid method 24 , smooth particle hydrodynamics 25 , largeeddy simulation 26 , phase field method 27 , direct numerical simulation 28 , and etc., were already utilized for RTI simula-tion. With consideration of stabilizing and destabilizing effects on RTI, viscosity, compressibility, surface tension, and etc., were investigated in previous simulations [29][30][31][32][33] .
Generally, there are three kinds of fluid dynamic models in instability researches: the macroscopic model, mesoscopic model, and microscopic model. The macroscopic model, such as Euler and Navier-Stokes (NS), is based on the continuity hypothesis. According to the Chapman-Enskog theory, the NS includes only the first-order thermodynamic non-equilibrium (TNE) effect. This kind of model has been longly and mostly applied in fluid instability research. However, in complex fluid systems with mid-scale structure and material interfaces, the mean free path (or mean relaxation time) cannot be ignored compared to characteristic scale (or characteristic time) of the system, which means that discrete (or nonequilibrium ) effect becomes pronounced and consequently challenges the physical rationality of NS model. The microscopic models such as molecular dynamics simulation are the most precise method in capturing information for flows, but they are restricted by the computational cost in spatiotemporal scales. Recently, the mesoscopic discrete Boltzmann model (DBM) was proposed to deal with such mid-scale structure behaviors and is capable of simulation of a relatively wider spectrum of spatiotemporal scale. Via selecting appropriate kinetic moments to keep values, it can be designed according to the extent of TNE or discreteness which is aimed to investigate. The DBM describes the system from a wider research perspective. A DBM is equivalent to a traditional macroscopic fluid dynamic model plus a coarse-grained model describing TNE effects 34 . Currently, DBM has made significant progress in various hydrodynamic instabilities including RTI, RMI, and KHI, which brought a series of new insights of abundant TNE behaviors [34][35][36][37][38][39][40] . Based on a multiple relaxation time DBM, besides investigating Prandtl number effects on RTI, Chen, Xu, and Zhang 37,38 uncovered the collaboration and competition roles between RMI and RTI and found high correlations between the nonuniformities of macroscopic quantity and nonequilibrium strength in the flow system. Zhang et al. 40 showed that via interfaces captured by appropriate TNE properties, the mixing of materials and that of internal energy resulted from hydrodynamic instability can be more conveniently probed. For example, the newly observed double-layer vortex has the correlation and difference of density mixing and temperature mixing/unification. Lai et al. 36 proposed a DBM model with acceleration term for RTI flow simulation and investigated the effect of compressibility on RTI, finding that the compressibility had a twostage effect that delays RTI growing at the initial stage and accelerates it at the later stage. Meanwhile, they presented an interface-tracking method through TNE characteristics.
Along with RTI flow, material mixing is a focus due to many applications 41 . The investigation of RTI mixing in complex flow systems requires a basic metric defining the mixing state. Several metrics were employed in previous researches. A coarse-grained indicator is the width of mixing layer defined as the distance between the amplitude of the spike and the bubble [42][43][44][45] , but the results and conclusions based on this metric may differ due to different definition of the edge of mixing layer 46 . To obtain a fine-grained description of mixing, Zhou et al. 47 proposed actual mixed mass which is decided by the local mass fraction of the two fluids and Ma et al. 48 applied this metric in an advanced RTI experiment. Recently, Lin et al. 39 studied the mixing state through mixing entropy which was also determined by mass fraction of each fluid based on a two-component DBM.
Based on a previous work of discrete Boltzmann investigation of RTI made by Lai et al. 36 , this paper intends to deepen the study from two aspects, the development of fluids interface and late-stage mixing. Correspondently, we will seek an interface-tracking method through embedding tracers with DBM and study the process of mixing in RTI. One-way coupled tracers will be introduced to provide a synchronous Lagrangian viewpoint to the Eulerian flow field, and these tracers will open up a new angle for RTI analysis. The rest of this paper is organized as follows. The numerical method and physical modeling are presented in Section §II. Simulation results are presented in Section §III along with discussions on RTI flow and mixing considering the effects of compressibility and viscosity in Subsection §III D and Subsection §III E respectively. Finally, conclusions are made in Section §IV.

II. DBM for RTI flow coupled with tracers
Regarding a complex flow system induced by RTI, though has been investigated for many years, it still faces two challenges. The first is how to model the process objectively. An accurate simulation is based on the effectiveness of the physical model and corresponding functions on the specific spatial and time scale. The error generated from physical modeling cannot be compensated for by any modification on algorithm. Secondly, it is a challenge to extract useful information from massive data generated from simulation. Only regarding these two challenges, it is possible to accurately model RTI flow and effectively obtain physics from the results.
A. The construction of DBM DBM is a coarse-grained physical model based on nonequilibrium statistical physics. It selects a set of kinetic properties decided by corresponding kinetic moments to study the system. Specifically, the deeper the non-equilibrium is, the more kinetic moments of distribution function f are needed. The initial several conserved kinetic moments(including density, momentum, and energy) are necessary to any nonequilibrium flows, whereas the non-conserved kinetic moments depend on the degree of non-equilibrium which is aimed to study. The Navier-Stokes equations are the case where the system deviates slightly from its thermodynamic equilibrium called the first-order non-equilibrium effect. In the case without external force, from the Boltzmann equation to a DBM, there are three steps: (i) simplification of the collision term, (ii) discretization of the particle velocity space, (iii) characterizing the nonequilibrium state and picking out meaningful nonequilibrium information. The first two steps are based on the properties of the system and seize the main contradiction. The third step is based on simulation results and seize the main contradiction. In the first two steps, the kinetic moments (or hydrodynamic quantities) used to describe the flow system must keep the same values before and after the simplification/discretization. In this way, the model is significantly simplified while the main features of described flow system are kept unchanged 34,49 .
To obtain a DBM for RTI flow, the acceleration effect of gravity should be taken into account. A DBM with acceleration effect was proposed by Lai et al. 36 , through approximately replacing the distribution function f in Boltzmann equation with its corresponding equilibrium distribution function f eq . The simplified BGK-Boltzmann equation 50 with a force term 51 is written as follow, where f = f (r, v,t) is the distribution function of fluid molecules and f eq is the equilibrium distribution function. r indicates spatial position, v is the fluid velocity, t indicates the time, a corresponds to external body force, R is the gas constant, T is the temperature, and τ is the relaxation time. Afterward, the continuum particle velocity space will be discretized. Through Chapman-Enskog expansion, a discrete Boltzmann model with at least 16 discrete velocities can cover fully two-dimensional compressible NSE on the description of fluid dynamics. According to previous researches, DBM provides credible result of flow simulation compared with classical solvers of NSE [52][53][54][55] .
Another feature of DBM is that it provides a quantitative indicator of local thermodynamic nonequilibrium effects by defining nonequilibrium moments ∆ * m,n as 34,49,56 , where M * m,n ( f ) represents kinetic central moments of f , defined as v * m v * · · · v * f dv * in continuum velocity space, in which used the relative velocity v * as v * = v − u. m indicates the number of velocity used in the moment and n indicates the tensor order respectively. In discrete velocity space, it only needs to replace the integration and continuum velocity with summation and discrete velocity.
The moments introduced are interrelated and constitute a relatively complete description of the system's nonequilibrium states and behaviors. The collaboration of TNE quantities given by DBM presents a rough but quantitative indicator for the specific nonequilibrium state of a fluid system. Each independent component of the TNE quantities describes the nonequilibrium state from its own perspective. For example, ∆ * 2 indicates viscous stress tensor and ∆ * 3,1 indicates heat flux. ∆ * 3 and ∆ * 4,2 are higher-order non-equilibrium indicators beyond the traditional NS model and indicate the flux of viscous stress and of heat flux, respectively.
In order to initiate a RTI flow, initial and boundary conditions are set according to Lai et al. 36 and Scagliarini et al. 57 . The acceleration directs from top to bottom and the denser fluid locates above the lighter fluid. Initially at the interface, a sinusoidal disturbance occurs with the form β c (x) = β 0 · cos (kx), where β is the relative coordinate of the interface, β 0 is the amplitude of perturbation, and k denotes the wavenumber related to wavelength as k = π/λ .

B. Tracers coupled in RTI flow
Tracers delineating fluid flow belong to a special kind of particles. Particles immersed in the fluid always have interaction with the surrounding phase but not necessarily have prominent influence on it. This is valid by assuming that the inertial effect of particles is extremely weak, and could be understood based on the indication of Stokes number (Stk = u 0 · t 0 /d 0 , where t 0 is the dynamic relaxation time of particle, u 0 is local flow velocity, and d 0 is the diameter of particle). When Stk 1, an immersed particle responds immediately to the surrounding change and moves closely according to local streamline. Therefore, the velocity of tracer is an indication of the surrounding fluid velocity. In order to obtain the tracer velocity from the surrounding fluid, the Dirac function is expected to be used according to Peskin 58 , where u p is the tracer velocity, and r k indicates the spatial position of the k-th tracer. The subscript p denotes tracers and k indicate the tracer number. δ is the Dirac function. The above equation is universally satisfied through the flow field. In the simulation, the flow field is usually discretized on Eulerian grids, so tracers may locate at the interval between grid nodes. Therefore, it is necessary to discretize continuum Dirac function through approximation, where approximate discrete Dirac function ψ replaces the original continuum one. Two-dimensionally, ψ equals to the product of two one-dimensional discrete Dirac functions ϕ, There are many choices for approximate discrete Dirac function, the one selected in this paper is a simplified form as Peskin 58 , which keeps high-order accuracy as well as easy to implement, The schematic of implementing a discrete Dirac function on Eulerian grid and comparison of weights between the simplified form, the continuum form, and a classical form 59 is shown in Figure 1.
After the velocity is determined, the trajectory of each tracer can be captured by motion equation as,  To numerically capture the trajectory of each tracer, a 4thorder Runge-Kutta scheme is utilized for the discretization of the particle motion equation 60 . At each midpoint during Runge-Kutta processing, the tracer velocity is refreshed at the midpoint location, which guarantees a high-order accuracy of the calculation.

III. Results and discussion
A numerical simulation of 2D compressible RTI flow is conducted based on the above modeling and method. Mesh size is set of 1024 × 256 to guarantee enough resolution and to save computational cost. About 250000 tracer particles are randomly distributed to the whole region of simulation. Initial parameters are selected as Table I shows. Periodic boundary conditon and non-flux boundary condition are set to horizontal and vertical boundaries respectively for both fluids and tracers. The interface of the two fluids also separates tracers into two kinds. Tracers above are marked with type-a and tracers below are marked with type-b. This enables to distinguish the mixing of tracers during the simulation. Additionally, tracers set on the interface are marked as type-c to record interfacial deformation and changes of physical quantities on the interface.

A. Tracing RTI evolution
A comparison between evolution patterns of tracers and fluid density contours is shown in Figure 2. Fluid density is utilized for traditional evolution patterns of RTI. In an incompressible immiscible fluids system, the density contour is clear to show distinct difference between fluids in RTI. However, in the miscible and compressible fluids system which is investigated in this paper, the density cannot provide a clear interface due to transition layer of density is formed at the interface [61][62][63][64] . However, by marking two types of tracers with different colors, the tracer distribution provides a clear boundary of two fluids during each time of RTI evolution. From At this stage, the density cannot provide an exact delineation of mixing. Nevertheless, the distinct difference between two kinds of particles is kept during RTI mixing as shown in Figure 2(h). Clearly shown by the particle distribution patterns, fluids are mixed in the late stages and tiny vortex structures are generated. These tiny flow structures are important incentive to further mixing and their contribution will be analyzed in the following sections of this paper. Noted that, due to the discrete nature of tracers, there are some unavoidable white-noisy points on the tracer distribution pattern, but if enough tracers are included, the resolution will not interfere with the observation of the characteristic structures. Compared with a two-component DBM utilized to distinguish RTI mixing by Lin et al. 39 , which nearly doubles the computational cost of a single-component simulation, the current tracer method is more convenient and efficient.
Previously, most researches focused on the evolution of the interface of two fluids, a separated analysis on the heavy and light fluids is seldom. The velocity of each tracer represents local flow velocity and the two kinds of tracers may have different behavior corresponding to the characteristics of heavy and light fluids. A further investigation of the velocity distribution of tracers is meaningful for understanding and distinguishing stages of RTI flow. Therefore, a statistic of tracer velocity is presented and the tracer velocity distribution (TVD) is shown in Figure 3. At t = 0.5, when the flow initiates, distribution of tracers concentrates at the original point where velocity equals zero, as shown by Figures 3(a1) and (a2), indicating that the disturbance from the interface has not propagated long distance in the system. The outlines of two distribution patterns both form ring shapes. Along the outline, the 3(a2)) is denser. During the movement, both heavy and light fluid encounters the hindrance from the other, leading to the compression of fluids at the interface. Gradually, the symmetry between the upper and lower halves of the two kinds of TVD patterns fully diminishes. At t = 1.5, a layered pattern is formed in both TVD, as shown by Figures 3(b1) and 3(b2). The inner structures of two TVDs are different. As shown in Figure 3(b1), a 'drip' develops from the origin and encircled by a 'crown', indicating that a small part of type-a tracers aligned with the spike breaks through the hindrance of surrounding light fluid during RTI. At the same time, two inside structures are formed vertically inside an 'envelope', see Figure 3(b2). As shown by histograms in Figures 3(a1), 3(a2), 3(b1), and 3(b2), both horizontal velocity (u px ) and vertical velocity (u py ) of type-a and type-b tracers distribute normally although with small deviation. At t = 2.5, two petunia-like patterns are shown in Figures 3(c1) and 3(c2). The outlines of two TVD patterns resemble petals with a long 'stamen' in the center, indicating that many tracers are with zero horizontal velocity but the vertical velocity of tracers spans for a large range. The spike maintains its flow direction when it penetrates the light fluid, so the 'stamen' in Figure 3(c1) is slender than that in Figure 3(c2). At t = 2.5, the distribution of u px of both type-a and type-b tracers centralize at original points and are close to normal distribution. However, a distinct doublepeak distribution is observed from u py histogram of type-b tracers. At the time, KHI develops and leads to the vortex motion of both heavy and light fluids. Two fluids roll together in the same vortex. A small swelling of type-a u py histogram indicates that the KHI vortex consists of a small portion of heavy fluids and the double-peak in type-b u py histogram shows that the light fluid constitutes mainly in KHI vortex. At t = 3.5, when there is no distinct spike and bubble in the RTI system, spiral patterns are observed in two TVDs, as shown by Figures 3(d1) and 3(d2). The u px histograms of both the type-a and type-b tracers cover a longer range, so tracers transport randomly in this stage without evident groups. The two peaks in u py histogram of type-b tracers disappear, and approximately normal distribution is formed. Distinct difference is observed between u py histogram of type-a and type-b tracers. It concludes that the mixing statues of heavy and light fluids are not synchronous and the mixing of RTI system is heterogeneous. Here presents the most interesting observation of TVDs, for further understanding, investigation of the dynamic process is still important to be made.

B. Interfacial tracking manifestation
The interface tracking technique presents a problem in complex fluid systems and the tracer method brings an effective measure in interface tracking. Initially, the number of 128 type-c particles are distributed at the intermittent interface. Then, through recording the positions of type-c particles in RTI evolution, we can easily delineate characteristics of the interface during RTI growing. Examples are shown in Figure  4.
Initially, a start-up slight perturbation exists at the interface. At t = 0.2, fingers of two fluids (spike and bubble) are formed at the interface and then penetrate each other. However, the developing rates of the two parts are not consistent. At t = 0.2, the spike crosses the position of y = −0.015 while the bubble just reaches the position of y = 0.015. Then, the spike grows still faster than the bubble, which is indicated by figures of late times as t = 0.4, 0.6, and 0.8. Due to acceleration, the interface continually deforms, while at t = 1.0, the velocity difference on the shear direction of the interface initiates ob- servable KHI, with two large vortices taking their shapes at each side of the spike, which is in accordance with the classical observation made by Kull 13 . The development of KHI is a critical step toward the opening of turbulent mixing 65 and it also accounts for the formation of mushroom-shaped flow structure during the RTI development 66 . Because the interface is stretched by further development of KHI, the limited number of tracers distribute more sparsely at a late time, so it is impossible to sketch the profile of the interface precisely as a decreased number of tracers locate at the front, shown by the interface profile at t = 1.4. After that, the tracer description for the interface is unacceptable. Therefore, it is important to design a tracking approach with appropriate number or concentration of tracer afore the RTI interfacial tracking.
Furthermore, as the main place for hydro-instability initiation and development, it is worthwhile to pay special attention to physical features of the interface, which may also be helpful for controlling instability in ICF 67,68 . Also, the design and improvement of all the interface tracking techniques need a better physical understanding of the interface. Fundamentally, thermodynamic nonequilibrium (TNE) is the driving force of fluid system evolution. The complexity of non-equilibrium behavior determines the diversity of the description of nonequilibrium degree. Besides, the TNE behaviors of the interface are indispensable content in the RTI system.
Fundamentally, Knudsen number, viscosity, heat conductivity, the macroscopic gradients of density, flow velocity, temperature, and pressure are commonly used to characterize the degree of non-equilibrium. They all describe the degree of nonequilibrium of the system from their respective perspectives. But they are also highly condensed and averaged coarse-grained descriptions of certain information. Much specific information about non-equilibrium states, the specific allocation of internal energy in different degrees of freedom, etc., through which they are invisible and cannot be directly studied. Therefore, in addition to the description mentioned above, we need a more detailed description of TNE. In this respect, the non-conservative kinetic moments of ( f − f eq ) (as introduced in Section II A) can be used for a more detailed description of the nonequilibrium behavior 34 .
By extracting the TNE information of type-c tracers through interpolating the TNE information of fluid field, we hence demonstrate the TNE property of interface. In the following analysis of interfacial TNE, the period of t = 0 to 1.2 is selected during which the application of type-c tracers is acceptable. The mean TNE of interface is then calculated as Equation (8), where Λ i represents the length of interface and N is the total number of type-c tracers.
From the averages over the interface of each TNE component (∆ * 2 , ∆ * 3 , ∆ * 3,1 , and ∆ * 4,2 ) in Figures 5(a)-(d), we observe that most curves have a peak value at bout t = 0.1. At this time, there is no obvious distortion at the interface, which means the driving potential is accumulated firstly and then activates the development of instability. Figure 5(a) shows the changes of ∆ * 2,xx , ∆ * 2,xy , and ∆ * 2,yy over time, which corresponds to viscous stress in dimension. During the development, curves of ∆ * 2,xx and ∆ * 2,yy are nearly symmetrical along TNE value of 0, indicating the homogeneous of normal pressure on the interface. During the development, ∆ * 2,xy fluctuates and gradually increases, indicating that the momentum flux exchange on the shear direction is slightly enhanced during interface deformation. In Figure 5 at about t = 0.67, delineating that the vertical energy flux is not dominant over the horizontal flux in the system. From dimensional analysis, ∆ * 4,2 corresponds to the flux of energy flux, which indicates a more complex physical process. It is found that two components of ∆ * 4,2 , ∆ * 42,xx and ∆ * 42,xy , reach their peak values at t = 0.1, and decline until t = 0.5 and 1.1, respectively. Also, both the two components rise during the following time.
In the phase space constituted by each independent component of TNE provided by DBM, any points in this phase space correspond to a specific TNE state. The distance from a specific state to the original point in the non-equilibrium phase space defines the extent of non-equilibrium roughly, see Equation (9).
where D Λ indicates interfacial TNE strength by measuring the distance between a specific point in phase space and the origin point. ∆ * 2 , ∆ * 3 , ∆ * 3,1 , and ∆ * 4,2 are different components of a TNE state. Figure 6 shows the evolution of average TNE strength on interface (D Λ ). The mean TNE strength has its maximum value at the initial stage as t = 0.1, and then gradually decreases. The reason is that TNE triggers the initial instability of the system at the start of the simulation, and then the system evolves spontaneously to a new equilibrium state, so there is a peak of the mean TNE strength at the beginning and gradual decreases subsequently.
The analysis is made out of our most present knowledge, but TNE phenomena are so fruitful that cannot be fully investigated, so we leave it for a future study. In sum, a specific non-equilibrium state requires a comprehensive consideration of various non-equilibrium components, and their collaboration is combined to give a more complete description of the flow system.

C. Delineation of RTI mixing
As well as delineation of complex flows, the tracers provide an approach for analyzing the mixing state of fluids. In order to quantify the extent of mixing, the local mixedness (χ p ) is defined based on the spatial information of tracers as, where c a and c b are local number density of type-a and type-b tracers respectively. When there is no type-a or type-b tracers locally , χ p equals to the minimum value 0, while if c a = c b , χ p equals to the maximum value 1.0. Thus, χ p ranges from 0 to 1.0, indicating from the none mixing to the fully mixing state.
To evaluate the mixing state of the whole field, the averaged mixedness (χ p ) is calculated according to Equation (11), Ω Ω χ p (x, y) dΩ (11) where Ω represents the whole field. Transferring information of discrete tracers into continuous number density (c a or c b ) is a key step. By dividing the simulation area into a limited number of statistical cells, we could calculate the number density of tracers in a cell. However, the cell should be selected with proper size in order to make an accurate demonstration of mixing (ds x and ds y are cell lengths on x and y directions respectively). Statistical cells with sizes of n times of the grid size (ds x = n · dx and ds y = n · dy) is employed, and n ranges from 1.0 to 8.0 is presented. The growth of χ p over time under different statistical cells is shown in Figure 7. All the curves show similar growing process of mixedness. Initially, the mixedness is nearly zero, while after a short period, the mixing accelerates and mixedness grows rapidly. Shown by results from different statistical cells that, the mixedness based on small cells (n = 1.0 for example) is relatively low during the whole growing process, and mixedness based on large cells (n = 8.0 for example) shows evident perturbation. A smooth growing curve can be achieved based on medium-sized statistical cells. In this case, we select n = 2.0 mainly for the following research.  Figure 8 shows the distribution of mixedness based on n = 2.0. Four representative time (t = 0.5, 1.5, 2.5, and 3.5) of RTI development are selected. It is found that the mixing mainly occurs at the contacting area of two fluids. As the interface being elongated and distorted in RTI flow, the contacting area increases, leading to further mixing.
Through integrating the mixedness (χ p (x, y)) along y direc- tion and then dividing the height of simulation region L y , the average mixedness in vertical direction ( χ p ) is, The evolution pattern of χ p based on n = 2.0 is shown in  Figure 9. Initially, χ p is almost zero and fluctuates along x axis. As RTI develops, two peaks of χ p appear along x axis. Comparing the appearing time and positions of the two mixedness peaks with RTI flow patterns in Figure 2, the two peaks coincide with the development of KHI, which occurs at the interface of the two fluids on both sides. The two peaks of the χ p gradually increase with time, indicating the enhancement of KHI over time. Finally, as χ p is almost even, KHI vortex disappears and the mixing is found on smaller scales.
The evolution of χ p inspires us that it provides a good indicator for KHI occurrence, with both the time and spatial range of KHI. Values of χ p also indicate the intensity of KHI during the development. χ p presents at least two advantages on capturing KHI. There is no need of visualizing the flow and manually differentiating the rotation region. The calculation of χ p is done separately from the main simulation, which means it is a post-process and can be done at any time with simulation data.
Additionally, the distribution of χ p based on two different statistical cells with size of n = 1.0 and n = 8.0 are provided by Figures 16(a) and 16(b) in appendix A. The two figures show similar results with Figure 9 and the details are discussed.

D. Compressibility effect on RTI mixing
Compressibility plays an important role in the growth of RTI as pointed by many researchers 13,69,70 . However, the effect of compressibility on RTI mixing has not been well analyzed. In this case, we intend to discuss how compressibility impacts RTI mixing. Firstly, a compression factor is defined through dimensional analysis 36,71 . The speed of sound (c s ) in ideal gas hypothesis is given as, where γ is the isentropic expansion factor. The Navier-Stokes equation coupled with ideal gas equation (p = ρRT ) then writes, where σ is the stress tensor expressed as σ = µ ∇u + (∇u) T − 2 3 (∇ · u) I . The sound speed c s and relaxation time τ are kept in the equation. The wavenumber k and the acceleration g are selected as basic units. Therefore, the typical length is 1/k, the typical time is (g · k) −1/2 , and the typical velocity is (g/k) 1/2 . Navier-Stokes equation is transferred to a dimensionless form, where the viscous term at the end of the equation writes, In the above equation, H c is the square of the ratio of the typical speed to the sound speed, defined as the compression factor, H τ is the ratio of relaxation time to typical time, It is noted that following the same dimensional analysis, BGK-Boltzmann equation can be transferred as, where H τ is in accord with the physical meaning of the Knudsen number, indicating the nonequilibrium of the flow system, so H τ is defined as Knudsen number.
In the following analysis on compressibility effect, the sound speed in upper fluid is chosen as c s . At the same time, in order to keep the nonequilibrium scale of different compressible systems unchange, the Knudsen number (H τ ) is kept constant as an initial condition, which is accomplished by adjusting different τ and g. H τ = 2.483 × 10 −4 and six different characteristic compression factors H c = 0.011596, 0.017393, 0.023191, 0.028989, 0.034787, and 0.040584, are selected. In the discussion, physical time is calculated according to a typical time t * = t/ [2π/ (g · k)] −1/2 , which is derived from initial perturbation on interface.
Firstly, the effect of compressibility on mixing is shown by Figure 10, which provides the average mixedness at different times. The mixedness increases with time under each compression factor in Figure 10(a). For a long-range from the initial time, curves of different compression factors show relatively close value, but with the increase of mixedness accelerates, curves of high compression show larger mixedness than low compression curves. Especially, the curve with the highest H c stands out at t = 8.0. In this case, high compressibility may enhance the mixedness at the late stage of RTI mixing.
However, the mixedness at the late stage is too large compared with initial values. Two ranges of t * at the early time are selected to help us better understand the mixing process, as shown by Figure 10(b) and (c). At the beginning process (see Figure 10(b)), there is an obvious decline of mixedness before t * = 0.5. According to previous analysis, the mixing starts from the interface, so the penetration of two fluids may be weakened at this stage. A tracer distribution pattern at t * = 0.5 shows that the local tracer density is small around the interface, indicating each fluid is retreating to avoid further contact at the interface. After t * = 0.5, due to the strong effect of acceleration, RTI continues to develop, but a declineascent pattern of mixedness is recorded.
Lai et al. 36 pointed that the effect of compressibility on the development of nonequilibrium in RTI is not monotonic destabilization 72 , but with two stages, which initially inhibits and later enhances the development. On RTI mixing, the effect of compressibility is also with two stages, as shown by the intersection of H max c and H min c curves in Figure 10(c). After the startup, the mixedness under high compressibility is lower than that of low compressibility, which is opposite to the relationship at the late stage. In this case, the curve of high compressibility catches up with the low compressibility curve at about t * = 3.9.
In order to understand the two stages of the effect of compressibility, Figure 11 demonstrates the tracer patterns of RTI flow, with three selected time points, t * = 2.0, 4.0, and 7.0, approximately located at the first stage, the transition stage, and the late stage of mixedness development; three compression factors, H c = 0.011596, 0.028989, and 0.040584, which are the minimum, a medium, and the maximum value in Figure  10. At t * = 2.0, see Figures 11(a1) to (a3), the low compressibility system shows spike growing faster as well as larger deformation. This difference in large structure (spike, bubble, and etc.) leads to the greater mixedness of high compressibility system. At the transition point, t * = 4.0, see Figures 11(b1) to (b3), although fluids penetrate each other, the flow pattern in low compressibility system is simpler, without highly expanded small structures. However, in high compressibility system, the spike is elongated and unstably deforms, leading to complex flow patterns. At the late stage (t * = 7.0), see Figures 11(c1) to (c3), distinct difference is found between three patterns from left to right, corresponding to low to high compressibility. Compared with a low compressibility system, the flow in a high compressibility system is turbulent, associated with many small drops and tiny vortices. Contrarily, there is no evident dissociation of fluids in a low compressibility system. In the late stage, turbulence enhances the mixing while a more disordered state corresponds to the greater mixedness.
Fitting different curves in Figure 10 with an exponential function, χ p = χ 0 + C 0 · exp (t * /t 0 ), where χ 0 , C 0 , and t 0 are fitting parameters. Mode of mixing in different compressibility systems is analyzed. In the fitting function, t 0 is the characteristic time scale of the mixing, and C 0 is a characteristic mixedness. Both parameters decline with the increase of compression factor H c , as shown by Figure 12. Higher compressibility would lead to a twofold effect, a higher mixing rate (1/t 0 ) and a lower characteristic mixedness (C 0 ), the two together govern the mixedness. In a short period of time after the start of RTI, mixedness is mainly controlled by C 0 ; thus, compressibility inhibits the mixing. As time increases, the exponential term (exp(t * /t 0 )) dominates the changes of mixedness, so compressibility enhances the mixing at the late stage.
Intuitively, fluids with low compressibility are "hard" and those with high compressibility are "soft". At the early stage, a "hard" one has relatively high penetrating ability and leads to the quick formation of large structures, which contributes to the mixedness. At the late stage, small structures tend to generate in the "soft" system, which then leads to high mixedness.

E. Viscous effect on RTI mixing
Due to in many cases RTI exists in an extreme physical condition with drastic changes 73 , the viscosity of fluid is not constant and would furtherly impact RTI development. Generally, viscosity is considered a stabilization factor on RTI at the start-up stage as well as the reacceleration stage, due to viscous force acts as resistance and inhibits the development of KHI in late-stage which then inhibits the RTI 38,74 . However, the role of viscosity is not simple because, after the reacceleration stage, viscosity could lead to complex flow patterns which contribute to the mixing of fluids 33 . Therefore, we intend to discuss the role viscosity plays in mixing.
Ahead of the discussion on viscous effect, it is noted that in the single-relaxation-time discrete Boltzmann model (SRT-DBM) with the ideal gas equation, the product of relaxation time τ and temperature T governs the viscosity, Due to the pressure is changeable during the RTI development, the initial pressure at the interface, p 0 = 1.0, is utilized to determine the viscosity. Different values of relaxation time τ are employed to alter the viscosity. Additionally, as an intrinsic physical parameter of fluid, the relaxation time τ links with Knudsen number H τ either, which indicates the nonequilibrium state of the fluid and has already been shown in equation (18). Figure 13 shows the average mixedness obtained from RTI Mixing under different viscosity, with µ ranges from 4.0 × 10 −5 to 12.0 × 10 −5 . In Figure 13(a), all curves are close to each other at the beginning, but at late time, curves of low viscosity show greater mixedness than that of high viscosity. The late-stage result is in good agreement with the com-mon understanding that viscosity acts as an inhibition effect on RTI development, which could be extended to RTI induced mixing here. When magnifying curves of average mixedness near the start of RTI flow, a similar decline-ascent relationship is found between mixedness and time, which is in accord with Figure 10(b). After this period, an intersection of different curves is found in Figure 13(b), indicating viscosity is also not a monotonic effect on mixing. From the result, at the initial stage, high viscosity enhances the mixing, but in the late stage, high viscosity inhibits the mixing. An interesting phenomenon at the late stage is that, for a specific time, the mixedness does not continue to decrease with the increase of viscosity, as shown by Figure 13(c). The relationship between mixedness and viscosity at three time points, t = 3.0, 3.5, and 4.0, are demonstrated in this figure. The mixedness slowly drops and tends to a constant value with an increment of viscosity. The saturation of mixedness with viscosity indicates that large viscosity can suppress the RTI mixing and the effect is diminishing with the further increase of viscosity. However, at different times, the saturation mixedness is different and mixedness still increases with time. on mixing at different stages. Three time points, t = 1.2, 2.2, and 3.0 are selected, which approximately locate at the initial stage, the transition stage, and the late state respectively. Three different viscosity, µ = 4 × 10 −5 , 7 × 10 −5 , and 12 × 10 −5 ) are employes in the delineation. At t = 1.2, according to Figure 14(a), the spike in the high viscosity system is integral and with more contact with the bubble in the vortex region, which lead to higher value of mixedness. In this stage, we could also attribute the mixing to large-scale development. At t = 2.2, the mixedness of low viscosity system starts to become larger than that of high viscosity system. We can observe the elongated structure of fluids in Figure 14(b1), but in Figure 14(b3), which shows high viscous fluids, the development is constrained. Significant differences are observed at t = 3.0 through Figure 14(c), when the mixedness of low viscosity system fully transcends that of high viscosity system. At the time, as more tiny flow structures break down from mainstream and drops dissociate from the main body in the low viscosity system, the mixedness of low viscosity system significantly increases. In contrast, in high viscosity system, large structures are maintained. The results of RTI mixing at the late stage are in accord with the conclusion made by Liang et al. 31 that relatively high Reynolds number (corresponding to low viscosity in this paper) leads to high growth rate of RTI and generation of dissociated tiny fluid structures. However, detailed results presented above can make compensation for this understanding, that initially, viscosity is a favorable effect for mixing.
A scale analysis based on optimal exponential fitting (χ p = χ 0 + C 0 · exp (t * /t 0 )) of average mixedness over viscosity is conducted, and Figure 15 shows the characteristic time scale (t 0 ) and characteristic mixedness (C 0 ). The curves of two parameters both can be separated into two parts, which is different from the compressibility. The time scale increases with viscosity initially and then becomes nearly constant with slight decrease. Characteristic mixedness shows similar relationship with viscosity. At late stage, the stabilization of two parameters leads to the saturation effect of mixedness. From the scale analysis, the effect of viscosity is not only attributed to the transition of dominant role between the time scale and characteristic mixedness, but also to the different stages of t 0 and C 0 respectively. The two-stage effect of viscosity is also caused by the separated development of large structure and small structure at two stages. A fluid system with high viscosity is considered "hard" and one with low viscosity is considered "soft". At the first stage, a "hard" one has a relatively higher penetrating ability and leads to higher mixedness, while a "soft" one has low penetrating ability as well as low mixedness. However, at the late stage, more tiny structures tend to generate in the "soft" system due to it is more delicate than the "hard" one, which leads to higher mixedness. Additionally, a "hard" enough system can successfully suppress the generation of tiny structures, which causes a saturation effect of viscosity on mixedness.

IV. Conclusions
In this paper, tracers are introduced into the DBM simulation of compressible RTI flow. The tracers enable observation of fine structures with clear interfaces in the late RTI mixing stage. In the velocity space of tracers, interesting distribution patterns are shown at different stages. These patterns bring quite dense information for the RTI research, which opens a new way for analyzing and accessing significantly deep insights into the flow system. With the aid of tracing the interface as well as interfacial physical quantities, various nonequilibrium behaviors on the interface of fluids are systematically investigated.
A mixedness χ p is defined by the local spatial distribution of tracers. With statistics, the distribution of vertically averaged mixedness χ p along the horizontal direction shows an interesting perspective to study the occurrence of KHI as well as evaluation of KHI intensity.
The effect of two important physical factors, compressibility and viscosity, are explored for RTI mixing. Both of them show two-stage effects. Specifically, high compressibility and low viscosity initially inhibits and later enhances the RTI mix-ing. The effect can be roughly understood via the following physical picture. In the initial period, the mixing is mainly induced by the evolution of large structures, while at the late period the mixing is mainly induced by the generation of small structures. Low compressibility and/or high viscosity make the system "hard" which favors the initial evolution of large structures. High compressibility and/or low viscosity make the system "soft" which favors the generation of small structures in the late stage.
The field averaged mixedness χ p increases nearly exponentially with time, which indicates the existence of a characteristic time scale t 0 . The time scale t 0 shows different dependences on compressibility and viscosity. It monotonically decreases with compressibility, while shows a two-stage dependence on viscosity. It firstly increases and then slightly decreases with viscosity, which means there is an appropriate viscosity for the system to get minimum mixedness. It is indicated that one can inhibit the mixedness via adjusting system viscosity. At the late stage, for a fixed time, the field averaged mixedness shows an exponential decrease with the viscosity, which indicates two points, (i) there exists a characteristic viscosity µ 0 and (ii) the mixedness will not show a meaningful decrease with viscosity after a certain value reaches.

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