{'fraction_non_alphanumeric': 0.03296866860162238, 'fraction_numerical': 0.008062315497010742, 'mean_word_length': 4.27505888510861, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Random Diffusion Model with Structure Corrections

The random diffusion model is a continuum model for a conserved scalar density field driven by diffusive dynamics where the bare diffusion coefficient is density dependent. We generalize the model from one with a sharp wavenumber cutoff to one with a more natural large-wavenumber cutoff. We investigate whether the features seen previously -- namely a slowing down of the system and the development of a prepeak in the dynamic structure factor at a wavenumber below the first structure peak -- survive in this model. A method for extracting information about a hidden prepeak in experimental data is presented.


I. INTRODUCTION
In this work we study a generalization of the random diffusion model (RDM) introduced by Mazenko in Ref. 1, henceforth referred to as RDMI. The RDM is a model for nonlinear diffusion in colloidal systems. It is formulated in terms of a conserved density field and, in its simplest form, can be chosen to have Gaussian static statistics. The model can be motivated as a continuum generalization of the facilitated spin models of glassy dynamics [2][3][4][5][6][7][8]. These models have a density dependent kinetic coefficient which leads to slowing down in dense systems. The RDM similarly has a density dependent bare diffusion coefficient. In RDMI, this bare diffusion coefficient led to a slowing down as the density was increased. The more realistic RDM presented here exhibits similar behavior.
There has been much speculation but relatively few solid results in establishing the existence of a modecoupling theory (MCT) ergodic-nonergodic (ENE) transition in the field-theoretic models of the liquid-glass transition. The RDM is a candidate for the simplest such model, however this point will not be a focus of this paper. In RDMI, the model was shown to undergo an ENE transition at one-loop order, but not at two-loop order. The RDM is part of a larger class of models with density dependent bare diffusion coefficients including those discussed by Dean [9], Kawasaki and Miyazima [10], and Miyasaki and Reichman [11].
In RDMI, the Fokker-Planck formalism was used to generate a self-consistent perturbation expansion for the memory function. This memory function could be inserted into the kinetic equation which describes the timeevolution of the dynamic correlation function, the physical observable of interest. In the simplest realization of this model, one has a single control parameter, g, which controls the density dependence. RDMI treated a coursegrained system in what was termed the structureless approximation; the static structure factor of the system was assumed constant up to a wavenumber cutoff Λ.
Solving the kinetic equation revealed the following: 1. As g increases, the system slows down.

2.
A new peak, termed the prepeak, develops at a wavenumber q 0 away from zero. The form of this peak can be fit to a Gaussian where the amplitude A decreases with time and the peak width 1/ √ B narrows with time. This form reveals a new growing length in the problem ( √ B) and fits of the the amplitude change smoothly from an exponential form at zero coupling to a power law form near a critical value. 3. Above this critical value, the system becomes unstable and the prepeak amplitude grows without bound. Such behavior is clearly unphysical and as such, the model breaks down.
In this paper, we generalize the random diffusion model to a more realistic form where wavenumber integrations are cut off naturally in the theory. This requires incorporating a more realistic static structure factor into the theory. Here, the integrals are naturally cut off by the large wavenumber decay of the direct correlation function. Our new model depends on two parameters; we have a coupling constant analogous to that seen in the structureless case, but our model also depends explicitly on density.
For this new model we find the following: 1. As the coupling and density increase, the system again slows down.
2. A prepeak forms, now located between q = 0 and the first structure peak. This peak and the peaks due to the static structure factor can be fit to Gaussians with decaying amplitude and narrowing width just as before and the amplitude again changes from an exponential to a power law form as the system approaches the critical crossover.
3. Above the critical values of density and coupling constant, we again find that the prepeak grows without bound. Whereas previously this separation between stable and unstable growth was characterized by a single number, our model has a separation given by a critical line in coupling constant-density parameter space.
This improved model therefore preserves the features seen in RDMI.
In this work, we again discuss in detail the slowing down of the system and the interesting feature of the dynamically generated prepeak. Though the development of this prepeak in both RDMI and this paper is surprising, it is worth pursuing. The random diffusion model is simple, yet incorporates nonlinearities known to be present and the perturbation is straightforward and without the "tricks" sometimes used to force things into mode-coupling form. It is dynamically generated instead of arising via the static structure factor and we find here that it is robust in the choice of said statics; what one may have worried was an artifact of a coarse approximation in RDMI is shown here to survive in the same form.
The final feature of the model, the unstable growth of the system above critical values of the model parameters, is a feature which reveals the breakdown of the model. In this region, the system likely wants to nucleate, but no terms are present in the model to provide this stabilization. We discuss some loose features of the instability and then propose the terms, motivated by density function theory, which may cutoff the growth.
Because our model uses a realistic structure factor, we can draw a closer connection to experimental results than was possible in RDMI. In the last section of the paper we present a method by which one can extract information about a possible prepeak hidden in experimental measurements of the dynamic structure factor. Though few, some experimental allusions to prepeaks exist and comparison is desired.

A. Introduction
Let us introduce the random diffusion model in the Fokker-Planck context. For a fundamental density field φ(x), we take an effective Hamiltonian quadratic in φ given by where δφ(x) = φ(x) − n and n = φ . We will study density time correlations through the equilibrium intermediate dynamic structure factor given by where φ(q) is the Fourier transform of the fundamental field δφ, averages are given by f (φ) = D(φ)W φ f (φ) with an equilibrium probability distribution of the usual form, W φ = e −βH φ /Z, andD φ is the adjoint Fokker-Planck operator defined below. We will also make frequent use of the static (equal-time) correlation function given bỹ which is related to the static susceptibility appearing in Eq.
(2) byC The adjoint Fokker-Planck operator takes the form where Γ φ is a transport matrix which incorporates the bare diffusion coefficient [12]. In RDMI, Γ φ was given by where the bare diffusion coefficient D(φ) was For this work, however, we choose a more general form By relaxing the constraining delta-function and introducing the general functions f 0 and f 1 , we will have considerably more freedom in regulating integrals which result from application of perturbation theory. For now, we leave the functions undetermined, but will later set them in such a way as to give the short-time sum rules correctly and give the large-wavenumber dependence of the vertex in the memory function as in mode-coupling theory. Our choice for the transport matrix (both the form of RDMI and the more general form here) is motivated as a course grained alternative to the microscopic form D(φ) = D 0 most often used. Density dependent terms are expected to play a role in the dynamics and as such, have been incorporated here.
Our model therefore depends on the functions f 0 , f 1 and χ and the constants D 0 and D 1 . The effects of extending the transport matrix to higher order through additional constants (D 2 , etc.) and functions (f 2 , etc.) will be addressed in a future paper.
The physical diffusion coefficient is given by where D(q) comes from the Fourier transform of the average of the transport matrix: Thus, we have where we have used the fact that averages over odd powers of δφ vanish.

B. Memory Function Formalism
In order to study the time evolution of the dynamic structure factor C(q, t), we organize our theory using the memory function technique [14]. In Laplace transformed space, the kinetic equation is given by where and where K(q, t) = K (s) (q) + K (d) (q, t) is the memory function decomposed into a static piece and a dynamic piece. Inverse Laplace transforming this equation gives the form ∂C(q, t) ∂t = iK (s) (q)C(q, t) First, the static memory function is determined (without approximation) by the equilibrium average Using our transport matrix given by Eq. (9), we have The dynamic part of the memory function is more complicated to compute and can be determined via a perturbation expansion. If we assume that the nonlinearities are small, we may expand in powers of D 1 , the coefficient of the density correction to the diffusion coefficient. A full derivation is given in RDMI, but at lowest order, the dynamic memory function is given by where V (q, k) is the vertex function. We define our vertex in light of our new choice for Γ φ as which has a simple dependence on the function f 1 . We choose this form to mimic the traditional mode-coupling theory form. (For more on traditional MCT, see, e.g. Refs. 15 and 16.) Now that we have defined this vertex, it is set to arbitrarily high order.

C. Structure corrections
Perturbation expansions inevitably run into wavenumber integrals such as the one we derived for the dynamic memory function, Eq. (18), which are potentially divergent. Care must be taken when forming and evaluating the integrals, and one usually must either develop a rationale for implementing a large-wavenumber cutoff or structure the vertex such that the integrals remain finite.
Constructing a vertex which enforces convergence can be difficult. For example, Ref. [17] studies the Dean-Kawasaki model [9,18] applicable to colloids. This work analyzes the theory in terms of time-reversal symmetry, but this approach ultimately leads to a vertex that is not of a standard MCT form and to divergent integrals. They propose that the solution may be a cumbersome resummation of higher-order diagrams to renormalize the vertex, but leave a full solution to future work.
For the random diffusion model, divergent integrals first appeared in RDMI. The approach in that work was to take a sharp cutoff at finite wavenumber Λ and (for most of the paper) to take the static susceptibility to be constant such that χ −1 (q) = r. This approach was called the structureless approximation because it corresponded to a coarse-grained system where the short-distance degrees of freedom (including the first peak in the static structure factor) had been integrated out [13]. While this simplified the problem and kept the integrals finite, one has reason to worry whether the results are strongly influenced by the nature of the cutoff and the structureless nature of the structure factor. It is therefore desirable to to impose a different method to regulate the integrals which does not raise such concerns.
To regulate the wavenumber integrals in this work, we now define our functions χ, f 0 , f 1 introduced earlier. First, let us restore the structure to the model by defining χ(q) in terms of the full static structure factor, S(q), as where and C D (q) is the direct correlation function which decays to zero at large wavenumber. For simplicity, we will use the solution to the Percus-Yevick approximation for hard-spheres in this paper. (See e.g. Refs. 19 and 20 for details.) In principle, one could take any realistic analytic approximation for S(q) or use experimental results.
Next, the function f 0 (q) appears in the static memory function K (s) (q) and the physical diffusion coefficient D p (q). Since we have no regularization constraints on either of these functions, let us take f 0 (q) = 1 which gives and Finally, the function f 1 (q) appears only in the vertex (and therefore in the dynamic memory function K (d) (q, t)). As discussed above, we want to choose f 1 (q) so that the memory function integral remains finite as we take the cutoff to infinity. Let us reason our way through the appropriate choice.
As it stands, the long wavenumber behavior of the vertex (beside the explicit q 2 dependence) is governed by χ −1 (q) which approaches unity at as q → ∞. Let us instead let the vertex go as the direct correlation function C D (q) which approaches zero in the same limit. Thus, we want f 1 (q) ∼ χ(q)C D (q). Normalizing f 1 (q) to unity at q = 0, we therefore have Note that if we now put f 1 (q) back into the bare diffusion coefficient (Eq. (9)) and take the large-q (short distance) limit, the density dependent D 1 term goes to zero returning the microscopic result D(φ) =D.
This now gives the vertex (25) This is the traditional MCT form which decays to zero as either q or k → ∞. Such behavior, we will see, is sufficient to keep the memory function finite without imposing a finite integral cutoff.
To simplify the kinetic equation let us now introduce dimensionless variables. The dimensionless wavenumber, time and density (packing fraction) are given by and η = (π/6)nσ 3 (28) where σ is the hard sphere diameter and the dimensionless correlation function is given by After inserting the static and kinetic equations in terms of dimensionless variables, we have where we've defined the simplified memory function and the coupling constant Let us make a few comments on this form. First, the kinetic equation now depends only on two parameters, the density η and the coupling constant R. Note that the temperature dependence has been absorbed into the dimensionless variables. The kinetic equation can be solved numerically to get f (Q, T ) and the general behavior as η and R are varied can be studied.
Second, the coupling constant R can be either positive or negative (as the coefficient D 1 can be either positive or negative). However, since R enters only through R 2 , the overall sign of D 1 is irrelevant. We choose for our purposes the constraint 0 ≤ R ≤ 1. One can imagine combinations involving negative, but large values of D 1 n such that |R| > 1, but since our perturbation expansion assumed small D 1 , we are already overly generous by studying R up to unity.
Finally, note that our form for ∂f (Q, T )/∂T diverges as η −1 as η → 0. While this may seem strange, this is just a manifestation of the characteristic time scale of the system slowing down with density. One can choose to rescale the dimensionless time as T → T ′ = T π/6η and the explicit divergence disappears. For various reasons, we have chosen to use the unscaled time, though if one were interested in low density systems, the scaled variable would be a more appropriate choice.

A. RDMI: A Review
Before we solve our model, let us review in more detail the results from RDMI. In that model, the kinetic equation (at lowest order) is given by and whereC is a constant and Λ is the wavenumber cutoff.
In the non-interacting case (g = 0), the equation can be solved analytically to get the simple form As the coupling g increases, the decay at large wavenumber values slows and a small peak develops. This peak (termed the prepeak) can be fit to a Gaussian of the form where the amplitude A has the time dependence while B satisfies where Q 0 , A 0 , B 0 , T A , T B , α, β and E are positive, time-independent fit parameters. A plot of f (Q, T ) as determined in RDMI for a coupling g > g * . Each line is equally spaced in time and the times are earliest to latest as the peak amplitude goes from smallest to largest. Wavenumber space is normalized such that the cutoff corresponds with Q=1. f (Q, T )=0 for larger Q.
As the system approaches the critical coupling g * , the peak amplitude goes to a completely algebraic decay with time such that E → 0. Above g * , the above analytic equations cease to hold; the peak no longer decays at long times, but instead grows without bound, rendering the system unstable. Examples of the behavior of this model can be seen in Fig. 1.

B. General Results
For our model, the evolution of the dynamic correlation function is given by Eq. (30). While similar to the form of RDMI, the incorporation of a realistic static structure factor leads to a slightly more complicated form. We can numerically solve the kinetic equation for a particular set of the parameters η and R. (The details of this solution are given in Appendix A.) The upper bound on the value of η is 0.74 which is the packing fraction of a close-packed solid. The upper limit for R is 1 as discussed at the end of the last section.
Before looking at the complete problem, let us first examine the simpler R = 0 case to see some of the features which are common to all the solutions of the kinetic equation. In this limit there is no dynamic feedback from the memory function and the kinetic equation can be solved analytically to give Thus, the correlation function decays exponentially with time and its wavenumber dependence is strongly determined by the static structure factor S(Q). This is the well known de Gennes narrowing form if we take a scaled time T ′ = T π/6η. In Fig. 2 we plot the behavior of f (Q, T ) for various values of η. We see that a number of peaks form which decay to zero with time. Beside the peak centered at Q=0, these peaks are centered at the same wavenumbers as the peaks of the static structure factor. Even with R = 0, we see that we have considerable structure in the theory.
We now turn to the case of R = 0. At small η, say η = 0.10 or 0.20, the behavior is nearly indistinguishable from the R = 0 case even as R is increased through all allowed values. As we increase to intermediate values of η, say η = 0.30, the situation becomes more interesting. In Fig. 3(a) we plot results for η = 0.30 and R = 0.51 and see a slight slowing down of the decay. That is, the peaks in the R = 0.51 case persist longer before decaying to zero than those in the R = 0 case.
If we stay at these intermediate densities and ratchet up to extremely large R, we see an interesting new feature; between the Q=0 peak and the first structure peak, a new small peak appears. With time, this peak also decays away to zero and disappears before the first structure peak. Because of its position in Q-space, we will call this peak the prepeak because it is located at a smaller Q value than the first structure peak. As we push the coupling up to the limit R = 1, we find that the prepeak persists for longer and longer times, but still decays to zero with all the other peaks. (Fig. 3(b).) Looking at the behavior of the second and higher structure peaks, we find little of interest. These higher order peaks mirror the behavior of the first structure peaktheir decay slows with increasing R -and we find no new features (e.g., higher order prepeaks). This remains true for most of the work presented in this paper and we will therefore restrict ourselves to discussion of the behavior at wavenumbers around and below the first structure peak unless otherwise noted.
Moving to η = 0.40, we again slowly increase the value of R. Initially, we clearly see the first structure peak and note that it decays to zero with time as in the R = 0 case. We plot f(Q,T) for R = 0.52 in Fig. 4(a) as an example. Increasing R a bit further to R = 0.56 however, we cross over to a new regime. Now, the first structure peak amplitude initially decreases, then slows down and finally begins to increase. (Fig. 4(b).) If we continue the solution to long times, we find that this growth has no bound and in fact accelerates. The prepeak (which is initially only weakly visible) grows faster than the first structure peak and the two peaks are of comparable amplitude only after the model and the numerical solution break down. We plot the amplitude of the prepeak and first structure peak in Fig. 5, including late times where the solution has become unphysical.
We may continue in this manner and look at higher values of η. We see a very general pattern appear. As R is increased from zero, there is initially a period where all peaks decay to zero with time. The plots look qualitatively like the R = 0 plots, however as R is increased, the peaks decay more slowly. We call this region of parameter space the stable region. Continuing, we pass a critical  value of R and see that the amplitudes of the peaks no longer decay to zero, but turn around and grow without bound. We call this (interesting, yet unphysical) region the unstable region and label the critical coupling R * . We find that as we increase η, the value of R * at which we cross from the stable to unstable regime decreases. As we continue increasing the value of R past R * , we find that the peak turns around its growth at earlier and earlier times.
The behavior of the prepeak is quite curious. As we increase η, we find that the prepeak becomes more difficult to see in the plots. Recall that at η = 0.30, we could clearly see the prepeak even though we were in the stable region. At high values of η we do not see the prepeak in the stable region and often do not see the prepeak in the unstable region until near the numerical breakdown when both the prepeak and first structure amplitudes become comparable. When both peaks are visible, we note that the first structure peak turns around and grows at an earlier time than the prepeak, but that the prepeak grows at a faster rate.
Though the prepeak is sometimes not visible from the raw data, we can extract some information about the prepeak's behavior. We believe that the behavior of the prepeak and the first structure peak, while quantitatively different, are related and that we do not lose information by monitoring just the latter. We discuss this in more detail in the sections which follow.

A. Stable-unstable crossover
We have seen that our solution has two distinct regimes: the stable region and the unstable region. We now wish to quantify the crossover from stable to unstable behavior. Because this crossover depends on both the coupling constant and the density, we have a critical function R * (η).
To determine this line, we follow the approach developed in RDMI. Let us first hold η fixed and vary R. If we are in the stable region, the amplitudes of the peak heights decrease monotonically with time. In the unstable region however, we may identify a time T min where the amplitude of a particular peak reaches its lowest point. (See, for example, Fig. 5.) As we approach R * from larger R values, we find that the value of T min approaches infinity as we get arbitrarily close to R * . We may therefore find T min for values of R close to but larger than R * and fit these values to a function diverging at R * : R * (η) which divides the stable and unstable regimes [21].
Our results for R * (η) are plotted in Fig. 6.
With a technique in place, one now may ask which peak to use in determining T min . We saw in the previous section that the prepeak and the first structure peak do not reach their minima at the same time and that their unstable growth is not equally paced. It is therefore unclear if the results from fitting one set of data will match the results from fitting the other. It is also unclear whether there may in fact be a region where we have unstable growth in one peak with stable decay in the other. A careful analysis, however, reveals that this is not the case. Values for R * determined from the prepeak amplitude closely match those determined from first structure peak amplitude and we are unable to find any region where the stability/instablilty of one peak does not match that of the other. Thus, we chose to use the first structure peak amplitude data for determining the critical function since it is the most clearly visible in all plots.
We can note several things from the plot in Fig. 6. First, we again see that for small to moderate values of η the solution is stable for all studied values of R. The critical line does not cross the R = 1 mark until some value between η = 0.30 and η = 0.35. Second, even though a transition to the unstable region at these moderate densities seems possible, one should be cautious about such a conclusion. Recall that we expect R to be a small perturbation parameter. By this reasoning, we may wish to restrict ourselves to discussion of such a transition at higher densities where R * is indeed small. Such higher densities are are those more typical of very dense liquids, solids and glasses.

B. Analytical Fits
We now see that the dynamic structure factor has the general form where f 0 (Q, T ) is the initial peak centered at Q = 0, f pp (Q, T ) is the prepeak and f n (Q, T ) are the peaks due to the structure factor, of which we are only concerned in the first (n = 1). Each peak can be modeled as a Gaussian such that where x may be 0, pp or 1 depending on the peak in question.
The peak centered at Q = 0 is the simplest to treat. For small Q, we may neglect the interaction term and we are left with the de Gennes form given earlier by Eq. (40): Comparing this to the generic Gaussian form of Eq. (43), we have A 0 = 1, Q 0 = 0 and A fit to the data shows excellent agreement. We find, for example, that for η = 0.55, R = 0.14, Eq. (45) predicts B 0 = 102.38T while a fit yields B 0 = 102.26T . Near the prepeak and first structure peak, the situation is more complicated. In the stable regime, each Gaussian peak has an amplitude which fits nicely to where the fit parameters vary with the coupling constant R. In the R → 0 limit, we know from Eq. (40) that the first structure peak decays exponentially (A 0 = 1, α = 0, E = T π/6ηS(Q 1 )) and that the prepeak is non-existant (A 0 = 0). In the R → R * limit, we find that both peaks tend toward pure power law decays (E → 0). The peak width can be nicely fit to where β is approximately 1 below the critical coupling, but decreases very rapidly as R → R * . The peak centers Q x drift slightly at early times, but quickly reach a fixed asymptotic value. Example fits for A and B can be seen in Fig. 7. In the unstable regime, peak amplitudes initially decrease then increase without bound and the peak widths narrow. At early times and away from T min , the amplitudes may again be fit by Eq. (46) and the widths by Eq. (47). However, the subsequent transition and growth pieces do not easily lend themselves to a universal fit, especially at late times when the growth accelerates and the numerical solution breaks down.
We find that these Gaussian peaks (in both the stable and unstable regimes) are tending toward delta functions as T → ∞; the peak widths are narrowing to zero and the amplitudes are time-dependent. This is the same behavior which was seen with the prepeak in RDMI.

C. Unstable Regime
The behavior of the model in the unstable regime is clearly unphysical. What is happening? At first, one might suspect that the system is undergoing a liquid to glass transition. However, as was shown in RDMI, an ergodic non-ergodic transition can be found for this model at lowest order, but is not sustained at higher order. It is therefore more likely that the system is instead simply freezing from a liquid to a solid. Since the model is for a system in equilibrium, as the density and coupling are increased, the system may wish to nucleate. Our model, however, is not equipped to handle such a transition and so the peaks grow and are not properly stabilized.
One solution would be to extend the model by choosing an effective Hamiltonian known to support freezing. According to density functional theory, we can choose the effective Hamiltonian to be which, after expanding the logarithm and rearranging, generates higher order terms compared to the Hamiltonian used in this work: and the higher order coefficients are given by v = −1/6βn 2 and u = 1/12βn 3 . A full treatment of this model with the above static perturbations will be given elsewhere.

D. Prepeak Behavior
We have seen that the prepeak is not always visible near and beyond the transition from stable to unstable. Let us explore the behavior and significance of the prepeak in greater depth.

Extraction of Prepeak Information
The prepeak is not always easily identifiable from the plots, but is hidden by other features in the data (e.g., the tails of the Q=0 peak and/or the first structure peak) and is of relatively small amplitude. To disentangle the prepeak from the other behavior, recall that we know the R = 0 solution exactly (Eq. (40)) and that the solution includes no prepeak and has a purely exponential decay for the structure factor peaks.
Let us divide out the R = 0 contribution to f (Q, T ), Though f (Q, T ) depends on R through an integrodifferential equation and cannot be separated into purely R = 0 and R = 0 contributions in either an additive or multiplicative fashion, the method has the potential to yield some information since we have seen that the late where the prime simply designates that while the form for the peak is the same, the amplitudes and widths differ [22]. As a test of the above proposed methods, let us look at the plot of f (Q, T ) for η = 0.45, R = 0.42 ( Fig. 8(a)). We see that we are in the unstable region since the first structure peak is growing with time. However, we see no sign of any prepeak. Let us divide out the η = 0.45, R = 0 values from our original f (Q, T ) to get h(Q, T ) ( Fig. 8(b)). Now we find a very clear prepeak at Q ≈ 3.40 which grows very rapidly in amplitude and narrows in width with time along with a smaller peak in the location of the first structure peak which also grows and narrows.
If we repeat this technique with other combinations of η and R we find that regardless of the parameters, a pre-peak centered near Q ≈ 3.40 can always be found. Also, whereas the prepeak and first structure peak centers exhibit some early time transient drift when the raw results are plotted, these new plots show fixed peak centers.
Let us now try to fit these peaks to Gaussians. The peak amplitudes can be fit to where α is positive for the first structure peak and negative for the prepeak. The width is given by where β is approximately 1 for both peaks.

Connection to Experimental and Simulation Results
While a prepeak developed in both the structureless approximation of RDMI and our own more realistic calculations here, there is little experimental evidence for the existence of such a peak in hard-sphere liquids. For molecular systems, a number of experiments reveal prepeaks below the first structure peak [23][24][25], but this work is likely unrelated to our model since the peaks can generally be explained through physical mechanisms unique to molecular systems (e.g. hydrogen bonding leading to clustering, molecular reconfigurations, etc.) and since these prepeaks generate static structure prepeaks instead of the dynamic structure we see in our model. Work with molecular dynamics simulations on the other hand, has shown some dynamically generated prepeaks [26][27][28][29], but again this work is on molecular systems and can be explained by uniquely molecular dyanamics, e.g. translation-rotation coupling.
While no evidence for prepeaks in hard-sphere or monatomic systems has been reported, we show here that at higher densities the prepeak was not readily visible in the dynamic structure factor except in the (non-physical) long-time breakdown. It is therefore conceivable that a prepeak similar to that found here may in fact exist in experimental results, but has simply not been detected and extracted.
How could a prepeak hidden in experimental data be extracted? In the previous section, we showed that hidden prepeaks appear if one plots the raw data divided by the theoretically determined zero-coupling data as in Eq. (51). This form, however, is not useful to an experimentalist, so let us rewrite it in terms of density n, temperature β and the physical diffusion coefficient, D p =D -parameters which would be known or could be determined -and the static and dynamic structure factors (S(q) and C(q, t) respectively) -which would be measured. We then have For such a method to yield useful information, extremely good data resolution for the dynamic structure factor at small wavenumbers in the prepeak region would be needed, as would small uncertainties on the quantities going into Eq. (55). Poor resolution and large uncertainties would easily wash out the small values that one is trying to extract.

V. CONCLUSIONS
We began with the random diffusion model introduced in RDMI which describes a system undergoing diffusive dynamics with a density-dependent bare diffusion coefficient. While the model had previously been studied in the structureless approximation where the short distance degrees of freedom had been integrated out, we wished to study the model from a more realistic view that included such features. By relaxing the delta function constraint of the transport matrix and introducing the smoothing functions f 0 (q) and f 1 (q), we were able to introduce a realistic structure factor back into the problem while still keeping the wavenumber integrals finite. The resulting kinetic equation then depended on two variable parameters, the dimensionless density η and the coupling constant R.
In our investigations, we came to the following conclusions: 1. When the coupling constant is increased, there is a significant slowing down of the system; exponential decay gives way to algebraic decay.
2. For large enough coupling there is a transition where the system goes from stable to unstable. All peaks grow without bound and the solution becomes unphysical. The critical coupling required to reach this transition decreases as the density increases.
3. It is possible this instability may be cut off by including terms in the model which allow for nucleation. We propose terms motivated by density functional theory and the results of this model extension will be pursued in another work.
4. Near but below the transition, a new peak, termed the prepeak, sometimes can be seen between the Q = 0 peak and the first structure peak. When it can be seen, its amplitude is less than that of the first structure peak. Above the transition, the prepeak always appears, though it may only be seen at late times when the peak amplitudes are growing without bound.
5. Near but below the transition, both the prepeak and first structure peak narrow to delta functions, but with algebraically decaying amplitude.
6. It is possible to isolate the prepeak by dividing out the R = 0 solution. In such a case we find the prepeak and structure peaks to be well separated with growing amplitudes and narrowing widths. This technique can separate out the prepeak even when it cannot be seen in the raw data and we find that the prepeak is always centered near Q ≈ 3.40 regardless of coupling or density.
7. By the above technique, one should be able to take an experimentally determined dynamic structure factor and isolate the prepeak using only experimentally measured or controlled parameters. Such an isolation might reveal previously undetected prepeak information and give indirect measurements of parameters of the system not usually accessible such as the the coefficients of the terms in the bare diffusion coefficient.
These results confirm that the features seen in RDMI were not due to the nature of the structureless approximation, but were indeed inherent features of the model. In addition, the model is now sufficiently realistic to invite comparison with experimental results. While the prepeak seen in our model has not been reported experimentally, we find that at reasonably accessible parameter values, it is often hidden; we provide a method by which a hidden prepeak might be extracted. The model, however, still lacks a sufficient mechanism to allow for freezing. Future work introducing static perturbations to the Hamiltonian may rectify this. where If we ignore the intricacies of g(Q, T ) for the moment, we can advance our solution to f (Q, T ) one time step at a time (in the normal Euler fashion) as using the initial conditions that f (Q i , T 0 ) = 1. The function g(Q, T ) includes a time integral over the memory kernel which can be decomposed into a similar differential equation, and can thus similarly be solved. Note now that since the time steps are not even, interpolation may be needed to find values of functions off the lattice. This will be addressed in detail below. At the next nested level, the memory function is itself an integral over wavenumber. In three dimensions, dK = K 2 dKd sin θdφ. Integrating over φ and expressing dot products as |Q · K| = Q 2 + K 2 − 2QK cos θ reduces the integral to two dimensions, which are again solved by the Euler method using interpolations.
Thus, to summarize, as we advance f (Q, T ) from one time step to the next, we must first advance the memory function array by one time step, then perform the time integration over it to get the function g(Q, T ) and finally use g(Q, T ) to fully advance f (Q, T ).

Adaptive Step Routine
We desire an adaptive time step size since we expect our solutions to sometimes offer smooth decays to zero and other times offer rapid growth and eventual instability. Let us describe the routine.
When choosing step sizes, there is always a tradeoff between maximizing computation speed (larger step sizes) and minimizing error (smaller step sizes). To quantify this trade off, let us imagine advancing a function y(x) from initial value y(x 0 ) up to y(x f ) in two different ways. In the first method, let us take two steps of size h 1 = (x f − x 0 )/2 and call the result Y 1 . In the second method, let us instead take one large step such that h 2 = 2h 1 and call that result Y 2 . The difference between these two results is then △Y = |Y 1 − Y 2 |. If one now specifies a desired difference between these values △Y 0 -which in some sense quantifies the trade off between accuracy and speed described above -then we can infer that the best choice for the next step is given by We see that if the actual difference is smaller than our desired difference, we can afford to increase the step size, whereas if the actual difference is larger, we decrease our step size. In our numerical integration, we specify an initial time step and tolerance △Y 0 . After each time step, the program determines where in Q-space the fastest growth/decay is occurring (where one assumes △Y will be the largest) and then computes the two test values described above at that point; the program integrates up from f (Q, T j−1 ) to f (Q, T j+1 ) first with two steps of the current size and then again with one step at twice that. From these points, the next step is determined.
Two limiting actions are taken to prevent run-away changes from occurring. First, a safety factor of 0.9 is added to the right hand side of Eq. (A5). This serves to underestimate the new step size so that the adjustments always err on the side of smaller error rather than faster speed. Second, the new step size is never allowed to be more than four times the previous step size. This prevents dramatic changes in step size.
The choices of these limiting values (0.9 and four times) as well as our choice for the tolerance △Y 0 = 10 −6 were chosen through rough trial and error.

Interpolation Algorithm
Integrations over the memory kernel require sampling the function at times not on the previously computed Q i , T j lattice. In such cases, a two-dimensional bilinear interpolation was used. If one wants to know the value of f (Q, T ) and knows the value at surrounding points (Q 1 , T 1 ), (Q 1 , T 2 ), (Q 2 , T 1 ) and (Q 2 , T 2 ) where Q 1 < Q < Q 2 and T 1 < T < T 2 , then we can estimate the value of f (Q, T ) as and u ≡ (T − T 1 )/(T 2 − T 1 ).