{'fraction_non_alphanumeric': 0.04448688423299444, 'fraction_numerical': 0.013930236476998259, 'mean_word_length': 4.408932097827738, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large-scale structure of randomly jammed particles

We numerically analyse the density field of three-dimensional randomly jammed packings of monodisperse soft frictionles spherical particles, paying special attention to fluctuations occurring at large lengthscales. We study in detail the two-point static structure factor at low wavevectors in Fourier space. We also analyse the nature of the density field in real space by studying the large-distance behavior of the two-point pair correlation function, of density fluctuations in subsystems of increasing sizes, and of the direct correlation function. We show that such real space analysis can be greatly improved by introducing a coarse-grained density field to disentangle genuine large-scale correlations from purely local effects. Our results confirm that both Fourier and real space signatures of vanishing density fluctuations at large scale are absent, indicating that randomly jammed packings are not hyperuniform. In addition, we establish that the pair correlation function displays a surprisingly complex structure at large distances, which is however not compatible with the long-range negative correlation of hyperuniform systems but fully compatible with an analytic form for the structure factor. This implies that the direct correlation function is short-ranged, as we also demonstrate directly. Our results reveal that density fluctuations in jammed packings do not follow the behavior expected for random hyperuniform materials, but display instead a more complex behavior.

In particular, it was conjectured that randomly jammed packings of spherical particles are hyperuniform [1]. More precisely, the static structure factor S(k) of jammed assemblies of hard spheres has been fitted in the low-k regime to a non-analytic functional form [2,15] S( k) = A + B| k|, where A and B are numerical constants. Because A is numerically found to be very small, the hyperuniform property that S( k → 0) = 0 follows. The hyperuniformity of jammed packings was first observed in monodisperse system [2,16], then also in binary mixtures [17,18], and even in polydisperse systems [18][19][20]. However, recent simulations have challenged this finding [21][22][23], and provided evidence that randomly jammed packings are not strictly hyperuniform, in the sense that S( k → 0) > 0. In addition to the vanishing wavevector limit, the nonanalytic linear-k behavior in Eq. (1) has several direct consequences at large but finite lengthscale. First, it implies that the pair correlation function g(r) converges to its asymptotic limit from negative values and as a power law [2], g(r) − 1 ∝ −r −4 , a behavior which has never been observed directly, to our knowledge. A second consequence is that density fluctuations display anomalous behavior at large scale, so that the variance of the fluctuations of the number density in a subsystem of size R scales (in three dimensions) as ∆D(R) ∼ R −4 , instead of the weaker R −3 scaling expected from the central limit theorem [1,3]. Deviations from R −3 behavior were recently reported in experimental work [24]. A third consequence is that the direct correlation function c(r) should become long-ranged and exhibit a power law decay [1], c(r) ∼ r −2 . Numerical evidence in favour of this behavior was recently published [25].
To conclude that randomly jammed sphere packings are disordered hyperuniform materials, one should ideally establish that Eq. (1) together with all its direct consequences are consistently observed in a single system. In particular, the power law decay of the pair correlation function is difficult to observe, because it is masked by an additional oscillatory behavior with a wavelength corresponding approximately to the particle diameter σ, arising from purely local correlations. These oscillations tend to obscure the power law decay at large distances [2]. In the same vein, it was noted that these local correlations may also affect the scaling of the density fluctuations ∆D(R) to produce an apparent anomalous scaling that could be unrelated to hyperuniformity [21].
In this work we simultaneously measure S(k), g(r), c(r), and ∆D(R) for random packings of spherical particles. We use large-scale simulations in three dimensions to determine the nature of density fluctuations at large distances. In addition, we introduce a coarse-grained density field to successfully disentangle genuine largescale correlations of the density field from more local effects, and this allows us to accurately determine the large-scale behavior of g(r) and ∆D(R) . Our results reveal a surprisingly complex pattern of density fluctuations at large-scale in jammed packings, but the results differ from the hyperuniform behavior in Eq. (1) on two grounds. We find that S( k → 0) does not vanish and that the non-analytic linear wavevector dependence does not consistently account for our large body of numerical results.
The outline of the paper is as follows. We explain our protocol to produce jammed packings in Sec. II. We then study respectively the static structure factor, the pair correlation function, the local fluctuations of the density and the direct correlation functions in Secs. III, IV, V, and VI. We conclude in Sec. VII.

II. PREPARATION OF JAMMED PACKINGS
We focus on the same model as in our previous studies [22,26]. We consider a system of monodisperse spheres in three dimensions, where spheres interact through a pairwise harmonic potential [27] is the Heaviside step function, r ij is the distance between particles i and j and σ is the particle diameter. Throughout this work, we use σ and ǫ as the units of length and energy.
We generate configurations of randomly packed spheres by the following protocol. We first prepare a random configuration of N spheres in a cubic box of dimension L at ϕ = 0.8, where ϕ = πN σ 3 6L 3 is the packing fraction. The number of particles is N = 512000 unless otherwise noted (see Appendix A for a discussion of finite size effects). Then we use the FIRE algorithm to minimize the potential energy and to find the mechanical equilibrium state [28]. The algorithm is efficient enough that large system sizes can be studied. Next we gradually decrease the packing fraction in small steps and minimize the energy after each step to obtain the mechanical equilibrium state at the desired packing fraction. Our strategy is to find energy minima for configurations above jamming, and to observe how the structure of these packings changes as density is varied. This approach differs from, and is numerically much simpler than, studies focusing on packings produced precisely at jamming, for which very precise algorithms need to be devised [23,29].
We repeat these calculations from independent initial random configurations and obtain a large number of random packings. The number of independent configurations is 240 for the lowest density and 960 for the highest density addressed. We carefully tested the convergence against the number of samples for each measured quantity, as detailed in Appendix B.
We report various correlation functions which are obtained by averaging over results obtained for each independently produced configuration at a given density. Note that these different configurations are characterized by distinct values of the critical density for the jamming transition, thus the distance to the jamming density fluctuates from one sample to another [30]. However because our system is sufficiently large, the fluctuations of these distances are negligible for all the parameters studied in this work, thus the sample average at each density can be taken safely. Indeed, the standard deviation of the pressure, which gives an estimate of the fluctuations of these distances is less than 1% of the average pressure, even at the lowest density addressed.

III. STATIC STRUCTURE FACTOR
Let us start by focusing on the static structure factor in Fourier space. The number density fields in real and reciprocal spaces are respectively defined as [31] ρ where R i is the coordinate of particle i. Using these fields, we define the density-density correlation functions where · · · = 1 L 3 d 3 x is the translational average and ρ = N L 3 is the number density. We denote the spherical averages of these correlation functions as G(r) and S(k). As in the case of the radial distribution function and the static structure factors in liquid states [31], G(r) and S(k) are related through a Fourier transform.
Because the static structure factor in this limit is the standard deviation of the number of particles, this property means that density fluctuations decay with increasing the lengthscale of observation more rapidly than the prediction obtained from the central limit theorem. Thus, hyperuniformity implies the existence of some sort of long-range correlations in the density field. The static structure factor of jammed packings generated by fast compressions of hard spheres was studied numerically, and is reported to behave as in Eq. (1) with A = 6.1 × 10 −4 and B = 3.4 × 10 −3 [2]. The observation that A is very small has led to the claim that jammed packings are hyperuniform. It is also noted that the linear wavevector dependence in Eq. (1) means that S(k) becomes a non-analytic function at k = 0, a behavior which is never observed in ordinary liquid states at thermal equilibrium [31]. We first compute the static structure factor S(k) in the following way. We compute the density field ρ( k) on the lattice points of the form k = (n 1 , n 2 , n 3 )∆k, where the n i 's are integers and ∆k = 2π L , and calculate S( k) =   Fig. 1. The pressure values at these state points are 0.029, 0.017, 0.0059, and 0.0026, thus these states are almost equally separated from each other in a logarithmic scale for the pressure [30]. Irrespective of the density, S(k) displays nearly linear behavior in the wide range of wave numbers at k ≤ 3, as previously reported [2,18]. However, deviations from this linear behavior take place at k ≈ 0.3 and S(k) becomes nearly constant or even increases slightly with decreasing k even further. Such deviations were recently noted [21][22][23] and they are more clearly demonstrated by measuring the second derivative S ′′ (k). As shown in the inset of Fig. 1, S ′′ (k) = 0 only holds for k > 0.5 but clear and robust deviations appear at lower wavevectors. Our findings show that these deviations are not due to the numerical difficulty to produce packings precisely at the jamming transition [23].
We have confirmed that the existence of this upturn at k ≈ 0.3 is numerically robust by changing both the system size and the number of independent configurations used to perform the ensemble average, as shown in Appendices A and B. With decreasing the density towards jamming, S(k) slightly increases at large k and slightly decreases at low k, but this dependence saturates at ϕ = 0.66, which is far above jamming. This indicates that the jamming criticality, which affects many physical quantities, does not dramatically affect the low-k behavior of S(k). This conclusion was reached in previous work [21,22], and originates from the fact that the jamming criticality is sensitive to detailed features of particle interactions near contact, whereas we focus here on large-scale density fluctuations.
The numerical results for S(k) can not be fitted using the non-analytic functional form S(k) = A + B| k| at k < 0.3. In the bottom panel of Fig Instead, the numerical results can be well fitted into the empirical analytic function This function converges to the constant value S ≈ S 0 when kξ ≪ 1, and behaves as a linear function S ≈ S 0 ξk/π when kξ ≫ 1 and it is thus well-suited to describe the saturation observed in the data near k ≈ 0.3. We display these functions for the values n = 2, 4 and 6 where S 0 is fixed to the value S 0 = 0.016 and the remaining fitting parameters are ξ = 7.32, 8.26 and 8.38 for n = 2, 4 and 6 respectively, are again determined by a least square fitting. The convergence of S(k) to the constant at k → 0 is well captured by both functions and the sharpness of the upturn at k ≈ 0.3 is better captured by the values n = 4 and 6. Note that the parameter ξ gives an estimate of the wavelength corresponding to the wavenumber where the upturn occurs. In summary, S(k) shows strong deviations from the proposed non-analytic linear behavior at low-k, and ex-hibits a sharp crossover which is well captured by an empirical analytic function given by Eq. (5).

IV. PAIR CORRELATION FUNCTION IN REAL SPACE
If S(k) has the non-analytic form S(k) = A + B| k|, then the real space density correlation function G(r) has the interesting property that its asymptotic decay becomes [2] Namely, there appears a long-range "negative" density correlation which does not have any characteristic lengthscale and is instead algebraic. Thus it is interesting to study the large-r behavior of G(r) to confirm whether such a correlation exists in jammed packings, because it represents the direct counterpart to the linear k dependence of the static structure factor. However numerical analysis of this putative power law behavior is very difficult because G(r) also has a strong oscillatory behavior which is caused by the short-range correlations of the density field. Physically, this is because the wellseparated wavevector regimes in the Fourier domain get entangled by the inverse Fourier transform. In order to overcome this problem, we introduce a coarse-grained density field and study its correlation function.

A. Coarse-grained density field
The coarse-grained density field ψ( x) is defined as where f ( x) is a Gaussian function which acts as a lowfrequency filter, and 1/ √ δ is a lengthscale controlling the filter width. Then, the coarse-grained density-density correlation function is defined as and so Q( r) can be seen as a coarse-grained version of G( r), and both functions should become equivalent when δ → ∞. To look into the large-r behavior, Q( r) is more suitable than G( r), because the Gaussian filter in Eq. (7) should eliminate the oscillatory behavior due to shortrange correlations, if δ is well-chosen. This statement becomes evident when we consider the relation between Q(r) and S(k) where When δ is small, f ( k) becomes extremely small at largek, thus the integral in Eq. (9) is not influenced by the first diffraction peak of S( k) which provides the oscillations in real space. In other words, local correlations are filtered out by the Gaussian filter which leaves the low-k behavior intact. We use δ = 0.3 (the filter width is 1/ √ δ ≈ 1.83) unless otherwise noted, which gives f ( k) 2 ≈ 10 −36 for | k| = 2π, where the first peak of S(k) is located. Note that the use of the Gaussian filter Eq. (7) is similar in spirit to the use of "definition I" in Ref. [21] although the Gaussian filter is presumably more efficient to suppress the effect of short-range correlations. Because Q(r) converges to 1 in the large-r limit, it is convenient to define which corresponds to the coarse-grained version of the total correlation function H( r) used in liquid state theory [31]. Hereafter, we denote P (r) the spherically averaged P ( r). Note that because f ( x) is properly normalized, the integration of P ( r) over the whole space gives S(0)/ρ. Before looking into the numerical results for jammed packings, we consider two simple examples of P (r). The first is the case when the static structure factor is constant S(k) = S(0). This simplification can apply for simple liquids because S(k) converges rapidly to the constant (compressibility) at low-k [31]. In this case, we can compute the integral in Eq. (9) which gives Namely, P (r) is characterized by a single Gaussian peak at r = 0 with a width given by 1/ √ δ. The oscillations in G(r) are perfectly suppressed in P (r) by the Gaussian filter, as expected. From Eq. (11), we observe that the integration of P (r) up to r = 1/ √ δ is enough to estimate S(0), showing that density fluctuations at the local scale are sufficient to recover the macroscopic limit in simple liquids.
The second example is the non-analytic function S(k) = A + B| k|, as in Eq. (1). We insert this function into Eq. (9) and evaluate the integral numerically, using A and B values from the fit shown in Fig. 1. The obtained P (r) is plotted in Fig. 2. We observe that P (r) not only has the positive peak at r = 0, it also has a negative peak near r ≈ 5. The bottom panel of Fig. 2 shows that the asymptotic decay of P (r) for r → ∞ is characterized by the power-law as expected also for G(r). Thus, in contrast to our first example, P (r) is now long-ranged. To recover the density fluctuations in the macroscopic limit, we must now integrate P (r) up to infinity to take into account the negative power law correlation at longer lengthscales.  Fig. 1. (b) A zoom on the negative correlation near r ≈ 5 of (a). The inset illustrates the r −4 power law decay of P (r) from negative values as r → ∞.

B. Numerical results
We now move to the numerical measurement of P (r) in our jammed packings. To this end, we evaluate P (r) in the following way. We discretize the simulation box into a set of lattice points x = (n 1 , n 2 , n 3 )∆x, compute the density field ψ( x) on each lattice point, and then calculate the correlation function P ( r) = 1 ρ 2 ψ( x)ψ( x + r) −1, followed by translational and spherical averages to finally get P (r). It should be noted that P (r) does not converge to 0 as r → ∞ in a finite system with a fixed number of particles (the same remark applies to G(r) in equilibrium systems [32]), and thus special care should be exercised about finite size effects, as explained in Appendix A.
The obtained P (r) is plotted in Fig. 3 for various packing fractions. The first important observation from that plot is that the oscillatory behavior found for G(r) with a short wavelength approximately given by σ, is very effi- The measured coarse-grained density correlation function P (r) at ϕ = 0.7, 0.68, 0.66, 0.653. using three different representations to emphasize the peak at r = 0 (a), the negative dip near r ≈ 5 (b), and the positive correlation near r ≈ 10 (c).
ciently suppressed for P (r), which allows us to carefully discuss the large-distance behavior of the pair correlation function in real space. This validates our idea to introduce a coarse-grained density field.
For all densities, P (r) has the peak at r = 0 which expresses the density fluctuations at the microscopic length- scale. We also observe that P (r) has a negative peak at r ≈ 5. This negative peak is very similar to the one of P (r) obtained from S(k) = A + B| k| (compare with Fig. 2), which affects density fluctuations in the macroscopic limit. However, a new feature emerges at larger distances, since the negative peak in the measured P (r) is followed by a positive correlation peak near r ≈ 10. This behavior is in sharp contrast with the non-analytic case S(k) = A + B| k|, where density correlations remain negative with a power law decay. It is also clear that these basic features of P (r) do not depend on the chosen density, and in particular on the distance to the critical density of jamming. Thus, the negative correlation at r ≈ 5 and the appearance of the positive correlation at r ≈ 10 are not associated to the jamming criticality but are instead robust features of jammed packings.
We look at the large-r behavior of P (r) more closely, and plot the quantity r 4 P (r) in Fig. 4. As expected from the structure of P (r) in Fig. 3, we find that r 4 P (r) is first positive, then becomes negative near r ≈ 5, and finally positive again for larger distances. In particular, we find that the large distance behavior of P (r) is qualitatively distinct from the non-analytic behavior expected for hyperuniform materials, Eq. (12), and it is actually even more complicated than expected.
To understand the origin of the second positive peak near r ≈ 10, we also plot in Fig. 4 the quantity r 4 P (r) obtained from analytic models of S(k) with various n, from Eq. (5). Although the analytic model with n = 2 does not show the second positive peak, the analytic model with n ≥ 4 reproduces the peak very well. Recall that the value of n corresponds in Fourier space to the sharpness of the crossover behavior near k ≈ 0.3, it is then clear that the second positive peak in P (r) or r 4 P (r) is the direct consequence of this low-k feature in S(k). Indeed, the second positive peak is located at a lengthscale comparable to value of ξ extracted from fitting S(k), which corresponds to the wavenumber of the crossover. Moreover, the expression for S(k) from the analytic model has complex singularities (branch points) at (kπ)/ξ = (−1) 1/n and for n > 2 these singularities correspond to an exponentially-decaying oscillatory behaviour at large distances. Note that the position of the second positive peak depends on the parameter δ in the coarse-graining in general, but it becomes independent of δ when √ δ is sufficiently large compared to the wavenumber of the kink in S(k). We have checked that our choice δ = 0.3 satisfies this condition.
We finally note that the second positive peak of r 4 P (r) does not decay rapidly with distance and seems to persist even at very large-r, at least up to r ≈ 20. Because the behavior of r 4 P (r) at r → ∞ is related to a linear wavevector of S(k) as k → 0, a positive value for r 4 P (r) at large-r is related to the upturn of S(k) at very low k, see the inset in Fig. 1. To check whether r 4 P (r) finally goes to zero at r → ∞ would require studying an even larger system size, which is beyond the scope of the present work.

V. DENSITY FLUCTUATIONS IN SUBSYSTEMS
Density fluctuations in jammed particle systems are often discussed in terms of the number density fluctuations within a spherical cavity [1,3,21,24]. The physical reason is that the suppression of density fluctuations expected for hyperuniform materials at low wavevector should correspond to an anomalous behavior of density fluctuations measured in subsystems of increasing sizes in real space [1,3]. Since we found above that S(k) in jammed packings is incompatible with hyperuniform behavior, mathematical consistency requires that local density fluctuations should display normal scaling, but this expectation contrasts with earlier reports [2,24].
The number density in a spherical cavity is defined as where R is the radius of a spherical cavity centered at position x. The variance of the fluctuations of the density is then defined as This variance is directly related to S(k) as follows where J n (x) is the Bessel function of order n. From the central limit theorem, we expect that this variance is proportional to R −3 , namely that it scales as the inverse volume of the cavity when R → ∞. However, in a hyperuniform system, it is known that this variance becomes proportional to R −4 , namely that it is dominated by the surface term. For jammed packings, this variance was measured and was shown to follow the surface term, which was considered as an evidence that jammed particle systems are hyperuniform [1,3]. However, it was pointed out that this quantity is considerably affected by the surface term due to the slow decay of the Bessel function in Eq. (15), and that extremely large values of R are required to eliminate this effect [21]. Physically, the reason is again that large wavevectors in Eq. (15) contribute significantly to the final result even though they originate from purely local correlation effects. To detect hyperuniformity, one should instead focus on low wavevectors corresponding to density fluctuations at large distances. It is therefore pertinent to the analyse local fluctuations of the coarse-grained density ψ( x), because this quantity will not be influenced by local correlations.
Following the above definitions for the number density ρ( x), we now define the coarse-grained density in a spherical cavity as and we define the variance of this density as This variance is also directly related to S(k) as Compared to Eq. (15), the term f (k) 2 now appears in the integral over wavevectors and its effect is to reduce the influence of large wavevectors and thus of short-range correlations over the measurement of local fluctuations of the density. For a truly hyperuniform materials, both ∆D(R) and ∆Ψ(R) would have the same R −4 anomalous scaling behavior. Plugging the simulation results for S(k) at ϕ = 0.66 into Eq. (15) and (18) and numerically evaluating the integrals, we calculate ∆D(R) and ∆Ψ(R) and plot them in Fig. 5. In the numerical integrations, we linearly extrapolated the simulation results for S(k) to k → 0. Therefore, this calculation gives reliable results for R ≤ L/2, where L = 74 is the dimension of the simulation box (N = 512000), but the reliability at larger R of course depends on the reliability of this extrapolation. The fluctuations of the number density ∆D(R) show the (apparent) hyperuniform scaling R −4 for a wide range of cavity sizes from R = 1 to 100. On the other hand, the fluctuations of the coarse-grained density ∆Ψ(R) show a normal scaling with inverse volume R −3 over almost the entire range of cavity sizes. Both functions seem to become equivalent at very large cavity sizes only, above R ≈ 200.
This result means that the apparent hyperuniform behavior of ∆D(R) mainly results from short-range correlations of the number density appearing in the integral (15), and is thus unrelated to a suppression of density fluctuations occurring at large scale. The normal R −3 scaling found for the coarse-grained density fluctuations is consistent with the behavior of the structure factor S(k) which does not show signs of hyperuniformity.
We also find that the integral in Eq. (15) is strongly affected by the large-k region of S(k) and a numerical integration up to k = 30 is required to obtain a converged result for the quantity ∆D(R) . This observation confirms our conclusion and is also consistent with the results reported in Ref. [21]. This means that future exploration of hyperuniformity in particle systems based on local fluctuations of the density should employ a coarse-grained density field in order to more directly detect large-scale effects and to decrease the influence of short-range correlations.

VI. DIRECT CORRELATION FUNCTION
We finally focus on the direct correlation function C(r). This function is defined through the Ornstein-Zernike equation [31] where H( r) = G( r) − 1 is the total correlation function. The convolution integral contained in this equation simplifies in the reciprocal space, and this relation can be rewritten as where C(k) is the Fourier transform of C(r). We use Eq. (20) to convert our numerical data for S(k) into C(k), and we then calculate C(r) by performing an inverse Fourier transform. Both S(k) obtained directly from the simulations and from the various models of S(k) are used. Note that the Fourier transform from C(k) to C(r) requires a special care because C(k) has a k −1 behavior at large k which stems from the divergence of the first peak in G(r) near jamming [30,33], and is thus a purely local effect again. To control the convergence of this numerical calculation, we multiply C(k) with a Gaussian window function, w(k) = e −(k/kw) 2 , before performing the numerical Fourier transform to ensure its convergence. We checked that all values k w ≥ 5 give essentially the same result for C(r) at distances r > 1. Thus, we fix k w = 100 in the following. In Fig. 6, we show C(r) obtained from direct simulations and from the various models for S(k). This function displays a sharp dip near r = 1 which reflects the physics at contact, as also seen in hard spheres [25]. There are two features at large distances to be observed in the direct correlation C(r) obtained from the simulations. First, C(r) follows a power law scaling |C(r)| ∼ r −2 over a modest range of distances r ≈ 1 − 6, but it then exhibits a sudden drop near r ≈ 8 where it changes sign. This power-law scaling was related in earlier work [25] to the linear behavior of S(k), because S(k) ∼ k directly implies from Eq. (20) that C(k) ∼ k −1 , and thus C(r) ∼ r −2 .
To understand the observed behavior of C(r) in more detail, it is instructive to use our various models for S(k) to obtain "synthetic" C(r) functions, which can then be r -2 n n n FIG. 6. Direct correlation function C(r) obtained from numerical simulations at ϕ = 0.66 and calculated from various models for S(k) from Eqs. (1,5). The analytic models with n ≥ 4 reproduce the measured behavior very well.
compared to the simulation results. The non-analytic model S(k) = A + B| k| indeed gives the expected powerlaw scaling r −2 of C(r) up to very large distances, and it shows no sign of a sudden drop. On the other hand, the analytic model with n ≥ 4 almost perfectly reproduces both the power-law scaling at very short distances r < 6 and also the drop near r ≈ 8. The amplitude of these functions then goes rapidly to zero at large distance and shows no sign of a power law behavior (not shown).
These observations imply that the drop of C(r) originates again from the saturation of S(k) observed at very low wavector, k ≈ 0.3, in Fig. 1. Note in addition that the drop in C(r) takes place at r ≈ 8, which is almost identical to ξ, the lengthscale corresponding to the wavenumber of crossover in S(k), and is also close to the second positive peak in P (r) discussed in Fig. 3. Our conclusion is then that a hyperuniform material should exhibit a r −2 power law scaling of the direct correlation function up to very large distances, whereas a non-hyperuniform material with an analytic structure factor would instead display a strong drop of C(r) at a finite distance. The second behavior is the one that is most consistent with our numerical results.
Comparing our results to an earlier report in Ref. [25], we conclude that our simulation results for C(r) are almost identical to the ones obtained for compressed hard spheres. This suggests that S(k) for such hard sphere packings is presumably equivalent to our own data, and presumably contains also a crossover to a nonhyperuniform behavior at low wavevector. Thus, our conclusions on the large-scale structure of jammed configurations may also apply for the configurations analyzed in Ref. [25], despite the difference in the preparation protocols. A more quantitative comparison between the two sets of results would be interesting.

VII. CONCLUSION
We have studied the nature of spatial correlations of the density field at large distances in random packings of spheres.
Our results show that a mathematical description of the density-density correlations in terms of a random hyperuniform structure does not describe jammed packings very well. Instead, density fluctuations do not vanish asymptotically at large distances, and they do not display an anomalous lengthscale-dependence when measured in subsystems. The real space counterparts of the proposed non-analytic linear wavevector dependence of the structure factor are not consistent with our numerical observations. We conclude therefore that density fluctuations display a complex behavior at large scale, but this behaviour appears analytic and not hyperuniform.
It would be interesting to understand the complex structure of such random packings from analytic theory, for instance by extending previous work [34,35] to also describe large-scale structural properties. Another interesting question is whether our findings affect the relevance of random sphere packings from the viewpoint of material science [10][11][12][13][14], or whether such details do not affect the transport properties of such materials. In that case, the structure of jammed packings would still be of considerable theoretical interest, as it forms a puzzle related to the structure of amorphous solids at large lengthscales. This question has been rarely addressed in the field of glassy materials, where only structural features from local to medium range are more typically discussed.
Because we focus on the large-scale structure, it is very important to keep finite size effects under control. In Fig. 7 we show that the low wavevector crossover behavior observed near k ≈ 0.3 for N = 512000 particles remains at the same position when using a larger system with N = 4048000 particles. The overall behavior of S(k) including the upturn at k ≈ 0.3 is unchanged, and thus it does not originate from any finite size artefact. Note that the upturn at even lower wavectors seems even more pronounced in the larger system. Similar data were obtained in two dimensions in Ref. [21].
In Fig. 8 we show the system size dependence of the pair correlation of the coarse-grained density field P N (r). We use the notation P (r) only for the one in the infinite system size limit and P N (r) is used for the one measured with a finite N . In contrast to the case of S(k), P N (r) exhibits a system size dependence. The large-r limit of P N (r) of the N = 512000 system is clearly negative, while that of the N = 4048000 system is closer to zero. Such behavior of the radial distribution function in thermal systems is a well-known finite size effect. We briefly review the method commonly-used to correct for this effect [32]. In the grand canonical ensemble, the density-density correlation function for an open system ρ( x)ρ( x + r) can be written as ρ( x)ρ( x + r) = N p(N ) ρ( x)ρ( x + r) N (A1) where ρ( x)ρ( x + r) N is the density-density correlation function for a closed system with N particles, and p(N ) is the probability that the open system contains N particles [31]. When N is sufficiently large, p(N ) becomes a function with a sharp peak at the most probable N , which we denote by N * . Therefore, by expanding ρ( x)ρ( x + r) N around N = N * , and taking the summation over N , we obtain where we introduce G(r) = ρ −2 ρ( x)ρ( x + r) and G N (r) = ρ −2 ρ( x)ρ( x + r) N as in Eq. (3). The same derivation can be directly applied for the coarse-grained density-density correlation function Q(r) as defined in Eq. (8). The derivation may hold even in our athermal system if one assumes that p(N ) is still characterized by a sharp peak, as in thermal systems, which is a reasonable hypothesis. Thus one can expect where Q N (r) is the value of Q(r) for a closed system with N particles. Because we focus on the large-r region, we can safely set Q N * (r) = 1 in the second derivative. This leads to the finite size correction for P (r): where P N (r) = Q N (r) − 1. To check the reliability of this analysis, we plotted P N (r) + S(0) N for N = 512000 and 4048000 in the bottom panel of Fig. 8. The results obtained from different system sizes become nearly identical, thus the finite size effect is corrected. We use P (r) obtained in this way in this work.  9. Dependence on the number of independent samples of (a) S(k) at ϕ = 0.7 and (b) r 4 P (r) at ϕ = 0.66. These results show that the low-k behavior of S(k) and the large-r behavior of P (r) are statistically significant.

Appendix B: Number of samples
In Fig. 9 we show that our measured results for S(k) and for P (r) are obtained for a sufficiently large number of independent samples. The low-k region of S(k) and the large-r region of P (r) have relatively large sampleto-sample fluctuations, but they eventually converge very well when the number of samples is larger than 480. This test confirms the robustness of our data set with respect to the chosen number of independent samples.