Scott Delbecq

Marc Budinger

Jérôme Piaton

A MORE INTEGRATED DESIGN APPROACH FOR EMBEDDED MECHATRONIC SYSTEMS: APPLICATION TO ELECTRICAL THRUST REVERSER ACTUATION SYSTEMS

Keywords: Sizing, Multidisciplinary System Design Optimization, Design Space Exploration, Electromechanical Actuator, Thrust Reverser

The design of embedded mechatronic systems like electromechanical actuation systems involves heterogeneous knowledge due to the interference of several engineering specializations and the multiple physical laws that govern their behaviour. This results in costly iterations during the design process and non-optimal solutions. Multidisciplinary System Design Optimization techniques provide theoretical foundations and computational tools for optimizing large and multidisciplinary systems. The approach taken uses these techniques within a proposed design and sizing methodology. It allows a holistic sizing of mechatronic engineering systems with emphasis placed on reusability and rapid decision making. This methodology has been applied in an industrial context to size an electrical thrust reverser actuation system. This more integrated design approach has resulted in significant assessments and insights in early trade-off studies for this complex aircraft subsystem.

Introduction

Environmental protection is a key driver for next generation aircrafts as drastic goals of Flightpath 2050 are asked by ACARE european council [START_REF] Kallas | Flightpath 2050 europe's vision for aviation[END_REF]. In order to reduce drastically fuel burn, emissions and noise, new aircraft concepts have been proposed. They investigate innovative propulsion systems [START_REF] Charier | System of compact contra-rotating propellers[END_REF], propulsion-airframe integration for boundary layer ingestion [START_REF] Gray | Modeling boundary layer ingestion using a coupled aeropropulsive analysis[END_REF], structural alloys, systems and equipment architectures. One research axis consists in electrifying propulsion, systems and equipment by betting on the advantages of electrical technologies in terms of integration, maintenance and power efficiency. More than ever, this more electrical aircraft course of action is driving the research and technology of the aeronautical industry and academics. The research axes that concerns non-propulsive systems can be separated in two main axes. The first axis investigates the concept of bleedless aircrafts which consists in replacing pneumatic systems by electrical systems. The second axis investigates the concept of hydraulicless aircrafts which consists in replacing hydraulic systems by electrical systems. These electrical systems are more eco-friendly than hydraulic systems due to the removal of irritating Skydrol fluid, more efficient and are also easier to integrate and maintain. The hydraulicless aircraft concept has lead to the electrification of actuation systems for flight controls [START_REF] Van Den Bossche | Model Based Design and Experimental Validation of Electro-mechanic Actuator Systems for a Nose Landing Gear[END_REF], landing gears [START_REF] Doberstein | Model Based Design and Experimental Validation of Electromechanic Actuator Systems for a Nose Landing Gear[END_REF] and thrust reversers.

The integration of these new technologies increases risks of unexpected and critical behaviours of systems leading to costly re-designs. These new electrical systems have to prove that they meet the severe reliability and availability requirements of commercial aircraft when operating in a harsh vibratory and thermal environment. Therefore, one of the challenges for developing a more electrical aircraft is to master the design and integration of the forthcoming mechatronic systems like actuation systems. Hence, emphasis must be placed on developing new design approaches and methodologies to design and size these new technologies that would ensure an optimized and successful integration during future aircraft developments. In parallel, Multidisciplinary Design Optimization (MDO) tools have significantly increased in sophistication in the recent decades. Thus, they have lead to many applications such as aircraft conceptual design (0D analytic models), aircraft trajectory optimization (0D-1D lumped parameter models) or high fidelity wing design (3D Finite Element Method models). This paper will focus on a developed designs and sizing methodology and its application to electrical thrust reverser actuation systems (ETRAS). For that, the paper is organized as follows. Section 2 presents the methodology and its associated framework. Section 3 outlines the architecture, design drivers and sizing scenarios of the ETRAS. Section 4 describes the modelling tasks achieved for this case study. Section 5 presents the design optimization/exploration formulation and the results obtained. Conclusions and future work are presented in Section 6.

Sizing methodology and associated framework

The design of mechatronic systems involves several stakeholders which have different type of engineering specialization knowledge. The knowledge is available in different engineering teams of the design office as shown in Fig. 1.

During system design the design drivers of the system are well known by the system architect which enables him to manage the other engineering specializations (electronics, electrical machine design, mechanical design, software...) design tasks. The other engineering specialization experts are well aware of technological limits of their components and therefore assess rapidly their performances.

Experienced system architects are able to make system design trade-offs without inter- Fig. 1 Different engineering teams and knowledge flows acting with other engineering specialization experts. However, for unconventional systems or in the absence of a multi-domain experienced system architects it is not possible. The proposed methodology permits better consideration for engineering specialization knowledge during system sizing. A methodology consists in a process, methods and a framework. This section is organized as so.

Process

A process can be defined as successive steps. The first step consists in the definition of requirements, design drivers and sizing scenarios of the system. Once sizing scenarios are defined, modelling needs are identified and modelling activities can be outlined. The second step is thus the implementation of the models requested for sizing. The third steps assembles elementary models into reusable groups that are manageable by a single individual. These groups correspond to component level sizing models. The fourth step assembles these component sizing models to form a complete system sizing model. In the fifth step, the system sizing model is analyzed and optimized. The final and sixth step explores the design space using the system sizing model. The next part summarizes these steps and outlines the associated methods. Reverser Actuation Systems 2.2 Methods

Step 1: Requirements, design drivers and sizing scenarios definition

The definition of requirements, design drivers and sizing scenarios has to be achieved through exchange of information between stakeholders during early development phase. The information consists in design knowledge represent needs and constraints of each of the stakeholders. These needs and constraints are translated into requirements and design drivers used to design and size the system. An efficient method to represent each of this knowledge and their interaction is through a system breakdown structure. Nodes represent components, their design drivers and the sizing scenarios. Using different shapes of nodes help distinguish them.

Step 2: Elementary computational model generation

In order to provide low computational cost models, the elementary computational model used are algebraic. Algebraic models can be generated using different methods. Analytic models can be obtained by deriving equations representing the laws of physics. Scaling laws are algebraic equations based on reference component [START_REF] Budinger | Estimation models for the preliminary design of electromechanical actuators[END_REF]. Linear regression of component data sheets provide the possibility to obtain algebraic model that compute component geometrical parameters (length, mass) with respect to performances (torque, speed). Surrogate modelling can be used to fit empirical data or simulation data. Surrogate modelling techniques that generate algebraic models using Finite Element Method (FEM) simulations are available [START_REF] Sanchez | Dimensional analysis and surrogate models for the thermal modeling of Multiphysics systems[END_REF]. Surrogate modelling can also be used for representing lumped parameter simulations.

Step 3: Elementary models assembly into reusable component sizing models

The assembly of elementary algebraic models enables to use acausal modelling for increasing reusability of models. The combination of graphbased methods and symbolic computation meth-ods provides techniques to implement an acausal environment. Such methods also permit to detect computation singularities. This step enables to generate flexible inputs/outputs component sizing models composed of elementary algebraic models whilst checking their solvability.

2.2.4

Step 4: Component sizing models assembly into a system sizing model

The assembly of component sizing models to form a total system sizing model is intricate because of interconnections between variables that come from the system layer. The method used to represent the system model is the N2 diagram [START_REF] Kapurch | NASA systems engineering handbook[END_REF]. It is a diagram that represents in the shape of a matrix the functional or physical interfaces between system elements. In addition, a hierarchical decomposition method is used to deal with the complexity of system containing a large number of elements.

2.2.5

Step 5: Analysis and optimization using the system sizing model Analysis and optimization using the system sizing model can be achieved using Multidisciplinary System Design Optimization techniques.

Gradient-based methods are used with symbolic differentiation for elementary models derivatives and Unified Derivatives Equations [START_REF] Martins | Review and unification of methods for computing derivatives of multidisciplinary computational models[END_REF] for total derivatives. This leads to low analysis and optimization times which is very interesting for rapid decision making during the sizing process.

2.2.6

Step 6: Design space exploration using the system sizing model Design space exploration using the system sizing model can be achieved using Multidisciplinary System Design Optimization techniques for multidisciplinary analysis and Design of Experiments. Visualization methods like parallel coordinates and scatter plots provide efficient features for design and sizing purposes.

Framework

Step 3 to 6 are supported by a in-house framework called BOA (Bind your models, Optimize your system, Accelerate your design) [START_REF] Delbecq | A framework for sizing embedded mechatronic systems during preliminary design[END_REF]. To match the process this framework is decomposed in four environments: Block Generation, Sizing Procedure, Design and Exploration. It was developed using Python programming language. The computational core relies on key scientific computing packages like Scipy, Numpy, Sympy, Net-workX, pyDOE and OpenMDAO [START_REF] Gray | OpenM-DAO: An open source framework for multidisciplinary analysis and optimization[END_REF]. The user interface was implemented using PyQt5. Fig. 2 illustrates the user interface of the sizing procedure environment.

Design variable block

Input block

Objective block Causal block Fig. 2 Sizing Procedure environment This framework has been used to achieve the sizing of aileron rotary [START_REF] Sanchez | Surrogate modeling technique for the conceptual and preliminary design of embedded actuation systems and components[END_REF] and linear [START_REF] Delbecq | A framework for the conceptual and preliminary design of embedded mechatronic systems[END_REF] electromechanical actuators and more recently a primary flight control actuation system [START_REF] Delbecq | Optimization of primary flight control actuation system using parametric sizing models of actuators, power electronics and structural analysis[END_REF]. For the first time, it is used to achieve a conceptual sizing of an ETRAS. This example is described in the following sections. The problem is nonexhaustive but includes the main challenges of ETRAS design and sizing.

Architecture, design drivers and sizing scenarios

The ETRAS is supplied by the AC voltage network of the aircraft. It has to deploy and stow the engine nacelle transcowl(s) during landing. The transcowl is part of the thrust reverser system and can have either two symmetrical transcowls (Cduct) or one same transcowl (O-duct). A O-duct transcowl is considered in this case study.

Architecture

The ETRAS includes an electrical power chain and an electromechanical power chain. The electrical power chain main components are an Autotransformer Rectifier Unit (ATRU), an inverter, a braking resistor and a housing. The electromechanical power chain main components are a brushless motor, flexshafts and linear actuators. This architecture is illustrated in Fig. 3. The ATRU includes an autotransformer, two rectifiers and two interphase inductors. An actuator includes a bevel gear, a thrust bearing, a ball screw and its tube (housing). We focus here on the sizing of components that have a significant effect on the total mass of the system. Hence, the sizing of the system includes the autotransformer, the braking resistor, the housing which acts as a heat sink for the braking resistor and inverter, the motor, flexshafts and actuators.

Actuator

Design drivers

The system design drivers considered for this study are the requirements such as fail-safe system, aerodynamic loads and the transcowl mass. The Aborted Take Off (ATO) deploy and stow performances are specified as a maximum time to achieve full stroke. Hence, position/speed mission profile for both sequences are considered as a system design driver.

Design drivers are given for each component. The main design driver of the autotransformer, the braking resistor and the housing is their maximum operating temperature. Since the ETRAS application is very short (typically a few seconds), only heat capacities of components are used to estimate there maximum temperature. The braking resistor is also selected on its peak power. The braking resistor resistance has also to be carefully adapted to DC bus voltage. The main design drivers of brushless motors for dynamic and transient applications are maximum torque, maximum speed, maximum temperature and inertia. For flexshafts their maximum/fatigue torque and maximum speed are technological limits whereas rotational stiffness and inertia parasitic characteristics. Bevel gears are selected on their reduction ratio and maximum/fatigue torque whereas inertia is a parasitic characteristic that has an effect on the sizing problem. Thrust bearings are selected using the maximum/fatigue axial force. The main design drivers of ball screws are their pitch, the maximum/fatigue axial force and axial stiffness.

Sizing scenarios

The first sizing scenario considered is the ATO deploy at cold temperatures (high friction in actuators and flexshafts). The torque and speed mission profile depends on the deploy position/speed profile, inertias, transcowl mass, aerodynamic loads and mechanical efficiency. It has an significant effect on motor maximum torque and speed and the autotransformer peak power. The aerodynamic load are position dependent. The deploy sequence begins with a phase where the motor drives the transcowl and follows with a phase where the motor brakes the transcowl and thus generates energy.

The second sizing scenario considered is the ATO deploy at hot temperatures (low friction in actuators and flexshafts). The high mechanical efficiency leads to more braking torque for the motor, higher peak power ad more energy to dissipate for the braking resistor.

The third sizing scenario taken into account is the stow at cold temperatures. Loads are smaller that the ATO deploy but there is no braking phase so the motor is always driving the transcowl. Hence, it does not determine peak torque or power for motor and autotransformer but generate important energy losses that they have to contain thanks to their intrinsic heat capacity.

The fourth and fifth sizing scenarios taken into account are jamming failure modes. As the system has to be fail-safe, all components shall resist to such event. The loads taken are for ATO deploy. The stress generated in mechanical components depends on inertias, stiffness, transcowl mass and motor torque/speed operating point at jamming. Fig. 4 outlines the architecture considered for the sizing problem, the different design drivers and sizing scenarios as well as their interactions and therefore the high complexity of the ETRAS design.

Sizing models

The sizing models used enable to compute the quantity of interest of sizing scenarios. Estimation models are also used to determine to component parameters involve in sizing scenarios models. Efforts are made to use algebraic models for optimization by introducing 0D models (scaling laws, analytic), surrogate models of 0D-1D (Dymola) and 3D models (FEM).

Sizing scenarios models

A surrogate model of lumped parameter model simulation is used to represent ATO deploy (cold and hot) and stow (cold). A 0D-1D lumped parameter model is achieved using Dymola Software. It computes the maximum power, torque, speed and energy at transcowl level during motor and generator (Deploy only) phases with respect to aerodynamic loads, total inertia of the system and a one degree of freedom trapezoidal speed profile. The trapezoidal speed profile is implemented so that the only degree of freedom is maximum speed at transcowl level. The model uses the inverse simulation features of Dymola as shown in Fig. 5.

A Latin Hypercube Sampling of 1000 sam- ples is achieved on total system reflected inertia and speed profile parameter. The simulations are executed using a co-simulation Functional Mockup Unit. Post-analysis computes the quantity of interest (energy, peak power...) for each simulation and then a third-order polynomial response surface surrogate is used to represent them with respect to total system inertia and speed profile parameter.

As the ETRAS operates during a very short time, components maximum temperature sizing scenario are computed with respect to the equivalent energy of losses and heat capacity.

Maximum speed of motor and flexshafts are checked using analytic models.

Jamming sizing scenarios models compute the peak torque/force applied on mechanical components with respect to inertias, stiffness, aerodynamic loads, motor torque and speed.

They use an analytic model (Equation 1) that assumes that the equivalent kinetic energy of spinning inertias is transformed into potential elastic energy of mechanical stiffness.

F jam = M eq K eq V 2 (1)
Where F jam is the jamming force, V the equivalent mass's (M eq ) speed and K eq the equivalent stiffness at the jamming point.

Estimation models

The autotransformer sizing is achieved using scaling laws and an existing design. The autotransformer is sized for constant magnetic saturation as it operates during a very short time. However, the thermal behaviour has to be checked. Hence, the heat capacity is also estimated using a scaling law. The braking resistor electrical resistance and thermal resistance are obtained using manufacturer analytic model. The housing heat capacity, which acts as a heat sink for the braking resistor and the inverter, is obtained with an analytic model.

The motor electromagnetic model is obtained using surrogate modelling of FEM simulations (Fig. 5). It considers the magnetic saturation of the core material. It computes the electromagnetic torque with respect to three high level variables, the stator diameter, the yoke thickness and current density in the windings. Models that compute Joules and iron losses are also implemented. Details of this model can be found in [START_REF] Sanchez | Dimensional analysis and surrogate models for the thermal modeling of Multiphysics systems[END_REF].

Fig. 6 Brushless motor electromagnetic FEM model using FEMM Flexshafts diameter, stiffness (Equation 2), inertia and mass are estimated using scaling laws with respect to peak torque.

K f s = G -1 l re f ( T f s T f s re f ) 4 3 L -1 f s (2)
Where K f s is the flexshaft stiffness, G l re f and T f s re f respectively the linear torsional deflection and peak torque of a reference flexshaft, T f s the flexshaft peak torque and L f s the length of the flexshaft. The bevel gear is sized using a scaling law that assumes the tooth contact stress to be the main effect on gear mass.

Ball screw parameters are estimated using scaling laws with respect to peak axial load. However, buckling load is check using a scaling law as well.

Thrust bearing scaling law assumes Hertz contact stress as the main design driver.

System sizing model

The system sizing model assembles all the sizing scenarios and parameter estimation models. The mission profile analysis evaluates the ATO deploy (cold and hot) and Stow (cold) with respect their respective speed profile parameter and total system inertia. It analyzes the more constraining variables in the different operating quadrants for the components sizing such as maximum speed, maximum load and energy. These sizing inputs are then distributed to the different component sizing bricks. The system level design variables are two speed profile parameters, spur gears reduction ratio, motor stator diameter, yoke thickness and current density. The system level outputs of the sizing model are autotransformer, braking resistor, housing and motor temperatures. Motor and flexshaft maximum speeds have also to be evaluated.

The sizing problem includes two multidisciplinary couplings. First, the total system inertia is required to compute mission profile analysis but is evaluated using the inertia of components estimated with respect to the outputs of mission profile analysis. The second multidisciplinary coupling comes from the jamming scenarios. The inertias and stiffness of components are required to achieve jamming load analysis but are evaluated using components sizing models that require the value of this jamming load. Hence, these multidisciplinary couplings have to be solved to provide a consistent design solution.

A Multidisciplinary System Design Optimization monolythic formulation is used to solve them. It consists in adding a normalized design variable and an inequality constraint for each coupling and the system total mass as an objective. For example, the first coupling can be described by:

F nom = f (J s , α) J s = g(F nom ) (3)
Where F nom is the nominal force from the mission profile analysis, α the ATO deploy speed profile parameter and J s the system total inertia. Is transformed into:

F nom = k os • h(α) J s = g(F nom ) F nom ≥ f (J s , α) (4)
Where h = f (J s = 0, α).

This formulation has a lighter computational cost compared to other monolythic formulations like Individual Discipline Feasible (IDF) or All At Once (AAO) [START_REF] Martins | Multidisciplinary Design Optimization: A Survey of Architectures[END_REF]. Nevertheless, the use of the inequality is only possible because decreasing the nominal force F nom decreases the total system mass and that is the objective.

Design optimization and exploration

Formulation

The representation of the system sizing model and optimization problem is given in Fig. 7 using the Extended Design Structure Matrix diagram [START_REF] Lambe | Extensions to the design structure matrix for the description of multidisciplinary design, analysis, and optimization processes[END_REF].

The sizing optimization implemented minimizes the total mass of the ETRAS M tot with respect to deploy and stow speed profile parameter α d , α s , motor diameter D mot , yoke thickness e y , current density J cur and consistency variables for mission profile load analysis k os F and jamming scenario k os jam . The constraints come from technological limits of components like maximum speed Ω * and maximum temperature Θ * . Two consistency constraints introduced previously are used to solve the two multidisciplinary couplings. The optimization problem formulation is the following: minimize M tot with respect to α d , α s , N red , D mot , e y , J cur , k os F , k os jam subject to

Ω mot -Ω mot max ≤ 0 Ω f s -Ω f s max ≤ 0 F dyn -F nom ≤ 0 F jam s -F nom mech ≤ 0 F jam bs -F nom mech ≤ 0 Θ at -Θ at max ≤ 0 Θ hous -Θ hous max ≤ 0 Θ br -Θ br max ≤ 0 Θ mot -Θ mot max ≤ 0 (5)
Performing this optimization provides the possibility to assess rapidly integration parameters such as mass and dimensions of the system in order to make system integration trade-offs.

For example such sizing code can assess the effect of changes of requirements like deploy time on the total system mass. In addition, the effect of changes of technologies like motor inertia can be evaluated rapidly.

Results

The study of the effect sizing scenarios on the global design is used to illustrate both the utilization of the methodology and the sizing problem. The sizing scenarios such as jamming and many other have to be included in the design in order to avoid non-compliant design solutions. Results show that these jamming scenarios have a significant impact on the mass of the actuators. As they are the main source of mass, it directly impacts the system total mass. Therefore, emphasis must be placed on implementing sizing codes that enable the consideration of a large number of sizing scenarios, design variables and design constraints.

Discussions

The proposed more integrated design approach uses a in-house framework. This framework provides the possibility to capitalize the knowledge through sizing models libraries and sizing projects libraries. This way sizing scenarios can be reused for other projects. This is a tremendous feature as generations of experts have started to retire and leave organizations with their expertise. The sizing scenarios of the ETRAS case study presented in this paper come from the main design drivers of the system. Nevertheless, the sizing of the ETRAS includes also scenarios that have not been taken into account such as geometric integration or different voltage/power supply configurations. However, this confirms that design sizing methodologies are mandatory to design such complex systems in a holistic and tightly coupled manner.

Conclusion and future work

The more electrical aircraft course of action has lead to the integration of new electrical technologies like electromechanical actuation systems. The introduction of such technologies requires to developed design and sizing methodologies in order to achieve successful developments. Such design and sizing methodology and its associated in-house framework have been outlined. It enables to integrate different engineering specialization knowledge to achieve a holistic sizing of the system. This permits to capitalize expertise knowledge throughout the design process that can be reused easily for forthcoming projects. In addition, the framework includes rapid optimization and exploration capabilities as emphasis is placed on rapid decision-making during early design phases. Then, this methodology was applied to the sizing of an ETRAS where a large number of disciplines and scenarios, various types of models (0D,1D,3D) and multiple couplings were involved as typical MDO applications.

This paper has shown that it is possible to integrate and capitalize different engineering specialization expertise in one same design and sizing environment. The ETRAS had not yet been treated in scientific literature. This paper is an introduction to the complexity of such system as only the main sizing scenarios were considered. However, this offers interesting perspectives for future work. The mission profile has a significant effect on the total system sizing. Hence, future work will include the optimization of both mission profile and motor control strategies.

Fig. 3

 3 Fig. 3 ETRAS power architecture

Fig. 4 Fig. 5

 45 Fig. 4 Non-exhaustive system breakdown graph of the ETRAS sizing problem with component design drivers (oval), system design drivers (octogone) and sizing scenarios (hexagone)

Fig. 8

 8 Fig. 8 Mass breakdown optimization results without (left) and with (right) jamming scenarios

Fig. 7

 7 Fig. 7 XDSM diagram for the ETRAS sizing optimization problem

A more integrated design approach for embedded mechatronic systems: Application to Electrical Thrust Reverser Actuation Systems

y , J (0) cur , k t,(0) osF , k Contact Author Email Address mailto: scott.delbecq@safrangroup.com mailto: marc.budinger@insa-toulouse.fr

Copyright Statement

The authors confirm that they, and/or their company or organization, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and distribution of this paper as part of the ICAS proceedings or as individual off-prints from the proceedings.


Ce Zheng

Malcolm Egan

Laurent Clavier

Gareth W Peters

Jean-Marie Gorce

Copula-Based Interference Models for IoT Wireless Networks

As the Internet of Things (IoT) is largely supported by wireless communication networks in unlicensed bands, there has been a proliferation of technologies that use a large variety of protocols. An ongoing challenge is how these networks can coexist given that they have different power levels, symbol periods, and access protocols. In this paper, we study the statistics of interference due to IoT networks that transmit small amounts of data. A key observation is that sets of active devices change rapidly, which leads to impulsive noise channels. Moreover, these devices operate on multiple partially overlapping resource blocks. As such, we characterize the joint distribution and propose a tractable model based on copulas. Using our copula model, we derive closed-form achievable rates. This provides a basis for resource allocation and network design for coexisting IoT networks.

I. INTRODUCTION

With the increasing scale of wireless network deployments for the Internet of Things (IoT), an ongoing challenge is to ensure that these networks can coexist. A key issue is that interference from a large number of devices, even if they operate at low power levels, can degrade the performance of other communication networks. This means that the interference statistics are difficult to characterize and has lead to a number of experimental studies on the interference in various contexts [START_REF] Lauridsen | Interference measurements in the European 868 MHz ISM band with focus on LoRa and SigFox[END_REF]- [START_REF] Vejlgaard | Interference impact on coverage and capacity for low power wide area IoT networks[END_REF].

One feature observed in IoT networks is the presence of impulsive interference, where high amplitude interference is significantly more likely than in Gaussian models. This behavior has been observed both in experimental studies [START_REF] Vejlgaard | Interference impact on coverage and capacity for low power wide area IoT networks[END_REF] and also in theoretical analysis [START_REF] Pinto | Communication in a Poisson field of interfererspart I: interference distribution and error probability[END_REF], [START_REF] Egan | Wireless communciation in dynamic interference[END_REF]. As a consequence, Gaussian models are often not appropriate and the interference statistics lie in a more general class of models.

Introducing non-Gaussian interference models implies that the interference statistics are not simply characterized by their mean and variance. This issue is amplified in settings where a number of frequency bands are used for transmissions. In these cases, the covariance matrix is not sufficient in order to characterize the joint interference statistics over multiple frequency bands. This is particularly evident when the frequency bands used by different users only partially overlap [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF], such as in non-orthogonal multiple access (NOMA) schemes [START_REF] Dai | Nonorthogonal multiple access for 5G: solutions, challenges, opportunities, and future research trends[END_REF] including sparse code multiple access (SCMA) [START_REF] Nikopour | Sparse code multiple access[END_REF].

Due to the rich dependence structure possible for the joint distribution of non-Gaussian random vectors, a key question is how it should be modeled. For the purposes of analysis, simulation and design, it is highly desirable that the class of distributions captures the realistic dependence between the interference on different bands, and admits tractable estimation and simulation methods.

A popular method in statistics for tractably modeling non-Gaussian multivariate distributions is based on copulas [START_REF] Nelson | An Introduction to Copulas[END_REF].

Copula models have been applied in wireless communications to model fading in MIMO systems [START_REF] Gholizadeh | On the capacity of MIMO correlated Nakagami-m fading channels using copula[END_REF], for blind source separation [START_REF] Kitchen | Copula techniques in wireless communications[END_REF], and modeling interference in ultra wideband or ad hoc networks [START_REF] Yan | Skew-t copula for dependence modelling of impulsive (alpha-stable) interference[END_REF], [START_REF] Soret | Modeling interference with alpha-stable distributions and copulae for receiver design in wireless communications[END_REF].

In the copula modeling approach, the joint distribution function of a random vector in R n , X = [X 1 , . . . , X n ] is given in the form

F (x 1 , . . . , x n ) = C(F 1 (x 1 ), . . . , F n (x n )), (1)
where

C : [0, 1] n → [0, 1]
is called a copula function, and F i , i = 1, . . . , n are the marginal distribution functions. When both the joint and marginal distributions admit density functions, the joint probability density function is then given in the form

p X (x 1 , . . . , x n ) = c(F 1 (x 1 ), . . . , F n (x n )) n i=1 p Xi (x i ). ( 2
)
That is, the joint probability density function decomposes into the product of the marginal densities and another function c : [0, 1] n → R + which captures dependence between the different components of X.

In this paper, we introduce copula models for interference in the context of IoT wireless networks. In particular, we adopt the class of t-copulas with α-stable marginal distributions. This choice of marginals is consistent with theoretical models for interference arising from Poisson spatial fields of interferers [START_REF] Pinto | Communication in a Poisson field of interfererspart I: interference distribution and error probability[END_REF], [START_REF] Ilow | Analytic alpha-stable noise modeling in a poisson field of interfers or scatterers[END_REF]- [START_REF] Gulati | Statistics of co-channel interference in a field of Poisson-Poisson clustered interferers[END_REF] that exhibits heavy-tails and consistent with experimental studies [START_REF] Vejlgaard | Interference impact on coverage and capacity for low power wide area IoT networks[END_REF].

To show that t-copulas are a good fit for the interference dependence structure, we study simulated data sets arising from a Poisson spatial field of interferers and a SCMA scheme. In our model, interferers randomly select frequency bands to transmit on, leading to partial overlaps. Our fits are shown to be close in terms of the Kullback-Leibler (KL) divergence. Since parameters from t-copulas can also be efficiently estimated and the resulting random vectors easily simulated, our proposed interference model satisfies several criteria for a desirable model.

In order to emphasize the tractability of working with the t-copula model, we also study achievable rates of a communication system with this class of interference. In particular, we consider a system where two devices communicate in a point-to-point fashion with interference modeled via the t-copula. This setting arises when the active interferer sets change rapidly; that is, on a symbol-by-symbol basis [START_REF] Pinto | Communication in a Poisson field of interfererspart I: interference distribution and error probability[END_REF]- [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF]. Such a setting is realistic when the point-to-point link exploits a different protocol to that of the interferers; for example on the 868 MHz band, where Zigbee coexists with the wide area networks of SigFox and LoRa. We derive limits on the resulting additive noise channel. In particular, achievable rates are derived, which form a tractable basis for further analysis and design of IoT networks. II. SYSTEM MODEL Consider a network of interfering devices, which is observed by a device at the origin. The locations of the interfering devices form a homogeneous Poisson point process with intensity λ. Each interfering device transmits over a subset of orthogonal frequency bands, B = {1, 2, . . . , K}.

At each time t, interfering devices independently transmit on a band k ∈ B, with probability p > 0 (p is same for each band). The probability p can be interpreted as proportional to the quantity of data each interfering device seeks to transmit. As such, the probability that an interfering device transmits on M out of K bands is K M p M (1-p) K-M . Due to independent thinning, the locations of the subset of transmitting devices on band k are also governed by a Poisson point process. We call Φ the set of interferers, a Poisson point process with intensity λ, and Φ k (t) the set of active interferers at time t on frequency band k, a Poisson point process with intensity λp.

The interference observed at the origin at time t on frequency band k is given by

z k (t) = j∈Φ k (t) r j (t) -η/2 h j,k (t)x j,k (t), (3)
where r j is the distance from device j ∈ Φ k (t) to the origin, η is the path loss exponent, h j,k (t) is a Rayleigh fading coefficient with h j,k (t) ∼ CN (0, 1), and x j,k (t) is the baseband emission. Suppose that in addition to the set of interferers, there is an interference-limited point-to-point commmunication system consisting of a device located at the origin and a transmitter located at a distance d from the origin. Hence at time t, the device at the origin observes the signal on the k-th frequency band given by

y k (t) = d -η/2 g k (t)w k (t) + z k (t), (4)
where g k (t) is the Rayleigh fading coefficent on the k-th band distributed according to g k (t) ∼ CN (0, 1), w k (t) is the baseband emission on band k and z k (t) is the interference in (3) on band k. In the following, we will focus on the interference limited regime and neglect the thermal noise arising at the receiver.

Motivated by the fact that IoT networks transmit noncontinguously small quantities of data and there exist protocols that have large variations in symbol periods, we assume that the interference at time t is independent of the interference at any other time t , corresponding to different symbols of the point-to-point communication system. The main consequence of this assumption is that the interference random vector z(t) = [z 1 (t), . . . , z K (t)] T is independent in t and hence the additive noise channel in ( 4) is memoryless and stationary. We remark that although z(t) is independent in time, in general the elements of z(t) corresponding to different frequency bands will be dependent.

Due to the additive nature of the channel in (4), it can be viewed as a stationary memoryless real-valued 2Kdimensional vector additive noise channel. Note that since the channel is both memoryless and stationary, we drop the time index and write the channel as

Y = AX + Z, (5)
where A is a real diagonal matrix, and Y, X and Z are all random vectors on the measurable space1 (R 2K , B(R 2K )).

The noise vector Z is given by

Z = [Re(z 1 ), Im(z 1 ), . . . , Re(z K ), Im(z K )] T . (6)
III. COPULA-BASED INTERFERENCE MODELS In this section, we develop our copula-based interference model and validate its use on a simulated data set obtained from the Poisson spatial field of interferers detailed in Section II. As α-stable models will play a key role in our analysis, we first recall their properties that will be used in developing our model.

A. α-Stable Distribution Preliminaries

The α-stable random variables have heavy-tailed probability density functions, which have been widely used to model impulsive signals [START_REF] Nikias | Signal processing with alpha-stable distributions and applications[END_REF], [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF]. The probability density function of an α-stable random variable is parameterized by four parameters: the exponent 0 < α ≤ 2; the scale parameter γ ∈ R + ; the skew parameter β ∈ [-1, 1]; and the shift parameter δ ∈ R. As such, a common notation for an α-stable random variable X is X ∼ S α (γ, β, δ). In the case β = δ = 0, X is said to be a symmetric α-stable random variable.

In general, α-stable random variables do not have closedform probability density functions. Instead, they are usually represented by their characteristic function, given by [19, Eq. 1.1.6]

E[e iθX ]= exp -γ α |θ| α (1-iβ(signθ) tan πα 2 )+iδθ , α = 1 exp -γ|θ|(1+iβ 2 π (signθ) log |θ|)+iδθ , α = 1 (7)
It is possible to extend the notion of an α-stable random variable to the multivariate setting. In general d-dimensional symmetric α-stable random vectors are represented via their characteristic function, given by [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF] E[e iθ•X ] = exp -

S d-1 d k=1 θ k s k α Γ(ds) , ( 8
)
where Γ is the unique symmetric measure on the ddimensional unit sphere S d-1 . A particular class of α-stable random vectors are an instance of the sub-Gaussian α-stable random vectors2 , defined as follows.

Definition 1. Any vector X distributed as

X = (A 1/2 G 1 , . . . , A 1/2 G d ), where A ∼ S α/2 ((cos πα/4) 2/α , 1, 0), (9)
and G = [G 1 , . . . , G d ] T ∼ N (0, σ 2 I) is called a sub-Gaussian α-stable random vector in R d with underlying Gaussian vector G.

Sub-Gaussian α-stable random vectors also play an important role in studying complex α-stable random variables; that is, a random variable with α-stable distributed real and imaginary components. In particular, the generalization of symmetric α-stable random variables to the complex case is known as the class of isotropic α-stable random variables, defined as follows.

Definition 2. Let Z 1 , Z 2 be two symmetric α-stable random variables. The complex α-stable random variable

Z = Z 1 + iZ 2 is isotropic if it satisfies the condition e iφ Z (d) = Z for any φ ∈ [0, 2π). ( 10
)
Due to the fact that baseband signals are typically complex, isotropic α-stable random variables will play an important role in the interference characterization.

The following proposition [19, Corollary 2.6.4] highlights the link between isotropic α-stable random variables and sub-Gaussian α-stable random vectors. Proposition 1. Let 0 < α < 2. A complex α-stable random variable Z = Z 1 + iZ 2 is isotropic if and only if there are two independent and identically distributed zero-mean Gaussian random variables G 1 , G 2 with variance σ 2 and a random variable A ∼ S α/2 ((cos(πα/4))

2/α , 1, 0) independent of (G 1 , G 2 ) T such that (Z 1 , Z 2 ) T = A 1/2 (G 1 , G 2 ) T . That is, (Z 1 , Z 2 )
T is a sub-Gaussian α-stable random vector.

We remark that isotropic complex α-stable random variables are closely related to sub-Gaussian random vectors as can be observed from a comparison with Definition 1. Moreover, unlike the isotropic (or circularly symmetric) Gaussian case (α = 2), isotropic α-stable random variables with α < 2 do not have independent real and imaginary components. This dependence arises from the characterization in Proposition 1 through the dependence of the α-stable random variable A in both the real and imaginary components.

B. Characterization of Interference Marginal Distributions

The marginal distributions of the interference are characterized in the following theorem.

Theorem 1. Suppose that h j,k (t)x j,k (t) is an isotropic complex random variable and

E[|Re(h j,k (t)x j,k (t))| 4/η |] < ∞. (11)
Then, z k (t) in (3) converges almost surely to an isotropic 4/ηstable random variable. Moreover, if the fading coefficients h j,k (t) are i.i.d. and baseband emissions x j,k (t) are also i.i.d. then the scale parameters of all real and imaginary components are equal, given by

σ y k (t) = πλpC -1 4 η E[|Re(h j,k (t)x j,k (t))| 4 η ] η 4 , (12)
where

C α = 1-α Γ(2-α) cos(πα/2) , if α = 1 2/π, if α = 1. (13)
Proof. We sketch the proof as follows. By the mapping theorem for homogeneous Poisson point processes, the process {r 2 j } is a one-dimensional Poisson point process with intensity πλp. By an application of the LePage series representation of symmetric α-stable random variables, it follows that z k (t) converges almost surely to

z k (t) = Z r + iZ i , (14)
where Z r , Z i are symmetric 4/η-stable random variables. Finally, z k (t) is isotropic due to the hypothesis that each h j,k (t)x j,k (t) is an isotropic complex random variable. By Definition 2, it then follows that z k (t) is an isotropic complex 4/η-stable random variable.

We remark that the condition in Theorem 1 is satisfied for most common models for fading and baseband emissions. In particular, if h j,k (t) is Rayleigh and x j,k (t) is Gaussian then the condition is satisfied.

C. Copula Interference Model

As shown in (1), copulas provide a parameteric representation of the joint distribution for a random vector. As a consequence of Sklar's theorem [START_REF] Nelson | An Introduction to Copulas[END_REF], any joint distribution function admits a copula C. This copula is also unique if the marginal distributions are continuous, which is the case for the α-stable marginals obtained in Section III-B.

Nevertheless, obtaining the exact copula is a challenging problem, both in its derivation and also for efficient parameter estimation and simulation. As such, we seek a copula that captures the joint dependence in the interference random vector Z in (6). This copula is required to:

1) well-approximate the true joint distribution of the interference random vector Z; 2) and have a closed-form representation or be computationally feasible to numerically evaluate. A good candidate for the copula to approximate Z can be obtained from standard parametric classes of copulas; namely, archimedean and t-copulas. As will be shown by a numerical study in Section III-D, the t-copula is particularly effective.

To this end, we recall the following definition of the t-copula [START_REF] Nelson | An Introduction to Copulas[END_REF]. Let F v be the distribution function of the univariate tdistribution, given by

F v (x) = x -∞ Γ( v+1 2 ) √ vπΓ( v 2 ) 1 + t 2 v -v+1 2 dt, (15)
parameterized by the degree of freedom v ∈ N. Moreover, the joint distribution F v,Σ of a n-dimensional multivariate tdistribution is given by

F v,Σ (x) = x1 -∞ • • • xn -∞ Γ( v+d 2 ) Γ( v 2 ) (πv) d |Σ| 1+ t T Σ -1 t v -v+d 2 dt, (16)
parameterized by the degree of freedom v ∈ N and the n × n correlation matrix Σ. The t copula is then defined as

C t v,Σ (u) = F v,Σ (F -1 v (u 1 ), . . . , F -1 v (u n ))). (17)

D. Numerical Results

We now study the quality of the t-copula model by fitting it to a data set obtained via Monte Carlo simulations. The Monte Carlo simulations implement the system model described in Section II. Recall that the system consists of a network of interfering devices. Each of these devices transmits on each band k = 1, ..., K, with a probability p > 0. The quantity pK can be interpreted as the average quantity of frequency resources used by each device per transmission.

Our goal is to obtain a statistical model of the interference that fits the simulated data well in terms of the KL divergence. In all the numerical results, the KL divergence between the copula model and the simulated data set is obtained using the Information Theoretical Estimators (ITE) toolbox [START_REF] Szabó | Information theoretical estimators toolbox[END_REF].

In the experiments, we compare three models all with αstable marginal distributions motivated by Theorem 1:

1) The t-copula model detailed in Section III-C.

2) The 2-dim. sub-Gaussian model consisting of independent two-dimensional sub-Gaussian random vectors. In this model, the 2K-dimensional random interference vector Z is decomposed into K two-dimensional random vectors (corresponding to the real and imaginary parts of each band). Each two-dimensional random vector is sub-Gaussian (see Definition 1), independent from each of the other K -1 two-dimensional random vectors. This model is exact when interfering devices only transmit on a single band.

3) The 2K sub-Gaussian model consisting of a 2Kdimensional sub-Gaussian random vector. This model corresponds to the scenario where all devices transmit on every band in B, i.e., p = 1 (see [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF]Theorem 3]). In Fig. 1 and2, we plot the KL divergence, denoted as D KL , of the three models with varying p in 2 and 8 bands, respectively. As expected, the 2-dim. sub-Gaussian model has a very low KL divergence for p ≈ 0 and the 2K sub-Gaussian model has a very low KL divergence for p ≈ 1. We also observe that for p > 0.7 (K = 2 bands) and p > 0.4 (K = 8 bands), the t-copula model performs better than or comparably to the other two models. This suggests that a useful model for the interference statistics is to adopt the 2-dim. sub-Gaussian model for small p and the t-copula model for larger p.

We also observe that as the number of bands K varies, the value of p giving the largest divergence between the t-copula model and the 2-dim. sub-Gaussian model (denoted by p * ) varies, illustrated in Table I. In particular, as K increases, p * decreases. This suggests that for larger numbers of bands K, the t-copula model is the best model for a wider range of p. Finally, note that as K increases, the KL divergence at p * very slowly increases. This suggests that using the 2-dim. sub-Gaussian model combined with the t-copula model has an error that is nearly independent of the dimension. We also make the same comparison with the Gaussian model in Table II., which shows that it is a poor model.

IV. INFORMATION THEORETIC LIMITS

In this section, we study achievable rates of a point-topoint communication system described in Section II with interference governed by the copula model in Section III.

A. Optimization Problem for the Capacity

Motivated by the model developed in Section III, the joint distribution of the interference in ( 3) is of the form

F Z (z 1,r , z 1,i , . . . , z K,r , z K,i )=C(F 1,r (z 1,r ),. . ., F K,i (z K,i )), ( 18
)
where C is a copula.

At present, no bounds on the capacity for power-constrained inputs in general α-stable noise channels are known [START_REF] Egan | Wireless communciation in dynamic interference[END_REF], [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF], [START_REF] Freitas | Capacity bounds for additive symmetric alpha-stable noise channels[END_REF]. As such, we follow [START_REF] Egan | Wireless communciation in dynamic interference[END_REF], [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF] and consider the constraints

E[|Re(X)| r ] ≤ c E[|Im(X)| r ] ≤ c |Re(X)| ≤ A |Im(X)| ≤ A, (19)
where 0 < r < α, c > 0 and A > 0. Note that although a power constraint is not explicitly imposed, the constraints on the amplitude imply that the power is finite. Following the standard procedure used to study Gaussian noise channels, we relax the amplitude constraints, which via the noisy channel coding theorem [START_REF] Han | Information-Spectrum Methods in Information Theory[END_REF] leads to the following optimization problem for the capacity

C = sup µ∈P(R 2K ) I(X; Y) subject to E µ [|X| r ] c, (20)
where P(R 2K ) corresponds to the set of probability measures on (R 2K , B(R 2K )).

B. Achievable Rates

The following theorem gives achievable rates for the channel in (5).

Theorem 2. An achievable rate for the channel in (4) subject to the constraints in ( 20) is given by

C ≥ R = 1 α log   1 +   |a| min r c C(r,α) σ Z   α    , (21)
where

C(r, α) = 2 r+1 Γ r+1 2 Γ(-r/α) α √ πΓ(-r/2) . ( 22
)
and |a| min = min{|a ii |, i = 1, . . . , 2K} with a ii the i-th diagonal element of A, σ Z is the scale parameter of each element of Z, given by [START_REF] Kitchen | Copula techniques in wireless communications[END_REF].

Proof. By Theorem 1, the scale parameters of each component of the noise vector Z in (5) are equal. Therefore, there exists a random vector W with α-stable marginals with unit scale parameters and constant σ Z ∈ R + given in [START_REF] Kitchen | Copula techniques in wireless communications[END_REF] such that

Z d = σ Z W. (23)
By [START_REF] Zolotarev | Mellin-stieltjes transforms in probability theory[END_REF], the fractional moment of a symmetric α-stable random variable V with scale parameter σ V is given by

E[|V | r ] = C(r, α)σ r V , (24)
where

C(r, α) = 2 r+1 Γ r+1 2 Γ(-r/α) α √ πΓ(-r/2) . (25)
As such, recalling that A is a diagonal matrix and choosing

X d = σ X A -1 X I
, where X I is a random vector with the same distribution as W,

σ X = |a| min r c C(r, α) (26)
and |a| min = min{|a ii |, i = 1, . . . , 2K} with a ii the i-th diagonal element of A. This ensures the constraint in ( 20) is satisfied.

By the stability property of scalar symmetric α-stable random variables,

Y:=(σ α X + σ α Z ) 1 α V =σV, (27)
where V is a random vector with symmetric α-stable marginals, each with scale parameter 1.

Let X 1 , X I,1 , W 1 be the first elements of X, X I , W, respectively. Then,

I(X; Y) = I(X; AX + Z) ≥ I X 1 ; σ a 1 X 1 σ + Z 1 σ = I X 1 ; σ σ X σ X I,1 + σ Z σ W 1 = log σ σ Z + h σ X σ X I,1 + σ Z σ W 1 -h(W 1 ) = log σ σ Z , (28)
where the inequality follows from the chain rule for mutual information [START_REF] Cover | Elements of Information Theory[END_REF]. The last equality follows from the fact that X I,1 and W 1 are symmetric α-stable random variables with unit scale parameters. Therefore, the scale parameter of [START_REF] Cover | Elements of Information Theory[END_REF] to (28) then yields the desired result.

σ X σ X I,1 + σ Z σ W 1 is 1 σ α X +σ α Z (σ α X + σ α Z ) 1 α = 1. Applying (

C. Numerical Results

Through the scale parameter σ Z , the achievable rates in Theorem 2 implicitly depend on the probability that a given device access a given band, and also the intensity λ of the Poisson point process governing the locations of the interferers. In particular, we observe that when all other model parameters-including α, η and the fractional moments of the fading and base-band emmisions-are constant, the achievable rate is proportional to log 1 + Q λp , where Q is a constant depending on the remaining fixed model parameters.

In Fig. 3, we plot the achievable rate in [START_REF] Freitas | Capacity bounds for additive symmetric alpha-stable noise channels[END_REF] with varying probability p a given user accesses a given band. Observe that the achievable rate decreases as p increases. The main reason for this is that as p increases, the scale parameter of the interference also increases. We remark that in the achievable rate in [START_REF] Freitas | Capacity bounds for additive symmetric alpha-stable noise channels[END_REF], the copula dependence structure of the interference does not explicitly appear. We note however that it does appear in the distribution of the input X. Achievable Rate (nats) Fig. 3. Plot of the lower bound achievable rate and area spectral efficiency V. CONCLUSIONS A basic challenge in communication over the ISM bands is that a number of large-scale IoT technologies have been proposed that have significantly different protocols. These differences lead to non-Gaussian interference and, in systems transmitting over multiple frequency bands, dependence between the interference on each band.

To provide tractable yet accurate models for the interference, we have proposed interference models based on the tcopula and α-stable marginals. This model has been validated on a simulated data set, which exhibits non-trivial dependence structures.

To illustrate the tractability of the copula interference model, we considered a point-to-point communication link in the presence of interference from an IoT system. We derived a closed-form expression for achievable rates. These results provide a basis to study resource allocation and network design, which we intend to explore in future work.

  Zheng, C. and L. Clavier are with Université de Lille 1, IEMN UMR CNRS 8520, IRCICA USR CNRS 3380, Villeneuve d'Ascq, France. M. Egan and J.-M. Gorce are with Univ. Lyon, INSA Lyon, INRIA, CITI, France. L. Clavier is also with IMT Lille Douai, France. G.W. Peters is with the School of Mathematical and Computer Sciences, Heriot-Watt University, UK.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Plot of the KL divergence between the simulated data set and three statistical models: t-copula model; 2-dim. sub-Gaussian model; and 2K sub-Gaussian model. In the figure, K = 2 frequency bands are considered.

TABLE I VALUE

 I OF p GIVING THE LARGEST KL DIVERGENCE BETWEEN THE t-COPULA MODEL AND 2-DIM. SUB-GAUSSIAN MODEL FOR DIFFERENT NUMBERS OF BANDS.

	K	2	3	4	5	6	7	8
	p *	0.692	0.634	0.576	0.519	0.469 0.425	0.388
	D KL	0.136	0.192	0.221	0.247	0.275 0.287	0.290

B(R

2K ) denotes the Borel σ-algebra of R 2K .

There exist also sub-Gaussian stable vectors allowing for more general dependence structure[START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF], but are not necessary for the purposes of this paper.

ACKNOWLEDGEMENTS

This work has been (partly) funded by the French National Agency for Research (ANR) under grant ANR-16-CE25-0001 -ARBURST.


Alice Bates

Alessandro Daducci

Emmanuel Caruyer

Multi-Dimensional Diffusion MRI Sampling Scheme: B-tensor Design and Accurate Signal Reconstruction

b-tensor encoding enables the separation of isotropic and anisotropic tensors. However, little consideration has been given as to how to design a b-tensor encoding sampling scheme. In this work, we propose the first 4D basis for representing the diffusion signal acquired with b-tensor encoding. We study the properties of the diffusion signal in this basis to give recommendations for optimally sampling the space of axisymmetric b-tensors. We show, using simulations, that the proposed sampling scheme enables accurate reconstruction of the diffusion signal by expansion in this basis using a clinically feasible number of samples.

Introduction

Multi-dimensional diffusion MRI (MD-dMRI) is a recent diffusion signal modelling and reconstruction framework where the diffusion signal is derived from the diffusion tensor distribution (DTD). Advanced magnetic gradient modulation schemes, resulting in b-tensor encoding, enable separation of isotropic and anisotropic tensors . Often diffusion tensors are assumed to be axisymmetric, so the DTD is 4D . While algorithms have been proposed for recovering the DTD from b-tensor encoded measurements of the diffusion signal , little consideration has been given to the b-tensor sampling scheme. Currently, b-tensors are empirically chosen, often from uniform sampling of b-tensor parameters , typically resulting in on the order of 1000 samples to recover the full DTD . To our knowledge, no b-tensor sampling scheme has been proposed for optimally sampling the MD-dMRI signal. We present the first 4D basis for representing the MD-dMRI signal. This enables us to study the properties of the signal in this basis to give recommendations for sampling the space of axisymmetric b-tensor, while enabling accurate reconstruction of the MD-dMRI signal.

MD-dMRI Signal Basis and b-tensor Sampling Scheme

is the diffusion signal acquired with b-tensor encoding, where is the DTD, are second order symmetric positive-definite diffusion and b-tensors respectively, and . For a discrete set of diffusion tensor populations, where is the proportion of tensor population and is the number of populations. The axisymmetric diffusion tensor can be parameterised as , where is a rotation operator. The b-tensor can be parametrised in terms of , Expanding , using the above parameterisation, gives a separable expression in , and the other three parameters, :

. This enables to be expanded in a separable and orthogonal 4D basis which is the product of a 1D basis for dimension and a 3D basis for and dimensions. As is a function of negative exponential of , we use an exponential modulated with a Laguerre polynomial, which together form an orthogonal basis over dimension . Similarly, we use the spherical Laguerre basis , a 3D orthogonal basis with an exponential weighting function in the radial direction, for dimensions, leading to expansion:

where is the spherical Laguerre basis, with and . are spherical harmonics of maximum degree , are scale factors and are maximum orders for and bases respectively. Coefficients are defined by inner product . Due to the separability of this basis, we can design a separable transform, with Gauss-Laguerre quadrature used to choose the weights , and sample locations for the and dimensions. We use the scheme to place the samples in the angular dimension. For a band-limited signal, this quadrature enables exact reconstruction and efficient sampling with the number of samples equal to number of coefficients.

(B) = exp(- ( ) exp(- ( ) exp(-2 ( ( (cos )), S S 0 ∑ d=1 N D w d b s D iso ) d b l D iso ) d b l D iso ) d D Δ ) d P 2 β d cos = cos Θ cos + sin Θ sin cos(Φ - ) β d θ d θ d ϕ d (B) S S0 b s , Θ b l Φ (B) S S 0 b s b s 7 8 , Θ, Φ b l ( , , Θ, Φ) = ( ) ( , Θ, Φ) S S 0 b s b l ∑ p=0 P ∑ n=0 N ∑ ℓ=0 L ∑ m=-ℓ ℓ c pnℓm Z p b s B nℓm b l ( , Θ, Φ) = ( ) (Θ, Φ) B nℓm b l X n b l Y m ℓ ( ) = exp( ) ( ), X n b l n! (n+2)! ζ l 3 ‾ ‾ ‾‾‾‾‾ √ -b l 2ζ l L 2 n b l ζ l ( ) = exp( ) ( ) Z p b s 1 ζ s √ -bs 2ζ s L p bs ζ s (Θ, Φ) Y m ℓ L , ζ l ζ s N , P b l b s c pnℓm = ⟨ ( , , Θ, Φ), ( ) ( , Θ, Φ)⟩ c pnℓm S S0 b s b l Z p b s B nℓm b l = ( (j)) ( (i)) ( (j), (i), Θ, Φ) (Θ, Φ) sin ΘdΘdΦ. c pnℓm ∑ j=0 J w j Z p b s ∑ i=0 I w i X n b l ∫ π Θ=0 ∫ 2π Φ=0 S S 0 b s b l Y m ℓ 06/

Methods

The Magic DIAMOND model is used to study maximum orders needed for expansion of the diffusion signal in the proposed basis and to evaluate reconstruction accuracy of the proposed sampling scheme. Magic DIAMOND , an extension of the DIAMOND model for axisymmetric tensor acquisitions, describes the DTD through a mixture of non-central matrix-variate Gamma distributions, where each diffusion compartment is represented through a single distribution with two shape parameters, that set axial and radial heterogeneity of a compartment, and a scale parameter. Compartments are used to represent diffusion molecules hindered and restricted by fascicles, and free water diffusion. We used a range of diffusivities and distribution parameters representative of those in the human brain.

Results and Discussion

As Magic DIAMOND model is linear, it is sufficient to study each compartment individually; the truncation order of the compartment with the highest truncation order is sufficient to represent all combinations of compartments. Fig. 1 shows the spherical harmonic band-limit as a function of . Fig. 2 shows the reconstruction error as a function of increasing truncation order for maximum and of . It was found that were sufficient for reconstruction error smaller than for all compartments, resulting in the sampling scheme (Fig. 3) with 192 samples, the same order as several state-of-the-art sampling schemes in q-space . Fig. 4 shows the mean reconstruction error evaluated at 1080 uniformly placed measurements for 640 realisations of Magic DIAMOND.

Conclusion

We present the first 4D basis and reconstruction algorithm for representing and reconstructing the axisymmetric b-tensor encoded diffusion signal. We study the properties of the signal to inform how many and which b-tensors should be used, thereby providing sampling recommendations. Although preliminary, the high reconstruction accuracy achieved on Magic DIAMOND model is promising for interpolation of b-tensors encoded measurements to enable algorithms for recovering the DTD to achieve higher accuracy and/or reduce the number of measurements.

Figures

  Figures

Fig. 1 :

 1 Fig.1: Maximum absolute reconstruction difference between the reconstructed signal and the ground truth as a function of the spherical harmonic bandlimit and value. required to accurately represent (maximum absolute reconstruction difference is above , as required by MD-dMRI Gaussian diffusion assumption ) the diffusion signal at different -values is shown by the black line.

Fig. 2 :

 2 Fig.2: Maximum absolute reconstruction difference between the reconstructed signal and the ground truth as a function of the truncation order (left) and the truncation order (right). The truncation order should be high enough so that the reconstruction error is smaller than (shown by the black line).

Fig. 3 :

 3 Fig.3: Proposed b-tensor sampling scheme: a) the b-tensors shown for a single orientation, and b) the orientations used in the proposed sampling scheme, projected onto a single sphere, samples on the inner most to outer most shell are shown in green, red (two shells with this sampling) and blue for each shell respectively. Locations where antipodal symmetry is used to infer the value of the signal are lighter in colour.

Fig. 4 :

 4 Fig.4: Mean absolute reconstruction error over 1060 uniformly placed samples on a Cartesian grid for 640 realisations of the magic DIAMOND model.

Acknowledgements

Alice Bates was supported by the Australian Research Council's Discovery Projects funding scheme (Project no. DP170101897).

Alessandro Daducci is supported by the Rita Levi Montalcini, MIUR for the recruitment of young researchers.


M Kachel

J Heymes

L Ammour

M Bautista

G Bertolone

A Dorokhov

S Fieux

F Gensolen

M Goffe

C Hu-Guo

F Lefebvre

F Pain

P Pangaud

L Pinot

P Gisquet

P Lanièce

C Morel

M.-A Verdier

M Winter

L Zimmer

J Baudot

Implantable CMOS Pixel Sensor for Positron Imaging in Rat Brain

Keywords: Semiconductor Detectors, CMOS Pixel Sensors, Beta Probe, Brain Imaging, Neuroimaging

Implantable CMOS pixel sensor for positron imaging in rat brain

Introduction

Neuroimaging results of anaesthetised or highly restrained animals suffers from undesired effects that are limiting their validity [START_REF] Gao | Time to wake up : Studying neu-447 rovascular coupling and brain-wide circuit function in the un-448 anesthetized animal[END_REF][START_REF] Alstrup | Anaesthesia for positron emis-451 sion tomography scanning of animal brains[END_REF]. In contrast, real time monitoring of molecular processes in brains of fully awake animals has a strong potential to improve the robustness of observations that couple neuroimaging and behavioural studies [3,[START_REF] Cherry | Functional whole-brain imaging in behaving rodents[END_REF].

Approaches based on PET imaging [START_REF] Patel | 458 Imaging dopamine release with Positron Emission Tomography 459 (PET) and 11C-raclopride in freely moving animals[END_REF][START_REF] Kyme | 462 Optimised Motion Tracking for Positron Emission Tomography 463 Studies of Brain Function in Awake Rats[END_REF][START_REF] Schulz | Integrating PET with behavioral neuroscience using 466 RatCAP tomography[END_REF] have only partly achieved this goal due to a remaining mechanical support, or an overall weight, or a limited space that does not guarantee the animal freedom of movement for behavioral studies.

Another approach consists in performing invasive measurements using an autonomous positron sensitive probe implanted in specific brain regions by stereotaxic surgery.

An advanced solution using a needle-shaped passive silicon sensor was developed with PIXSIC [START_REF] Märk | A wireless beta-472 microprobe based on pixelated silicon for in vivo brain studies 473 in freely moving rats[END_REF][START_REF] Godart | PIXSIC : A Pixellated Beta-Microprobe 479 for Kinetic Measurements of Radiotracers on Awake and Freely 480 Moving Small Animals[END_REF][START_REF] Balasse | PIXSIC, a Pixelated 484 β+-Sensitive Probe for Radiopharmacological Investigations in 485 Rat Brain : Binding Studies with [18F]MPPF[END_REF][START_REF] Balasse | PIXSIC : 489 A Wireless Intracerebral Radiosensitive Probe in Freely Moving 490 Rats[END_REF]. Figure 1 illustrates the full PIXSIC concept. The miniaturisation of 19 the sensor signal treatment within the PICPUS ASIC [START_REF] Godart | PIXSIC : A Pixellated Beta-Microprobe 479 for Kinetic Measurements of Radiotracers on Awake and Freely 480 Moving Small Animals[END_REF] 20 allows wireless data transmission through a relative light 21 system carried on a backpack. PIXSIC actually releases 22 physical hindrances to the rodent movements. While in-23 vasiveness of this technique is a major point of concern, 24 no impact from the small-size probe direct implantation 25 was observed. However, the PIXSIC setup suffers from two 26 main limitations. The needle probe implanted in the brain is composed of 10 passive pixels, each directly wire-bonded to one channel of the readout chip : PICPUS. This chip embeds amplifiers, filter, digitisation, and memories, and is enclosed in dental cement above the skull. A micro-controller steers and gathers data from the readout chip. It is located on a backpack carried by the rodent during experiments. It also embeds active and passive circuits for the polarisation of the sensor and readout chip, and a battery for full autonomy.

it was not really satisfactory for both the constraints.

The usage of monolithic active pixels sensors (MAPS)

was proposed [START_REF] Ammour | MAPSSIC, a Novel CMOS Intra-Cerebral β+ 492 Probe for Deep Brain Imaging in Awake and Freely-Moving 493 Rat : a Monte-Carlo Study[END_REF] to improve the implanted probe concept.

MAPS can embed an amplification and digitisation stage

in each pixel, basically replacing the PICPUS chip and suppressing a priori sensitivity to electromagnetic environment. Since the amplification allows a small primary signal, the sensitive thickness can be thinner than the overall sensor thickness. Hence a depth of typically a few tens of micrometers generate enough signal for efficient detection and does not need to extent to the full sensor thickness also limiting its sensitivity to γ-rays. The latter can be set optimally with respect to invasiveness and stiffness.

In addition, the probe integration is simplified due to the absence of an additional circuit on the rodent head, as depicted in figure 2.

However, the usage of an implanted active sensor can only be effective under the condition that the local temperature is marginally affected by the probe operation (typically ∆T <1 • C). Consequently, a key design feature is the power dissipated by the sensor.

Following the ideas developed in [START_REF] Ammour | MAPSSIC, a Novel CMOS Intra-Cerebral β+ 492 Probe for Deep Brain Imaging in Awake and Freely-Moving 493 Rat : a Monte-Carlo Study[END_REF] The monolithic active pixel sensor IMIC is a needle shaped probe composed of 16 × 128 pixels of 30 × 50 µm 2 , which provide binary information whether a pixel has been hit or not. IMIC is wire bonded to a PCB acting as an interface to the connector and is enclosed in cement on the head of the rodent. The chip is steered and read out only through digital signals by using a micro-controller located on the backpack together with a battery, which is embedding wireless transmission with the data acquisition PC.

probe were 200 × 500 µm 2 , with no segmentations but along its length. Such a segmentation proved to be granular enough to map subcortical brain activity [START_REF] Godart | PIXSIC : A Pixellated Beta-Microprobe 479 for Kinetic Measurements of Radiotracers on Awake and Freely 480 Moving Small Animals[END_REF] and therefore should be matched quite easily by MAPS, which can feature much higher granularity.

According to first positron count rate observations with the PIXSIC probe [START_REF] Godart | PIXSIC : A Pixellated Beta-Microprobe 479 for Kinetic Measurements of Radiotracers on Awake and Freely 480 Moving Small Animals[END_REF], an average positron flux of about 100 hits/sensor/s is expected over the required dimensions.

This level of positron flux corresponds to at least three orders of magnitude less than standard achievements of CMOS pixel sensors used for charge particle detection.

However, a stringent constraint arises from operation in a living animal brain.

As a matter of fact, active architectures involve power consumption and thus heat dissipation, which could locally overheat brain tissues in contact with the sensor. To avoid behavioural impact or worst, tissue damage, from excess heat dissipation, finite element computations were performed in order to set a limit on the dissipated power. Due to the simplistic nature of the above calculation, the power of 1 mW per sensor is taken as an indicative limit.

139

The actual guideline chosen for the sensor design follows The produced sensor layout is shown on the top of figure 3.

Pixel design

To satisfy simultaneously the main requirements on the pixel functionalities and performances (digital output, low data-throughput, low detection threshold and low power), a modified version of the ALPIDE (ALice PIxel DEtector) [START_REF] Mager | the Monolithic Active Pixel Sensor for the 503 ALICE ITS upgrade[END_REF] circuit developed for the ALICE experiment is implemented. The initial compact pixel architecture amplifies, shapes (with clipping to prevent saturation), and discriminates the signal collected by the diode using only 12 transistors, as depicted on the bottom of figure 3).

Pixels integrate any incoming charge between two consecutive readouts. Since the binary output of the front-end indicating the presence of a hit is only available for a few tens of µs, a 1-bit memory was added per pixel to keep this information until the pixel readout time. The memory is reset after the pixel is read out, so that if a second particle hit the same pixel before readout, the 1-bit memory cannot account for it. The matrix readout strategy is detailed in the next sub-section. Pixels integrate particle hits during the time t wait between two readout cycles, plus the readout time itself :

t integration = t readout + t wait [s]. (2)
For each readout cycle, data are provided to the microcontroller as a vector of 2048 bits to be sent further away through a wireless protocol. In view of limiting the data rate defined by equation 3 and hence the overall power consumption of the MAPSSIC setup, it is better to maximise the waiting time t wait during two readout cycles. However, avoiding pile-up of hits during the same integration period requires to minimise t wait .

Data rate = 2048 t integration [bits/s] (

The foreseen average hit-rate of about 100 hits/matrix/s in our experiments, leads to an optimisation for the integration time of around 50 millisecond. In such conditions the data rate stays at a few tens of kilo-bits per second while keeping the pile-up probability at the per mil level. The main characteristics of the IMIC sensor are sum-271 marised in table 1.

272 DACs driven by the SPI protocol, the bonded chips were confirmed to be fully functional. All results reported below were obtained with these sensors at room temperature and under relative obscurity.

The measured dissipated power by the whole sensor is 160 µW, close to the value expected from simulations in table 1 considering the low absolute value for a complete chip.

Three types of sources have been used to characterise the IMIC sensor. Single frames from each source type are displayed in figures 5a-5d. Soft X-rays (5.9 and 6.4 keV rays) from the 55 Fe source (figure 5a) generate small clusters of mainly two pixels, while the β particles, either from β-( 90 Sr : figure 5b) or from β+ ( 18 F : figure 5c) yield clusters about ten times larger in average. This difference comes from the relatively long ionising path of β-rays in silicon. In contrast, X-rays undergo point-like conversion through photoelectric effect.

Images from 511 keV γ-rays were obtained by placing the 18 F source far away (see sub-section 4.3) from the sensor. The cluster size observed for these γ-rays is similar to the one for positrons (see figure 5d). This indicates that the cluster size information cannot be used to identify (and suppress) the unwanted γ-rays hits on an individual basis information.

To demonstrate that the sensor internal readout performs as expected, an image of a metallic washer was acquired using the 55 Fe source (as shown in figure 5e). The ring is approximately 500 µm thick with an external diameter of 5 mm and an opening diameter of 2 mm. The ring structure can easily be recognized on the image as it completely stops the soft X-rays.

Validation of long integration times

As explained in section 3, the modest awaited activities of the radiotracers allow running the IMIC sensor with a long integration time in the order of 100 ms, which in turn reduces the data flow and simplifies the whole system. This readout strategy is made possible by the in-pixel memory and was tested by measuring a range of activities at integration times varying from 20 ms to 1 s.

For this test, the IMIC sensor was indirectly illuminated with a 90 Sr β source through an aluminium shield. The activity was modulated using three different aluminium thicknesses : 1, 3 and 7 mm. Additional measurements were performed without the source in order to monitor dark count rate. For each combinations of shield and integration time, 1000 frames were acquired with the same detection threshold. Results are presented in figure 6. The mean entries per integration time were converted into a mean activity measured over the whole matrix :

A(t int ) = entries t int [hits/matrix/s]. (4)
The saturation limit corresponds to the situation where 333 every pixel over the matrix has been fired. This limit de-334 pends on the integration time through the simple relation : The analysis of the dark count map revealed that a few 337 pixels (denominated noisy pixels in the following) domi-338 nate the production of fake hits, which explains the dark 339 count rate evolution with the integration time. As there 340 are one to two noisy pixels systematically firing during 341 one integration period, it follows that the dark count rate 342 simply decreases as 1/t integration along with equation 4.

343

The low level for the dark count rate means that the IMIC 344 sensor can measure robustly activities down to approxima-345 tely 10 hits/matrix /s.

346

With the 90 Sr source attenuated by only 1 mm of alumi-347 num, the mean activity follows the saturation limit curve.

348

At t int = 10 ms, the mean activity seen by the sensor is 349 approximately 60000 hits /matrix/s.

350

Using thicker shielding, the mean activity on the sen-351 sor remains below the saturation limit and is constant for

FirstFigure 1 :

 1 Figure1: Schematic of the PIXSIC setup. The needle probe implanted in the brain is composed of 10 passive pixels, each directly wire-bonded to one channel of the readout chip : PICPUS. This chip embeds amplifiers, filter, digitisation, and memories, and is enclosed in dental cement above the skull. A micro-controller steers and gathers data from the readout chip. It is located on a backpack carried by the rodent during experiments. It also embeds active and passive circuits for the polarisation of the sensor and readout chip, and a battery for full autonomy.

  along with the requirement to minimize the power consumption, a first MAPS based on the CMOS pixel sensors technology has been developed. The chip is named IMIC (Imageur Moléculaire Intra-Cérébral ) and is part of the MAPSSIC (MAPS Sonde Intra-Crânienne) project to develop and exploit the first MAPS-based intracerebral positron probe. This paper first describes in section 2 the requirements for the active sensor to be implanted in the brain. The design and operation features of IMIC, optimized for low count rate and power dissipation, are detailed in section 3. 76 Section 4 provides first characterization results, focussing 77 on functional validation and sensitivity with laboratory 78 sources.

79 2 .

 2 Sensor Requirements 80 Neuroimaging exploits a variety of radioactive isotopes, 81 typically 18 F, 15 O and 11 O, which emit β+ at energies 82 from 0 to 1.7 MeV. The probe ability to measure activi-83 ties depends obviously on what energy range the sensor 84 is sensitive to. Consequently an important figure of merit 85 is the detection threshold, which should be as low as pos-86 sible (this requirement is an outcome of the study in [12]) 87 and depends on the signal-over-noise ratio of the sensing 88 element. The threshold value will be discussed in the next 89 section 3 in view of the sensitive thickness and noise level 90 of the designed sensors. 91 As previously mentionned, the usage of MAPS for this 92 kind of probe intends to both enhance the overall robust-93 ness and simplify the system beyond the implanted part of 94 the probe. Thus, we have included into the requirements 95 the digital conversion of the signal inside the pixels hence 96 facilitating further data transmission. It is also expected 97 that the readout architecture of the sensor optimises data 98 throughput to the micro-controller embedded on the ani-99 mal backpack with respect to the measured activity.

Figure 2 :

 2 Figure2: Schematic of the MAPSSIC setup. The monolithic active pixel sensor IMIC is a needle shaped probe composed of 16 × 128 pixels of 30 × 50 µm 2 , which provide binary information whether a pixel has been hit or not. IMIC is wire bonded to a PCB acting as an interface to the connector and is enclosed in cement on the head of the rodent. The chip is steered and read out only through digital signals by using a micro-controller located on the backpack together with a battery, which is embedding wireless transmission with the data acquisition PC.

  The simulated geometry includes a volume of silicon representing the implanted part of the probe (10 mm × 500 µm × 600 µm) in a spheric medium representing the brain at a constant temperature of 37 • C. The highest temperature gradient at the direct probe contact stays below the required 1 • C for a dissipated power below 4 mW. Considering the possible implantation of multiple probes composed of 2 to 4 sensors, we have set a limit of 1 mW dissipated per sensor.

140 1 .

 1 the requirement to achieve the lowest possible power dis-Technology choice and overall geometry144We have chosen to design the sensor in the same ima-145 ger CMOS process with 180 nm feature size as used for 146 various charged particle detection applications[START_REF] Baudot | Optimisation of CMOS pixel sensors for high 498 performance vertexing and tracking[END_REF]. The 147 sensitive volume consist of a 18 µm thick low-doped epi-148 taxial layer, which has demonstrated 100 % detection ef-149 ficiency for minimum ionising particles for various pixel 150 sizes below 50 µm [13]. Such realization guarantee a priori 151 excellent detection performances for β+ within the energy 152 range of interest. 153 In addition, the thin detection layer presents two key ad-154 vantages. First, it mitigates the sensitivity to the anni-155 hilation photons as recommended in [12]. Secondly, it al-156 lows thinning the sensor down to thicknesses from 100 to 157 300 µm, thus achieving an optimum between invasiveness 158 and mechanical stiffness.

159 160Figure 3 :

 1593 Figure 3: Top : layout of the 12000 µm long and 610 µm wide IMIC sensor. Bottom left : magnified layout of a single pixel underlying the location of its main functional components. Bottom right : Pixel circuit schematic derived from the ALPIDE sensor [14].

199 200

 199 Simulations of the pixel architecture using Cadence201 Virtuoso Analog Design Environment have evaluated the 202 equivalent noise charge (ENC) at the collection node be-203 low 10 electrons, hence allowing to set discriminators for 204 the digitisation at an equivalent level of ≈ 100 e -. The 205 corresponding energy threshold for the detection equals 206 ≈ 360 eV, well below the maximum probability value of 207 energy loss in 18 µm of silicon by β+ with energies within simulations predict that the overall front-211 end dissipates about 55 nW per pixel (or 113 µW for the 212 whole matrix), well below the upper limit. The schematic 213 of the pixel and its layout are shown in figure 3.214 215

Figure 4

 4 Figure 4 depicts the principle of operation of the IMIC readout.

3. 4 .

 4 A fully programmable needle shaped sensor In order to minimize the complexity of the external circuitry, IMIC limits the number of required connections by using a command decoder controlled through Serial Peripheral Interface (SPI) protocol. This controller also allows the parametrisation of on-chip digital-to-analog converters (DACs) used for the polarisation of the front-end. The steering signals (and the data) are sent independently between the micro-controller located on the rodent backpack and the sensor. The micro-controller connects to the exter-258 nal world through a wireless protocol. The overall system 259 description will be the subject of a subsequent publica-260 tion. 261 The digital block driving IMIC is located on the sensor 262 edge near the few connection pads located outside of the 263 brain, opposite to the pixel matrix deeply implanted in the 264 tissues. This block adds a power consumption of approxi-265 mately 2 µW. The simulated total IMIC power consump-266 tion amounts to approximately 115 µW. This dissipation 267 does not depend on the activity measured but is entirely 268 driven by the number of pixels populating the matrix.269 270

Figure 5 :

 5 Figure 5: Images acquired with the IMIC sensor. Single frame aquired with a 55 Fe soft-X-ray source (5a), a 90 Sr β source (5b), and a 18 F positron source (5c). A single annihilation photon from the 18 F positron source hit (5d). Image composed of 1000 accumulated frames of a metallic spacer acquired with 55 Fe (5e).A few IMIC sensors were bonded to a printed circuit

Figure 6 :

 6 Figure6: Mean activity seen by the sensor for integration times between 20 ms and 1 s of the IMIC sensor dark counts, and with various activities of 90 Sr β-rays modulated by using various aluminum shielding thicknesses.

352

  every integration time. Thus, in the ballpark of expected 353 average activities (100 hits/matrix/s), integration time up 354 to 1 s can be chosen without information loss and with a 355 safety margin against activity fluctuation of one order of 356 magnitude. This read-out condition allows low data rates, 357 and limits the power consumption of the overall system.

358 4 . 3 .

 43 Preliminary sensitivity study 359 Measurements using 18 F radioisotope, produced at the 360 CYRCé cyclotron at IPHC, in a non-sealed aqueous so-361 lution were performed in view of studying the differential 362 sensitivity to prompt β+ and annihilation photons.

363A 378 379TimeFigure 8 :

 3788 Figure 8: Left : Evolution with time of the pixel count per 100 ms frame duration for the three source configurations described in the text. Right : Average pixel count per frame over five minutes data acquisition for the same configurations.

  is γ when only the distant tank was used.The counting rates for the three source configurations 389 are displayed in figure8b. Within uncertainties, the ave-390 rage counting rate with both sources ( 18 F+γ) or only with 391 the closer tank ( 18 F) are equivalent, they amount respec-392 tively to 56.0 ± 0.7 and 57.9 ± 0.8 pixels/frame. This 393 observation is coherent with the very small counting rate 394 (0.13 ± 0.04 pixels/frame) measured when only the far-395 away source (γ) is used.396The experiment illustrates qualitatively the very small re-397 lative sensitivity of the thin IMIC sensor to γ-rays with 398 respect to β+.

Following [ 12 ]

 12 , we have designed a monolithic active 401 pixel sensor to equip the next generation of intracerebral 402 positron probes for deep brain imaging adapted to fully 403 conscious and freely moving rodents. The IMIC chip fea-404 tures 2048 active pixels over a sensitive area of 480 × 6400 µm 2 . 405 The in-pixel amplification and 1-hit binary memory allows 406 both a low detection threshold (400 eV for the energy de-407 posited) and fully digital matrix readout. This front-end 408 architecture benefits from an extremely low power, addres-409 sing the main challenge of heat dissipation within the im-410 plantation tissues. Finally, being fabricated on a 18 µm 411 thick sensitive layer, the sensor is mostly transparent to 412 γ-rays. 413 The first IMIC sample characterisations confirmed the 414 successful operation of the sensor as a β+ activity loca-415 tor. The measured power dissipation per sensor reached 416 160 µW. This performance sets the current state-of-the-417 art for a probe built on the CMOS pixel technology, and 418 should allow the implantation of a few probes in a single 419 animal with acceptable overheating of the local brain tis-420 sues. Foreseen in-vivo experiments will be conducted to 421 investigate the full compliance of the IMIC sensor with 422 the application requirements.

A

  second generation sensor can benefit from a fully 425 depleted thin sensitive layer, which would allow to increase 426 the pixel size, hence decrease the total number of pixels 427 and leading to a reduction of the total dissipated power. 428 Currently, pairs of IMIC sensors are assembled back-to-429 back to form the implantable part of the final MAPSSIC 430 probe with a total volume of 560 × 610 × 12000 µm 3 . The 431 rest of the construction process for the whole probe and 432 its tests in realistic environment will be the subject of a 433 forthcoming publication.

  507

  Readout strategy for IMIC. When a pixel detects a hit, the event is memorised until rolling-shutter readout. A frame is composed of an adjustable waiting time and the readout time. During the readout, the data of each pixel is sent outside the sensor, bit after bit. The data of the memory is reset after readout to allow the detection and memorisation of a new hit.

		2000	Hit
			Memorization
		1800	Hit Readout
			Waiting Time
		1600	Readout Time Rolling shutter Readout
			Frame
		1400
		1200
	Pixel	1000
		800
		600
		400
		200
		0 0	0.5	1	1.5	2	2.5	3	3.5	4
						Period
	Figure 4:
										216	3.3. Readout strategy
										217	IMIC readout mechanism is kept simple for the sake
										218	of space and power dissipation, which does not depend on
										219	the hit rate, and relies on the rolling-shutter principle. A
										220	token is propagated through the matrix in a shift register
										221	at a fixed clock rate to select one line. The 16-bits word
										222	formed by all the pixel memories in the line is serialised
										223	towards the output of the sensor at a higher clock rate
										224	(at least 16 times), which defines f clock . This operation is
										225	repeated 128 times to gather data from the entire pixel
										226	array, thus the time t readout needed for this full readout
										227	cycle amounts to
										t readout =	2048 f clock	[s].	(1)

Table 1 :

 1 Summary of the simulated characteristics of the IMIC sensor.

	Feature	Value
	Sensor size [µm 3 ]	610 × 12000 × 280
	Pixel pitch [µm 2 ]	30 × 50
	Sensitive depth [µm]	18
	Det. threshold [eV deposited]	≈ 400
	Pixel power consumption [nW/pixel] ≈ 55
	Matrix size	16 × 128
	Sensitive area [µm 2 ]	480 × 6400
	Total power consumption [µW]	≈ 115
	4. IMIC characterisation

442

The authors thank the CYRCé team at IPHC for produ-443 cing and making available the 18 F solutions.

444