Boundary Reflections of Rolling Waves in Cubic Anisotropic Material

Rolling waves have unconventional circular polarizations enabled by the equal-speed propagation of longitudinal and transverse waves in elastic solids. They can transport non-paraxial intrinsic (i.e. spin) mechanical angular momentum in the media. In this work, we analyze the rolling wave reflections and their effects on the non-paraxial spins in a cubic elastic half-space with an elastically supported boundary. Reflected waves from both normal and general oblique incidences are investigated. We show that, by adjusting the stiffness of the elastic boundary, we can precisely control the spin properties of the reflected waves, paving the way towards a broad category of spin manipulation techniques for bulk elastic waves.


INTRODUCTION
Recent advances in mechanical and acoustic metamaterials have made it possible to manipulate the general intrinsic spin angular momentum (Holanda et al., 2018;Zhu et al., 2018;An et al., 2020;Burns et al., 2020;Rückriegel et al., 2020;Shi and Yang, 2020) carried by elastic waves in the bulk (Long et al., 2018;Wang et al., 2018;Shi et al., 2019;Toftul et al., 2019;Long et al., 2020). Our previous work  demonstrated that, when longitudinal and transverse elastic waves propagate at the same speed, rolling waves may emerge, carrying a spin that is orthogonal to, and in general not parallel to, the wave vector. This is only possible in anisotropic materials that satisfies certain equal-wave-speed criteria . As the simplest anisotropic material characterized by only three independent elastic constants, cubic materials have been investigated in many studies (Stoneley, 1955;Borcherds, 1973). Thomas (1966) analyzed the longitudinal wave velocity at different directions using orientation components, and Mielnicki (1972) studied both longitudinal and transverse waves using approximation method. The bounds of elastic wave speeds (Zuo and Hjelmstad, 1997), surface wave polarizations (Chadwick and Wilson, 1992), and supersonic surface waves (Every, 2015) have all been the topics of an active line of research. Benefited from the fast development in additive manufacturing technology, cubic symmetric metamaterials can now provide us with a broader range of possible elastic constants (Favre et al., 2018;Tancogne-Dejean and Mohr, 2018;Lohmuller et al., 2019) even to the higher-order elasticity effects, such as the compression torsion metamaterial (Frenzel et al., 2017;Wang and Liu, 2020), hyper-stress tensor with second-gradient elasticity (Weeger, 2021), and chiral-micropolar metamaterials (Chen et al., 2021).
Although the stable propagation of rolling waves with non-paraxial spins have been demonstrated , effects of boundary reflections on the wave spin angular momenta have yet to be investigated. The elastic wave reflection process is often both rich in physics and useful in applications (Schoenberg, 1980;Bedford and Drumheller, 1994;Vavrycuk and Psencík, 1998;Lee et al., 2017;Chen et al., 2019). However, studies on reflections from elastically supported boundaries of cubic materials are lacking. This further motivates us to analyze the reflection of rolling waves in a cubic material incident either normally or obliquely at such boundaries. The present work fills the gap in understanding and predicting rolling reflections and offers a new technique to manipulate the spin of elastic bulk waves.
In this paper, we first review the general governing equations for plane waves in an anisotropic medium. Next, the equal-wavespeed criteria in cubic materials are derived. Furthermore, we provide a general formulation for the spin angular momentum of bulk elastic waves. Lastly, the normal and oblique incident rolling waves are analyzed with different boundary conditions, where both spin-preserving and spin-flipping phenomena may emerge.

PLANE WAVES IN ANISOTROPIC MEDIUM
With consideration of long-wavelength limit (i.e. quasi-static limit) only, we focus on the plane wave with wavelength much larger than any micro-structure size, which implies the dispersive behavior from geometry is negligible. The general bulk elastic wave equation without body force is, where σ is the stress tensor and u is the displacement vector. The general solution of plane waves is where u 0 denotes wave amplitude, and k the wave vector. We consider a plane wave propagating along an arbitrary direction specified by the unit vector, k̂ l 1 e 1 + l 2 e 2 + l 3 e 3 , where l 1 , l 2 , l 3 are direction cosines, and e 1 , e 2 , e 3 are the standard orthonormal basis vectors of Cartesian coordinates. Then, we have where k is the wavenumber. The time derivative of eq. 2 gives us, With the definition of a matrix L, L l 1 0 0 0 l 3 l 2 0 l 2 0 l 3 0 l 1 0 0 l 3 l 2 l 1 0 We have the spatial derivative in Voigt notation (Carcione, 2007), Combining the above equations, we arrive at where is the Kelvin-Christoffel matrix.

CONSTITUTIVE RELATION OF CUBIC MATERIALS
Next, we focus our discussions on cubic materials characterized by the cubic crystallographic point groups. As illustrated in Figure 1, the symmetry conditions include four axes of threefold rotational symmetry that can be identified with the body diagonals of a cubic unit cell (Bückmann et al., 2014;Dirrenberger et al., 2013;Authier, 2003). Rather counterintuitively, this symmetry requirement does not include any axis of four-fold (90 degree) rotational symmetry (Ashcroft and Mermin, 1976, p. 121). The stiffness matrix of cubic material is given as (Authier, 2003;Norris, 2006;Zhang et al., 2020), For waves propagating along the z-direction, the direction cosines are l 1 0, l 2 0, l 3 1.
Then, the Kelvin-Christoffel matrix components become, Therefore, the governing eq. 11 becomes, We note that u 3 represents the longitudinal wave, while u 1 and u 2 represent the transverse waves. Therefore, from eq. 16, we obtain the longitudinal and transverse wave speeds as, Thus, for waves traveling in the z-direction of a cubic material, in order to have a rolling wave with non-praxial spin polarization , we need the equal-speed criterion: By symmetry, for the waves propagating in x− or y−direction, the equal-speed criterion is exactly the same as eq. 18.

INTRINSIC SPIN ANGULAR MOMENTUM OF BULK ELASTIC WAVES
To calculate the spin density of the bulk elastic waves, we consider the displacement field of a general plane wave u Importantly, here m, n and l are complex-valued, so that they contain the information not only about amplitudes but also about the relative phase differences among displacement components of the traveling wave. The spin angular momentum density, as a real-valued vector, can be calculated as (Berry, 2009;Long et al., 2018;Burns et al., 2020): where (·)* denotes complex conjugation, and the spin-1 operator is a third-order tensor defined as Hence, the spin density for a general traveling wave is We note that the above formulation makes physical sense only if all displacement components propagate at the same wave speed. Particularly for cubic materials, the equal-speed criterion eq. 18 needs to be true. In this case, we can have stable propagation of spins pointing at any direction by adjusting the complex amplitudes m, n and l. This can be realized by simply changing the relative phase differences. Hence, the direction of the spin vector s can be completely independent from the direction of the wave vector k.

NORMAL INCIDENCE AND REFLECTION
We now investigate normal reflections of a rolling wave at a general elastic boundary. This serves as an example non-paraxial spin manipulation. Considering a rolling wave along the zdirection normally incident on a flat surface, we have the instantaneous wave displacement fields at time t 0 as u I ũ I exp(ikz) and u R ũ R exp(−ikz) with where the superscripts, I and R, denote the incident and reflected waves, respectively. For the displacement with small amplitude, the linear strain terms are dominant. In this case, by neglecting higher-order strain terms, the strains can be calculated from displacements by where the comma ",″ in subscripts denotes the derivative operation. The reflection occurs at z 0. So we have e ikz e −ikz 1, and the strain vector in Voigt notation becomes The cubic constitutive relations are It is easy to compute the stress components for incident and reflection waves, respectively: For a elastically supported cubic half-space, the boundary conditions are where (K x , K y , K z ) are components of the distributed stiffness per unit area (Zhang et al., 2017) representing a general elastic foundation supporting the solid surface, and the total stress and displacement are u 0 Substituting the displacement and stress components into elastic boundary conditions, we get Solving the above system of equations, we obtain Because of the equal-speed requirement of rolling wave (C 11 C 44 ), the amplitude of reflection wave becomes, with the generalized amplitude ratio of normal reflection given as We can also write in the form of a reflection matrix (Zoeppritz, 1919;Ursin and Tjäland, 1996): with the reflection matrix R̂ diag(R x , R y , R z ). This complexvalued non-dimensional parameter R j is the key connection between the incident and reflected waves and deserves further analysis. We note that |R j | 1 is consistent with the fact that all wave energy is reflected. Hence, we have Here, the phase angle ϕ of R j corresponds to the phase change during the reflection, and it dependents on the boundary-bulk stiffness ratio K j /C 11 and the wave number k, as shown in Figure 2. Therefore, by adjusting the elastic stiffness at the boundary, one can manipulate the spin of the reflected waves.
Although all properties of both the bulk and the reflection surface are assumed to be independent of the incidence wave frequency, here the reflection phase change can be frequencydependent, as the angular wave number k appears in the ratio. The emergence of frequency dependency can be intuitively FIGURE 2 | The complex-valued amplitude ratio of normal reflection, R j , is determined by the boundary-bulk stiffness ratio, K j /C 11 , and wave number k.
Frontiers in Mechanical Engineering | www.frontiersin.org August 2021 | Volume 7 | Article 704192 explained by the role of wavelength during reflection. We note from eq. 25 that k first appears in the strain calculations since, for a given wave displacement amplitude, the strains in the propagation direction, ε zj , are actually inversely proportional to the wavelength: Longer wavelength gives rise to a smaller strain and vice versa. Consequently, the stresses, σ zj , and the force acting on the boundary springs are wavelength-dependent as well.
If the boundary is supported by an elastic foundation with finite stiffness per area, K j , we have the following: At the low-frequency and long-wavelength limit, the force per area acting on the boundary, |σ zj |∝ C 11 k → 0. So, the boundary hardly moves, and the incidence wave effectively "sees" a rigid surface; At the high-frequency and short-wavelength limit, the force per area acting on the boundary, |σ zj |∝ C 11 k → ∞. So, the boundary moves a lot, and the incidence wave effectively "sees" a free surface. Next, we consider some special cases and focus on the in-xzplane rolling waves with n I n R 0. We note that the elastic boundary condition degenerates into traction-free boundary condition (Neumann type) when the boundary stiffness K j → 0. Then, the amplitudes of reflected wave become Similarly, the elastic boundary condition degenerates into the fixed rigid boundary condition (Dirichlet type) when the stiffness K j → ∞. Then, the amplitudes of reflected wave become In both cases above, by the definition of spin given in eq. 22, we can conclude s R s I , as shown in Figures 3B,C. So, the spin is unaffected by the reflection process. In addition, for the fixed-free hybrid boundary (K x → ∞ and K z → 0, Figure 3D), we have Similarly, for the free-fixed hybrid boundary (K x → 0 and K z → ∞, Figure 3E), we have Thus, both hybrid boundaries will flip the spin for any incident rolling wave.
Summarizing different cases for the normal incidence, we note that, for K x → 0 and K z → 0, the boundary becomes traction-free (Neumann type) and we haveũ R ũ I with no phase change. For K x → ∞ and K z → ∞, the boundary becomes rigid (Dirichlet type) and we have an out-of-phase reflected wave withũ R −ũ I .
In both cases, we have s R s I , so the spin is preserved. In contrast, for hybrid boundaries (K x → ∞, K z → 0) and (K x → 0, K z → ∞), we get s R −s I , so the spin is flipped due to the difference in phase change between the longitudinal and transverse displacement components during the reflection process. Figure 3 illustrates the results using the example of an incident rolling wave carrying a non-paraxial spin of s I y −1. The reflection process is spin-preserving in both rigid and free boundaries, while being spin-flipping for both hybrid freefixed and hybrid fixed-free boundary conditions.

THE OBLIQUE INCIDENCE AND REFLECTION AT 45 DEGREE
Next we investigate the scenario in which the boundary plane is oblique to the incident rolling wave at the 45°angle. We focus on the incident rolling wave propagating towards the positive z−direction, As shown in Figure 4A, the incidence angle here is α I 45°. Importantly, we note that, due to material anisotropy, the FIGURE 3 | (A) Rolling wave reflections with normal incidence at a general elastic boundary. In all sub-plots, the incidence wave carries a spin angular momentum, s I y −1. Here s I and s R indicate the spin density of incidence wave and reflection wave, respectively. They can be calculated using the eq. 22. The reflection process is (B) spin-preserving at the fixed rigid boundary, (C) spin-preserving at the free boundary, (D) spin-flipping at the fixed-free boundary, and (E) spin-flipping at the free-fixed boundary.
Frontiers in Mechanical Engineering | www.frontiersin.org August 2021 | Volume 7 | Article 704192 reflection angle(s) will not necessarily be equal to the incidence angle. In fact, a single incidence wave, in general, may result in multiple reflected waves towards different directions (Achenbach, 2012;Graff, 2012). Hence, the reflection angle(s) need to be determined by a detailed analysis on boundary conditions and material properties. In this case, we can write the reflection wave as a general form in xz-plane, where θ is the angle of wave vector with respect to the x-axis. It follows that the reflection angle can be expressed as α R θ + 45°, θ ∈ (−45°, + 45°). Next, we look into the general conditions for rolling waves propagating in the xz-plane. The normalized wave vector, as given in eq. 3, is now characterized by cosine directions l 2 0 and l 2 1 + l 2 3 1. The element expressions of matrix Γ become, Γ 11 C 11 l 2 1 + C 44 l 2 3 , Γ 22 C 44 , Γ 33 C 44 l 2 1 + C 11 l 2 3 , Γ 12 0, Γ 13 (C 12 + C 44 )l 1 l 3 , Γ 23 0.
For any linear material to be statically stable, it needs a positive strain energy function, and hence a positive-definite stiffness matrix (Ting, 1996), which leads to the following conditions: C 11 > 0 and − C 11 < C 12 < C 11 (52) Then, in order to determine the reflection angle(s) in this case, we apply Snell's law for reflection Here, the incidence angle is α I 45°. Combining it with ω I ω R and v I C 11 /ρ , we arrive at 1 2C 11 √ sin(α R ) C 11 ± |C 12 + C 11 |l 1 l 3 sin(θ + 45°) C 11 ± |C 12 + C 11 |l 1 l 3 with the cosine directions, l 1 cos(θ), l 3 sin(θ).
For θ ∈ (−45°, + 45°) with the positive-definite conditions of eq. 52, the only solution to above equation is θ 0 (57) which implies α R 45°. Therefore, there can be one reflection wave vector only, and the reflected wave will propagate along x-direction at the speed v R v I . Hence, we can rewrite the reflection wave as, For brevity, denoting the wave number as k I k R k, we obtain The corresponding strain wave amplitudes are And the stress wave amplitudes are To analyze the normal and transverse elastic boundary condition, we introduce a local coordinate system x′ Ψx specified by the transformation matrix, where the φ 45°as shown in Figure 4. The general elastic boundary conditions in the local coordinate system are where (K x′ , K y′ , K z′ ) are components of distributed stiffness per unit area, and are components of total stress and total displacement on the boundary surface. We next perform coordinate transforms using In the local coordinates, x′, we obtain σ′ I ik 2 Similarly, the local stress and displacement of reflection wave are, Substituting the stress and displacement components eqs 67 -eqs 70 into the boundary conditions eqs 64,65, we obtain ik 2 [(C 11 + C 12 )l I − 2C 11 m I + (C 11 + C 12 )m R − 2C 11 l R ] Then, we can solve the system of eqs 71-73 for the complexvalued amplitude components of the reflected wave: where κ + ik 2 √ (C 11 + C 12 ), κ − ik 2 √ (C 11 − C 12 ).
We can also write in the form of a reflection matrix (Zoeppritz, 1919;Ursin and Tjäland, 1996).

FIGURE 4 | (A)
Illustration of rolling wave reflections at 45°incidence. In all sub-plots, the incidence wave carries a spin angular momentum, s I y +1. Here s I and s R indicate the spin density of incidence wave and reflection wave, respectively. They can be calculated using the eq. 22. The reflection process is (B) spin-preserving at the fixed rigid boundary, (C) spin-preserving at the free boundary, (D) spin-flipping at the fixed-free boundary, and (E) spin-flipping at the free-fixed boundary. Note that the free boundary in (C) has the specially effect of giving rise to an additional linearly polarized reflected wave. Although the spin is still exactly preserved, the overall reflected wave has an elliptical rolling polarization. where with denominator D K x′ (K z′ + κ + ) + (K z′ + 2C 11 ik 2 √ )(κ − + K x′ ). We note that R ml ≠ R lm , which means that the complex-valued reflection matrix is not Hermitian. A closer scrutiny reveals that the amplitudes of R ml and R lm are always the same. The two complexvalued reflective coefficients only differ in their phase angles.
We also note that the term R nn for the out-of-xz-plane wave is independent and very similar to the normal incidence scenario, as shown in eqs 38,40 and Figure 2. All other terms in eq. 79 are illustrated in Figure 5 to show how they change with the wave number, k, and the boundary stiffness ratio, K z′ /K x′ . We observe in Figure 5B that both R ml and R lm shows a phase jump of 180°. This is because of two reasons: 1) The small wave number k 0.001 leads to the fact that both R ml and R lm are almost purely imaginary; and 2) At the boundary stiffness ratio of K z′ /K x′ 3, values transition from a negative imaginary number to a positive imaginary number. This results in the 180°phase jump.
Clearly, the out-of-xz-plane wave n R only depends on n I : n R −n I for K y′ ∞.
The four extreme cases for in-xz-plane wave polarizations are: Case I ( Figure 4B), K x′ → ∞, K z′ → ∞: Case II ( Figure 4C), K x′ → 0, K z′ → 0: with C r C 12 /C 11 . Here, the second term in l R is associated with l I , and m R l I . So, this term does not have any contribution to the circularly rolling polarization, and it gives rise to an additional linearly polarized wave. Consequently, the spin is still exactly preserved, but the overall reflected wave has a elliptically rolling polarization.
Case III ( Figure 4D), K x′ → ∞, K z′ → 0: Case IV ( Figure 4E), K x′ → 0, K z′ → ∞: As illustrated in Figure 4C, the results in eq. 83 for Case II are special and remarkable. The reflected wave from the free boundary is different from all other cases. However, the reflected spin density is indeed still the same to that of the incidence wave. The reason is that the additional term in l R does not contribute to the circularly rolling polarization. To further investigate this especially interesting case, we calculate the polarization of reflected waves with different ratios of C r C 12 /C 11 and plot the result in Figure 6. Note that the ratios of C r −1 in Figure 6A and C r 1 in Figure 6E are actually extreme scenarios not attainable in statically stable materials due the positive definite conditions listed in eq. 52.

THE OBLIQUE INCIDENCE AND REFLECTION AT ARBITRARY ANGLE
Next, we look into the situation with an arbitrary incidence angle. From Snell's law eq. 53, we can obtain sin(α I ) C 11 √ sin(α R ) C 11 ± (C 12 + C 11 )l 1 l 3 where l 1 cos(θ), l 3 sin(θ). Noticed from Figure 4A that the incidence angle is identical to the rotation angle, i.e., α I φ and FIGURE 6 | In the oblique reflection at 45°, elliptically rolling polarizations of the reflected waves from the free boundary may appear due to different material stiffness constants C 11 and C 12 . An incident spin which is the same of that shown (A) is assumed throughout. (A) also represents the purely circularly rolling polarization for C r −1. (B-F) show the elliptically rolling polarizations for the cases of C r −0.5, C r 0, C r 0.5 and C r 1, respectively. Note that all polarizations are scaled for the purpose of clearly demonstrating differences in elliptical eccentricity.
then the reflection angle α R π/2 − φ + θ. The reflection angle α R is the solution of the transcendental equation depending on the incidence angle α I and the material constant ratio C r C 12 /C 11 : where "±" gives rise to two different reflected wave modes. In general, those two modes have two different reflection angles. Next, we numerically solve eq. 87 for the reflection angle α R and illustrate the results in Figure 7. The points P 1 and P 2 in the plots indicate the cases of the normal incidence and the 45°i ncidence, respectively. Both cases have been thoroughly discussed in previous sections. The next critical point P 3 is where the reflection angle α R → π/2 for Mode 1. Plugging α R π/2 into eq. 87, we obtain: (1 + C r ) 2 sin 6 (α I 3 ) + sin 2 (α I 3 ) − 1 0.
Thus, we arrive at the analytical solution of the above equation, We know that, with Δ > 0, the cubic function always has one real root and two complex conjugate roots. In this case, we can only obtain one critical angle α I 3 . Therefore, we conclude that, for incidence angle larger than the critical angle α I 3 , the reflection wave of Mode 1 will no longer exist.
Here, from eq. 52 we have the constraint −1 < C r C 12 /C 11 < 1. Hence, we investigate all values of C r within this range and find the following four representative scenarios for Mode 2 of the reflected waves: • For 0 < C r < 1, as shown in Figures 7A,B, we have one local maximum (point P 4 ) and no local minimum for Mode 2. The reflection angle α R will increase together with the incidence angle α I until it reaches to the maximum at point P 4 . Then, the anomalous trend appears, and the reflection angle α R will decrease if we further increase the incidence angle α I . Consequently, for Mode 2, we may have two different incidence angles with the same reflection angle. • For −0.149 < C r ≤ 0, as shown in Figures 7C,D, we have one local maximum (point P 4 ) and one local minimum (point P 5 ). Remarkably, as shown in the inset of Figure 7D, we may have three different incidence angles with the same reflection angle for Mode 2. Also, we note that, with the decreasing value of C r , point P 4 and point P 5 move closer to each other. • For C r −0.149, as shown in Figure 7E, we have no local maximum, no local minimum, but a saddle point (point P 6 ). This saddle point P 6 is the result of point P 4 and point P 5 merging together. • For −1 < C r < −0.149452 as shown in Figure 7F, we have a monotonic curve for Mode 2. The reflection angle α R always increases with the increasing incidence angle α I . And this is the normal phenomenon generally understood.
In the first two scenarios listed above, we have the anomalous trend phenomenon: The reflection angle α R does not have a monotonic relationship with the incidence angle α I . This may provide us with a new technique to control the direction of reflected waves.
Furthermore, we can take the derivative of both sides of eq. 87 with respect to the incidence angle α I and obtain the changing rate of reflection angle α R of Mode 2, dα R dα I 1 − 1+Cr 2 sin(2α R + 2α I ) sin(2α I ) − (1 + C r )sin 2 (α I ) cos(2α R + 2α I ) sin(2α R ) + (1 + C r ) sin 2 (α I )cos(2α R + 2α I ) We note that this derivative of α R with respect to the α I is the root of another transcendental equation and rarely has closed-form solution. Thus, we solve it numerically and plot the results in Figure 8. The same four represented scenarios discussed above can be easily categorized here as different number of roots of

CONCLUSION
Cubic materials are anisotropic elastic materials that satisfy the symmetry requirements illustrated in Figure 1, and their stiffness matrix is given by eq. 13. If the equal-speed criteria for longitudinal and transverse waves, as given in eq. 18, also holds in a cubic material, then we can have stable propagation of any spin angular momentum carried by bulk elastic waves. In this paper, we focus on rolling waves, which have a spin vector orthogonal to its wave vector, in cubic materials. We investigate the reflections of rolling waves incident either normally or obliquely at an elastically supported flat boundary. Behaviors of spin-preserving, spin-flipping, and general arbitrary phase shifts are shown by adjusting the boundary conditions. This work provides a new approach to manipulating bulk elastic wave spins. Several other key findings are summarized below: 1) For the normal incidence of rolling waves, the reflection is spin-preserving at both the fixed rigid boundary and the free boundary, while being spin-flipping at both the fixed-free and free-fixed boundaries. Most generally, we can even introduce the arbitrary phase shifts between reflection and incidence waves by adjusting the elastic boundary conditions. 2) For oblique incidence at 45°, we prove that, contrary to the general anisotropic case, there will be only 1 reflected wave with a unique reflection angle of 45°. 3) With the 45°incidence, the reflection is spin-preserving at a fixed rigid boundary, while being spin-flipping at both fixed-free and free-fixed boundaries. The free boundary is a uniquely special case here: While the reflection is still spin-preserving, depending on the cubic material constant ratio, C 12 /C 11 , it will also introduce an additional linearly polarized part. Hence, an incidence wave with circularly rolling polarization will result in a reflected wave with elliptically rolling polarization. 4) For an arbitrary incidence angle, we showed that each incidence wave will in general result in two reflected wave modes with different reflection angles. One of the mode (Mode 1 in Figure 7) exists only when the incidence angle is less than a critical value α I 3 . Also, The other mode (Mode 2 in Figure 7) may show anomalous trend. 5) We developed a new way to manipulate the elastic wave spin.
Further, the present study can be extended to more general cases, i.e., arbitrary incidence angle for reflection and refraction problem, the reflection matrix properties of arbitrary incidence angle, visco-elastically supported boundary condition, or breaking the limit of statically stable materials by using metamaterials.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
PZ is the primary contributor to this study and this manuscript. PW is the secondary contributor to this study and this manuscript.

FUNDING
This work was supported by start-up funds of the Dept. Mechanical Engineering at Univ. Utah.