Cyril Cano

Éric Chassande-Mottin

Nicolas Le Bihan

FAST AND ACCURATE GRAVITATIONAL-WAVE MODELLING WITH PRINCIPAL COMPONENT REGRESSION

Keywords:

Inference from gravitational-wave observations relies on the availability of accurate theoretical waveform models to compare with the data. This paper considers the rapid generation of surrogate time-domain waveforms consistent with the gravitational-wave signature of the merger of spin-aligned binary black holes. Building on previous works, a machinelearning model is proposed that allows for highly-accurate waveform regression from a set of examples. An improvement of about an order of magnitude in accuracy with respect to the state of the art is demonstrated, along with a significant speed up in computing time with respect to the reference generation software tools.

MOTIVATIONS

The first observation of gravitational waves (GW) by the LIGO/Virgo collaboration in 2015 [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF] marked the advent of gravitational-wave astronomy. Since this date about 50 events have been detected [START_REF] Abbott | GWTC-2: Compact binary coalescences observed by LIGO and Virgo during the first half of the 3rd observing run[END_REF] associated to the merger of compact star binaries, primarily binary black holes, the focus of this work.

Theoretical waveform models (or templates) are used to detect the gravitational-wave signals using matched filtering [START_REF] Abbott | GWTC-2: Compact binary coalescences observed by LIGO and Virgo during the first half of the 3rd observing run[END_REF], or to infer the astrophysical parameters of their source using Bayesian samplers [START_REF] Thrane | An introduction to bayesian inference in gravitational-wave astronomy: Parameter estimation, model selection, and hierarchical models[END_REF]. For both those tasks, a large number of template waveforms (about 10 5 to 10 6 ) are required to be computed to cover the relevant parameter space.

Waveform models are deduced from the resolution of the source dynamics which is a difficult relativistic problem (see [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF] for a recent review). The evaluation of recent and accurate waveform models is computationally expensive and actually dominates the computational budget for parameter inference.

In the next decade, LIGO and Virgo detectors are expected to conduct at least two major observing runs with improved sensitivity, leading to a large increase in the number of detected signals. The analysis of those future observations calls for numerically efficient, yet accurate waveform generators.

The waveform morphology being reasonably smooth and slowly varying with respect to its parameters (i.e., the binary component masses and their spins), it is possible to fit a signal model based on a generic amplitude and phase evolution from a set of example waveforms. We aim at reproducing noisefree waveform models with high fidelity thanks to a fast learning algorithm. The targetted accuracy is 10 -3 in mismatch, a figure-of-merit defined later in Sec. 3.5. This has been successfully realized using reduced-order modelling (based on tensor spline fitting) [START_REF] Pürrer | Frequency domain reduced order model of aligned-spin effective-one-body waveforms with generic mass-ratios and spins[END_REF][START_REF] Lackey | Surrogate model for an alignedspin effective one body waveform model of binary neutron star inspirals using Gaussian process regression[END_REF] and more standard machine learning approach such as mixture-of-experts regression [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] and artificial neural networks [START_REF] Khan | Gravitational-wave surrogate models powered by artificial neural networks[END_REF][START_REF] Chua | Reduced-order modeling with artificial neurons for gravitational-wave inference[END_REF].

In this work, we propose a model with reduced complexity based on principal component regression, that is able to improve the overall regression accuracy by about an order of magnitude. It can be implemented using off-the-shelf algorithms from Scikit-learn software library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. These performances are notably achieved thanks to a different choice of features.

GRAVITATIONAL-WAVE WAVEFORMS

Gravitational waves are oscillations of the space-time curvature propagating through space [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF]. So far, the only sources detected are coalescences of compact star binaries, composed of either black holes and/or neutron stars. Here we focus on "spin-aligned" BBH where the spins S i ∈ R 3 for i = 1, 2 of both objects are aligned with the orbital angular momentum L i.e., the z axis normal to the orbital plane. In this case, the parameter space reduces to four parameters {m 1 , m 2 , S 1z , S 2z } [START_REF] Bohé | Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors[END_REF] with m i for i = 1, 2 the mass of each black hole. The gravitational-wave signal is also determined by extrinsic parameters such as the luminosity distance r in megaparsecs (Mpc), i.e. the distance between the source and the observer, and the direction of the line of sight parametrized by the two angles (ι, ϕ), namely the inclination of the source and the initial phase, see Fig. 1.

General relativity predicts the existence of two polarization modes for gravitational waves, the plus and cross polarizations, denoted h + (t) and h × (t). It is customary to express the two polarizations as a single complex-valued time series h(t) = h + (t) -ih × (t) referred to as gravitationalwave strain.

The computation of h(t) needs to resolve the source dynamics. There are no exact close-form solution for this relativistic problem, however accurate approximations are available. In the spin-aligned case and for binaries with low mass ratio q (defined as q = m 1 /m 2 and with q ≥ 1) and low total mass M = m 1 + m 2 , the following expression [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF] can be used for the strain

h(t) ∝ M r a(t) 1 + cos 2 ι 2 cos(Φ(t) -2ϕ) -i cos ι sin(Φ(t) -2ϕ)] (1)
where a(t) and Φ(t) stand for the amplitude and phase of the dominant mode in the spin-weighted spherical harmonics expansion of h(t) (see [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF] for details). General relativity's scale invariance implies h(t; m 1 , m 2 ) = h(λt; λm 1 , λm 2 ) [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF]. As a consequence, the total mass M * = m 1 + m 2 can be fixed to a fiducial value 1 and in the case of spin-aligned BBH, waveforms can be described using a reduced number of three intrinsic parameters, e.g., {q, S 1z , S 2z }. The extrinsic parameters r * , ι and ϕ correspond to simple scaling or phase factors that can be applied a posteriori in the waveform computation. We thus fix r * = 1 Mpc, and (ι * , ϕ * ) = (0, 0) in the sequel, with no loss of generality. With such a setting, the GW signal h(t) ∝ a(t)e -iΦ(t) appears as an amplitude and frequency modulated signal, often referred to as "chirp" for short, with 1 Here, we used M * = 20 M with M the solar mass. a(t) and Φ(t) its amplitude and phase parameters exhibiting non-oscillatory and smooth behaviours (see Fig. 2).

A range of models allows the generation of approximated yet accurate waveforms for data analysis purposes [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF]. Here, we use SEOBNRv4 [START_REF] Bohé | Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors[END_REF] based on the effective-one-body formalism (EOB). Fig. 2 shows an example of the expected gravitational-wave signal from a BBH merger.

WAVEFORM GENERATION MODEL

Following [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] we now introduce a generative model able to regress waveforms given θ = {m 1 , m 2 , S 1z , S 2z } the set of features that collects the input astrophysical parameters. The proposed regressor is a principal component regressor [START_REF] Hastie | The elements of statistical learning[END_REF], that consists in regressing the PCA coefficients of the attributes from polynomial combining of the features. The waveform generation essentially performs the reverse process: obtain estimates for the PCA coefficients from the regressor, that are inverted to compute attributes, which in turn are mapped to waveforms.

Mapping to attributes

Over its entire duration, the chirp signal h(t) goes through a sequence of three phases associated to different dynamical regimes: inspiral when the two black holes are far apart, merger when they are close and "plunge" onto each other and ringdown when the merger remnant settles down to equilibrium [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF]. The amplitude and phase evolve over different timescales during these phases. To capture their variations with a uniform accuracy over the entire waveform duration, amplitude and phase are discretized in time with a varying sampling resolution.

Time grid goes from t start = -20 s to t end = 0.006 with t = 0 at maximum amplitude. Concretely, 4096 sample points are computed over the non-uniform time grid t = sign (t) |t| 1 α where α = 0.35 [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF]. The sampling rate is reduced in the inspiral phase where both the amplitude and phase evolutions are slow, while it is increased when getting closer to the merger where the amplitude reaches its maximum, and where the phase evolves rapidly. This resampling has two benefits: data reduction during inspiral2 , and high accuracy near the merger.

Data reduction

The training set is reduced through a principal component analysis (PCA) of each waveform attribute after an alignment to 0 at t = t start for PCA efficency purposes 3 . The minimal number of principal components (PC) required to reach maximal capabilities of the model is determined by crossvalidation (see Sec. 4.2).

Regression

Schmidt et al. [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] use a mixture of experts (MoE) regressor to infer the PCA coefficients from the set of features {q, χ 1z , χ 2z } where χ i = S i /m 2 i for i = 1, 2 are the dimensionless spins. MoE is an ensemble learning method, based on a weighted sum of linear regressors called "experts" [START_REF] Hastie | The elements of statistical learning[END_REF].

In this work, a single linear regressor is preferred and applied to a different set of features, expanded with polynomial combinations up to a pre-determined order. This maximum order is computed to maximize a score, see Sec. 4.2.

Choosing the initial set of features has a significant impact on the final score and regression accuracy. To leading order, the amplitude and phase evolution are known [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF] to depend on the chirp mass M = (m 1 m 2 ) 3/5 /M 1/5 , the mass ratio q and the effective spin χ eff = (qχ 1z + χ 2z )/(1 + q). It is natural to think that those physically motivated parameters are good candidates to fit the data.

Systematic tests with subsets made of different feature combinations were performed and the final model use the feature set leading to the best score (see Sec. 4.2.2).

Waveform generation

Once the model is trained, the usage goes as follows: the PC coefficients of the attributes are predicted for the parameters of the desired compact binary. The attribute time series are deduced from the predicted PC coefficients and the waveform amplitude and phase are interpolated from the attributes on the (uniform) time grid used for data analysis. The GW polarizations h + (t) and h × (t) are finally computed by applying (1) with the requested distance, inclination and initial phase.

Scoring

The regression accuracy is evaluated by a specific metric called mismatch or unfaithfulness [START_REF] Chatziioannou | Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals[END_REF]. The mismatch between two waveforms h, g ∈ C N is defined as:

mismatch (h, g) = min τ ∈R 1 - | h τ , g | h τ g (2)
where h τ (t) := h(t-τ ) and with the scalar product:

f, g = h(f )g * (f ) S(f )
df defined in the Fourier domain and with f being the frequency variable. This metric is a loss function (smaller is better) and it is phase-shift and time-shift invariant. The definition in (2) allows for a frequency-dependent weighting, usually fixed to the GW detector noise power spectrum density S(f ). Here, we assume a flat noise curve S(f ) = 1, which leads to a conservative constraint as the regression should then be equally accurate at all frequencies.

Errors in the waveform approximation lead to systematic errors in the astrophysical parameters estimates obtained from the observations. Those systematic errors from mis-modeling should be smaller than the statistical errors (due to the presence of noise in the observations). This principle leads to the rule of thumb (see e.g., App. G of [START_REF] Chatziioannou | Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals[END_REF]) stating that the mismatch should be < N/(2 SNR 2 ) where N is the effective number of intrinsic parameters and SNR is the signal-to-noise ratio. Schmidt et al. [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] achieve a median mismatch value of 5 × 10 -4 , with tails going to 10 -1 in the worst case. In the spin aligned case with N = 3 effective parameters, this corresponds to an applicability range that goes up to SNR = 54 (3 in the worst case).

RESULTS

Training and testing datasets

A dataset of 4000 randomly distributed BBH waveforms was simulated with the SEOBNRv4 model from the LALSimulation software library [START_REF]LIGO Algorithm Library -LALSuite[END_REF]. The mass ratio q was randomly drawn uniformly in the range [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF]20] and spins χ 1z , χ 2z uniformly over [-0.8, 0.95]. Such dataset is comparable to the one used in [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF].

For the subsequent analysis, the dataset is splitted into a training and a testing set that corresponds to 80 % and 20 % of the main dataset respectively. All reported results are obtained with the testing set.

Hyperparameter tuning

Number of PC after truncation

The overall modelling error is essentially due to the truncation of the PC expansion and the regression for each attribute. We first evaluate the former in Fig. 4 (circles and boxes). The phase approximation plays a dominant role. For instance, the truncation to only one PC leads to a mismatch of 5 × 10 -5 for the amplitude and 7 × 10 -1 for the phase.

Fig. 4 shows that the overall error after regression (crosses and plus signs) stabilizes at 6 PCs for the phase with a mismatch median score of about 10 -5 . This is the retained number of PCs after truncation as this indicates the regressor fails to estimate higher order PC coefficients. For simplicity, the PCA expansion is truncated to the same number of PCs for both the amplitude and phase. Fig. 3. Median mismatch vs number of principal components retained after PCA truncation. The mismatch is computed after fitting the amplitude only (phase is exact) for circles and plus signs, and the phase only for boxes and crosses. Circles and boxes represent mismatches solely due to the PCA truncation whereas crosses and plus signs represent mismatches due to the overall model. The dashed line indicates the selected PC truncation level at 6 PC.

Feature selection

Different feature combinations were tested and compared through their mismatch median value and dispersion. Promising sets of features were first pre-selected based on their r 2 score [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] obtained for the regression of the first PC of the phase, which essentially governs the overall performance (see Sec. 4.2.1). This procedure was performed for subsets that collect up to 6 parameters from the pool that includes q, M, χ eff and χ iz , m i and 1/m i for i = 1, 2.

We found that about twenty feature sets lead to the best median mismatch of order 10 -5 . Among those sets, the two following {χ 2z , χ eff , M} and {χ 1z , χ 2z , q, m 2 } are of particular interest. The first has the nice property to possess three features only, similarly to the number of intrinsic physical parameters. However, the mismatch distribution obtained from the testing set has much smaller variance with the second feature set. Using this feature set it is possible to produce many of the terms that appear in the so-called post-Newtonian expansion that approximates the signal phase before the final merger [START_REF] Schmidt | Gravitational waves from binary black hole mergers: Modelling and observations[END_REF] (e.g., the component spins S 1z , S 2z before normalization, or the reduced mass m 1 m 2 /(m 1 + m 2 )), thus explaining the good performances. Interestingly feature sets such as {q, χ 1z , χ 2z } (used in [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF]) or {q, χ eff , M} (a "natural" candidate composed of the physically motivated parameters) do not perform as well, with a median mismatch of 3 × 10 -2 and 1 × 10 -1 respectively.

Maximum order of polynomial features

The feature set {q, m 2 , χ 2z , χ eff } is expanded by polynomial combinations of the initial features up to a maximum order [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The minimum order necessary to achieve best performances was determined by cross-validation. We found that a 7 th order polynomial feature expansion is sufficient to minimize the median mismatch.

Accuracy and runtime benchmarking

Fig. 5 shows the mismatch distribution obtained with the testing set using {χ 1z , χ 2z , q, m 2 } as feature set. The median mismatch is 1.8 × 10 -5 (average is 6.8 × 10 -5 ) and the 5% and 95% percentiles are 2.2 × 10 -6 and 1.6 × 10 -4 , resp. The worst case mismatch is ∼ 10 -2 . As shown in Fig. 5 we were not able to reproduce the performances given [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] with the Schmidt et al. model. Instead we found a median mismatch of 2.3 × 10 -3 on the testing set.

No matter which reference is considered ( [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF] or our reproduction here), we obtain an improvement by at least an order of magnitude compared to the Schmidt et al. model. This makes the approximation error much smaller than the numerical and modelling error intrinsic to the SEOBNRv4 model [START_REF] Bohé | Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors[END_REF], uniformly over the considered parameter space. Following the rule given in Sec. 3.5, this corresponds to an applicability range that goes up to SNR ≈ 288 (12 in the worst case) which covers the loudest SNR expected during the up-coming LIGO and Virgo science runs.

We evaluated the speed up factor defined as the ratio between the runtimes required to compute the original waveform model 4 (using the LALSimulation software library [START_REF]LIGO Algorithm Library -LALSuite[END_REF]) to that of the proposed regression model. To perform this measurement, a set of 500 waveforms was generated by drawing the total mass uniformly between 40 and 100 M with the other properties as described in Sec. 4.1. The median speed up is measured to be ∼ 10 2 (comparable to [START_REF] Schmidt | Machine learning gravitational waves from binary black hole mergers[END_REF]). This shows that the proposed method can greatly accelerate parameter estimation with Bayesian samplers [START_REF] Thrane | An introduction to bayesian inference in gravitational-wave astronomy: Parameter estimation, model selection, and hierarchical models[END_REF]. The runtime of the generation with the regression model is dominated by the interpolation step from the non-uniform time grid of attributes to the uniform time grid of the user. The speed up factor could thus be further improved by dropping the non-uniform time grid (and thus the need for interpolation). Preliminary tests indicate that the same accuracy level can be obtained without the non-uniform time grid.

PERSPECTIVES

This work paves the way to further extensions relevant to gravitational-wave astronomy, such as the inclusion of higher/subdominant modes other than ( = 2, m = ±2). Another important extension is the case of precessing binaries (with arbitrary, misaligned spins) leading to arbitrarily polarized waveforms (not necessarily circularly polarized as it is the case here). In this more complex case with more input parameters and a larger waveform variability, fast and highly accurate modelling is still an open problem.

While the method presented here is evaluated in the context of gravitational-wave astronomy, it applies in principle to any area that involves chirp signals (see e.g., [?] for examples). One of the lessons learned from this study is that pragmatic regression approaches based on off-the-shelf tools may provide better results than more sophisticated approaches. A judicious choice of features appears critical to optimize the approximation accuracy.

2 rFig. 1 .

 21 Fig.1. Physical parameters of a compact star binary -Intrinsic parameters {m 1 , m 2 , S 1 , S 2 } (component masses and spins) and extrinsic parameters such as the direction of the line of sight (ι, ϕ) and the distance r. For spin-aligned binaries component spins are aligned with the orbital angular momentum J so the intrinsic parameters reduce to {m 1 , m 2 , S 1z , S 2z }.

Fig. 2 .

 2 Fig. 2. Top: example of gravitational-wave waveform obtained with the SEOBNRv4 model for a spin-aligned BBH. Bottom: associated amplitude and phase waveform attributes obtained from the waveform following Sec. 3.1.

Fig. 4 .

 4 Fig. 4. Mismatch histograms obtained with the proposed principal component regression model (orange) and that of Schmidt et al. (purple) [7].

This allows to easily process waveforms sampled at 16384 Hz on a standard laptop.

[START_REF] Thrane | An introduction to bayesian inference in gravitational-wave astronomy: Parameter estimation, model selection, and hierarchical models[END_REF] The amplitude and phase offsets subtracted by the alignment procedure can be fitted and added back at the generation stage to produce the full waveform. Though we don't detail this part here but this can be done with good accuracy with the same regressor.


Côme Olivier

Gaëlle Poignand

Matthieu Malléjac

Vicente Romero-García

Guillaume Penelet

Aurélien Merkel

Daniel Torrent

Jensen Li

Johan Christensen

Jean-Philippe Groby

Non-reciprocal and even Willis couplings in periodic thermoacoustic amplifiers

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Swapping source and receiver has long provided identical frequency responses due to reciprocity features, thus preventing the possible tuning of the transmission coefficient in opposite directions. In acoustics and elasticity, the reciprocity of wave propagation can be broken via three means: by using spatiotemporal dependent material properties 1 , by combining nonlinear properties with an asymmetry 2-4 , or by introducing an external bias 5,6 .

Non-reciprocal devices and their applications have recently been reviewed in Ref. [7]. Among the various possibilities offered by the introduction of an external bias in acoustic systems, thermoacoustics, although effective 8,9 , is often neglected. Thermoacoustics is concerned with the thermal interaction between an acoustic wave and the surrounding waveguide wall [START_REF] Swift | Thermoacoustics: A unifying perspective for some engines and refrigerators[END_REF][START_REF] Garrett | [END_REF] .

Although research in thermoacoustics has mostly been oriented towards the development of thermodynamic machines and the understanding of the phenomena controlling their behavior -on the one hand, the thermoacoustic prime-movers (or engines) 12,13 ; on the other hand, thermoacoustic heat-pumps and refrigerators 14,15 -more fundamental research has relied on the thermoacoustic amplification or damping of acoustic waves in recent years. Quasiperiodic and chaotic oscillations 16,17 or synchronization phenomena involved in thermally driven autonomous oscillators 18,19 , thermoacoustic shock waves 20,21 and solitary waves 22 , thermoacoustic diode 8 or a PT -symmetric system 23 have been the focus of a growing interest. Thermoacoustics still appears to be overlooked as far as non-reciprocal systems are concerned.

Meanwhile, Willis materials have received an increasing interest since the publication of the seminal work 24 , thanks to their analogy with bi-isotropic electromagnetic metamaterials [START_REF] Lindell | Electromagnetic waves in chiral and bi-isotropic media[END_REF] . The Willis coupling parameters couple the potential and kinetic energies in the acoustic conservation relations, therefore enhancing the ability to control waves in these metamaterials compared to other materials that do not exhibit such coupling. Three categories of Willis coupling parameters have been highlighted: the even coupling related to the structure asymmetry, the odd coupling related to non-local effect (although its existence is questioned for one-dimensional systems), and the non-reciprocal coupling [START_REF] Sieck | [END_REF] . Only a few non-reciprocal systems 1, [27][28][29] have been modeled and analyzed as Willis materials so far. Interestingly, non-reciprocal Willis coupling was found to be non zero already in the first-order Taylor expansion of the solution, while even Willis coupling was found to be non zero in the second-order Taylor expansion of the solution. None of these articles derived closed form expressions of the coupling parameters and investigated the possible induced asymmetry by the non-reciprocity.

In this article, we derive the closed form expressions of the effective properties of a one-dimensional periodic arrangement of thermoacoustic amplifier following Ref. [30]. The procedure relies on the Padé's approximation of the total transfer matrix that links the state vectors at both sides of the unit cell and directly provides both even and non-reciprocal Willis coupling terms. The analysis of the effective parameters as well as the dispersion relation show that the system possesses a coalescence point in the k-space 31 that resembles a PT -broken phase and yields to zero-group-velocity for given temperature differences.

II. GENERAL STATEMENT

A. Wave propagation in a duct submitted to a temperature gradient Assuming an implicit time dependence e -iωt and following Refs. [10] and [32], the equations that describe the acoustic wave propagation in a duct of cross-sectional area S in the presence of a longitudinal temperature gradient ∂T m /∂x are:

∂p ∂x = iωρ m Sφ(1 -f ν ) V , ∂V ∂x = iωSφ γp m [1 + (γ -1)f κ ]p + f κ -f ν (1 -f ν )(1 -Pr) 1 T m ∂T m ∂x V , (1)
where p is the acoustic pressure, V = SφV is the volume flow rate, with V the particle velocity and φ the porosity across the waveguide, ρ m , p m , and T m are the mean-state values (across the section) of the density, pressure (atmospheric pressure), and temperature respectively, f ν and f κ are complex frequency-dependent functions, which account respectively for viscous and thermal losses 10 (see also Appendix A), γ is the specific heat ratio, and Pr is the Prandtl number. Introducing the state vector W =< p, V > T , where T is the transpose operator, these equations can thus be cast in the matrix form

∂ ∂x W =   0 iωρ(x) iωC(x) G (x)   W = A(x)W, (2)
where ρ(x) is the effective density, C(x) is the effective compressibility, i.e., the inverse of the effective bulk modulus K(x), and G (x) is a 'gain' term depending on the temperature gradient ∂T m /∂x, that represents the thermoacoustic amplification/attenuation. The choice of the volume flow rate as second component of the state vector instead of the particle velocity is motivated by the continuity conditions (continuity of pressure and volume flow rate) at interfaces between different materials and ducts of different cross-sections. Note that the effective bulk modulus and density are thus both divided by the waveguide cross-section because of the use of the volume flow rate as second component of the state vector. Thus, the effective bulk modulus and density are bulk modulus per unit area (given in Pa.m -2 ) and density per unit area (given in kg.m -5 ). Note in addition that the effective bulk modulus and density are spatially-dependent in an implicit way (but obviously frequency-dependent), notably because the fluid constants that appear in both functions f ν and f κ depend on the temperature distribution T m (x). The propagation matrix A depending on x and not commuting for different values of x, the solution of the system Eq. (2), which relates the state vector at a position l to that at a position 0 via W(l) = M l W(0) involves a matricant M l that is evaluated through Peano series expansion 33

M l =Id + l 0 A(x)dx + l 0 A(x) x 0 A(ζ)dζ dx + • • • . (3)
Of particular interest is the first-order expansion, which reads as

M (1) l = Id + l 0 A(x)dx =   1 iω ρl iω Cl 1 + Ḡ l   , (4)
where ρ, C, and Ḡ are the mean values of ρ(x), C(x), and G (x) over the length l, i.e.,

• = 1 l l 0 •dx. This first-order iteration corresponds to a first-order Taylor expansion of exp(Al), when A does not depend on x. Introducing k = ω ρ C, the determinant of M

(1) l reads in the first-order expansion det M

(1) l

= 1 + Ḡ l + O ( kl) 2 , (5)
which clearly shows that the system is not reciprocal in the presence of the 'gain' term, i.e., when ∂T m /∂x = 0. The term Ḡ l is a marker of the non-reciprocity that makes det M

(1) l = 1. When Ḡ = 0 (i.e., ∂T m /∂x = 0), the system falls back to reciprocal with det M

(1) l = 1. The term G = 1 + Ḡ l is thus called non-reciprocity term in the following. In reality, Eq. ( 4) relies on two assumptions and is O ( kl) 2 , i.e., the system length is much smaller than the wavelength, but also O ( Ḡ l) 2 , i.e., the 'gain' per unit length is weak.

Effectively, the dimensionless term Ḡ l mainly depends on the temperature gradient and is thus driven independently of kl. Note that Eq. ( 5) turns out to be O ( kl) but O ( Ḡ l) 2 .

The determinant does not depend on kl and only depends on Ḡ l, thus emphasizing the central role of Ḡ l in the system non-reciprocity.

B. Calculation of the effective properties from the knowledge of the transfer matrix of the unit cell

We assume a one-dimensional non-reciprocal and asymmetric system composed of a d-

periodic
where ρ e is the effective density, C e is the effective compressiblity, χ a e is the even Willis coupling, and χ nr e is the non-reciprocal Willis coupling respectively equal to

ρ e = -4i ωd t 12 1 + Tr(T) + det(T)
, C e = -4i ωd

t 21 1 + Tr(T) + det(T) , χ a e = -2i ωd t 11 -t 22 1 + Tr(T) + det(T) , χ nr e = -2i ωd det(T) -1 1 + Tr(T) + det(T) . (8)
When considering reciprocal systems, det(T) = 1 and thus χ nr e = 0. The expressions of the effective parameters then reduce to those provided in Ref. [30]. In addition, for symmetric systems t 11 = t 22 and as a consequence χ a e = 0.

III. EXPRESSION OF THE EFFECTIVE PROPERTIES OF A PERIODIC THER-MOACOUSTIC AMPLIFIER

We consider a d-periodic thermoacoustic amplifier, the unit cell of which consists of five elements in series of respective lengths l j , j = 1, • • • , 5, such that d = 5 j=1 l j , see Fig. 1(a). The l 2 -long porous material, where the thermoacoustic effect takes place (customarily called stack), is subjected to a temperature gradient applied with the help of an ambient and a hot heat exchangers (denoted AHX and HHX) of respective lengths l 1 and l 3 . The heat exchangers impose a temperature difference ∆T = T H -T C , between the hot T H and the ambient T C temperatures. The HHX is followed by an empty tube of total length l 4 + l 5 . The temperature distribution loops back to the temperature of the AHX along the l 4long portion, which is therefore called the thermal buffer tube (TBT). In this study, the resulting temperature distribution is considered one-dimensional, homogeneous in the heat exchangers, and varies linearly in the stack and the TBT, as depicted on Fig. 1(a). The temperature gradient ∂T m /∂x is thus positive along the stack, while it is negative along the TBT, here also with reference to Fig. 1(a). No temperature gradient is applied on the elements 1, 3, and 5. These elements are thus reciprocal. The state vector at the end of the unit cell W(d) is related to that at the other end W(0) via the full transfer matrix T which is the multiplication of the five matrices modeling the propagation along the five elements of the unit cell

W(d) = TW(0) = T l 5 M
(1)

l 4 T l 3 M
(1)

l 2 T l 1 W(0), (9)
where T l = exp (Al) is the usual transfer matrix along a distance l with A the propagation matrix of the element, the first-order Taylor expansion of which takes the form

T l =   cos(kl) iZ sin(kl) i sin(kl)/Z cos(kl)   ≈   1 iωρl iωCl 1   + O (kl) 2 , ( 10
)
where k is the wavenumber and Z is the impedance (divided by the duct cross-sectional area in our case). Equation (10) formally corresponds to Eq. ( 4) in the absence of 'gain' term. Viscothermal losses are accounted for in each element, see Appendix A. Considering the first-order Taylor expansion of the total transfer matrix elements, we end up with

ρ e = 2 d [1 + G 2 G 4 ] (ρ 1 l 1 + ρ2 l 2 + [ρ 3 l 3 + ρ4 l 4 + ρ 5 l 5 G 4 )] G 2 ) , C e = 2 d [1 + G 2 G 4 ] C 5 l 5 + C4 l 4 + C 3 l 3 + C2 l 2 + C 1 l 1 G 2 G 4 , χ a e = -i [1 -G 2 G 4 ] dω [1 + G 2 G 4 ] , χ nr e = i [1 -G 2 G 4 ] dω [1 + G 2 G 4 ] , (11)
where all the effective parameters are in O ((kd) 2 ).

This result calls for several comments.

First, all the effective parameters are impacted by the non-reciprocal terms, notably via a common denominator 1 + G 2 G 4 . This denominator arises from the assumption that both Ḡ2 and Ḡ4 are O ((kd)). This implies that ρ e S/ρ 0 < 1 and K e S/γP 0 < 1 at low frequencies, where ρ 0 and P 0 are respectively the ambient (air) density and atmospheric pressure. Please note that Eqs. ( 11) also rely on the assumption that Eq. ( 4) is

O ( Ḡ d) 2 . Thus, G 2 G 4 = 1 + Ḡ2 l 2 + Ḡ4 l 4 + O ( Ḡ d) 2 rigorously.
Second, the thermoacoustic amplifier appears as a Willis material and exhibits both even and non-reciprocal coupling parameters at the first-order Taylor expansion. This constitutes a major difference with respect to laminate or resonant Willis materials, which require second-order Taylor expansion to exhibit Willis coupling [START_REF] Sieck | [END_REF]30 .

Third, χ a e = -χ nr e , which implies that the top left element of the effective propagation matrix is null in the first-order expansion. The propagation matrix thus strongly resembles Eq. (2). This suggests that the thermoacoustic fundamental equations, Eq. (1), are a specific form of those of a non-reciprocal Willis material. Together with the second comment, we conclude that the even Willis coupling arises from the asymmetry of the temperature gradient and appears in reaction to the non-reciprocity. This becomes clear when both temperature gradients are of opposite sign so as to lead to Ḡ2 l 2 = -Ḡ4 l 4 . In this case,

G 2 G 4 = 1+O ( Ḡ d) 2
and thus χ a e = -χ nr e = 0 in the first-order expansion, while the system material organization is asymmetric. Both Willis coupling parameters only depend on the temperature gradient and not on the material asymmetry in the first-order expansion.

Fourth, when Ḡ2 and Ḡ4 are both equal to zero, the effective parameters collapse to those of a usual laminate fluid in the first-order expansion, with χ a e = χ nr e = 0. Note that the Willis even coupling term does not vanish in that case when second-order expansions are used because the unit cell involves different materials and is asymmetric 30 .

Fifth, the form of both ρ e and C e exhibits different cumulative impacts of the temperature gradients. While the temperature gradient impacts the following density terms, i.e., the density terms following a non-reciprocal element are multiplied by the product of all the preceding non-reciprocity terms, it impacts the preceding compressibility terms, i.e., the effective compressibility terms preceding a non-reciprocal element are multiplied by all the following non-reciprocity terms. The effect of the temperature gradient on the density derives directly from its flow source nature. Its effect on the compressibility is more complex and can be analyzed in the way that the temperature gradient acts as a volume flow amplifier rather than a pure flow source, with the volume flow rate depending on the previous compressibilities.

Sixth

IV. EXPERIMENTAL AND NUMERICAL VALIDATIONS OF THE EFFECTIVE PARAMETERS AND DISCUSSION

The experimental set-up consists of a waveguide of rectangular cross-section (2a × 2b = 1 cm × 5 cm) containing either a single thermoacoustic cell or 15 identical thermoacoustic cells of length d = 19.5 cm in both cases, as depicted in Fig. 1(a-c). Only the results for the single-unit cell are discussed in this section, while those for the 15 identical cells are discussed in Appendix D. The stack is a cordierite substrate with a pore density of 400 CPSI, of length l 2 = 15 mm. Its pores are square of half-width r s = 585 µm and its porosity is φ = 0.85.

The ambient heat exchanger is made of a section of 1/8 " aluminum honeycomb (porosity φ = 0.945), with l 1 = 15 mm, used to ensure a uniform temperature distribution over the cross-section. The hot heat exchanger is a section of the stack material of length l 3 = 3 mm supporting a resistive wire providing sufficient heating power to maintain temperature up to

T H = 373 K.
The TBT is assumed to be of length l 4 = 23 mm. The thermoacoustic system is placed between two waveguides of the same cross section for the measurement of the scattering matrix by means of the two-source method [START_REF] Munjal | [END_REF] . As illustrated in Fig. 1, each of these measurement sections consists of a straight 1.2 m duct, on which a pair of 1/4 " microphones are flush mounted, a moving coil loudspeaker is side mounted in a small enclosure, and ends with an anechoic termination. Each microphone pair is separated by a distance of 65 cm, allowing measurements in the 10 Hz to 200 Hz range, provided that they have been precisely calibrated (both in amplitude and phase) to limit the errors on the measured and calculated quantities 35 . A relative calibration is therefore realized using a homemade small cavity coupler. The low frequency anechoic terminations attached to the measurement sections consist of a larger cavity (1.5 cm × 7 cm) with an adjustable length closed by stainless steel wire mesh screens 36 . The reflection coefficient of such terminations has an amplitude lower than 5% in the frequency range of interest, which is a sufficiently low value for the needs of this study. The sources are driven at amplitudes that generate sufficiently low acoustic pressures to ensure a linear system behavior during all measurements. -open circles-for ∆T = 0 K (blue curves), ∆T = 22 K (cyan curves), ∆T = 42 K (green curves), ∆T = 75 K (orange curves), and ∆T = 100 K (red curves). Note that the direct calculation with the Transfer Matrix Method is conducted by discretizing the stack and the TBT (subjected to temperature gradients) in piecewise constant elements. All curves are found in good agreement. Slight discrepancies are nevertheless noticed, notably for the normalized bulk modulus (see Fig. 2(e)), which are attributed to the assumption of the linear variation of the temperature distribution, but also to the measurement difficulties at very low frequency. As pointed out in the previous section, the two Willis couplings are almost purely imaginary at low frequencies and are of equal modulus but opposite sign. Interestingly, the Willis coupling parameters do not vanish at low frequency as it is the case for laminate structures or detuned Helmholtz resonators 30 . When ∆T = 0 K, χ nr e = χ a e = χ nr num = 0, but χ a num = 0 although very small. This comes from the fact that χ a num captures the second-order expansion element which accounts for the material asymmetry 30 . While the numerically evaluated effective parameters are in good agreement with the measured parameters whatever the temperature gradient, the parameters given in Eq. ( 11) slightly deviate from both numerical and experimental results with increased temperature differences, mainly at low frequencies. This is particularly visible on the Willis coupling and is due to the use of the first-order Peano series expansion, for which the 'gain' term is approximated by its mean value. Effectively, the only use of the mean value can become too harsh an approximation when the temperature gradient becomes large. This is a direct translation of the assumption that Eq. ( 4) is O ( Ḡ l) 2 . Nevertheless, the overall trend is correctly captured. This is also visible in Figs. 3(a-d), which depict the absolute values of the reflection and transmission coefficients for each incidence, i.e., that where amplification occurs T + and R + and that where attenuation occurs T -and R -. Note that the transfer matrix is calculated via Padé's approximation when A e is used, thus validating the derivation of the expressions given in Eqs. (11) The results as calculated with the expressions given in Eq. ( 11) are depicted with solid curves, as numerically evaluated from the Transfer Matrix method are depicted with dashed curves, and as measured are depicted with the curves marked with circles.

The real and imaginary parts of the dispersion relations are depicted in Figs. 2(i-j) respectively. The higher the temperature difference, the more asymmetric k ± e , which translates the non-reciprocal nature of the configuration. The trajectories of both k ± e are modified when a temperature gradient is applied. In particular, both k ± e are shifted in the negative Im(k)-space at low frequencies. Thus, k + e is amplified while k - e is attenuated. When Im(k + e ) crosses back to the positive half Im(k)-space, i.e., no amplification is observed anymore, |T + | moves from values larger than 1 to values lower than 1 because the structure is almost acoustically transparent. This is highlighted by the magenta arrows, Figs. 2(j) and 3(a), when ∆T = 75 K, for example. A coalescence point in the k-space, where both k ± e are equal, is reached for a temperature difference around ∆T = 75 K as pointed out by the black arrows in Figs. 2(i,j). A zero-group-velocity point is reached at this frequency, where both modes are amplified, i.e., Re (k ± e ) > 0, while Im (k ± e ) < 0 at this frequency. This is due to the almost purely imaginary nature of both Willis coupling parameters as explained in Section III. This zero-group-velocity point is purely due to the even and non reciprocal Willis coupling parameters and does not translate into the effective compressibility and density, as assumed in Ref. [9]. In practice, the even Willis coupling parameter is usually purely imaginary at low frequency, see Ref. [30] for example, and thus a band gap can be opened at zero frequency. Nevertheless, the even Willis coupling is usually non zero in the second-order Taylor expansion for purely asymmetric laminated or resonant structures and usually vanishes notably in the case of Helmholtz resonators. This bandgap is thus rarely or never noticed. The story is completely different in non-reciprocal systems, in which the even coupling is non zero already in the first-order expansion and does not necessarily vanish at low frequencies. It is also accompanied by the non-reciprocal Willis coupling that shifts the mode in the Im(k)-space. This behavior was already noticed at much higher frequencies and at a higher pass band in Ref. [27]. This coalescence point in the k-space does not seem to translate into a particular behavior of the scattering matrix, although the reflection and transmission coefficients present inflection points at a frequency that is close. Interestingly, the k - e branch presents infinite group velocity points for lower frequencies and higher temperature differences, which find translation neither in the effective parameters, nor in the scattering elements.

V. CONCLUSION

A periodic one-dimensional thermoacoustic amplifier is analyzed as a Willis material.

The closed form expressions of the effective properties are derived from Padé's approxi-mation of the total transfer matrix that links the state vectors at both sides of the unit cell. A first-order Taylor expansion of the transfer matrix elements is sufficient to derive the effective properties notably both the even and non-reciprocal Willis couplings. The odd Willis coupling is absent, while the even Willis coupling is found to be only related to the asymmetry of 'gain', i.e., the temperature gradient that is applied to the unit cell, already in the first-order expansion. Even and non-reciprocal Willis couplings are found of equal modulus but opposite sign, which suggests that the even coupling appears as a counter reaction to the non-reciprocal coupling and that the thermoacoustics fundamental equations describe specific forms of Willis materials. We show that both Willis couplings are almost purely imaginary at low frequencies, thus enabling a zero-group-velocity point and more importantly the opening of an amplification band at vanishing frequency. The effective parameters and scattering properties of a single unit and a 15 units system are validated against experimental results. These results pave the way for a better physical understanding of Willis couplings in non-reciprocal systems, for easing the engineering application of Willis materials, and for various applications of such systems to further control the acoustic waves at very low frequencies. Concurrently, this article also aims at promoting thermoacoustics as an excellent mean to design non-reciprocal systems. thermal couplings between the oscillating gas and the surrounding solid walls [START_REF] Swift | Thermoacoustics: A unifying perspective for some engines and refrigerators[END_REF]37 and are defined for various geometries of the waveguide channels (i.e., the waveguide itself in the case of an empty duct such as the TBT, or the pores of a porous material, as in the case of the stack or the heat exchangers). Considering δ ν,κ the frequency-dependent viscous and thermal acoustic boundary layer thicknesses given respectively by δ κ = 2κ/ω and δ ν = √ Prδ κ , classical results for cylindrical channels of radius r (adequate approximation for the honeycomb heat exchangers) yield

f ν,κ = 2 (i + 1)r/δ ν,κ J 1 ((i + 1)r/δ ν,κ ) J 0 ((i + 1)r/δ ν,κ ) , ( A1
)
where J n is the n-th order Bessel function of the first kind and for rectangle channels of dimension 2a×2b (such is the case of the empty waveguide and for the stack square channels where a = b = r s )

f ν,κ = 8 π 2 2 m n F ν,κ m,n , (A2)
with m, n = 1, 3, 5, . . . and where

F ν,κ m,n = m 2 n 2 1 + i π 2 2 δ ν,κ r 2 b 2 m 2 + a 2 n 2 (a + b) 2 -1 , (A3)
with r = 2ab/(a + b) the characteristic length of the channel. The temperature dependence of the fluid (in this instance air) is accounted for through the variations of its viscosity and thermal conductivity following Ref. [38].

matrix of which is given by Eq. ( 7), Σ ± = ik ± = iωχ nr ± iω (χ a ) 2 + ρ/K = iωχ nr ± iωσ, and

V -1 = 1 √ 2   1/Kσ (σ -χ a )/σ -1/Kσ (σ + χ a )/σ   . (B2)
From the expression of Σ ± , the link with anisotropic fluid 39 becomes clear via the terms iωχ nr . Introducing Z ± = K(σ ± χ a ), left-multiplying Eq. (B1) by V -1 , and expressing the state vectors in terms of the scattering coefficients R + , T + , R -, and T -(see Fig. 1(b)), leads to the 2 following systems of equations

  1 Z - -1 Z +     T + T + /Z 0   =   e Σ + L Z -e Σ + L -e Σ -d Z + e Σ -L     1 + R + (1 -R + )/Z 0   ,   1 Z - -1 Z +     R -+ 1 (R --1)/Z 0   =   e Σ + L Z -e Σ + L -e Σ -L Z + e Σ -L     T - -T -/Z 0   , (B3)
where Z 0 is the impedance (divided by the duct cross-sectional area in our case) of the surrounding fluid. Introducing r ± = (Z 0 + Z ± )/(Z 0 -Z ± ), these equations can be inverted to yield

r + = -(R + R --T + T -+ 1) ± (R + R --T + T -+ 1) 2 -4R + R - 2R + , r -= -(R + R --T + T -+ 1) ± (R + R --T + T -+ 1) 2 -4R + R - 2R - , e 2iωσL = r -(R + r + + 1) R + + r + , e 2iωχ nr L = (R -+ r + ) (R -r -+ 1) (T -) 2 r + . ( B4
)
When T + = T -= T , this set of equations collapses to that derived in 30 for reciprocal system. From these equations, χ a exp , χ nr exp , K exp , and ρ exp are subsequently recovered. Special attention must be paid to the choice of the sign in the first two equations and to folding arising when L = N d, N ∈ N.

clearly visible around 60 Hz, 120 Hz and 180 Hz. Note that the transfer matrix is this time calculated via exp (A e L) when A e is used, because L is no longer small with respect to the wavelength. The real and imaginary parts of the dispersion relations are depicted in Figs.

4(m-n) respectively. Slight discrepancies are again noticed for each subfigure, notably for the normalized bulk modulus and for large ∆T , which are attributed to the assumption of the linear variation temperature distribution, but also to the very low frequency range of the measurements. Comments are similar to those already reported for the case of a single cell, notably with the presence of the coalescence point. Note that the amplification and attenuation are more pronounced and more visible on T + and T -this time and that the amplitude of the reflection coefficients are not small this time. This result validates the effective properties derived in Eqs. (11) and assumption of one-dimensional system. Real and imaginary parts of the dispersion relations (m-n). The results depicted with solid curves are calculated with the expressions given in Eq. ( 11), the dashed curves represent the numerical evaluation with the Tranfer Matrix method, and the circles depict the measurements. * jean-philippe.groby@univ-lemans.fr

  FIG. 1. (Color online) (a) Scaled representation of the unit cell together with the temperature distribution along the unit cell. The length of the experimental cell is d = 195 mm . (b) Photography of the elements of the thermoacoustic cell and (c) schematic view of the experimental setup for measuring the scattering matrix of the thermoacoustic amplifier.

  , both the even and non-reciprocal coupling parameters are almost purely imaginary and non-zero at low frequency for a thermoacoustic amplifier. The dispersion relation thus exhibits at first glance a band gap at zero frequency that is shifted along the imaginary wavenumber axis, i.e., k ± e = -iω|χ nr e | ± ω -|χ a e | 2 + ρ e C e , see Appendix B.

  FIG. 2. (Color online) Real and imaginary parts of the even Willis coupling parameter (a,b), normalized density (c,d), normalized bulk modulus (e,f), and non-reciprocal Willis coupling parameter (g,h) for ∆T = 0 K (blue curves), ∆T = 22 K (cyan curves), ∆T = 42 K (green curves), ∆T = 75 K (orange curves), ∆T = 100 K (red curves). Real and imaginary parts of the dispersion relations (i-j). The inset of the subfigure (i) shows a zoom on the real part of the dispersion relation at low frequency. The results depicted with solid curves are calculated with the expressions given in Eq. (11), the dashed curves represent the numerical evaluation with the Tranfer Matrix method, and the circles depict the measurements.

Figures 2

 2 Figures2(a-h) depict the real and imaginary parts of the even and non-reciprocal Willis coupling parameters, normalized density and bulk modulus as evaluated from the expressions

FIG. 3 .

 3 FIG. 3. (Color online) Transmission and reflection coefficients along the amplification, T + and R + (a,c), and attenuation, T -and R -(b,d), directions, for ∆T = 0 K (blue curves), ∆T = 22 K (cyan curves), ∆T = 42 K (green curves), ∆T = 75 K (orange curves), ∆T = 100 K (red curves).

FIG. 4 .

 4 FIG. 4. (Color online) Real and imaginary parts of the even Willis coupling parameter (a,b), normalized density (c,d), normalized bulk modulus (e,f), and non-reciprocal Willis coupling parameter (g,h) for ∆T = 0 K (blue curves), ∆T = 22 K (cyan curves), ∆T = 42 K (green curves), ∆T = 45 K (orange curves), ∆T = 100 K (red curves). Corresponding transmission and reflection coefficients along the amplification, T + and R + (i,k), and attenuation, T -and R -(j,l), directions.

1

  H. Nassar, X. Xu, A. Norris, and G. Huang, J. Mech. Phys. Solids 101, 10 (2017).2 B.-I. Popa and S. A. Cummer, Nat. Commun. 5, 3398 (2014).

  repetition of a unit cell of respective propagation matrix A e . The system is homogenized and A e does not depend on x. The state vectors at both sides of the unit cell W(d) and

	W(0) are related to each other through the 2 × 2 transfer matrix T = exp (A e d) of elements
	t ij , (i, j) ∈ [1, 2]. Following Ref. [30], the propagation matrix is correctly approximated by
	the inversion of the first-order Padé's approximation of the matrix exponential,
	A e ≈	2 d	(T + Id) -1 (T -Id)
	≈	2 d	1 Tr(T) + det(T) + 1	  t 11 -t 22 + det(T) -1 2t 21	 t 22 -t 11 + det(T) -1 2t 12  ,	(6)
	with Tr(T) = t 11 +t 22 . Equation (6) directly provides the elements of a non-reciprocal Willis
	material
			A e = iω	  χ a e + χ nr e C e	ρ e e + χ nr -χ a e	  ,

ACKNOWLEDGMENTS C.O., G.P. and G.P would like to acknowledge the support of the ANR SelfiXs project (ANR-18-CE92-0001). J.-P.G., M.M., and V.R.-G. would like to acknowledge the support of the ANR-RGC METARoom project (ANR-18-CE08-0021). J.L. would like to acknowledge the support of the Research Grants Council in Hong Kong (Grant No. 16302218). D.T. acknowledges financial support through the "Ramón y Cajal" fellowship under grant number RYC-2016-21188.

Appendix A: Expression of the effective properties in presence of viscothermal losses

The viscothermal losses are accounted for in the different elements composing the unit cell through the complex functions f ν and f κ . These functions describe the viscous and The two state vectors at both sides of the L-thick sample, whose propagation matrix is assumed to be independent of x by definition, are related by

where diag is the diagonal matrix, Σ ± are the eigenvalues of A e and V the corresponding eigenvector matrix. For an asymmetric and non-reciprocal Willis material, the constitutive Once the total transfer matrix T is calculated, it is directly assimilated to exp(A e d),

A e being homogenized. From Eq. (B1), it is clear that the eigenvectors of T and A e are identical and that the exponential of the eigenvalues of A e d are the eigenvalues of T. We immediately end up with

where Λ ± are the eigenvalues of T and V the associated eigenvector matrix. The four quantities χ nr num , χ a num , ρ num , and K num are subsequently evaluated.