The topological characteristics of waves in elastic structures are determined by the geometric phase of waves and, more specifically, by the Berry phase, as a characterization of the global vibrational behavior of the system. A computational procedure for the numerical determination of the geometrical phase characteristics of a general elastic structure is introduced: the spectral analysis of amplitudes and phases method. Molecular dynamics simulation is employed to computationally generate the band structure, traveling modes' amplitudes and phases, and subsequently the Berry phase associated with each band of periodic superlattices. In an innovative procedure, the phase information is used to selectively excite a particular mode in the band structure. It is shown analytically and numerically, in the case of one-dimensional elastic superlattices composed of various numbers of masses and spring stiffness, how the Berry phase varies as a function of the spatial arrangement of the springs. A symmetry condition on the arrangement of springs is established, which leads to bands with Berry phase taking the values of 0 or π. Finally, it is shown how the Berry phase may vary upon application of unitary operations that mathematically describe transformations of the structural arrangement of masses and springs within the unit cells.

## I. INTRODUCTION

The paradigm of the plane wave has served as the foundation of our understanding of sound and elastic waves. The four wave characteristics of frequency (ω), wave vector (*k*), amplitude (*A*), and phase (*ϕ*) undergird this paradigm. For the last 20 years or so, the manipulation of spectral and refractive properties of elastic waves using host materials that exploit frequency and wave vectors has led to significant advances in the fields of phononic crystals and acoustic metamaterials (Deymier, 2013). Elastic waves' spectral properties include the formation of stop bands caused by either Bragg-like scattering or resonant processes, and the capacity to achieve narrow band spectral filtering by introducing structural defects in the wave supporting material. Negative refraction, zero-angle refraction, and other unusual refractive properties arise from the manipulation of the full range of the dispersion relations for elastic waves over both frequency and wave vector through changes in the structure.

The amplitude and phase characteristics of elastic waves have garnered more attention recently. When sound waves propagate in media that exhibit symmetry breaking conditions, they may exhibit amplitudes with a geometric phase leading to non-conventional topologies (Deymier and Runge, 2017) and, under certain conditions, to non-reciprocal wave propagation. Examples for the breaking of symmetries, including time-reversal symmetry, chiral symmetry, and particle-hole symmetry, have been discussed elsewhere (Deymier and Runge, 2017). An alternative realization of intrinsic parity symmetry breaking, which is comprised of a one-dimensional (1D) harmonic crystal with masses attached to a rigid substrate through harmonic springs, has been shown to possess a spin-like topology that can be described by a Dirac-like equation (Deymier *et al.*, 2014, 2015). Complementing intrinsic symmetry breaking are extrinsic topological phononic structures, which have been created using a periodic spatial modulation of the stiffness of a 1D elastic medium such that its directed temporal evolution breaks both time-reversal and parity symmetries (Croënne *et al.*, 2017; Deymier *et al.*, 2017; Deymier and Runge, 2016; Nassar *et al.*, 2017b; Nassar *et al.*, 2017a; Swinteck *et al.*, 2015; Trainiti and Ruzzene, 2016). As the phase of these amplitudes plays a central role in the exploration of the topological characteristic of elastic structures, we present here a computational tool for the calculation of these phases for a general elastic structure. While in the current manuscript we test the computational tool on simple mass and spring models, the method can be applied to continuous systems, including rods, beams, and plates. The method is employed for a number of examples that may intrinsically break parity symmetry to assess the geometric phase of these elastic systems. Though trivial and non-trivial topologies are typically associated with integral multiples of π, the examples considered here demonstrate the full range of geometric phases attainable showing the flexibility of the method.

In Sec. II, we illustrate the computational tool for the spectral analysis of amplitudes and phases (SAAP) of elastic waves in the 1D case. In Sec. III, we address the calculation of the Berry connection and the Berry phase as a characterization of the global vibrational behavior of an elastic system and its topological character. In Sec. IV, we present results of SAAP calculations for a variety of 1D mass-spring systems with varying spring constant arrangements, discuss them in comparison with analytical results, and also analyze the behavior of the Berry connection and Berry phase with the change of Ansatz and origin of the system of reference. Finally, we draw conclusions in Sec. V.

## II. SAAP OF ELASTIC WAVES

Normal modes are synchronous time-periodic motions during which all coordinates of the system vibrate in a synchronous fashion, reaching their maximum and minimum values at the same time. Consider an elastic system composed of masses connected by harmonic springs, when we apply a set of initial conditions in the form of random displacements and/or velocities of the masses, after an initial transient response, a steady state is reached. In the case of a 1D periodic system composed of an arrangement of $Nc$ identical unit cells, each containing some number of masses $Nm$, we define $un,Ni(t)$ as the displacement from equilibrium as a function of time *t* of the *n*th mass in the *N _{i}*th unit cell of the chain. In general, multiple frequencies can be supported for a single wave number

*k*, we denote the frequencies as $\omega j(k)$, which refers to the

*j*th lowest frequency the system supports for wave number

*k*. The

*j*labels the bands in the elastic band structure of the system. Therefore, the complex amplitude $An,j(k)$ of the oscillation of the

*n*th mass on the

*j*th band is defined as the following projection of $un,Ni(t)$ onto a plane wave with wave vector

*k*and frequency $\omega j(k)$:

where $\tau 0$ is the total time over which the elastic waves are sampled, and *L* is the unit cell length. The phase of the complex amplitude is given by

where

Note that the phase of $z$ is defined on the interval $\u2212\pi $ to $\pi $.

The outlined procedure is general for any displacement pattern. For the sake of illustration and testing, this paper focuses on the application of this procedure to discrete 1D harmonic elastic systems. To calculate the traveling modes' amplitudes and phases, we introduce a procedure that is illustrated in the workflow in Fig. 1. The dispersion relation is first calculated by finding the unknown frequencies $[\omega j(k)]$ for a given wave number $k$. Employing molecular dynamics (MD) simulation with the chosen initial conditions $un,Ni(0)=cos(kNiL)$ and $u\u0307n,Ni(0)=0$ and Born-von Karman periodic boundary conditions for which $eikLNc=1$, the elastic equations of motion are solved with the Runge-Kutta (4,5) formula, which engages all the springs. In the first Brillouin zone, $kL$ is limited to the interval $\u2212\pi $ to $\pi $ with a spacing of $2\pi /Nc$. After each MD run for a specific wave number, the frequency spectrum is found by the temporal integral of the displacement $un,Ni(t)$,

Equation (3) is like a Fourier transform whose peak positions define the frequencies for the given value of *k*, labeled $\omega j(k)$. The elastic band structure is fully generated by applying this procedure to all available values of *k* for the finite system within the first Brillouin zone.

Once the band structure $\omega j(k)$ is known, we calculate the elastic wave amplitudes and phases using another set of MD simulations with new initial conditions $un,Ni(0)=cos(kNiL)$ and $u\u0307n,Ni\u2009(0)=\u2212\omega j(k)sin(kNiL)$. This new set of initial conditions now sets the values of the velocity to those prescribed by the computed band structure. We emphasize here that one has to use specific traveling elastic wave initial conditions instead of random initial conditions to obtain the phases since random initial conditions add or subtract a constant value to the phase. This is a distinction between the SAAP method and other eigenvalue/vector approaches. After each set of MD simulations, we project the calculated displacement $un,Ni(t)$ onto plane waves to calculate the complex amplitude of the *n*th mass, $An,j(k)$, following Eq. (1). The phase of the displacement is calculated using Eq. (2). We have implemented SAAP in MATLAB (The MathWorks, Natick, MA) and tested it for several systems. As we have seen above, in the SAAP method Newton's equations of motion are solved as a function of time via a MD approach, and the solutions are projected onto plane waves, which allows the extraction of the amplitudes, including their phases.

## III. BERRY CONNECTION AND BERRY OR ZAK PHASE

A further characterization of the behavior of the amplitudes is realized by the Berry connection for a particular band of the band structure. When the system contains a finite number of unit cells, $\u2009NC$, the Brillouin zone is discretized and the Berry connection is given by (Resta, 2000)

where $Nm$ is the number of masses per unit cell, $\Delta kL=2\pi /NC$, and $A\u0303n,j(k)$ is the normalized complex amplitude for the *n*th mass in a unit cell on band *j*, and is calculated as

The $A\u0303n,j(k)$ with $n\u2208[1,Nm]$ represent the complex *n*th components of a unit vector in a $Nm$-dimensional complex space. This normalized amplitude unit vector (n.b., the normalization is for the full vector) can be defined as

This amplitude unit vector evolves along some parametric curve as *k* is varied. The Berry connection characterizes the variation in orientation of the unit vector along some path in the complex space of amplitudes parametrized by *k.* Summing the Berry connection for the possible $k$ values of some specific band over the closed path in $k$-space defined by the first Brillouin zone gives the Berry or Zak phase (Berry, 1984; Zak, 1989). The geometric phase that characterizes the property of bulk bands in periodic systems is known as the Zak phase (Zak, 1989), whereas Berry introduced the general concept of geometric phase earlier (Berry, 1984). Several recent works have been done to correlate the bulk properties of a crystal through the geometrical phases of the bulk bands (Atala *et al.*, 2013; Chen *et al.*, 2018; Pal and Ruzzene, 2017; Wang *et al.*, 2015; Wang *et al.*, 2018; Xiao *et al.*, 2014). The importance of finding the topological structures through geometric phases is that according to the bulk edge correspondence principle, at an interface between topologically different crystals, a localized state forms (Trainiti and Ruzzene, 2016). Topological interface states have prospects for realization of disorder-robust one-way transport of information. The Berry's phase $\varphi B$ of each band is then defined as (Resta, 2000)

In Eq. (7), Im takes the imaginary part of its argument. For infinite systems or continuous *k*, a differential form of the above expression is used (Berry, 1984). The evolution of the amplitude unit vector in the *N _{m}*-dimensional space parametrized by the wave number

*k*generates a manifold. The Berry phase is the net phase accumulated by the amplitude unit vector over the entire manifold for a closed path. That is to say, the Berry phase characterizes the topology of the manifold.

Notice that the calculation of the Berry phase is conditioned by the existence of a closed path in the *k*-space with the amplitudes or projections $A\u0303n,j(k)$ being periodic in reciprocal space, but this periodicity may be affected by the Ansatz used to describe the dynamics of the system. In fact, the spectrum $[\omega j(k)]$ is invariant under unitary transformations, and therefore to calculate it there is a freedom of choosing a particular Ansatz. In what follows we will use a compact Ansatz for computational convenience, but also illustrate how using a general Ansatz will impact the Berry phase. We show in the Appendix that these two Ansätze possess the same Berry connection and Berry phase.

## IV. RESULTS AND DISCUSSION

We proceed to do the spectral analysis of elastic wave amplitudes and phases of 1D mass-spring systems with identical masses but variable spatial modulations of the spring stiffness. The unit cell length parameter *L* is given by the number of masses per unit cell $Nm$ times the inter-mass spacing $a$, that is, $L=Nma$. Again, the system contains $Nc$ unit cells of $Nm$ masses each.

### A. Two mass unit cell

We first consider a system of two identical masses per unit cell with spring constants modulation for which the analytical solutions of the eigen-problem are available and can be compared with numerical results obtained via SAAP. The equations of motion for masses 1 and 2 of the $Ni$ th unit cell interacting with spring constants $\beta 1$ and $\beta 2$ in an infinite 1D system are (see Fig. 2)

In Eq. (8), $m$ is the mass, $\beta 1$ and $\beta 2$ are the force constants of the spring, and $t$ is time. The dots denote differentiation with respect to time.

#### 1. Analytical solution

We seek traveling wave solutions of Eq. (8) using the following compact Ansatz:

which, as anticipated in connection to the Berry phase calculation, have amplitudes $An$ that are periodic in reciprocal space. That is also in contrast to the Ansatz of the general form $un,Ni(t)=An\u2032eikxneikNiLei\omega t$, where $An\u2032eikxn$ is periodic in reciprocal space. We note that both Ansätze are related by a unitary transformation that yields the same spectrum, Berry connection, and Berry phases (as is shown in the Appendix), however, the compact Ansatz of Eq. (9) is computationally more convenient.

with $(\beta 1+\beta 2\u2009e\u2212ikL)=\gamma $ and $(\beta 1+\beta 2\u2009\u2212m\omega 2)=\alpha $. Taking the complex conjugate of Eq. (10) and exchanging rows and columns, we have

Further, from Eq. (10) we obtain the dispersion relation

where the first sign corresponds to the acoustic branch and the second sign corresponds to the optical branch in the band structure. Choosing these solutions for the frequencies, we have for the amplitudes

Taking into account Eq. (12) leads to the following selection for the amplitudes of the acoustic branch:

since for the acoustic branch $A1,1/A2,1=\gamma /\gamma *.$ For the optical branch

which satisfies Eq. (14) taking the minus sign.

In Fig. 3, we plot the real and imaginary components of the normalized complex amplitude, $A\u0303n,j(k)$, where the normalization constant is chosen as in Eq. (5). Figure 3(a) shows the real and imaginary components of the normalized amplitudes for both masses and both the acoustic and optical branches when the ratio of the spring constants is $\beta 1/\beta 2=0.5$. Note that, for our compact Ansatz, the real component of the amplitude is an even function of *k* and the imaginary component is an odd function, the change in sign of the imaginary component leads to a discontinuity in the optical branch amplitudes. Figure 3(b) displays the same information about the real and imaginary components of the amplitude when the spring constant ratio is $\beta 1/\beta 2=2$. Again, we note the discontinuity in the imaginary component of the optical branch amplitude. From Figs. 3(a) and 3(b), we observe that although the real components of the complex amplitudes are equal for masses 1 and 2, the magnitude of the imaginary components are equal but of opposite signs. As a result, from Figs. 4(a) and 4(b) it can be seen that each mass attains a phase value both in the acoustic and optical branches to sustain the traveling wave, and the magnitudes of the phases of masses 1 and 2 are the same but of opposite signs. We note that this relationship likely arises due to the symmetry of the system as mass 1 is connected with springs $\beta 1$ and $\beta 2$ to its right and left, respectively, and mass 2 is connected with springs $\beta 2$ and $\beta 1$ to its right and left, respectively. Further, the phase value is an odd function of the wave number, i.e., $\varphi n,j(k)=\u2212\varphi n,j(\u2212k)$.

In Fig. 5, we plot a representation of the manifold spanned by the real and imaginary parts of the amplitude unit vector, $A\u0303\u2192$, within the Brillouin zone. At each *k* point, mapped onto the angle to form a ring, we have an arrow pointing in the $A\u0303\u2192$ direction (thin red arrow: $k$ varies from $\u2212\pi /L$ to $0$; thick blue arrow: $k$ varies from $0$ to $\pi /L$). We take $A\u03032,j(k)$ along the normal to the ring plane and $A\u03031,j(k)$ along the radius with the positive direction pointing away from the center of the ring. As $A1=A2*$, the $A\u0303\u2192=(A\u03031,j(k),A\u03032,j(k))$ vectors form a 45° angle with the normal to the ring plane, and hence there exists only twists but not a change in the direction of $A\u0303\u2192$. Moreover, the amplitude unit vector in Fig. 3(a) ($\beta 1/\beta 2=0.5$) generates a manifold with the imaginary part taking the form of a closed ribbon with a single twist at *k *=* *0 [see Fig. 5(b)]. The total accumulated geometric phase over the Brillouin zone is therefore $\pi $. In contrast, the amplitude unit vector in Fig. 3(b) ($\beta 1/\beta 2=2$) generates a manifold with the imaginary part taking the form of a closed ribbon with two twists at *k *=* *0 and at the boundaries of the Brillouin zone [see Fig. 5(a)]. The total accumulated geometric phase over the Brillouin zone is therefore zero. These two manifolds differ in their topology by one twist.

We now characterize the behavior of the amplitudes for each band as the wave number spans the Brillouin zone. Using Eq. (4), defining the Berry connection, we have

$BCj(k)$ given by Eq. (17) is purely real with either positive or negative values and therefore a Berry phase of either 0 or $\pi $. As an example, in Fig. 6 we plot the Berry connection for the two mass unit cell with $\beta 1/\beta 2=2$ and $\beta 1/\beta 2=0.5,$ for both the acoustic [$BC1(k)$] and optical [$BC2(k)$] branches. As we see, the Berry connections are all real. Using Eq. (7), we find that the Berry phase is $0$ for each branch when $\beta 1/\beta 2=2$ and $\pi $ when $\beta 1/\beta 2=0.5,$ as anticipated. Therefore, the topology of the elastic waves (Fig. 5) supported by these configurations are characterized by dissimilar Berry phases.

One important point to realize is that the two mass unit cell has inversion symmetry, and we have found that its Berry phase is always either 0 or $\pi $. It will be shown later in Sec. IV B 1 that in the case of unit cells with more than two masses, inversion symmetry (as imposed by the values of the spring constants) also leads to a multiple of $\pi $ (either 0 or $\pi $) Berry phase. However, multiple of $\pi $ Berry phase values, including zero, are also conditioned by the selection of the origin of the unit cell. That is similar to what Zak found for electron systems with inversion symmetry, except that for electrons there is no dependency of the Berry phase on the choice of origin (Zak, 1989), i.e., for centrosymmetric crystals, when the origin is at a center of inversion, the only allowed values are either 0 or $\pi $ (mod $2\pi $) for electrons.

The arrangements of the spring constants $\beta 1/\beta 2=2$ to $\beta 1/\beta 2=0.5$ can be related to each other by a change in origin of the unit cell from mass 1 to mass 2, i.e., by shifting the origin to the left by an inter-mass spacing, as shown in Fig. 7.

Prior to translating the origin, the dynamical problem is described by Eq. (10), which can be written in the compact form $M\u2194kA\u2192k=0,$ where $A\u2192k=[A1,A2]T$. The dynamical matrix is defined by

After the shift in origin, the dynamical equations take the form $M\u2194k\u2032\u2009A\u2192\u2032k=0,$ where

and $A\u2192\u2032k=[A\u20321,A\u20322]T$ with $An\u2032$ being the amplitudes of the compact Ansatz $un,Ni(t)=An\u2032eikNiL+i\omega t;\u2009\u2009n=1,2$.

The general and compact Ansätze are related by the unitary transformation $A\u2192\u2032k=S\u2194kA\u2192k$ with

and correspondingly the dynamical matrices $M\u2194k$ and $M\u2194\u2032k$ transform as $M\u2194k=S\u2194k\u2020M\u2194\u2032k\u2009S\u2194k.$ The unitary matrix has the following properties: $S\u2194\u2212k=S\u2194k*$ and

The Berry connection of the system after the shift of the origin is given by

Using Eqs. (20) and (21) we find $S\u2194k\u2020\Delta S\u2194k=\Delta k(000\u2212iL)$. Therefore, Eq. (22) yields

where we have used the fact that for the system with a two mass unit cell, $|A1|2=|A2|2=1/2$. After summing over the first Brillouin zone we get that the Berry phase $(\varphi B)$ transforms upon the shift in origin according to

It is clear that if the Berry phase for the original system is $0$ (or $\pi $) then the Berry phase for the system with the shifted origin will be $\pi $ (or $\pi +\pi =2\pi =0,$ since the Berry phase is defined mod $2\pi $). These findings are consistent with the Berry phase values previously obtained [see Eq. (17) and Fig. 6].

#### 2. Numerical study

In order to numerically calculate the complex amplitude $An,j(k)$ of each mass in the unit cell, we first calculate the frequency spectrum $\omega j(k)$ for corresponding wave number $k$ using the initial conditions as $un,Ni(0)=cos(kNiL\u2009)$ and $u\u0307n,Ni(0)=0$ (as mentioned in Sec. II).

Figure 8(a) shows the power spectrum of the amplitudes for two different wave numbers for the two mass system with $\beta 1=2\beta 2.$ For each specific $k$ value we clearly obtain two frequencies for the acoustic branch [$j=1$; smaller $\omega j(k)$] and the optical branch [$j=2$; larger $\omega j(k)$]. Figure 8(b) shows the resultant dispersion relations obtained numerically. As is well known, Fig. 8 depicts the presence of a band gap due to spatially varying stiffness. The width of the gap depends on the relative values of $\beta 1$ and $\beta 2$ (when $\beta 1=\beta 2$ there is no gap). In Fig. 8(b), we also compare the dispersion relations obtained theoretically using Eq. (13) (*, asterisk). The results are in excellent agreement. The dynamical trajectories generated by the MD simulation can also be analyzed within the framework of the spectral energy density (sed) method (Thomas *et al.*, 2010) for analyzing the elastic band structure of superlattices. However, the sed method is computationally very expensive (because one must scan for all possible frequency values) and does not capture the information on the phases of the amplitudes, whereas using the SAAP method we can easily calculate the elastic band structure since the frequency spectrum $\omega j(k)$ for corresponding wave number $k$ is known, while also extracting the phases.

Numerical simulations of Eqs. (1), (2), and (8) for the two mass unit cell produce the normalized complex amplitudes and phases predicted theoretically, validating our SAAP procedure and numerical implementation.

Let us now focus on two different unit cells with $\beta 1/\beta 2=1/2$ and $\beta 1/\beta 2=1/8$ spatial variation of the stiffness. In Fig. 9, we plot the normalized complex amplitudes both for $\beta 1/\beta 2=1/2$ and $\beta 1/\beta 2=1/8$, and observe that both the real and imaginary values at the boundary of the Brillouin zone (i.e., at $kL=\xb1\pi $) remain constant regardless of the spatial variation of the stiffness, though the magnitude of the real and imaginary values are different for different wave numbers. This indicates that regardless of the spatial variation of the elastic stiffness, if $\beta 1/\beta 2<1$, the Berry phase under these conditions is invariant and takes on the value $\pi $. Similarly, if $\beta 1/\beta 2>1$, the Berry phase is always $0$.

Another important aspect of calculating the phase value is that the phase can be used to selectively excite the acoustic or optical branch. As shown before in Fig. 8(a), for a particular wave number both acoustic and optical branches are excited. However, if we introduce the phase value as part of the initial condition, we only excite one of the branches depending on the chosen phase. For example, if we use the initial condition corresponding to the acoustic branch, i.e., $un,Ni(0)=cos[kNiL+\varphi nAc]$, $u\u0307n,Ni(0)=\u2212\omega 1(k)sin[kNiL+\varphi nAc]$, then we are able to excite only the acoustic branch as shown in Fig. 10, and vice versa.

We note that for a superlattice of two identical masses per unit cell with no spring constants modulation, i.e., $\beta 1=\beta 2=\beta 0$, there is no bandgap and the band structure consists of two folded bands. The dispersion relations for the two bands are $\omega 1=2\beta 0/m|\u2009sin\u2009(kL/2)|$ and $\omega 2=2\beta 0/m2\u2212sin2(kL/2)$, where $kL=[\u2212\pi ,\pi ]$, and the subscript $j$ indicates the branch order. The complex amplitudes are

Hence, the equation for Berry connection is

which is purely real with positive value and therefore a trivial Berry phase of $0$.

### B. System with $Nm$ unit cell

After validating the SAAP method by calculating the elastic wave band structure and the phases for the two mass unit cell, both analytically and numerically, we now explore it further by analyzing more complicated systems with $Nm>2$ with spatial variation of linear stiffness defined by $\beta 1$, $\beta 2,\u2009\u2026,\beta Nm$.

#### 1. Analytical study

The equations of motion for $Nm$ masses in unit cell $Ni$ interacting with spring constants $\beta 1,\beta 2,\u2026,\beta Nm$ forming an infinite 1D chain are

From Eq. (25) we note that the system preserves time reversal symmetry, i.e., $t\u2192\u2212t$; however, by choosing the values of $\beta n,Ni$, the inversion symmetry can be broken. Hence, we will study the traveling waves in this infinite system of masses under different conditions, including broken inversion symmetry.

The traveling wave solutions of Eq. (25) of the compact form $un,Ni(t)=AneikNiLei\omega t;n=1,2,\u2026,\u2009Nm$ lead to the system of equations

where $A\u2192k=[A1,A2,\u2026,ANm]T$,

and $\alpha pq=\u2212\beta p\u2009\u2212\beta q+m\omega 2$. Taking the complex conjugate of Eq. (26) and rearranging the rows and columns, we have

where $A\u0303\u2192k=[ANm*,ANm\u22121*,\u2026,A1*]T$ and

By comparing Eq. (27) with Eq. (26) we see that if $\beta i=\beta Nm\u2212i$, then $M\u2194k=M\u0303\u2194k$ and hence we have $Ai=ANm\u2212i+1*;$$i=1,2,..,Nm/2$. Using Eq. (4) of the Berry connection, the discrete representation of the Berry connection for the system with $Nm$ mass unit cell is

which is purely real with either positive or negative value. Therefore, the Berry phase should be either $0$ or $\pi $ if $\beta i=\beta Nm\u2212i\u2009,\u2009\u2009i=1,2,..,Nm/2$. However, if $\beta i\u2260\u2009\beta Nm\u2212i$, then comparison between Eqs. (27) and (26) does not lead to an explicit relation between the amplitudes of the traveling waves. Therefore, the Berry phase can take values, including those different from $0$ or $\pi $, as was seen from Eqs. (4) and (7). Note that the conditions leading to a Berry phase $0$ or $\pi $ define systems with inversion symmetry.

#### 2. Numerical study

We consider elastic superlattices containing three and four masses per unit cell with different combinations of variations in the spatial stiffness. The dynamical trajectories generated by the MD simulation are analyzed using the SAAP method.

##### a. Three mass unit cell.

From the general analysis, we see that if $\beta i=\beta Nm\u2212i,$ where $i=1,2,\u2026,\u2009Nm/2$, the Berry phase is always either 0 or $\pi $. An arrangement for which the three mass unit cell will have either 0 or $\pi $ Berry phase is that $\beta 1$ equal $\u2009\beta 2$. In the numerical simulation of the system with a three mass unit cell, we first explore two contrasting combinations of elastic stiffness coefficients: (i) $\beta 1=\beta 2\u2009,\beta 3>\beta 1$, and (ii) $\beta 1=\beta 2\u2009,\beta 3<\beta 1$. We show subsequently that when the Berry phase of case (i) is $\pi $, then the Berry phase of case (ii) is $0$.

In Figs. 11 and 12, we plot the normalized complex amplitude and phase (in radians) for three mass unit cell with $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$ and $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$. From Figs. 11(a) and 12(a), we observe that although the real components of the complex amplitudes are equal for masses 1 and 3, the magnitude of the imaginary components are equal but of opposite signs. As a result, from Figs. 11(b) and 12(b) we note that the magnitude of the phases of masses 1 and 3 are the same but opposite in sign. This relationship likely arises due to the symmetry of the system as mass 1 is connected with springs $\beta 1$ and $\beta 3$ to its right and left, respectively, and mass 3 is connected with springs $\beta 3$ and $\beta 1$ to its right and left, respectively.

Figures 13 and 14 show the Berry connection values and manifolds both for $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$ and $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$ cases. Using Eq. (4) of Berry connection, it is clear from the Berry connection plot that the Berry phase value is $\pi ,0,\pi $ for acoustic, first optical, and second optical branches, respectively, for $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$ [Fig. 13(a)], and $0$ for all branches in the case $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$[Fig. 14(a)]. This likely marks a topological transition point where a gap closing and reopening process could be seen if $\beta 3$ is tuned continuously from $0.5\beta 1$ to $2\beta 1$. A previous study has shown this sort of change in the Zak phase of bands with system parameters in a superlattice composed of pipes with changing diameters (Xiao *et al.*, 2015). Figures 13(b) and 14(b) show real and imaginary representations of the manifold generated by the evolution of the amplitude unit vector, $A\u0303\u2192$, along the Brillouin zone. We take $A\u03033,j(k)$ along the normal to the plane of the *k-*space ring, $A\u03032,j(k)$ along the radius, and $A\u03031,j(k)$ along the tangent. The $A1=A3*$ constraint leads to the $Re(A\u0303\u2192)$ vectors forming a 45° angle with the normal to the ring plane, and hence there are only twists but not a change in the direction of $A\u0303\u2192$. Similar to two mass unit cell, the three mass unit cell's amplitude unit vector in Fig. 11(a) (for $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$) and Fig. 12(a) (for $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$) generates a manifold with the imaginary part taking the form of a closed ribbon with either: (i) a single twist at *k *=* *0 [see Fig. 13(b)] or (ii) two twists, one at *k *=* *0 and the other at the boundaries of the Brillouin zone [see Fig. 14(b)]. The total accumulated geometric phase over the Brillouin zone is therefore either $\pi $ (for $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$) or 0 (for $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$).

Now we focus on the case when $\beta i\u2260\beta Nm\u2212i,$ e.g., $\beta 1=2\beta 2\u2009,\beta 2=\beta 3$. Under the conditions described in Sec. IV B 1, we can conclude that the Berry phase need not be restricted to either $0$ or $\pi $. Figure 15 shows the normalized complex amplitudes and phases, and Fig. 16 shows the Berry connection and manifolds for such a case. From Figs. 15(a) and 15(b) we observe that the amplitudes, and hence the phases, of each mass do not appear to have any explicit relation among them. As a result, as is seen from the Berry connection plot [Fig. 16(a)], the Berry phase for each branch does not take on values that are multiples of $\pi $. However, though the Berry phase for any individual band may not be a multiple of $\pi $, we find that the sum of the Berry phases over all bands is an integer multiple of $2\pi $ (similar to Rudner *et al.*, 2016), i.e., $\u22120.83+1.09\u22120.25\u22450(2\pi )$. Furthermore, in contrast to the cases of $\beta 1=\beta 2\u2009,\beta 3=2\beta 1$ or $\beta 1=\beta 2\u2009,\beta 3=0.5\beta 1$, for three mass unit cell with $\beta 2=\beta 3\u2009,\beta 1=2\beta 2$, the amplitude $Re(A\u0303\u2192)$ vectors do not always form a 45° angle with the normal to the Brillouin zone ring plane. Hence, in addition to twists, there is also a change in the direction of $A\u0303\u2192$ [see Fig. 16(b)]. Therefore, although the amplitude unit vector generates a manifold with the imaginary part taking the form of a closed ribbon with two twists at *k *=* *0 and at the boundaries of the Brillouin zone [see Fig. 16(b)], the real part now also forms a closed ribbon with a warp (partial twist). This manifold results in a Berry phase with an arbitrary value.

The two arrangements of the spring constants $(\beta 1,\beta 2,\beta 3)=(1,1,2)$ and $(\beta 1,\beta 2,\beta 3)=(2,1,1)$ are related by a shift in origin. Starting with the three mass unit cell with the sequence of spring constants $(\beta 1,\beta 2,\beta 3)=(1,1,2)$, if one moves the origin from mass 1 in one unit cell to mass 3 in a neighboring unit cell to the left of the original cell, one obtains a unit cell with the sequence $(\beta 1,\beta 2,\beta 3)=(2,1,1)$. This transformation is represented mathematically by the unitary matrix

The Berry connection for the $(\beta 1,\beta 2,\beta 3)=(2,1,1)$ unit cell, $BC(2,1,1)$, relates to that of the original unit cell, $BC(1,1,2)$, through the relationship: $BC(2,1,1)\u2243BC(1,1,2)\u2212iL\Delta k|A1|2$. Finally, after summing over the Brillouin zone, the Berry phase becomes, $\varphi B,(2,1,1)\u2243\varphi B,(1,1,2)+L\Delta k\u2211k=1Nc|A1(k)|2.$ Starting with the (1,1,2) system and after calculating numerically $\u2211k=1Nc|A1(k)|2$ for all three branches, we can predict the Berry phase calculated numerically using the SAAP method for the (2,1,1) system reported in Fig. 16(a).

##### b. Four mass unit cells.

The analytical study above has shown that if $\beta i=\beta Nm\u2212i,$ where $i=1,2,\u2026,\u2009Nm/2$, the Berry phase is always either 0 or $\pi $. Hence, in the numerical calculations for the four mass unit cell, we use both $\beta 1=\beta 3$ and $\beta 1\u2260\beta 3$ cases. Calculated values of the Berry phase for a representative number of four mass unit cell are shown in Table I. From Table I, we see that depending on the variations in the spatial stiffness, different topologies can be seen based on the different values of the Berry phase.

$(\beta 1,\beta 2,\beta 3,\beta 4)$ . | Berry phase . | ||||
---|---|---|---|---|---|

Acoustic branch . | First optical branch . | Second optical branch . | Third optical branch . | Summation . | |

(1,1,1,1/2) | 0 | 0 | 0 | 0 | 0 |

(1/2,1,1,1) | 1.41 | 2.15 | 1.92 | 0.69 | $6.17\u22452\pi $ |

(1,1/2,1,1) | $\pi $ | $\pi $ | $\pi $ | $\pi $ | $4\pi $ |

(1,1,1/2,1) | −1.41 | −2.15 | −1.92 | −0.69 | $\u22126.17\u2245\u22122\pi $ |

$(\beta 1,\beta 2,\beta 3,\beta 4)$ . | Berry phase . | ||||
---|---|---|---|---|---|

Acoustic branch . | First optical branch . | Second optical branch . | Third optical branch . | Summation . | |

(1,1,1,1/2) | 0 | 0 | 0 | 0 | 0 |

(1/2,1,1,1) | 1.41 | 2.15 | 1.92 | 0.69 | $6.17\u22452\pi $ |

(1,1/2,1,1) | $\pi $ | $\pi $ | $\pi $ | $\pi $ | $4\pi $ |

(1,1,1/2,1) | −1.41 | −2.15 | −1.92 | −0.69 | $\u22126.17\u2245\u22122\pi $ |

The evolution of the arrangement of the spring constants starting from $(\beta 1,\beta 2,\beta 3,\beta 4)=(1,1,1,1/2)$ to (1,1,1/2,1) can again be understood as an effect of a change in origin. Similar to three mass unit cell, if we move the origin from mass 1 to mass 4 by shifting the origin to the left by an inter-mass spacing, we obtain $(\beta 1,\beta 2,\beta 3,\beta 4)=(1/2,1,1,1)$ configuration unit cell. Such a change of origin is described by the unitary matrix $S\u2192k$,

Upon this translation, the Berry connection becomes $BC(1/2,1,1,1)\u2243BC(1,1,1,1/2)\u2212iL\Delta k|A1|2$, and Berry phase transforms as $\varphi B,(1/2,1,1,1)\u2243\varphi B,(1,1,1,1/2)+L\Delta k\u2211k=1Nc|A1(k)|2$. For $(\beta 1,\beta 2,\beta 3,\beta 4)=(1,1,1,1/2)$ spring arrangement, we find numerically $|A1(k)|2=5.3,\u2009\u20098.5,\u2009\u20097.4,$ and $2.6$ for acoustic, first optical, second optical, and third optical braches, respectively. With $L\Delta k=\pi /12$ and $\varphi B,(1,1,1,1/2)=0$ for all branches (as shown in Table I, row 1), we find $\varphi B,(1/2,1,1,1)\u22431.39,\u2009\u20092.23,\u2009\u20091.94,$ and $0.68$ for the acoustic, first optical, second optical, and third optical branches, respectively. These predicted values match the numerical values reported in Table I, row 2, up to numerical errors.

Now, if we move the origin from mass 1 to mass 3 by shifting the origin to the left by two inter-mass spacings, we obtain the arrangement $(\beta 1,\beta 2,\beta 3,\beta 4)=(1,1/2,1,1)$. The unitary matrix corresponding to this translation is

The Berry connections are related via the relationship $BC(1,1/2,1,1)\u2243BC(1,1,1,1/2)\u2212iL\Delta k(|A3|2+|A4|2).$ From the analytical study we observed that if $\beta 1=\beta 3$, then $A1=A4*$ and $A2=A3*$, therefore $|A3|2+|A4|2=1/2$. Hence, the Berry phase becomes $\varphi B,(1,1/2,1,1)\u2243\varphi B,(1,1,1,1/2)+\pi .$ Since $\varphi B,(1,1,1,1/2)=0$ for all branches, hence, $\varphi B,(1,1/2,1,1)$ should be $\pi $ for all branches, as is verified numerically in Table I, row 3.

Finally, translating the origin from mass 1 to mass 2 by shifting the origin to the left by three inter-mass spacings, we obtain the $(\beta 1,\beta 2,\beta 3,\beta 4)=(1,1,1/2,1)$ unit cell. The unitary matrix is now

The Berry connection transforms according to $BC(1,1,1/2,1)\u2243BC(1,1,1,1/2)\u2212iL\Delta k(|A2|2+|A3|2+|A4|2)\u2243BC(1,1,1,1/2)\u2212iL\Delta k+iL\Delta k|A1|2$. Hence, the relation between the Berry phases becomes $\varphi B,(1,1,1/2,1)\u2243\varphi B,(1,1,1,1/2)\u2212L\Delta k\u2211k=1Nc|A1(k)|2=\u2212\varphi B,(1/2,1,1,1)$. Therefore, the Berry phase values for the $(1/2,1,1,1)$ unit cell and the $(1,1,1/2,1)$ unit cell should be the same but opposite in sign, as is shown in Table I, rows 2 and 4. Finally, in the case of rows 2 and 4, since $\beta 1\u2260\beta 3$, we obtain Berry phase values which are not multiples of $\pi $, however, the summation of Berry phases equals a multiple of $2\pi $.

## V. CONCLUSIONS

In this paper, we have introduced the SAAP method for the calculation of the amplitudes and phases of elastic waves in general elastic structures. Although tested here for 1D superlattice systems, it is easily generalizable to two-dimensional (2 D) and three-dimensional (3 D) systems. The method entails the use of MD twice with differing initial conditions to allow for the extraction of the band structure, complex amplitudes, and the Berry phase of all the bands in the elastic band structure. We have tested and applied the method to 1D harmonic chains with periodically repeating unit cells containing $Nm=2,\u2009\u20093,$ and 4 elastically coupled masses. The masses are coupled via elastic springs of various stiffness $\u2009\beta i$. We focused on the dependence of the Berry phase of the various elastic bands upon variation of the arrangement and symmetry of the spring stiffness within the unit cell of the superlattice. Various configurations of the unit cell lead to different Berry phases. From our analytical and numerical study, we found that there is inversion symmetry for $Nm=2$ mass unit cell and $Nm>2$ if $\beta i=\beta Nm\u2212i,$ where $i=1,2,\u2026,\u2009Nm/2$. In those cases, the Berry phase is either 0 or $\pi $. We also observe that for $Nm>2$ if $\beta i\u2260\beta Nn\u2212i$, the Berry phase of each band is no longer a multiple of $\pi $, however, the summation of the Berry phase over all bands is an integer multiple of $2\pi $. We also show that unit cells that may be related by a translation of the origin may possess different values of the Berry phase. However, those phases can be related through straightforward operations involving unitary transformations. The eigenvalues of the dynamical problem, i.e., the band structure, are independent of these unitary transformations. All possible normalized eigenvectors related by unitary transformations form a hypersphere in the $Nm$-dimensional complex space of the amplitudes. The eigenvectors of a given choice of Ansatz describe a *k*-parametric trajectory on top of the hypersphere. These eigenvectors generate a manifold, which may or may not be conserved under unitary transformations. In the case of the general and compact Ansätze, the manifold is conserved, which explains the Berry connection and Berry phase invariance. In contrast, Ansätze connected by a unitary transformation representing a change of the origin of coordinates in real space, define different trajectories on the hypersphere as parametric functions of *k*, yielding different Berry connections and Berry phases.

In summary, the SAAP method is a useful tool to explore the topological characteristics of elastic waves in elastic structures, including continuum systems such as rods, beams, and plates. Indeed, topological elastic structures are the subject of very active exploration currently and the ability to calculate the phases of various band structures for arbitrary configuration will accelerate the pace at which these explorations can be undertaken. The phase of amplitudes for elastic systems is intimately related to topological characteristics, including non-reciprocal wave propagation, pseudospins, and phonons that possess fermionic and “imaginary masses” (Deymier *et al.*, 2018). We anticipate that the SAAP method will serve as a practical tool for the computational investigation of the phase characteristics of general elastic structures. In fact, we will apply it next to the investigation of systems with spatio-temporal modulations in the stiffness as well as continuous rods.

## ACKNOWLEDGMENTS

This work was supported by National Science Foundation Award Emerging Frontiers in Research and Innovation (EFRI) No. 1640860. We also would like to thank Nicholas Boechler at University of California at San Diego for useful suggestions.

### APPENDIX

##### 1. Ansätze invariance of the Berry connection and Berry phase

The traveling wave solutions of Eq. (25) for the Ansatz of the general form $un,Ni=An\u2032eikxneikNiLei\omega t;\u2009\u2009n=1,2,\u2026,\u2009Nm$ lead to the system of equations

where $A\u2192k\u2032=[A\u20321,A\u20322,\u2026,ANm\u2032]T$ and

Here $xpq=xp\u2009\u2212xq$ and $\alpha pq=\u2212\beta p\u2009\u2212\beta q+m\omega 2$ as defined before. The transformation linking the amplitudes, and therefore Eqs. (26) and (A1), is given by the unitary matrix

as it verifies the identities $S\u2194k\u2020M\u2194k\u2032S\u2194k=M\u2194k$ and $A\u2192\u2032k=S\u2194kA\u2192k$.

Based on the two Ansätze, we can either use the amplitudes $An$ of the *compact Ansatz* $un,Ni(t)=AneikNiLei\omega t$ or $An\u2032eikxn$ of the *general Ansatz* $un,Ni=An\u2032eikxneikNiLei\omega t$ to calculate the Berry connection and the Berry phase. If we use $An\u2032eikxn$ or in the general form $S\u2194k\u2020A\u2192\u2032k,$ $S\u2194k\u2020=diagonal(eikx1,eikx2,\u2026,eikxNm),$ to calculate the Berry connection, then in virtue of the unitary character of $S\u2194k$, we have

Therefore, the Berry connections and hence Berry phases are the same regardless of the choice of Ansatz.

## References

*Springer Series in Solid-State Sciences*(

*Springer Series in Solid-State Sciences*(