Houda Ferradi

Rémi Géraud

Diana Maimut

David Naccache

Hang Zhou

Backtracking-Assisted Multiplication

This paper describes a new multiplication algorithm, particularly suited to lightweight microprocessors when one of the operands is known in advance. The method uses backtracking to find a multiplicationfriendly encoding of the operand known in advance. A 68hc05 microprocessor implementation shows that the new algorithm indeed yields a twofold speed improvement over classical multiplication for 128-byte numbers.

Introduction

A number of applications require performing long multiplications in performance-restricted environments. Indeed, low-end devices such as the 68hc05 or the 80c51 microprocessors have a very limited instructionset, very limited memory, and operations such as multiplication are rather slow: a mul instruction typically claims 10 to 20 cycles.

General multiplication has been studied extensively, and there exist algorithms with very good asymptotic complexity such as the Schönhage-Strassen algorithm [START_REF] Schönhage | Schnelle Multiplikation grosser Zahlen[END_REF] which runs in time O(n log n log log n) or the more recent Fürer algorithm [START_REF] Fürer | Faster integer multiplication[END_REF], some variants of which achieve the slightly better O(2 3 log * n n log n) complexity [START_REF] Harvey | Even faster integer multiplication[END_REF]. Such algorithms are interesting when dealing with extremely large integers, where these asymptotics prove faster than more naive approaches.

In many cryptographic contexts however, multiplication is performed between a variable and a pre-determined constant:

-During Diffie-Hellman key exchange [START_REF] Diffie | New directions in cryptography[END_REF] or El-Gamal [START_REF] El Gamal | A public key cryptosystem and a signature scheme based on discrete logarithms[END_REF] a constant g must be repeatedly multiplied by itself to compute g x mod p. -The essential computational effort of a Fiat-Shamir prover [START_REF] Feige | Zero knowledge proofs of identity[END_REF][START_REF] Feige | Zero-knowledge proofs of identity[END_REF] is the multiplication of a subset of fixed keys (denoted s i in [START_REF] Feige | Zero knowledge proofs of identity[END_REF]). -A number of modular reduction algorithms use as a building-block multiplications (in N) by a constant depending on the modulus. This is for instance the case of Barrett's algorithm [START_REF] Barrett | Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor[END_REF] or Montgomery's algorithm [START_REF] Montgomery | Modular multiplication without trial division[END_REF].

The main strategy to exploit the fact that one operand is constant consists in finding a decomposition of the multiplication into simpler operations (additions, subtractions, bitshifts) that are hardware-friendly [START_REF] Bernstein | Multiplication by integer constants[END_REF]. The problem of finding the decomposition with the least number of operations is known as "single constant multiplication" (SCM) problem. SCM ∈ NP-complete as shown in [START_REF] Cappello | Some complexity issues in digital signal processing[END_REF], even if fairly good approaches exist [START_REF] Avižienis | Signed-digit number representations for fast parallel arithmetic[END_REF][START_REF] Dempster | Constant integer multiplication using minimum adders[END_REF][START_REF] Dempster | Use of minimum-adder multiplier blocks in FIR digital filters[END_REF][START_REF] Wu | Closed-form expression for the average weight of signed-digit representations[END_REF] for small numbers. For larger numbers, performance is unsatisfactory unless the constant operand has a predetermined format allowing for ad hoc simplifications.

In this paper, we propose a completely different approach: the constant operand is encoded in a computation-friendly way, which makes multiplication faster. This encoding is based on linear relationships detected amongst the constant's digits (or, more generally, subwords), and can be performed offline in a reasonable time for 1024-bit numbers and 8-bit microprocessors. We use a graph-based backtracking algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF] to discover these linear relationships, using recursion to keep the encoder as short and simple as possible.

Multiplication Algorithms

We now provide a short overview of popular multiplication methods. This summary will serve as a baseline to evaluate the new algorithm's performance.

Multiplication algorithms usually fall in two broad categories: general divide-and-conquer algorithms such as Toom-Cook [START_REF] Cook | On the minimum computation time of functions[END_REF][START_REF] Toom | The Complexity of a Scheme of Functional Elements Realizing the Multiplication of Integers[END_REF] and Karatsuba [START_REF] Karastuba | Multiplication of many-digital numbers by automatic computers[END_REF]; and the generation of integer multiplications by compilers, where one of the arguments is statically known. We are interested in the case where small-scale optimizations such as Bernstein's [START_REF] Bernstein | Multiplication by integer constants[END_REF] are impractical, but general purpose multiplication algorithms à la Toom-Cook are not yet interesting.

Throughout the paper we will assume unsigned integers, and denote by w the word size (typically, w = 8), a i , b i and r i the binary digits of a, b and r respectively:

a = n-1 i=0 2 wi a i , b = n-1 i=0 2 wi b i , and r = a × b = 2n-1 i=0 2 wi r i .

Textbook Multiplication

A direct way to implement long multiplication consists in extending textbook multiplication to several words. This is often done by using a mad1 routine.

A mad routine takes as input four n-bit words {x, y, c, ρ}, and returns the two n-bit words c , ρ such that

2 n c + ρ = x × y + c + ρ. We write {c , ρ } ← mad(x, y, c, ρ).
If such a routine is available then multiplication can be performed in n 2 mad calls using Algorithm 1. The MIRACL big number library [START_REF]Certivox: The MIRACL big number library[END_REF] provides such a functionality. This approach is unsatisfactory: it performs more computation than often needed. Assuming a constant-time mad instruction, Algorithm 1 runs in time O(n 2 ).

Algorithm 1: Mad-based computation of r = a × b. Input: a, b ∈ N. Output: r ∈ N such that r = a × b. 1 for i ← 0 to 2n -1 do 2 ri ← 0 3 end for 4 for i ← 0 to n -1 do

Karatsuba's Algorithm

Karatsuba [START_REF] Karastuba | Multiplication of many-digital numbers by automatic computers[END_REF] proposed an ingenious divide-and-conquer multiplication algorithm, where the operands a and b are split as follows:

r = a × b = (2 L ā + a) × (2 L b + b),
where typically L = nw/2. Instead of computing a multiplication between long integers, Karatsuba performs multiplications between shorter integers, and (virtually costless) multiplication by powers of 2. Karatsuba's algorithm is described in Algorithm 2.

Algorithm 2: Karatsuba's algorithm to compute r = a × b.

Input: a, b ∈ Z. Output: r ∈ Z such that r = a × b. 1 u = ā × b 2 v = a × b 3 w = (ā + a)( b + b) -u -v 4 r = 2 2L × u + 2 L × w + v 5 return r
This approach is much faster than naive multiplication -on which it still relies for multiplication between short integers -and runs2 in Θ(n log 2 3 ).

Bernstein's Multiplication Algorithm

When one of the operands is constant, different ways to optimize multiplication exist. Bernstein [START_REF] Bernstein | Multiplication by integer constants[END_REF] provides a branch-and-bound algorithm based on a cost function.

The minimal cost, and an associated sequence, are found by exploring a tree, possibly using memoization to avoid redundant searches. More elaborate pruning heuristics exist to further speedup searching. The minimal cost path produces a list of operations which provide the result of multiplication.

Because of its exponential complexity, Bernstein's algorithm is quickly overwhelmed when dealing with large integers. It is however often implemented by compilers for short (32 to 64-bit) constants.

The Proposed Algorithm

Intuitive Idea

The proposed algorithm works with an alternative representation of the constant operand a. Namely, we wish to express some a i as a linear com-bination of other a j s with small coefficients. It is then easy to reconstruct the whole multiplication b × a from the values of the b × a j only.

The more linear combinations we can find, the less multiplications we need to perform. Our algorithm therefore tries to find the longest sequence of linear relationships between the digits of a. We call this sequence's length the coverage of a.

Yet another performance parameter is the number of registers used by the multiplier. Ideally at any point in time two registers holding intermediate values should be used. This is not always possible and depends on the digits of a.

As an example, consider the set of relations of Table 1. All words are expressed ultimately in terms of the values of a 3 and a 7 . In Table 1, we express a as a subset of words A ∈ {a 0 , . . . , a n-1 } and build a sparse table U where U i,j ∈ {-1, 0, 1, 2, =}, which encodes linear relationships between individual words. During multiplication, U describes how the different a i can be derived from each other. Hence it suffices to compute b × a 3 and b × a 7 to infer all other b × a i by long integer additions. Note that the algorithm only needs to allocate three (n + 1)-word registers reg1, reg2 and reg3 to store intermediate results.

The values allowed in U can easily be extended to include more complex relationships (larger coefficients, more variables, etc.) but this immediately impacts the algorithm's performance. Indeed, the corresponding search graph has correspondingly many more branches at each node.

Operations can be performed without overflowing (i.e. so that results fit in a word), or modulo the word size. In the latter case, it is necessary to subtract b w from the result, where w is the word size, to obtain the correct result. This incurs some additional cost.

Backtracking Algorithm

Algorithm 3: macro Step(u, w).

1 (p d+1,0 , p d+1,1 ) ← (u, w) 2 Backtrack(d + 1)
Algorithm 4: macro EncodeDep(c, opcode).

1 if c < 256 then 2 if vc = False then 3 (vc, p d,2 , p d,3 ) ← (True, c, opcode) 4
Step(a, b)

5

Step(a, c)

6
Step(b, c) Linear combinations amongst words of a are found by backtracking [START_REF] Knuth | The Art of Computer Programming[END_REF], the of which is given in Algorithm 5. Our implementation focuses on linear dependencies amongst 8-bit words, as our main recommendation for applying the proposed multiplication algorithm is exactly an 8-bit microprocessor. We take advantage of recursion and macro expansion (see Algorithms 3 to 5) to achieve a more compact code. In this implementation, p encodes the current depth's three registers of Table 1 as well as the current operation. With suitable listing, Algorithm 6 outputs a set of values being related, along with the corresponding relation. The dependencies that we take into account in our C code (given in Appendix A) don't go beyond depth 2. Thus, the corresponding operations are C = {+, -, ×2}. We also add these operations performed modulo 256, to obtain more solutions. The alternative to this approach is to consider a bigger depth, which naturally leads to more possibilities.

Our program takes as an input an integer p that represents the percentage of a being covered (i.e. the coverage is p/100 times the length of the a). In a typical lightweight scenario, a 128-byte number is involved in the multiplication process. Our software attempts to backtrack over a coverage-related number of values out of 256. It follows immediately that at most a 50% coverage would be required for performing such a multiplication (as byte collisions are likely to happen).

The program takes as parameter the list of bytes of a. If some bytes appear multiple times, it is not necessary to re-generate each of them individually: generation is performed once, and the value is cloned and dispatched where needed.

Note that if precomputation takes too long, the list of a i can be partitioned into several sub-lists on which backtrackings are run This would entail as many initial multiplications by the online multiplier but still yield appreciable speed-ups. With the encoding of a generated by Algorithm 6, it is now possible to implement multiplication efficiently.

Multiplication Algorithm Algorithm 7: Virtual Machine

Input: b, instr = (opcode, i, j, t, p) k , R Output: r 1 r ← 0 2 foreach (i, v) ∈ R do 3 reg[i] ← v × b 4 PlaceAt(v, reg[i])
To that end we make use of a specific-purpose multiplication virtual machine (VM) described in Algorithm 7. The VM is provided with instructions of the form opcode t, i, j, p that are extracted offline from U . Here, opcode is the operation to perform, i and j are the indices of the operands, t is the index of the result, and p ← w × t is the position in r where to place the result, w being the word size. The value of p is pre-computed offline to allow for a more efficient implementation.

We store the result in a 2n-byte register initialized with zero. We also make use of a long addition procedure PlaceAt(p, i) which "places" the contents of the (n + 1)-byte register reg[i] at position p in r. PlaceAt performs the addition of register reg[i] starting from an offset p in r, propagating the carry as needed.

Finally, we assume that the list R = (i, v) k of root nodes (position and value) of U is provided.

After executing all the operations described in U , Algorithm 7 has computed r ← a × b.

Remark 1 (Karatsuba Multiplication).

Using the notations of Algorithm 2 one can see that in settings where a is a constant, the numbers u, v, w all result from the multiplication of b, b and b + b (which are variable) by ā, a and ā + a (which are constant).

Hence our approach can independently be combined with Karatsuba's algorithm to yield further improvements.

Performance

The algorithm has an offline step (backtracking) and an online step (multiplication), which are implemented on different devices.

The offline step is naturally the longest; its performance is heavily dependent on the digit combination operations allowed and on how many numbers are being dealt with. More precisely, results are near-instant when dealing with 64 individual bytes and operations {+, -, ×2}. It takes much longer if operations modulo 256 are considered as well, but this gives a better coverage of a, hence better online results. That being said, modulo 256 operations are slightly less efficient than operations over the integers ( 1.5 more costly), since they require a subtraction of b afterwards.

Table 2 provides comparative performance data for a multiplication by the processed constant π2 1024 . Backtracking this constant took 85 days on an Altix UV1000 cluster. As a final remark, note that one can also reverse the idea and generate a key by which multiplication is easy. This can be done by progressively picking VM operations until an operand (key) with sufficient entropy is obtained. While this is not equivalent to randomly selecting keys, the authors conjecture that, in practice, the existence of linear relations 3 between key bytes should not significantly weaken public-key implementations.

5 c ← 0 6 for j ← 0 to n -1 do 7 {c

 67 , ri+j} ← mad(ai, bj, c, ri+j)

Algorithm 5 :

 5 macro Backtrack(d).1 if d > d max then 2 (d max , p max ) ← (d, p) 3 end if 4 (a, b) ← (p d,0 , p d,1 ) 5 for opcode ∈ C do 6 EncodeDep(opcode(a, b), opcode) 7 end for

Algorithm 6 : 14 Backtrack(0) 15 ( 23 (

 6141523 Main Backtracking Program. Input: A = N -1 i=0 256 i Ai. Output: Ui,j 1 // Initialization 2 for i = 0 to 255 do 3 vi ← True 4 end for 5 for i = 0 to N -1 do 6 vA i ← False 7 end for 8 d max ← -1 9 // Backtracking 10 for i = 0 to 255 do 11 for j = i + 1 to 255 do 12 if vi = vj = False then 13 (p0,0, p0,1, vi, vj) ← (i, j, True, True) vi, vj) ← (False, False) 19 // U -matrix reconstruction 20 Ui,j ← 0 21 for i = 0 to 255 do 22 (in1, in2, out, opcode) ← p max i Ui,in 1 , Ui,in 2 , Ui,out) ← (1, 1, opcode) 24 end for 25 return Ui,j

Table 1 .

 1 An example showing how linear relationships between individual words are encoded and interpreted.

	step a0 a1 a2 a3 a4 a5 a6 a7 a8 a9	meaning	reg1 reg2 reg3
	0	1	=	1	a5 ← a3 + a7	a5	a3	a7
	1		= 2		a4 ← a5 + a5	a5	a4	a7
	2	=	-1	1	a0 ← a7 -a4	a0	a4	a7
	3		2	= a9 ← a4 + a4	a0	a4	a9
	4	1 =		-1 a1 ← a0 -a9	a0	a1	a9
	5	1		= 1 a8 ← a1 + a9	a8	a1	a9
	6	1 =		1	a2 ← a1 + a8	a8	a1	a2
	7	2	=		a6 ← a2 + a2	a8	a6	a2

9 end foreach 10 return r

	5 end foreach
	6 foreach (opcode, i, j, t, p) ∈ instr do
	7	reg[t] ← opcode(reg[i], reg[j])
	8	PlaceAt(p, t)

Table 2 .

 2 Performance on a 68hc05 clocked at 5 MHz

		Time	ram	Code Size
	Usual Algorithm 188 ms 395 bytes 1.1 kilobytes
	New Algorithm	72 ms 663 bytes 1.7 kilobytes

An acronym standing for "Multiply Add Divide"

When repeated recursively.

A Source Code

# include < stdio .h > # include < string .h > # include < math .h > # include < stdlib .h > int False = 0; int True = 255;


Antoine Marteau

Innocent Niyonzima

Gerard Meunier

Janne Ruuskanen

Nicolas Galopin

Paavo Rasilo

Olivier Chadebec

Gérard Meunier

Magnetic Field Upscaling and B-Conforming Magnetoquasistatic Multiscale Formulation

published or not. The documents may come

L'archive ouverte

This paper proposes two novel methods for upscaling the macroscopic magnetic field strength from the local solutions of Bconforming magnetoquasistatic multiscale problems. Unlike the volume average method classically used, these methods yield accurate values of the macroscale magnetic field for problems with strong locally-confined eddy currents which enables B-conforming multiscale formulations of eddy currents problems at higher frequencies.

Index Terms-Eddy currents, heterogeneous multiscale method, homogenization, nonhomogeneous media

I. INTRODUCTION

M derived in Section III, and the validation on a 2D eddy currents problem is presented in Section IV.

II. MULTISCALE MODELING

Let σ and ν be the electrical conductivity and magnetic reluctivity. We are interested in solving the following Maxwell's equations [START_REF] Feddi | Homogenization technique for Maxwell equations in periodic structures[END_REF] on an open domain Ω ⊂ R d in the magnetoquasistatic regime: curl h = j, div b = 0, curl e = -∂ t b, j = σe, h = νb, [START_REF] Feddi | Homogenization technique for Maxwell equations in periodic structures[END_REF] where h, b, j and e are the magnetic field strength, magnetic flux density, electric current density and electric field, respectively. The conducting domain is denoted by Ω c , and the nonconducting one by Ω N = Ω \ Ω c . The dimension d is 2 or 3.

A. Fine-scale model on periodic geometry

Let us consider a periodic medium denoted by Ω ϵ included in Ω, that is a tiling obtained by translation of a periodic unit cell Y = Y c ∪ Y N , where Y c and Y N are its conducting and nonconducting regions. Without loss of generality, the cell Y is chosen as a rectangular hexahedron. In this paper, we make the assumption that the conducting region of the cell Y c is simply connected and does not intersect the boundaries of the periodic cell, i.e. ∂Y ∩ Y c = ∅, which is tantamount to an insulated periodic cell.

We use a usual B-conforming formulation [START_REF] Kameari | Calculation of Transient 3D Eddy Current Using Edge-Elements[END_REF]:

         curl ν curl a + σ∂ t a = 0 in Ω ϵ , (2a) a × n = a × n on ∂Ω ϵ , ( 2b
) curl ν 0 curl a = j s in Ω \ Ω ϵ , (2c) a × n = 0 on ∂Ω, (2d
) where the modified magnetic vector potentials a verifies b = curl a and e = -∂ t a, and has to be gauged in non conducting parts of Ω. Here a known current source j s and magnetic reluctivity of the vacuum ν 0 are assumed in the non homogenized domain Ω \ Ω ϵ . ULTISCALE methods have been extensively studied in the electromagnetic community thanks to their easy parallelization and ability to speed up numerical simulations of problems involving heterogeneous materials such as laminated cores, stranded inductors and soft magnetic composites. These methods make possible the simulation of 3D nonlinear eddy currents problems with composite materials, unlike classical methods such as the finite element method (FEM) which are computationally too expensive for these problems.

Multiscale methods generally necessitate computation of macroscopic quantities from finer scale quantities, also called upscaling. Previous works [START_REF] Feddi | Homogenization technique for Maxwell equations in periodic structures[END_REF] - [START_REF] Sabariego | Nonlinear Computational Homogenization Method for the Evaluation of Eddy Currents in Soft Magnetic Composites[END_REF] use B-conforming formulation that necessitate the upscaling of magnetic field strength, which is usually done using the volume average of the finer scale magnetic field h. In the early work from El Feddi et. al. [START_REF] Feddi | Homogenization technique for Maxwell equations in periodic structures[END_REF], upscaling of magnetic field using volume average is done at very high frequency. In Bottauscio et. al. [START_REF] Bottauscio | Nonlinear Homogenization Technique for Saturable Soft Magnetic Composites[END_REF] and Niyonzima et. al. [START_REF] Sabariego | Nonlinear Computational Homogenization Method for the Evaluation of Eddy Currents in Soft Magnetic Composites[END_REF], the multiscale solution is validated at low and medium frequencies. But Meunier et. al. [START_REF] Meunier | Homogenization for Periodical Electromagnetic Structure: Which Formulation?[END_REF] warned that the upscaling of the magnetic and the electric fields h and e c an unfortunately be non trivial for eddy currents problems. Indeed, a simple volume average of the finescale magnetic field is not compatible with the Maxwell-Ampère law at macroscopic scale in presence of strong locally confined eddy currents. In their paper, the authors proposed that Bconforming formulation should be avoided to prevent upscaling h.

This paper proposes two novel methods for the upscaling of the magnetic field strength that are always valid. These methods only depend on the solutions of the mesoscale problem and the topology of its conducting region. The paper is outlined as follow, Section II recalls the multiscale model for the Bconforming formulation. The proposed upscaling methods are

B. The two scale model

In our multiscale approach, the reference problem (1) is replaced by the macroscale problem and many mesoscale problems both governed by equations similar to [START_REF] Bottauscio | Nonlinear Homogenization Technique for Saturable Soft Magnetic Composites[END_REF] and solved in an iterative scheme. Mesoscale fields are defined on the periodic cell Y and admit the following decomposition [START_REF] Meunier | Homogenization for Periodical Electromagnetic Structure: Which Formulation?[END_REF]:

h(y) = H + J × y d -1 + h c (y), b(y) = B + b c (y), e(y) = E -∂ t B × y d -1 + e c (y), j(y) = J + j c (y), (3)
where capital letters denote the macroscale fields assumed constant on Y , the subscript • c is used to denote correction fields assumed periodic on the cell Y and y represents the local coordinate vector in the cell. Periodicity relates to tangential components of e c and h c and to normal components of b c and j c on ∂Y . Here J is null as Y is insulated.

The mesoscale B-conforming formulation is derived using ( 1) and (3) as follow:

curl ν (B + curl a c ) + σ ∂ t a c + B × y d -1 = 0, (4)
where the unknown potential a c has a periodic tangential component across ∂Y , verifies b c = curl a c and has to be gauged in Y N . The magnetic flux density, magnetic field strength and electric field are

b = B + curl a c , h = νb and e = -∂ t a c -∂ t B × y d-1 .
To define the macroscopic problem, let Ω H be the homogenized domain replacing Ω ϵ . The macroscopic vector potential A is the solution of (2b)-(2d) in Ω \ Ω H , but the equation in Ω H must be defined for macroscopic fields. Since the current density is zero at macroscale, Ω H is not conducting. Hence macroscopic Maxwell-Ampère's law reads

curl H = 0 in Ω H , ( 5
)
where H is the upscaled magnetic field which will be computed from the mesoscale solution obtained by imposing the source B to the cell problem. In general, H is a nonlinear function of B, and ( 5) can be solved iteratively, e.g. with :

curl

H(B k ) + ∂H ∂B k (B k+1 -B k ) = 0 in Ω H , (6)
where

B k+1 = curl (A k+1
). The Jacobian ∂H ∂B k can be computed using finite differences, its columns are given by:

∂H ∂B i k = 1 δB H(B k + δBe i ) -H(B k ) for i = 1, . . . d ,
where e i are the standard basis vectors, δB is a scalar perturbation and k is the nonlinear iteration index. This Jacobian computation requires solving d + 1 problems on each cell. The definitions of macroscopic fields are needed to define the decomposition [START_REF] Sabariego | Nonlinear Computational Homogenization Method for the Evaluation of Eddy Currents in Soft Magnetic Composites[END_REF]. The usual definition of macroscopic fields uses the volume average <> Y of the mesoscopic ones:

U =< u > Y = 1 |Y | Y u for a field u on Y,
where |Y | is the volume of Y . This definition always works for B and J , i.e.

B = < b > Y , J = < j > Y , (7)
because the surface average of the flux of b and j through the faces of ∂Y respectively sum up to volume average of b and j in Y (because div b = 0 and div j = 0), meaning that definition ( 7) is compatible with the continuity of B • n and J • n on ∂Ω H . Combining (3) and ( 7)

yields < b c > Y = 0, which is verified in (4) because Y b c = Y curl a c = ∂Y a c × n = 0,
where the integral on ∂Y vanishes due to anti-periodicity of a c × n.

The magnetic field h can also be upscaled using < h > Y when eddy currents in Y are weak or non existent. Indeed, some previous work could obtain correct results using < h > Y for dynamical problems with eddy currents (see e.g. [START_REF] Feddi | Homogenization technique for Maxwell equations in periodic structures[END_REF], [START_REF] Bottauscio | Nonlinear Homogenization Technique for Saturable Soft Magnetic Composites[END_REF] and [START_REF] Sabariego | Nonlinear Computational Homogenization Method for the Evaluation of Eddy Currents in Soft Magnetic Composites[END_REF]), although the upscaling of h cannot generally be done in this way, as will be shown later. The problem comes from the fact that eddy currents impact the average of h without impacting the average of j, leading the macroscopic Maxwell-Ampère's law to be wrong.

III. MAGNETIC FIELD UPSCALING

Although all the geometries and numerical tests presented in this paper are done with an insulated cell, the material of this section is derived in the general setting of a conducting cell, but still with the assumption that Y c is simply connected. In this section, we demonstrate how to properly upscale H by removing the part of h that generates the local eddy currents in j. More precisely, we define a potential t 0 such that curl t 0 = j e such that j e are the local eddy currents, and such that t 0 = 0 when there are no eddy currents (j e = 0). Then H can be defined as

H =< h -t 0 > Y . ( 8
)
The rest of the current j M = jj e is the flow of the macroscopic current in the cell.

A. A setting to study currents in the cell

The splitting j = j e +j M will be defined with a Helmholtz-Hodge decomposition (HHD) [START_REF] Bhatia | The Helmholtz-Hodge Decomposition-A Survey[END_REF]- [START_REF] Rundong | 3D Hodge Decompositions of Edge-and Face-Based Vector Fields[END_REF]. Let Γ c = ∂Y c be the boundary of the cell conducting domain, Γ out c = Γ c ∩ ∂Y the boundary of Y through which the current can flow and Γ in c = Γ c \ Γ out c the interface between the conducting and non conducting domain in Y . The current density verifies div j = 0, (j • n)n is periodic through Γ out c , and the current cannot flow to the nonconducting domain, i.e. j • n| Γ in c = 0. We denote the space for j by H #,out 0,in (div; Y c ), where 0,in is the usual notation for null trace on Γ in c and where #,out denotes the trace periodicity on Γ out c . It is hard to study j in H #,out 0,in (div; Y c ) because it has mixed boundary conditions, and neither Γ in c nor Γ out c are closed and simply connected when the cell is not insulated. We thus introduce a simpler equivalent setting.

Let Y c be the periodized Y c , that is Y c but with the opposite components of Γ out c identified. The boundary ∂Y c is Γ in c periodized, so that ∂Y c is a closed surface while Γ in c is not. There is a natural one to one correspondence between fields on Y c and the fields on Y c with periodic trace on Γ out c . In particular H #,out 0,in (div; Y c ) is isomorphic to H 0 (div; Y c ). We use symbols j, t for the fields defined on Y c corresponding to the Γ out c -periodic fields j, t defined on Y c .

B. Helmholtz-Hodge decomposition

In general, any vector field j on Y c can be orthogonally decomposed as follows with a HHD (see eq. ( 7) from [7] 1 ):

j = grad φ + curl t0 + jt in Y c ( 9
)
where φ ∈ H(grad; Y c ), t0 ∈ H 0 (curl; Y c ) and jt ∈ H 1 (Y c ) is a tangential harmonic field, i.e. it is curl-and divergencefree and jt •n| ∂Yc = 0. H 1 (Y c ) is the first co-homology space, its dimension β 1 is the finite number of loops that Y c forms, and its elements (e.g. jt ) are harmonic fields looping in its loops [START_REF] Rundong | 3D Hodge Decompositions of Edge-and Face-Based Vector Fields[END_REF]. In addition, the flux of t0 through each connected component of ∂Y c must be 0 to nullify its normal harmonic components.

The general decomposition ( 9) is simpler for j with

grad φ = 0 in Y c . (10)
Indeed, by definition, all three j, jt and curl t0 are divergence free and tangential to ∂Y c . Applying the divergence operator on (9) yields div (grad φ) = 0. Similarly, taking the trace of (9) yields grad φ • n| ∂Yc = 0. Then

∥grad φ∥ 2 = Yc grad φ • grad φ = - Yc div(grad φ) φ + ∂Yc (grad φ • n) φ = 0.
At this point, we can define the decomposition of j on Y c . Let t 0 ∈ H #,out 0,in (curl; Y c ) be the field corresponding to t0 , j M the one corresponding to jt and define j e = curl t 0 . We have the unique orthogonal decomposition

j = j e + j M in Y c . ( 11
)
Moreover, the decomposition verifies < j e > Yc = 0, and

< j M > Y = J (12) (j e = 0 =⇒ < t 0 > Y = 0). ( 13
)
Proof of (12):

Yc j e = Yc curl t 0 = ∂Yc t 0 × n = Γ in c t 0 × n + Γ out c
t 0 × n = 0 because both integrals vanish due to t 0 × n nullity on Γ in c and anti-periodicity on Γ out c . Then 1 The correspondence between [START_REF] Rundong | 3D Hodge Decompositions of Edge-and Face-Based Vector Fields[END_REF] notations and ours is :

Yc = M, ∂Yc = ∂M, j = ω 2 , t0 = α 1 n , jt = h 2 n , curl = d 1
and -grad = δ 3 . As h 2 t and η 2 can both be written as codifferential δ 3 , we merged them with δ 3 β 3 t in our grad φ term. The decomposition ( 9) is still orthogonal in 2D because hn and ht +η are always orthogonal (see comments about Friedrichs decomposition).

< j M > Y = < j > Y -< j e > Y = < j > Y = J .
Proof of (13): when curl t0 = 0 in Y c , there exist ψ0 ∈ H 0 (grad; Y c ) such that t0 = grad ψ0 (because t0 contains no harmonic part). Now coming back in Y c , one gets t 0 = grad ψ 0 where ψ 0 ∈ H #,out 0,in (grad; Y c ). Finally,

Yc t 0 = Yc grad ψ 0 = Γ in c ψ 0 n + Γ out c ψ 0 n = 0,
where both integrals vanish due to the nullity of ψ 0 n on Γ in c and its anti-periodicity on Γ out c . To sum up, the decomposition (11) is orthogonal and unique and its properties have a desirable physical meaning. The "macroscopic current flow" j M flows through ∂Y , averages to J and has zero curl. On the other hand, the eddy currents j e have zero average and contain all "curling" currents (curl j = curl j e ). The potential t 0 is defined to represent the part of h that creates j e . Indeed curl (h -t 0 ) = jj e = j M and the average of t 0 vanishes when eddy currents vanish (13).

C. Macroscopic magnetic field computation

As stated before [START_REF] Geuzaine | Gmsh[END_REF], we propose to define the macroscopic field as H =< ht 0 > Y . This way, the local currents j e no longer impact H, which is consistent with the macroscopic Maxwell-Ampère's law which relates to J but not j e .

The potential t 0 can be extracted from j by solving the projection curl (jcurl t 0 ) = 0 in Y c . When j is the solution of ( 4) and the cell is insulated, the FEM weak form reads:

Yc σ ∂ t a c + B × y d -1 curl t ′ 0 + Yc curl t 0 curl t ′ 0 = 0 (14) for all t ′ 0 , with t 0 , t ′ 0 ∈ H 0 (curl; Y c ).
Here the periodicity of t 0 is not applied because Γ out c = ∅. And nothing is needed to ensure that t 0 cannot contain a normal harmonic part, because Γ in c is connected. A decomposition example is shown Fig. 1. Alternatively, when Y c is insulated, it is possible to compute H in another way. Introducing < u > ∂Y ∥ the tangential average of a vector field u over ∂Y , it is defined by:

< u > ∂Y ∥ = i=x,y,z 1 |∂Y i | ∂Y i u t , (15)
where ∂Y i are faces of ∂Y with normal n = e i and

u t = (n × u) × n. One can prove that < u > ∂Y ∥ = < u > Y whenever u t is Y -periodic and curl u = 0 on Y . In an insulated cell, curl (h -t 0 ) = j M = 0, thus H =< h -t 0 > Y =< h -t 0 > ∂Y ∥ =< h > ∂Y ∥ , (16)
where the last equality is due to t 0 | ∂Y = 0. So the simple boundary integration < h > ∂Y ∥ can be used to upscale H without doing the FEM computation of t 0 . But this alternative is only available for insulated cell, and might lead to more numerical error than using the volume average [START_REF] Geuzaine | Gmsh[END_REF].

IV. VALIDATION The validation is done on a 2D insulated cell of width 100 µm with a conducting disks of radius 40 µm (see Fig. 4). If not mentioned otherwise, we consider f = 100 MHz, ν = ν 0 everywhere, σ = 1.01×10 7 S/m in the disk and σ = 0 outside.

A. Tests on a single cell problem

A first validation is done using the H-conforming formulation of the cell problem with the imposed source H(t) known at every time step t. Figure 2 shows the validation that the proposed methods work, unlike the volume average. Time (ns) To explain why previous use of < h > Y as the upscaling formula were successful, the plots on Fig. 3 show that the error committed is correlated with the percentage of the Joule losses in the total electromagnetic power of the cell. This percentage depends on the physical parameters, geometry, and especially on the frequency. Using H = < h > Y below f = 500 kHz would work fine to stay below 1% of error. In general, simulation of magnetic cores with small Joule losses are unlikely to require our approach.

H (A/m) H < h > Y < h -t0 > Y < h > Y ∥

B. Multiscale model validation

A full 2D multiscale problem is solved with FEM Bconforming formulation at macro and meso scales. The studied periodic geometry Ω ϵ is a square of 10 by 10 cells. A source field is generated by an external inductor coil surrounding Ω ϵ and fed by a sinusoidal current. The coil has a radius of 40 µm and is 670 µm away from Ω ϵ (Fig. 4

top).

Heterogeneous Multiscale Method (HMM) [START_REF] Sabariego | Nonlinear Computational Homogenization Method for the Evaluation of Eddy Currents in Soft Magnetic Composites[END_REF] is used for the multiscale problem, i.e. a cell problem is associated to each element of the macro mesh, which do not have to correspond to the real geometry tiling. The macro and meso problems are solved alternatively. The reference mesh takes into account the real geometry of Ω ϵ . Problems are meshed with Gmsh [START_REF] Geuzaine | Gmsh[END_REF] and solved using GetDP [START_REF] Dular | GetDP[END_REF].

Using < h > ∂Y ∥ , HMM yields the correct value of the total Joule losses over one period with 1.34% of error compared to the reference, versus 90.2% with < h > Y (Fig. 5). V. CONCLUSION This paper presented the ins and outs of magnetic field homogenization on geometries with insulated periodic cells. Firstly by demonstrating the limitation of the volume average and its dependency on the eddy currents power ratio in the total electromagnetic power. Secondly by proposing a general definition of the homogenized magnetic field which is mathematically well defined and comprehensive with physical insights. And finally by validating the proposed definitions on multiscale solutions.

As future perspectives, the method should be validated on a 3D geometry, and with nonzero macroscopic current. It can also be extended to geometries with multi-connected cell conducting domains. Finally, the same approach could be used for the electric field upscaling.

Fig. 1 .

 1 Fig.1. Example of HHD in a 2D insulated conducting sphere at high frequency, from left to right : j = j e (two viewpoints), t 0 , h and ht 0 .

Fig. 2 .

 2 Fig. 2. Comparison of the magnetic field upscaling strategies, over one period of H, which is sinusoidal. The proposed methods < ht 0 > Y (8) and < h > Y ∥ (16) match H, but the volume average < h > Y doesn't.

Fig. 3 .

 3 Fig. 3. Log-log error contour maps in the meso cell exited by a sinusoidal B source, at different frequencies and at different sphere relative permeabilities. Left : percentage of difference between T 0 ∥H(t)-< h(t) > Y ∥dt and T 0 ∥< h(t) > Y ∥dt. Right : percentage of the Joule losses power percentage in the total electromagnetic power.

Fig. 4 .

 4 Fig. 4. Geometry for multiscale test with the two inductor cuts on the sides (top), meshes of corners of Ω ϵ and Ω H (left & middle), one cell mesh (right).

Fig. 5 .

 5 Fig. 5. Joule losses power comparison in the multiscale simulation.


Enrico Di Russo

Lorenzo Rigutti
email: lorenzo.rigutti@univ-rouen.fr

Enrico Di Russo
email: enrico.dirusso@unipd.it

Correlative atom probe tomography and optical spectroscopy: an original gateway to materials science and nanoscale physics

Keywords: Atom Probe Tomography, Photoluminescence, Optical Spectroscopy, Correlative Microscopy, Semiconductors, Heterostructures

in 2018. His research interests included the development of the "photonic atom probe" and the study of dissociation reactions occurring during atom probe tomography (APT). Afterward, he worked as postdoc at CEA-LETI, developing a unique expertise to study optoelectronic devices combining APT with other microscopy techniques such as electron holography and cathodoluminescence. In 2020, he moved toward the University of Padova to develop an innovative method to produce direct band-gap GeSn alloys by pulsed laser melting.

Introduction

The use of ultrafast laser pulses, introduced since 2005, has opened up the possibility of the analysis of non-metallic materials by Atom Probe Tomography (APT), with important feedback in many fields of materials science (physical metallurgy, study of oxides, etc.) and nanoscience (spin valves, tunnel junctions, heterostructures, etc.), where it plays a key role in better understanding the links microstructure-physical properties. At the same time, multi-technique, or correlative approaches have been developed. These allow not only to reach an ever more advanced level in the interpretation of the link between structural and functional characteristics of materials, but also to achieve a deeper understanding of laser-matter interaction phenomena under high field. In this article, we will shortly review the particular case of the correlation of optical spectroscopy techniques such as photoluminescence (PL) spectroscopy with APT. A wide class of materials can be studied in order to understand the relationship between optical properties and the 3D distribution of the atomic constituents. Among the optical properties we can mention the faculty of absorbing or emitting light upon excitation and how these responses depend on wavelength and polarization. These properties are crucial in the case of nanoscale systems constituting the building blocks of light-emitting or photovoltaic devices, but also of single photon sources for quantum cryptography. Objects such as quantum wells, quantum dots, nanoparticles embedded in insulating matrices, but also extended or point-like defects can be the object of such studies. Optical properties not only depend on the 3D distribution of the atomic species as it can be retrieved from APT datasets. Other structural factors, such as crystal symmetry or the degree of structural order, may also play a role, and a complete picture can be found in these cases through the complementary correlation of other techniques such as TEM.

However, optical properties may also be influenced by environmental factors, such as temperature, presence of electric or strain fields, and by the specific morphology of the nanoscale system under study. This also opens the way to the study of the physics of APT in an original way, as the instrumental development of an atom probe capable of recording the PL signal during APT analysis demonstrates.

We will thus show how the relationship between optics and APT may on one hand disclose original opportunities in materials science, but also how it may reach far beyond this, towards a sort of nanoscale laboratory capable of generating new horizons of understanding and beauty.

Correlative methods with optical spectroscopy

A possible classification of correlative methods based on optical spectroscopy and APT is the following:

Statistical correlation. This is the standard method applied to the interpretation of optical and structural properties of nanoscale light emitters via APT. The optical properties are studied on a macroscopic portion of the sample, (e.g., a thin film or a large bunch of nanoscale objects). Micro-PL may allow for selecting a microscopic portion of the sample and extract its signal. The APT analysis is conducted on a nm-scale -object extracted from the sample, without a strict relationship with the region probed by optical spectroscopy. A significant number of studies based on statistical correlation can be cited, for instance the studies on (In,Ga)N/GaN quantum wells performed by the teams led by R. Oliver and G.W. D. Smith/A. Cerezo/M. Moody. These works have assessed for instance the randomness of the In distribution within quantum wells grown on the polar c-plane of GaN and characterized the behavior of interfaces [START_REF] Galtrey | Three-dimensional atom probe studies of an InxGa1-xN∕GaN multiple quantum well structure: Assessment of possible indium clustering[END_REF][START_REF] Galtrey | Compositional inhomogeneity of a high-efficiency InxGa1-xN based multiple quantum well ultraviolet emitter studied by three dimensional atom probe[END_REF]. These are important elements for the interpretation of the carrier localization mechanisms -and thus for the explanation of the radiative efficiency -within these quantum confined emitters that constitute the active region of standard LEDs [START_REF] Weisbuch | Disorder effects in nitride semiconductors: impact on fundamental and device properties[END_REF]. Other statistical studies have pointed out that the (In,Ga)N behavior is different in quantum wells grown on non-polar planes, on which In clustering becomes possible [START_REF] Tang | Indium clustering in a-plane InGaN quantum wells as evidenced by atom probe tomography[END_REF][START_REF] Tang | Insight into the impact of atomic-and nano-scale indium distributions on the optical properties of InGaN/GaN quantum well structures grown on m-plane freestanding GaN substrates[END_REF]. However, there are many other examples of statistical correlation in the literature, not only on wide bandgap semiconductors [START_REF] Jeon | Alloy Fluctuations Act as Quantum Dot-like Emitters in GaAs-AlGaAs Core-Shell Nanowires[END_REF], but also on semiconductor nanoparticles within insulating matrices [START_REF] Roussel | Atomic characterization of Si nanoclusters embedded in SiO 2 by atom probe tomography[END_REF][START_REF] Beainy | Structural and optical study of Ce segregation in Ce-doped SiO1.5 thin films[END_REF][START_REF] Beainy | Direct Insight into Ce-Silicates/Si-Nanoclusters Snowman-Like Janus Nanoparticles Formation in Ce-Doped SiOx Thin Layers[END_REF] and on III-V nanowires [START_REF] Mancini | Three-dimensional nanoscale study of Al segregation and quantum dot formation in GaAs/AlGaAs core-shell nanowires[END_REF].

Sequential correlation on a single nanoscale object.

In this advanced approach, illustrated within the scheme of Fig. 1-(a), the optical and the APT analyses are carried out on the same individual nanoscale object. Intermediate TEM-micro-Photoluminescence (μPL) correlation is also possible and may add significant information which is not accessible by APT only (for instance, crystal symmetries, presence of extended defects such as stacking faults or dislocations, etc.) or improve the information obtained by APT, for instance by complementary electron tomography (ET). However, we can mention only rare cases in which an individual sample survived three different types of analyses [START_REF] Rigutti | Correlation of Microphotoluminescence Spectroscopy, Scanning Transmission Electron Microscopy, and Atom Probe Tomography on a Single Nano-object Containing an InGaN/GaN Multiquantum Well System[END_REF][START_REF] Mancini | Carrier Localization in GaN/AlN Quantum Dots As Revealed by Three-Dimensional Multimicroscopy[END_REF]. It is crucial to notice that sequential methods are based on the quite restrictive hypothesis that an APT sample is optically active. This requirement is far from obvious, as typical preparation methods require the use of focused ion beam (FIB), which introduces significant damage to the sample in terms of surface amorphization and swarms of point defects within up to several tens of nm from the surface. It is thus important that light emitters within APT samples are characterized by a strong localization of the carriers involved in the radiative transition: only in this way their wavefunctions avoid being influenced by the damaged regions close to the surfaces, which would otherwise yield an efficient source of non-radiative recombination paths. Up to now, III-N and II-O crystals have shown to be suitable for this approach.

In situ correlation via the Photonic Atom Probe. The last type of correlation is achieved through a dedicated instrument that combines structural analysis and optical spectroscopy capabilities, similarly to what is implemented in the case of Cathodoluminescence (CL) spectroscopy in an electron microscope (SEM or TEM). In the case of laser-assisted APT, it is a simple consideration that a laser pulse exciting a suitable sample at a proper wavelength and intensity can trigger ion evaporation in a Tomographic Atomic Probe and, at the same time, excite electron-hole pairs that can recombine radiatively and generate a PL signal. An instrument exploiting both ionic and photonic signals has been recently developed at GPM in Rouen, and is referred to as the Photonic Atom Probe [START_REF] Houard | A photonic atom probe coupling 3D atomic scale analysis with in situ photoluminescence spectroscopy[END_REF], whose principle is schematically illustrated in Fig. 1(b). This instrument can not only yield new performances in terms of correlative optical spectroscopy and atomic scale microscopy, but also represent a new tool for the study of surface and bulk light-matter interaction phenomena at the nanoscale, with or without an intense applied electric field.

From a Micro-Wire to a Nano-Tip

In order to correlate the 3D information provided by APT with the local PL signal, a bulk sample can be progressively milled by FIB until an optical active field emission tip is formed [START_REF] Rigutti | Correlation of Microphotoluminescence Spectroscopy, Scanning Transmission Electron Microscopy, and Atom Probe Tomography on a Single Nano-object Containing an InGaN/GaN Multiquantum Well System[END_REF][START_REF] Mancini | Carrier Localization in GaN/AlN Quantum Dots As Revealed by Three-Dimensional Multimicroscopy[END_REF][START_REF] Houard | A photonic atom probe coupling 3D atomic scale analysis with in situ photoluminescence spectroscopy[END_REF][START_REF] Mancini | Multi-microscopy study of the influence of stacking faults and three-dimensional In distribution on the optical properties of m-plane InGaN quantum wells grown on microwire sidewalls[END_REF]. This process leads to substantial variations in the luminescence spectrum, reducing the emission intensity and gradually revealing more and more details of the PL fine structure. This approach was successfully attempted to study the emission of a (In,Ga)N multi-quantum wells (QWs) system grown on the lateral and upper surface of GaN micro-wires (fig. 2-(a)). The PL spectra of a micro-wire is represented in the top part of fig. 2-(b). Even if QWs typically present well defined transition energies that are strictly linked to their morphology, composition and homogeneity, the observed spectrum presents broad and continuous peaks in a wide spectral region between 2.9 -3.3 eV. This suggests the presence of composition and/or morphology inhomogeneities in the QWs. Then, this spectrum was compared with those associated with a 500 nm diameter nano-cylinder and with a ~100 nm diameter needle-shaped APT tip containing both only upper part of the multi-QWs system, that are represented in fig. 2-(b).

Interestingly, discrete emission lines appear in the nano-cylinder PL spectrum. Similar lines are even better resolved in the case of the nano-tip, revealing the presence of several 5 meV large peaks that can be interpreted in terms of local carrier recombination within single QWs. Such a huge difference between the emission spectrum of the whole micro-wire and that associated to the resulting nano-tip is due to the very little volume that is probed with the µPL laser in the latter case. The nano-tip was then analyzed by APT. The 3D reconstruction depicted in fig. 2-(c) shows the position of In atoms, revealing the QWs, but also the formation of vertical In-rich planes associated to the presence of staking faults. Finally, the correlation between the structural and optical properties of the QWs was established though an effective mass calculation. The electron-holes wavefunctions (fig. 2-(f)), as well as the transition energies, were calculated on the bandgap landscapes (fig. 2-(e)) obtained on the basis of 3D chemical maps by including or neglecting the presence of In-rich stacking faults. Results as summarized in fig. 2-(f). The carrier localization in correspondence of the In-rich zones and of the stacking faults is also associated with a decrease in transition energies estimated to be about 15 meV. This result is in good agreement with the peak energy distribution observed in fig. 2-(b). Although it is not possible to strictly correlate each single QW with the relative emission line in the PL spectrum, this study confirms the possibility to statistically correlate the optical emission of a nano-tip with its APT reconstruction. [11] American Chemical Society and [START_REF] Mancini | Multi-microscopy study of the influence of stacking faults and three-dimensional In distribution on the optical properties of m-plane InGaN quantum wells grown on microwire sidewalls[END_REF] American Institute of Physics).

The Multi-Microscopy Approach

A more advanced statistical study was attempted performing a combination of sequential µPL, ET, and APT on a set of several tens of samples [START_REF] Mancini | Carrier Localization in GaN/AlN Quantum Dots As Revealed by Three-Dimensional Multimicroscopy[END_REF]. The system designated for this study consists in a stack of GaN layers containing truncated hexagonal pyramidal quantum dots (QDs) arranged in columns.

These structures are embedded in a matrix of AlN, acting as barrier leading to the confinement of both electron and hole wavefunctions within the QDs. Because the confinement energy depends on both dimension and shape of the QDs, the 3D information provided by techniques such as APT or ET is crucial to predict the transition energy spectrum emitted by a specimen. In particular, ET is the ideal technique to determine the dimensions of 3D nano-objects embedded in a tip with sub-nanometer precision and it could also provide information to improve the accuracy of APT reconstructions. However, morphology variations at a scale of less than 1 nm are typically associated to reconstruction artifacts.

The 3D reconstruction of Ga atom positions provided by APT is illustrated in fig. 3-(a). Here, QDs appear as high-density Ga-rich regions with blurred outlines along the plane perpendicular to the tip evaporation direction, due to limitations of APT lateral resolution associated to low-rage diffusion phenomena and local aberration effects. These problems are accentuated in correspondence of the GaN/AlN interfaces, due to the difference of the materials evaporation fields, as revealed by the composition map reported in fig. 3-(b). It should be noted that in this sample Ga/Al interdiffusion phenomena cannot occur during the growth process due to thermodynamic reasons. The vertical resolution of APT, in contrast, corresponds to an atomic monolayer (ML), that it is approximately equal to about 0.25 nm. Three different approaches were then pursued in order to statistically correlate the QDs morphology with their recombination energy.

The 3D composition maps of the QDs extracted from the APT reconstruction reveal interface monolayer fluctuations at the bottom of the QDs, as indicated by the red arrows in fig. 3-(c, top). The formation of atomic steps on the barrier surface could be ascribed to the formation of extended defects (e.g., dislocations) that propagate through the heterostructure. However, it is also possible to describe the QDs through three parameters only as truncated pyramids, exploiting the data obtained by ET (fig. 3-(c), middle). This simple description may then be enriched by adding the bottom interface fluctuations suggested by APT. These fluctuations consist in either one or two monolayers, with 1.5 nm diameter cylindrical shape (fig. 3-(c), bottom).

Comparing the emission spectra calculated based on APT data only (fig. 3-(e)) with the experimental PL peak distribution (fig. 3-(d)), the expected transitions energies are overestimated. In fact, lateral resolution limitations of APT artificially leads to the mixture of Al and Ga in correspondence of the QDs lateral interfaces, increasing the material bandgap within QDs. Contrariwise, the transition energy distribution expected for pyramidal QDs is quite close to the emission energies observed in PL measurements, as displayed in fig. 3-(f). In addition, the introduction of interface fluctuations has the effects of localize hole states (fig. 3-(c, bottom)), lowering the average transition energy and providing a narrower energy distribution. Thus, the best statistical agreement with this experimental distribution is obtained considering pyramidal QDs with the geometrical parameters extracted from ET and bottom interfaces affected by monolayer fluctuations, as suggested by APT.

Photonic Atom Probe and Super-resolution Photoluminescence

The possibility of measuring PL spectra simultaneously with the APT analysis, i.e., during specimen evaporation, is the specific characteristic of the Photonic Atom Probe. This opens up new methods and insights in the correlative analysis of optically active nanoscale systems. The PAP can be seen as the APT counterpart of cathodoluminescence in a scanning (transmission) electron microscope: it upgrades the instrument with optical spectroscopic capabilities, and it offers possibilities that go beyond the separated application of optical spectroscopy and structural/chemical/morphological study: as in many other systems, the whole is more than the sum of the parts. PAP has obviously very specific features. An example of this is discussed in the following. The PAP analysis of a set of 5 ZnO QWs embedded inside a (Mg,Zn)O alloy grown on a ZnO substrate is synthetically displayed in Fig. 4. These QWs have features that were particularly suitable for testing the new instrument: first, they have different thicknesses, which should result in different PL emission energies; second, they are quite closely packed with each other, with a mutual distance of around 20 nm. The set of PL acquired during the APT analysis is reported in Fig. 4-(a). The series begins with the spectral index 1 and finishes with 53. Three main spectral regions may be identified as (i) the emission of the ZnO substrate, around 3.35 eV, (ii) the emission of the QWs, between 3.35 and 3.55 eV, and (iii) the emission from the (Mg,Zn)O alloy, which constitutes a band peaked at 3.67 eV. The basic concept for the interpretation of these data is that the signal originating from a given emitter disappears when the emitter is no longer there because it is evaporated. Disappearing spectral features are highlighted by the red arrows in the graph. This quite obvious state of things allows attributing the different spectral signatures to the different sample features that may be retrieved through the APT reconstruction. However, it is also possible to plot the PL data as differential spectra, as shown in Fig. 4-(b). A differential spectrum being simply the difference between two consecutive spectra, this series of data emphasizes when main variations in the spectral intensity or energy occur. In this way, disappearing features in the PL series of spectra become peaks in the series of differential spectra. Such peaks (corresponding to the red arrows in Fig. 4-(a)) are marked within blue circles in Fig. 4-(b). The correlation of PL with APT may be straightforwardly carried out by analyzing the 3D reconstructions, in the form displayed in either Fig. 4-(c) or (d). As a guide for the reader, the position of the evaporation front during selected PL spectra acquisitions is underlined by red curves. The selected indices i = 7, 11, 16, 24 correspond to the spectra in which the PL signals marked by the red arrows in Fig. 4-(a) disappear. It can be easily recognized that this sets an immediate correspondence between the disappeared spectrum and the QW which is about to be evaporated (as a matter of fact, QW PL signals have the tendency to disappear before the QW is completely evaporated, but this is expected because the structure of the QW is strongly affected even before its complete evaporation). Furthermore, differential spectra can be exploited in order to encode the intensities of different spectral intervals into the positions of the atoms, as it is displayed in fig. 4-(c). Here it becomes immediately visible that the 3.42 eV signal (coded in green) is generated by QW5, while the 3.49 eV signal (in red) from QW4.

In any case, the technique allows singling out the signal of an individual quantum confined system and to discriminate it among the neighboring emitters, which are as close as 20 nm. This spatial resolution (which is relative to the APT axis of analysis only) is over one order of magnitude smaller than the laser wavelength, which approximately corresponds to the standard spatial resolution of µPL as set by the diffraction limit. The PAP allows thus for super-resolution PL, which is a nice and (a posteriori) straightforward consequence of collecting a light signal from a slowly and gently evaporating tip. The cross-sectional composition map of Fig. 4-(d) exhibits quite clearly the alloy distribution, which can be studied by statistical means, and the presence of interface roughness. It is well known that these features may be quite effective in the localization of carriers within the QW (especially of holes, localizing over the scale of only several nm). Such features are considered when calculations are performed in order to understand whether the behavior of the light emitter is expected within the most common models. For instance, we adopted the quite simplified but robust approach of conducting effective mass 6 bands calculations, and we found a good agreement between the experimental PL energies and those calculated on the basis of the morphological and compositional characteristics of the system as obtained by APT. These results demonstrate that the PAP can be profitably applied to materials science, in particular to the correlative study of optically active materials. Other studies have already been carried out in (In,Ga)N/GaN QDs systems, and further studies are planned. Instrumental upgrades have already resulted in the possibility of longer and more reliable analyses. In the meantime, it is necessary to explore its limits in terms of application to materials and systems (quantum emitters, extended or point defects, etc.).

A nanoscale laboratory for the study of light-matter-field interaction

Yet, the PAP being a microscope equipped with optical spectroscopy capabilities, its vocation reaches beyond the mere application to materials science studies. The PAP (as the APT itself, but differently) is also a nanoscale laboratory, in the sense that it makes it possible to play with light, field and matter in order to induce new or just so far unobserved phenomena, and to understand something more about their physics. One example of this is represented by the opto-mechanical studies of field-emission tips as a function of the applied field or of the evolution of their shape. Again, the PAP represents a quite original possibility of studying the optical properties of a nanoscale object which changes its shape over time. As it is known, the application of a high voltage to a tip results in a stress building up at the apex. Many APT are aware of that, especially when their sample fractures because it did not withstand mechanical stress. The relationship between field F and stress 𝜎 is quite simple, 𝜎 = 1 2 ⁄ 𝜀 0 𝐹 2 , which is by the way the same as in any metallic surface, including that of a planar capacitor.

Of course, we can wonder whether this formula also holds in the case of semiconducting tips. We know from past works that the electric field is supposed to be expulsed from field-emission tips under high voltage [START_REF] Silaeva | Do Dielectric Nanostructures Turn Metallic in High-Electric dc Fields?[END_REF], but the measurement of stress in a nanoscale object is never an obvious experimental problem. This problem was considered in the PAP experiment reported in Fig. 5 [START_REF] Dalapati | Situ Spectroscopic Study of the Optomechanical Properties of Evaporating Field Ion Emitters[END_REF]. In this case, we

5-(d).

As the apex stress was not a priori known, we let it span an interval within which the best fitting value of 𝜎 =1.25 GPa could be found. This value corresponds to an apex average surface field F = 16

V/nm, which is consistent with the microscopic field evaluated through the charge state ratio statistics, Fµ = 20 V/nm. The whole picture is also consistent with the stress-field relationship which is indirectly reported in the left panel of the graph: the PL shift exhibits indeed a quadratic dependence on the applied voltage. These results also confirm our past piezospectroscopic studies of color centers in diamond nanoneedles, which were performed in a proto-PAP not yet equipped with APT detectors [START_REF] Rigutti | Optical Contactless Measurement of Electric Field-Induced Tensile Stress in Diamond Nanoscale Needles[END_REF]. They also provide an independent way of assessing surface electric fields, which play an important role in the field evaporation of compound materials [START_REF] Mancini | Composition of Wide Bandgap Semiconductor Materials and Nanostructures Measured by Atom Probe Tomography and Its Dependence on the Surface Electric Field[END_REF][START_REF] Russo | Composition Metrology of Ternary Semiconductor Alloys Analyzed by Atom Probe Tomography[END_REF].

Further phenomena are currently being studied. The dataset of Fig. 5-(b) indicates that not only there is an energy shift, but also an intensity variation which is far from obvious: why does the QW PL intensity increase during the first phases of the evaporation? The mechanism leading to this intensity variation is possibly to several factors: (i) the modification of the light emission pattern through the morphology variation, which could couple more or less radiation into the spectrometer (ii) the modification of charge generation, also related to the evolution of the absorption of light during evaporation and (iii) the evolution of the transport of charge carriers within the system -also related to the morphology change, to the stress or to the progressive erosion of regions with a high density of defects at the surfaces.

Perspectives

Currently, correlative APT and optical spectroscopy displays an intriguing set of possibilities that could be applied in materials science or in nanoscale physics. Even if in the opinion of the authors there is very little potential for its widespread application because of the complexity of the issued datasets, it may represent an original way to tackle specific materials science problems. For instance, the study of coupled QDs, the study of parametric dependences of spectral signatures on applied fields, all of this correlated with the sub-nanoscale resolution 3D distribution of atomic species. On the other hand, its "nanoscale lab" vocation discloses the possibility of new experiences targeting the transport of charge carriers in APT specimens, the implementation of optically measured strain-stress curves, maybe also a deeper understanding of fracture events (they produce an electrical arc which could be recorded by the PL spectrometer) and of the chemical processes occurring at evaporating surfaces, with significant feedback for field emission physics, APT metrology and correlative microscopy.

Figure 1

 1 Figure 1 Advanced correlative methods based on photoluminescence spectroscopy and atom probe: (a) The sequential correlation is based on the analysis of a single individual nanoscale sample by µPL first, then by APT. An intermediate TEM step is also possible, in order to access the structural properties otherwise not accessible by APT. (b) The in-situ correlation has recently become possible thanks to the instrumental development of the Photonic Atom Probe, in which PL spectroscopy may be carried out during APT analysis.

Figure 2 (

 2 Figure 2 (a) Micro-wire (insert) and field-emission tip containing a set of (In,Ga)N/GaN quantum wells; (b) µPL spectra issued from microwires (top), from a 500 nm diameter cylinder fabricated by FIB and from an APT tip (bottom); (c) 3D reconstruction of the position of In atoms within the multi-QW system; (d) Composition maps from selected QWs; (e) Bandgap landscape (left) and calculated electron/hole wavefunctions based on APT reconstructions with (top) and without (bottom) the effect of a stacking fault (Adapted with permission from ref.

Figure 3

 3 Figure 3 Multi-microscopic analysis of a set of GaN/AlN quantum dots: (a) The position of the Ga atoms in a portion of APT reconstructed volume; (b) correspondent composition map (thin slice); (c) composition maps adopted for the calculation of the electron/hole wavefunction and energies. Top: extracted from the APT maps. Middle: simplified pyramid geometry based on ET data. Bottom: simplified geometry + ML fluctuation; (d) Statistical distribution of PL energies recorded from 15 tips; (e) calculated PL energies based on APT data; (f) calculated PL energies based on simplified geometries including or not ML fluctuations; (g) average calculated energies (Adapted with permission from ref. [12] American Chemical Society).

Figure 4

 4 Figure 4 In-situ Photonic Atom Probe analysis of a system of 5 ZnO/(Mg,Zn)O quantum wells: (a) Complete set of collected µPL spectra, the red arrows indicating the approximate location of the disappearance of the PL peaks of consecutive QWs; (b) Differential spectra recorded in the sub-region indicated by the blue rectangle in (a). The optical signatures of over 5 quantum wells are highlighted by the blue circles; (d) 3D reconstruction of the positions of Mg atoms; the color of the atoms is color-coded according to the PL intensities in the spectral intervals indicated at the top and in the graph in (b); (e) 2D-distribution of the Mg II-site fraction calculated over a 5 nm thick slice close to the tip axis. The red dashed lines indicate the approximate position of the specimen surface during the acquisition of the i-th spectrum (Adapted with permission from ref. [15] American Chemical Society).

  analyzed a ZnO/(Mg,Zn)O double QW system, of which only one QW is optically active (in the APT reconstruction of Fig.5-(a), the optically active QW is the thinner one on the left, closer to the ZnO substrate). The series of PL spectra is shown in Fig.5-(b): the QW signal is clearly appearing emitting at an energy of around 3.45 eV, while the other two signals can be easily identified (see also Fig.(4)) as the emissions of the ZnO substrate and of the (Mg,Zn)O barrier. The PL energy of the QW gradually shifts over time. Its shift evolves with the progression of the evaporation, as show in the right panel of the graph in Fig.5-(d). The left panel of the same graph also indicates a dependence on the application of the voltage before the evaporation. Quite clearly, this energy shift is related to the stress induced by the electric field. It is far from straightforward to evaluate the stress at the QW position, but it is fortunately possible. First, one has to understand how the apex stress is transmitted to the rest of the structure down the tip shaft -which can be performed once one knows the morphology of the tip and the differences between the elastic constants of the materials composing the structure; then, from this stress the strain tensor can be retrieved (through the elastic constants). Finally, the strain tensor allows for the calculation of the bandgap shift via the set of deformation potentials. Two examples of calculations of locally resolved bandgap shifts are reproduced in Fig.5-(c) for two different sample morphologies, i.e., at the evaporation beginning and close to the evaporation of the emitting QW. The bandgap shift represents per se a good estimation of the PL shift, which is reported in the graph of Fig.

Figure 5

 5 Figure 5 Study of stress-induced PL shift in a ZnO/(Mg,Zn)O QW within a Photonic Atom Probe. (a) The reconstructed position of the Mg atoms. The emitting QW is the one at roughly 50 nm from the ZnO substrate / (Mg,Zn)O interface (set at 0); (b) Complete set of PL spectra acquired during the APT acquisition; (c) calculation of the distribution of the stress-induced bandgap shift within the tip at two different evaporation stage (beginning, top, and evaporation of QW1, bottom); (c) Comparison between the experimental and a set of calculated trends of the QW PL energy as a function of the applied voltage before evaporation (left panel) and as a function of the apex position during the evaporation. The best fitting curve is obtained for a stress at the apex equal to 1.25 GPa (Adapted with permission from ref. [17] American Physical Society).

Acknowledgements

The authors acknowledge funding from French National Research Agency (ANR) in the framework of the projects EMC3 Labex AQURATE, EMC3 Labex IDEPOP and ANR-13-JS10-0001-01 TAPOTER and cofunded in the framework of RIN IFROST, EMC3 Labex IDEPOP and CPER BRIDGE projects by European Union with European Regional Development Fund (ERDF) and by Region Normandie.