Velocity statistics inside coherent vortices generated by the inverse cascade of 2d turbulence
aa r X i v : . [ phy s i c s . f l u - dyn ] S e p Velocity statistics inside coherent vortices generated by the inverse cascade of d turbulence I.V.Kolokolov and V.V.Lebedev
Landau Institute for Theoretical Physics, RAS,142432, Ak. Semenova 1-A,Chernogolovka, Moscow region, Russia;NRU Higher School of Economics,101000, Myasnitskaya 20, Moscow, Russia.
We analyze velocity fluctuations inside coherent vortices generated as a result of the inverse cascade in thetwo-dimensional ( d ) turbulence in a finite box. As we demonstrated in [1], the universal velocity profile, es-tablished in [2], corresponds to the passive regime of flow fluctuations. The property enables one to calculatecorrelation functions of the velocity fluctuations in the universal region. We present results of the calcula-tions that demonstrate a non-trivial scaling of the structure function. In addition the calculations reveal stronganisotropy of the structure function. PACS numbers: 68.55.J-, 68.35.Ct, 68.65.-k
I. INTRODUCTION
Effects of the counteraction of (relatively fast) turbulencefluctuations with a coherent (relatively slow) flow is one of thecentral problems of turbulence theory [3]. Usually the fluidenergy is transferred from the slow large-scale flow to turbu-lent pulsations [4]. However, in some cases the energy cango from small-scale fluctuations to the large-scale ones thatcan lead to formation of a non-trivial mean flow [5]. Evenbasic problems such as to determine at which mean velocityturbulent fluctuations are sustained is still object of intenseinvestigations [6]. There is still no consistent theory for themean (coherent) profile coexisting with turbulent fluctuations,so that even the celebrated logarithmic law for the turbulentboundary layer is a subject of controversy [7]. Here, we con-sider an important case: two-dimensional ( d ) turbulence in arestricted box where large-scale coherent structures are gener-ated from small-scale fluctuations excited by pumping. Thisprocess occurs because in two dimensions the non-linear hy-drodynamic interaction favors the energy transfer to largerscales [8–10].Already, the first experiments on d turbulence [11] haveshown that in a finite box with small bottom friction, the en-ergy transfer to large scales leads to the formation of coherentvortices. First numerical simulations [12–14] also show ap-pearing coherent vortices in d turbulence. Subsequent morepronounced numerical simulations [15] and experiments [16]demonstrated that these vortices have well-defined and repro-ducible mean velocity (vorticity) profiles. This profile is quiteisotropic with a power-law radial decay of vorticity inside thecoherent vortex. The profile in that region depends neitheron the boundary conditions (no-slip in experiments, periodicin numerics) nor on the type of forcing (random in numericsversus parametric excitation or electromagnetic forcing in ex-periments). The same flow profile is formed both in the statis-tically stationary case where the mean flow level is stabilizedby the bottom friction and in the case where the average flowis still not stabilized and increases as time runs.In the paper [2] results of intensive simulations of d turbu-lence were reported, they demonstrated that the vortex polarvelocity profile is flat in some interval of distances from the vortex center, that we call the universal interval. The meanvorticity in the interval is inversely proportional to the dis-tance r from the vortex center. In the same paper a theoreticalscheme based on conservation laws and on symmetry argu-ments was proposed that explains the flat velocity profile. Thescheme predicts the value of the polar velocity U = p ǫ/α (where ǫ is the energy production rate and α is the bottomfriction coefficient), that is in excellent agreement with thenumerics [2].In the work [1] we performed an analytical investigationof the coherent vortex in the universal interval. As a result,we established that the flat velocity profile corresponds tothe passive regime of the flow fluctuations where their self-interaction can be neglected. The passive regime admits con-sistent analytical calculations that confirm validity of the value U = p ǫ/α for the polar velocity. Besides, we found ex-pressions for the viscous core radius of the vortex and for theborder of the universal region where the flat velocity profile isrealized. The results reported in the work [1] explain why noflat velocity profile was observed in early simulations [12–14]and imply that at some conditions a large number of coherentvortices could appear instead of a few vortices in numerics[2, 15] and experiment [16].In this paper we examine the spatial structure of the flowfluctuations. The passive nature of the fluctuations admits adetailed analytical analysis. We find the pair correlation func-tions of the velocity fluctuations in the universal interval atscales less than the distance r from the vortex center and largerthan the pumping length. There the correlation function pos-sesses a definite scaling, the scaling is strongly anisotropic.The structure function of the velocity in the range is a linearfunction of the separation between the points. If the dissipa-tion is strong enough, then it can restrict this region of thelinear profile from above. At the end of the paper we discussapplicability conditions of the results and possible extensionsof our scheme. II. GENERAL RELATIONS
We consider the case where d turbulence is excited in afinite box of size L by an external forcing. It is assumed tobe a random quantity with homogeneous in time and spacestatistical properties. We assume also that correlation func-tions of the pumping force are isotropic. The main object ofour investigation is the stationary (in the statistical sense) tur-bulent state caused by such forcing. To excite turbulence theforcing should be stronger than dissipation related both to thebottom friction and the viscosity at the pumping scale. Thatimplies that the characteristic velocity gradient of the fluctua-tions produced by the forcing should be much larger than theflow damping at the pumping scale. The velocity gradient isestimated as ǫ / k / f , where ǫ is the energy flux (energy pro-duction rate per unit mass) and k f is the characteristic wavevector of the pumping force. Thus we arrive at the inequalities ǫ / k / f ≫ α, νk f . (1)Here α is the bottom friction coefficient and ν is the kinematicviscosity coefficient, therefore νk f is the viscous dampingrate at the pumping scale k − f . In simulations, hyperviscosityis often used. In the case the inequalities (1) are still obliga-tory for exciting turbulence, where νk f has to be substitutedby the hyperviscous damping rate at the pumping scale k − f .If the inequalities (1) are satisfied then turbulence is ex-cited in the box and random pulsations of different scalesare formed due to non-linear hydrodynamic interaction. Thepumped energy flows to larger scales whereas the pumped en-strophy flows to smaller scales [8–10]. Thus two cascadesare formed: the energy cascade (inverse cascade) realized atscales larger than the forcing scale k − f and the enstrophy cas-cade realized at scales smaller than the forcing scale k − f . Inan unbound d system the inverse energy cascade is termi-nated by the bottom friction at the scale L α = ǫ / α − / , (2)where a balance between the energy flux ǫ and the bottomfriction is achieved. The enstrophy cascade is terminated byviscosity (or hyperviscosity) [5].In a finite box the above two-cascade picture is realized ifthe box size L is larger than L α . Here we consider the oppo-site case L < L α . Then the energy, transferred by the nonlin-earity to the box size L by the inverse cascade, is accumulatedthere giving rise to a mean (coherent) flow. We analyze thestatistically stationary case where the mean flow is alreadyformed and stabilized by the bottom friction. To describe theflow, we use the Reynolds decomposition, that is the flow ve-locity is presented as the sum V + v where V is the velocityof the coherent flow and v represents velocity fluctuations onthe background of the coherent flow. Let us stress that V is anaverage over time, it possesses a complicated spatial structure.As numerics and experiment show, the coherent flow con-tains some vortices separated by a hyperbolic flow. The char-acteristic velocity V of the coherent motion can be estimatedas V ∼ p ǫ/α . The estimate is a consequence of the energy balance: in the stationary case the incoming energy rate ǫ isequal to the bottom friction rate. The characteristic mean vor-ticity in the hyperbolic region is estimated as Ω ∼ L − p ǫ/α .However, inside the coherent vortices the mean vorticity Ω ismuch higher than the estimate [2, 15, 16]. The maximal valueof the mean vorticity Ω is achieved in the viscous core of thevortex. The radius of the core can be estimated as ( ν/α ) / [1]. III. COHERENT VORTEX
Here we examine the flow inside the coherent vortex. Weattach the origin of our reference system to the vortex centerthat is determined as the point of maximum vorticity. Thedefinition corresponds to the procedures used in the works[2, 15, 16] to establish the mean vortex profile. The posi-tion of the vortex center fluctuates, in the laboratory experi-ments it fluctuates near a fixed position determined by the cellgeometry. For the periodic setup (used in the numerics) thevortex center can shift essentially from its initial position, andonly the average relative position of the vortices is fixed. Thereference system is not inertial, and the velocity of the vor-tex center is subtracted from the flow velocity in the system.However, the flow vorticity in the reference system coincideswith one in the laboratory reference system.As it was established experimentally and numerically [2,15, 16], in the chosen reference system the mean flow pos-sesses the axial symmetry. Such flow can be characterizedby the polar velocity U dependent on the distance r fromthe vortex center. Then the mean vorticity is calculated as Ω = ∂ r U + U/r . To obtain an equation for the profile U ( r ) ,one has to use the complete Navier-Stokes equation. Assum-ing that the average pumping force is zero one finds after av-eraging the Reynolds equation [17]. Outside the viscous corewhere the viscous term is irrelevant we arrive at αU = − (cid:18) ∂ r + 2 r (cid:19) h uv i , (3)where v and u are radial and polar components of the velocityfluctuations and angular brackets mean time averaging.To analyze the flow fluctuations inside the vortex, it is con-venient to start from the equation for the fluctuating vorticity ̟∂ t ̟ +( U/r ) ∂ ϕ ̟ + v∂ r Ω+ ∇ ( v ̟ − h v ̟ i ) = φ − ˆΓ ̟, (4)that is obtained from the same Navier-Stokes equation. Here ϕ is polar angle, φ is curl of the pumping force, v is fluctu-ating velocity, and the operator ˆΓ presents dissipation includ-ing some terms. Among the terms are the bottom friction α and the viscosity term, − ν ∇ . For the case of hyperviscositythe last contribution to ˆΓ is substituted by ( − p +1 ν p ( ∇ ) p where p is an integer number. An additional contribution to ˆΓ is related to the non-linear interaction of the fluctuations.Though the interaction is weak, it could be larger than α and − ν ∇ because of the smallness of the contributions.After solving the equation (4) one can restore the fluctuat-ing velocity from the relation ̟ = ∂ r v + v/r − ∂ ϕ u/r andthe incompressibility condition ∂ r u + u/r + ∂ ϕ v/r = 0 . Thescheme enables one to avoid calculation of the pressure, thatis related to the velocity by a non-local relation. However,at restoring the velocity from the vorticity we still encounternon-local expressions. IV. UNIVERSAL INTERVAL
Further we consider the region outside the vortex corewhere the coherent velocity gradient is large enough,
U/r ≫ ǫ / k / f . (5)In this case fluctuations in the interval of scales between thepumping scale k − f and the radius r are strongly suppressedby the coherent flow. The inequality (5) means that the meanvelocity gradient U/r is larger than the gradient of the ve-locity fluctuations in the region at all scales larger than k − f .Therefore the passive regime is realized there, that is the self-interaction of the velocity fluctuations is weak. The intervalof scales outside the vortex core where the inequality (5) issatisfied will be called further the universal interval of scales..Moreover, the passive regime is realized for scales smallerthan the pumping scale k − f . Indeed, in the direct cascade thevelocity gradients can be estimated as ǫ / k / f , upto logarith-mic factors weakly dependent on scale, see [18–20]. There-fore the inequality (5) means domination of the coherent ve-locity gradient in the interval of scales where the direct cas-cade would realize. The passive regime can be consistentlyanalyzed. Then one neglects the nonlinear in the velocity fluc-tuations term in Eq. (4) staying with a linear equation for thevorticity fluctuation ̟ . The equation enables one to express ̟ in terms of the pumping φ and then to calculate correlationfunctions of ̟ via the correlation functions of φ .Further we focus on the case where the pumping φ is shortcorrelated in time and has Gaussian statistics. Direct calcula-tions [1] show that in this case h uv i = ǫ/ Σ , (6)where Σ is the local shear rate of the coherent flow Σ = r∂ r ( U/r ) = ∂ r U − U/r. (7)The expression (6) is derived at the condition Σ ≫ Γ f , where Γ f is the damping of the velocity fluctuations at the pumpingscale. Validity of the condition is guaranteed by the inequali-ties (1,5). Some additional condition νk f ≫ α is needed forvalidity of the expression (6), the inequality is assumed to besatisfied in our scheme. (Note that the inequality is satisfiedin numerics [2].) The opposite case needs some additionalanalysis that is out of scope of our work.Substituting the expression (6) into Eq. (3) one finds a so-lution U = p ǫ/α, Σ = − U/r, (8) for the mean profile. Thus we arrive at the flat profile of thepolar velocity found in Ref. [2] and confirmed analytically inRef. [1]. It is characteristic of the universal region.The left-hand side of the inequality (5) diminishes as r grows. Therefore it is broken at some r ∼ R u . Substitutingthe expression (8) into Eq. (5) one obtains R u = L / α k − / f = ǫ / α − / k − / f . (9)Note that R u can be larger or smaller than the box size L ,depending on the system parameters. The case R u > L is,probably, characteristic of the numerics [15] and the experi-ments [16], then the passive regime is realized everywhere inthe box. In contrast, in numerics [2] the universal region isrelatively small, R u < L , and is well separated from the outerregion, that is not completely passive. V. VORTICITY FLUCTUATIONS
Since the flow fluctuations inside the universal region arepassive we can use the linearized version of the equation (4) ∂ t ̟ + ( U/r ) ∂ ϕ ̟ + v∂ r Ω + ˆΓ ̟ = φ, (10)Since the pumping is assumed to be short correlated in time,its statistics is determined by the pair correlation function h φ ( t, k ) φ ( t ′ , k ′ ) i = 2(2 π ) ǫδ ( k + k ′ ) δ ( t − t ′ ) k χ ( k ) , (11)for the space Fourier transform of φ . The function χ ( k ) has aprofile with the characteristic pumping wave vector k f and isnormalized: Z d k (2 π ) χ ( k ) = 1 . (12)Then ǫ is the energy (per unit mass, per unit time) pumped tothe system, that is the energy flux.We analyze the fluctuations near a radius r = r with scalesmuch smaller than the radius. Then the shear approximationfor the mean velocity can be used. We pass to the referencesystem rotating with the angular velocity Ω( r ) and expand allterms in the equation (10) in x = r − r and x = r ϕ . Weassume that the parameter ( k f r ) − is small. Then the term v∂ r Ω in Eq. (10) can be discarded and we end up with thefollowing equation: ∂ t ̟ + Σ x ∂ ̟ + ˆΓ ̟ = φ. (13)Let us rewrite the evolution equation (13) for the spatialFourier components of the vorticity ̟ k : ∂ t ̟ ( k ) − Σ k ∂ k ̟ ( k ) + Γ( k ) ̟ ( k ) = φ ( t, k ) . (14)Solving the evolution equation (14) one obtains a formal so-lution ̟ ( t, k ) = Z t dτ φ [ τ, k + Σ( t − τ ) k , k ] (15) × exp − t Z τ dτ ′ Γ (cid:20)q ( k + Σ( t − τ ′ ) k ) + k (cid:21) , Now we can find the simultaneous pair vorticity correlationfunction for the Fourier transform from Eq. (11) h ̟ ( t, k ) ̟ ( t, k ′ ) i = 2(2 π ) ǫδ ( k + k ′ ) ∞ Z dτ q χ ( q ) exp (cid:20) − Z τ dτ ′ Γ( q ′ ) (cid:21) . (16)Here q = ( k + Σ τ k , k ) , (17)and q ′ differs from q by a substitution τ → τ ′ . The factor δ ( k + k ′ ) in the expression (16) reflects space homogeneitythat is not destroyed by a shear flow.Let us consider scales larger than k − f , that is wave vectors k ≪ k f . At the condition the main contribution to the integral(16) is gained from times τ ∼ k f / (Σ k ) . In the case | k | ≫ Γ f k f / Σ the dissipation is irrelevant and the last exponentialfactor in Eq. (16) can be substituted by unity. Here, as above, Γ f is the flow damping at the pumping scale. Passing then tothe integration over the wave vector (17), one obtains h ̟ ( t, k ) ̟ ( t, k ′ ) i = δ ( k + k ′ ) 2(2 π ) ǫq f Σ | k | (18) q f = ∞ Z dq q χ ( q ) . (19)Here, we replaced the lower integration limit | k | in the inte-gral (19) by zero, since the integral is gained at q ∼ k f ≫ k .The wave vector q f is of the order of the inverse pumpinglength.There can exist an interval of the wave vectors r − < | k | < Γ f k f / Σ where the dissipation is relevant. Then theexponential factor in Eq. (16) is relevant. Therefore the ex-pression (18) should be corrected by an additional small factor exp( − A ) , A ∼ Γ f k f / (Σ | k | ) . Thus the vorticity correlationsare strongly suppressed in the region of wave vectors. VI. VELOCITY CORRELATION FUNCTIONS
Knowing the vorticity correlation function, one can calcu-late the velocity correlation functions using the relation v α ( k ) = iǫ αβ k β k ̟ ( k ) , (20)valid for the Fourier transforms. If k f ≫ | k | ≫ Γ f k f / Σ and k f ≫ | k | , then one finds from Eqs. (18,20) h v ( k ) v ( k ′ ) i = 2(2 π ) δ ( k + k ′ ) q f ǫ Σ | k | k , (21) h u ( k ) u ( k ′ ) i = 2(2 π ) δ ( k + k ′ ) q f ǫ Σ k k | k | , (22) h v ( k ) u ( k ′ ) i = − π ) δ ( k + k ′ ) q f ǫ Σ k k k | k | . (23) If r − < | k | < Γ f k f / Σ then the expressions are stronglysuppressed due to dissipation.It follows from the expressions (21,23) that the averages h v i and h u i are determined by the infrared integrals. There-fore h v i , h u i ∼ k f ǫ Σ r if Γ f k f r ≪ Σ , (24) h v i , h u i ∼ ǫ Γ f if Γ f k f r ≫ Σ . (25)The average h uv i needs an additional analysis [1]. It shows,that the quantity is gained at small scales and is determined bythe expression (6).A special problem is calculation of h u i where u is zeroangular harmonics of the fluctuating polar velocity. It is ac-counted by absence in the equation for u an advection termrelated to the average flow. Therefore the quantity h u i is de-termined solely by the damping. Strictly speaking, calculationof h u i is outside of our shear approximation. However, ourlogic can be easily extended to the case to obtain h u i ∼ ǫk f r Γ f (26)An explanation of the expression is based on the expression h u i = Z dϕ π h u ( r ) u ( r ) i , where the points r and r are separated by the same distance r from the vortex center and ϕ is the angle between the vec-tors r and r . The factor ǫ/ Γ f is the contribution to h u i caused by the pumping that is effective if ϕ . ( k f r ) − . Thecontribution (26) should be taken into account besides (24),the latter is related to the sum of non-zero angular harmonics.In the case (25) the contribution (26) is small in comparisonwith one related to non-zero harmonics.The average h u i was calculated previously in the paper[22] where the contribution related to the pumping was ig-nored and the non-linear effects were taken into account in-stead. The approach is correct outside the universal region, at r > R u where R u is determined by the expression (9). At theborder, where r ∼ R u , our estimate (26) coincides with oneobtained in Ref. [22].It is worth to characterize scales where the expressions (21-23) are correct by the velocity structure functions. One finds S ( x , x ) = h [ v ( x , x ) − v (0 , i = 2 q f ǫ Σ π Z d k | k | k (cid:0) − e ik x + ik x (cid:1) , (27)correct if k − f ≪ | x , x | ≪ r, k − f Σ / Γ f . Infrared diver-gence in the integral (27) can be regularized by substituting k → k + µ , where µ ∼ /r . A result of the integration canbe expressed via the function J ( z ) = ∞ Z dq e − z q + µ ≈ π µ + z [Γ f − µz )] . (28)Particularly, one finds S = 2 q f ǫ Σ Re (cid:20) π µ − J + x ∂ J (cid:21) , (29)where J = J ( x − ix ) , Calculating the expression (29), onefinds S ≈ q f ǫ Σ (cid:20) | x | + x arctan (cid:18) x | x | (cid:19)(cid:21) . (30)Analogous expressions can be found for other componentsof the structure function: S = (cid:10) [ u ( x , x ) − u (0 , (cid:11) ≈ q f ǫ Σ (cid:20) x arctan (cid:18) x | x | (cid:19) − | x | ln (cid:18) µ q x + x (cid:19)(cid:21) , (31)and S = h [ v ( x , x ) − v (0 , u ( x , x ) − u (0 , i≈ − q f ǫ Σ x arctan (cid:18) x | x | (cid:19) . (32)In the region | x | , | x | ≫ k − f Σ / Γ f the pair correlation func-tions are strongly suppressed and the structure functions aredominating by the single-point averages. VII. DISCUSSION
We analyzed correlations of the velocity fluctuations insidea coherent vortex generated as a result of the inverse cascadein a finite d cell. Our attention was concentrated on the uni-versal region inside the vertex where the mean velocity has theflat profile. We analyzed the fluctuations on a distance r fromthe vortex core and with scales less than r . The amplitude ofthe velocity fluctuations grows as the scale grow as in the tra-ditional inverse cascade. However, the expressions (30,31,32)demonstrate linear profile, that is different from the / powerlaw in the traditional inverse cascade. Let us stress also that in our case the fluctuations are strongly anisotropic. Notealso that at some conditions viscous dissipation can come intogame, that leads to suppressing the fluctuations at the largestscale (below r ).We performed our calculations in the reference systemwhere origin is attached to the vortex center and rotating withthe angular velocity Ω dependent on the radius r and coincid-ing with angular velocity of the mean flow at the distance r .In this reference system the correlation time of the pumpingattached to the bottom of the cell cannot be larger than Ω − .That justifies our approach (where the pumping is assumedto be short correlated in time) since the angular velocity Ω is the largest characteristic rate in the universal region. Notealso, that use of the rotating reference system implies an im-plicit angular averaging of the correlation functions (besidesthe time averaging).The universal region is restricted from above by the radius(9). At larger distances from the vortex center the flow fluctu-ations are not completely passive, and our scheme is, strictlyspeaking, incorrect. In this case the traditional inverse cascadeis realized on scales smaller than l , where l ∼ ǫ / Σ − / isdetermined by the balance between the effective shear rate Σ of the mean flow and the characteristic velocity gradient in theinverse cascade. However, the flow fluctuations are passiveat scales larger than l . That is the region where our schemeis applicable. And the only difference is that the role of thepumping length is played just by the scale l .Probably, our results can be extended for some types ofthree-dimensional turbulent flows. Note, as an example, theturbulence excited at the fluid surface [23, 24] where the in-verse cascade is observed. It is a subject of future investiga-tions. Acknowledgments
We thank valuable discussions with G. Boffetta and G.Falkovich. The work is supported by RScF grant 14-22-00259. [1] I. V. Kolokolov and V. V. Lebedev, Phys. Rev. E , 033104(2016).[2] J. Laurie, G. Boffetta, G. Falkovich, I. Kolokolov, and V. Lebe-dev, Phys. Rev. Lett., , 254593 (2014).[3] A. A. Townsend, The Structure of Turbulent Shear Flow, Cam-bridge University Press, Cambridge, 1976.[4] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov, Cam-bridge University Press, Cambridge, England, 1995.[5] G. Boffetta and R. E. Ecke, Ann. Rev. Fluid Mech. , 427(2012).[6] K. Avila et al. , Science , 192 (2011).[7] M. Buschmann and M. Gad-el-Hak, AIAA , 4 (2003).[8] R. H. Kraichnan, Phys. Fluids, , 1417 (1967).[9] C. E. Leith, Phys. Fluids, , 671 (1968). [10] G. K. Batchelor, Phys. Fluids, , 233 (1969).[11] J. Sommeria, J. Fluid Mech. , 139 (1986).[12] L. M. Smith and V. Yakhot, Phys. Rev. Lett. , 352 (1993);[13] L. M. Smith and V. Yakhot, J. Fluid Mech. , 115 (1994).[14] V. Borue, Phys. Rev. Lett. , 1475 (1994).[15] M. Chertkov, C. Connaughton, I. Kolokolov, and V. Lebedev,Phys. Rev. Lett , 084501 (2007).[16] H. Xia, M. Shats, and G. Falkovich, Phys. Fluids , 125101(2009).[17] A.S.Monin and A.M.Yaglom, Statistical Fluid Mechanics, vol.1: Mechanics of Turbulence, Edited by J.L.Limley, Dover Pub-lications, Inc., Mineola, New York, 1971, Chapter 3.[18] R. H. Kraichnan, J. Fluid Mech., , 525 (1971); , 155(1975). [19] G. Falkovich and V. Lebedev, Phys. Rev. E , R1800 (1994); , 3883 (1994).[20] G.Falkovich and V.Lebedev, Phys. Rev. E , 045301(R)(2011).[21] I.V. Kolokolov and V.V. Lebedev, Pisma v ZhETF, , 181(2015) [JETP Letters, , 164 (2015)]. [22] G.Falkovich, Proc. R. Soc. A , 20160287 (2016).[23] A. von Kameke, F. Huhn, G. Fernandez-Garsia, A. P. Munuzuri,and V. Perez-Munuzuri, Phys. Rev. Lett. , 074502 (2011).[24] N. Francois, Y. Xia, H. Punzmann, S. Ramsden, and M. Shats,Phys. Rev. X,4