Navigating the Hilbert space of nonseparable elastic states in arrays of periodically coupled one-dimensional waveguides

A planar array of three one-dimensional elastic waveguides mutually coupled periodically along their length and driven externally is shown theoretically and numerically to support nonseparable superpositions of states. These states are the product of Bloch waves describing the elastic displacement along the waveguides and spatial modes representing the displacement across the array of waveguides. For a system composed of finite length waveguides, the frequency, relative amplitude, and phase of the external drivers can be employed to selectively excite specific groups of discrete product modes. The periodicity of the coupling is used to fold bands enabling superpositions of states that span the complete Hilbert space of product states. We show that we can realize a transformation from one type of nonseparable superposition to another one that is analogous to a nontrivial quantum gate. This transformation is also interpreted as the complex conjugation operator in the space of the complex amplitudes of individual waveguides. © 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/5.0014259., s


I. INTRODUCTION
Superposition of states and entanglement are at the core of today's second quantum revolution. 1 Entangled superpositions of states possess two distinct attributes, namely, nonlocality and nonseparability. Nonlocality is a unique feature of quantum mechanics, but nonseparability is not. The notion of "classical entanglement," i.e., local nonseparable superposition of states, has received a significant amount of attention, theoretically and experimentally, in the area of optics. [2][3][4][5][6][7][8] However, to date, much less attention has been paid to other classical waves such as elastic waves; yet, remarkable new behaviors of sound, analogous to quantum physics, such as the notions of elastic pseudospin [9][10][11][12][13][14][15][16][17][18][19] and Zak/Berry phase, [20][21][22][23][24][25][26][27][28] are emerging. Recently, elastic nonseparable states, 29 analogous to "classically entangled" states, were observed in externally driven systems composed of parallel arrays of one-dimensional elastic waveguides coupled elastically and uniformly along their length. [30][31][32] In these theoretical and experimental studies, the classically nonseparable states are Bell states constructed as a superposition of elastic waves, each being a product of a plane wave part and a spatial eigen mode. The plane wave part describes the elastic wave along the length of the waveguides, and the spatial eigen mode characterizes the amplitude and phase difference between waveguides, i.e., across the array of waveguides. In the case of a finite number of finite length waveguides, the plane wave and spatial parts are both discrete. It is therefore possible to produce elastic Bell states, whose amplitude can be represented using the usual ket notation of quantum mechanics in the form A|S 1 ⟩|k 1 ⟩ + D|S 2 ⟩|k 2 ⟩, where |Si⟩ are discrete spatial states and |kj⟩ are discrete plane wave states. This superposition that cannot be factored as a single product is thus nonseparable. When the constitutive materials are dissipative, A and D are complex. By tuning A and D through the amplitude, phase, and/or frequency of the external drivers, one is able to experimentally navigate a portion of that Bell state Hilbert space. 30,31 However, in these previous studies, the system did not allow varying the spatial modes and plane wave states independently. This system could not be used to explore the complete Hilbert space of product elastic states because it does not support the general superposition A|S 1 ⟩|k 1 ⟩ + B|S 1 ⟩|k 2 ⟩ + C|S 2 ⟩|k 1 ⟩ + D|S 2 ⟩|k 2 ⟩. Manipulation of such a general state is highly desirable as it would allow for the creation of Bell states of the form B|S 1 ⟩|k 2 ⟩ + C|S 2 ⟩|k 1 ⟩ by tuning the complex coefficients A and D to zero. This limitation inherent to the previously studied system can be overcome by considering an externally driven system composed of a parallel array of one-dimensional elastic waveguides coupled elastically and nonuniformly along their length. In the case of waveguides coupled periodically along their length, elastic states are now products of a spatial part and a Bloch wave part. It is the objective of the present study to show that band folding within a finite Brillouin zone due to periodic coupling between waveguides offers the possibility of selecting spatial and Bloch wave states independently. Using periodically coupled one-dimensional waveguides, one can therefore navigate a more sizeable portion of the Hilbert space of nonseparable spatial mode/Bloch wave product states. Navigation of the Hilbert space of product states is achieved by stimulating externally each waveguide in the array with a driver. The product states can be tuned by varying the relative amplitude and phase of each driver as well as the frequency of the drivers.
Sections II-V introduce the theory of externally driven arrays of one-dimensional waveguides coupled periodically along their length. We present the mathematical formulation behind achieving nonseparable superpositions of states, which are the product of plane waves and spatial modes. We also show how these states can be manipulated by varying the drivers' characteristics. We demonstrate the possibility of achieving transformations between different nonseparable superpositions of product states reminiscent of nontrivial quantum gates. In Sec. VI, we introduce a mass-spring model of a parallel array of discrete waveguides coupled elastically and periodically along their length. This model is used to validate and illustrate the theory. We calculate the band structure of this model. Periodicity-induced band folding and finite length of the waveguides lead to modes that can be excited externally to realize superpositions of states that are the product of plane waves and spatial modes that span the complete available Hilbert space of product states. A section of that Hilbert space is explored by tuning the frequency and the relative amplitudes of the external drivers stimulating the individual waveguides. We show that transformation from one type of nonseparable superposition to another one in that space can be interpreted as the complex conjugation operator in the space of the complex amplitudes of individual waveguides.

II. ARRAY OF PARALLEL WAVEGUIDES COUPLED PERIODICALLY ALONG THEIR LENGTH
We consider three one-dimensional elastic waveguides coupled elastically along their length (see Fig. 1). In the long wavelength limit, the wave equation takes the form The parameter β is proportional to the speed of sound in the waveguides.
x represents the position along the waveguides. I 3×3 is the identity matrix. M 3×3 is the matrix describing the elastic coupling between the three waveguides. The 3 × 1 vector U = (U 1 , U 2 , U 3 ) represents the displacement in waveguides 1, 2, and 3, respectively. We now assume that the coupling is periodic along the direction x. For instance, in the case of a planar array of three waveguides, the coupling matrix takes the form where and α(x) is a periodic function with the period P. The parameter α is proportional to the stiffness of the coupling. We write the function α(x) as a Fourier series with the reciprocal lattice number G = m 2π P , where m is an integer. We seek solutions in the form of Bloch waves, namely, where g is also a reciprocal lattice number. Inserting Eqs.
(2)-(4) into (1) yields We multiply Eq. (5) by e −i(k ′ +g ′ )x and integrate over one period In obtaining Eq. (6), we have used the relation ∫ dxe i(a−b)x = δ a,b , where δ a,b is the usual delta function.
We now define λn and En with n = 1, 2, 3 as the eigen values and eigen vectors of the matrix N 3×3 . That is, we can write N 3×3 En = λnI 3×3 En. These eigen vectors are equivalent to spatial states across the parallel array of waveguides. Since En form a complete orthonormal basis, we can write the displacement vector as FIG. 1. Schematic representation of three one-dimensional elastic waveguides coupled periodically along their length with a period P (stiffness of the coupling is represented schematically by different shades of gray).
Using Eq. (7), Eq. (6) becomes (here, we have dropped the prime for convenience) By factoring I 3×3 En, Eq. (8) reduces to This represents an infinite set of coupled equations. These equations can be solved by truncating the summation over G to a finite number of terms. For the sake of simplicity, we will examine the case when G is limited to G = 0 and − 2π P . This corresponds to periodic coupling realized as the sum of a uniform coupling parameter and one single sinusoidal function. This approximation allows for mathematical tractability while still capturing the effect of band folding due to the periodicity.
If g = 0, Eq. (9) simplifies to The dispersion relation is obtained by setting the determinant of the matrix in Eq. (11) equal to zero, (12) Equation (12) is rewritten as For a given ω, and n, two wave numbers correspond to the G = 0 dispersion curves. The other two wave numbers are associated with the dispersion curves translated by the reciprocal number G = − 2π P (see Fig. 2), i.e., dispersion curves translated into the second Brillouin zone to the left of the first Brillouin zone. Equivalently, these latter P and their corresponding discrete states. The first Brillouin zone is delimited by vertical lines at − π P and π P . The driving frequencies ωD and ω ′ D are represented as a horizontal dashed line. The wave numbers k2 and k3 and corresponding characteristic frequencies ω2, ω3, ω ′ 2, and ω ′ 3 are labeled. The gaps that may have formed at k = − π P are not represented for the sake of simplifying the illustration. The length of the finite waveguides is L = 5P.
wave numbers correspond to dispersion branches, which fold within the first Brillouin zone. If more G vectors were kept in the calculation, the band structure would include dispersion curves translated not only in the second Brillouin zone but also in third, fourth, and higher order Brillouin zones. The extension of this translated dispersion relation would be equivalent to multiple folding of the bands in the first Brillouin zone. However, a single band folding is sufficient here to illustrate the condition for observing a nonseparable superposition of states spanning the complete Hilbert space of product states.
It is worth analyzing Eq. (14) in a few simple cases. In the absence of periodicity in the coupling between the waveguides, α 2π / P = 0 and η 2 n = 0. Equation (12) reduces to the two dispersion relations ω 2 = β 2 k 2 +λnα 0 and ω 2 = β 2 (k + 2π P ) 2 +λnα 0 . These are the dispersion relations for waveguides coupled uniformly along their length in a supercell of length P. Note that in that case at k = 0 and k = − 2π P , the lowest positive and real frequencies are given by ω 2 = λnα 0 .
For the system with periodic coupling, if we insert the condition ω 2 = λnα 0 (i.e., γn = 0) in Eq. (14), we find four nonzero solutions for the wave number kn(ω) = − π P ± √ ( π P ) 2 ± ηn. This indicates that in the vicinity of the wave number origin, the bands of the periodically coupled systems have been lowered compared to those of the uniformly coupled system. Finally, at the edge of the Brillouin zone, k = − π P , Eq. (14) gives the following condition: (( π P ) 2 + γn) At the edge of the Brillouin zone, the degeneracy of the folded bands is lifted by opening a gap. For instance, the positive frequency at k = − π P is two valued.

III. DRIVEN ARRAY OF ELASTIC WAVEGUIDES COUPLED PERIODICALLY ALONG THEIR LENGTH
Let us now drive the system externally at some position x = 0. Equation (1) is replaced by where ωD is the frequency of the driver. The 3 × 1 vector, F, is expressed in the En basis: F = ∑ n FnEn. The Fn's are therefore defined as the dot product F ⋅ En between the two 3 × 1 vectors. Again, we seek solutions in the form of Bloch waves with the driver's frequency, The index D of uD,n indicates that this is the amplitude of the driven system. Following the algebraic steps of Sec. II, we obtain the set of equations In obtaining the right-hand side of Eq. (17), we have used Similarly to Sec. II, Eq. (17) represents an infinite set of coupled equations. Again, limiting the reciprocal lattice numbers to G = 0 and − 2π P , and g = 0 and 2π P , we obtain the set of two linear equations The driven resonant amplitudes are therefore given by where If we employ Eq. (12), we can express the denominator, Δn, as where ωn and ωn,P are the characteristic frequencies corresponding to the nth eigen mode of the coupling matrix with the wave numbers k and k + 2π P [see Eq. (14)]. If the system exhibits some damping, taking the form η ∂ ∂t , where η is a damping coefficient, in Eq. (15), one needs to replace ω 2 D by ω 2 D − iηωD. Importantly, in that case, the amplitudes ( uD,n(k) uD,n(k + 2π P ) ) are complex quantities. Complex amplitudes are important to establish the analogy between nonseparable elastic states and locally "entangled" quantum states described in Sec. IV.

IV. NONSEPARABLE ELASTIC STATES SUPPORTED BY THE PERIODICALLY COUPLED WAVEGUIDES
Here, we demonstrate that the displacement field of the driven system can be written as nonseparable superposition of states. Solving Eq. (17) by truncating the sum over the reciprocal lattice number beyond what was done in Sec. III would result in complex driven amplitudes uD,n(k + g). These complex amplitudes are proportional to Fn, so we rewrite them as uD,n(k + g) = FnũD,n(k + g).
In the case of the three waveguides arranged in a planar array, the coupling matrix N 3×3 admits three eigen vectors and three eigen values. One of the eigen vectors for which λ 1 = 0, namely, , does not involve transfer of energy between the waveguides via the coupling. This trivial case is equivalent to three independent waveguides. The other two eigen modes of the coupling matrix with eigen values λ 2 = 1 and λ 3 = 3 are We now suppose that we are driving the system with F 1 = 0 (to eliminate the E 1 mode) and F 2 ≠ F 3 ≠ 0. The displacement field of Eq. (16) becomes We note that in Eq. (23), we use a discrete summation over the wave number k. This condition can be realized with finite length waveguides obeying periodic boundary conditions. This system supports a discrete set of standing waves. Furthermore, let us suppose that ωD is close to two characteristic frequencies ω 2 and ω 3 corresponding to the second and third eigen modes of the coupling matrix with the wave numbers k 2 and k 3 inside the first Brillouin zone and k 2 + g and k 3 + g inside the other Brillouin zones, respectively (see Fig. 2). The amplitudesũ D,2 (k 2 + g) andũ D,3 (k 3 + g) are expected to dominate, and the sum over k in Eq. (23) can be approximated as The choice of the driving frequency ωD determines how many waves with different k will be composing the superposition of states given by Eq. (23).
If we change the driving frequency from ω d to ω ′ d such that now the largest amplitudeũ D ′ ,2 occurs at the characteristic frequencies ω ′ 2 with the wave number k 3 or k 3 + g and the largest amplitudeũ D ′ ,3 occurs at ω ′ 3 and k 2 or k 2 + g (Fig. 2), then the displacement field takes the new form Equation (26) can be expressed in terms of the product basis, We can simplify these expressions further by considering that the measurements of the displacement fields given by Eqs. (25) and (26) are performed at a distance x = QP, where Q is an integer, from the driver. In that case, ϕ 1,P = E 2 e ik 2 QP , ϕ 2,P = E 2 e ik 3 QP , ϕ 3,P = E 3 e ik 2 QP , and ϕ 4,P = E 3 e ik 3 QP . The product basis becomes independent of the reciprocal lattice number, g, since e igQP = 1. The displacement fields are then reduced to The amplitudes of UD(QP, t) and U D ′ (QP, t) take the form of nonseparable linear combinations of product states. Equations (28a) and (28b) are isomorphic to Bell states. By controlling the driving frequency, one is able to navigate the Hilbert space H E,k by transforming a state of the form Without the periodicity of the coupling, that is, without the folding of the bands within the first Brillouin zone (or, equivalently, the translation of bands within the second Brillouin zone), the state B|E 2 ⟩|k 3 ⟩ + C|E 3 ⟩|k 2 ⟩ would not be accessible to observation. Navigation of the Hilbert space enables transformation of states analogous to quantum gates. For instance, we can also define the basis {ϕ 1,P , ϕ 2,P , ϕ 3,P , ϕ 4,P } in terms of the 4 × 1 vectors In that representation, the two displacement fields become Changing the driving frequency from ωD to ω ′ D is therefore equivalent to applying the following transformation matrix to the amplitude vectors: The most significant aspect of this transformation matrix is that it contains non-zero off-diagonal terms of the two blocks. The physical effects of this transformation are directly measurable in terms of amplitudes and relative phases of the waveguide displacement at the location x = QP.
For instance, one measures the three components (U 1,D (QP, t), U 2,D (QP, t), U 3,D (QP, t)) and expresses them as (G 1 , G 2 e iϕ 12 , G 3 e iϕ 13 ), where ϕ 12 and ϕ 13 are the phases of the displacement of waveguides 2 and 3 relative to waveguide 1. If the amplitude of the first waveguide is normalized to 1, the four quantities G 2 , G 3 , ϕ 12 , and ϕ 13 completely characterize the nonseparable state given by Eq. (29a). Similarly, we can characterize the nonseparable state U D ′ (QP, t) with two other relative amplitudes and two relative phases G ′ 2 , G ′ 3 , ϕ ′ 12 , and ϕ ′ 13 . We also note that the complex quantities A 2,2 and A 3,3 can be tuned via the control of the driver amplitudes F 2 and F 3 and their relative phase. The same is true for the quantities A 2,3 and A 3,2 . Indeed, rewriting F 2 = |F 2 | and F 3 = |F 3 |e iθ , where θ is the phase difference between F 2 and F 3 , we can add θ to the phase difference between ∑ gũ D,2 (k 2 + g) and ∑ gũ D,3 (k 3 + g) when driving at ωD. A similar phase difference can be added between ∑ gũD ′ ,2 (k 3 + g) and ∑ gũD ′ ,3 (k 2 + g) when driving at ω ′ D . Control of the frequency of the drivers and their relative phase provides significant flexibility to manipulate the coefficients of the transformation matrix given by Eq. (30). The significance of the transformation matrix in Eq. (30) is its nonzero off-diagonal terms. With appropriate tuning of the complex amplitudes via the external drivers or of the geometry of an array of waveguides or physical characteristics of the constitutive materials, one may be able to achieve elastic generalized analogs of nontrivial quantum gates such as the Pauli X gate.

V. NAVIGATING THE COMPLETE HILBERT SPACE, H E,k
To navigate the space, H E,k , we need to excite four amplitudes simultaneously with a single driving frequency, to obtain a displacement field of the form UD(QP, t) = (B 2,2 ϕ 1,P + B 2,3 ϕ 2,P + B 3,2 ϕ 3,P + B 3,3 ϕ 4,P )e iω d t . (31) One can excite such a state in the vicinity of a crossing point between the bands corresponding to the coupling modes, E 2 and E 3 . The crossing occurs between the E 2 band with g = 2π P and the E 3 band with g = 0. Using Eq. (14), such a crossing occurs at k 2 (ω) = k 3 (ω). The continuous bands in Fig. 2 exhibit such a crossing point. In the case of finite length waveguides with periodic boundary conditions, and driving the system with a frequency near the crossing point frequency, Eq. (23) can be simplified into the following form: To obtain Eq. (23), we are assuming that only four classes of points in the vicinity of the crossing point frequency contribute significantly to the amplitude of the displacement field. The resonant frequency associated with these points is the closest to the driving frequency. Again, performing measurements at x = QP, where Q is an integer, from the driver, using e igQP = 1, and defining the basis ϕ 1,P = E 2 e ik 2 QP , ϕ 2,P = E 2 e i(− π P )QP , ϕ 3,P = E 3 e i(− π P )QP , and ϕ 4,P = E 3 e ik 2 QP , Eq. (32) becomes This equation can be rewritten in the desirable form of Eq. (31) by defining B 2,2 = [∑ g F 2ũD,2 (k 2 + g)], B 2,3 = [∑ g F 2ũD,2 (− π P + g)], B 3,2 = [∑ g F 3ũD,3 (− π P + g)], and B 3,3 = [∑ g F 3ũD,3 (k 2 + g)]. In the bracket notation, this state takes the general form A|E 2 ⟩|k 2 ⟩ + B|E 2 ⟩|k 3 ⟩ + C|E 3 ⟩|k 2 ⟩ + D|E 3 ⟩|k 3 ⟩. Such a state would not be accessible without the folding of the bands resulting from the periodicity of the coupling.
The product Hilbert space can be explored by changing the frequency to values above or below the frequency of the crossing point. This effectively changes the complex coefficients A, B, C, D. The magnitude and relative phase of F 2 and F 3 offer additional degrees of freedom to manipulate the complex amplitudes of the superposition of states.
By tuning the magnitude of the coupling stiffness as well as the length of the finite waveguides (with periodic boundary conditions), one can create more general states. For instance, in Fig. 3, we are schematically representing g = 0 and 2π P for two E 2 and E 3 bands in the case of finite length waveguides L = 9P. In this particular case, we have four classes of modes a, b, c, and d flanking a crossing point Here, the length of the finite waveguides is L = 9P. The first Brillouin zone is delimited by vertical lines at − π P and π P . The driving frequency ωD is represented as a horizontal dashed line. The wave numbers k2 and k3 correspond to the E3 (a) and E2 (c) modes and E2 (b) and E3 (d) modes, respectively. The gaps that may have formed at k = − π P are not represented for the sake of simplifying the illustration.
within the first Brillouin zone. Driving the system at the frequency ωD shown in the figure, the measured displacement field at QP will take the form UD(QP, t) = {[∑ g F 2ũD,2 (k 2 + g)]ϕ 1,P + [∑ g F 2ũD,2 (k 3 + g)]ϕ 2,P + [∑ g F 3ũD,3 (k 2 + g)]ϕ 3,P + [∑ g F 3ũD,3 (k 3 + g)]ϕ 4,P }e iω D t , with the product basis ϕ 1,P = E 2 e ik 2 QP , ϕ 2,P = E 2 e ik 3 QP , ϕ 3,P = E 3 e ik 2 QP , and ϕ 4,P = E 3 e ik 3 QP . By driving the system at a frequency, ω t . This is also a nonseparable state. The transition from one of these nonseparable states to the other one is solely allowed by the band folding due to the periodicity of the coupling.
We note that Eq. (34) becomes separable if we can find a driving frequency such that

ARTICLE scitation.org/journal/adv
Then, Eq. (34) can be written in the product form The relative amplitudes and phase of the three components (U 1,D (QP, t), U 2,D (QP, t), U 3,D (QP, t)) expressed as (G 1 , G 2 e iϕ 12 , G 3 e iϕ 13 ) are simply dependent on the factor F 2 E 2 + F 3 E 3 . The factor (D 2 e ik 2 QP + D 3 e ik 3 QP ) only adds a global phase to each waveguide. The phases of the displacement of waveguides 2 and 3 relative to waveguide 1 parallel the relative phase of the external drivers. This is a signature of separability of the elastic states. Nonseparable states will exhibit relative phases between waveguides, which differ from the relative phase of the drivers. It is therefore possible with this simple measurement to differentiate nonseparable from separable superpositions of elastic states.
Finally, we also recall that in the model system composed of three coupled waveguides, there are only two non-trivial discrete eigen vectors describing the spatial states, namely, E 2 and E 3 . This represents, effectively, a two-level system. Furthermore, due to the finite length of the waveguides, the plane wave states are discrete. We have chosen driving frequencies such that only two plane wave states may contribute the most to the displacement field. The system behaves like a two-level system in plane wave states. The Hilbert space of this restricted system is therefore a 2 × 2 dimensional space. The periodicity of the coupling between waveguides, leading to band folding, enables the complete exploration of this Hilbert space. Increasing the number of coupled waveguides will increase the number of non-trivial spatial eigen modes. For instance, in the case of four coupled waveguides, the spatial states represent a three-level system. The spatial states of N coupled waveguides would span an N − 1 dimensional space. By increasing the length of the waveguides, one increases the number of discrete plane wave states within a given range of frequencies. It is therefore possible to extend the space of dominant plane wave states beyond the two investigated here. For instance, one can consider a driving frequency that would excite three plane wave states with significant amplitudes. In that case, the plane wave states would span a three-dimensional space corresponding to a three-level system. Generalizing to n plane wave states corresponding to an n-level system, an elastic system of N waveguides could support superpositions of states in an (N − 1) × n Hilbert space of product states.

VI. COMPUTER MODEL AND SIMULATION
To validate the theory in Sec. V, we develop a numerical model based on coupled mass-spring waveguides that can be used to shed light on the complex behavior of nonseparable states and hence on the complex Hilbert space, H E,k . For this, we consider a set of 1D discrete elastic waveguides composed of chains of identical masses connected to each other via harmonic springs. We further assume that the masses of each chain are constrained to move only in the horizontal direction. This limits the model to displacements with a single polarization as was the case in Sec. V. Dissipative effects in the medium are modeled by linear viscous damping elements.
For a set of three coupled mass-spring chains with a total of Nm identical masses in each chain, the discrete elastic equations of motion are given by where un, vn, and wn are the displacements of the nth mass of waveguides 1, 2, and 3, respectively. These quantities form the spatially discrete components of the 3 × 1 vector U used in Secs. II-V for a continuous system. m is the mass, and the viscous damping coefficient η models the dissipation. The term knn describes the coupling constant of the nearest-neighbor interaction along a waveguide, and kc describes the stiffness of the springs that couple the masses between the waveguides organized in a planar array. The spring constant kc has a sinusoidal modulation so that where k c,0 is the unmodulated spring constant and Δkck c,0 is the modulation amplitude. Here, we take kc to be the same for all coupled chains. The superposition of a constant coupling stiffness and a simple sinusoidal variation parallels our choice of reciprocal numbers limited to the two values G = 0 and − 2π P in Secs. II-V.
We use Born-von Karman periodic boundary conditions for which e ikN c P = 1; Nc is the total number of unit cells along the waveguides. As mentioned above, in the first Brillouin zone, kP is limited to the interval −π to π with a spacing of 2π/Nc. The set of Eqs. (36a)-(36c) is the discrete equivalent of Eq. (1). A vectorized fourth order Runge-Kutta time integration scheme is used to numerically calculate the dynamics of the system. To computationally generate the band structure, the spectral analysis of amplitudes and phases method 24 is employed, which we have recently developed and entails the use of molecular dynamics simulation. The band structure in the absence of damping is reported in Fig. 4. The specific values of the model physical constants are presented in the figure caption. From Fig. 4, it is clear that the E 1 mode does not involve transfer of energy between the waveguides via the coupling and hence does not fold in the Brillouin zone. On the other hand, near the edges of the first Brillouin zone, the bands corresponding to the E 2 and E 3 spatial eigen modes fold and hence intersect. Because of the finite length of the waveguides, four discrete modes (two closed and two open circle modes) near these intersections form two pairs of modes with the same wave numbers. This configuration is isomorphic to that schematically presented in Fig. 3.
Each mode represented by a single closed or open circle in the band structure of Fig. 4 corresponds to a separable state in the Hilbert space H E,k , i.e., states expressible as products of eigen modes spatially varying across the array of waveguides and Bloch wave degrees of freedom along the waveguides. From Fig. 4, we identify different isofrequency values that correspond to nearly overlapping E 2 and E 3 modes. For instance, the inset of Fig. 4 shows that driving the model at a frequency of 0.622 or 0.6 will lead to non-separable superpositions of states with characteristics similar to those of the theoretical superpositions U . To excite such separable or nonseparable superpositions of states numerically, the coupled elastic system is driven with the external force ⃗ Fe iω d t ; ⃗ F = (1 − μ)E 2 + μE 3 . Moreover, by tuning the parameter μ, the contribution of the two eigen modes (E 2 , E 3 ) is numerically manipulated. Finally, to drive the system, the prescribed base periodic excitation F u (t) is applied to the first mass of each chain, where (F u 0 , F v 0 , F w 0 ) and (ϕ u 0 , ϕ v 0 , ϕ w 0 ) are the amplitudes and phases of the applied force to the chain (1, 2, 3). Figure 5 shows the variations of the maximum amplitudes of the last mass of each waveguide Ci; C 1 = max(uN m (t)), C 2 = max(vN m (t)), C 3 = max(wN m (t)) and the phase difference (ϕij) between the transmissions for each pair of waveguides as a function of μ, where ϕ 12 = 180 π cos −1 ( |u Nm (t)||w Nm (t)| ). For a pure eigen mode, the phase difference between the transmissions for each pair of chains is ϕ E 2 12 = ϕ E 2 23 = π/2 and ϕ E 2 13 = π for E 2 and ϕ E 3 12 = ϕ E 3 23 = π and ϕ E 3 13 = 0 for E 3 . From Fig. 5, we indeed see that manipulation of the parameter μ can be used to tune the eigen mode superposition, since for μ ≠ 0 and μ ≠ 1, the output mode is a linear combination of spatial eigen modes E 2 , E 3 with corresponding k-labeled Bloch waves. In Fig. 5, we drive the system at different frequencies, between ω (c,d) D = 0.6 that is in the vicinity of the two modes c and d and ω (a,b) D = 0.622 that is in the vicinity of the two modes a and b. In particular, for ωD = 0.618, we note that at μ = 0, the spatial character across the planar array of waveguides that is measured is a pure E 2 state. For 0 < μ < 0.2, the measured spatial state is a superposition of modes E 2 , E 3 . At μ = 0.2, the output mode is the well-defined superposition (1 −1 0). For μ > 0.2, the measured spatial output suggests a state very close to E 3 until μ = 1, where the measured spatial character is purely E 3 . Figure 6 shows the variations of the phase difference as a function of the driving frequency ω d and for different values of μ. We observe that some of the phase differences between pairs of waveguides change sign between the frequencies 0.6165 and 0.618. To shed light on these numerical results, we write an approximation to the amplitude of the displacement field, UD, by considering only the contribution of the modes a, b, c, d within the first Brillouin zone, where the complex amplitudes resulting from the four modes are written as simple resonances with damping. k 2 (= ka,c) is the wave number of modes a and c. The wave number k 3 (= k b,d ) corresponds to modes b and d.
We now consider a special case by choosing ωD between ωc and ω d such that −ω 2 D + ω 2 c = −δ and −ω 2 D + ω 2 d = δ. Note that δ is small. In that case, only two terms dominate. The approximate displacement field has now the components When choosing another special case by driving the system at a frequency ω ′ D that lies between ωa and ω b such that −ω ′2 D + ω 2 a = −δ and −ω ′2 D + ω 2 b = δ, the components of the approximate field amplitude become The ratio of the amplitude measured at x = QP along the first waveguide relative to that measured along the second waveguide is given by

ARTICLE scitation.org/journal/adv
where the phase difference ϕ 12 = ϕ 1 − ϕ 2 . The ratio of amplitudes between waveguides 1 and 2 for the driving frequency ω ′ D is obtained as Since the frequencies of the modes a, b, c, d are within a few percent of each other, the driving frequencies ωD and ω ′ D necessary to transform the displacement field from UD to U D ′ are also within a few percent of each other. Therefore, Eq. (42) approximates the complex conjugate of Eq. (41). Subsequently, one has ϕ ′ 12 ∼ −ϕ 12 . This is what is observed in Fig. 6. This observation can also be made for the ratios 32 and , that is, ϕ ′ 32 ∼ −ϕ 32 that automatically applies to ϕ 23 . However, with these simple cases, the approximate complex conjugation may not be true for ϕ 13 in light of Eqs. (40a) and (40c). We recall that for a given set of driving parameters F 2 and F 3 , i.e., for a given value of μ, raising of the driving frequency from ωD to ω ′ D executes the matrix operation of Eq. (30) on the displacement field amplitudes expressed in the H E,k basis. Complex conjugation is the signature of this operation in the space of the complex amplitudes of individual waveguides (C 1 e iϕ 1 , C 2 e iϕ 2 , C 3 e iϕ 3 ) or of the complex amplitudes relative to the first waveguide (G 1 , G 2 e iϕ 12 , G 3 e iϕ 13 ).
To quantify the level of control the model allows, on the extent of the Hilbert space that can be explored, we calculate the coefficients in the transformation matrix of Eq. (30). For the sake of completeness, we rewrite Eq. (28) by including mode E 1 , although we anticipate its contribution to be negligible in light of the way the system is stimulated. At the driving frequencies ω D = 0.6 and ω D' = 0.622, we get UD(QP, t) = (A 1,1 E 1 e ik 1 x + A 2,2 E 2 e ik 2 x + A 3,3 E 3 e ik 3 x )e iω D t , (43) In writing Eqs. (43) and (44), we have used the fact that at the end of each waveguide, i.e., at x = QP = Nc where the measurements are performed, e igx = 1. Furthermore, in the mass-spring system with periodic boundary conditions, the wave numbers are integer multiples of 2π/Nc, and hence, e ik 1 x = e ik 2 x = e ik 3 x = 1. Simplifying the equations above, we obtain The measured amplitudes as a function of the driving parameter μ on the left of Eqs. (45) and (46) are used to calculate the complex coefficients A 1,1 , A 2,2 , A 3,3 , A 2,3 , and A 3,2 [as shown in Fig. 7(a)]. We find, as anticipated, that A 1,1 [shown by an asterisk in Fig. 7(a)] is zero. With the frequency ωD and for μ = 0, A 3,3 = 0 and the spatial state is a pure mode E 2 . For μ = 1, the spatial state is a pure mode E 3 with A 2,2 = 0. Similarly, with the frequency ω D ′ and for μ = 0, A 3,2 = 0 and the spatial state is also a pure mode E 2 . For μ = 1, the spatial state is a pure mode E 3 with A 2,3 = 0. Following Ref.

ARTICLE scitation.org/journal/adv
we can express the entropy of "entanglement" for the two driving frequencies We note that when |A 2,2 | = |A 3,3 |, S = ln 2. The superposition of states is maximally "entangled." Maximum entanglement is also achieved for the driving frequency ω D ′ when |A 2,3 | = |A 3,2 |, and S ′ = ln 2. We can calculate the ratios of complex amplitude coefficients, A 2,3 /A 2,2 and A 3,2 /A 3,3 , constituting the off-diagonal terms in the non-trivial transformation matrix (quantum-like gate) of Eq. (30). The modulus and the corresponding phase of these ratios are reported in Fig. 7(b). The moduli remain constant for all values of μ. However, from Fig. 7(b), it is clear that we can navigate the phases between 0 and 3π/10. Therefore, we have liberty to manipulate the drivers such that one can tune the block transformation matrix between the two nonseparable states of Eqs. (43) and (44). For instance, when μ → 0, the first block of the transformation matrix becomes proportional to the Pauli X matrix: ( 0 A 2,3 /A 2,2 A 2,3 /A 2,2 0 ) ∝ ( 0 1 1 0 ), while the second block is proportional to the nontrivial matrix: ( 0 A 3,2 /A 3,3 A 3,2 /A 3,3 0 ) ∝ ( 0 cos 3π/10 + i sin 3π/10 cos 3π/10 + i sin 3π/10 0 ).
Finally, we remark that the case discussed above is only one possible example of manipulation of states within the Hilbert space of product states and exploration of that space can be enlarged by choosing other input parameters for the drivers such as changing the relative phase of the drivers.

VII. CONCLUSIONS
A planar array of three one-dimensional elastic waveguides mutually coupled periodically along their length can be driven externally to support nonseparable superpositions of states. These states are the product of Bloch waves characterizing the elastic displacement along the waveguides and spatial modes representing the displacement across the array of waveguides. For a system composed of finite length waveguides, the frequency of external drivers can be employed to excite specific discrete modes. The periodicity induced band folding realizes superpositions of states that span the complete available Hilbert space of product states. By tuning the frequency and the relative amplitudes of the external drivers stimulating individual waveguides, we can navigate a broad segment of the Hilbert space of product states. Transformation from one type of nonseparable superposition to another one in the product Hilbert space is analogous to a generalized Pauli gate. This transformation can also be interpreted as the complex conjugation operator in the space of the complex amplitudes of individual waveguides. Using similar elastic waveguides, although not periodic, we demonstrated Hadamard gate operations on the coherent superposition analogous to the single qubit gate operation in quantum computing. 17 Finally, nonseparable or "classically entangled" states of elastic waves offer the advantage of stability over entangled states of true quantum systems. Nonseparable superpositions of elastic waves are robust against decoherence and will not require operating at cryogenic temperatures to maintain the delicate balance of the superpositions. Nonseparable superpositions of elastic waves do not suffer from the phenomenon of wave function collapse upon measurement. A coherent superposition of quantum states collapses into a pure state upon measurement. Multiple statistical measurements are, therefore, necessary to obtain information on the original superposition. Moreover, the present demonstration of (a) "classically entangled" superposition of elastic states that span their complete Hilbert space and (b) a transformation from one type of nonseparable superposition to another one analogous to a nontrivial quantum gate is a significant step toward realizing the potential of elastic waves in capturing quantum-like phenomena previously associated with optics. Recently, the use of classical light with entangled degrees of freedom has found applications in quantum information 33 and metrology. 34,35 The present study suggests that the same sort of applications can be realized with acoustic systems. The periodicity of the coupling between elastic waveguides was the determining factor for widening the Hilbert space accessible to elastic superposition of states. This coupling can be easily manipulated by choices of materials and fabrication. Moreover, due to the flexibility of the elastic system, the coupling can easily be tailored to be linear, periodic, and nonlinear (with different types or degrees of nonlinearity). 36,37 AUTHORS' CONTRIBUTIONS All authors contributed equally to this work.

ACKNOWLEDGMENTS
We acknowledge financial support from the W.M. Keck Foundation.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.