Sergio Castillo-Perez

Joaquin Garcia Alfaro

Joaquin Garcia-Alfaro
email: joaquin.garcia-alfaro@acm.org

Spyware-based menaces against web applications

Keywords: Network Security, Privacy, Data Protection, Countermeasures

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The use of the web paradigm in all kinds of business models and organizations is becoming truly pervasive. Indeed, its use is becoming an emerging strategy in all kinds of application software companies [START_REF] Cary | A viable solution to enterprise development and systems integration: a case study of web services implementation[END_REF]. It allows the design of fully interactive applications which can potentially be used by thousands of users around the world. The existence of new technologies for the improvement of traditional web features allows software engineers to conceive new services and spaces which are not longer restricted to specific operating systems. Traditional information systems related to education, health care, banking, or even emergency response, can highly benefit with the use of this technology.

The current complexity of the web paradigm has however a direct impact on the security of web browsers and more precisely in the treatment of its resources. Attacks against browsers can compromise the security and privacy of web users. This can have serious consequences given the pervasive presence of malicious software, such as Spyware. Spyware can indeed be secretly installed on Web browsers to steal sensitive data, such as user identities, passwords, and financial data [START_REF] Hu | Is Spyware an Internet Nuisance or Public Menace?[END_REF]. Web browsers must therefore include, in addition to the expected value offered to their users, reliable mechanisms to ensure security and privacy of its users. We survey in this paper some Spyware-based techniques that can be used by malicious entities to violate web users' privacy. We present a sample scenario that shows how the privacy of a user accessing web services, may be violated by Spyware associated to the Web browser. We then discuss some defense mechanisms that might be used to reduce the risk posed by the Spyware threat.

Paper organization -The remainder of this paper has been organized as follows. Section 2 presents the Spyware threat and develops our motivation scenario. Section 3 overviews some defense mechanisms to reduce the risk of the Spyware threat. Section 4 closes the paper.

Spyware

Spyware is a type of malicious software (malware) [START_REF] Skoudis | Malware: Fighting Malicious Code[END_REF] which is installed secretly on the client side to monitor users' behaviors and/or steal sensitive information, such as passwords and credit card numbers. Spyware uses Internet communications to send the stolen information towards malicious databases without the user knowledge or consent. This captured information can eventually be sold with malicious purposes, such as electronic fraud [START_REF] Garcia-Alfaro | Protecting On-line Casinos against Fraudulent Player Drop-out[END_REF], identity theft, spamming, etc. Moreover, Spyware is responsible of system misbehavior, as the appearance of navigation windows with unwanted advertising, web browser hijacking, slow Internet connections speed and backdoor installation to partially control the users' system. Spyware is often complemented by evasion mechanisms, often referred in the literature as rootkits, to ease the process of evading its presence from network and system operators, as well as to avoid being eliminated from the infected system. These mechanisms include the use of third party techniques for hiding processes, files and network connections; and the use of encryption and anti-debugging techniques. The use of these evasion mechanisms, and the corresponding change of basic system functionality features and applications, is the main reason of the system misbehavior aforementioned.

The term Spyware must not be confused with the term virus or worm. Though all three terms refer to malware, they are quite different regarding their operation modes, especially their infection mode. The main difference is that Spyware does not usually self-replicate. In other words, a system affected by Spyware does not spread the infection to other computers. is generally spread onto a system disguised as a different application that the user intends to install. In this case, the infection depends on a direct interaction required by the user. The use of non-certified applications, system updates, or plug-ins, often downloaded from malicious web sites or from malicious peer-to-peer systems, is the most common way of Spyware infection. Alternative ways of Spyware infection also include the exploitation of system application vulnerabilities. For instance, the exploitation of web browsers vulnerabilities during user browsing to force the system to download and install the Spyware. We analyze more in detail the infection process in the sequel.

Infection Process

Traditional Spyware infection mechanisms include the following two strategies: (1) infection processes started by victim users and (2) self-infection processes automatically started through existing system vulnerabilities. In the first case, malicious code contained inside a trusted piece of software gets installed into the system under user's consent. The camouflage of malicious code is often associated with illicit software, such as cracking of commercial applications or altered versions of peer-to-peer (P2P) applications and/or web browser add-ons. The use of social engineering is a perfect example of how manipulating security-unaware user into performing actions, such as opening infected emails or visiting websites associated with the installation of the illicit software. In any case, we consider that in this first category, the infection process requires user's tacit approval.

In the second case, the infection mechanism exploits existing vulnerabilities in system applications. We include in this category the subsequent infection after a user agrees to visit a malicious website or opens a document which exploits either the web browser or the visualization application. The difference now is that the infection does not require user's approval. The existence of a vulnerability here is enough to trigger the infection process by executing the malicious code hidden in the malicious web site or docu-ment. In this case, the installation of the Spyware into the system goes unnoticed by the user. This strategy is often combined with more sophisticated attacks, such as Crosssite scripting attacks (XSS for short) and/or DNS Spoofing redirections. Once the infection process has been successfully performed, the Spyware is ready to start capturing and profiling both user and system data. Different strategies can be then used by the Spyware to implement the subsequent gathering and evasion processes. We describe in the sequel the use of user space vs. kernel space strategies.

User Space Execution -In this case, the malware only has regular user privileges. Today's most common Spyware belongs to this category, since its development is straightforward -even unskilled attackers can do it based on automatic tools. Due its simplicity and the fact that such Spyware acts under user privileges, the detection and disinfection of this category of Spyware are fairly trivial. The strategies employed by this class of malware, in order to intercept information and to evade detection/removal, is mainly based on diverting the execution flow of the affected application, ceding control to a certain code of the malware. Two different techniques might be used: (1) modification of the memory addresses that are mapped to the shared libraries functions used by the affected applications; and (2) exploitation of the add-on/plug-in mechanisms provided by the affected applications. Both techniques would allow the Spyware to get complete control of the application flow, and so to implement the interception and evasion routines.

Kernel Space Execution -The previous technique, although it allows to infect user space memory areas, does not enable the Spyware to further propagate the infection throughout the system. Given that most of the affected applications running in user space are almost always stored in system space areas, the Spyware would need administration privileges to keep spreading the infection everywhere the system. The execution of the Spyware in Kernel space provides with such a functionality. Specifically, the Spyware code is now loaded and executed within the memory space assigned for operating system routines. The execution within this memory space area provides to the malware higher resistance to be detected and/or removed from the system. Indeed, the Spyware code is now running with administration privileges. This implies higher complexity in the Spyware code, which makes it less common. The common way to manage the execution of malware at kernel space relies on the modification of memory addresses associated to the operating system calls (syscalls for short).

Web Browser Infection Scenario

We show in this subsection a hypothetical scenario where the privacy of a given user, who is accessing web services, may be violated by Spyware associated to the browser. The technique used by the Spyware of our scenario to intercept user's data is to exploit the add-on mechanisms of certain programs -in our case, the web browser. Most of today's web browsers allow this option in order to ease the adding of new enhancements by third party programmers. To do so, the main application's routine delegates the execution flow to the new functions and grants them the permission to control certain events. For instance, the main routine may grant the access to both read and modify the Document Object Model (DOM). If so, these complements may control the web browser application during the loading of a web page, the processing of uploaded/downloaded files, etc. Hence, the way that this pluggable technology works can be exploited to put in place the necessary Spyware routines.

Let us suppose that a malicious user wants to obtain valid credit card numbers with the aim of carrying out electronic fraud. Let us also assume that the malicious user finds out a given vulnerability associated to a particular version of a web browser that allows the execution of arbitrary code. Analyzing the vulnerability, our malicious user would prepare various web servers containing the necessary code that exploits the vulnerability and would spread the necessary links (e.g., by using social engineering) [START_REF] Garcia-Alfaro | Protecting On-line Casinos against Fraudulent Player Drop-out[END_REF]. The visit of the victims to those malicious web sites leads to the Spyware installation by exploiting the vulnerability. The objective of the Spyware programmed by the malicious user is to capture and redirect towards a secret place those credit card numbers introduced by the victims during their electronic web-based sessions. Hence, the aforementioned pluggable techniques will be used by the installed Spyware.

To guarantee the success of the previous scenario, the attacker must persuade a large number of victims (i.e., users of the vulnerable web browser) to visit the malicious web sites. The more victims the attacker would successfully fool, the higher the attacker's profit. The time factor would also be decisive in this process since, as soon as the browser vulnerability is reported and fixed, the likelihood of success for the attack will drastically be reduced. With the objective of speeding up the infection process, several other strategies such as DNS poisoning, Spamming, and cross-site scripting attacks, are going to be used by the attacker to lead as much victims as possible towards the malicious web sites.

After the infection, the Spyware -now seen as a browser add-on -remains active, but silent, while waiting for the matching of data associated to e-commerce transactions (i.e., to capture data such as credit card numbers, expiration numbers, passwords, etc.). All the information captured by the Spyware is eventually sent, once postprocessed and protected, to the attacker's database. Let us notice that although the web browser can use cryptographic operations to ensure the authenticity, integrity and confidentiality of the information exchanged with e-commerce servers (e.g., use of SSL/TLS communications), this does not protect the data from being stolen by the Spyware which is locally executed at the client side.

If we analyze the use of SSL/TLS in the TCP/IP model, we can see that the information remains protected between the application layer (HTTPS) and the transport layer (TCP). However, since the Spyware that infected the browser is executed at the application layer (i.e., HTTPS), it intercepts the information before being protected by SSL/TLS. Therefore, it can store and modify such information before it is sent to the following layers of the TCP/IP model. We should also notice that the Spyware does not raise either any suspicious alert associated to the use of SSL/TLS certificates, since the verification process is not affected by the infection. Similarly, the use of verification codes, such as CVC2, CVV or CID, does not help either to detect or prevent the malicious activities associated to the Spyware. These codes are also provided by the users and certainly captured by the Spyware without any difficulties.

Prevention Techniques

We can consider three main forms of prevention to avoid the infection of a system by Spyware. The first one, more ambitious, is to encourage users to improve their computer hygiene. This includes the enhancement of their behavioral patterns, such as keeping their operating systems and associated components updated, to avoid downloading software from untrusted sources, and to ignore emails or attachments coming from unknown senders.

The second form of prevention, much more technical, can be defined as a protection design pattern. The objective is to avoid the infection of a system or to reduce the damage of the infection. This includes hierarchical protection domains, often called protection rings. The result is a layered structure of trust in the operating system. This structure of trust is generally complemented by specific hardware that allows a protected separation between trusted and untrusted processes. It provides, moreover, several levels of access to the resources of the system. Rings, or domains, are arranged in a hierarchical manner: from most privileged (most trusted) to least privileged (least trusted). This protection mechanism is specially relevant to reduce the risk of Spyware targeting the kernel level of an operating system. It increments the technical difficulties to solve to get control of the system services. Attacks would target now some kind of low-level hardware vulnerabilities [START_REF] Embleton | SMM Rootkits: a New Breed of OS Independent Malware[END_REF].

The third defense consists on using automatic applications to increase the odds of detecting and isolating incoming Spyware. Depending on the detection strategy, this automatic software can be classified in two main categories: syntactic signatures based analysis and semantic-aware analysis. Syntactic signatures based analysis intends to collect as much information about suspicious applications as possible. It avoids moreover the execution of suspicious files. Then, using this characterization, it determines whether the file can be considered as Spyware. The main detection mechanism basically compares the suspicious file towards a database of characterizing signatures, generally seen as a sequence of low-level instructions. This pattern-based signature database is generated by analyzing previously known Spyware samples and, as a consequence, it must be regularly updated by the users. As Moser et al. points in [START_REF] Moser | Limits of Static Analysis for Malware Detection[END_REF], there are some evasion schemes based on some kind of obfuscation, such as polymorphism or metamorphism, that can be used by the Spyware authors to remain undetectable. To do so, some kind of code mutation alter the files without changing its behavior. So, every infection can be performed by a different image, defeating pattern-based signatures detection.

Semantic-aware analysis [START_REF] Egele | Dynamic Spyware Analysis[END_REF], on the other hand, can be implemented by modeling the interaction of the suspicious files and the system environment. This strategy can be seen as a way to overcome the deficiencies of the previous strategy, since any mutation of the executable image does not affect the final behaviour of Spyware. This can be understood if we consider that syntactic signatures does not take into account the semantic of the instructions. Thus, this detection mechanism is also sometimes known as behavior-based detection. According to [START_REF] Jacob | Behavioral Detection of Malware: From a Survey Towards an Established Taxonomy[END_REF], we can classify these strategies in two categories depending how data collection to model the behaviour is extracted, that is static or dynamic. In the last years, several behavioural detection strategies has been developed based on different algorithms, such as graph isomorphism, model checking or heuristics. We encourage the reader to consult [START_REF] Jacob | Behavioral Detection of Malware: From a Survey Towards an Established Taxonomy[END_REF] for a taxonomy details. The combination of static and dynamic strategies can be merged to enhance the detection process [START_REF] Kirda | Behavior-Based Spyware Detection[END_REF]. However, these techniques are not completely effective for the detection of today's and future generations of Spyware. First, because it is quixotic to believe that we can build sets of signatures that are complete enough to detect all kinds of Spyware [START_REF] Lee | Investigating Factors Affecting the Adoption of Anti-spyware Systems[END_REF]; and second, because the strategy of determining whether a program exhibits a specific malicious behavior independent of a suspicious file is a undecidable problem [START_REF] Christodorescu | Semantics-Aware Malware Detection[END_REF].

Conclusions

We introduced in this paper the Spyware menace for compromising the security and privacy of web browser resources. This menace can especially harm web-based applications, such as online education, health care, and electronic banking. We surveyed some Spyware-based techniques that can be used by malicious entities to infect a system; and briefly introduced some prevention techniques. We plan to analyze more in deep these techniques in a forthcoming version of this paper.

Acknowledgments -We acknowledge the financial support received from the Spanish Ministry of Science and Innovation and the FEDER funds (grants TSI2006-03481, TSI2007-65406-C03-03 E-AEGIS, and CONSOLIDER-INGENIO 2010 CSD2007-00004 ARES).


Alexandre Graell
email: alexandre.graell@telecom-bretagne.eu

Eirik Rosnes

Stopping Set Analysis of Repeat Multiple-Accumulate Codes

In this work, we consider a stopping set analysis of repeat multiple-accumulate (RMA) code ensembles formed by the serial concatenation of a repetition code with multiple accumulators. The RMA codes are assumed to be iteratively decoded in a constituent code oriented fashion using maximum a posteriori erasure correction in the constituent codes. We give stopping set enumerators for RMA code ensembles and show that their stopping distance hmin, defined as the size of the smallest nonempty stopping set, asymptotically grows linearly with the block length. Thus, the RMA code ensembles are good for the binary erasure channel. Furthermore, it is shown that, contrary to the asymptotic minimum distance dmin, whose growth rate coefficient increases with the number of accumulate codes, the hmin growth rate coefficient diminishes with the number of accumulators. We also consider random puncturing and show that for sufficiently high code rates, the asymptotic hmin does not grow linearly with the block length, contrary to the asymptotic dmin, whose growth rate coefficient approaches the Gilbert-Varshamov bound as the rate increases. Finally, we give iterative decoding thresholds to show the convergence properties.

I. INTRODUCTION

The performance of a concatenated code in the error floor region for an additive white Gaussian noise channel is dominated by the distance spectrum of the code and, in particular, by its minimum distance d min . Turbo codes [START_REF] Berrou | Near Shannon limit error-correcting coding and decoding: Turbo-codes[END_REF] and single serially concatenated codes [START_REF] Benedetto | Serial concatenation of interleaved codes: Performance analysis, design, and iterative decoding[END_REF] are asymptotically bad codes, in the sense that their minimum distance does not asymptotically grow linearly with the block length. Recently, it was shown that the minimum distance d min of very simple serially concatenated code ensembles built from the concatenation of a repeat code with two or more accumulators increases as the number of accumulators increase. In fact, in the limit of infinitely many accumulators, the d min approaches the Gilbert-Varshamov bound (GVB) [START_REF] Pfister | The serial concatenation of rate-1 codes through uniform random interleavers[END_REF]. Furthermore, in [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF][START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF], it was shown that, indeed, repeat multiple-accumulate (RMA) code ensembles are asymptotically good, in the sense that their d min asymptotically grows linearly with the block length. A method for the calculation of a lower bound on the growth rate coefficient was also given in [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF].

For the binary erasure channel (BEC) with iterative decoding, a similar role to that of the distance spectrum with A. Graell i Amat was funded by a Marie Curie Intra-European Fellowship within the 6th European Community Framework Programme. E. Rosnes was funded by the Norwegian Research Council under Grants 174982 and 183316.

C 0 C 1 C L K w qw N N N h h 1 π 1 π L Fig. 1.
Encoder structure for RMA codes.

maximum-likelihood decoding is played by the stopping set distribution. Stopping sets for iterative belief-propagation (BP) decoding of low-density parity-check (LDPC) codes were introduced in [START_REF] Di | Finite-length analysis of low-density parity-check codes on the binary erasure channel[END_REF]. The concept of stopping sets was later extended in [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF] to turbo decoding by introducing turbo stopping sets. In this work, we extend the asymptotic minimum distance analysis in [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF][START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF] by considering stopping sets for iterative constituent code oriented decoding using maximum a posteriori (MAP) erasure correction in the constituent codes of RMA codes. We derive a closed-form expression for the subcode input-output support size enumerating function (SIOSEF) [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF] of an accumulate code and obtain stopping set enumerators for RMA code ensembles. We then analyze the asymptotic behavior of the stopping set distribution and show that RMA code ensembles exhibit an asymptotic stopping distance h min , defined as the size of the smallest nonempty stopping set, that grows linearly with the block length. Therefore, RMA codes are good for the BEC. However, contrary to the asymptotic d min , whose growth rate coefficient increases with the number of accumulate codes, the asymptotic h min growth rate coefficient diminishes with the number of accumulators. We also consider random puncturing of RMA codes and show that for sufficiently high code rates, the asymptotic h min does not grow linearly with the block length, contrary to the asymptotic d min , whose growth rate coefficient approaches the GVB as the rate increases [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF]. Finally, we consider extrinsic information transfer (EXIT) charts [START_REF] Brink | Convergence behaviour of iteratively decoded parallel concatenated codes[END_REF] to analyze the convergence properties of RMA codes. We will show that repeat accumulate-accumulate (RAA) codes are good codes for the BEC, while increasing the number of accumulators decreases both the stopping distance and the convergence threshold.

II. ITERATIVE CONSTITUENT CODE ORIENTED DECODING

AND STOPPING SETS The encoder structure of RMA codes is depicted in Fig. 1. It is a serial concatenation of a (Kq, K) repetition code C 0 with input block length K and rate 1/q with L ≥ 1 identical rate-1, memory-one, accumulators C l , l = 1, . . . , L, with generator polynomials g(D) = 1/(1+D), through interleavers π 1 , . . . , π L . The overall nominal code rate (avoiding termination) is denoted by R = K/N = 1/q, where N denotes the output block length. Higher rates may be obtained by puncturing the output of the most inner accumulator C L .

In this paper, we consider stopping sets for iterative constituent code oriented decoding using MAP constituent decoders (turbo decoding) over the BEC of RMA codes. With iterative constituent code oriented decoding we mean a decoding strategy that iterates between the constituent decoders. MAP decoding of the constituent codes can be implemented efficiently on a trellis representation of the constituent code [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF]. Thus, the complexity of iterative constituent code oriented decoding of RMA codes is linear in the block length N when the number of constituent code activations is independent of N .

A. Stopping Sets for RMA Codes

We will now give the formal definition of a stopping set for RMA codes, adapted from the definition in [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF] for turbo stopping sets. The generalization to the case with puncturing is straightforward. In the definition, we need the concept of support set of a binary vector and of a binary linear code. The support set χ(x) of a binary vector x = (x 1 , . . . , x N ) (of length N ) is the set of nonzero coordinates. As an example, with x = (0, 1, 1, 0, 1), χ(x) = {2, 3, 5}. The support set χ(C) of a binary linear code C is the union of the support sets of each codeword in C. Finally, an interleaver will be regarded as a mapping from the set of coordinates of its input sequence to the set of coordinates of its output sequence. Definition 1. Let C RMA denote a given RMA code with interleavers π 1 , . . . , π L . A set S = S(π 1 , . . . , π L ) ⊆ {1, . . . , N } of the coordinates of the output sequence (or codeword) is a stopping set if and only if there exist L + 1 linear subcodes Ĉl ⊆ C l ⊆ {0, 1} N , l = 0, . . . , L, with support sets χ( Ĉl ) such that

χ( ĈL ) = S and π l (χ( Ĉl-1 )) = χ( Ĉl ), l = 1, . . . , L.
The size of a stopping set S is its cardinality.

In Definition 1 we used the fact that the mapping between the input support set and the output support set is an identity mapping for rate R = 1 encoders. Note also that Definition 1 does not exclude the empty set. Thus, the empty set is formally a stopping set of size zero. We stress the fact that the concept of stopping sets for RMA codes, as defined above, is conceptually different from the traditional concept of stopping sets used in connection with iterative BP decoding of LDPC codes, but it reduces to the traditional concept of stopping sets for Tanner graphs when the constituent codes are single parity-check codes. We can prove the following theorem. Theorem 1. Let C RMA denote a RMA code that we use to transmit information over the BEC. The received vectors are decoded iteratively using constituent code oriented decoding using MAP erasure correction in the constituent codes until either the codeword has been recovered, or the decoder fails to progress further. Then the set of erased positions that remain when the decoder stops is equal to the unique maximum-size stopping set that is contained in the (initial) set of erased positions.

From Theorem 1 it follows that an important parameter for code performance is the stopping distance h min .

III. SUPPORT SETS AND STOPPING SET ENUMERATORS

Let C denote an (N, K) binary linear code. Partition all the subcodes of C of dimension d, d = 0, . . . , K, into equivalence classes based on their support sets. In particular, all subcodes within a specific subcode class are required to have the same support set, but the subcodes may have different dimensions. We define the SIOSEF [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF] of C as

A C (W, H) = K w=0 N h=0 a C w,h W w H h
where W and H are dummy variables, and a C w,h is the number of subcode classes of C of input support set size w and output support set size h. In the rest of the paper, with a slight abuse of notation, we will refer interchangeably to both A C (W, H) and a C w,h as the SIOSEF of a code C.

A. Stopping Set Enumerators for RMA Code Ensembles

Benedetto et al. introduced in [START_REF] Benedetto | Serial concatenation of interleaved codes: Performance analysis, design, and iterative decoding[END_REF] the concept of uniform interleaver to obtain average weight enumerators for concatenated code ensembles from the weight enumerators of the constituent encoders. The same concept can be used to obtain stopping set enumerators of concatenated code ensembles by combining the SIOSEFs of the constituent encoders. Let sC w,h be the ensemble-average input-output stopping set size enumerating function (IOSSEF) of the code ensemble C with input and output block length K and N , respectively, denoting the average number of stopping sets of input size w and output size h over C. Also, denote by sC h = K w=0 sC w,h the ensembleaverage stopping set size enumerating function (SSEF) of the code ensemble C, giving the average number of stopping sets of size h over C. Using the concept of uniform interleaver [START_REF] Benedetto | Serial concatenation of interleaved codes: Performance analysis, design, and iterative decoding[END_REF], the IOSSEF of a RMA code ensemble C RMA can be written as [START_REF] Benedetto | Serial concatenation of interleaved codes: Performance analysis, design, and iterative decoding[END_REF][START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF]

sCRMA w,h = N h1=0 • • • N hL-1=0 a C0 w,qw a C1 qw,h1 N qw L-1 l=2 a C l h l-1 ,h l N h l-1 a CL hL-1,h N hL-1 = N h1=0 • • • N hL-1=0 sCRMA w,h1,...,hL-1,h (1)
where sCRMA w,h1,...,hL-1,h is called the conditional support size enumerating function of C RMA .

The evaluation of (1) requires the computation of the SIOSEFs of the constituent encoders, which is addressed below.

K w ⌊ qw 2 ⌋ d1=1 N -h1 d1 h1-1 d1-1 h1-d1 qw-2d1 N qw L-1 l=2 j h l-1 2 k d l =1 N -h l d l h l -1 d l -1 h l -d l h l-1 -2d l N h l-1 • j h L-1 2 k dL=1 N -h dL h-1 dL-1 h-dL hL-1-2dL N hL-1

B. SIOSEFs for Memory-One Encoders and the Repetition Code

Theorem 2. The SIOSEF for rate-1, memory-one, convolutional encoders with generator polynomials g(D) = 1/(1+D) and g(D) = 1 + D that are terminated to the zero state at the end of the trellis and with input and output block length N can be given in closed form as

a 1 1+D w,h = a 1+D h,w = ⌊ w 2 ⌋ d=1 N -h d h -1 d -1 h -d w -2d (2)
for positive input sizes w. Also, a 1 1+D 0,0 = a 1+D 0,0 = 1. A detailed proof can be found in [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF]. Theorem 3. The SIOSEF for the (Kq, K) repetition code C 0 with input block length K can be given in closed form as

a C0 w,qw = K w . (3)
Using ( 2) and ( 3) in (1), we get the expression at the top of the page for the conditional support size enumerating function (with w > 0) of RMA code ensembles.

IV. FINITE-LENGTH ANALYSIS OF STOPPING DISTANCE

The ensemble-average SSEF can be used to lower bound h min for finite block lengths, similarly to the case of the d min [START_REF] Pfister | The serial concatenation of rate-1 codes through uniform random interleavers[END_REF][START_REF] Gallager | Low-Density Parity-Check Codes[END_REF]. In particular, the following bound holds.

Lemma 1. The probability that a code chosen randomly from the ensemble C with ensemble-average SSEF sC

h has stopping distance h min < is upper-bounded as

Pr(h min < ) ≤ -1 h=1 sC h . (4)
Lemma 1 can be used to obtain a probabilistic lower bound on the stopping distance of a code ensemble. In particular, if we set Pr(h min < ) = ε, where ε is any positive value between 0 and 1, we would expect that at least a fraction 1 -ε of the codes in the ensemble have a stopping distance h min of at least . In Fig. 2, we plot the probabilistic lower bound on h min from Lemma 1 for RMA codes with repeat factor q = 4, and L = 2, 3, and 4, as a function of the code length N . The bounds were obtained by setting ε = 0.5 in (4), i.e., at least half of the codes in the ensemble have stopping distance at least equal to the value indicated by the curves. All codes appear to have linear h min growth rate. The best growth rate is obtained for L = 2, while increasing the number of accumulate codes decreases the growth rate. For instance, the bound on h min for the RAA code ensemble indicates a stopping distance of at least 132 for block length N = 1000 bits. The value is reduced to 122 and to 80 for the repeat triple-accumulate (RAAA) and the repeat quadrupleaccumulate (RAAAA) code ensembles, respectively.

In the following, we obtain the asymptotic expression for the stopping set enumerator of RMA code ensembles, and we show that, indeed, their h min asymptotically grows linearly with the block length.

V. ASYMPTOTIC ANALYSIS OF STOPPING DISTANCE

We define the asymptotic stopping set size spectral shape function as [START_REF] Gallager | Low-Density Parity-Check Codes[END_REF] r s (ρ) = lim

N -→∞ sup 1 N ln sC ⌊ρN ⌋
where sup(•) denotes the supremum of its argument, ρ = h N is the normalized stopping set size, and N is the block length. If there exists some abscissa ρ 0 > 0 such that sup ρ≤ρ * r s (ρ) < 0 ∀ρ * < ρ 0 , and r s (ρ) > 0 for some ρ > ρ 0 , then it can be shown (using Lemma 1 for example) that, with high probability, the stopping distance of most codes in the ensemble grows linearly with the block length N , with growth rate coefficient of at least ρ 0 . On the other hand, if r s (ρ) is strictly zero in the range (0, ρ 0 ), it cannot be proved directly whether h min grows linearly with the block length or not. In [START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF], it was shown that the spectral shape function of RMA codes, for the codeword case, exhibits this behavior, i.e., it is zero in the range (0, ρ 0 ) and positive for some ρ > ρ 0 , where ρ means here the normalized output weight. By combining the asymptotic spectral shapes with the use of bounding techniques, it was proved in [START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF]Theorem 11] that the minimum distance of RMA codes indeed grows linearly with the block length with growth rate coefficient of at least ρ 0 . We remark that in the rest of the paper, with a slight abuse of language, we sometimes refer to ρ 0 as the exact value of the asymptotic growth rate coefficient. However, strictly speaking, ρ 0 is only a lower bound on it. Now, by using Stirling's approximation for the binomial coefficient n k n→∞ -→ e nH(k/n) where H(•) is the binary entropy function with natural logarithms, and the fact that w, h 1 , . . . , h L-1 , d 1 , . . . , d L , and h can all be assumed to be of the same order as N [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF], sCRMA w,h1,...,hL-1,h can be written as sCRMA w,h1,...,hL-1,h = d1,...,dL

exp {f (α, β 1 , . . . , β L-1 , γ 1 , . . . , γ L , ρ) N + o(N )}
when N -→ ∞, where α = w K is the normalized input stopping set size, ρ = h N is the normalized output stopping set size, β l = h l N is the normalized output support set size of code C l , γ l = d l N , and the function f (•) is given by

f (β 0 , β 1 , . . . , β L-1 , γ 1 , . . . , γ L , ρ) = H(β 0 ) q - L l=1 H (β l-1 ) + L l=1 (1 -β l )H γ l 1 -β l + L l=1 β l H γ l β l + L l=1 (β l -γ l )H β l-1 -2γ l β l -γ l ( 5
)
where for conciseness we defined β 0 = α and β L = ρ. Finally, the stopping set size spectral shape function for RMA code ensembles can be written as [4]

r CRMA s (ρ) = sup 0<β0,...,βL-1≤1 0<γ1,...,γL≤1 f (β 0 , β 1 , . . . , β L-1 , γ 1 , . . . , γ L , ρ) (6)
To analyze the asymptotic stopping distance behavior of RMA code ensembles, we must solve the optimization problem in (5)- [START_REF] Di | Finite-length analysis of low-density parity-check codes on the binary erasure channel[END_REF]. Details can be found in [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF]. The numerical evaluation of ( 5)-( 6) is shown in Fig. 3 for RAA code ensembles with q = 2, 3, 4, 5, and 6. We observe that the stopping set size spectral shape function for the rate R = 1/2 RAA code ensemble is strictly positive, meaning that the ensemble is bad for the BEC. For 3 ≤ q ≤ 6, the function r CRMA s (ρ) is zero in the range (0, ρ 0 ) and positive for some ρ > ρ 0 . In this case, we cannot conclude directly whether h min grows linearly with the block length or not. Define the function

ψ(u, ρ) = sup max(0,u-ρ)≤γ≤min(ρ,1-ρ,u/2) -H(u) + ρH γ ρ + (1 -ρ)H γ 1 -ρ + (ρ -γ)H u -2γ ρ -γ .
We can prove the following theorem (see [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF] for details), extending the results in [START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF] to the stopping distance case.

Theorem 4. Define ρ 0 = max{ρ * ∈ [0, 1/2) : r CRMA s (ρ) = 0 ∀ρ ≤ ρ * }. Assuming that lim u-→0 ψ(u,ρ) u < 0 ∀ρ < ρ 0 , then ∀ρ * > 0 lim N -→∞ Pr (h min ≤ (ρ 0 -ρ * )N ) = 0
for L ≥ 3 and q ≥ 2, and L = 2 and q ≥ 3. Thus, if ρ 0 > 0 and r CRMA Fig. 3. Asymptotic stopping set size spectral shape function for the RAA code ensemble with q = 2, 3, 4, 5, and 6. We remark that it can be verified that the assumption in Theorem 4 always holds for the numerical values of ρ 0 that we have found. The exact values of ρ 0 are given in Table I (except for q = 6). For comparison purposes we also give the asymptotic growth rate coefficient of the d min computed in [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF]. As expected, the asymptotic growth rate coefficient of h min is smaller than for d min . We can now prove the following theorem (see [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF] for details).

TABLE I ASYMPTOTIC h min GROWTH RATE COEFFICIENTS FOR RMA CODE ENSEMBLES. q = 2 q = 3 q = 4 q = 5 ρ0 (hmin) (L = 2) N/A 0.
Theorem 5. The typical h min of RMA code ensembles for L ≥ 3 and q ≥ 2 grows linearly with block length.

In Table I, we also report the asymptotic growth rate coefficient ρ 0 for RAAA and RAAAA code ensembles. Interestingly, contrary to the asymptotic d min growth rate coefficient, which increases with the number of accumulators and tends to approach the GVB [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF][START_REF] Fagnani | Spectra and minimum distances of repeat multiple accumulate codes[END_REF], the asymptotic growth rate coefficient of h min decreases with the number of accumulators concatenated in series. The results are in agreement with the finite-length analysis in the previous section.

On the other hand, the h min of RA code ensembles does not grow linearly with the block length (see [START_REF] Graell I Amat | Good concatenated code ensembles for the binary erasure channel[END_REF] for details).

A. RMA Codes with Random Puncturing

In this section, we consider high rate RMA code ensembles obtained by randomly puncturing the output of the most inner accumulator C L of a RMA code ensemble. Denote by λ (0 ≤ λ ≤ 1) the puncturing permeability rate, i.e., the fraction of bits surviving after puncturing, and by R ′ = R/λ = 1/(λq)

a C punct. w,h ′ = N h=h ′ a C w,h h h ′ N -h λN -h ′ N λN . (7)
Using Stirling's approximation in [START_REF] Rosnes | Turbo decoding on the binary erasure channel: Finite-length analysis and turbo stopping sets[END_REF] and coupling it with ( 5) and ( 6), the stopping set size spectral shape function of a punctured RMA code ensemble r

CRMA-punct. s (ρ ′ ), where ρ ′ = h ′ λN , can be derived.
The values of ρ 0 corresponding to r CRMA-punct. s (ρ ′ ) are given in Table II for L = 2, 3, and 4 mother RMA code ensembles with q = 4 for several nominal code rates R ′ . Asymptotic h min linear growth can be guaranteed for some rates R ′ > 1/q. However, it is interesting to note that the asymptotic linear growth rate property breaks down for heavy puncturing of the mother code ensemble. This behavior gets worse as the number of encoding stages increases. These results are also supported by the finite-length analysis. In Fig. 2, we plot the probabilistic lower bound on h min from Lemma 1 for punctured q = 4 RMA code ensembles with λ = 3/4 and L = 2 and 3. The results are in agreement with the asymptotic analysis; for L = 2 and permeability rate λ = 3/4 linear growth rate is guaranteed. However, applying the same puncturing to the L = 3 RMA code ensemble breaks down the linear growth rate property. Note that these results are in contrast with the results in [START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF], where it was observed that the asymptotic normalized d min gets closer to the GVB for higher rates with random puncturing.

VI. EXIT CHARTS ANALYSIS

In this section, we address iterative constituent code oriented decoding of RMA code ensembles on the BEC by using EXIT charts analysis [START_REF] Brink | Convergence behaviour of iteratively decoded parallel concatenated codes[END_REF] to estimate the convergence thresholds. Note that for the BEC the EXIT functions of the repeat and the accumulate codes can be given in closed form [START_REF] Shi | Exact EXIT functions for convolutional codes over the binary erasure channel[END_REF]. The convergence thresholds for RMA code ensembles for q = 3, . . . , 6 are given in Table III. Note that RAAA and RAAAA code ensembles show very poor thresholds, which make them impractical. From the EXIT charts analysis and the asymptotic analysis in Section V, it arises that double serially concatenated code ensembles are good ensembles for the BEC, since they provide both high asymptotic h min growth rates and good convergence behavior, while adding more encoding stages penalizes both the asymptotic h min growth rate and the convergence threshold.

VII. CONCLUSION

In this paper, we considered RMA code ensembles on the BEC. By deriving a closed-form expression for the SIOSEF of the accumulate code, we obtained stopping set enumerators for RMA code ensembles. We then analyzed the asymptotic behavior of the stopping distance and showed that RMA code ensembles exhibit a stopping distance that grows linearly with the block length. Contrary to the minimum distance, whose growth rate coefficient increases with the number of accumulators, the growth rate of h min diminishes with the number of accumulators. We also considered puncturing of the RMA code ensembles and showed that linear growth can be obtained for rates larger than 1/q. However, the linear growth property breaks down for heavy puncturing.

Fig. 2 .

 2 Fig.2. Probabilistic lower bound on the stopping distance h min versus block length N for RMA codes with q = 4 and L = 2, 3, and 4.

s

  (ρ) ≥ 0 ∀ρ, then almost all codes in the ensemble

TABLE II ASYMPTOTIC

 II h min GROWTH RATE COEFFICIENTS FOR PUNCTURED RMA CODE ENSEMBLES WITH q = 4 AND NOMINAL CODE RATE R ′ = 1/(λq).

	R ′	1/4	0.28	0.29	3/10	5/16	0.33	1/3	11/30	2/5	0.41	0.42	0.43
	λ	1	25/28	25/29	5/6	4/5	25/33	3/4	15/22	5/8	25/41	25/42	25/43
	ρ0 (L = 2)	0.1289	0.1192	0.1142	0.1077	0.0977	0.0819	0.0788	0.0474	0.0188	0.0112	0.0045	N/A
	ρ0 (L = 3)	0.1194	0.0694	0.0528	0.0373	0.0198	0.0004	N/A	N/A	N/A	N/A	N/A	N/A
	ρ0 (L = 4)	0.0784	0.0112	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A

TABLE III CONVERGENCE

 III THRESHOLDS FOR RMA CODE ENSEMBLES. rate of the punctured RMA code. The SIOSEF of a randomly punctured code C punct. with input support set size w, output support set size h before puncturing, and output support set size h ′ after puncturing is given by[START_REF] Kliewer | New results on the minimum distance of repeat multiple accumulate codes[END_REF]

		q = 3	q = 4	q = 5	q = 6
	RAA	0.4965	0.5422	0.5719	0.5935
	RAAA	0.3259	0.3531	0.3718	0.3860
	RAAAA	0.1957	0.2105	0.2209	0.2290
	the nominal code


Sandrick Le Maguer

Nestor Pena

Michel Ney

Matched absorbing medium techniques for full-wave TLM simulations of microwave and millimeter-wave components

Keywords: Circuit hyperfrequence, Electromagnetisme, Milieu absorbant, Modelisation, Methode numerique, Ligne transmission, Methode _ matricielle, Ligne microruban, Ligne coplanaire : Adapta llon impedance

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

This paper investigates the absorbing layer tech niques for three-dimensional transmission-line matrix (3D-TIM) simulations. Applications are the full-wave analysis of structures such as microwave and millime ter-wave circuits open to free space. Two approaches are presented, namely the matched layer (ML) and the perfectly matched layer (PML) technique which is used for general cases. For the ML technique, which is effi cient for matching guides operating in a quasi-TEM condition, a comparison of the dif f erent TIM condensed node is carried out. Concerning the more general case of the PML technique, the theoretical development of a new general three-dimensional (3D) TIM condensed cell is presented. Results pertaining to the characterization of planar circuits are compared with other solutions or measurements. It is found that a proposed unified approach, which includes the TIM simulation for mat ched layer media, yields very good performance in terms of the absorbing boundary condition (ABC). Cet article presente des techniques des milieux absor bants pour la modelisation par la mithode de la matrice des lignes de transmission tridimensionnelle ( 3D-TIM). Les applications sont !'analyse e!ectronu1gnitique de cir cuits microondes et millimetriques ouverts sur l'espace libre. Les deux approches discutees utilisent !es milieux adaptes (ML) et !es couches parfaitement adaptees (PML), ces demieres etant utilisees pour le cas general. La tech nique ML, qui est ejficace pour !'adaptation de guides operant en modes quasi statiques, est itudiee dans le cas de dijferents types de cellules TIM condensees. Pour la simulation des milieux PML applicables aux cas plus generaux, le developpement theorique d'une nouvelle cellule TLM tridimensionnelle est presente. Dans le cadre de caracterisation de discontinuites en technolo gie plaquee, des resultats de simulations sont compares avec des solutions trouvees par d' autres mithodes ou la me sure. On trouve que !' approche unifiee pro po see donne de tres bonnes performances en termes de condi tions absorbantes aux limites (ABC).

I. INTRODUCTION

The lLM method has been widely used for field com putations in the time domain. It is a general method that can be applied to a wide class of problems, including non linear and/or non-stationary cases [l]- [START_REF] Hoefer (y{ | The transmission line matrix method theory and applications[END_REF]. However, like any method whose formulation leads to volume dis cretization, absorbing boundary conditions (ABC) must be used to limit the computational domain when open structures are investigated. MRese structures can be, for instance, open to free space (antenna) or of infinite extend (e.g., waveguide discontinuity problems). Gene rally, ABC should have absorb waves over a wide fre-quency range and at any incident angle. They are usually classified into two groups : those based on boundary operators and those using absorbing layer techniques. For 1LM, the most popular boundary operators are based on one-way equation solutions [START_REF] Eswarappa (ch | Absorbing boundary condi tions for time-domain TLM and FDTD analysis of electromagnetic structures[END_REF][START_REF] Morente | Absorbing boun dary conditions for the TLM method[END_REF] or extrapolation techniques [START_REF] Chen | Absorbing and connecting boundary conditions for the TLM method[END_REF]. MReir algorithm is applied to the limiting boundary where field values at future time are computed by linear combination of field values in the computatio nal domain at different locations and previous time steps. Algorithms are rather simple to apply but these techniques yield numerical instabilities under certain conditions and regularization procedures must be applied to obtain stable results [START_REF] Eswarappa (ch | Absorbing boundary condi tions for time-domain TLM and FDTD analysis of electromagnetic structures[END_REF]. In addition, such ABC should be located at relatively large distances from scat terers as they cannot absorb evanescent waves.

Concerning absorbing layer techniques, equations governing fields in the lossy layer, terminated by certain arbitrary boundary conditions, must be solved by a nume rical scheme. Therefore, the computational domain must be extended, unlike the boundary operator approach. For instance, the so-called matched layer (ML) technique consists in surrounding the computational domain with a layer of fictitious lossy medium with electric and magne tic conductivities which are chosen accordingly. This technique yields excellent results for any frequency when normal incidence is approximately achieved [START_REF] Le Maguer | Etude des couches adap tees pour le calcul des parametres S de structures planaires en TLM 3D. NUMELE C[END_REF]. This condition arises, for instance, for matched load simula tions of guides operating under quasi-TEM conditions. In addition, the ML technique has the great advantage of being inherently stable and the layer medium can be directly modeled by the standard symmetrical condensed node (scN) without extra programming.

The extension of the ML technique to operate for arbitrary incidence yields to the perfectly matched layer (PML) technique, which was first proposed by Berenger for finite-diM¤erence time-domain (FDm) computations [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]- [START_REF] Berenger | A perfectly matched layer for fre e-space simu lation in finite-difference computer codes[END_REF]. The layer surrounding the computational domain is perfectly matched in the sense that electroma gnetic waves, whether traveling or evanescent, theoreti cally penetrate into it with no reflection and refraction in all cases of incidence and at arbitrary frequency. The implementation of the PML technique in the FDTD scheme is straightforward as the layer equations are directly written in a difference form using the standard Yee's algorithm [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Concerning the PML-M}M case, two different approaches can be considered:

•interfacing the RM-computational network with FDTD mesh simulating the PML medium (non-uni fied or split algorithm) [START_REF] Eswarappa | Implementation of Beren ger absorbing boundary conditions in TLM by interfacing FDTD perfectly matched layers[END_REF]-[10]. •using a TLM node in the PML medium (unified algo rithm) [11 ].

The second approach has the potential advantage of avoiding field spatial interpolation which is required to interface TLM and FDm algorithms. Indeed, it is well known that field components are not defined at the same location in both cells. One objective of this paper is to present a new M}M node that can simulate PML media with the objective to use a unified algorithm for PML used as ABC. To develop this new node, the PML field equations are rewritten such as to establish Maxwell's equations pertaining to an anisotropic lossy medium with sources and characterized by real diagonal electric and magnetic conductivity tensors. Magnetic and elec tric current source densities are controlled by field spa tial time variations [l l ]- [START_REF] Zhao | A general approach for the development of unsplit-field time-domain implementations of perfectly matched layer for FDTD grid truncation[END_REF]. Then, using a rigorous field formulation [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF], equations are approximated while keeping 1LM formalism, i.e., in terms of incident and reflected voltages that travel in the network of inter connected nodes which simulates the PML medium.

II. THE TLM METHOD

The 1LM method is a time-domain numerical tech nique that models an electromagnetic field in a compu tational domain by means of voltage impulses that travel in a network of interconnected transmissions lines. More fundamentally, the 1LM network is able to model Maxwell-Ampere (MA) and Maxwell-Faraday (MF) laws in presence of electric 1; and magnetic �' current den sities in lossy media [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF]. These two equations in local form are given respectively by:

o l v /\ Tl = £ �7 + 8e "F + -;: � /\ all /\ -rl n (2) Y' M| MĊ =-µ a t -MU m H -]m
in which £, µ, 8 e and 8 m are real diagonal tensors of permittivity, permeability, electric and magnetic conduc tivities, respectively, given by:

(3) £ = E 0 diag(E u , E YY' E zz ), il = µ 0 diag( µu , µYY' µzz ), (4) 8e = diag( MU ex' (Jey' cre ) 8 m = diag( MUm x ' Mmmy' Mmmz ), and £ 0 , µ 0 are vacuum permittivity and permeability, respectively. The computational domain is then filled with a Car tesian mesh of cells with dimension L1 x , L1 Y , L1 z . In each cell, Fand dtangential components are sampled at the center of every face as well as total field components at the center, as illustrated by Figure 1 a and 1 b, respecti vely. There is a total of 15 samples for each field which can be combined to generate eighteen incident and reflected voltages which are contained in vectors ctand F, respectively. In accordance with the numbering used for the symmetrical condensed 1LM-node (scN) [l], electric field components and voltages are numbered correspondingly, as illustrated by Figure 2a. For ins tance, samples Ei and Hf are combined in the follo wing way [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF]: [START_REF] Chen | Absorbing and connecting boundary conditions for the TLM method[END_REF] { al } = _l (L1 E'"!:.Z '1.ffi)

All expressions concerning incident and reflected vol tages are included in Appendix A. Now comes the time-evolution of the TLM algorithm that computes the fields in the structure. First of all, sup posing that at time (n -1/2) � t the incident voltage vec tor ct is known at every node, then, the reflected voltage 1f from every node can be computed via S at time (n + 1/2) � t as :

(7) 't< n + 11 2) = S {(<11 -11 2 )
in which S element expression can be found in [lJ for the lossless case. Finally, reflected voltages become incident voltages to the neighbouring nodes for the next iteration.

Although the procedure described by [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] is very simple, it is relatively costly in terms of computer expenditure and an accelerated procedure is used ins in Appendix B. Then, reflected voltages b at time (n + l/2) � t are determined from the previously compu ted {P, 7/n} and the incident voltage 71 which is known at time (n -1/2)Mĭ via expressions written in Appendix B. Hence, the matrix Sis not used, which reduces the number of operations per iteration [START_REF] Tr | Efficient Computation Algorithms for TLM[END_REF]. In lLM, six degrees of freedom are available to consider alternate cells : six stubs or six link-line impe dances or any combination thereof. With the objective to reduce the number of stubs and increase the maxi mum time-step value, Scaramuzza and Lowery [START_REF] Scaramuzza | Hybrid symmetrical conden sed node for the TLM method[END_REF] pro posed the hybrid scN (HSCN) lLM node that uses only three stubs but needs three different link-line characte ristic impedances. Several years later, a new super scN (sscN) lLM node that has no stubs was proposed by Trenkic et al. [17]. Constitutive parameters and/or paral lelpiped cell dimensions are included in six diM¥erent values of the link-line characteristic impedance. The drawback of these new cells is the interface conditions that must be applied between two different media. However, the reduction in terms of computer costs makes them potentially attractive.

Finally, boundaries are accounted for by using a proper reflection coefficient in the adjacent lLM arm. The value of the reflection coeM §cients should translate the condition on tangential field components on the boundary. For instance, it is -1 and + 1 for perfectly conducting and a magnetic wall, respectively. For other types of physical boundary conditions such as those that yield frequency dependent reflection coefficients, the procedure becomes more complex, yet possible to implement in the lLM scheme (eg., [START_REF] Ney | Modeling of highly conducting boundaries with TLM[END_REF]).

III. ABSORBING LAYER THEORY

111.1.

Matched layer (ML)

The ML technique consists in enclosing the computa tional domain (for instance with parameters Er and µ r ) by a lossy passive medium whose constitutive parame ters should fulfill the following relation [START_REF] Berenger | Calcul de la diffraction a !'aide d'une methode aux differences finies[END_REF] :

�=� = S (8) E h µ r µO
where cr and cr are the electric and magnetic conductie m vities of the ML medium, respectively. This lossy medium can be shown to have an intrinsic impedance that corresponds to :

As a result, it is perfectly matched for waves that impinge at normal incidence for any frequency. The apparent reflection coefficient (ARC) is evaluated by considering the wave returning from the back of the layer terminated by an electric wall and expressed in dB by: -�20 (10) ARC= 20 log (p = e c ) where o is the layer thickness and c the velocity of light in a medium with intrinsic impedance Z m .

It has been found that purely numerical reflections exist for a wave propagating through two media with different conductivities but same intrinsic impedance. The level of reflection is proportional to the conducti vity difference at the interface. To circumvent this pro blem, a tapered loss profile is used in the layer, at a price of increasing its thickness o. This approach allows a small difference in conductivities across the interface between the computational domain and the layer, while keeping the same ARC which is now expressed in dB as :

_i. J,s1 �i <lz) (11) 20 log (p = exp c z1
for a wave propagating in the (Oz) direction and where z1 is the coordinate of the interface in the xy-plane. For the numerical implementation, the conductivity spatial distribution is stepwise approximated such that each cell has constant conductivity value. Thus, the layer can be directly simulated by standard loaded scN or other more recent nodes such as the HSCN or the sscN.

111.2.

Perfectly matched layer (PML)

The PML technique is an extension of the ML and theoretically allows a perfect match for any case of inci dence. To achieve good absorption in any direction, Berenger [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] proposed to split all field components into two sub-terms. For instance, the electric field x-compo nent is written as E x = Ex y + E, z -Substituting �nto (1) and ( 2) for sourceless situation, leads to a set ot twelve equations that can be written in a generic form as [START_REF] Berenger | A perfectly matched layer for fre e-space simu lation in finite-difference computer codes[END_REF] :

(12) (13)
where (s,ri,sl E {(x, y, z),(y, z, x), (z, x, vl} and (cr r, cr* r ) which are fictitious electric and magnetic 11•� 11.s . h Id conductivities, respectively. These parameters s ou satisfy [START_REF] Berenger | A perfectly matched layer for fre e-space simu lation in finite-difference computer codes[END_REF]. Originally, the time evolution of every suh term is computed in the PML medium using the classical FDTD Yee's scheme. Thus, since these sub-terms are defined at the same location as the corresponding total component, it is straightforward to transfer the algo rithm through the interface between the PML medium and the FDTD computational domain. Hence, a fully uni fied algorithm can be achieved when FDTD is used.

One can show that a grouped form of (12) and (13) allows the interpretation of a PML medium as a classical lossy medium with controlled current source densities

[ll]- [START_REF] Zhao | A general approach for the development of unsplit-field time-domain implementations of perfectly matched layer for FDTD grid truncation[END_REF]. For instance, with E x = E" + E x z' (12) can be written in terms of (E x , E9) and (E � , E xz ) :

Hence, one retrieves (1) for E x , with cre x equal to cr 2 and cr y in ( 12) and (13), respectively. In addition, the electric current source density J e x can be identified as (cr y -cr 2 )E xy and (cr 2

-cr y )E x z in (12) and ( 13), res pectively. Also, one may notice that J e x is controlled by they-spatial variation of H 2 in (12), and by the z spatial variation of H Y in (13).

According to the above observations, one can consi der the extension of the scN to PML medium simulation, if the controlled source evolution can be accounted for by adding appropriate voltages to the basic SCN. This concept was successfully applied to the two-dimensio nal nM node [11].

IV. IMPLEMENTATION IN THE TLM ALGORITHM

IV.1. ML implementation

Generally, S-parameter computation requires two simulations. One has the purpose of generating the field solution without the discontinuity (reference or thru line). The other simulation performed with the disconti nuity inserted, generates the total field solution. This allows the separation of incident and reflected voltages on the guiding structure that feeds the discontinuity under investigation, from which S-parameter can be computed. A major advantage of time-domain simula tions is that S-parameters are calculated on a wide fre quency band by exciting the structure with a narrow impulse (generally gaussian signal). Then, a Fourier transform of the time-domain response is performed to obtain the frequency characterization of the structure.

For the above applications, several precautions must be taken. First of all, input and output ports need to be matched (infinite line) over the frequency range of inter est. Thus, ABC must be placed behind the excitation plane and where a matched load would be connected. Then, points at which the time responses are recorded and should be located far enough from the discontinuity and the excitation plane in such a way that only the dominant mode prevails.

Quasi-TEM mode in microstrip or coplanar wave guide propagates nearly along the longitudinal direction. Thus, these modes impinge almost normally on the plane where ABC is enforced. These features fulfill the basic assumption of using the ML technique. For planar structures, metal strips and substrates are simply exten ded through the absorbing layer. Then, the mesh is sim ply terminated in the layer by the intrinsic impedance of the adjacent medium whether it is air or dielectric mate rial (an example is shown in Figure 3). It has been found that better matching is achieved by terminating the layer this way than by a conducting wall as it is usually done in FDTD simulations. Indeed, by terminating the compu tational domain directly using such conditions, return loss better than 30 dB is achieved. In fact, residual reflections are mainly produced at the interface due to the slight deviation of quasi-TEM modes from normal incidence condition. At the other limits of the computa tional domain (except if a ground plane is present) the trivial condition, which consists in terminating the nM mesh by the intrinsic impedance of the adjacent medium, is enforced. If these conditions are applied at locations suM¨ciently remote from the planar structure where the energy is confined, results are not signifi cantly perturbed. For instance, locating the limit about six to eight cells above the substrate has been found to yield satisM£ctory results for all practical situations. It is important to mention that using ML instead would not be practical since evanescent waves (reactive fields) predo minantly prevail above and on the side of the structure. The last condition to be aware of, pertains to the case of non-homogenous matched layers. A condition obtained from Snell-Descarte's law must be taken into account to balance absorption in air and dielectric [START_REF] Bahr | Application of the PML absorbing condition to the FDTD analysis of microwave circuits[END_REF]: [START_REF] Scaramuzza | Hybrid symmetrical conden sed node for the TLM method[END_REF] cr l E l cr 2 E 2 where cr1, E1 and cr 2 , E 2 are adjacent media electric losses and permittivity, respectively. Note that magnetic losses remain the same in all media when dielectric materials only are involved. Finally, losses in the absor bing layer are directly included in the TLM node scatte ring matrix. In the case of cubic cells, this is done by adding loss (infinitely long) stubs whose characteristic conductance and resistance is given by g e = Z 0 cr e dl and rm= crmM/Z 0 , respectively [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF]. ;:, n = [ -cr; Hry -cr; H H or (20) 7?i = [ cr z E <z cr z E yz O f . 7!1 = [ cr ;Hxz cr; H v z O f where T denotes transposed vector. Examining ( 17) to [START_REF] Bahr | Application of the PML absorbing condition to the FDTD analysis of microwave circuits[END_REF] shows that there are different ways to group PML equations. Each yields a different numerical scheme that has specific dispersion characteristics. It has been found that the algorithm implemented in TLM and expressed by [START_REF] Berenger | Calcul de la diffraction a !'aide d'une methode aux differences finies[END_REF] and (18.a) yields the best performance in terms of numerical stability and, consequently, is used in the pre sent paper. For instance, the sub-term E xz involved in [START_REF] Herring | Accurate modelling of zero thickness septa with the symmetrical condensed node[END_REF] needs two samples on faces perpendicular to z-direction, i.e., at the same locations as E x samples with subscript 2 and 9 (see Fig. 1 ). This is the reason why one uses the nota tion (E�2, E�2). Likewise, E yz samples involves in (23) are written (E� z , En, H . samples involves in (24) are written (JF, 5 Y, H 7 xy ), and, finally, H samples involved in yx yx yx (25) are written (H 6, H 1. 0 ) .

Proceeding further with the sub-term E xz , [START_REF] Herring | Accurate modelling of zero thickness septa with the symmetrical condensed node[END_REF] is numerically integrated between time (n -1/2)dt and ndt as well as between the spatial limits (k -1/2)Mċz and (k + 1/2)t,.z, yielding:

(26) t,. E" -!._ (d e z + t,. ,;• xz x xz 2 This yields (30) the first expression of the second group of the algorithm described in the next section. A spatial integration procedure between Ki1z and (k + 1/2)Az, allows one to establish a relation similar to (29) from which b 2 0 can be then computed. Finally, the same pro cedure is applied to all sub-terms yielding the generic PML-11.M cell which comprises 18 +8 = 26 voltages as illustrated in Figure 2b.

At this point, it is important to stress that the 11.M node developed above uses only set of conductivities (i.e. cr z , O"! * 0 ). Some auMÈors refer to uniaxial PML media.

For simulating layers that completely surround the com putational domain, it is customary to use all conductivi ties at corners, where three uniaxial PML (i.e., perpendicu lar to x, y, z axes, respectively) meet. A more general 11.M node that accounts for all conductivities can be develo ped following exactly the above procedure. Such a node would include a total of 30 voltages. However, the pro blem of using multiaxial PML at corners is still under investigation and, consequently, not dealt with here.

IV.2.2 PML-TLM algorithm

As for the standard MňM node, the PML-11.M algorithm uses the accelarated procedure briefly described before.

According to [START_REF] Trenkic | Theory of the symmetrical super-condensed node for the TLM method[END_REF], the only modification pertains to accounting for the source terms at the cell center and at time nM´t, by using typical expression such as ( 27 Finally, to complete the procedure, the definition of vol tages that model controlled sources is given in Appendix C.

IV.2.3 Boundary conditions

The algorithm for both lLM and TLM-PML nodes are now established. The last point pertains to determining incident voltages at nodes adjacent to the interface bet ween the lLM and lLM-PML network. As far as the first eighteen voltages are concerned, they are treated by the basic lLM scheme [START_REF] Hoefer (y{ | The transmission line matrix method theory and applications[END_REF]. For the remaining voltages, let us first consider the interface perpendicular to the z-axis (Fig. 4.a). In the lLM domain (region 1), there are two reflected voltages from the adjacent nodes, say b� l and b�

) , that impinge on the interface at time (n + l/2)ill .

On the other side in the PML medium (region 2), there are four voltages b� J , b�J , b� ) and b�i reflected form the PML-lLM nodes adjacent ot that interface (see Fig. 2).

Concerning (b� l , b� > ) and (b� l , b� \ they are simply transfered across the interface for the next iteration [START_REF] Hoefer (y{ | The transmission line matrix method theory and applications[END_REF].

On the other hand, b�J and b�i are not transfered through the interface as the lLM node in region 1 has only 18 voltages. However, a�J and a�i incident ot lLM PML node must be computed. First of all, note that the knowledge of (b� ) , b� ) ) and (b� ) , b� » allows one to compute the field components tangential to the inter face. Then, the magnetic field components H

(2J = b ( 2 J -C z � H !n + 11 2) a 2 0 2 0 x o y y = b (2) -c (b (2) -b (l) ) 2 0 x 9 2 c/2) = b� 2) + C Z � H ( n + 11 2) 22 L2 xOx x = b (2) + c (b (l) -b< 2 >) 22 x 4 8
Now, let us consider the case where the PML-lLM network is terminated by the vacuum intrinsic impe dance Z 0 in the plane perpendicular to x-axis. As the characteristic impedance of the lLM link lines is also Z 0 , such termination means zero reflection in the link lines (note that the global reflection is not zero). At time (n + 1!2) � t b3, b 6 and b 2 5 impinge on the boundary (see Fig. 2b). Thus, a 3 and a6 are zero and one can easily show that from the definition (A. 1) in Appendix A, the component Ez is given by: (47 The above procedures can he easily applied to other reflection coefficient values when a 1LM-P\1L algorithm is used.

) A E (n + 11 2) = / u z z ')6,
V.

RESULTS

V.1. Matched layer results

V.1.1.

Influence of the conductivity profile

The reflection coefficient is computed by using the SCN for different profiles (linear, geometric, . .. ) ARC and geometries (microstrip and coplanar). In addition, sub strate with low and high permittivity is considered. The parabolic profile has been found to yield the best perfor mances:

(49)

a (z) = a max (i-)2
However, as mentioned before, the discretization of the profile leads to an abrupt change in the conductivity value near the interM£ce that provokes numerical reflec tions. An alternate profile, called semi-parabolic profile, has a factor five times lower than the parabolic profile, which reduces the conductivity step value near the inter face and, consequently, the numerical reflections. Howe ver, since the penetrating wave must be suM¦ciently atte nuated to obtain the same ARC, the conductivity profile has a much larger conductivity factor at the end of the layer. Simulation results shown below are performed with a 15-Al layer thickness, where Al is the space step in the direction of propagation. This choice was found to be the best compromise between reflection and memory requirement.

Return loss extraction for microstrip and co-planar structures is performed using the geometries shown in Figure 5 good spatial sampling must be used to describe the slot for co-planar line. If this condition is not fulfilled, the quasi-TEM mode is not accurately described, which results in a higher reflection level from the matched layer interface. It was found that three cells for descri bing the slots has produced sufficiently good results. Using a coarser mesh, return loss may get 20 dB lower in the examples presented below. Finally, excitation consisted of a Gaussian pulse whose variance was selec ted according to the desired frequency range characteri zation. Sources were located between the ground plane and the strip while they were across the slots in the co planar lines (even mode). First of all, simulations were performed on planar structures shown in Figure 5, wiMÈ Er = 2. 2. For the microstrip case (Fig. 6a), the optimum layer corresponded to a 40 dB ARC (parabolic profile) and a computed reflection coefficient (CRC) which is always lower than -52 dB. A slight decrease in absorption (still good), is obtained in the co-planar case illustrated in Figure 6b ( -30 dB ARC with -50 dB CRC). For both cases, there is some substantial improvement as compared to stan dard matched boundaries. Finally, as mentioned before, less reflection than the ARC value occurs because of the matched loads used at the back of the absorbing layer. Also, one can observe that the CRC increases as the ARC decreases. This can be explained by the fact that the adjacent conductivity ratio at the interface increases with ARC (see (11) with (49)).

Results have shown that the substrate permittivity did not influence the absorbing perMXrmance of the ML up to E r = 6. However, as the permittivity increases further, the absorption perMXrmance slightly deteriorates. For ins tance, on the frequency range (0 -100 GHz), results pertaining to two different substrates (Er = 9. 9 and E r = 12. 8) are reported in Figure 7. Optimal layer is obtained for a -30 dB ARC, yielding a -46 dB CRC. Results (not shown) are very similar for the co-planar line. Note that one can lower low-frequency reflection by decreasing the ARC down to -100 dB (see Fig. 7). In this case, the CRC is lower than -50 dB on the 0-30 GHz band. This result can be explained by the fact that quasi-TEM modes are closer to normal incidence in the low MZequency range. On the other hand, semi-parabolic profile optimi zation may be performed on the high MZequency range. Results compared with parabolic profile are given in Figure 8. A -50 dB CRC is obtained on 30 to 100 GHz band, using a -100 dB ARC.

--• FIG. 7. -Matched layer ere (parabolic profile) for several arc values in the microstrip line case, with high permittivity substrate (E, = 9. 9).

... '

Fa cteur de reflexion calcu/e du milieu ML (projil parabolique) pour diff erent es valeurs du fa cteur de reflexion apparent ARC dans le cas du microruban avec subsrrat a fo rte penniuivite relative (9, 9).

From the above results, when dealing with high per mittivity substrates, one can improve the absorption by using parabolic profile in the low frequency range and semi-parabolic profile in the high MZequency range. However, results have shown that improvement with semi-parabolic profile is no longer efficient for co-pla nar guides. Finally, one should mention that the perfor mance of ML for planar circuit's characterization is suffi cient as a reflection level below -46 dB is very good for S parameter extraction. Fac/eur de reflexion calrn/e du milieu ML (projil parabolique el semi-parabolique) dans le cas du microruban avec subs/rat a Jorie penni//ivite relative (9.9).

V. 1.2. Absorption perfonnance for different

TLM mode simulations

The various TMĻM scN nodes have diM¤erent dispersion characteristics and, in addition, some of them require additional interface condition between different media [START_REF] Ney | Modeling of highly conducting boundaries with TLM[END_REF]. Figure 9 shows the comparison between the various SCN in the case of the microstrip line, with two different substrate permittivities (E r = 2. 2 and 9. 9). It can be observed that the overall performance of the sscN is not as good as for other nodes, especially for high permittivity (Fig. 9a). These differences could be explai ned by the difference of the sscN dispersion characteris tic as compared to the HSCN and the scN. MRis is confir med by Figure 9b in which the performance shown by the HSCN and the scN tends to be identical, likewise their dispersion characteristics for low permittivity values. Further investigations are necessary to find explanations on this issue.

V. 1.3. Planar circuits S-parameter com putation

The S-parameters are computed with the basic SCN for the microstrip line loaded by an open stub shown in Figure 10. Results are compared with measurements in Figure 11. Excellent agreement is found in both magni tude and phase. Note that a 2 % resonance frequency shift MXr s ll is observed. This phenomenon is well known in TI .M modeling and is due to coarseness error as explained further. The resonant frequency shift was 2. 1 % and 2. 4 % for the HSCN and sscN, respectively. They confirm the slightly lower accuracy produced by the sscN. However, if one considers the very close per formance obtained by the three nodes, one should men tion the suhstantial memory gain of the HSCN and sscN as compared to the scN (17 % and 33 %, respectively). To conclude from the ahove results, the HSCN seems to be the best compromise in terms of memory and accuracy. Discontinuite en microruban sur un substrat de hauteur 254 µm et de pennittivite relative 2,2.

V.2. PML results

V.2.1. Conductivity profile tests

The proposed model is tested for validation in the case of wide-band matched load simulation for TE 1 0mode propagation in a empty rectangular waveguide (WR-28 with a = 32 mm = 36 111). This canonical case ,� is very relevant to testing absorbing conditions. Indeed, all cases of incidence can be achieved (except near to normal), including evanescent modes. The guide is exci ted hy a current sheet [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF] which is located near th PML interface (Fig. 12). The conductivity profile used in the PML is parabolic with crmax = 25 Sim. The time excita tion is a Gaussian pulse to achieve wideband characteri zation centered on the dominant mode operating range.

Output points from which time sequences are extracted are P 1 and P 2 (fig. 12). A FFT procedure is applied to obtain the frequency characteristics. Results shown in Figure 12 yield a reflection coeM¦cient which is lower than -58 dB over the 25-40 GHz frequency range. This shows some significant improvement as compared to reflection levels reported with other ABC used in ll.M [START_REF] Eswarappa (ch | Absorbing boundary condi tions for time-domain TLM and FDTD analysis of electromagnetic structures[END_REF]. 12. -Relle<:tion coefficient from the PML medium used as matched loads in a WR-28 re<:tangular waveguide (umax = 25 Sim. L = 60 !:J.J, 8 = 20 !:J.J, a= 36 !:J.J, P1 and P 2 Iocated at one !:J.J from the PML interface).

Facteur de reflexion pour la propagation du mode dominant dans un guide d' 011des tennine par des PML ( uma.' = 25 Sim, L = 60 !:J.J. 8 = 20 !:J.J, a = 36 !:J.I. p I et Pz places a un !:J.J de !'interface).

V.2.2. S-parameters extraction

The next canonical case is the wide-band characteri zation of a capacitive iris in a WR-28 rectangular wave guide in which PML technique described above is used for matching. The magnitude of S 11 is shown in Figure 13 in which FDID [10] and the proposed unified TI.M tech nique are compared with a closed-form solution used as reference [START_REF] Marcuvitz | Wa veguide handbook[END_REF]. Some slight discrepancies can be obser ved with the reference for both numerical methods, in particular the TI.M. This is due to the presence of sharp edges near the iris where field gradients are high and for -Computed reflection coefficient (cRc) for a capacitive iris in a WR-28 rectangular waveguide terminated by PML (a = 66 tJ.J for the fine mes h).

Facteur de reflexion calcule po11r 1m iri.1• capacitif da11.1• wz g11ide d'ondes wr-28 tennine par des pm! (a = 66 Ill pour le mail/age fin).

which the TI.M produces some additional error as com pared to FDID [START_REF] Herring | Accurate modelling of zero thickness septa with the symmetrical condensed node[END_REF]. In order to reduce this error, an extrapolation technique is applied. From two successive simulations using finer mesh size, results for the zero mesh-size limit are extrapolated. This technique yields excellent results as shown in Figure 13. It is important to notice that PML interface is located close to the iris where evanescent higher-order modes prevail. It seems to indicate that PML can absorb evanescent modes. Finally, some cases of instabilities for iris simulations were observed when fine mesh was used. This problem could be resolved by decreasing the time-step. This resulted in an increase of the number of iterations and, thus, of the computation time. The source of instabilities is currently under investigation.

The next case for validation pertains to planar cir cuits shown in Figures 14 and15. In both cases, the PML are used as matched loads in the direction of propaga tion. In other sides, both PLM and PML-11.. M networks are terminated by Z0. Results for S-parameters are presented in Figure 14, where they are compared with those obtai ned with the FDID-PML unified algorighm and measure ments [START_REF] Shee | Applica tion of the three-dimensional finite-difference time-domain method to the analysis of planar microstrip circuit s[END_REF] . Some good agreement can be observed up to 15 GHz. Beyond that frequency, the comparison should be considered with caution as the authors men tioned that the measurements were not reliable.

Finally, the example of the patch antenna (Fig. 15) is used to test the ability of the PML technique when the layer is located in the very near field (evanescent waves). For this purpose, the PML interface is located within two cells from the antenna edges perpendicular to the microstrip line. Results pertaining to the input reM©ection coefficient are shown in (Fig. 15). Compari son with PML interface located 10 cells from the antenna showed only a slight diM¥erence of 0. 3 dB at the first resonance. One may conclude that PML can be located very near the structure without affecting the accuracy -Input reflection coefficient of a patch antenna fed by a microstrip line (substrate : h = 0. 794 mm, £, = 2,2) computed with PML used as absorbing boundary conditions (ABC).

Fa cteur de reflexion ii /'entree d'une antenne plaque alimentee par un micromban (substrat : h = 0,794 mm, £, = 2,2) et avec des PML comme conditions aux limites ab.wrbantes.

that pertains to input impedance computation. Conse quently, the computational domain can be reduced which produces some significant reduction of the com putational cost, as compared to other ABC techniques.

VI. CONCLUSION

Simulation of absording media with a unified lLM algorithm was presented. The matched layer technique was found to be very attractive in applications that involve wide-band matched loads of guides that operate in quasi-MŇM situation. In this case, the absorbing medium can be simulated directly with the basic MºM algorithms (scN, HSCN and sscN) without modification. When quasi normal incidence assumption is no longer valid, the per fe ctly matched layer technique (PM) must be rather used. A new lLM cell capable of simulating uniaxial PML media was introduced. It is an extension of the classical lLM scN and, thus ensures the continuity of the algorithm across the lLM-PML network interface. The concept of sources controlled by local fields was used and translated in terms of additional voltages (stubs) connected to the basic scN. An accelerating procedure was applied for reducing the number of operations per iteration and boun dary treatment MYr the new PML-lLM cell presented.

Simulation results in the case of planar circuits such as microstrip and co-planar lines, showed that ML tech nique yields 45 dB return loss, at least, over the fre quency range up to 100 GHz and relative permittivity substrate up to 6. For a higher value of Er, return loss decreases but is still at a good level for practical pur poses. All MºM node simulations of ML showed numeri cal stability. Finally, the lLM sscN node is less perfor mant than other cells but requires less memory than the other nodes.

The new lLM-PML cell was found to be eM¨cient MYr simulation of general matched load. In addition, it was found that the PM was able to be located where evanes cent waves prevailed which was not the case for other absorbing techniques currently used in lLM. In certain solutions, numerical instabilities occurred but could be circumvented by reducing the time step at the price of increasing the number of iterations. The source of this problem is currently under investigation.

Manuscrit reru le 15 octo bre 1997 (B 1 ) A £11 A { C ( + + + + 11y )<11 -11 2 ) computed above in the reflected voltage vector /}z c where (S, 'fl, s.) E {(x, y, z), (y, z, x), (z , x, y)} . Finally, one should notice that the time step Lit is selected as the maximum value for which normalized immittances are real positive in every cell. This can be interpreted as the numerical stability criterion for the TLM method.

  Key words : Microwave circuit : Electromagnetism : Absorbing medmn ; Modelmg ; Numencal method : Transmission line : Matrix method : Microstrip line ; Coplanar line : Impedance matching TECHNIQUE DES COUCHES ABSORBANTES ADAPTEES, APPLIQUEES A L' ANALYSE RIGOUREUSE DES COMPOSANTS MICROONDES ET EN BANDE MILLIMETRIQUE, PAR LA METHODE TLM Resume :

b 2 xFIG. 1 .

 21 FIG. 1. -Field samples for the basic TLM (SCN) cell : (a) electric field. (b) magnetic field. E chantillonnage des champs dans la cellule de base TI.M ( SCN) :(a) champ electrique, (b) champ magnetique.

  tead. First, field components CE\ lf n } at the center of the cell are computed at time n � t. They depend on the (b) Fm. 2. -Incident voltages in the TLM cell; (a) basic scN, (b) new PML-TLM cell. Tensions incidentes dans la cellule TLI/: (a) cellule de base SCiV. (bi nouvelle cellule PML-TIN. source terms I' and I' via expressions that are included e m � .

  FIG. 3. -Matched layer configuration for microstrip geometry.Configuration de la couche adaptee (ML) pour une geometrie de tvpe microntban.

IV. 2 .

 2 The PML•TLM cell Let us consider the interface between the computational domain (E r , µ r ) and the PML medium (E r , µ r , cr x , cr _ ;, cr y , cr;, cr z , cr;) perpendicular to the z-axis. In this case, there is theoretically no reflection at the interface if cr = cr = cr * = cr * = 0 and cr cr * satisfy. By setting x y x y z ' z cr y = 0 into (12) and (13 ), one can notice an electric loss and current source terms in (14.a). The time evolution of the latter, given by (14.b), is controlled by the magne tic field component normal to the interface. On the other hand, there is only a source term in (15.a) controlled by the spatial variation of magnetic field components tan gential to the interface (see (15.b)). By applying the same procedure to other field components, one can derive from (1) and (2) the following relations: (17 .a) (17.b) ( 18.a) (18.b) v A Fl= ££ a !! + 8 7t + I' r 0 at en en v A Fl= £ £ a!! + I' r 0 -ac-et v M| 7t = -µ µ a ll -8 Fl-I' r 0 at m n m n in which subscript e, m pertain to electric and magne tic, respectively and n, t to normal and tangential com ponents to the interface, respectively. Also , 8 en = diag(cr z , crz, 0),8m n = diag(cr ;, cr;. 0) and source vec tors are given by: (19) I' -[-E -E O ] T en -cr z xy cr z xy '

IV. 2

 2 .1. Numerical scheme One can notice that 8 e = 0 and crmn= diag(cr ;,cr;,o) present in (17.b) and (18.a) can easily be added in the classical algorithm (see Appendix A), by setting A x = AY = A z = 1, B z = 1 and B x = BY = B in (B. 1) to (B. 6), where: 4 /\ 1 cr z *s cr s (21) B = --, where , R = ---= Z 0 _ z _ , and in which conductivities are deduced from (8). According to relations (B. 1) to (B. 6) in Appendix B, the source terms must be determined at the cell center and time ndt . However, they are easily computed as they are proportional to the sub-terms E xz ' E yz ' H ry, H yx which are governed by the following equations : (22) aE xz + z Co cr z E = -Co z aHY at o T x z E , aE yz + z Co cr z E = Co z aH X at O E r yz E , O az z aHry = _ c 0 � 0 at µ r ay Z o aH vx = � aE at µ r ax In order to model these sub-terms, one adds two samples per sub-term to the cell illustrated by Figure 1.

+

  Q. /j._ 1;n 4 . \"""' x z /\ /\ where C can be found in Appendix B, G = R is given x by (21 ), and Tis a time-delay operator defined by TU" = [f 11-1 >. Samples are taken at time (n+ 112)dt unless otherwise indicated. Thus. the operation TV denotes the quantity U at time (n-1/2)M´t. Hence, (26) yields the value of E: z as : (27) d r;n � (a + a ) (n -li2 J x' " "' xz = 4+ R 19 2 0 in which : (28) { ;j = t (L1 x f2 z :t Cl 0 L1�) and { a 2 0} = _!_ (L1 pz +. C 7 A ffi) b 2 0 2 x 9 .t'"'()'-' y 9 The last step to evaluate the contribution of Exz into � at the cell center, is to insert (27) into (B. 1 ). This yields (30) which is the first algorithm expression des cribed in the next section. Other relations are derived in a similar fashion. To completely account for the influence of the sub term E x z ' remains the evaluation of the evolution of b 19 and b 2 0 from incident voltages and field components at the center of the cell. To achieve this, one again inte grates numerically (22) but, this time, from time (n-1/2)MĬ and (n + 1/2)i1t and spatially from (k-1/2)Az and kAz (from the face to the cell center). This yields: /\ + Si_ L1fi n = + C./; i1 Ff' 4 . x z . O yy In order to establish a new expression that contains b 19 at time (n + 1/2)i1t as a function of a 19 , a 2 0 at time (n-1/2)i1t and the magnetic field component H Y (�) at the cell center, (27) and (28) are introduced into (29).

�P 4 25 24 ( 35 )

 42435 ) and new relations for reM©ected voltages {b 1 9 , ... , b 26 }. Hence, relations (B. 1) to (B. 6) see Appendix B) become:(30) i1 P = { C (a + a + a + a + Y a )<n -11 2J x x x I 12 2 9 sx 13 ( a + a )<n -112) } 4 + G 19 2 0i1 P = { C (a + a + a + a + Y a ) <n -11 2 ) = { C (a + a + a + a + Y a )<n -11 2 J z z z 5 7 6 I 0 sz 15 Z oi1P� = s{ D x ( -a 4 + a 8 + a 5a 7 + z a )<n -11 2) + G( a + a 7 )<n -112) } sx 16 4 23 _4 Z oi1 Y H� = s{ D Y (a 2a 9 -a6 + a 1 0 + z a )<n -11 2) + G ( a + a )<n -11 2)} sy 17 Zoi1 z H � = { D z (-a 1 + a 12 + a 3a 11 /\ /\ + z a ) <n -1 / 2)} sz 18 where Band G = R are given by (21). Then, {b 1 9 , ... , b 26 } are computed using (B. 10) in Appendix B and {b 19 , ... , b 26 ) from the expressions below: (36) b (n + l / 2 ) = -C Z A un + B (n -11 2 ) n + 1 / 2) = + C Z A H n + B (n -l / 2) + 11 2) = -C Z d Ff' + B (n -

  + l / 2 ) = +D d p + (n -112) 24 x x z a 24 (42) (43) b (n + 1 / 2) = + D d p + (n -1 / 2) 25 y x z C l25 b (n + 112) = _ D d p + (nl / 2)

y y 9 2 ,

 2 x and H Y at time (n + 112) are given by (see (A. 1) in Appendix A): (44) z � Mė t.n + 1 / 2) = b ( l )b ( 2! n + 11 2) = b ( 2) _ b ( l J 0 Then, a�J and a�i can be computed from their definition (see (C. 1) in Appendix C):

Region 2 :FIG. 4 .

 24 FIG. 4. -Voltage at PML boundaries (a) interface between TLM and PML-TLM networks (b) at the limit of the PML-TLM network terminated by Z 0 Tensions incidente. 1• ii la fr ontiere du reseai1 PML-TIM: (a) inter fa ce entre le reseau TIM classique et le reseau PML-TLM (b) reseau PML-TIM te1min par Z0

FIG. 5 .

 5 FIG. 5. -Geometry used for ML absorption performance : (a) microstrip line. (b) coplanar line.Geometries utilisees pour le calcul du fa creur de reflexion des ML.(a) microruban, (b) ligne coplanaire.

FIG. 6 .

 6 FIG. 6. -Matched layer ere (parabolic profile) for several arc values (Er = 2. 2) : (a) microstrip line. (b) coplanar line. Facteur de reflexion calcu/e du milieu ML ii profil parabolique pour dif f erenr.1• valeurs de ARC (Er = 2, 2) : (a) microruban. (b) ligne coplanaire.

Fm. 8 .

 8 -Matched layer CRC (parabolic and semi-parabolic pro file) for microstrip with high relative permittivity substrate (9. 9).

FIG. 9 .

 9 FIG. 9. -Matched layer reflection comparaison of the three TLM nodes : (a) microstrip E , = 9. 9, ARC= -100 dB. (b) microstrip E , = 2,2, ARC = -40 dB. Comparaison des valeurs du f acteur de reflexion calcule pour le milieu Ml entre Les trois cellule. 1• TIM: (a) microruban, E , = 9. 9, ARC = -100 dB. (b) microruban, E , = 2,2, ARC = -40 dB.

  FIG. 10. -Microstrip discontinuity on a 254-µm thick substrate with E , = 2,2.

  FIG. 11. -S11 scattering coefficient of the discontinuity illustrated in Figure 10, (6.1 = 124 µm, 6.t = I. 97. 10• 1 3 s, 2500 iterations) with SCN-ML matched load simulation : (a) magnitude, (b) phase. Parametres de repartition S11 de la discontinuite de la figure JO ca/cutes avec des charges adaptees simu/ees par un milieu ML et la cellule SCN (6.1 = 124 µm , 6.t = 1,97. 10-13 s, 2500 iterations): (a) module. (b) phase.

  FIG.12. -Relle<:tion coefficient from the PML medium used as matched loads in a WR-28 re<:tangular waveguide (umax = 25 Sim. L = 60 !:J.J, 8 = 20 !:J.J, a= 36 !:J.J, P1 and P 2 Iocated at one !:J.J

  FIG.13.-Computed reflection coefficient (cRc) for a capacitive iris in a WR-28 rectangular waveguide terminated by PML (a = 66

FIG. 14

 14 FIG.14.-,\ '-parameter' of a micro;;trip filter (�uh;;trate : h = 0. 794 mm, E, = 1.2) computed with PML matched load simulation.Parmnetre.1• S d'un filtre micromban (.1•11b.1•trat : h = 0. 794 mm. E, = 2.1) te11ni11e par desPML.

  FIG.15. -Input reflection coefficient of a patch antenna fed by a

  Figure 1 [14] :

x a l ll 1 2 ll 2 ll9 sx ll1 3 (

 3 B. 5) Z 0 Li Y H; =B Y { D/a 2 -a 9 -a6 + a 1 0 + Z a )(n -l !l) _ _!_ __!___ L'1 l }sy 17 4 µ yy y my (B.6) Z 0 L'1::.H: = B z { D z (-a 1 + a 12 + a 3 -a ll" 1 s } + z a )(n -1 / 2 ) ___ L'1 I''Lit is the time step, s = 2c 0 Lit with c 0 the free space velocity of light. Factors C� and D� are given by:

  Stub

f 5 � = 4 (

 54 admittances (Y5� andZ,�) are given by : �� -1) and Z, � = 4 ( � -1) with s {x, y, z} "'where the normalized resistance R� = -s --"' • Finally, sources er � . r �) corresp MÖ ond Z o µ ;; e"::> m� to a spatial average value [14] at time n .1 t which can be interpreted as value at the cell center. The next step consists of combining the field values whose expression is given by (A. 1) in Appendix A. Finally, one can show [14] that reflected voltages at time (n 112) Lit can be computed as MYllows : (B. 10) -;--) b (n + 11 2) = � b n _ � c a P where � a P = {a 12 • ag, a ll , ag, C l 7, a 1 0• as, Cl 4 , a 2 , a 6, a 3 , a l ' a l 3 ' aw a l5 ' a w a l7' a ls}r

  Appendix C: Definition of additional voltages for the PML medium Incident (dPM L) and reflected z;: M L) voltage vectors added to the scN for the modeling of PML media are given below : ; z + Ct Z 0 Li y H� Li e z -C Z Li H Y x 9 x 0 y 9 .1 y E� z -C y Z 0 L'1 x H; Li £ YZ + C Z .1 H x y 8 y 0 x 8 Za Li x Hj+ D x L'1 z E� Z Li H xy -D .1 E z 0 x 7 x z 7 yx z Z 0 Li y H 6 -D y Li2 E 6 Z 0 Li y H�� +D y Li z EJ 0 Li E' z -C Z Li H Y x 2 x 0 y 2 Li x E� z + C x Z0 .1 y H� Li £ YZ + C Z .1 H x y 4 y 0 x 4 yz C Li H x .1 y E sy Z a x 8 Z Li H'<Y-D .1 E z 0 x 5 x z 5 Z0 Li x Hj + D x Li z E; Z Li H yx + D .1 E z 0 y 6 y z 6 yx Z Z 0 Li y H 1 0 -D y Li ft10

y x � (/ 26 4 Appendix B: Accelerated procedure for (lf) calcu lation [START_REF] Pel | A general formulation of a three-dimensio nal TLM condensed node with the modeling of electric and magnetic losses and current sources[END_REF] First of all, the field components at the cell center are computed at time nAt from (7t) known at time (n -1 \2) A t and adding the source contribution at time n d t: