Pierre Guiraud

Stefano Giordano

Olivier Bou Matar

Philippe Pernod

Raphael Lardat

Thermoacoustic wave generation in multilayered thermophones with cylindrical and spherical geometries

A thermoacoustic sound generation model, based on the classical balance equations of the continuum mechanics, is here developed for the cylindrical and the spherical thermoacoustic wave generation. In both geometries, the model considers an arbitrary multilayered structure, where each layer can be fluid or solid and it is characterized by the fully coupled thermo-visco-acoustic response. It means that the viscous behavior and the thermal conduction are considered in each layer. The model is based on a unified representation of cylindrical or spherical thermoacoustic waves, which is valid for both fluid and solid phases. Thanks to the continuity of temperature, particle velocity, normal stress and heat flux between adjacent layers, the model can be implemented by means of a versatile matrix approach, allowing flexible analysis and design of cylindrical or spherical thermophones. Any thermoacoustic variable can be determined at any position, any frequency and for any input power. The results are compared with the models already existing in the literature and the underlying physics is thoroughly discussed. The analysis is focused on the better understanding of the thermoacoustic generation with application to the state of the art of the thermophone technology.

I. INTRODUCTION

Traditional ways for sound generation involve the use of an electroacoustic transducer that is either composed of a piezoelectric material, or of a speaker consisting of a membrane, a coil and a magnet. Those electroacoustic transducers are widely adopted from everyday sound reproduction to scientific use for ultrasonic sound generation. Additionally, when we talk about the ultrasonic region of the spectrum we must mention other technologies used for electroacoustic transduction. For instance, the nondestructive evaluation and testing (NDE and NDT) of materials commonly involve the use of ultrasonic waves. The standard transduction techniques make use of the laser generation of ultrasound, [START_REF] Davies | Laser-generated ultrasound: its properties, mechanisms and multifarious applications[END_REF] or the acoustic wave generation by microwaves. [START_REF] Hosten | Acoustic wave generation by microwaves and applications to nondestructive evaluation[END_REF] In these cases, the acoustic field is typically generated in the form of pulses or bursts. An alternative and more recent method is based on capacitive micromachined ultrasonic transducers (CMUTs). [START_REF] Salim | Capacitive Micromachined Ultrasonic Transducers: Technology and Application[END_REF] This technique uses an electrostatic mechanism: the acoustic wave is generated utilizing the attraction force between the plates of a condenser. They are fabricated by means of the silicon technology, using standard processes of microelectronics, and exhibit a larger bandwidth than the piezoelectric devices together with a lower noise level.

However, dynamic loudspeakers, piezoelectric transducers and CMUTs are all based on mechanical vibrations and thus are inherently resonant. Today, no alternative for wideband sound generation is currently available. In the early XXth century, Arnold and Crandall started theorizing a phenomenon observed in the eighteen hundreds, where samples with low thermal capacity and high thermal conductiva) Electronic mail: stefano.giordano@univ-lille.fr (Author to whom correspondence should be addressed) AC generator ( ω 2 )

Heat wave (ω):

x FIG. 1. Schematic of the thermoacoustic sound generation through the thermophone principle. An electric current of frequency ω/2 is applied to the active layer, which is therefore subjected to an heating process at frequency ω (Joule effect). Heat is transmitted to the adjacent fluid, which is then subjected to a periodic sequence of compressions and dilations. This phenomenon leads to the generation of an acoustic wave at frequency ω that freely propagates outside the thermal layer.

ity produced sound when supplied with an alternative electrical current. [START_REF] Arnold | The Thermophone as a Precision Source of Sound[END_REF][START_REF] Weisendanger | The Thermophone[END_REF] This kind of samples are called thermophone. When supplied with an alternative electric power, the temperature of the sample will follow the electric oscillations (Joule effect) and the air in the vicinity (in the so-called thermal layer) will compress and dilate because of the exchanged heat, thus generating an acoustic pressure wave (see Fig. 1). Since no resonating part are involved in the process, the sound generation is wideband.

If the interest in this technology was largely overshadowed by the rise of traditional speakers, the modern possibilities to design and produce complex micro-and nano-materials allowed the fabrication of more efficient materials from the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458 thermoacoustic point of view. This point has rekindled the interest in thermophones theory and applications over the past 20 years. [START_REF] Shinoda | Thermally induced ultrasonic emission from porous silicon[END_REF] Most thermophones are carbon, [START_REF] Tian | Graphene-on-paper sound source devices[END_REF][START_REF] Fei | Low-voltage driven graphene foam thermoacoustic speaker[END_REF] or metal based. [START_REF] Niskanen | Suspended metal wire array as a thermoacoustic sound source[END_REF][START_REF] Tian | Flexible, ultrathin, and transparent sound-emitting devices using silver nanowires film[END_REF][START_REF] Dutta | Gold nanowire thermophones[END_REF] They can be in the geometry of a thin film, [START_REF] He | Transparent, flexible, ultrathin sound source devices using Indium Tin oxide films[END_REF][START_REF] Tsai | Fabrication of efficient thermoacoustic device with an interdigitated-like electrode on indium tin oxide glass[END_REF][START_REF] Tao | Flexible graphene sound device based on laser reduced graphene[END_REF] of forest of nanotubes or wires, [START_REF] Niskanen | Suspended metal wire array as a thermoacoustic sound source[END_REF][START_REF] Suzuki | Study of carbon-nanotube web thermoacoustic loud speakers[END_REF][START_REF] Aliev | Underwater sound generation using carbon nanotube projectors[END_REF] or in porous, foam like, 3D geometries. [START_REF] Aliev | Alternative nanostructures for thermophones[END_REF][START_REF] Lee | Analysis of acoustical performance of Bi-layer graphene and graphene-foam-based thermoacoustic sound generating devices[END_REF][START_REF] Guiraud | Two temperature model for thermoacoustic sound generation in thick porous thermophones[END_REF][START_REF] Ngoh | Experimental characterization of three-dimensional Graphene's thermoacoustic response and its theoretical modelling[END_REF] Alongside those experimental tryouts, many theoretical models have been developed and improved over the years. For instance, in 2008, Xiao et al. improved the Arnold and Crandall model by adding the heat stored in the thermophone using heat flow balance equation. [START_REF] Xiao | Flexible, stretchable, transparent carbon nanotube thin film loudspeakers[END_REF] Similarly, in 2013, Daschewski et al. added the influence of a substrate (or backing) to better describe the free standing thermophone response, and also considered the effect of the viscous dissipation in the propagation medium. [START_REF] Daschewski | Physics of thermo-acoustic sound generation[END_REF][START_REF] Daschewski | Influence of thermodynamic properties of a thermo-acoustic emitter on the efficiency of thermal airborne ultrasound generation[END_REF] Successively, in 2017, La Torraca et al. added conduction and convection contributions to the model, and also took the thickness of the substrate into account with a more refined technique. [START_REF] La Torraca | Physical modeling and characterization of thermo-acoustic loudspeakers made of silver nano-wire films[END_REF] These approaches are called piston based models. They are valid in many cases but may show limitations depending on the configuration under investigation (especially at high frequencies). Another method used for modeling the thermoacoustic effect is based on the classical conservation laws of continuum mechanics applied to the propagation media (typically without considering the viscosity). This system of equations was firstly elaborated by MacDonald et al. to study the photoacoustic effect in 1978. [START_REF] Mcdonald | Generalized theory of the photoacoustic effect[END_REF] More recently, the same set of equations has been adopted also for the thermoacoustics analysis. This coupled set of balance equations was solved in 2010 by Hu et al. for a thermophone placed onto a substrate, generating sound in a perfect gas. [START_REF] Hanping | Model for thermoacoustic emission from solids[END_REF] This model was validated against Shinoda experimental data, [START_REF] Shinoda | Thermally induced ultrasonic emission from porous silicon[END_REF] and was later adapted to evaluate the far field pressure and to consider the influence of the so-called Heat Capacity per Unit Area (HCPUA). [START_REF] Tian | Graphene-on-paper sound source devices[END_REF] This parameter is defined as the product of three factors: the density of the thermophone generating layer (kg/m 3 ), the specific heat of the same layer (J/(kg•K)), and the generating layer thickness (m). The HCPUA must be low enough to have an efficient thermoacoustic radiation. The refined methodology of Hu et al. has been used in 2012 to study the influence of the main thermophone parameters on its wideband frequency response. [START_REF] Hu | Wideband flat frequency response of thermo-acoustic emission[END_REF] Successively, in 2012, Hu et al. modified this approach to investigate the thermophones with spherical geometry (or acoustic monopoles). [START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF] This approach was generalized in 2014 to arbitrary sources using an arbitrary distribution of acoustic monopole, [START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF] and in 2017 to represent a cylindrical geometry. [START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] Lastly, in 2019, Mao et al. added the influence of viscosity to the Hu et al. equations. [START_REF] Mao | Thermo-acoustic radiation of free-standing nano-thin film in viscous fluids[END_REF] Furthermore, based on previous investigations, [START_REF] Lim | Theory of suspended carbon nanotube thinfilm as a thermal-acoustic source[END_REF][START_REF] Tong | Gap separation effect on thermoacoustic wave generation by heated suspended CNT nano-thinfilm[END_REF] a generalized theory was developed by Liu et al. in 2018 to describe point source, line source and line array thin film thermophones in free field and in half space, also with an air gap between the substrate and the generating layer. [START_REF] Liu | Thermo-acoustics generated by periodically heated thin line array[END_REF] Due to the rather low efficiency of thermophones, it is of prime interest to optimize the thermophone acoustic radiation, especially for loudspeaker applications. One direction considered to perform this optimization makes use of the flexibility of the fabrication processes. Indeed, modern methods can be exploited to design and produce optimal shapes of thermoacoutic generators. However, despite the exploration of several geometries as thermophone samples (films, nanotube forests, sponges, and so on), they are often represented by thin sheets in the modeling activity. In fact, most of the developed theoretical models are based on planar generators and plane waves propagation. Only some investigations, as briefly discussed above, generalized the plane wave model to spherical or cylindrical symmetries. However, due to the application oriented purpose of those models, simplifications have been made and only specific configurations have been studied. For instance, concerning the spherical model, the objective was more about using acoustic superposition theory to explore complex thermophone geometries than about understanding all the possibilities offered by the spherical symmetry itself (for example with multilayered systems). Furthermore, a refined investigation of the near field generation of the cylindrical geometry could improve the understanding of the thermoacoustic mechanisms in single wires, in nanotubes or even in the tubular branches of porous or foam like thermophones. If the plane wave analysis is useful in many cases at a macroscopic scale, the geometry is typically more complex at a microscopic scale and the analysis of cylindrical or spherical waves needs to be considered to better understand the near field generation.

In this work a multilayer model for cylindrical and spherical wave generation is elaborated. It is based on the conservation equations of the continuum mechanics and takes into account the fully coupled thermo-visco-acoustic behavior of each fluid or solid layer. This methodology was recently proposed for plane wave generation, [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] and is here generalized to the cylindrical and the spherical geometries. With respect to previous approaches, the following improvements can be mentioned: (i) a versatile matrix approach, based on the continuity of the main variables across the interfaces between adjacent layers, allows the consideration of a completely arbitrary multilayered structure composed of both fluid and solid layers; (ii) the main dissipative phenomena, namely the thermal conduction and the viscous behavior, are fully considered in each layer of the system; (iii) since the thermo-visco-acoustic wave propagation is introduced in each layer, the acoustic wave propagation in solids is then considered, and allows the study of possible resonances and antiresonances of the structure. The approach presented here will be compared with some previously existing models and a complete analysis of the effects of the main important physical and geometrical parameters is performed.

In Sections II and III, the physical and mathematical formalism adopted for studying the cylindrical and spherical geometries is introduced for both fluid and solid phases. Moreover, the matrix methodology to develop the model for a multilayered structure is discussed. Then, in Section IV, this technique is applied to specific thermophone configurations. In particular, in Section IV A, some results are compared with the model from the literature, and, in Section IV B, different thermophone structures (full thermophone, hollow thermophone and hollow thermophone with substrate) are investigated in more detail and their defining parameters will be modified to better assess their influence. In particular, both This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

frequency and position dependent analyses will be performed.

II. MULTILAYER MODEL FOR CYLINDRICAL WAVES

It can been proven that the linearized conservation equations in any fluid with viscosity and thermal conduction can be written as [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF][START_REF] Landau | Fluid Mechanics[END_REF]

     ρ 0, f B f ∂ p ∂t -ρ 0, f α T ∂ T ∂t + ρ 0, f ∇ • v = 0, ρ 0, f ∂ v ∂t = -∇p + µ∇ 2 v + (λ + µ) ∇( ∇ • v), ρ 0, f C p ∂ T ∂t -α T T 0 ∂ p ∂t = κ∇ 2 T. (1)
where the pressure p (Pa), the temperature variation T (K) and the particle velocity vector v (m/s) are the main variables depending on time t (s) and space r (m). Moreover, ρ 0, f is the density (kg/m 3 ), B f the bulk modulus (Pa) , α T the coefficient of volumetric expansion (1/K), λ and µ the first and second viscosity coefficients (Pa•s), C p the specific heat at constant pressure (J/(kg•K)), T 0 the ambient temperature (K) and, finally, κ the thermal conductivity (W/(m•K)). The subscript f in most parameters means fluid. All these parameters will be considered as constants in each cylindrical layer, and they can only vary from one layer to another. It is worth noticing that in some practical applications it could be useful to consider the dependence of some parameters on the temperature. However, the introduction of this dependence transforms the system from linear to nonlinear, eliminating the possibility of doing the mathematical analysis developed in this article. Only more time-consuming numerical methods (finite elements or similar techniques) can consider such kind of dependence.

The proposed system of differential equations describes the heat transfer and the mechanical wave propagation while taking into account dissipative phenomena introduced through the viscous behavior of the fluid. Equation ( 1) is obtained by a combination of the linearized classical conservation laws (of mass, linear momentum and energy) of continuum mechanics and the linearized constitutive equations for a fluid [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF][START_REF] Landau | Fluid Mechanics[END_REF] . This approach is justified by the fact that we are here considering small variations of the relevant quantities p, v, T around given equilibrium values.

Under the hypotheses of small deformations and small temperature variations around equilibrium, a similar system of equations can be stated for an isotropic and homogeneous solid as [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF][START_REF] Landau | Theory of Elasticity[END_REF]

     ρ 0,s C V,s ∂ T ∂t = κ s ∇ 2 T -α T,s B 0 ∂ ∂t ∇ • uT 0 + Q 0 , ρ 0,s ∂ 2 u ∂t 2 = (λ 0 + µ 0 ) ∇( ∇ • u) + µ 0 ∇ 2 u +(λ s + µ s ) ∇( ∇ • v) + µ s ∇ 2 v -α T B 0 ∇T + b, (2)
with main variables given by the particle displacement vector u (m), the particle velocity vector v = ∂ u/∂t (m/s), and the temperature variation T (K). We also introduced the Lamé elastic coefficients λ 0 and µ 0 (Pa), the bulk modulus B 0 = λ 0 + (2/3)µ 0 (Pa), the coefficient of volumetric expansion α T,s (1/K), the first and second viscosity coefficients λ s and µ s (Pa•s), the specific heat at constant volume C V,s (J/kg•K), the externally applied body forces b (N) and the supplied thermal power density Q 0 (W/m 3 ). The subscript s in most parameters means solid. In this case, dealing with a solid material, we considered only the conservation of energy and linear momentum since the mass density is considered constant within the linearized formalism.

Those conservation equations, for both fluids and solids, have been solved for a one-dimensional model, assuming plane wave generation and propagation, in recent literature. [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] Considering a sytem with cylindrical symmetry, we shall now solve those equations assuming a cylindrical wave generation, thus generalizing previous results. [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] Explicit analytical solutions for the variables T , p, v and q, the heat flux, are obtained in a single material, and then a transfer matrix method is introduced to investigate the thermoacoustic behavior of systems with an arbitrary distribution of fluid and solid cylindrical layers. This will allow to create a multilayer model, which we will later use to investigate near field sound generation in real life microscopic thermophone structures.

The cylindrical thermophone is supposed to be infinitely long in the z direction. Moreover, all physical properties are supposed to be independent of the angle θ on the xy plane due to the cylindrical symmetry. The main variables describing the thermoacoustic wave generation are therefore solely dependent on the radius r. We analyze the effects of these assumptions on the main equations for both fluid and solid layers in the next sections.

A. Thermoacoustic cylindrical waves in a fluid

Using the general equations for a fluid given in Eq.(1), a cylindrical geometry is now investigated as shown in Fig. 2. The pressure, the temperature variation and the particle velocity can be written as

     p = p(r), T = T (r), v = r r v(r), (3)
with r = x e x + y e y and || r|| = r = x 2 + y 2 (see Fig. 2 for the definition of coordinates). In order to write Eq.(1) in cylindrical coordinates, each component of those balance equations has to be rewritten using the prescriptions of Eq.( 3). This leads to

     iω 1 B f p -iωα T T + 1 r ∂ ∂ r (rv) = 0, iωρ 0, f v = -∂ p ∂ r + (λ + 2µ)(-1 r 2 v + 1 r ∂ v ∂ r + ∂ 2 v ∂ r 2 ), iωρ 0, f C p T -iωα T T 0 p = κ( 1 r ∂ T ∂ r + ∂ 2 T ∂ r 2 ), (4)
with 1 r ∂ ∂ r (rv) = 1 r (v + r ∂ v ∂ r ).
Here, we also substituted the time derivative ∂ /∂t with iω by assuming a monochromatic sinusoidal permanent regime with angular frequency ω.

The aim is to write the pressure p and the particle velocity v as a function of the temperature T and find the differential equation that T solves. To this aim Eq.(4a) can be rewritten as

p = α T B f T - B f iω ( 1 r v + ∂ v ∂ r ).
(5) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Here, the hydrostatic cylindrical pressure p does not take the viscosity into account and so, in order to have a continuity relation between different adjacent media (solid, fluid), we need to define p, which represents the normal component of the total stress tensor T = -p Î + TV in cylindrical coordinates (with negative sign). Knowing the viscous stress TV being defined as

PLEASE CITE THIS ARTICLE AS

TV = V ∂ ε ∂t , (6)
with V the viscous tensor

V Ŝ i j = V i jkh S kh = 2µS i j + λ S kk δ i j , (7)
p is written in the form

p = -T n • n = p Î n • n -TV n • n = p -2µ ∂ v ∂ r -λ 1 r ∂ ∂ r (rv) = α T B f T - B f iω + λ 1 r ∂ ∂ r (rv) -2µ ∂ v ∂ r , (8)
where n = r/r. The variable p here introduced will be useful to develop the transfer matrix method for multilayered structures. Anyway, by deriving Eq.( 5), and injecting this result in Eq.(4b), we have

iωρ 0, f v = -α T B f ∂ T ∂ r + (λ + 2µ + B f iω ) ∂ ∂ r ( 1 r ∂ ∂ r (rv)). (9)
Similarly, by injecting Eq.( 5) in Eq.(4c) and using the thermodynamic relation 35

ρ 0, f (C p -C V ) = T 0 α 2 T B f , (10)
we have after straightforward elaborations

1 r ∂ ∂ r (rv) = κ α T T 0 B f ( 1 r ∂ T ∂ r + ∂ 2 T ∂ r 2 ) - iωρ 0, f C V α T B f T 0 T. ( 11
)
Injecting now Eq.( 11) into Eq.( 9) allows us to write v as a function of T only

v = - 1 iωρ 0, f α T B f + (λ + 2µ + B f iω ) iωρ 0, f C V α T B f T 0 ∂ T ∂ r + (λ + 2µ + B f iω ) κ iωρ 0, f α T T 0 B f ∂ ∂ r ( 1 r ∂ T ∂ r + ∂ 2 T ∂ r 2 ). ( 12
)
To obtain a pure equation for T we apply the operator

1 r ∂ ∂ r [r(•)]
to Eq.( 9) and we combine the result with Eq.( 11). Hence, by introducing the notations

T ii = DT = 1 r ∂ ∂ r (r ∂ T ∂ r ), T iv = (T ii ) ii = DDT = 1 r ∂ ∂ r r ∂ ∂ r [ 1 r ∂ ∂ r (r ∂ T ∂ r )] , (13)
a fourth order differential equation for T can be obtained

0 = (λ + 2µ + B f iω )κT iv + -(λ + 2µ + B f iω )iωρ 0, f C V + α 2 T B 2 f T 0 + iωρ 0, f κ T ii -ω 2 ρ 2 0, f C V T, (14)
and it is similar to the equation for plane waves but with different definition of the derivatives of T . [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] Following exactly the same procedure previously introduced for plane waves and using the notations 35

             γ = C p C V , C 0 = B f ρ 0, f γ, [m/s] l κ = C 0 κ B f C p , [m] l V = λ +2µ ρ 0, f C 0 , [m] (15)
Eq.( 14) can be finally written as

l κ (l V -i C 0 ωγ )T iv -[1 + iω C 0 (l κ + l V )]T ii - ω 2 C 2 0 T = 0. ( 16
)
The four solutions of the algebraic equation associated to Eq.( 16) are written in the form ±θ ac, f and ±θ th, f as

   θ ac, f = iω C 0 1 -1 2 iω C 0 l V -1 2 l κ iω C 0 (1 -1 γ ) , θ th, f = iωγ C 0 l κ 1 + 1 2 iω C 0 l κ (1 -1 γ ) + l V (1 -γ) . (17)
They are associated to progressive and regressive acoustical waves and thermal waves, respectively, propagating in the fluid phase. The details of their calculation is found in Appendix B 1. The general form of the solution for the fourth order differential equation in Eq.( 16) is obtained in AppendixA.

It can be written as

T (r) =AH
0 (iθ th, f r), (18)
with A, B,C, D integration constants and H

i (z) and H

(2)

i (z) being the cylindrical Hankel function of order i and of first and second kind. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

B. Thermoacoustic cylindrical waves in a solid

The cylindrical geometry will now be investigated for a thermo-visco-elastic solid medium. Based on Eq.( 2) and using Eq.( 3), the main equations in a solid become

ρ 0,s C V,s iωT = κ s 1 r ∂ ∂ r (r ∂ T ∂ r ) -α T,s B 0 T 0 1 r ∂ ∂ r (rv) + Q 0 , ρ 0,s iωv = λ 0 +2µ 0 iω + λ s + 2µ s ∂ ∂ r 1 r ∂ ∂ r (rv) -α T,s B 0 ∂ T ∂ r , (19) with
v = 1 iωρ 0,s 1 α T,s B 0 T 0 λ 0 + 2µ 0 iω + λ s + 2µ s × κ s ∂ ∂ r ( 1 r ∂ ∂ r (r ∂ T ∂ r )) -ρ 0,s C V,s iω ∂ T ∂ r -α T,s B 0 ∂ T ∂ r . (20)
Then, by applying to Eq.( 20) the operator 1

r ∂ ∂ r [r(•)],
and using again Eq.(19a), a differential equation for T is eventually found as

ρ 0,s iωQ 0 =(λ s + 2µ s + λ 0 + 2µ 0 iω )κ s T iv -(λ s + 2µ s + λ 0 + 2µ 0 iω )iωρ 0,s C V,s + α 2 T,s B 2 0 T 0 + iωρ 0,s κ s T ii -ω 2 ρ 2 0,s C V,s T, (21)
using the notation for derivatives introduced in Eq. (13). Apart from the supplied input power Q 0 , Eq.( 21) is similar to Eq.( 14) for fluid with only B f changed to λ 0 + 2µ 0 for the solid. Henceforth, even if the solutions θ of the associated algebraic equation will be numerically different between fluid and solid, the whole process to find T in the fluid, resulting in Eq.( 18) (see Appendix A), is still valid for a solid. In order to add the effect of the supplied energy to the final solution for T , we use the assumption that Q 0 is spatially uniform, and we easily find the particular solution,

T Q 0 = Q 0 ρ 0,s C v,s iω , ( 22
)
which can be used in the following general solution

T (r) =AH (1) 0 (iθ ac,s r) + BH (2)
0 (iθ ac,s r)+ CH

0 (iθ th,s r) + T Q 0 . ( 23
)
We start now the investigation of the pressure behavior in the solid case. By definition, a force d f applied on a surface dS within a solid can be written using the stress tensor T as d f = T ndS, where n is the unit vector perpendicular to dS. [START_REF] Landau | Theory of Elasticity[END_REF] Therefore, the elastic pressure p in the direction n can be calculated as

p = -Te n • n = -T e i j X i X j r 2 , (24)
with n = r r and X i , i = 1, 2, 3, are the coordinates of r in cartesian coordinates. Here, Te is the elastic part of the stress tensor. It is important to note that a plane strain hypothesis has been used in this analysis. This means that the stress is acting perpendicularly to the length of the cylinder leading to X 3 = 0 and ∂ ∂ X 3 = 0. With ε the infinitesimal deformation tensor de-

fined as 37 ε = 1 2 ( ∇ u + ∇ u T ), (25)
the elastic stress tensor elements are written 37

T e i j = 2µ 0 ε i j + λ 0 ε kk δ i j -α T (λ 0 + 2 3 µ 0 )δ i j T. ( 26
)
with the displacement vector u defined as u = r r u(r). The three elastic tensor elements necessary for the pressure calculation are the following

T e 11 = 2µ 0 ∂ u 1 ∂ X 1 + λ 0 ε kk -α T,s B 0 T, (27)
T e 12 = µ 0 ∂ u 1 ∂ X 2 + ∂ u 2 ∂ X 1 , (28)
T e 22 = 2µ 0 ∂ u 2 ∂ X 2 + λ 0 ε kk -α T,s B 0 T. (29)
The partial derivatives of u i are

∂ u i ∂ X j = X i X j r 2 ∂ u ∂ r + r 2 δ i j -X i X j r 3 u, (30)
and therefore we obtain the pressure after straightforward calculations

p = -T e 11 X 2 1 r 2 -2T e 12 X 1 X 2 r 2 -T e 22 X 2 2 r 2 = α T,s B 0 T - 2µ 0 iω ∂ v ∂ r - λ 0 iω 1 r ∂ ∂ r (rv). (31)
Here, the pressure p is only defined with the elastic stress tensor and so, in order to have a stress continuity between different media (solid, fluid), we need to define p as (using Eqs.( 6) and ( 7)

) p = -T n • n = -Te n • n -TV n • n = p -2µ s ∂ v ∂ r -λ s 1 r ∂ ∂ r (rv) = α T,s B 0 T -2( µ 0 iω + µ s ) ∂ v ∂ r + ( λ 0 iω + λ s ) 1 r ∂ ∂ r (rv) , (32)
where the total stress tensor composed of elastic and viscous parts is considered. The definition of p within the solid layer is crucial to impose the stress continuity at the interface in multilayered systems, as discussed in the next Section. Due to the similarity of the differential equation for T in a fluid and in a solid, the temperature variation has been here This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. written as in the fluid in Eq.( 23). The differences are found in the wavenumbers θ ac and θ th and in the presence of the term depending on the local heating. In order to explicit the differences in wavenumbers between solid and liquid media, a first order solution of the algebraic equation associated to Eq.( 21) is determined (details in Appendix B 2) as follows

PLEASE CITE THIS ARTICLE AS

                         θ ac,s = iω C 0,s (1 -1 γ s + C L C 0,s ) -1/2 × 1 -l V,s iω 2C 0,s (1 -1 γ s + C L C 0,s ) -1 -l κ,s iω 2C 0,s × (1 -1 γ s + C L C 0,s ) -1 -C L C 0,s (1 -1 γ s + C L C 0,s ) -2 , θ th,s = ( iω l κ,s C L ) 1/2 (1 -1 γ s + C L C 0,s ) 1/2 1 + iω 2C 0,s × l κ,s (1 -1 γ s + C L C 0,s ) -1 -(1 -1 γ s + C L C 0,s ) -2 C L C 0,s + l V,s (1 -1 γ s + C L C 0,s ) -1 - C 0,s C L . ( 33
)
We note that Eq.( 33) is consistent with Eq.( 17) if C L is replaced by C 0,s γ s . This similitude is coming from the comparison of the coefficient of the second degree in both cases.

C. Transfer matrix method with N cylindrical layers

Having solved the equations for cylindrical thermoacoustic radiations in both fluid and solid phases, we will now investigate the behavior of a multilayered cylindrical thermophone. As seen in Fig. 3, a system with N layers defined by concentric cylinders is considered. The external layer, at the extremity (layer N), is supposed to be semi infinite, and each layer j ∈ [1, N -1] potentially has a volumetric input source Q 0, j (Q 0,N = 0). The thickness of each circular crown j is r jr j-1 with r j the external radius of the j-th concentric cylinder. The type of material of each layer (fluid or solid) is not explicitly specified since the material differences are taken into account in the inherent parameters of their defining matrix, as discussed below.

To build the multilayered model, we first define the vector heat flux in each circular crown as

q = -κ ∇T = -κ r r ∂ T ∂ r , (34)
and his scalar version q = -κ ∂ T ∂ r . We also define T Q 0 as in Eq.( 22) and the functions F and G as

F (η) = α T B f - B f iω + λ + 2µ (L 1 η 2 + L 2 η 4 ), (35)
G (η) = L 1 + L 2 η 2 . ( 36
)
Here, the coefficients L 1 and L 2 are defined through Eq.( 12) for a fluid and Eq.( 20) for a solid, when these expressions are written in the form v = L 1

∂ T ∂ r + L 2 ∂ ∂ r ( 1 r ∂ ∂ r (r ∂ T ∂ r )).
We also remark that when Eq.( 35) is written for a solid layer, one must substitute B f with λ 0 + 2µ 0 and all parameters must correspond to those of the solid material. In order to rewrite the final expressions of velocity, pressure and heat flux in terms of Hankel functions, we use the following derivative formulas and recurrence properties 38

C ν (z) = -C ν+1 (z) + ν z C ν (z), (37)
C ν (z) = C ν-1 (z) - ν z C ν (z), (38)
C ν-1 (z) + C ν+1 (z) = 2ν z C ν (z), (39)
with C ν any cylindrical function of order ν. The general form of the parameters p, v, q, T in an arbitrary layer can be found to be

                                                 p(r) = A F (θ ac )H (1) 0 (iθ ac r) -2µi θ ac r G (θ ac )H (1) 1 (iθ ac r) +B F (θ ac )H (2) 0 (iθ ac r) -2µi θ ac r G (θ ac )H (2) 1 (iθ ac r) +C F (θ th )H (1) 0 (iθ th r) -2µi θ th r G (θ th )H (1) 1 (iθ th r) +D F (θ th )H (2) 0 (iθ th r) -2µi θ th r G (θ th )H (2) 1 (iθ th r) +α T B f T Q 0 , v(r) = -Aiθ ac G (θ ac )H (1) 1 (iθ ac r) -Biθ ac G (θ ac )H (2) 1 (iθ ac r) -Ciθ th G (θ th )H (1) 1 (iθ th r) -Diθ th G (θ th )H (2) 1 (iθ th, f r), q(r) = Aiθ ac κH (1) 1 (iθ ac r) + Biθ ac κH (2) 1 (iθ ac r) +Ciθ th κH (1) 1 (iθ th r) + Diθ th κH (2) 1 (iθ th r), T (r) = AH (1) 0 (iθ ac r) + BH (2) 0 (iθ ac r) +CH (1) 0 (iθ th r) + DH (2) 0 (iθ th r) + T Q 0 .
(40) These equations apply for both a fluid and a solid layer provided that the appropriate changes in parameters are applied. Additionally, when we deal with a solid layer we have to replace B f by λ 0 + 2µ 0 in F and G , as previously discussed. Moreover, for the solid layers, in the expression for p we have to substitute µ with µ + µ 0 iω as follows p(r) =A F (θ ac )H

(1)

0 (iθ ac r) -2(µ + µ 0 iω )i θ ac r G (θ ac )H (1)
1 (iθ ac r) + ...

We are able to rewrite Eq.( 40) in matrix form for the layer j This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

as follows    p v q T    j = W j    A j B j C j D j    +    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j    , (42)
with

T Q 0 ,s = Q 0, j
ρ 0, j C V, j iω . The general form of the matrix W j is then written as

W j =            F (θ ac )H (1) 0 (iθ ac r)+ F (θ ac )H (2) 0 (iθ ac r)+ -2iµ θ ac r G (θ ac )H (1)
1 (iθ ac r) -2iµ θ ac r G (θ ac )H

(2)

1 (iθ ac r) -iθ ac G (θ ac )H (1) 1 (iθ ac r) -iθ ac G (θ ac )H (2)
1 (iθ ac r)

iθ ac κH (1)
1 (iθ ac r) iθ ac κH

(2)

1 (iθ ac r) H (1)
0 (iθ ac r) H

(2) 0 (iθ ac r)

F (θ th )H (1) 0 (iθ th r)+ F (θ th )H (2) 0 (iθ th r)+ -2iµ θ th r G (θ th )H (1) 1 (iθ th r) -2iµ θ th r G (θ th )H (2) 1 (iθ th r) -iθ th G (θ th )H (1) 1 (iθ th r) -iθ th G (θ th )H (2)
1 (iθ th r)

iθ th κH (1)
1 (iθ th r) iθ th κH

(2)

1 (iθ th r) H (1)
0 (iθ th r) H

(2)

0 (iθ th r)            , (43)
where, of course, all parameters correspond to the j-th layer.

Here, as before, when we deal with a solid layer, we must substitute B f with λ 0 + 2µ 0 and we have to consider the correct values for θ ac and θ th .

As shown in Fig. 3, r = 0 corresponds to the center of the cylinder and due to the cylindrical symmetry only positive values of r are considered. In order to avoid any singularity at r = 0, we define in the first layer the conditions

A 1 = B 1 and C 1 = D 1 .
This assumption is based on the property H

(1)

1 (z) ∼ -H (2)
1 (z) when z → 0. [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] The last layer N is considered to be semi infinite so only progressive waves are assumed and we therefore impose A N = C N = 0 (absence of regressive waves). In addition to those boundary conditions, the continuity of the parameters p, v, q and T is assumed between two adjacent layers for obvious physical reasons. Of course, it means that we do not take into account imperfect interfaces. [START_REF] Pavanello | Generalized interface models for transport phenomena: Unusual scale effects in composite nanomaterials[END_REF][START_REF] Pavanello | How imperfect interfaces affect the nonlinear transport properties in composite nanomaterials[END_REF] This continuity is automatically verified if we adopt the matrix form given in Eq. (42). By imposing the continuity of the variables at interfaces between layers we found

   A 1 A 1 C 1 C 1    = M 0    0 B N 0 D N    + ∑ n∈[1,N-1] M Q 0,n    α T,n B f ,n T Q 0 ,n 0 0 T Q 0 ,n    , ( 44
)
where n is summed over layers with a volumetric source and

M 0 = W -1 1 (r 1 ) N-1 ∏ j=2 W j (r j-1 ) W -1 j (r j ) W N (r N-1 ), (45)
with

                 M Q 0,1 = -W -1 1 (r 1 ), M Q 0,n = W -1 1 (r 1 ) ∏ n-1 j=2 W j (r j-1 ) W -1 j (r j ) × 1 -W n (r n-1 )W -1 n (r n ) , ∀n ∈ [2, N -1], M Q 0,N = W -1 1 (r 1 ) ∏ N-1 j=2 W j (r j-1 )W -1 j (r j ) . (46)
Eq.( 44) can be rewritten in a more explicit form as

   1 0 -M 0 (1, 2) -M 0 (1, 4) 1 0 -M 0 (2, 2) -M 0 (2, 4) 0 1 -M 0 (3, 2) -M 0 (3, 4) 0 1 -M 0 (4, 2) -M 0 (4, 4)       A 1 C 1 B N D N    =                   ∑ n M Q 0,n (1, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (1, 4)T Q 0 ,n ∑ n M Q 0,n (2, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (2, 4)T Q 0 ,n ∑ n M Q 0,n (3, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (3, 4)T Q 0 ,n ∑ n M Q 0,n (4, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (4, 4)T Q 0 ,n                   . ( 47
)
It is important to remark that Eq.( 44) or Eq.(47) represents a system of linear equations giving A 1 , C 1 , B N , and D N . The determination of these four parameters allows the calculation of any coefficient (ABCD) j using either

   A j B j C j D j    = W -1 j (r j )   Wj+1(r j )    A j+1 B j+1 C j+1 D j+1    +    α T, j+1 B f , j+1 T Q 0 , j+1 0 0 T Q 0 , j+1    -    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j       , (48) with j ∈ [1, N -1], or    A j B j C j D j    = W -1 j (r j-1 )   Wj-1(r j-1 )    A j-1 B j-1 C j-1 D j-1    +    α T, j-1 B f , j-1 T Q 0 , j-1 0 0 T Q 0 , j-1    -    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j       , (49)
with j ∈ [2, N]. In Eq.(48), we calculate recursively the coefficients of the layers moving from the last external layer to the center of the structure. Alternatively, in Eq.(49), we calculate recursively the coefficients of the layers going from the central This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

layer to the external one of the structure. The two procedures are theoretically equivalent and, depending on the system under consideration, it is possible to choose the better one from the numerical point of view. Finally p, v, q, T can be found for any value of r by means of Eq.( 42). This methodology allows for studying the thermoacoustic wave generation and propagation in arbitrarily complex systems with cylindrical geometry.

III. MULTILAYER MODEL FOR SPHERICAL WAVES

In section II the balance equations for fluids, Eq.( 1), and the balance equations for solids, Eq.( 2), were solved for thermoacoustic cylindrical waves. Using the same procedure and starting equations, the problem will be now solved under the assumption of spherical waves.

The system of equations will be solved first in a fluid and then in a solid medium using a spherical coordinate system with the hypothesis that all physical variables depend on the radius r while they are independent of azimuth and elevation angles. The solutions obtained describe therefore thermoacoustic spherical waves in a medium with thermal conduction and viscous dissipation. As before, the general solutions obtained are used to develop a multilayered spherical model able to represent spherical thermophones of technological interest.

Due to the similarities between the Sections II and III, it is advised for a first read to skip this section if the analysis of results is of more interest to the reader.

A. Thermoacoustic spherical waves in a fluid

Using the general equations for a fluid given in Eq.( 1), thermoacoustic spherical waves are now investigated in systems as shown in Fig. 4 with spherical symmetry. The main variables appearing in Eq.( 1) now can be written as

     p = p(r), T = T (r), v = r r v(r), (50)
with r = x e x + y e y + z e z , || r|| = r = x 2 + y 2 + z 2 and x, y, z, the adopted coordinates (see Fig. 4). In order to write Eq.( 1) in spherical coordinates, with only the radial abscissa r, each variable of those equations has to be rewritten using the assumptions in Eq.( 50). Hence, we have that Eq.( 1) becomes

       iω 1 B f p -iωα T T + 1 r ∂ ∂ r (rv) = 0, iωρ 0, f v = -∂ p ∂ r + (λ + 2µ) ∂ ∂ r 1 r 2 ∂ ∂ r (r 2 v) , iωρ 0, f C p T -iωα T T 0 p = κ 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ). ( 51
)
The aim is to write the pressure p and the particle velocity v as a function of the temperature variation T and find the differential equation that T solves, similarly as before. To this aim, Eq.(51a) can be rewritten in the same form as Eq.( 5) As previously discussed, the pressure p does not take the viscosity into account and so, in order to have a continuity relation between different adjacent media (solid, fluid), we need to define p, which represent the normal component of the total stress tensor T = -p Î + TV in spherical coordinates (with negative sign). We obtain

p = α T B f T - B f iω 1 r 2 ∂ ∂ r (r 2 v). (52)
p = -T n • n = p Î n • n -TV n. n = p -2µ ∂ v ∂ r -λ 1 r 2 ∂ ∂ r (r 2 v) = α T B f T - B f iω + λ 1 r 2 ∂ ∂ r (r 2 v) -2µ ∂ v ∂ r . (53)
This function must be continuous across the interfaces between adjacent materials. This point will be exploited later on. By deriving Eq.( 52) with respect to r, and injecting the result in Eq.(51b), we obtain

iωρ 0, f v = -α T B f ∂ T ∂ r + (λ + 2µ + B f iω ) ∂ ∂ r 1 r 2 ∂ ∂ r (r 2 v) .
(54) In the same way, by substituting Eq.(52) in Eq.(51c), we have

1 r ∂ ∂ r (rv) = κ α T T 0 B f 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ) - iωρ 0, f C V α T B f T 0 T, (55)
having used the thermodynamic relation stated in Eq. (10). Finally, by injecting Eq.(55) into Eq.( 54), we can write v as a function of

T only v = - 1 iωρ 0, f α T B f + (λ + 2µ + B f iω ) iωρ 0, f C V α T B f T 0 ∂ T ∂ r + (λ + 2µ + B f iω ) κ iωρ 0, f α T T 0 B f ∂ ∂ r 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ) . (56)
To obtain a pure equation for the temperature T , we apply now the operator 1

r 2 ∂ ∂ r [r 2 (•)]
to Eq.( 54) and then we combine the result with Eq.(55). So, we can redefine (with respect to the cylindrical geometry) the notations for the differential operators as follows

T ii = DT = 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ), T iv = (T ii ) ii = DDT = 1 r 2 ∂ ∂ r r 2 ∂ ∂ r [ 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r )] ,
(57) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

and the fourth order differential equation for T can be written as

0 =(λ + 2µ + B f iω )κT iv -(λ + 2µ + B f iω )iωρ 0, f C V + α 2 T B 2 f T 0 + iωρ 0, f κ T ii -ω 2 ρ 2 0, f C V T, (58)
which is similar to Eq.( 14) but, importantly, with a different definition of the differential operators. Adopting the same notations introduced studying the thermoacoustic cylindrical waves in a fluid, Eq.( 58) can be rewritten as

l κ (l V + C 0 iωγ )T iv -[1 + iω C 0 (l κ + l V )]T ii - ω 2 C 2 0 T = 0, ( 59
)
with the parameters defined in Eq.( 15) and the differential operators in Eq.( 57). The four solutions of the associated algebraic equation are written of the form ±θ ac, f and ±θ th, f . They assume exactly the same values defined as in Eq.( 17). Moreover, it can be easily shown that if T is written in the form T (r) = e ±θ i r r then

DT = θ 2 i T, DDT = θ 4 i T. (60)
where the operator D is defined in Eq.( 57) and θ i being a solution of the associated fourth degree algebraic equation. Then, T (r) = e ±θ i r r is a solution for Eq.( 59). Hence, the final general solution of Eq.( 58) is eventually found as

T (r) = A e -θ ac, f r r + B e θ ac, f r r +C e -θ th, f r r + D e θ th, f r r , (61)
with A, B,C, D integration constants.

B. Thermoacoustic spherical waves in a solid

The spherical geometry will now be investigated in a solid phase. Based on Eqs.( 2) and (50), the balance equations in a solid become

ρ 0,s C V,s iωT = κ s 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ) -α T,s B 0 T 0 1 r 2 ∂ ∂ r (r 2 v) + Q 0 , ρ 0,s iωv = λ 0 +2µ 0 iω + λ s + 2µ s ∂ ∂ r 1 r 2 ∂ ∂ r (r 2 v) -α T,s B 0 ∂ T ∂ r , (
v = - 1 iωρ 0,s C V,s iωρ 0,s α T,s B 0 T 0 λ 0 + 2µ 0 iω + λ s + 2µ s +α T,s B 0 ∂ T ∂ r + κ s iωρ 0,s α T,s B 0 T 0 × λ 0 + 2µ 0 iω + λ s + 2µ s ∂ ∂ r ( 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r )). ( 63
)
By multiplying Eq.( 63) by 1

r 2 ∂ ∂ r [r 2 (•)
] and using Eq.(62a), we get, using the notation introduced in Eq.(57), a differential equation for T

(λ s + 2µ s + λ 0 + 2µ 0 iω )κ s T iv -(λ s + 2µ s + λ 0 + 2µ 0 iω )iωρ 0,s C V,s + α 2 T,s B 2 0 T 0 + iωρ 0,s κ s T ii -ω 2 ρ 2 0,s C V,s T = ρ 0,s iωQ 0 . (64)
Apart from the supplied input power Q 0 , Eq.( 64) is similar to Eq.( 58) for fluids, with only B f changed to λ 0 + 2µ 0 for the solid. Henceforth, the solutions of the associated algebraic equation are given in Eq. (33). Moreover, the solution for T obtained in Eq.( 61) is still valid if we add the effect of the supplied energy defined in Eq. (22). Therefore, we finally get the general solution

T (r) = A e -θ ac,s r r +B e θ ac,s r r +C e -θ th,s r r +D e θ th,s r r +T Q 0 , (65)
with A, B,C, D integration constants.

As before, we analyze the pressure in the solid phase with spherical geometry. By definition, the pressure is given by

p = -Te n • n = -T e i j X i X j r 2 , (66)
with n = X r . Using strain and stress tensors defined in Eqs.( 25) and ( 26), respectively, and the assumption u = X r u(r), the elastic tensor can be written [START_REF] Landau | Theory of Elasticity[END_REF]

Te

=          2µ 0 ∂ u 1 ∂ X 1 + λ 0 ε kk -α T,s B 0 T µ 0 ∂ u 1 ∂ X 2 + ∂ u 2 ∂ X 1 µ 0 ∂ u 1 ∂ X 3 + ∂ u 3 ∂ X 1 µ 0 ∂ u 1 ∂ X 2 + ∂ u 2 ∂ X 1 2µ 0 ∂ u 2 ∂ X 2 + λ 0 ε kk -α T,s B 0 T µ 0 ∂ u 2 ∂ X 3 + ∂ u 3 ∂ X 2 µ 0 ∂ u 1 ∂ X 3 + ∂ u 3 ∂ X 1 µ 0 ∂ u 2 ∂ X 3 + ∂ u 3 ∂ X 2 2µ 0 ∂ u 3 ∂ X 3 + λ 0 ε kk -α T,s B 0 T          . (67)
The partial derivatives of u can be calculated as

∂ u i ∂ X j = X i X j r 2 ∂ u ∂ r + r 2 δ i j -X i X j r 3 u, (68)
and therefore

T e i j X i X j r 2 = 2µ 0 r 2 (X i X j ) 2 r 2 ∂ u ∂ r - (r 2 δ i j -X i X j )X i X j r 3 u + (λ 0 ε kk -α T,s B 0 T ) X i X j δ i j r2 , (69)
with

ε kk = Tr(ε) = 2 r u + ∂ u ∂ r .
This finally leads to

p = α T,s B 0 T - 2µ 0 iω ∂ v ∂ r - λ 0 iω 1 r 2 ∂ ∂ r (r 2 v). (70)
In this expression, the pressure p is only defined through the elastic stress tensor and therefore, in order to have a stress continuity relation between different adjacent media (solid, fluid), This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. we need to also define p. It is based on the total stress tensor obtained by adding elastic and viscous components

PLEASE CITE THIS ARTICLE AS

p = -T n • n = -Te n • n -TV n • n = p -2µ s ∂ v ∂ r -λ s 1 r 2 ∂ ∂ r (r 2 v) = α T,s B 0 T -2( µ 0 iω + µ s ) ∂ v ∂ r + ( λ 0 iω + λ s ) 1 r 2 ∂ ∂ r (r 2 v) . (71)
To conclude, the system of equations for a solid given in Eq.(62) has been solved in a spherical geometry. The general solution for the temperature T has been obtained, and the other variables can be calculated on the base of this quantity. Lastly, by definition of the elastic stress tensor, the pressure was determined in the solid with spherical symmetry. The viscous stress tensor was then added to define the normal component of the total stress tensor, useful to have a continuous quantities at interfaces.

C. Transfer matrix method with N spherical layers

A flexible model with multiple concentric spherical layers, as shown in Fig. 5, will be now elaborated using the solutions obtained in fluid and solid media. We consider a system with N spherical layers. The external one is supposed to be semiinfinite and each internal layer, j ∈ [1, N -1], potentially has a volumetric source Q 0, j (Q 0,N = 0). The thickness of the j-th shell is given by r jr j-1 , with r j being the external radius of the j-th layer. The type of material of each layer (fluid or solid) is not explicitly specified since materials are taken into account in the inherent parameters of their defining matrix, as explained below. In each spherical shell, the thermal flux q is defined as

q = -κ ∇T = -κ r r ∂ T ∂ r , (72)
and its scalar version is q = -κ ∂ T ∂ r . As before, the functions F and G are defined as

F (η) =α T B f - B f iω + λ + 2µ (L 1 η 2 + L 2 η 4 ), (73) G (η) =L 1 + L 2 η 2 . ( 74
)
Here, the coefficients L 1 and L 2 are defined through Eq.( 56) for a fluid and Eq.( 63) for a solid, when these expressions are written in the form v = L 1

∂ T ∂ r + L 2 ∂ ∂ r ( 1 r 2 ∂ ∂ r (r 2 ∂ T ∂ r ))
. We also remark that when Eq.( 73) is written for a solid layer, one must substitute B f with λ 0 + 2µ 0 and all parameters must correspond to those of the solid material. The general form of the variables p, v, q, T in a given shell can be directly found to be

                                                   p(r) = A F (θ ac ) -2µ 2 r ( 1 r + θ ac )G (θ ac ) e -θacr r +B F (θ ac ) -2µ 2 r ( 1 r -θ ac )G (θ ac ) e θacr r +C F (θ th ) -2µ 2 r ( 1 r + θ th )G (θ th ) e -θ th r r +D F (θ th ) -2µ 2 r ( 1 r -θ th )G (θ th ) e θ th r r +α T B f T Q 0 , v(r) = -A( 1 r + θ ac )G (θ ac ) e -θacr r -B( 1 r -θ ac )G (θ ac ) e θacr r -C( 1 r + θ th )G (θ ac ) e -θ th r r -D( 1 r -θ th )G (θ th ) e θ th r r , q(r) = A( 1 r + θ ac )κ e -θacr
Those equations apply for both fluid and solid layers with appropriate changes of parameters. In particular, we have to replace B f with λ 0 + 2µ 0 in F and G , and we have to replace µ with µ + µ 0 iω in the expression for p, which become for solid layer

p(r) = A F (θ ac ) -2(µ + µ 0 iω ) 2 r ( 1 r + θ ac )G (θ ac ) e -θ ac r r +... (76)
We can rewrite Eq.(75) in matrix form in the j-th layer as

   p v q T    j = W j    A j B j C j D j    +    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j    , (77)
with T Q 0 , j = Q 0, j ρ 0, j C V, j iω . We also defined W j (r) = W can be written in an ar-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458 bitrary layer as

W (a) =          F (θ ac )+ F (θ ac )+ -2µ 2 r ( 1 r + θ ac )G (θ ac ) -2µ 2 r ( 1 r -θ ac )G (θ ac ) -( 1 r + θ ac )G (θ ac ) -( 1 r -θ ac )G (θ ac ) ( 1 r + θ ac )κ ( 1 r -θ ac )κ 1 1 ... ... F (θ th )+ F (θ th )+ -2µ 2 r ( 1 r + θ th )G (θ th ) -2µ 2 r ( 1 r -θ th )G (θ th ) -( 1 r + θ th )G (θ th ) -( 1 r -θ th )G (θ th ) ( 1 r + θ th )κ ( 1 r -θ th )κ 1 1          , and
W (b) =      e -θacr r 0 0 0 0 e θacr r 0 0 0 0 e -θ th r r 0 0 0 0 e θ th r r      , (78)
where all parameters must correspond to the layer under consideration. Here, as before, when we deal with a solid layer, we must substitute B f with λ 0 + 2µ 0 and we have to consider the correct values for θ ac and θ th . All variables within our system are defined for r > 0. In order to avoid any singularity at r = 0, we impose in the first layer the conditions A 1 = -B 1 and C 1 = -D 1 , coming from the simple property e z ∼ 1 + z when z → 0. The last N-th layer is consider to be semi infinite and only progressive waves are therefore assumed. It means that we have to impose B N = D N = 0 (absence of regressive waves). By imposing the continuity of the main four variables between adjacent layers we directly find

   A 1 -A 1 C 1 -C 1    = M 0    A N 0 C N 0    + ∑ n∈[1,N-1] M Q 0,n    α T,n B f ,n T Q 0 ,n 0 0 T Q 0 ,n    , (79)
with n summed over all layers with a volumetric source and where

M 0 = W -1 1 (r 1 ) N-1 ∏ j=2 W j (r j-1 ) W -1 j (r j ) W N (r N-1 ), (80)
and

                 M Q 0,1 = -W -1 1 (r 1 ), M Q 0,n = W -1 1 (r 1 ) ∏ n-1 j=2 W j (r j-1 )W -1 j (r j ) × 1 -W n (r n-1 )W -1 n (r n ) , ∀n ∈ [2, N -1], M Q 0,N = W -1 1 (r 1 ) ∏ N-1 j=2 W j (r j-1 )W -1 j (r j ) . (81)
Of course, as discussed before, this is true only if we neglect any interface imperfection. [START_REF] Pavanello | Generalized interface models for transport phenomena: Unusual scale effects in composite nanomaterials[END_REF][START_REF] Pavanello | How imperfect interfaces affect the nonlinear transport properties in composite nanomaterials[END_REF] We remark that Eq.( 79) can be rewritten in a more explicit form as

   1 0 -M 0 (1, 1) -M 0 (1, 3) -1 0 -M 0 (2, 1) -M 0 (2, 3) 0 1 -M 0 (3, 1) -M 0 (3, 3) 0 -1 -M 0 (4, 1) -M 0 (4, 3)       A 1 C 1 A N C N    =                   ∑ n M Q 0,n (1, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (1, 4)T Q 0 ,n ∑ n M Q 0,n (2, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (2, 4)T Q 0 ,n ∑ n M Q 0,n (3, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (3, 4)T Q 0 ,n ∑ n M Q 0,n (4, 1)α T,n B f ,n T Q 0 ,n + ∑ n M Q 0,n (4, 4)T Q 0 ,n                   . ( 82
)
Then, we observe that Eq.(79) or Eq.(82) represents a system of four linear equations giving the values of the parameters A 1 , C 1 , A N , and C N . This allows the calculation of any coefficient (ABCD) j using either the decreasing recursion

   A j B j C j D j    =W -1 j (r j )   Wj+1(r j )    A j+1 B j+1 C j+1 D j+1    +    α T, j+1 B f , j+1 T Q 0 , j+1 0 0 T Q 0 , j+1    -    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j       , (83) ∀ j ∈ [1, N -1], or the increasing recursion    A j B j C j D j    =W -1 j (r j-1 )   Wj-1(r j-1 )    A j-1 B j-1 C j-1 D j-1    +    α T, j-1 B f , j-1 T Q 0 , j-1 0 0 T Q 0 , j-1    -    α T, j B f , j T Q 0 , j 0 0 T Q 0 , j       , (84)
∀ j ∈ [2, N]. Finally p, v, q, T at any point, i.e. for any r, of each layer can be found with Eq.(77).

As before, this methodology allows for studying the thermoacoustic wave generation and propagation in arbitrarily complex systems with spherical geometry.

IV. ANALYSIS OF THERMOPHONE STRUCTURES GENERATING CYLINDRICAL AND SPHERICAL WAVES

In this Section, the cylindrical and spherical models, elaborated in Sections II and III, are used to analyze different multilayered thermophone structures generating cylindrical and This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458 The results for cylindrical and spherical models are here compared and it will be shown that they display a rather similar behavior. The cylindrical and spherical geometries are non conventional shapes for thermophones at the macroscopique scale. Nonetheless, those geometries are worth being investigated for different reasons. For instance, the cylindrical shape analysis could be applied to the microstructure of porous thermophones, where the foam struc-ture is often constituted of hollow nano-cylinders connecting the pores. This is the case, for example, for carbon foams or three-dimensional graphene. [START_REF] Ngoh | Experimental characterization of three-dimensional Graphene's thermoacoustic response and its theoretical modelling[END_REF] Additionally, the spherical geometry is of interest since, if the radius of the sphere is very small, it would correspond to a thermoacoustic monopole. As such, acoustical theory relative to monopoles could be applied assuming only a different input pressure generated from the monopoles. This would allow for easy comparison with regular sound generator and facilitate thermoacoustic design. In this Section, will focus on three different cases, which can be explained by means of Fig. 6: • a full cylinder or sphere (R b = 0) composed of the thermophone material;

ρ [kg• m -3 ] C p [J•kg -1 • K -1 ] C v [J•kg -1 • K -1 ] B f [Pa] α T [K -1 ] Gas,
λ [Pa•s] µ [Pa•s] λ 0 [Pa] µ 0 [Pa] κ [W•K -1 •m -1 ] Gas,
• a hollow cylinder or sphere (0 < R b < R s ) with air in the core region;

• a hollow cylinder or sphere (0 < R b < R s ) with a different material (the so-called substrate) in the core region.

In all cases, the external propagating medium is typically air or water even if other materials can be envisaged. Nevertheless, for the sake of definiteness, our results are obtained and shown only with air as external medium. The full thermo-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

phone structure, without core, will allow to estimate the radiation of an infinitely thin and long cylindrical thermophone as well as a monopole. The two other cases are motivated by the chemical vapor deposition (CVD) method of fabrication of porous graphene thermophone [START_REF] Ngoh | Experimental characterization of three-dimensional Graphene's thermoacoustic response and its theoretical modelling[END_REF][START_REF] Loeblein | Configurable three-dimensional boron nitride-carbon architecture and its tunable electronic behavior with stable thermal performances[END_REF] . More specifically, in graphene foams, the graphene sheets are grown on a nickel template, which is then removed leaving the foam composed of hollow branches. The case with substrate at the core represents a branch of foam before removal of the template while the hollow case represents branches after nickel removal. The parameters adopted for thermophone material, air, water and substrate are those shown in Table I and Table II. The parameters describing the substrate correspond to the titanium, [START_REF] Tong | Gas-Filled Encapsulated Thermal-Acoustic Transducer[END_REF] while the parameters describing the thermophone active layer correspond to the three-dimensional carbon. [START_REF] Loeblein | Configurable three-dimensional boron nitride-carbon architecture and its tunable electronic behavior with stable thermal performances[END_REF] Moreover, Table III allows for an overview of the presented plots and their varying parameters. In all plots, the recording distance from the thermophone surface was fixed at ∆r 0 = 50 × 10 -6 m (acoustic near field). In order to facilitate the comparison between different plots, the input power is normalized at P in = 1W, and the following procedure has been applied. For the cylindrical geometry we can write

Q 0 = P in /V , V = πL(R 2 s -R 2 b
) and A = 2πR s L (L being the length of the cylinder and A its external surface). Similarly for the spherical geometry we have that

Q 0 = P in /V , V = (4/3)π(R 2 s -R 2 b ) and A = 4πR 2
s . Therefore, with fixed radii and P in = 1W, we can easily calculate the effective cylinder length L and the thermal power density in each configuration. Here, the thermophone external surface A was assumed to be the same for both cylindrical and spherical models. The sound pressure level (SPL) is defined as SPL=20log 10 (p rms /p re f ) where p rms is the root mean square pressure | p | / √ 2 (with p the complex pressure previously defined) and p re f is the reference sound pressure being, by definition, 20µPa in air and 1µPa in water.

A. Comparison with models from the literature

In order to validate our multilayer models, they must first be confronted to models from the literature. It was seen in Section I that very few models actually investigated cylindrical or spherical thermoacoustic sound generation. Moreover, not many thermophones have a cylindrical or a spherical macroscopic geometry and therefore it is not possible to draw a direct comparison between experimental data and our models for the moment. Nevertheless, one investigation concerning those geometries has been made by Hu et al., where they extended their plane wave geometry model 26 to spherical [START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF][START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF] and then cylindrical [START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] waves. Those models are based on the conservation equations of continuum mechanics. Alternatively, a piston based model for cylindrical shapes is also found in the Tong et al. work. [START_REF] Tong | Theory and modeling of cylindrical thermo-acoustic transduction[END_REF] This technique was then applied to model thinline arrays. [START_REF] Liu | Thermo-acoustics generated by periodically heated thin line array[END_REF] All these models have been implemented and compared with our approach in the next discussion.

Cylindrical geometry

In Fig. 7, a full cylinder radiating in free field is considered and our results are compared to the Tong model [START_REF] Tong | Theory and modeling of cylindrical thermo-acoustic transduction[END_REF] and to the Yin and Hu model [START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] . In his model, Tong and colleagues consider a uniformly heated cylindrical conductor and therefore they neglect the influence of the heat capacity of the thermophone. This is justified by the fact that in their future work they investigate thinline arrays whose length are greatly larger than the radius of cylinder. They only defined the radius in order to take into account the input power density (see Eq.( 5) of Ref. [43]). In our model, we also assumed a cylinder length much larger than the radius to simplify the mathematical problem. Nonetheless, the radius is of critical importance, not only for the input power density but also for the HCPUA calculation and for mechanical resonance determination. The final pressure equation for the Tong model is written as p g = DH (1) 0 (k g r) + p * g with D a constant determined in Eq.( 12) of Ref. [43] and p * g a particular solution found in Eq.( 11) of Ref. [43].

Regarding Yin and Hu model, it takes into account the thermal wave inside the solid but assumes that the input power is only distributed at the external surface of the thermophone. In fact, all the models developed by Hu and colleagues (for both planar, cylindrical and spherical thermoacoustic wave generation) [START_REF] Hanping | Model for thermoacoustic emission from solids[END_REF][START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF][START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF][START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] take into consideration the electric and then thermal energy entering the system distributed only on the external surface of the active layer of the thermophone. It means that they consider the material covered by a conductive 2D film that provides the energy to the system. They adopted this particular geometry since in their first planar analysis they wanted to reproduce the experimental results obtained by Shinoda et al., [START_REF] Shinoda | Thermally induced ultrasonic emission from porous silicon[END_REF] actually using this kind of configuration in the experimental setup. Therefore, to compare these results with our thermophone structure where the energy is distributed uniformly within the whole volume of the active layer, we conveniently modified the models by Hu and colleagues. This modification simply removes the surface energy distribution and replaces it with a volume energy distribution. We performed this operation for the case of planar waves in Ref. [35], and we apply here the same procedure for cylindrical and spherical acoustic generation. For the sake of brevity, we omitted all the mathematical details of this modification but all developments can be found in Ref. [44]. Anyway, the resulting pressure equation has been found as

p g = Q 0 iωρ s C p,s κ s M 0 σ s d a d t k g K 1 (k g R s )K 0 (σ g r) -σ g K 1 (σ g R s )K 0 (k g r) d a k g K 1 (k g R s ) κ s M 0 σ s K 0 (σ g R s ) -κ g σ g K 1 (σ g R s ) + -d t σ g K 1 (σ g R s ) κ s M 0 σ s K 0 (k g R s ) -κ g k g K 1 (k g R s ) =p g,sur f ace Q 0 iωρ s C p,s κ s σ s M 0 , (85)
for volume density power. Here, p g,sur f ace is the original result of the Yin and Hu model with a surface density power. [START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] Moreover, the quantities involved in Eq.( 85) are defined in the original publication. [START_REF] Yin | Analysis of cylindrical thermo-acoustic radiation in gas[END_REF] Here, this equation does not take into account a substrate. In Fig. 7, panel a), one can find the comparison of the frequency spectrum obtained through our model for a full cylindrical thermophone with the literature results discussed above. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. It is seen that there is a good agreement between the original Yin et Hu model with the presented model up to 20kHz. There is then an increase of SPL in the Yin and Hu original model that reaches a plateau between 1 and 8MHz before decreasing in the very high frequency range. This overestimation in the high frequency range is due to the neglect of the thermal capacity of the thermophone. Indeed, the power being supplied only on its surface, the energy stored by means of the heat capacity is strongly underestimated. On the other hand, the modified Yin and Hu model, with volume density power, is coherent with our model up to around 10MHz. After 10MHz both models exhibit a decreasing behavior of SPL, although at a higher rate in our case. This is due to the way the wavenumbers are calculated in the Yin and Hu model. Similarly as for plane wave, [START_REF] Hanping | Model for thermoacoustic emission from solids[END_REF] Yin and Hu used a zero order approximation for the thermal and acoustical wavenumbers with respect to viscosity and thermal conduction. If the first order approximations defined in Eq.( 17) are used in Yin and Hu model, the consistency between both models is observed up to around one hundred MHz, as seen in Fig. 7,panel b). This proves that, in the high frequency range, it is necessary to take losses into account in the wavenumbers calculation. More specifically, thermal conduction losses are crucial to get a correct high frequency spectrum.

PLEASE CITE THIS ARTICLE AS

Lastly, a small deviation can be seen between our model and the modified Yin and Hu model around 1kHz. This is due to the chosen special functions, solutions of the temperature differential equation. In Eq.( 18) the progressive and regressive waves are defined through Hankel functions of the first and second kind, whereas Yin and Hu actually used modified Bessel function. A better agreement between the spectra can be found in this frequency range if the modified Bessel functions are replaced by Hankel functions.

Regarding the Tong model in Fig. 7, panel a), the 20dB/dec increase in the low frequency range fits with our model. How-ever, due to the lack of consideration of the wave propagation within the thermophone region as well as an assumed uniform heat distribution, there is an overestimation of SPL above 1kHz. In a similar fashion as for the surface density model of Yin and Hu, the thermal capacity of the thermophone is neglected in the Tong model, making this approach only valid for very thin thermophones or at very low frequencies. Notably, a zero order approximation is used for the wavenumbers determination. However, since the model is oversimplified, when a first order approximation is used, the spectrum remains overestimated, as seen in Fig. 7, panel b).

It is worth discerning that there are no acoustic resonances due to size and shape of the thermophone in neither literature models. Only our multilayer model takes into account the acoustic propagation inside the solid, allowing us to determine antiresonances caused by disruptive interference. Models in the literature neglected this influence since it was assumed to be in a frequency range higher than the one of interest. The first antiresonance observed in our model is seen around 1MHz, which is much higher than the human hearing range, making the assumption above reasonable in most of the literature. Antiresonances at higher frequencies are also observed and display variable intensities due to the chosen sampling step in the frequency axis. All antiresonances should be of similar amplitude, as can be seen by introducing a very small frequency sampling step in the analysis. We can confirm that all the notches observed in Fig. 7 (and the following ones) are real antiresonances of the system (and not numerical artifacts) since we checked the conditioning number of the matrices defining the model and we observed a regular working of the numerical procedure in this range of frequencies. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

Spherical geometry

In Fig. 8, panel a), a full sphere radiating in free field is considered and our results are compared to the model of Hu et al. (2014). [START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF] Moreover, in Fig. 8, panel b), a hollow sphere with core (substrate) is considered and our results are compared with the model of Hu et al. (2012). [START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF] As already discussed for plane and cylindrical waves, the Hu models assume a surface density energy input and therefore they must be modified in order to be comparable with our approach. Assuming a volume density power input, the modified pressure equation for a full sphere becomes

p g = R s r d a d t k + g Z (e -σ g ∆r - σ + g k + g e -k g ∆r ) Q 0 iωρ s C p,s κ s σ - s * = p g,sur f ace Q 0 iωρ s C p,s κ s σ - s * , (86)
and the one for a sphere with core (substrate) is

p g = d a d t k + g Z (e -σ g ∆r - σ + g k + g e -σ g ∆r ) × κ s Q 0 iωρ s C p,s r R s σ - s * + R b κ b σ - b * κ b σ - b * -κ s σ - s σ s (1 -M)e σ s ∆R = p * g,sur f ace r κ s Q 0 iωρ s C p,s R s σ - s * + R b κ b σ - b * κ b σ - b * -κ s σ - s σ s (1 -M)e σ s ∆R . ( 87
)
It is easily seen that Eq.( 86) is equivalent to Eq.(87) if we impose M = 1. The form of these expressions is similar to the modified Hu equation used for the plane wave case in a recent paper. [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] For all parameters used in Eqs.( 86) and (87), the reader is asked to refer to the original publications. [START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF][START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF] Similarly to Fig. 7 a), Fig. 8 a) (full sphere) displays a change of behavior between the original Hu model with surface density power and the modified one with volume density power. On the one hand, the original version keeps increasing after 20kHz, before reaching a plateau between 1MHz to 80MHz and then decreases. This is once again due to the neglect of the heat capacity of the thermophone and to the assumption of a surface density power. On the other hand, from 10kHz to 10MHz there is a good agreement between the modified model and our multilayer one. Above 10MHz the same issue as for plane and cylindrical waves is observed, where the zero order approximation of the wavenumbers in Hu models prevent us from observing the same decrease. A perfect match is found if the wavenumbers are changed accordingly (not shown here for the sake of brevity).

Below 10kHz, in Fig. 8, panel a), the observed spectrum slopes in both Hu models and ours are different. This interesting point is explained by the adopted distance at which we calculate the physical quantities of interest (e.g., the SPL). In Fig. 8, the pressure is calculated at ∆r 0 = 50µm distance. However, based on Eq.(17b) the thermal layer can be estimated with the approximated equation

L th = 2 C 0 l k 2ωγ = 2α T
ω and then L th at 3kHz is about 50µm in air. It is known that L th is longer for lower frequency and therefore our recording distance is not large enough to have an accurate evaluation of the acoustic pressure. In this low frequency range, we are at a distance ∆r 0 inside the thermal layer in which the thermoacoustic process occurs leading to altered results. The properties of the acoustic wave should be calculated sufficiently far from the thermal layer. This scenario can be observed only in our model where all layers are fully studied with all relevant physical variables. In this particular case, the Hu model does not implement the thermal layer and therefore the acoustic wave is fully accessible just outside the generating layer. Anyway, if the position at which we cal-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. culate the SPL is chosen further away (larger ∆r 0 ), there is a better agreement between the models also at low frequencies. Now regarding Fig. 8, panel b) (hollow sphere with core), similar observations can be made between the original and modified Hu models at high frequencies. In fact, above 20kHz only little differences can be observed between the spectra of the full sphere and the hollow sphere with substrate. Below 20kHz, due to the small size of the substrate (R b = R s /3) thermal losses will occur where the thermal layer is longer than the radius of the thermophone. As observed for plane waves, [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] energy is lost in the substrate and results in the decrease in the calculated SPL. It is also seen that numerical instabilities can appear at low frequencies in our model, making Hu et al. model more robust in this specific case.

PLEASE CITE THIS ARTICLE AS

Finally, like for plane and cylindrical waves, in both plots of Fig. 8 only our model display high frequency antiresonances. This is due to the fact that this is the only model with a full implementation of the acoustic field within the solid layers.

B. Analysis of thermophone geometries

The validity of our multilayer model for cylindrical and spherical geometries has now been proven by the comparison with other models from the literature. As previously explained, it is important to investigate those geometries as foam like thermophones have cylindrical branches at a microscopic scale. Additionally, acoustic theory based on monopoles (generating spherical wave) is very useful and can be easily applied for thermophones if a monopole thermoacoustic response is known. In the following, the three configurations of interest previously defined will be investigated in more detail and their defining parameters will be modified to better assess their influence. In particular, both frequency and position dependent analyses will be performed. along the z direction is closer to the observation point than the surface of the sphere. Even though diffraction is not taken into account in this analysis, the inherent shape of the cylindrical thermophone induces more efficient low frequency radiations. Finally, Fig. 9 show that computational limitations are more easily reached at low frequency for spherical geometry and for thin thermophones. Fig. 10 displays the spatial behavior of the temperature variation, particle velocity and heat flux within a full cylinder and a full sphere for comparison. We considered a frequency of 10kHz and we plotted these quantities from the center of the thermophone to twice the recording distance used for the frequency spectrum analysis. It is seen that both cylinder and sphere display very similar profiles. Notably, we remark that at the center of the thermophone the temperature exhibits a maximum, while the heat flux and the velocity are zero. At this frequency the thermal layer length is 26µm and can be seen in all figures but mostly in the velocity plot.

Hollow cylindrical or spherical thermophone

We show in Fig. 11 the frequency response of thermophones shaped as hollow cylinders and hollow spheres, with a fixed external radius R s and a variable internal radius R b . Air is assumed at the core of those thermophone structures. Fig. 11, panel a), show that for a small internal radius, relatively to the external one, the frequency spectrum is similar to the one of a full sphere. With increasing R b , the heat capacity of the thermophone is reduced because of its thinner thickness. This leads to an increase in high frequency SPL, while the low frequency range is not sensibly modified. Another observable change is the slight modification of the antiresonance frequency for large internal radius. New resonances due to the new geometry are observed above 1MHz. However, those changes are not as significant as the modifications induced by changing the external radius previously observed for a full thermophone. Those changes in resonances would not have This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. been observable without our multilayer model. We also observe that the spectra for R b = 0.1R s and R b = 0.01R s are very similar since very small cavities do not modify the heat capacity of the system and therefore their effects is not discernible. Indeed, for R b approaching zero we obtain again the response of the full thermophone. In a similar fashion as for a full thermophone (Section IV B 1), the behavior of a hollow sphere (Fig. 11, panel b) displays a similar behavior as a hollow cylinder (Fig. 11, panel a). The same analysis can then be applied.

PLEASE CITE THIS ARTICLE AS

Finally, Fig. 12 shows the spatial behavior of the main variables at 10kHz for this new geometry. We adopted R b = R s /3. Similar behaviors are observed for both the hollow cylinder and the hollow sphere. The temperature exhibits its maximum value in the thermophone material at the internal radius R b and the thermal layer is observed on the outer boundaries of R s . The heat flux and the velocity are also zero at r = 0. It is then not possible to observe the thermal layer inside the inner radius as waves cancel each other out.

Hollow cylindrical or spherical thermophone with substrate

In Fig. 13, we present the results for hollow cylinder and hollow sphere thermophones, where a given substrate is assumed in the internal core. The frequency spectra for a cylindrical and for a spherical shapes are found in Fig. 13, panel a), and 13, panel b), respectively. Here we considered a varying inner radius of the core. Two main behaviors are observable in Fig. 13, panel a).

Firstly, below 10kHz, it is seen that for larger values of the core radius more losses are observed and the resulting SPL is consequently reduced. In this case, a considerable quantity of heat is dissipated in the substrate and this phenomenon is even more pronounced at low frequencies with a larger thermal layer. Conversely, the cylinder with the smallest core has the highest SPL as it is closer to a full cylinder.

Secondly, above 10kHz, since the thickness of the thermophone is reduced by increasing R b , its HCPUA is reduced as well and this trend increases the SPL at high frequencies (as previously seen for a hollow thermophone). Furthermore, when the thickness of the thermophone R s -R b starts to be similar to the thermal layer length, the influence of the substrate diminishes, as already observed for plane waves. [START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] This means that the cylinder with the thickest substrate is much more efficient above 1MHz. Summing up, we underline that thin substrates are recommended at low frequencies and thick substrates are convenient at high frequencies.

Additionally, the antiresonance wavelength of the thermophone is directly proportional to its size but also to its stiffness. Therefore, having a substrate at the core changes the overall or effective stiffness of the thermophone shifting the first antiresonance to a higher frequency range.

Similar observations can be made for Fig. 13, panel b). Again, we observe that the numerical model for the spherical case can be unstable at low frequencies, especially for thin core substrate.

The results in Fig. 14 are complementary to those of Fig. 13 as they consider the same structure but with a fixed ratio between inner and outer radius (R s = 2R b ) and a variable outer radius. Across the whole frequency range, the thinner thermophone has overall the best performance. Its low HCPUA, due to its small size, improves the high frequency radiation. As the same time, its thin core with the substrate also limits the losses at low frequencies. At last, the shift toward the high frequency range of the first antiresonance of the solid for a thicker radius is also clearly distinguishable. It is however reminded that smaller samples have a higher input power density, which theoretically helps the observation and analysis but might create real experimental issues both in terms of fabrication and heat density management.

The spatial behavior of temperature variation, particle velocity and heat flux for this configuration can be found in This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

FIG.
12. Temperature variation, particle velocity, and heat flux calculated through the multilayer model in hollow cylinder (continuous blue curves) and hollow sphere (dashed red curves) thermophones. All curves show the absolute value of the corresponding complex quantities. While the air is represented by the yellow regions, the green region corresponds to the thermophone material. Here, we adopted f=10kHz.

Fig. 15. We can observe that the temperature drops inside the substrate. The heat flux entering the substrate is larger at R b than at R s , as opposed to Fig. 12 where the heat flux was at its maximum in correspondence to the outer region of the thermophone material (where the thermoacoustic generation occurs). This leads to a smaller particle velocity and thus a reduced SPL with respect to Fig. 12.

V. CONCLUSIONS

We introduced in this work a thermoacoustic model able to describe the acoustic wave generation in multilayered thermophones with cylindrical or spherical geometry. One important point is that we introduced the thermal conduction and the viscous behavior in all layers, which are therefore characterized by the fully coupled thermo-visco-acoustic constitutive response. From the thermodynamic point of view, these two contributions describe the dissipative behavior of the system. up to a dozen layers. The calculation time may increase with a larger number of layers but it remains extremely competitive compared to numerical methods based on finite elements or other computational techniques. This efficiency is based on the use of linearized equations for wave propagation. Of course, if we want to include nonlinear phenomena or temperature dependent parameters, then we are forced to switch to more expensive numerical techniques.

In Section IV A, our multilayer cylindrical and spherical models were confronted to models from the literature. Actually, very few models currently exist for cylindrical or spherical geometry of thermophones but our model proved to be coherent with those currently available. In particular, it has been proven to be able to fit different configurations (full thermophone, hollow thermophone and hollow thermophone with substrate), while being the only model considering a refined estimation of the wavenumbers as well as the mechanical resonances and/or antiresonances generated by the solid layers. Indeed, we fully implemented in our model the elastic wave propagation within the solid layers, suitably coupled with the thermal and viscous behavior. We underline that the multilayer model sometimes can exhibit numerical instabilities due to the computation of the inverse of not well-conditioned matrices. Those limitations have been observed for low frequencies and particularly in spherical geometries.

In Section IV B, the influence of the geometrical parameters of the cylindrical or spherical thermophones on the thermoacoustic generation was analyzed. Firstly, it was seen that cylindrical and spherical systems have displayed similar results and the qualitative physical analysis is valid in both cases. Importantly, similarly to the one-dimensional case with plane waves, the thickness of the thermophone directly influ-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. and spherical (dashed red curves) thermophones with a substrate as internal core. All curves show the absolute value of the corresponding complex quantities. Three different regions can be observed: the substrate (red region), thermophone (green region), and the propagating air medium (yellow region). Here, we adopted f=10kHz.

PLEASE CITE THIS ARTICLE AS

ences the HCPUA of the device. Thinner thermophone are then privileged for high frequency radiation, although high power density has to be expected and one has to be careful not to damage real samples. The use of a substrate at the core was proven to lead to thermal losses as expected. To optimize this configuration, thin substrates are recommended at low frequencies and thick substrates are convenient at high frequencies.

Concerning the applications, the spherical model can be useful to obtain the response of the thermoacoustic monopole, which can be integrated in more complex simulation systems. Moreover, the cylindrical model can be used to study nanoor micro-structured thermophones. In this case, the thickness of layers has once again proven to be a decisive factor for HCPUA tuning. In particular, carbon foam thermophone are characterized by hollow cylindrical branches in the porous microstructure. Interestingly, the growth of branches by CVD can adapt their thickness by controlling the fabrication process and growth time. This is useful to optimize the thermoacoustic response. Finally, in the CVD based samples production we can choose whether to keep the substrate (typically a nickel matrix) to improve the mechanical stability, or to remove the matrix to increase the thermoacoustic efficiency at high frequencies, at the risk of high input power density. Also intermediate solutions can be envisaged.

Through these general solutions, wavenumbers θ ac, f and θ th, f associated to acoustical and thermal waves respectively can be approximated to the first order in l V and l k by

θ ac, f = ± iω C 0 1 - 1 2 iω C 0 l V - 1 2 iω C 0 l k (1 - 1 γ ) , (B5)
θ th, f = ± iωγ C 0 l k 1 + 1 2 iω C 0 l k (1 - 1
γ ) + 1 2 iω C 0 l V (1 -γ) .
(B6)

Wavenumbers in solid media

In order to explicit the differences in wavenumbers between solid and liquid media, a first order approximation solution of Eq.( 21) will be determined for the associated homogeneous fourth order differential equation of T . Eq.( 21) becomes 0 =(λ s + 2µ s + λ 0 + 2µ 0 iω )κ s θ 4 -(λ s + 2µ s + λ 0 + 2µ 0 iω )iωρ 0,s C V,s + α 2 T,s B 2 0 T 0 + iωρ 0,s κ s θ 2

-ω 2 ρ 2 0,s C V,s . (B7)

It will be seen that this equation is also the one associated to spherical wave radiation following Eq.(64) in the following Section III B. Using similar notations as in Eq.( 15) we define (B11) Using Eq.(B11), Eq.(B10), the first order Taylor's series approximation and a first order approximation on l κ,s and l V,s (weak conduction and weak viscosity)

                   γ s = C p,s C V,
                           1 ± ac b 2 = 1 ± 1 2 l κ,s C L iω -ω 2 C 2 0,s (1 -1 γ s + C L C 0,s ) -2 , -c b = iω C 0,s (1 -1 γ s + C L C 0,s ) -1/2 × 1 + iω C 0,s (l κ,s + l V,s )(1 -1 γ s + C L C 0,s ) -1/2 , -b a = ( iω l κ,s C L ) 1/2 (1 -1 2 iω C L l V,s ) ×(1 -1 γ s + C L C 0,s ) 1/2 × 1 + 1 2 iω C 0,s (l κ,s + l V,s )(1 -1 γ s + C L C 0,s ) -1 .
(B12) This leads to AC generator ( ω 2 )

                         θ ac,s = iω C 0,s (1 -1 γ s + C L C 0,s ) -1/2 × 1 -l V,s iω 2C 0,s (1 -1 γ s + C L C 0,s ) -1 -l κ,s iω 2C 0,s × (1 -1 γ s + C L C 0,s ) -1 -C L C 0,s (1 -1 γ s + C L C 0,s ) -2 , θ th,s = ( iω l κ,s C L ) 1/2 (1 -1 γ s + C L C 0,s ) 1/2 1 + iω 2C 0,s × l κ,s (1 -1 γ s + C L C 0,s ) -1 -(1 -1 γ s + C L C 0,s ) -2 C L C 0,s + l V,s (1 -1 γ s + C L C 0,s ) -1 - C 0,s C L . ( B13
Heat wave (ω):

x

FIG. 2 .

 2 FIG. 2. Example of a cylindrical thermophone composed of three cylindrical regions and two interfaces. The external region represents the propagating medium of the emitted cylindrical thermoacoustic wave.

  ac, f r) + BH

  ac, f r) +CH (1) 0 (iθ th, f r) + DH

  ∂ ∂t = iω, Q 0 (the heat density energy supplied) spatially uniform, b = 0 and knowing that ∂ u ∂t = v. When injecting 1 r ∂ ∂ r (rv) from Eq.(19a) into Eq.(19b), v can be written as a function of T alone as

  th,s r) + DH

FIG. 3 .

 3 FIG. 3. Schematic of a generalized multilayered thermophone model with N cylindrical layers.

FIG. 4 .

 4 FIG. 4. Example of a spherical thermophone composed of three spherical regions and two interfaces. The external region represents the propagating medium of the emitted spherical thermoacoustic wave.

r 2 ∂∂

 2 62) with ∂ ∂t = iω, Q 0 (the heat density energy supplied) spatially uniform, b = 0 and knowing that ∂ u ∂t = v. When injecting 1 r (r 2 v) from Eq.(62a) into Eq.(62b), v can be written as a function of T alone as

FIG. 5 .

 5 FIG. 5. Schematic of a generalized multilayered thermophone model with N spherical layers.

r+ B( 1 r

 1 θ ac )κ e θacr r +C( 1 r + θ th )κ e -θ th r r + D( 1 rθ th )κ e θ th r r , T (r) = A e -θacr r + B e θacrr +C e -θ th r r + D e θ th r r + T Q 0 .

FIG. 6 .

 6 FIG. 6. Schematic of the investigated models for cylindrical and spherical wave generation. The core of the thermophone (air or substrate) in light yellow has a radius R b and the external radius of the thermophone layer (light green) is R s

FIG. 7 .

 7 FIG. 7. Frequency spectrum of a full cylinder radiating in free field. The presented multilayer model is compared to the Tong model 43 and to the Yin and Hu model 30 (original, with surface energy density, and modified, with volumetric energy density). The wavenumbers for Tong model and for Yin and Hu model are calculated through a zero order approximation as in the original literature in panel a) and corrected with the first order approximation (see Eq.(17)) in panel b).

FIG. 8 .

 8 FIG.8. Frequency spectrum of a full spherical thermophone in panel a), and of a hollow spherical thermophone with substrate at the core in panel b), both radiating in free field. The proposed multilayer model is compared to the Hu model for the spherical geometry with and without substrate.[START_REF] Hu | Analysis of spherical thermoacoustic radiation in gas[END_REF][START_REF] Hu | Solution for acoustic field of thermoacoustic emission from arbitrary source[END_REF]

FIG. 9 .

 9 FIG.9. Frequency spectrum of a full cylinder (a) and a full sphere (b), with a variable radius R s .

FIG. 10 .

 10 FIG.10. Temperature variation, particle velocity, and heat flux calculated by the proposed multilayer model for full cylindrical (continuous blue curves) and full spherical (dashed red curves) thermophones. All curves show the absolute value of the corresponding complex quantities. The thermophone material (green region) is on the left and the propagating air medium (yellow region) is on the right. Here, we used f=10kHz.

FIG. 11 .

 11 FIG. 11. Frequency spectrum of a hollow cylinder in panel a), and of a hollow sphere in panel b), with a fixed external radius R s and a variable internal radius R b .

FIG. 13 .

 13 FIG. 13. Frequency spectrum of a cylinder in panel a), and a sphere in panel b), with a substrate as internal core. The external radius R s is fixed and the internal radius R b is variable.

FIG. 15 .

 15 FIG.15. Temperature variation, particle velocity, and heat flux calculated with the multilayer model for cylindrical (continuous blue curves) and spherical (dashed red curves) thermophones with a substrate as internal core. All curves show the absolute value of the corresponding complex quantities. Three different regions can be observed: the substrate (red region), thermophone (green region), and the propagating air medium (yellow region). Here, we adopted f=10kHz.

2 ac b 2 2 ac b 2 ,

 2222 s , C 0,s = B 0 ρ 0,s γ s , [m/s] C L = λ 0 +2µ 0 ρ 0,s C 0,s , [m/s] l κ,s = C 0,s κ B 0 C p,s , [m] l V,s = λ s +2µ s ρC 0,s , [m].(B8)Following the same thought process as for fluids, Eq.(B7) can be written l κ,s l V,s + l κ,s of the asymptotic solutions for the equations of the form aθ 4 + bθ 2 + c = 0 has been proven to be   θ ac,s = ± -c b 1 + 1 , θ th,s = ± -b a 1 -1(B10) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5b = -(1 -1 γ s + C L C 0,s ) 1 + iω C 0,s (l κ,s + l V,s )(1 -1 γ s + C L C 0,s ) -1 , c = -ω 2

  air (g) 1.20 9.96 × 10 2 7.17 × 10 2 1.01 × 10 5 3.33 × 10 -3 Thermophone (s) 0.03 × 10 3 2.38 × 10 2 2.38 × 10 2 2.78 × 10 5 0.6 × 10 -6 Substrate (b) 4.5 × 10 3 5.23 × 10 2 5.23 × 10 2 1.14 × 10 11 9.00 × 10 -6 Fluid, water (w) 9.99 × 10 2 4.43 × 10 3 4.17 × 10 3 2.15 × 10 9 3.03 × 10 -4

TABLE I .

 I Parameters describing the physical behavior of the materials constituting the thermophone systems investigated.

TABLE II .

 II Other parameters describing the physical behavior of the materials constituting the thermophone systems investigated.

Full cylindrical or spherical thermophoneIn Fig.9, panel a) and panel b), we show the frequency response of a full cylinder and a full sphere, respectively, by considering a variable radius R s .In Fig.9, panel a), reducing the radius of the cylinder increases the high frequency SPL. In fact, a 20dB/dec profile is observed from 100Hz to almost 1MHz for the smaller radius. This increase could once again be attributed to the smaller HCPUA induced by the reduced radius. On the other hand, having a large radius seems to create a more flat frequency response in this frequency range, which reminds us of the profile observed for plane waves.[START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] This is explained by the fact that if the cylinder radius becomes infinitely large, the cylinder can be approximated by a plane. Therefore, in this case, cylindrical and plane waves become similar. It is also seen that the first observable antiresonance is at a lower frequency for thicker thermophones. For the smallest radius R s = 0.5µm, the antiresonance is even higher than 100MHz (i.e., out of the range of the plot). This behavior was expected since the wavelengths of those mechanical resonances are directly proportional to the size of the thermophone. This phenomenon was also observed for plane waves.[START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] Lastly, it is important to remember that those plots are normalized with total power. In other words, 1W was applied to all simulated cylinder but for smaller cylinder the actual power density was higher due to the reduced volume.In Fig.9, panel b), we observe the spectra of full spheres with variable radius. We can make the same analysis developed for Fig.9, panel a), since results display a very similar behavior. In order to facilitate the comparison, as previously declared, the same total power and the same external surface have been imposed for both cylindrical and spherical geometries. We observe similar SPL levels above 1MHz for both shapes. Conversely, below 1MHz, the frequency spectrum profile of the sphere display lower SPL than the cylinder. This could be explained by the fact that the surface of the cylinder This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

Moreover, for both the cylindrical and the spherical geometry, we considered layers that can be either fluid or solid. It is interesting to remark that, starting from the classical balance equations of the continuum mechanics for fluids and solids, we obtained a unified formalism for describing cylindrical and spherical thermo-visco-acoustic waves. This unified picture has been exploited to elaborate the multilayer model, which has been implemented for an arbitrary number of layer and for an arbitrary disposition of fluid or solid phases. This approach allows therefore the analysis and the design of complex thermophone structures with cylindrical or spherical geometry. This contribution is complementary to what already exists in the literature with regard to stratified one-dimensional thermophone structures.[START_REF] Guiraud | Multilayer modeling of thermoacoustic sound generation for thermophone analysis and design[END_REF] From the computational point of view, the implementation of the model is very efficient since it is based on matrices (see Sections II C and III C) and the calculation time is of the order of a few seconds for systems with This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

ACKNOWLEDGMENT

This work was supported by the French DGA/MRIS (Direction Générale de l'Armement/Mission pour la Recherche et l'Innovation Scientifique), by the ANRT (Association Nationale Recherche Technologie) and is part of a project developed in collaboration with CINTRA, UMI 3288 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

Appendix A: Fourth order differential equation in cylindrical coordinates

The four solutions of the algebraic equation aθ [START_REF] Arnold | The Thermophone as a Precision Source of Sound[END_REF] bθ 2c = 0, associated to Eq.( 16), are written in the form ±θ 1 and ±θ 2 . The fourth order differential equation given in Eq.( 16) can so be written as

with D defined as in Eq.( 13), b = a(θ 2 1 +θ 2 2 ) and c = -aθ 2 1 θ 2 2 . We define Ψ as

Leaving aside the factor a, Eq.(A1) becomes

whose solution is

with A, B integration constants and H

i (z) and H

i (z) being the cylindrical Hankel function of order i and of first and second kind. The equation for T becomes

This non homogeneous differential equation gives the general solution of the temperature variation as

with ps being a particular solution of Eq.(A5) and C, D other two integration constants. In order to find the final solution of T , the particular solution of Eq.(A5) can be determined by finding separate solutions for the term AH

(1) 0 (iθ 1 r) and BH

(2) 0 (iθ 1 r), which then can be added to form the final particular solution.

We are looking for a particular solution of the form

The substitution of Eq.(A7) into the second order inhomogeneous equation of this kind leads to a system of equations of the form da dr H

(1)

Based on standard properties of Hankel functions, [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] we have

Eq.(A8) can be rewritten in matrix form as

.

(A10) By calculating the inverse matrix and by using the Wronskian relation [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]

H

(1)

we found the solutions

0 (iθ 1 r)H

(2)

0 (iθ 1 r)H

(1)

These integrals can be performed by means of the following property [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] zC ν (gz

with C ν and D ν being any cylinder function of order ν, not necessarily distinct, and g, h being different constants. This leads to

(A15) Using Eq.(A7) and the Wronskian relation in Eq.(A11), we have the first particular solution T 1

We now investigate the second part of the particular solution.

• Case 2:

0 (iθ 1 r). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0039458

We are looking for a particular solution of the form

Beside the change of the constant from A to B and the fact that we have now a second kind Hankel function in the inhomogeneous term, the calculation can be repeated in the same way.

We eventually find

A final solution of Eq.( 14) for T is finally found as

0 (iθ ac, f r) +CH ( 1)

with

Appendix B: Wavenumbers calculation

Wavenumbers in fluid media

In order to find simplified explicit expressions for θ ac, f and θ th, f , solutions of the algebraic equation associated to Eq.( 14), we make the assumption of weak conduction and weak viscosity. Hence, the asymptotic solutions of ax 4 + bx 2 + c = 0 for small values of a are investigated. If a = 0, the first couple of solutions is given by x 2 = -c b , which means x = ± -c b . So, we search ε such that x = ± -c b + εa, in order to obtain the solutions to the first order in the parameter a. These solutions correspond to a regular perturbation.

However, the second couple of solutions diverges to infinity when a → 0, thus corresponding to a singular perturbation.