Finite elements computational modeling of coupled elastic waveguides

The theoretical study of one-dimensional-infinite systems of elastically coupled parallel waveguides has established the existence of band structures with pseudo-spin characteristics. Those systems, which are named f-bits, have been shown to exhibit a spinor character associated with directional degrees of freedom, which makes them potential quantum mechanical analogs. The realization of such systems is challenged by the three-dimensional and finite nature of physical elastic waveguides. We address this problem, and with it the design of f-bits in general, by developing finite elements models based on COMSOL Multiphysics®. We model systems of one or more coupled finite length Al rods. The analysis of their dispersion relations, transmission spectra, and amplitudes establishes their f-bit character. For three coupled finite length Al rods, the elastic field is associated with wavefunctions, tensor products of a spinor part related to the directional degrees of freedom, and an orbital angular momentum part representing the phase of the coupled waveguides. We demonstrate the possibility of creating non-separable states between these degrees of freedom. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5127207


I. INTRODUCTION
The investigation of harmonic vibrational properties of materials across spatiotemporal scales is going through a renascence, driven by the similarity of the underlying mathematical structure as compared with that of quantum wave mechanics. 1 A perspective showing that phenomena such as, for example, topological effects, non-reciprocal propagation of waves, coherence of wave amplitudes, and nonseparability are not unique to the quantum realm, [2][3][4][5][6][7] therefore opening the possibility of finding similar applications to those of quantum mechanical nature, but without the drawbacks of the fragile coherence and statistical constraints of the quantum field.
However, perhaps one of the most exciting ideas is to find quantum analogs to those phenomena of quantum coherence and entanglement or nonseparability, driving the foreseen second quantum revolution. 8 This would mean the potential for the development of quantum computing analogs based on classical physics. In fact, coherent waves can be formed in classical physics in general, and the concept of classical nonseparability has been established in classical optics. [3][4][5][9][10][11][12][13][14][15][16][17][18][19] Furthermore, coherence and nonseparability of elastic vibrations and sound waves have been analytically predicted. 1,6 However, the investigation of these potential quantum analog phenomena relies on simplified theoretical models that are unidimensional and infinite in extension, raising questions about the practical realizations of such phenomena, due to their intrinsic three-dimensional and finite nature. The theoretical analysis of feasibility and design of realistic quantum analog devices is then underpinned on the solution of three-dimensional equations, which are usually out of reach of analytical methods, forcing us to resort to numerical procedures.
Pseudo-spin characteristics in acoustic 20,21 and elastic 22-26 systems of topological nature have been investigated. Here, we focus on f-bit systems, 1,27 that is, systems with elastic displacement fields akin to wavefunctions, which exhibit a spinor character associated with directional degrees of freedom. In Ref. 27, we experimentally demonstrated the realization of a f-bit in an acoustic waveguide. In the case of coupled f-bit systems, the elastic field is also described by the eigenvectors of the coupling matrix, which are, in turn, analogs to the Orbital Angular Momentum (OAM) eigenvectors in relativistic quantum mechanics. The elastic spinor and OAMs form a tensor product basis whose linear combinations describe waves in the elastic field and provide the framework for the analysis of nonseparability. 7,28 We describe modeling techniques and their results, which have been successfully used for the exploration and design of f-bit systems, 29 made of elastically coupled parallel waveguides in three dimensions and of finite lengths. The techniques are based on finite elements computational models and have been implemented in COMSOL Multiphysics®. We proceed in Sec. II with the analysis of practical feasibility and the computational design of f-bit systems, starting with a review of the one-dimensional-infinite models defining it and continuing with the introduction of the three-dimensional-finite model of a f-bit based on a single Al rod. Afterward, in Sec. III, we proceed to investigate two f-bit systems, providing first a discussion of the corresponding one-dimensional-infinite models and the notion of nonseparability, followed by the analysis of two f-bit systems, made of three coupled Al rods of finite length, including the investigation of non-separable vibrational modes. Conclusions are drawn in Sec. IV.

A. Simple model
The simplest f-bit 1 is made by an infinite linear chain of masses, coupled by massless springs, with each mass also attached to its own side spring. The other end of the side springs is itself attached to fixed positions or very heavy masses. The equation of motion governing that system is where u n is the displacement of the nth mass in the chain, k 0 is the elastic constant of the coupling springs, and k I is the elastic constant of the side springs. The normal mode solutions are given by with A being the maximum amplitude of the oscillations, a the distance between masses, k the wave vector, and ω the angular frequency. Such solutions lead to the dispersion relation which has a cutoff frequency at k ¼ 0 of ffiffiffiffiffiffiffiffiffiffi k I =m p .

B. Klein-Gordon equation of a waveguide with side elastic coupling
Taking the small k (long-wave) limit for the equation of the discrete f-bit [Eq. (1)], we have the Klein-Gordon (KG) like equation, as the equation of motion for the displacement u(x, t) of a one-dimensional and infinity string along the x axis. The KG equation is used to describe relativistic particles in quantum mechanics.
The normal modes of the KG equation has the form which is straightforward to prove that it yields the dispersion relation proving that, as in the case of the discrete system, the dispersion branch is parabolic with a cutoff frequency at k ¼ 0 of ffiffiffiffiffiffiffiffiffiffi k I =m p . As for the amplitude A, it can take any values in principle, but for the eigensolutions, it is taken as a normalization factor. In any case, A is independent of k and ω.

C. Dirac equation of a waveguide with side elastic coupling
The KG equation [Eq. (4)] admits a Dirac-like factorization into two sets, one for each sign of the α term, of first-order differential equations, where Ψ 2Â1 (x, t) is a two-component vector, each component being a solution of the KG equation. Dirac factorized the KG equation and established the theoretical foundation for the positron. 30 The mathematical formulation of the Dirac equation naturally reveals the spin of relativistic particles. The normal modes of the Dirac system of equations [Eq. (7)] would be which leads to the same dispersion relation as that of the KG equation [Eq. (8)] and to the amplitudes, with s 0 being the normalization constant. Notice that the amplitude vector has components that depend on k and ω in a relative way, and they correspond to the spinor components in Dirac's relativistic quantum mechanics.
D. Three-dimensional one f-bit system

General approach
Real systems are neither infinite nor one-dimensional, but we use the results for the discrete and continuum f-bits to help the exploration of the practical realization of such systems.
The specific systems of interest are single or coupled Al rods, with the coupling provided by epoxy glue and with or without fixed boundary sections. For a three-dimensional and finite system of that nature, we turn to solve numerically the three-dimensional version of the driven elastic wave equations, for homogeneous and isotropic solid objects labeled with m, given by where U m (r, t) is the elastic displacement vector field as a function of the position r at time t for object m (waveguides, coupling, or glue regions), ρ m is the mass density of object m, C m is its elastic constant matrix, f m is the external force acting on it, and η is a damping coefficient (introduced mainly to avoid singularities in the solutions). We need to specify also the boundary conditions of the surfaces of all the objects. Models of the systems of interest were created in COMSOL Multiphysics®, and its mechanical module is used to solve Eq. (10). The Al elastic properties used for the simulations are Young's modulus (E) of 60 GPa, the Poisson ratio (μ) of 0.333, and the density (ρ) of 2660 kg=m 3 . The elastic properties were determined experimentally using ultrasound waves and the density from mass and volume measurements. For the epoxy, we use the values E ¼ 4:08 GPa, μ ¼ 0:378, and ρ ¼ 1142 kg=m 3 . 31 In all cases, tests of convergence of eigen-frequencies with COMSOL Multiphysics® grid size were performed. The COMSOL Multiphysics® grids usually providing convergence are those predefined as extra-fine, which yields a number of grid points in the order of millions.

Eigen solutions
The plane wave ansatz reduces the space of solutions to a specific subset function of wave vectors and frequencies, with the k-vectors defined by boundary conditions. However, given the complexity of the boundaries in the models, it is advisable to resort to solve the three-dimensional equations in the (r, ω) space instead of the (k, ω) space. For that, we assume eigenmodes of the form which leads to the eigenmode equation We solve this equation numerically for given geometries, properties, and boundary conditions via finite elements using COMSOL Multiphysics®.

Driven 3D finite systems
We restrict ourselves to forces of the form where F m is the force applied on the object m at the surface S m d , which is a section of the area at one end of the rods of around 25% of the total area. The solutions in the form of Eq. (11) allow one to write the equations of motion as which are then solved under the same displacement conditions at the surface joining objects and free boundary conditions in the rest of the surfaces that are neither joining objects nor being driven.

Measure of transmission
As a measure of the transmission, we calculate the root mean squared (RMS) of the z-component of the displacement (u z or u 3 ), over an area on the free end faces of the rod of the same magnitude as the driven area S d , to mimic the RMS of the voltage generated by transducers in the experiments. We then generate plots of driving frequencies (ω or f ) vs RMS(u 3 ) to analyze the potential f-bit characteristics. This approach has been implemented in COMSOL Multiphysics® and also the classification of modes described next.

Classification of vibrational modes
It is also useful for the investigation to classify the vibrational modes into longitudinal, transversal, or mixed according to some measure. With that goal, we first define and calculate the velocity V of the elastic field displacements assuming that its shape is constant in time. For that, we choose an arbitrary value of the displacement field, let us say c, and follow the time evolution r c (t) of such isodisplacement through the equation taking the time derivative, we get the equations for the velocity Journal of Applied Physics ARTICLE scitation.org/journal/jap where q and p label the Cartesian components, is the q-component of the displacement velocity at position r, and @up(r, t) @xq is the q component of the gradient of the p-component of the displacement. Equation (17) forms a system of three linear equations and three unknowns, which are the components of the velocity V p (r).
Using Eq. (17), we then define and calculate the RMS of the projection of the velocity on the displacement direction, RMS(cos(θ)), along the longitudinal central line of the rods as RMS(cos(θ)) ; This quantity is finally used to classify the modes according to Notice that in the case of modes with mainly a longitudinal character, the velocity will be parallel to the displacement field along the central line, and therefore, the RMS(cos(θ)) values will be close to one. In the case of transversal modes, the velocity will be perpendicular to the displacement field, and therefore, the RMS(cos(θ)) values will be close to zero. Mixed modes and modes of other nature (e.g., torsion) will give RMS(cos(θ)) values farther from zero or one.

Long rod with free and fixed boundaries
The simplest approximation to a real system close to a f-bit should be a long rod with a fixed boundary condition on all the rod surfaces except at the end faces. We will refer to such a surface as a "cylindrical surface" and to the boundary condition as a "fixed cylindrical surface." The large length compared to the diameter resembles the one-dimensionality, while the fixed boundary condition mimics the effect of the fixed side springs.
For comparison, and validation of the approach, we first calculate the modes of an Al rod of diameter D ¼ 0:01 m and length L ¼ 1 m with a free boundary condition. Figure 1 shows the results for that system. Figure 1(a) provides the classification of the vibrational modes according to the values RMS(cos(θ)) and the response RMS(u 3 ). It shows that the system supports a large variety of vibrational modes all across the range of frequencies, signaling the presence of several dispersion branches. Nevertheless, the agreement between the longitudinal mode classification and the mode shapes shown in Figs. 1(c) and 1(d) points to the validity of the mode classification approach. From the classification of mode results, we plot the frequencies of the longitudinal ones in Fig. 1(b) as a function of k z ¼ n=2L, with n being the number of nodes, that is, the dispersion relation along the k z or axial direction. As expected, it is a straight line with a slope or speed of sound of 5007 m/s, in good agreement with our measured experimental value.
Fixing or immobilizing the entire cylindrical surface of the Al rod, with the same geometrical parameters (D and L) as for the rod with a free boundary condition, leads to the results shown in Fig. 2. That is, as anticipated in Sec. III B, the system shows only the characteristic f-bit like behavior of the dispersion relation [ Fig. 2(b), blue squares]. Also, comparing Figs. 1(a) and 2(a), it is clear that the fixed cylindrical surface induces a separation of the longitudinal and mixed normal modes in frequency and eliminates the purely transversal modes. However, at frequencies higher than 320 kHz, there are mixed modes and modes that may have other characters showing up attesting the three-dimensional nature of the rod.
It is also interesting to explore the behavior of the dispersion relation with changes in the diameter of the rod. Figure 3 shows the result of decreasing the diameter by a factor of two in contrast with increasing the length by the same factor. It shows that the reduction in the diameter shifts the entire dispersion relation by a factor of around two higher in frequency, conserving the shape of the curve. That is, the diameter of the rod also controls the frequency cutoff by modifying the effective value of the stiffness of the side coupling media α. This is an important design element for the f-bits.

Long rod with two fixed boundary lines
Fixing all the cylindrical surfaces of the rod is harder to implement in practice than just fixing lines on the cylindrical surface along the rod axis (Fig. 4). In experiments, the fixed two line condition can be achieved by sandwiching the rod between heavy steel plates and applying pressure on the plates, which will create the fixed boundary condition along the lines of the Hertzian contact between the rods and the plates. 27 More in tune with material dimensions available for experiments, 27 we examine next the vibrational modes and transmission properties of an Al rod of length L ¼ 0:6096 m and diameter D ¼ 12:7 mm, with two lines fixed on the cylindrical surface and opposite to each other across the diameter (Fig. 4). From the analysis of the eigen and driven mode nodes and shapes, we extract the longitudinal f-bit frequencies as a function of the number of nodes. Using those frequencies and knowing that k z ¼ n=2L, n being the number of nodes, due to the free boundary condition at the end of the rod, we reconstructed the f-bit dispersion branch. The main results of such an examination are shown in Fig. 5.
As in the case of a fixed cylindrical boundary, the acoustic linear branch is also converted to a f-bit branch [ Fig. 5(a)] with a cutoff frequency of around 54 kHz, which can also be seen in the transmission spectrum in Fig. 5(b). The mode shapes conserved a longitudinal main character as can be noticed from Figs. 5(c) and 5(d).

III. REALIZATION OF TWO f-BITS USING FINITE ELASTIC COMPONENTS
In Secs. I and II, we have theoretically established, using simple models, that f-bits and coupled f-bits can be realized by fixing boundary lines on single and coupled waveguides. We have numerically verified using a three-dimensional finite model the theoretical results for the single rod as a waveguide. Next, we extend our approach to the design of two f-bit three-dimensional systems taking the form of multiple coupled waveguides.

Equation (4) is the Klein-Gordon (KG) equation for a f-bit,
which we proceed to extend to more than one coupled waveguide, including damping. The system is a set of N parallel unidimensional waveguides coupled along sides to neighbor waveguides by a massless elastic unidimensional media. The equations of motion governing such a system are given in a compact form by with u(x, t) is a vector whose components are the displacement fields along x for each waveguide, β is the speed of sound in the waveguide, η is the damping term divided by the density, α represents the strength of the coupling between the waveguides divided by the density, M is the N Â N coupling matrix, and f (x, t) is the force acting on each rod at position x and time t divided by the density. Equation (20) has the mathematical form of a generalized KG equation, 1 if the damping coefficient η is set to zero. Its normal modes have the usual ansatz

Journal of Applied Physics
ARTICLE scitation.org/journal/jap which reduces Eq. (20) from its differential form to the algebraic eigen-problem for the vector amplitude A, where I is the N Â N unitary matrix. The matrix M is real and symmetric, which means that it has N real eigenvalues λ m , each corresponding to an eigenvectorê m , whose components are all real. Projecting Eq. (22) along the eigenvectorê m and taking into account that it is an eigenvector of M with the eigenvalue λ m , we get which has solutions only if ω takes the values as a function of k, Therefore, the dispersion relations [Eq. (24)] of this system have a quadratic character with a cutoff frequency ffiffiffiffiffi ffi λ m p α if λ m = 0, providing conditions for having multiple f-bit branches.

B. Generalized Dirac equation for multiple coupled waveguides
The generalized KG equation, Eq. (20) with η ¼ 0, admits a Dirac factorization, 1,6 which generates the two equations where I NÂN and I 2NÂ2N are antidiagonal matrices with unit elements, Ψ 2NÂ1 is a 2N dimensional vector that represents the modes of vibration of the N waveguides projected in the two possible directions of propagation along the x-axis, ffiffiffiffiffiffiffiffiffiffiffiffi ffi M NÂN p is the square root of the coupling matrix, and σ 1 and σ 2 are the Pauli matrices

Journal of Applied Physics
ARTICLE scitation.org/journal/jap C. One-dimensional infinite systems with two f-bit branches and nonseparability Systems with more than one f-bit branch open the possibility of realizing acoustic nonseparable states. In the case of three elastic unidimensional parallel waveguides coupled along their lengths, the coupling matrix has the form and choosing the components of Ψ 6Â1 as the plane waves along the basis of the eigenvectors e m,3Â1 of the matrix ffiffiffiffiffiffiffiffiffiffi ffi M 3Â3 p , the amplitudes have the form where s m,2Â1 is given by with s 0 being a complex constant and e m,3Â1 are the eigenvectors of the coupling matrix M, (30) each corresponding to the eigenvalues λ 1 ¼ 0, λ 2 ¼ 1, and λ 3 ¼ 3 and the dispersion relations Equation (28) has the same form as a Dirac spinor of relativistic quantum mechanics, for which the amplitudes are defined by the product of the orbital angular momentum (OAM) eigenvectors and the spinor vector amplitudes. In this case of coupled long waveguides, the OAM part is the eigenvectorsê m of the coupling matrix M. The spinors amount to the two possible amplitudes for the longitudinal elastic waves s m,2Â1 . From there, we name theê m eigenvectors as OAM vectors, or just OAMs, and the s m,2Â1 two-vectors as spinor or spinor-like terms. The elastic spinor components depend on the eigen-frequencies and k vectors and their ratio is fixed, which gives their a coherent character in the sense that if we know one of the spinor components, then we also know the second one.
Equation (25) is linear; hence, a linear combination of solutions of the form given by Eq. (27) is also a solutions. A linear superposition for a given frequency ω ¼ ω 2 ¼ ω 3 of the form will be nonseparable given that the spinor and spatial plane waves cannot be extracted as a common factor, as it is possible in the case where both frequencies and k vectors are the same, yielding a tensor product of a linear combination of the coupling matrix eigenvectors (or the OAM part) with the spinor-like part. So far, we have examined idealized one-dimensional infinite systems. We turn our attention next to the design of real systems exhibiting two f-bit branches and nonseparable states.
D. Three-dimensional two f-bit system

Normal modes, dispersion curves, and mode shapes
The three coupled Al rod system with free boundaries [ Fig. 6(a)] has a large number of eigenmodes and dispersion branches, which make the analysis difficult. To overcome the difficulties, we compare eigenmode shapes and driven-mode shapes, which are generated by driving the three rod system at one of their same side ends, with a harmonic force, along the theoretical OAMs. Observing the eigenand driven-mode shapes, we select those that are mainly longitudinal in character for each OAM branch [ Fig. 6(b)]. In some cases, at low and high k values, it was not possible to separate the mode shapes. Examples of selected mode shapes for each OAM branch are shown in Figs. 6(c)-6(e).
As expected, Fig. 6(b) shows the presence of one linear dispersion branch and two f-bit branches. The example of mode shapes generated by a color contour plot of the module of the elastic field vector clearly exhibits the OAM dependency and their mainly longitudinal character. In the case of theê 1 ¼ (1, 1, 1) FIG. 7. Displacement fields (mode shapes) of the system made of three glued Al rods described in Fig. 6(a) driven according to Eq. (33) as a function of the C 2 constant.

Journal of Applied Physics
branch at f ¼ 7:3 kHz, Fig. 6(e) demonstrates that there are two node regions across the three rods (dark blue) and they are vibrating longitudinally along their axes with the same amplitudes. For theê 1 ¼ (1, 0, À1) branch at f ¼ 21 kHz, Fig. 6(d) depicts the characteristic low excitation expected for the middle rod as well as four nodes going across the rods, attesting the longitudinal character of the oscillations. Similarly, for theê 1 ¼ (1, À2, 1) branch at f ¼ 31 kHz [Fig. 6(c)], we see that the central oscillations are wider as compared with the side rods and the presence of four nodes across the rods.

Observation of nonseparability in a two f-bit system
Let us turn our attention to the equally important question of the realization of a nonseparable linear combination of eigenmode vibrations. We start by choosing a frequency at which the modes for the (1, 0, À1) and (1, À2, 1) OAM branches can be excited, but the (1, 1, 1) cannot. From the dispersion relations in Fig. 6(b), we determined that such a frequency is at 32 960 Hz. Then, we drive the system with an external harmonic force of frequency 32 960 Hz, applied to all three same side ends of the rods. Such a force amplitude in the OAM basis has the form where C 2 and C 3 are constants and the 15.6 factor accounts for differences in the strength of the excitation. Such a factor was found by driving the system first with C 2 ¼ 0 and C 2 ¼ 1 and scaling the (1, 0, À1) OAM until the maximum displacements for both cases become the same. The module of the displacement vector fields for each value of the C 2 constant is shown in Fig. 7. We see how we go from the (1,À2,1) mode shape for C 2 ¼ 0 to the (1,0,À1) mode shape for C 2 ¼ 1. Any of the states in between are non-separable because they correspond to linear combinations of waves having different wave vectors [Eq. (32)] and OAM.
The footprint of the OAM effect is seen in the displacement field relation of the three rods. The side rods displace less than the central one for the (1, À2, 1) OAM (C 2 ¼ 0) and it transitions with increasing C 2 to a state where the central rod does not vibrate, but the side rods do vibrate with the same maximum amplitude of the displacement.
The spinor signature is in the interference pattern created by the forward and backward waves generated by the driving force and the free boundary condition on each rod. The signature is modulated by the OAMs so that even when the spinors themselves are the same for each rod, the interference pattern changes.

IV. CONCLUSIONS
The theoretical investigation of one-dimensional-infinite models has predicted the existence of f-bit systems, which have a parabolic band structure with a finite cutoff frequency in the long wavelength limit. The f-bits usually are described by KG-like equations, which admits a Dirac-like factorization, leading to amplitudes of the displacement field with pseudo-spin characteristics. In systems with more than one f-bit branch, such a dependency allows generating nonseparable states. However, physical or practical realization of such systems relies on the persistence of such properties in three-dimensional and finite systems and that is not guaranteed, in principle. Therefore, to design and investigate the feasibility of practical realization of f-bit systems, we have developed finite elements computational models based on COMSOL Multiphysics®, which allow us to solve the corresponding wave equations numerically as well as to drive the systems as needed.
We have applied the devised approach to extract dispersion relations, displacement fields, and transmission spectra of rods of finite lengths made of Al, either acting alone or coupled to other rods by epoxy glue along their longitudinal directions. We found that single rod systems with fixed cylindrical surfaces, or fixed boundary lines on that surface, have only a dispersion relation branch and it is of a f-bit character. Therefore, in these three-dimensional-finite systems with fixed boundaries, we observe the general theoretical tendency, predicted by one-dimensional-infinite models of single or coupled rods with fixed boundary lines, of behaving like pure f-bit systems.
We also found that nonseparability is a characteristic also present in three-dimensional-finite systems with at least two f-bit branches. More specifically, by driving a system of three coupled Al rods, we demonstrate that we are able to create elastic displacement fields that are nonseparable as a tensor product of an OAM part and a spinor-like term.
Hence, the f-bit characteristics predicted by one-dimensionalinfinite theoretical models persist in three-dimensional-finite systems according to the computational models, which, therefore, confirm the feasibility of the physical realization of f-bit systems and provide insights for their design. The approach and its results have been successfully used for the experimental design of such systems and could be adapted to the investigation of other quantum analog phenomena.