Analytical modeling of one-dimensional resonant asymmetric and reciprocal acoustic structures as Willis materials

As building blocks of acoustic metamaterials, resonant scatterers have demonstrated their ability to modulate the effective fluid parameters, which subsequently possess extreme properties such as negative bulk modulus or negative mass density. Promising applications have been shown such as extraordinary absorption, focusing, and abnormal refraction for instance. However, acoustic waves can be further controlled in Willis materials by harnessing the coupling parameters. In this work, we derive the closed forms of the effective parameters from the transfer matrix in three asymmetric and reciprocal one-dimensional resonant configurations and exhibit the differences in terms of coupling coefficients. The way in which Willis coupling occurs in spatially asymmetric unit cells is highlighted. In addition, the analysis shows the absence of odd Willis coupling for reciprocal configurations. These effective parameters are validated against experimental and numerical results in the three configurations. This article paves the way of a novel physical understanding and engineering use of Willis acoustic materials.


Introduction
Since the seminal work of Willis in the 80's [1], the eponymous materials have received an increasing attention, because of their analogy with bi-isotropic electromagnetic metamaterials [2]. The Willis coupling parameters couple the potential and kinetic energy in the acoustic conservation relations, therefore enhancing the ability to control waves in metamaterials compared to other materials that do not exhibit such coupling. These parameters have thus been employed to design PT symmetric [3], wave front shaping [4,5], or non-reciprocal [6] systems. Willis coupling arises from chiral inhomogeneities [7], asymmetric inhomogeneities, nonlocal effects, and nonreciprocal biases [8]. Although most of the works to date have focused on the experimental evidence [4,9,10], physical origins [11], calculation [12][13][14], and enhancement [15] of Willis coupling, only a few have focused on deriving a closed form of these parameters. The present article aims at filling this gap and therefore at easing its physical understanding and engineering use. Effectively, it turns out that various systems rely on asymmetric meta-atoms for perfect absorption in transmission problem [16], non-hermicity of the acoustic waves [17], PT symmetry [18], or more generally most of the double negative one-dimensional devices [11,19,20].
We focus on three resonant, asymmetric, and reciprocal one-dimensional unit cells and derive closed forms of the corresponding Willis coupling parameters, exhibiting different forms depending on the nature of the asymmetry. By analyzing these Willis parameters, we show that Willis coupling in resonant systems arises from second order Taylor expansion originating differently from multilayer systems [13,21,22] and exhibit the dipolar feature of the coupling via arm terms. In addition, the reciprocity condition directly provides even coupling and vanishing odd coupling [11], thus inducing doubt on the existence of odd coupling in one-dimensional reciprocal acoustic systems. While even Willis coupling parameters appear with opposite sign in the propagation matrix, odd Willis couplings appear with identical sign in the propagation matrix for reciprocal structures.
The article is organized as follows. The general procedure for the derivation of the effective parameters is detailed in section 2. It relies on the Padé's approximation of the total transfer matrix, which links the state vectors at both sides of a unit cell. The procedure is applied in section 3 to three different resonant asymmetric and reciprocal one-dimensional unit cells. While the first two unit cells are composed of detuned resonators either in parallel or in series of a duct, the third one combines resonators in parallel and in series in the duct, see figure 1. Two of these unit cells have already been studied as Willis materials [9,11], but the closed forms of the coupling parameters have not been provided. The derived forms clearly unveil dipolar feature of Willis coupling, as well as the differences between resonant and non-resonant asymmetries and nature of the resonant asymmetry. The effective parameters are validated against experimental and numerical (derived from [13]) results for each unit cell in section 4. The order of Taylor expansion required to correctly model an asymmetric resonant unit cell is questioned. Finally, concluding remarks and perspectives are provided. Additional information are also given in appendices.

Derivation of the Willis coupling parameters
We consider a one-dimensional asymmetric and reciprocal material composed of the periodic repetition of a unit cell of length d. The pressure and particle velocity (alternatively the flow in a duct) form a state vector W T =< p, V >, where T is the transpose operator. This state vector satisfies the following matrix equation which directly arises from the mass conservation and constitutive law: Assuming an implicit time dependence e −iωt , the matrix A reads as for Willis type materials, where ρ, K and χ are respectively the density, bulk modulus and Willis coupling parameter. Please note that a usual isotropic and symmetric effective fluid implies χ = 0. The state vectors at both sides of the unit cell, W(d) and W(0), are thus linked by where the term exp (Ad) is the matrix exponential of Ad also known as the transfer matrix T. Among the different ways of evaluating the matrix exponential [23] is the Padé's approximation. The transfer matrix T is thus approximated in the long wavelength regime, i.e. when the wavelength λ is much larger than the period d, by where I is the identity matrix. Assuming the transfer matrix T between the state vectors at both sides of the unit cell being known, equation (4) can be inverted and the expression of the constitutive matrix A becomes The Padé's approximation is of particular interest compared to other approximations such as the Taylor's expansion (see appendix A), because it allows us to account for the reciprocity of the material [24]. This material property imposes det (T) = 1, i.e. t 11 t 22 − t 12 t 21 = 1, where t ij are the elements of the matrix T. Accounting for the reciprocity, equation (5) can be rewritten in the form from which we immediately see that the diagonal terms are of opposite signs and we can identify , by comparison with equation (2). Density, bulk modulus and Willis coupling parameter can then be directly calculated via the elements of the transfer matrix that links the state vectors at both sides of the unit cell. Please note that the odd Willis coupling introduced in reference [11] is completely canceled from equation (6) only because of the reciprocal condition and that only even coupling remains. The effective parameters that are derived following this procedure are indicated by a subscript e in the following.

Explicit effective parameters of Willis materials composed of resonant elements
The total transfer matrix T is evaluated thanks to the knowledge of different elementary transfer matrices that are presented in appendix B. A circular duct of radius r and length d much smaller than the wavelength is considered. Only plane waves propagate in this duct, i.e. the frequencies are lower than that of the first cut-off of the duct and possible evanescent coupling between the unit cell elements is neglected. The reduced density and bulk modulus (or alternatively the wavenumber and reduced impedance) in the duct areρ andK (or alternativelyk = k andZ). Three asymmetric and reciprocal unit cells are considered and are represented in figure 1: (i) when the duct is loaded by two detuned Helmholtz resonators (HR) located at different positions (figure 1(a)) of impedancesZ (j) HR , j = 1, 2, (ii) when two detuned plates are clamped in the duct (figure 1(b)) of impedancesZ (j) p , j = 1, 2, and (iii) when the duct is loaded by an HR and a plate is clamped in it at a different position (figure 1(c)). Each resonant element possesses its own local dynamics which is assumed different from that of the duct, i.e. co-dynamic regime is assumed [25]. This is an important difference with respect to a laminated two-material unit cell for example, for which it is clear that a second order Taylor expansion of the total matrix elements in equation (7) is required to exhibit Willis coupling terms as demonstrated in appendix C. The order of the resonant element impedances becomes unclear notably around the resonance as explained in appendix B. Nevertheless, 1/Z HR (Lorentzian function) andZ p (inverse of a Lorentzian function) are assumed to vary like O (kd) 2 to ensure reciprocity of the configuration in the considered frequency range, see for example the discussion of appendix D in the presence of a single HR. Note that this relies more on a frequency analysis rather than on a purely kd analysis.

Unit cell composed of a duct loaded by two detuned Helmholtz resonators
We first consider a unit cell of length d composed of a straight circular duct loaded by two detuned HR of respective reduced impedancesZ (1) HR andZ (2) HR , assumed to be point-like resonators, and located at l (1) and l (1) + l (2) such that d = l (1) + l (2) + l (3) , as represented figure 1(a). Both inverse impedances 1 Z (1) HR and 1 Z (2) HR are assumed to be O (kd) 2 terms. This configuration is formally that considered in reference [3], where the loading quarter-wavelength resonators are replaced by HR. The total transfer matrix reads as where the transfer matrices accounting for the propagation along each length and for the loading HR are given in equations (B.1) and (B.2). Assuming kd λ, in such a way that kl (1) = ζ (1) kd, kl (2) = ζ (2) kd, and kl (3) = ζ (3) kd, with ζ (j) = l (j) /d, j = 1, 2, 3, are also much smaller than the wavelength, and making use of equation (7) leads to 3 , To get a grip on these equations, we compare them to those in the presence of a single HR presented in appendix D. The Willis coupling term appears as the sum of the two Willis coupling terms arising from each HR. Each of them clearly exhibits a momentum introduced by the resonator as testified by the moment arm term l (3) + l (2) − l (1) for the first HR and the moment arm term l (3) − l (1) − l (2) for the second HR. The Willis coupling parameter vanishes when the symmetry is introduced, for example, when the two HR are identical and when l (3) = l (1) . It appears in the form of an effective density divided by the impedance of the HR. The density term, second line of equation (9), is not the subject of any specific remark and simply reads as the effective density in the duct. The bulk modulus is the sum of the contribution of the effective properties of the duct in the absence of the loading resonators, i.e. 1/K, and of each HR, i.e. −i/dωZ (j) HR , j = 1, 2. The presence of the HR causes the effective bulk modulus to become negative for frequencies around their resonances [26,27]. All in all, the effective parameters appear as the sum of those of each elements, i.e. those of each segments and those of the HR, without any particular coupling. When a subperiodicity can be introduced, i.e. when the two HR are identical and when l 2 = d/2, the effective parameters reduce to those of a unit cell of length d/2 consisting in a single, possibly uncentered, HR with Willis parameter exhibiting a moment arm term equal to l (3) − l (1) .

Unit cell composed of a duct with detuned clamped plates
We now consider a unit cell of length d, as depicted in figure 1(b), composed of a straight circular duct in which two detuned plates, assumed to be point-like resonators, of respective reduced impedancesZ (1) p and Z (2) p are clamped at l (1) and l (1) + l (2) such that d = l (1) + l (2) + l (3) . The clamped plates (CP) simply replace the HR when compared to the previous configuration, but this time the impedances are in series instead of being in parallel. This configuration is formally that considered in reference [9]. The total transfer matrix reads as where the transfer matrices accounting for the propagation along each length and for the CP are given in equations (B.1) and (B.4). Proceeding similarly as in the previous section, equation (7) together with the second order Taylor expansion of the total transfer matrix elements provides 3 , The Willis coupling term appears as the sum of the Willis terms associated with each plate (see appendix D). The momentum seems opposite to that imposed by the HR with an arm term l (1) − l (2) − l (3) for the first CP and an arm term l (1) + l (2) − l (3) for the second CP. It appears in the form of an effective compressibility multiplied by the impedance of the plate. Again, this term vanishes when the symmetry is introduced as for example, when the two CP are identical and l (3) = l (1) . The density is the sum of the contribution of the effective density of the duct in the absence of the CP, i.e.ρ, and of each CP, i.e. −iZ p (j) /dω, j = 1, 2. The presence of the CP causes the density to be negative for frequencies lower than their resonances [28]. Finally, the effective bulk modulus is that in the absence of the clamped plate and is thus not the subject of any specific remark. The effective parameters again appear as the sum of those of each elements, i.e. those of each segments and those of the CP, without any particular coupling. When a plane of symmetry can be introduced, i.e. when the two CP are identical and when l 2 = d/2, the effective parameters reduce to those of a unit cell of length d/2 consisting in a single, possibly uncentered, CP with Willis parameter exhibiting a moment arm term equal to l (1) − l (3) .

Unit cell composed of a duct with a clamped plate and loaded by a detuned Helmholtz resonator
We finally consider a unit cell of length d composed of a straight circular duct in which a plate, assumed to be a point-like resonator, of reduced impedanceZ p is clamped at l (1) . The duct is in addition loaded by an HR, also assumed to be a point-like resonator, of reduced impedanceZ HR and located at l (1) + l (2) such that d = l (1) + l (2) + l (3) . A picture of the configuration is shown in figure 1(c). In other words, the first HR is replaced by a CP compared to the configuration studied in section 3.1 or the second CP is replaced by an HR compared to the configuration studied in section 3.2. Beyond the expected double negative effective property [29], this asymmetric configuration combines impedances in series and in parallel and is similar to that studied in reference [11]. The total transfer matrix reads as where the transfer matrices accounting for the propagation along each length, for the CP, and for the HR are given in equations (B.1), (B.4), and (B.2) respectively. Again, equation (7) together with the second order Taylor expansion of the total transfer matrix elements provides Contrary to the configurations studied in both previous sections, the Willis coupling term is not only the sum of the Willis terms associated with the presence of the CP and of the HR (see appendix D), but also exhibits a coupling between the CP and the HR thanks to the term iZ p /2 dωZ HR . Note that the moment arms are introduced by the first two terms, that related to the presence of the CP and that related to the presence HR, while the coupling term does not present moment arm term. In addition, the Willis coupling parameter never seems to vanish, because the configuration is structurally asymmetric. The density is the sum of the contribution of the effective density of the duct in the absence of the CP and of the HR, and of the CP, i.e. −iZ p /dω. In a similar way, the bulk modulus is the sum of the contribution of the effective bulk modulus of the duct in the absence of the CP and of the HR, and of the HR, i.e. −i/dωZ HR . Note that obviously no symmetry plane can be introduced for this configuration.

Experimental validation of the effective parameters and discussion
All experiments are conducted in a duct of radius r = 2.5 cm, see figures 1(a)-(c). The experimental set-up consists of a 4 microphone measurement system with a pair of microphones upstream and a second pair of microphones downstream of the sample. The microphones that compose each pair are separated by a distance of 2.5 cm. A step signal from 100 Hz to 1000 Hz with a step of 1 Hz is delivered by a loudspeaker at one end of the duct and an anechoic termination is mounted at the opposite end. Temperature, humidity and atmospheric pressure are recorded for each experiment. The transfer function between the loudspeaker and each microphone is recorded by an NI USB-4431 acquisition card driven by the INTAC software. Each sample is measured in both direct and reverse orientations in order to form an overdetermined system based on the scattering matrix [3,30] as explained in appendix E to solve for T, R + , and R − , i.e. the transmission, the direct orientation reflection, and the reverse orientation reflection coefficients. These coefficients are also calculated by the transfer matrix method (TMM) using the total transfer matrix relying on the elementary matrices and the Pade's approximation of transfer matrix relying on the evaluated effective parameters I − A e d/2 −1 I + A e d/2 , equation (4), whereZ is the reduced impedance of the surrounding medium. From the measured direct and inverse orientation reflection and transmission coefficients, the effective parameters are reconstructed following the procedure described in appendix F. The effective properties can also be directly evaluated from the total transfer matrix as explained in appendix G, which appears as a numerical version of the procedure described in reference [13]. Figures 2(a), (e) and (i) depict respectively the absolute values of the two reflection and transmission coefficients, calculated by the TMM with the total transfer matrix and by the Pade's approximation using the effective properties derived in section 3 for the three configurations considered in the present article. A single unit cell is measured for each configuration. This measurement is entirely sufficient, because one-dimensional structures are assumed. The first configuration (see figure 1(a)) consists in two detuned HR, the dimensions of the cavities and necks of which are l (1) c = 8 cm, l (2) c = 4.2 cm and r (1) c = r (2) c = 2.15 cm, and l (1) n = l (2) n = 2 cm and r (2) n = r (1) n = 3 mm, separated by a distance l (2) = 5 cm. The resonant frequencies of both HR are thus f (1) HR ≈ 165 Hz and f (2) HR ≈ 230 Hz. The two remaining dimensions l (1) and l (3) are chosen identical, i.e. l (1) = l (3) = 1 cm, such that d = 7 cm. The second configuration (see figure 1(b)) consists in two CP separated by a distance l (2) = 5 cm + h (1) p + h (2) p /2. The first CP is a plastic shim plate of thickness h (1) p = 254 μm and material properties ρ (1) p = 1400 kg m −3 , ν (1) p = 0.41, and E (1) p = 4.6 (1 − 0.03i) GPa already used in reference [31]. The second CP is a poroelastic plate of thickness h (2) p = 3.1 mm and material properties ρ (2) p = 40 + 380i/ √ ω kg m −1 , ν (2) p = 0.1, and E (2) p = 470 − 0.007iω kPa already used in reference [32]. The resonant frequencies of both CP are thus f (1) p ≈ 380 Hz and f (2) p ≈ 255 Hz and the two remaining dimensions l (1) and l (3) are chosen almost identical, i.e. l (1) = 5 mm + h (1) p /2 and l (3) = 5 mm + h (2) p /2, such that d = 6 cm + h (1) p + h (2) p . The last configuration (see figure 1(c)) is composed of the poroelastic plate and a loading HR, the cavity and neck dimensions of which are l c = 5 cm and r c = 2.15 cm, and l n = 2 cm and r n = 2 mm, separated by a distance l (2) = 5 cm + h p /2. The resonance frequency of the CP is still f p = 255 Hz and that of the HR is f HR = 140 Hz and the two remaining dimensions l (1) and l (3) are again chosen almost identical, i.e. l (1) = 1 cm + h p /2 and l (3) = 1 cm, such that d = 7 cm + h p . The lengths of the unit cells are thus much smaller than the wavelengths over the frequency range considered in each configuration. The neck was manufactured by fused deposition modeling to fit in an initial radius of 1 cm. The neck also presents some rugosity inherent to this rapid manufacturing technique, which can influence the viscothermal losses in these resonator narrow elements [33]. The configurations were designed to present a reflection coefficient that vanishes at a specific frequency while the reflection coefficient of the reverse sample is different from zero at this frequency. The calculated coefficients are found in very good agreement with the experiments, which therefore valid the derived effective properties and prove the asymmetry of the acoustic response of each configuration. Some discrepancies are visible, notably for the third configuration at very low frequencies. These discrepancies are attributed to measurement complexity at low frequencies and to possible evanescent coupling between the resonant elements that is not accounted for in our modeling. At high frequency, the calculated transmission coefficient based on the effective parameters starts to deviate from its value for the second configuration because the scale separation is not ensured anymore and the Pade's approximation is not valid anymore, see figure 2(e). In addition, slight shifts between the coefficients calculated with the total transfer matrix and with the effective parameters can also be notice at low frequencies and around the resonances, see figures 2(e) and (j). These discrepancies are attributed to the assumption we made on the dependency of the resonant element impedances, that are supposed to vary like kd. Around the resonances, this assumption is not fully valid and the effective properties would require additional terms to be Taylor expanded to the second order.
The normalized effective properties, i.e. χ e S, ρ e /ρ, and K e S/γP 0 , are respectively depicted in figures 2((b), (f), (j)), ((c), (g), (k)) and ((d), (h), (l)) for each configuration. The numerical and experimental effective properties match those evaluated in section 3. Generally speaking, the real parts of the effective densities (CP) and bulk moduli (HR) are negative in stop bands and follow regular trends [26][27][28][29]. Some discrepancies are visible for the third configuration again due to the fact that the phenomenon are encountered at very low frequencies, but also to the assumption made on the order of the resonant element impedances. The amplitude of the Willis parameter is usually much smaller than the other ones, partially due to the fact that this parameter is only normalized by the cross-sectional area of the main duct S, in the absence of a normalization value for this parameter. Note that χ e S vanishes away from the resonances in case of detuned HR, while it tends to infinity at low frequency when a plate is involved, see figures 2(f) and (j). A similar remark is made concerning the effective density ρ e /ρ. This is due to the presence of the termZ p in both χ e and ρ e , equations (11) and (13), which varies like 1/ω at low frequencies. This also blurs the required order of Taylor expansion to derive these effective parameters even away from the resonances in case of unit cells involving plates as it can be seen in figures 2(e) and (i). This divergence is nevertheless physical, because a plate should be a rigid wall at low frequency. Note that to reach the request accuracy near to resonance, non-local effective medium are probably requested, which are outside the scope of the current paper. The dispersion relation is also plotted for the configuration comprising a straight duct with a CP and loaded by a detuned HR in figures 2(m) and (n). The grey areas highlight the frequency range where the real parts of both the effective density and bulk modulus are negative. The viscothermal losses and Willis coupling induce a slight shift in the negative index part. The dispersion relation for the other two configurations (not shown here) are found in good agreement.

Conclusion
We derive closed forms of the effective properties of Willis materials thanks to a Pade's approximation of the total transfer matrix linking the state vectors at both sides of an asymmetric and reciprocal one-dimensional unit cell. We primarily show that the reciprocal condition leads to the unique existence of the even Willis coupling. Similarly to the case of laminated structures, which are usually in nonresonant acoustic structures, second order Taylor expansion of the transfer matrix elements is sufficient to exhibit Willis coupling parameters. Nevertheless, this result relies on a strong assumption on the resonant element behavior and its veracity is questioned around the resonance. Higher order Taylor expansion of the full transfer matrix terms involving the resonant elements may therefore be more suitable. The dipolar feature of the Willis coupling is clearly evidenced, because the coupling parameter presents moment arm terms. This also links Willis material to higher order strain gradient theories [34]. We also evidenced different types of coupling terms due to the asymmetry, either absent when the unit cell consists in detuned identical type resonators or due to a physical asymmetry when the unit cell involves different types of resonators. Beyond the derivation of these closed-forms, we also show that various asymmetric structures can be modeled as Willis materials. Each effective parameter (effective density, effective bulk modulus, and effective Willis coupling parameter) is validated against experimental and numerical results, showing the robustness of the derivation method. This work paves the way for the engineering use of Willis materials, from perfect absorption devices in transmission problems, to double negative structures, but also to metaporous materials and liquid foams [35][36][37].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Padé's approximation versus Taylor's expansion
In the long-wavelength regime, the first order Taylor's expansion of the transfer matrix reads as which can be inverted to provide This expression does not ensure the reciprocity of the material and the diagonal terms are not opposite each other. Taylor's expansion, at least of first order, is not a sufficiently robust tool to derive effective parameters, when compared to Padé's approximation.

Appendix B. Elementary transfer matrices
This section details the different elementary transfer matrices that are used to evaluate the total transfer matrix T. In particular, the impedance of a side branch HR to a duct of radius r, assumed to be a point-like resonator, is [38]Z HR = −iZ n 1 − k nZ n Z c δ tan (k c l c ) −Z n Z c tan (k n l n ) tan (k c l c ) tan (k n l n ) − k nZ n Z c δ tan (k n l n ) tan (k c l c ) +Z n Z c tan (k c l c ) , where k n andZ n are respectively the wavenumber and reduced impedance of the neck, k c andZ c are respectively the wavenumber and reduced impedance of the cavity, and δ = 0.82(1 − 1.35r n /r c + 0.31(r n /r c ) 3 )r n + 0.82(1 − 0.235r n /r − 1.32(r n /r) 2 + 1.54(r n /r) 3 − 0.86(r n /r) 4 )r n is the correction length, with r n and r c the radius of the neck and of the cavity respectively. In the absence of losses and correction length, the latter impedance reduces tō Z HR = −iZ 1 − k 2 l n S c l c /S n k (S n l n + S c l c ) at low frequencies, from which it follows that 1/Z HR varies like ω when ω → 0, but like a Lorentzian around the resonance frequency.
where the expression of T l (2) and of T l (1) are given equation (B.1). Assuming kd λ, in such a way that k (1) l (1) = ζ (1) kd and k (2) l (2) = ζ (2) kd, with ζ (j) = l (j) /d, j = 1, 2, are also much smaller than the wavelength and making use of equation (7) leads to if we only Taylor expand the elements of the transfer matrix to the first order. In the opposite, it leads to if we Taylor expand the elements of the transfer matrix to the second order. Of particular interest is the fact that the eigenvalues of A e , i.e. the wavevectors, are identical for both orders of expansion (assuming respectively first and second order expansion) and equal to k e = ± iω √ K e /ρ e , while the eigenvectors exhibit Z ± e = K e (k e ± χ e ) at the second order against Z ± e = K e k e at the first order. In the lossless case, the Willis coupling parameter is purely imaginary.

Appendix D. Derivation of the Willis parameters in case of a unit cell presenting a single resonator
This section is motivated by nonlocal aspects of the Willis coupling and the fact that structures are often bounded in practice. Moreover, it helps in a better understanding of coupling terms that arise in case of more complex asymmetric unit cells.

Appendix D.1. Derivation of the Willis parameters in case of a unit cell presenting a single Helmholtz resonator
We consider a unit cell composed of a straight duct of length d and radius r loaded by an HR of reduced impedanceZ HR and located at l (1) such that d = l (1) + l (2) . The total transfer matrix reads as where the transfer matrix accounting for the propagation along each length and the HR are provided in equations (B.1) and (B.2). Assuming kd λ, in such a way that kl (1) = ζ (1) kd and kl (2) = ζ (2) kd, with ζ (j) = l (j) /d, j = 1, 2, are also much smaller than the wavelength, but also that 1/Z HR is a O (kd) 2 term (because it varies like ω at low frequencies), and making use of equation (7) leads to At first glance, the effective parameters appear as the sum of those of a laminated two-material unit cell, equation (C.3) and additional terms related to the presence of the HR. The main tube is identical along the segments 1 and 2, thusK (1) =K (2) andρ (1) =ρ (2) making the Willis parameter associated with the two materials to collapse. The effective bulk modulus and density then read as those proposed in [26,27]. The presence of the HR affects the bulk modulus, which can become negative around the HR resonance. The Willis parameter only accounts the possible phase shift, when a material is bounded by bounds that do not coincide with the natural unit cell bounds, i.e. when the HR cannot be centered in the unit cell. Please note that the Willis parameter simply vanishes when l (1) = l (2) = d/2. Nevertheless, the term l (2) − l (1) translates a momentum within the unit cell. While ρ e and 1/K e are also O (kd) 2 terms, χ e is a O (kd) 3 term, when 1/Z HR is assumed to be a O (kd) 2 term. To prove that 1/Z HR should be a O (kd) 2 term (at least in the low frequency), we adopt a reductio ad absurdum. Let us assume a centered HR (i.e. the Willis parameter vanishes) and that 1/Z HR is a O (kd) term for example. The equation (D.2) are thus derived from a first order Taylor expansion of the elements of the total transfer matrix. Under these assumptions, det (T) ≈ 1 − iωρd/Z HR + O (kd) 2 = 1 does not satisfy the reciprocity condition. Thus, 1/Z HR is assumed to be a O (kd) 2 term. Note that this relies more on a frequency analysis rather than on a kd analysis. We can therefore use the power of ω instead of the power of kd to express the order of Taylor expansion in the following.
Turning the sample over and repeating the procedure (second set of measurements are mark with), gives another set of equations

(E.2)
These two sets of equations form an overdetermined system to solve for R + , R − , and T.

Appendix F. Recovery procedure of the effective parameters from the measured reflection and transmission coefficients
The two state vectors at both sides of the sample are related by where diag is the diagonal matrix, Σ ± are the eigenvalues of A and V the corresponding eigenvector matrix. For a Willis material, the constitutive matrix of which is given in equation (2), Σ ± = ±iω χ 2 + ρ/K = ±iωσ, Introducing r ± = Z + Z ± / Z − Z ± , where Z ± = K (σ ± χ), the reflection coefficients at the interface between a semi-infinite background medium (of impedance Z) and the Willis material in the direct and reverse orientations, these equations can be inverted to yield The sign of the first two equations are chosen to satisfy the passivity condition and χ e , ρ e , and K e are subsequently evaluated. These equations are slightly different from those given in [9] and extend the method proposed in [41,42], which was already used in [43].

Appendix G. Direct numerical calculation of the effective properties from the total transfer matrix
Once the total transfer matrix T is calculated, it is directly assimilated to exp A num