Laurent Gourvès
email: laurent.gourves@dauphine.fr

Jérôme Monnot

Fanny Pascual
email: fanny.pascual@lip6.fr

Daniel Vanderpooten

Bi-objective matchings with the triangle inequality

Keywords: Bi-objective optimization, approximation algorithm, matching

HAL is

Introduction

For many optimization problems, there is a need to deal with several aspects of a solution, at the same time. In multi-objective optimization, a problem has a set of feasible solutions S and k > 1 objective functions f i : S → IR, i ∈ {1, . . . , k}. For each f i , it is also specified if the function should be minimized or maximized. Suppose every f i should be maximized and let opt i denote max s∈S f i (s), for all i ∈ {1, . . . , k}. For a given instance of a multi-objective problem, the ideal point (opt 1 , . . . , opt k ) is the image of a (not necessarily feasible) solution reaching optimality on all objective functions, at the same time.

In practice, it is unlikely that the ideal point is the image of a feasible solution. We thus need to resort to approximation. A feasible solution s is said to be a ρ-approximation of the ideal point (or, in short, a ρ-approximate solution), for ρ ∈ (0, 1], if f i (s) ≥ ρ opt i , i ∈ {1, . . . , k}.

The existence of a feasible ρ-approximate solution is not guaranteed for every possible ratio ρ ∈ (0, 1]. Therefore, it is challenging to identify, for a given problem, the best approximation factor under which the existence of a feasible ρ-approximate solution is always guaranteed. It is also relevant to have a constructive approach: what is the largest ratio ρ ∈ (0, 1] such that a feasible ρ-approximation of the ideal point can be computed in polynomial time?

In this article, we focus on matching problems, which are central in combinatorial optimization problems. In particular, we deal with a bi-objective matching problem. The set of feasible solutions is every subset of pairwise non-adjacent edges. Each edge e has a non-negative weight w(e) and a non-negative length ℓ(e), where both objectives satisfy the triangle inequality.

After a short review of the literature related to our problem (section 2), we give in section 3 a formal definition of the model, together with a motivation for restricting ourselves to two objectives satisfying the triangle inequality. Section 4 is devoted to the computational complexity of the bi-objective matching problem on complete graphs. In section 5 we give an approximation algorithm. The result is a general method which, combined with a technical lemma, provides an approximation algorithm with performance guarantee 1/3.

In section 6 we reuse the general method but we substitute the technical lemma with computational results, to give an insight about possible improvements of the previous performance guarantee of 1/3. Some upper bounds are given in section 7. We end this article with concluding remarks.

Related work

This article falls in the field of multi-objective combinatorial optimization (see [START_REF] Ehrgott | Multicriteria Optimization[END_REF] Chapter 8 in particular). A prominent solution concept in multi-objective optimization is the Pareto set P (a.k.a. Pareto curve [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF]). For every feasible solution s, P contains a feasible solution s * such that the image of s * weakly dominates1 the image of s. However, the notion of Pareto set is sometimes problematic since, for some instances, P must contain a number of solutions that is not polynomial in the input size. Therefore, no algorithm can compute P in polynomial time. Another issue comes from computational complexity because the multi-objective version of many polynomial single-objective problems is NP-hard (this is for example the case for the minimum cost spanning tree problem [START_REF] Aggarwal | Minimal spanning tree subject to a side constraint[END_REF] and for the minimum assignment problem [START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF]). We can remedy to these difficulties with the help of approximation. If the maximization is retained for every f i , then we say that s ′ is a ρ-approximation of s, with ρ ∈ (0, 1], when f i (s ′ ) ≥ ρ f i (s) for all i ∈ {1, . . . , k}. An approximate Pareto set P ρ contains a ρ-approximation for every feasible solution.

Interestingly, the notion of approximate Pareto set helps to circumvent the two aforementioned difficulties. For a constant number of objectives, Papadimitriou and Yannakakis [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF] have shown that for every ρ ∈ (0, 1], there always exists an approximate Pareto set P ρ such that |P ρ | is polynomial in both the input size and (1ρ) -1 . Moreover, for the multi-objective version of some paradigmatic problems in combinatorial optimization (e.g. minimum cost spanning tree), P ρ can be computed in polynomial time.

Usually, no explicit upper bound on |P ρ | is given. However, we can be interested in fixing the size of an approximate Pareto set and try to provide the best performance guarantee. Because P ρ may contain an overwhelming number of solutions, limiting its size is useful (in a decision process for example). When the size of an approximate Pareto set is upper bounded, we cannot guarantee that every approximation factor can be achieved. Intuitively, the larger the size of the approximation set, the more accurate it can be. In fact, it is not difficult to see that a ρ-approximation of the ideal point is just a ρ-approximate Pareto set of size 1. Therefore, we interchangeably talk about ρ-approximation of the ideal point and ρ-approximate Pareto set (with a single solution).

Approximation of the ideal point has been done for the multi-objective versions of many well-studied optimization problems. In [START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF][START_REF] Rasala | Existence theorems, lower bounds and algorithms for scheduling to meet two objectives[END_REF][START_REF] Angel | A note on scheduling to meet two min-sum objectives[END_REF], the authors study scheduling problems with two classical objective functions (e.g. makespan, sum of completion times, etc). They show the existence (and computation) of schedules that are approximately optimal for the objective functions, at the same time.

A 1/3-approximation of the ideal point exists for the bi-objective maximum spanning tree problem [START_REF] Escoffier | Fair solutions for some multiagent optimization problems[END_REF]. In the same article, the authors give an approximation of two other multi-objective problems: max sat and cut-complement. In [START_REF] Gourvès | Approximate tradeoffs on weighted labeled matroids[END_REF], the 1/3-approximation presented in [START_REF] Escoffier | Fair solutions for some multiagent optimization problems[END_REF] has been generalized to a bi-objective problem on simple matroids.

A bi-objective version of max cut is studied in [START_REF] Angel | Approximation algorithms for the bi-criteria weighted max-cut problem[END_REF]: there always exists a 1/2-approximation of the ideal point whereas a 0.439-approximation can be computed in polynomial time. These results were recently generalized to the maximization of two sub-modular, symmetric and nonnegative functions2 [START_REF] Bazgan | Approximation with a fixed number of solutions of some multiobjective maximization problems[END_REF].

The ideal point approach has also been used on the bi-objective version of the traveling salesman problem [START_REF] Manthey | Deterministic algorithms for multi-criteria max-tsp[END_REF][START_REF] Manthey | On approximating multicriteria tsp[END_REF][START_REF] Bazgan | Single approximation for the biobjective max tsp[END_REF]. For the maximization version of the problem, the best ratios for the approximation of the ideal point come from [START_REF] Bazgan | Single approximation for the biobjective max tsp[END_REF]: 5/12ε if the two objectives satisfy the triangle inequality, 3/8ε if only one objective satisfies the triangle inequality, and 1+2 √ 2 14

ε without the triangle inequality.

Actually, the bi-objective maximum matching (where the triangle inequality is satisfied) has already been studied in [START_REF] Bazgan | Approximation with a fixed number of solutions of some multiobjective maximization problems[END_REF]. This article provides a general approximation method for a class of biobjective problems. An application of the general method to the bi-objective maximum matching yields a 1/6-approximation of the ideal point. In addition, the problem is shown NPhard when the instance is a collection of complete graphs in which the triangle inequality is satisfied [START_REF] Bazgan | Approximation with a fixed number of solutions of some multiobjective maximization problems[END_REF]. The case of a complete graph was left open. In the present article, we improve on the approximation result since a 1/3-approximation is provided (Corollary 1). Moreover the complexity of the problem (previously left open) is shown on complete graphs (Theorem 1).

It is noteworthy that bi-objective versions of the assignment problem (i.e. matchings in bipartite graphs) have been studied previously. Exact methods for the minimization version can be found in [START_REF] Ulungu | The two-phases method: An efficient procedure to solve bi-objective combinatorial optimization problems[END_REF] and [START_REF] Przybylski | Two phase algorithms for the bi-objective assignment problem[END_REF].

To conclude, let us mention a popular approach in multi-objective combinatorial optimization. It consists in turning all but one objective functions in budget constraints. Given a vector (B 2 , B 3 , . . . , B k ) ∈ R k-1 , a solution s ∈ S is said to be optimal if it optimizes f 1 (s) under the budget constraints f j (s) ≥ B j (or f j (s) ≤ B j ), j = 2, . . . , k. To mention a few examples, this method is used for some multi-objective network-design problems [START_REF] Ravi | Many birds with one stone: multi-objective approximation algorithms[END_REF] or shortest path [START_REF] Hassin | Approximation schemes for the restricted shortest path problem[END_REF]. The matching problem, and generalizations to matroids and independence systems are studied in [START_REF] Berger | Budgeted matching and budgeted matroid intersection via the gasoline puzzle[END_REF][START_REF] Grandoni | New approaches to multi-objective optimization[END_REF].

Model

Let G = (V, E) be a complete graph. Each edge e ∈ E has a non-negative weight w(e) and a nonnegative length ℓ(e). We assume that both w and ℓ satisfy the triangle inequality 3 . The weight and length of a set of edges M ⊂ E are defined as w(M ) = e∈M w(e) and ℓ(M ) = e∈M ℓ(e), respectively. Both w(M ) and ℓ(M ) should be maximized. Since it is unlikely that a matching reaches optimality for both objectives, our aim is to find a single matching which constitutes the best approximation of the ideal point.

Note that we have to restrict to a bi-objective matching where the two edge values fulfill the triangle inequality, otherwise there is not always a single matching which constitutes an ε-approximate Pareto set with any given positive ratio ε:

-If the triangle inequality does not hold then at least 3 solutions are required to constitute an approximate Pareto set for any given positive ratio ε. See the instance on Figure 1 which contains 3 matchings with values (1 + ε, 0), (0, 1 + ε) and (ε, ε), respectively. No solution is an ε-approximation of any other one, so the three solutions are required to provide an ε-approximate Pareto set. -If only one objective fulfills the triangle inequality, then at least 2 solutions are required to constitute an ε-approximate Pareto set with any given positive ratio ε: consider a triangle with values (0, 1), (1, 0), and (1, 0) (see Figure 2). -If we consider 3 objectives which all fulfill the triangle inequality, then at least 3 solutions are required to constitute an ε-approximate Pareto set with any given positive ratio ε: consider a triangle with values (0, 1, 1), (1, 0, 1), and (1, 1, 0) (see Figure 3). (0,1) (0,1,1)

( 0 , 1 + ε ) (ε, ε) ( 1 + ε , 0 )
(1,0) (1,0)
(1,1,0) (1,0,1)

Complexity

We first investigate the complexity of the decision version of the biobjective maximum matching problem. This corresponds to the problem RestrictedMaxMatching defined as follows. Let (x 1 , x 2 ) ∈ N 2 , does there exist a matching M such that the weight of M is at least x 1 , and the length of M is at least x 2 ? Theorem 1. RestrictedMaxMatching is NP-complete, even in a graph in which both objectives fulfill the triangle inequality.

Proof. We will show that Partition polynomially reduces to RestrictedMaxMatching. We recall the Partition problem: an instance of Partition is a set E = {e 1 , . . . , e n } of n elements, and a cost function c : E → N such that e∈E c(e) = 2B and 0 < c(e) < B. The question is: is there a subset

E ′ ⊆ E such that e∈E ′ c(e) = e∈E\E ′ c(e) = B?
The corresponding instance of RestrictedMaxMatching is as follows: we fix x 1 = x 2 = (4n+1)B, and we consider a complete graph G composed of 4n vertices : n vertices {a 1 , . . . , a n }, n vertices {b 1 , . . . , b n }, n vertices {c 1 , . . . , c n }, and n vertices {d 1 , . . . , d n }. For each i ∈ {1, . . . , n}, the cost of edge (a i , b i ) is equal to (2B + c(e i ), 2B), and the cost of edge (b i , c i ) is equal to (2B, 2B + c(e i )). All the other edges have a cost equal to (2B, 2B). This graph satisfies the triangle inequality since on each objective the cost of each edge is at least 2B and at most 3B.

Let us show that there is a feasible solution of RestrictedMaxMatching on this instance if and only if (E, c) is a yes instance of Partition. If (E, c) is a yes instance of Partition (i.e. there is a subset

E ′ ⊆ E such that e∈E ′ c(e) = e∈E\E ′ c(e) = B
), then there is a solution M of the corresponding instance of the RestrictedMaxMatching problem: for each i ∈ {1, . . . , n}, if e i ∈ E \ E ′ then the edges (a i , b i ) and (c i , d i ) are in M , and else, if

e i ∈ E ′ then the edges (b i , c i ) and (a i , d i ) are in M (no other edge is in M ). Since e∈E ′ c(e) = e∈E\E ′ c(e) = B, the weight of M is then x 1 = 2n(2B) + B = (4n + 1)B, and the length of M is x 2 = 2n(2B) + B = (4n + 1)B.
Conversely, let us now show that if there is a feasible solution M of RestrictedMax-Matching then (E, c) is a yes instance of Partition. We assume that M is a perfect matching, and is thus made of 2n edges; it is possible because the graph is complete, it contains an even number of nodes and any matching can be completed to become perfect without decreasing its value because all the edges have positive values for each objective. The weight of M is |M |2B + (ai,bi)∈M c(e i ). Since the weight of M is at least x 1 = (4n + 1)B, we have

(ai,bi)∈M c(e i ) ≥ B. Let E 1 = {e i ∈ E|(a i , b i ) ∈ M }. We have: ei∈E1 c(e i ) ≥ B. Likewise, the length of M is |M |2B + (bi,ci)∈M c(e i ). Since the length of M is at least x 2 = (4n+1)B, we have (bi,ci)∈M c(e i ) ≥ B. Let E 2 = {e i ∈ E|(b i , c i ) ∈ M }. We have: ei∈E2 c(e i ) ≥ B. Note that for each i ∈ {1, . . . , n}, at most one edge in {(a i , b i ), (b i , c i )} belongs to M , since M is a matching. Thus E 1 ∩E 2 = ∅. Since e∈E c(e) = 2B, we have ei∈E1 c(e i ) = ei∈E2 c(e i ) = B,
which is a solution to the (E, c) instance of Partition.

⊓ ⊔

Approximation

This section provides a general method, described in Theorem 2, which, combined with some results for small graphs, described in Lemma 2, provides an approximation algorithm with performance guarantee 1/3.

Theorem 2. Let k be a fixed integer satisfying k ≥ 2. Let ρ be a real satisfying 0 < ρ ≤ 1.

Suppose that for every q ≤ 2k and every complete graph K q of q nodes, there always exists a matching F q satisfying w(F q ) ≥ ρw(M q w ) and ℓ(F q ) ≥ ρℓ(M q ℓ ), where M q w (resp. M q ℓ ) denote a maximum weight (resp. length) matching of K q . Then, there exists a single matching which constitutes a ρ k-1 k -approximate Pareto set for the biobjective Maximum Matching problem. Proof. Let k ≥ 2. Build two maximum matchings of K n with respect to w and ℓ respectively. Denote them by M w and M ℓ . We can suppose wlog that both M w and M ℓ are maximal with respect to the inclusion and then

|M w | = |M ℓ |. The subgraph G ′ = (V, M w ∪ M ℓ
) is made of p connected components denoted by C i for i ∈ {1, . . . , p}. One C i may be one isolated node but every other connected component contains at least one edge. Each C i is a cycle of even size or a path. Note that there is at most one path with at least two edges in M w ∪M ℓ (because the graph is complete and we can assume that M w are M ℓ are of maximum size) and this path contains an even number of edges. Each path containing a single edge must be in M w ∩ M ℓ (because M w and M ℓ are of maximum size: at least one endpoint of each edge should be adjacent to one edge of M ℓ -and this is the same thing for M w ).

Consider a connected component C i with at least two edges and denote by n i its number of vertices. Thus

C i contains n i edges if C i is a cycle and n i -1 if C i is a path.
-If n i ≤ 2k, then by hypothesis, there exists in K ni (the subgraph induced by the vertices of C i ) a matching F ni such that:

w(F ni ) ≥ ρw(M ni w ) ≥ ρw(C i ∩ M w ) and ℓ(F ni ) ≥ ρℓ(M ni ℓ ) ≥ ρℓ(C i ∩ M ℓ ) (1)
since M ni w (resp. M ni ℓ ) is an optimum weight (resp. length) matching of K ni . -Now, suppose that n i ≥ 2k +1. The edges of C i ∩M w are labeled with numbers in {1, . . . , k} as follows. If C i is a cycle, then it is explored from an arbitrarily chosen vertex that we denote v 0 . If C i is a path, then it is explored from an endpoint denoted by v 0 . During the exploration, the first edge of C i ∩ M w is labeled 1, the second edge of C i ∩ M w is labeled 2, . . . , the k + 1 edge of C i ∩ M w is labeled 1, and so on. Formally, the j-th edge of C i ∩ M w is labeled j -⌊ j-1 k ⌋k. See Figure 4 for an illustration for a cycle of length 12 and k = 3. For j ∈ {1, . . . , k}, L j (C i ∩ M w ) are the edges of C i ∩ M w with label j. Assume that j w is the index such that the matching of label j w is the lightest, i.e., w(L jw (C i ∩ M w )) = min j=1,...,k {w(L j (C i ∩ M w ))}. Thus, we deduce that

w(C i ∩ M w \ L jw (C i ∩ M w )) ≥ k -1 k w(C i ∩ M w ) (2)
Moreover, since L jw (C i ∩ M ℓ ) = ∅ (the edges which have been labelled belong to C i ∩ M w and in a connected component C i with at least two edges

(C i ∩ M w ) ∩ (C i ∩ M ℓ ) = ∅), we have: ℓ(C i ∩ M ℓ \ L jw (C i ∩ M ℓ )) = ℓ(C i ∩ M ℓ ) (3)
Observe that if C i is a cycle, then all the paths in C i \ L jw (C i ∩ M w ) except one have a length 2(k -1) + 1 = 2k -1 and the last path P i 0 (containing v 0 ) has a length at most equal to 2(2(k -1)) + 1 = 4k -3. See Figure 5 for an illustration of the worst case with n i = 10 and k = 3; here edge e corresponds to matching L 3 (C i ∩ M w ). Let us focus on the path P i 0 when C i is a cycle and assume that the length of this path is at least 2k (actually, it is at least 2k + 1 because P i 0 is of odd length). This path contains at most 4k -3 edges and by construction the two end edges are in M ℓ . We label the edges of P i 0 ∩ M ℓ with numbers in {1, . . . , k} as previously starting from an end edge with label 1. See Figure 6 for an illustration with n i = 10 and k = 3; Using a similar approach, L j (P i 0 ∩ M ℓ ) denotes the set of edges of P i 0 ∩ M ℓ with label j for j ∈ {1, . . . , k}. Assume that j ℓ is the index such that the matching of label j ℓ is the lightest, i.e., ℓ(L j ℓ (P i 0 ∩ M ℓ )) = min j=1,...,k {ℓ(L j (P i 0 ∩ M ℓ ))}. Thus, we obtain:

ℓ(P i 0 ∩ M ℓ \ L j ℓ (P i 0 ∩ M ℓ )) ≥ k -1 k ℓ(P i 0 ∩ M ℓ ) (4)
w(P i 0 ∩ M w \ L j ℓ (P i 0 ∩ M ℓ )) = w(P i 0 ∩ M w ) (5)
If C i is a path, then all the paths in C i \ L jw (C i ∩ M w ) contain at most 2k -1 edges. In this case, we set P i 0 = ∅.

In conclusion, we have the following property: each path P of C i \ (L jw (C i ∩ M w ) ∪ L j ℓ (P i 0 ∩ M ℓ )) is contained in a complete subgraph K q ′ where q ′ ≤ 2k (in particular, it is the case of P i 0 when C i is a cycle with length at most 2k -1). By hypothesis there is a matching F q ′ of K q ′ such that:

w(F q ′ ) ≥ ρw(M q ′ w ) ≥ ρw(P ∩ M w ) and ℓ(F q ′ ) ≥ ρℓ(M q ′ ℓ ) ≥ ρℓ(P ∩ M ℓ ) (6)
since M q ′ w (resp. M q ′ ℓ ) is optimal for w (resp. ℓ). The union of these matchings F q ′ , denoted by F i , is a matching of G induced by the vertices of C i . Moreover, using Inequalities (2), (3), ( 4), [START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF], and (6) we get:

w(F i ) ≥ ρw(C i \ (L jw (C i ∩ M w ) ∪ L j ℓ (P i 0 ∩ M ℓ )) ∩ M w ) ≥ ρ k -1 k w(C i ∩ M w ) (7) ℓ(F i ) ≥ ρℓ(C i \ (L jw (C i ∩ M w ) ∪ L j ℓ (P i 0 ∩ M ℓ )) ∩ M w ) ≥ ρ k -1 k ℓ(C i ∩ M ℓ ) (8)
By summing up previous inequalities (1), or [START_REF] Angel | A note on scheduling to meet two min-sum objectives[END_REF], or (8) according to the case and by adding the connected components made of a single edge, we get a matching F such that w(F

) ≥ ρ k-1 k i w(C i ∩ M w ) = ρ k-1 k w(M w ) and ℓ(F ) ≥ ρ k-1 k i ℓ(C i ∩ M ℓ ) = ρ k-1 k ℓ(M ℓ ). ⊓ ⊔
We show now the existence of 1 2 -approximations for graphs of small size. For this purpose, the following simple result will be useful.

Lemma 1. Let G = (V, E) be a weighted complete graph with |V | ≥ 4 and weight function h satisfying the triangle inequality. For any four vertices i, j, k, l, at least one of the following two inequalities holds :

h(i, k) + h(j, l) ≥ 1 2 (h(i, j) + h(k, l)) h(i, l) + h(j, k) ≥ 1 2 (h(i, j) + h(k, l))
Proof. Using the triangle inequality, we have h(i, k) + h(j, k) ≥ h(i, j), h(i, l) + h(j, l) ≥ h(i, j), h(i, k) + h(i, l) ≥ h(k, l), and h(j, k) + h(j, l) ≥ h(k, l). Summing up these inequalities we get h(i, k) + h(j, l) + h(i, l) + h(j, k) ≥ h(i, j) + h(k, l). The result follows. ⊓ ⊔ Lemma 2. For complete graphs with at most 6 vertices, with values on edges w and ℓ satisfying the triangle inequality, there exists a single matching M such that w(M ) ≥ 1 2 w(M w ) and ℓ(M ) ≥ 1 2 ℓ(M ℓ ), where M w and M ℓ are optimal matchings for the weight and the length, respectively.

Proof. Let n ∈ {3, 4, 5} since the case n = 2 is trivial. Build two optimal matchings on K n with respect to w and ℓ denoted respectively by M w and M ℓ . We assume wlog that the subgraph G ′ induced by M w ∪M ℓ is connected since otherwise we deal with connected components separately. Thus, the different cases are depicted in Figure 7. Let us show now that for each of these cases, there exists a matching M satisfying w(M ) ≥ 1 2 w(M w ) and ℓ(M ) ≥ 1 2 ℓ(M ℓ ).

-Case n = 3. G ′ is a path. We have -Case n = 4. G ′ is a cycle. We have M w = {(1, 2), (3, 4)} and M ℓ = {(1, 4), (2, 3)}. Using Lemma 1 for w, at least one of the following two inequalities holds :

M w = {(1, 2)} and M ℓ = {(2, 3)}. If ℓ(1, 2) ≥ ℓ(2, 3)/2 then M = {(1, 2)}. If w(2, 3) ≥ w(1, 2)/2 then M = {(1, 2)}. Otherwise, M = {(1, 3)} since,
w(1, 4) + w(2, 3) ≥ 1 2 (w(1, 2) + w(3, 4)) (9) w(1, 3) + w(2, 4) ≥ 1 2 (w(1, 2) + w(3, 4)) (10)
Using Lemma 1 for ℓ, at least one of the following two inequalities holds :

ℓ(1, 2) + ℓ(3, 4) ≥ 1 2 (ℓ(1, 4) + ℓ(2, 3)) (11) ℓ(1, 3) + ℓ(2, 4) ≥ 1 2 (ℓ(1, 4) + ℓ(2, 3)) (12)
If ( 9) holds then M = {(1, 4), (2, 3)}. If [START_REF] Bazgan | Approximation with a fixed number of solutions of some multiobjective maximization problems[END_REF] holds then M = {(1, 2), (3, 4)}. Finally, if neither ( 9) nor ( 11) hold then both ( 10) and ( 12) must be true and we have M = {(1, 3), (2, 4)}.

-Case n = 5. G ′ is a path. We have

M w = {(1, 2), (3, 4)} and M ℓ = {(2, 3), (4, 5)}. If w(1, 2) + w(4, 5) ≥ 1 2 (w(1, 2) + w(3, 4)) and ℓ(1, 2) + ℓ(4, 5) ≥ 1 2 (ℓ(2, 3) + ℓ(4, 5
)), we have M = {(1, 2), (4, 5)}. Otherwise, we have:

• either w(1, 2)+w(4, 5) < 1 2 (w(1, 2)+w [START_REF] Aggarwal | Minimal spanning tree subject to a side constraint[END_REF][START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF]). This implies w(3, 4) ≥ 1 2 (w(1, 2)+w(3, 4)) and w(3, 5) ≥ 1 2 (w(1, 2) + w(3, 4)), since otherwise the converse inequalities, summed up with the initial inequality, would lead respectively to w(1, 2) + w(4, 5) + w(3, 4) < w(1, 2) + w(3, 4) and w(1, 2) + w(4, 5) + w(3, 5) < w(1, 2) + w(3, 4) which are clearly wrong, using the triangle inequality in the second case. Moreover, using Lemma 1 for ℓ, we have either ℓ(2, 5)

+ ℓ(3, 4) ≥ 1 2 (ℓ(2, 3) + ℓ(4, 5)), in which case M = {(2, 5), (3, 4)}, or ℓ(2, 4) + ℓ(3, 5) ≥ 1 2 (ℓ(2, 3) + ℓ(4, 5)), in which case M = {(2, 4), (3, 5)}. • or ℓ(1, 2) + ℓ(4, 5) < 1 2 (ℓ(2, 3) + ℓ(4, 5)). This implies ℓ(2, 3) ≥ 1 2 (ℓ(2, 3) + ℓ(4, 5)) and ℓ(1, 3) ≥ 1 2 (ℓ(2, 3) + ℓ(4, 5
)), since otherwise the converse inequalities, summed up with the initial inequality, would lead respectively to ℓ(1, 2)+ℓ(4, 5)+ℓ(2, 3) < ℓ(2, 3)+ℓ(4, 5) and ℓ(1, 2) + ℓ(4, 5) + ℓ(1, 3) < ℓ(2, 3) + ℓ(4, 5), which are clearly wrong, using the triangle inequality in the second case. Moreover, using Lemma 1 for w, we have either w(1, 4) + w(2, 3)

≥ 1 2 (w(1, 2) + w(3, 4), in which case M = {(1, 4), (2, 3)}, or w(1
, 3) + w(2, 4) ≥ 1 2 (w(1, 2) + w(3, 4), in which case M = {(1, 3), (2, 4)}.
-Case n = 6. G ′ is a cycle. We have

M w = {(1, 2
}. We first show that there exists i ∈ {1, 2, 3} such that ℓ(M i ) ≥ 1 2 ℓ(M ℓ ). By contradiction, assuming that ℓ(M i ) < 1 2 ℓ(M ℓ ) for all i ∈ {1, 2, 3}, we multiply by 2 each inequality and sum them up, which gives 2ℓ [START_REF] Ehrgott | Multicriteria Optimization[END_REF][START_REF] Rasala | Existence theorems, lower bounds and algorithms for scheduling to meet two objectives[END_REF] [START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF][START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF]. Applying the triangle inequality to each term between brackets gives values respectively larger than or equal to ℓ(1, 6), ℓ(2, 3), and ℓ(4, 5), contradicting the inequality. In the following, we assume wlog, after possibly renumbering the vertices, that ℓ(M 1 ) ≥

+ (ℓ(1, 2) + ℓ(2, 5) + ℓ(5, 6)) + 2ℓ(2, 3) + (ℓ(1, 2) + ℓ(1, 4) + ℓ(3, 4)) + 2ℓ(4, 5) + (ℓ(3, 4) + ℓ(3, 6) + ℓ(5, 6)) + ℓ(1, 4) + ℓ(2, 5) + ℓ(3, 6) < 3ℓ(1, 6) + 3ℓ(2, 3) + 3ℓ
1 2 ℓ(M ℓ ). If w(M 1 ) ≥ 1 2 w(M w ) then M = M 1 . If w(M 1 ) < 1 2 w(M w ) then, by Lemma 1, we get w(M 4 ) ≥ 1 2 w(M w ) where M 4 = {(1, 2), (3, 5), (4, 6)}. If ℓ(M 4 ) ≥ 1 2 ℓ(M ℓ ) then M = M 4 . Thus, we assume ℓ(M 4 ) < 1 2 ℓ(M ℓ ).
Consider now the two additional matchings M 5 = {(1, 5), (2, 6), (3, 4)} and M 6 = {(1, 3), (2, 4), (5, 6)}. We show that the four following conditions hold.

there exists i ∈ {4, 5, 6} such that ℓ(M i ) ≥ 1 2 ℓ(M ℓ ) (13)
there exists i ∈ {2, 3, 4} such that ℓ(M i ) ≥ 1 2 ℓ(M ℓ ) (14)
there exists i ∈ {2, 4, 6} such that ℓ(

M i ) ≥ 1 2 ℓ(M ℓ ) (15)
there exists i ∈ {3, 4, 5} such that ℓ(

M i ) ≥ 1 2 ℓ(M ℓ ) (16)
Regarding [START_REF] Manthey | On approximating multicriteria tsp[END_REF], the real condition is: there exists i ∈ {4, 5, 6, w} such that ℓ(M i ) ≥ 1 2 ℓ(M ℓ ). Obviously, if i = w then M = M w , and therefore we can assume [START_REF] Manthey | On approximating multicriteria tsp[END_REF], provided that the real condition is true. Assume by contradiction that ℓ(M i ) < 1 2 ℓ(M ℓ ) for all i ∈ {4, 5, 6, w}. Then, adding up these four inequalities, we obtain (ℓ(1, 2) + ℓ(2, 6)) + (ℓ(1, 5) + ℓ(5, 6)) + (ℓ(1, 2) + ℓ(1, 3))+(ℓ(2, 4)+ℓ [START_REF] Aggarwal | Minimal spanning tree subject to a side constraint[END_REF][START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF])+(ℓ(3, 4)+ℓ [START_REF] Aggarwal | Minimal spanning tree subject to a side constraint[END_REF][START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF])+(ℓ(4, 6)+ℓ [START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF][START_REF] Rasala | Existence theorems, lower bounds and algorithms for scheduling to meet two objectives[END_REF]) < 2ℓ(1, 6)+2ℓ(2, 3)+2ℓ [START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF][START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF]. Applying the triangle inequality to each term between brackets leads to contradict the inequality. Regarding [START_REF] Bazgan | Single approximation for the biobjective max tsp[END_REF], we proceed again by contradiction, assuming that ℓ(M i ) < 1 2 ℓ(M ℓ ) for all i ∈ {2, 3, 4}. Multiplying by two the inequality for i = 4 and adding the two other inequalities, we obtain ℓ(1, 6) + (ℓ(1, 4) + ℓ(4, 6)) + ℓ(2, 3) + (ℓ(2, 5) + ℓ(3, 5)) + (ℓ(3, 4) + ℓ(3, 5)) + (ℓ(4, 6) + ℓ(5, 6)) + 2ℓ(1, 2) < 2ℓ(1, 6) + 2ℓ(2, 3) + 2ℓ [START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF][START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF]. Applying the triangle inequality to each term between brackets leads to contradict the inequality. Conditions ( 15) and ( 16) are established similarly (each time multiplying by 2 the inequality for i = 4).

Since ℓ(M 4 ) < 1 2 ℓ(M ℓ ), by [START_REF] Manthey | On approximating multicriteria tsp[END_REF] we have either ℓ(M 5 ) ≥ 1 2 ℓ(M ℓ ) or ℓ(M 6 ) ≥ 1 2 ℓ(M ℓ ). We conclude the proof by examining these two cases:

• Consider the case ℓ(M 5 ) ≥ 1 2 ℓ(M ℓ ). If w(M 5 ) ≥ 1 2 w(M w ) then M = M 5 . Otherwise, if w(M 5 ) < 1 2 w(M w ) then by Lemma 1, we get w(M 2 ) ≥ 1 2 w(M w ). If ℓ(M 2 ) ≥ 1 2 ℓ(M ℓ ) then M = M 2 . If ℓ(M 2 ) < 1
2 ℓ(M ℓ ) then, by ( 14) and ( 15) we get respectively ℓ(M 3 ) ≥

1 2 ℓ(M ℓ ) and ℓ(M 6 ) ≥ 1 2 ℓ(M ℓ ). If w(M 3 ) ≥ 1 2 w(M w ) then M = M 3 . If w(M 3 ) < 1
2 w(M w ) then, by Lemma 1, we have w(M 6 ) ≥ 1 2 w(M w ) and thus M = M 6 . • Consider finally the case ℓ(M 6 ) ≥ 1 2 ℓ(M ℓ ). If w(M 6 ) ≥ 1 2 w(M w ) then M = M 6 . Otherwise, if w(M 6 ) < 1 2 w(M w ) then by Lemma 1, we get w(M 3 ) ≥ 1 2 w(M w ). If ℓ(M 3 ) ≥ 1 2 ℓ(M ℓ ) then M = M 3 . If ℓ(M 3 ) < 1 2 ℓ(M ℓ ) then, by ( 14) and ( 16) we get respectively ℓ(M 2 ) ≥

A computer-aided approach

Lemma 2 establishes the existence of 1/2-approximate matchings for complete graphs with at most n = 6 nodes. It seems difficult to establish a general result, and our proofs become more and more tedious as n grows. In order to test the existence of a single 1/2-approximate matching for larger values of n, we propose a mixed 0-1 linear programming formulation. Actually, this formulation provides, for any value of n, the best possible approximation ratio ρ * .

To introduce this formulation, we need to consider, for a given n:

the set M n of all maximal (perfect or nearly perfect) matchings which can be defined on a complete graph of size n. Note that the number of such maximal matchings is n!! when n is odd and (n -1)!! when n is even4 . -the set I ∆ n of all possible instances corresponding to a complete valued graph G = (V, E) with |V | = n where each edge e ∈ E has two values w(e) and ℓ(e) which satisfy the triangle inequality. We assume wlog that 0 ≤ w(e) ≤ 1 and 0 ≤ ℓ(e) ≤ 1 for each edge e ∈ E. Moreover, we impose that M w and M ℓ , which are optimal matchings for w and ℓ respectively, are such that w(M w ) = ℓ(M ℓ ) = 1. If these conditions are not met, then compute w(M w ) and ℓ(M ℓ ), and for each edge e, replace w(e) by w(e)/w(M w ) and ℓ(e) by ℓ(e)/ℓ(M ℓ ). This transformation has no incidence on the approximation guarantee because it is a ratio. We also assume the worst case situation where the subgraph induced by M w ∪ M ℓ is connected since otherwise we deal with connected components separately. Therefore, M ∪M ℓ is either a cycle or a path which alternates edges of M w and M ℓ and covers all vertices (see Figure 7 for the cases n = 3, 4, 5, 6).

Our formulation aims at identifying an instance I ∈ I ∆ n minimizing a variable ρ such that for all matchings M ∈ M n we have w(M ) ≤ ρ or ℓ(M ) ≤ ρ. Let ρ * n be the optimal value for variable ρ. Therefore, for any ρ < ρ (with at least one equality).

The mixed 0-1 linear program corresponding to the previous formulation is given by:

(P n )                                        min ρ s.t. w ij ≤ w ik + w kj and ℓ ij ≤ ℓ ik + ℓ kj , ∀(i, j, k) ∈ V 3 (17)
e∈Mw w e = 1 and e∈M ℓ ℓ e = 1 [START_REF] Hassin | Approximation schemes for the restricted shortest path problem[END_REF] e∈M w e ≤ 1 and

e∈M ℓ e ≤ 1 ∀M ∈ M n ( 19
) e∈M w e ≤ ρ + z M , ∀M ∈ M n ( 20
) e∈M ℓ e ≤ ρ + 1 -z M , ∀M ∈ M n (21) 0 ≤ w ij ≤ 1 and 0 ≤ ℓ ij ≤ 1, ∀(i, j) ∈ E, i < j (22) z M ∈ {0, 1}, ∀M ∈ M n (23) 0 ≤ ρ ≤ 1 (24)
where continuous decision variables w ij and ℓ ij , defined in (22), represent an instance satisfying the triangle inequality through constraint [START_REF] Ravi | Many birds with one stone: multi-objective approximation algorithms[END_REF], and 0-1 variables, defined in (23), are used to model the disjunctive constraints w(M ) ≤ ρ or ℓ(M ) ≤ ρ imposed by ( 20) and (21), whereas constraints ( 18) -( 19) define matchings M w and M ℓ as optimal matchings.

Owing to the quickly growing cardinality of M n , we were able to solve (P n ) for n up to 10, using CPLEX 12.6. The results are reported in Table 1. It is worth noting that these values coincide, up to the last (15th) computed decimal digit, with simple fractions larger than or equal to 1/2 (see last row of Table 1 1. The ratios obtained with the computational approach for 3 ≤ n ≤ 10.

These computational results cannot be considered as totally valid owing to the limited precision of the computations. They suggest, however, the existence of 1/2-approximate matchings for complete graphs satisfying the triangle inequality, up to 10 nodes. Admitting this and using again Theorem 2, this improves the result of Corollary 1 about the approximation ratio from 1/3 to 2/5. More generally, these results tend to confirm that the best approximation could be 1/2.

Some upper bounds

In this section, we describe instances which give upper bounds on the ratio ρ such that a ρ-approximate matching is guaranteed to exist. In particular, we show that there is no ( 12 + ε)approximate algorithm, for all ε > 0, and that there is is no ( 23 + ε)-approximate solution if we restrict ourselves to large graphs.

Let us define a family of instances containing n = 4k nodes for any positive integer k (see Figure 8 for an illustration with k = 2), such that each instance does not admit any ( 12 + ε)approximate matching, for all ε > 0. The nodes are denoted by {1, . . . , 2k} ∪ {1 ′ , . . . , 2k ′ }. The graph is complete and the edges are partitioned as follows :

-E w = {{i, i ′ }|i ∈ {1, . . . , 2k} and i ′ ∈ {1 ′ , . . . , 2k ′ }}. Each edge of E w has value (1/2k, 1/4k).

-E l = {{i, i + 1}|i ∈ {1, . . . , 2k} and i is odd}. Each edge of E l has value (0, 1/2k).

-E l ′ = {{i ′ , i ′ + 1}|i ′ ∈ {1 ′ , . . . , 2k ′ } and i ′ is odd}. Each edge of E l ′ has value (0, 1/2k). -E 0 the set containing the remaining edges. Each edge of E 0 has value (0, 1/4k).

We can take the 2k edges (i, i ′ ) and have a matching of maximal weight. We can take the 2k edges with value (0, 1/2k) and have a matching with maximal length. The maximal weight and length are both equal to 1. It is not difficult to see that the triangle inequality is satisfied on both coordinates.

A Pareto optimal matching consists in picking x ≤ k edges of value (0, 1/2k) in E l , x edges of value (0, 1/2k) in E l ′ and 2(kx) edges of value (1/2k, 1/4k) in E w . The resulting matching has weight k-x k and length k+x 2k . If k/3 is integral, then x = k/3 gives min{ k-x k , k+x 2k } = 2/3; we have an upper bound of 2/3. If k/3 is not integral, then taking x = 0 when k < 3 and x = ⌈k/3⌉ when k ≥ 3 allows us to maximize min{ k-x k , k+x 2k }. This gives an upper bound of max{1/2, 2/3 -1/k}.

We

Fig. 1 .

 1 Fig. 1. At least 3 solutions are required if the edge valuations do not fulfill the triangle inequality.

Fig. 2 .

 2 Fig. 2. At least 2 solutions are required if one edge valuation does not fulfill the triangle inequality.

Fig. 3 .

 3 Fig. 3. At least 3 solutions are required with 3 edge valuations fulfilling the triangle inequality.

Fig. 4 .

 4 Fig. 4. The labels of the edges when ni = 12 and k = 3. Bold edges belong to Mw and dashed edges belong to M ℓ .

Fig. 5 .

 5 Fig. 5. On the left, the labels of the edges of Ci when ni = 10 and k = 3 (Bold edges belong to Mw and dashed edges belong to M ℓ ). On the right, the path P0 when the edges with label 3 for matching L3(Ci ∩ Mw) have been deleted.

Fig. 6 .

 6 Fig.6. The labels of the edges of P i 0 ∩ M ℓ when ni = 10 and k = 3 (Bold edges belong to Mw and dashed edges belong to M ℓ ) for the example described in Figure5.

Fig. 7 .

 7 Fig. 7. The four cases n = 3, 4, 5, 6. Bold edges belong to Mw and dashed edges belong to M ℓ .

  by the triangle inequality on ℓ and w, ℓ(1, 2) < ℓ(2, 3)/2 and w(2, 3) < w(1, 2)/2 give ℓ(1, 3) ≥ ℓ(2, 3)/2 and w(1, 3) ≥ w(1, 2)/2, respectively.

  ),[START_REF] Aggarwal | Minimal spanning tree subject to a side constraint[END_REF][START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF],[START_REF] Stein | On the existence of schedules that are near-optimal for both makespan and total weighted completion time[END_REF][START_REF] Rasala | Existence theorems, lower bounds and algorithms for scheduling to meet two objectives[END_REF]} and M ℓ = {(1, 6), (2, 3), (4, 5)}. Consider the three following matchings M 1 = {(1, 2), (3, 6), (4, 5)}, M 2 = {(1, 6), (2, 5), (3, 4)}, and M 3 = {(1, 4), (

1 2 ℓCorollary 1 .

 21 (M ℓ ) and ℓ(M 5 ) ≥ 1 2 ℓ(M ℓ ). If w(M 2 ) ≥ 1 2 w(M w ) then M = M 2 . If w(M 2 ) < 1 2 w(M w ) then, by Lemma 1, we have w(M 5 ) ≥ 1 2 w(M w ) and thus M = M 5 . ⊓ ⊔ We can build in polynomial time a single matching which constitutes a 1/3approximate Pareto set for the biobjective Maximum Matching problem satisfying the triangle inequality on each objective. Proof. By using Theorem 2 with ρ = 1/2, k = 3 and Lemma 2 the result follows.⊓ ⊔

  * n , for any instance I ∈ I ∆ n , there exists M ∈ M n such that w(M ) > ρ and ℓ(M ) > ρ. Moreover, ρ * n , corresponding to a feasible solution, is such that w(M ) ≤ ρ * n or ℓ(M ) ≤ ρ * n for all matchings M ∈ M n . It follows that ρ * n is the largest value such that, for any instance I ∈ I ∆

n , there exists M ∈ M n such that w(M ) ≥ ρ * n and ℓ(M ) ≥ ρ * n

  ).Table

	n	3	4	5	6	7	8	9	10
	ρ * n	0.5	0.5	0.5 0.5556 0.5714 0.5 0.5714 0.6
		1/2	1/2	1/2	5/9	4/7	1/2	4/7	3/5

Table 2 .

 2 get the upper bounds (UB) provided in table 2 when the number of nodes n is equal to 4k. Upper bounds on the best approximation ratio for a tradeoff when n = 4k.

	k	1	2	3	4	5	6	7
	n	4	8	12	16	20	24	28
	UB 1/2	1/2	2/3	1/2	1/2	2/3 11/21

For k objectives functions to be maximized, s * weakly dominates s when fi(s * ) ≥ fi(s) for all i ∈ {1, . . . , k}.

This bi-objective problem generalizes max cut.

For every triple of edges a, b and c forming a triangle, we have that w(a) + w(b) ≥ w(c) and ℓ(a) + ℓ(b) ≥ ℓ(c).

The double factorial r!! is the factorial r! restricted to numbers of the same parity as r, e.g.

5!! = 5 * 3 * 1 and 8!! = 8 * 6 * 4 * 2.

) (0, 1 8 )

Fig. 8. An 8-node instance where no (1/2 + ǫ)-approximate matching exists. Dotted edges have value (0, 1 8 ), bold edges have value (0, 1 4 ) and dashed edges have value ( 1 4 , 1 8 ).

Concluding remarks

In this article we were interested in determining, and computing, an approximate tradeoff for the bi-objective matching problem. We focused on the case where each objective satisfies the triangle inequality. We proposed a general method which, combined with results stating the existence of a good tradeoff in graphs of small size, implies the existence and computation of a good tradeoff in a graph of any size. For the existence of a good tradeoff in graphs of small size, we followed two approaches. The first one is analytical but it is limited to graphs containing at most 6 nodes. We believe that such an analysis can be conducted for slightly larger graphs. The second approach is a computer-aided determination of an approximation ratio ρ * n , such that any complete graph on n nodes must admit a ρ * n -approximate matching. The values manipulated in the computational approach have a limited precision so we cannot pretend that the given ratios hold in any case. However, they provide an insight in what ρ * n should be equal to. More generally, we conjecture that a 1/2-approximate matching exists in any bi-objective instance. We followed a constructive path but an interesting future direction is to consider this question, but only under an existential approach.

Our last observation deals with how ρ * n depends on the number of nodes n. From the results given in sections 6 and 7, we observe that ρ * n is surprisingly non-monotone with n.


Pedro Rodrigues
email: pedro.rodrigues@centraliens.net

Benedetta Franzelli

Ronan Vicquelin

Olivier Gicquel

Nasser Darabiha

Unsteady dynamics of PAH and soot particles in laminar counterflow di↵usion flames

Keywords: Soot, PAH, Laminar flame dynamics, Soot sectional model, Particle size distribution

Unsteady dynamics of PAH and soot particles in laminar counterflow diffusion flames

Introduction

Due to incomplete combustion, soot emissions have e↵ects on both human health and environment. Soot emissions are also considered as an important cause of global warming [1]. Consequently, important e↵orts are made both experimentally and numerically [2,3,4] to understand soot production mechanisms in order to control their emission.

Most of the combustion facilities are characterized by high Reynolds number flames where turbulent eddies are expected. The local strain rate usually fluctuates in a wide amplitude range and with random fluctuation frequencies [5]. These turbulent eddies are also responsible for variable length scale recirculation zones, introducing a wide range of residence times for soot particles, strong intermittency and dynamics features in soot production [6,7].

One of the most popular approaches used to simulate turbulent non-premixed flames is the flamelet approach, based on a quasi-steady response of the flame characteristics to the local strain rate fluctuations [8,9].

In the optic of applying such models to numerical simulations of turbulent flames, the response of soot to strain rate fluctuations can be investigated by looking at unsteady laminar counterflow di↵usion flames [10,11,12]. Specifically to soot context, previous experimental works have been performed in a di↵usion laminar flame by introducing sinusoidal velocity variations at both opposed nozzles [13,14] . They showed that soot production response to these fluctuations was phase-lagged and damped when increasing the oscillation frequency. A particular hierarchical behavior was observed: soot volume fraction response is more phaselagged and damped compared to soot precursors response, which are also more phase-lagged and damped than the temperature response [15]. Cuoci et al. [16] numerically investigated these flames with good prediction of unsteadiness soot dynamics, confirming the experimental observations. Nevertheless, a lack of knowledge remained on the origin of soot response to unsteady strain fluctuations. Moreover, when computing counterflow di↵usion flames with unsteady velocities at the nozzle exits, a phase lag exists between the global strain rate and the local strain rate [16], increasing the complexity of the phenomena.

The objective of the present work is to characterize the response of soot to strain rate oscillations and to identify the physical phenomena underlying the phase lag and damping observed in soot production. In order to avoid the phase lag between the global and the local strain rate, a strain-imposed formulation is considered in this work and unsteadiness is introduced by varying the imposed flame strain rate a(t) with time for a given pulsation !, an initial strain rate A 0 and fluctuation amplitude a 1 :

a(t) = A 0 + a 1 sin(!t) = A 0 ⇥ 1 + ↵sin(2⇡ f t) ⇤ . (1)
Both analytical and numerical approaches are considered in this paper to study the evolution of the soot precursors and of the particle size distribution (PSD) with the strain rate a(t).

The paper is organized as follows. First, an analytical model is proposed in Section 2 in the limit of a linear behavior, i.e. small oscillation amplitudes. This model predicts the unsteady response on the basis of steady flame results. Then, soot production in unsteady laminar flames is numerically studied using a detailed sectional model. The modeling strategy is introduced in Section 3. The flame response is then investigated for the configuration described in Section 4.1. The unsteady behavior is analyzed in Section 4.2 for di↵erent frequencies at small amplitude in terms of global quantities and PSD. Analytical results will be compared to the numerical ones in Section 4.3 to prove their validity. The causes of phase lag and damping in soot production will then be identified by combining information from numerical and analytical results. Finally, numerical simulations at high amplitudes are analyzed in Section 4.4 to completely characterize the soot response to unsteady strain rate oscillations and to discuss the limits of the analytical model.

Analytical model for pulsed sooted flames

In order to investigate the response of soot production to strain rate fluctuations, a linearized analytical model is developed in the following to predict the response of the maximum of a flame variable ✓ to strain rate oscillations at a given pulsation !. The complex form of the fluctuating strain rate a 1 (t) = a(t) A 0 is denoted by â1 (!) = ↵A 0 e i(!t+⇡/2) . The corresponding response of the maximum value of

✓, namely ✓ max (t) = ✓ max 0 + ✓ max 1 (t) with ✓ max 1 (t) = ✓ max 1 (!)sin(!t ' ✓ max (!)), is represented by the com- plex number ✓max 1 (!) = ✓ max 1 (!)e i(!t+⇡/2 ' ✓ max (!))
. This response is fully characterized by the transfer function T ✓ max (!) = ✓max 1 (!)/â 1 (!). Starting from the previous works [10,11,12], the transfer function is split into two terms: the transfer function T finite,✓ unst (!), introducing an equivalent steady strain rate A ✓ seen by the quantity ✓, and the transfer function T ✓ max |A ✓ steady (!), describing the response of ✓ max to the equivalent steady strain rate A ✓ .

Equivalent steady strain rate

Following [10,11,12], under the assumption of infinitely fast chemistry, the unsteady flame acts at each time t as an equivalent steady counterflow flame at constant strain rate equal to the instantaneous strain rate A(t) verifying:

dA dt = 2A 2 (t) + 2A(t)a(t). (2)
Assuming a linear response of A(t) = A 0 + A 1 (t) with a(t) = A 0 + a 1 (t), i.e. small fluctuations of a(t) around A 0 , the transfer function T inf unst (!) between Â1 (!) and â1 (!) in the case of infinitely fast chemistry is given by:

T inf unst (!) = 1 1 + j!/(2A 0 ) . (3)
When finite-rate chemistry is considered, the equivalent strain rate A ✓ for a given variable ✓ (T or Y k ) is given by [10]:

@A ✓ @t = A ✓ (t) A(t) A ✓ (t) ⌦✓ (t) (d✓/dA) ⌦ = A(t) A ✓ (t) ✓ (t) , (4)
where ⌦✓ = !max T /(⇢c p ) for ✓ = T and ⌦✓ = W k !max k /⇢ for ✓ = Y k . ⇢ and c p are evaluated at the position where !✓ is maximum. (d✓/dA) ⌦ represents the variation of ✓ where !✓ is maximum with a steady strain rate. ✓ (t) is defined as ✓ (t) = (d✓/dA) ⌦ A ✓ (t)/ ⌦✓ (t).

To find the linearized response of A ✓ (t), A ✓ (t) and ✓ (t) are written as:

A ✓ (t) = A 0 + A ✓ 1 (t) and ✓ (t) = ✓ 0 + ✓ 1 (t)
, where A 0 , ✓ 0 , are the values of respectively A ✓ (t), and ✓ (t) for the initial steady flame. By linearizing Eq. ( 4), one obtains:

@A ✓ 1 @t = A 1 (t) A ✓ 1 (t) 1 ✓ 0 . (5)
Combining the Fourier transform of Eq. ( 5) and Eq.

(3), the following transfer function T finite,✓ unst (!) between Â✓ 1 (!) and â1 (!) is obtained:

T finite,✓ unst (!) = Â✓ 1 (!) â1 (!) = 1 1 + j! ✓ 0 1 1 + j!/(2A 0 ) . (6)

Steady response of the maximum value

Once the equivalent steady strain rate A ✓ (t) is known, the flame response can be analyzed by looking at steady conditions. For ✓ 2 {T, Y k }, it has been observed that in the neighborhood of a given strain rate A 0 , the dependency of ✓ max for a steady flame with strain rate A is given by [17]:

✓ max (A)/✓ max (A 0 ) = (A/A 0 ) p ✓ (7)
with p ✓ a characteristic constant. Linearizing Eq. ( 7) with

✓ max (t) = ✓ max 0 + ✓ max 1 (t) gives ✓ max 1 (t) = p ✓ ✓ max 0 A 0 A ✓ (t)
. However, the forthcoming comparison with the detailed computation demonstrates the requirement to introduce an additional delay in the response to the unsteady strain rate oscillations. Linking this delay to the chemical time seems particularly relevant for soot precursors and particles, whose chemistry is mainly sequential so that all the reactions necessary for the formation have to respond before getting the response of ✓ max . This delay is then assumed to be equal to the characteristic chemical time scale ⌧ ✓ of the quantity of interest ✓ defined in Appendix A: ✓ max (t) reacts then at the equivalent strain rate A ✓ (t ⌧ ✓ ). The validity of this hypothesis will be verified in Section 4.3. The response of ✓ max 1 (t) is therefore expressed as:

✓ max 1 (t) = p ✓ ✓ max 0 A ✓ (t ⌧ ✓ )/A 0 ) T ✓ max |A ✓ steady (!) = p ✓ ✓ max 0 e j!⌧ ✓ /A 0 ( 8
)
where T ✓ max |A ✓ steady (!) = ✓max 1 (!)/ Â✓ 1 (!) represents the transfer function between ✓max 1 (!) and Â✓ 1 (!).

Transfer function T ✓ max (!)

From the definitions of ✓ 0 and ⌧ ✓ , ✓ 0 can be rewritten as ✓ 0 = ⌧ ✓ ✓ with ✓ = (d✓/dA) ⌦ • (A 0 /✓ max (A 0 )), a dimensionless parameter characterizing the steady response of the quantity ✓ to strain rate. Then, by combining Eqs. ( 6) and ( 8), the transfer function T ✓ max (!) between ✓max 1 (!) and â1 (!) is given by T ✓ max (!) = T ✓ max |A ✓ steady (!)T finite,✓ unst (!). Gain and phase lag of ✓ max are expressed respectively by

G ✓ max (!) = 20log 10 (|T ✓ max (!)| / |T ✓ max (! = 0)|) and ' ✓ max (!)
, where:

8 > > > > > > > > > > > < > > > > > > > > > > > : |T ✓ max (!)| = a 1 p ✓ ✓ max 0 A 0 1 p 1 + ⌘(!) 2 1 r 1 + ✓ 2⌘(!) ✓ Da ✓ ◆ 2 ' ✓ max (!) = tan 1 (⌘(!)) + tan 1 (2⌘(!) ✓ /Da ✓ ) + 2⌘(!)/Da ✓ (9)
with

⌘(!) = !/(2A 0 ) = ⇡ f /A 0 the Stokes number and Da ✓ = A 1 0 ⌧ 1
✓ the Damköhler number associated with ✓. From Eq. ( 9), it can be deduced that three nondimensional parameters are responsible for the phase lag and damping of the response of ✓:

• ⌘ compares the characteristic time associated to the strain rate A 0 to the imposed frequency f and is responsible for the filtering of the flow structure. The damping response of all the quantities increases when ⌘ increases.

• Da ✓ is directly responsible for the phase lag and damping response due to the low chemical time scale of the analyzed quantity ✓. The lower Da ✓ is, the more the response of ✓ max is phase-lagged and damped.

• ✓ , which represents the steady response of the quantity ✓ to strain rate, also contributes to the damping response of ✓ with unsteady strain fluctuations. For high values of ✓ , the damping response will be high.

This identified behavior is valid for all the quantities but is more significant in the case of species with large chemical characteristic time scales (Da ✓ ⌧ 1), which is the case of soot precursors and particles. Equation (9) allows the predictions of the unsteady response of ✓ from information on steady flames. Its validity will discussed, in particular for soot production, in Sec. 4.3.

Detailed modelisation of soot production

Parallely to the asymptotic analysis, the behavior of pulsed laminar di↵usion flames is investigated numeri-cally. In order to obtain an accurate numerical prediction of soot and its precursors, detailed models for both gas and solid phase are considered in this work and described below.

Sectional method for solid phase

The soot population is evaluated by using a sectional method. Each section i represents particles with a volume between V min i and V max i . The soot mass fraction Y s,i of the i th section is given by the following transport equation:

@⇢Y s,i @t + r • (⇢(u + v T )Y s,i ) = r • (⇢D s,i r(Y s,i )) + ⇢ s Qs,i (10)
where ⇢ is the gas phase density, u is the gas velocity, v T is the thermophoretic velocity of the particles given in [18], D s,i is the di↵usion coe cient of particles of the i th section defined in [19], ⇢ s is the constant soot density (chosen equal to ⇢ s = 1860 kg/m 3 ) and Qs,i is the production rate of the soot volume fraction for the i th section accounting for nucleation, condensation, surface growth, oxidation and coagulation.

The models used to close this production rates are based on those used by Karkar el al. [20]. Several improvements have been made in the present work and are presented below.

Nucleation corresponds to the formation of the smallest solid particles. Here, the coalescence of two dimers is considered for the formation of these particles [21]. Condensation is considered as the coalescence of a dimer at a soot particle surface.

Dimers are formed from the collision of two polycyclic aromatic hydrocarbons (PAH). According to [22], only the collision of PAH with four-aromatic rings and more is considered. Here, seven PAH have been considered, from the pyrene (A4) and up to the coronen (A7). A quasi-steady-state hypothesis is considered between their production from the gaseous phase and their consumption by nucleation and condensation [20].

Soot surface reactions are responsible for both soot particles surface growth and oxidation. These phenomena are described through the HACA-RC mechanism [23]. The oxidation reaction by OH has been updated based on recent experimental results of [24].

According to previous works [21], big soot particles can not be considered as spherical particles. A soot particle of a given volume V and surface S is considered as a fractal aggregate of n p = S 3 / ⇣ 36⇡V 2

⌘

primary spherical particles with a diameter d p = 6V/S . For each soot particle, S is estimated by fitting numerical results from [21] as a function of 3 . V 1 denotes the volume from which a soot particle is no longer considered as spherical. S C 2 and V C 2 are respectively the surface and volume of a spherical molecule composed of two carbon atoms.

V by S /S C 2 = V/V C 2 2/3 for V < V 1 and S /S C 2 = V/V C 2 ✓(V)/3 for V > V 1 , with ✓(V) = 2 + 0.175log 10 (V/V 1 ) and V 1 = 320 nm
Particle nucleation and condensation as well as coagulation are described through the Smoluchowski equation [START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen[END_REF]. This equation is expressed as a function of the collision diameter d c of the soot particles, calculated as in [21] as a function of d p , n p and their fractal dimension D f (chosen equal to D f = 1.8).

Gaseous phase description, radiation model and solving strategy

The detailed kinetic scheme KM2, due to [22], has been considered in this study. It involves 202 species and 1 351 reactions, and has been validated for the estimation of PAH up to coronen.

A radiative source term has been added in the energy transport equation. It is computed at each point as a function of the absorbed and the emitted radiative powers [START_REF] Soufiani | [END_REF], which are expressed as a function of the corresponding detailed intensities. Radiative properties are expressed via a narrow-band model [START_REF] Soufiani | [END_REF]. For the gaseous species, CO 2 , H 2 O and CO are considered as the main contributors in radiation fluxes. For soot particles, the absorption coe cient  ⌫,soot = 5.5⌫ f V is estimated as a function of the soot volume fraction f V and the wave number ⌫ [27].

The above models as well as Eqs. (10) are introduced into a system of 1-D equations describing a counterflow configuration based on self similar approximation and imposed variable strain rate or injection velocity [28]. The coupled gas and soot sections transport equations are solved using the REGATH package [29].

Validation test cases

The proposed modeling strategy has been first validated on soot prediction for:

• PSD description (not shown): numerical results show a fair agreement with the experimental data [30] obtained with the BSS (Burner Stabilized Stagnation) flame technique for laminar premixed ethylene flames.

• Response to steady strain rate (Fig. 1): the response of soot volume fraction as a function of the global steady strain rate has been reproduced numerically on the steady counterflow di↵usion flame experimentally investigated by Decroix et al. [13] for three fuels. • Response to unsteady strain rate: an unsteady counterflow propane/air di↵usion flame has been investigated by imposing an oscillating velocity at the injection as experimentally done in [13,14].

A phase lag of 125 between the minimum average soot volume fraction and maximum imposed velocity was experimentally observed for A 0 = 15s 1 , ↵ = 30% and f = !/2⇡ = 25Hz. With the numerical computation, a phase lag of 114 was obtained.

The good agreement of the numerical results with experimental data confirms the validity of the retained modeling strategy.

4. Detailed simulations of soot production in unsteady laminar di↵usion flames at imposed strain rate

Numerical configuration

A counterflow propane/air di↵usion flame is considered here by varying the imposed strain rate a(t) from an initial flame at A 0 = 60 s 1 . Ten frequencies and three amplitudes have been considered. Pure propane and pure air stream, both at 294 K are supplied through the two opposed nozzles at a distance L = 12.7 mm, discretized with more than 400 points. For each studied frequency, ten signal periods were computed. Once the permanent regime was attained, the response of each variable was studied in terms of gain and phase lag.

PAH and soot particles response

Results for small strain rate fluctuations (↵ = 10%) are considered here.

Figure 2 (left) presents the unsteady response of the soot maximum volume fraction and pyrene (A4) maximum mass fraction (the smallest considered soot precursor) to the unsteady imposed strain rate during two oscillating cycles. Quantities have been normalized with their respective steady values at the lowest and highest strain rates for three frequencies. The higher the frequency, the more f max V and Y max A4 fluctuations are dumped and phase-lagged. Looking at the results in the a-space (Fig. 2, right) enables a clear comparison with the quasi-steady solution (grey line). A quasi-steady response is observed at low frequency ( f = 0.1Hz), while for higher frequencies, solutions step aside from the steady results. The temporal evolution of the PSD is also studied here by looking at the four instants A.,B.,C.,D. of Fig. 2 separated by 90 in one pulsation period. Results for three frequencies are presented in Fig. 2 together with the quasi-steady state at the spatial position x soot , where soot volume fraction is maximum, close to the stagnation point. At each time, the PSD shows a bi-modal nature with one peak for small particles (generated by nucleation) and another for large aggregates (due to condensation and coagulation). In the quasi-steady case, from point A. to C. the characteristic flow time decreases (since a(t) increases), so that particles have less time to coagulate. The position of the aggregates peak translates then towards smaller diameter values: the higher the strain rate, the smaller are the aggregates composing the soot population. Inversely, from point C. to A., the strain rate decreases, particles have the time to coagulate and bigger aggregates populate the PSD.

The unsteady PSDs follow such dynamics, but their response is a↵ected by the phase-lag already observed on the global f V . Indeed, at f = 5Hz, the PSD responds in a quasi-steady way, whereas the phase-lag e↵ect is more and more evident on PSD for higher frequencies. The response of the PSD is also more and more damped so that at f = 60Hz only small PSD fluctuations are observed between the four instants. For high oscillation frequencies, the PSD is observed to not oscillate anymore since the oscillations are completely damped.

Comparison with analytical predictions

In order to understand the processes governing the PSD evolution, results for the di↵erent sections are now investigated.

Figure 4 presents the response in terms of gain and phase lag of maximum temperature, Y max A2 , Y max A4 , maximum soot mass fraction of two sections (sections 12 and 16, whose mean diameter are indicated in Table 1) and f max V . The response of precursors and soot is more phase-lagged and damped than temperature. Moreover, phase-lag and damping increases with their size (not shown for all precursors). Big particles are the main contributions to soot volume fraction, so that f V response is mainly governed by the last soot sections.

A good agreement is obtained between the numerical results (lines) and the analytical model (symbols) described in Section 2.2. Discrepancies are mainly observed at high frequencies but the hierarchical behavior between temperature, soot precursors and soot sections is well predicted. This confirms that soot dynamics are mainly governed by the three parameters identified with the analytical model. In particular, soot response is mainly due to its slow chemistry compared to the flame. To identify the main physical processes contributing to such a long chemical time, the characteristic time scales for nucleation (⌧ nu ), condensation (⌧ cond ), surface growth (⌧ sg ) and coagulation (⌧ coag ) have been estimated for di↵erent sections from the steady flame at A 0 following Appendix A. Table 1 presents these characteristic time scales normalized by the flame time scale (⌧ T = 0.31 ms) for five soot sections. All the characteristic time scales increase with the soot particle size, in particular for ⌧ cond and ⌧ coag which depend on the collisions rate. The particle number density of the last sections being smaller than for small particles sections, the number of particles available for collision is lower so that the characteristic time scales of collisional phenomena increases with the particle size.

The long characteristic time scale of f max V , governing the phase-lag and damping of the unsteady response, is then mainly due to condensation and coagulation phenomena of the biggest particles. represents the mean collisional diameter of a soot particle in the i th section.

Section d mean c,i (nm) ⌧ nu ⌧ T ⌧ cond ⌧ T ⌧ sg ⌧ T ⌧ coag ⌧ T 1 1.2 3. 4
In order to study the validity of the assumption on the induced delay time due to slow chemistry presented in Section 2.2, Fig. 5 presents the obtained numerical phase lag due to this delay time as a function of the ex-pected one. This phase lag is obtained by substracting to the obtained numerical phase lag ' num ✓ max (!) the theoretical phase lag of the equivalent steady strain rate A ✓ of the variable of interest ✓. Good results are obtained for all the species and temperature, which confirms the representativity of the chosen variable (⌧ ✓ ). Nevertheless, some discrepancies exist and future investigations are still needed in order to define in a more precise way this delay time.

Numerical results at high amplitudes

In order to study the soot dynamics at higher amplitudes, computations have been performed for amplitudes ↵ of 30% and 60%. Table 2 compares the numerical results for the phase lag and amplitude gain of Y max A4 and f max V for three amplitudes and three frequencies. The phase lag increases with the frequency in a similar way for all the amplitudes. The gain remains almost the same for ↵ = 10% and ↵ = 30%, but decreases for ↵ = 60%. The di↵erence of the numerical behavior between ↵ = 30% and ↵ = 60% highlights the non-linear e↵ects for such amplitudes, which cannot be described by the linear analytical model whose predictions do not depend on the perturbation amplitude.

Conclusions

Response of sooting propane-air counterflow di↵usion flames to imposed strain rate harmonic oscillations were numerically investigated with a detailed description for the gas and the solid phases. The unsteady behavior of soot particles and precursors production, as well as the PSD evolution, were studied both analytically and numerically. It has been observed that the higher the oscillation frequency is, the more PAHs and soot particles fluctuations are damped and phase-lagged so that unsteady solutions are farther and farther away from the quasi-steady state. The phase-lag and damping increase with the size of PAHs and soot particles.

f ✓ ↵ = 10% ↵ = 30% ↵ = 60% (Hz) G ✓ max ' ✓ max G ✓ max ' ✓ max G ✓ max ' ✓ max 30 A4 - 4
An analytical model has been proposed to predict the observed phase lags and dampings assuming a linear behavior. Three non-dimensional parameters (⌘, Da ✓ and ✓ ) govern the unsteady response. Soot particles are characterized by long time scales mainly due to condensation and coagulation phenomena. Indeed, compared to the gas species, their dynamics, particularly the additional identified phase lag, are mainly governed by the Da ✓ parameter.

Therefore, models developped for numerical simulations of soot production in turbulent flames have to correctly reproduce these observed features in order to represent unsteady behaviors such as soot intermittency. On the one hand, these behaviors highlight the limits of flamelet regime assumption based on quasi-steady hypothesis, implying major complexities in modeling for turbulent calculations. In this sense, the presented results support the need for specific techniques [31] to account for PAHs response to unsteady strain rate fluctuations. On the other hand, the reduced models have to provide a good prediction of ⌘, Da ✓ and ✓ for PAHs and soot. As an example, representative soot precursors have to be chosen in terms of these three parameters in order to obtain the good unsteady behavior of soot production: large precursors dynamics (such as pyrene and coronen) have to be reproduced. The proposed analytical model will be very useful for the development of models that reproduce the dynamics of soot and their precursors in turbulent flames.

AppendixA. Chemical characteristic time scales

Species and flame characteristic time scales.

To study the unsteady response of each chemical species separately, the characteristic time scale ⌧ k for the k th species is defined as [32]:

⌧ k = [X k ] max / !max k = (⇢Y k ) max / ⇣ W k !max k ⌘ (A.1)
where [X k ], Y k , W k and !k are the molar concentration, the mass fraction, the molecular weight and the molar production rate of the k th species.

In the same way, the characteristic time scale ⌧ T of the flame can be defined as ⌧ T = ⇣ ⇢c p T ⌘ max / !max T , with c p the mixture mass specific heat capacity and !T the heat release rate.

Soot sections characteristic time scales.

As for the species, a time scale ⌧ s,i for the soot particles in the i th section can be defined as:

⌧ s,i = (⇢Y s,i ) max / ⇣ ⇢ s ( Qs,i ) max ⌘ (A.2)
By perturbing each volume fraction production rate relative to each phenomenon (ph) by a small value (typically 1%), the characteristic time scales ⌧ nu,i for nucleation, ⌧ cond,i for condensation, ⌧ sg,i for surface growth, ⌧ ox,i for oxidation and ⌧ coag,i for coagulation can be expressed as:

⌧ ph,i = ⇥ (⇢Y s,i /⇢ s ) max ⇤ / h ( Qph,i ) max i (A.3)
where h ( Qph,i ) max i corresponds to the variation of the peak volume fraction production rate of the phenomenon for the i th section.

Figure 1 .

 1 Figure 1. Maximum soot volume fraction f max V as a function of the global strain rate. Comparison between the present numerical (lines)and experimental (symbols) data from[13].

Figure 2 .

 2 Figure 2. Normalized response of soot maximum volume fraction ( f max V ), pyrene maximum mass fraction (Y max A4 ) to the unsteady imposed strain rate (a(t)).

Figure 3 .

 3 Figure 3. Unsteady variations of the PSD at f max V position.

Figure 4 .

 4 Figure 4. Comparison between analytical model predictions (lines) and numerical results (symbols) of amplitude gain and phase lag for maximum temperature, naphtalene (A2) and pyrene (A4) maximum mass fractions, maximum mass fractions of the 12 th , the 16 th soot sections and maximum soot volume fraction. Analytical results for T max and Y max A2 are superposed.

Figure 5 .

 5 Figure 5. Numerical results of the induced phase lag due to the delay time imposed by slow chemistry as a function of the expected ones.

Table 1 .

 1 Comparison of normalized characteristic time scales of nucleation, condensation, surface growth and coagulation. d mean

	3.2 3.6	-

Table 2 .

 2 Numerical analysis of the impact of the strain fluctuation amplitude (↵) on the pyrene maximum mass fraction and soot maximum volume fraction gain (G in dB) and phase lag (' in deg).

	30 60 60 120 A4 f V A4 f V 120 f V	54 146 87 -14 269 -14 287 -4 53 -4 -4 158 -13 -5 -9 -8 91 -10 -24 -16 126 -16 134 -18 -33 460 -33 524 -41	56 142 87 256 130 497