Donatello Conte

Giuliano Grossi
email: grossi@di.unimi.it

Raffaella Lanzarotti

Jianyi Lin

Alessandro Petrini

Analysis of a parallel MCMC algorithm for graph coloring with nearly uniform balancing

We propose the analysis of a scalable parallel MCMC algorithm for graph coloring aimed at balancing the color class sizes, provided that a suitable number of colors is made available. Firstly, it is shown that the Markov chain converges to the target distribution by repeatedly sampling from suitable proposed distributions over the neighboring colors of each node, independently and hence in parallel manner. We prove that the number of conflicts in the improper colorings genereted thoughout the iterations of the algorithm rapidly converges in probability to 0. As for the balancing, given to the complexity of the distributions involved, we propose a qualitative analysis about the balancing level achieved. Based on a collection of multinoulli distributions arising from the color occurrences within every node neighborhood, we provide some evidence about the character of the final color balancing, which results to be nearly uniform over the color classes. Some numerical simulations on big social graphs confirm the fast convergence and the balancing trend, which is validated through a statistical hypothesis test eventually.

Introduction

Graph theory is a branch of discrete mathematics which also has cross-disciplinary aspects with computer science as well as other branches of natural and physical sciences. It plays significant roles in modelling real-world problems and exhibits numerous applications in Pattern Recognition, Operation Research, Chemistry, Physics and Engineering disciplines among others.

Within this theory, the graph coloring problem is one of the oldest and among the most popular constraint satisfaction problems. Basically, it consists in finding an assignment of colors to the vertices of a given graph such that no two adjacent vertices share the same color. Graph coloring is extensively used in many Pattern Recognition applications: social networks problem such as Community Identification in Dynamic Social Networks [START_REF] Tantipathananandh | A framework for community identification in dynamic social networks[END_REF], summarization of social networks messages [START_REF] Mosa | Graph coloring and aco based summarization for social networks[END_REF], Improving Friends Matching in Social Networks [START_REF] Murad | Improving friends matching in social networks using graph coloring[END_REF], and for Collective Spammer Detection [START_REF] Fakhraei | Collective spammer detection in evolving multi-relational social networks[END_REF]. It can also be used for recording medical or biometric images [START_REF] Lupas ¸cu | Semi-automatic registration of retinal images based on line matching approach[END_REF][START_REF] Cuculo | Using sparse coding for landmark localization in facial expressions[END_REF] and for finding good resource allocation scheme for device-to-device (D2D) communications [START_REF] Tsolkas | A graph-coloring secondary resource allocation for d2d communications in lte networks[END_REF][START_REF] Comi | Hardware-accelerated high-resolution video coding in virtual network functions[END_REF] in modern wireless communication systems [START_REF] Janis | Device-to-device communication underlaying cellular communications systems[END_REF].

Moreover, recast balancing as graph coloring problem frequently appears in parallel scientific computing to identify equitable subsets of independent tasks that allow for a better utilization of hardware resources, essentially by overcoming inefficiency due to small color classes [START_REF] Lu | Algorithms for balanced graph colorings with applications in parallel computing[END_REF].

A common characteristic of these problems is that the graphs have very large size, thus requiring a speed up of the traditional greedy sequential coloring heuristics [START_REF] Coleman | Estimation of sparse Jacobian matrices and graph coloring problems[END_REF] typically by introducing parallelization techniques. For instance, in [START_REF] Deveci | Parallel graph coloring for manycore architectures[END_REF] and [START_REF] Chen | Efficient and high-quality sparse graph coloring on gpus[END_REF] the authors focus on thread scalability for hundreds or thousands of threads, showing that it is possible to achieve both good performance and high quality with massive parallelism. Thus, thread scalability in parallel algorithms, namely the ability of hard-ware and software to deliver greater computational power when the amount of resources is increased, is crucial and widely exploited in this work.

In the literature, parallel graph coloring problem has been tackled by several approaches but, at the best of our knowledge, very few of them address the problem of balancing the color classes in parallel manner. One category of them is based on searching for a maximal independent set of vertices on a progressively shrunk graph and the concurrent coloring of the vertices in the found independent set. Often the independent set itself is computed in parallel using some variant of the Luby's algorithm [START_REF] Luby | A simple parallel algorithm for the maximal independent set problem[END_REF]. Examples of such approaches are [START_REF] Jones | A parallel graph coloring heuristic[END_REF]; [START_REF] Gjertsen | Parallel heuristics for improved, balanced graph colorings[END_REF]. Another category includes methods that color as many vertices as possible concurrently, tentatively tolerating potential conflicts, while detecting and solving conflicts afterwards (e.g. [START_REF] Boman | A scalable parallel graph coloring algorithm for distributed memory computers[END_REF]). Despite these solutions are effective in producing a proper coloring, generally minimizing the number of colors, they produce highly skewed color classes, undesirable for many applications, such as parallel job scheduling, that requires balancing among the classes. At the other extreme, one could search for a coloring being equitable, that is a coloring that guarantees that the sizes of any two color classes differ by at most one [START_REF] Meyer | Equitable coloring[END_REF]). This constraint is computationally very demanding and somehow too stringent for practical applications; moreover class size constraints are known to usually raise computational hardness when performing partitioning [START_REF] Bertoni | Size constrained distance clustering: separation properties and some complexity results[END_REF]. Balanced coloring relaxes the equitable constraint requiring that any two color class sizes differ at most by an integer l greater than 1. Few approaches have been proposed to tackle Balanced graph coloring (e.g. [START_REF] Robert | Parallel heuristics for improved, balanced graph colorings[END_REF]; [START_REF] Lu | Balanced coloring for parallel computing applications[END_REF]). However, the limit of these methods is still that they are intrinsically sequential thus not scalable, becoming unfeasible on large graph. A promising direction of research on graph coloring concerns the Markov Chain Monte Carlo (MCMC) methods that allow sampling from non analytic complex distributions. The idea is to define an ergodic Markov chain whose steady state distribution is defined over the set of colorings we wish to sample from. Within the framework of graph coloring using Markov chains several contributions have been proposed. In [START_REF] Jerrum | A very simple algorithm for estimating the number of k-colourings of a low-degree graph[END_REF] a simple sequential solution based on the Glauber dynamics has been adopted. The Glauber dynamics produces a Markov chain on a proper coloring where at each step a random vertex v is recolored, choosing a color uniformly at random from the permissible ones.

To deal with large graphs, in [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF] we presented an algorithm based on MCMC method producing balanced graph coloring in a parallel way. Moreover, in that preliminary work we showed the effectiveness of this stochastic approach through experiments on random graphs. In this paper we provide some advances outlined in the following points:

• we analyze the convergence the proposed MCMC algorithm on the basis of the stochastic rules used in the different stages of the overall sampling process;

• we provide an intuitive and non-rigorous analysis of the balancing mechanism, whereas we assess the quality of the final balancing through a suitable asymptotic hypothesis test of statistical fitting;

• we have conducted experiments on social graphs (extracted from real data) in order to show that both convergence and balancing properties do not depend on "regular" topologies as those exhibited by random graphs.

As for experiments, to achieve significant speedups we leverage on modern many-core GPUs architectures (NVIDIA, 2019). Numerical results also show that the number of threads used by the parallel algorithm scales well with the graph sizes, even if compared with the GPU exploitation by the greedy strategy.

The remainder of the paper is organized as follows: Section 2 recall the basic principles of the MCMC algorithm described in [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF]; the convergence analysis is given in Section 3 while in Section 4 we show numerical results on some real graphs of considerable size. Section 5 concludes the paper.

A parallel MCMC algorithm for graph coloring

In this section we briefly outline the parallel algorithm for graph coloring based on the MCMC sampling technique presented in [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF].

Notations

We will consider a simple undirected graph G = V, E with n = |V| vertices and the set [k] = {1, . . . , k} of colors used to label the vertices. A k-coloring c : V → [k], also represented as vector c = (c 1 , . . . , c n ) ∈ [k] n , is called proper if adjacent vertices receive different colors, otherwise it is termed improper. It is well-known that, if ∆(G) is the maximum degree of G, k = ∆(G) + 1 colors are sufficient to properly color the graph by a sequential greedy algorithm. For a given coloring c, let N(v) denote the neighborhood of node v in G, and c N(v) ⊆ [k] be the set of colors occupied by vertices N(v) and cN(v) its complement; let us denote the respective cardinalities with h v (c) = |c N(v) | and hv (c) = |c N(v) |. Given c, an edge uv ∈ E with c u = c v is a conflict and #(c) : [k] n → N counts the number of conflicts. We will also consider the absolute frequency of the color j in c:

f j (c) = |{u ∈ V : c u = j}|.
Hereafter we will use lowercase letters, e.g. c, c , c * , for given colorings and uppercase for random colorings, e.g. C, C , C * . For example the probability of C = c given C = c will be denoted

P(C = c | C = c) or P(c | c) for short.

Markov Chain Monte Carlo for sampling colorings

This novel parallel Markov Chain Monte Carlo (MCMC) technique relies on a 1 st -order ergodic Markov chain (C i ) ∞ i=1 visiting a sequence of (possibly improper) k-colorings c ∈ [k] n and whose stationary distribution π strongly depends on the set of conflicts involved in c. As usual, we consider the Gibbs distribution as target stationary distribution for the Markov chain:

π(c) = e -β#(c) Z(β) , with Z(β) = c ∈[k] n e -β#(c ) . (1)
The role of parameter β in (1) aims to penalize improper colorings, leading π toward the uniform distribution over the proper colorings only with exponential decrease (as β increases). It is known from MCMC theory that the latter distribution is asymptotically approached in the sampling process when the chain is suitably constructed using the well-established Metropolis-Hastings algorithm [START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF]. This construction requires the specification of a proposal probability encapsulating the acceptance ratio and a transition probability for the chain.

As for the transition probabilities of the Markov chain, given a coloring C = c we sample the successive coloring C * in two phases acting according to a typical "rejection sampling" scheme [START_REF] Voss | An introduction to statistical computing: a simulation-based approach[END_REF]: first a candidate coloring C is generated according to a suitable proposal probability r(c, c ) := P(c | c), then the proposal C is accepted effectively as successive coloring C * according to the acceptance ratio α(c, c ) := min π(c )r(c ,c) π(c)r(c,c ) , 1 :

P(c | c) :=        α(c, c ), c c 1 -α(c, c ), otherwise .
The proposal coloring C is sampled with probability r(c, c ) as follows. Each node v ∈ V is drawn independently and with identical distribution P(c v | c) of colors so that the overall proposal probability is

r(c, c ) = v∈V P(c v | c). (2)
Notice that, in the construction of the acceptance ratio also the backward probability r(c , c) is required, hence r(c, c ) is called forward probability.

The choice of the node proposal probability P(c v | c) is a key step and is hence detailed distinctly in the following subsection. It is also important to observe from the computational viewpoint that the independent drawing of all c v , v ∈ V, allows for the generation of the new coloring in a parallel manner.

Proposal distribution for new colorings

Here we specify the algorithm for the proposal distribution procedure so that the stochastic evaluations follow consequently from the analysis of the color generation. First, the behavior of the algorithm splits into two cases based on the old coloring c. When there is some conflict locally for v, namely c v ∈ c N(v) , the new proposed color C v for v shall be redrawn with the aim of reducing the possible conflicts. We draw it from the free colors cN(v) following a nearly uniform distribution of C v = j given c:

η v ( j, c) =        1-ε h v (c) k-h v (c) , if j ∈ cN(v) ε, if j ∈ c N(v) . (3)
The rationale behind such definition is that we want to generate with high probability, a color equally likely among those free, in order to aim at the balancing objective of the method. Nevertheless, we keep a negligible chance ε > 0 to pick a color that is not free, in order to widen the search space.

As for the case of no conflict for v, i.e. c v ∈ cN(v) , it is desirable to keep nearly surely the old color c v to facilitate the convergence of the algorithm, or otherwise pick another color with a small chance ε. Hence, in the case of no conflict, for the node v the proposal color C v = j given c is distributed as

ζ v ( j, c) =        1 -ε (k -1), if j = c v ε, if j c v . (4)
With the above definitions we derive the following conditional distribution for proposal color (also called forward):

P(c v | c) =        η v (c v , c), if h v (c) < k, c v ∈ c N(v) ζ v (c v , c), otherwise.
(5)

The backward probabilities

r(c , c) = P(C = c | C = c ) = v∈V P(C v = c v | C = c )
can be obtained by symmetrical reasoning, i.e. exchanging the role of c and c in the calculations outlined above. This allows to compute then the acceptance ratio α(c, c ).

The main procedural steps of the MCMC algorithm described above are sketched in Algorithm 1.

Algorithm 1: Parallel MCMC Graph Coloring Input : Graph G = V, E with n = |V|; Number k of colors; Gibbs parameter β 1 Output: Random proper coloring C ∈ [k] n 1 C ← random initial coloring ∈ [k] n while #(C) > 0 do foreach v ∈ V in parallel do 2 Calculate C N(v) 3 h v (C) ← |C N(v) | 4 P(c v | C) ← Compute according to (5) 5 C v ← Generate with distribution P(c v | C) 6 C ← Proposed coloring (C 1 , C 2 , ..., C n ) 7 r(C, C ) ← Compute forward probability 8 r(C , C) ← Compute backward probability 9 α(C, C ) ← min r(C ,C) r(C,C ) e -β(#(C )-#(C)) , 1 10 Accept C ← C with probability α(C, C )

Algorithm analysis

We now deal with the study of the behavior of the highly scalable parallel MCMC algorithm for graph coloring proposed in [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF], and outlined in the previous section. First, we claim that in the Markov chain, the color distribution stochastically converges to proper colorings by repeatedly sampling a pool of easier proposal distributions built upon the sets of neighboring node colors. Then, we also develop a discussion about the quality of balancing achieved, the latter being the main goal of this modelling.

Convergence analysis

As for the convergence properties of Algorithm 1, due to the intrinsic randomness in drawing colors we cannot guarantee that the number of conflicts #C strictly decreases at each iteration. Therefore we give a characterization of the convergence in probabilistic terms, studying a slight variant of Algorithm 1 where we bring the parameters to the extreme values: acceptance ratio α(C, C ) = 1 and ε = 0.

Lemma 1. Given current coloring C, let #W and #W * be number of nodes in G having some conflict at, respectively, any current and successive iteration of Algorithm 1, where we set acceptance ratio α(C, C ) = 1 and ε = 0. For any integer r > 1 and provided a number of colors k ≥ ∆(G) + r, the following expectation inequality holds:

E[#W * | C] ≤ ρ#W for some 0 < ρ = 1 -(1 -1/r) ∆(G) < 1.
Proof. Preliminarily, we introduce the convenient notation of indicator expression 1(cond) that takes value 1 if the condition cond is true, and 0 otherwise. The number of conflicts in the new coloring C * can be written as #

C * = uv∈E 1(C * u = C * v ), while the number of local conflicts is # v C * = u:uv∈E 1(C * u = C * v )
; similar definitions can be given for #C and #C v . Given the current coloring C we partition V into nodes with conflicts, W = {u ∈ V : # u C > 0}, and conflict-free nodes W = V W.

First, notice that the number of local conficts at any v ∈ W is simply

# v C * = u∈W:uv∈E 1(C * u = C * v ) due to u∈ W:uv∈E 1(C * u = C * v ) = 0, since when u is not in W it maintains the old color C * u = C u C * v . Hence, in the summation in # v C * each term 1(C * u = C * v
) is a Bernoullian variable with parameter

p uv := | CN(u) ∩ CN(v) | | CN(u) | • | CN(v) | ≤ min | CN(u) |, | CN(v) | | CN(u) | • | CN(v) | ≤ 1 r , conditioned on C, i.e. 1(C * u = C * v ) | C ∼ Bernoulli(p uv ).
Moreover, these terms are stochastically independent since the drawing of C * u for each neighbor u of v is done autonomously from the other neighbors. Given C, the variable # v C * is hence a sum of independent Bernoullian variables having parameters p uv with uv ∈ E, u ∈ W, namely a so-called Poisson-Binomial random variable:

# v C * | C ∼ PoissonBinomial({p uv : uv ∈ E, u ∈ W})
We define the Bernoulli variable B * v := 1(# v C * > 0) indicating whether v has some conflict in the new coloring. Thanks to the observation above, its distribution conditioned on C is easily determined as

(B * v | C) = (1(# v C * > 0) | C) ∼ Bernoulli(1-Π u∈W:uv∈E (1-p uv ))
when v ∈ W, since Π u∈W:uv∈E (1p uv ) is the probability that v does not have any conflict with the neighbors in the new coloring C * . Clearly,

B * v = 0 if v W. The number #W * = #{u ∈ V : # u C * > 0} = v∈W B *
v of conflicting nodes in the new coloring C * thus satisfies the following relationships

E[#W * | C] = E        v∈W B * v | C        = v∈W E[B * v | C] = v∈W        1 - u∈W:uv∈E (1 -p uv )        ≤ v∈W        1 - u∈W:uv∈E 1 - 1 r        ≤ v∈W        1 -1 - 1 r deg v        ≤        1 -1 - 1 r ∆(G)        v∈W B v =        1 -1 - 1 r ∆(G)        #W = ρ#W where 0 < ρ := 1 -(1 -1/r) ∆(G) < 1. That is, the claim E[#W * | C] ≤ ρ#W holds.
Directly following the previous result, we have the following convergence.

Theorem 1. The number of conflicts #C converges in probability to 0, i.e. lim t→∞ P(#(C t ) < δ) = 1 ∀δ > 0.

Proof. From Lemma 1, using the monotonicity property of the expected value we have

E[E[#W * | C]] ≤ ρE[#W].

By the Tower Rule of the conditional expectation (E[E[#W

* | C]] = E[#W * ]) we can give the bound E[#W * ] ≤ ρE[#W].
The coloring C * is the one successive to C; hence building the sequence of successive colorings C t , t = 0, 1, 2, ... by means of the algorithm above, one can guarantee E[#W t ] ≤ ρ t E[#W 0 ] for any t. Hence, applying lim t→∞ to both side, we have the convergence in expectation of the number of conflicting nodes to 0:

lim t→∞ E[#W t ] = lim t→∞ ρ t E[#W 0 ] = 0.
Now, it is well known that convergence in expectation (namely, L 1 -convergence in the proper probability space (Ω, F , P)) implies the convergence in probability, i.e. lim t→∞ P(#W t < δ) = 1 for any δ > 0.

Since the number #C t of conflicts is easily bounded above by 1 2 ∆(G)#W t , it converges in probability to 0 as well:

lim t→∞ P(#C t < δ) = 1
for any δ > 0.

Qualitative analysis of the balancing

The balancing of color classes is an important aspect but it is difficult to prove it rigorously, due to the complex role that the distributions arising from the model play in the evolution of the algorithm towards attaining the final coloring. Nevertheless, in this section we sketch a qualitative non-rigorous analysis of this process in order to provide some intuition and rationale about the nearly uniform character of final color frequencies, regardless of the topology of the graphs at hand.

To carry out this analysis, we resort to the (n, k)-Poisson Multinomial Distribution (PMD) [START_REF] Daskalakis | On the structure, covering, and learning of poisson multinomial distributions[END_REF] which is the distribution of the sum of n independent random variables supported on the set B k = {e 1 , . . . , e k } of standard basis vectors in R k , representing the k colors provided to the algorithm. We replace the representation of the color set [k] with B k for technical reasons in this section only.

At each iteration step t of the Algorithm 1, each node v ∈ V randomly draws a color with probability distribution (5). Let F (t) v ⊆ B k be the set of free colors, at step t, for node v whose color is represented by the random variable C (t) v ∈ B k , and let

p (t) v = p (t) v,1 , . . . , p (t)
v,k be the distribution of the new color for v in B k . Clearly, due to the random drawing triggered by ( 5) we have

p (t) v, j ≈        1/|F (t) v |, if e j ∈ F (t) v 0, if e j F (t) v (6)
in case of conflict for v (eq. ( 3)), or

p (t) v, j ≈        1, if e j = C (t) v 0, otherwise , ( 7
)
when v is conflict-free (eq. ( 4)). This means that the process of assigning a new color C (t+1) v ∈ B k to node v, can be expressed by a multinoulli distribution with parameter p (t) v , i.e.,

C (t+1) v ∼ Multinoulli(p (t) v ).
Notice that each node in the color C (t+1) v is drawn independently from each other.

Let X (t) = v∈V C (t) v be the number of occurrences of all colors in the graph coloring at time t. In particular, for each j ∈ [k], X (t) j represents the frequency of color e j within the coloring C (t) . Hence, the random vector X (t) is the sum of independent but not identically distributed multinoulli random variables, which follows a (n, k)-PMD, where n is the number of vertices and k the number of colors. Let ξ(x) = P(X (t) = x) denote the probability mass function of X (t) , where x is a k-dimensional vector with non-negative integer entries summing up to n. Clearly, in order to have the desired balancing, it is crucial to assign high probability ξ(x) to the elements x with equalized entries. In the following, we provide a qualitative and non-rigorous analysis of the dynamics of color assignments and its distribution evolution in order to support the experimental evidence that the color classes are quite balanced among each other. Let us investigate the following steps in the coloring process generated by Algorithm 1.

Step 1. Assume that in the first step the colors are assigned independently and uniformly at random at each node. As a consequence, the set of free colors F (1) v of v comes up uniformly among all subsets of size F (1) v . Let us observe how a new coloring is generated by focusing on a color e j through all the n independent trials, each one corresponding to a node v. To consider the overall chances e j has, we have to focus in turn on the probability sets P j = p (1) v, j : v ∈ V . In particular, for large n and under the above independence assumptions, it results that p (1) v, j is the same on average for all v due to the uniformness of F (1) v , so that e j is assigned nearly the same chance to be selected in each trial. In other words, the averages over all nodes P j should be approximately 1/k for large graph size n.

Step t. Assuming that at time t -1 the (in general improper) coloring C (t-1) is quite balanced in terms of color frequencies, we can conclude that also the new coloring C (t) could achieve a good balancing. This should happen thanks to the following two facts. First, some colors are definitely fixed since they are randomly chosen in some previous step, as described by the (7) and this contribute to uniformly spread the colors in the neighborhoods of those nodes that have to make a choice. Second, those nodes that have to update the color, are in the same condition described in Step 1, thus they repeat the same independent drawings.

Figure 1 provides an empirical demonstration of the balancing trend exhibited by the MCMC algorithm and described in the above steps for a social graph (labeled ca-coauthors-dblp) introduced in the next section on numerical simulations. In fact, the balancing quality in each repetition (light blue lines) and on average (red line) is maintained quite constant over all iterations (Step t) apart from the first iteration (Step 1) that requires an adjustment as a consequence of the random initial coloring.

Numerical simulations

In [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF] we offered a preliminary assessment of a novel parallel stochastic technique for recovering balanced colorings on big random graphs. Here we extend the evaluation to a small sample of real-world instances, where graphs describe relations among individuals, to identify differences and analogies in behavior between synthetic and real-world data. A second purpose of these experiments is to empirically validate both algorithm convergence and balancing properties in the light of the analysis delineated in the previous section.

Hence, as in [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF], we tested our GPU implementation of the presented strategy, hereafter called MCMC-GPU, against a non-parallel implementation (called MCMC-CPU) of the same algorithm which runs on a standard CPU. The purpose of this comparison is to evaluate the scalability and speed-up of the parallel approach taking the sequential implementation as baseline reference. Also, we tested the MCMC-GPU implementation against a fast and fully parallel greedy coloring strategy inspired by [START_REF] Luby | A simple parallel algorithm for the maximal independent set problem[END_REF] work, called Luby-GPU. This simple coloring strategy provides a proper coloring in very short time, at the cost of a very high class unbalance. All the algorithms were implemented 1 in C++ and both MCMC-GPU and Luby-GPU leverages the NVidia CUDA programming paradigm (NVIDIA ( 2019)). All the host code is singlethreaded and it is meant to be run on a single machine, hence no other libraries (such as pthread, C++11 threads, OpenMP, MPI, etc...) are required. Both CPU and GPU code is compiled with NVidia nvcc v10.1, using the underlying GNU GCC v5.4.0 as C++ compiler.

For running the experiment we used a workstation with 2× Intel Xeon CPU E5-2620 v3 @ 2.40GHz, 64GB of RAM, Linux 16.04 and an NVidia GTX980 GPU featuring 2880 cores and 4GB of VRAM. MCMC-CPU runs are single-threaded, hence were run on a single core of the host CPU. For Luby-GPU and MCMC-GPU, on the host the processes are singlethreaded as well, and parallelization occurs only on the GPU side with the 2880 cores.

The three coloring approaches have been tested on four realworld social graphs taken from the [START_REF] Rossi | An interactive data repository with visual analytics[END_REF] repository. They greatly vary in terms of size (number of nodes and edges) and average/maximum degree, thus providing a small but relatively diverse array of use cases. Their relevant features are summarized in Table 1. The maximum number of possible colors to be used is a parameter in our algorithm and it has to be fixed in advance. We fixed this parameter in such a way that the results are stable for each trial. Table 1 shows the values of this parameter for each graph used in the experiment.

For measuring the overall balancing quality of a coloring c, we used the following quantity, we call unbalancing index, which is closely related to the classical root mean square deviation of the color class sizes {n j : j = 1, . . . , k}:

Γ n,k (c) =         1 k k j=1 n j - n k 2         1/2 . ( 8
)
1 The software used for the experiments is freely available on the GitHub repository: https://github.com/phuselab/MCMC_Colorer Naturally, a perfectly balanced coloring satisfies Γ n,k (c) = 0.

Figure 2 top shows the average unbalancing index over 35 repetitions for each coloring algorithm applied to the 4 sample graphs. Giving to its greedy design, Luby-GPU achieved a worse balancing performance, being largely outperformed by MCMC strategy. For the latter we set a number of colors a bit greater than the one obtained by Luby (for details see Table 1). ca-hollywood-2009). We also remark that in our experiments, MCMC-GPU produces a proper coloring for each graph in less time than Luby-GPU does, showing that sometimes it is even faster than the greedy approach.

For statistically assessing the balancing, we test the uniformness of color distribution through a hypothesis test starting from the balancing index (8). Given a coloring c produced by MCMC-GPU, we consider the statistic defined as

K 2 = k Γ n,k (c) 2 n = k j=1 (n j -n/k) 2 n/k .
In fact, this expression is well known in statistical inference as the Pearson's chi-squared statistic for goodness-of-fit test (Mood et al., 1973, §IX.5.2), that is used for assessing whether the observed categorical outcomes of the experiment have frequencies (n j ) following a hypothesized categorical distribution, namely uniform distribution in our case (i.e. color class sizes n/k). With this so-called null hypothesis H 0 , it was proved using the Central Limit Theorem and a suitable geometric transformation (see [START_REF] Cochran | The χ2 test of goodness of fit[END_REF]) that such statistic K 2 has asymptotically a χ 2 (k -1) distribution (chi-squared distribution with k -1 degrees of freedom) as n → ∞: hence, χ 2 (k -1) is a good approximation of the true distribution for large number n of nodes, as in our experiments.

The key point is that the statistic K 2 is large and the corresponding p-value is small, when the produced color class sizes are unbalanced, and vice-versa.

Results of the test are reported in Table 2 for 35 runs: notice that the range of rather large p-values leads to accept the H 0 hypothesis stating the color balancing in all analyzed graphs.

Conclusions and discussion

We have studied the convergence properties of a scalable parallel algorithm for graph coloring problem based on Markov Chain Monte Carlo techniques [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF]. We also have sketched a qualitative analysis of the achieved balancing, remarking that it is as important as difficult to formally derive an analytic form for the distribution on the final colorings, due to the complex dynamics among the applied probability distributions that the model brings into play during its evolution. It would be also worth it to study MCMC models that incorporate in the target distribution not only the term promoting proper colorings (as in the current one), but also new specific terms helping in recovering suitable families of colorings not necessarily equalized in the color classes. For example, in case of application of graph coloring to the deployment of independent tasks performed concurrently on a pool of processors (typically within datacenters or on distributed architectures), in order to gain computational efficiency the target distribution should look for colorings having color classes of size multiple of the number of processors.

In [START_REF] Conte | A Parallel MCMC Algorithm for the Balanced Graph Coloring Problem[END_REF] the effectiveness of this stochastic approach in finding balanced colorings has been shown through experiments on random graphs. Here we have extended the algorithm assessment on some kind of social graphs just to show that both convergence and balancing properties are preserved although the topologies of the graphs considered are not as "regular" as those exhibited by random graphs. As for the experimental analysis, further studies on realistic applications where graph topology and size (especially big graphs) are taken into account should be conducted. For instance, it is important to stress the heuristic ability on graphs where edge density is not uniform and the presence of quasi-clique subsets interferes with the needs of MCMC model to locally find low conflict level with high probability.

Moreover, being an heuristic technique based on a number of parameters that affect the behavior and ultimately its ability to find highly balanced solutions, a both theoretical and experimental deep analysis focused on some of the most relevant parameters is worth a great attention. For instance, it is the case of the number of colors required as input to the MCMC algorithm, which gives rise to a clear trade off between convergence speed and balancing quality, as emerged since the early experiments. There are also meta-parameters related to the MCMC technique itself that affect the results, such as the parameter β appearing in the formulation of the target stationary distribution for the Markov chain (Eq. 1). Furthermore, besides its effect in penalizing improper colorings and leading π toward the uniform distribution over the proper ones, β or some function of it would also have a relevant role in defining new interesting metrics that can promote colorings that capture graph properties other than the balancing, such as those associated to planar graphs.

Figure 1 :

 1 Figure 1: Color balancing quality of a sample graph (the social graph ca-coauthors-dblp) during the iterations of the parallel MCMC algorithm for the first 10 states (colorings) visited by the Markov chain.

Figure 2 :

 2 Figure 2: Average unbalancing index (upper plot) of final colorings achieved by the three algorithms for the graphs of Table 1 and related average execution times (lower plot), expressed in seconds. Results are on 35 repetitions (standard deviation as error bar) and plotted in logarithmic scale.

Figure 2

 2 Figure2bottom shows the average computing times over the trials. Being inherently sequential, MCMC-CPU is the slowest and a direct comparison with MCMC-GPU results in an average speed-up ranging from 28× (for sc-pwtk) up to 69× (for ca-hollywood-2009). We also remark that in our experiments, MCMC-GPU produces a proper coloring for each graph in less time than Luby-GPU does, showing that sometimes it is even faster than the greedy approach.For statistically assessing the balancing, we test the uniformness of color distribution through a hypothesis test starting

Table 1 :

 1 Features of the graphs used in the experiments: number of nodes and edges (columns 2-3), maximum and average degree (columns 4-5), average number of colors found by Luby-GPU (column 6) and the parameter of maximum number of colors to be used, set for MCMC-GPU (and MCMC-CPU) (last column).

	Graph	Nodes	Edges	Max deg	Avg deg	N. colors by Luby	N. colors by MCMC-GPU
	sc-shipsec5	179.1K	4.4M	75	24	33	50
	sc-pwtk	218K	11.4M	179	51	43	85
	ca-coauthors-dblp	540K	30M	3.3K	56	337	460
	ca-hollywood-2009	1.1M	56M	11.5K	105	2209	2400

Table 2 :

 2 Test for fitting uniform color distribution on the result of 35 runs with the analyzed graphs. Ranges of observed Pearson's chi-squared test-statistics and corresponding p-values are shown.

	Graph	K 2 statistic range	p-value range
	sc-shipsec5	[19.51, 50.52]	[0.413, 0.999]
	sc-pwtk	[31.14, 64.57]	[0.940, 1.000]
	ca-coauthors-dblp	[368.09, 473.42]	[0.310, 0.998]
	ca-hollywood-2009 [2205.43, 2446.97] [0.242, 0.997]


Hugo Merbouche

Isabella Boventer
email: isabella.boventer@cnrs-thales.fr

Victor Haspot

Stéphane Fusil

Vincent Garcia

Diane Gouéré

Cécile Carrétéro

Aymeric Vecchiola

Romain Lebrun

Paolo Bortolotti

Laurent Vila

Manuel Bibes

Agnès Barthélémy

Abdelmadjid Anane

Voltage-Controlled Reconfigurable Magnonic Crystal at the Submicron Scale

Keywords: magnonics, functional oxides, reconfigurable magnonic crystal, voltage control, spin waves, frequency filtering

Multiferroics offer an elegant means to implement voltage-control and on the fly reconfigurability in microscopic, nanoscaled systems based on ferromagnetic materials. These properties are particularly interesting for the field of magnonics, where spin waves are used to perform advanced logical or analogue functions. Recently, the emergence of nano-magnonics is expected to eventually lead to the large-scale integration of magnonic devices. However, a compact voltage-controlled, on demand reconfigurable magnonic system has yet to be shown. Here, we introduce the combination of multiferroics with ferromagnets in a fully epitaxial heterostructure to achieve such voltagecontrolled and reconfigurable magnonic systems. Imprinting a remnant electrical polarization in thin multiferroic BiFeO 3 with a periodicity of 500 nm yields a modulation of the effective magnetic field in the micron-scale, ferromagnetic La 2/3 Sr 1/3 MnO 3 magnonic waveguide. We evidence the magneto-electrical coupling by characterizing the spin wave propagation spectrum in this artificial, voltage induced, magnonic crystal and demonstrate the occurrence of a robust magnonic bandgap with > 20 dB rejection.

Spin waves, the collective spin excitations in a magnetically ordered material and their associated quanta, magnons, are promising candidates to be used as information carriers for the next generation of information processing devices. [1][2][3] Among all beyond-CMOS (complementary metal-oxide semiconductor) based approaches, magnonics combines room-temperature operation at GHz frequencies, nano-scale integration, digital and analog architectures. [4][START_REF] Schneider | Realization of Spin -Wave Logic Gates[END_REF][START_REF] Zografos | Exchange -Driven Magnetic Logic[END_REF][START_REF] Ustinov | Multifunctional Nonlinear Magnonic Devices for Microwave Signal Processing[END_REF] Magnons can interact directly or indirectly with frequencies from the microwave to the optical regime. Thus, magnonics enables the coupling to other paradigms such as quantum computing and opto-electronics rendering it a rich playground to test innovative ideas to overcome the power consumption bottleneck of conventional CMOS-based electronics. [START_REF] Tabuchi | Coherent Coupling between a Ferromagnetic Magnon and a Superconducting Qubit[END_REF][START_REF] Lachance-Quirion | Hybrid Quantum Systems Based on Magnonics[END_REF][START_REF] Lachance-Quirion | Entanglement-Based Single-Shot Detection of a Single Magnon with a Superconducting Qubit[END_REF][START_REF] Kusminskiy | Coupled Spin-Light Dynamics in Cavity Optomagnonics[END_REF][START_REF] Parvini | Antiferromagnetic Cavity Optomagnonics[END_REF][13][START_REF] Chumak | Fundamentals of Magnon-Based Computing[END_REF] However, the compatibility with state of the art CMOS based technologies will still be required for any practical magnonic device. Therefore, signal transducing and processing should be advantageously performed through voltage rather than current driven signals.

Correspondingly, piezo-electric transducers and modulators have been proposed. [START_REF] Duflou | Micromagnetic Simulations of Magnetoelastic Spin Wave Excitation in Scaled Magnetic Waveguides[END_REF][START_REF] Sadovnikov | Straintronics: Reconfigurable Spin-Wave Routing in Strain-Controlled Bilateral Magnetic Stripes[END_REF] In wave-based devices, signal processing most often relies on artificial crystals. Magnonic artificial crystals (MACs) have therefore been extensively studied. [START_REF] Chumak | Magnonic Crystals for Data Processing[END_REF][START_REF] Qin | Low-Loss YIG-Based Magnonic Crystals with Large Tunable Bandgaps[END_REF][START_REF] Frey | Reflection-Less Width-Modulated Magnonic Crystal[END_REF][START_REF] Choudhury | Voltage Controlled On-Demand Magnonic Nanochannels[END_REF][START_REF] Merbouche | Frequency Filtering with a Magnonic Crystal Based on Nanometer-Thick Yttrium Iron Garnet Films[END_REF] Achieving reconfigurability has however proven to be difficult. Indeed, most practical realisations of MACs are based on a periodic structuration which prohibits any reconfiguration of the frequency response. To date, reconfigurability has only been demonstrated in a few cases relying on the magnetic stray field modulation, albeit the spatial period of the MAC is fixed by construction. [START_REF] Chumak | A Current-Controlled, Dynamic Magnonic Crystal[END_REF][START_REF] Topp | Making a Reconfigurable Artificial Crystal by Ordering Bistable Magnetic Nanowires[END_REF][START_REF] Krawczyk | Review and Prospects of Magnonic Crystals and Devices with Reprogrammable Band Structure[END_REF] Recently, theoretical propositions based on micromagnetic simulations have shown the interest of voltage-controlled reconfigurable MAC. [START_REF] Wang | Voltage-Controlled Nanoscale Reconfigurable Magnonic Crystal[END_REF] Such devices based on the switching of multiferroic domain states would simultaneously allow for electric field control and low-power operation. [START_REF] Hu | Non-Volatile Ferroelectric Switching of Ferromagnetic Resonance in NiFe/PLZT Multiferroic Thin Film Heterostructures[END_REF][START_REF] Rana | Towards Magnonic Devices Based on Voltage-Controlled Magnetic Anisotropy[END_REF] Here, we present a fully reconfigurable, voltage-controlled device which passively couples magnonics to the magneto-electric degrees of freedom in a multiferroic/ferromagnetic heterostructure (Fig. 1) and operates at 300 K. Regarding the multiferroic, we employ a 23 nm thick layer of BiFeO 3 (BFO). This material is one of the most studied multiferroics and exhibits both a sizeable magneto-electric coupling and a piezo-electric effect of 40 pm/V. [START_REF] Ramesh | Multiferroics: Progress and Prospects in Thin Films[END_REF][START_REF] Sando | BiFeO 3 Epitaxial Thin Films and Devices: Past, Present and Future[END_REF][START_REF] Park | Structure and Spin Dynamics of Multiferroic BiFeO 3[END_REF][START_REF] Garcia | Artificial Multiferroic Heterostructures for an Electric Control of Magnetic Properties[END_REF][START_REF] Spaldin | Multiferroics: Past, Present, and Future[END_REF][START_REF] Spaldin | Advances in Magnetoelectric Multiferroics[END_REF] For the ferromagnetic layer which serves as the spin wave medium, we utilize a 13 nm thick layer of La 2/3 Sr 1/3 MnO 3 (LSMO). [START_REF] Hemberger | Magnetic, and Electrical Properties of Single-Crystalline La 1-x Sr x MnO 3 (0.4 < x < 0.85)[END_REF] Based on this structure, we create our reconfigurable and tuneable crystal by using the piezo-force microscopy (PFM) technique to generate a uniformly or periodically polarized ferroelectric pattern in the BFO. [START_REF] Alexe | Nanoscale Characterisation of Ferroelectric Materials : Scanning Probe Microscopy Approach[END_REF][START_REF] Gruverman | Piezoresponse Force Microscopy and Nanoferroic Phenomena[END_REF] In turn that pat-tern modulates the spin wave propagation in the LSMO, and allows to obtain spin wave filtering at the magnonic bandgap of the crystal. We identify the creation of the bandgap -the hallmark of a magnonic crystal -by a careful comparison of the transmitted spin wave spectra whilst the adjacent ferroelectric domain state is changed between a uniform down (virgin) or up or periodically down-up ferroelectric polarization. That comparison (cf. Fig. 2) allows us to unambiguously attribute the observed frequency filtering (cf. Fig. 3) in the amplitude representation at 3.54 GHz to a bandgap in the magnonic medium due to the periodic domain configuration of BFO. In our study we show the modulation of the spin wave propagation by applying a voltage to the BFO layer which switches the ferroelectric polarization states. Our conceptual realization of such voltage control allows for other perspectives to go beyond PFM based writing. For instance, the theoretical propositions for all-solid state devices discussed by the work of Q. Wang et al. [START_REF] Wang | Voltage-Controlled Nanoscale Reconfigurable Magnonic Crystal[END_REF] can be adapted to the BFO polarisation control using an array of voltage-gates. Writing or erasing the MAC could be simply performed using voltage spikes on a desired set of electrodes.

RESULTS AND DISCUSSION

The device layout and experimental setup for the creation of a reconfigurable magnonic crystal with our multiferroic/ferromagnetic BFO/LSMO heterostructure are shown in Fig. 1.

In order to maximize the influence of the BFO interfacial polarization state we minimized the thickness of the LSMO layer to 13 nm where the growth conditions are optimized such that the dynamical magnetic properties are improved. For instance, we obtained a value of ∝ 6 • 10 -3 for the Gilbert damping parameter which is comparable to that of thicker films. [START_REF] Madon | Room Temperature Ferromagnetic Resonance in Hetero-Epitaxial BTO -BFO/LSMO Magnetoelectric Composite[END_REF] Further details on growth and the nanofabrication process are given in the methods section and section 1 of the supporting information (cf. Fig. S1). We excite and detect spin waves in the Damon-Eshbach configuration [START_REF] Damon | Magnetostatic Modes of a Ferromagnetic Slab[END_REF] by means of propagative spin wave spectroscopy [START_REF] Vlaminck | Spin-Wave Transduction at the Submicrometer Scale: Experiment and Modeling[END_REF] (cf. Fig. 1 (a)) and deduce the transmission inductance L 21 parameter from the recorded data (cf. Supporting information). In order to both write and read the ferroelectric domain configuration of BFO, we employ piezo-force microscopy (PFM) to define the out-of-plane polarization state between the microwave antennae within an area of ≈ 2×2 µm (cf. Fig. 1 x y z (d)). The number of the writing steps starting from the virgin state is labelled as P n,m , where n denotes the writing cycle defined as the sequence of virgin, up, periodic domain polarization states and m the specific configuration in that writing cycle. The antenna arms exhibit a periodicity of 500 nm. Hence, the antennae predominantly excite spin waves around K ≈ 6 rad/µm corresponding to a maximum of the antenna's response spectrum (cf. Fig. 1 (b)). Therefore, the present study focuses on a periodic ferroelectric domain pattern with the same periodicity of 500 nm. The resulting Bragg wavevector of the periodic ferroelectric pattern then corresponds to spin waves of the wave vector of 6 rad/µm, where the down and up domains have a width of 300 nm and 200 nm, respectively. Matching the periods of the antennae and the ferroelectric domains allows to set the bandgap at the most efficiently excited wavevector. Consequently, such a matching optimizes the coupling between the antennae and supports the unambiguous interpretation of the results. Please note, that at room temperature spin wave excitations in BFO occur in the sub-THz 40 band and cannot be coherently excited here. Thus, any observation of a spin wave transmission signal can be fully attributed to the spin wave excitation in the LSMO layer. state is compared to the up state. Apart from small shifts in resonance frequency, there is no change between both uniform up and down ferroelectric domain states in the BFO indicating that the PFM writing process does not affect the LSMO properties. We attribute the small frequency difference to a slight misalignment of the sample with respect to the magnetic field between different PFM writing cycles. Furthermore, since we used less averaging cycles for the S-parameters during the measurement, the noise floor of the virgin state is higher than for the uniform up (cf. Fig. 2 b) state. Hence, by the comparison of the two different uniform ferroelectric polarization states before writing a periodic structure (cf. Fig. 2 a-b) and after rewriting back to a uniform state (cf. Fig. 2 c-d), we can reliably and reproducibly switch between the various ferroelectric polarization domain states P n,m which are also temporally stable.

In order to model the dispersion relation (cf. Fig. 2 (c), red line) for our system, we utilize the approach developed by B.A. Kalinikos and A. N. Slavin [START_REF] Kalinikos | Theory of Dipole-Exchange Spin Wave Spectrum for Ferromagnetic Films with Mixed Exchange Boundary Conditions[END_REF][START_REF] Kalinikos | The Dipole-Exchange Spin Wave Spectrum for Anisotropic Ferromagnetic Films with Mixed Exchange Boundary Conditions[END_REF] by considering both an inplane (ip) and out-of-plane (oop) (uniaxial) anisotropy for spin waves in the Damon-Eshbach configuration which reads as:

ω 2 n = (Ω n +ω M (1-P nn )+ω M • N zz ) • (Ω n +ω M • P nn • sin 2 (φ-φ n )+ω M • N xx ) (1)
where

Ω n = ω Hn +ω H ip u • cos 2 (φ-φ 0 ) with ω H ip u = γH u for the in-plane (ip) anisotropy field and ω Hn = γ((H ext +H demag )+M s • L 2 ex • K 2 n ) (H ext : external field, H demag : demagnetization field, L ex : exchange length, K n = k 2
x +k 2 n , wave vector of spin wave), ω M = γM s , φ n = arctan kn kx (k n = nπ w , w: width of waveguide), φ the in-plane angle between the saturation magnetization (here along the y-axis) and the waveguide axis (x-axis), where H oop u denotes the out-of-plane anisotropy field, N xx = -H ip u /M s sin 2 (φ-φ 0 ) and the matrix elements P nn = 1-( 1Knt (1-e -Knt )) (t: ferromagnetic layer thickness). [START_REF] Kalinikos | Theory of Dipole-Exchange Spin Wave Spectrum for Ferromagnetic Films with Mixed Exchange Boundary Conditions[END_REF] The modeled dispersion is shown (red solid line) (H demag = -13.7 Oe, k x = 6 rad/µm, L ex = 5.89 nm, Furthermore, the low value of α G allows us to extend the distance between the emitting and receiving antennae for the spin wave transmission signal to the micrometre range. In a similar fashion to Ref. [START_REF] Vlaminck | Spin-Wave Transduction at the Submicrometer Scale: Experiment and Modeling[END_REF] we can both theoretically and experimentally estimate the spin wave attenuation length of our system. Theoretically, for k x ≈ 6 rad/µm, group velocity v g ≈ 350 m/s estimated from our dispersion relation and ω/2π = 3.54 GHz, we obtain L th att = v g τ ≈ 2.1 µm, where τ is the lifetime defined as 1/τ = α G (ω Hn + ω M 2 (1+N xx +N zz )). We experimentally estimate the attenuation length from the decay of the spin wave amplitude, i.e. 2|∆L 21 |/∆L 11 = e (-D eff /Latt) , where D ef f = D+D 0 corresponds to the effective distance between the antennae, with the antenna spacing denoted as D and as D 0 the width of the antennae arms, respectively. Then, we obtain L exp att ≈ 2 µm very close to the theoretically expected value which validates the overall modeling approach.

N zz = -H oop u /M s , (a) (b) (c) (d)
H ip u =
As a result, the study of the spin wave dispersion spectrum allows now to directly attribute changes in the spin wave transmission signal for periodically altered ferroelectric domains to a periodic modulation, i.e. the emergence of a bandgap (cf. Fig. 3). In Fig. 3 (a)-(b), we show the transmitted spin wave's envelope (amplitude, (a)) and phase ((b)) in the LSMO layer from the measurement for the periodic arrangement of up and down polarized ferroelectric domains in the BFO layer (green) and compare it both to the uniformly up polarized state (red) and reconfigured down polarized state (blue). At ω/2π = 3.54 GHz (cf. Fig. 3 (a)), the spin wave amplitude drops by ≈ 20 dB compared to the amplitude in the uniformly up polarized case. Further, at that frequency, we also observe an additional phase jump in the phase spectrum for the case of the periodic domain pattern. The occurrence of such an additional phase jump was not observed for the uniform ferroelectric domains of the BFO layer, where both amplitude and phase show no discontinuity over the measured The amplitudes are normalized by their mean value. The frequency filtering from the magnonic crystal emerges at ≈ 3.54 GHz where the spin wave amplitude drops by 20 dB for the periodic P 0,↑↓ state compared to the uniform P 0,↑ state. The gap vanishes for reversing the FE domain state back to a uniform down state. (b) Display of the wrapped phase data for the uniform P 0,↑ , periodic P 0,↑↓ and uniform P 1,↓ state. The magnonic bandgap corresponds to the additional phase jump at the bandgap frequency which vanishes again after resetting the FE domain to uniform down (P 1,↓ ). The data is shown for a magnetic field of 415 Oe. Measurements are conducted at 300 K. frequency window. In particular, the gap and the corresponding phase discontinuity due to the change in group velocity vanish by the following reconfiguration of the system back to the uniform P 1,↓ FE polarization state. Furthermore, both the frequency of the gap and the phase jump are field dependent. As illustrated in Fig. 4 we observe a shift to The propagative spin wave spectroscopy used here does not provide direct access to the wavevector of the excited spin wave. However, it is possible to infer the wavevector through the phase φ of the dynamical magnetization m. It depends on the wavevector via m ∝ e iωt-kx(ω) • x , i.e. is described as a plane wave, which allows us to assess discontinuities in the dispersion of a MAC. The distance from the antenna exciting the spin wave which propagates to the second, receiving antenna is D eff . As a result, the accumulated phase between x = 0 and x = D eff corresponds to φ(ω) = -k x (ω)D eff . Then, the frequency dependence of wavevector can be derived from the unwrapped phase φ (typically the VNA measures the wrapped phase) by k x (ω) = -φ/D eff +K 0 , where K 0 is an integration constant that is set by the phase of the excitation field.

However, by introducing a spacial periodicity, we expect the dispersion relation to split in different branches separated by a magnonic band gap. Interestingly in the vicinity of k x = K B (the MAC's Bragg wavevector) the phase relation established previously no longer holds due the branches degeneracy. Indeed, the reflections induced by the periodic domains strongly interfere near the band gap (Bragg conditions). Standing spin-waves are formed, thus the plane wave ansatz does not hold anymore. Additionally, the position of those standing waves can abruptly change the coupling with the antennas and the value D eff when they move from the lower to the upper dispersion branch. Hence, the step-like phase jump for the MAC (P 0,↑↓ state) is reminiscent of an additional large phase accumulation, probably much larger than 2π as the 20dB drop in the transmitted amplitude indicates. Measuring the value of this phase-jump would require a time-resolved measurement scheme. Our observations demonstrate that a change in the ferroelectric state unquestionably modulates the spin wave dispersion of the ferromagnetic layer.

However, the underlying governing physical mechanism that couples the polarisation state of BFO and LSMO is still to be identified. The most probable mechanism would imply the coupling of the ferroelectric state to the antiferromagnetic order in BFO that induces a modulation of an the interfacial exchange field. As a matter of fact BFO/LSMO heterostructures are known to exhibit magneto-electric control of the exchange bias field (up to 160 Oe at 10 K) but this field vanishes when the temperature is increased on macroscopic devices towards room temperature. [START_REF] Wu | Reversible Electric Control of Exchange Bias in a Multiferroic Field-Effect Device[END_REF][START_REF] Yi | Tailoring Magnetoelectric Coupling in BiFeO 3 /La 0.7 Sr 0.3 MnO 3 Heterostructure through the Interface Engineering[END_REF] However, given the small sizes of the written ferroelectric domains (200 nm and 300 nm) in our system, it is possible that there is a remanent exchange coupling between BFO and LSMO that is sufficient in amplitude to generate a MAC and enables the observed frequency filtering. Correspondingly, as well as a thickness dependent study of the LSMO layer in addition to the modulation of the ferroelectric polarization of the BFO and a temperature dependent study on the role of the exchange bias could give further insight into the relevance of these effects and the role of the interface for these types of oxide heterostructures.

CONCLUSION

In summary, we implemented a tuneable MAC filtering scheme based on BFO/LSMO thin film heterostructures. Our approach allows us to reliably and reversibly filter certain frequencies in the spin wave dispersion, analogous to magnonic crystals, by controlling the ferroelectric domain state order of the BFO layer. Once written in a certain domain pattern., no further energy input is required in the MAC filtering, which results in minimal power consumption for the operation of our device. Furthermore, the writing voltage does not exceed 5 V. Thus, the introduction of multiferroics, reconfigurable and low power voltage controlled materials into microscopic magnonic systems represents a key step towards large scale integrated magnonic systems.

METHODS

Growth

Both LSMO and BFO are perovskites with well-matched unit cell parameters. They are epitaxially grown on an NGO substrate using pulsed laser deposition (PLD). [START_REF] Béa | Tunnel Magnetoresistance and Robust Room Temperature Exchange Bias with Multiferroic BiFeO 3 Epitaxial Thin Films[END_REF] Please note, that in order to reduce leakage effects of the film between the electrode and the PFM tip, the BFO layer is not pure BFO but has been doped with 5 % of Manganese (Mn).

Regarding the growth conditions, the BFO/ LSMO //NGO (001) is grown by pulsed laser ablation under oxygen-rich growth conditions of 800 • C and 0.2 mbar for the LSMO film, and 650 • C and 6 • 10 -3 mbar for the BFMO film. The frequency of the laser is 2 Hz and 5 Hz, respectively. Annealing at 650 • C is carried after the growth of the bilayer for 30 minutes in an environment rich in oxygenates (200 mbar).

Device fabrication

For the fabrication of our device we employed a combination of laser lithography and electron ion beam machining for the smallest structures such as the microwave antennae.

As for the waveguide's dimensions it exhibits a width of 2 µm and length of 90 µm. A thin (d AlO 3 = 10 nm) alumina oxide (AlO 3 ) layer is sandwiched between the BFO and a titanium (d Ti = 20 nm) layer. The titanium layer enhances the adhesion of the gold for the addition of the microwave antennae for the propagative spin wave spectroscopy. These antennae (Au, d = 130 nm) exhibit a meander-like structure with a periodicity of 500 nm with several periods narrowing the wavevector band of excited spin waves.

METHODS BFO/LSMO film characterization

In order to characterize the surface quality of both components of our heterostructure, we performed reflection high-energy electron diffraction (RHEED) spectroscopy. The spectroscopic data is shown in Fig. S1 for LSMO on the NGO substrate only ((a)) and for the complete heterostructure with BFO ((b)). Both the RHEED patterns are well separated from each other, with an elongation in the vertical direction indicating a flat surfaces. Please note that, the uniaxial in-plane anistropy H ip u in the LSMO originates from our orthorhombic NGO substrate with a (001)/K(110) orientation which induces an anisotropic strain, where K denotes the anisotropy plane. For this substrate, the LSMO grows with an in-plane angle of 45 • to the [100] orthorombic direction in NGO in (001) pseudocubic configuration. 1 As a result, the BFO/LSMO waveguides exhibit a [110] uniaxial anistropy with respect to the LSMO crystallographic direction, yielding an angle φ 0 of 45 • between the direction of the easy axis and the waveguide axis (x-axis) which needs to be included in the dispersion relation. In addition, we include an out-of-plane anisotropy H oop u (along the z-axis) due to the magneto-elastic anisotropy induced by the slight lattice mismatch between LSMO and the NGO substrate leading to a compressive strain for our thin LSMO film. 2 The magnetic field is along the y-axis direction, in the DE configuration.

MEASURED AND DEDUCED PARAMETERS

Specifically, we measure the device's properties by probe-tip ferromagnetic resonance (FMR) with a vector network analyser (VNA) by recording the scattering matrix S for a two-port system. 3 Afterwards, for the discussion of spin wave propagation and modulation in LSMO, similar to Vlaminck et al., 4 we calculate the corresponding inductance matrix L

=      L 11 L 12 L 21 L 22      (1) = Z 0 δ      α+S 21 • S 12 2S 21 2S 12 β+S 21 • S 12      (2)
, where δ = ω 0 In order to determine the saturation magnetization, the Gilbert damping parameter and the contribution to the linewidth from inhomoegenous broading denoted by ∆H 0 , we performed ferromagnetic resonance (FMR) measurements on our pure films before any patterning. In Fig. S2 (a), we show the frequency dependence of the peak-to-peak linewidth ∆H pp shows an agreement between the experimental and analytically derived numbers. Please note, that due to the required in situ growth of LSMO and subsequently of BFO during the pulsed laser deposition process, only the FMR of the complete heterostructure could be measured. Thus, the obtained values are effective values of the heterostructure. The Curie Temperature has been also measured by resistivity measurements and was found to be T c ≈ 390 K.

Fig. 1 .

 1 Fig. 1. (a) Experimental setup of our BFO/LSMO heterostructure shown for the periodic configuration to obtain a reconfigurable magnonic crystal on the sub-micron scale. The spin wave is excited in Damon-Eshbach configuration and the ferroelectric domain state is written via piezo-force microscopy (PFM). The inductance matrix for a two-port system using vector network analysis is measured via probe-tip ferromagnetic resonance at room temperature (T = 300 K) (b) Analytical calculation of the antenna's response function with maximum at wavevector K ≈ 6 rad/µm. (c) Image of the radio frequency antennae with spacing of 500 nm and 2 µm distance between each arm and antenna, respectively. (d) Corresponding PFM imaging of the ferroelectric state after one domain writing cycle of a down, up and periodic configuration (P n,m , n: writing cycle, m: domain configuration)

Fig. 2 (

 2 Fig. 2 (a) and (b) shows the L 21 inductance signal of the spin wave for the two uniform up (red) and down polarized (blue) ferroelectric domain states in the BFO layer for the amplitude (in linear scale) and the imaginary part, respectively. The initial (virgin) down

Fig. 2 .

 2 Fig. 2. (a) Spin wave transmission measurement showing the inductance amplitude |L 21 | normalized by the mean value for the uniform domain states before writing the periodic state and subsequent reconfiguration. The insets show the corresponding domain state. (b) Imaginary part of the spin wave signal of L 21 . In (a) and (b) the P 0,↓ (blue) state is compared to the uniform P 0,↑ (red) state. Apart from a small frequency shift (≈ 15 MHz) and small change in the amplitude the spin wave propagation is reversible between the uniform domain states of the BFO. (c) Colormap of (|L 21 |)(H ext , ω) for the reconfigured P 1,↓ state after switching from the up state to the periodic configuration and back to the uniform down state. The red curve is the theoretical ω(H ext , k x = K B )(cf. equation. (1)). (d) Im(L 21 ) of the P 1,↓ state at 413 Oe. All measurements were conducted at T = 300 K between 3.3-3.9 GHz. The field independent background is attributed to residual radio frequency standing waves after a VNA short-open-load-through calibration.

  100 Oe, H oop u = -100 Oe, H offset = -100 Oe an offset taking additional contributions to the effective field and uncertainties in the determination of the saturation magnetization into account) together with the recorded dependence of ω n=1 (H) for the lowest mode number n = 1 and Im (L 21 ) at 413 Oe in Fig. 2 c-d for the P 1,↓ -BFO domain state. By employing ferromagnetic resonance, we find 4πM eff = (3323±23) G and a Gilbert damping α G as low as 6 • 10 -3 of our 13 nm thick LSMO layer (Supporting information section 2 (cf. Fig. S2)).

Fig. 3 .

 3 Fig. 3. (a) Amplitude representation of the reconfigurable magnonic crystal in the BFO/LSMO heterostructure. The transmitted spin wave amplitudes |L 21 | are shown in logarithmic scale for the uniform P 0,↑ (red), periodic P 0,↑↓ (green) and uniform P 1,↓ (blue) state (cf. inset) and compared.

Fig. 4 .

 4 Fig.4. Field dependence of the magnonic bandgap shown for magnetic field values from 405-458 Oe). The total amplitude varies due to the remaining, field independent background in the system. Hence, the amplitude of the frequency filtering changes. This variation is also present in the uniform domain state

  (a) Clear Rheed pattern of solely LSMO on NGO and (b) Rheed pattern after the addition of the BFO layer from our raw data as: L

( 1 -FERROMAGNETICS 2 .

 12 S 11 )(1-S 22 )-S 12 • S 21 ] with the resonance frequency ω 0 , α = (1+S 11 )(1-S 22 ), β = (1+S 22 )(1-S 11 ) and Z 0 is the impedance. Magnetic characterization of the unpatterned BFO/LSMO heterostructure measured with field-modulated ferromagnetic resonance at room temperature. (a) Frequency dependence of the linewidth for the determination of the Gilbert damping parameter (α G = (5.95±0.3) • 10 -3 ) and the intrinsic linewidth at zero applied external field of the specific film (∆H 0 = (3.3±0.9) Oe). (b) Ferromagnetic resonance frequency as a function of the applied static magnetic field and determination of the effective magnetization (4πM eff = (3323±23) G).

  including a linear fit to determine the damping. As mentioned in the main text, the effective Gilbert damping of our heterostructure made of a BFO/LSMO//NGO stack is determined to α G ≈ 6 • 10 -3 . In addition, the dispersion at T = 300 K, i.e. ω FMR (H res ) including a fit of the Kittel law for thin films is shown in Fig.S2 (b). The fit yields 4πM eff = (3323±23) G. A comparison of that value to the Probe-Tip FMR measurements of the frequency dispersion of the transmitted spin wave and using the analytical formula (cf. Equation (3) in main part)

ACKNOWLEDGEMENTS

This work has been supported from the European Union's Horizon 2020 research and innovation program within the FET-OPEN project CHIRON under grant agreement No. 801055. VH is grateful for support from the project OISO under grant agreement ANR-17-CE24-0026 from the Agènce National de la Recherche.

Propagative spin wave spectroscopy

The spin wave spectra are obtained via probe-tip ferromagnetic resonance. The resulting transmission and reflection scattering matrix S ij is recorded by means of network analysis employing a two-port (i and j) vector network analyzer setup.

SUPPORTING INFORMATION

We have provided supporting information for our work which is available online. The contents are:

• Material's characterizations including a RHEED diffraction pattern.

• Information on the measured scattering parameters and deduced inductance matrix