Christelle Peyron

Claire Rampon
email: claire.rampon@univ-tlse3.fr

New neurons never sleep (ou) Young neurons tickle memory during REM sleep

Memory formation is a dynamic process and sleep is part of it. Consolidation of memories relies on finely orchestrated brain activities occurring during the post-learning sleep period. In this issue of Neuron, Kumar and colleagues provide evidence that the activity of adult-born hippocampal neurons, during REM sleep, is critical for the consolidation of episodic memory.

Sleep is thought to be critical for the consolidation of encoded events, a process wherein temporary and labile encoded representations form stable and long-lasting memories. The sleeping brain

generates oscillatory patterns of activity that allow information processing. In mammals, sleep occurs in cycles of rapid-eye-movement (REM) sleep and non-REM sleep, two states that are highly different based their electrophysiologic characteristics, their neurobiology and their functions. Non-REM sleep is characterized by high amplitude slow-waves (0-4Hz) rhythm and the presence of sharp-wave-ripples while an oscillatory activity of 6-8 Hz called theta rhythm is predominant during REM sleep. Theta oscillation is observed in local field potential recordings from cortical structures including the hippocampus. Accumulating evidence indicate that REM sleep plays a key role in memory consolidation, and reducing hippocampal theta rhythms during REM sleep results in the impairment of hippocampo-dependent memories [START_REF] Boyce | Causal evidence for the role of REM sleep theta rhythm in contextual memory consolidation[END_REF]. It indicates that theta oscillation may provide a favorable state for shaping memory. However, the underlying mechanisms and cell populations involved remain to be discovered.

In the 60's, [START_REF] Altman | Autoradiographic Investigation of the Distribution and Utilization of Intraventricularly Injected Adenine-3h, Uracil-3h and Thymidine-3h in the Brains of Cats[END_REF] revealed for the first time the existence of newly born neurons in the adult brain, a unique form of brain plasticity. Since this pioneering work, a wealth of evidence has accumulated indicating that neurons generated in the dentate gyrus throughout adulthood are causally involved in several hippocampal functions and that they support the expression of elaborated cognitive functions including the formation of episodic memories. Between 4 and 7 weeks of age in mouse, these adult-born neurons go through a critical period characterized by a low activation threshold and a high excitability, allowing them to respond to even weak inputs [START_REF] Mongiat | Unique processing during a period of high excitation/inhibition balance in adult-born neurons[END_REF]. Owing to these particular characteristics, the maturing new neurons are expected to play a critical role in sensory and memory processes.

In this issue of Neuron, Kumar and colleagues (2020) deployed a sophisticated combination of transgenic and imaging strategies to overcome the scarcity and sparse activity of new neurons, with the aim of investigating their contribution to memory consolidation during sleep. They used calcium neural imaging in freely moving mice equipped with a miniaturized microendoscope to examine new neurons' activity. By crossing two transgenic lines of mice, they generated offspring in which they managed to record the activity of young adult-born granule neurons no older than 4 weeks, i.e. their critical time period, by probing calcium events. Consistent with a previous report [START_REF] Danielson | Distinct Contribution of Adult-Born Hippocampal Granule Cells to Context Encoding[END_REF], Kumar and colleagues first observed that adult-born granule neurons of untrained mice exhibit a sparse activity during wakefulness which further slows down during sleep periods. Then, the authors subjected mice to contextual fear conditioning, a protocol eliciting a simple form of learning, in order to explore the activity of adult-born hippocampal neurons during memory consolidation after learning. Within this task, mice have to learn the association between a new environmental context and the delivery of electric foot shocks, this learning involving the hippocampus. When exposed to this context again, the animals display freezing, a typical fear behavior, routinely used as a proxy to score memory.

A single trial of aversive conditioning is sufficient for learning, which conveniently allows to distinguish the acquisition from the consolidation phase of memory. Moreover, this procedure generates a robust and long-lasting memory whose kinetic is well described [START_REF] Izquierdo | Fear Memory[END_REF]. The authors report that the activity of adult-born hippocampal neurons is significantly reduced during the periods of REM sleep, in mice that later expressed a strong fear memory in comparison to control animals with a weak memory. Interestingly, this effect is only observed during REM sleep, which suggests a statedependent phenomenon. These findings raise a number of questions including why does the very sparse activity of adult-born neurons during REM sleep matter so much to memory consolidation ?

What are the underlying biological mechanisms ?

In search of a comprehensive and dynamic picture of adult-born neuron's activity across the entire memory process, the authors took up the challenge of monitoring the activity of individual adult-born neurons throughout the learning, consolidation and retrieval phases of associative fear memory.

Interestingly, a majority of the adult-born neurons that were active during learning also participated to subsequent activity recorded during REM sleep episodes of memory consolidation. How were these neurons selected during the learning task and how were they subsequently identified to be reactivated during memory consolidation, has yet to be deciphered. Work by [START_REF] Trouche | Recruitment of adult-generated neurons into functional hippocampal networks contributes to updating and strengthening of spatial memory[END_REF] proposed that the new neurons active during the learning phase could be tagged to promote their preferential recruitment when a similar experience is encountered, suggesting reactivation during retrieval.

Further work shall determine whether the reactivation of the same neuronal assembly during REM sleep consolidation relies on such neuronal tagging and what the nature of this tag is. However, in Kumar et al, young neurons remained inactive during memory retrieval. Given that the adult-born neurons contribute to learning and consolidation but not to memory retrieval, it is tempting to speculate that they do not conserve memory, but rather participate in building a downstream neuronal network in charge of storing memory.

To establish the causal contribution of new neuron activity during REM sleep to memory consolidation, Kumar and colleagues used optogenetic tools to silence adult-born neurons during REM sleep. They found that theta oscillation at the surface of the skull above the hippocampus or, directly in the dentate gyrus was not modified. New neurons are few and their activity is sparse. The authors estimated that they managed to silence 467 events over the 6hrs of experimental recording. Such activity is relatively reduced compared to the global activity that can be recorded in the dentate gyrus during REM sleep, possibly explaining the lack of effects on theta oscillations. The monitoring of calcium events confirmed that the neurons were activated but it provided no information regarding the actual pattern of neuronal activity, due to a limited temporal resolution and sensitivity of the technology. The timing of neuronal activity, with respect to local theta oscillation, is critical for the modulation of synaptic strength [START_REF] Orr | Hippocampal synaptic plasticity is modulated by theta rhythm in the fascia dentata of adult and aged freely behaving rats[END_REF] but the firing mode of the new neurons is not yet known. Do they fire at the peak, the dip, the ascendant or descendant phase of theta oscillation? The timing of their firing relative to the theta oscillation is probably critical to their effectiveness in promoting memory retention. Challenging experiments will thus be required in the near future to answer these questions and further understand why young neurons and theta rhythm are conducive to memory stabilization.

Memory formation involves synaptic changes both on mature and immature hippocampal neurons.

For instance, learning of contextual fear accelerates the formation of dendritic spines on adult-born neurons, while memory consolidation strengthens their spine connectivity [START_REF] Petsophonsakul | Memory formation orchestrates the wiring of adult-born hippocampal neurons into brain circuits[END_REF]. Here, Kumar et al show that although sparse, the activity of new neurons during REM sleep controls the expression of genes typically involved in synaptic plasticity. The authors therefore investigated whether the new neurons active during REM sleep underwent structural changes.

Silencing adult-born neurons during REM sleep-dependent memory consolidation had no impact on the increase in their spine density triggered by learning but led to the production of elongated spine necks. It was reported that upon synaptic stimulation, long-necked spines can rapidly be mobilized and their neck shorten to increase synaptic strength [START_REF] Araya | Activity-dependent dendritic spine neck changes are correlated with synaptic strength[END_REF]. This mechanism may be triggered during REM sleep to specifically strengthen and stabilize synapses of the immature new neurons that are recruited upon learning.

Last, the authors provide the ultimate argument that adult-born neurons, of a critical age and with an activity finely-tuned to REM sleep contribute to memory consolidation. Indeed, adult-born neurons were silenced at any time but REM sleep, at different maturation stages (2, 4 and 10 weeks of age) or activated at random frequency during the REM sleep post-training period. It was only when the silencing or the random activation of young adult-born neurons (4 weeks-old) occurred during REM sleep, that fear memory was impaired. These data demonstrate that the contribution of adult-born neurons to memory consolidation happens only when the brain is in a specific state (REM sleep), their activity is tune to theta oscillations and while they display low activation threshold and high excitability. Future experiments should unravel how REM sleep provides an appropriate context for adult-born neurons to deliver information allowing the embedding of a memory trace in downstream circuits although their activity, as currently recorded by calcium imaging, is sparse.

Acknowledgements: This work is supported by the French Centre national de la recherche scientifique.

The authors have no conflict to declare.


Rémi Dorget

Thierry Lubin

Sabrina Ayat

Jean Lévêque

3-D Semi-analytical model of a superconducting axial flux modulation machine

Keywords: Axial flux machine, high-temperature superconductivity, Semi-analytical model

This article presents a 3-D semi-analytical model of a partially superconducting axial flux modulation machine for the computation of the flux density distribution, the torque and the Eddy-current losses in the armature. The modulation of the superconducting coil magnetic field by the superconducting bulks is calculated by solving the Maxwell's equations in scalar potential with the variable separation method. Using Fourier series analysis and the convolution theorem, the bulk region is modelled by a domain with a spatially variable magnetic permeability. This variable permeability enables considering a diamagnetic bulk as well as a more complex magnetic behavior using the Bean model. The results of the proposed model are compared with 3-D finite elements models and good agreement is obtained for a significantly reduced computation time. The proposed model appears to be a suitable tool for optimisation and sizing processes.

I. INTRODUCTION

The use of high temperature superconductors (HTS) in electrical machines is seen as a potential solution to enable aircraft electric propulsion [START_REF] Boll | A holistic system approach for short range passenger aircraft with cryogenic propulsion system[END_REF]- [START_REF] Berg | HTS Electrical System for a Distributed Propulsion Aircraft[END_REF]. To this end, HTS machines must reach a high specific power [START_REF] Engineering | Commercial Aircraft Propulsion and Energy Systems Research: Reducing Global Carbon Emissions[END_REF], recent studies have shown that power to mass ratios over 20 kW/kg are achievable for MW-class generators and motors using liquid hydrogen as coolant [START_REF] Filipenko | Concept design of a high power superconducting generator for future hybrid-electric aircraft[END_REF]- [START_REF] Corduan | Design of superconducting AC machines for hybridelectric aircraft[END_REF]. Among the technologies considered, the flux modulation machine is a promising candidate as shown by several projects and articles on its radial flux form [START_REF] Masson | Scaling Up of HTS Motor Based on Trapped Flux and Flux Concentration for Large Aircraft Propulsion[END_REF]- [START_REF] Malé | Analytical calculation of the flux density distribution in a superconducting reluctance machine with HTS bulks rotor[END_REF] and axial flux form [START_REF] Masson | Design of HTS Axial Flux Motor for Aircraft Propulsion[END_REF], [START_REF] Colle | Design of a superconducting machine and its cooling system for an aeronautics application[END_REF] including a recently built 50 kW prototype [START_REF] Colle | Test of a Flux Modulation Superconducting Machine for Aircraft[END_REF].

This last configuration, presented in Fig. 1, is composed of a static HTS coil fed with DC current generating a magnetic field in the machine. The rotor consists of several HTS bulks (one per pole pair) used in the superconducting mixed state. When the HTS bulks are subjected to the coil magnetic field, an electrical current is induced in the material to screen the magnetic field and the flux lines will concentrate in the apertures between the bulks. In the air gap, this flux concentration will generate a spatially variable magnetic field which generates an electromotive force in the two three-phases conventionnal armatures located on either side of the rotor when it rotates.

Due to the large magnetic flux that can be generated with HTS, no ferromagnetic material is required to reach a consequent magnetic loading, thus making this air-cored topology very light. The HTS coil is a static component, hence the machine magnetisation can be controlled from the coil without any slip rings. The absence of slip ring is a major advantage of this structure for reliability reason.

In order to properly exploit the potential of the flux modulation topology, an optimisation of the design must be performed and requires a precise and fast model. Numerical modelling using the finite element method (FEM) might be the most precise solution. However, the high computational cost of numerical method limits the number of parameter explored during the sizing process. Therefore, analytical and semianalytical modelling using the separation of variables method and Fourier series analysis can be used for the electromagnetic design of electrical machines [START_REF] Corduan | Topology Comparison of Superconducting AC Machines for Hybrid Electric Aircraft[END_REF], [START_REF] Devillers | A review of subdomain modeling techniques in electrical machines: Performances and applications[END_REF]. For common radial-flux structures such as induction machines, permanent magnet or wound rotor synchronous machines, cylindrical 2-D models are very efficient and have been extensively developed in the literature [START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF]- [START_REF] Lubin | Improved Analytical Model for Surface-Mounted PM Motors Considering Slotting Effects and Armature Reaction[END_REF]. These topologies can also be modelled efficiently in 3-D [START_REF] Meessen | General Formulation of Fringing Fields in 3-D Cylindrical Structures Using Fourier Analysis[END_REF]. Nevertheless, axial flux machines do not have any symmetries that allow them to be accurately modelled in 2-D. Good results can be obtained under some assumptions to consider the end effects [START_REF] Lubin | Development of a 2-D Analytical Model for the Electromagnetic Computation of Axial-Field Magnetic Gears[END_REF], though 3-D calculation is more suitable for that kind of problems. For the simulation of axial flux machines in 3-D there is typically two solutions, the first one being the use of a cylindrical frame which will lead to the most precise results but requires to manipulate Bessel functions [START_REF] Lubin | Improved 3-D Analytical Model for Axial-Flux Eddy-Current Couplings With Curvature Effects[END_REF]. A second option would be the use of a cartesian coordinate system by linearisation at the average radius symplifying the model by neglecting the curvature effects [START_REF] Dolisy | A New Analytical Torque Formula for Axial Field Permanent Magnets Coupling[END_REF], [START_REF] Lubin | 3-D Analytical Model for Axial-Flux Eddy-Current Couplings and Brakes Under Steady-State Conditions[END_REF]. These assumptions will be described in more details further in this article.

Moreover, flux modulation inductors have been analysed with analytical models, the radial flux topology has been solved in [START_REF] Malé | Analytical calculation of the flux density distribution in a superconducting reluctance machine with HTS bulks rotor[END_REF] while the axial flux topology with disk shaped bulks has been solved in [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF]. In either case, the problem was adressed by superposing two 2-D models. The solution has been obtained by making the product of the magnetic field created by the superconducting coil(s) by a modulation function calculated by considering the bulks as perfectly diamagnetic domains. Both of these models have shown excellent results, but the perfect diamagnetism hypothesis is not necessarily verified in practice because it assumes an infinite or at least very high critical current density in the bulk [START_REF] Colle | Test of a Flux Modulation Superconducting Machine for Aircraft[END_REF], [START_REF]Relationship between local microstructure and superconducting properties of commercial YBa2Cu3O7-$\updelta$ bulk[END_REF].

More recently, various research groups have developed models including subdomains with a spatially variable permeability, thus reducing the number of domains required to simulate machine with slots or variable reluctances [START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF]- [START_REF] Zhang | Analytical Prediction of Halbach Array Permanent Magnet Machines Considering Finite Tooth Permeability[END_REF]. These models are based on the convolution therorem that allows considering complex space distribution of the magnetic permeability in the subdomains. Thanks to this method, it is possible to simply take into account the saturation of the magnetic circuit [START_REF] Djelloul-Khedda | Nonlinear Analytical Prediction of Magnetic Field and Electromagnetic Performances in Switched Reluctance Machines[END_REF], [START_REF] Zhao | Analytical Modeling and Comparison of Two Consequent-Pole Magnetic-Geared Machines for Hybrid Electric Vehicles[END_REF] and numerous radial-flux topologies have been modelled in recent years [START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF], [START_REF] Zhang | Analytical Prediction of Halbach Array Permanent Magnet Machines Considering Finite Tooth Permeability[END_REF].

A 3-D model of the axial flux modulation machine shown in Fig. 1 is presented in this paper. In section 2, the electromagnetic problem as well as the required assumptions are presented. In section 3, the magnetic field distribution of the inductor is calculated with the variable separation method, Fourier series analysis and the convolution theorem for perfectly diamagnetic bulks. The section 4 presents the calculation of the machine power from the electromotive force as well as the eddy current losses in the armature. The results and the accuracy of this model will be discussed in the section 5. Finally, the model is improved in section 6 using the Bean model to consider the penetration of the magnetic field in the HTS bulks.

II. PROBLEM DESCRIPTION

The axial flux modulation inductor with ring segment shaped bulks has an inherent 3-D geometry and cannot be simply described in 2-D as was done for the axial flux inductor with round shaped bulk using two axisymmetric models [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF]. Yet, some symmetries can still be exploited to simplify the problem. Indeed, due to the periodicity along azimuthal direction and to the symmetry in the axial direction, only half of a pole pair need to be considered. Thus, the geometrical problem to be solved according to the symmetries is presented in Fig. 2 with the geometrical parameters. The inner and outer radius of the bulk are Ri and Re, respectively. The inner radius of the coil is Rb and its thickness Wb. The cylindrical frame is centred at the middle of the machine so Z1 represent half of the bulk thickness and Zb half of the coil length. Finally, a bulk has an opening angle of 𝛼2𝜋/𝑝 with p being the number of pole pairs and α is the bulk pole-arc to pole-pitch ratio.

To simplify the resolution, the flux density distribution produced by the inductor will be obtained by separately calculating the magnetic field created by the coil and the deflection of the flux lines due to the bulks. The first problem is a 2-D axisymmetric model of the coil alone, which has already been described in [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF] and will not be detailed in this paper. The second problem is a 3-D model of a bulk placed in a homogeneous external magnetic field from which a modulation function will be extracted. This modulation function represents the ability of the bulk to divert flux lines. Similar methods have already been applied in 2-D for various flux modulation machines [START_REF] Malé | Analytical calculation of the flux density distribution in a superconducting reluctance machine with HTS bulks rotor[END_REF], [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF]. The calculation of the modulation function will be described in detail in this article. Ultimately, the actual flux density in the air gap is obtained by multiplying the magnetic field created by the coil by the modulation function. Considering the geometry, the more intuitive way to calculate the modulation function would be in a cylindrical coordinate system (r,θ,z). However, it would lead to a complex system of differential equations to solve and requires the use of the Bessel functions. It has been observed in several references [START_REF] Dolisy | A New Analytical Torque Formula for Axial Field Permanent Magnets Coupling[END_REF], [START_REF] Lubin | 3-D Analytical Model for Axial-Flux Eddy-Current Couplings and Brakes Under Steady-State Conditions[END_REF], [START_REF] De La Barriere | 3-D Formal Resolution of Maxwell Equations for the Computation of the No-Load Flux in an Axial Flux Permanent-Magnet Synchronous Machine[END_REF] that the curvature effects can be neglected by slightly changing the geometry from a cylindrical problem to a cartesian problem using a linearisation at the mean radius 𝑅 = (𝑅 + 𝑅 )/2. Thus, the new geometry is presented in Fig. 3, the y-coordinate represents the radial direction while the x-coordinate represents the azimuthal direction. The x-coordinate is related to the θ-coordinate by 𝑥 = 𝑅 𝜃, the y-coordinate is related to the r-coordinate by 𝑦 = 𝑟 -𝑅 and the z-coordinate is identical for both frames. The height of the rectangle shaped bulk is 𝐷 = 𝑅 -𝑅 and its length is 𝛼𝜏, 𝜏 = 𝑅 being the length of a pole pair. The impact of these assumptions on the model precision will be discussed in section 5.

III. . 3-D SEMI-ANALYTICAL MODEL

The cross-sectional views of the cartesian problem for the modulation function are shown in Fig. 4 and Fig. 5 in planes (xz) and (yz) respectively. The external homogeneous magnetic field will be modelled here by a permanent magnet placed on a ferromagnetic plate with an assumed infinite permeability creating a magnetic field along the z-axis which is modulated by the HTS shield. The iron plate and the magnet must be located far enough from the bulk surface otherwise the plate disturbs the magnetic field distribution and the result would not correspond to the modulation of a homogeneous field. This distance has been set at 𝑍 = 5𝑅 to avoid such issues. The x and z components of the modulation function will be obtained by dividing the x and z components of the flux density in the presence of the bulk with the z-component of the flux density in the absence of the bulk. The ferromagnetic plate and the magnet surface are at a distance Z3 and Z2 from the bulk centre respectively. The magnet remanent magnetization is written Brem and the magnetization 𝑀 ⃗ is in the opposite direction to the z-axis, 𝑀 ⃗ = -𝐵 𝜇 ⁄ 𝑢 ⃗. It should be noted that the uniform applied magnetic field could have been modelled by applying a Dirichlet condition on the scalar potential on the boundary at Z2. However, it has been decided to model it through a magnet for the sake of understanding.

In this section, the bulk is considered as an almost perfect diamagnetic material with a relative permeability 𝜇 = 10 while the other parts of the problem have a relative permeability equal to 1. The whole geometry is divided into three regions:

The HTS bulk region (region I)

The air gap region (region II)

The magnet region (region III) Since there is no current in any of the regions, 𝐻 ⃗ , the magnetic field in region i, will be expressed in term of a magnetic scalar potential Φ i , which is defined as:

𝐻 ⃗ = -∇ ⃗ × 𝛷 (𝑖 = 𝐼, 𝐼𝐼, 𝐼𝐼𝐼) (1)
A. Boundary conditions Solving this problem requires the definitions of boundary conditions. The boundaries on the x-direction and y-direction will define the Sturm-Liouville problem eigenvalues and define the general solution form while the boundary and interface conditions in the z-direction will be used to calculate the unknown coefficients.

Due to the bulks periodicity in the x-direction, the magnetic problem presents even boundary conditions in the x-direction. Thus, the domains are limited to 𝑥 ∈ [-𝜏/2 ; 𝜏/2], and the boundary conditions are:

𝛷 (-𝜏/2, 𝑦, 𝑧) = 𝛷 (𝜏/2, 𝑦, 𝑧) (2)
The conditions in the y-direction appear less naturally. Indeed, far enough from its surface, the bulk has no more influence on the magnetic field. Hence, we define the limits of the domain in the y-direction as -𝐿/2 and 𝐿/2 on which the boundary conditions are imposed as: 𝛷 (𝑥, -𝐿/2, 𝑧) = 𝛷 (𝑥, 𝐿/2, 𝑧) = 0

(3)

The value of L, which is the domain size, must be chosen carefully. If the domain limits are too close to the bulk, the boundaries will have a significant impact on the result and will lead to errors. However, if L is too large, too many harmonics would be required to solve the problem, thus increasing the computational cost of the model. For the rest of the article, the value 𝐿 = 3𝐷 has been chosen after trials and errors.

Eventually, in the z-direction, the ferromagnetic plate and the axial symmetry in 𝑧 = 𝑍 and 𝑧 = 0, respectively, lead to the orthogonality of the magnetic field for both cases that expressed in scalar potential give: 𝛷 (𝑥, 𝑦, 0) = 0 (4)

𝛷 (𝑥, 𝑦, 𝑍 ) = 0 (5)

B. Separation of variables and Fourier series representation

The magnetic scalar potential solution in each domain is obtained using the separation of variables. Considering the boundary conditions in the x and y directions, Φ i is expressed in terms of 2-D complex Fourier series:

𝛷 (𝑥, 𝑦, 𝑧) = 𝜑 , (𝑧)𝑒 𝑒 (6)
The harmonics ranks are denoted by n and k for the x and y directions respectively. N represents the highest spatial harmonic considered in both directions. Ideally, the value of N should be infinite. In practice, it is selected to satisfy the balance between accuracy and computational cost. The impact of N on the computation time and model precision will be discussed further. Similarly the magnetic field 𝐻 ⃗ components are also expressed as Fourier series:

𝐻 ⃗ = 𝐻 𝑢 ⃗ + 𝐻 𝑢 ⃗ + 𝐻 𝑢 ⃗ (7

C. Convolution theorem

The magnetic permeability is homogeneous in the regions II and III and equal to 𝜇 . However in the region I, the permeability is equal to 𝜇 𝜇 with 𝜇 = 10 at the bulk location and μ0 elsewhere as shown in Fig. 6. The permability in the region I is inhomogeneous and must be decomposed in Fourier series as well:

𝜇 (𝑥, 𝑦) = 𝜇̂ , 𝑒 𝑒 (15)
Where the coefficients 𝜇̂ , are expressed as:

𝜇̂ , = ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ 𝜇 1 + 𝛼 𝐷 𝐿 (𝜇 -1) 𝑓𝑜𝑟 𝑛 = 𝑘 = 0 , 𝜇 𝐷 𝜋𝑛𝐿 (𝜇 -1) sin(𝑛𝜋𝛼), 𝑓𝑜𝑟 𝑛 ≠ 0 & 𝑘 = 0 𝜇 𝛼 𝜋𝑘 (𝜇 -1) 𝑠𝑖𝑛 𝑘𝜋 𝐷 𝐿 𝑓𝑜𝑟 𝑛 = 0 & 𝑘 ≠ 0 , 𝜇 𝜋𝑛𝑘 (𝜇 -1) sin(𝑛𝜋𝛼) sin 𝑘𝜋 𝐷 𝐿 , 𝑓𝑜𝑟 𝑛 ≠ 0 & 𝑘 ≠ 0 (16)
It should be noted here that the bulk relative permeability 𝜇 cannot be equal to zero as the resolution of the magnetostatic problem will involve the inverse of the permeability. Therefore, the permeability value must be small enough to behave almost as a perfect diamagnetic material. Conversely, a too small value of 𝜇 would cause numerical problems. Thus, the value 𝜇 = 10 has been chosen after trials.

The constitutive relation between 𝐵 ⃗ and 𝐻 ⃗ in nonmagnetized regions is:

𝐵 ⃗ = 𝜇 (𝑥, 𝑦)𝐻 ⃗ (17)
The right part of ( 17) being the product of two Fourier series decomposition, the convolution theorem can be applied allowing to express the coefficients of each component of 𝐵 ⃗ as follows:

𝐵 , (𝑧) = 𝜇̂ , 𝐻 , (𝑧) (18) 𝐵 , (𝑧) = 𝜇̂ , 𝐻 , (𝑧) (19) 𝐵 , (𝑧) = 𝜇̂ , 𝐻 , (𝑧) (20)
It can be seen from equations ( 18) to ( 20) that each harmonic of 𝐵 ⃗ depends on all harmonics of 𝐻 ⃗ , whereas for homogeneous permeability regions, only harmonics of the same rank are interdependent. Hence, the application of Maxwell's equations in this region leads to a system of differential equations to solve.

D. Matrix form

In order to solve the system of differential equations, ( 18), [START_REF] Meessen | General Formulation of Fringing Fields in 3-D Cylindrical Structures Using Fourier Analysis[END_REF] and [START_REF] Lubin | Development of a 2-D Analytical Model for the Electromagnetic Computation of Axial-Field Magnetic Gears[END_REF] are rewritten in matrix form:

𝑩 𝒙 = 𝝁 𝑯 (21)
𝑩 𝒚 = 𝝁 𝑯 (22)
𝑩 𝒛 = 𝝁 𝑯 (23)
In matrix form, every quantity 𝑋(𝑥, 𝑦, 𝑧) expressed as Fourier series (𝐵 , 𝐵 , 𝐵 , 𝐻 , 𝐻 , 𝐻 and 𝛷) with coefficients 𝑋 , (𝑧) can be represented in matrix form by gathering the coefficient in a column vector X as follows:

𝑿 = 𝑋 = 𝑋 , (𝑧) (24)
The relations between the vector element index m and the harmonic ranks n and k are:

𝑘 = 𝑓𝑖𝑥 𝑚 -1 2𝑁 + 1 -𝑁 (25)
𝑛 = 𝑚𝑜𝑑(𝑚 -1,2𝑁 + 1) -𝑁 (26)
The functions fix(u) and mod(u,v) are the rounding toward zero of u and the remainder of the division of u and v, respectively.

The matrix 𝝁 is the permeability convolution matrix and is defined according to [START_REF] Lubin | Improved Analytical Model for Surface-Mounted PM Motors Considering Slotting Effects and Armature Reaction[END_REF], [START_REF] Meessen | General Formulation of Fringing Fields in 3-D Cylindrical Structures Using Fourier Analysis[END_REF] and [START_REF] Lubin | Development of a 2-D Analytical Model for the Electromagnetic Computation of Axial-Field Magnetic Gears[END_REF]:

𝝁 = 𝜇 , = 𝜇̂ , (27)
The relations between the matrix element indexes m, l and the harmonic ranks n, k, n', k' are:

𝑘 = 𝑓𝑖𝑥 𝑚 -1 2𝑁 + 1 -𝑁 (28)
𝑛 = 𝑚𝑜𝑑(𝑚 -1,2𝑁 + 1) -𝑁 (29)
𝑘′ = 𝑓𝑖𝑥 𝑙 -1 2𝑁 + 1 -𝑁 (30)
𝑛′ = 𝑚𝑜𝑑(𝑙 -1,2𝑁 + 1) -𝑁 (31)
In the regions II and III, the matrix form will also be used although the permeability convolution matrixes 𝝁 and 𝝁 are much simpler:

𝝁 = 𝝁 = 𝜇 𝑰 (𝟐𝑵 𝟏) 𝟐 (32)
With 𝑰 (𝟐𝑵 𝟏) 𝟐 the (2𝑁 + 1) 𝑏𝑦 (2𝑁 + 1) identity matrix.

E. Magnetostatic field equation

Using the matrix form described above, the relation (1) becomes:

𝑯 𝒙 = -𝑗 2𝜋 𝜏 𝑵 𝒙 𝜱 𝒊 (33)
𝑯 𝒚 = -𝑗 2𝜋 𝐿 𝑵 𝒚 𝜱 𝒊 (34)
𝑯 𝒛 = - 𝜕𝜱 𝒊 𝜕𝑧 (35)
With the matrices Nx and Ny being defined as:

𝑵 𝒙 = 𝑁 , = 𝑚𝑜𝑑(𝑚 -1,2𝑁 + 1) -𝑁 𝑓𝑜𝑟 𝑚 = 𝑙 0 𝑓𝑜𝑟 𝑚 ≠ 𝑙 (36)
𝑵 𝒚 = 𝑁 , = 𝑓𝑖𝑥 𝑚 -1 2𝑁 + 1 -𝑁 𝑓𝑜𝑟 𝑚 = 𝑙 0 𝑓𝑜𝑟 𝑚 ≠ 𝑙 (37)
The magnetostatic field equation in scalar potential is then derived from the Gauss law for magnetism:

∇ • 𝐵 ⃗ = 0 (38)
Substituting ( 21), ( 22) and ( 23) into (38) leads to the following equation since the divergence of the magnetization is equal to zero in each region:

𝝁 𝒊 𝑯 𝒛 + 𝑗 𝑵 𝒙 𝝁 𝒊 𝑯 + 𝑗 𝑵 𝒚 𝝁 𝒊 𝑯 = 0 (39)
Finally, the magnetostatic equation to solve in scalar potential is obtained by substituting ( 33), ( 34) and ( 35) in [START_REF] Ainslie | Modelling and comparison of trapped fields in (RE)BCO bulk superconductors for activation using pulsed field magnetization[END_REF]:

𝜕 𝜱 𝒊 𝜕𝑧 -(𝑽 ) 𝜱 𝒊 = 0 (40)
Where in the region I:

𝑽 = (𝝁 𝑰 ) 𝑵 𝒙 𝝁 𝑰 𝑵 𝒙 + (𝝁 𝑰 ) 𝑵 𝒚 𝝁 𝑰 𝑵 𝒚 (41)
Whereas in the region II and III, the matrix V is diagonal:

𝑽 = 𝑽 = 2𝜋 𝜏 (𝑵 𝒙 ) 𝟐 + 2𝜋 𝐿 𝑵 𝒚 𝟐 (42)
F. General solution to the magnetostatic equation In each region i, the general solution of ( 40) is:

𝜱 𝒊 = 𝑒 𝑽 𝑪 𝟏 𝒊 + 𝑒 𝑽 𝑪 𝟐 𝒊 (43)
Where 𝑪 𝟏 𝒊 and 𝑪 𝟐 𝒊 are column vectors with unknown coefficients of each region. In order to evaluate the matrix exponentials involved in (43), eigendecomposition is applied to 𝑒 𝑽 :

𝑒 𝑽 = 𝑷 𝒊 𝑒 𝑫 (𝑷 𝒊 ) 𝟏 (44)
Where 𝑫 and 𝑷 𝒊 are the diagonal eigenvalue matrix and the eigenvector matrix of 𝑽 , respectively. As a result, the general solution (43) becomes:

𝜱 𝒊 = 𝑷 𝒊 𝑒 𝑫 (𝑷 𝒊 ) 𝟏 𝑪 𝟏 𝒊 + 𝑷 𝒊 𝑒 𝑫 (𝑷 𝒊 ) 𝟏 𝑪 𝟐 𝒊 (45)
As 𝑽 and 𝑽 are diagonal matrices, their eigenvector matrices are equal to the identity matrix 𝑰 (𝟐𝑵 𝟏) 𝟐 and 𝑽 , 𝑽 , 𝑫 and 𝑫 are equal.

G. Calculation of the coefficients

The determination of the six vectors of coefficients of the general solution will require the resolution of a system of linear equations and the inversion of a 6(2𝑁 + 1) 𝑏𝑦 6(2𝑁 + 1) matrix [START_REF] Sprangers | Magnetic Saturation in Semi-Analytical Harmonic Modeling for Electric Machine Analysis[END_REF]. Thus, for 𝑁 = 10, the solution will involve the inversion of a 2646 by 2646 matrix which would have an important computational cost. Nevertheless, it is possible to drastically reduce the size of the linear system by rewriting (45) according to the boundary conditions in the z-direction.

In matrix form, the boundary conditions ( 4) and ( 5) as well as the conditions between each domains lead to the following equations:

𝜱 𝑰 | = 𝟎 (46) 𝜱 𝑰 | = 𝜱 𝑰𝑰 | (47) 𝝁 𝑰 𝜕𝜱 𝑰 𝜕𝑧 = 𝜇 𝜕𝜱 𝑰𝑰 𝜕𝑧 (48) 𝜱 𝑰𝑰 | = 𝜱 𝑰𝑰𝑰 | (49) 𝜕𝜱 𝑰𝑰𝑰 𝜕𝑧 - 𝜕𝜱 𝑰𝑰 𝜕𝑧 = 𝑴 (50) 𝜱 𝑰𝑰𝑰 | = 𝟎 ( 51
)
Where 𝑴 is the vector representing the magnetization in the region III expressed as:

𝑴 = 𝑀 = -𝐵 𝜇 , 𝑚 = 2𝑁 + 2𝑁 + 1 0 , 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (52)
In the region I, according to the conditions (46) and (47) and to the superposition theorem, the general solution can be rewritten as:

𝜱 𝑰 = 𝑷 sinh (𝑫 𝑧) sinh (𝑫 𝑍 ) (𝑷 𝑰 ) 𝟏 𝒂 (53)
In the region II, according to the conditions (47) and (48):

𝜱 𝑰𝑰 = sinh (𝑫 (𝑧 -𝑍 )) sinh (𝑫 (𝑍 -𝑍 )) 𝒂 + sinh (𝑫 (𝑧 -𝑍 )) sinh (𝑫 (𝑍 -𝑍 )) 𝒃 (54)
Eventually, in the region III, according to (49) and (51):

𝜱 𝑰𝑰𝑰 = sinh (𝑫 (𝑧 -𝑍 )) sinh (𝑫 (𝑍 -𝑍 )) 𝒃 (55)
The function 𝑠𝑖𝑛ℎ(𝑿) is the hyperbolic sine of the matrix 𝑿 and defined similarly to the hyperbolic sinus for scalar values. The function hyperbolic cosine of the matrix 𝑿 𝑐𝑜𝑠ℎ(𝑿) can also be defined.

Thanks to the reformulations of the general solutions, the relations (53), ( 54) and (55) show only 4 coefficients vectors to be determined, 𝒂 , 𝒂 , 𝒃 and 𝒃 . Furthermore, substituting (53) and ( 54) into (47) results in the relation:

𝒂 = 𝒂 (56)
In a similar way, (54) and ( 55) into (49) lead to:

𝒃 = 𝒃 (57)
As a result, there are only two coefficients vectors remaining. Ultimately, the conditions (48) and ( 50) are used to calculate 𝒂 and 𝒃 . These two relations are expressed as the following system:

𝑬 𝑬 𝑬 𝑬 𝟐𝟐 𝒂 𝒃 = 𝟎 (𝑫 ) 𝑴 (58)
Where the component of the matrix E are:

𝑬 = 𝝁 𝑰 𝑷 𝑫 cosh (𝑫 𝑍 ) sinh (𝑫 𝑍 ) (𝑷 𝑰 ) 𝟏
The linear system (58) has 2(2𝑁 + 1) unknown coefficients and requires, for 𝑁 = 10, the inversion of only a 882 by 882 matrix. It should be noted that the zero harmonics cause the matrix E to be non-invertible since the solution (43) is only valid for non-zero harmonics. To overcome this issue, a very small value is given to the zero harmonics such as 10 -5 instead of absolute zero [START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF].

H. Calculation of the flux modulation function

The flux modulation function is to be calculated in the region II as it corresponds to the location of the airgap and the armature. In order to compute the torque, only the axial component of the magnetic flux density is required while the eddy current losses in the windings involve the tangential component in the θ-direction as well. Thus, in the cartesian frame, only the x-component and z-component of the magnetic field in the region II have to be calculated. As the magnetic scalar potential is now known in each region, the vectors 𝑩 𝒙 and 𝑩 𝒛 are obtained by:

𝑩 𝒙 𝑰𝑰 = -𝜇 𝑗 2𝜋 𝜏 𝑵 𝒙 𝜱 𝑰𝑰 (63) 𝑩 𝒛 𝑰𝑰 = -𝜇 𝜕𝜱 𝑰𝑰 𝜕𝑧 (64)
The transition from the matrix form to the spatial form is made using the relations:

Where B0 is obtained from the Ampere law:

𝐵 = - 𝐵 (𝑍 -𝑍 ) 𝑍 (69)

I. Magnetic field of the whole inductor

Once the flux modulation function is determined in the cartesian frame, the transition to the cylindrical frame is achieved by using the following relations between the coordinates:

𝑟 = 𝑅 + 𝑦 (70) 𝜃 = 𝑥 𝑅 (71)
In the cylindrical frame, the z-component of the modulation function is identical to that of the cartesian frame 𝐹𝑀 (𝑥, 𝑦, 𝑧) and the θ-component is equal to the x-component 𝐹𝑀 (𝑥, 𝑦, 𝑧).

Ultimately, the magnetic flux density generated by the whole inductor is obtained from the following relations:

𝐵 (𝑟, 𝜃, 𝑧) = 𝐵 (𝑟, 𝑧) 𝐹𝑀 (𝑥, 𝑦, 𝑧) (72)
𝐵 (𝑟, 𝜃, 𝑧) = 𝐵 (𝑟, 𝑧) × 𝐹𝑀 (𝑥, 𝑦, 𝑧)

Where 𝐵 and 𝐵 are the z-component and the θ-component respectively of the magnetic flux density produced by the inductor while 𝐵 is the magnetic flux density in the z-direction generated by the HTS-coil alone. 𝐵 is obtained using the 2-D axisymmetric model presented in [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF]. IV. ELECTROMAGNETIC POWER AND EDDY-CURRENT LOSSES CALCULATION In a three-phase electrical machine, the electromagnetic torque 𝑇 can be calculated from:

𝑇 = 3𝐸 𝐼 𝛺 ( 74
)
Where Ω is the rotational speed , Irms is the rms value of the armature electrical current and E1 is the rms value of the fundamental component of the electromotive force which is obtained from:

𝐸 = 𝑝𝛺𝜆 √2 ( 75
)
Where λ1 is the phase flux linkage amplitude due to the fundamental spatial harmonic of the inductor. In the present case, the air cored stator complexify the calculation of 𝜆 . Indeed, the inhomogeneity of the magnetic field in the armature region means that each turn of the windings has a different electromotive force. The expression used to account for the magnetic field inhomogeneities and calculate λ1 is:

𝜆 = 𝑁 𝐾 𝐾 𝐿 × 𝐾 (𝑟)𝐵 (𝑟, 𝑧) cos(𝑝𝜃) 𝑟𝑑𝑟𝑑𝜃𝑑𝑧 (76)
The geometric parameters involved in (76) are detailled in Fig. 7. 𝑅 and 𝑅 are respectively the inner and outer radius of the inner turn of a coil composing the armature. 𝑊 is the width of a coil conductor bundle in the θ-direction, 𝐿 the thickness of a coil along the z-direction, 𝑊 the thickness of the coil head, 𝜃 the opening angle of a stator coil and 𝑁 the number of turns per phase. Several kinds of windings can be modelled from (76) through the winding coefficients 𝐾 and 𝐾 . 𝐾 is the distribution factor while 𝐾 is the pitch factor whose calculation is detailed for several types of axial flux air cored windings in [START_REF] Kamper | Analysis and Performance of Axial Flux Permanent-Magnet Machine With Air-Cored Nonoverlapping Concentrated Stator Windings[END_REF]. In the following, only the case of a concentrated winding in which 𝑊 = 𝑊 will be detailed but conventional overlapping stator windings can be considered as well.

The value of 𝐾 can be classicaly calculated for differents windings according to [START_REF] Bianchi | Use of the star of slots in designing fractional-slot single-layer synchronous motors[END_REF]. In the case of concentrated windings, 𝐾 is expressed as:

𝐾 = 2 𝑝𝜃 sin 𝑝 2 (𝜃 -𝜃 ) 𝑠𝑖𝑛 𝑝𝜃 2 (77)
Where the angle 𝜃 is defined as:

𝜃 = 𝑊 𝑅 (78)
The fundamental spatial harmonic of the inductor 𝐵 (𝑟, 𝑧) is obtained from the previous section inductor model as follows:

𝐵 (𝑟, 𝑧) = 𝐵 (𝑟, 𝑧) × 𝐵 , (𝑧) 𝐵 𝑒 + 𝐵 , (𝑧) 𝐵 𝑒 (79)
The integrals over r and θ in (76) allow the calculation of the flux under a pole. Since the flux value changes in the z-direction, the third integral over z is made to calculate the average of this flux in the armature.

The function 𝐾 (𝑟) accounts for the coil heads influence on the flux linkage. Indeed, the outermost region of the windings enclose more flux than the innermost region. Thus, 𝐾 (𝑟) is represented in Fig. 8. The influence of the coil heads depends on the machine arrangement, of which the two main ones are shown in Fig. 9. On the one hand, in the case of Fig. 9 (a) where the coil heads are not covered by the bulks, their impact on the flux linkage is small and not considering 𝐾 (𝑟) would not lead to a significant error on the torque. On the other hand, for Fig. 9 (b) where the coil heads are covered, omitting 𝐾 (𝑟) could lead to an underestimation of the torque depending on the machine geometry. In any case, the computation of the triple integral of (76) is done numerically.

Another important consequence of the air cored armature topology is that the copper conductors are directly subjected to the magnetic field of the inductor. Hence, the eddy-current losses in the windings are higher than in iron cored armatures and the calculation of these losses is an important part of this machine sizing process.

The computation of the losses is done by neglecting the magnetic field generated by the armature currents and the associated losses as the main contribution to these losses comes from the inductor. The calculation is done for each turn by using the formula for the losses in a conductor subjected to an external magnetic field which can be found in [START_REF] Tourkhani | Accurate analytical model of winding losses in round Litz wire windings[END_REF]. The latter method for axial flux machines with an air cored armature is described in detail in [START_REF] Wang | Calculation of eddy current loss in axial field permanent-magnet machine with coreless stator[END_REF].

V. COMPARISON WITH 3-D FEM SIMULATIONS

In this section, the precision and efficiency of the proposed method will be assessed through the comparison of the three following models:



The semi-analytical model (SAM) in which the inductor magnetic field is obtained from the model proposed in section 3 while the torque and the eddy current losses are calculated from the armature model of section 4. This model is implemented on the software MATLAB 2019.



A semi-FEM model (SFEM) in which the inductor magnetic field is obtained from a 3-D stationary model in A-formulation with the software COMSOL Multiphysics while the torque and the eddy current losses are calculated from the armature model of section 4.



A complete FEM model (CFEM) in which the inductor and the armature are simulated on COMSOL with the torque being calculated using the Laplace force on the windings. The eddy current losses are not computed for this model. Fig. 10 shows the mesh and simulated geometry of the CFEM (a) and SFEM (b). As the winding layout represented is concentrated, the armature periodicity is different from the rotor periodicity. Therefore, the only symmetry plane is (𝑟𝜃) in the CFEM whereas the SFEM model only need to consider 1 4 ⁄ of a bulk and the corresponding portion of the coil.

This study will identify the consequences of the different assumptions from which the SAM is derived. The different models will be tested with the parameters listed in Table . 1.

A. Flux density distribution in the air gap Fig. 11 and Fig. 12 show the axial and azimuthal components of the flux density in the air gap at the mean radius obtained with the SAM for 5 and 20 harmonics and by the numerical models (SFEM and CFEM). For 𝑁 = 20, a good agreement is observed with the numerical simulations while 𝑁 = 5 also gives a reasonable estimate of the flux density distribution but does not allow predicting precisely the maximum value of the tangential component.

Fig. 13 and Fig. 14 show the axial component of the flux density in the air gap around the bulk outer and inner radii respectively. In both cases, the effects of the linearisation at the mean radius can be seen as the analytical curves slightly differ from the numerical curves. Nonetheless, the accuracy of the results seems to be more than satisfactory. Additionally, Fig. 15 shows the flux density norm in the winding region computed by the CFEM. It appears that the flux density is higher on the top portion of the armature compared to the lower part.

B. Precision on the torque and the eddy current losses

Observing the flux density distribution in the air gap is not sufficient to properly assess the model precision. For this purpose, the torque and the eddy current losses are more relevant values. Indeed, the torque indicates the overall accuracy of the fundamental over the whole armature area while the eddy current losses inform on the harmonic's precision.

First of all, Fig. 16 shows the computational time of the different models as a function of the number of harmonics.

Various elements can be observed. Firstly, the computation time of the SFEM increases a little with the number of harmonics while the CFEM time is constant. This is because the spectral analysis is only made for the eddy current losses computation in the SFEM while they are not calculated in the CFEM. In addition, the computation time of the SFEM is more important than CFEM because of this spectral analysis. Subsequently, the calculation time of the SAM rises exponentially with N as it strongly increases the matrices sizes. However, the resolution of the SAM remains below the numerical methods for 𝑁 < 20 and is of the order of a second for 𝑁 < 10. Thus, as long as it is accurate enough for 10 harmonics, the proposed model is significantly faster than numerical methods.

In order to ascertain the latter condition for the torque, Fig. 17 shows the relative error in absolute value of the SAM compared to the SFEM and the CFEM as a function of N. The CFEM is here taken as a reference since it simulates the exact geometry of the inductor and armature. Thus, the comparison of the SAM and the CFEM provides the overall accuracy of the model while the comparison of the SAM and the SFEM only gives the inductor precision. The inductor model appears to be precise even for 𝑁 = 5 as the error between the SAM and the SFEM is below 4 %. Additionally, for 𝑁 = 30, the error between the SAM and the SFEM drops to 0.15 %. The overall accuracy of the model is also good as the maximum error is below 7 %. Eventually, it should be noted that the error between the SFEM and the CFEM is not influenced by N as it corresponds to the armature model error and accounts for 3.2 % of the overall error. Moreover, Table . 2 compares the torque calculated by the different models for the machine of Table . 1 with and without the covering of the coil heads as presented on Fig. 9. This table illustrate the impact of the coil heads on the torque as was described in section 4. When the coil heads are not covered, neglecting their impact, as done in the SFEM (b), lead to an additional 4 % error. In contrast, when the coil heads are covered, the additional error made by the SFEM (b) is 15.5 %. This table also shows the significant performance improvement that can be achieved by covering the coil heads. Eventually, Fig. 18 presents the normalized eddy current losses obtained by the SAM and the SFEM as a function of the number of harmonics. For this analysis, the reference value is represented by the losses computed by the SFEM with 𝑁 = 30 as it considers the exact magnetic field with a large number of harmonics. The first thing to observe here is the inaccuracy of the SAM for 𝑁 = 5. Indeed, the value calculated with the SAM is 27 % lower than the reference value but the spread decreases rapidly as the number of harmonics increase. Hence, while 5 harmonics are sufficient to get an accurate idea of the torque, taking less than 10 harmonics leads to a significant error on the eddy current losses. Nevertheless, it is important to distinguish between the two sources of error here. Firstly, the eddy current losses calculation requires a sufficient number of harmonics since high ranks harmonics still generate losses that are not negligible. This source of error is present even if the magnetic field distribution is completely accurate and corresponds to the error made by the SFEM. Thus, above N = 20, this error is almost null while it accounts for 18.4 % of the total error at 𝑁 = 5. Subsequently, the other source of error is the inaccuracy of the flux density distribution obtained from the SAM. As all harmonics are interdependent, a large number of harmonics is required to calculate their value accurately. It can therefore be deduced that this source of error is responsible for 8.6 % of the total error at 𝑁 = 5.

C. Influence of the geometrical parameters

According to the previous results, the SAM with 𝑁 = 10 would appear to be a good compromise between accuracy and calculation speed. However, the geometrical parameters presented in Table . 1 are the result of an optimisation from previous works using the SFEM model. Therefore, it is important to investigate the sensitivity of the model's accuracy with the variation of the geometric parameters. For this purpose, the following parameters will be explored: Additionally, the study will be performed for the cases represented in Fig. 9 (a) and Fig. 9 (b) with the coil heads covered or not by the bulks. Within this context, the aspect ratio of the bulk is defined as follows:

𝐴𝑠𝑝𝑒𝑐𝑡 𝑟𝑎𝑡𝑖𝑜 = 𝑅 -𝑅 𝛼𝜏 (80)
Since all the parameters explored affect the aspect ratio, the latter seems to be a relevant value to understand the results. Relative error on the torque Fig. 19 shows the relative error on the torque of the SAM with 𝑁 = 10 compared to the SFEM function of the aspect ratio. Several elements can be observed. Firstly, there is a clear correlation between the relative error and the aspect ratio as linearisation at the mean radius has less impact on bulks with low aspect ratios. For the same reason, the model is less precise for machines with coil heads covered. Secondly, it can be observed that the SAM has a strong tendency to underestimate the torque. Lastly, the SAM accuracy remains good even for large aspect ratios as the maximum error is inferior to 6 %. It should also be noted that, as for today, the use of bulks with very small or very large aspect ratios is unlikely since ring segment shaped bulks are mostly manufactured by machining disk or square shaped bulks [START_REF] Dorget | Review on the Use of Superconducting Bulks for Magnetic Screening in Electrical Machines for Aircraft Applications[END_REF]. Thus, HTS bulks with an aspect ratio close to 1 are easier to produce. Fig. 20 shows the relative error on the eddy current losses of the SAM with 𝑁 = 10 compared to SFEM function of the aspect ratio. The same observations as in Fig. 19 can be made. Indeed, there is a correlation between the error and the aspect ratio and the SAM tends to underestimate the losses. However, it is clear that the aspect ratio is not the only parameter influencing the error. Additionally, the average error is much higher for the losses than for the torque as was observed for the parameters of Table . 1. Indeed, as previously said, every harmonics of the SAM model are interdependent, which means that to obtain an accurate value of 10 harmonics, many more must be calculated.

VI. FLUX PENETRATION IN THE BULK

The method described above as well as the numerical simulation are based on the main assumption that the bulk are perfectly diamagnetic which corresponds to type II superconductors with an infinite critical current density. However, even the most efficient superconductors does not have a perfect diamagnetic behaviour [START_REF] Dorget | Review on the Use of Superconducting Bulks for Magnetic Screening in Electrical Machines for Aircraft Applications[END_REF]. Therefore, it is important for the model to approach the bulk magnetic behaviour in order to have a more accurate idea of the performance of a flux modulation machine. The most common method to finely simulate this behaviour would be a time-dependant H-formulation considering the E-J law [START_REF] Ainslie | Modelling of bulk superconductor magnetization[END_REF]- [START_REF] Ainslie | Modelling and comparison of trapped fields in (RE)BCO bulk superconductors for activation using pulsed field magnetization[END_REF]. However, the computation in 3-D takes several minutes making it difficult to integrate into a sizing study. Alternatively, the SAM can be improved using the Bean model [START_REF] Bean | Magnetization of High-Field Superconductors[END_REF] to account for the magnetic penetration in the bulk.

𝑑 = 𝐻 𝐽 (81)
In the SAM, complex permeability distributions can be considered. Hence, the penetration of the magnetic field is considered by a linear decline of the permeability with a penetration depth calculated from (81). The value of 𝐻 is computed using the coil axisymmetric model [START_REF] Colle | Analytical Model for the Magnetic Field Distribution in a Flux Modulation Superconducting Machine[END_REF] giving the following relation:

𝐻 = 𝐵 (𝑅 , 0) 𝜇 (82)
Thus, to include this feature in the SAM, only the matrix 𝝁 has to be changed. To do so, the Fourier series coefficients 𝜇̂ , of the permeability need to be recalculated. Fig. 22 shows the permeability distribution of the region I for an incomplete penetration of the bulk. Since the permeability is symmetric along x and y, only a quarter of the region I needs to be considered for the Fourier series coefficients calculation. The bulk centre is still diamagnetic but the permeability in the three sectors at the edges is expressed as:

In sector (a):

𝜇 (𝑥, 𝑦) = 𝜇 𝑑 𝑦 - 𝐷 2 + 𝑑 (83)
In sector (b):

𝜇 (𝑥, 𝑦) = 𝜇 𝑑 𝑥 - 𝛼𝜏 2 + 𝑑 (84)
In sector (c):

𝜇 (𝑥, 𝑦) = 𝜇 1 - 𝑥 - 𝛼𝜏 2 𝑦 - 𝐷 2 𝑑 ( 85
) 0 0.5 1 1.5 2 2.5 3
Bulk aspect ratio

Care should be taken as Fig. 22 represents only the incomplete penetration case. Indeed, three other cases should be considered:

The penetration is complete in the y-direction only (𝑑 > and 𝑑 < ) and only sectors (a) and (c) are present in the bulk.

The penetration is complete in the x-direction only (𝑑 < and 𝑑 > ) and only sectors (b) and (c) are present in the bulk.

The penetration is complete in both directions (𝑑 > and 𝑑 > ) and the sector (c) recovers the whole bulk.

Once again, 10 𝜇 must be added to the 𝜇̂ , coefficient such that the matrix 𝝁 is invertible. In order to evaluate the relevancy of the improved SAM (ISAM), another reference model is made on COMSOL (HFEM). This model is a 3-D time dependant simulation of the inductor using an H-formulation of the Maxwell equations. The superconducting bulk is modelled through an E-J law with a constant 𝐽 [START_REF] Berger | Influence of Temperature and/or Field Dependences of the $E-J$ Power Law on Trapped Magnetic Field in Bulk YBaCuO[END_REF] and the power law index equal to 20. The coil is supplied with a current ramp up to its nominal value in 60 s and then the current is kept constant for 60 s to reach the steady state at which the torque is computed using the armature model of section 4.

Fig. 23 shows the current density distribution normalized to 𝐽 in the bulk calculated by the HFEM. It should be noted that the symmetries exploited here are identical as in the SFEM, only 1 4

⁄ is therefore simulated. Fig. 24 shows the axial component of the flux density in the air gap obtained with the SAM and ISAM for N = 20 compared to the HFEM for 𝐽 = 1000 𝐴/𝑚𝑚 . It can be observed that the HTS bulk behaviour as predicted by the HFEM is not completely diamagnetic for this 𝐽 . The modulation predicted by the ISAM differs from that of a diamagnetic bulk, though there is still a difference between the penetration calculated by the HFEM and the ISAM. Fig. 24 is not sufficient to properly assess the relevancy of the ISAM as the behaviour still seems to be closer to perfect diamagnetism than to the actual behaviour as calculated by the HFEM. For this purpose, Fig. 25 shows the electromagnetic torque calculated by the HFEM and the ISAM for N = 10 function of the bulk critical current density. The torque value is normalized to the diamagnetic torque as the latter represents the maximum achievable torque. In both cases, a transition is observed from complete penetration at low 𝐽 to close to diamagnetism behaviour at high 𝐽 and both models predict similar behaviours. Yet, the transition does not occur at the same 𝐽 values and there is an important difference between 100 𝐴/𝑚𝑚 and 1000 𝐴/𝑚𝑚 which leads to an important error on the torque. Above the latter threshold, both curves are converging to the diamagnetic behaviour. Thus, the application of the Bean model in the ISAM allows to predict the global trend in magnetic behavior but cannot accurately calculate the torque for bulk with a small 𝐽 . This can be explained by the geometrical assumptions of the Bean model which can be applied for infinitely long bulks. In the present case, the bulk thickness in the z-direction is the smallest dimension of the bulk, hence the end effects are very important and cannot be accounted for as the convolution method only allows the permeability to be variable in the x and y directions. Indeed, as shown on Fig. 23, altough the current density is located at the bulk edges as predicted by the Bean model, the latter predict that its value is either 𝐽 or 0 which is not exactly the case. Additionnaly, the currents also penetrate the bulk through the top and bottom surfaces of the bulk (i.e. surfaces along the (𝑟𝜃) plane) which is not taken into account in the HFEM. Nevertheless, the ISAM model is still more accurate than the SAM model and approaches the HFEM results for high 𝐽 . It should be noted that the materials most likely to be used for magnetic shielding are RE-Ba-Cu-O which have 𝐽 above 3000 𝐴/𝑚𝑚 at temperatures close to liquid hydrogen [START_REF]Relationship between local microstructure and superconducting properties of commercial YBa2Cu3O7-$\updelta$ bulk[END_REF] which places them within the area of validity of the ISAM. Lastly, the ISAM requires around a second of computation time for 𝑁 = 10 while the HFEM requires between 15 minutes and 3 hours to be solved depending on the 𝐽 value which makes it impossible to be integrated in an optimisation routine. Moreover, this high computationnal cost make is difficult to use a thin mesh which explain the rough aspect of the current density map of Fig. 23.

Therefore, the ISAM model appears to be the best compromise between precision and calculation time but must be used carefully as it is suitable only for bulks with a behaviour close to the diamagnetism.

VII. CONCLUSION

In this paper, a 3-D semi-analytical model of a superconducting axial flux modulation machine has been developed to calculate the torque and the eddy-current losses. The inductor model has been made using separation of variables method, Fourier series analysis and the application of the convolution theorem. This last tool allowed to consider domains with a variable permeability in two directions and enabled a resolution in 3-D using only 3 domains and therefore a drastically reduced computation time. The torque was derived from the electromotive force whereas the eddy-current losses were calculated from well known analytical formulas and the inductor flux density distribution. Eventually, as the variable permeability domains permit considering complex permeability distribution, the model has been improved using the Bean model to consider the magnetic penetration in the HTS bulk.

Comparisons with the FEM simulations have shown the ability of the proposed model to quickly provide a precise result when the bulks are considered diamagnetic. Indeed, for 10 harmonics, the torque and the eddy-current losses of the test machine can be predicted with less than 5 % error for the torque and 10 % error for the losses in less than 2 seconds. When the critical current density of the bulk is considered, the model has shown its ability to predict the trend in magnetic behaviour but does not accurately simulate bulks with low current densities. Thus, the validity area of this model still only comprises bulks with high critical current density and appears to be the best compromise between computation time and precision. In any case, the performances of the magnetic screens is the cornerstone of high specific power flux modulation machines. Thus, for aicraft applications, only large HTS bulks with large and homogeneous critical current density are relevant and likely to be used.

The model presented in this article is a powerful tool for optimisation and sizing processes. It allows efficient design for prototypes as well as extrapolations and theoretical studies on the potential of this machine technology for aeronautical applications. Additionnaly, the modelling method in 3-D can be adapted to other kind of axial flux machines with region with a variable permeability such as axial flux permanent magnet machines or switch reluctance couplings.

Fig. 1 .

 1 Fig. 1. Exploded view of the active elements in a partially superconducting axial flux modulation machine.

Fig. 2 .

 2 Fig. 2. Geometrical parameters of the inductor modelled parts according to the symmetries.

Fig. 3 .

 3 Fig. 3. Representation of the equivalent bulk geometry in the cartesian frame.

Fig. 4 .

 4 Fig. 4. Sectional view of the cartesian problem in the (xz) plane with the different regions and boundary conditions.

Fig. 5 .

 5 Fig. 5. Sectional view of the cartesian problem in the (yz) plane with the different regions and boundary conditions.

Fig. 6 .

 6 Fig. 6. Magnetic permeability distribution in the region I over the (xy) plane. Thanks to the Fourier series representation, only the z-dependant coefficients of the scalar potential 𝜑 , (𝑧), the magnetic field components 𝐻 , (𝑧), 𝐻 , (𝑧), 𝐻 , (𝑧) and the flux density components 𝐵 , (𝑧), 𝐵 , (𝑧), 𝐵 , (𝑧)

  cosh (𝑫 (𝑍 -𝑍 )) sinh (𝑫 (𝑍 -𝑍 )) -cosh (𝑫 (𝑍 -𝑍 )) sinh (𝑫 (𝑍 -𝑍 ))

  the function 𝑑𝑖𝑎𝑔(𝑿) returns a column vector composed of the elements of the diagonal of the matrix 𝑿 and 𝑿 is the transpose of matrix 𝑿. Eventually, x-component and the z-component of the flux modulation function are the ratio of 𝐵 (𝑥, 𝑦, 𝑧) and 𝐵 (𝑥, 𝑦, 𝑧) respectively and the magnetic flux density generated by the magnet without the bulk B0: 𝐹𝑀 (𝑥, 𝑦, 𝑧) = 𝐵 (𝑥, 𝑦, 𝑧) 𝐵 (67) 𝐹𝑀 (𝑥, 𝑦, 𝑧) = 𝐵 (𝑥, 𝑦, 𝑧) 𝐵

Fig. 7 .

 7 Fig. 7. Geometrical parameter of an armature coil.

Fig. 8 .

 8 Fig. 8. Function 𝐾 (𝑟) accounting for the coil head influence on the flux linkage.

Fig. 9 .

 9 Fig. 9. Representation of the two main arrangements of the armature in relation to the bulks. (a) The bulks cover only the purely radial part of the armature giving the relations: 𝑅 = 𝑅 and 𝑅 = 𝑅 . (b) The bulks also cover the coil heads giving the relations: 𝑅 = 𝑅 + 𝑊 and 𝑅 = 𝑅 -𝑊 .

Fig. 10 .

 10 Fig. 10. Geometry and mesh of the CFEM (a) and SFEM (b). The grey regions are not simulated thanks to the symmetries.

Fig. 11 .

 11 Fig. 11. Axial component of the inductor flux density at the average radius, 3 mm from the bulk surface (𝑟 = 𝑅 and 𝑧 = 𝑍 + 3 𝑚𝑚).

Fig. 13 .

 13 Fig. 13. Axial component of the inductor flux density around the outer bulk radius, 3 mm from the bulk surface (𝑟 = 𝑅 -10 𝑚𝑚 and 𝑧 = 𝑍 + 3 𝑚𝑚).

Fig. 12 .

 12 Fig. 12. Azimuthal component of the inductor flux density at the average radius, 3 mm from the bulk surface (𝑟 = 𝑅 and 𝑧 = 𝑍 + 3 𝑚𝑚).

Fig. 14 .

 14 Fig. 14. Axial component of the inductor flux density around the inner bulk radius, 3 mm from the bulk surface (𝑟 = 𝑅 + 10 𝑚𝑚 and 𝑧 = 𝑍 + 3 𝑚𝑚).

Fig. 15 .

 15 Fig. 15. Flux density norm distribution over the winding region from the CFEM model

Fig. 16 .

 16 Fig. 16. Computational time of the different models versus the number of harmonics.

Fig. 17 .

 17 Fig. 17. Relative error in absolute value on the torque of the SAM compared to the SFEM and the CFEM.

Fig. 18 .

 18 Fig. 18. Normalized eddy current losses calculated by the SAM and the SFEM versus the number of harmonics.



  𝑅 in the range [100 mm : 250 mm].  The ratio 𝑅 𝑅 in the range [0.45 : 0.7].  𝛼 in the range [0.3 : 0.7].

Fig. 19 .

 19 Fig. 19. Relative error on the torque of the SAM with N = 10 compared to the SFEM for the explored parameters.

Fig. 20 .

 20 Fig. 20. Relative error on the eddy current losses of the SAM with N = 10 compared to the SFEM for the explored parameters.

Fig. 21 .

 21 Fig. 21. Implementation principle of the Bean model in the SAM model. The penetration of the magnetic field is modelled as a variation of the permeability instead of a current density.

Fig. 21

 21 Fig. 21 presents the implementation principle of the Bean model in the SAM model. The Bean model considers a linear penetration of the magnetic field at the edges of the superconductor. The penetration depth of the magnetic field 𝑑 is linked to the bulk critical current 𝐽 and the applied magnetic field 𝐻 through the following relation:

Fig. 22 .

 22 Fig. 22. Magnetic permeability distribution in a quarter of the region I for an incomplete penetration.

Fig. 23 .

 23 Fig. 23. Normalized current distribution in 1 4⁄ the bulk at its centre calculated by the HFEM for 𝐽 = 1000 𝐴/𝑚𝑚 : (a) corresponds to the (𝑟𝑧) cross section highlighted in red; (b) corresponds to the (𝑟𝜃) cross section highlighted in cyan.

Fig. 24 .

 24 Fig. 24. Axial component of the inductor flux density at the average radius, 3 mm from the bulk surface (𝑟 = 𝑅 and 𝑧 = 𝑍 + 3 𝑚𝑚) for 𝐽 = 1000 𝐴/𝑚𝑚 and N = 20.

Fig. 25 .

 25 Fig. 25. Electromagnetic torque versus the bulk critical current density for the ISAM (N = 10) and the HFEM. The torque value is normalized to the diamagnetic torque.

TABLE . 1

 . . PARAMETERS OF THE STUDIED MACHINE.

	Symbol Description	Value
	𝑅	Bulk outer radius	148 mm
	𝑅	Bulk inner radius	85 mm
	𝑅	Bulk average radius	116.5 mm
	𝑅	HTS coil inner radius	176 mm
	𝑅	Outer radius of the armature coil inner turn	148 mm
	𝑅	Inner radius of the armature coil inner turn	85 mm
	𝑝	Number of pole pairs	5
	𝛼	Bulk pole-arc to pole-pitch ratio	0.5
	𝑊	Armature coil thickness	18 mm
	𝐿	Armature coil length	20 mm
	𝑊	HTS coil thickness	8.8 mm
	𝑍	HTS coil half length	53.5 mm
	𝑍	Bulk half thickness	10 mm
	𝐸	Air-gap length	3 mm
	𝑑	Winding strand diameter	0.25 mm
	𝐾	Winding filling factor	0.7
	𝑁 𝐼	Ampere-turns per phase	60.5 kA
		HTS coil current density	613 A/mm 2
	𝛺	Mechanical rotational speed	5000 rpm

TABLE . 2

 . . TORQUE CALCULATED BY EACH MODEL FOR THE CASES WHERE THE COIL HEADS ARE COVERED AND UNCOVERED. (A) NORMAL SFEM. (B) SFEM MODEL WITH THE FLUX LINKAGE COMPUTED ONLY ON THE PURELY RADIAL PART OF THE WINDINGS.

	Model	Coil heads uncovered	Coil heads covered
	CFEM	1366 Nm	2206 Nm
	SAM N = 10	1301 Nm	2050 Nm
	SFEM (a)	1322 Nm	2114 Nm
	SFEM (b)	1254 Nm	1711 Nm

Computational time (s)

ACKNOWLEDGEMENTS

The authors would like to thank the Direction Générale de l'Armement (DGA), the Agence de l'Innovation de Défense (AID), and the Agence Nationale de la Recherche (ANR). This work was supported by the project FROST under Grant ANR-19-ASMA-0001k.

Normalized torque