Jet-induced jammed states of granular jet impacts

The impacts of granular jets for both frictional and frictionless grains in two dimensions are numerically investigated. A dense flow with a dead zone emerges during the impact. From our two-dimensional simulation, we evaluate the equations of state and the con- stitutive equations of the flow. The asymptotic divergences of pressure and shear stress similar to the situation near the jamming transition appear for the frictionless case, while their exponents are smaller than those of the sheared granular systems, and are close to the extrapolation from the kinetic theoretical regime. In a similar manner to the jam- ming for frictional grains, the critical density decreases as the friction constant of grains increases. For bi-disperse systems, the effective friction constant defined as the ratio of shear stress to normal stress, monotonically increases from near zero, as the strain rate increases. On the other hand, the effective friction constant has two metastable branches for mono-disperse systems because of the coexistence of a crystallized state and a liquid state.

The impacts of granular jets for both frictional and frictionless grains in two dimensions are numerically investigated. A dense flow with a dead zone emerges during the impact. From our two-dimensional simulation, we evaluate the equations of state and the constitutive equations of the flow. The asymptotic divergences of pressure and shear stress similar to the situation near the jamming transition appear for the frictionless case, while their exponents are smaller than those of the sheared granular systems, and are close to the extrapolation from the kinetic theoretical regime. In a similar manner to the jamming for frictional grains, the critical density decreases as the friction constant of grains increases. For bi-disperse systems, the effective friction constant defined as the ratio of shear stress to normal stress, monotonically increases from near zero, as the strain rate increases. On the other hand, the effective friction constant has two metastable branches for mono-disperse systems because of the coexistence of a crystallized state and a liquid state.

Introduction
Non-equilibrium phenomena induced by impacts have been extensively studied in various contexts, such as nuclear reactions [1][2][3], nanotechnology [4,5], water-bells [6,7] and granular flows [8][9][10][11][12][13][14][15][16][17][18]. Crater morphology is studied via an impact process of a free-falling water drop or a grain onto a granular layer [9,10], while a sinking grain produces a sand jet [11]. The impact of a granular jet on a target produces a sheet-like scattered pattern or a cone-like pattern, depending on the ratio of the target diameter and the jet diameter [8], which is also found in water-bell experiments with low surface tensions [6,7]. Cheng et al. suggested that the fluid state of a granular jet after an impact is similar to the Quark Gluon Plasma(QGP), which behaves as a perfect fluid through their experiment [8]. Recently, we reported that the shear viscosity during the impact is well described by the kinetic theory of the granular gas [19][20][21][22][23][24], though the small shear stress observed in the experiment is reproduced through our three-dimensional (3D) simulation [12,13]. Because the shear viscosity, at least, for 3D is not anomalous, the correspondence between a granular flow and QGP would be superficial.
To discuss the fluid state of granular jets, we need to know the details of rheology of moderate dense granular flows. A typical situation of the study for a dense granular flow is the flow on an inclined plane [25][26][27]. Bagnold proposed the constitutive equation for dense granular flows that the shear stress is proportional to the square of the shear rate [25], so called Bagnold's scaling, which has been verified experimentally [26] and numerically [27,28] under several conditions such as the flow down an inclined plane. Dense granular flows, however, have more variety of rheological constitutive equations for flows on inclined planes [29][30][31][32][33]. Conventional one would be the constitutive equation presented by Jop and coworkers [30], where the effective friction constant, defined as the ratio of shear stress to pressure, saturates from a static value at zero shear rate to a maximum value as the shear rate increases. The power-law friction law, which is also different from Bagnold's scaling for dense granular flows, is proposed via extensive simulations [34][35][36][37].
The aim of this paper is to investigate the rheological properties for two-dimensional (2D) granular jet impacts. Although some previous numerical studies on granular jets used 2D simulations to reproduce 3D experiments for the computational efficiency [14][15][16], it is unclear whether the rheological properties in 2D granular jets are qualitatively the same as those in 3D. Therefore, to clarify the qualitative difference between 2D and 3D granular jets is necessary. Because grains are easily packed through the impact in 2D, the system would be near the jammed state. Thus, we can investigate rheological properties of very dense granular fluids after the impact of granular jet flow, which cannot be achieved by 3D simulations and experiments. As a result, correlated flows appear in 2D granular jets, while uncorrelated flows characterized by the granular kinetic theory is realized in 3D jets. There are another advantage for the visualization to use 2D system even for experiments to know detailed properties of particles in granular jets, such as contact networks (force chains) and the effect of crystallization for mono-disperse case. We also stress that it is easy to perform 2D or one layer experiments for granular jets.
In this paper, we perform 2D simulations for the granular jet in terms of the discrete element method (DEM) [54]. This paper is complementary to the previous 2D DEM study [16], and hard core simulations supplemented by the simulation of a perfect fluid model [15]. Indeed, although Huang et al. reported that the relevant role of the contact stress in a 2D granular jet, they were not interested in the critical behavior of jammed grains induced by the jet. Guttenberg suggested that the friction constant does not play a significant role, at least, in the scattering angle [14], while the effects on the jammed state induced by jets have not been studied in his paper. This paper is organized as follows: After the introduction of our numerical model in Sec. 2, we analyze the profile of the local stress tensor, the area fraction and the granular temperature. We also discuss the rheology of the granular jets for the frictionless case in 2D to compare their behavior with the jamming transition for a bi-disperse frictionless case. The effect of the friction constant is discussed in Sec. 4. In Sec. 5, our numerical results for a mono-disperse case are shown and the paper is concluded in Sec. 6. In the Appendix A, we comment on the artificial burst-like flow in 2D, which appears in the case of large µ p for soft 2/20 grains. In the Appendix B, we discuss the effect of the inhomogeneity of the temperature to Balgnold's scaling in terms of the method of Green's function. Fig. 1 (i) A typical snapshot of the simulation for the frictionless case withφ 0 = 0.90. Blue particles and green particles denote grains with diameter 0.8d and d, respectively and red particles are wall-particles. (ii) The corresponding contact forces among grains are visualized as black colored arrows. The average coordination number Z ≃ 0.526 and 71.5% of particles are not in contact in the region 0 < x ≤ 10 and |y| < R tar .

Model
We adopt DEM to simulate the jet [54]. We mainly focus on bi-dispersed soft core particles of the diameter d and 0.8d with the same mass m to avoid the crystallization. When the particle i at the position r i and the particle j at r j are in contact, the normal force F n ij is given by ij ≡ k n (R i + R j − r ij ) and F (vis) ij ≡ −η n (g ij ·r ij ), where r ij ≡ |r i − r j | and g ij ≡ v i − v j with the velocity v i and the radius R i of the particle i. The tangential force is given by , where the sign function is defined to be sgn(x) = 1 for x ≥ 0 and sgn(x) = −1 for otherwise,F t ij ≡ k t δ t ij − η tδ t ij with the tangential overlap δ t ij and the tangential component of relative velocityδ t ij between i th and j th particles. We examine the value of µ p from µ p = 0.2 to 1.0. Here, we adopt parameters k n = 4.98 × 10 2 mu 2 0 /d 2 , η n = 2.88u 0 /d, with the incident velocity u 0 for the frictionless case and µ p = 0.2. The value µ p = 0.2 is close to the experimental value for nylon spheres [55]. We use k n = 1.99 × 10 3 mu 2 0 /d 2 , η n = 5.75u 0 /d for µ p = 0.4, and k n = 7.96 × 10 2 mu 2 0 /d 2 , η n = 10.15u 0 /d for µ p = 1.0. These sets of parameters imply that the duration times are, respectively, t c = 0.10d/u 0 for the frictionless case and µ p = 0.2, t c = 0.05d/u 0 for µ p = 0.4 and t c = 0.01d/u 0 for µ p = 1.0, the restitution coefficient for a normal impact is unchanged e = 0.75 for the frictionless case and for all µ p . The reason why we adopt these parameters for large µ p is that many overlaps among grains lead to the artificial burst-like flow, if we adopt the identical t c to frictionless case, as is shown in the 3/20 Appendix A. For the tangential parameters, we choose k t = 0.2k n , η t = 0.5η n . We adopt the second-order Adams-Bashforth method for the time integration of Newton's equation with the time interval ∆t = 0.02t c .
An initial configuration is generated as follows: We prepare a triangular lattice with distance between grains 1.1d and remove particles randomly to reach the desired density. We control the initial area fraction φ 0 /φ ini ≡φ 0 before the impact as 0.30 ≤φ 0 ≤ 0.90 with the initial area fraction before the removal φ ini = 0.612, 0.780 for the bi-disperse and the monodisperse case, respectively, and 8,000 particles are used. We average numerical data over the time 180.0 ≤ tu 0 /d < 300.0 after the impact. The initial granular temperature, which represents the fluctuation of particle's motion, is zero. The wall consists of particles in one layer with the same diameter d and the same mass m, which are connected to each other and with their own initial positions via the spring and the dashpot with the spring constant k w = 10.0mu 2 0 /d 2 and the dashpot constant η w = 5.0η n , respectively. A typical snapshot of our simulation and that of the contact force network are shown in Fig. 1 (i) and (ii), respectively. Blue, green and red particles denote grains with diameter 0.8d and d, and wall-particles, respectively in Fig. 1 (i) and all of the corresponding contact force network among grains are visualized as black colored arrows in Fig. 1 (ii). It is easily found that the contact force network emerges during the impact. It should be noted that the average coordination number Z ≡ i =j Θ(R i + R j − r ij )/N ≃ 0.526 and 71.5% of particles are not in contact in the region 0 < x ≤ 10 and |y| < R tar , where Θ(x) and N represent the Heaviside function and the number of particles in the region.
We evaluate physical quantities near the wall in two regions: 0 < x ≤ 5d and 5d < x ≤ 10d, where we call (a) and (b) layers in the followings, respectively. We use R jet /d = 15.0 and R tar /R jet = 2.2 with the jet radius R jet . We adopt the Cartesian coordinate, where y = 0 is chosen to be the jet axis, and divide the calculation region into the y direction y = −5∆y, −4∆y, · · · , 0, · · · , 5∆y, with ∆y ≡ R tar /5. Then we estimate physical quantities in the corresponding mesh region with k∆y < y < (k + 1)∆y (k = −5, −4, · · · , 4). Numerical data are averaged over ten initial configurations with the sameφ 0 and error bars in figures denote their variance.
We calculate the stress tensor as in Ref. [56]. The microscopic definition of the stress tensor at r is given by where i and j are indices of particles, µ, ν = x, y, the contact force between i th and j th particles F ij µ and i denotes the summation over the particles denoted by i located at r. A is the are of each mesh at r and u iµ (r) = v i µ −v µ (r) with the mean velocityv µ (r) in the mesh at r.

Rheology of Granular Jets for the frictionless case
In this section, our numerical results of granular jet, for 2D frictionless cases are presented. The results for frictional grains will be reported in Sec. 4

Existence of the dead zone and the profile of the are fraction
Chicago group suggested the existence of the dead zone near the target, where the motion of the grains is frozen, [14,15]. Ellowitz et al. suggested that the dead zone exists in the sense that the velocity of grains are close to zero in Ref. [15]. However, as is shown in our previous paper [13], although the velocity of grains at the center is small, the fluctuation of the velocity, i.e. the granular temperature T g , defined by T g ≡ i∈c mu 2 i /DN c with the number of grains N c in the mesh c and the spatial dimensions D, is the largest at the center in 3D (D = 3).
On the other hand, we verify the existence of the actual frozen layer (a) i.e. T g ≃ 0. The fluctuation of the grain velocity in (a) layer is suppressed, while the motion is not frozen in (b) layer for 2D granular jets (D = 2). The numerical data for T g in 2D for the frictionless case are shown in Fig. 2 (i). T g is the smallest at the center y ≃ 0 in (a) layer, which cannot be found in our previous 3D study (see Fig. 2 in Ref. [13]), while T g is the largest at y ≃ 0 in (b) layer. In very recent paper by Chicago group, it is suggested that the dead zone also exists in 3D experiment by introducing the effective temperature whose definition is not explicitly written [17]. 1 The profile of the packing fractions divided by φ max with φ max ≡ π/(2 √ 3) ≃ 0.907 in 2D are shown in Fig. 2 (ii) for the frictionless case. In 3D, the packing fraction divided by φ 3D max ≃ 0.740 ranges within 0.30 < φ/φ 3D max < 0.75. Compared with 3D, grains in 2D are well packed: 0.79 < φ/φ max < 0.94. Note that φ in (b) layer is almost independent of the position, while φ in (a) layer strongly depends on the position.

Profile of the stress tensor
The profiles of the stress tensor for (a) and (b) layers of frictionless grains are shown in Fig. 3 (i) and (ii), respectively. We stress that there exists a large normal stress difference between σ xx and σ yy in each layer as in 3D case [12].
The profiles of the stress tensor in (a) and (b) layer forφ 0 = 0.90 are shown in (i) and (ii), respectively. There exist large normal stress differences between σ xx and σ yy , and the shear stress is much smaller than the normal stress in (b) layer, though it is not so small in (a) layer.
Ellowitz et al. suggested that the profile of the velocity and the pressure for the granular jet are reproducible from the simulation of a perfect fluid [15] but our result may not support their claim. Indeed, the shear stress looks small but finite. Moreover, the large normal stress difference exists in both layers, which does not exist in the perfect fluid. We should note that they have not discussed the stress tensor itself in details, though they reproduce some similar feature through their hard core simulation. In addition, Huang et al. indicated the relevant role of the contact stress in their DEM simulation, which may be an indirect objection to the perfect fluidity of the jet flow [16].

6/20
We follow the analysis in Ref. [34]. Here, we introduce two dimensionless numbers consisting of pressure: I T ≡ T g /P d 2 and I s ≡ D xy m/P with pressure P ≡ (σ xx + σ yy )/2. We plot numerical data on φ vs I T plane and φ vs I s plane, in Fig. 4 (i) and (ii), respectively. Comparing φ in (a) with (b) layers against the identical I T , φ in (b) layer has a little larger value than φ in (a) layer at the same I T , while all φ against I s are collapsed on a universal curve (Fig. 4 (ii)) We can fit the data by the equations with constants φ T , a T , α T , φ s , a s and α s . Fitting parameters are determined simultaneously by using Levenberg-Marquardt algorithm [57]. The obtained equations of states are written as where the comparison of Eqs. (2) and (3) with numerical data for the frictionless case are shown in the main figure of Fig. 4 (i) and (ii), respectively. From Eqs. (2) and (3) which suggest the pressure diverging at φ T or φ s , the granular particles are well packed with the fraction sufficiently close to the jamming point.   (2) and (3). The insets denote numerical data for (i) log 10 I T vs log 10 |φ T − φ| and (ii) log 10 I s vs log 10 |φ s − φ| to examine how good the fitting results are.
The obtained parameters from our simulation are φ T = 0.858 ± 0.006, a T = 0.980 ± 0.1, α T = 1.15 ± 0.1, φ s = 0.834 ± 0.001, a s = 3.94 ± 0.5 and α s = 1.36 ± 0.05, where the error originates from the fitting. We also plot log 10 I T vs log 10 |φ T − φ| and 7/20 log 10 I s vs log 10 |φ s − φ| and the corresponding slope α T /2, α s /2 in the inset of Fig. 4 (i) and (ii), respectively, to examine how good our fitting results are, by using obtained critical densities φ T and φ s . Note that the conventional jamming point φ J ≃ 0.8425 at which the pressure diverges is located between φ s and φ T and close to φ s [51]. The asymptotic divergences of pressure for the frictionless case are described as In a conventional picture based on the extrapolation of the kinetic theory, the divergence of the pressure is expected to originate from the divergence of the radial distribution function i.e.
as φ → φ c , with the radial distribution function g(φ) the critical density φ c = 0.82 and the freezing density φ f = 0.69 [20]. In our case, the data are not far from P d 2 /T g ∼ P/mD 2 xy ∼ (φ c − φ) −1 expected from the conventional view based on the extrapolation of the kinetic theory, where T g ∼ md 2 D 2 xy is assumed. On the other hand, Hatano demonstrated an elegant scaling law in the vicinity of φ J , where the corresponding exponents are estimated as α s = 2.8 and α T = 1.7 from his data of the jamming transition [48]. Otsuki and Hayakawa showed the phenomenological explanation of the critical behavior near φ J and they predicted α s = 4.0 and α T = 2.0 [44,47]. It should be stressed that the critical scaling of the jamming transition is analyzed in the D xy → 0 limit, and the critical exponents strongly depend on the choice of the jamming point. Because the strain rate cannot be controlled in our setup, the jamming point is not clearly defined. Moreover, there are no data above the jamming transition in which the residual stress exists. Thus, our obtained exponents are smaller than those of the jamming transition for sheared granular particles. We note that the data for I T in (a) and (b) layers are separated, due to the difference of the profile of T g .
From Eqs. (4) and (5), T g and D xy are expected to satisfy The validity of Eq. (9) is verified in Fig. 5, which can be independent check of the scaling laws (4) and (5). From Fig. 5, Eq. (9) well reproduces the data for φ < φ s . Numerical data for φ ≃ 0.84 around the center, deviates from Eq. (9), which may result from the existence of the source point. The velocity field at the center is singular, compared with other regions. Actually, the similar deviation of the numerical data at the center from the theory can be found in our previous 3D study, in terms of the pressure and the shear viscosity [12]. The relation T g ∝ md 2 D 2 xy , which is equivalent to the Bagnold's scaling, is known to be derived from the energy balance equation for dense granular flow in the case that the heat flux can be negligible [28]. The dotted line in Fig. 5 represents the curve which can be derived 8/20 from Ref. [19] by taking the frictionless limit, which is written as Because there exists the unique critical density φ c for f σ (φ) and f Tg (φ), the conventional curve does not exhibit the critical behavior: as φ → φ c , while md 2 D xy /T g → 0 as φ → φ s in our setup, due to the two critical densities φ s < φ T . It should be noted that the functional form of f σ , f Although there exist the inhomogeneity of T g , as in the dead zone near the target. This is because we can generalize the discussion of Bagnold's scaling, at least, if the inhomogeneity is small (see Appendix B). Indeed, our case satisfies the condition that the gradient of the thermal velocity 2T g /m is much smaller than that of the velocity field. Because the momentum transfer plays major roles in the energy balance equation, where the only relevant time scale would be the shear rate. The detail analysis for the inhomogeneity of T g is shown in Appendix B.  Fig. 5 The numerical data for md 2 D xy /T g with Eq. (9) are compared for φ < φ s . The deviated data exist around φ ≃ 0.84 in (a) layer, which may result from the existence of the source point at y ≃ 0. We also plot f σ (φ)/f Tg (φ) as the dashed line for comparison. 9/20

Constitutive equation
3.4.1. Effective friction constant. Let us discuss I s dependence of effective friction constant µ * ≡ −σ xy /P to obtain the constitutive equation. Numerical data for the frictionless case are shown in Fig. 6. The behavior of µ * is conventionally described as where µ * starts from a static value of µ s at zero shear rate and converges to a limiting value of µ max at high I s . We obtain µ s = 0.0153 ± 0.009, µ max = 0.521 ± 0.06 and I 0 = 0.0820 ± 0.02 for the frictionless case by fitting. Thus, µ * (I s ) can be fitted by the conventional relation (14), which is denoted by a solid line in Fig. 6. It should be stressed that µ s is close to zero.  Fig. 6 Numerical data for the bi-dispersed frictionless grains are plotted on µ * vs I s plane. Red and blue points denote data for (a) and (b) layer for severalφ 0 , respectively. All points are fitted into the phenomenological equation in Eqs. (14) and (15), where we cannot judge which equation is better from the data.

Some researchers proposed a different constitutive equation called the power-law friction
which also well reproduces numerical data [32][33][34][35][36][37], where β ranges from 0.28 to 1.0, depending on the dimension, microscopic parameters and the friction constant of grains. Numerical data can be fitted by Eq. (15) within error bars, where we obtain b = 1.18 ± 0.1 and β = 0.592 ± 0.03, assuming µ s = 0 for the frictionless case. The fitting result of Eq. (15) is denoted by a dotted line in Fig. 6. As can be seen in Fig. 6, there is no significant difference between Eqs. (14) and (15). We, of course, cannot discuss the superiority of one of frictional laws from our simulation. We should stress that the fitting for both Eqs. (14) and (15) leads to almost zero µ s . This implies that the residual stress is negligible in the granular fluid after the jet impact.

3.4.2.
The asymptotic divergence of the shear stress. Let us discuss the asymptotic divergence of the shear stress: 10/20  Fig. 7 Equation (17), which denotes the divergence of the shear stress, is independently checked for the bi-disperse frictionless case. Red and blue points denote data for (a) and (b) layer for severalφ 0 , respectively. Equation (17) well reproduces numerical results. Because σ xy itself in (b) layer are small, the error bars in (b) layer are larger than those in (a) layer.
with an exponent β s . By using the divergence of the pressure (5) and the power-law friction µ * ∝ I β s , we obtain the constitutive equation for σ xy − σ xy mD 2 which is checked independently against the numerical data (Fig. 7). The exponent is estimated to be β s = (1 − β/2)α s ≃ 0.96. We also plot log 10 |φ s − φ| vs log 10 | − σ xy /mD 2 xy | with obtained φ s and the slope −β s in the inset of Fig. 7. The numerical data for large φ are deviated from the theoretical curve, due to the small shear stress and shear rate at y ≃ 0. Because we use the power-law friction β > 0, the divergence of shear stress may be slightly weaker than that of P . Here, Bagnold's scaling is still satisfied even in the vicinity of the "jammed" density, which is in contrast to the actual jamming transition [38][39][40][41][42][43][44][45][46][47][48][49][50][51].
Let us compare our observed critical behavior of the shear stress with the case of the jamming transitions, in details as well as that of the extrapolation of the kinetic theory. When we adopt the extrapolation of the kinetic theory, we have −σ xy /mD 2 xy is used and dimensionless shear viscosity η * ≡ η/η 0 is introduced with η 0 ≡ mT g /4πd 2 and shear viscosity η ≡ −σ xy /D xy . The extrapolation from the kinetic regime by Garcia-Rojo et al. predicts that σ xy diverges at density different from P d 2 /T g and β s = 1.0 [49]. Therefore, our analysis based on the power-law friction Eq. (15) predicts the results similar to Garcia-Rojo et al [49]: β s = 1.0 and φ s < φ T . On the other hand, the exponents for the divergence of −σ xy /mD 2 xy at the jamming transition for sheared granular materials are estimated to be β s ≃ 2.6 from data in Ref. [48] and β s = 4.0 in Ref. [44]. Thus, our corresponding exponent β s = 0.96 is much smaller than those of the jamming transition for sheared granular systems and rather close to the result of the kinetic theory. Otsuki et al. [44] studied the difference of soft core jamming and the asymptotic divergence of hard core systems. Then, they confermed the 11/20 exponent β s can only deviate from 1.0 in very narrow critical region, in which the soft core effect becomes relevant. Although our system has high density, the number of particles in contacts is still not large. Therefore, we may regard the granular fluid after the impact as a hard core fluid.

Effect of the friction constant
We examine how the fluid state depends on the friction constant from the simulation for µ p = 0.2, 0.4 and 1.0. It is noteworthy that the separation between (a) and (b) layer exists for larger µ p , even on φ vs I s plane. The results for µ p = 0.2 and µ p = 1.0 are shown in Fig.  8 (i) and (ii), respectively. Because there are two branches on φ vs I s plane, we adopt Eq. (3) to fit the data in (a) or (b) layer, separately. Figure 9 denotes the critical densities and the exponents for each µ p , where φ s in both (a) and (b) layer slightly decrease as µ p increases, and α s in (a) layer increases as µ p increases, while it decreases in (b) layer. We note that the decrease of our critical densities φ s both in (a) and (b) layer are gentler than φ L of the jamming transition for sheared granular systems [39].  Fig. 8 Fitting results of Eq. (3) for µ p = 0.2 (i) and µ p = 1.0 (ii). As µ p becomes larger, data for (a) and (b) layer deviates from each other.
In contrast, the friction law is little affected by the friction of grains. The results for the friction law are shown in Fig. 10 (i) for µ p = 0.2 and (ii) for µ p = 1.0, where the numerical data can be fitted by both Eqs. (14) and (15). We stress that µ * monotonically increases from near zero, as I s increases, even for large µ p .

Result for the mono-disperse case
Here, we discuss the impact of granular jets in 2D for the mono-disperse case. A typical snapshot zoomed near the target is shown in Fig. 11, where grains are crystalized near the wall. The black solid lines in Fig. 11 (i) are drawn by hand to clarify the grain boundary between the crystallized region and the disordered region, where the boundary becomes a slip line. We also visualize all of the corresponding contact force network in Fig. 11 (ii).
The notable difference of the mono-disperse cases from the bi-disperse cases appears in the friction law. We plot µ * (I s ) for the mono-disperse case of frictionless grains, µ p = 0.2 and µ p = 1.0 in Fig. 12 Fig. 9 The µ p dependence of the critical density φ s and the exponent α s in (i) and (ii), respectively. φ s in both (a) and (b) layer slightly decrease as µ p increases. The purple solid line denotes the corresponding critical density of jamming for frictional granular particles φ L [39]. The exponents in (a) layer increases as µ p increases, while they decrease in (b) layer.  Fig. 10 Numerical data for µ p = 0.2 (i) and µ p = 1.0 (ii) can be fitted into Eqs. (14) and (15) within error bars, where we cannot judge which equations are better. The friction law is little affected by the friction of grains. It should be noted that µ * monotonically increases from near zero, as the increment of I s , even for large µ p .
(b) layer cannot be fitted by a single curve, unlike the bi-disperse case. Judging from the snapshot (Fig. 11), grains, at least, in (a) layer are partially crystallized. Therefore, it is reasonable that the response of the crystallized region is different from that in disordered regions in (b) layer. The behavior of µ * (I s ) in (a) layer, which are observed in both frictional and frictionless cases, can be understood as follows. Because of the crystallization, a grain is trapped in a crystallized region. However, as I s increases, the grain can escape from the crystallized region. Thus, µ * decreases as I s increases.
The macroscopic friction µ * (I s ) for frictionless grains are different from that for frictional grains in (b) layer. The most remarkable difference between the frictionless and the frictional cases is the existence of peak of µ * at a small I s for the frictionless case, while there is no such 13/20 Fig. 11 Typical snapshot of the simulation for the frictionless and mono-disperse grains case withφ 0 = 0.90 near the target. Green particles denote mobile grains in the granular jet and red particles are wall-particles (i). The black solid lines are drawn by hand to clarify the grain boundary between the crystallized region and the disordered region, where the boundary becomes a slip line. Crystallization into a triangular lattice can be seen near the region enclosed by the black lines. All of the corresponding contact forces between grains are visualized as black colored arrows in (ii). a peak for frictional cases. Because a frictional grain can roll over grains, grains easily form a cluster. Therefore, the boundary between such clusters becomes a slip line. Thus, µ * (I s ) would be constant as I s becomes smaller. On the other hand, because a frictionless grain can neither roll over them nor slip, it is trapped in the crystallized region even for the large I s . Thus, µ * (I s ) for frictional and frictionless cases exhibit different behaviors in (b) layer. However, we should stress that there exist two metastable branches for both frictionless and frictional cases.

Discussion and Conclusion
We have performed two-dimensional simulations for the impact of a granular jet and discussed its rheology. We confirmed the existence of the dead zone, as is reported in Ref. [15] at least in (a) layer, unlike our previous three-dimensional cases [12,13]. There exists large normal stress difference, which has not been reported previously. The shear stress is much smaller than the normal stress, at least in (b) layer. We need to solve the inconsistency in (a) layer with Ref. [15].
We have analyzed the rheology of frictionless grains after the jet impact. We found that the pressure and the shear stress diverge with exponents similar to the extrapolations from the kinetic regime, and their exponents are smaller than those of the jamming transition for sheared granular systems. We adopted the power-law friction for µ * Eq. (15) to obtain the critical exponent for σ xy /mD 2 xy . The discrepancy between our case and the jamming transition for sheared granular systems would originates from (i) our system cannot reach the true jamming transition and (ii) the uncontrollability of D xy in our setup. The jamming point for a sheared system φ J is located between φ s < φ J < φ T and is close to φ s Our analysis 14 Fig. 12 Numerical data for the frictionless case (i), µ p = 0.2 (ii) and µ p = 1.0 (iii) are plotted on µ * vs I s plane. Red and blue points denote data for (a) and (b) layer for several φ 0 , respectively. Unlike bi-disperse cases, µ * (I s ) for (a) and (b) cannot be fitted into a single curve. µ * (I s ) in (i)-(iii) show similar behavior in (a) layer. However, in (b) layer, because frictionless grains cannot roll over the crystallized region, µ * (I s ) for the frictionless case shows similar dependence on I s to that in (a) layer, while the corresponding µ * (I s ) for frictional cases do not.
based on the power-law friction is consistent with that by Garcia-Rojo et al [49], where σ xy diverges at the density different from P d 2 /T g [49].
The effects of the friction of grains µ p have been discussed. Although Guttenberg [14] suggested that µ p does not play a significant role, at least, in the scattering angle via the approximate hard-sphere method [58], we found that the existence of the friction affects rheology of granular fluids after the impact. The separation between (a) and (b) layer appears for larger µ p , even on φ vs I s plane. The critical fraction φ s decreases as µ p increases, which is similar to the behavior of critical fraction of jammed frictional grains φ L . The corresponding exponent α s increases (decreases) as the increment of µ p in (a) layer ((b) layer).
The effective friction constant µ * (I s ) for the mono-disperse case has two branches because of the coexistence of the crystallized state and a liquid state. On the other hand, µ * (I s ) for the bi-disperse case can be described by known constitutive equations for dense granular flow [30][31][32][33][34][35][36][37].
Finally, let us comment on the rheological model proposed in a recent paper of Chicago group [17]. It is suggested that the granular fluid after the impact may be described by the plastic flow without the viscous stress and with the isotropic pressure. This suggestion is interesting, but our data may not support their suggestion. In fact, our data suggest the existence of viscous term (the shear stress depends on the location), the pressure is anisotropic, and no evidence of the existence of the residual stress as is shown in Fig. 6.
Appendix A. On artificial burst-like flows for the large µ p case In this appendix, we comment on the artificial burst-like flow in 2D, which appears in the case of large µ p with softer grains than those in the text. After the impact of a jet composed of softer grains with large µ p , the burst occurs when a grain slips, because large tangential force can be accumulated before the slip of a grain. In Fig. A1, we show the time evolutions of T g at −∆y < y < 0 in (a) layer for (i) the frictionless case and µ p = 0.2, and (ii) µ p = 1.0 with several stiffness, where T g for the frictionless case and µ p = 0.2 reaches the small steady values, while T g raise many times after the impact for µ p = 1.0 with large t c , due to the slip events. As t c becomes smaller, the burst-like flows are suppressed. Thus, we use harder grains for large µ p . Though there are a few small raises of T g for the frictionless case, they are out of our averaging time. Appendix B. Inhomogeneity of T g Here, let us discuss the effect of the inhomogeneity of T g to Bagnold's scaling. We demonstrate that T g ∼ md 2 D 2 xy may be valid because the gradient of 2T g /m is much smaller than that of the velocity field in our setup.
The numerical data for the profile of 2T g /m,v x and |v y | are shown in Fig. B1(i)(ii). Red empty (i) and blue filled (ii) points denote the data in (a) and (b) layer, respectively. The corresponding triangle, square and circle points are the data for 2T g /m,v x and |v y |, respectively in Fig. B1. Although 2T g /m is the smallest at the center y ≃ 0 in (a) layer, the inhomogeneity of T g is much smaller than that ofv x and |v y |. In particular, it is notable that |v y | linearly increases as |y| → R tar from zero.
The relation T g ∼ md 2 D 2 xy can be derived from the energy balance equation in the case that the heat flux can be negligible: where Γ(φ, T g ) = (1 − e 2 )Γ(φ)T respectively. The corresponding triangle, square and circle points are the data for 2T g /m, v x and |v y |, respectively. Although 2T g /m is the smallest at the center y ≃ 0 in (a) layer, the inhomogeneity of T g is smaller than that ofv x and |v y |. |v y | linearly increases as y → R tar from zero.
(B13) Here, Green's function, which satisfies the condition G(x = 0, y|x ′ , y ′ ) = 0 can be con- structed as G(x, y|x ′ , y ′ ) = G free (x, y|x ′ , y ′ ) − G free (x, y| − x ′ , y ′ ) (B14) We numerically calculate the integral in Eq. (B13) as δT (x i , y j ) ≃ k,l ∆x∆y d 2 G(x k , y l |x i , y j )T (x k , y l ) − l ∆yδT C (y l ) ∂G ∂x (x k , y l |x i , y j ) (B16) 18/20 with ∆x ≡ 5d = |x i − x i+1 |, ∆y = |y i − y i+1 | and compared with numerical data for δT (x, y) = T g (x, y) − T 0 in Fig. B2. Here, the value of T 0 is the average of T g (x = 3∆x, y) along y. The inhomogeneity of δT appeared in the dead zone rapidly decreases as x → ∞ both for simulation data (i) and Eq. (B13) (ii). We note that the characteristic length for the inhomogeneity is estimated as 1/kd ≃ 0.696 < ∆x. Interestingly, the analytic result in terms of Green's function well agrees with that of our simulation in the inhomogeneous region for x < 1/k. This result means that the heat flux from the boundary would not play an important role and the granular temperature is determined through the local shear rate. Furthermore, as discussed in Sec. 3.3, the numerical data for md 2 D xy /T g near (x, y) ∼ (0, 0) deviates from Eq. (9), due to the singularity of the source line. Thus, although there exist the inhomogeneity of the temperature, T g is locally determined through shear rate, i.e. T g ∼ md 2 D 2 xy is still valid, because the gradient of 2T g /m is much smaller than that of the velocity field in our setup.