Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

HTML conversions sometimes display errors due to content that did not convert correctly from the source. This paper uses the following packages that are not yet supported by the HTML conversion tool. Feedback on these issues are not necessary; they are known and are being worked on.

  • failed: xpatch
  • failed: xpatch

Authors: achieve the best HTML results from your LaTeX submissions by following these best practices.

License: arXiv.org perpetual non-exclusive license
arXiv:2312.12875v1 [cond-mat.dis-nn] 20 Dec 2023

, .

Dynamic Screening by Plasticity in Amorphous Solids

H. George E. Hentschel Dept. of Physics, Emory University, Atlanta Ga. 30322    Anna Pomyalov Dept. of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel    Itamar Procaccia Dept. of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel Sino-Europe Complex Science Center, School of Mathematics, North University of China, Shanxi, Taiyuan 030051, China    Oran Szachter Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem, Israel 9190
(December 20, 2023)
Abstract

In recent work it was shown that elasticity theory can break down in amorphous solids subjected to nonuniform static loads. The elastic fields are screened by geometric dipoles; these stem from gradients of the quadrupole field associated with plastic responses. Here we study the dynamical responses induced by oscillatory loads. The required modification to classical elasticity is described. Exact solutions for the displacement field in circular geometry are presented, demonstrating that dipole screening results in essential departures from the expected predictions of classical elasticity theory. Numerical simulations are conducted to validate the theoretical predictions and to delineate their range of validity.

I Introduction

The response of amorphous solids to oscillatory strain is a well studied subject, mainly in protocols employing simple shear, cf. for example [1, 2, 3]. While important and indicating a lot of interesting physics, systems under oscillatory shear do not succumb easily to analytic scrutiny, making the comparison of analytic predictions to experiments or simulations quite difficult to accomplish. In this work we turn to oscillatory forcing by another protocol, of an oscillatory inflation of an inner circular boundary, creating non-uniform oscillatory strain on an amorphous solids that is contained between this and a much larger outer circular boundary that is maintained stationary.

The reason for this choice is two-fold. First, in a series of recent works it was shown that classical elasticity theory needs to be reconsidered and amended to describe correctly the mechanical response of amorphous solids subjected to time-independent non-uniform strains [4, 5, 6, 7, 8, 9, 10]. A new phase is formed in which the elastic fields are being screened by emergent geometric dipoles in the resulting displacement field. The transition from a quadrupolar screening in the solid phase to a dipolar screening is reminiscent of the structural transition in 2D crystals from hexagonal solids to hexatics, except that here the transition is in the emergent displacement field and the structure of the solid is always amorphous [9]. Our theory provides a classical field theory for describing the mechanical state of a deformed amorphous solids, and it predicts anomalous behavior that is observed in both numerical simulations and experiments. In this paper we extend the theory to a new direction, to describe mechanical responses of amorphous solids to dynamical loading.

The second reason for choosing the present geometry is that it allows considerable analytic progress. We consider the effects of inertia and plastic responses in amorphous matter. For the sake of obtaining analytic solutions we will focus on an amorphous solid contained in an annulus of outer radius routsubscript𝑟outr_{\rm out}italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT and inner radius rinsubscript𝑟inr_{\rm in}italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT, such that the inner radius oscillates with a fixed frequency ΩΩ\Omegaroman_Ω. Initially, before oscillations begin, the material will be brought to mechanical equilibrium in which the resultant force on each particle vanishes. Once oscillations start, we will be interested in the time-dependent displacement field, evaluated with respect to the initial equilibrated positions. We will show that classical elasticity fails to predict correctly this displacement field, and that an appropriate theory requires taking into account the screening introduced by quadrupolar and dipolar charges that form due to plastic responses.

The structure of the paper is as follows: in Sect. II we review the dynamics as expected for our configuration from the solution of the equations of motion dictated by classical elasticity. We then solve analytically for the displacement field of a purely elastic medium which is subjected to oscillatory inflation of an inner boundary. We find that the dynamics is already non-trivial, exhibiting interesting features. Section III introduces quadruple and then dipole screening, leading to predictions of a rich array of expected solutions in which new length scales emerge spontaneously, breaking down elasticity theory. In Sect. IV we describe simulation results to test the prediction of the theory. The comparison of theory to simulations calls for some careful considerations. First, it is important that the dynamics will describe oscillations around a well defined mechanical equilibrium state. Second, we need to deal with the issue of dissipation. In the numerics we introduce dissipation by damping terms in the particle collisions. In the theory there are dissipative terms proportional the rate of change of displacement fields. This requires careful discussion. Finally, nonlinear effects should be kept at bay. When all these are considered, we find indeed dynamical responses that are in accord with the novel theory. Section V offers conclusion and a discussion of the road ahead.

II Dynamics of amorphous solids in two dimensions

II.1 The displacement field in purely elastic solids

We start by reviewing the classical approach to dynamics and dissipation á la Landau & Lifshitz [11], and then we develop the basic ideas that are called for to accommodate the physics of amorphous solids.

In a purely elastic medium of mass density ρ𝜌\rhoitalic_ρ, denote the stress field by 𝝈𝝈{\bm{\sigma}}bold_italic_σ, and the displacement field by 𝒅𝒅{\bm{d}}bold_italic_d. The equation of equilibrium is given by

𝝈=0.𝝈0\nabla\cdot{\bm{\sigma}}=0\ .∇ ⋅ bold_italic_σ = 0 . (1)

i.e. the force per unit volume is zero in equilibrium . Once we are not in mechanical equilibrium the force does not vanish but is equal to the acceleration times the mass per unit volume

ρ𝒅¨=𝝈.𝜌¨𝒅𝝈\rho\ddot{{\bm{d}}}=\nabla\cdot{\bm{\sigma}}\ .italic_ρ over¨ start_ARG bold_italic_d end_ARG = ∇ ⋅ bold_italic_σ . (2)

For an isotropic body, one can write 𝝈𝝈\nabla\cdot{\bm{\sigma}}∇ ⋅ bold_italic_σ in terms of the displacement to get [11]

ρ𝐝¨=μ2𝐝+(μ+λ)(𝐝),𝜌¨𝐝𝜇superscript2𝐝𝜇𝜆𝐝\rho\ddot{\mathbf{d}}=\mu\nabla^{2}\mathbf{d}+(\mu+\lambda)\mathbf{\nabla}% \left(\nabla\cdot\mathbf{d}\right)\ ,italic_ρ over¨ start_ARG bold_d end_ARG = italic_μ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_d + ( italic_μ + italic_λ ) ∇ ( ∇ ⋅ bold_d ) , (3)

where μ𝜇\muitalic_μ and λ𝜆\lambdaitalic_λ are Lamé coefficients.

II.1.1 Dynamical response to oscillations of an inner circle

Elasticity theory holds equally well in two and in three dimensions. The theory described below can be easily extended to three dimension, and see for example [8]. In this paper the simulations are performed in two dimensions for computational efficiency, and therefore we specialize the equations to the case of an oscillating inner circle of initial radius r0subscript𝑟0r_{0}italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and time dependent radius rin(t)subscript𝑟in𝑡r_{\rm in}(t)italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t ), in a system bounded by a fixed, rigid outer circular boundary of radius routsubscript𝑟outr_{\rm out}italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT. The simulations were done according to the following protocol. At the beginning of each periodic driving cycle, the inner circle was inflated to rin(t=0)=r0+Δsubscript𝑟in𝑡0subscript𝑟0Δr_{\rm in}(t=0)=r_{0}+\Deltaitalic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t = 0 ) = italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + roman_Δ, and the system was equilibrated. Then we periodically inflate and deflate the central disk according to

rin(t)=rin(t=0)δsin(Ωt)subscript𝑟in𝑡subscript𝑟in𝑡0𝛿Ω𝑡r_{\rm in}(t)=r_{\rm in}(t=0)-\delta\sin(\Omega t)\,italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t ) = italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t = 0 ) - italic_δ roman_sin ( roman_Ω italic_t ) (4)

where Δ,δ,ΩΔ𝛿Ω\Delta,\delta,\Omegaroman_Δ , italic_δ , roman_Ω are parameters. After a given number of oscillations, the system is again allowed to equilibrate. Then the procedure is repeated for ten times. The reason for this initial set up is to guarantee that the subsequent oscillations are around a well defined equilibrated state and the resulting dynamics is reproducible. It is important that the oscillating displacement field is measured with respect to a well-defined equilibrated state. In order to compare the numerical responses to the offered theory, we calculate and analyze the angle-averaged radial displacement.

Clearly, in a realistic amorphous solid such an oscillatory inflation exerts work, and without dissipation the system will not reach a stationary state. Thus in the numerical simulations presented below, cf. Sect. IV, we will add a small dissipative term to the dynamics of the disks that comprise the amorphous solids. This results in the system reaching an oscillatory steady state. At this point we continue to focus on an ideal elastic medium to solve Eq. (3) as it stands, without adding dissipation that in an ideal elastic solid does not exist.

Refer to caption
Figure 1: Examples of the function fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) according to the purely elastic theory as shown in Eq. (13), for different values of the driving frequency ω=1.5,5,7.5𝜔1.557.5\omega=1.5,5,7.5italic_ω = 1.5 , 5 , 7.5 ( dot-dashed green, dashed red and solid blue lines, respectively).

In the presence of circular symmetry we can focus on the radial component of the displacement field, i.e

dr(r,t)𝒅(r,t)𝒓/r.subscript𝑑𝑟𝑟𝑡𝒅𝑟𝑡𝒓𝑟d_{r}(r,t)\equiv{\bm{d}}(r,t)\cdot{\bm{r}}/r\ .italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r , italic_t ) ≡ bold_italic_d ( italic_r , italic_t ) ⋅ bold_italic_r / italic_r . (5)

The equation for drsubscript𝑑𝑟d_{r}italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT becomes

ρd¨r=2μ+λr2[r2dr′′+rdrdr].𝜌subscript¨𝑑𝑟2𝜇𝜆superscript𝑟2delimited-[]superscript𝑟2superscriptsubscript𝑑𝑟′′𝑟superscriptsubscript𝑑𝑟subscript𝑑𝑟\rho\ddot{d}_{r}=\frac{2\mu+\lambda}{r^{2}}[r^{2}d_{r}^{{}^{\prime\prime}}+rd_% {r}^{\prime}-d_{r}]\ .italic_ρ over¨ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT = divide start_ARG 2 italic_μ + italic_λ end_ARG start_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT start_POSTSUPERSCRIPT start_FLOATSUPERSCRIPT ′ ′ end_FLOATSUPERSCRIPT end_POSTSUPERSCRIPT + italic_r italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT - italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ] . (6)

This equation of motion is solved with the boundary conditions

dr(rin,t)=rin(t)rin(t=0)=δsinΩt,subscript𝑑𝑟subscript𝑟in𝑡subscript𝑟in𝑡subscript𝑟in𝑡0𝛿Ω𝑡\displaystyle d_{r}(r_{\rm in},t)=r_{\rm in}(t)-r_{\rm in}(t=0)=-\delta\sin% \Omega t\ ,italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT , italic_t ) = italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t ) - italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t = 0 ) = - italic_δ roman_sin roman_Ω italic_t ,
dr(rout)=0.subscript𝑑𝑟subscript𝑟out0\displaystyle d_{r}(r_{\rm out})=0\ .italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) = 0 . (7)

The displacement field is measured from the equilibrated configuration obtained after the expansion of the inner boundary to r0+Δsubscript𝑟0Δr_{0}+\Deltaitalic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + roman_Δ. Accordingly we seek a solution of the form

dr(r,t)=fΩ(r)sin(Ωt),subscript𝑑𝑟𝑟𝑡subscript𝑓Ω𝑟Ω𝑡d_{r}(r,t)=f_{\Omega}(r)\sin(\Omega t)\ ,italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r , italic_t ) = italic_f start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_r ) roman_sin ( roman_Ω italic_t ) , (8)

where

r2fΩ′′(r)+rfΩ(r)+fΩ(r)(ρr2Ω22μ+λ1)=0,superscript𝑟2superscriptsubscript𝑓Ω′′𝑟𝑟superscriptsubscript𝑓Ω𝑟subscript𝑓Ω𝑟𝜌superscript𝑟2superscriptΩ22𝜇𝜆10r^{2}f_{\Omega}^{\prime\prime}(r)+rf_{\Omega}^{\prime}(r)+f_{\Omega}(r)\left(% \frac{\rho r^{2}\Omega^{2}}{2\mu+\lambda}-1\right)=0\ ,italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ ′ end_POSTSUPERSCRIPT ( italic_r ) + italic_r italic_f start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_r ) + italic_f start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_r ) ( divide start_ARG italic_ρ italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ + italic_λ end_ARG - 1 ) = 0 , (9)

The dynamic response solves the Bessel equation (9). Interestingly, this equation is of the same form as that found in anomalous static anomalous elasticity cf. Ref. [4]. Here the frequency related expression acts as dynamical screening. The solution of Eq. 9 together with boundary conditions (II.1.1) in terms of the first order Bessel J1(rω~)subscript𝐽1𝑟~𝜔J_{1}(r\tilde{\omega})italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r over~ start_ARG italic_ω end_ARG ) and Neumann Y1(rω~)subscript𝑌1𝑟~𝜔Y_{1}(r\tilde{\omega})italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r over~ start_ARG italic_ω end_ARG ) functions is

fω~(r)=δ[Y1(rω~)J1(routω~)J1(rω~)Y1(routω~)][Y1(rinω~)J1(routω~)J1(rinω~)Y1(routω~)]subscript𝑓~𝜔𝑟𝛿delimited-[]subscript𝑌1𝑟~𝜔subscript𝐽1subscript𝑟out~𝜔subscript𝐽1𝑟~𝜔subscript𝑌1subscript𝑟out~𝜔delimited-[]subscript𝑌1subscript𝑟in~𝜔subscript𝐽1subscript𝑟out~𝜔subscript𝐽1subscript𝑟in~𝜔subscript𝑌1subscript𝑟out~𝜔f_{\tilde{\omega}}(r)=-\delta\frac{\left[Y_{1}\left(r\tilde{\omega}\right)J_{1% }\left(r_{\rm out}\tilde{\omega}\right)-J_{1}\left(r\tilde{\omega}\right)Y_{1}% \left(r_{\rm out}\tilde{\omega}\right)\right]}{\left[Y_{1}\left(r_{\rm in}% \tilde{\omega}\right)J_{1}\left(r_{\rm out}\tilde{\omega}\right)-J_{1}\left(r_% {\rm in}\tilde{\omega}\right)Y_{1}\left(r_{\rm out}\tilde{\omega}\right)\right]}italic_f start_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG end_POSTSUBSCRIPT ( italic_r ) = - italic_δ divide start_ARG [ italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r over~ start_ARG italic_ω end_ARG ) italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r over~ start_ARG italic_ω end_ARG ) italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) ] end_ARG start_ARG [ italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG ) ] end_ARG (10)

where ω~Ωρ/(2μ+λ)~𝜔Ω𝜌2𝜇𝜆\tilde{\omega}\equiv\Omega\sqrt{\rho/(2\mu+\lambda)}over~ start_ARG italic_ω end_ARG ≡ roman_Ω square-root start_ARG italic_ρ / ( 2 italic_μ + italic_λ ) end_ARG

An important feature of this solution is that fω~subscript𝑓~𝜔f_{\tilde{\omega}}italic_f start_POSTSUBSCRIPT over~ start_ARG italic_ω end_ARG end_POSTSUBSCRIPT is fully determined by the geometry of the system, the prescribed frequency, i.e. by rin,rout,Ωsubscript𝑟insubscript𝑟outΩr_{\rm in},r_{\rm out},\Omegaitalic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT , italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT , roman_Ω and the mechanical properties of the media μ𝜇\muitalic_μ and λ𝜆\lambdaitalic_λ. In particular, the number of nodes of the oscillating Bessel functions is determined by the prescribed frequency ΩΩ\Omegaroman_Ω. In order to achieve presentation of the results which does not depend on the system size and material properties, we rewrite the equation and its solutions in dimensionless variables χ𝜒\chiitalic_χ and τ𝜏\tauitalic_τ. This is done using a characteristic length the outer radius routsubscript𝑟outr_{\rm out}italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT and characteristic time tcrout/cdsubscript𝑡𝑐subscript𝑟outsubscript𝑐𝑑t_{c}\equiv r_{\rm out}/c_{d}italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ≡ italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT / italic_c start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT where cdsubscript𝑐𝑑c_{d}italic_c start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT is the dilatational speed

cd(2μ+λ)/ρ.subscript𝑐𝑑2𝜇𝜆𝜌c_{d}\equiv\sqrt{(2\mu+\lambda)/\rho}\ .italic_c start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT ≡ square-root start_ARG ( 2 italic_μ + italic_λ ) / italic_ρ end_ARG . (11)

With this,

χr/rout,τt/tc.formulae-sequence𝜒𝑟subscript𝑟out𝜏𝑡subscript𝑡𝑐\chi\equiv r/r_{\rm out},\quad\tau\equiv t/t_{c}\ .italic_χ ≡ italic_r / italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT , italic_τ ≡ italic_t / italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT . (12)

In dimensionless units Eq. (10) reads

fω(χ)=δχY1(χω)J1(χoutω)J1(χω)Y1(χoutω)Y1(χinω)J1(χoutω)J1(χinω)Y1(χoutω)subscript𝑓𝜔𝜒subscript𝛿𝜒subscript𝑌1𝜒𝜔subscript𝐽1subscript𝜒out𝜔subscript𝐽1𝜒𝜔subscript𝑌1subscript𝜒out𝜔subscript𝑌1subscript𝜒in𝜔subscript𝐽1subscript𝜒out𝜔subscript𝐽1subscript𝜒in𝜔subscript𝑌1subscript𝜒out𝜔f_{\omega}(\chi)=-\delta_{\chi}\frac{Y_{1}\left(\chi\omega\right)J_{1}\left(% \chi_{\rm out}\omega\right)-J_{1}\left(\chi\omega\right)Y_{1}\left(\chi_{\rm out% }\omega\right)}{Y_{1}\left(\chi_{\rm in}\omega\right)J_{1}\left(\chi_{\rm out}% \omega\right)-J_{1}\left(\chi_{\rm in}\omega\right)Y_{1}\left(\chi_{\rm out}% \omega\right)}italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) = - italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT divide start_ARG italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ italic_ω ) italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT italic_ω ) - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ italic_ω ) italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT italic_ω ) end_ARG start_ARG italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT italic_ω ) italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT italic_ω ) - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT italic_ω ) italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT italic_ω ) end_ARG (13)

Note that in this result δχ=δ/routsubscript𝛿𝜒𝛿subscript𝑟out\delta_{\chi}=\delta/r_{\rm out}italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT = italic_δ / italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT, χout=1subscript𝜒out1\chi_{\rm out}=1italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT = 1, χin=rin/routsubscript𝜒insubscript𝑟insubscript𝑟out\chi_{\rm in}=r_{\rm in}/r_{\rm out}italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT = italic_r start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT, and the frequency ω𝜔\omegaitalic_ω is dimensionless.

To get familiarity with the type of functions involved we show in Fig. 1 a few examples of fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) as given by Eq. (13). Later we will learn how dissipation and dipole screening change this function in realistic situations.

II.1.2 Adding dissipation

In reality, in any experimental or model system, dissipation will be crucial to balance the work done by an oscillating pulsar at the center of the circular system. In the simulations, we employ disks that interact via Hertzian repulsive forces (see Appendix A), but we introduce dissipation by adding to the particle dynamics a dumping force proportional to their velocity 𝒗isubscript𝒗𝑖{\bm{v}}_{i}bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, 𝑭iγ~𝒗isubscript𝑭𝑖~𝛾subscript𝒗𝑖{\bm{F}}_{i}\equiv-\tilde{\gamma}{\bm{v}}_{i}bold_italic_F start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ≡ - over~ start_ARG italic_γ end_ARG bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, which is applied to all the disks, to remove any excess energy. In the macroscopic modeling we need therefore to take into account an effective viscous term, that is added to Eq.(3):

ρ𝐝¨+γ(r)𝒅˙=μ2𝐝+(μ+λ)(𝐝).𝜌¨𝐝𝛾𝑟˙𝒅𝜇superscript2𝐝𝜇𝜆𝐝\rho\ddot{\mathbf{d}}+\gamma(r)\dot{{\bm{d}}}=\mu\nabla^{2}\mathbf{d}+(\mu+% \lambda)\mathbf{\nabla}\left(\nabla\cdot\mathbf{d}\right).italic_ρ over¨ start_ARG bold_d end_ARG + italic_γ ( italic_r ) over˙ start_ARG bold_italic_d end_ARG = italic_μ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_d + ( italic_μ + italic_λ ) ∇ ( ∇ ⋅ bold_d ) . (14)

At this point the function γ(r)𝛾𝑟\gamma(r)italic_γ ( italic_r ) is kept unspecified, later it will be chosen to accommodate the results of the numerical simulations. For any choice of γ(r)𝛾𝑟\gamma(r)italic_γ ( italic_r ) the appearance of 𝒅˙˙𝒅\dot{{\bm{d}}}over˙ start_ARG bold_italic_d end_ARG will now mix sine and cosine functions. In dimensionless units we seek a solution for the radial component in the form

dχ(χ,τ)=fω(χ)sin(ωτ)+gω(χ)cos(ωτ).subscript𝑑𝜒𝜒𝜏subscript𝑓𝜔𝜒𝜔𝜏subscript𝑔𝜔𝜒𝜔𝜏d_{\chi}(\chi,\tau)=f_{\omega}(\chi)\sin(\omega\tau)+g_{\omega}(\chi)\cos(% \omega\tau)\ .italic_d start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT ( italic_χ , italic_τ ) = italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) roman_sin ( italic_ω italic_τ ) + italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) roman_cos ( italic_ω italic_τ ) . (15)

Next we solve for the functions fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) by substitution the ansatz Eq. (15) into the boundary value equation Eq. (14), matching together terms in cosωτ𝜔𝜏\cos{\omega\tau}roman_cos italic_ω italic_τ and sinωτ𝜔𝜏\sin{\omega\tau}roman_sin italic_ω italic_τ. We find

ω2gωγ(χ)tcρωfω=1χ2[χ2d2gωdχ2+χdgωdχgω]ω2fω+γ(χ)tcρωgω=1χ2[r2d2fωdχ2+χdfωdχfω].superscript𝜔2subscript𝑔𝜔𝛾𝜒subscript𝑡𝑐𝜌𝜔subscript𝑓𝜔1superscript𝜒2delimited-[]superscript𝜒2superscript𝑑2subscript𝑔𝜔𝑑superscript𝜒2𝜒𝑑subscript𝑔𝜔𝑑𝜒subscript𝑔𝜔superscript𝜔2subscript𝑓𝜔𝛾𝜒subscript𝑡𝑐𝜌𝜔subscript𝑔𝜔1superscript𝜒2delimited-[]superscript𝑟2superscript𝑑2subscript𝑓𝜔𝑑superscript𝜒2𝜒𝑑subscript𝑓𝜔𝑑𝜒subscript𝑓𝜔\displaystyle\begin{split}-\omega^{2}g_{\omega}-\gamma(\chi)\frac{t_{c}}{\rho}% \omega~{}f_{\omega}&=\frac{1}{\chi^{2}}[\chi^{2}\frac{d^{2}g_{\omega}}{d\chi^{% 2}}+\chi\frac{dg_{\omega}}{d\chi}-g_{\omega}]\\ -\omega^{2}f_{\omega}+\gamma(\chi)\frac{t_{c}}{\rho}\omega~{}g_{\omega}&=\frac% {1}{\chi^{2}}[r^{2}\frac{d^{2}f_{\omega}}{d\chi^{2}}+\chi\frac{df_{\omega}}{d% \chi}-f_{\omega}]\ .\end{split}start_ROW start_CELL - italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT - italic_γ ( italic_χ ) divide start_ARG italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ end_ARG italic_ω italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_CELL start_CELL = divide start_ARG 1 end_ARG start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_ARG start_ARG italic_d italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_χ divide start_ARG italic_d italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_ARG start_ARG italic_d italic_χ end_ARG - italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ] end_CELL end_ROW start_ROW start_CELL - italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT + italic_γ ( italic_χ ) divide start_ARG italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ end_ARG italic_ω italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_CELL start_CELL = divide start_ARG 1 end_ARG start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_ARG start_ARG italic_d italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_χ divide start_ARG italic_d italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT end_ARG start_ARG italic_d italic_χ end_ARG - italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ] . end_CELL end_ROW (16)

These two coupled real equations can be solved by introducing one complex variable

z(χ)=gω(χ)+ifω(χ).𝑧𝜒subscript𝑔𝜔𝜒𝑖subscript𝑓𝜔𝜒z(\chi)=g_{\omega}(\chi)+if_{\omega}(\chi)\ .italic_z ( italic_χ ) = italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) + italic_i italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) . (17)

Combining Eqs. (16) according to Eq. (17) we find

ω2ziγ(χ)tcρωz=1χ2[χ2d2zdχ2+χdzdχz].superscript𝜔2𝑧𝑖𝛾𝜒subscript𝑡𝑐𝜌𝜔𝑧1superscript𝜒2delimited-[]superscript𝜒2superscript𝑑2𝑧𝑑superscript𝜒2𝜒𝑑𝑧𝑑𝜒𝑧-\omega^{2}z-i\gamma(\chi)\frac{t_{c}}{\rho}\omega z=\frac{1}{\chi^{2}}[\chi^{% 2}\frac{d^{2}z}{d\chi^{2}}+\chi\frac{dz}{d\chi}-z]\ .- italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_z - italic_i italic_γ ( italic_χ ) divide start_ARG italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ end_ARG italic_ω italic_z = divide start_ARG 1 end_ARG start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_z end_ARG start_ARG italic_d italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_χ divide start_ARG italic_d italic_z end_ARG start_ARG italic_d italic_χ end_ARG - italic_z ] . (18)

We can rewrite Eq. (18) as

χ2d2zdχ2+χdzdχ+z(1+χ2[ω2+iγ(χ)tcρω])=0.superscript𝜒2superscript𝑑2𝑧𝑑superscript𝜒2𝜒𝑑𝑧𝑑𝜒𝑧1superscript𝜒2delimited-[]superscript𝜔2𝑖𝛾𝜒subscript𝑡𝑐𝜌𝜔0\chi^{2}\frac{d^{2}z}{d\chi^{2}}+\chi\frac{dz}{d\chi}+z\Big{(}-1+\chi^{2}[% \omega^{2}+i\gamma(\chi)\frac{t_{c}}{\rho}\omega]\Big{)}=0\ .italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_z end_ARG start_ARG italic_d italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_χ divide start_ARG italic_d italic_z end_ARG start_ARG italic_d italic_χ end_ARG + italic_z ( - 1 + italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_γ ( italic_χ ) divide start_ARG italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ end_ARG italic_ω ] ) = 0 . (19)

The solutions of Eq. (19) are a combination of a first order Bessel J1[ζ(χ)]subscript𝐽1delimited-[]𝜁𝜒J_{1}[\zeta(\chi)]italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ ) ] and Neumann functions Y1[ζ(χ)]subscript𝑌1delimited-[]𝜁𝜒Y_{1}[\zeta(\chi)]italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ ) ] of the complex variable ζ(χ)=χ[ω2+iγ(χ)tcρω]𝜁𝜒𝜒delimited-[]superscript𝜔2𝑖𝛾𝜒subscript𝑡𝑐𝜌𝜔\zeta(\chi)=\chi~{}\sqrt{[\omega^{2}+i\gamma(\chi)\frac{t_{c}}{\rho}\omega]}italic_ζ ( italic_χ ) = italic_χ square-root start_ARG [ italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_γ ( italic_χ ) divide start_ARG italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ end_ARG italic_ω ] end_ARG.

Thus we can write

z(χ)=GJ1(ζ)+HY1(ζ).𝑧𝜒𝐺subscript𝐽1𝜁𝐻subscript𝑌1𝜁z(\chi)=GJ_{1}(\zeta)+HY_{1}(\zeta)\ .italic_z ( italic_χ ) = italic_G italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_ζ ) + italic_H italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_ζ ) . (20)

To determine the coefficients G𝐺Gitalic_G and H𝐻Hitalic_H we need to fit the boundary conditions. They are fω(χout)=gω(χout)=0subscript𝑓𝜔subscript𝜒outsubscript𝑔𝜔subscript𝜒out0f_{\omega}(\chi_{\rm out})=g_{\omega}(\chi_{\rm out})=0italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) = italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) = 0 at the outer boundary. At the oscillating inner boundary fω(χin)=δχsubscript𝑓𝜔subscript𝜒insubscript𝛿𝜒f_{\omega}(\chi_{\rm in})=-\delta_{\chi}italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) = - italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT, gω(χin)=0subscript𝑔𝜔subscript𝜒in0g_{\omega}(\chi_{\rm in})=0italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) = 0. Accordingly, z(χin)=iδχ𝑧subscript𝜒in𝑖subscript𝛿𝜒z(\chi_{\rm in})=-i\delta_{\chi}italic_z ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) = - italic_i italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT, and z(χout)=0𝑧subscript𝜒out0z(\chi_{\rm out})=0italic_z ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) = 0. Using these boundary conditions, we define

D=J1[ζ(χin)]Y1[ζ(χout)]J1[ζ(χout)]Y1[ζ(χin)].𝐷subscript𝐽1delimited-[]𝜁subscript𝜒insubscript𝑌1delimited-[]𝜁subscript𝜒outsubscript𝐽1delimited-[]𝜁subscript𝜒outsubscript𝑌1delimited-[]𝜁subscript𝜒inD=J_{1}[\zeta(\chi_{\rm in})]Y_{1}[\zeta(\chi_{\rm out})]-J_{1}[\zeta(\chi_{% \rm out})]Y_{1}[\zeta(\chi_{\rm in})]\ .italic_D = italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) ] italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) ] . (21)

The coefficients in Eq. (20) are computed to be

G=iδY1[ζ(χout)]/D,H=iδJ1[ζ(χout)]/D.formulae-sequence𝐺𝑖𝛿subscript𝑌1delimited-[]𝜁subscript𝜒out𝐷𝐻𝑖𝛿subscript𝐽1delimited-[]𝜁subscript𝜒out𝐷G=-i\delta Y_{1}[\zeta(\chi_{\rm out})]/D\,,\quad H=i\delta J_{1}[\zeta(\chi_{% \rm out})]/D\ .italic_G = - italic_i italic_δ italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] / italic_D , italic_H = italic_i italic_δ italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_ζ ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] / italic_D . (22)

Finally, fωsubscript𝑓𝜔f_{\omega}italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT and gωsubscript𝑔𝜔g_{\omega}italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT of Eq. (15) are determined as (z)𝑧\Re(z)roman_ℜ ( italic_z ) and (z)𝑧\Im(z)roman_ℑ ( italic_z ).

These exact solutions are expected to be relevant as long as dipole screening is absent. To assess the effects of the latter we now add quadrupole and dipole contributions to the theoretical discussion.

III Anomalous Dynamics due to plastic Screening

In this section we derive the equations of motion in the presence of quadrupolar and dipolar screening. This derivation follows the ideas presented in Ref. [4] for the static problem, supplemented with inertial and dissipative contributions as required. In each case we need to write down the appropriate Lagrangian which reflects the interactions that are taken into account. We denote the quadrupolar tensor field as Qαβsuperscript𝑄𝛼𝛽Q^{\alpha\beta}italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT. Physically, this is the field associated with the eigenstrains of the Eshelby quadrupoles that are formed by plastic responses. The dipolar field PαβQαβsuperscript𝑃𝛼subscript𝛽superscript𝑄𝛼𝛽P^{\alpha}\equiv\partial_{\beta}Q^{\alpha\beta}italic_P start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT ≡ ∂ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT is simply the gradient of this density. The displacement field is the vector field dαsubscript𝑑𝛼d_{\alpha}italic_d start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT and the strain tensor uαβ12(αdβ+βdα)subscript𝑢𝛼𝛽12subscript𝛼subscript𝑑𝛽subscript𝛽subscript𝑑𝛼u_{\alpha\beta}\equiv\tfrac{1}{2}(\partial_{\alpha}d_{\beta}+\partial_{\beta}d% _{\alpha})italic_u start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ≡ divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( ∂ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT + ∂ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ). The Euclidean metric is denoted gμνsubscript𝑔𝜇𝜈g_{\mu\nu}italic_g start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT.

In all cases we start with the Euler-Lagrange equations in the standard form [12], with ΦΦ\Phiroman_Φ playing the role of a fundamental field, which in our case can be the displacement or the quadrupole density.

Φν((νΦ))=0.Φsubscript𝜈subscript𝜈Φ0\frac{\partial\mathcal{L}}{\partial\Phi}-\partial_{\nu}\left(\frac{\partial% \mathcal{L}}{\partial\left(\partial_{\nu}\Phi\right)}\right)=0\ .divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ roman_Φ end_ARG - ∂ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ ( ∂ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT roman_Φ ) end_ARG ) = 0 . (23)

III.1 Quadrupole Screening

To derive the equations of motion under quadrupole screening we assume that the density of quadrupoles is sufficiently low, or that the gradients of their density are negligible. Thus to quadratic order we can write the Lagrangian for an elastic medium with quadrupoles as

\displaystyle\mathcal{L}caligraphic_L =ρ2gμνd˙μd˙ν12σμνuμν12ΛαβγδQαβQγδabsent𝜌2superscript𝑔𝜇𝜈subscript˙𝑑𝜇subscript˙𝑑𝜈12superscript𝜎𝜇𝜈subscript𝑢𝜇𝜈12subscriptΛ𝛼𝛽𝛾𝛿superscript𝑄𝛼𝛽superscript𝑄𝛾𝛿\displaystyle=\frac{\rho}{2}g^{\mu\nu}\dot{d}_{\mu}\dot{d}_{\nu}-\frac{1}{2}% \sigma^{\mu\nu}u_{\mu\nu}-\frac{1}{2}\Lambda_{\alpha\beta\gamma\delta}Q^{% \alpha\beta}Q^{\gamma\delta}= divide start_ARG italic_ρ end_ARG start_ARG 2 end_ARG italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Λ start_POSTSUBSCRIPT italic_α italic_β italic_γ italic_δ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT (24)
ΓγδαβαdβQγδsubscriptsuperscriptΓ𝛼𝛽𝛾𝛿subscript𝛼subscript𝑑𝛽superscript𝑄𝛾𝛿\displaystyle-\Gamma^{\alpha\beta}_{\gamma\delta}\partial_{\alpha}d_{\beta}Q^{% \gamma\delta}- roman_Γ start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT

Computing the derivative with respect to the quadrupole field we find

Qσκsuperscript𝑄𝜎𝜅\displaystyle\frac{\partial\mathcal{L}}{\partial Q^{\sigma\kappa}}divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ italic_Q start_POSTSUPERSCRIPT italic_σ italic_κ end_POSTSUPERSCRIPT end_ARG =\displaystyle== ΛαβγδQαβδσγδκδΓγδαβαdβδσγδκδsubscriptΛ𝛼𝛽𝛾𝛿superscript𝑄𝛼𝛽superscriptsubscript𝛿𝜎𝛾superscriptsubscript𝛿𝜅𝛿subscriptsuperscriptΓ𝛼𝛽𝛾𝛿subscript𝛼subscript𝑑𝛽superscriptsubscript𝛿𝜎𝛾superscriptsubscript𝛿𝜅𝛿\displaystyle-\Lambda_{\alpha\beta\gamma\delta}Q^{\alpha\beta}\delta_{\sigma}^% {\gamma}\delta_{\kappa}^{\delta}-\Gamma^{\alpha\beta}_{\gamma\delta}\partial_{% \alpha}d_{\beta}\delta_{\sigma}^{\gamma}\delta_{\kappa}^{\delta}- roman_Λ start_POSTSUBSCRIPT italic_α italic_β italic_γ italic_δ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT italic_δ start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_γ end_POSTSUPERSCRIPT italic_δ start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_δ end_POSTSUPERSCRIPT - roman_Γ start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_δ start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_γ end_POSTSUPERSCRIPT italic_δ start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_δ end_POSTSUPERSCRIPT (25)
=\displaystyle== ΛαβσκQαβΓσκαβuαβ=0.subscriptΛ𝛼𝛽𝜎𝜅superscript𝑄𝛼𝛽subscriptsuperscriptΓ𝛼𝛽𝜎𝜅subscript𝑢𝛼𝛽0\displaystyle-\Lambda_{\alpha\beta\sigma\kappa}Q^{\alpha\beta}-\Gamma^{\alpha% \beta}_{\sigma\kappa}u_{\alpha\beta}=0\ .- roman_Λ start_POSTSUBSCRIPT italic_α italic_β italic_σ italic_κ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT - roman_Γ start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_σ italic_κ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = 0 .

At this point we do not have dipoles in the Lagrangian, or

(νQσκ)=0.subscript𝜈superscript𝑄𝜎𝜅0\frac{\partial\mathcal{L}}{\partial\left(\partial_{\nu}Q^{\sigma\kappa}\right)% }=0\ .divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ ( ∂ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_σ italic_κ end_POSTSUPERSCRIPT ) end_ARG = 0 . (26)

Returning to Eq. (25) we rewrite it in the form

Qαβ=ΛαβσκΓσκγδuγδ=Γ~αβγδuγδsuperscript𝑄𝛼𝛽superscriptΛ𝛼𝛽𝜎𝜅subscriptsuperscriptΓ𝛾𝛿𝜎𝜅subscript𝑢𝛾𝛿superscript~Γ𝛼𝛽𝛾𝛿subscript𝑢𝛾𝛿Q^{\alpha\beta}=-\Lambda^{\alpha\beta\sigma\kappa}\Gamma^{\gamma\delta}_{% \sigma\kappa}u_{\gamma\delta}=\tilde{\Gamma}^{\alpha\beta\gamma\delta}u_{% \gamma\delta}italic_Q start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT = - roman_Λ start_POSTSUPERSCRIPT italic_α italic_β italic_σ italic_κ end_POSTSUPERSCRIPT roman_Γ start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_σ italic_κ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT = over~ start_ARG roman_Γ end_ARG start_POSTSUPERSCRIPT italic_α italic_β italic_γ italic_δ end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT (27)

Here we denoted the inverse of ΛαβσκsubscriptΛ𝛼𝛽𝜎𝜅\Lambda_{\alpha\beta\sigma\kappa}roman_Λ start_POSTSUBSCRIPT italic_α italic_β italic_σ italic_κ end_POSTSUBSCRIPT as ΛαβσκsuperscriptΛ𝛼𝛽𝜎𝜅\Lambda^{\alpha\beta\sigma\kappa}roman_Λ start_POSTSUPERSCRIPT italic_α italic_β italic_σ italic_κ end_POSTSUPERSCRIPT and defined a new (presently unknown) tensor of coefficients 𝚪~bold-~𝚪{\bm{\tilde{\Gamma}}}overbold_~ start_ARG bold_Γ end_ARG. Eq. (27) is important, showing that the quadrupolar tensor is not some mysterious entity, but that it is induced by the strain field in the system, with the quartic tensor 𝚪~bold-~𝚪{\bm{\tilde{\Gamma}}}overbold_~ start_ARG bold_Γ end_ARG providing the link.

Next we consider the derivative with respect to ρdσsubscript𝜌subscript𝑑𝜎\partial_{\rho}d_{\sigma}∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT:

(ρdσ)=ρgμνd˙μδtρδνσσμνδμρδνσΓγδαβδαρδβσQγδsubscript𝜌subscript𝑑𝜎𝜌superscript𝑔𝜇𝜈subscript˙𝑑𝜇subscriptsuperscript𝛿𝜌𝑡subscriptsuperscript𝛿𝜎𝜈superscript𝜎𝜇𝜈subscriptsuperscript𝛿𝜌𝜇subscriptsuperscript𝛿𝜎𝜈subscriptsuperscriptΓ𝛼𝛽𝛾𝛿subscriptsuperscript𝛿𝜌𝛼subscriptsuperscript𝛿𝜎𝛽superscript𝑄𝛾𝛿\frac{\partial\mathcal{L}}{\partial\left(\partial_{\rho}d_{\sigma}\right)}=% \rho g^{\mu\nu}\dot{d}_{\mu}\delta^{\rho}_{t}\delta^{\sigma}_{\nu}-\sigma^{\mu% \nu}\delta^{\rho}_{\mu}\delta^{\sigma}_{\nu}-\Gamma^{\alpha\beta}_{\gamma% \delta}\delta^{\rho}_{\alpha}\delta^{\sigma}_{\beta}Q^{\gamma\delta}divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ ( ∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ) end_ARG = italic_ρ italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_ρ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT italic_δ start_POSTSUPERSCRIPT italic_ρ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - roman_Γ start_POSTSUPERSCRIPT italic_α italic_β end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_ρ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT (28)

Using Eq. (23) we find

ρgμσd¨μρ(σρσ+ΓγδρσQγδ)=0𝜌superscript𝑔𝜇𝜎subscript¨𝑑𝜇subscript𝜌superscript𝜎𝜌𝜎subscriptsuperscriptΓ𝜌𝜎𝛾𝛿superscript𝑄𝛾𝛿0\rho g^{\mu\sigma}\ddot{d}_{\mu}-\partial_{\rho}\left(\sigma^{\rho\sigma}+% \Gamma^{\rho\sigma}_{\gamma\delta}Q^{\gamma\delta}\right)=0italic_ρ italic_g start_POSTSUPERSCRIPT italic_μ italic_σ end_POSTSUPERSCRIPT over¨ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT - ∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT ( italic_σ start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT + roman_Γ start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT ) = 0 (29)

At this point we can define a normalized stress field and write

ρd¨σ=ρσ~ρσ,σ~ρσσρσ+ΓγδρσQγδ=ξσρσ,formulae-sequence𝜌superscript¨𝑑𝜎subscript𝜌superscript~𝜎𝜌𝜎superscript~𝜎𝜌𝜎superscript𝜎𝜌𝜎subscriptsuperscriptΓ𝜌𝜎𝛾𝛿superscript𝑄𝛾𝛿𝜉superscript𝜎𝜌𝜎\rho\ddot{d}^{\sigma}=\partial_{\rho}\tilde{\sigma}^{\rho\sigma}\ ,\quad\tilde% {\sigma}^{\rho\sigma}\equiv\sigma^{\rho\sigma}+\Gamma^{\rho\sigma}_{\gamma% \delta}Q^{\gamma\delta}=\xi\sigma^{\rho\sigma},italic_ρ over¨ start_ARG italic_d end_ARG start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT = ∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT over~ start_ARG italic_σ end_ARG start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT , over~ start_ARG italic_σ end_ARG start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT ≡ italic_σ start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT + roman_Γ start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_γ italic_δ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_δ end_POSTSUPERSCRIPT = italic_ξ italic_σ start_POSTSUPERSCRIPT italic_ρ italic_σ end_POSTSUPERSCRIPT , (30)

where the last equality follows from Eq. (27) and homogeneity and isotropy, with ξ𝜉\xiitalic_ξ being a scalar number. The upshot of this calculation is that the dynamics of the displacement field is unchanged in form compared to Eq. (2), except for a renormalization of the elastic moduli. Thus in the present case the introduction of the dissipative term and the analytic solutions for the displacement field follow verbatim the theory presented in the previous section.

III.2 Dipole Screening

The situation changes qualitatively with dipole screening, the equations of motion change their form. To derive these equations we will assert that the renormalization due to quadrupole screening is already included in our Lagrangian, in the form of a renormalized stress tensor (that will be again denoted as σμνsuperscript𝜎𝜇𝜈\sigma^{\mu\nu}italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT). To quadratic order in the fields we therefore write

\displaystyle\mathcal{L}caligraphic_L =\displaystyle== ρ2gμνd˙μd˙ν12σμνuμν12ΛαβPαPβΓβαdαPβ𝜌2superscript𝑔𝜇𝜈subscript˙𝑑𝜇subscript˙𝑑𝜈12superscript𝜎𝜇𝜈subscript𝑢𝜇𝜈12subscriptΛ𝛼𝛽superscript𝑃𝛼superscript𝑃𝛽subscriptsuperscriptΓ𝛼𝛽subscript𝑑𝛼superscript𝑃𝛽\displaystyle\frac{\rho}{2}g^{\mu\nu}\dot{d}_{\mu}\dot{d}_{\nu}-\frac{1}{2}% \sigma^{\mu\nu}u_{\mu\nu}-\frac{1}{2}\Lambda_{\alpha\beta}P^{\alpha}P^{\beta}-% \Gamma^{\alpha}_{\beta}d_{\alpha}P^{\beta}divide start_ARG italic_ρ end_ARG start_ARG 2 end_ARG italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Λ start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT italic_P start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT - roman_Γ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT (31)
=\displaystyle== ρ2gμνd˙μd˙ν12σμνuμν12ΛαβγQγαδQδβ𝜌2superscript𝑔𝜇𝜈subscript˙𝑑𝜇subscript˙𝑑𝜈12superscript𝜎𝜇𝜈subscript𝑢𝜇𝜈12subscriptΛ𝛼𝛽subscript𝛾superscript𝑄𝛾𝛼subscript𝛿superscript𝑄𝛿𝛽\displaystyle\frac{\rho}{2}g^{\mu\nu}\dot{d}_{\mu}\dot{d}_{\nu}-\frac{1}{2}% \sigma^{\mu\nu}u_{\mu\nu}-\frac{1}{2}\Lambda_{\alpha\beta}\partial_{\gamma}Q^{% \gamma\alpha}\partial_{\delta}Q^{\delta\beta}divide start_ARG italic_ρ end_ARG start_ARG 2 end_ARG italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Λ start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_γ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_α end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_δ italic_β end_POSTSUPERSCRIPT
\displaystyle-- ΓβαdαγQγβ.subscriptsuperscriptΓ𝛼𝛽subscript𝑑𝛼subscript𝛾superscript𝑄𝛾𝛽\displaystyle\Gamma^{\alpha}_{\beta}d_{\alpha}\partial_{\gamma}Q^{\gamma\beta}\ .roman_Γ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_γ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_β end_POSTSUPERSCRIPT .

Computing the derivative with respect to the dipole field we find

ρ(ρQσκ)=σ[ΛακPα+Γκαdα]=0.subscript𝜌subscript𝜌superscript𝑄𝜎𝜅subscript𝜎delimited-[]subscriptΛ𝛼𝜅superscript𝑃𝛼subscriptsuperscriptΓ𝛼𝜅subscript𝑑𝛼0{\partial_{\rho}}\frac{\partial\mathcal{L}}{\partial\left(\partial_{\rho}Q^{% \sigma\kappa}\right)}=-\partial_{\sigma}[\Lambda_{\alpha\kappa}P^{\alpha}+% \Gamma^{\alpha}_{\kappa}d_{\alpha}]=0\ .∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ ( ∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_σ italic_κ end_POSTSUPERSCRIPT ) end_ARG = - ∂ start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT [ roman_Λ start_POSTSUBSCRIPT italic_α italic_κ end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT + roman_Γ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ] = 0 . (32)

The expression in the square brackets is a constant, that can be taken as zero using the translational invariance of the displacement field. From this point onward we lose the gauge freedom of the displacement field since we chose a gauge. As before, we use this equation to express the dipole field in terms of the fundamental displacement field:

Pα=ΛακΓκβdβ.superscript𝑃𝛼superscriptΛ𝛼𝜅subscriptsuperscriptΓ𝛽𝜅subscript𝑑𝛽P^{\alpha}=-\Lambda^{\alpha\kappa}\Gamma^{\beta}_{\kappa}d_{\beta}\ .italic_P start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT = - roman_Λ start_POSTSUPERSCRIPT italic_α italic_κ end_POSTSUPERSCRIPT roman_Γ start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT . (33)

Computing the derivatives with respect to the displacement field and its derivative we find

dσ=ΓβαγQγβδασ=ΓβαPβδασ,subscript𝑑𝜎subscriptsuperscriptΓ𝛼𝛽subscript𝛾superscript𝑄𝛾𝛽subscriptsuperscript𝛿𝜎𝛼subscriptsuperscriptΓ𝛼𝛽superscript𝑃𝛽subscriptsuperscript𝛿𝜎𝛼\frac{\partial\mathcal{L}}{\partial d_{\sigma}}=-\Gamma^{\alpha}_{\beta}% \partial_{\gamma}Q^{\gamma\beta}\delta^{\sigma}_{\alpha}=-\Gamma^{\alpha}_{% \beta}P^{\beta}\delta^{\sigma}_{\alpha}\ ,divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ italic_d start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT end_ARG = - roman_Γ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_γ end_POSTSUBSCRIPT italic_Q start_POSTSUPERSCRIPT italic_γ italic_β end_POSTSUPERSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT = - roman_Γ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT , (34)
(ρdσ)=ρgμνd˙μδtρδνσσμνδμρδνσ.subscript𝜌subscript𝑑𝜎𝜌superscript𝑔𝜇𝜈subscript˙𝑑𝜇subscriptsuperscript𝛿𝜌𝑡subscriptsuperscript𝛿𝜎𝜈superscript𝜎𝜇𝜈subscriptsuperscript𝛿𝜌𝜇subscriptsuperscript𝛿𝜎𝜈\frac{\partial\mathcal{L}}{\partial\left(\partial_{\rho}d_{\sigma}\right)}=% \rho g^{\mu\nu}\dot{d}_{\mu}\delta^{\rho}_{t}\delta^{\sigma}_{\nu}-\sigma^{\mu% \nu}\delta^{\rho}_{\mu}\delta^{\sigma}_{\nu}\ .divide start_ARG ∂ caligraphic_L end_ARG start_ARG ∂ ( ∂ start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ) end_ARG = italic_ρ italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT over˙ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_ρ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - italic_σ start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT italic_δ start_POSTSUPERSCRIPT italic_ρ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_δ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT . (35)

Combining these two equations together we write

ΓβσPβgμσd¨μ+μσμσ=0.subscriptsuperscriptΓ𝜎𝛽superscript𝑃𝛽superscript𝑔𝜇𝜎subscript¨𝑑𝜇subscript𝜇superscript𝜎𝜇𝜎0-\Gamma^{\sigma}_{\beta}P^{\beta}-g^{\mu\sigma}\ddot{d}_{\mu}+\partial_{\mu}% \sigma^{\mu\sigma}=0\ .- roman_Γ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT - italic_g start_POSTSUPERSCRIPT italic_μ italic_σ end_POSTSUPERSCRIPT over¨ start_ARG italic_d end_ARG start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT + ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_σ start_POSTSUPERSCRIPT italic_μ italic_σ end_POSTSUPERSCRIPT = 0 . (36)

Inverting for d¨μsuperscript¨𝑑𝜇\ddot{d}^{\mu}over¨ start_ARG italic_d end_ARG start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT we find the modified equation of motion

d¨μ=μσμσΓβσΛβκΓκγdγsuperscript¨𝑑𝜇subscript𝜇superscript𝜎𝜇𝜎subscriptsuperscriptΓ𝜎𝛽superscriptΛ𝛽𝜅subscriptsuperscriptΓ𝛾𝜅subscript𝑑𝛾\ddot{d}^{\mu}=\partial_{\mu}\sigma^{\mu\sigma}-\Gamma^{\sigma}_{\beta}\Lambda% ^{\beta\kappa}\Gamma^{\gamma}_{\kappa}d_{\gamma}over¨ start_ARG italic_d end_ARG start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_σ start_POSTSUPERSCRIPT italic_μ italic_σ end_POSTSUPERSCRIPT - roman_Γ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT roman_Λ start_POSTSUPERSCRIPT italic_β italic_κ end_POSTSUPERSCRIPT roman_Γ start_POSTSUPERSCRIPT italic_γ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT italic_d start_POSTSUBSCRIPT italic_γ end_POSTSUBSCRIPT (37)

In a homogeneous and isotropic medium we find

ΓβσΛβκΓκγ=κ2gαγ.subscriptsuperscriptΓ𝜎𝛽superscriptΛ𝛽𝜅subscriptsuperscriptΓ𝛾𝜅superscript𝜅2superscript𝑔𝛼𝛾\Gamma^{\sigma}_{\beta}\Lambda^{\beta\kappa}\Gamma^{\gamma}_{\kappa}=\kappa^{2% }g^{\alpha\gamma}\;.roman_Γ start_POSTSUPERSCRIPT italic_σ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT roman_Λ start_POSTSUPERSCRIPT italic_β italic_κ end_POSTSUPERSCRIPT roman_Γ start_POSTSUPERSCRIPT italic_γ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT = italic_κ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_g start_POSTSUPERSCRIPT italic_α italic_γ end_POSTSUPERSCRIPT . (38)

We note that Eq. (37) breaks translational symmetry together with the introduction of a length scale, since κ𝜅\kappaitalic_κ has the dimension of an inverse scale. This is further discussed at great length (and excuse the pun) in Sect. V below.

III.3 Solutions of Anomalous Dynamic dipole screening

As done before, we write σ𝜎\nabla\cdot\sigma∇ ⋅ italic_σ in terms of the displacement, to get the explicit equation in the case of dipole screening. As before, we add the dissipative term to the resulting equation and end up with

ρ𝒅¨+γ(r)𝒅˙=μ2𝐝+(μ+λ)(𝐝)+κ2𝐝.𝜌¨𝒅𝛾𝑟˙𝒅𝜇superscript2𝐝𝜇𝜆𝐝superscript𝜅2𝐝\rho\ddot{{\bm{d}}}+\gamma(r)\dot{{\bm{d}}}=\mu\nabla^{2}{\bf d}+(\mu+\lambda)% \nabla(\nabla\cdot{\bf d})+\kappa^{2}{\bf d}\,.italic_ρ over¨ start_ARG bold_italic_d end_ARG + italic_γ ( italic_r ) over˙ start_ARG bold_italic_d end_ARG = italic_μ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_d + ( italic_μ + italic_λ ) ∇ ( ∇ ⋅ bold_d ) + italic_κ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_d . (39)

We note that the addition of the screening term κ2𝒅superscript𝜅2𝒅\kappa^{2}{\bm{d}}italic_κ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_italic_d does not change the linearity of the equation, and therefore we can use the solution presented in Subsect. II.1.2 with very little modification. Referring again to a solution in the form of Eq. (15), to solve for fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) we only need to change the variable ζ𝜁\zetaitalic_ζ to a new variable ζ~~𝜁\tilde{\zeta}over~ start_ARG italic_ζ end_ARG where

ζ~(χ)χκ~2+ω2+iγ~(χ)ω,~𝜁𝜒𝜒superscript~𝜅2superscript𝜔2𝑖~𝛾𝜒𝜔\tilde{\zeta}(\chi)\equiv\chi\sqrt{\tilde{\kappa}^{2}+\omega^{2}+i\tilde{% \gamma}(\chi)~{}\omega}\ ,over~ start_ARG italic_ζ end_ARG ( italic_χ ) ≡ italic_χ square-root start_ARG over~ start_ARG italic_κ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i over~ start_ARG italic_γ end_ARG ( italic_χ ) italic_ω end_ARG , (40)

where κ~=κtc/ρ,γ~=γ(χ)tc/ρformulae-sequence~𝜅𝜅subscript𝑡𝑐𝜌~𝛾𝛾𝜒subscript𝑡𝑐𝜌\tilde{\kappa}=\kappa~{}t_{c}/\sqrt{\rho},\tilde{\gamma}=\gamma(\chi)~{}t_{c}/\rhoover~ start_ARG italic_κ end_ARG = italic_κ italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT / square-root start_ARG italic_ρ end_ARG , over~ start_ARG italic_γ end_ARG = italic_γ ( italic_χ ) italic_t start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT / italic_ρ. The Eqs. (21) and (22) are unchanged except that the everywhere ζ𝜁\zetaitalic_ζ is replaced by ζ~~𝜁\tilde{\zeta}over~ start_ARG italic_ζ end_ARG.

To provide a feeling to the nature of the solutions of Eq. (39) we show in Figs. 2 and 3 the functions fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) for various frequencies and different values of the screening parameter κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG. In all cases the dissipation function γ~(χ)=γ0[χ/χin(t=0)]ϵ~𝛾𝜒subscript𝛾0superscriptdelimited-[]𝜒subscript𝜒in𝑡0italic-ϵ\tilde{\gamma}(\chi)=\gamma_{0}[\chi/\chi_{\rm in}(t=0)]^{\epsilon}over~ start_ARG italic_γ end_ARG ( italic_χ ) = italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT [ italic_χ / italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_t = 0 ) ] start_POSTSUPERSCRIPT italic_ϵ end_POSTSUPERSCRIPT with γ0=1.0subscript𝛾01.0\gamma_{0}=1.0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 1.0 and ϵ=0.5italic-ϵ0.5\epsilon=0.5italic_ϵ = 0.5. The choice of this value of ϵitalic-ϵ\epsilonitalic_ϵ has been made purely by comparison to simulations. But one can rationalize it by noting that the geometric spreading of 2-dimensional waves originating from a point source is expected to decay like r0.5superscript𝑟0.5r^{-0.5}italic_r start_POSTSUPERSCRIPT - 0.5 end_POSTSUPERSCRIPT [13]. On the other hand, since the number of collisions between disks is increasing like r𝑟ritalic_r, in balance one can expect that a dissipation that increases like r0.5superscript𝑟0.5r^{0.5}italic_r start_POSTSUPERSCRIPT 0.5 end_POSTSUPERSCRIPT should suffice to keep the energy budget finite.

Refer to caption
Figure 2: (a) and (b) The functions gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) as predicted with the screening theory, for fixed driving frequency ω=1.5𝜔1.5\omega=1.5italic_ω = 1.5, and fixed dissipation, γ0=1.0subscript𝛾01.0\gamma_{0}=1.0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 1.0 and ϵ=0.5italic-ϵ0.5\epsilon=0.5italic_ϵ = 0.5, for three values of the screening parameter κ~=3,5~𝜅35\tilde{\kappa}=3,5over~ start_ARG italic_κ end_ARG = 3 , 5, and 7777 (dot-dashed green, dashed red and solid blue line , respectively. The boundary value δχ=0.3subscript𝛿𝜒0.3\delta_{\chi}=0.3italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT = 0.3.
Refer to caption
Figure 3: (a) and (b) The functions gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) with the same dissipation parameters as in Fig. 2 and the screening parameter κ=7𝜅7\kappa=7italic_κ = 7, for three values of the frequency, ω=1.5𝜔1.5\omega=1.5italic_ω = 1.5 (solid dark blue line– the same as in Fig. 2), 5555 (dashed light blue line), and 7.57.57.57.5 (dot-dashed cyan line).
Refer to caption
Figure 4: (a) Comparison of the measured functions gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) to the predictions of the screening theory, taking account of dissipation. The frequencies are ω=1.90𝜔1.90\omega=1.90italic_ω = 1.90 (green circles), ω=4.44𝜔4.44\omega=4.44italic_ω = 4.44 (red squares) and ω=6.34𝜔6.34\omega=6.34italic_ω = 6.34 (blue diamonds). The respective fitting parameters are γ0=0.75,κ~=1.8formulae-sequencesubscript𝛾00.75~𝜅1.8\gamma_{0}=0.75,\tilde{\kappa}=1.8italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0.75 , over~ start_ARG italic_κ end_ARG = 1.8 (green line), γ0=0.9,κ=4.5formulae-sequencesubscript𝛾00.9𝜅4.5\gamma_{0}=0.9,\kappa=4.5italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0.9 , italic_κ = 4.5 (orange line) and γ0=1,κ=6.5formulae-sequencesubscript𝛾01𝜅6.5\gamma_{0}=1,\kappa=6.5italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 1 , italic_κ = 6.5 light blue line). The value of the displacement at the boundary δ1subscript𝛿1\delta_{1}italic_δ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT differ from the amplitude of the oscillations δχsubscript𝛿𝜒\delta_{\chi}italic_δ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT due to finite size of the angular-averaging band. For the shown examples the values of δ1subscript𝛿1\delta_{1}italic_δ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and δ2subscript𝛿2\delta_{2}italic_δ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT (cf. Eq. (41)) correspond to the values of fωsubscript𝑓𝜔f_{\omega}italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPTand gωsubscript𝑔𝜔g_{\omega}italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT, respectively, at the first available point of χ𝜒\chiitalic_χ.
Refer to caption
Figure 5: The frequency dependence of (a) the anomalous screening parameter κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG and (b) the prefactor of the damping term γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. The values of κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG and γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT were averaged over ten different configurations. The error-bars correspond to the standard deviation from the mean values. The lines serve to guide the eye only.

IV Numerical Simulations

IV.1 Setup

The numerical setup consists of a 2D concentric circular enclosure with an inner circle of radius r0subscript𝑟0r_{0}italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, filled with a bi-disperse mixture of frictionless disks. The disks interact through normal Herzian forces only [see Appendix A]. The interaction with the outer wall is also Hertzian, with the interaction coefficient equal to the interaction coefficient Knsubscript𝐾𝑛K_{n}italic_K start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT between disks. Additional viscous dumping force 𝑭i=γ¯𝒗isubscript𝑭𝑖¯𝛾subscript𝒗𝑖{\bm{F}}_{i}=-\bar{\gamma}{\bm{v}}_{i}bold_italic_F start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = - over¯ start_ARG italic_γ end_ARG bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is applied to all discs to remove excess of energy from the system. Here γ¯¯𝛾\bar{\gamma}over¯ start_ARG italic_γ end_ARG is the microscopic damping coefficient (to be distinguished from the macroscopic γ(r)𝛾𝑟\gamma(r)italic_γ ( italic_r ) of Eq. (14)), and 𝒗isubscript𝒗𝑖{\bm{v}}_{i}bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the velocity of the i𝑖iitalic_i’th disc. The simulations are carried out using LAMMPS [14] for the dynamics of the granular system. The driving is realized by periodic inflation of the radius of the inner enclosure wall which is placed at the center of the system.

The granular system was prepared at a desired packing fraction and pressure, starting with the creation of a randomly distributed set of N=10,000𝑁10000N=10,000italic_N = 10 , 000 particles, half of which has a radius ra=0.5subscript𝑟𝑎0.5r_{a}=0.5italic_r start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = 0.5 and and the other half with rb=0.7subscript𝑟𝑏0.7r_{b}=0.7italic_r start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT = 0.7, with r0=4subscript𝑟04r_{0}=4italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 4. Next the system is equilibrated by energy minimization. Then the radius of the enclosing circle routsubscript𝑟outr_{\rm out}italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT is decreased by small decrements from the initial radius to a value corresponding to the largest desired packing fraction ϕitalic-ϕ\phiitalic_ϕ according to rout=(ra2+rb2)N/2ϕ+r02subscript𝑟outsubscriptsuperscript𝑟2𝑎subscriptsuperscript𝑟2𝑏𝑁2italic-ϕsuperscriptsubscript𝑟02r_{\rm out}=\sqrt{(r^{2}_{a}+r^{2}_{b})N/2\phi+r_{0}^{2}}italic_r start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT = square-root start_ARG ( italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT + italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ) italic_N / 2 italic_ϕ + italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG. At each step, the system is again equilibrated. A set of configurations for various pressures p𝑝pitalic_p, or equivalently packing fractions ϕitalic-ϕ\phiitalic_ϕ, was obtained.

Having a system at mechanical equilibrium we can start oscillating the radius of the inner circle. The procedure followed was explained in Subsect. II.1.1. As stated there, to see reproducible results it is important to guarantee that the oscillations take place around a true stationary state, in the sense that one is guaranteed that stopping the oscillations at any point and performing energy minimization would result in the same mechanical equilibrium independent of when and where the oscillatory driving is stopped.

Another lesson from the numerical experiments is that for the present system of Hertzian disks one is limited in the frequencies of driving that can be simulated. Increasing the frequency too much results in the creation of a hole around the oscillating boundary - the system does not have enough time to return to contact with the boundary when this boundary recedes back too rapidly. We therefore limited the simulations, and the comparison with the theory, to relatively low frequencies as we will see next. It is likely that in systems in which there are attractive forces between the constituents the creation of the hole can be eliminated, and higher frequencies could be studied. On the other hand, it is possible that for any system higher frequencies would trigger nonlinear interactions, requiring additional scrutiny of the equations of motion. At this point we leave these interesting questions to future study.

IV.2 Results for the dominant radial mode

To compare the theory with the numerical experiments, we note that the measured functions gω(χ),fω(χ)subscript𝑔𝜔𝜒subscript𝑓𝜔𝜒g_{\omega}(\chi),f_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) , italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) are calculated by angle-averaging over a band of finite width and therefore we amend the boundary conditions at the inner boundary as

fω(χin)=δ1,gω(χin)=δ2,formulae-sequencesubscript𝑓𝜔subscript𝜒insubscript𝛿1subscript𝑔𝜔subscript𝜒insubscript𝛿2f_{\omega}(\chi_{\rm in})=-\delta_{1}\ ,\quad g_{\omega}(\chi_{\rm in})=\delta% _{2}\ ,italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) = - italic_δ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) = italic_δ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , (41)

giving

G𝐺\displaystyle Gitalic_G =\displaystyle== (δ2iδ1)Y1[ζ~(χout)]/D,subscript𝛿2𝑖subscript𝛿1subscript𝑌1delimited-[]~𝜁subscript𝜒out𝐷\displaystyle(\delta_{2}-i\delta_{1})Y_{1}[\tilde{\zeta}(\chi_{\rm out})]/D\,,( italic_δ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_i italic_δ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] / italic_D ,
H𝐻\displaystyle\ Hitalic_H =\displaystyle== (δ2+iδ1)J1[ζ~(χout)]/D,subscript𝛿2𝑖subscript𝛿1subscript𝐽1delimited-[]~𝜁subscript𝜒out𝐷\displaystyle(-\delta_{2}+i\delta_{1})J_{1}[\tilde{\zeta}(\chi_{\rm out})]/D\,,( - italic_δ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_i italic_δ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] / italic_D , (42)
D𝐷\displaystyle Ditalic_D =\displaystyle== J1[ζ~(χin)]Y1[ζ~(χout)]J1[ζ~(χout)]Y1[ζ~(χin)],subscript𝐽1delimited-[]~𝜁subscript𝜒insubscript𝑌1delimited-[]~𝜁subscript𝜒outsubscript𝐽1delimited-[]~𝜁subscript𝜒outsubscript𝑌1delimited-[]~𝜁subscript𝜒in\displaystyle J_{1}[\tilde{\zeta}(\chi_{\rm in})]Y_{1}[\tilde{\zeta}(\chi_{\rm out% })]-J_{1}[\tilde{\zeta}(\chi_{\rm out})]Y_{1}[\tilde{\zeta}(\chi_{\rm in})]\,,italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) ] italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] - italic_J start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ) ] italic_Y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ over~ start_ARG italic_ζ end_ARG ( italic_χ start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ) ] ,

and ζ~(χ)~𝜁𝜒\tilde{\zeta}(\chi)over~ start_ARG italic_ζ end_ARG ( italic_χ ) is defined by Eq.(40).

In Fig. 4 we present comparisons between the measured functions gω(χ)subscript𝑔𝜔𝜒g_{\omega}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and fω(χ)subscript𝑓𝜔𝜒f_{\omega}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT ( italic_χ ) and the prediction of the theory with dissipation and screening for a number of driving frequencies. The parameters κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG and γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT that were used to fit the analytical solution defined by equations (15),(20),(40), and (IV.2), are typical for our system. The quality of the fits is typical, as long as the frequency of oscillations is not too high, we find very consistent agreement between the analytical theory and the simulations. As can be seen in Appendix B, the transverse components of the displacement field (as well as the nonlinear contributions) are much smaller than the radial ones, and thus we do not refer to these components in this paper.

The dependence of the inverse scale κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG and the dissipative parameter γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT on the frequency of oscillations appears linear. In Fig. 5 we plot an ensemble-averaged values (over ten independent configuration ) of the κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG and γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. The observed linearly is a reflection of the linear regime of the system’s dynamics. It is likely that for higher frequencies, with nonlinear effects becoming relevant, also this observed linearity will be lost.

Our results indicate that as the frequency decreases, κ~~𝜅\tilde{\kappa}over~ start_ARG italic_κ end_ARG goes to zero. This means that the low frequency limit does not coincide with a single inflation of the inner circle [4]. The training of our system by repeated oscillation until a stable steady state is obtained removes the anomalous responses to a single inflation step. The system becomes quasi-elastic and not static-anomalous.

V Conclusions and the road ahead

It is important to stress that consistent and reproducible results for the measured displacement field under oscillatory forcing depend crucially on insisting upon the existence of true stationary state around which the oscillations take place. As explained, the requirement is that stopping the dynamics at any point in time, and applying energy minimization algorithms result in the same equilibrated configuration up to some severe tolerance requirements. This requires repeated training of the system until such a stable steady state configuration is attained. In terms of the energy picture such a stable steady state corresponds to a sufficiently deep local minimum around which the oscillations take place.

A second requirement for the successful correspondence between the analytic theory and the simulations is a synchronization between the oscillatory driving and the system’s dynamics. In the present case, when the driving frequency exceeds 2π2𝜋2\pi2 italic_π in dimensionless units, the system reaction cannot follow the driving, entering a dynamical regime that cannot be treated with the present linear theory. Needless to say, this regime is of interest, but is beyond the scope of the present paper.

Finally, the results of this study indicate very strongly that linear classical elasticity fails to provide a valid description of the response of amorphous solids to oscillatory driving, in much the same way as it fails to describe the response to static driving as has been shown in recent work. It is therefore worthwhile to study in the near future also the effect of screening on oscillatory driving by simple or pure shear, since these modes of forcing have attracted a lot of recent interest due to the presence of shear banding and fracture. Such a future direction will automatically necessitate a nonlinear extension of the screening theory, an effort that is under present active research.

Refer to caption
Figure 6: The azimuthal contribution to the displacement field (a) gωϕ(χ)superscriptsubscript𝑔𝜔italic-ϕ𝜒g_{\omega}^{\phi}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ϕ end_POSTSUPERSCRIPT ( italic_χ ) and (b) fωϕ(χ)superscriptsubscript𝑓𝜔italic-ϕ𝜒f_{\omega}^{\phi}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ϕ end_POSTSUPERSCRIPT ( italic_χ ) at the driving frequency (blue circles) are much smaller that the dominant radial ones gωχ(χ)superscriptsubscript𝑔𝜔𝜒𝜒g_{\omega}^{\chi}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) and fωχ(χ)superscriptsubscript𝑓𝜔𝜒𝜒f_{\omega}^{\chi}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) (green squares). The radial functions are the same is in Fig.4 with ω=6.34𝜔6.34\omega=6.34italic_ω = 6.34.
Refer to caption
Figure 7: The nonlinear contribution to the displacement field (a)g2ωχ(χ)superscriptsubscript𝑔2𝜔𝜒𝜒g_{2\omega}^{\chi}(\chi)italic_g start_POSTSUBSCRIPT 2 italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) and (b) f2ωχ(χ)superscriptsubscript𝑓2𝜔𝜒𝜒f_{2\omega}^{\chi}(\chi)italic_f start_POSTSUBSCRIPT 2 italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) at the double frequency 2ω2𝜔2\omega2 italic_ω (orange triangles) are much smaller that the dominant radial ones gωχ(χ)superscriptsubscript𝑔𝜔𝜒𝜒g_{\omega}^{\chi}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) and fωχ(χ)superscriptsubscript𝑓𝜔𝜒𝜒f_{\omega}^{\chi}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) (blue circles). The radial functions are the same as in Fig.4 with ω=6.34𝜔6.34\omega=6.34italic_ω = 6.34.

Appendix A Hertzian interaction

The disks interact via normal Hertzian forces

F=δ~rarbra+rb[Knδ𝐧ijmeffγn𝐯𝐧].𝐹~𝛿subscript𝑟𝑎subscript𝑟𝑏subscript𝑟𝑎subscript𝑟𝑏delimited-[]subscript𝐾𝑛𝛿subscript𝐧𝑖𝑗subscript𝑚effsubscript𝛾𝑛subscript𝐯𝐧F=\sqrt{\tilde{\delta}}\sqrt{\frac{r_{a}r_{b}}{r_{a}+r_{b}}}[K_{n}\delta{\bf n% }_{ij}-m_{\rm eff}\gamma_{n}\bf{v}_{n}]\,.italic_F = square-root start_ARG over~ start_ARG italic_δ end_ARG end_ARG square-root start_ARG divide start_ARG italic_r start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_r start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT end_ARG start_ARG italic_r start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT end_ARG end_ARG [ italic_K start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT italic_δ bold_n start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT - italic_m start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT bold_v start_POSTSUBSCRIPT bold_n end_POSTSUBSCRIPT ] . (43)

Here δ~~𝛿\tilde{\delta}over~ start_ARG italic_δ end_ARG is the overlap distance of two discs, Knsubscript𝐾𝑛K_{n}italic_K start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT is the elastic constant of the normal contact, γnsubscript𝛾𝑛\gamma_{n}italic_γ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT is the viscoelactic constant for normal contact, and meff=mamb/(ma+mb)subscript𝑚effsubscript𝑚𝑎subscript𝑚𝑏subscript𝑚𝑎subscript𝑚𝑏m_{\rm eff}=m_{a}m_{b}/(m_{a}+m_{b})italic_m start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT = italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT / ( italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT + italic_m start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ) is the effective mass. In our simulations all mj=1subscript𝑚𝑗1m_{j}=1italic_m start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = 1, Kn=2×104,γn=500formulae-sequencesubscript𝐾𝑛2superscript104subscript𝛾𝑛500K_{n}=2\times 10^{4},\gamma_{n}=500italic_K start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT = 2 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT , italic_γ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT = 500. Additional viscous dumping force 𝑭iv=ν~𝒗isubscriptsuperscript𝑭𝑣𝑖~𝜈subscript𝒗𝑖{\bm{F}}^{v}_{i}=-\tilde{\nu}{\bm{v}}_{i}bold_italic_F start_POSTSUPERSCRIPT italic_v end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = - over~ start_ARG italic_ν end_ARG bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is applied to all discs to remove excess of energy from the system. Here ν~=0.9~𝜈0.9\tilde{\nu}=0.9over~ start_ARG italic_ν end_ARG = 0.9 is the damping coefficient and visubscript𝑣𝑖v_{i}italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the velocity of the disc. The simulations are carried out using LAMMPS [14] for the dynamics of the granular system using LJ units.

Appendix B The dominance of the driving mode over nonlinear and azimuthal contributions

Here we present the evidence that the driving mode gωχ,fωχsuperscriptsubscript𝑔𝜔𝜒superscriptsubscript𝑓𝜔𝜒g_{\omega}^{\chi},f_{\omega}^{\chi}italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT , italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT is dominant in the considered range of parameters. Since the original displacement field 𝒅(𝒓)𝒅𝒓{\bm{d}}({\bm{r}})bold_italic_d ( bold_italic_r ) has both radial dr(r)subscript𝑑𝑟𝑟d_{r}(r)italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r ) and azimuthal components defined as 𝒅ϕ(𝒓)𝒅(𝒓)dr(r)r^superscript𝒅italic-ϕ𝒓𝒅𝒓subscript𝑑𝑟𝑟^𝑟{\bm{d}}^{\phi}({\bm{r}})\equiv{\bm{d}}({\bm{r}})-d_{r}(r)\hat{r}bold_italic_d start_POSTSUPERSCRIPT italic_ϕ end_POSTSUPERSCRIPT ( bold_italic_r ) ≡ bold_italic_d ( bold_italic_r ) - italic_d start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( italic_r ) over^ start_ARG italic_r end_ARG. In principle we can solve our equations for both components, but it turns out that the azimuthal component is small compared to the radial one and we did not consider it for the present conditions. In Fig. 6 we compare the radial contributions at the largest driving frequency in our range with the azimuthal contribution at the same frequency, gωϕ(χ)superscriptsubscript𝑔𝜔italic-ϕ𝜒g_{\omega}^{\phi}(\chi)italic_g start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ϕ end_POSTSUPERSCRIPT ( italic_χ ) and fωϕ(χ)superscriptsubscript𝑓𝜔italic-ϕ𝜒f_{\omega}^{\phi}(\chi)italic_f start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ϕ end_POSTSUPERSCRIPT ( italic_χ ). This justifies the present neglect of this component. In addition, we confirm that the radial components of the lowest order nonlinear contribution with double frequency, g2ωχ(χ)superscriptsubscript𝑔2𝜔𝜒𝜒g_{2\omega}^{\chi}(\chi)italic_g start_POSTSUBSCRIPT 2 italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ) and f2ωχ(χ)superscriptsubscript𝑓2𝜔𝜒𝜒f_{2\omega}^{\chi}(\chi)italic_f start_POSTSUBSCRIPT 2 italic_ω end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ end_POSTSUPERSCRIPT ( italic_χ ), are also negligible at our range of frequencies, as is shown in Fig. 7.

Acknowledgements.
We thank Michael Moshe for useful discussions at the early stages of this project. This work has been supported by the joint grant between the Israel Science Foundation and the National Science Foundation of China, and by the Minerva Foundation, Munich, Germany.

References

  • [1] Davide Fiocco, Giuseppe Foffi, and Srikanth Sastry. Oscillatory athermal quasistatic deformation of a model glass. Phys. Rev. E, 88:020301, Aug 2013.
  • [2] Nikolai V. Priezjev. Shear band formation in amorphous materials under oscillatory shear deformation. Metals, 10(3), 2020.
  • [3] Ido Regev, Ido Attia, Karin Dahmen, Srikanth Sastry, and Muhittin Mungan. Topology of the energy landscape of sheared amorphous solids and the irreversibility transition. Phys. Rev. E, 103:062614, Jun 2021.
  • [4] Anaël Lemaître, Chandana Mondal, Michael Moshe, Itamar Procaccia, Saikat Roy, and Keren Screiber-Re’em. Anomalous elasticity and plastic screening in amorphous solids. Phys. Rev. E, 104:024904, Aug 2021.
  • [5] Bhanu Prasad Bhowmik, Michael Moshe, and Itamar Procaccia. Direct measurement of dipoles in anomalous elasticity of amorphous solids. Phys. Rev. E, 105:L043001, Apr 2022.
  • [6] Avanish Kumar, Michael Moshe, Itamar Procaccia, and Murari Singh. Anomalous elasticity in classical glass formers. Phys. Rev. E, 106:015001, Jul 2022.
  • [7] Chandana Mondal, Michael Moshe, Itamar Procaccia, Saikat Roy, Jin Shang, and Jie Zhang. Experimental and numerical verification of anomalous screening theory in granular matter. Chaos, Solitons and Fractals, 164:112609, 2022.
  • [8] Harish Charan, Michael Moshe, and Itamar Procaccia. Anomalous elasticity and emergent dipole screening in three-dimensional amorphous solids. Phys. Rev. E, 107:055005, May 2023.
  • [9] Yuliang Jin, Itamar Procaccia, and Tuhin Samanta. An intermediate phase between jammed and un-jammed amorphous solids, 2023.
  • [10] Chandana Mondal, Michael Moshe, Itamar Procaccia, and Saikat Roy. Dipole screening in pure shear strain protocols of amorphous solids, 2023.
  • [11] Lev Davidovich Landau and Eugin M Lifshitz. Course of Theoretical Physics Vol 7: Theory of Elasticity. Pergamon press, 1959.
  • [12] Sean M Carroll. Spacetime and geometry. Cambridge University Press, 2019.
  • [13] Sebastiano Foti, Carlo G Lai, Glenn J Rix, and Claudio Strobbia. Surface wave methods for near-surface site characterization. CRC press, 2014.
  • [14] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comp. Phys. Comm., 271:108171, 2022.