Low-frequency tunable topological interface states in soft phononic crystal cylinders

Topological phononic crystals have attracted intensive attention due to their peculiar topologically protected interface or edge states. Their operating frequency, however, is generally fixed once designed and fabricated. Here, we propose to overcome this limitation by utilizing soft topological phononic crystals. In particular, we design a simple one-dimensional periodic system of soft cylindrical waveguides to realize mechanically tunable topological interface states for longitudinal waves. To this end, we employ the nonlinear elasticity theory and its linearized incremental version to fully account for both geometric and material nonlinearities of the system. We derive the dispersion relation for small-amplitude longitudinal motions superimposed on the finitely deformed state. In addition, our analytical results provide information about the corresponding Bloch wave modes, displacement field distributions, and signal transmission coefficients for finite cylindrical waveguides with identical or various unit-cell topologies. Our numerical results illustrate the low-frequency topological interface state occurring at the interface between two topologically distinct soft phononic cylinders. Moreover, we show that the corresponding frequency in the overlapped band gap can be continuously adjusted by an external force. This analytical result is also validated by the finite element simulations. Finally, we provide the topological phase diagrams to demonstrate the tunable position and existence condition of the topological interface states when tuning the external loading. The low-frequency tunable topological interface states with remarkable field enhancement may find a wide range of potential applications such as tunable energy harvesters, low-pass filters and high-sensitivity detectors for biomedical applications.


Introduction
Phononic crystals (PCs) have attracted intensive attention thanks to their outstanding properties in the manipulation of acoustic/elastic waves. The intrinsically artificial periodic composites can give rise to the wave band gap (BG) -a special state, where acoustic/elastic waves are prohibited within a certain frequency range. This prominent characteristic of PCs can be attributed to the Bragg scattering (Kushwaha et al., 1993), local resonance (Liu et al., 2000) and inertial amplification (Yilmaz et al., 2007). The unique BG character and strong dispersive properties in passbands may produce anomalous wave behaviors such as acoustic/elastic wave filtering, focusing, directional propagation, negative refraction and cloaking (Hussein et al., 2014;Ma and Sheng, 2016;Cummer et al., 2016;Slesarenko et al., 2018;Krushynska et al., 2018).
Recently, the topic of topological acoustic or mechanical PCs has emerged to offer exciting opportunities for designing materials with broadband one-way transport properties. These peculiar material systems can generate topologically protected unidirectional, backscatteringimmune interface or edge states (Zhang et al., 2018;Ma et al., 2019). The theoretical development has motivated some experimental efforts in realization of the topological PCs (Süsstrunk and Huber, 2015;Xiao et al., 2015;Peng et al., 2016;Lu et al., 2017;Yu et al., 2018;Yin et al., 2018;Li et al., 2018;Ding et al., 2019), indicating the feasibility of the concept for potential applications in wave filters, energy harvesters, acoustic rectifiers, vibration isolators, acoustic imaging, and bio-sensors. The topological characteristics of PC band structures can be characterized by their topological invariants such as the Berry phase (Zak, 1989) for two-dimensional (2D) systems or Zak phase (Atala et al., 2013) for one-dimensional (1D) media. Since topological PCs possess the global properties of band structures, their nontrivial topological states are extremely robust to the defects and boundary effects (Xiao et al., 2014;Wang et al., 2015).
Topological PCs can be categorized into the following three classes. The first approachan analogue of the quantum Hall effect -is to break the time-reversal symmetry of the systems and realize the topologically protected edge states through introducing gyroscopic inertial effects (Wang et al., 2015), external flow fields (Chen and Wu, 2016) or time-modulated materials (Chen et al., 2019a). This type of topological PCs -referred to as the acoustic/mechanical Chern insulators -has been experimentally verified by Ding et al. (2019). This method, however, is still challenging to be applied in real-life situations, due to the complicated implementation of external uniform motion in the lattice and the inherent dynamic instabilities and noise in a moving medium. The second strategy exploits the pseudospindependent edge states, and breaks the spatial-inversion symmetry of the time-reversal invariant topological PCs, also known as quantum spin Hall topological insulators Xia et al., 2017;Chen et al., 2018b). These topological PCs are intrinsically based on the spin-orbit coupling mechanism -an analogue to the quantum spin Hall effect. In principle, the quantum spin Hall topological insulators have a double Dirac degeneracy in band structures and the topological phase inversion appears at the double Dirac point, where the pseudospin state forms the topological edge state (Deng et al., 2019). Depending on the appropriate polarization excitation, the quantum spin Hall topological insulators support robust forward or backward edge states due to their time-reversal symmetry (Yang et al., 2018;Yu et al., 2018). The last method is based on the quantum valley Hall effect, which provides a pair of valley vortex states with opposite chirality (Lu et al., 2017;Pal and Ruzzene, 2017;Zhu et al., 2018;Wang et al., 2020b). The quantum valley Hall PCs also possess the topologically protected one-way edge states along the interface of two domains with different valley vortex states. A systematic comparison between acoustic topological states based on valley Hall and quantum spin Hall effects can be found in the recent work by Deng et al. (2019). In the 1D phononic systems, a combination of PCs with different topological properties may result in topological interface states within the overlapped BGs; see, for example, the observations and predictions in discrete spring-mass , water wave (Yang et al., 2016), acoustic (Xiao et al., 2015) and elastic (Yin et al., 2018) systems. For more detailed discussion of the recent progress in the topological acoustic/mechanical sys-tems, interested readers are referred to the comprehensive review articles by Zhang et al. (2018) and Ma et al. (2019).
A major limitation of passive topological PCs is that their operating frequency range of topological states is fixed and extremely narrow (usually corresponding to a single frequency of transmission peak in the BG). To realize a wider operating frequency range of topologically protected states, several methods have been proposed to design actively tunable topological PCs. By adjusting the airflow velocity field and unit-cell geometric size, the time-reversal symmetry of a 2D PC was broken to tune the BG topology and realize a tunable topological acoustic Chern insulator (Chen and Wu, 2016). The intelligent magnetoelastic materials were introduced by Feng et al. (2019) into the topological system to realize magnetically tunable topological interface states for Lamb waves in 1D PC slabs. The periodic electric boundary conditions were exploited by Zhou et al. (2020b) to generate actively tunable topologically protected interface mode in a 1D piezoelectric rod system. Wang et al. (2020a) investigated the topological interface mode in a 1D granular PC composed of two sub-lattices, which can be tuned by varying the pre-compression between the spheres. Although all the above-mentioned works demonstrate the tunability of topologically protected edge/interface states, they operate in the high-frequency scenarios.
Soft PCs offer both the low-frequency operating ranges and high tunability by external stimuli such as pre-stretch (Bertoldi et al., 2017;Galich et al., , 2018Wu et al., 2018a;Li et al., 2019;Parnell and De Pascalis, 2019;Gao et al., 2019), electric Wu et al., 2017Wu et al., , 2018bMao et al., 2019;Chen et al., 2020;Wang et al., 2020c;Zhu et al., 2020) or magnetic field (Karami Mohammadi et al., 2019). These abilities motivate the exploration of low-frequency tunable topological states in soft PCs. By changing the filling ratio and tuning the mechanical load, the dynamically tunable topological interface state was experimentally observed by Li et al. (2018) in a circular-hole soft PC plate made of two domains with different topological properties. The design of a 2D quantum valley Hall PC was presented by Liu and Semperlotti (2018), where the topological states at the domain interface are triggered by geometric nonlinear effects due to the applied strain. These two works, however, neglect changes in material stiffness induced by the pre-stretch. The influence of the essential material nonlinearity for soft matter was considered by Nguyen et al. (2019) in the context of soft topological PCs. They designed a 2D quantum valley Hall PC consisting of soft annular cylinders embedded in an elastic matrix, and utilized the prestretch and inflation to actively tune the frequencies of topologically protected edge states. Zhou et al. (2020a) designed a soft membrane-type PC for the voltage-controlled quantum valley Hall effect in a dielectric elastomeric membrane with sprayed metallic particles. Recently, Huang et al. (2020) examined a 1D soft periodic system composed of topologically different plates and realized tunable topological interface states by applying deformation. Nevertheless, the geometric structures and loading ways in the aforementioned soft topological PCs are relatively complex, and it is not easy to combine two differently deformed domains while keeping a smooth interface.
It is well known that the Bragg band gaps (BGs) are usually produced by the material periodicity, geometric periodicity, and periodic boundary conditions (Hussein et al., 2014). With appropriate geometric and/or material design of metamaterials, their wave propagation and attenuation behaviors (such as transmission and reflection) could be optimized (Derın et al., 2019;Abdulkarim et al., 2020;Ozturk et al., 2020a,b). Inspired by previous works (Xiao et al., 2015;Yin et al., 2018), here we design a 1D soft phononic crystal cylinder (PCC) composed of step-wise sub-cylinders to realize low-frequency tunable topological interface states under the application of mechanical load. We fully account for both the geometric and material nonlinearities in our theoretical and numerical analyses. The proposed soft topological PCC is a single-phase material structure allowing to induce different deformation states in its base elements while preserving a smooth interface between the base elements. Due to the low stiffness of soft materials, their operating frequency range is much lower than that of hard materials with the same structure. We focus on (i) tunability of the topological interface states (in the low-frequency range) by an applied axial force, and (ii) influence of the strain-stiffening effect on the tunability and existence of the topological interface states. To this end, we derive analytically the dispersion relations and acoustic characteristics for small-amplitude longitudinal waves propagating in the finitely deformed PCC. This information is complemented by our numerical calculations including finite element (FE) simulations, elucidating the relations between the morphology, applied loading and material nonlinearity effects on the band structures, transmission characteristics and topological phase diagrams. This paper is organized as follows. The theoretical background on nonlinear elasticity theory and its associated incremental theory (Ogden, 1997) is summarized in Sec. 2. The nonlinear static response of the proposed soft PCC with alternating cross-sections is analyzed in Sec. 3. Section 4 describes the derivations of the dispersion relation, transmission coefficient and displacement field of a finite PCC waveguide with various periodic unit cells. Numerical calculations are described in Sec. 5. For a mixed finite neo-Hookean PCC waveguide, the frequency of topological interface states is lowered monotonically by the increasing axial force. However, for a Gent PCC waveguide, the axial force affects the topological interface state frequency in a non-monotonous way that an increase in the axial force leads to the continuous decrease of frequency to a minimum value, and then the frequency is increased reversely by a further increase of the axial force. Section 6 concludes the work with a summary and discussion. Some mathematical derivations and FE simulation procedures are provided in Appendices A-C.

Nonlinear elasticity
We consider a deformable continuous body that occupies the undeformed reference configuration B r in the Euclidian space with the boundary ∂B r and the outward unit normal N. An arbitrary material point labelled as X in the undeformed configuration is identified by the position vector X. Subjected to a mechanical loading, the body deforms and moves to the deformed or current configuration B t with the boundary ∂B t and the outward unit normal n t , such that the point X occupies a new position x = χ(X, t) at time t in B t , where an invertible vector function χ is defined for all points in B r . The deformation gradient tensor is defined as F = ∂x/∂X = Gradχ, where 'Grad' is the gradient operator with respect to B r . The components of the deformation gradient tensor are F iα = ∂x i /∂X α , where Roman and Greek indices are associated with B t and B r , respectively. The local measure of the volume change is denoted by J = det F > 0.
In the absence of body force, the equilibrium equations can be written in Eulerian and Lagrangian forms, respectively, as where τ = J −1 FT is the Cauchy stress tensor and T is the nominal stress tensor. Here 'div' and 'Div' denote the divergence operators relative to B t and B r , respectively. Note that the nominal stress tensor is the transpose of the first Piola-Kirchhoff stress tensor and that both of them are non-symmetric two-point tensors like the deformation gradient tensor. In index notation, the equilibrium equations (1) read τ ij,i = 0 and T αi,α = 0, where the Einstein summation convention is used.
Consider a compressible hyperelastic material described in terms of its strain energy density function Ω (F) (per unit reference volume) such that or in index notation with the angular momentum conservation τ ij = τ ji . Alternatively, the strain energy density function Ω can be expressed in terms of the principal stretches, i.e., Ω = Ω(λ 1 , λ 2 , λ 3 ) with J = λ 1 λ 2 λ 3 (Ogden, 1997). Thus, referring to the principal axes of τ , the corresponding principal Cauchy stresses τ i (i = 1, 2, 3) are expressed as The mechanical boundary conditions on ∂B t can be written in Eulerian form as where t a is the applied mechanical traction vector per unit area of ∂B t .

Incremental motions superimposed on finitely deformed state
A time-dependent infinitesimal incremental motionẋ (X, t) is superimposed on a finitely deformed configuration B 0 (with the boundary ∂B 0 and the outward unit normal n). Here, the incremental quantities are represented by a superposed dot. The incremental equation of motion in the updated Lagrangian form is whereṪ 0 = J −1 FṪ is the push-forward counterpart of the Lagrangian incremental stress tensorṪ, u =ẋ (X, t) is the incremental displacement vector, and ρ = ρ 0 J −1 is the current mass density in B t , with ρ 0 denoting the mass density in the reference configuration B r . The subscript 0 indicates the resulting push-forward quantities and the subscript t following a comma represents the material time derivative. The linearized incremental constitutive law for a compressible hyperelastic material iṡ where H = gradu is the incremental displacement gradient tensor; 'grad' is the gradient operator with respect to B 0 . The fourth-order instantaneous elasticity tensor A 0 is represented in component form by in which A indicates the referential elasticity tensor with its components given by Following Ogden (1997) and referring to the principal axes of τ , the non-zero components of A 0 for compressible isotropic hyperelastic materials can be expressed, in terms of the three principal stretches λ i , as where Ω i = ∂Ω/∂λ i and Ω ij = ∂ 2 Ω/∂λ i ∂λ j . In the updated Lagrangian form, the incremental mechanical boundary conditions, which are to be satisfied on ∂B 0 , are written asṪ where the superscript (·) T signifies the usual transpose of a tensor andṫ A 0 is the updated Lagrangian incremental mechanical traction vector per unit area of ∂B 0 .

Nonlinear deformation of a soft PCC
Consider a single phase hyperelastic PCC structure with periodically varying crosssectional areas as shown in Fig. 1(a). Each unit cell has two wider sub-cylinders 1 and 3 of length L (1) /2 and inner radius R (1) 0 , sandwiching a narrower sub-cylinder 2 of length L (2) and inner radius R (2) 0 . Here and thereafter, the superscript (·) (p) denotes the physical quantities of the sub-cylinder p (p = 1, 2, 3). The wider sub-cylinders p = 1 and 3 have identical geometric sizes (i.e., R 0 and L (1) = L (3) ). In the undeformed configuration, the unit-cell length is L = L (1) + L (2) along the Z direction. Note that the 1D PCC with inversion symmetry has two inversion centers and without loss of generality, we assign the origin to be at the center of the wider sub-cylinder (see Fig. 1). By varying the initial length fraction of sub-cylinder 1 or 2, the Bragg BG could exhibit the evolutionary process of open, close and reopen. The topological transition point where the BG closes is mechanically tunable, which is our main goal and is shown in Sec. 5.
As shown in Fig. 1(b), under the application of tensile axial force, the length of the deformed PCC becomes longer and its lateral size becomes smaller, respectively, than those of the undeformed PCC. The length fraction of the deformed sub-cylinder is different from that of the undeformed one. We note that, due to the geometric inhomogeneity, complex local deformations can develop near the interfaces between the wider and narrower subcylinders when subjected to an axial force F N (Parnell and De Pascalis, 2019). These local deformations, however, are only confined in small regions in the vicinity of the interfaces and barely affect the response of topological interface states in the low-frequency regime of interest here. Therefore, the nonlinear deformation can be approximately assumed uniform in the theoretical model. As we shall show, this assumption has been validated by the FE simulations (see Sec. 5.2). The deformed configurations of the soft PCC and its unit cell are shown in Fig. 1(b). The uniform axisymmetric deformations can be described in two cylindrical coordinate systems (R, Θ, Z) and (r, θ, z) as follows: where λ 1 and λ 3 are the principal stretches along the radial and axial directions, respectively. Thus, the geometric sizes of each sub-cylinder become where λ (p) 1 and λ (p) 3 represent the principal stretches of sub-cylinder p, and r (p) 0 , l (1) /2 = l (3) /2 and l (2) are the radius of sub-cylinder p and the lengths of the wider and narrower subcylinders in the deformed state, respectively. In addition, l = l (1) + l (2) is the length of the deformed unit cell. Therefore, the deformation gradient tensor of sub-cylinder p is expressed as 3 ], which will be determined by the initial boundary conditions on the lateral surface.
In order to analyze the longitudinal wave propagation in the soft PCC, the compressible Gent material model (Gent, 1996) is adopted to characterize the hyperelastic cylinder, which is described as where µ and Λ are the shear modulus and the first Lamé's parameter in the undeformed configuration, I denotes the first strain invariant and J (p) = (λ 3 . The bulk modulus is then calculated as K = Λ+2µ/3. The parameter J m is the dimensionless Gent constant used to characterize the strain-stiffening effect of the PCC. Recall that the compressible neo-Hookean model is recovered from Eq. (15) when J m → ∞. It should be emphasized that the material properties (i.e., µ, Λ, K and J m ) are the same for the three sub-cylinders, but I (p) 1 and J (p) are different because of various cross-sections. In virtue of Eqs. (5) and (15), we obtain the principal Cauchy stress components for sub-cylinder p as Considering the axial force F N applied along the z direction as well as the traction-free boundary condition on the lateral surface r = r where is the area of the deformed cross-section of sub-cylinder p. Therefore, the nonlinear algebraic equations (16) and (17) can be utilized to completely determine the principal stretch ratios λ 3 ) once the axial force F N is prescribed.

Analysis of incremental longitudinal wave propagation
After obtaining the nonlinear axisymmetric deformations in Sec. 3, the solutions of the superimposed incremental longitudinal waves in an initially deformed PCC are derived in Sec. 4.1. The transfer matrix method (Yeh et al., 1977) in conjunction with the Bloch-Floquet theorem (Kittel, 2005) is then employed in Sec. 4.2 to derive the dispersion relation of incremental wave motions in an infinite PCC, which in turn determines the displacement mode shape of unit cell in Sec. 4.3. Furthermore, the transmission coefficient of a finite PCC with identical unit cells is provided in Sec. 4.4. For a finite cylindrical waveguide consisting of two types of PCCs with different unit cells, we derive its transmission coefficient and displacement distribution in Secs. 4.5 and 4.6, respectively.

Wave solutions of incremental motions
For each sub-cylinder p, the incremental constitutive law (8) for the superimposed longitudinal waves can be expressed in component form aṡ where the non-zero components of the instantaneous elasticity tensor A 0 for the compressible Gent model characterized by Eq. (15) may be derived from Eq. (11) as As is well-accepted in the classical rod theory, the assumption of 1D stress state (Wu et al., 2018b) is made in the following derivation, which results inṪ where where A 2 is the effective elastic stiffness. Thus, Eq. (22) is the reduced incremental constitutive relation in the updated Lagrangian form.
Here, by introducing the incremental axial displacement w, we can rewrite H 33 as H (p) 33 = dw (p) /dz. Due to the postulation of 1D stress state as well as the applied axial force, the incremental equation (7) of motion is simplified, if ignoring the lateral inertial effect, aṡ Note that taking the lateral inertial effect into account will lead to the Love rod theory (Graff, 1991), which exceeds the scope of the present study. Substituting Eq. (22) into Eq. (23) and which is the incremental wave equation for the superimposed longitudinal motions, where all physical fields depend on z and t only.
Consequently, for the harmonic time-dependency e −iωt with ω being the angular frequency, the incremental axial displacement in sub-cylinder p of the nth unit cell (see Fig. 1(b)) can be written as where w n are the undetermined complex coefficients denoting the amplitudes of incident and reflected waves, respectively, and 0 /ρ (p) ) represents the axial wave number in sub-cylinder p. Owing to the difference in cross-sections for various sub-cylinders, it is appropriate to choose the incremental axial force Q (p) n rather than the incremental stressṪ (p) 033n to be continuous at the interfaces delimiting the sub-cylinders. Inserting Eq. (25) into Eq. (22) and multiplying the resultant expression by the deformed cross-sectional area leads to the incremental axial force as where Q is the axial force amplitude.

Dispersion relation for an infinite PCC
For simplicity, the state vector in sub-cylinder p of the nth unit cell (see Fig. 1 The state vectors are not independent of each other and can be connected through the interfacial continuity conditions between different sub-cylinders, which are expressed as

Substituting Eqs. (25)-(28) into Eq. (29) and noting s
0 , we rewrite the displacement and force continuity conditions (29) as where . Through some mathematical manipulations, the transfer relation in Eq. (30) can be expressed as where M is the 2×2 unit-cell transfer matrix that relates the state vector in one sub-cylinder of a unit cell to that in the same sub-cylinder of the adjacent unit cell, and its components are Note that M is a unimodular matrix (Yeh et al., 1977), and we have The deformed PCC is still periodic along the axial direction. Based on the Bloch-Floquet theorem (Kittel, 2005), the relation of state vectors of the same sub-cylinder in adjacent unit cells takes the form as follows: where q is the Bloch wave number. It follows from Eqs. (31) and (34) that the state vector of the Bloch wave satisfies the following eigenvalue problem: As a necessary condition for nontrivial solutions, the determinant of the coefficient matrix of Eq. (35) vanishes, which yields the following dispersion relation as Thus, Eq. (36) determines the relation between q and ω (i.e., the band structure) for incremental longitudinal waves.

Displacement mode shape of the deformed unit cell for Bloch waves
To interpret the topological characteristics of band structures, it is feasible to examine the symmetry properties of mode shapes of passband/BG edge states (Xiao et al., 2014). Thus, we will provide the derivation of displacement mode shapes in the deformed unit cell in this subsection. Corresponding to the eigenvalue e −iql , the eigenvector of the transfer matrix for sub-cylinder 3 of the first unit cell is obtained from Eq. (35) as where S 1 ] T denotes the state vector made of the coefficients of incident and reflected waves in sub-cylinder 3. Making use of the interfacial continuity conditions (29) 2,3 for the displacement in the first unit cell, we have where S (1) 1 ] T and S (2) 1 ] T are the state vectors in sub-cylinders 1 and 2 of the first unit cell.
Once the state vectors S (p) 1 are determined from Eqs. (37) and (38), we obtain the displacement distributions in the deformed unit cell for incremental Bloch waves as

Transmission coefficient of a finite PCC with identical unit cells
Now consider a finite PCC with N identical deformed unit cells arranged in the axial direction (see Fig. 2(a)). Inserting Eq. (28) into Eq. (31), using Eq. (33) and performing matrix transfer N times, we have where a are the amplitude coefficients of incident and reflected waves at the incident side, respectively, and M t is the global transfer matrix.
To calculate the transmission spectrum in the finite hyperelastic PCC, we set the reflection coefficient at the output side to be zero (i.e., b (3) N = 0). As a result, the wave coefficient where M tij are the components of the global transfer matrix. The transmission coefficient t N defined as the absolute square of a 0 is then calculated as 4.5. Transmission coefficient of a finite waveguide with two types of different unit cells In order to investigate the existence of topological interface states, the transmission behaviors of a finite cylindrical waveguide composed of two types of different deformed unit cells (i.e., N unit cells of S1-type and M unit cells of S2-type arranged consecutively in the axial direction) are considered in this subsection, and its schematic diagram is shown in Fig. 2(b). Based on Eq. (40), the transfer relation of N unit cells of S1-type is expressed as where M (1) t is the transfer matrix of N unit cells of S1-type. Analogous to Eqs. (25) and (27), the axial displacement and force amplitudes in the mth unit cell for the S2-type PCC are where z ′ = z − N l − l ′ , (·) ′ stands for the related parameters and physical quantities of the S2-type unit cell, and m is chosen to vary from 0 to M − 1 for the purpose of illustration (see Fig. 2(b)). Referring to Eqs. (31) and (33), we can obtain the transfer relation between the state vectors for sub-cylinder 3 in two adjacent S2-type unit cells as Similar to Eq. (43), the transfer relation of the last M − 1 unit cells of S2-type is where M (2) t indicates the transfer matrix for the last M − 1 S2-type unit cells. Furthermore, the interfacial continuity condition between the two different PCCs and those in the first S2-type unit cell (m = 0) are written as Using Eqs. (25), (27) and (44), Eq. (47) becomes (48) Through some mathematical manipulations, Eq. (48) is rewritten as where M int is the 2 × 2 interfacial transfer matrix between the two different PCCs and its components are omitted here due to their redundancy. Combining Eqs.
where K = M is the 2 × 2 global transfer matrix. Provided that the reflection coefficient at the output side is equal to zero (i.e., b where K ij are the components of K. Note that Eq. (51) has the same form as Eq. (42).
4.6. Displacement field of a finite waveguide with two types of different unit cells In this subsection, we will derive the incremental displacement distribution of a finite deformed waveguide consisting of N unit cells of S1-type and M unit cells of S2-type. In view of the second formula of Eq. (50) Fig. 2 which, combined with the transfer relation (31) between two adjacent unit cells, yields the state vector of sub-cylinder 3 in the nth S1-type unit cell as Utilizing the interfacial continuity conditions (29) 2,3 for the displacement in the nth S1-type unit cell, we obtain the state vectors S (1) n and S (2) n of sub-cylinders 1 and 2 in terms of S n . Similarly, the state vector of sub-cylinder 3 in the mth S2-type unit cell is written as where S ′(3) 0 can be achieved from Eqs. (49) and (53). Analogous to the S1-type unit cell, the state vectors S

Results and discussion
This section will elucidate the tunable effects of mechanical load on the topological interface state of longitudinal waves propagating in the hyperelastic PCC characterized by the neo-Hookean and Gent models, respectively. As described in Sec. 3, the developed theoretical model assumes a uniform nonlinear deformation when subjected to an axial force and neglects the locally nonuniform deformation near the interfaces of different sub-cylinders. Therefore, the effectiveness of this hypothesis will be validated in Subsec. 5.2.2 by performing the FE simulations based on the commercial software package ABAQUS.
In the following numerical calculations, the undeformed unit cell shown in Fig. 1(a) has the radius R (1) 0 = 0.5 cm and the length L (1) = φ 0 L for sub-cylinder 1, along with the corresponding parameters R (2) 0 = 0.4 cm and L (2) = (1 − φ 0 ) L for sub-cylinder 2, where the total length is L = 10 cm and φ 0 is the initial length fraction of sub-cylinder 1. Note that the radii and the total lengths of the undeformed unit cells are the same for the two base elements of the PCC waveguide, which ensures a smooth interface between the two PCC elements after deformation. In addition, the hyperelastic PCC is characterized by the commercial product Zhermarck Elite Double 32 made of silicon rubber  with its initial density, shear modulus and first Lamé's parameter given as ρ 0 = 1040 kg/m 3 , µ = 0.444 MPa and Λ = 22.2 MPa, respectively. The dimensionless axial force is defined as F N = F N /µS (2) , where S (2) is the initial cross-sectional area of sub-cylinder 2. We define the normalized Bloch wave number as q = ql/2π ranging from −0.5 to 0.5 within the first Brillouin zone (Kittel, 2005). The ordinary frequency f measured in Hz is given by f = ω/(2π).

Nonlinear static deformation
First, we examine the nonlinear axisymmetric deformation of the soft PCC under the action of axial force. The analytical results are calculated from the nonlinear algebraic equations (16) and (17). Fig. 3(a) shows the variations of axial (λ 3 ) and radial (λ (2) 1 ) stretch ratios of sub-cylinder 2 with the dimensionless axial force F N for both the neo-Hookean and Gent (J m = 20) models. Indeed, the results for the neo-Hookean model are also recovered by those for the Gent model with a large enough value of J m (e.g., J m = 2000). Clearly, when the axial tensile force is applied, the length of sub-cylinder 2 increases and its radius shrinks owing to the Poisson's effect. The induced stretch variations of the two material models overlap for F N ≤ 1.25. Nevertheless, the difference between the axial stretches of the neo-Hookean and Gent models becomes more evident with an increase in the axial force when F N > 1.25. Specifically, at the same axial force level, the Gent PCC with the strainstiffening effect experiences a smaller deformation as compared to the neo-Hookean PCC. Similar results can be obtained for sub-cylinder 1 in the unit cell. Intuitively, the sub-cylinder with a smaller cross-section will be elongated more than the counterpart with a larger one in the unit cell. For completeness, we plot in Fig. 3(b) the dependence of the overall stretch ratio, defined as on the initial length fraction φ 0 for sub-cylinder 1 with a larger cross-section. The results are shown for the neo-Hookean and Gent PCCs subjected to different axial force levels, namely, F N = 1, 2 and 4. As expected, the PCC with a smaller φ 0 (i.e., a smaller length fraction of sub-cylinder 1 with a larger radius) develops larger deformations when subjected to an identical axial force. As the axial force level is increased, this variation trend strengthens irrespective of the material model. However, this variation trend of the Gent PCC weakens for a large axial force activating the strain-stiffening effect, compared to that of the neo-Hookean PCC (such as F N = 4 in Fig. 3(b)). Note that we have selected the sub-cylinders with a relatively small difference in their initial cross-sectional areas to diminish the influence of the inhomogeneous deformation at the interface, especially for a larger axial force F N . Next we will examine the tunable topological interface states propagating in the neo-Hookean and Gent PCC waveguides separately, to illustrate the effects of geometric and material nonlinearities.  By varying the initial length fraction φ 0 of sub-cylinder 1, the second BG for F N = 0 closes at the center (q = 0) of the first Brillouin zone with a Dirac cone formed by accidental degeneracy (see Fig. 4(b)). The Dirac cone (where the two bulk bands have linear dispersion) occurs at φ 0 = 0.5 in the absence of an axial force. The degeneracy is broken for any configuration such that φ 0 = 0.5, and the Dirac cone will be opened to form the second BG. This is illustrated, if we decrease φ 0 from 0.5 to 0.35 and increase φ 0 from 0.5 to 0.65, in Figs. 4(a) and (c) for the unloaded case. Therefore, the second BG of the soft PCC without axial force exhibits the evolutionary process of open, close and reopen when adjusting the geometric parameter φ 0 . Remarkably, as we shall show below, this transition corresponds to the switching of topological characteristics. Similar BG inversion process is observed for the axially loaded PCC as demonstrated in Figs. 4(d)-(f). A topological transition point (i.e., the point for two bands to cross) is also obtained in Fig. 4(e) for the case of F N = 1, where the Dirac cone appears at a different length fraction φ 0 ≃ 0.524 due to the nonlinear deformation of unit cell.
It should be pointed out that the topological characteristics of a BG is completely determined by the summation of the Zak phases of all the bulk bands below this gap (Zak, 1985;Xiao et al., 2014Xiao et al., , 2015. Originating from electronic systems (Kohn, 1959;Zak, 1989), the so-called Zak phase is a special type of Berry phase, which is a topological invariant characterizing the topological property of bulk bands in 1D periodic systems. The Zak phase for the jth isolated band of the 1D PCC system is defined as (Xiao et al., 2015;Yin et al., 2018) where z is the deformed axial coordinate, r denotes the position in the cross-section plane, l is the length of the deformed unit cell, ρ and c are the current mass density and longitudinal wave velocity in the deformed configuration, respectively, and W j,q (z, r) = w j,q (z, r) e −iqz represents the periodic in-cell part of the normalized Bloch displacement eigenfunction w j,q (z, r) in the jth band with Bloch wave number q. The factor 1/(2ρc 2 ) is the weight function for the elastic system. Thus, to distinguish the topological properties of different PCC configurations, it is necessary to obtain the Zak phase value. Here, we make use of a discretized form (A.3) to numerically calculate the Zak phase (the detailed derivation of the numerical procedure is presented in Appendix A). Given the inversion symmetry of the 1D PCC with respect to its central cross-sectional plane, the calculated Zak phase value is quantized at either 0 or π (Zak, 1989). This quantization also holds true for the PCC under the action of an external mechanical load. The calculated Zak phase values are marked in magenta on the corresponding bands in Fig. 4. Moreover, the BG topological property can be characterized by the BG sign ς, which is related to the Zak phase and given by (Xiao et al., 2014) sgn where the integer h is the number of band crossing points below the jth BG. The second BGs with ς > 0 and ς < 0 are represented by the yellow and green stripes in Fig. 4, respectively, to demonstrate different BG topological properties. For simplicity, the two PCC configurations with φ 0 = 0.35 and φ 0 = 0.65 will be referred to as S1-configuration and S2-configuration, respectively. While the Zak phase of the first band is the same for S1-and S2-configurations in the absence of an axial force (see Figs. 4(a) and (c)), the Zak phase of the second bulk band varies from 0 to π, which exhibits a topological phase transition. According to Eq. (57), the topological properties of the second BG for S1-and S2-configurations are thus different in spite of the overlapped BG frequencies (see Figs. 4(a) and (c)). This indicates the existence of a topological state in this BG at the interface separating S1-and S2-configurations. This topological interface state will be discussed in the next Sec. 5.2.2. The phenomenon is also observed in PCC subjected to an axial loading (see Figs. 4(d)-(f) for F N = 1). In this case, however, the BG frequencies for the loaded PCC are lower than those for the unloaded PCC (compare Figs. 4(d) and (f) with Figs. 4(a) and (c)). For example, the second BG frequency range of S2-configuration is from 336 Hz to 376 Hz for F N = 0, while for F N = 1 the corresponding BG frequency varies from 263 Hz to 297 Hz. Thus, the applied axial force can tune the BG position and the frequency of the topological transition point.
In addition to the direct calculation of Zak phase and BG sign from Eqs. (56) and (57), the symmetry analysis method of the edge states at the two Brillouin zone symmetry points can also be employed to verify the topological phase transition and to identify the BG topological property (Kohn, 1959;Zak, 1989;Xiao et al., 2014). To perform the symmetry analysis, we make use of Eq. (39) and calculate the absolute value of the displacement field w(z) for the six band-edge states A − F (indicated with cross symbols in Figs. 4(d) and (f)). Fig. 5 shows the dependence of |w(z)| on the normalized axial coordinate z * = z/l in the deformed neo-Hookean unit cell with F N = 1. For a bulk band, the Zak phase is θ Zak j = 0 if the two edge states at the Brillouin zone center and boundary possess the same symmetric property. Otherwise, the Zak phase value should be θ Zak j = π. For S1-configuration, the edge states A and B of the second band exhibit the symmetric distributions with respect to the unit-cell center (i.e., even eigenmodes associated with a nonzero displacement amplitude at the unit-cell center) (see Figs. 5(a) and (b)), and hence the corresponding Zak phase is 0. For S2-configuration, however, the edge states D and E of the second band are symmetric and antisymmetric (odd eigenmode related to the zero displacement amplitude at the unit-cell center) (see Figs. 5(d) and (e)), and its Zak phase is π. Therefore, the Zak phase of the second band is altered after the band crossing, which is in full agreement with the topological phase transition shown in Figs. 4(d) and (f). In addition, Figs. 5(b) and (f) show the symmetric displacement fields with respect to the center of unit cell for edge states B and F , while the displacement fields in Figs. 5(c) and (e) are antisymmetric for edge states C and E. Thus, we can observe the eigenmode switching of the two edge states across the second BG, which characterizes the topological band inversion. Furthermore, if two states at the lower or upper edges of the overlapped BG possess different symmetries, the sign ς of this BG is opposite and thus an interface state exists inside the BG (Xiao et al., 2014). Here, the lower edge states B and E in Figs. 5(b) and (e) for S1-and S2-configurations are symmetric and antisymmetric, respectively. This indicates the different signs ς of the second BG and various topological properties, as shown in Figs. 4(d) and (f).

Transmission spectra and displacement distributions
Here, we make use of Eq. (42) to calculate the transmission spectra (i.e., the transmission coefficient t N versus the frequency f ) of a finite-size neo-Hookean PCC. Figs. 6(a) and (b) show the results for the finite PCC composed of 10 identical S2-type unit cells subjected to the axial forces F N = 0 and 1, respectively. The transmission spectra based on our theoretical analysis agree well with the corresponding band structures (compare Figs. 6(a) and (b) with Figs. 4(c) and (f)). The transmission coefficient in the BG range approaches zero, indicating that the wave propagation is prohibited; but in the passing bands, the elastic wave is allowed to propagate in the finite PCC as the transmission ratio almost equals one. The spectrum does not show any peaks in the BGs for the finite-size PCC with identical unit cells. In addition, the transmission spectra have been calculated independently by means of the finite element (FE) simulations. The details about the FE simulation procedures are described in Appendix B. Figs. 6(c) and (d) shows the FE results for the finite-size PCC under the action of axial forces F N = 0 and 1, respectively. Here, the attenuation intensity T (dB) is defined as T = 20 log(A output /A input ), where A input and A output are the average displacement amplitudes of the input and output signals. The comparison of the theoretical and numerical results demonstrates a good agreement between these intendant methods. There are, however, slight differences in the form of transmission spectra due to various calculation expressions. The FE calculations show that the transmission coefficient reaches a dip within the BG range, which stands for a strong attenuation of the wave propagation; however, the attenuation intensity is larger than zero in the passing bands, implying that the output signal can be obviously detected. Particularly, the frequency ranges of the second BG for F N = 0 are (316 Hz, 383 Hz) and (307 Hz, 385 Hz) for the theoretical results and FE simulations, respectively. For F N = 1, the corresponding results based on the theoretical model and FE simulation are (248 Hz, 302 Hz) and (240 Hz, 299 Hz), respectively. We note, however, that the difference in the central frequencies between the theoretical and FE predictions does not exceed 2%.
Recall that the topological phase transition exists and the topological properties of the second BG are different for S1-and S2-configurations (see Sec. 5.2.1). Thus, a topological interface state should appear in this overlapped BG at the interface delimiting the S1-and S2-configurations. Motivated by this prediction, we design a mixed finite-size neo-Hookean waveguide consisting of 5 S1-type and 5 S2-type unit cells. We utilize the theoretical formula (51) to calculate the corresponding transmission spectra of the mixed-type PCC waveguide. Fig. 7 shows the theoretical ((a) and (b)) and FE simulation ((c) and (d)) results for the mixed soft PCC waveguide. The results are illustrated for the unloaded PCC (Figs. 7(a) and (c)), and for the PCC subjected to the axial force F N = 1 (Figs. 7(b) and (d)). We observe from Fig. 7 that the sharp transmission peaks emerge in the overlapped BG frequency range for both unloaded and loaded cases. In the unloaded mixed PCC waveguide, the peak frequencies of theoretical prediction and FE simulation are 355 Hz and 352 Hz, respectively. In the loaded case with F N = 1, the corresponding peak frequencies shift down to 276.5 Hz (theory) and 275 Hz (simulation). Thus, the application of axial force tunes the position of the transmission peak. We note that the peak frequencies predicted by the FE simulations and theoretical model are almost identical, which further validates the effectiveness of our theoretical assumption. Figure 8 depicts the spatial distributions of the displacement modes for the mixed finitesize neo-Hookean waveguide at the transmission peak frequencies corresponding to the unloaded and loaded (F N = 1) states. Figs. 8(a) and (b) show theoretical results for the absolute value |w(z)| of displacement distributions at the peak frequencies 355 Hz (F N = 0) and 276.5 Hz (F N = 1). The normalized axial coordinate is defined as z * = z/l, where l is the length of the deformed unit cell of S1-configuration. Through the FE simulations, Figs. 8(c) and (d) display the displacement mode shapes at the FE peak frequencies 352 Hz (F N = 0) and 275 Hz (F N = 1). We observe that the displacement field is localized at the interface between the two PCC elements, and it decays dramatically towards the ends of the mixed waveguide. This is an evident sign of the interface state as a result of the topological conflict of the distinct states. In particular, the displacement amplitude at the interface is 6 times more than the input signal (see Figs. 8(a) and (b)).
It is worth noting that the topological interface state is different from the concept of resonant mode (Li et al., 2017;Chen et al., 2019b). The resonant mode is greatly affected by the boundary conditions, local resonance and excitation location, whereas the topological interface state method -based on the topological property conflict -provides a robust mechanism against the wave propagation direction or boundary conditions (Yin et al., 2018). To validate that the occurrence of transmission peaks is ascribed to the topological interface state, we have reversed the input and output ends, and then calculated the transmission spectra as well as the displacement distributions for the unloaded and loaded (F N = 1) cases. The results (omitted here) indicate that the topological interface state is still observed at the same transmission peak frequencies as those in Fig. 7.

Topological phase diagram
Next, to analyze the effect of axial force on topological interface states, we examine the topological phase diagram. Fig. 9 shows the two edge-state frequencies of the second BG as functions of the geometric parameter φ 0 of the neo-Hookean PCC for different levels of the applied axial force. The variations of the second BG frequencies with φ 0 are calculated by setting q = 0 in Eq. (36). Fig. 9 illustrates that an increase in the axial force results in a lower frequency of the topological transition point (or band crossing point). However, there is only a slight change in φ 0 where the band crossing happens. In particular, for F N = 0, 0.5, 1 and 1.5, four topological transition points occur at φ 0 = 0.5, 0.514, 0.524 and 0.527, with their topological transition frequencies being 355 Hz, 310.6 Hz, 275.9 Hz and 251.5 Hz, respectively. Note that for the neo-Hookean PCC subjected to a fixed positive tensile loading, the BG central frequency increases with an increase in φ 0 . For example, for PCC under the action of F N = 1, the BG central frequency increases from 272.5 Hz at φ 0 = 0.35 to 280.4 Hz at φ 0 = 0.65. This is due to the fact that for a smaller φ 0 , the thinner sub-cylinder 2 occupies more in the unit cell. In the neo-Hookean PCC with smaller φ 0 the unit cell elongates more, thus increasing the overall geometric size of the structure (recall Fig. 3(b)). Moreover, according to Ding et al. (2006), the natural frequency of cylindrical structure decreases with an increase in the slenderness or length-to-radius ratio. Therefore, the BG central frequency shifts down towards a lower frequency for PCC with a smaller φ 0 . We note, however, that this trend is reversed if a compressive axial force is applied (F N < 0); namely, the second BG central frequency decreases with an increase in φ 0 . The four groups of solid curves represent the two edge-state frequencies of the second BG as functions of the initial length fraction φ0 for four different axial forces. The topological phase curve (dash-dotted line) connects the topological transition points corresponding to different axial forces, and the arrow denotes the direction in which the axial force increases. The yellow and green filled regions indicate the BG signs with ς > 0 and ς < 0 respectively, which are labelled at two sides of the topological phase curve.
Following Xiao et al. (2014), the topological transition point can be obtained analytically in Appendix C. The dash-dotted line shown in Fig. 9 is referred to as the topological phase curve for different axial force levels ranging from −0.3 to 5, which is determined from the formulae in Appendix C. We see from Fig. 9 that all the topological transition points formed by the close of the second BG are connected by the topological phase curve, which provides the frequencies of transmission peaks in Figs. 7(a) and (b) for axial forces F N = 0 and 1. Moreover, the topological phase curve divides the topological phase diagram into two regions with different topological properties, the BG signs of which are indicated in Fig. 9. According to the topological phase diagram, we can design conveniently tunable topological interface states in a soft PCC. For example, we can construct a topological waveguide composed of two types of PCCs with φ 0 = 0.35 and 0.65. As the axial force F N increases from 0 to 2, the frequency of topological interface state is tuned from 355 Hz to 236.3 Hz. It should be emphasized that if the BG has no common frequency range for the chosen geometric size and axial force, there will be no topological interface state.
Thus, the band structure and topological interface state can be actively tuned towards a lower frequency by applying the axial force in the neo-Hookean waveguide. Its PCC elements should be properly selected to have different topological properties, so that the existence of low-frequency tunable topological interface states is guaranteed.

Influence of the strain-stiffening effect
Here, we consider a PCC made out of nonlinear material with a strong strain-stiffening behavior. In particular, we employ the Gent material model and analyze the stiffening effect on the tunable topological interface states. When the axial force is not large enough (e.g., F N ≤ 1.25), the band structure and transmission behavior of the Gent model resemble those of the neo-Hookean model (recall Figs. 4 and 6) since the hyperelastic material has not reached the stiffening stage (see Fig. 3) for the given locking parameter J m = 20. Therefore, in Fig. 10, we display the topological transition of band structures and the transmission spectrum for a large axial force F N = 2. In Figs. 10(a) and (c), the calculated Zak phase 0 or π is indicated on the related bulk band and the second BGs are marked by yellow (ς > 0) and green (ς < 0) stripes with different topological properties. Similar to the previous observations in the neo-Hookean PCC (recall Fig. 4), the band inversion can be obtained for the Gent PCC. Here, the second BG closes as φ 0 increases from 0.35 to 0.494, and it reopens with a further increase in φ 0 . Other illustrations are not repeated for brevity. Nevertheless, in view of the triggered strain-stiffening effect, the position of the second BG for the Gent PCC is higher than that of the neo-Hookean PCC for F N > 1.25. Particularly, for the S1configuration PCC subjected to the axial force F N = 2, the frequency limits of the second BG vary from (225 Hz, 254 Hz) for the neo-Hookean model to (256 Hz,286 Hz) for the Gent model. Next, we examine the wave characteristics of a mixed Gent PCC waveguide subjected to the axial force F N = 2. Fig. 11(a) demonstrates the transmission spectrum of the waveguide consisting of 5 S1-type unit cells (φ 0 = 0.35) and 5 S2-type unit cells (φ 0 = 0.65). The transmission spectrum is calculated with the help of Eq. (51). The transmission peak emerges in the second overlapped BG, with the peak frequency 271.7 Hz for this case (see Fig. 11(a)). The corresponding displacement field distribution at 271.7 Hz is shown in Fig. 11(b). It is confirmed that the displacement is mainly confined in vicinity of the interface (with the amplitude nearly 6 times over the input signal) and attenuates rapidly towards the ends of the hyperelastic waveguide. To further explore the influence of strain-stiffening effect on the topological interface state, we plot the corresponding topological phase diagram in Fig. 12. The results are depicted for the Gent PCC (J m = 20) with four topological transition variations at F N = 0, 1, 2 and 4. The topological phase curve calculated according to the theoretical formulae (see Appendix C) is also included in Fig. 12. The topological phase curve and the band inversion curves for the Gent PCC are almost the same as those for the corresponding neo-Hookean PCC (compare with Fig. 9) for the range of loadings not reaching the stiffening stage (i.e., when F N ≤ 1.25). In particular, the frequency of the topological transition point continuously decreases with an increase in axial force level in the range of F N ≤ 1.25. However, a further increase in the axial force leads to different variation trend in the topological transition point frequency: after decreasing to the lowest value 270.8 Hz at around F N = 1.8, the frequency starts to increase conversely and rapidly (see Fig. 12). This is a unique feature of the nonlinear PCC with a strong strain-stiffening effect.
Remarkably, the tensile force -if large enough -can cause a reverse trend in the second BG central frequency versus φ 0 in Gent PCC. For example, in the Gent PCC subjected to F N = 2, the BG central frequency decreases with an increase in φ 0 . This reverse trend is even more prominent for the Gent PCC subjected to F N = 4 (see the corresponding curves in Fig. 12). Note that this BG central frequency trend reversion is not observed in the neo-Hookean PCC (see Fig. 9). This is again an important manifestation of the strain-stiffening effect that begins to prevail over the deformation-induced geometric change upon achieving a certain level of the applied tensile force.
Interestingly, we find from Fig. 12 that the length fraction φ 0 corresponding to the topological transition point has a more obvious variation with the increasing axial force applied to the Gent PCC, compared with the neo-Hookean case (see Fig. 9). This phenomenon combined with the frequency non-monotonous change results in a peculiar state for morphologies enclosed by the topological phase curve. In this morphology domain (denoted by the violet region in Fig. 12), the BG sign can switch depending on the applied axial force. For example, for F N = 1, this special area is on the left of the topological phase curve with ς > 0, while for F N = 4, the violet area turns to be on the right of the topological phase curve with ς < 0. Therefore, for F N ≤ 1.8, ς > 0; but ς < 0 for F N > 1.8.
In addition, when the axial force exceeds F N = 4, the BG frequency ranges for PCC elements (whose unit-cell length fractions φ 0 take values from different sides of the topological . The four groups of solid curves represent the two edge-state frequencies of the second BG as functions of φ0 for four different axial forces. The topological transition points are located on the topological phase curve (dash-dotted curve), and the arrow denotes the direction in which the axial force increases. The yellow and green filled regions indicate the BG signs with ς > 0 and ς < 0 respectively, which are labelled at two sides of the topological phase curve. The violet region marks the morphologies for which the topological property switches depending on the axial force.
transition point) do not overlap (see Fig. 12), and hence the topological interface state does not exist for any combination of the geometric parameters. For example, for Gent PCC subjected to the axial force F N = 4.1, the topological transition point is φ 0 = 0.444 and the corresponding frequency is 333.8 Hz. In vicinity of this topological transition point, the BG frequency limits for φ 0 = 0.44 are (333.8 Hz, 334.4 Hz), while those for φ 0 = 0.45 become (332.9 Hz, 333.8 Hz); thus, there is no overlapped BG frequency.
Furthermore, for a mixed finite Gent PCC waveguide composed of 5 S1-type unit cells (φ 0 = 0.35) and 5 S2-type unit cells (φ 0 = 0.65), the topological interface state exists until the axial force reaches the level F N ≈ 3.8. When the applied axial force exceeds this level, the frequency limits of the second BG have no overlapped part, and although the topological properties of the two PCC elements are different, the topological interface state cannot be activated.

Conclusions
We studied a class of 1D soft PCs possessing the topologically protected interface states. The large-deformation ability of soft waveguides combined with the material stiffening effect is exploited to tune the topologically protected states. In particular, we illustrate this concept based on the example of the 1D waveguide composed of two types of soft PCCs with different topological characteristics. Here, the topological interface state for longitudinal waves is tunable by application of an external axial force. First, we utilized the nonlinear elasticity theory combined with the assumption of uniform deformations to determine the nonlinear static response of PCCs under the action of an axial force. Next, the dispersion relation for small-amplitude longitudinal waves, the unit-cell mode shape as well as the displacement field distribution and signal transmission coefficient of the finite-size cylindrical waveguide were derived analytically. Finally, the theoretical predictions along with the numerical calculations were analyzed to elucidate how the external loading and the material stiffening affect the frequency tunability and the existence of topological interface states. Our main observations are summarized below: (1) The BG inversion process (i.e., the BG open, close and reopen process) accompanied by the topological phase transition can be realized by altering the initial PCC geometric parameter and be tuned by adjusting the axial force.
(2) For the neo-Hookean waveguides, the axial tensile force lowers monotonically the frequency of topological interface states owing to the generated elongation in the whole system.
(3) For the Gent waveguides, the frequency of topological interface states varies with the axial force in a non-monotonous way; this is a result of the competition between the geometric change and the material strain-stiffening effect.
(4) In reference to the topological phase diagram, the tunable position and existence condition of topological interface states are clearly demonstrated when changing the axial force in a properly pre-designed system.
Our results -based on the example of 1D soft PCC waveguides -indicate the possibility to realize on-demand tunability of the topological interface states. The present study provides guidelines for further design of actively tunable topological wave devices operating at the low-frequency range. These systems may find a wide range of potential applications such as tunable energy harvesters, low-pass filters and high-sensitivity biomedical detectors.
It should be emphasized that the tunable topological interface or edge states in 2D metamaterial systems could be achieved by means of the electromechanical biasing fields Zhou et al., 2020a), which is an interesting topic to be addressed for further applications in the future.
By discretizing q and noting the relation ln (x − 1 + 1) → x − 1 in the limit of x → 1, Eq. (A.1) with ∆q → 0 can be rewritten in a discretized form as which, after some simplifications, yields Consequently, after obtaining the Bloch displacement distribution w j,q of the deformed unit cell in Subsec 4.3, we can exploit the relation W j,q = w j,q e −iqz and Eq. (A.3) to calculate the Zak phase numerically.

Appendix B. FE simulations
To validate our theoretical model, numerical simulations are conducted by using Abaqus, an FE analysis and solver software. The simulations to calculate the transmission spectra are implemented for the finite PCC structure, including its geometric parameters and material properties of silicon rubber (Zhermarck Elite Double 32) . Here, we adopt the neo-Hookean hyperelastic model and establish an axisymmetric structural model with a fine mesh of 8-node hybrid elements (i.e., element type CAX8H in Abaqus).
In order to understand the longitudinal wave propagation behaviors in a pre-deformed structure, the static analysis (i.e., pre-stretching the structure) and frequency domain analysis (i.e., wave propagation in the structure) are performed consecutively in Abaqus: Step 1 (Static analysis): The boundary conditions in accordance with the theoretical model are applied to both sides of the structure to simulate its extension procedure. The resultant displacement and stress fields in the deformed structure are recorded and saved.
Step 2 (Frequency domain analysis): The displacement and stress fields calculated in the static analysis are imported to the structural model as the initial deformed state. For calculating the transmission spectra, a sinusoidal axial-displacement excitation over a frequency range of interest is imposed to one input side of the finite-size PCC structure and its average displacement amplitude is calculated as the input signal A input . Additionally, the average displacement amplitude at output side is collected as the output signal A output . Thus, the attenuation intensity T (dB) is defined as T (dB) = 20 log (A output /A input ), which describes the elastic wave transmission behavior.

Appendix C. Conditions for the band crossing
Following the argument of Xiao et al. (2014) for 1D photonic crystals, this appendix will provide the conditions for two bands to cross for a soft PCC subjected to an axial force.
In this work, we focus on the second BG. In order to make the second BG close, we have γ = m 1 = m 2 = 1 and hence α = 1 is a rational number. We can take two steps to obtain the position of the band crossing point (or topological transition point). First, for a given axial force F N , the mechanical biasing field state and the parameters c (p) are determined. Second, using the relation α = 1 and the condition that the length of the undeformed unit cell is kept fixed (i.e., L = L (1) + L (2) ), we can obtain the corresponding geometric sizes L (1) and L (2) (or initial length fraction φ 0 = L (1) /L) and then calculate the frequency ω c of topological transition points from Eq. (C.3).