Esteban P Busso

L Wright

H E Evans

L N Mccartney

S R J Saunders

S Osgerby

J Nunn

A physics-based life prediction methodology

Keywords: TBC, Finite Element Modelling, Life Prediction Model, MCrAlY

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The ability to predict the lifetime of thermal barrier coating (TBC) systems is a high priority for gas turbine users. Typical TBC systems are made of a top layer of yttria stabilized zirconia (YSZ -the insulating layer) and a metallic low pressure plasma sprayed (LPPS) bond coat. During exposure at high temperatures, the metallic bond coat forms a thermally grown oxide (TGO) that consists predominantly of alumina. The formation of the TGO plays a crucial role in the failure of the TBC. The usual mechanism of TBC failure is associated with spallation at or close to the TGO interfaces within either the YSZ or the bond coat. A typical crack pattern in an EB-PVD MCrAlY bond coat after 1000 o C exposure for 700 h is shown in Fig. 1(a). In the micrograph of Fig. 1(b), possible initiation sites are also shown of interfacial microcracks between the TGO and the bond coat of an EB-PVD TBC system exposed to 1000 o C for 600 h.

Approaches to predict the mechanical integrity of TBCs have generally concentrated on the development of life prediction models based on empirical/phenomenological fatigue life relationships which link the life of the TBC to the summation of damage due to mechanical straining and oxidation [START_REF] Padture | Thermal barrier coatings for gasturbine engine applications[END_REF]. Delamination models, which have been considered in some theoretical depth, deal with edge cracking and buckling [START_REF] He | Simulation of stresses and delamination in a plasma-sprayed thermal barrier system upon thermal cycling[END_REF]. However, the latter are widely perceived only as a possible but not critical part of the failure process since an extensive delamination zone at the base of the ceramic top coat is necessary for crack propagation and buckling to occur. The sub-critical growth of such a zone is generally the slow, life-determining stage in TBC systems. However, the conditions under which this could occur have received little attention.

In this work, a novel mechanistic approach is proposed to predict the life of TBC systems. It relies on a statistical treatment of the TBC's morphological characteristics and on a continuum mechanics framework to quantify the local TBC stresses responsible for sub-critical crack nucleation and growth.

Two key aspects characterise this life prediction methodology. Firstly, it incorporates the complex interaction between interfacial and microstructural features, local oxidation mechanisms and time-dependent processes (i.e. oxidation kinetics, diffusion, creep and sintering of the ceramic material). Secondly, its calibration is based on TBC damage inferred from non-destructive measurements and on the numerically predicted local TBC stresses responsible for the initiation of such damage. The nondestructive measurements are based on a combination of fluorescence measurements using piezo-spectroscopy and high resolution thermal images. Thus, the methodology requires TBC damage information generated by non-destructive measurements, a suitably defined failure criterion based on the dominant failure mechanisms, morphological characteristics of the TGO obtained from microstructural observations, and a numerically generated database containing the local predicted TBC stresses.

The emphasis of the paper will be on the different stages involved in the proposed life prediction methodology and on how to integrate the experimental and numerical data. Moreover, the mechanistic study of local TBC stresses will be summarised and, when necessary, reference to related published work will be given where further details about the approach can be found. The full validation of the methodology and an investigation into the effects of significant thermal cycling are to be reported in a future publication. The paper is structured as follows. In Section 2, an overview of the TBC damage measurements and overall mechanistic approach is first given. Representative numerical results of the TBC stress analyses are then presented in Section 3. Section 4 describes the different stages involved in the TBC life prediction methodology and highlights its predictive capabilities and potential for further development. Finally, Section 5 presents some general conclusions. 2 Mechanistic Approach To Model TBC Stresses

Mechanistic Considerations

Experimental work (e.g. [START_REF] Padture | Thermal barrier coatings for gasturbine engine applications[END_REF][START_REF] He | Simulation of stresses and delamination in a plasma-sprayed thermal barrier system upon thermal cycling[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part I: Model Formulation[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF][START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF]) has shown that the failure modes of EB-PVD TBCs are very dependent on the properties of the bond coat, the YSZ, and the TGO growth rate. In addition, the stress conditions in and around the TGO layer are significantly affected by the aspect ratio of the wavy TGO morphology. One of the most commonly observed TBC failure mechanisms is the spallation of the top YSZ coating driven by the high interfacial tractions that develop at the TGO interfaces upon cooling (e.g. see [START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part I: Model Formulation[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF][START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF]). Experimental evidence of interfacial damage is shown in Fig. 1(b), where the arrows indicate possible initiation sites of interfacial cracks between the TGO and the bond coat. The driving force for the nucleation of such types of defects is intrinsically linked to the local stress conditions, as illustrated in Fig. 2. Here, a schematic diagram of a TBC system shows the different relevant stresses and tractions, which need to be considered when studying the failure of the rough undulating interfaces often found in TBC systems.

The out-of-plane stress component, parallel to the X 2 axis in Fig. 2, is responsible for the formation of microcracks lying in the X 1 -X 3 plane, which are generally seen within the YSZ near the TGO interface. The interfacial TGO tractions, normal to either the top coat-TGO or the bond coat-TGO interfaces, drive the nucleation of interfacial microcracks and propagation of interfacial defects. Finally, the maximum principal stress is generally identified as the main 'driving force' for the nucleation of cohesive microcracks within the TGO.

Figure 2

Diagram illustrating the structure of the TBC system and some of the stresses that are used in determining the failure of the system.

TBC Damage Measurements

The non-destructive measurements of TBC damage were based on a combination of high resolution thermal images and fluorescence measurements using Raman spectroscopy. The latter relies on the fluorescence of trace amounts of Cr ions in the TGO and on the fact that spectral positions depend on the magnitude of the stresses in the sampled TGO volume. Thus, stress changes associated with microcracking / debonding of the alumina layer can be readily detected. The stress state can be estimated from the peak shift of the spectra with respect to the unstressed sapphire reference. In the most general case, the frequency change is related to the stress components via the relation

∆ν = Π σ + Π σ + Π σ (1)
where ij Π are piezo-spectroscopic coefficients for the system and σ ij are the stress components. By assuming a random distribution of the crystallographic texture in the TGO, isotropic properties apply so that 11 22 33

Π = Π = Π = Π leading to ( ) 11 22 33 ∆ν = Π σ + σ + σ , ( 2
)
implying that what is measured is the hydrostatic stress, an invariant of the stress field. Alternatively, if an equi-biaxial plane stress state is assumed, then the frequency shift can be related to just the biaxial stress component (e.g., see Selculk and Atkinson [START_REF] Selcuk | Analysis of the Cr3+ luminescence spectra from thermally grown oxide in TBCs[END_REF]). However, such a plane stress assumption is particularly suited to rather flat oxide layers which, as will be discussed later on in the text, are not generally seen in TBC systems. Such an approximation will not be correct for rough TGOs, as the waviness of the TGO gives rise to non-zero stress components in the direction normal to the TBC surface. to define the boundary between the 'damaged' and 'undamaged' parts of the TBC. The required damaged region is defined as a region initially under high compressive stress that moves to a lower compressive or tensile stress state such that the mean stress is above the defined threshold value. The results obtained by piezospectroscopy can be corroborated metallographically by examining sections taken through the TBC sample in order to define a relationship between observed cracking/damage and the measured mean stress in the damaged region. Alternatively, thermography may also provide the required information (see for example [START_REF] Nunn | Application of thermography in the evaluation of early signs of failure of thermal barrier coatings[END_REF]), as shown by the thermography measurements of Fig. 3(d) conducted on the same specimens as those used for the piezo-spectroscopy measurements. The thermographs shown in Fig. 3(d) were obtained at room temperature when radiating the surface of the specimen with low level infra-red. The apparent temperature difference between the black and grey pixels shown on the specimen in this figure is of the order of 1°C, and the background is approximately 3°C cooler than the specimen. A comparison between Figs. 3(c) and (d) reveals good agreement between the stress-relieved areas with the coolest ones from the thermal images, giving credence to the understanding that these are indeed areas where damage has occurred.

Mechanistic Approach to Determine the Local TBC Stresses

The mechanistic approach followed here is based on the work of Busso et al. [START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part I: Model Formulation[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF][START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF] and accounts for the effects of YSZ sintering, bond coat creep and TGO morphology on the local stresses responsible for interfacial cracking. The oxidation-induced stresses that control the TBC degradation and eventual failure, such as the components illustrated in Fig. 2, are calculated using a constitutive framework coupled with the diffusion of oxidant species through a multi-phase solid. The calculations are performed using the ABAQUS finite element (FE) analysis code. The overall thickness of the oxide layer is taken to be a phenomenological function of temperature and time. The proposed internal state variable constitutive model accounts for the phase transformation associated with the oxidising phase, and incorporates the effect that the local volumetric expansion of the newly formed oxide has on the generation of inelastic volumetric strains and stresses.

The strains associated with the oxidation-induced phase transformation of the bond coat into the internal primary oxide were calculated in several stages:

1.

The total thickness h of the oxide layer at the current time was first computed from a phenomenological oxide growth kinetics relation in terms of temperature and accumulated time, t, at the maximum temperature in the cycle, θ max [START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF],

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ θ - θ = max 0 1 1 exp R O q R Q t A h (3)
where A 0 is a proportionality constant, q a growth exponent, Q 0 an apparent activation energy, θ R a reference temperature, and R the universal gas constant. For the TBC system of interest in this work, the calibrated oxidation parameters are A 0 = 10 -6 m/s q , θ R = 2424 K, Q 0 = 766,900 J/mol and q = 0.332. Thus the units of h and max θ are m and K, respectively.

2.

The total oxide thickness was then subdivided into an internally grown part, of thickness, h in , and an externally grown component of thickness,

h out = h -h in ,
with their ratio given by h in /h out = 1 /(PBR -1). Here PBR = 1.28 is the Pilling-Bedworth ratio corresponding to the primary oxidation reaction.

3.

The volumetric strains associated with the oxidation reaction, whose mean value is given by e T = log e (PBR / 3), were applied to the row of elements in the FE model undergoing internal oxidation. Note that this expansion is irreversible and generally anisotropic, with the anisotropy of the transformation determined experimentally. The in-plane normal strains are given by ε X1 = ε X3 = 3 e T r/(1 + 2r), and the out-of-plane normal strain is ε X2 = 3 e T /(1 + 2r), where r = 0.005/0.435 = 0.0115 (see [START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF] for further details).

The constitutive model for the oxidising bond coat material was implemented into the FE method as a user-defined material subroutine and used to update the stresses at each material point within the unoxidised and oxidising bond coat in the finite element models of the TBC, as described in the next Section. For further details about the approach, refer to [START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part I: Model Formulation[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF][START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF].

The temperature dependence of the elastic material properties associated with the bond coat and the TGO are shown in Table A1 of Appendix A. The table also includes the yield stress of the bond coat, used in an incremental model of perfect plasticity, which is separate from the power law creep model used for the bond coat and the TGO (see Appendix A). The YSZ top coat is modelled as a transversely isotropic elastic solid, with properties dependent on the degree of sintering (parameterised by sintering time) and temperature. The YSZ elastic moduli, taken from [START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF], are given in Appendix A.

Identification of an Appropriate Criterion for TBC failure

It is expected that, for a given oxidation temperature, a TBC system of the type studied here would fail by the time the oxide reaches a critical thickness. The definition of the critical TGO thickness can be made from the type of data presented in Fig. 4, which shows the effect of oxidising temperature on the time to TBC spallation for MCrAlY type coatings and for PtAl based TBCs. The data points in the figure taken from the open literature [START_REF] Evans | Delamination Processes in Thermal Barrier Coating Systems[END_REF] are the closed symbols. The open symbols represent the data obtained as part of this work. In each case, circles and diamonds are MCrAlY coatings and triangles are PtAl. As in the original data set, there does not appear to be a systematic difference between the two sets although the current data set on the aluminides is limited. The striking feature of these data is the large spread (a factor of ~10) in the times to spallation at any given temperature, even for the present samples which were nominally identical and tested in the same facility. This large scatter may be due to an uncontrolled variation in the bond coat surface roughness [START_REF] Evans | Delamination Processes in Thermal Barrier Coating Systems[END_REF] since this is expected to produce a corresponding variation in the stresses developed within the TBC system. An objective of the present FE analysis is to evaluate this possibility. The results from the current data, in this regard, are typical of those found in the literature. They also show similar failure times to those reported elsewhere so that the two sets of data appear mutually consistent. The continuous line in Fig. 5 represents the time at a given temperature to reach an alumina TGO thickness of 3 µm, calculated using Eq. 3. It can be seen that this prediction at the 3µm TGO thickness level forms an approximate lower bound to the data set with the possible exception of tests at 1000 o C where they provide an underestimate.

Identification of Representative Periodic Unit Cells and FE Model

The morphological characteristics of the TGO interface regions of the EB-PVD MCrAlY bond coating, see Figs. 1 and2, can be represented by the series of undulations shown in Fig. 5. Here, the random nature of the interface morphology can be idealised by a periodic representative profile such as the unit cell shown in Fig. 5. Two interface roughness parameters characterise the unit cell, namely (i) the half-period a and (ii) the amplitude b. Note that whilst b/a constitutes the most important roughness parameter, b is also required to account for the fact that its value relative to that of the TGO thickness will affect the local stresses. Typical TGO-YSZ profiles were identified from SEM micrographs of the interfacial regions for the EB-PVD TBC /LPPS bond coat system of interest and analysed to obtain the statistical distributions of the geometric parameters b and b/a. In order to obtain the representative spread of 'b', 'a' and 'b/a' ratios of the undulations of the bond coat surface, a specimen was cross sectioned and polished to reveal the interface. Over 100 measurements were made on a total of 28 scaled micrographs. Typical mean and standard deviation roughness values for the TBC system were found to be 2.39 µm and 1.62 µm for b and 0.32 and 0.12 for b/a. Three sets of interfacial roughness parameters (b, b/a) were chosen for the stress analysis calculations from the measured roughness distributions, so as to be as close as possible to the mean value ± one standard deviation for each interfacial parameter. This yielded, The TBC is assumed to be subjected to a generic temperature cycle consisting of three stages: a heating transient from a minimum θ min to a maximum, θ max , at a constant rate, followed by a dwell period of length, t hold , and a cooling transient to θ min at the same rate as in the heating transient. In the work conducted here, θ min is always assumed to be 25 ˚C, and θ max is set to typical TGO temperatures found during service in land-based gas turbines, namely either 1000 ºC or 1100 ºC. Based on the experimental evidence shown in Fig. 4 and the discussions in Section 2.4, the holding time t hold was chosen such that the oxide layer would have reached the critical thickness of 3 µm after one thermal cycle.

As previously discussed, the aim of the FE modelling work was to carry out a set of parametric calculations, the results of which are required as inputs for the prediction of coating life, as will be described in Section 4. The parametric variables in the numerical study were the maximum cycle temperature, θ max , and the interfacial parameters b/a and b (here assumed to be perfectly correlated). Finite element models of the periodic unit cells for the three representative TGO-YSZ profiles were then developed. The models of the TGO-YSZ interfacial regions were constructed from the single geometric unit of the form shown in Fig. 5. The surface of the TBC was assumed to be traction free, and symmetry and periodic boundary conditions were assumed on the planes X 1 = 0 and X 1 = a, respectively (see Fig. 5). In this way, the interface is effectively modelled as an infinite series of undulations of period 2a and amplitude b. The initial thicknesses of the bond coat and YSZ top coat were taken to be 50 and 125 µm, respectively. The temperature within the TBC system was assumed to be uniform.

The 2D model of the TBC system relied on generalised plane strain assumptions (i.e. strain along X 3 , see Fig. 5, is uniform but not necessarily zero). Furthermore, the TBC system was assumed to be attached to a bulk substrate whose deformation is determined solely by its thermal expansion, that is, it is elastically rigid so that its Young's and shear moduli are effectively assumed to be infinite. This implies that the in-plane displacements of the X 1 and X 3 boundaries of the representative unit cell model are controlled by the expansion and contraction of the substrate. Note that the shape of the curved parts of the periodic unit cell is chosen to be consistent with typical TGO interface profiles observed in SEM micrographs such as that shown in Fig. 1. Here, the interface morphology was constructed by linking straight line segments and circular arcs in order to yield a curve with a continuous first derivative.

Stresses Obtained from Parametric FE Studies

Some representative results from the finite element analyses are presented and discussed in this Section. Note that the TGO thickness is 3 µm in all the plots. Similar results were obtained for oxidation at 1100 o C, but contours at this higher temperature are not shown in this paper for brevity. It can be seen that the extent and magnitude of the most highly stressed areas increases with interface roughness, and that these regions are located in the YSZ top coat around the valleys. It is also worth noticing that both the TGO and the YSZ regions tend to be generally under some compression at 25 o C for the rougher TBCs. Thus these results show that, for smoother interfaces, peak tensile stresses develop within the TGO but, for rougher interfaces (such as those shown in Fig. 1), significant tensile stresses occur at the bond coat / TGO apex.

Figures 7(a) to (c) show the contour plots of the out-of-plane stress component σ 22 when the TBC was cooled down to 25 ºC after (a) 61 h, (b) 167 h and (b) 241 h oxidation at 1000 ºC. At this temperature, the corresponding oxide thicknesses are 1.91, 2.65 and 3.00 µm, respectively. The initial interfacial geometry used was b=2.47 µm, b/a = 0.323. It can be seen that the most highly stressed regions are in the YSZ around the valleys, and that the magnitudes of the out-of-plane stresses in these regions increase rapidly with oxide thickness. Results at 1100 ºC show similar trends. In general, a comparison between the results obtained at different temperatures shows clearly that the TBC stresses do not depend just on the oxide thickness but also on the kinetics of oxidation. Whilst the stress contours calculated at different temperatures share some features, there are significant differences. A summary of the results obtained from the parametric finite element analyses at both 1000 and 1100 ºC is shown in Appendix B. The predicted effects can be seen of the geometric parameters and peak cycle temperature on the stress maxima relevant to the different TBC failure modes at room temperature.

Methodology for TBC Life Prediction

The proposed TBC lifetime prediction methodology can be considered as a threestage approach which requires the following inputs:

A well characterised thermal cycle (defined by θ max , θ min ,and t hold, and a specified heating and cooling rate, see Section 2.5).

The accumulated time defined as t acc = t hold N where N is the number of thermal cycles. The dominant failure mechanism for the TBC coating system and loading conditions of interest (as discussed in Sections 2.1 and 2.4).

A statistical distribution of the TBC's interface roughness (defined by the ratio b/a and amplitude b) obtained from actual TBC micrographs. Experimental measurements of TBC damage accumulation with thermal loading (see Section 2.

2). A full description of the material properties of all constituents of the TBC (see Appendix A).

A numerically generated database of the relevant local TBC stresses obtained from a parametric study such as that presented in Section 3 (see also refs. [START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part I: Model Formulation[END_REF][START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF][START_REF] Busso | A Mechanistic Study of Microcracking in Transversely Isotropic Ceramic-Metal Systems[END_REF]). A PC-based software tool, using Excel spreadsheets, has been developed to carry out some of the calculations required at the different stages of the life prediction methodology, which are described next. Fig. 8 illustrates the methodology as a flow chart.

Note that all interpolations mentioned in this Section are linear interpolations, unless specified otherwise. Furthermore, here the times have been converted into equivalent thicknesses of TGO and the interpolation has then been carried out with respect to TGO thickness. This enables the interpolation to be carried out more accurately due to the fact that the stresses are likely to depend more linearly on TGO thickness than on accumulated time. This Section illustrates the methodology using real measurement data and the FE results given in Appendix B. It should be noted that whilst the measured and modelled TBC systems are similar, they are not identical. No experimental data were available for the modelled system. The differences between the systems make some aspects of the software tool less reliable: in particular, the example shown contains more extrapolation of FE results than is desirable due to the long time over which the damage measurements were taken. It is, however, emphasised that the purpose of this paper is to describe the methodology, and not to validate the approach by comparing predictions with experimental results. For validation of the methodology, the FEA runs would need to be repeated using model properties that correspond more closely with those of the measured system.

Stage 1: Quantification of TBC Damage

This stage deals with the way in which the information about TBC damage (see Section 2.2) is integrated into the life prediction methodology. Two critical assumptions are made in this stage. Firstly, the total damaged area within the TBC is assumed to be directly associated with areas exhibiting detectable stress relief, which constitutes a measure of damage as shown in Fig. 3. The measured fraction of damaged area can be quantified by a generic dimensionless scalar variable or "damage parameter", D, and expressed as a function of accumulated time at the oxidising temperature, t acc , for given values of θ max and θ min , so that:

) t , , ( D D acc max min θ θ = ,
where

N t t hold acc = , ( 5
)
and N is the total number of thermal cycles. Thus Eq. 5 relates the accumulated time to the proportion of areas sampled that exhibited damage, which can be derived from measured stress maps such as those shown in Fig. 3(c). The measured t acc vs. D curve for the case where θ max = 1000 o C is shown in Fig. 9. Note that the TBC failed after approximately 700 h oxidation, failure in this case being taken as damage of more than 22%.

Thermal cycle characterisation

Damage accumulation measurements

Interface roughness distributions

Step 1: Quantification of TBC damage Link the damage measurements at each time to a critical interface roughness for a given thermal cycle.

Curve of critical interface roughness versus accumulated time for a given thermal cycle.

Step 2: Interpolation of finite element results

Interpolate between the FE results curves to get a critical stress value for the damage mechanism at each accumulated time

Material properties

Results of parametric FE study

Thermal cycle characterisation

Failure mechanism

Curve of critical stress versus accumulated time for a given thermal cycle.

Step 3: Life prediction Using the curve generated in step 2, calculate an estimated failure time and hence an estimated remaining life for the coating.

Failure mechanism Prediction of remaining coating life

Figure 9 Measured fraction of damaged area against time, derived from measured stress maps such as those shown in Fig. 3 for the case when θ max = 1000 o C.

Secondly, this stage assumes that all of the measured damage occurs only at or near interfacial regions having values of b/a that are larger than some critical value r c . Hence, through the damage measurement D, a time t acc can be associated with a critical value r c . This definition of r c can be written as

( ) ∫ ∞ = c r dx x f D , ( 6
)
where f is the probability density function of b/a. Figure 10 illustrates Eq. 6 graphically, using a fitting of the data obtained from a statistical treatment of the TGO's morphological characteristics. The area under the probability density function of A more convenient way of determining the critical value r c from a given measured value of D is by replotting the information contained in Fig. 10 as the corresponding cumulative probability function. This is shown in Fig. 11, which can now be used directly to obtain the critical b/a values corresponding to every measurement of D that is available, making it possible to assemble a curve of r c against t acc for the given thermal history. It is worth pointing out that the dependence of Fig. 11 on temperature is only through the use of t acc since the relationship between D and b/a is only a property of the coating and not of the thermal cycle. The symbols in Fig. 11 show the cumulative distribution of the measured data, and the line without symbols shows the normal distribution with the same mean and variance, i.e. the integral of that shown in Fig. 10. It is important to note the differences between using the statistical normal distribution shown in Fig. 10, and the measurement data shown in Fig. 11. It is clear from Fig. 10 that the use of the normal distribution means that there is a non-zero probability of obtaining a negative value of r c . This is physically meaningless since the definitions of b and a are lengths rather than directed vectors. This problem could be removed if a different distribution (such as the log-normal) was used to describe b/a.

The output from this stage is the relationship between the critical interfacial roughness ratio r c and accumulated oxidation time t acc for a given set of θ max and θ min temperatures. This can be obtained by combining the information presented in Figs. 9 and 11. The resulting t acc vs. r c relation for the thermal cycle with θ min = 25 o C and θ max = 1000 o C is shown in Fig. 12. Thus, this curve gives the minimum value of roughness of TGO regions that will develop damage (i.e. cracks) after the TBC has been exposed to the oxidising temperature (1000 o C) for different lengths of (accumulated) time. If measurements of D generated during different thermal cycles are available, then this step can be repeated many times for other values of θ max and θ min . If enough curves are generated in this way, then interpolation between the damage versus time curves will enable the expected response at other values of θ max and θ min to be estimated.

Figure 11

Determination of the critical value of b/a from its cumulative distribution function. The line with symbols shows the measurement data, and the smooth line without symbols is the cumulative normal distribution with the same mean and variance as the measured data.

Figure 12

Critical value of b/a, r c , as a function of accumulated time for a given thermal history.

Stage 2: Interpolation and Extrapolation of the Finite Element Results

This stage combines the r c vs. t acc dependencies identified in the previous stage (Fig. 12) with the results of the parametric finite element calculations (Section 3 and The life prediction methodology proposed here was numerically implemented and accounts for the following options for the stress interpolation, depending on which stress component controls TBC failure: (i) the peak maximum principal stress within the TGO, (ii) the peak out-of-plane stress in the YSZ top coat, and (iii) the peak interfacial TGO tractions. To illustrate the methodology, and in the absence of conclusive evidence about the dominant failure mode of this TBC system, the results presented here use the peak normal traction, σ NT , on the TGO-bond coat interface (Fig. 13(a)) and the peak maximum principal stress within the TGO, σ MPS (Fig. 13(b)). As mentioned in Section 2.1, the first of these stresses drives the nucleation and propagation of interfacial cracks, and the second drives nucleation and propagation of cohesive cracks within the TGO.

The interpolation process for a chosen thermal profile proceeds as follows:

Step 1:

Inputs: Accumulated time t acc .

Step: using Eq. Output: values of h c and r c .

Step 2: Inputs: FE stress results.

Step: For each curve of peak stress vs. accumulated time obtained by FEA, use Eq. 3 and the temperature θ at which the results were calculated, to transform accumulated time into TGO thickness.

Output: Curves of peak stress vs. h for a set of values of θ and b/a.

Step 3: Inputs: Curves of peak stress vs. h from previous step, value of h c from Step 1.

Step: Interpolate the data on each curve to obtain the peak stress at a TGO thickness of h c for the values of θ and b/a corresponding to that curve.

Output: Values of peak stress at a TGO thickness h c for a set of values of θ and b/a (values lie on a 2-d surface).

Step 4: Inputs: Peak stress values from previous step, value of θ max corresponding to thermal profile associated with measured damage value D.

Step: Interpolate between the values from Step 3 with respect to the temperatures used for FE, to obtain the peak stress at a TGO thickness h c and a temperature θ max for a range of values of b/a.

Output: Values of peak stress at a TGO thickness h c and a temperature θ max for a set of values of b/a (values generate a curve).

Step 5: Inputs: Peak stress vs. b/a curve from step 4, value of r c from Step 1.

Step: Interpolate the curve from Step 5 with respect to b/a to obtain the value of the peak stress corresponding to a TGO thickness of h c , a temperature θ max , and a value of b/a equal to r c . This value is σ C .

Output: A value of σ C at a time t acc and a temperature θ max .

These five steps can be repeated for different values of t acc to build up a curve of accumulated time versus critical stress for a given thermal profile. The resulting critical stress versus accumulated time curves for the 25-1000 °C temperature history are shown in Figs. 13(a) and (b). Note that the critical stress σ C is that required to initiate a local crack for the given value of r c at the given time t acc , and is not the stress required to cause total failure at that time.

As was mentioned in Section 2.5, the FE models were run until the TGO thickness reached 3 µm, since it is expected that the coating system is likely to fail by this point. Hence the FE results are only available for TGO thicknesses less than or equal to 3 µm. If some of the values of t acc correspond to TGO thicknesses greater than 3 µm, then the software carries out an extrapolation of the FE results. In some cases this may lead to physically unlikely stress values, and so results must be interpreted with care.

Stage 3: Life Prediction

In this third and last stage the remaining life of the coating is estimated. The failure criterion used in this stage is based on the understanding that global TBC failure occurs by a cleavage-type mechanism, that is, it is stress, or more strictly speaking energy, driven. The implied assumption here is that failure is initiated from a microdefect of a given size, requiring a specific stress for the operation of the cleavage mechanism. Under this assumption, for a given temperature cycle, the failure of any given region of the TBC will occur when the local 'driving force' (i.e. stress component which controls TBC failure, denoted as σ c in Section 4.2) reaches a global failure value, σ f , independently of all other factors. It is possible to estimate σ f from measurements made on TBCs that have failed, by taking σ f to be the value of σ c corresponding to the last measurement made before failure (see ref. [START_REF] Busso | A Mechanistic Study of Oxidation-Induced Degradation in a Plasma-Sprayed Thermal Barrier Coating System. Part II: Life Prediction Model[END_REF] for an example). The global failure value σ f is a property of the material system being examined and is effectively a material property.

The inputs required in this stage are, for each temperature history of interest, the global failure stress σ f , the current value of accumulated time, t 0 , and the curve of σ c versus t acc generated as the output of the previous stage. The procedure to predict the remaining lifetime is: (i) find the accumulated time from the curve of σ c versus t acc that corresponds to the critical stress for global failure σ f as shown in Fig. 13. This time is the failure time, t fail .

(ii) subtract t 0 from t fail to obtain the remaining life of the coating.

If the value of t 0 is not known, it can be estimated from a damage measurement by interpolating the curve of damage versus accumulated time (shown in Fig. 9). It is possible that more than one value of t fail is obtained from the interpolation in (i).

It should be noted that, in principle, an alternative approach could be followed based on a measured mean stress value rather than the measured damage used so far. It will also rely on an interpolation of the finite element results to estimate an effective roughness ratio (r m ) from the measured stress, σ m , actual time t acc and the θ max of interest. However, if the stress is measured using piezo-spectroscopy, the measured stress would, in the most general case, be limited to the mean triaxiality, as discussed in Section 2.2, rather than the component responsible for failure. Thus specific failure modes within the TBC would not be able to be addressed using a measured mean stress approach.

Conclusions

A mechanism-based life prediction methodology for thermal barrier coating systems has been proposed. It relies on a combination of information about the TBC's morphological characteristics, accumulated TBC damage inferred from nondestructive fluorescence measurements, and numerically predicted local TBC stresses responsible for the initiation of such damage. It incorporates the complex interaction between interfacial and microstructural features, local oxidation mechanisms and time-dependent processes. The methodology is applied to predict the life of an EB-PVD TBC with an MCrAlY bond coat. The results of parametric finite element studies using periodic unit cell techniques revealed the magnitudes of the local TBC stresses known to lead to the failure of this type of coating. The maximum TGO stresses responsible for microcrack nucleation were found to increase with oxidation time and TGO roughness.

The proposed novel procedure to assess TBC failure has been shown to provide a physically-based yet technologically feasible method to predict the onset of TBC damage in-situ.

Appendix A : TBC Properties

A.1 Elasto-visco-plastic properties of the bond coat and oxide The steady state creep model for the bond coat and TGO is given by

n cr R Q A σ θ ε exp ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = & (A.1)
where θ is the temperature in K, and σ ~ is the uniaxial equivalent Mises stress. The values of the material properties A, n and Q for the bond coat and TGO are shown in Table A2.

= = = + = (A.2)
The values of E 2 , E p /E 2 , and ν pp are given in Tables A1, A2 and A3, respectively.

Figure 1

 1 Figure 1 SEM micrographs showing (a) a typical microstructure of a TBC system (2466 MCrAlY bond coat, EB-PVD YSZ, CMSX4 substrate) after 700 h at 1000 o C, and (b) possible initiation sites (see arrows) of interfacial microcracks between the TGO and the bond coat of a TBC system (bond coat: 2412 MCrAlY, EB-PVD YSZ, substrate: CMSX4) exposed to 1000 o C for 600 h.

Figure 3

 3 Figure 3 Typical damaged area maps measured in the current TBC system (2466 MCrAlY bond coat, EBPVD YSZ on a CMSX4 substrate): (a) photograph of a coated 20 mm diameter test sample after 700 h test at 1000 0 C, (b) typical SEM micrograph of a failing TBC system (mean TGO thickness is 3.3 µm), (c) stress maps measured by piezo-spectroscopy and (d) corresponding thermographs after 400, 600 and 700 hours exposure at 1000 0 C. (The sampling points in (c) were taken at sites in an array of cell dimensions 1x1 mm.) Typical results from the experimental measurements of TBC damage accumulation with thermal loading, conducted in this study to underpin the development of the life prediction methodology, are shown in Figs. 3(a) to (d).Figure 3(a) shows a photograph of a TBC coated superalloy sample oxidised for 700 h at 1000 o C. Delaminated areas around the edge of the specimen can be clearly seen. Figure 3(b) shows an SEM micrograph of a typical local TGO crack, which appears to have propagated along the TGO-bond coat interface. In Fig. 3(c), the results of TBC damage measured by piezo-spectroscopy are shown at three different oxidation times.Here, the grey area represents a compressive stress of around 4 GPa and black regions represent 0 GPa. Thus, the red areas are measured regions where the stresses have been relaxed by damage (or microcracks). The identification of such regions requires the specification of a threshold value for the measured mean stress

Figure 3 (

 3 a) shows a photograph of a TBC coated superalloy sample oxidised for 700 h at 1000 o C. Delaminated areas around the edge of the specimen can be clearly seen. Figure3(b) shows an SEM micrograph of a typical local TGO crack, which appears to have propagated along the TGO-bond coat interface. In Fig.3(c), the results of TBC damage measured by piezo-spectroscopy are shown at three different oxidation times. Here, the grey area represents a compressive stress of around 4 GPa and black regions represent 0 GPa. Thus, the red areas are measured regions where the stresses have been relaxed by damage (or microcracks). The identification of such regions requires the specification of a threshold value for the measured mean stress 12µm 20 mm

Figure 4

 4 Figure 4 Experimental data showing the effect of oxidising temperature on the time to TBC spallation for MCrAlY type coatings and for PtAl based TBCs (triangles). The line represents the time at a given temperature to reach a TGO thickness of 3 µm, calculated using Eq. 3.

Figure 5

 5 Figure 5 Illustration of the surface roughness geometric parameters a and b and a typical finite element mesh of a periodic unit cell used in the modelling.

3 •

 3 b = 0.79 µm, b/a = 0.20, • b = 2.47 µm, b/a = 0.32, and (4) • b = 4.54 µm, b/a = 0.48.

Figures 6 (

 6 Figures 6 (a) to (c) shows the contour plots of the out-of-plane stress component (σ 22 ) when the TBC was cooled down to 25 o C after 241 h oxidation at 1000 o C. The results are shown for TBCs with the three different TGO morphologies specified in Eq. 4.Note that the TGO thickness is 3 µm in all the plots. Similar results were obtained for oxidation at 1100 o C, but contours at this higher temperature are not shown in this paper for brevity. It can be seen that the extent and magnitude of the most highly stressed areas increases with interface roughness, and that these regions are located in the YSZ top coat around the valleys. It is also worth noticing that both the TGO and the YSZ regions tend to be generally under some compression at 25 o C for the rougher TBCs. Thus these results show that, for

Figure 6

 6 Figure 6 Three typical out-of-plane stress contour patterns at 25 o C for different geometric parameters with θ max = 1000 ºC: (a) b =0.79 µm, b/a = 0.20, (b) b= 2.47 µm, b/a = 0.32 and (c) b = 4.54 µm, b/a = 0.48. In all cases, the TGO thickness is 3 µm.

Figure 7

 7 Contour plots of the out-of-plane stress component σ 22 after cooling to room temperature from 1000 ºC for (a) 61 h, (b) 167 h, and (c) 241 h oxidation. The corresponding TGO thicknesses are 1.91, 2.65 and 3.00 µm, respectively. (Also, b = 2.47 µm, b/a = 0.323).

Figure 8 :

 8 Figure 8: Flow chart of the life prediction methodology.

  Fig. 10 corresponding to a measured value of D when time is equal to t acc , enables the estimation of the critical roughness value r c .

Figure 10

 10 Figure 10 Critical value r c as determined from the statistical normal distribution of the b/a values.

Figure 13

 13 for each thermal cycle to obtain a relationship between the accumulated time t acc and the critical value of the local stress responsible for local microcrack initiation and growth within the TBC, denoted by σ C . This stress is regarded as critical because it is associated with the measured values of the damage parameter D via its dependence on r c . The results of the parametric analyses are in the form of curves of the peak values of the stress component which controls TBC failure versus accumulated time for given values of b/a and θ max . By interpolating between these curves with respect to b/a, θ max , (assuming θ min remains always the same), and accumulated time (or TGO thickness), the value of the stress σ C at a time t acc for the corresponding critical value r c and temperatures θ max and θ min can be found. Illustration of interpolation to obtain TBC failure time, t fail , and remaining lifetime for a coating whose failure is controlled by (a) the peak normal traction on the TGO-bond coat interface and (b) peak maximum principal stress within the TGO.

  3 for h in terms of t acc and θ max , and values of r c vs. t acc as in Fig 12, determine the values of TGO thickness (h c ) and r c for the time t acc and the temperature θ max .

Table A1 :

 A1 Elastic-plastic material properties of the bond coat and oxide

	Temperature		Bond coat			TGO
	(°C)	E (GPa)	ν	α (×10 6 , K -1 )	σ Y (MPa)	E (GPa)	ν	α (×10 6 , K -1 )
	20	200	0.3	13.6	426	400	0.23	8.00
	200	190	0.3	14.2	412	390	0.23	8.20
	400	175	0.31	14.6	396	380	0.24	8.40
	600	160	0.31	15.2	362	370	0.24	8.70
	800	145	0.32	16.1	284	355	0.25	9.00
	1000	120	0.33	17.2	202	325	0.25	9.30
	1100	110	0.33	17.6	114	320	0.25	9.60

Table A2 :

 A2 Creep properties of the bond coat and oxide Effect of sintering time and temperature on the transversely anisotropic elastic moduli of YSZ The elastic properties of the transversely isotropic YSZ are here fully defined by E 2 , E p /E 2 (function of sintering time only), and ν pp (function of temperature only), where E p is the isotropic in-plane Young's modulus, E 2 is the out-of-plane Young's modulus, and ν pp is the Poisson's ratio linking the in-plane directions. The other elastic constants are given by

														Bond coat				TGO
							n								4.1					2.3
				Q (kJ/ (kg mol))				263×10 3			424×10 3
				A (1/(s MPa -n ))					13.3				6805
	A.2 ( G p 1 2	E	p ν	pp	)	,	ν	2	p	ν	p	2	E E	p 2	,	G	2	G	p	E E	p 2	,	ν	p	2	ν	pp	E E	p 2	.

Table A3 :

 A3 E 2 (GPa) as a function of temperature and sintering time

	Time			Temperature (°C)
	(hours)	25	200	400	600	800	1000	1200
	0	183	175	168	163	159	155	153
	20	196	187	179	174	169	166	164
	60	203	194	186	180	176	172	170
	100	205	196	188	182	178	174	171
	200	207	198	190	184	179	176	173
	400	209	199	191	185	180	177	174
	600	209	200	192	186	181	177	175

Table A4 :

 A4 E p /E 2 as a function of sintering time

	Time (hours)	0	20	60	100	200	400	600
	E p /E 2	0.01	0.169	0.258	0.281	0.303	0.318	0.325

Table A5 :

 A5 ν pp as a function of temperature

	Temperature	25	200	400	600	800	1000	1200
	(°C)
	ν pp	0.1	0.1	0.1	0.11	0.11	0.12	0.12

Appendix B: Summary of the predicted effects of the geometric parameters and peak cycle temperature on the stress maxima relevant to the different TBC failure modes at room temperature.


Anna Maria Vegni
email: annamaria.vegni@uniroma3.it

Valeria Loscrí
email: valeria.loscri@inria.fr

Chirality effects on channel modeling for THz-band wireless communications in LoS/NLoS propagation

Keywords: Chirality effects, full-wave model, THz channel modeling, nanoscale communications

Chirality effects on channel modeling for THz-band wireless communications in LoS/NLoS propagation

Anna Maria Vegni, Valeria Loscrì

Introduction

The increasing demand of wireless communication systems based on largescale use of nanodevices according to the Internet of Things idea has enhanced the attention on accurate and more precise channel models working in the Terahertz frequency band. This was a consequence of new sensing components to be used in a large variety of applications, including e.g. health and life sciences.

Nowadays, the Terahertz band [START_REF] Piesiewicz | Short-Range Ultra-Broadband Terahertz Communications: Concepts and Perspectives[END_REF] has been identified in the range [0.06, 10] THz, and represents one of the most promising spectrum bands to enable ultra-highspeed communications. The use of this frequency band is envisioned to address the spectrum scarcity and capacity limitations of current wireless systems, thus providing a plethora of applications, like ultra-high-speed wireless communications.

In order to address these needs for an optimal wireless communication networks in the THz band, it is necessary to develop new channel models working in THz band, that include specific phenomenons that are more evident when working in optic and almost-optic frequencies [START_REF] Kurner | Towards THz Communications -Status in Research, Standardization and Regulation[END_REF]. Researchers are investigating specific Terahertz spectrum peculiarities, like the very high molecular absorption loss or the very high reflection loss. Indeed, differently from traditional lower frequency bands, where the propagation is mainly influenced by the spreading loss only, the physical mechanisms in a THz-band wireless transmission are a very high molecular absorption loss, and spreading loss. Also, the molecular absorption attenuates the propagated signals in THz-band wireless communications. These phenomenons are specific for THz band and cause a very high and frequencyselective path loss for Line-of-Sight (LoS) links, while high reflections occur in Non-Line-of-Sight (NLoS) propagation due to the type of material and the roughness factor of reflecting surfaces [START_REF] Kurner | Towards THz Communications -Status in Research, Standardization and Regulation[END_REF][START_REF] Moldovan | LOS and NLOS channel modeling for terahertz wireless communication with scattered rays[END_REF][START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF].

From all these considerations, it follows that existing channel models used for lower frequency bands cannot be re-used for THz communications, and then the need to design novel channel models for THz band is a challenge. In this paper, we omit channel models for lower frequency (i.e., ultra-wideband), but investigate recent works in channel modeling for THz band.

Many works have already investigated all these specific phenomenons occurring in the THz band, and consequently novel channel models for THz band have been presented so far. One of the main channel models for communications in THz band is that presented by Jornet and Akyildiz in [START_REF] Jornet | Channel Modeling and Capacity Analysis for Electromagnetic Wireless Nanonetworks in the Terahertz Band, Wireless Communications[END_REF], which includes both molecular absorption and free space loss. In [START_REF] Llatser | Time-and Frequency-Domain Analysis of Molecular Absorption in Short-Range Terahertz Communications, Antennas and Wireless Propagation Letters[END_REF], Llatser et al. investigate the influence of molecular absorption in a short-range THz communication model on both frequency and time domain. Piro et al. [START_REF] Piro | Terahertz Communications in Human Tissues at the Nanoscale for Healthcare Applications[END_REF] address radio propagation models for nanoscale wireless communications in the field of nanomedicine. Finally, Boronin et al. [START_REF] Boronin | A molecular noise model for THz channels[END_REF] present a model for describing the phenomena of molecular noise occurring in THz band. This concept has been broadly surveyed in [START_REF] Kokkoniemi | A discussion on molecular absorption noise in the terahertz band[END_REF] by distinguishing from time and frequency dependence.

In [START_REF] Moldovan | LOS and NLOS channel modeling for terahertz wireless communication with scattered rays[END_REF] Moldovan et al. present a deterministic channel model working in THz band, both in case of LoS and NLoS propagation, and by assuming molecular absorption and spreading loss, together with a ray-tracing approach. In [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF] Han Another feature to be considered is the scattering effect. Piesiewicz et al. [START_REF] Piesiewicz | Scattering Analysis for the Modeling of THz Communication Systems, Antennas and Propagation[END_REF] study scattering properties of typical building materials, for channel modeling in THz band. The authors observed that in realistic environments, with rough surfaces, the modeling of the propagation effects should take into account the reduction of the reflected power in the specular direction due to diffuse scattering, even if not extremely large. Also, robust scattering models are needed and can be implemented in ray-tracing algorithms to enhance the propagation modeling. Finally, Kokkoniemi et al. [START_REF] Kokkoniemi | Frequency and Time Domain Channel Models for Nanonetworks in Terahertz Band, Antennas and Propagation[END_REF] consider the scattering effect into the channel modeling in the THz band. They theoretically analyze multiple scattering effects, referring to a sequence of scattering events from small particles, such as aerosols.

It is shown that the small particle scattering can result into significant additional loss that needs to be taken into account with the loss depending on the density and size distribution of the particles.

As can be noticed from previous works, the channel modeling in THz band is addressed through transfer functions that consider specific effects, like molecular absorption loss, spreading loss, reflected, diffracted, and scattered rays. However, the electromagnetic chirality effect and the particular kind of the medium (also in terms of power excitation, with respect to the intensity of the electromagnetic field) are a few aspects that are often neglected in classical modeling of transmission channels. In fact, related works on channel modeling in the THz band consider the behavior of the channel by assuming molecular absorption loss function, the signal attenuation due to free space, and the particle scattering effect. All those aspects are assumed through appropriate functions (e.g. the Friis's function, the directivity of the transmitting antenna) added to the amplitude of the transmitted electromagnetic field. These features can be addressed in terms of constitutive relations of the electromagnetic medium, where the propagation of the electromagnetic signal occurs. On the other side, in this paper, we aim to consider a punctual setup of the issue, based on the determination of whole structure of the electromagnetic field that propagates in the channel, and then, a more appropriate definition of the constitutive relations of the medium should result important for this aim.

To summarize, the main goals addressed in this paper are (i) the full-wave propagation model in a non-conventional material (namely, a chiral metamaterial), and (ii) how the chirality effects affect the channel transfer functions on LoS and NLoS propagation. Chiral metamaterials, namely metamaterials that do not present any mirror symmetry plane, are materials where electric (magnetic) field generates also magnetic (electric) polarization [START_REF] Sihvola | Metamaterials in electromagnetics[END_REF]. They represent a very interesting category of metamaterials, since they offer great possibilities in the control of polarization such as optical activity. Specifically, they allow the polarization rotation of a linearly polarized wave.

As outlined in [START_REF] Kenanakis | Controlling THz and far-IR waves with chiral and bianisotropic metamaterials[END_REF], another aspect really interesting of chiral metamaterials is the exploitation of chirality property in order to reduce the attractive Casimir force, which is detrimental in micro-and nano-electromechanical systems. It has been discovered that chiral metamaterials are very good to realize the repulsive Casimir force, and this force is strictly dependent on the strength of the chirality. In particular, the repulsive force leads to nanolevitations and opens up many opportunities such as frictionless operation of nanomotors. Furthermore, chiral metamaterials can provide negative refractive index, and this type of materials has gained a lot of attention in the last few years, since they can be used in all the applications that require "perfect lenses", by overcoming diffraction limit. This paper is organized as follows. In Section 2 we introduce a full-wave propagation model of the electromagnetic field that propagates in the THz band. Starting from the concept of electromagnetic chirality [START_REF] Lindman | Om en genom ett isotropt system av spiralformiga resonatorer alstrad rotationspolarisation av elektromagnetiska vagorna, Ofversigt af Finska Vetenskaps-Societetens forhandlingar[END_REF][START_REF] Engheta | Chiral Materials and Chiral Electrodynamics: Background & Basic Physical Principles[END_REF][START_REF] Lakhtakia | Recent contributions to classical electromagnetic theory of chiral media: what next?[END_REF][START_REF] Lindell | Electromagnetic waves in chiral and bi-isotropic media[END_REF], this approach considers the chiral effects following the change in the propagation velocity and in the refractive index, due to the chiral impurities inside the propagation medium.

These effects are evaluated also in the case when the considered medium exhibits a Giant Optical Activity (GOA) [START_REF] Fang | Dual-band Terahertz Chiral Metamaterial with Giant Optical Activity and Negative Refractive Index based on Cross-wire Structure[END_REF][START_REF] Zhou | Terahertz chiral metamaterials with giant and dynamically tunable optical activity[END_REF]. Section 3 describes the equivalent channel model by considering the chirality-affected channel transfer functions, both for LoS and NLoS propagation. In Section 4 the performances have been assessed in terms of path loss in cases of LoS and NLoS propagation mechanisms, for biological and indoor applications. Finally, conclusions are drawn at the end of the paper.

Full-Wave Propagation in THz Band

This section is devoted to the design of a full-wave propagation model for THz band. Starting from the classic harmonic macroscopic Maxwell's equations, we consider the electromagnetic propagation inside a generic complex material, under the assumption that it is a linear and chiral medium.

Maxwell's equations represent the differential relations between the electric displacement D and the electric field E, as well as the magnetic displacement B and the magnetic field H, combined together with the constitutive relations characterizing the medium where the electromagnetic field exists according to assigned electric and magnetic sources.

For real-world materials, the constitutive relations are rarely simple, except approximately, and usually determined by experiment. For materials without polarisation and magnetisation effects (i.e., vacuum), the constitutive relations are:

B = µ 0 H, D = ε 0 E, (1)
where ε 0 is the electric permittivity, and µ 0 is the magnetic permeability.

For homogeneous materials, the electric permittivity is ε = ε 0 ε r , where ε r is the relative electric permittivity. Analogously, the magnetic permeability changes as µ = µ 0 µ r , where µ r is the relative magnetic permeability. These quantities are constant throughout the material, while in inhomogeneous materials they depend on space location within the material. For isotropic materials, ε and µ are scalar quantities, while for anisotropic materials they are tensors e.g., in crystal structures.

In this formulation unconventional materials (i.e. metamaterials) are also included, even if our considerations will be focused essentially on traditional dielectrics. Materials can be also generally dispersive, so that ε and µ depend on the frequency of any incident electro-magnetic (EM) waves. Even more generally, in the case of non-linear materials B and D are not necessarily proportional to E, and H, respectively. In general, D and B depend on both E and H, on location and time.

Now, let us consider a time-harmonic generic linear material, where chiral (magneto/electric-optical) effects are included in the following constitutive relations:

   B = ξ • E + µ • H D = ε • E + ζ • H (2)
where the symbol • represents the scalar product operator. In (2), ξ, µ, ε, and ζ are specific tensor quantities of the material.

From (2) we observe the chirality property through the dependence of (i) E in B, and (ii) H in D. Furthermore, the displacement field existing inside the material is generated by an excitation expressed in terms of intensity of the incident electro-magnetic field. Therefore, the material under consideration is a linear chiral medium.

We remind that the chiral effects mean that:

• an electric field applied on the material provides not only an electric induction, but also a magnetic displacement,

• a magnetic field applied on the material provides not only a magnetic induction, but also an electric displacement, unlike other non-chiral materials.

We assume this specific material with the aim of assessing the performance (i.e., transmission of the electromagnetic power) of an optical channel. This performance is linked to the channel transfer function. Also, this formulation allows a more rigorous assessment of the channel properties, that is verified for all the frequency range (i.e., from 0 to infinity).

For this aim, in our proposed approach we need to define the expressions of the EM field generated by a transmitter that is assumed as an impressed source distribution in a given volume V , characterized by assigned electric and magnetic vector functions J i,e (r) and M i,m (r) 1 , respectively, where r is the observation vector. In the following, the spatial dependence r will be omitted.

The computation of the EM field in the channel is carried out through the solutions of linear differential equations, arising from the following Maxwell's equations:

   ∇ × E = -jωB -M i,m ∇ × H = jωD + J i,e (3)
where Kong's vector ∇ in Cartesian space has been used i.e.,

∇ = ∂ x x + ∂ y ŷ + ∂ z ẑ = ∂ ∂x x + ∂ ∂y ŷ + ∂ ∂z ẑ (4)
The system in (3) can be written as in the following form:

   ∇ + jωξ • E = -jωµ • H -M i,m ∇ -jωζ • H = jωε • E + J i,e (5)
where Kong's operator ∇ is defined as a 3 × 3 matrix:

∇ =      0 -∂ z ∂ y ∂ z 0 -∂ x -∂ y ∂ x 0      (6)
and the solving differential equations for E and H are given, respectively:

   ∇ -jωζ • µ -1 • ∇ + jωξ -ω 2 ε • E = f e (J i,e , M i,m ) ∇ + jωξ • ε -1 • ∇ -jωζ -ω 2 µ • H = f m (J i,e , M i,m ) (7)
Equation ( 7) derives from well-known Maxwell's equations [21]. In this paper, we apply classical electromagnetic theory to the THz band, by assuming a signal generated by a linear polarized electric field propagating in a complex medium described by ξ, µ, ε, and ζ. Thus, from the authors' best knowledge, the following equations represent novel results, useful to derive the electromagnetic field in a complex medium. However, the proposed equations can be appropriately applied to both traditional and non-conventional (non-traditional) materials, under the assumption of specific material constitutive relations.

The solutions of the problem can be obtained via ( 7), together with the divergence equations. Following this approach, we obtain the expression of the total electromagnetic field, but in order to give a simplified formulation for every regular observation point P, that is ∀P / ∈ V , we pose the condition of electric linear polarization (e.g., x). Under this hypothesis, we get:

   E = E 1 x H = -(ρ R1 • ξ R1 E 1 ) x + j ρ 23 ω ∂E 1 ∂z -ρ R2 • ξ R1 E 1 ŷ + j ρ 33 ω ∂E 1 ∂y + ρ R3 • ξ R1 E 1 ẑ (8)
where ρ R1 , ρ R2 , and ρ R3 are the row vectors of ρ = µ -1 tensor, ξ R1 is the first row vector of ξ tensor, and ρ 23 , ρ 33 are elements of ρ tensor. Notice that E satisfies the following partial differential equation:

f 1 (r) E 1 +f 2 (r) ∂E 1 ∂x +f 3 (r) ∂E 1 ∂y +f 4 (r) ∂E 1 ∂z +f 5 (r) ∂ 2 E 1 ∂y 2 +f 6 (r) ∂ 2 E 1 ∂z 2 = 0, (9)
where functions f i (r) with i = {1, 2, . . . , 6} depend on the elements of ξ, µ, ε, and ζ tensors, and assume complicated expressions. As an instance, we present the f 1 (r) term under the condition of spatial homogeneity of ξ, µ, ε, and ζ tensors:

f 1 (r) = (ρ 23 + ρ 32 ) det Ω 1,1 + ω 2 v 1 • ξ C1 -ε 11 det Ω 1,2 , (10)
with

Ω 1,1 =      -ρ 32 ρ 33 ω 2 ζ R2 • ρ C1 ξ 11 + ζ R2 • ρ C2 ξ 21 + ζ R2 • ρ C3 ξ 31 -ε 21 ρ 22 -ρ 23 ω 2 ζ R3 • ρ C1 ξ 11 + ζ R3 • ρ C2 ξ 21 + ζ R3 • ρ C3 ξ 31 -ε 31 ζ R1 • ρ C2 -ζ R1 • ρ C3 0      , (11)
Ω 1,2 =      ρ 13 ρ 32 -ρ 33 -ρ 12 ρ 22 -ρ 23 ζ R2 • ρ C2 -ζ R3 • ρ C3 ζ R1 • ρ C2 -ζ R1 • ρ C3      , (12)
v 1 = ζ R1 • ρ C1 x + ζ R1 • ρ C2 ŷ + ζ R1 • ρ C3 ẑ, (13)
where the indexes C1, C2, and C3 represent the column vectors.

According to this approach, the Poynting's vector defined as:

P = E × H * 2 , (14)
can be split into two components i.e., P = P y ŷ + P z ẑ, whose expressions are

   P y = -E 1 2 j ρ 33 ω ∂E 1 ∂y + ρ R3 • ξ R1 E 1 * P z = E 1 2 j ρ 23 ω ∂E 1 ∂z -ρ R2 • ξ R1 E 1 * (15)
From ( 15), we can easily distinguish that Poynting's vector is formed by two components that are orthogonal to each other (i.e., P y and P z ), and to the electric polarization axis (i.e., x). Therefore, if we refer this solution to the case of a transmission channel, where a giant optical activity could be activated, we can deduce that one of these terms is responsible of LoS power (W LoS ), that is the power at the receiver traveling along the direct (LoS) path, and the other one contributes to the NLoS power (W NLoS ), that is the power at the receiver traveling along the indirect (NLoS) path, according to the receiver location. Obviously, similar considerations hold for a ŷ/ẑ-axis electric polarization.

In the following, we attribute the presence of the indirect path as due at reflection/scattering/diffraction centers located at z = h. Therefore, the reflection/scattering/diffraction characteristics of the transmissive channel can be evaluated through the specific knowledge of the local planar geometry associated to these centers. In this paper, we consider the effects of the reflection generated by these centers allocated in the generic S points, while diffraction coefficient requires the knowledge of the canonical geometry of the diffraction points. This consideration limits the cases of real geometry that can be examined, and imposes a statistical approach to the problem more than a deterministic one we are dealing with. Finally, scattering effects can be inserted either at the level of constitutive relations or added in dB as a global term to the power flow evaluation.

By assuming Figure 1, where the transmitter and the receiver are aligned along the ŷ-axis, we can derive that the LoS power (i.e., W LoS ) is determined by P y , while the NLoS power (i.e., W NLoS ) is contributed both by P y , and P z .

ŷ ẑ Tx Rx θ π 2 -θ h P y ŷ P z ẑ ẑ d 1 d 2 β i β r β i β r Q ≡ 0,η,h (
So, at the receiver the W LoS can be written as

W LoS = S R n • P y ŷdS = = 1 2 S R j ρ * 33 ω E 1 ∂E * 1 ∂y -ρ * R3 • ξ * R1 E 1 E * 1 dS, (16)
being n the normal vector to the S R receiving surface i.e., n = -ŷ.

We observe that ( 16) is a general formula that can be applied in the case of conventional and/or unconventional materials. In this way it can be applied either to biological propagation media, as well as indoor wireless communications.

Moreover, ( 16) is expressed in terms of the electric field and depends on the ele-ments of the constitutive tensors. The explicit representation of the electric field can be obtained after integration of ( 9), as it will be indicated below.

On the other side, in order to obtain the W NLoS , we need to evaluate the indirect path followed by the Poynting's vector through the generic reflection/ scattering/ diffraction point S laying in the z = h plane, to reach the receiving point.

We distinguish these points S into two classes i.e., (i) Q points belonging to the straight line (x = 0, z = h), and (ii) all the remaining R points in the z = h plane. In this paper, we refer to the first class Q points, as shown in Figure 1.

Let us consider the effect of the local mirror reflection in points Q in z = h. In this generic reflection point, the dielectric discontinuities in the permittivity give a local specular reflection coefficient that acts on the two incident Poynting's vectors From Figure 1, it is easy to verify that for Q points belonging to the z = hplane, the incident/reflected β angle, related to the incident Poynting's vector, is connected to the elevation angle θ by the following simple relation:

P (i) components i.e., P (i) = P y (i) ŷ + P z (i) ẑ,
cos θ cos(θ -2β) = d 2 d 1 , (17)
where d 1 and d 2 are the distances of these points from the transmitter and the receiver, respectively. Notice that for R points, a similar formula for the NLoS ŷ-axis Poynting's vector component can be established as well.

Let us now determine P (r) in Q points. With reference to Figure 2 (a), we can

Q z = h γ π 2 -γ P z i ( ) P y i ( ) 2γ π 2 -2γ Q z = h γ P y,
γ n Q ( ) n Q ( ) P y r ( ) (a) (b)
Figure 2: Schematic of the reference geometry for P (r) evaluation, in the case of (a)

P (i)
y , and (b)

P (i) z components.
write:

P (r) Q = P y (r) Q ŷ + P z (r) Q ẑ = = -P y,y (r) + P z,y (r)
Q

ŷ -P y,z (r) + P z,z (r)
Q ẑ = = P y NLoS Q ŷ + P z NLoS Q ẑ, (18)
and from the Snell's law, we obtain:

P y (r) Q = -P y (i) Q cos (2γ) + P z (i) Q sin (2γ) c (γ) = P y NLoS Q , (19)
and

P z (r) Q = -P y (i) Q sin (2γ) + P z (i) Q cos (2γ) c ⊥ (γ) = P z NLoS Q , ( 20
)
where c (γ) and c ⊥ (γ) are the reflection coefficients on the local plane reflecting surface centered on Q point for parallel and orthogonal incidence, respectively.

From ( 19) and ( 20), we observe that the components along ŷ and ẑ of the incident Poynting's vector contribute to the mirror reflection in point Q. The new reflection components i.e., P

y,y , P

y,z , P

z,y , and P (r) z,z , propagate towards the receiver with the same propagation characteristics of the incident components, as the propagation medium is the same of the direct path (LoS path). Therefore, such components will arrive at the receiving end with an amplitude proportional to that of the incident waves, as influenced from the reflection in point Q, but with the same propagation characteristics.

To determine W NLoS for the generic Q point, we observe that if the receiving surface is (-ŷ)-oriented, only P

y,y and P (r)

z,y will contribute, while W LoS is determined via P y (see Figure 1). These two quantities can be written implicitly as functions of the electric field, but to give an explicit form of these, the partial differential equation in [START_REF] Kokkoniemi | A discussion on molecular absorption noise in the terahertz band[END_REF] satisfied by the electric field must be solved. This can be done by writing the electric component as a factorized expression that is

E 1 = e 1 (x) e 2 (y) e 3 (z) . (21)
Under this hypothesis, the expression of W LoS in ( 16) becomes:

W LoS = j e 2 (d) 2ω de * 2 (y) dy y=d N 1 - |e 2 (d)| 2 2 N 2 , (22)
where N 1 and N 2 are two scalar quantities i.e.,

       N 1 = x 2 x 1 z 2 z 1 ρ * 33 |e 1 (x)| 2 |e 3 (z)| 2 dxdz N 2 = x 2 x 1 z 2 z 1 ρ * R3 • ξ * R1 |e 1 (x)| 2 |e 3 (z)| 2 dxdz (23)
where x 1 , x 2 , z 1 , and z 2 are the extension limits of S R .

Finally, to determine W NLoS at the receiving end, we observe that it can be written as the sum of infinite individual contributions of Q and R reflecting points P NLoS y,k , but can be approximated to N T element with k = {1, 2, . . . , N T }, over S R :

W NLoS = N T h=1 S R n • P h (r) dS = N T h=1 S R n • P NLoS y,h ŷ + P NLoS z,h ẑ dS. ( 24
)
To give a simplified example, we can consider the case of a single reflecting Q point (i.e., h = 1), by assuming again that the x-polarized electric field is written in a factorized form. Following this hypothesis, we obtain:

W NLoS h=1 = e 2 (d) x 2 x 1 z 2 z 1    j ω ρ * 33 de * 2 (y) dy y=d + ρ * 23 de * 2 (y) dz y=d + -ρ * R3 + ρ * R2 • ξ * R1 e * 2 (d)    |e 1 (x)| 2 |e 3 (z)| 2 dxdz. (25)
Finally, from ( 16) and (25), the total power W , due to the contribution of a single generic reflection point, at the receiver is computed as

W | h=1 = W LoS + W NLoS h=1 .
We can conclude that all the remarks regarding the specific features of chiral metamaterials permit us to characterize them as perfect candidates for nanosystem applications realization, and more specifically, they represent the ideal candidates for operations in the THz and also higher frequencies. In fact, in comparison to natural materials, chiral metamaterials show strong response to the THz radiation that shows a very great technological potential in several sectors such as imaging, sensing and also communications. Also, materials exhibiting a giant tunable optical activity show a very great potential in the achievement of active THz polarization components (such for example tunable polarizers, polarization filters, etc.). As outlined by Kenanakis et al. in [START_REF] Kenanakis | Optically controllable THz chiral metamaterials[END_REF], giant optical activity is a very important property particularly in the THz regime where there is a lack of optical components.

Equivalent Channel Model

In this section, we will present a channel model for THz band, based on ray tracing technique, and that considers the chirality effects of the EM field by taking advantage of the results in Section 2. Indeed, we firstly derived the received electric power in LoS and NLoS propagation in a chiral metamaterial from Maxwell's equations for the specific material. In this section we keep on exploiting Maxwell's equations to obtain the chirality dependence of the channel transfer functions in LoS and NLoS scenarios. More details on how the chirality affects the channel transfer functions are provided in the Appendix at the end of the paper.

The use of ray tracing techniques for channel modeling in THz band has been already adopted in [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF], where the authors present a multi-ray propagation model consisting of LoS, reflected, scattered, and diffracted paths. Indeed, due to the very short wavelength in the THz band, the ray tracing approach allows accurate modeling. We refer to the work in [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF], and consider the chirality effects on the channel impulse response, affecting both the LoS and NLoS ray propagation in the simplified case of spatial homogeneous chirality. The channel model will be the combination of many individual narrow sub-bands with a flat-band response. In the i-th frequency sub-band, the narrow-band channel impulse response is expressed as a superposition of N i rays, among which the n-th ray experiences frequency-dependent attenuation α i,n that is also affected by chirality effects of the propagation medium.

In the case of stationary environment, with a fixed transmitter and receiver, the channel response of the multi-ray model is given as

h i (τ ) = N i n=1 α i,n δ (τ -τ n ), (26)
where τ and τ n are the propagation delays for the direct ray, and the n-th NLoS ray, respectively.

Analogously in [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF], by assuming a homogeneous chirality-affected LoS propagation, and

N (i)
Ref reflected rays, the multi-ray channel model in the i-th frequency sub-band can be described as:

h i (τ ) = α i,LoS δ (τ -τ LoS ) I LoS + N (i) Ref l=1 α (Ref) i,l δ τ -τ (Ref) l , ( 27
)
where I LoS is the indicator function that is equal to 1 or 0 for the presence of LoS path or not, respectively. For the LoS path, α LoS and τ LoS refer to the attenuation and the delay, respectively. For the l-th reflected path, α

(Ref) i,l
, and τ The chirality-affected LoS channel transfer function can be expressed as:

H LoS (f ) = H Abs (f ) H Spr (f ) e -j2πf τ LoS , ( 28
)
where H Abs is the transfer function due to the molecular absorption loss that accounts for the attenuation that part of the wave energy is converted into internal kinetic energy to the molecules in the propagation medium [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF], i.e.

H Abs = e -1 2 k(f )d , (29)
while H Spr is the transfer function due to the spreading loss, when the chirality is active i.e.,

H Spr = ν Chir 4πf d , (30)
by assuming d as the distance between the transmitter and the receiver. In (30), ν Chir is the propagation velocity of the EM field in a chiral homogeneous isotropic medium i.e.,

ν Chir = c n Chir , (31)
where n Chir is the refractive index, and being c = 1/ √ ε 0 µ 0 , we can write

n Chir = µ r ε r + ξ 2 cr . ( 32
)
where ξ cr is the relative chirality parameter, and ε r (µ r ) is the relative electric (magnetic) permittivity (permeability), respectively. From (32), the expression in (31) becomes

ν Chir = c ε r µ r + ξ 2 cr . (33)
Finally, it follows that τ LoS = d/ν Chir is the time-of-arrival of the direct ray.

For the NLoS case, by considering d 1 as the distance between the transmitter and a generic reflecting point Q, and d 2 as the distance between Q and the receiver, the frequency-dependent transfer function of the NLoS ray propagation affected by the homogeneous chirality effect, i.e. H NLoS , is given by [START_REF] Han | Multi-Ray Channel Modeling and Wideband Characterization for Wireless Communications in the Terahertz Band, Wireless Communications[END_REF] H

NLoS (f ) = ν Chir 4πf (d 1 + d 2 ) e -j2πf τ Ref -1 2 k(f )(d 1 +d 2 ) • R (f ) , (34)
where

τ NLoS = τ LoS + (d 1 + d 2 -d)/ν
Chir is the time-of-arrival of the reflected ray, and R(f ) is the rough surface reflection loss of EM waves at THz band. It depends on the material, the shape and the roughness of the surface on which EM field has been reflected.

According to the Kirchhoff scattering theory, the reflection coefficient for a rough surface can be obtained by multiplying the smooth surface reflection coefficient derived from the Fresnel equations, Ψ(f ), with the Rayleigh roughness factor ς(f ), as

R(f ) = Ψ(f )ς(f ), (35)
where Ψ is chosen according to the polarization of the incident Poynting's vector, and can be approximated as:

Ψ (f ) ≈ -exp -2 cos (β i ) n 2 Chir -1 . ( 36
)
and the Rayleigh roughness factor is defined as:

ς (f ) = exp - 8π 2 f 2 σ 2 cos 2 (β i ) ν 2 Chir , (37)
with σ as the rough surface height standard deviation, commonly assumed as Gaussian-distributed.

By combining the aforementioned models in (28) and ( 34), the channel model in the i-th frequency sub-band in (27) becomes as:

h i (τ ) = ν Chir 4πf i d e -1 2 k(f i )d δ (τ -τ LoS ) I LoS + + N (i) Ref l=1 ν Chir 4πf i (d 1 +d 2 ) e -1 2 k(f i )(d 1 +d 2 ) e - 2 cos(β i ) √ n 2 t,Chir -1 e - 8π 2 f 2 i σ 2 cos 2 (β i ) ν 2 Chir •δ τ -τ (Ref) l . (38)
Now, we are able to compute the total path loss, i.e. A [dB], both in the case of LoS and NLoS propagation. In the first case, it is contributed by the spreading path loss (i.e., A s ), and the absorption path loss (i.e., A a ), and it can be expressed as:

A LoS = A s + A a = 20log 10 4πd λ Chir + 10γdlog 10 e, (39)
where λ Chir = ν Chir /f is the wavelength in the considered medium affected by the homogeneous chirality, and γ is the absorption coefficient measuring the amount of absorption loss of the EM field in the medium.

On the other hand, from (34) we can derive the total path loss for NLoS propagation i.e., A NLoS , as follows:

A NLoS = 20log 10 ν Chir 4πf (d 1 +d 2 ) + +20log 10 e -1 2 α(d 1 +d 2 ) e - 2 cos(β i ) √ n 2 Chir -1 e - 8π 2 f 2 σ 2 cos 2 (β i ) ν 2 Chir . (40)

Simulation Results

In this section, we distinguish simulation results for LoS and NLoS propagation models, in the case of (i) biological (for very small distances, in the order of millimeters), and (ii) indoor (for small distances, in the order of meters) wireless communications, and w.r.t. the chirality effects on these two types of propagation models. Since in this paper our aim is basically addressed to the path loss behavior

Scenario Parameters

LoS In NLoS propagation case, we have a similar behavior as those observed in LoS propagation. Figure 5 regarding the propagation in bio-media shows small variations for different values of ε r , with levels of attenuation higher than those exhibited in the LoS propagation scenario. In fact, in Figure 5 (a) we observe a path loss ranging around [30, 90] dB, while values increase for higher distances as Finally, we investigate the path loss for LoS and NLoS scenarios in media that exhibit high values of ξ Chir , that is when the environment where the electromagnetic field is propagating is a chiral unconventional material (Chiral Meta-Material). In Figure 7 we observe a reduction of the attenuation in the band [0.06, 10] THz w.r.t the remaining (higher) THz band, but with higher levels of path loss. Furthermore, it is observed [START_REF] Fang | Dual-band Terahertz Chiral Metamaterial with Giant Optical Activity and Negative Refractive Index based on Cross-wire Structure[END_REF][START_REF] Zhou | Terahertz chiral metamaterials with giant and dynamically tunable optical activity[END_REF][START_REF] Zhao | Chiral Metamaterials: retrieval of the effective parameters with and without substrate[END_REF][START_REF] Kenanakis | Flexible Chiral Metamaterials in the Terahertz Regime: a Comparative Study of Various Designs[END_REF] that the GOA produces also a growth in the effective values of ε r and µ r . This difference in the values of ξ Chir is due to the higher electromagnetic interactions that arise when the doping impurities (cross-wire structures) are used in GOA materials.

ε r = [1, 3, 12] µ r = 1 ξ Chir = 1 γ = 0.6 NLoS ε r = [1, 3, 12] µ r = 1 ξ Chir = 1 β 1 = π/6, σ = 10
In the case of GOA scenario and indoor LoS propagation (Figure 8), we ob- scenario with GOA effect, we notice low values of path loss for lower frequencies, while increasing trends are observed with higher frequency range, both for the biological and indoor cases (see Figure 9 and Figure 10, respectively). In the typical range of Terahertz nanocommunications i.e., [0.06, 10] THz, NLoS path loss with GOA has low and almost flat values, thus representing an efficient transmission scenario. For higher frequencies, the behavior is strongly dependent on the relative electrical permittivity, and then on the specific features of the transmission material.

To summarize, we can conclude that chirality effects act differently on the LoS path and on the NLoS path. Specifically, low values affect the LoS attenuation with a progressive increase for higher frequency range up to 60 THz, and this effect is more consistent with the growth in the values of ε r . Higher values are also experienced for higher distances (indoor scenario). For a better understanding we obtain

kk • E 0 -ω (ξ -ζ) k × E 0 + ω 2 (µε -ξζ -k • k) E 0 = 0, (43)
with k • k = k 2 ,
where k is the complex amplitude of the propagation vector k.

In order to find the dispersion relation, we need to compute k • E 0 , and then we have to consider the divergence equations for both the electrical and magnetic displacement D 0 and B 0 , respectively. Those equations arise to the following linear homogeneous system of two equations for the unknown quantities k • E 0 , and k • H 0 :

   k • B 0 = ξk • E 0 + µk • H 0 = 0 k • D 0 = εk • E 0 + ζk • H 0 = 0 (44)
We notice that the system in (44) provides (i) ∞

It is evident that only trivial solutions are acceptable. Indeed, Eq. ( 46) is never that is

k = ± µ 0 µ r ε 0 ε r + ξ 2 0 ξ 2 r = ± µ 0 ε 0 (ε r µ r + ξ 2 r ) = ± √ µ 0 ε 0 ε r µ r + ξ 2 r , (52)
being ξ 0 = √ µ 0 ε 0 . It follows that the propagation vector will be:

k = ± √ µ 0 ε 0 ε r µ r + ξ 2 r k. (53)
Finally, the amount ν Chir is the relative refractive index for a reciprocal material (except SNG), whose expression is:

ν Chir = ε r µ r + ξ 2 r . (54)

  et al. deal with a multi-ray channel model that considers not only LoS and NLoS propagation effects, but also reflected, scattered and diffracted paths. This work represents an accurate channel model, where the use of ray tracing technique is exploited. Because of the very short wavelength in the THz band, THz waves propagate quasi-optically and can be modeled highly accurate with a ray optical approach. Ray tracing has already enabled good capabilities at low Terahertz frequency i.e., at 0.3 THz. In [10], Sheikh et al. use a ray-tracing model to simulate the multipath wave propagation in a realistic office environment, and include the frequency dependent atmospheric attenuation.

) Q d PFigure 1 :

 d1 Figure 1: Schematic of the reference geometry, where the transmitter (Tx) and the receiver (Rx) lay along the ŷ-axis with the generic reflecting point Q belonging to {x = 0, z = h}.

  and is depending on the specific location of point Q and the reflecting local plane. This local specular reflection coefficient generates the reflected Poynting's vector P (r) = P y (r) ŷ + P z (r) ẑ = P y NLoS ŷ + P z NLoS ẑ, in the NLoS path, contributing to W NLoS .

  attenuation and the delay, respectively.Through the use of Wiener-Khinchin theorem, the attenuations and delays in the i-th frequency sub-band can be expressed by the module of transfer functions H LoS , and H NLoS , respectively as LoS, and NLoS (i.e., reflected) propagation paths.

in

  LoS and NLoS propagation (assuming reflected paths only) in a chiral affected medium, in the simulation results, we focus only on LoS and NLoS transfer functions with chirality effects, assuming a flat behavior for the molecular absorption loss (i.e., frequency independent behavior). Thus, the following results represent the behavior of transmitted signals propagating along LoS and NLoS directions in a chirality-affected channel, by omitting the frequency-dependent molecular absorption effect.In this framework, we discuss the influence of low values of chirality parameter i.e., ξ Chir = 1, separately from the higher values occurring in the GOA medium i.e., ξ Chir = 10.Simulation results have been obtained under the assumption of homogeneous chirality and with a linear electric polarization. Standard values of the chiral parameter ξ Chir have been chosen varying in the range 0 ≤ |ξ Chir | ≤ 1. However higher values for ξ Chir have been considered[START_REF] Fang | Dual-band Terahertz Chiral Metamaterial with Giant Optical Activity and Negative Refractive Index based on Cross-wire Structure[END_REF][START_REF] Zhou | Terahertz chiral metamaterials with giant and dynamically tunable optical activity[END_REF][START_REF] Zhao | Chiral Metamaterials: retrieval of the effective parameters with and without substrate[END_REF][START_REF] Kenanakis | Flexible Chiral Metamaterials in the Terahertz Regime: a Comparative Study of Various Designs[END_REF] to evaluate the effects of this parameter in the particular media such as the chiral complex materials where a giant optical activity with strong circular dichroism takes place. In this special media, |ξ Chir | can be extended up to |ξ Chir | ≈ 5 ÷ 10.

- 6 GOAε

 6 r = [40, 60, 80] µ r = 40 ξ Chir = 10 (ii) d 1 = d 2 = [1, 10] mm (m), for biological (indoor) LoS and NLoS scenarios,respectively. In GOA scenario, the values of ε r have been chosen in accordance with[START_REF] Fang | Dual-band Terahertz Chiral Metamaterial with Giant Optical Activity and Negative Refractive Index based on Cross-wire Structure[END_REF]-[START_REF] Kenanakis | Flexible Chiral Metamaterials in the Terahertz Regime: a Comparative Study of Various Designs[END_REF]. In the following, we describe the simulation results for LoS, NLoS and GOA scenarios.

Figure 3 Figure 3 :

 33 Figure 3 depicts the total path loss for LoS propagation in a biological environment in the frequency range [0.06, 60] THz, and for different values of electrical permittivity ε r corresponding to the case of (i) vacuum (ε r = 1), (ii) calcium (ε r = 3), and (iii) silicon (ε r = 12). We observe that low values in ξ Chir , i.e. ξ Chir = 1, contribute to a lower rate of attenuation (i.e., < 60 -80 dB) in the lower THz band [0.06, 10] THz, for observation points laying in the range of millimeters (see Figure 3 (a) and (b)). Graphs are plotted for same different valuesin the relative electric permittivity ε r , and can be of some interest for a biological propagation medium. We also observe that path loss increases for higher values

depicted in Figure 5 Figure 4 :

 54 Figure 4: Total path loss for LoS propagation in an indoor environment with a chirality-affected channel (ξ Chir = 1), for (a) d = 1 m, and (b) d = 10 m.

Figure 5 :

 5 Figure 5: Total path loss for NLoS propagation in a biological environment with a chiralityaffected channel (ξ Chir = 1), for (a) d 1 = d 2 = 1 mm, and (b) d 1 = d 2 = 10 mm distance.

Figure 6 :

 6 Figure 6: Total path loss for NLoS propagation in an indoor environment with a chirality-affected channel (ξ Chir = 1), for (a) d 1 = d 2 = 1 m, and (b) d 1 = d 2 = 10 m distance.

Figure 7 :

 7 Figure 7: Total path loss for LoS propagation in a biological environment with GOA effect (ξ Chir = 10), for (a) d = 1 mm, and (b) d = 10 mm. ε r and µ r effective values have been chosen in accordance with [19]-[24].

Figure 8 :

 8 Figure 8: Total path loss for LoS propagation in an indoor environment with GOA effect (ξ Chir = 10), for (a) d = 1 m, and (b) d = 10 m. ε r and µ r effective values have been chosen in accordance with [19]-[24].

Figure 9 :

 9 Figure 9: Total path loss for NLoS propagation in a biologic environment with GOA effect (ξ Chir = 10), for (a) d 1 = d 2 = 1 mm, and (b) d 1 = d 2 = 10 mm. ε r and µ r effective values have been chosen in accordance with [19]-[24].

Figure 10 :

 10 Figure 10: Total path loss for NLoS propagation in an indoor environment with GOA effect (ξ Chir = 10), for (a) d 1 = d 2 = 1 m, and (b) d 1 = d 2 = 10 m. ε r and µ r effective values have been chosen in accordance with [19]-[24].

1

 1 autosolutions for k • E 0 = 0, k • H 0 = 0, if and only if µε -ξζ = 0,(45)and (ii) trivial solutions k • E 0 = 0, k • H 0 = 0, if and only if µε -ξζ = 0.

Table 1

 1

	collects the parameters used in the simulation results in the cases of
	LoS, NLoS, and GOA scenarios. Distance values are (i) d = [1, 10] mm (m), and

Table 1 :

 1 Parameters used in the simulations for LoS, NLoS, and GOA cases. Common parameters are f = [0.06, 60] THz, ε 0 = 8.85 × 10 -12 , µ 0 = 4π10 -7 .

The index i means impressed source distribution, e and m mean electric and magnetic, respectively.

Conclusions

In this paper, we have presented a theoretical approach to evaluate the chirality effects on the electromagnetic propagation in a real complex medium. A full-wave model has been considered and the electromagnetic field has been studied in the THz frequency band with particular emphasis on the determination of the W LoS and W NLoS power terms.

Chirality has been shown to be responsible of low attenuation in LoS propagation, while it increases for NLoS scenario due to the the presence of the reflected rays in two types of environment considered (i.e., biological and indoor applications). Finally, in the case of GOA both in LoS and NLoS, the path loss has a flat behavior for lower frequencies (i.e., < 10 THz).

Appendix

The Maxwell's equations for the harmonic non-conventional medium defined under the hypothesis of plane wave become:

By posing the magnetic field as

verified since

and by assuming a reciprocal material i.e.,

we get ξ r > 0, and ζ r > 0. So, we can deduce that (46) becomes

and this is never null except for Single NeGative (SNG) materials, which are not considered in this paper.

By assuming a reciprocal material, (46) is not verified, and only the trivial solutions are acceptable. It follows that k • E 0 = 0, k • H 0 = 0 and then, based on previous hypothesis, the dispersion relation becomes:

The values of k derive from (50) as: