Nikolai Romashchenko

Benjamin Linard

Eric Rivals

Fabio Pardi

ON COMPUTATIONAL BIOLOGY AND BIOINFORMATICS 1 Computing Phylo-k-mers

Keywords: phylogenetic placement, k-mers, enumeration algorithms, empirical runtimes, phylogeny ✦

Finding the correct position of new sequences within an established phylogenetic tree is an increasingly relevant problem in evolutionary bioinformatics and metagenomics. Recently, alignment-free approaches for this task have been proposed. One such approach is based on the concept of phylogenetically-informative k-mers or phylo-k-mers for short. In practice, phylo-k-mers are inferred from a set of related reference sequences and are equipped with scores expressing the probability of their appearance in different locations within the input reference phylogeny. Computing phylo-k-mers, however, represents a computational bottleneck to their applicability in real-world problems such as the phylogenetic analysis of metabarcoding reads and the detection of novel recombinant viruses.

Here we consider the problem of phylo-k-mer computation: how can we efficiently find all k-mers whose probability lies above a given threshold for a given tree node? We describe and analyze algorithms for this problem, relying on branch-and-bound and divide-and-conquer techniques. We exploit the redundancy of adjacent windows of the alignment to save on computation. Besides computational complexity analyses, we provide an empirical evaluation of the relative performance of their implementations on simulated and real-world data. The divide-and-conquer algorithms are found to surpass the branch-and-bound approach, especially when many phylo-k-mers are found.

INTRODUCTION

A LIGNMENT-FREE approaches in bioinformatics are mo- tivated by the fact that sequence alignment is a complex task requiring the use of memory and time-consuming algorithms. Moreover, alignments are potentially inaccurate, sensitive to sequencing errors, and difficult to apply to genomes with permuted structures [START_REF] Zielezinski | Alignment-free sequence comparison: benefits, applications, and tools[END_REF]. Many alignmentfree methods for solving various problems in bioinformatics (e.g., de novo assembly, genome comparison, read correction, read clustering) rely on the decomposition of a sequence into its constituent k-mers, that is, its overlapping substrings of length k.

Recently, we and other authors proposed a probabilistic extension of the notion of k-mers [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], [START_REF] Scholz | Rapid screening and detection of inter-type viral recombinants using phylo-k-mers[END_REF]. In this development, many more k-mers are inferred from a set of reference sequences beyond the ones that are actually within those sequences. This inference aims at predicting k-mers that may be present in relatives of the reference sequences (e.g., within their ancestors or within "cousin" sequences). Moreover, for any given location in the phylogeny of the reference sequences, one can estimate the probability of observing any given k-mer, meaning that probability scores can be assigned to the inferred k-mers. Key to this inference are probabilistic models of sequence evolution, which rely on a phylogenetic tree describing the evolutionary history of the reference sequences. The inferred k-mers are intended to be informative about the phylogenetic origin of newlyobserved sequences containing them. For these reasons, they are called phylo-k-mers.

• N. Romashchenko, B. Linard, E. Rivals, F. Pardi were with MAB, LIRMM, CNRS, University of Montpellier, 34095 Montpellier, France. B. Linard was also with UR MIAT, F-31320, INRAE, University of Toulouse, 31320 Castanet-Tolosan, France. E-mails: nromashchenko@lirmm.fr, benjamin.linard@inrae.fr, rivals@lirmm.fr, pardi@lirmm.fr Every phylo-k-mer w is associated with scores describing how probable w is to appear at a predefined set of nodes in the reference phylogeny (more detail in the Preliminaries). These scores can be used to determine the likely phylogenetic origin of any given query sequence while avoiding the need to align the query to the reference sequences. This task is also known as phylogenetic placement, which is an increasingly relevant problem in phylogenetic analysis and metagenomics [START_REF] Czech | Metagenomic analysis using phylogenetic placement -a review of the first decade[END_REF]. While phylo-k-mers were initially developed for the phylogenetic placement of metabarcoding reads [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], they were also applied to the detection and analysis of virus recombinants composed of fragments from different viral types of the same species [START_REF] Scholz | Rapid screening and detection of inter-type viral recombinants using phylo-k-mers[END_REF].

The main bottleneck of this technique lies in the very large number of phylo-k-mers, which comes from the fact that we need to consider up to 4 k k-mers for DNA and 20 k for protein sequences. Although we can reduce this number by only considering phylo-k-mers with probability scores above a certain threshold, the thresholds that have been used in practical applications are typically low. Thus, finding phylo-k-mers remains computationally challenging. While previous works only considered the accuracy and speed of sequence classification based on already computed phylo-k-mers [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], [START_REF] Scholz | Rapid screening and detection of inter-type viral recombinants using phylo-k-mers[END_REF], here we focus on algorithms for computing phylo-k-mers.

In the following, we consider a number of algorithms for this problem. While one of these algorithms has already been described to some degree in the literature (e.g., [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF], [START_REF] Pizzi | Finding significant matches of Position Weight Matrices in linear time[END_REF], [START_REF] Salmela | Algorithms for weighted matching[END_REF]), the others are novel. We analyze the complexities of all the presented algorithms and compare their running times over simulated and real-world datasets. Both theoretical analyses and empirical evaluations show that the new algorithms can be significant improvements over the existing one, especially when a large number of phylo-kmers must be output.

Related works

A problem similar to phylo-k-mer computation arises in the context of sequence motifs, precisely of Position-Specific Scoring Matrices (PSSMs), also known as Position Weight Matrices (PWMs) or weighted patterns. PSSMs represent DNA and protein sequence motifs (e.g., transcription factor binding sites) as a matrix of probabilities for each nucleotide, or amino acid, at each position in the motif. An important problem is to find significant matches of such weighted patterns in collections of genome-sized sequences. In existing algorithmic solutions to this problem, one of the preliminary steps is to enumerate all possible motif instances that reach a threshold score for a given PSSM. This step is similar to the problem of phylo-k-mer computation, with some important differences that we discuss below. Previous literature showed that the tree of all prefixes of motifs with high-enough scores could be efficiently explored in a depth-first [START_REF] Salmela | Algorithms for weighted matching[END_REF], [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF] and breadth-first manner [START_REF] Pizzi | Fast search algorithms for position specific scoring matrices[END_REF], [START_REF] Pizzi | Finding significant matches of Position Weight Matrices in linear time[END_REF].

However, in the context of phylo-k-mers, the computation is more challenging: the PSSM-based approaches only involve a single execution per profile, and the number of profiles to process is usually in the hundreds [START_REF] Kulakovskiy | HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP-Seq analysis[END_REF], [START_REF] Fornes | JASPAR 2020: update of the open-access database of transcription factor binding profiles[END_REF]; on the other hand, computing phylo-k-mers may well require processing millions of matrices, as it must process each of the k-wide sub-matrices of several input matrices originating from different parts of the reference phylogeny. Another difference is that, for phylo-k-mers, score threshold values are typically much lower than for PSSM matching, meaning that a larger fraction of the possible k-mers can reach the threshold. Finally, phylo-k-mer computation assumes processing matrices related to each other, both because kwide sub-matrices overlap and because of the phylogenetic relatedness of the input matrices. We exploit the overlap between sub-matrices to improve the running time of phylok-mer computation.

PRELIMINARIES

Notation

Let Σ be a finite ordered alphabet of cardinality σ. We consider strings (or sequences) over alphabet Σ. Let k be a positive integer. Let Σ k denote the set of all possible strings of length k over Σ. Given a string s, the length of s is denoted by |s|. For any two integers 1 ≤ i ≤ j ≤ |s|, s i denotes the i th character of s, and the substring of s starting in position i and ending at position j is denoted by s i . . . s j . A substring s i . . . s j is a prefix of s if i = 1, and a suffix of s if j = |s|. For a set X, |X| denotes the number of elements in X.

We consider matrices whose rows are indexed by characters of Σ and whose columns are indexed as the positions of a multiple alignment. A column stores the probability of occurrences of each possible character at that position. Hence, we term such matrices probability matrices since the values of a column sum to one. For a σ × m probability matrix P , P α,j denotes the element on row α (with α ∈ Σ) and column j of P (with 1 ≤ j ≤ m); the same element is denoted by P ij if α is the i-th element of Σ. For two integers i, j such that 1 ≤ i ≤ j ≤ m, P [i : j] denotes the matrix P restricted to columns from i to j included.

v 1 u 1 u 3 a b v 3 v 4 c v 2 u 2 u 4
Fig. 1: A toy reference tree (solid lines) with three leaves a, b, c (filled squares), which correspond to the (observed) reference sequences, for which a multiple alignment is given as input. To this reference tree, we add the nodes in

V ′ = {u 1 , u 2 , u 3 , u 4 , v 1 , v 2 , v 3 , v 4
} (filled circles), representing unobserved relatives of a, b, c. Some of these nodes represent ancestral sequences (u 1 , u 2 , u 3 , u 4 ), while some others represent "cousin" sequences (v 1 , v 2 , v 3 , v 4 ) related to the reference tree via newly added edges (dashed lines).

For each of these nodes, we can obtain probability matrices {P ui }, {P vi }, on the basis of the input alignment and of the reference tree. These matrices are the input of the phylo-kmer computation problem.

Phylo-k-mers at a glance

Consider a multiple alignment of reference sequences and a phylogenetic tree T = (V, E) describing the evolutionary history leading up to the reference sequences. We add to T a set of nodes V ′ , representing sequences that are unknown relatives of the reference sequences. (See Figure 1 for an example.) Let m be the number of columns (sites) in the alignment. For each node u ∈ V ′ , we compute a σ × m probability matrix P u describing the probability at u of any character in Σ, at any site in the alignment, conditional to the sequences observed at the leaves of T (i.e., the aligned reference sequences). P u can be derived from the tree likelihood conditional to the character at u by applying Bayes' theorem, which is standard in phylogenetics (see, e.g., section 4.4.2.1 in [START_REF] Yang | Computational molecular evolution[END_REF]). Then, the complexity of computing all matrices P u is equal to that of computing conditional tree likelihoods across all tree nodes, which for a constant-size alphabet can be done in O(|V ∪ V ′ | • m) time [START_REF] Bryant | Likelihood calculation in molecular phylogenetics[END_REF] with Felsenstein's algorithm [START_REF] Felsenstein | Evolutionary trees from DNA sequences: a maximum likelihood approach[END_REF].

Given P u , we can then define a probability score S u (w) associated to any given k-mer w and to the node u. See Definition 1 below for a definition of S u (w) (where the superscript is dropped for simplicity). Informally, S u (w) approximates the probability of w to appear in a sequence positioned at node u, based on the chosen model of sequence evolution and on the sequences at the leaves of T . We call the pair (w, S u (w)) a phylo-k-mer.

The interest of phylo-k-mers is that finding the nodes u that maximize the product of S u (w) over all k-mers in a query sequence provides a good estimate of the evolutionary origin of the query [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], [START_REF] Scholz | Rapid screening and detection of inter-type viral recombinants using phylo-k-mers[END_REF]. Moreover, this can be computed without aligning the query to the reference sequences, making this approach very scalable to large numbers of queries. For a detailed treatment of phylo-k-mers, see [START_REF] Romashchenko | Computing informative k-mers for phylogenetic placement[END_REF]. While the matrix P u and score function S u are relative to a particular node u, in the following, we assume that node u is fixed and therefore omit this dependency. We simply write P and S.

The problem of phylo-k-mer computation

Here, we study the problem of enumerating k-mers and their scores relative to a probability matrix P and a threshold score value ε ∈ [0, 1). P contains probabilities P α,j of observing different characters α ∈ Σ at every site j of the multiple alignment. Starting from an alignment site j, or position j, we can calculate the score of a k-mer w = w 1 w 2 . . . w k for this position by taking the product of corresponding probabilities: S(w, j) = P w1,j • P w2,j+1 • . . . • P w k ,j+k-1 . We say that w obtains the score of S(w, j) at position j. We only consider k-mers that obtain scores greater than ε for at least one position. For such a k-mer w, we say that w reaches the threshold at position j if S(w, j) > ε. The final score S(w) is the maximum of S(w, j) obtained among all positions. Definition 1 formalizes this problem.

ALGORITHMS

Phylo-k-mer computation has been implemented in RAP-PAS [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF] but has not been described explicitly. Here, we describe an algorithm similar (but not equivalent) to the one of RAPPAS and present new algorithms for this problem. Unless otherwise stated, all described algorithms approach the problem window-by-window: given a window W = P [j : j + k -1] of k consecutive columns in P , we list all k-mers that reach the threshold for the window, as well as their scores. Let Z be the set of such k-mers for the window W . If w ∈ Z, we call w alive in the window, and we call it dead otherwise. Then, we can obtain the solution for P by simply taking the union of sets Z for every window and setting the score of each k-mer to the maximum score obtained across all windows.

In the analysis of the algorithms, we adopt the word-RAM model of computation. It assumes operating on words of size b and performing arithmetic and bitwise operations in constant time [START_REF] Hagerup | Sorting and searching on the word RAM[END_REF]. Also, we assume that the alphabet size σ is constant. Finally, we assume that any k-mer can be represented with a constant number of machine words, implying b = Θ(log σ k ). Those assumptions imply that we can operate on binary representations of k-mers (e.g., writing a k-mer to memory) in constant time.

Branch-and-bound

Given a window W , a naïve algorithm, implemented in RAPPAS, iterates over possible prefixes in a depth-first manner. For a prefix p = w 1 . . . w l with a score l j=1 W wj ,j > ε, it expands p by one character and checks whether the score of the expanded prefix also reaches the threshold. As soon as a prefix obtains a score ≤ ε, it is rejected. This algorithm can be naturally improved with the lookahead bound technique (introduced in [START_REF] Wu | Fast probabilistic analysis of sequence function using scoring matrices[END_REF], also used in [START_REF] Beckstette | Fast index based algorithms and software for matching position specific scoring matrices[END_REF], [START_REF] Korhonen | MOODS: fast search for position weight matrix matches in DNA sequences[END_REF], [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF]). Consider a lookahead bound array L of p ← 2 ⌈log 2 σ⌉ p + i -1 ▷ Update the binary representation of p and its score 3: BRANCHANDBOUND(i, 1, 0, 1); 20: end for 21: return Z elements L j = k h=j+1 max a∈Σ W a,h giving maximum possible scores achieved in W by suffixes of different lengths. Then, a prefix p = w 1 . . . w l of length l can be rejected if l j=1 W wj ,j ≤ ε/L l . By analogy with k-mers, we call p alive if its score reaches ε/L l , and dead otherwise. Note that a prefix is alive if and only if it is the prefix of an alive k-mer, i.e., an element of Z.

s ← s • W ij ; 4: if s ≤ ε/L j then ▷ Lookahead
Figure 2 gives the pseudocode of the recursive depthfirst branch-and-bound algorithm. Similar algorithms were described for preprocessing PSSMs in depth-first [START_REF] Salmela | Algorithms for weighted matching[END_REF], [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF] and breadth-first [START_REF] Pizzi | Fast search algorithms for position specific scoring matrices[END_REF], [START_REF] Pizzi | Finding significant matches of Position Weight Matrices in linear time[END_REF] manners. In some cases (e.g., [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF]), the columns of the PSSM were ordered by conservation to facilitate early rejection of prefixes. This idea can easily be adapted for phylo-k-mer computation by ordering the columns in each window by the entropy of the probability distribution that they define. However, in practice, we did not find this to be worth the computational overhead it involves (see section B in Appendix). Proof. Let us consider the call tree of the algorithm where every tree node of depth j corresponds to considering a prefix of length j. We call a node alive if it corresponds to an alive prefix, and dead otherwise. Let ξ j A and ξ j D be the numbers of visited nodes of depth j that are alive and dead, respectively. Trivially, ξ k A = |Z|. Note that every alive prefix of length j -1 is extended into at least one alive prefix of length j, implying that ξ j-1

A ≤ ξ j A . Therefore, ξ 1 A ≤ ξ 2 A ≤ • • • ≤ ξ k-1 A ≤ ξ k A , and k j=1 ξ j A ≤ kξ k A = k|Z|. Now, let us count dead nodes: ξ j D < σξ j-1
A , and since

ξ j-1 A ≤ ξ j A , then ξ j D < σξ j A . Therefore, k j=1 ξ j D < k j=1 σξ j A = σ k j=1 ξ j A ≤ σk|Z|.
Finally, the total number of visited nodes is k j=1 (ξ j A + ξ j D ) < k|Z| + σk|Z| = (σ + 1)k|Z| = O(k|Z|), assuming that σ is a constant. We visit every node in constant time by virtue of the word-RAM model. Besides that, it takes Θ(σk) to precompute L. Then, the total time complexity is O(σk + k|Z|) = O(k|Z|).

Theorem 1 shows the worst-case complexity of the branch-and-bound algorithm to be O(k • |Z|). However, the algorithm achieves optimal best-case complexity: consider ε = 0 and W consisting of strictly positive probabilities, for which |Z| = σ k . The algorithm visits k j=0 σ j = (σ k+1 -1)/(σ -1) = Θ(σ k ) = Θ(|Z|) nodes; including preprocessing time, it takes Θ(k + |Z|) = Θ(|Z|) time in the best case. Finally, we note that it is possible to construct examples for which |Z| = Θ(k c ) for a small constant c, and branch-and-bound runs in

Θ(k c+1 ) = Θ(k • |Z|) time,
showing that the upper bound in Theorem 1 is tight in these cases. See Appendix A for one such example.

Divide-and-conquer

We present a new algorithm for the problem of phylo-k-mer computation. It applies the divide-and-conquer technique to compute scores of prefixes and suffixes for a given window W of size k. It also relies on a score bounding technique similar to the one discussed above. Consider the array {max a∈Σ W a,j : j = 1 . . . k} giving maximum score values for every column. Then, let M be a data structure answering range product queries M (j 1 : j 2 ) in constant time (see Appendix E):

M (j 1 : j 2 ) = j2 l=j1 max a∈Σ W a,l
We start with constructing M for W , which can be done in time linear in the size of W . Trivially, if M (1 : k) ≤ ε, we return an empty list as no k-mer is alive in the window. Otherwise, we split W into two subwindows of sizes ⌊k/2⌋ and ⌈k/2⌉. We compute L, defined as the list of ⌊k/2⌋-mers that reach the score of ε l = ε/M (⌊k/2⌋ + 1 : k) in the left subwindow. Similarly, we compute R, the list of ⌈k/2⌉-mers that reach the score of ε r = ε/M (1 : ⌊k/2⌋) in the right subwindow. Note that every ⌊k/2⌋-mer in L must be a prefix of at least one alive k-mer, and every ⌈k/2⌉-mer in R is a suffix of an alive k-mer. The procedure described above is applied recursively to every subwindow until, at the bottom of the recursion, we process a column j and select 1-mers reaching the score of ε/ k l=1,l̸ =j max a∈Σ W a,l . We combine the results of the recursive calls as follows: if |L| < |R|, swap them; sort R (the smaller of the two lists) by score. Finally, for every l ∈ L, consider the elements r ∈ R in descending order of scores; include the sequence obtained by concatenating l and r in the output, while the concatenated sequences are alive. Figure 3 gives the pseudocode of this algorithm. Proof. First note that both σ ⌊h/2⌋ and σ ⌈h/2⌉ are Θ(σ h/2 ), so we drop all notation expressing integer rounding. Z ← empty list; swapped ← f alse

3: if h = 1 then 4: return {(i -1, W i,j ) : W i,j > ε ′ for i ← 1 . . . σ } 5:
end if 6:

ε l ← ε ′ / M (j + ⌊h/2⌋ : j + h -1) 7:
ε r ← ε ′ / M (j : j + ⌊h/2⌋ -1)

8:

L ← DC(j, ⌊h/2⌋, ε l )

9:
R ← DC(j + ⌊h/2⌋, ⌈h/2⌉, ε r ) We consider a generic recursion call acting on a subwindow of size h ≤ k and analyse its runtime complexity. The total runtime will be obtained by summing runtimes over all recursion calls. Lines 2-12, except the calls to DC(), take constant time. Line 13 involves sorting a list of (h/2)mers, whose size is trivially at most σ h/2 . This can be done in no more than c

•σ h/2 log σ h/2 = c ′ •h/2•σ h/2 = O(h•σ h/2
) time (for some positive constants c, c ′ , and assuming σ is constant). We then move on to the complexity of the loops at lines 14-24 to complete the analysis of a single recursion call.

Let φ h denote the number of alive h-mers for a recursive call acting on a window of size h, φ h ≤ σ h . Note that every (h/2)-mer in L gives rise to at least one alive h-mer, and leads to considering at most one dead h-mer. Thus, at most one dead h-mer is considered per alive h-mer. Then, the total number of h-mers considered (dead and alive) by the loops is Θ(φ h ). In other words, lines 16-22 are executed Θ(φ h ) times, each of which takes constant time under the assumptions of the word-RAM model. In conclusion, overall, a single recursion call takes O(φ h + h • σ h/2 ) time.

The top-level recursion call takes O(|Z| + k • σ k/2 ) time, which is the claimed asymptotic complexity. We now show that the total time spent in deeper recursion calls is dominated by this complexity. Trivially, the runtime for a single recursion call of depth d ≥ 1 is bounded above by

c 0 • (|Z| + k • σ k/2 ) c • σ k/2 c • σ k/4 . . . c • σ n . . . 2 1 0 . . . . . . . . . . . . c • σ k/2 . . . . . . . . . c • σ k/4 . . . . . . c • σ = k • c • σ . . . 4 • c • σ k/4 2 • c • σ k/2 O(k • σ k/2 ) + + + + • • • + + + = = = + + • • •
(k • σ k/2 ). c • σ h = c • σ k/2 d
, for some positive constant c (because

φ h ≤ σ h and h • σ h/2 ≤ c ′ • σ h , for c ′ ≥ 2).
σ i = σ k/2+1 -1 σ -1 = O(σ k/2 ).
Theorem 2 gives an upper bound for running time of the algorithm as a function of the output size. Intuitively, the algorithm achieves linear complexity in output size for |Z| sufficiently large. This can be illustrated by the same example as for branch-and-bound: if ε = 0 for W of positive values, then all σ h h-mers are alive for every recursive call. It is then easy to see that the top call runs in Θ(σ k ) time, while all other calls take Θ(kσ k/2 ) in aggregate, giving a total runtime of complexity Θ(|Z|) = Θ(σ k ).

Divide-and-conquer with Chained Windows

While the problem of computing phylo-k-mers (Definition 1) is defined for a σ × m matrix containing many σ × k windows, the algorithms described above only consider one window at a time. Thus, they ignore an important property of the sequence of windows of P : two adjacent windows share (k-1) identical columns, meaning that some computation is redundant. Based on this observation, we suggest an improvement to the divide-and-conquer algorithm, illustrated in Figure 5.

We explain the idea for specific input and later will show how to generalize it to any input. Let k be an even value, and let the matrix P be such that max a P a,j is constant, ∀j ∈ {1, . . . , m}. Then, local thresholds ε l and ε r for a fixed recursion level are equal and constant for all windows. Consider a window W at position j, for which we recursively process its right subwindow W [k/2 + 1 : k], obtaining the list R of alive (k/2)-mers and their scores > ε r . Then, the list R is identical to the list L + of alive (k/2)-mers for the left subwindow W + [1 : k/2] of another window W + starting at position (j + k/2): it corresponds to the same range of columns (see Figure 5a) and is computed for the same threshold. Naturally, we can reuse R to compute the phylok-mers of W + . This allows us to make only one top-level recursive call for W + instead of two (Figure 5a). We iterate over windows with a step of k/2, always keeping the list R of the preceding window for the next one. A sequence of windows at a distance of k/2 from each other is called a chain of windows. We need to process k/2 such chains starting at positions 1, 2, . . . , k/2 to cover all windows of P . Figure 5b illustrates this idea. Note that we still have to make both recursive calls for the first window of every chain.

The described example relies on the assumption that the threshold ε r computed for W is equal to the threshold ε l computed for W + , which from here onwards we call ε + l , to distinguish it from the threshold for the left subwindow of W . This allowed us to assume R = L + , where L + is the list of alive (k/2)-mers for the left subwindow of W + . Of course, ε r and ε + l are generally not equal, meaning that R ̸ = L + . However, it is easy to see that one of these lists is always contained in the other: if ε + l < ε r , then R is a subset of L + , and vice versa otherwise. To be sure not to lose any alive (k/2)-mer for the subwindow shared by W and W + , we then compute the list of (k/2)-mers that reach min(ε r , ε + l ). This list equals R ∪ L + , the largest of R and L + .

The problem now becomes how to retrieve R from

W W L W R W + L W + R W + W - R W - (a)
... R ∪ L + , when computing alive k-mers for W , and how to retrieve L + from R ∪ L + , when computing alive k-mers for W + . This can be achieved as follows: rearrange R ∪ L + to separate all its elements that have a score greater than the pivot value of max(ε r , ε + l ) (corresponding to the (k/2)-mers that are in the smaller of R and L + ) from those that have a score less or equal to max(ε r , ε + l ) (corresponding to the (k/2)-mers that are only in the larger of R and L + ). Once the rearrangement around the pivot is performed, retrieving R and L + from their union is trivial.

Figure 6 presents the pseudocode of this algorithm for even values of k, where the PARTITION algorithm of quicksort [START_REF] Cormen | Introduction to algorithms[END_REF] is used to rearrange R ∪ L + , using max(ε r , ε + l ) as pivot. Note that the algorithm substitutes the top level of the recursion, and uses the divide-and-conquer from subsection 3.2 for deeper recursive calls. The CHAIN function iterates over windows of the chain starting at position j. We assume that the data structure for range product queries is precomputed beforehand. (k/2)-mers for the two subwindows are combined in a way similar to the one described in subsection 3.2.

Finally note that the Chained Windows technique above can also be adapted to the case of odd k, by splitting every window into three subwindows of sizes ⌊k/2⌋, 1, and ⌊k/2⌋ respectively, meaning that chains will now contain windows that are ⌈k/2⌉ sites apart from each other. We also note that the technique could in theory be adapted at every recursion level, so that only a single call to DC(j, h, ε ′ ) is performed for each valid pair (j, h), with ε ′ set to the minimum value across all possible sub-windows from which the call to DC(j, h, ε ′ ) could be executed. We leave a more thorough investigation of this idea for future work.

EXPERIMENTS

We implemented the described algorithms (https://github. com/nromashchenko/pk-algs) as part of the new software IPK (https://github.com/phylo42/IPK) and ran them on simulated and real-world data, using an Intel(R) Xeon(R) Fig. 6: Divide-and-conquer with Chained Windows for even k. This pseudocode is executed on the entire probability matrix P , and not on a window-by-window basis. See Definition 1 for the meaning of global parameters P and ε.

1: for all j ← 1 . . . k/2 do 2:

L ← empty list ▷ for the previous, current, and next window of the chain for all (W -, W , W + ) ∈ CHAIN(P , j) do ▷ "look behind" and "look ahead" score thresholds 4:

ε - r ← ε/M W -(1 : k/2) if W -̸ = nil else ε 5: ε + l ← ε/M W + (k/2 + 1 : k) if W + ̸ = nil else ε 6: Z W , L ← DCCW(L, ε - r , ε + l ) 7:
end for 8: end for 9: return all lists Z W 10:

11: procedure DCCW(L, ε - r , ε + l ) 12:
Z ← empty list; swapped ← f alse 13:

ε l = ε/M (k/2 + 1 : k) ▷ Local thresholds 14: ε r = ε/M (1 : k/2) 15: 16: if L is empty then ▷ If W is the first window 17: L ← DC(1, k/2, ε l ) 18:
end if

19: R ← DC(k/2 + 1, k/2, min(ε r , ε + l )) 20: n l ← PARTITION(L, ε l ) if ε - r < ε l else |L| ▷ Find
for all (l, s l ) ∈ L[1 : n l ] do 29: for all (r, s r ) ∈ R[1 : n r ] do 30: if s l • s r ≤ ε then 31: break 32: end if 33: x ← r • 2 ⌈log 2 σ⌉⌊k/2⌋ + l if swapped else l • 2 ⌈log 2 σ⌉⌈k/2⌉ + r 34: Add (x, s l • s r ) to Z 35:
end for 36: end for 37: return Z, (L if swapped else R) 38: end procedure W-2133 CPU @ 3.60GHz (8Mb cache size) machine with 62 Gb RAM (running under Linux 5.4.0-109-generic) and GCC 9.4.0. We measured the wall-clock time spent by every algorithm to process every window of the input matrices, and the peak memory consumption while processing all matrices.

In the first experiment, we generated a thousand random matrices of one thousand positions each, as follows. Every a ∈ {A, C, G, T } for every position gets a random score from the uniform distribution over [0, 1]. Then, every column is normalized so that its values sum up to one. Note that this means that the algorithms are tested over about one million windows of size k.

In the real-world experiments, we take benchmark datasets previously used in other studies related to phylogenetic placement. Each dataset specifies a reference alignment and a reference tree. We infer two P u matrices per edge of the reference tree, as it is typically done for phylogenetic placement applications [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF] (see also Figure 1).

The first real-world dataset, neotrop [START_REF] Mahé | Parasites dominate hyperdiverse soil protist communities in neotropical rainforests[END_REF], consists of 512 eukaryotic 18S rRNA sequences of 2.8 Kbp length, resulting in 2042 matrices of size 4 × 2817 (≈ 5.7M k-wide submatrices in total). The second real-world dataset, D155 [START_REF] Linard | Rapid alignment-free phylogenetic identification of metagenomic sequences[END_REF], consists of 155 complete Hepatitis C Virus (HCV) genome sequences, of 9.5 Kbp length, resulting in 614 matrices of size 4 × 9552 (≈ 5.9M k-wide sub-matrices in total). We calculate the P u matrices using RAXML-NG [START_REF] Kozlov | RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference[END_REF].

We focus on threshold values of ε = (1.5/4) k as this is the default in RAPPAS as well as in IPK. Thus, the threshold value does not depend on the input matrix, contrary to commonly used dynamic thresholds for PSSM based on pvalues. However, it depends on the length of the k-mers computed. We run algorithms for k equal to 6, 8, 10, and 12, which are common values for processing DNA datasets for RAPPAS (whose default value of k for DNA is 8).

Running time per window as a function of the number of alive k-mers

Figure 7 shows the mean running time per window of the three algorithms we have presented here: branch-andbound (BB), divide-and-conquer (DC), and divide-andconquer with Chained Windows (DCCW), plotted against the number of alive phylo-k-mers in the window, for k = 10 and ε = (1.5/4) k (see Appendix F for other threshold values). Note that many different windows may correspond to a single value of the x-axis. Each point in Figure 7 shows the average time over all windows that happened to have the same number of alive k-mers. Both axes are in log-scale. From left to right, Figure 7 shows the plot for simulated data (Random dataset), for neotrop and for D155 datasets.

First, let us observe the relative performance of the three algorithms. In experiments both on simulated and realworld data, BB (red points) showed a better running time for phylo-k-mer-poor windows (|Z| < 25) than DC (green points). However, BB showed the worst running time for most values of |Z|. Let us now compare DC (green points) against DCCW (blue points). For nearly all values of |Z|, DCCW showed better or similar mean running time compared to DC. For real-world datasets, the gain in running time for DCCW is higher for phylo-k-mer-poor windows than for phylo-k-mer-rich windows.

As for the dependence of mean processing times on |Z|, note that if we keep k constant (as done in Figure 7), the time complexity of BB is Θ(|Z|) (because of Theorem 1, and because every element of Z is part of the output). The linear dependence of BB (red points) on |Z| is somewhat more visible in the random dataset than in the real-world datasets. As for the two divide-and-conquer algorithms, for low values of |Z|, the runtime seems to be dominated by a term that is constant in |Z|, which is consistent with the analysis provided in Theorem 2.

Interestingly, we remark a strong spread of the points for very high values of |Z| (extreme right of each panel in Figure 7), which mostly affects BB. This is due to the fact that for very large values of |Z|, only a few windows contribute to the computation of the mean processing time. For this reason, the computed means have an increasingly large variance. If we exclude large values of |Z|, a large number of windows contribute to the computation of the mean processing time for most other parts of the plot. To check this, Figure 10 in Appendix C plots the number of windows contributing to each value of |Z|. The fact that this phenomenon appears to affect BB more than the other two algorithms suggests a higher variance of its running time.

Figure 10 also allows us to appreciate the difference between the simulated and the real-world datasets. Compared to the simulated dataset, the real-world datasets (especially D155) contain an over-representation of windows contributing with a large number of alive k-mers. Despite these differences, the three panels in Figure 7 are fairly similar. Unsurprisingly, the plot for the random dataset offers a somewhat less noisy version of the other two plots.

Running time over all windows

From Figure 7, we can see that the relative performance of the algorithms is dependent on the number of alive k-mers in it. In Figure 8, we look at the overall performance of the algorithms per dataset, averaging processing times over all windows in a single dataset. This has the effect of naturally weighting the contribution of phylo-k-mer-rich and -poor windows according to their frequency. Figure 8 shows the mean running times for different values of k, which also allows us to examine their dependence on k.

With the possible exception of DC for k = 6, we note that the divide-and-conquer algorithms are faster than BB across most experiments. The speed-up of DCCW over BB varies from about 1.4x (for k = 6) to between 4.4x and 5.2x for k = 12. In all three datasets, the advantage of the two versions of divide-and-conquer for phylo-k-mer-rich regions appears to outweigh by far any potential disadvantage for phylo-kmer-poor regions. As for the dependence on k, the roughly linear plot confirms the exponential dependence of running times on k (as |Z| is typically an exponential function of k).

Memory consumption

Memory consumption of the three algorithms is very close in practice. We provide measurements and discuss them in Appendix (see section D, Table 1).

CONCLUSION AND FUTURE WORK

We have described the problem of phylo-k-mer computation and algorithms for solving it. In particular, we designed two algorithms based on the divide-and-conquer approach, one of which exploits the redundancy of adjacent probability matrix windows of the input alignment. To the best of our knowledge, these two algorithms are novel, even when considering a problem similar to phylo-k-mer computation arising in the literature about motif searches. Experiments on simulated and real-world data suggest that the new algorithms perform better than the previously known branchand-bound algorithm in terms of running time, especially when a large number of phylo-k-mers must be output. The algorithmic results presented here, paired with an effective C++ implementation, made it possible to improve the running times of RAPPAS (written in Java) by two orders of magnitude [START_REF] Romashchenko | Computing informative k-mers for phylogenetic placement[END_REF]. It makes it practical for the successor of RAPPAS (manuscript in preparation) to use parameter values that were hardly feasible before, e.g., values of k > 10. Note that all the required preprocessing steps (construction of the references and the computation of the P u matrices) are independent of k, so phylo-k-mer computation from P u is indeed the bottleneck here.

One direction for future research is to exploit the relationships between nodes encoded in the phylogeny. Let u, v be tree nodes that are closely located in the reference tree (e.g., in terms of the length of the path separating them). Then, the corresponding probability matrices P u , P v can also be expected to be close to each other in terms of probability values, potentially giving rise to similar sets of phylo-k-mers. The question is thus to find an algorithm that, instead of computing from scratch the phylo-k-mers for node v, updates the list of phylo-k-mers computed for node u by taking advantage of the proximity of P u and P v .

program under the Marie Skłodowska-Curie grant agreement № 956229. FP is funded by the ANR via the project Fish-Predict. to avoid rearranging the same (k -1) columns in memory when switching from one k-sized window to the next one. Figure 9 presents the mean time spent by the branchand-bound (BB) and its modified version (BB sorted) to perform phylo-k-mer computation for different values of k on the random and the neotrop datasets. We could not obtain any significant gain in speed by sorting columns, unlike that in [START_REF] Martin | Fast and accurate genomescale identification of DNA-binding sites[END_REF]. This is likely due to the fact that, for phylo-k-mer computation, we typically deal with a larger proportion of alive k-mers in a window than in PSSM applications. In other words, it shows the distribution of window "phylok-mer richness" for different datasets. Note that for realworld datasets (neotrop, D155), distributions have multiple peaks in the mid-range, unlike the one of the random dataset. This is due to the biological nature of these datasets: for this data, the probability distributions of the columns of the corresponding matrices P u are indeed almost guaranteed to be not uniform. This figure highlights that the impact of more phylo-k-mer-rich windows on the overall computation time is higher for real-world datasets than for simulated data.

APPENDIX C PHYLO-k-MER RICHNESS

APPENDIX D MEMORY CONSUMPTION

To evaluate and compare the memory requirements of the presented algorithms, we measured the peak RAM TABLE 1: Peak memory consumption (maximum resident set size in Megabytes) of the process performing the computation of phylo-k-mers for all input matrices of a given dataset using each of the presented algorithms. Every value is the average of measurements for three independent runs. consumption as follows. For every algorithm, we ran an individual process that performed reading input data for a given dataset (or simulating input data) and phylo-k-mer computation (for k = 10 and the default threshold value) for all windows of all input matrices. We measured the maximal resident size reached in the process's lifetime using GNU time. We ran every process three times to average the measurements.

The resulting values (shown in Table 1) are virtually identical for different algorithms. While BB showed the best numbers in all experiments, the degradation of DC's and DCCW's memory consumption is under 0.03% compared to BB. This can be explained by the fact that, for all algorithms, memory consumption is dominated by the size of the input and output. For the input, we keep all matrices P u in memory to optimize the overall computation for speed regardless of which algorithm is used. The output is accumulated across multiple windows of P u , as required by Definition 1.

Definition 1 (P

 1 Phylo-k-mer Computation). Input: An integer k > 1; a σ × m probability matrix P ; a threshold value ε ∈ [0, 1). Output: All pairs (w, S(w)) | w ∈ Σ k : S(w) > ε , where S(w) := wj ,l+j-1 .

Fig. 2 :

 2 Fig. 2: Depth-first branch-and-bound algorithm 1: procedure BRANCHANDBOUND(i, j, p, s) ▷ Considers extending prefix p of score s by the i-th character of the alphabet 2:

Theorem 1 .

 1 Depth-first branch-and-bound runs in O(k • |Z|) time for one window of k columns.

Theorem 2 .

 2 Divide-and-conquer runs in O(kσ k/2 + |Z|) time for one window of k columns.

Fig. 3 :

 3 Fig.3: Divide-and-conquer algorithm 1: procedure DC(j, h, ε ′ ) ▷ Lists all the h-mers reaching the score of ε ′ in a window starting at site j 2:

Fig. 4 :

 4 Fig. 4: Illustration of the work required by each recursion call of the algorithm in Figure 3. If we exclude the root, the sum for the remaining nodes is O(k • σ k/2 ).

  Now note that there are at most 2 d calls of depth d, which gives us a total of O( log k d=1 2 d σ k/2 d ) for all depths (excluding the top-level call; see the rightmost column in Figure 4). Because 2 d ≤ k, we can write this as O(k • S), where S = σ k/2 + σ k/4 + σ k/8 + • • • + σ 2 + σ. To conclude the proof, note that

Fig. 5 :

 5 Fig. 5: Illustrations for the divide-and-conquer algorithm with Chained Windows for even k. (a) For k = 4 and σ = 4, two windows W and W + at a distance of (k/2) = 2 from each other share (k/2) = 2 columns (gray area). Thus, the (k/2)-mers alive in W R or W + L can be computed with a single recursive call. (b) An example of three chains of windows for k = 6 and σ = 2. The arrows indicate the starting positions of the different windows within the same chain. The curly braces indicate the first two windows of the first chain. In this example, all possible windows are covered with three chains.

3 D155Fig. 7 :

 37 Fig. 7: Average time in to process a window of the alignment plotted against the number of phylo-k-mers alive for k = 10 and ε = (1.5/4) k , for the three algorithms considered here: branch-and-bound (BB), divide-and-conquer (DC), and divide-and-conquer with Chained Windows (DCCW). Both axes are in log-scale.

2 D155Fig. 8 :

 28 Fig. 8: Time (in microseconds, log-scale) to process a window for different values of k, averaged across all windows encountered in a single dataset.

2 NeotropFig. 9 :

 29 Fig. 9: Mean time (in microseconds, log-scale) to process a window for different values of k by branch-and-bound on original data (BB) and on windows with sorted columns (BB sorted).

5 D155Fig. 10 :

 510 Fig. 10: Number of windows (y-axis) that have |Z| alive kmers (x-axis) for the three datasets used in experiments.

Figure 10

 10 Figure 10 plots the numbers of windows having different numbers of alive k-mers for k = 10 and ε = (1.5/4) k .In other words, it shows the distribution of window "phylok-mer richness" for different datasets. Note that for realworld datasets (neotrop, D155), distributions have multiple peaks in the mid-range, unlike the one of the random dataset. This is due to the biological nature of these datasets: for this data, the probability distributions of the columns of the corresponding matrices P u are indeed almost guaranteed to be not uniform. This figure highlights that the impact of more phylo-k-mer-rich windows on the overall computation time is higher for real-world datasets than for simulated data.

Fig. 11 :

 11 Fig. 11: Average time in microseconds to process a window of the alignment plotted against the number of phylo-k-mers alive for k = 10 and different values thresholds: (a) ε = (1/4) k . (b) ε = (2/4) k . Both axes are in log-scale.

ACKNOWLEDGMENTS

We thank Bastien Cazaux for his comments and the valuable discussion of the algorithms and the French Ministry of Research (MNERT) for a fellowship to NR. ER is funded by the ALPACA project, which has received funding from the European Union's Horizon 2020 research and innovation

Nikolai Romashchenko is a research engineer at the CNRS (French National Centre for Scientific Research) and the LIRMM (Laboratory of Computer Science, Robotics and Microelectronics of Montpellier). He received a doctoral degree in Computer Science at the University of Montpellier in 2021. His scientific interests are algorithms for bioinformatics, metagenomics, and phylogenetics, with a particular focus on alignment-free methods.

Benjamin

APPENDIX A WORST-CASE PERFORMANCE OF DEPTH-FIRST BRANCH-AND-BOUND

Here, we show a case where |Z| = Θ(k c ) for a small constant c, and depth-first branch-and-bound runs in Θ(k c+1 ) = Θ(k • |Z|). Consider the instances of the phylo-kmer computation problem with the following form: suppose the alphabet is binary and that all the columns of P are identical, with P 0,j = p > 1/2, and P 1,j = 1-p < 1/2. Since we are only interested in the behavior of the algorithm on a single window, we can assume P has exactly k columns. The score of any binary sequence w ∈ {0, 1} k is given by:

where h(w) is the number of 1s in w (or equivalently the Hamming distance between w and 0 k ). Note that S(w) is strictly decreasing in h(w). Now suppose that we set ε = S(1 c+1 0 k-c-1 ) = p k-c-1 (1p) c+1 , for some constant c. (Note that since c is constant and k is not, we can assume c ≪ k.) Then a kmer w is alive if and only if h(w) ≤ c, i.e., it has at most c 1s. Because of this,

Let us now consider the set of k-mers with h(w) = c + 1, i.e., whose number of 1s is exactly c + 1. There are exactly k c+1 = Θ(k c+1 ) such k-mers. We now prove that each of these k-mers has a different dead prefix that is visited by the algorithm: Let w be such that h(w) = c + 1 and let p w be the maximal alive prefix of w, ending with the character preceding the last 1 in w. Because p w is an alive prefix, it is visited by the algorithm, as well as its dead extension p w 1 (also a prefix of w), which however is immediately recognized as dead, as it cannot be extended in any alive kmer. Thus each of the Θ(k c+1 ) with h(w) = c+1 has a dead prefix p w 1 visited by the algorithm, and moreover all the prefixes p w 1 obtained in this way are clearly different, as p w uniquely determines w.

Because the total number of visited dead prefixes for this example is bound below by a function in Θ(k c+1 ), the running time of depth-first branch-and-bound is Ω(k c+1 ) = Ω(k • |Z|). Combining this result with the statement of Theorem 1, we obtain that on this example depth-first branchand-bound runs in Θ(k c+1 ) = Θ(k • |Z|) time.

APPENDIX B BRANCH-AND-BOUND WITH SORTED COLUMNS

We compared the performance of the branch-and-bound algorithm as described in the main text against its modified version. Namely, we sorted columns of the matrix by their entropy before the phylo-k-mer computation. Thus, branchand-bound visits the columns not in their original order but in descending order of their informativeness. To reduce the computational overhead, we did not rearrange the columns in memory but kept their original order while giving the algorithm the appropriate order to visit them. The latter is

APPENDIX E RANGE PRODUCT QUERIES

Range product queries over column-maximum values can be solved using a classic technique for range sum queries over an array. First, we store maximum values of every column in an array W max of length k. Then, we compute the array

The answer to M (j 1 :

APPENDIX F ALTERNATIVE THRESHOLDS

We include measurements of the running times of the presented algorithms for the values of thresholds different from ε = (1.5/4) k for the fixed value of k = 10. Figure 11 presents these measurements.

With the exception of Figure 11(a) for the random dataset, these plots confirm the observations made in the main text. In particular, BB only outperforms DC for phylo-k-mer-poor windows, and DCCW generally has the lowest running time. The plot for the random dataset and ε = (1/4) k is special: here the distribution of |Z| for different windows is very concentrated around its mean.


Olivier Tosoni
email: olivier.tosoni@cea.fr

Lukas Schäfer
email: lukas.schaefer@tu-darmstadt.de

High-coercivity copper-rich Nd-Fe-B magnets by powder bed fusion using laser beam method

Keywords: Additive Manufacturing, Powder Bed Fusion, Nd-Fe-B, Permanent Magnets, Coercivity

 40-µm mean-sized Nd-Fe-B powders with tailored compositions were developed for PBF-LB  4 wt.% Dy alloy composition close to sintered with 1% Cu content yields highest coercivity  Fine-grained equiaxed microstructure with almost no α-Fe can be obtained  Annealed PBF-LB-built Nd-Fe-B samples reach Br = 0.62 T and Hcj = 1790 kA.m -1

Introduction

Additive manufacturing (AM) of rare-earth (RE) permanent magnets (PM) is a relatively new research topic, aiming at enabling resource efficient production of RE-PM's with complex geometry and specific stray field distribution. RE-PM's and in particular Nd-Fe-B magnets play an important role in the energy transition and electrification of mobility, their demand will grow dramatically in the next decades [START_REF] Gauss | Rare Earth Magnets and Motors: A European Call for Action. A report by the Rare Earth Magnets and Motors Cluster of the European Raw Materials Alliance[END_REF]. Due to the concentration of the RE industry in a limited number of countries and with the growing demand driven by magnet applications required to combat climate change, the markets are prone to supply shortages, sharp price fluctuations and unwanted geopolitical dependencies, which are detrimental for the industry and more generally, the transition to a net-zero CO2 emission scenario. Besides the search for materials containing less RE and the enhancement of extrinsic properties of the existing hard magnetic materials [START_REF] Gutfleisch | Magnetic Materials and Devices for the 21st Century: Stronger, Lighter, and More Energy Efficient[END_REF] [3], net-shape processes like powder injection molding or additive manufacturing appear as a way to use magnets more efficiently without wasting material in the machining step to final geometry.

AM is by nature a near-net-shape manufacturing process generating minimal waste, and it enables much higher geometrical freedom than conventional magnet manufacturing processes of pressing and sintering or injection molding, which enable the design and manufacture of optimized magnetic topologies. However, research on AM of functional materials such as Nd-Fe-B, which has a very specific microstructure and tailored grain boundaries required to transfer the intrinsic magnetic properties of the alloy into the bulk magnet, has been only initiated in the past couple of years, as the focus of AM materials research has been previously on structural materials.

The most straightforward way to additively manufacture Nd-Fe-B magnets, which has been applied by several authors, is to use a polymer-based AM method with commercially available coercive Nd-Fe-B powder mixed with PA12 as binder to make bonded magnets [START_REF] Huber | 3D Printing of Polymer-Bonded Rare-Earth Magnets With a Variable Magnetic Compound Fraction for a Predefined Stray Field[END_REF] [START_REF] Mascheroni | Additive Manufacturing of Bonded Nd-Fe-B -Effect of Process Parameters on Magnetic Properties[END_REF]. Other authors propose a method to form an epoxy-based ink with the Nd-Fe-B powder [START_REF] Compton | Direct-write 3D printing of NdFeB bonded magnets[END_REF]. Due to the absence of grain alignment (texture), the magnetic properties remained low (remanence below 0.45 T).

Gandha et al. [START_REF] Gandha | Additive manufacturing of highly dense anisotropic Nd-Fe-B bonded magnets[END_REF] used material extrusion based big area additive manufacturing (BAAM) to produce anisotropic bonded Nd-Fe-B magnets (70 vol.% MQA-38-14 Nd-Fe-B powder, 30 vol.% nylon-12) with properties of Hcj = 1030 kA.m -1 , Br = 0.94 T and BHmax = 140 kJ.m -3 . Grain alignment was induced by applying a post-process heat treatment in a magnetic field at temperatures high enough to soften the nylon-12 matrix and allow for the Nd-Fe-B particles to align with the field. The same group then proposed a mathematical analysis using a multiphysic modelling of the alignment process in a powder/polymer mixture under magnetic field in order to help predict the magnetic properties in a functionalized AM machine [START_REF] Sarkar | In-situ magnetic alignment model for additive manufacturing of anisotropic bonded magnets[END_REF]. Although the magnetic properties are promising as they are comparable to injected anisotropic bonded magnets [START_REF] Mitarai | Development of Compound for Anisotropic Bonded Nd Magnets Using d-HDDR Magnet Powder[END_REF], the major downside of this approach is the poor dimensional accuracy and surface finish to the finished magnets.

To get the maximum properties from Nd-Fe-B magnets, it is desirable to manufacture fully dense and textured magnets (such as sintered magnets are) without any binder (such as polymer bonded magnets). From all available AM technologies, powder bed fusion using laser beam (PBF-LB) is best suited for producing fully dense metallic parts with high geometrical freedom, good dimensional accuracy and reasonable surface finish [START_REF] Zhang | Additive Manufacturing of Metallic Materials: A Review[END_REF] [START_REF] Gruber | Comparison of dimensional accuracy and tolerances of powder bed based and nozzle based additive manufacturing processes[END_REF]. Furthermore, the highly selective and localized laser melting in the PBF-LB process enables the formation of fine-grained microstructures and is well suited for processing alloys where rapid solidification is of benefit [START_REF] Liu | Dislocation network in additive manufactured steel breaks strength-ductility trade-off[END_REF], as with Nd-Fe-B. It must be noted that post process heat treatments are often required after PBF-LB to reach optimal microstructure and properties [START_REF] Goll | Additive Manufacturing of Bulk Nanocrystalline FeNdB Based[END_REF].

An early work on the PBF-LB of RE-transition metal compounds is the one from Moore et al. in 2013 where magnetocaloric heat exchangers with specific internal porosity or microchannels were produced from La(Fe,Si)13 [START_REF] Moore | Selective laser melting of La(Fe,Co,Si)13 geometries for magnetic refrigeration[END_REF]. Jacimovic et al. [START_REF] Jaćimović | Net Shape 3D Printed NdFeB Permanent Magnet[END_REF] were first to report fully dense (with no binder matrix) Nd-Fe-B magnets processed by PBF-LB. They used a commercially available gas atomized nanocrystalline RE-lean Nd(Pr-Zr-Ti-Co)-Fe-B MQP-S-11-9 spherical powder from Magnequench and claimed properties of Hcj = 695 kA.m -1 , Br = 0.59 T, and BHmax = 45 kJ.m -3 in the as-built magnets. According to the authors, the rapid solidification during the process enabled the re-formation of a nanocrystalline microstructure similar to the one of the starting powder, and the realization of hard magnetic properties. Several other groups addressed the synthesis of Nd-Fe-B magnets by PBF-LB [START_REF] Urban | Influences of process parameters on Rare Earth Magnets produced by Laser Beam Melting[END_REF] [17] [START_REF] Skalon | Influence of Melt-Pool Stability in 3D Printing of NdFeB Magnets on Density and Magnetic Properties[END_REF]. Among them, Bittner et al. [19] [20] optimized the laser melting parameters to enhance the properties of this alloy further to Hcj = 921 kA.m -1 , Br = 0.63 T and BHmax = 63 kJ.m -3 . Huber et al. [START_REF] Huber | Coercivity enhancement of selective laser sintered NdFeB magnets by grain boundary infiltration[END_REF] used grain boundary infiltration with low melting point eutectic alloys to further enhance the coercivity of PBF-LB processed MQP-S powder to 1211 kA.m -1 when infiltrated with Nd50Tb20Cu30. The disadvantage of this method is the additional processing step and the corresponding decrease in volume fraction of the 2-14-1 phase and hence remanence (0.39 T in their case). Volegov et al. [START_REF] Volegov | Additive manufacturing of heavy rare earth free high-coercivity permanent magnets[END_REF] used a powder mixture of 80% MQP-B powder with 20% (NdPr)3Cu0.25Co0.75 powder to realize in-situ and in one single step the infiltration during the building of the parts, and obtained an even higher room-temperature coercivity of about 1250 kA.m -1 .

Goll et al. [START_REF] Goll | Additive Manufacturing of Bulk Nanocrystalline FeNdB Based[END_REF] took the RE-lean MQP-S-11-9 composition used by multiple research groups and added excess RE to synthesize stoichiometric and over-stoichiometric powder compositions to be manufactured with PBF-LB. With 18 at.% RE and 2-step post-process annealing (600°C for 10 min, then 500°C for 60 min), properties of Hcj = 925 kA.m -1 , Br = 0.58 T, BHmax = 62.3 kJ.m -3 were achieved. Laser parameters leading to extremely rapid solidification during PBF-LB, combined with post-process annealing was found to be important in realizing high coercivity microstructure (formation and wetting of the RE-rich grain boundary phase).

An alternative approach by using RE-rich Pr-Fe-Cu-B material for PBF-LB was investigated by Schäfer et al. [START_REF] Schäfer | Qualification of Pr-Fe-Cu-B Alloys for the Additive Manufacturing of Permanent Magnets by alloy design[END_REF]. The authors show that subsequent heat treatments can significantly increase the coercivity, which is accompanied by the formation of an intermetallic Pr6Fe13Cu secondary phase. Goll et al. [START_REF] Goll | Additive Manufacturing of Textured FePrCuB[END_REF] further demonstrated the possibility to produce a textured microstructure in such Fe73.8-Pr20.5-Cu2.0-B3.7 alloy magnets by using a non-rotating laser scanning strategy. This resulted in notable grain alignment along the laser scanning direction and increase from 0.5 T (isotropic structure) to 0.67 T in remanence.

In this work, a Cu-rich Nd-Fe-B alloy and powder were developed, with a composition otherwise close to Dy-doped sintered magnets for PBF-LB processing. This work was carried out in the framework of the EIT Raw Materials funded KIC project 3DREMAG, whose main objective was the development, upscaling and demonstration of an alternative Nd-Fe-B powder for AM. Therefore, all the processing steps for the production of the powder have been voluntarily kept as close as possible to the established production routes of Nd-Fe-B magnets and a demonstration of scalability was realized (pre-industrial 63-kg batch of powder). After laser parameter optimization and post-process annealing, properties of Br = 0.62 T, Hcj = 1790 kA.m -1 and BHmax = 65 kJ.m -3 were obtained. This represents the highest ever reported coercivity for additively manufactured Nd-Fe-B; remanence and energy product remain in the range of other PBF-LB reports and lower than reports with additive manufacturing processes allowing grain alignment.

Materials and methods

Methodology

For this work, three experimental steps, which are in short single track scanning, suction casting, and the AM of simple shaped magnets, have been used to obtain additively manufactured magnets with hard magnetic properties.

The single track scanning (laser melting along one track) was carried out on strip cast flakes fixed to the building platform. Different chemical compositions typical for sintered Nd-Fe-B magnets have been used to compare the influence of stoichiometry on the melt pool geometry. This method has been used in various studies to study the laser-material interaction in the context of PBF-LB alloy research [START_REF] Yadroitsev | Single track formation in selective laser melting of metal powders[END_REF] [26] [START_REF] Aboulkhair | On the formation of AlSi10Mg single tracks and layers in selective laser melting: Microstructure and nano-mechanical properties[END_REF] [28] [START_REF] Fan | Selective laser melting of alumina: A single track study[END_REF]. It allows to study the relationship between the most influencing laser melting parameters (speed and power) and the alloy composition in terms of melt pool size/shape, solidification microstructure and defect formation such as cracking in the case of brittle materials like Nd-Fe-B before the production of powder batches. This reduces the lead-time of alloy development for PBF-LB considerably during the first iterations with the given alloy.

The second step is suction casting experiments to simply test the composition of the alloy in laboratory scale. The suction casting process is used to test compositional modifications of strip cast alloys and the impact on their magnetic properties after rapid solidification. The cooling rate is estimated to reach orders of 10 3 K.s -1 based on previous works [START_REF] Kozieł | Estimation Of Cooling Rates In Suction Casting And Copper-Mould Casting Processes[END_REF] [31] [START_REF] Pawlik | Investigation of the cooling rate in the suction casting process[END_REF]. Although it is not possible to simulate the local, complex melting and solidification behavior during PBF-LB, the suction casting is able to indicate trends based on the influence of compositional variations on the intrinsic magnetic rapidly solidified microstructure and the resulting extrinsic magnetic properties (coercivity and remanence). Additional heat treatments can be easily applied to simulate the post process heat treatment for the AM parts.

These preliminary experiments allowed to identify an Nd-Fe-B alloy composition suitable to the PBF-LB process by performing dedicated small-scale experiments. The last step was therefore the larger scale (several kg) production of powder with this composition to additively manufacture parts with varying parameters to optimize the building conditions and understand the relationships between microstructure and macroscopic properties.

Experimental

Single-track scanning was conducted on 200-300 µm thick strip cast flake substrates having different chemical compositions typical for sintered Nd-Fe-B magnets. The six compositions that were used are summarized in the Table 1. Especially the varying proportion of the REE Nd and Dy as well as the total fraction of the REE and the boron level in the composition are known to strongly influence the magnetic properties. The flakes were mounted on the building plate as shown in Figure 1a. Their surfaces were wiped with Ethanol before processing to remove any dirt and grease. The laser melting experiments were performed using an SLM 125 HL PBF-LB system from SLM Solutions GmbH, having a IPG-YLR-400-SM fiber laser producing 400 W maximum nominal output power at 1070 nm wavelength. The x-y scanning on the processing plane is performed with a galvanometer scanner (intelliSCAN 25 from SCANLAB GmbH). After all the optics, the laser spot size is 80 ± 6 µm at focal plane with a Gaussian intensity distribution. Argon 5.0 was used as the shielding gas in the processing chamber, the atmosphere having oxygen content of < 0.1% during laser melting. Before the actual single tracks (as indicated with lines in Figure 1a), a square area of 10 x 10 mm 2 was laser-melted to have a more representative substrate microstructure to the actual situation of PBF-LB, where the previous laser-melted layer acts as the substrate. This has been shown to affect the single track studies [START_REF] Antikainen | Single-Track Laser Scanning as a Method for Evaluating Printability: The Effect of Substrate Heat Treatment on Melt Pool Geometry and Cracking in Medium Carbon Tool Steel[END_REF]. 21 single tracks (6 mm length with 0.5 mm spacing) for the actual analysis were laser melted on each substrate, perpendicular to the scanning of the previously laser-melted surface. The laser power was varied between 50 and 200 W and scanning speed between 500 and 2000 mm.s -1 in the single scan tracks. The test setup of the building platform is shown in Figure 1a for three mounted strip cast flakes. After laser melting, the flakes were cross-sectioned perpendicular to the melt tracks, mirror polished and etched with Nital to reveal the melt pool boundaries and microstructure. Optical microscope images of the melt pool cross-sections were acquired using ZEISS Axio Observer Inverted Microscope at a magnification of 50x. The ImageJ software was used to measure the melt pool dimensions from the images. After that, around 0.5-1.0 mm of the material was ground, polished, etched and imaged again following the same procedure to have two cross-sections of each melt track analyzed for the presence of cracks. It is however acknowledged that this is not an absolute analysis of crack presence and some cracks that do not continue along the melt track for a considerable distance will remain unnoticed.

Ref

The suction casting experiments have been performed in a Buehler MAM-1 Arc Melter. The sample chamber was purged and flushed with Argon several times to prevent strong oxidation of the sample material. In the Cu cooling plate, a Cu mold is inserted with a 4 mm times 0.5 mm wide cavity. Attached to the Cu mold is a secondary vacuum chamber which induces the suction of the liquid metal into the Cu mold; the depth of the cavity is several centimeters, allowing all the material to be solidified before it reaches the opposite edge. The obtained suction cast plates were annealed at 1075°C for 30 minutes, then at 750°C for 1 h and finally at 500°C for 1 h, with a cooling speed of 300°C.h -1 between the plateaus. The temperatures and holding times for the heat treatment follows typical post processing treatment of Nd-Fe-B permanent magnets produced by powder metallurgical routes.

To prepare PBF-LB experiments, an 18 kg batch of alloy with the nominal composition Nd12.1Pr0.8Dy1.5Febal.B6.3Co1.24Al1Cu1Ga0.2 was cast using a strip casting setup from ULVAC. The molten materials were induction heated up to 1450°C, then poured through two tundishes on a cooled rotating copper wheel having a linear speed 1.2 m.s -1 , to obtain ribbons with a 300 µm mean thickness. The ribbons were then submitted to a heat treatment at 1050°C under vacuum for 2 hours to grow the grains prior to hydrogen decrepitation (HD), and jet milling. This additional step allows avoiding the generation of excessively large amounts of fine particles due to the typical Nd-rich spacing of Nd-Fe-B strip cast ribbons being by one order of magnitude smaller than the targeted granulometry for AM; fine particles are known to decrease the flowability and spreadability of powders on the building platform. Hydrogen decrepitation was performed under 800 mbar H2 pressure at 250°C: these conditions allow to reduce the amount of absorbed hydrogen compared to standard hydrogenation conditions (a few bar at room temperature) in order to adapt to the coarser granulometry needed for PBF-LB compared to sintering (medium size around 40 µm vs. 5 µm). Partial dehydrogenation was obtained by exposing the decrepitated material to 550°C under primary vacuum for 3 hours. Jet milling was performed with a 3-nozzle opposed fluidized bed AFG100 setup from Hosokawa using a 3.5 bar N2 pressure and 2.5 mm diameter nozzles. The granulometric parameters of the obtained 17 kg powder batch, measured through an in-line analysis with a laser diffraction granulometer from Malvern Analytics, are the following: d10 = 18 µm, d50 = 42 µm and d90 = 85 µm. A micrography of the powder and its granulometric distribution are shown in Figure 2, confirming the very low amount of fine particles (below 10 µm). Eventually, the powder was sieved (mesh = 50 µm) in order to remove bigger particles before being used for PBF-LB experiments, with a yield above 99%. This yield is higher than one would expect according to the granulometry of the powder, which is explained by the observed irregularity in the particles' shapes.

Finally, 10 mm diameter and 5 mm high cylindrical samples were produced by laser powder bed fusion under inert atmosphere (Ar 5.0, O2 content in the building chamber below 1000 ppm) using a ProX200 AM machine from 3DSystems. The arrangement of the AM cylindrical samples on the building platform is shown in Figure 1b. The laser wavelength was λ = 1064 nm, the spot size was equal to 70 µm; laser power, laser speed and hatching distance were varied. The scanning pattern consisted of parallel lines, scanned in only one direction (no power during the return), with a rotation of 90° between two successive layers. Stainless steel substrates were used without any pre-heating.

A few experiments were conducted with cubic instead of cylindrical samples, and some of them without rotation of the scanning lines between successive layers.

The non-spherical powder has a relatively poor spreading behavior compared to commercial spherical powders for other alloys, but the narrow granulometry and especially the low proportion of fine particles much below 10 µm allowed to process it in an equipment where the powder is not conveyed through narrow nozzles that might get clogged. Still, the spreading was improved thanks to surface preparation (roughening) of the substrate, allowing to increase the adherence of the powder to the substrate and avoid depletion. After the construction, parts were removed from the substrate by a sharp knock at the base, and their circular faces were manually rectified by polishing to make them parallel, yielding pieces of roughly 4 mm height for magnetic characterization. Annealing of the AM magnets was performed under secondary vacuum at temperatures between 400°C and 700°C, with a cooling rate of 30°C.min -1 after each plateau, here again following typical post processing treatment of Nd-Fe-B permanent magnets produced by powder metallurgical routes. Density measurements were performed with the Archimedes method in pure ethanol on as-build and polished samples. J-H curves were obtained with an AMH-300-P hysteresisgraph from Laboratorio Elettrofisico. Scanning Electron Microscopy (SEM) was carried out with a XL30 tungsten filament microscope from Philips and a MERLIN Field Emission Gun (FEG) microscope from Zeiss, both equipped with a back-scattered electron (BSE) detector for chemical contrast. Prior to microstructural characterization, samples were embedded in an epoxy resin and their surface was mirror polished with a diamond solution down to 1 µm roughness. The measurement of the oxygen content was performed by the Instrumental Gas Analysis method with a commercial equipment (EMGA-830 by Horiba).

(a) (b)

Results and discussion

Screening of alloys by single track melting experiments

The single track screening parameters were chosen based on prior experience of suitable energy density limits in PBF-LB of Nd-Fe-B. Figure 3a shows the cracks observed for all the studied laser parameter and alloy combinations and Figure 3b the number of cracks versus the RE content of the alloy. We have further marked in Figure 3a whether the cracking was observed in the heat-affected zone (HAZ), in the weld (WELD) or whether cracks were seen on both (BOTH), or if no cracks were detected (-). It can be seen that both type of cracks were observed in most of the alloys. As already reported [START_REF] Skalon | Influence of Melt-Pool Stability in 3D Printing of NdFeB Magnets on Density and Magnetic Properties[END_REF] [34], cracking is always present to some extent, and increases with the energy density.

The only alloy composition where cracks in the weld metal were not observed with any of the studied parameter combinations (which, as explained in the previous section, does not mean there was no cracking at all) was the Alloy E. It can be further seen that this alloy composition has the most parameter combinations that do not result in cracking in the heat-affected zone either (i.e. this alloy has the widest processing window).

Figure 4 shows the etched cross-sections of two single tracks on the Alloy E composition, produced with high (Figure 4a, 0.15 J.mm -1 ) and low (Figure 4b, 0.06 J.mm -1 ) linear energy input, resulting in significantly different melt pool sizes. With high energy input (Figure 4a), large melt pool is formed, with a relatively slow cooling rate resulting in mainly dendritic solidification microstructure. Dendrites as long as 65 µm grow from the edges towards the center of the melt pool. 89% of the total melt pool area as measured from the cross-section has a dendritic structure, the remaining being equiaxed in the center. With low energy input (and hence small melt pool, Figure 4b), only a small region (11% of total melt area) at the edges of the melt pool solidifies dendritically (5-6 µm long) as seen in Figure 4b, while the majority of the melt pool solidifies in an equiaxed grain structure due to extremely rapid cooling rate. It can be also seen that the Alloy E composition has a distinctively different behavior from the other studied compositions, namely that the effect of the laser parameters on the formation of dendritic or equiaxed structure is most pronounced.

Figure 5a shows the melt pool penetration depth and Figure 5b the width of the molten track as measured from the etched cross-sections for all the studied compositions as a function of laser line energy. Here again, it can be seen that the Alloy E composition distinguishes from the other ones by having the largest melt pools, which is related to thermal conductivity but can also be affected by the metastable melting point of the material [START_REF] Fan | Selective laser melting of alumina: A single track study[END_REF]. This may explain the larger share of dendrites observed in the melt pools compared to the other compositions, as a larger melt pool would result in slower solidification rate and hence dendrite growth is promoted.

As stated before, the aim of this screening is not an optimization of the mechanical properties of the magnets, nor the complete suppression of cracks, but a comparative evaluation of the behavior of different Nd-Fe-B alloys in order to identify alloys that should be more prone to processing by PBF-LB. Cracking needs to be minimized but a stronger emphasis is made on magnetic properties in this work, provided that samples are mechanically stable enough for magnetic measurements. At this stage, Alloy E distinguishes from the others and seems a good candidate. Original data is provided in the Supplementary Material section (Table 3) 3.2 Confirmation of the composition choice and investigation of the copper addition by suction casting Three of the aforementioned alloys have been tested by suction casting, namely the alloys D, E, F, which differ mainly by their dysprosium content: no Dy for the Alloy D, 3.7 wt.% in the alloy E, and 10.1 wt.% in the Alloy F. The suction casting reproduces a form of rapid solidification whereby cooling rates of 10 3 K.s -1 can be achieved [START_REF] Kozieł | Estimation Of Cooling Rates In Suction Casting And Copper-Mould Casting Processes[END_REF] [31] [START_REF] Pawlik | Investigation of the cooling rate in the suction casting process[END_REF]. Although this technique cannot replace the optimization process and the local, complex melting and solidification processes, as well as the high cooling rates of up to 10 6 K.s -1 during PBF-LB [35] [36] [37] [START_REF] Buchbinder | High Power Selective Laser Melting (HP SLM) of Aluminum Parts[END_REF], these experiments only serve as a first test to assess the potential of the alloys and the copper addition regarding their magnetic properties, especially their coercivity.

Speed

The Dy content is known to have a positive impact on the coercivity of a magnet [START_REF] Sagawa | Nd-Fe-B permanent magnet materials[END_REF]. Yet, contrarily to expectations, the Alloy F with highest Dy content shows the lowest coercivity after the suction casting experiment Figure 6a. It must be mentioned that the unexpected low coercivity for Alloy F with high Dy content can be explained by its microstructure which shows a smaller fraction of RE-rich grain boundary phase than the other alloys, as depicted in the SEM-BSE images (Figure 6b,c,d). Increasing the RE content has a beneficial effect on the coercivity by forming an intergranular matrix and separating the hard magnetic (Nd,Dy)2Fe14B grains. The fraction of RE rich phase in Alloy F appears to be undeniably lower in comparison to the other alloys under investigation, which explains the low coercivity despite the high Dy content. The Alloy E with a moderate Dy addition provides the highest coercivity and is used for further compositional modifications.

Starting with the selected composition, the composition was varied by increasing the Cu content by 0.4 wt.%, 0.6 wt.% and 0.8 wt.%. It is known that Cu improves the wettability of the hard magnetic phase by forming low melting eutectic with REE and, moreover, decreases the grain size in Nd-Fe-B based magnets [START_REF] Fidler | Overview of Nd-Fe-B magnets and coercivity[END_REF] [41] [42] [START_REF] Li | The role of Cu addition in the coercivity enhancement of sintered Nd-Fe-B permanent magnets[END_REF]. Pieces of the strip cast material were weighed and Cu was added systematically by arc melting the mixture 2 times, followed by another melting and suction casting. The magnetic hysteresis was then measured after the annealing cycle described in the Experimental section.

Beside the decreasing saturation magnetization, the magnetic measurements show the highest coercivity value for a Cu addition of 0.8 wt.%, as shown in Figure 7a. This measured coercivity of ~0.5 T is remarkably high for cast material and emphasizes the important effect of Cu on the magnetic properties of Nd-Fe-B magnets. The microstructure shows small grains of (Nd,Dy)2Fe14B particles which are well separated by a RE rich grain boundary phase (Figure 7b). Similar to sintered Nd-Fe-B magnets, the coercivity can be explained by the fine grain size (less than 10 µm) and magnetic decoupling of the (Nd,Dy)2Fe14B grains.

Additive manufacturing of magnets by PBF-LB

Following these preliminary experiments, additively manufactured batches were realized by PBF-LB as described before, using the Alloy E composition with an additional 0.8 wt.% Cu, yielding the aforementioned molar stoichiometry Nd12.1Pr0.8Dy1.5Febal.B6.3Co1.24Al1Cu1Ga0.2. In these experiments, the layer thickness was kept constant but the laser speed, laser power and hatching distance were varied. Even though the level of cracking could be reduced after the study presented in Section 3.1, all samples present some cracks. This is reported in many publications on PBF-LB of Nd-Fe-B alloys [START_REF] Jaćimović | Net Shape 3D Printed NdFeB Permanent Magnet[END_REF] [21] [13] [START_REF] Goll | Additive Manufacturing of Textured FePrCuB[END_REF] [44] [START_REF] Bittner | Microstructure and magnetic properties of Nd-Fe-B permanent magnets produced by laser powder bed fusion[END_REF] and it is due to the brittle behavior of the intermetallic material [START_REF] Bittner | Laser powder bed fusion of Nd-Fe-B permanent magnets[END_REF]. The analysis of the J-H loops of the additively manufactured magnets reveals surprisingly large differences in the magnetic properties depending on the operating parameters. Figure 8a shows how the coercivity of the as-built magnets evolves with the laser speed for several values of laser speed: for speed values below a certain threshold, the coercivity is very low (< 200 kA.m -1 ), whereas above the threshold, the coercivity increases dramatically to reach values as high as 1200 kA.m -1 . The threshold value increases in proportion to the laser power, which suggests that the relevant parameter may be the volumetric energy density (VED) rather that speed or laser power by itself. We define the volumetric energy density as 𝑉𝐸𝐷 = 𝑃 𝑣. 𝑒. ℎ where P stands for the laser power, v for the laser speed, e for the layer thickness and h for the hatching distance. The relevance and the limits of VED as a synthetic parameter has been discussed in literature [START_REF] Bertoli | On the limitations of Volumetric Energy Density as a design parameter[END_REF]; in particular, VED does not capture the complexity of the physics of PBF-LB and cannot be considered as an absolute indicator. Still, within certain boundaries, it may remain an interesting hint to target a parameter window that is likely to yield good material properties. Figure 8b displays the coercivity and the relative density versus the VED for about 30 samples synthesized with moderate laser power (100 W < P < 250 W) and constant layer thickness (e = 30 µm), with varying laser power, laser speed and hatching distance. Relative density is consistently satisfying (> 94%) for VED values above 37 J.mm -3 , whereas coercivity is high (> 800 kA.m -1 ) for VED values below 46 J.mm -3 , thus defining a preferential parameter window. For example, for a given laser power P, the scanning speed must be above P/(46 J.mm -3 .e.h) to obtain high coercivity. Figure 11 shows demagnetization curves of samples built with different VED values. It should be mentioned that samples with VED below 30 J.mm -3 were not successfully built and do not appear in the data and that samples with coercivity above 1 T (800 kA.m -1 ) all have a remanence comprised between 0.55 T and 0.60 T. Because no special means were applied to obtain grain alignment (texture), the maximum remanence value that one can expect is half of the saturation magnetization of the alloy. The alloy that was used for these experiments has a high RE content and high copper content; this decreases to the amount of Fe available to form the magnetic phase Nd2Fe14B, which has a saturation magnetization of 1.6 T. With this composition, only about 80 % of the volume can consist of magnetic phase, yielding a saturation magnetization of about 1.3 T. Therefore, the remanence values that were obtained with these high-coercivity samples are close to the physical limit, which is about 0.65 T for an isotropic structure. On the contrary, samples with low coercivity (Hcj << Br/µ0) were measured with artificially low remanence, which is due to the self-demagnetizing of the samples between the pulsed 6 T magnetization and the close-loop J-H measurement. 4)

The comparative analysis of the microstructure of low-coercivity and high-coercivity samples reveals significant differences. In the samples with low coercivity (Figure 9a), α-Fe dendrites are numerous and clearly visible all over the sample, gathered along curved lines corresponding to the bottom of the melting pools. In the vicinity of these α-Fe dendrites chains, large zone with columnar magnetic phase grains can be found almost systematically. The remaining zones present an equiaxed finegrained microstructure where the typical grain size is below 1 µm. This microstructure type is known to provide high magnetic hardness and has already been observed in PBF-LB AM magnets [START_REF] Goll | Additive Manufacturing of Bulk Nanocrystalline FeNdB Based[END_REF]. However, when its proportion in the magnet volume is too low, and when soft magnetic iron is ubiquitous, the magnetic hardness is lost at the macroscopic scale and the coercivity remains very low. On the contrary, samples that have a high coercivity in the as-built state contain almost no α-Fe dendrites and most of the volume is occupied by equiaxed fine grains (Figure 9b). These differences in the microstructure are due to the cooling rate corresponding to different VED values: if the VED is just sufficient to melt the last powder layer with only shallow re-melting of previous layers, the quenching of the alloy after the laser melting is very fast, which favors the emergence of a submicronic equiaxed microstructure (Figure 9c) and impedes the formation of α-Fe and the dendritic growth of the 2:14:1 phase, as shown in previous works [START_REF] Goll | Additive Manufacturing of Bulk Nanocrystalline FeNdB Based[END_REF] [START_REF] Jaćimović | Net Shape 3D Printed NdFeB Permanent Magnet[END_REF]. On the contrary, if the VED is higher, the laser provides excess energy that causes a bigger melt pool and deeper re-melting, and therefore a reduction of the thermal gradient and cooling rate, thus enabling the columnar growth of the 2:14:1 phase and giving sufficient time for iron dendrites to form and dendrites to grow. These trends can be to some extent influenced by compositional changes of the alloy, and this work shows that a 0.8wt.% copper addition (i.e. about 5-8 times more copper than in standard sintered Nd-Fe-B compositions) is helpful to promote the small-grained equiaxed microstructure. The oxygen content of the additively manufactured samples is about 1400 ppm, whereas the oxygen content of the as-cast ribbons is about 200 ppm. This moderate oxygen uptake, comparable to the sintering process, could be achieved thanks to the low amount of fine particles in the powder combined to the good tightness of the PBF-LB equipment and transfer vessels.

Annealing of the AM magnets was carried out after preliminary differential scanning calorimetry experiments, which allowed to estimate that a liquid phase appears at about 440°C for this composition (Figure 10). The best annealing condition generally comprises a last plateau at a temperature slightly above the temperature of apparition of the first liquid phase of the alloy. The annealing study was focused on samples that already had a high coercivity in the as-built state. Figure 11 shows demagnetization curves of as-built and annealed samples. Remanence does not exceed 0.65 T, which is the maximum value expected for Nd-Fe-B alloys without grain alignment. The absence of texture was confirmed by magnetic measurements done on other samples which had a cubic shape, and which we measured along the three axes, obtaining exactly the same J(H) curves in the three directions. In parallel, some attempts to obtain texture with laser patterns without rotation of the scanning direction between successive layers were carried out unsuccessfully, i.e. the remanence remained below 0.65 T.

The plateau temperature providing the larger coercivity enhancement appeared to be between 470°C and 490°C, with a coercivity as high as 1765 kA.m -1 . However, this condition does not allow to remove the few α-Fe dendrites present in the as-built state, and no obvious modification of the microstructure can be pointed out with standard SEM microscopy. The dissolution of α-Fe dendrites can be achieved with an annealing plateau at 600°C, but its duration has to remain very short in order to avoid grain growth. For this reason, a 10-min plateau at 600°C followed by a second plateau at lower temperature has proved to produce good results [START_REF] Goll | Additive Manufacturing of Bulk Nanocrystalline FeNdB Based[END_REF] and was applied in this work. Eventually, the better magnetic properties were obtained with a 2-stage annealing, with a first 10-min plateau at 600°C followed by a 1-hour plateau at 470°C, with a cooling speed equal to 30°C.min -1 after each plateau. Figure 9d displays the microstructure of this sample with the optimal annealing cycle, showing the complete removal of α-Fe dendrites and the absence of noticeable grain growth.

Conclusion

In this study, we have presented a method for screening the alloy composition of Nd-Fe-B magnets for PBF-LB using single track scanning and suction casting, to ultimately produce magnets using PBF-LB from non-spherical jet milled powder. Using this approach, a copper-rich Nd-Fe-B composition was synthesized by a process close to the industrial standards (strip casting, annealing, hydrogen decrepitation, nitrogen jet milling), and high coercivity isotropic Nd-Fe-B magnets realized from this powder using PBF-LB. Magnetic properties (after annealing) of Br = 0.62 T, Hcj = 1790 kA.m -1 and BHmax = 65 kJ.m -3 were obtained, showing the possibility to produce magnets with PBF-LB having high coercivity comparable to sintered magnets. This represents the highest reported coercivity for additively manufactured Nd-Fe-B magnets as of today ( Figure 12). A relatively low energy input (VED = 38-46 J.mm -3 in this study) during PBF-LB processing was identified necessary to realize high coercivity, as low energy input leads to a small melt pool size, extremely rapid cooling and hence a fine, equiaxed microstructure without dendrite growth or excessive α-Fe is formed. Further studies should focus on inducing grain alignment to PBF-LB processed Nd-Fe-B magnets to reach higher remanence and energy product.

Figure 1 :

 1 Figure 1: Views of the experimental building platform setup. (a) SLM 125 setup used for single laser track experiments on three flakes with marked examples for the tracks. (b) ProX200 setup for the AM of cylindrical parts with the powder of the selected alloy.

Figure 2 :

 2 Figure 2: (a) SEM picture of the Nd-Fe-B powder produced in this work and used for PBF-LB experiments in the ProX200 machine; (b) Volumetric particle size distribution of this powder, with d10 = 18 µm, d50 = 42 µm, and d90 = 85 µm.

Figure 3 :Figure 4 :

 34 Figure 3: Cracking tendency for the studied alloy compositions (a) as a function of the laser parameters:-corresponds to no cracks, HAZ to the presence of cracks in the heat-affected zone only; MELT in the melt pool zone; and BOTH in both zones; N/A corresponds to one isolated case where the laser vector missed the sample, so that no crack data was available.

Figure 5 :

 5 Figure 5: Size of the melt pools measured from the cross-sections (a) penetration depth (b) width of the melt pool.Original data is provided in the Supplementary Material section (Table3)

Figure 6 :Figure 7 :

 67 Figure 6: (a) J-H loops of suction cast alloys of compositions D, E, F without copper addition (b, c, d) Microstructure of suction cast alloys of, respectively, composition D with no Dy, composition E with 3.7 wt.% Dy, and composition F with 10.1 wt.% Dy

Figure 8 :

 8 Figure 8: (a) Coercivity of the as-built magnets as a function of laser speed for various laser power values: hatching distance is 70 µm and layer thickness 30 µm (b) Coercivity and relative density of the AM magnets as a function of VED; in this series, laser speed, laser power and hatching distance have been varied, the layer thickness is 30 µm. Red dashed lines indicate the upper and lower limit of the VED window leading to both sufficient density (>94% of theoretical) and high coercivity. Original data is provided in the Supplementary Material section(Table 4)

Figure 9 :

 9 Figure 9: Microstructures of dense (95%) PBF-LB processed magnets (a) with high energy density and low coercivity; some zones containing numerous α-Fe dendrites are marked by blue circles (b) with low energy density and high coercivity (c) same sample as (b) with higher magnification (d) with low energy density and high coercivity, after a two-step 600°C + 470°C thermal annealing

Figure 10 :

 10 Figure 10: Differential Scanning Calorimetry curve of a PBF-LB processed Nd-Fe-B sample with copper-rich composition, revealing the apparition of liquid phase at temperatures at 441°C (lower than standard Nd-Fe-B alloy).

Figure 11 :

 11 Figure 11: Demagnetization curves of PBF-LB processed magnets with various volumetric energy densities (in grey), and with thermal annealing or not (in color). In the legend, Pxxx stands for the laser power in W, Vxxx stands for the laser velocity in mm.s -1 , and Hxx stands for the hatching distance in microns.

Figure 12 :

 12 Figure 12: Results of this work compared to state-of-the-art Nd-Fe-B magnets synthetized by PBF-LB [15] [17] [18] [19] [20] [21] [13] [24]

Table 1 :

 1 Chemical

	.	Dy	Nd	Pr	Al	B	Co	Cu	Ga	Ti	Zr	Fe	Total REE
	Alloy A 0.51 23.3 7.38 0.19 1.04 0.94 0.10 0.16 0.05 0.02 Bal.	31.18
	Alloy B 3.03 21.4 6.53 0.19 1.07 1.04 0.18 0.15 0.10 0.03 Bal.	30.95
	Alloy C -	25.4 7.97 0.28 0.95 0.53 0.14 0.11 0.05 0.00 Bal.	33.37
	Alloy D -	22.8 7.5	0.23 0.89 1.99 0.16 0.22 -	-	Bal.	30.30
	Alloy E 3.7	21.7 6.1	0.12 1.03 1.10 0.11 0.20 -	-	Bal.	31.50
	Alloy F	10.1 15.1 5.0	<0.02 0.96 1.98 0.17 0.21 -	-	Bal.	30.20

compositions (wt.%) of the studied strip cast flakes.

Table 2

 2 summarizes the best magnetic properties of the PBF-LB processed magnets in this work.

	Best

magnetic properties of PBF-LB processed magnets with 0.96wt.% Cu content built with 38<VED<46 J.mm -3 Br [T] Hcj [kA.m -1 ] BHmax [kJ.m -3 ]

	As-built	0.58	1210	48.2
	After 1h-470°C annealing	0.61	1761	60.0
	After 10min-600°C + 1h-470°C annealing	0.62	1790	65.0

Table 2 :

 2 Magnetic properties achieved in PBF-LB 0.96wt.% Cu PBF-LB processed magnets with 38 < VED < 46 J.mm-3

Table 4 :

 4 Data relative to the Figure 8 -Intrinsic coercivity, remanence, energy product and relative density as a function of process parameters

	Scanning speed V [mm.s -1 ]	Laser power P [W]	Hatching distance H [µm]	Volumetric Energy Density VED [J.mm -3 ]	Intrinsic coercivity Hcj [kA.m -1 ]	Remanence Br [T]	Energy product BHmax [J.mm -3 ]	Relative density [%]
	250	50	70	95	193	0.45	16.9	97.4
	500	50	70	48	1115	0.59	48.9	96.7
	500	100	70	95	89	0.36	7.2	98.1
	750	100	70	63	404	0.51	26.1	97.5
	750	100	90	49	470	0.56	24.1	96.9
	1250	100	55	48	951	0.48	31.3	96.8
	1000	100	70	48	1096	0.58	46.9	94.8
	1250	100	70	38	1210	0.58	46.2	96.8
	1750	100	55	35	1063	0.54	36.7	97.4
	1500	100	70	32	1197	0.58	46.1	98.5
	750	150	55	121	51	0.25	3.1	98.2
	750	150	70	95	179	0.44	15.3	95.8
	1250	150	55	73	229	0.47	17.7	93.9
	750	150	90	74	121	0.42	6.7	96.9
	1000	150	70	71	75	0.35	5.0	97.4
	1250	150	70	57	733	0.56	38.4	97.7
	1750	150	55	52	692	0.54	34.6	96.9
	1500	150	70	48	1050	0.59	46.6	95.9
	1250	150	90	44	1219	0.57	45.6	96.5
	1750	150	70	41	1185	0.58	46.0	96.5
	2000	150	70	36	1083	0.53	37.3	98.3
	1750	150	90	32	970	0.52	30.2	94.6
	750	200	70	127	45	0.19	2.1	94.6
	1000	200	70	95	69	0.31	4.3	90.8
	1250	200	70	76	183	0.44	14.0	90.2
	1500	200	70	63	611	0.54	31.5	97.2
	1750	200	70	54	915	0.55	36.5	94.8
	2000	200	70	48	945	0.57	43.0	96.6
	1250	250	70	95	83	0.31	4.0	94.2
	1500	250	70	79	134	0.40	10.7	97.3
	1750	250	70	68	330	0.49	19.6	97.8
	2000	250	70	60	620	0.52	28.5	96.3

Acknowledgment

The authors would like to express special gratitude to Mr. Bastien Marchais (CEA) for this help in the PBF-LB experiments and to Jürgen Gassmann from the Fraunhofer IWKS Institute Hanau for supplying the Nd-Fe-B alloys. TUDa also thanks the DFG Collaborative Research Centre/Transregio 270 HoMMage Project ID No. 405553726, TRR 270 (projects A01, A10). This work was conducted under the EIT Raw Materials funded 3DREMAG project (project number 19242).

CRediT authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary material

Width of the laser track [µm]