Stéphanie Chaillat

Huy Duong Bui

Resolution of linear viscoelastic equations in the frequency domain using real Helmholtz boundary integral equations

Keywords: Computational solid mechanics, Viscoelasticity, Boundary integral equations Mécanique des solides numérique, Viscoélasticité, Équations intégrales de frontière

Boundary integral equations are well suitable for the analysis of seismic waves propagation in unbounded domains. Formulations in elastodynamics are well developed. In contrast, for the dynamic analysis of viscoelastic media, there are very seldom formulations by boundary integral equations. In this Note, we propose a new and simple formulation of time harmonic viscoelasticity with the Zener model, which reduces to classical elastodynamics if a compatibility condition is satisfied by boundary conditions. Intermediate variables which satisfy the classical elastodynamic equations are introduced. It makes it possible to utilize existing numerical tools of time harmonic elastodynamics.

Résumé

Résolution des équations de la viscoélasticité linéaire en régime harmonique à l'aide des équations intégrales de frontière de Helmholtz réelles. Les équations intégrales de frontière sont bien adaptées à l'étude de la propagation des ondes sismiques notamment pour les domaines non bornés. Les formulations pour l'élastodynamique sont très développées. Par contre, pour la dynamique des corps viscoélastiques modélisant le sol, peu de formulations par équations intégrales de frontière existent. On propose dans cette Note une nouvelle formulation simple de la viscoélasticité en domaine fréquentiel par un modèle de Zener qui se ramène à la formulation classique de l'élastodynamique, si les conditions à la frontière respectent une condition de compatibilité. Des variables intermédiaires qui respectent les équations de l'élastodynamique classique sont introduites. Il est alors possible de réutiliser les outils numériques existant pour l'élastodynamique harmonique. Pour citer cet article : S. Chaillat, H.D. Bui C. R. Mecanique 335 (2007).

Version française abrégée

Les équations intégrales de frontière sont bien adaptées à l'étude de la propagation des ondes sismiques dans les domaines non bornés. Les formulations en élastodynamique sont bien développées [START_REF] Semblat | Seismic site effects in deep alluvial basin : numerical analysis by boundary element method[END_REF]. D'un autre côté, il y a très peu de formulations équations intégrales de frontière pour l'étude dynamique des milieux viscoélastiques. La méthode classique consiste à utiliser des coefficients de Lamé complexes [START_REF] Dominguez | Boundary Elements in Dynamics International Series on Computational Engineering[END_REF]. Dans cette Note, on propose une nouvelle méthode simple pour la viscoélasticité en domaine fréquentiel, qui se ramène à l'élastodynamique classique par un simple changement de variables. Le principal avantage de cette méthode est de conserver des constantes élastiques réelles et donc de permettre de réutiliser les outils numériques existant pour l'élastodynamique en domain fréquentiel [START_REF] Bonnet | Boundary integral equation method for solids and fluids[END_REF].

Le modèle rhéologique utilisé dans cette Note est un modèle de Zener (Fig. 1). L'idée principale est de chercher le déplacement sous la forme découplée mais restrictive [START_REF] Semblat | Seismic site effects in deep alluvial basin : numerical analysis by boundary element method[END_REF]. Ensuite, les notations introduites dans [START_REF] Goriacheva | Contact problem of rolling of a viscoelastic cylinder on a base of the same material[END_REF] sont utilisées [START_REF] Schanz | Wave Propagation in Viscoelastic and Poroelastic Continua[END_REF]. Ces variables intermédiaires vérifient les lois de l'élasticité. En introduisant l'angle ψ tel que tan ψ = αω, on montre (3) que les variables u * et u sont déphasées de ψ mais ont la même pulsation ω.

On note de la même manière l'angle φ tel que tan φ = βω. En utilisant le fait que les variables u * et σ * doivent être en phase, on montre que u et σ doivent être déphasées de θ = ψφ (4). De plus, pour que u * et σ * puissent satisfaire l'équation de la dynamique [START_REF] Gaul | Dynamics of viscoelastic solides treated by boundary element approaches in time domain[END_REF], θ doit être petit. Ceci revient à ajouter l'hypothèse ω ≪ ω 1 = 1/|α -β|.

Supposons que les données u d sur ∂Ω 1 et σ d .n sur ∂Ω 2 soient compatibles sur la frontière (afin de pouvoir noter u sous la forme (1), donc déphasées de θ = ψφ). Par exemple, les données suivantes sont compatibles :

(i) u d = 0 sur ∂Ω 1 et σ d .n = 0 sur ∂Ω 2 (ii) u d = 0 sur ∂Ω 1 et σ d .n = 0 sur ∂Ω 2
Ainsi les variables u * et σ * qui s'en déduisent sont en phase. La solution du problème visco-élastique linéaire est alors obtenue en résolvant l'équation intégrale classique de l'élastodynamique pour u * et σ * . Connaissant la solution u * (resp. σ * ), on peut calculer sans difficulté u (resp. σ). En effet, les solutions du problème élastodynamique u * et σ * ont respectivement une amplitude égale à v(x) cos ψ et w(x) cos φ (d'après (3) et ( 4)). Pour connaître l'amplitude de u (resp. σ), il faut donc multiplier l'amplitude de u * (resp. σ * ) par cos ψ (resp. cos φ) où ψ = tan -1 αω et φ = tan -1 βω.

En pratique, la résolution numérique de cette formulation peut être effectuée en utilisant les travaux effectués en élastodynamique [START_REF] Chaillat | A Fast Multipole Method formulation for 3D elastodynamics in the frequency domain C[END_REF].

Introduction

The main advantage of boundary element method (BEM) is that only the domain boundary is discretized. As a result, the method is well suitable for the study of problems in unbounded domains. So, the boundary integral formulation of linear elasticity is used to study seismic wave propagation [START_REF] Semblat | Seismic site effects in deep alluvial basin : numerical analysis by boundary element method[END_REF].

But, the ideal model of a linear elastic soil is not adapted in a lot of cases. It is necessary to take into account the soil damping factor and so to use a formulation for dynamic viscoelasticity.

In time domain, various methods have been proposed for BEM formulation of viscoelasticity [START_REF] Schanz | Wave Propagation in Viscoelastic and Poroelastic Continua[END_REF]. They are sorted into three kinds. The first formulation is developed by applying the elastic-viscoelastic cor-respondence principle. An integral transform (according to time) is performed on the boundary integral equation of elastodynamics. Generally, Laplace transform is considered, for example in the works of Rizzo [START_REF] Rizzo | An application of the correspondence principle of linear viscoelasticity theory[END_REF] or Kusama [START_REF] Kusama | Boundary element method applied to linear viscoelastic analysis[END_REF]. Then, the viscoelastic fundamental solutions are obtained by applying the elasticviscoelastic correspondence principle to elastodynamic solutions. Various works deal with the reduction of the last step: the back transform to time domain.

The second class of methods uses the fundamental solutions of elastodynamics in time domain. Once the convolution with time shape functions is integrated analytically, the equation is transformed in Laplace domain where the correspondence principle is applied. Then, a numerical inverse transformation is required [START_REF] Gaul | Dynamics of viscoelastic solides treated by boundary element approaches in time domain[END_REF] to lead to a time domain boundary element formulation. Various works use the "convolution quadrature method" developed by Lubich [START_REF] Lubich | Convolution quadrature and discretized operational calculus[END_REF] to evaluate the convolution [START_REF] Schanz | A new visco-and elastodynamic time domain boundary element formulation[END_REF].

The last class of methods directly requires the knowledge of the viscoelastic fundamental solutions in time domain. Using differential systems of Kelvin and Boltzmann models, Mesquita [START_REF] Mesquita | A simple Kelvin and Boltzmann viscoelastic analysis of three-dimensional solids by boundary element method Engineering analysis with boundary elements[END_REF] determines integral formulations adapted to each model. With those methods, only for the simplest viscoelastic models, the fundamental solutions are available analytically and one does not consider works in elastodynamics.

In frequency domain, the usual method is to replace the Lamé's constants by complex values [START_REF] Dominguez | Boundary Elements in Dynamics International Series on Computational Engineering[END_REF].

In this Note, a simple method to formulate the boundary integral equations for viscoelasticity, with a Zener model, is presented. This method, based on the introduction of new intermediate variables, reuses the classical formulation of elastodynamics and presents the advantage to keep real values of Lamé's constants. To the authors knowledge, a formulation similar to ours does not exist in the literature.

Rheological model

Various rheological models exist to model the viscoelastic behavior of a material [START_REF] Flügge | Viscoelasticity[END_REF] In this model, if ǫ = 0, so the dashpot has not effect, it is called a "delayed modulus". The model is equivalent to two springs connected in series. On the contrary, if ǫ = ∞, the dashpot does not have time to react, it is called an "instantaneous modulus".

Formulation of the problem

In the following, usual typeface letters denote scalar quantities while boldface letters denote vectors, matrices or tensors. The partial derivate is denoted using a comma ( ∂f ∂x = f ,x ).

The study is made in the frequency domain. The main idea is to not define the displacement in the classical form u(x, t) = u(x)e iωt (ω denoting the circular frequency) but under the restictive condition:

u(x, t) = v(x) cos ωt (1)
where v(x) is a real. As a result, the variables x and t are uncoupled. The three-dimensional generalization of the Zener constitutive law has been proposed by I. Goriacheva in [START_REF] Goriacheva | Contact problem of rolling of a viscoelastic cylinder on a base of the same material[END_REF] and is used in the present Note, denoting σ the stress tensor and ǫ the strain tensor:

           σ * = σ + β ∂σ ∂t ǫ * = ǫ + α ∂ǫ ∂t u * = u + α ∂u ∂t (2)
The coefficients α and β (α > β) are determined by:

α = η k 1 , β = η k 0 + k 1
The tensors σ * and ǫ * are linked by the constitutive equation of isotropic linear elasticity, with the Lamé's coefficients λ and µ of the delayed modulus, σ * = Lǫ * .

With notation (1), we note that u(x, t) and u(x, t) are dephased of π/2. Using notation (2), we obtain u * (x, t) = v(x)[cos ωt-αω sin ωt]. Noting the angle ψ such that tan ψ = αω (0 ≤ ψ < π/2, i.e. cos ψ = 0), it follows:

u * (x, t) = v(x) cos ψ cos(ωt + ψ) (3)
It is now clear that the variables u * and u are dephased of ψ but have the same circular frequency ω.

Then, noting σ(x, t) = w(x) cos(ωt + θ), we obtain σ * (x, t) = w(x)[cos(ωt + θ)βω sin(ωt + θ)]. Defining in the same way that for u, the angle φ such that tan φ = βω (0 ≤ φ < π/2, i.e. cos φ = 0), it follows:

σ * (x, t) = w(x) cos φ cos(ωt + φ + θ) (4)
But, the variables σ * and u * are known to satisfy the linear elasticity equations, as a result they have to be in phase. It follows that ψ = φ + θ. Finally, σ and u have to be dephased of θ = ψφ.

Boundary integral formulation

The boundary integral equation method for this formulation of viscoelasticity is now defined using the well-known method for elastodynamics. In fact, the main advantage of this formulation is that only a simple change of variables is introduced.

Boundary integral equation. The quantities σ * and u * which are linked by the elastic law, are now shown to satisfy the dynamic equation. σ and u can satisfy the dynamic equation divσρü ∼ = 0 if and only if σ and u are almost in phase, that is to say if and only if the angle θ is small. The difference introduced by the dephasing between σ and u, in the dynamic equation, is proportional to ρω 2 θ. It can be easily proved that θ ∼ = |α -β|ω. As a result, ρω 2 θ is proportional to ρω 3 |α -β| and we remark that we have divσρü ∼ = 0 to order O(ω 3 ). It follows that:

divσ * -ρü * ≡ (divσ -ρü) + β(div σ -ρ ... u) + (β -α)ρ ... u ∼ = (β -α)ρ ... u (5)
Using the definition of u (1), we obtain that divσ *ρü * ∼ = 0 and that the difference introduced in this dynamic equation is equally proportional to ρω 3 |α -β| if and only if:

|α -β|ω ≪ 1 that is to say |θ| ≪ 1. (6)
As a result, if ω is much less than the limit frequency ω 1 = 1/|α -β|, the notation (1) is compatible with elastodynamics.

Since |α -β| is proportional to the viscosity coefficient η (for a given set of elastic constant k 0 and k 1 ), the lower the coefficient η, the higher the limit frequency ω 1 (in soil mechanics η is small so that ω 1 is very large). Thus the hypothesis of "low" frequency ω ≪ ω 1 (including the quasi-static case ω = 0), which we suppose in the following, is satisfed in soil mechanics, for a large frequency range. This leads to the boundary integral formulation which is the same that for elastodynamics:

divσ * + ρω 2 u * = 0, σ * = Lǫ * .
This formulation is now recalled [START_REF] Bonnet | Boundary integral equation method for solids and fluids[END_REF]. The stress vector T n u on a plane of normal n is defined by the operator (λ and µ representing the Lamé's constants):

T n = 2µ ∂ ∂n + λn.div + µn ∧ rot.
Noting Ω the region of space occupied by an elastic solid with isotropic constitutive properties, the displacement u at an interior point x ∈ Ω is given by:

u k (x) = ∂Ω [(T n u(y)) i U k i (x, y; ω) -u i (y)(T n U k (x, y; ω)) i ]dS y (x ∈ Ω), (7)
where U k i (x, y; ω) denotes the i-th component of the elastodynamic fundamental solution, in the frequency domain, for an infinite space.

When x ∈ ∂Ω, a singularity occurs in y = x. With the help of a well-documented limiting process, equation ( 7) yields the integral equation:

c ik (x)u i (x) = (P.V.) ∂Ω [(T n u(y)) i U k i (x, y; ω) -u i (y)(T n U k (x, y; ω)) i ]dS y (x ∈ ∂Ω), (8)
where (P.V.) indicates a Cauchy principal value (CPV) singular integral and the free term c ik (x) is equal to 0.5δ ik in the usual case where ∂Ω is smooth at x.

Boundary conditions. We consider a domain Ω, of boundary ∂Ω on which mixed but independent boundary conditions are imposed (∂Ω = ∂Ω 1 + ∂Ω 2 and ∂Ω 1 ∂Ω 2 = ∅). On ∂Ω 1 , the imposed displacement is noted as in (1):

u d (x, t) = u d (x) cos ωt,
on ∂Ω 2 , the stress vector is written in the same manner:

t d (x, t) = w d (x) cos ω(t + θ),
Assume that the data u d on ∂Ω 1 and σ d .n on ∂Ω 2 are compatible (if u can be defined in form (1) so the variables are dephased of θ = ψφ). For example, the following data are compatible:

(i) u d = 0 (circular frequency ω) on ∂Ω 1 and σ d .n = 0 on ∂Ω 2 (ii) u d = 0 on ∂Ω 1 and σ d .n = 0 (circular frequency ω) on ∂Ω 2 As a result, the intermediate variables u * and σ * , are necessary in phase. The problem in u * and σ * is solved using the well-known boundary integral formulation of elastodynamics in the frequency domain. Having the solution u * (resp. σ * ), u (resp. σ) is easily computed. Indeed, the solutions of the elastodynamic problem u * and σ * have respectively an amplitude equal to v(x) cos ψ and w(x) cos φ ((3) and ( 4)).

So, to compute the amplitude of u (resp.σ), one only has to multiply the amplitude of u * (resp.σ * ) by cos ψ (resp. cos φ) where ψ = tan -1 αω and φ = tan -1 βω.

Conclusion

A new and simple formulation of time harmonic viscoelasticity (including the quasi-static case) have been presented. Hence, it has been shown that this problem reduces to a classical elastic problem by a simple change of variables if the boundary conditions respect a restrictive condition. It makes possible to reuse existing numerical tools of time harmonic elastodynamics. The speed up of the computation can be done using the Fast multipole method [START_REF] Chaillat | A Fast Multipole Method formulation for 3D elastodynamics in the frequency domain C[END_REF].

1 k 0 Figure 1 .

 101 Figure 1. Zener model

Email addresses: chaillat@lms.polytechnique.fr (corresponding author) (Stéphanie Chaillat), hdb@lms.polytechnique.fr (Huy Duong Bui).


Zhongrui Chen

Gabriel Canard

Cloé Azarias

Denis Jacquemin

Olivier Siri

Straightforward Metal-Free Synthesis of an Azacalix[6]arene Forming a Host-Guest Complex with Fullerene C 60

Contrasting with the well mastered azacalix [4]arene, the synthesis and characterization of azacalix[6]arene derivatives remains a challenge for chemists. In this framework, bulky n-octylamino groups appended on the 1,3-metadiaminobenzene 2b can lead to the concomittant and straightforward one-pot formation of the azacalix[6]arene 4b and of the thermodynamically favoured azacalix[4]arene 3b through simple consecutive nucleophlic aromatic substitutions. The conformations of 3b and 4b were studied in solution and in the solid state using their molecular structures, their NMR and UV-visible spectra combined with first principle TD-DFT calculations. Preliminay solution experiments show that 4b can form a 1:1 host-guest complex with C60.

Introduction

Since their first complete characterization in the early 1970s, calixarenes have become the subject of an intense research and have been applied in various fields such as catalysis and specific delivery of drugs.

1 Among the different ways to tune the cavity sizes and shapes of these macrocycles, one can focus on the connection between aryl units by heteroatoms thereby forming heterocalixarenes.

2 Compared to thiacalixarenes incorporating sufur-bridging atoms that were intensively studied, 3 the research devoted to the syntheses and applications of an emerging class of heterocalixarenes, namely, azacalixarenes is still in its infancy. 4 This class of calixarenes are of peculiar interest because the connecting nitrogen atoms may bring supplementary functionalities to the macrocycle. For instance, the nitrogen atoms can form supramolecular bonds and their alkylation was used to produce derivatives with an improved solubility and/or featuring an inherent chirality. 5 Moreover, electron rich polyazacyclophanes can be used as macrocyclic scaffolds for hole-and spin-containing molecular-based electronics.

6

The synthesis of azacalixarenes relies on two principal strategies. The first one is the use of successive Buchwald-Hartwig palladium based aryl aminations that can produce azacalixarenes containing from three to more than ten building aryl or heteroaryl units. 4,7 Another catalyst-free synthetic preparation of such azacyclophanes can be achieved using successive Nucleophilic Aromatic Substitutions (SNAr). 4,8 This strategy was show to produce only azacalix [4]arenes when the starting reactants are simple aryl units connected through their meta-positions. 4,8 In this case, the experimental conditions are probably responsible of the exclusive formation of the thermodynamically favoured meta-linked azacalix [4]arenes. Interestingly, the importance of using a kinetic vs. thermodynamic control during successive SNAr was illustrated for oxacalixarenes when studying the competing formation of poly(m-phenylene oxides) vs. m-cyclooligomers starting from the condensation of meta-dihydroxyand metadihalido-aryl units 9 which often affords oxacalix [4]arenes rather selectively over other larger macrocycles. 10 Besides the use of non-equilibrating conditions, 11 the synthesis of large moxacalix[n]arenes (n > 4) by convenient SNAr strategies can be achieved following fragment coupling approaches 12 or by introducing bulky groups in the core of the cyclophanes.

11b,13

To the best of our knowledge, such metal-free syntheses from simple aryl units are not yet described when considering large azacalix[n]arenes (n > 4) where the expansion of the cavity could bring, for example, unique recognition properties 14 or a peculiar interaction with G-quadruplexes.

15

The present work illustrates how two peripheral octyl chains appended on a 2,4-positions in tetraminobenzene can lead to the formation of an azacalix [6]arene by following a SNAr strategy using simple aryl units. This unprecedented azacyclophane is easily separated by chromatography from the large amounts of the thermodynamically favoured four member macrocycle that is concomitantly produced. A structural investigation of their respective conformation is conducted in the solid state and in solution using single crystal X-ray diffraction, NMR studies in connection with (timedependent) density functional theory [(TD)-DFT] calculations.

+ NH 2 H 2 N NHR RHN F F NO 2 O 2 N 1 2a : R = H 2b : R = C 8 H 17 HN HN NH NH NH HN NO 2 O 2 N NHR RHN O 2 N NO 2 RHN NHR NO 2 NO 2 NHR NHR + H c H d H a H b H a H b H c H d
These first-principles calculations are combined with preliminary experiments to enlighten the impact of the macrocycle expansion towards the complexation of [60]fullerene.

Results and discussion

Concomitant formations of the azacalix[6]arene 4b and the azacalix[4]arene 3b.

A few years ago, we reported the synthesis of the tetraaminoazacalix[4]arene 3a produced by the SNAr between the electron-deficient aryl 1 and 1,2,4,5-tetraminobenzene 2a in MeCN in the presence of an excess of i Pr 2 NEt acting as a nonnucleophilic base (Scheme 1). 16 This one step strategy can be used to prepare N-substituted derivatives such as the tetra-noctylamino-azacalix[4]arene 3b starting from the N-substituted tetra-aminobenzene 2b even if a two step procedure was shown to be more efficient (Scheme 1). 17 The condensation of Assuming that this pre-organisation would be rather disfavoured if N-substituents are added to the nucleophile 2 and knowing that the production of larger macrocycles is usually increased in more diluted conditions, a first experiment was conducted within the same conditions but using a twice diluted concentration of 2b (Table 1, entry 1). If the cyclic tetramer 3b appeared to be still the major reaction product, we could isolate a significant amount of the hexamer 4b with a yield of 3.4% and with a simple chromatographic purification facilitated by the quite different polarities of the two macrocycles on silica (see the experimental section). The highest yield of 4b reached 4.5% when pursuing the dilution of the reaction mixture (Table 2, entry 2) while, as expected, more concentrated conditions leads to a higher 3b/4b ratio (Table 1, entries 3-4). We underline that highly ethanolconcentrated mixtures produce only linear oligomers and polymers because this solvent, by weakening the intramolecular hydrogen bonds, prevent the pre-organisation of oligomers towards the cyclisation (Table 1, entry 5). Solid state conformations of 3b and 4b.

Despite the growth of numerous single crystals from different concentrated solutions of 3b and 4b, their molecular structures were all suffering from partial disorder affecting the long alkyl moieties and the co-crystallized solvent molecules. Nevertheless, we are still confident that the corresponding structural solutions and data refinements are of high enough qualities to allow us a first qualitative examination of their respective conformations in the solid state.

The crystal structure of 3b features a typical 1,3-alternate conformation that is commonly observed in the structure of all the azacalix[4]arenes built on two equivalents of 1 (Figure 1). 8c,8f,16b This conformation, is also found by DFT as the most stable (vide infra), and is resulting from two stabilizing factors. The first one is the sp 2 character of the nitrogen bridging atoms that are conjugated with their adjacent electronwithdrawing dinitrobenzene rings illustrated by the short corresponding NC bon lengths averaging 1.35 Å. The second one is the formation of four strong and intramolecular N-H•••O 2 N hydrogen bonds (average distance of 2.03 Å) that involve the acidic and bridging N-H moieties (Figure 1). Interestingly, this molecular arrangement induces the location of the two aromatic protons H a of the dinitroaryl units (Scheme 1) in the anisotropic shielding cones of the two nearly parallel oriented (interplane angle of 4°) tetra-aminobenzene separated by 4.74 Å. Very few examples of single crystal X-ray structures of azacalix[6]arenes were reported previously. In 2006, Vale et al. described the structure of an aza[6]metacyclophane bearing six p-tert-butylphenyl groups in a 1,3,5-alternate conformation with an approximate S 6 point group of symmetry.

7d Two years later, an analogous derivative bearing six supplementary methoxy groups was shown to adopt, in the solid state, a definite 1,2,3-alternate conformation stabilized by intraannular hydrogen bonds. 18 In 2002 and 2005, Yamamoto et al.

reported the crystal structure of an azacalix[6]heteroarene adopting a triangular shape in which each of the six pyridine rings were alternatively inclined upward and downward to the macrocyclic plane due to the electrostatic repulsion between the lone pairs of the N-pyridine and of the N-bridging atoms.

7b,19

Two molecules featuring the same conformation are contained in the asymmetric unit of the crystal structure of 4b (Figure 2). At first glance, 4b adopts an 1, 2,4,5 alternate (u,d,d,u,d,d) conformation 1a,20 which no remarkable symmetry element. Because of the two different building blocks, each kind of fragments is found in the two aryl groups that are orientated upward. The molecular cavity is rendered inaccessible to solvent molecules by the upward tetra-amino benzene unit. This latter aryl group is probably folding inside the cavity following two weak CH•••N hydrogen bonds that cannot be precisely determined regarding the low resolution of the structure. As we show below, theoretical calculations deliver a very similar conformation as the most stable amongst those tested. Interestingly, if this conformation was kept frozen in solution, one would expect two different sets of signals for each kind of aromatic protons. For example, as in the structure of 3b, the aromatic proton H a (Scheme 1) of the upward dinitroaryl moiety is in the anisotropic shielding cones of the adjacent downward orientated and nearly parallel tetraaminobenzene moieties while the two other equivalent H a protons are only affected by one shielding cone.

Conformations of 3b and 4b in solution.

A first insight of the conformations of 3b and 4b in solutions is given by their 1 H NMR spectra recorded in CDCl 3 at room temperature (Figure 3). The azacalix[4]arene 3b has a single and symmetrical 1,3-alternate conformer in which the intraannular aromatic protons H a are located inside the anisotropic shielding cone of the two adjacent aryl rings. This particular effect is put into evidence by the corresponding singlet at a high field chemical shift of = 5.49 ppm. As stated above, this conformation is stabilized by four NH•••O 2 N hydrogen bonds that are illustrated by the down-field signal of the bridging NH signals at  = 8.85 ppm. The assignments of the NMR signals of 3b are also consistent with theoretical calculations (see the ESI). As for 3b, the 1 H HMR spectrum of the larger macrocycle 4b features only one signal for each aromatic and NH protons (Figure 3). This result underlines that 4b either adopts a symmetrical favourite conformation in solution or that this spectrum is resulting from different equivalent or non-equivalent structures that are rapidly interconverting on the NMR time scale.

The UV-Visible spectrum of a diluted dichloromethane solution of 3b (Figure 4) is characterized by an intense band in the UV region (340 nm) together with a shoulder lying at lower energy (400 nm). These bands originate from π-π* transitions of the electronically independent diamino-dinitro-benzene units as shown previously by TD-DFT calculations performed on similar derivatives.

8h A noticeable consequence of the macrocycle's expansion is evidenced by the corresponding absorption spectrum of 4b (Figure 4) displaying two bands, though having maxima located at equal energies, that have a different intensity ratio. Moreover, the ring expansion induces the broadening of the visible band that extends up to 600 nm. This effect can be attributed to an effective π-conjugation between two consecutive aryl units in 4b which is, in contrast, interrupted by the restricted and "frozen" 1,3 alternate conformation of 3b.

21

We have searched the conformational space for both 3b and 4b using DFT calculations. For the former, DFT predicts that the 1,3-alternate conformation is the most stable, the following one being at least more than 10 kcal.mol -1 less stable (see the ESI for details). Interestingly, despite the fact that a greater conformational flexibility would be expected for 4b, DFT provides exactly the same conclusion as for 3b: the conformer corresponding to the XRD structure is much more stable than the others. Together with the experimental NMR spectra, this indicates that there is probably an interconversion of 4b between equivalent structures alike the one displayed in Figure 2. This analysis is also supported by a comparison between theoretical and experimental NMR chemical shifts (see the ESI). When comparing the experimental and theoretical UV/Vis spectra of 3b (Figure S4 in the ESI,) one observes that, despite a blue-shift of the bands (of 0.25 eV and 0.27 eV for the shortand long-wavelength band, respectively), theory provides a shape that matches well the one of Figure 4. Indeed, though the relative intensity of the long-wavelength maxima is somehow underestimated, the theoretical separation between the two bands (0.55 eV) is in excellent agreement with the experimental value of 0.58 eV. The orbital composition of each relevant electronic transition, as well as the orbital shapes of 3b can be found in the ESI. One can conclude from these data that the excitations responsible for both absorption bands of 3b mainly imply molecular orbitals centered on the phenyl rings bearing the nitro groups, though there is small contribution from the phenyl rings bearing the amino groups for the longest wavelength band. In Figures S6, we compare the experimental and theoretical spectra of 4b. Even though the shape of the theoretical convoluted spectra does not perfectly the experimental curve, theory provides accurate trend Indeed, when going from 3b to 4b, one notes: i) no important changes of the intense absorption at ca. theoretical and measured spectra; ii) the spectrum of extends to larger wavelengths which also fits observation. Interestingly by analysing, the additional long-wavelength band specific to notices that it corresponds to a charge-transfer band, as the key occupied MOs are localized on the amino (donor) whereas their unoccupied counterparts are c on the nitro-phenyl moieties (see Table S7 andFigure Previous works have focused on the recognition properties of large azacalix[n]heteroarenes (n=5-10) towards guests such as fullerene derivatives.

14 Consequently, we ran preliminary solution experiments to study if host-guest complexes are formed between C 60 and the azacalixarenes tetrachloroethane solutions at a millimolar scale (0.4 mM), no interaction occurs between 3b and increasing amounts of C as illustrated by the corresponding UV-Visible spectra that result from a simple addition of the absorption spectra of each component (Figure S8 in the SI). On the contrary, when the same experiment was conducted using 4b, a clear modification of the shape and intensity of the visible band was detected and reveal the probable formation of a 1:1 com through the presence of two isosbestic points (Figure S

We have performed DFT calculations on the possible complexes between C 60 and both 3b and 4b. detailed in the ESI. For 3b-C 60 the most stable aggregate that could be identified presents an interaction free energy as small as ca. -3 kcal.mol -1

. In contrast, for 4b-C 60 , we could found two complexes presenting interaction energies of compared to the corresponding 4b free conformer 5 and S16). In these two compounds, the C with the nitro-phenyl moieties. Interestingly, we note that this interaction takes place with two rings for the most stable "free" conformer found in XRD (left-hand-si three rings for a less stable conformer (right Figure 5). For this latter case, the gain of interaction energy is however too limited to compensate the less stable character of the free conformer (10 kcal.mol -1 , see Table the former is foreseen to be present in solution.

For the most stable 4b-C 60 aggregate, we have determined the UV/Vis spectra with TD-DFT, using CAM functional known to provide a sound description of long charge-transfer excited-states. Compared to the free TD-DFT calculations revealed several additional transitions of weak intensity in the 450-500 nm domain, which fits the experimental findings (Figure S9). Amongst all these TD , we compare the experimental and . Even though the shape of the perfectly match with the experimental curve, theory provides accurate trends.

one notes: i) no important ca. 300 nm in both the spectrum of 4b also fits the experimental analysing, the nature of this wavelength band specific to 4b, one clearly transfer band, as the key occupied MOs are localized on the amino-phenyl rings counterparts are centered phenyl moieties (see Table S7 and Figure S7). Such , due to the orthogonal nature Previous works have focused on the recognition properties 10) towards guests such Consequently, we ran preliminary guest complexes are and the azacalixarenes 3b and 4b. In chloroethane solutions at a millimolar scale (0.4 mM), no and increasing amounts of C 60 Visible spectra that he absorption spectra of each in the SI). On the contrary, when the , a clear modification of the shape and intensity of the visible band was detected and reveal the probable formation of a 1:1 complexation through the presence of two isosbestic points (Figure S9).

We have performed DFT calculations on the possible 4b. The results are the most stable aggregate that could be identified presents an interaction free energy as small , we could found two complexes presenting interaction energies of ca. -9 kcal.mol . Interestingly, we note that this o rings for the most stable side in Figure 5) or for a less stable conformer (right-hand-side in For this latter case, the gain of interaction energy is however too limited to compensate the less stable character , see Table S2), and only the former is foreseen to be present in solution.

aggregate, we have determined DFT, using CAM-B3LYP, a functional known to provide a sound description of long-range states. Compared to the free 4b, these DFT calculations revealed several additional transitions of 500 nm domain, which fits the . Amongst all these TD-DFT Compound 2b was prepared according to literature procedure. 17 All reagents and solvents anhydrate were purchased from Alfa Aesar or Sigma Aldrich and were used without further purification. Column chromatography was performed using silica gel (60-120 mesh). Analytica chromatography (TLC) was performed on precoated silica gel 60 F254 (0.5 mm) aluminium plate. Visualization of the spots on TLC plates was achieved by exposure to UV light. NMR spectra were recorded on a Bruker AC250 ECS400 NMR spectrometer. Chemical shifts are reported in delta (δ) units, expressed in parts per million using the residual protonated solvent signal as an internal standard (For proton: CDCl 3 , 7.26 ppm; Acetone-d 6 , 2.05 ppm. For ppm; Acetone-d 6 , 30.8 ppm. The multiplicity of designated by the following abbreviations: s, singulet; br s, broad singulet; d, doublet; t, triplet; br t, broad triplet; m, multiplet. Coupling constants, J High-resolution mass spectrometry (HRMS) performed on a SYNAPT G2 HDMS (Waters). Absorption spectra were recorded using 1 mm vertical transitions, the most intense takes place at 474 nm (f=0.02) and the largest MO contribution to it is an HOMO to se orbitals are displayed in Figure 6 and it is crystal clear that they imply a charge-transfer from the (HOMO) to the C 60 (LUMO).

complexes identified by DFT calculations.

ight) of the most stable 4b-C60 aggregate.

was prepared according to literature All reagents and solvents anhydrate were fa Aesar or Sigma Aldrich and were used without further purification. Column chromatography was 120 mesh). Analytical thin layer chromatography (TLC) was performed on precoated silica gel-60 F254 (0.5 mm) aluminium plate. Visualization of the spots on TLC plates was achieved by exposure to UV light. 1 H and 13 C NMR spectra were recorded on a Bruker AC250 or on a JEOL . Chemical shifts are reported in delta (δ) units, expressed in parts per million using the residual protonated solvent signal as an internal standard (For proton: , 2.05 ppm. For 13 C: CDCl 3 , 77.0 , 30.8 ppm. The multiplicity of the signals is designated by the following abbreviations: s, singulet; br s, broad singulet ; d, doublet; t, triplet; br t, broad triplet; m, multiplet. Coupling constants, J, are reported in Hertz (Hz).

tion mass spectrometry (HRMS) analysis was performed on a SYNAPT G2 HDMS (Waters). Absorption 1 mm and 10 mm quartz cells in a

Varian Cary 50 UV-vis spectrophotometer and spectroscopic grade solvents.

Synthesis of compounds 3b and 4b. To a solution of 1,5difluoro-2,4-dinitrobenzene 1 (m= 135 mg, 0.66 mmol, 0.95 equiv.) in anhydrous MeCN (7 mL), was added 2b (m= 300 mg, 0.69 mmol, 1.0 equiv.). The flask was closed by a septum, and cooled in a water-ice bath. The system was degassed 3 times by pump-argon cycling. Under argon degassed DIPEA (720 µL, 4.13 mmol, 6.0 equiv.) was added dropwise by a syringe. After stirring at 0 °C for 3 hrs, the solution was stirred at room temperature for 5 days. The resulting solid in suspension was filtered and rinsed with acetone (100 mL) affording 3b (m= 127 mg, 121 µmol) as an orange powder. J HH = 7.0 Hz, 12H), 1.71-1.17 (m, 72H), 0.88 (t, 3 J HH = 6.8 Hz, 18H) ppm. 13C NMR (63 MHz, CDCl 3 ): δ 146. 7, 143.7, 129.5, 125.1, 121.6, 112.4, 95.7, 95.0, 44.2, 31.8, 29.5, 29.4, 29.2, 27.1, 22.6, 14.0 Crystallography.

For 3b intensity data of single crystal X-ray diffraction were collected on a Bruker-Nonius KappaCCD diffractometer using MoKα radiation (λ = 0.71073 Å). Data collection was performed with COLLECT, 22 cell refinement and data reduction with DENZO/SCALEPACK. [START_REF] Otwinowski | Macromolecular Crystallography[END_REF] For 4b intensity data of single crystal X-ray diffraction were collected with a Rigaku Oxford Diffraction SperNova diffractometer using Cu-Kα radiation (λ = 1.54184 Å). The structures of 3b and 4b were solved with SIR92 [START_REF] Altomare | [END_REF] and SHELXS 25 respectively while SHELXL2013 25 used for full matrix least squares refinement. The hydrogen atoms were introduced at geometrical positions and their U iso parameters were fixed to 1.2U eq (parent atom) for the aromatic or CH 2 carbon atoms and to 1.5U eq (parent atom) for the methyl groups. CCDC-1525759 (3b) and CCDC-1525760 (4b) contain the supplementary crystallographic data. These data can be obtained free of charge from The Cambridge Crystallographic Data Centre via www.ccdc.cam.ac.uk/data_request:cif.

Computational details.

All (TD-)DFT calculations have been performed using the Gaussian 09.D01 program. [START_REF] Frisch | Gaussian 09 Revision D.01[END_REF] We have tightened both the selfconsistent field (10 -10 a.u.) and geometry optimization (10 -5 a.u.) convergence thresholds and use an accurate (99,590) pruned DFT integration grid. The macrocycles have been modeled by replacing the C 8 H 17 alkyl chains with Methyl groups. All our (TD-)DFT calculations relied on the PBE0 hybrid functional 27 using empirical dispersion corrections, as described by the so-called D3-BJ model. 28 CAM-B3LYP 29 was used to model the spectrum of the 4b-C 60 complex as this functional is adequate for long-range charge-transfer phenomena. Following a basis set combination approach proposed previously, 30 we used the 6-31G(d) atomic basis set for determining the geometrical and vibrational parameters whereas the TD-DFT transition energies have been computed with 6-31+G (d,p) [6-31G(d) for the 4b-C 60 complex]. In computing the conformation and interaction energies, single point calculations using the highly extended def2-QZVP basis set have been performed as well so to avoid the need of correcting for BSSE effects. The nature of all stationary points was confirmed by analytical Hessian calculations that returned 0 (minima) imaginary vibrational modes. These calculations gave access to estimates of the free energies (G). The NMR spectra were computed using the cc-pVTZ basis set. We note that the TMS reference signal was computed with the same protocol. Environmental effects [here, (CHCl 2 ) 2 for transitions and chloroform for NMR] have been accounted for electronic transition calculations, using polarizable continuum model (PCM), 31 as implemented in Gaussian09, [START_REF] Frisch | Gaussian 09 Revision D.01[END_REF] applying the linearresponse regime for excited-states. During the calculation of complexation and conformational energies, the solvent has also been introduced for single point calculation on the gas phase geometry using the PCM or SMD model.

32

Conclusions

Successive metal-free nucleophilic aromatic substitutions between an activated meta-dihalobenzene and a metadiaminobenzene bearing two n-octylamino groups can lead to the concomitant and straightforward one-pot formation of the azacalix[6]arene 4b and of the thermodynamically favoured azacalix[4]arene 3b. Multiple analyses including single crystal X-ray diffraction, NMR and UV-Visible spectroscopy were combined with by TD-DFT calculations and put into evidence that both macrocycles feature a single conformer in solution which is the one observed in the crystalline solid state. If the corresponding 1,3-alternate conformer of the azacalix[4]arene 3b was expected, the 1,2,4,5 alternate conformation adopted by the larger macrocycle 4b has been rarely observed previously. This conformation remains flexible in solution and allows the formation of stable host-guest complexes between the azacalix[6]arene 4b and C60 which barely interacts with the smaller macrocycle 3b.

1

  and 2a affords a single macrocycle 3a irrespective of the experimental conditions with a very high yield reaching 92% when the reactant concentrations in MeCN is ca. 0.01 M. The macrocyclisation of the [2+2] product is probably favoured by intramolecular N-H•••O 2 N hydrogen bonding interactions involving the NH bridges that pre-organise the intermediate oligomers towards the cyclisation.

  8a

Fig 1 .

 1 Fig 1. Two views of the single-crystal X-ray structure of 3b (N-octyl chains are omitted for clarity). The dotted lines show the intramolecular NH•••O2N hydrogen bonds.

Fig 2 .

 2 Fig 2. Two views of the single-crystal X-ray structure of 4b (N-octyl chains and solvent molecules are omitted for clarity, the hydrogen atoms locations were not fitted).

Fig 3 .

 3 Fig 3. Partial 1 H NMR spectra of 3b and 4b recorded in CDCl3 at room temperature (asterisks denote residual solvent peaks; the full spectra are given in the SI).

Fig 4 .

 4 Fig 4. UV-Visible spectra of diluted dichloromethane solutions of 3b (dashed line) and 4b (solid line).

  S transfer was not possible in 3b, due to the orthogonal nature of the phenyl planes in that compound.Complexation of C 60 by 3b and 4b.

  In these two compounds, the C 60 is interacting

Fig 5 .

 5 Fig 5. View of the two most stable 4b-C60 complexes identified by DFT calculations.

Fig 6 .

 6 Fig 6. HOMO (left) and LUMO (right) of the most stable

  Synthesis of azacalix[4]arenes 3a-b and of the azacalix[6]arene 4b.

		O 2 N	NO 2
	RHN	HN	NH	NHR
	RHN	HN	NH	NHR
		O 2 N	NO 2
		3a : R = H (92%)
		3b : R = C 8 H 17 (33-65%)

Scheme 1 i Pr 2 NEt Solvent 4a : R = H (0%) 4b : R = C 8 H 17 (0-5%)

Table 1

 1 Screening of reaction conditions (using 1 equivalent of 2b).

	Entry	1 a	Solvent	[2b]	i Pr2NEt b	Temp.	3b/4b
				(M)		(°C)	Yields (%)
	1	0.95	MeCN	0.02	6	i) O °C (2 h)	32.4/3.4
						ii) rt. (12 h)
						iii) reflux (5 h)
	2	0.95	MeCN	0.03	6	i) O °C (3 h)	44.4/4.5
						ii) rt. (120 h)
	3	0.77	MeCN	0.11	10	i) O °C (4 h)	60.5/2.6
						ii) rt. (72 h)

a Equivalents of 1. b Equivalents of base.

Acknowledgements

This work was supported by the Centre National de la Recherche Scientifique, the Ministère de la Recherche et des Nouvelles Technologies (PhD grant to Z. C.) and the Agence Nationale de la Recherche in the frame of the EMA project (PhD grant of C. A.). D.J. thanks the LUMOMAT project for support. This works used computational resources of the IDRIS/CINES and of the CCIPL.

Notes and references


Romain De Monchy

François Destrempes

Ratan K Saha

Guy Cloutier

Emilie Franceschini
email: franceschini@lma.cnrs-mrs.fr

Coherent and incoherent ultrasound backscatter from cell aggregates

Keywords: numbers: 43.80.Cs, 43.80.Qf, 43.35.Bf quantitative ultrasound, structure factor, aggregate, effective medium theory, ultrasound tissue characterization

I. INTRODUCTION

Ultrasound scattering from biological tissues exhibits a frequency dependence which is related to the structure of the insonified medium. For determining the tissue microstructure, one approach consists in fitting the measured frequency-dependent backscatter coefficient (BSC) to a theoretical BSC derived using an appropriate theoretical scattering model. Two theoretical scattering models are often used in the field of quantitative ultrasound (QUS) imaging. The first model, named the spherical Gaussian model, describes a tissue as a random inhomogeneous continuum with impedance fluctuations. 1 The second model considers tissues as randomly distributed discrete scatterers with an impedance differing from a homogeneous background medium. 2 In both models, the scatterers are assumed to be independently and uniformly randomly distributed, which may correspond to the case of a low scatterer concentration. However, biological tissues often have high cellular concentrations and complex structures that invalidate the assumption of those models. 3,[START_REF] Franceschini | Structure factor model for understanding the measured backscatter coefficients from concentrated cell pellet biophantoms[END_REF] For example, some tumors have densely packed cells and eventually aggregating cells (such as breast sarcoma [START_REF] Oelze | Examination of cancer in mouse models using high-frequency quantitative ultrasound[END_REF] ); and blood contains a high volume fraction of red blood cells (between 30 and 50%) that form aggregates having rouleaux structures in normal blood or clumps in some pathologies, as in diabetes mellitus. [START_REF] Schmid-Schonbein | Red cell aggregation in blood flow. I. New methods of quantification[END_REF][START_REF] Schmid-Schonbein | Erythrocyte aggregation: causes, consequences and methods of assessment[END_REF] Therefore, theoretical scattering models have been proposed to include high cellular concentration together with aggregating cells, and to estimate scatterer structures. [START_REF] Savery | A point process approach to assess the frequency dependence of ultrasound backscattering by aggregating red blood cells[END_REF][START_REF] Yu | Experimental ultrasound characterization of red blood cell aggregation using the structure factor size estimator[END_REF] A scattering model called the Effective Medium Theory combined with the Structure Factor Model (EMTSFM) was recently proposed. [START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF] In this model, the aggregates of cells are viewed as individual scatterers which have effective properties determined by the acoustical characteristics and concentration of cells within aggregates. The approximation of cell aggregates as homogeneous effective scatterers is combined with the structure factor model to consider the concentrated medium; i.e., to consider the interference effects caused by the correlations between the spatial positions of effective scatterers. The EMTSFM allows characterization of the radius and of the compactness of cell aggregates, as shown in a previous simulation study. [START_REF] Franceschini | Comparison of three scattering models for blood characterization[END_REF] The goal of the current study was to further develop the EMTSFM and to extend its validity into a larger frequency range. In that aim, the Effective Medium Theory (EMT) was modified to decompose the differential backscattering cross section of a single cell aggregate into coherent and incoherent components as proposed by Morse and Ingard. [START_REF] Morse | Theoretical Acoustics[END_REF] The coherent component corresponds to the average wave emerging from the effective scatterer, and the incoherent component represents the fluctuation of the scattering wave around its average within the effective scatterer. Note that only the coherent component was taken into account in our previous works. [START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF][START_REF] Franceschini | Comparison of three scattering models for blood characterization[END_REF] A new theoretical expression for the incoherent component based on the structure factor is proposed. This theoretical improvement was assessed using three dimensional (3D) computer simulations of ultrasound backscatter from aggregating cells. The new incoherent component based on the structure factor was compared to the incoherent component based on a Gaussian direct correlation function proposed by Morse and Ingard. [START_REF] Morse | Theoretical Acoustics[END_REF]

II. THEORY

In this section, the theoretical expression for the ultrasound backscattering response from one aggregate of cells is proposed. Then, the theoretical BSC expression for an ensemble of identical aggregates that was developed in Ref. 10 is briefly recalled.

A. Differential backscattering cross-section σag of one aggregate of cells

In the sequel, we consider identical cells distributed randomly inside an aggregate. It is assumed that the incident wavelength is large compared to the cell size. Consequently, the cell shape can be approximated by a sphere of radius a having an equivalent volume V c = 4 3 πa 3 . Cells are described in terms of their mass density ρ c and compressibility κ c , and their surrounding medium is characterized by its mass density ρ 0 and compressibility κ 0 . An aggregate of cells, denoted V ag , is assumed to be spherical as it occurs in some tissues (e.g., red blood cell aggregates in pathological cases [START_REF] Schmid-Schonbein | Red cell aggregation in blood flow. I. New methods of quantification[END_REF][START_REF] Schmid-Schonbein | Erythrocyte aggregation: causes, consequences and methods of assessment[END_REF] or breast sarcoma tumors [START_REF] Oelze | Examination of cancer in mouse models using high-frequency quantitative ultrasound[END_REF] ).

The EMT model assumes that an aggregate of cells can be treated as an effective homogeneous spherical scatterer of radius r ag . As previously studied by Morse and Ingard, [START_REF] Morse | Theoretical Acoustics[END_REF] the differential backscattering cross-section σ ag of an aggregate can be decomposed into two components:

σ ag = σ ag,coh + σ ag,inc , (1)
where the coherent component σ ag,coh corresponds to the average wave emerging from the effective scatterer, and the incoherent component σ ag,inc describes the fluctuation of the scattering wave around its average within the effective scatterer, as detailed below.

The scattering amplitude

Using Born and far-field approximations, the backscattered wave pressure is expressed in terms of the backscattered amplitude Φ ag ; i.e., p s (r) = e ikr /r Φ ag (k), with 13

Φ ag (k) = k 2 4π Vag γ κ (r 0 ) -γ ρ (r 0 ) e 2ik n0•r0 d 3 r 0 ,
where k is the wavenumber, r 0 is the position in the three dimensional space, n 0 is the incident wave direction, γ κ (r 0 ) = (κ(r 0 )-κ 0 )/κ 0 and γ ρ (r 0 ) = (ρ(r 0 )-ρ 0 )/ρ(r 0 ) are the fractional variations in compressibility and mass density, respectively. One has

γ κ (r 0 ) -γ ρ (r 0 ) = γ κ -γ ρ (if r 0 is inside the cells) 0 (if r 0 is outside the cells),
(2) where γ κ = (κ c -κ 0 )/κ 0 and γ ρ = (ρ c -ρ 0 )/ρ c . Assuming N identical cells inside an aggregate with centers located at positions r j , j = 1, ..., N , and considering the change of variable r 0 = r 0 -r j , the scattering amplitude can be expressed as

Φ ag (k) = k 2 (γ κ -γ ρ ) 4π V c V -1 c Va e 2ik n0•r 0 d 3 r 0 N j=1 e 2ik n0•rj ,
(3) where V a is a sphere of radius a centered at the origin (see Appendix A for details). The integral term in Eq.

(3) is equal to

V -1 c Va e 2ik n0•r 0 d 3 r 0 = 3 sin(2ka) -2ka cos(2ka) (2ka) 3 = F 0 (k; a), (4)
and F 0 (k; a) 2 = F (k; a) is the spherical form factor. [START_REF] Insana | Acoustic scattering theory applied to soft biological tissues[END_REF] From there, the scattering amplitude resulting from one aggregate is expressed as

Φ ag (k) = k 2 (γ κ -γ ρ ) 4π V c F 0 (k; a) N j=1 e 2ik n0•rj . ( 5
)

The coherent component

The coherent component of the differential backscattering cross-section of the effective scatterer is obtained from the coherent scattering amplitude as

σ ag,coh (k) = | Φ ag (k) | 2 ,
where • represents the expected value of a random variable. One computes:

Φ ag (k) = k 2 (γ κ -γ ρ ) 4π V c F 0 (k; a) N j=1 e 2ik n0•rj = k 2 (γ κ -γ ρ ) 4π φ i V ag F 0 (k; a) 1 N N j=1 e 2ik n0•rj , (6)
where V ag = 4 3 πr 3 ag is the volume of an aggregate of radius r ag , and φ i = N V c /V ag is the aggregate compactness. Next, we define

ρ ag = φ i ρ c + (1 -φ i )ρ 0 κ ag = [φ i /κ c + (1 -φ i )/κ 0 ] -1 (7)
which are obtained from the average values of 1/ρ and κ, respectively, within the aggregate, as mentioned in Eq. (8.2.23) of Morse & Ingard. [START_REF] Morse | Theoretical Acoustics[END_REF] Denoting γ κ,ag = (κ ag -κ 0 )/κ 0 and γ ρ,ag = (ρ ag -ρ 0 )/ρ ag , one computes (γ κ -γ ρ )φ i = γ κ,ag -γ ρ,ag . Based on these definitions, it follows that the coherent scattering amplitude of the aggregate is expressed as:

Φ ag (k) = k 2 (γ κ,ag -γ ρ,ag ) 4π V ag F 0 (k; a) 1 N N j=1
e 2ik n0•rj .

(8) Note that in this derivation, no approximation was used, but rather algebraic manipulations and definitions.

The expression F coh (k) = F 0 (k; a)1/N N j=1 e 2ik n0•rj corresponds to the Fourier transform of the function 1/N χ a (r -r j ) that describes a system of spheres, where χ a (r -r j ) = 1 inside a sphere V a of radius a centered at r j , and 0 elsewhere. Since it is not straightforward to compute the expression F coh (k), we approximate this quantity with the expression F 0 (k, r ) given in Eq. ( 4) for an equivalent full sphere of radius r . The full sphere choice was motivated here by the postulated spherical distribution of cells in aggregates. In order to match the low frequency approximations of the scattering amplitudes from the system of spheres and from the equivalent full sphere, their gyration radii R g should be equal. The gyration radius of a system of N spheres of radius a is equal to

R g = 3 5 a 2 + 1 N N j=1 |r j | 2 , (9)
where r j is the position vector of the jth cell with respect to the aggregate center, as given in Eq. ( 6) of Saha et al. [START_REF] Saha | Assessment of accuracy of the structure-factor-size-estimator method in determining red blood cell aggregate size from ultrasound spectral backscatter coefficient[END_REF] So, the radius r of the equivalent full sphere must be expressed as: r = 5 3 R g , as demonstrated in Appendix B. Henceforth, from Eq. ( 8), the coherent scattering amplitude resulting from one aggregate is approximated as

Φ ag (k) ≈ k 2 (γ κ,ag -γ ρ,ag ) 4π V ag F 0 (k; 5 3 R g ). ( 10
)
All together, the coherent component of the differential backscattering cross-section of the effective sphere is modeled as:

σ ag,coh (k) ≈ k 4 (γ κ,ag -γ ρ,ag ) 2 16π 2 V 2 ag F (k; 5 3 R g ). ( 11
)
One should note that in the former work, 10 the spherical form factor F given in Eq. ( 11), was computed by considering the external radius of the aggregate r ag (instead of 5/3R g in the present work) . Comparison between the former and the present coherent form factors will be presented in the subsection IV.A.

The incoherent component

Inside the aggregate, each cell produces a scattered wave as a sphere of radius a, with density ρ c and compressibility κ c , in a surrounding medium of density ρ ag and compressibility κ ag throughout the aggregate. Therefore, the wavenumber inside the aggregate is modified as [START_REF] Morse | Theoretical Acoustics[END_REF] k ag = k √ ρ ag κ ag / √ ρ 0 κ 0 . Based on Eq. ( 5), the incoherent component of the differential backscattering cross-section of the effective scatterer is equal to

σ ag,inc (k) = |Φ ag (k ag ) -Φ ag (k ag ) | 2 (12) = k 4 ag (γ κ -γ ρ ) 2 16π 2 V 2 c N F (k ag ; a)S 0 (k ag ; a, φ i ), (13)
where

S 0 (k; a, φ i ) = 1 N N j=1 e 2ik n0•rj - N j=1 e 2ik n0•rj 2 (14
) is the structure factor for a collection of N randomly distributed identical spheres of radius a and of concentration φ i within an aggregate. Assuming that the aggregate is sufficiently large, Eq. ( 14) corresponds to Eq. (10) in Ref. 15. This latter structure factor was analytically calculated as established by Wertheim, [START_REF] Wertheim | Exact solution of the Percus-Yevick integral equation for hard spheres[END_REF] which is denoted here S(k ag ; a, φ i ). One thus obtains the expression

σ ag,inc (k) ≈ k 4 ag (γ κ -γ ρ ) 2 16π 2 V ag V c φ i F (k ag ; a)S(k ag ; a, φ i ).
(15) In the sequel, we compute the differential backscattering cross-section σ ag = σ ag,coh + σ ag,inc of an aggregate using Eqs. (1), ( 11) and (15).

Let us mention that one can derive an alternative expression of the incoherent component of the differential backscattering cross-section based on a Gaussian direct correlation function σ ag,inc,G as proposed by Morse and Ingard. [START_REF] Morse | Theoretical Acoustics[END_REF] The proposed computation of σ ag,inc,G is a modified version of the incoherent component of Ref. 12 and is expressed as follows:

σ ag,inc,G (k) ≈ k 4 ag (γ κ -γ ρ ) 2 16π 2 V ag φ i (1-φ i ) 8 √ 2π 3/2 3 d 5 k 2 ag e -d 2 kag , ( 16
)
where d is the correlation distance related to the cell radius a as d = 2 1/6

3 1/3 π 1/6 a ≈ 0.643092×a, as proposed by Insana and Brown (see Eq. ( 81) p.107 in Ref. 13). Details can be found in Appendix C. The differential backscattering cross-section of an aggregate based on Eqs. (1), ( 11) and ( 16) is denoted σ ag,G in the sequel.

B. Backscatter coefficient BSC EMTSFM from an ensemble of identical aggregates

The BSC from an ensemble of identical aggregates is obtained by summing the contributions from individual effective scatterers and modeling the effective sphere interaction by a statistical mechanics structure factor S ag as follows: [START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF][START_REF] Franceschini | Comparison of three scattering models for blood characterization[END_REF]

BSC EMTSFM (k) = m ag σ ag (k)S ag (k; r ag , φ ag ), (17)
where m ag is the number density of effective spheres, which is related to the volume fraction of effective spheres φ ag as m ag = φ ag /V ag and σ ag is computed using Eqs.

(1), ( 11) and (15). In Eq. ( 17), S ag is the structure factor for a collection of randomly distributed identical effective spheres of radius r ag and of concentration φ ag , analytically computed as established by Wertheim. [START_REF] Wertheim | Exact solution of the Percus-Yevick integral equation for hard spheres[END_REF] The volume fraction of effective scatterers is equal to the volume fraction of cells in the tissue φ divided by the aggregate compactness φ i : φ ag = φ/φ i .

III. 3D SIMULATION METHOD

The 3D simulation study was conducted based on the Structure Factor Model (SFM). The SFM is an ultrasound scattering model largely used to predict the frequency dependence of the BSC from aggregated cells in the field of ultrasound tissue characterization. [START_REF] Franceschini | Structure factor model for understanding the measured backscatter coefficients from concentrated cell pellet biophantoms[END_REF][START_REF] Savery | A point process approach to assess the frequency dependence of ultrasound backscattering by aggregating red blood cells[END_REF][START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF][START_REF] Saha | Assessment of accuracy of the structure-factor-size-estimator method in determining red blood cell aggregate size from ultrasound spectral backscatter coefficient[END_REF][START_REF] Franceschini | Modeling of ultrasound backscattering by aggregating red blood cells[END_REF] The SFM consists of summing contributions from individual cells and modeling cell interaction by a statistical mechanics structure factor, whatever the complexity of the cells spatial distribution such as cell aggregates. The structure factor computation for a complex particle spatial distribution was described in details in Section 6.3.1 of Ref. 17. The current section focuses on the procedure for obtaining cell aggregates in the computer simulations, and then on the simulation of the differential backscattering cross section of one cell aggregate, denoted σ ag,sim , and of the BSC sim from aggregated cells based on the SFM.

A. Simulation of spatial distribution of cells within one aggregate

We describe here how the cell distribution inside one aggregate was computed. The cell particle radius a was set to 2.75 µm for all simulations, which corresponds to the red blood cell size usually used in blood computer simulations. [START_REF] Franceschini | Comparison of three scattering models for blood characterization[END_REF] We specified the aggregate radius r ag defined as the radius of the external envelope (Fig. 1) and the aggregate compactness φ i that fixes the number of cells N within the aggregate. N cells were uniformly randomly distributed such that cells within a radius of r ag could overlap, and the total number of overlapping pairs was counted. Then, the system was able to evolve by moving randomly selected cells to another position than the initial. If the number of overlapping pairs of the new system equaled or was lower than in the previous system, the displacement was accepted and the procedure was repeated until no overlapping was detected. This process allowed to reach aggregate compactness φ i up to 40% in 3D, whereas a method using a random sequential absorption [START_REF] Hinrichsen | Geometry of random sequential adsorption[END_REF] would give compactness near 30%. Figure 1 illustrates spatial arrangements of cells within a single aggregate for two aggregate compactnesses of 10% and 40% with an identical aggregate radius r ag /a = 7. The radius of gyration was computed using Eq. ( 9). For the aggregate example given in Fig. 1, a change of compactness from 10% to 40% induces a slight increase of the radius of gyration from R g = 12.56 µm to R g = 14.02 µm.

B. Computation of a simulated σ ag,sim curve

The simulated differential backscattering cross section of one cell aggregate σ ag,sim , was computed using the SFM [START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF][START_REF] Franceschini | Comparison of three scattering models for blood characterization[END_REF] as follows

σ ag,sim (k) = N σ b (k)S sim (k), ( 18
)
where σ b is the differential backscattering cross section of a single cell given by the fluid-filled sphere expression [START_REF] Morse | Theoretical Acoustics[END_REF][START_REF] Insana | Acoustic scattering theory applied to soft biological tissues[END_REF]

σ b (k) = k 4 V 2 c γ 2 z 4π 2 F (k; a), (19)
and S sim is the simulated structure factor representing the spatial positioning of the cells inside the aggregate and defined as

S sim (k) = 1 N N j=1 e -2ikn0•rj 2 (20)
where the cell positions r j are given using the procedure described in section III.A. The simulated structure factor S sim was determined from the 3D Fourier transform of the spatial distribution of cells. 17 A mean σ ag,sim curve was computed by averaging over 1000 different realizations for averaging purposes.

C. Computation of a simulated BSC from several aggregates

Random distributions of non-overlapping aggregates were computed within the simulated cubical volume 480 3 µm 3 by specifying the volume fraction φ of cells of radius a, the aggregate radius r ag and the aggregate compactness φ i . The spatial distribution of cells inside each aggregate was performed using the procedure described in section III.A. The simulated cubical volume was periodized; i.e., interactions between aggregates were determined under periodic boundary conditions, in order to remove the edge effects. The BSC from aggregated cells was then computed using the SFM as follows

BSC sim (k) = φ V c σ b (k)S sim (k), (21)
where S sim is the simulated structure factor representing the spatial positioning of the cells inside the simulated tissue. For each distribution of cells, the 3-D Fourier Transform of the cell spatial positioning was computed to obtain the simulated structure factor. [START_REF] Franceschini | Modeling of ultrasound backscattering by aggregating red blood cells[END_REF] . The mean BSC sim curve was obtained by averaging over 100 realizations for averaging purposes. Only 100 realizations were used here (against 1000 realizations in the case of the σ ag,sim computation) because one simulated volume involve several thousands of cells (against a few dozen of cells per aggregate for the σ ag,sim computation).

IV. RESULTS

A. Comparison of the simulated and theoretical σag

The frequency-dependent differential backscattering cross section of one cell aggregate σ ag,sim computed with the SFM is given in Fig. 2 for a single aggregate of radius r ag /a = 6 and of compactness φ i = 30%. Also shown in Fig. 2(a) (and Fig. 2(b), respectively) are the theoretical differential backscattering cross section considering only the coherent component (and considering both coherent and incoherent components, respectively). In Fig. 2(a), the dashed line (or the solid line, respectively) represents the theoretical σ ag,coh given in Eq. (11) using F (r ag ) (or F ( 5/3R g ), respectively). The first peaks of σ ag,sim and σ ag,coh using F ( 5/3R g ) occur at the same frequency (approximately 22.1 MHz). On the other hand, the first peak of σ ag,coh using F (r ag ) does not match perfectly the simulation and occurs at a lower frequency around 20.5 MHz.

In Fig.

2(b), the solid line represents σ ag = σ ag,coh + σ ag,inc using Eqs. (1), ( 11) and ( 15) and the dashed line represents σ ag,G = σ ag,coh + σ ag,inc,G using Eqs. (1), (11) and (16). The consideration of the incoherent component allows to better match the σ ag,sim behavior at higher frequencies, whatever the considered expression for the incoherent component σ ag,inc or σ ag,inc,G . Indeed, the σ ag,coh curve shows very deep dips (Fig. 2(a)), whereas both σ ag and σ ag,G curves match much better the dip behavior of the σ ag,sim curve.

Figure 3 represents the simulated σ ag,sim (symbols) for a radius r ag /a = 6 and aggregate compactnesses φ i = 10, 20, 30 and 40%. One can note that the σ ag,sim curves have more pronounced frequency dips when the aggregate compactness increases. For the lower aggregate compactness, only the first peak is well pronounced and the first dip is smooth, whereas the first two peaks and the first dip are clearly enhanced for the higher aggregate compactnesses. This behavior suggests that for compact aggregates, the incident wave tends to be scattered as if the aggregate was a well defined sphere, whereas for the lower aggregate compactness, the boundaries of the aggregates are not well defined. The occurence of the first frequency peak also shifts toward higher frequencies from 22.6 MHz to 25.8 MHz when the aggregate compactness decreases from 40% to 10%. Also represented in Fig. 3 are the theoretical σ ag (solid line) and σ ag,G (dashed line). Both models are very accurate at φ i = 30%. The σ ag,G curve matches better the simulation results for the highest aggregate compactness φ i = 40%, whereas the σ ag curve matches better the two lower compactnesses φ i = 10% and 20%.

Figure 3 also allows to compare the general σ ag behavior over all frequencies, however, in order to evaluate a workable frequency bandwidth for each model, the product kr ag for which the relative error was less than 20% was determined, as done previously in our twodimensional study (see Fig. 7 of Franceschini et al. [START_REF] Franceschini | Forward problem study of an effective medium model for ultrasound blood characterization[END_REF] ). First, we assess the improvement of the kr ag limit for the σ ag,coh expression using F 5/3R g (see the solid curves in Fig. (4)). In comparison with the former modeling using F (r ag ), the use of F 5/3R g increases the kr ag limit from 1.07 to 1.76 on average for φ i = 10%, and from 1.47 to 1.91 on average for φ i = 40% (the limit being computed as the mean kr ag for r ag /a values varying from 5 to 9). Secondly, we focus on the improvement of the kr ag limit considering the incoherent component expressions (see the dashed curves in Fig. (4)). For the lowest aggregate compactness of 10%, consideration of the incoherent component based on the structure factor σ ag,inc allowed to increase the kr ag limit from 1.76 to 2.12, whereas the incoherent component based on the Gaussian direct correlation function σ ag,inc,G did not allow to improve the kr ag limit. For the highest aggregate compactness of 40%, there was no significant difference in the kr ag limit whether the incoherent component was considered or not: it means that the kr ag limit increase is only due to the modification of F 5/3R g in the coherent component. To conclude, the overall best results are obtained by taking into account (1) the modification of F 5/3R g in the coherent component, and

(2) the incoherent component based on the structure factor σ ag,inc . Thus, in the sequel, for the results presented later, we will use the differential backscattering cross section σ ag computed with Eqs. (1), ( 11) and (15).

B. Comparison of the simulated and theoretical BSC

Figure 5 (a) shows simulated and theoretical frequency dependent BSCs for three aggregate sizes r ag /a = 5, 7 and 9, a fixed aggregate compactness φ i = 40% and a fixed cell volume fraction φ = 16%. The symbols represent the BSC sim computed with the SFM. The solid and dashed lines depicts the BSC EMTSFM using the proposed model of Eqs. (1), ( 11), (15), and (17). The theoretical BSC EMTSFM matches well the simulated BSC sim , especially for the amplitude and the frequency occurence of the first two peaks and dips, whatever the aggregate size considered in this work. In Fig. 5(b) are shown simulated and theoretical BSCs for three aggregate compactnesses φ i = 10, 25 and 40%, at a lower φ = 4%, and a fixed aggregate radius of r ag /a = 7. Note that for the study of the aggregate compactness variation, we limit the total cell volume fraction to φ = 4% and the lower aggregate compactness to φ i = 10% (corresponding to a value φ ag = 40%), because the maximum volume fraction of aggregates achievable in computer simulation is approximately 40% using the method described section III.C.

The influence of the structure factor of the effective medium S ag on the BSC curves was also studied. In that aim, we observe the spectral slope (i.e. the linear slope of the BSC as a function of frequency in a loglog scale) of the BSC sim curves before the first peak occurrence. The spectral slope increases above 4 for the BSC sim curves with an aggregate concentration φ ag = φ/φ i = 40% presented in Fig. 5(a), and for the BSC sim curve with the clustering condition (φ = 4%, φ i = 10%, φ ag = φ/φ i = 40%) presented in Fig. 5(b). On the other hand the BSC sim curve shows a fourth power frequencydependence (i.e. a spectral slope of 4) with the clustering condition (φ = 4%, φ i = 40%, φ ag = φ/φ i = 10%) presented in Fig. 5(b). So, the more the aggregate concentration φ ag increases, the greater the influence of S ag becomes observable in the BSC frequency dependence.

V. DISCUSSION

A. On the use of simple spatial distribution of cells

The method used in this study to obtain the cell spatial distribution was not based on a generic physical model of interactions between cells. Rather, we used a simple and fast method to generate hundreds of simulated media with a large scale ratio between the cell size (2.75 µm) and the whole medium size (480 3 µm 3 ). The purpose was to build a controlled medium containing non-overlapping spherical aggregates, all aggregates having the same radius and compactness with a unique cell spatial distribution defined by a structure factor. The main advantage of this method was the possibility to have various aggregate compactnesses with the same size of aggregates, which allowed to demonstrate the role of the incoherent scattering component on the BSC frequency dependence at high frequencies. Indeed, the σ ag,sim curves had well pronounced peaks and dips for the higher aggregate compactnesses, whereas the first peak was well enhanced and the first dip was smooth for the lower aggregate compactness (see Fig. 3). Theoretical predictions of the σ ag based on the structure factor for the incoherent component proposed in Eq. ( 15) agreed well with the simulations and validated the proposed theoretical modeling.

In the 3-D computer simulations of this study, we had to limit the total volume fraction of cells to a maximum of 16%. Indeed, the procedure we chose to distribute the cells within aggregates allowed reaching a maximum aggregate compactness φ i,max of 40%. The upper limit of 40% is easily understandable because placing nondeformable identical spheres randomly enough inside a larger sphere without overlapping is not straightforward. In comparison, a standard method using a random sequential absorption [START_REF] Hinrichsen | Geometry of random sequential adsorption[END_REF] would give a compactness close to 30%. The maximum aggregate volume fraction φ ag,max was thus also fixed to 40%. As a consequence, the maximum value of the total volume fraction φ max was limited to φ max = φ i,max φ ag,max = (40%) 2 = 16%.

B. The benefit of considering the incoherent component

A new theoretical modeling of the differential backscattering cross section from a single cell aggregate was developed in this work and consisted in taking into account a coherent component and an incoherent component in the expression of σ ag . This modeling was compared to numerical simulations based on the SFM. Concerning the coherent component, a slight modification of the σ ag,coh expression using F 5/3R g was proposed and allowed the EMT to better match the frequency occurrences of peaks and dips observed in simulations. Concerning the incoherent component, two expressions of σ ag,inc based on the structure factor and on a Gaussian direct correlation function proposed by Morse and Ingard 12 were compared. The numerical study demonstrated the superiority of the formulation using the structure factor over that with the Gaussian direct correlation function, regarding the kr ag limit study of Fig. 4. The consideration of the incoherent component σ ag,inc based on the structure factor allowed to approximate the simulations satisfactorily for an average kr ag limit around 2, against an average kr ag limit comprised between 1.07 and 1.47 with the former model considering only the coherent component.

It is interesting to observe that the simulated σ ag,sim curves showed a fourth power frequency dependence (Rayleigh scattering) before the first peak occurrence (see Figures 2 and3). When the aggregate compactness varied (Fig. 3), the σ ag,sim frequency dependence differed mainly at high frequencies after the first peak. To better understand this high frequency behavior, we analyzed the theoretical respective influence of each scattering component (coherent and incoherent) on the σ ag behavior, while varying the aggregate compactness. According to Eqs. ( 11) and ( 15), the σ ag frequency dependence is mainly determined by the spherical form factor from an effective sphere F 5/3R g that intervenes in the coherent component, and the structure factor S(a, φ i ) that intervenes in the incoherent component. As shown in Fig. 6(a), an increase of the compactness induces a decrease of the structure factor amplitude S(k; a, φ i ), and a slight shift toward lower frequencies of the form factor F (k; r ag ) peaks and dips. In Fig. 6(b) (or (c), respectively) are plotted the theoretical backscattering cross sections (solid lines), the corresponding coherent component σ ag,coh (dashed lines) and incoherent component σ ag,inc (dotted lines) as defined in Eq. (1), for the same aggregating condition with r ag /a = 6 and φ i = 10% (or φ i = 40%, respectively). The relative influence of the incoherent component is more important for lower aggregate compactnesses, giving rise to smoother peaks and dips. Also, the frequency shift of the peaks and dips observed in the simulations can be explained theoretically by a dependance of the form factor F 5/3R g on the aggregate compactness. Indeed, the gyration radius value depends on the aggregate compactness as mentioned in section III.A. In the framework of an inverse problem approach, it would be useful to approximate the gyration radius as a function of the radius r ag of the aggregate envelope and the aggregate compactness φ i . Based on the 3D building of isotropic aggregates using the method described in section III.A, the relationship between R g , r ag and φ i was empirically determined as follows: R g ≈ 3/5 0.96 × r ag + 7.06 × 10 -6 φ i -2.71 × 10 -6 .

(22) This approximation, obtained from a linear regression analysis, shows less than 5% relative error for all tested radii r ag /a (varying from 4 to 9) and compactnesses (varying from 5% to 40%). The estimated R-squared value and p-value were 0.996 and 5.39 × 10 -101 respec-tively.

VI. CONCLUSION

The EMT was further developed to decompose the differential backscattering cross section of a single cell aggregate into coherent and incoherent components. The coherent component, corresponding to the average wave emerging from the effective scatterer, was approximated by the scattering of an effective fluid sphere, whose radius depends on the gyration radius of the cell aggregate. An important contribution of this new EMT is the taking into account of the incoherent component based on the structure factor, which allowed to approximate the computer simulations satisfactorily for a product kr ag up to 2.

The EMTSFM assumes that all cells are aggregated and that aggregates are identical and isotropic. Therefore, the BSC behavior obtained in all simulations showed a pronounced first peak. In experimental conditions when insonifying aggregated red blood cells 9 or clusters of tumor cells 5 , the BSC behavior was smoother and the peaks were less pronounced. The reason behind this might be that real tissues contain several sizes of aggregates, and since the location of BSC peaks are different for different aggregate populations, a relatively smoother BSC curve can be obtained. Therefore, future improvements should consider incorporating polydispersity in aggregate size and compactness to provide an optimal EMTSFM for the inversion of experimental data.

We then compute based on Eq. ( 23)

Vag γ κ (r 0 )-γ ρ (r 0 ) e 2ik n0•r0 d 3 r 0 = N j=1 (γ κ -γ ρ ) Vag χ(r 0 -r j )e 2ik n0•r0 d 3 r 0 .
Next, considering the change of variable r 0 = r 0 -r j , one obtains Vag χ(r 0 -r j )e 2ik n0•r0 d 3 r 0 = Va e 2ik n0•(r 0 +rj ) d 3 r 0 .

Lastly, since e 2ik n0•(r 0 +rj ) = e 2ik n0•r 0 e 2ik n0•rj , one concludes that

Vag γ κ (r 0 ) -γ ρ (r 0 ) e 2ik n0•r0 d 3 r 0 = (γ κ -γ ρ ) Va e 2ik n0•r 0 d 3 r 0 N j=1 e 2ik n0•rj = (γ κ -γ ρ )V c V -1 c Va e 2ik n0•r 0 d 3 r 0 N j=1 e 2ik n0•rj . Appendix B
This appendix gives the computation steps to obtain the low frequency approximation of F coh (k) = F 0 (k; a)1/N N j=1 e 2ik n0•rj . Firstly, by assuming centrosymetric aggregates, one has

1 N N j=1 e 2ik n0•rj = 1 N N j=1 cos(2k n 0 • r j ) = 1 N N j=1 sin(2k |r j |) 2k |r j | ≈ 1 - 2 3 k 2 1 N N j=1 |r j | 2 ,
where the last approximation corresponds to the second order Taylor expansion of the expression. An equivalent approach can be found in Guinier & Fournet 19 pp. 7 and 8. A second order Taylor expansion of F coh (k) thus gives:

F coh (k) ≈ 1 - 2 3 k 2   3 5 a 2 + 1 N N j=1 |r j | 2   = 1 - 2 3 k 2 R 2 g ,
by using the definition of R g given in Eq. ( 9), and the Taylor expansion of F 0 (k; a)(≈ 1 -2 5 k 2 a 2 ). Lastly, since for the equivalent full sphere of radius r , F 0 (k; r ) ≈ 1 -2 5 k 2 r 2 , we need r = (5/3)R g in order to match the low frequency approximations of the form factors F coh (k) and F 0 (k; r ).

Appendix C

This appendix gives the computation steps to obtain the incoherent component of the differential backscattering cross-section σ ag,inc,G based on a Gaussian direct correlation function as proposed by Morse and Ingard. [START_REF] Morse | Theoretical Acoustics[END_REF] Firstly, one has γ κ,ag (r 0 ) -γ ρ,ag (r 0 ) = χ(r 0 )(γ κ -γ ρ ), where χ is the characteristic function of the cells: χ(r 0 ) = 1 inside cells and χ(r 0 ) = 0, otherwise. Moreover, we use the following approximation γ κ,ag (r 0 ) -γ ρ,ag (r 0 ) ≈ φ i (γ κ -γ ρ ).

From there, one obtains the expression

Φ ag,inc (k) = Φ ag (k ag ) -Φ ag (k ag ) ≈ k 2 ag 4π Vag (γ κ -γ ρ )(χ(r 0 ) -φ i )e 2ikag n0•r0 d 3 r 0 .
Note that the function (γ κ -γ ρ )(χ(r 0 ) -φ i ) corresponds to the function δ modified from Ref. 12 (Eq. (8.2.24)). Therefore, the incoherent component of the differential backscattering cross-section can be computed using the direct correlation function Υ of the compressibility and density fluctuations as follows:

σ ag,inc (k) = |Φ ag,inc (k)| 2 ≈ k 4 ag (γ κ -γ ρ ) 2 16π 2 V ag Vag Υ(∆r)e 2ikag n0•∆r d 3 ∆r, (24)
where

Υ(∆r) = 1 V ag Vag (χ(r 0 ) -φ i )(χ(r 0 + ∆r) -φ i ) d 3 r 0 .
The function Υ depends only on ∆r = |∆r| and must have the following properties: the mean value of Υ(∆r) is zero and lim |∆r|→∞ Υ(∆r) = 0. Under the postulate of a Gaussian model, the expression of the direct correlation function Υ is given by Morse & Ingard 12 (Eq. (8.2.26)):

Υ(∆r) ≈ Υ(0) 1 - ∆r 2 3d 2 e -∆r 2 /(2d 2 ) ,
where Υ(0) = φ i (1-φ i ), and the correlation distance d is related to the cell radius a as d = 2 1/6 3 1/3 π 1/6 a ≈ 0.643092 × a as defined by Insana and Brown 13 (Eq. (81) p.107). Performing the integral in Eq. ( 24), one then obtains the Gaussian incoherent component expressed in Eq. ( 16). 11) and ( 15), (or the theoretical σ ag,G using Eqs. (1), ( 11) and ( 16), respectively). 11) and ( 15), (or the theoretical σag,G using Eqs. (1), ( 11) and ( 16), respectively).

FIGURE CAPTIONS Figure 1 .

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 .

 2 Figure 2. (Color online) Differential backscattering cross section of one cell aggregate of radius r ag /a = 6 and of compactness φ i = 0.3. The symbols represent the σ ag,sim computation. (a) The dashed line (or the solid line, respectively) represents the theoretical coherent component σ ag,coh given in Eq. (11) using F (r ag ) (or F 5/3R g , respectively). (b) The solid line (or the dashed line, respectively) represents the theoretical σ ag(or σ ag,G , respectively), both using F 5/3R g in the coherent component.

Figure 3 .

 3 Figure 3. (Color online) Differential backscattering cross sections of one cell aggregate for different aggregate compactnesses. The symbols represents the σ ag,sim computation. The solid line (dashed line, respectively) represents the theoretical σ ag using Eqs. (1), (11) and (15), (or the theoretical σ ag,G using Eqs. (1), (11) and (16), respectively).

Figure 4 .

 4 Figure 4. (Color online) Averaged kr ag limits, as a function of the aggregate size ranging for r ag /a =5 to 9, for three compactnesses φ i = 10% (a), φ i = 25% (b) and φ i = 40% (c). The solid lines correspond to the limits of σ ag,coh computed using F 5/3R g (crosses), and using F (r ag ) (squares). The dotted lines correspond to the limits of σ ag (diamonds) and σ ag,G (circles), both computed using the modification of the radius in F 5/3R g . Note that the curves marked with crosses and circles are superimposed in Figs. (a) and (b).

Figure 5 .

 5 Figure 5. (Color online) Frequency dependent backscatter coefficients for various aggregate sizes and compactnesses.Symbols represent simulation results and lines represent theoretical BSCs computed using Eq. (17). (a) The volume fraction occupied by the cells in the medium is φ = 16% and the aggregate compactness is fixed to φ i = 40%. Radii of aggregates vary from r ag /a = 5 to r ag /a = 9. (b) The volume fraction occupied by the cells in the medium is φ = 4% and the aggregate size is fixed to r ag /a = 7. Aggregate

Figure 6 .FIG

 6 Figure 6. (Color online) (a) Theoretical form factor F k; 5/3R g and structure factor S(k; a, φ i ) computed as a function of the frequency for two compactnesses φ = 10% and φ i = 40% and aggregate size r ag /a = 6. (b) and (c) Theoretical differential backscattering cross sections of one cell aggregate σ ag , its coherent component σ ag,coh and incoherent component σ ag,inc , as a function of frequency.

Frequency

  FIG. 2. (Color online) Differential backscattering cross section of one cell aggregate of radius rag/a = 6 and of compactness φi = 0.3. The symbols represent the σag,sim computation. (a) The dashed line (or the solid line, respectively) represents the theoretical coherent component σ ag,coh given in Eq. (11) using F (rag) (or F 5/3Rg , respectively). (b) The solid line (or the dashed line, respectively) represents the theoretical σag (or σag,G, respectively), both using F 5/3Rg in the coherent component.

Frequency

FIG. 4 . 16 Frequency

 416 FIG. 4. (Color online) Averaged krag limits, as a function of the aggregate size ranging for rag/a =5 to 9, for three compactnesses φi = 10% (a), φi = 25% (b) and φi = 40% (c). The solid lines correspond to the limits of σ ag,coh computed using F 5/3Rg (crosses), and using F (rag) (squares). The dotted lines correspond to the limits of σag (diamonds) and σag,G (circles), both computed using the modification of the radius in F 5/3Rg . Note that the curves marked with crosses and circles are superimposed in Figs. (a) and (b).

  ,sim (r ag /a=6, i =10%) ag,sim (r ag /a=6, i =30%) ag,sim (r ag /a=6, i =20%) ag,sim (r ag /a=6, i =40%) FIG. 3. (Color Differential backscattering cross sections of one cell aggregate for different aggregate compactnesses. The symbols represents the σag,sim computation. The solid line (dashed line, respectively) represents the theoretical σag using Eqs. (1), (

			10	-14
	Differential backscattering cross section	of one cell aggregate (dB.mm -1 .sr -1 )	10	-16	ag ag,G	ag ag,G

ag

ACKNOWLEDGMENT This research was supported by the Canadian Institutes of Health Research (grant #MOP-84358), the CNRS, the CSIR New Delhi, as well as the Labex MEC (ANR-10-LABX-0092) and the A * MIDEX project (ANR-11-IDEX-0001-02) funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR). F. Destrempes was an invited researcher at the Laboratoire de Mécanique et d'Acoustique LMA, CNRS, UPR 7051, Aix-Marseille University, Centrale Marseille, F-13402 Marseille Cedex 20, France, in July 2013 to work on this project.

Appendix A

The model underlying Eq. ( 2) is equivalent to the following expression for the tissue function

where χ(r 0 ) = 1 inside a sphere V a of radius a centered at the origin, and 0 elsewhere.