Nathanaël Barrot
email: nathanael.barrot@lamsade.dauphine.fr

Jérôme Lang
email: jerome.lang@lamsade.dauphine.fr

Conditional and Sequential Approval Voting on Combinatorial Domains

Several methods exist for making collective decisions on a set of variables when voters possibly have preferential dependencies. None is based on approval voting. We define a family of rules for approval-based voting on combinatorial domains, where voters cast conditional approval ballots, allowing them to approve values of a variable conditionally on the values of other variables. We study three such rules. The first two generalize simple multiwinner approval voting and minimax approval voting. The third one is an approval-based version of sequential voting on combinatorial domains. We study some properties of these rules, and compare their outcomes.

Introduction

Collective decisions on combinatorial domains cover a large variety of common situations, such as finding a set of dates for a series of meetings, electing a committee of representatives, adopting or rejecting each of a series of yes-no questions about common facilities to be built, or finding a common set of movies for a group to watch together.

Several methods have been proposed and studied for making such collective decisions. An important source of problems is when voters have preferential dependencies between variables: someone may be willing to attend the second meeting on a Friday provided that the first meeting is not already on a Friday; someone may want Ann to be elected at the department council only if both Betty and Charles are not elected too (they all belong to the same research group and would have too much joint power); I may want to watch one comedy and one drama, but not two films of the same genre.

It has been argued many times that making decisions independently on each variable cannot take such preferential dependencies into account. If we want to deal with them, we have to allow voters to express these dependencies. Several classes of methods have been proposed and studied (see [Lang and Xia, 2016] for a recent survey). Direct methods work by asking voters a one-shot report of their preferences, expressed in some language that allows preferential dependencies to be (at least partly) specified; sequential methods proceed by eliciting the voters' preferences on a first (fixed) variable, decide its value, broadcast it, then eliciting the voters' preferences on a second variable conditionally on the value chosen for the first one, and so on.

Most of the existing methods for voting in combinatorial domains assume that voters express, implicitly or explicitly, rankings over the values of variables, or else cardinal utility functions. Both may induce a cognitive burden to the voters, as well as an increased complexity in communication and computation. On the other hand, in simple domains, approval voting is a well-studied method for elections, where each voter may approve of any number of candidates, and the candidate with the largest number of approvals wins. Approval voting has nice properties, and is easy to use and understand [Brams and Fishburn, 2005].

Approval voting can be extended to combinatorial domains, and especially to committee elections. The simplest way, if we want to select k candidates, consists in selecting those with the largest k approvals. Other ways of extending approval voting to multiple winner elections have been considered and are surveyed in [Kilgour, 2010]. Similar methods can be used for multiple referenda, as argued in [START_REF] Amanatidis | [END_REF], due to the structural similarity of committee elections and multiple referenda: in both cases, a decision has to be made over a set of possible subsets (of candidates in one case, of binary issues in the other case) subject to feasibility constraints (such as the number of winners in committee elections, or budget constraints in multiple referenda). However, the methods we know of for multiwinner approval voting do not allow voters to express conditional preferences; they must vote as if they had separable preferences.

We propose a family of methods for reconciliating the simplicity of approval voting with the possibility for voters to express preferential dependencies. In Section 3 we define conditional approval ballots where voters approve some values of variables conditionally on the value of some other variables. In Section 4 we define two first rules based on conditional approval ballots, where the winning alternative minimizes respectively the maximum or the sum, over all voters, of the number of 'disagreements', and we show that these rules generalize simple (or 'minisum') approval voting and minimax approval voting [Brams et al., 2007]. In Section 5 we define sequential conditional approval voting: given a fixed order of variables, the voters approve some of their values, one variable after the other. We conclude in Section 6.

Background and related work

In a committee (or multi-winner) election, there is a set of candidates C, a set of voters N , each casting a vote (or a ballot), and a number k of winners (or sometimes a more flexible constraint, or even no constraint at all such as in Hallof-Fame elections). In multiwinner approval ballotting, each voter casts an approval ballot where she approves as many candidates as she wants (in some variants, there are also constraints on the number of approvals on a ballot).

In simple multiwinner approval voting, the winners are simply the k candidates approved most often. Simple approval voting can be criticized for failing to guarantee a sufficient level of fairness or representativeness, and a number of other rules that map a collection of approval ballots to a set of winners have been studied; see [Kilgour, 2010] for a survey. One prominent such rule is minimax approval voting, where the elected committee minimizes the maximum, over all voters, of the number of 'disagreements' (i.e., formally, the Hamming distance) between the voter's ballot and the chosen committee [Brams et al., 2007]. Unlike simple approval voting, minimax approval voting is NP-hard [START_REF] Moti | [END_REF] and manipulable, but it can be efficiently approximated [START_REF] Legrand | [END_REF][START_REF] Caragiannis | Approximation algorithms and mechanism design for minimax approval voting[END_REF][START_REF] Byrka | [END_REF]; see [Misra et al., 2015] for its parameterized complexity. [START_REF] Amanatidis | [END_REF] generalize these two rules to other aggregation functions, and [START_REF] Baumeister | [END_REF] extend them beyond approval ballots.

Other multiwinner approval voting rules, such as proportional or satisfaction approval voting [Brams and Kilgour, 2014], have been recently studied from the point of view of computation [Aziz et al., 2015b] or of their properties [Aziz et al., 2015a]. A number of works address the computation of full proportional representation with approval ballots, two good representatives being [Procaccia et al., 2008] and [START_REF] Skowron | [END_REF]. The computational aspects of strategic behaviour for single-and multi-winner approval voting are adressed in [START_REF][END_REF] and in [START_REF] Baumeister | [END_REF]. The properties of multiwinner voting rules are studied more generally in [START_REF] Elkind | [END_REF].

Combinatorial domains are sets of alternatives consisting of the Cartesian product (or sometimes, a subset of it) of finite domain values corresponding to issues, variables, attributes, seats, or even individuals (in the case of committee elections). The main difficulty in voting in combinatorial domains is the presence of preferential dependencies. To deal with them, some approaches make use of compact representations, such as [Rossi et al., 2004;[START_REF] Conitzer | [END_REF][START_REF] Li | [END_REF]. Some others make use of sequential voting, variable after variable, such as [START_REF] Lang | [END_REF]dalla Pozza et al., 2011;Airiau et al., 2011]. None of them makes use of approval ballots. See [Lang and Xia, 2016] for a recent overview of voting in combinatorial domains.

Conditional approval ballots

Let X = {X 1 , . . . , X p } be a set of variables, each of them associated with a finite value domain D i . Let D = D 1 × . . . × D p , and D * ⊆ D be a set of feasible alternatives (by default, D * = D). A group of voters have to decide of a common outcome in D * . When X i is binary, we will usually take D i = {x i , x i }. For J ⊆ {1, . . . , p} we note D J = × j∈J D j , and for d = (d 1 , . . . , d p ) ∈ D, d J denotes the tuple d projected to the indices in J.

For instance, in a designated-post committee election [Benoît and Kornhauser, 2010], X is a set of different seats to be filled and their domains are the respective sets of candidates for each position (an alternative is feasible only if it does not assign the same candidate to more than one seat). In a standard committee election where the size of the committee is fixed to k, X can be identified with the set of candidates, with binary domains {elected, not elected}, and an alternative is feasible iff k variables have value 'elected'. In a multiple referendum [START_REF] Brams | [END_REF]], X is a set of binary issues.

As a running example, a group of friends have to decide on a common menu. The variables are X 1 (main dish), with D 1 = {m, f, v} (meat, fish, vegetarian dish), and X 2 (drink), with D 2 = {r, w, b} (red wine, white wine, beer).

The most common way of using approval ballotting in such domains consists in asking voters to approve values for each variable separately. This does not require much communication, but it is not possible for voters to report conditional preferences. For instance, a voter may approve {m, f } ⊆ D 1 and {r, w} ⊆ D 2 but cannot express that he approves the red wine only if the main dish collectively chosen is meat.

One could also think of asking voters to specify all combinations of values they approve (possibly in a compact way by means of a propositional formula). Here, preferential dependencies can be expressed, but the communication is costly, and moreover, a problem is that once a voter obtains a value he does not like for a variable, his opinion about the other variables does not count. Assume for instance that a voter hates fish and beer, and wants red wine with meat. It is reasonable to expect that he approves the set of menus {mr, vw, vr}, and disapproves mw, mb, f r, f w, f b, vb. However, if the collective decision turns out to give the value fish to X 1 , our voter won't have a chance to express that he still does not want to drink beer.

A trade-off between these two (extreme) methods consists in allowing voters to approve sets of values for each variable, conditioned on the value of some other variables. Definition 1. A conditional approval (CA) ballot over variables X 1 , . . . , X p with domains D 1 , . . . , D p is a pair

B = G, {A i | i = 1, . . . , p}
where G is a directed graph over {X 1 , . . . , X p }, and for each i, A i is a set of conditional approval statements {u :

a i (u)|u ∈ D P ar G (Xi) }, where a i (u) ⊆ D i , and P ar G (X i ) is the set of parents of X i in G. B ↓Xi|u = a i (u) denotes the projection of B over X i conditioned by u. We say that B is acyclic if G is acyclic. B is a separable CA ballot if G has no edges.
Most often, we will alleviate the notation and only write the CA statements (and not the graph) when we specify a CA ballot. Also, we omit some curly brackets when it cannot lead to any confusion (in particular, for singletons). Definition 2. Given an CA ballot B and an assignment d

= (d 1 , . . . , d p ) ∈ D 1 × . . . × D p , we say that d disagrees with B on variable X i if d i / ∈ B ↓Xi|d U , where U = P ar G (X i ).
The disagreement of d with respect to B, denoted by δ(d, B), is the number of variables X i on which d disagrees with B.

Example 1. Consider a voter expressing the CA ballot

B = [X 1 → X 2 ], {{f, v}; m : r; f : w; v : {r, w}}
This voter has a preference for not having meat: B ↓X1 = {f, v}; but in case the collective decision is to have meat, she still has her word to say about the drink, and approves (only) the red wine: B ↓X2 | m = {r}. If the collective decision is the vegetarian dish, then she approves both the red and the white wine: B ↓X2 | v = {r, w}. Now, d = mr disagrees with B on X 1 but not X 2 , therefore δ(d, B) = 1 (she is not satisfied with the choice of the dish, but is satisfied by the choice of the drink given that the dish is meat). Likewise, we have δ(f w, B) = δ(vr, B) = δ(vw, B) = 0 (she is satisfied by the choice of the dish and the choice of the drink conditional on the value of the dish),

δ(f r, B) = δ(f b, B) = δ(vb, B) = 1, and δ(mw, B) = δ(mb, B) = 2.
Equivalently, p -δ(•, B) can be seen as a utility function which counts the number of conditional preferences satisfied. Thus, even though voters are using approval-like ballots, the preferences they express are not dichotomous.

Observation 1. If B is acyclic and a i (u) = ∅ for all i ≤ p and u ∈ D P ar G (Xi) , then δ(d, B) = 0 for some d ∈ D. This, however, is no longer true if B is not acyclic: let for instance B = {x 1 : x 2 , x 1 : x 2 , x 2 : x 1 , x 2 : x 1 }, then for every d ∈ D we have δ(d, B) = 1.
Conditional approval ballots have a clear resemblance with CP-nets [Boutilier et al., 2004]. For the sake of simplicity, assume that variables are binary and let N (B) the CP-net obtained by replacing all conditional approval sets u : a i (u) by u : a i (u) a i (u), with This can be seen on this example:

x i = x i .
X = {X 1 , X 2 , X 3 }; B = {x 1 , x 1 : x 2 , x 1 : x 2 , x 1 x 2 : x 3 , x 1 x 2 : x 3 , x 1 x 2 : x 3 , x 1 x 2 : x 3 }. We have x 1 x 2 x 3 N (B) x 1 x 2 x 3 , because N (B) contains x 1 x 1 , however x 1 x 2 x 3 ≺ B x 1 x 2 x 3 , because δ(x 1 x 2 x 3 , B) = 1 < δ(x 1 x 2 x 3 , B) = 2. Note that δ(•, B) is a function from D to [0, p],
but not all functions from D to [0, p] can be expressed by a CA ballot; this can be seen by a simple counting argument. Therefore, assuming that a voter's preferences are representable by a CA ballot is a domain restriction.

4 Conditional minisum and minimax Definition 3. An n-voter conditional approval profile is a collection P = B 1 , . . . , B n of CA ballots. It is separable if every B i is a separable CA ballot. A CA irresolute rule is a function mapping any CA profile to a nonempty subset of D * .

Resolute rules can be obtained by combining an irresolute rule with a tie-breaking mechanism: if F is an irresolute rule and T is a tie-breaking mechanism induced by a fixed priority ranking over D * , then F T is the composition of F by T . (Note that T is a relation over an exponentially large set, and thus should be represented in some compact way; this has no impact in our results.) We now define two specific rules.

Conditional minisum

Definition 4. Given a CA profile P = B 1 , . . . , B n , the conditional minisum rule outputs the outcomes that minimize the total number of disagreements over all voters:

CondM iniSum(P ) = argmin d∈D * n i=1 δ(d, B i )
Example 2. Let D * = D = {m, f, v} × {r, w, b} and let P be the following 19-voter CA profile; the number on the top row denotes the number of voters; for instance, 5 voters express the CA ballot {{m, f }, m : r, f : w, v : {r, w}}. Note that the CA ballots of the second column are separable.

5 4 4 3 3 {m, f } m {m, f, v} {f, v} v m : r m : b m : {r, w} m : r m : b f : w f : b f : w f : w f : b v : {r, w} v : b v : {r, b} v : {r, w} v : {r, w, b}
The table below gives the disagreement values for each group of voters, and the total disagreement value. If P is separable, then the outcome of CondM iniSum(P ) can be obtained simply by decomposing the vote into a series of standard approval votes, variable by variable. Let C be a set of candidates; a (standard) approval profile is a collection of approval ballots (A 1 , . . . , A n ) where A i ⊆ C for each i; the (standard) approval rule App maps any approval profile to a winner (or a set of tied winners), and is defined by

mr mw mb f r f w f b vr vw vb 5 0 1 1 1 0 1 1 1 2 4 1 1 0 2 2 1 2 2 1 4 0 0 1 1 0 1 0 1 0 3 1 2 2 1 0 1 0 0 1 3 2 2 1 2 2 1 0 0 0 total
App(A 1 , . . . , A n ) = argmax c∈C |{i : c ∈ A i }|.
Observation 3. If P = B 1 , . . . , B n be a separable CA profile. Then CondM iniSum(P ) = App(P ↓X1 ) × . . . × App(P ↓Xp ), where P ↓Xi = (B ↓Xi 1 , . . . , B ↓Xi n ). In particular, if all variables are binary, we are in the context of a multiple referendum without preferential dependencies between issues, and where the decision is made issue by issue according to majority.

On the other hand, classical committee elections correspond to a collection of such separable conditional approval ballots, where binary variables correspond to candidates, together with a constraint on the cardinality of the committee.

Observation 4. Let D = {yes, no} p (plus possibly a cardinality constraint), and P = B 1 , . . . , B n a separable CA profile. Then CondM iniSum(P ) coincides with the output of simple multiwinner approval voting (also called minisum in [Brams et al., 2007]).

Conditional minimax

Conditional minisum has a utilitarianistic flavor: minimizing the sum of disagreements corresponds to maximizing social welfare. We may instead take an egalitarian point of view and minimize the maximum disagreement: Definition 5. Given a CA profile P = B 1 , . . . , B n , the conditional minimax rule outputs the outcomes that minimize the maximum number of disagreements over all voters :

CondM iniM ax(P ) = argmin d∈D * n max i=1 δ(d, B i )
Example 2, continued. The only outcome that has a disagreement at most 1 (in fact, exactly 1) with every conditional ballot is f b, therefore CondM iniM ax(P ) = {f b}: every agent is either happy with the fish, or unhappy with the fish, but in that case, happy with the beer given that the dish is fish.

A similar result as Observation 4 holds for conditional minimax: Observation 5. Under the same assumptions as in Observation 4, CondM iniM ax(P ) is the output of minimax approval voting [Brams et al., 2007]).

Computation

Computing a winning committee for minimax approval voting is NP-hard [START_REF] Moti | [END_REF], which of course carries over to conditional minimax approval voting. On the other hand, the polynomial time computation of winning committees for simple approval voting does not carry over to conditional minisum (even without feasibility constraints): Proposition 1. Given a CA profile P and an integer k, deciding whether there exists d ∈ D such that i δ(d, B i ) ≤ k is NP-complete, even for binary variables and, for each voter, an acyclic dependency graph with maximal indegree 1.

Proof. The problem is clearly in NP. The hardness proof is based on a reduction from MAX2SAT. Consider an instance I of MAX2SAT, defined by a set of variables V = {x 1 , . . . , x p }, a set of 2-clauses C = {C 1 , . . . , C n }, and an integer t. We create a set of binary variables X = {X 1 , . . . , X p } and, for each clause C i = l j ∨ l k , where l j ∈ {x j , x j } and l k ∈ {x k , x k }, and j < k, a voter i whose CA ballot B i is defined as follows (with the convention l j = x j when l j = x j ):

• the dependency graph has a single edge X j -→ X k .

• the CA statements are {x j , x j }, l j : {x k , x k }, l j : {l k }, and {x q , x q } for each q = k. First, assume that there exists an assignment of the variables V satisfying at least t clauses. From that assignment, we define an alternative d as follows: for all q ≤ p, d q = x q if x q is assigned to true, and d q = x q if x q is assigned to false. Then, we study the number of disagreements between B i and d. There are two cases to consider:

• If C i is satisfied, then either l j is true or l k is true, and from the definition of B i , we get δ(d, B i ) = 0.

• If C i is not satisfied, then both l j and l k are false, and from the definition of B i , we get δ(d, B i ) = 1.

Thus, since the assignment satisfies at least k clauses, we have

n i=1 δ(d, B i ) ≤ p -t.
Conversely, if there exists a d such that n i=1 δ(d, B i ) ≤ p -t, then by a similar line of reasoning as above, we can construct an assignment that satisfies at least t clauses.

On the positive side, there is a model-preserving translation from winner determination for conditional minisum approval voting into maximum satisfiability. We omit the formal details of the general construction but only give an example: the CA ballot B = {x, x : y, x : {y, y}} is translated into the set of clauses C B = {x, x ∨ y, x ∨ y ∨ y} (equivalent to {x, x ∨ y}); then a CA profile P = B 1 , . . . , B n is translated into a multi-set of clauses C P = C B1 ∪ . . . ∪ C Bn . For all d ∈ D, we have i δ(d, B i ) = α if and only if the number of clauses that are not satisfied in C P is equal to α. 1This means that conditional minisum approval voting can be solved by off-the-shelf MAXSAT solvers. Moreover, in the case of binary variables with acyclic dependency graphs, it is easy to check that this transformation preserves the differential approximation ratio of 4.34/(m + 4.34) obtained for MAXSAT in [START_REF] Escoffier | [END_REF], where m represents the number of clauses.2

Properties

We now study some of the properties satisfied by conditional minisum and conditional minimax. Because our rules have CA ballots as input, some of the standard properties have to be adapted. Clearly, our rules satisfy anonymity (the outcome is independent from the identity of voters), and their irresolute versions satisfy value neutrality (the outcome is unchanged after any renaming of the values of some variable) and variable neutrality (the outcome is unchanged after any renaming of the variables). Now we define four other important properties. For two profiles P = (P 1 , . . . , P n ) and Q = (Q 1 , . . . , Q m ), we note P +Q = (P 1 , . . . , P n , Q 1 , . . . , Q m ). Definition 6. Let F be a CA irresolute rule.

• F satisfies reinforcement (resp. weak reinforcement) if for any two profiles P, Q, if

F (P ) ∩ F (Q) = ∅ then F (P + Q) = F (P ) ∩ F (Q) (resp. ⊇ F (P ) ∩ F (Q)).
• F satisfies monotonicity if for any profile

P = (B 1 , . . . , B n ), if d = (x 1 , . . . , x p ) ∈ F (P ) and B i
is obtained from B i by adding x j to the set of values of X j approved conditioned on d P ar G i (Xj ) , then

d ∈ F (B 1 , . . . , B i-1 , B i , B i+1 , . . . , B n ).
The following properties, participation and strategyproofness, are defined only for resolute rules.3 Moreover, they require a voters' preference i to be induced from her CA ballot. We say that a voter with CA ballot B i has δ-induced preferences if for all x, y ∈ D * , x i y if and only if δ(x, B) ≤ δ(y, B). δ-induced preferences generalize Hamming-induced preferences [START_REF] Lainé | [END_REF]. Definition 7. Let F T the resolute version of some irresolute rule F , for some tie-breaking mechanism T . F T satisfies

• F T satisfies δ-participation if for any CA profile P = (B 1 , . . . , B n ) and CA ballot B n+1 , we have δ(F T (P + {B n+1 }), B n+1 ) ≤ δ(F T (P ), B n+1 ).
• F T satisfies δ-strategyproofness if for any CA profile

P = (B 1 , . . . , B n ) and CA ballot B i , δ(F T (B 1 ,. . .,B i-1 ,B i ,B i+1 ,. . .,B n ), B i ) ≥ δ(F T (P ), B i ).
We have the following results for conditional minisum and minimax. All proofs are simple (we omit them). Proposition 2.

• CondM inSum satisfies reinforcement.

• CondM inM ax does not satisfy reinforcement, but satisfies weak reinforcement.

• CondM inSum and CondM inM ax satisfy monotonicity.

• for any T , CondM inSum T and CondM inM ax T satisfy δ-participation.

Since minimax approval voting [START_REF] Caragiannis | Approximation algorithms and mechanism design for minimax approval voting[END_REF] is manipulable, this is a fortiori the case for CondM inM ax T . On the other hand, simple approval voting (or 'minisum') is strategyproof; however, as soon as preferences are conditional, strategyproofness is lost: Proposition 3. CondM inSum T is δ-manipulable, even for two binary variables.

Proof. Without loss of generality, assume T favors x 1 x 2 over x 1 x 2 . Consider the two-voter CA profile consisting of one ballot {x 1 , x 1 : x 2 , x 1 : x 2 } and one ballot {x 1 , x 1 : x 2 , x 1 : x 2 }. The outcome is x 1 x 2 . Now, if the second voter casts the ballot {x 1 , x 1 : x 2 , x 1 : x 2 } instead, then the outcome is x 1 x 2 , which is her preferred alternative under the assumption that her preferences are δ-induced.

Sequential conditional approval voting

For the sake of simplicity, in this section we assume that D * = D (handling feasibility constraints is possible but makes the definitions more complicated).

Sequential conditional approval voting (SCAV) is an approval-based version of sequential voting in combinatorial domains [START_REF] Lang | [END_REF]. Its intuitive principle is that decisions are taken variable by variable, following a fixed order; for each variable, the winning value is the one that has the maximal approval score given the values of the variables that come before it. The SCAV rule is defined only for CA profiles whose dependency graph is compatible with a given order O of the variables. In the rest of this section, O = X 1 > . . . > X p is a fixed order over X .

Definition 8 (O-legal conditional approval profile). A CA ballot is O-legal if its dependency graph contains no edge from X j to X i with i < j.

A CA profile is O-legal if it is composed of O-legal CA ballots.
Definition 9. The sequential conditional approval voting rule (SCAV) is the rule defined on O-legal CA profiles as follows. Let App Ti be the standard approval voting rule on D i , together with a tie-breaking mechanism T i , and T = (T 1 , . . . , T p ). For each O-legal P = (B 1 , . . . , B n ), let

P ↓Xi|x * 1 ,...,x * i-1 = (B ↓Xi|x * 1 ,...,x * i-1 1 , . . . , B ↓Xi|x * 1 ,...,x * i-1 n ); then SCAV T (P ) = (x * 1 , . . . , x * p ), where • x * 1 = App T1 (P ↓X1 ); • for each i = 2, . . . , p, x * i = App Ti (P ↓Xi|x * 1 ,...,x * i-1
). Similarly, we define the SCAV irresolute rule: SCAV (P ) is the set of all (x 1 , . . . , x p ) such that x 1 ∈ App(P ↓X1 ) and for all i = 2, . . . , p, x i ∈ App(P ↓Xi|x * 1 ,...,x * i-1 ). 45 A cheap protocol for SCAV T is composed of p rounds: at round i, it elicits only the voters' approval ballots

B ↓Xi|x * 1 ,...,x * i-1 j . Thus, the communication complexity of SCAV T is O(pn max i |D i |).
Example 2 (continued). Let O = X 1 > X 2 . The approval scores for m, f and v are respectively 13, 12 and 10: the (unique) selected value for X 1 is m. Given X 1 = m, the approval scores for r, w and b are respectively 12, 4 and 7: the selected value for X 2 is r. The outcome is mr.

For irresolute SCAV, some winning alternative can be computed in polynomial time, and determining whether a given alternative is a winner is also polynomial-time computable: at each round, it suffices to check that x i ∈ App(P ↓Xi|x * 1 ,...,x * i-1 ). Also, winner determination is polynomial for resolute SCAV (with a polynomial-time computable tie-breaking mechanism).

Recall that winner determination for conditional minisum is NP-hard, even for an O-legal profile. Moreover, both conditional minisum and SCAV coincide with simple approval voting when restricted to separable profiles. Thus, we may wonder how good an approximation SCAV is to conditional minisum. The approximation ratio is different whether we measure the quality of an alternative by its number of disagreements i δ(•, B i ) (to be minimized) or its number of agreements i (p -δ(•, B i )) (to be maximized):

4 Note that writing

B ↓X i |x * 1 ,...,x * i-1 j
is a slight abuse of notation, since O-legality implies only that P arG j (Xi) ⊆ {X1, . . . , Xi-1}.

Proposition 4. Consider a combinatorial domain with p variables X 1 , . . . , X p with |D i | = α i for all i.

• The largest possible ratio, over all O-legal CA profiles, between the disagreement score of the SCAV winner and the disagreement score of the conditional minisum winner, is

p i=1 (1 -1 αi )/(1 -1 α1 ).
• The largest possible ratio, over all O-legal CA profiles, between the agreement score of the conditional minisum winner and the agreement score of the SCAV winner, is

( 1 α1 + p -1)/ p i=1 1
αi . Proof sketch. Consider n CA ballots, with the dependency graph whose edges are {X i → X j | 1 ≤ i < j ≤ p}. First, we look for an upper bound of the disagreement score of an alternative that wins for SCAV. Assume that (x 1 1 x 2 1 . . . x p 1 ) is the winning alternative for SCAV. Then, for (x 1 1 x 2 1 . . . x p 1 ) to be winning for SCAV, x 1 1 , respectively x 1 1 : x 2 1 , . . . , x 1 1 x 2 1 . . . : x p 1 , has to be approved by at least n/α 1 voters, respectively n/α 2 , . . . , n/α p voters. It leads to a maximal disagreement score of n • p i=1 1 -1/α i . Now, we look for a lower bound for the disagreement score of a winning alternative for conditional minisum, under the condition that it is winning for SCAV. Assume that (x 1 2 x 2 2 . . . x p 2 ) is a winning alternative for conditional minisum. One way to avoid this alternative to be a SCAV winner is by eliminating it at the first round of the sequential vote. It implies that x 1 2 has to be approved by less than n α1 voters, and, in this case, x 1 2 : x 2 2 , . . . , x 1 2 x 2 2 . . . : x p 2 can be approved by any number of voters, and then, the disagreement score for (x 1 2 x 2 2 . . . x p 2 ) without being a winner for SCAV is at least n -n α1 . There exist other ways to avoid (x 1 2 x 2 2 . . . x p 2 ) to be a winning alternative for SCAV, but they lead to disagreement scores that are at least as large as that one. The proof for the bound relative to the agreement score is similar.

Thus, we get a worst-case ratio of ( p i=1 1 -1/α i )/(1 -1/α 1 ). This ratio is reached on a CA profile for which there is a perfect split between all values of X 1 (with each voter approving only one value), and then, conditionally on the value chosen for X 1 (by tie-breaking), there is also a perfect split for each other variable, while for some of the other values of X 1 , there is a perfect agreement between all voters on all other variables. The same profile shows that the bound relative to the number of agreements is reached too.

In particular, when all domains have the same cardinality α, then these ratios become p and (1 + (p -1)α)/p, and in the case of binary domains, p and 2 -1/p. Corollary 1. If the quality of a solution is measured by the number of disagreements, then SCAV is a

p i=1 (1 - 1 αi )/(1 -1 α1 )-approximation of CondM iniSum. If it is measured by the number of disagreements, then it is a ( 1 α1 + p -1)/ p i=1 1
αi -approximation of CondM iniSum. Proposition 5. SCAV satisfies anonymity, neutrality, reinforcement and monotonicity. For any tie-breaking mechanism T , SCAV T does not satisfy δ-participation.

Proof. Anonymity and neutrality are obvious. Reinforcement and monotonicity are easily proven by induction on the variables, and using the fact that standard approval voting satisfies them. For participation, consider the two-voter CA profile P = B 1 , B 2 with B 1 = {x 1 , x 2 , x 3 } and B 2 = {x 1 , x 2 , x 3 }. Assume without loss of generality that T favors x 1 x 2 x 3 over x 1 x 2 x 3 , thus F T (P ) = x 1 x 2 x 3 . Now, consider a third voter with preferences δ-induced by B 3 = {x 1 , x 1 : x 2 , x 1 : x 2 , x 1 : x 3 , x 1 : x 3 }. If she votes, the outcome is x 1 x 2 x 3 , which has only one agreement with her ballot, while x 1 x 2 x 3 had two. Proposition 6. For any tie-breaking mechanism T , SCAV T is not δ-strategyproof for p ≥ 3 (even for binary variables), and is strategyproof for two variables.

Proof. For p = 3, consider the profile in the proof of Proposition 5. If the third voter expresses his sincere ballot {x 1 , x 1 : x 2 , x 1 : x 2 , x 1 : x 3 , x 1 : x 3 }, the outcome is x 1 x 2 x 3 . If he expresses instead the ballot {x 1 , x 1 : x 2 , x 1 : x 2 , x 1 : x 3 , x 1 : x 3 } the the outcome is x 1 x 2 x 3 , which he prefers.

For p = 2, let X = {X 1 , X 2 }; suppose that voter i with sincere ballot B has a manipulation B ; let x * be the outcome for her sincere vote and x + the outcome after manipulation. Since B is a manipulation and i has δ-induced preferences, we have δ(x + , B ) < δ(x * , B). This implies (1) δ(x + , B ) < 2 and (2) δ(x * , B) > 0. Once the value of X 1 has been fixed, we are left with a standard approval vote on X 2 and i cannot obtain a better outcome than the sincere outcome; therefore, the manipulation has to bear on X 1 and has to change the value of X 1 , from, without loss of generality, x 1 to x 1 , with i initially approving x 1 but not x 1 . But this implies that δ(x + , B ) ≥ 1 and δ(x * , B) ≤ 1, which together with (1-2) implies δ(x + , B ) = δ(x * , B) = 1, therefore i, who has δ-induced preferences, cannot prefer x + to x * .

Conclusion

We have generalized approval voting to combinatorial domains and to nonseparable preferences. Our rules are natural generalizations of simple (or 'minisum') and minimax approval voting. They can be applied for all types of voting on combinatorial domains. Conditional minisum and minimax are computationally hard, but conditional minisum can be solved by MaxSAT solvers. Sequential conditional AV is easy to compute. Our rules (especially conditional minisum) satisfy a number of important properties.

There are other rules for multiwinner elections using approval ballots, where the satisfaction of a voter is not simply the Hamming distance to his preferred outcome (see [Kilgour, 2010] for a survey, and [Aziz et al., 2015b] for their computation). These rules could be adapted to CA ballots in the same way as we did for simple and minimax approval voting. The rules defined in [START_REF] Amanatidis | [END_REF], where the dissatisfaction degrees of the voters is aggregated by ordered weighted averages intermediate between max and sum, could also be generalized to CA ballots.

  The semantics of B and N (B) are different: if N (B) is the partial order over D induced from N (B), and B the order defined by d B d if δ(d, B) < δ(d , B), then B and N (B) are incomparable: Observation 2. There exists B and d, d such that d B d and d N (B) d.

  13 21 18 26 14 19 13 17 17 and we have CondM iniSum(P ) = {mr, vr}.

The translation works also for nonbinary variables, but is more complex, as there is one propositional variable per value di ∈ Di, and a feasibility constraint ensures that each variable takes exactly one value. More generally, when D * is a strict subset of D, the feasibility constraints are expressed as propositional formulas, and are in sufficiently many copies so as to be 'protected'.

Given an instance I of a combinatorial optimization problem, the differential approximation ratio measures "the relative position of the value of an approximated solution in the interval between the value of a worst feasible solution of I, and the value of a best solution of I"[START_REF] Escoffier | [END_REF].

For irresolute rules, we would require to first define a preference extension principle (see, e.g.,[START_REF] Brandt | [END_REF] for strategyproofness and[Brandl et al., 2015] for participation).

Note that the restriction of SCAV to separable profiles coincides with simple approval voting.

Acknowledgement

We thank Bernard Ries for helpful discussions. This research has been supported by the project ANR-14-CE24-0007-01 CoCoRICo-CoDec.


Nicolas Crouseilles

Shi Jin
email: sjin@wisc.edu

Mohammed Lemou

Liu Liu
email: lliu@math.wisc.edu

Jin ‡ Shi

Nonlinear Geometric Optics Based Multiscale Stochastic Galerkin Methods for Highly Oscillatory Transport Equations with Random Inputs

Keywords: highly oscillatory PDEs, nonlinear geometric optics, asymptotic preserving, uncertainty quantification, generalized polynomial chaos, stochastic Galerkin method, surface hopping

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Computational high frequency waves is challenging since one needs to numerically resolve the small wave length which is often prohibitively expensive [START_REF] Engquist | Computational high frequency wave propagation[END_REF]. Recently, a nonlinear geometric optics (NGO) based numerical method was introduced in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] for a class of highly oscillatory transport equations which allows the use of mesh size and time step independent of the wave length. The transport equations that were solved by this method are deterministic and are relevant to semiclassical approximations to quantum dynamics with band-crossings [START_REF] Chai | Semi-classical models for the Schrödinger equation with periodic potentials and band crossings[END_REF][START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF][START_REF] Morandi | Effective classical Liouville-like evolution equation for the quantum phasespace dynamics[END_REF][START_REF] Morandi | Wigner model for Klein tunneling in graphene[END_REF], a general non-adiabatic quantum mechanics phenomenon that can be found in a variety of important physics and chemistry problems, such as chemical reaction [START_REF] Tully | Trajectory surface hopping approach to nonadiabatic molecular collisions: the reaction of h + with d2[END_REF][START_REF] Tully | Molecular dynamics with electonic transitions[END_REF], Bose-Einstein condensation [START_REF] Wu | Nonlinear Landau-Zener tunneling[END_REF], and graphene [START_REF] Neto | The electronic properties of graphene[END_REF][START_REF] Fefferman | Honeycomb lattice potentials and Dirac points[END_REF].

Since these semiclassical models, and many other quantum models, use potential matrices that are obtained in empirical or ad hoc ways, they inevitably contain uncertainties. Uncertainty could arise in the potential, boundary or initial data, and forcing terms. In this paper, we are interested in the random uncertainty of the potential energy surface, and the goal is to develop efficient computational methods to compute the propagation of the uncertainty. To deal with uncertainty in general, the so-called gPC-SG methods (combination of generalized polynomial chaos (gPC) approximation with stochastic Galerkin (SG) projections) are known to be efficient for a wide range of partial differential equations with random uncertainties (see for example [START_REF] Ghanem | Stochastic finite elements: A spectral approach[END_REF][START_REF] Gunzburger | Stochastic finite element methods for partial differential equations with random input data[END_REF][START_REF] Le Maître | Spectral methods for uncertainty quantification, with applications to computational fluid dynamics[END_REF][START_REF] Xiu | Numerical Methods for Stochastic Computations[END_REF][START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]). However, direct application of such methods to highly-oscillatory problems with uncertainty becomes computationally expensive if one wants to correctly capture the effect of these oscillations, since one needs to resolve numerically the oscillations. In this paper, we are interested in extending the nonlinear geometric optics (NGO)-based method, which was developed in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] for a class of highly-oscillatory deterministic problems, to highly-oscillatory problems with uncertainty, in the framework of gPC-SG methods. In particular, we will develop here gPC-SG methods that not only allow the mesh size and time step, but also the order of the gPC approximation, to be independent of the small wave length.

We develop the gPC-SG methods for two model equations: a scalar model transport equation and a semiclassical surface hopping model developed in [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF]. Both models were studied in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] in the deterministic case. Here we assume random coefficients in these models that correspond to the band-gap in the non-adiabatic quantum dynamics, and the band gap could become small to model the so-called avoided crossing in which the quantum transition between bands is significant. With the random uncertainty, we first show that such a method is oscillatory for higher order derivatives in the random space, thus requires a frequency dependent discretization in the random space, which becomes prohibitively expensive to compute. We then modify this method by introducing a new "time" variable based on the phase, which is shown to be non-oscillatory in the random space, based on which we develop a gPC-SG method that can capture the pointwise solution with frequency-independent time step, mesh size as well as the degree of polynomial chaos. This method is then extended to the semiclassical model of surface hopping with random band gap, with the same numerical property.

The paper is organized as follows. In section 2, using a scalar equation with uncertain random coefficient, we develop the gPC-SG method by either solving the equation directly, or using the NGO approach. We prove that in both cases the method needs to use wave frequency dependent gPC order (although for the latter method such a dependence is one order more milder). Based on a theoretical result of [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], we introduce a new time variable using the phase, which will be shown to be non-oscillatory also in the random space, allowing us the develop a new NGO-based method which are capable of obtaining accurate pointwise solution without resolving the oscillations by any of the numerical parameters. Numerical examples will verify the theoretical property as well as the aforementioned numerical prop-erties. This method is then extended to a semiclassical surface hopping model in section 3, with the same numerical properties which are demonstrated numerically in section 4. We conclude the paper in section 5.

The one dimensional scalar equation with random inputs

In this section we consider a one dimensional problem with random inputs. Here the unknown is u(t, x, z) ∈ C, where x ∈ Ω ⊂ R, t ≥ 0, and z ∈ Iz ⊂ R n , n ≥ 1, is the random variable with a prescribed probability density function π(z). We consider the following scalar model for u:

∂tu + c(x)∂xu + r(u) = ia(x, z) ε u, u(0, x, z) = uin(x, z), (2.1)
where a(x, z), c(x), r(u) and uin(x, z) are all given functions, with a ≥ 0. We assume that Ω is a bounded interval and that periodic boundary conditions are considered for the space variable x.

A gPC-SG framework

We briefly describe the SG method. Let P n P be the space of the n-variate polynomials of degree less than or equal to P , P ≥ 1, and recall that

dim(P n P ) = card{k ∈ N n , |k| ≤ P } = n + P P := K,
where we have denoted k = (k1, . . . , kn) and |k| = k1 + • • • + kn. We consider the following inner product

f, g π = Iz f (z)g(z)π(z)dz, ∀f, g ∈ L 2 (π(z)dz),
where L 2 (π(z)dz) is the usual weighted Lebesgue space, and its associated norm is

f 2 L 2 (π(z)dz) = Iz |f (z)| 2 π(z)dz.
Consider a corresponding orthonormal basis {ψ k (z)} k∈N n ,|k|≤P of the space P n P , where the degree of ψ k is deg(ψ k ) = |k|. In particular

ψ k , ψ l π = Iz ψ k (z)ψ l (z)π(z)dz = δ kl , |k|, |l| ≤ P,
where δ kl is the Kronecker symbol. The commonly used pairs of {ψ k (z)} and π(z) include Hermite-Gaussian, Legendre-uniform, Laguerre-Gamma, etc. Since the family {ψ k (z)} |k|≤P has K elements we introduce its renumbering family ( ψ1, ψ2,

• • • ψK ), that is {ψ k (z)} k∈N n ,|k|≤P = ( ψ1, ψ2, • • • ψK ), deg( ψj) ≤ deg( ψj+1).
The SG method consists in seeking the solution to (2.1) as a projection onto the space P n P , that is

u(t, x, z) ≈ |k|≤P u k (t, x)ψ k (z) = K j=1 ũj(t, x) ψj(z), (2.2)
with

u k = u, ψ k π , ũj = u, ψj π , |k| ≤ P, j = 1, 2, ..., K.
From this approximation one can easily compute statistical moments, such as the mean and standard deviation, as

E(u) ≈ u (0,••• ,0) = ũ1, SD(u) ≈   1≤|k|≤P |u k | 2   1/2 = K j=2 |ũj| 2 1/2 . (2.3)

The direct gPC-SG method

We first introduce the gPC-SG method for the direct method (hereafter called gPC-SG-D)-solving (2.1) directly. For simplicity of illustration, assume c(x) does not depend on z. The case when c depends on z can be easily incorporated into the gPC-SG framework [START_REF] Gottlieb | Galerkin method for wave equations with uncertain coefficients[END_REF]. The gPC-SG solution for u is computed through the projection formula (2.2) as follows.

Denote the gPC coefficients

#» u (t, x) = (ũ1(t, x) • • • , ũK (t, x)) T .
Inserting the approximated solution (2.2) to (2.1) and conducting the standard Galerkin projection, one gets

∂t #» u + c(x)∂x #» u + #» γ (u) = i ε A(x) #» u , #» u (0, x) = Iz uin(x, z) ψ(z)π(z)dz , (2.4)
with

ψ(z) = ( ψ1(z), ψ2(z), • • • , ψK (z)) T , u(t, x, z) = K j=1 ũj(t, x) ψj(z). (2.5)
The j-th component of the vector

#» γ = (γ1, γ2, • • • , γK ) is given by γj(u) = Iz r(u) ψj(z)π(z)dz , for j = 1, • • • , K, (2.6)
and the symmetric, non-negative definite matrix A ∈ RK×K (recall that a ≥ 0) is defined by

Asj(x) = Iz a(x, z) ψs(z) ψj(z)π(z)dz , for s, j = 1, • • • , K. (2.7)
To solve numerically (2.4), we use a simple time splitting. To this aim, we fix a time step ∆t > 0 and set tn = n∆t for n ∈ N. As usual, we denote by

#» u n (x) = (ũ n 1 (x), ũn 2 (x), • • • , ũn K (x)
) an approximation of #» u (tn, x). We split system (2.4) in three parts which can be solved as follows.

• Oscillatory part

The oscillatory part

∂t #» u = i ε A(x) #» u ,
is solved exactly in time, for each x, to get

#» u n+1 (x) = exp i ε A(x)∆t #» u n (x) .

• Nonlinear part

To solve the nonlinear part ∂t #» u + #» γ (u) = 0, we choose to use the forward Euler method, that is

#» u n+1 = #» u n -∆t #» γ (u n ) . (2.8) To compute #» γ (u n ), we first compute u n (x, z) = K s=1
ũn s (x) ψs(z) using the Gauss-quadrature rule (where the quadrature points are chosen as the roots of the orthogonal polynomials determined by the distribution π of the random variables, see [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]) to get #» γ (u) given by (2.6) as follows. Define the Gauss-quadrature points as z (l) and the corresponding weights ω l , l = 1, • • • , Ng. The value of r(u) evaluated at x, t and at the Gauss-quadrature point z (l) is given by r(u(z

(l) )) = r K j=1 ũj(t, x) ψj(z (l) ) . By (2.6), one has γj(u) = Ng l=1 r(u(z (l) )) ψj(z (l) )ω l .
which is used in (2.8) to get #» u n+1 .

• Transport part

To solve the transport part

∂t #» u + c(x)∂x #» u = 0,
we use a pseudo-spectral method in space and a three-stage Runge-Kutta method in time which was introduced in [START_REF] Choi | The reconstruction of upwind fluxes for conservation laws: Its behavior in dynamic and steady state calculations[END_REF]. Denote T (u n ) = -c(x)F -1 (iξF(u n )), where F and F -1 are the (discrete) Fourier and inverse Fourier Transforms respectively, and ξ is the Fourier space variable. Then ) ).

#» u n,(1) = #» u n + 1 2 ∆tT ( #» u n ), #» u n,(2) = #» u n + 1 2 ∆tT ( #» u n,(1) ), (2.9)
#» u n+1 = #» u n + ∆tT ( #» u n,( 2
Notice that the spectral approximation is central type finite difference approximation with purely imaginary spectrum, thus one needs to use ODE solvers that have a stability region that contains part of the purely imaginary axis [START_REF] Choi | The reconstruction of upwind fluxes for conservation laws: Its behavior in dynamic and steady state calculations[END_REF]. The scheme (2.9) is such an ODE solver where the stability region takes the largest part of the imaginary axis among three-stage ODE solvers, and is of second order accuracy in time.

The NGO-based gPC-SG method

We first review the NGO based method introduced in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] for the deterministic onedimensional scalar equation,

∂tu + c(x)∂xu + r(u) = ia(x) ε u, u(0, x) = uin(x), (2.10)
where the functions uin, a, c and r are given and a(x) ≥ 0. Periodic boundary conditions in space are considered. For the sake of simplicity, we assume that the initial data uin does not depend on ε, namely, it is non-oscillatory. We focus on the case when r is nonlinear. We recall how nonlinear geometric optics (NGO) is utilized to solve problem (2.10) in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF]. Introduce a profile function U (t, x, τ ) which depends on the 2π-periodic variable τ , and a phase S(t, x) such that

U (t, x, S(t, x)/ε) = u(t, x), (2.11)
with u solving (2.10). Inserting this ansatz into (2.10) gives

∂tU + c(x)∂xU + 1 ε [∂tS + c(x)∂xS] ∂τ U + r(U ) = ia(x) ε U.
Due to the periodicity constraint on U , the following equation on the phase S should be imposed (see [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] for details), ∂tS + c(x)∂xS = a(x), S(0, x) = 0, (2.12)

then the equation for U is given by

∂tU + c(x)∂xU + r(U ) = - a(x) ε (∂τ U -iU ), U (0, x, 0) = uin(x). (2.13)
For convenience, we write (2.13) in terms of V = e -iτ U and get

∂tV + c(x)∂xV + e -iτ r(e iτ V ) = - a(x) ε ∂τ V, V (0, x, 0) = uin(x). (2.14)
One needs initial data V (0, x, τ ) for all τ to solve equation (2.14). The only requirement we have to ensure is V (0, x, 0) = uin(x). One critical idea is to use initial data for V so that the solution to (2.14) remains bounded uniformly in ε up to certain order of derivatives.

Introduce the operators L and Π,

Lg = ∂τ g, Πg = 1 2π 2π 0 g(τ )dτ, and
L -1 g = (I -Π) τ 0 g(σ)dσ = τ 0 g(σ)dσ + 1 2π 2π 0 σg(σ)dσ.
Following the work of [START_REF] Chartier | Uniformly accurate numerical schemes for highly-oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly-oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], a Chapman-Enskog expansion is used to give a suitable initial condition given by

V (0, x, τ ) = uin(x) + ε a(x) [G(0, uin) -G(τ, uin)] , (2.15)
with G(τ, uin) = L -1 (I -Π) e -iτ r(e iτ uin) .

(2.16)

It was proved in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] that the solution V to (2.14) and (2.15) has bounded (uniformly in ε) derivatives in both x and t up to second order. Thus it allows the construction of a scheme with a uniform accuracy with respect to ε. We recall in the sequel the scheme introduced in To solve the equation (2.12) for S, the pseudo-spectral method in space and 4-th order Runge-Kutta method in time are used. In other words, the following ODE system on Sj(t) ≈ S(t, j∆x):

∂tSj + c(xj)F -1 (iξF(S n ))j = a(xj), Sj(0) = 0, (2.17)
is solved by the 4-th order Runge-Kutta method. To recover the original solution u at the final time tN , we use (2.11) at t = tN and x = xj. Since S N j /ε does not necessarily coincide with a grid point τ l at the final time tN = N ∆t, thus one can perform a trigonometric interpolation u(tN , xj) ≈ V N j (τ = S N j /ε).

(2.18)

Note that a higher order method is necessary here to solve S since we need to construct the quantity S(t, x)/ε, and the error of S is divided by ε.

To solve (2.14) for V , we start with V 0 j given by (2.15) at x = xj. Then, the scheme for V reads (assume c(x) > 0)

V n+1 j -V n j ∆t + c(xj) V n j -V n j-1 ∆x + e -iτ r(e iτ V n j ) = - a(xj) ε ∂τ V n+1 j . (2.19)
Here the right hand side is discretized implicitly due to its numerical stiffness and the τ variable is discretized using the Fourier transform.

The gPC-SG-N1 method

Our aim now is to extend the NGO method presented in the previous section to solve a highly oscillatory problem with uncertainty, namely (2.1). By [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], and similarly to (2.12) and (2.14), the NGO-based method solves S from ∂tS + c(x)∂xS = a(x, z), S(0, x, z) = 0, (2.20) and solves V from ∂tV + c(x)∂xV + e -iτ r(e iτ V ) = -a(x, z)

ε ∂τ V, (2.21)
with initial data

V (0, x, τ, z) = uin(x, z) + ε a(x, z) [G(0, uin, z) -G(τ, uin, z)] (2.22)
with G(τ, uin, z) = L -1 (I -Π) e -iτ r(e iτ uin) ,

where uin(x, z) and a(x, z) now depend on z. Below we detail the gPC approach for this system.

The gPC formulation for S

One inserts the gPC-SG ansatz

S(t, x, z) ≈ |k|≤P S k (t, x)ψ k (z) = K j=1 Sj(t, x) ψj(z)
into (2.20) and conducts the Galerkin projection, to get

∂t #» S + c(x)∂x #» S = #» R(x), #» S (0, x) = 0, (2.24)
where the gPC coefficients #» S is defined by

#» S (t, x) = ( S1(t, x), • • • , SK (t, x)) T ,
and

#» R(x) = Iz a(x, z) ψ(z)π(z)dz. (2.25)
We recall the notation: ψ = ( ψ1, ψ2, • • • , ψK ).

The gPC formulation for V

We insert the gPC-SG ansatz

V (t, x, τ, z) ≈ |k|≤P V k (t, x, τ )ψ k (z) = K j=1 Ṽj(t, x, τ ) ψj(z)
into (2.21) and conduct the Galerkin projection to get

∂t #» V + c(x)∂x #» V + e -iτ #» γ (e iτ V ) = - 1 ε A(x)∂τ #» V , (2.26)
where A is the matrix given by (2.7) and

#» V (t, x, τ ) = ( Ṽ1(t, x, τ ), • • • , ṼK (t, x, τ )) T .
The j-th component γj of the vector

#» γ is γj(e iτ V ) = Iz r(e iτ V ) ψj(z)π(z)dz, for j = 1, • • • , K, (2.27)
which is computed by the Gauss-quadrature formula.

The initial data

Here we consider that a(x, z) depends on the random variable z. The gPC approximation of the initial condition (2.22) becomes

#» V (0, x, τ ) = T + εY , (2.28)
where the vectors T ∈ R K and Y ∈ R K are given by

T (x) = Iz uin(x, z) ψ(z)π(z)dz, Y (x, τ ) = Iz 1 a(x, z) [G(0, uin, z) -G(τ, uin, z)] ψ(z)π(z)dz,
with ψ(z) given in (2.5).

The fully discrete gPC-SG-N1 method

Now we discretize (2.20), (2.21) in time as in the deterministic case. The phase vector #» S is obtained by applying a pseudo-spectral method in space, which leads to the following semi-discretized equation

∂t #» S j + c(xj)F -1 (iξF( #» S ))j = #» R(xj), #» S j (0) = 0 , ∀j = 0, • • • Nx, (2.29)
where #» R(x) is given by (2.25). Here, #» S j (t) is an approximation of #» S (t, xj) and the same notation is used for other quantities. Then a 4-th order Runge-Kutta method is applied to march in time.

To solve the equation on #» V , knowing #» S , we start by the initial value #» V 0 j given by (2.28) at grid points x = xj. Then the system is advanced in time by a simple time-splitting algorithm as explained for the deterministic case in section 2.2. The nonlinear part (the #» γ term) and the transport part are treated in the same way as in section 2.2. For the oscillatory part, we use the backward Euler method in time:

# » Vj n+1 - # » Vj n ∆t = - 1 ε A(xj)∂τ # » Vj n+1 , (2.30)
where

# » Vj n (τ ) is an approximation of #» V (n∆t, xj, τ
) and A is given by (2.7). Let #» V n j (ζ) be the Fourier transform of # » Vj n (τ ) in the periodic variable τ , where ζ is the Fourier variable. The spectral method is used to discretize the τ -derivative, then (2.30) becomes (removing the dependency in ζ)

# » Vj n+1 = I + iζ ∆t ε A(xj) -1 # » Vj n . (2.31)
The updated values #» V n+1 j are then obtained by the inverse Fourier transform. Note that since A is non-negative definite, all its eigenvalues are real and non-negative. Thus the matrix I + iζ ∆t ε A has non-zero eigenvalues, thus is invertible. In the sequel, we refer this method as the gPC-SG-N1 method.

The relation between the gPC order K and ε

The main advantage of the NGO-based method over the direct method is that the former allows one to use ∆t and ∆x independent of ε, thus when ε is small, one can still use relatively larger ∆t and ∆x to get accurate solution at the mesh points. In this section, we prove that this remarkable property is not true for the random approximation, namely one can not use the gPC order K independent of ε either for the gPC-SG-D or gPC-SG-N1.

• The gPC-SG-D

To show the gPC-SG-D method does not satisfy the uniform accuracy property, we consider the simple case where the functions c and a are independent of x, and where z ∈ R:

∂tu + c∂xu = ia(z) ε u .
By the method of characteristics, the analytic solution is given by

u(t, x, z) = e ia(z) ε t uin(x -ct, z),
where uin is the initial data. Taking the q-th order derivative in z, one sees that

∂ q z u = O(1/ε q ).
Assume uin ∈ H q z (π(z)dz), where H q z is the Sobolev space in z defined as

||u|| H q z = q l=0 ||∂ l z u|| L 2 (π(z)dz) 1 2
.

By the standard approximation theory of orthogonal polynomials (see for example [START_REF] Shen | Spectral methods: Algorithms, analysis and applications[END_REF]), one has

||u - K j=1 ũj ψj||L2(π(z)dz) ≤ C ε q K q ,
where C is a constant depending on t, the derivatives of a(z) and uin with respect to z. This shows that one needs to take K inversely proportional to ε in order to ensure an accurate approximation in z.

• The gPC-SG-N1

We consider the gPC-SG-N1 method and show that the degree of the polynomial approximation needs to be chosen dependent on ε to ensure a spectral accuracy. To this aim, we analyze the projection error on the function V , solution to (2.21),

∂tV + c(x)∂xV + e -iτ r(e iτ V ) = - a(x, z) ε ∂τ V . (2.32)
From (2.32), clearly, the l-th derivative in z of V satisfies basically the same equation as the l-th derivative in x of V (except one involves the derivatives of a in z while the other involves the derivatives of a in x). Therefore, as long as a and uin have the same regularity in x and z, V will have the same regularity in x and z. Since with the choice of initial data (2.22), we only gaurantee bounded (in terms of ε) x derivatives of V up to the second order, therefore, we only have up to second bounded derivatives in z for V . This is clearly not enough for SG methods (nor for the stochastic collocation method, to be introduced later) which can achieve high order-up to spectral-accuracy if there is sufficient regularity in the solution.

A new NGO-based gPC-SG method (gPC-SG-N2)

In the following, we show how to construct a numerical method for which the derivatives of V with respect to z and x are bounded uniformly in ε at any order. To this aim, let us consider the following change of function

V (t, x, z, τ ) = W (S(t, x, z), x, z, τ ), (2.33)
where S solves (2.20). Then, W is the solution to the following problem:

∂sW + c a(x, z) ∂xW + 1 a(x, z) e -iτ r(e iτ W ) = - 1 ε ∂τ W, W (0, x, z, τ ) = V (0, x, z, τ ) = uin(x, z) + ε a(x, z) G(0, uin, z) -G(τ, uin, z) , with G(τ, uin, z) = L -1 (I -Π)[e -iτ r(e iτ uin(x, z))] .
(2.34)

We then have the following result for W .

Proposition 2.1 Let W be the solution of (2.34) on [0, T ], T > 0, with periodic boundary condition in x and τ . Then, up to the second order derivative in s, and arbitrary order in x and z derivatives of

W are bounded uniformly in ε ∈ [0, 1], that is, ∃C > 0 independent of ε such that, ∀s ∈ [0, T ] ∂ p s ∂ q x ∂ r z W (s) L ∞ τ,x (L 2 (π(z)dz)) ≤ C, for p = 0, 1, 2, and q, r ∈ N.
This proposition can be deduced from Lemma 2.6 of [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] in which we replace x by (x, z). In fact, the W -equation (2.34) was introduced in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] to prove the uniform in ε regularity and convergence of the NGO-based method. Note in particular that the high order derivatives in x and z are uniformly bounded at any order with respect to ε although the initial data for W is only prepared at the first order in ε. This guarantees the uniform accuracy in z and x, that is in particular

||W - K j=1 Wj ψj||L2(π(z)dz) ≤ C K q , ∀q positive integer,
where C is a positive constant which is independent of ε, and q depends only on the regularity in z of the initial data uin(x, z). Thus one can choose K independent of ε.

To summarize, to obtain a spectral accuracy in z at any order uniformly in ε, one has to solve equation (2.34) satisfied by W (instead of the equation on V ) and then recover V by formula (2.33). Note that the mapping t → S(t, x, z) is a increasing function in t for any fixed (x, z). Of course, this needs an interpolation in time to get back to the t variable for V .

The gPC-SG-N2 method

The above discussion motivates us to design the following new scheme, denoted by gPC-SG-N2 in the sequel.

The gPC system for solving the equation for S is the same as in section 2.3.1. For the equation on W , we insert the gPC-SG ansatz

W (s, x, τ, z) ≈ |k|≤P W k (s, x, τ )ψ k (z) = K j=1 Wj(s, x, τ ) ψj(z)
into (2.34) and conduct the Galerkin projection to get

∂s # » W + c(x)A (x)∂x # » W + e -iτ # » γ (e iτ W ) = - 1 ε ∂τ # » W , (2.35)
where A is the symmetric and positive-definite matrix given by

A ij (x) = Iz 1 a(x, z) ψi(z)ψj(z)π(z)dz, and # » W (s, x, τ ) = ( W1(s, x, τ ), • • • , WK (s, x, τ )) T .
The j-th component γ j of the vector

# » γ is , γ j (e iτ W ) = Iz 1 a(x, z) r(e iτ W ) ψj(z)π(z)dz, for j = 1, • • • , K, (2.36)
which is computed by the Gauss-quadrature formula. The initial data for W is the same as for V : W (0, x, τ, z) = V (0, x, τ, z), which is shown in section 2.3.1.

The fully discrete gPC-SG-N2 method

Let the final time be T . Define the grid points in space xj, j = 0, 1, • • • , Nx. As discussed in (2.18), to recover the solution u, one needs a trigonometric interpolation on

τ = S(T, xj, z (l) ) ε , for each quadrature points z (l) , l = 1, • • • , Ns, u(T, xj, z (l) ) ≈ V (T, xj, τ = S(T, xj, z (l) ) ε , z (l) ),
then one can compute the mean and standard deviation of u at (T, xj) in the random space,

E(u)(T, xj) = Ns l=1 u(T, xj, z (l) )ω l , SD(u)(T, xj) = Ns l=1 |u(T, xj, z (l) )| 2 ω l -(E(u)(T, xj)) 2 1/2
, where ω l are the corresponding quadrature weights. After obtaining the values of S at (T, xj, z (l) ) for all j, l, let S = max j,l S(T, xj, z (l) ) .

We choose s ∈ [0, S ] to be the time domain for computing W . Define the time discretization of W as s l , for l = 0, 1, • • • , L, with time step ∆s and the total number of time steps

L = S ∆s + 1.
To solve the equation on # » W , knowing #» S , we start by the initial value # » W (0, x, τ ) = #» V (0, x, τ ) given by (2.28) at x = xj. Then the system is advanced in time by a simple time-splitting algorithm as explained for the deterministic case in section 2.2. The nonlinear part (the #» γ term) and the transport part are treated in a similar way as in section 2.2. For the oscillatory part, we use the backward Euler method in time,

# » W n+1 j - # » W n j ∆s = - 1 ε ∂τ # » W n+1 j , (2.37)
where

# » W n j is an approximation of # » W (n∆s, xj). Let # » Ŵ n j (ζ) be the Fourier transform of # »
Wj n in the periodic variable τ , where ζ is the Fourier variable. The spectral method is used to discretize the τ -derivative, then (2.37) becomes

# » Ŵj n+1 = I + iζ ∆s ε I -1 # » Ŵj n .
Note that the matrix I + iζ ∆s ε I has non-zero eigenvalues, thus is invertible.

The interpolation step

To find the values of

V at time T , since V (T, xj, τ k , z (l) ) = W (S(T, xj, z (l) ), xj, τ k , z (l) ), one uses linear interpolation to find W at (S(T, xj, z (l) ), xj, τ k , z (l) ), for l = 1, • • • , Ns.
We search for the interval [sn l , sn l +1] that contains S(T, xj, z (l) ). Denote S(T, xj, z (l) ) = s † , sn l = s (1) , sn l +1 = s (2) , W (sn l ) = W (1) , W (sn l +1) = W (2) , then by the Lagrange interpolation formula,

W (s † ) = W (1) s † -s (2)
s (1) -s (2) + W (2) s † -s (1) s (2) -s (1) ,

which gives V (T, xj, τ k , z (l) ) = W (s † , xj, τ k , z (l) ).
Remark 2.2 Compare to gPC-SG-N1, the transport step is more costly since one has to multiply by the matrix A . If a = O(1), then the CFL condition for the transport equation requires just ∆t = O(∆x) thus independent of ε. However, for avoided crossing, typically a(x, z) = O( √ ε) for some region of x, thus the eigenvalue of A are of O(1/ √ ε), thus an explicit time discretization of the transport step will require ∆t = O( √ ε∆x), further increasing the computational cost. One may use implicit time discretization in the transport step. This will be studied in our future research.

Numerical tests for the 1D scalar equation

For numerical comparison, we also use the stochastic collocation (SC) method [START_REF] Gunzburger | Stochastic finite element methods for partial differential equations with random input data[END_REF][START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF]. Let {z (j) } Nc j=1 ⊂ Iz be the set of collocation nodes, Nc the number of samples. For each fixed individual sample z (j) , j = 1, . . . , Nc, one applies the deterministic solver to the deterministic equations as in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], obtains the solution ensemble uj(t, x) = u(t, x, z (j) ), then adopts an interpolation approach to construct a gPC approximation, such as

u(t, x, z) = Nc j=1 uj(t, x)lj(z),
where lj(z) depends on the construction method. The Lagrange interpolation method is used here by choosing lj(z (i) ) = δij. Depending on whether the direct or NGO-based method is used for the deterministic method, we will have the gPC-SC-D and gPC-SG-N1, gPC-SG-N2 methods.

Example 2.1 We consider the numerical example in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], with a involving a 1D random variable z that follows a uniform distribution,

r(u) = u 2 /(u 2 + 2|u| 2 ), c(x) = cos 2 (x), a(x, z) = (3/2 + cos(2x))(1 + 0.5z) > 0,
and the non-oscillatory initial data given by

u(0, x) = 1 + 1 2 cos(2x) + i 1 + 1 2 sin(2x) , x ∈ [-π/2, π/2].
In the following tests, for all the reference solutions obtained from gPC-SC-D method, Nc = 64, Ng = 32 quadrature points are used. In gPC-SG-N1 and gPC-SG-N2 methods, Ns stands for the number of quadrature points used in the final reconstruction step (to get u from V or W ).

Figures 1 and2 demonstrate how ε affects the choices of K for gPC-SG-D. Here the gPC-SC-D, computed with very small ∆t and ∆x, is used to obtain the reference solutions. For gPC-SG-D, we also use very small ∆x and ∆t, in order to concentrate on effect of ε on K. We output the solution at t = 0.1. In Figure 1, for ε = 1 and 0.1, one can use K = 4. However, for ε = 0.01, neither K = 4 nor K = 8 gives the correct solutions, especially around the center of the domain. One starts to see satisfactory numerical solutions when K = 10. For ε = 0.005 one needs to use K = 20, see Figure 2.

Figures 3 and4 compare mean of the solutions of gPC-SG-N1 and gPC-SG-N2 for ε = 0.01, and ε = 0.005 respectively. We take small ε and use gPC-SC-D to fully resolve the oscillations (to serve as a reference solution). One can observe from Figure 4 that gPC-SG-N2 is able to capture correctly the solution (its mean and standard deviation) at mesh points with ∆x, ∆t much larger than ε and K = 4. For the same grids, gPC-SG-N1 is a bit off near x = 0. If one uses a much finer spatial size ∆x = π/1000, gPC-SG-N1 becomes more accurate but still is a bit off near x = 0, as shown in the second row of Figure 3.

Some comparisons between the gPC-SG-N1 and gPC-SG-N2 methods are shown in Figure 5 in the case of ε = 3 × 10 -3 . We first plot the mean of real part of u (zoomed in solutions) with K = 4, Ns = 144 in gPC-SG-N1 and gPC-SG-N2 and observes that gPC-SG-N2 can capture the correct mean and standard deviation, while gPC-SG-N1 can not. Using K = 4, Ns = 800 in gPC-SG-N1 still does not give the good result, and further increasing the gPC order (K = 16) enables the mean and standard deviation of Re(u) to capture the oscillations. Thus gPC-SG-N2 does not require larger K or Ns when ε is small, whereas larger K and Ns are needed in gPC-SG-N1.

A semiclassical surface hopping model with random inputs

We consider a semiclassical surface hopping model in one dimension in space and momentum, with random inputs,

         ∂tf + + p∂xf + -∂x(U + E)∂pf + = bi f i + b i f i , ∂tf -+ p∂xf --∂x(U -E)∂pf -= -bi f i -b i f i , ∂tf i + p∂xf i + ∂xU ∂pf i = -i 2E ε f i + b i (f --f + ) + (b + -b -)f i , (3.1)
where (f

+ (t, x, p, z), f -(t, x, p, z), f i (t, x, p, z)) ∈ R+ × R+ × C, (t, x, p) ∈ R+ × R × R, and b ± ∈ C, b i ∈ C, E = E(x, z) ∈ R are given functions depending on (x, p, z).
This system is equipped with initial conditions f ± (0, x, p) = f ± in (x, p) and f i (0, x, p) = f i in (x, p).

We also define the zeroth moments of f ±,i as

ρ ±,i (t, x, z) = R f ±,i (t, x, p, z)dp .
Here ρ ± denote the densities of particles in the two bands. This model, introduced in [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF], was a semiclassical approximation to the nucleonic Schrödinger system that arises from the Born-Oppenheimer approximation with non-adiabatic phenomenon, in which particles can "hop" from one potential energy surface to the other ones. Here f ± stand for the particle density distributions in the two energy surfaces respectively, while f i is the off-diagonal entry of the semiclassical Wigner matrix. The right hand side of the system describes the interband transition between different potential energy surfaces with the energy gap ∆E = 2E , which can be random in reality. In particular, the so-called avoided crossing corresponds to where the minimum energy gap of order √ ε.

The direct method

In this part and the following sections, we assume that the function E(x, z) is given and random, whereas no randomness is introduced to the initial function.

Introduce the following matrix

A = diag(-∂x(U + E), ∂x(U -E), ∂xU, ∂xU ), (3.2)
The specific form of b i depends on the potential matrix in the nucleonic Schrödinger equation.

Here we consider b i = b i (x, p) ∈ R, b ± = 0, which correspond to two specific potential matrices considered in [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF]. Let

B =      0 0 2b i 0 0 0 -2b i 0 -b i b i 0 2E/ε 0 0 -2E/ε 0     
, and write f = (f + , f -, g := Re(f i ), h := Im(f i )) ∈ R 4 , then (3.1) can be written as

∂tf + p∂xf + A ∂pf = Bf, (3.3)
which will be solved using simple operator splitting method:

1. Solve ∂tf + p∂xf = 0 using spectral method in space and exact integration in time,

2.

Solve ∂tf + A ∂pf = 0 using spectral method in velocity and exact integration in time,

3.

Solve ∂tf = Bf to be specified below.

The first two steps are the same for both the direct and the NGO-based method (to be introduced in the next section 3.2). For step 3, since the matrix B now involves randomness, instead of solving the whole system exactly in time, we use a Crank-Nicolson method in time coupled with a substitution method to save computations. More precisely, we have

f +,n+1 -f +,n ∆t = 1 2 2b i g n + g n+1 , (3.4) f -,n+1 -f -,n ∆t = - 1 2 2b i g n + g n+1 , (3.5) g n+1 -g n ∆t = 1 2 -b i f +,n + b i f -,n + 2E ε h n -b i f +,n+1 + b i f -,n+1 + 2E ε h n+1 , (3.6) h n+1 -h n ∆t = - 1 2 2E ε g n + g n+1 . (3.7)
We used classical notations for the time approximation, for example f +,n ≈ f + (tn). From (3.4), one solves f +,n+1 in terms of g n+1 , and then substitute it into (3.6) to get

g n+1 , 1 + (b i ) 2 (∆t) 2 + E 2 ε 2 (∆t) 2 g n+1 = g n + ∆t -b i f -,n + b i f -,n + 2E ε h n -(b i ) 2 ∆tg n - E 2 ε 2 ∆tg n . ( 3.8)
Then the other unknowns f ±,n+1 and h n+1 can be obtained from (3.4), (3.6) and (3.7).

The NGO based method

We first review the NGO-based method introduced in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], wherein deterministic E = E(x) and deterministic initial data is considered.

We introduce the augmented unknowns (F ± , G, H)(t, x, p, τ ) satisfying f ± (t, x, p) = F ± (t, x, p, S(t, x, p)/ε), f i (t, x, p) = e iS(t,x,p)/ε (G + iH)(t, x, p, S(t, x, p)/ε).

(3.9) The phase function S(t, x, p) is designed to follow the main oscillations in the model. Assume periodicity in S/ε, S solves

∂tS + p∂xS + ∂xU ∂pS = 2E, S(0, x, p) = 0. (3.10) Then F ± , G and H solve                      ∂tF + + p∂xF + -∂x(U + E)∂pF + = - E + ε ∂τ F + + 2b i (G cos τ + H sin τ ), ∂tF -+ p∂xF --∂x(U -E)∂pF -= - E - ε ∂τ F --2b i (G cos τ + H sin τ ), ∂tG + p∂xG + ∂xU ∂pG = - 2E ε ∂τ G + b i (F --F + ) cos τ, ∂tH + p∂xH + ∂xU ∂pH = - 2E ε ∂τ H + b i (F --F + ) sin τ, (3.11)
where we denoted

E ± (t, x, p) = 2E -∂x(2U ± E)∂pS . (3.12)
The well-prepared initial conditions are given by (see [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF])

F + (0, x, p, τ ) = f + in -i ε E + bi f i in (1 -e -iτ ) -b i f i in (1 -e iτ ) , (3.13)
F -(0, x, p, τ ) = f - in + i ε E - bi f i in (1 -e -iτ ) -b i f i in (1 -e iτ ) , (3.14)
G(0, x, p, τ ) = Re(f i in ) - ε 2E b i (f + in -f - in ) sin τ, (3.15)
H(0, x, p, τ ) = Im(f i in ) + ε 2E b i (f + in -f - in )(cos τ -1), (3.16)
where f ± in , f i in are the initial data of (3.1). One can also develop the gPC-SG-D and gPC-SG-N1 methods as in section 2. Since our focus is on gPC-SG-N2, we will just put gPC-SG-N1 and gPC-SG-D in the appendix for future reference.

The gPC-SG-N2 for the system with random inputs

Similar to the discussion in section 2.4, which allows the gPC order independent of ε, we implement the gPC-SG-N2 scheme. First, one focuses on the approximation of the phase equation (3.10). We insert the approximated solution

S(t, x, p, z) ≈ K j=1 Sj(t, x, p) ψj(z)
into (3.10) and conduct the Galerkin projection. We denote by #» S (t, x, p) = ( S1(t, x, p), • • • , SK (t, x, p)) T , the vector of unknown satisfying

∂t #» S + p∂x #» S + ∂xU ∂p #» S = 2 Iz E(x, z) ψ(z)π(z)dz, #» S (0, x, p) = 0, (3.17)
with ψ(z) = ( ψ1(z), • • • , ψK (z)) defined in (2.5). A 4-th order Runge-Kutta in time and pseudo-spectral method in space and velocity is used to solve #» S . Define the augmented functions

F + (t, x, p, τ, z) = W1(S(t, x, p, z), x, p, τ, z), F -(t, x, p, τ, z) = W2(S(t, x, p, z), x, p, τ, z), G(t, x, p, τ, z) = W3(S(t, x, p, z), x, p, τ, z), H(t, x, p, τ, z) = W4(S(t, x, p, z), x, p, τ, z) .
Then

∂tF + = ∂sW1∂tS, ∂xF + = ∂sW1∂xS + ∂xW1, ∂pF + = ∂sW1∂pS + ∂pW1 . Denote E ± (t, x, p, z) = 2E(x, z) -∂x(2U (x) ± E(x, z))∂pS(t, x, p, z) .
By (3.10) and (3.11), we have

                         ∂sW1 + p Ẽ + (s, x, p, z) ∂xW1 - ∂x(U + E) Ẽ + (s, x, p, z) ∂pW1 = - 1 ε ∂τ W1 + 2b i Ẽ + (s, x, p, z) (W3 cos τ + W4 sin τ ), ∂sW2 + p Ẽ -(s, x, p, z) ∂xW2 - ∂x(U -E) Ẽ -(s, x, p, z) ∂pW2 = - 1 ε ∂τ W2 - 2b i Ẽ -(s, x, p, z) (W3 cos τ + W4 sin τ ), ∂sW3 + p 2E(x, z) ∂xW3 + ∂xU 2E(x, z) ∂pW3 = - 1 ε ∂τ W3 + b i 2E(x, z) (W2 -W1) cos τ, ∂sW4 + p 2E(x, z) ∂xW4 + ∂xU 2E(x, z) ∂pW4 = - 1 ε ∂τ W4 + b i 2E(x, z) (W2 -W1) sin τ. (3.18)
How to obtain Ẽ ± (s, x, p, z) will be introduced in next subsection 3.2.1. Now we focus on the approximation of the profiles (W1, W2, W3, W4). We insert the Galerkin approximation of each component

(W1, W2, W3, W4)(s, x, p, z) = K j=1
(( W1)j, ( W2)j, ( W3)j, ( W4)j)(s, x, p) ψj(z), into (3.18) and conduct the gPC Galerkin projection. Denote the associated gPC coefficients vector of (

# » W 1, # » W 2, # » W 3, # » W 4) as previously, for instance for # » W 1, # » W 1 = (( W1)1, • • • , ( W1)K ) T . (3.19)
Then, we get the gPC system

                   ∂s # » W 1 + pJ ∂x # » W 1 -C∂p # » W 1 = - 1 ε ∂τ # » W 1 + 2b i J ( # » W 3 cos τ + # » W 4 sin τ ), ∂s # » W 2 + pL∂x # » W 2 + C∂p # » W 2 = - 1 ε ∂τ # » W 2 -2b i L( # » W 3 cos τ + # » W 4 sin τ ), ∂s # » W 3 + pH∂x # » W 3 = - 1 ε ∂τ # » W 3 + b i H( # » W 2 - # » W 1) cos τ, ∂s # » W 4 + pH∂x # » W 4 = - 1 ε ∂τ # » W 4 + b i H( # » W 2 - # » W 1) sin τ, (3.20)
where matrices J , L, C, H are given by

Jmn(s, x, p) = Iz 1 Ẽ + (s, x, p, z) ψm(z)ψn(z)π(z)dz, Lmn(s, x, p) = Iz 1 Ẽ -(s, x, p, z) ψm(z)ψn(z)π(z)dz, Cmn(s, x, p) = Iz ∂xE(x, z) Ẽ + (s, x, p, z) ψm(z)ψn(z)π(z)dz, Hmn(x) = Iz 1 2E(x, z) ψm(z)ψn(z)π(z)dz.
The initial conditions of

# » W 1, # » W 2, # » W 3, # » W 4 are # » W 1(0, x, p, τ ) = #» F + (0, x, p, τ ), # » W 2(0, x, p, τ ) = #» F -(0, x, p, τ ), # » W 3(0, x, p, τ ) = #» G(0, x, p, τ ), # » W 4(0, x, p, τ ) = #» H(0, x, p, τ ).
According to (3.13)- (3.16), by the gPC-SG method, the initial conditions for gPC coefficient vectors

#» F + , #» F -, #» G and #» H are given by #» F + (0, x, p, τ ) = f + in #» Q1 -iε bi f i in (1 -e -iτ ) -b i f i in (1 -e iτ ) #» Q2, (3.21)
#» F -(0, x, p, τ ) = f - in #» Q1 + iε bi f i in (1 -e -iτ ) -b i f i in (1 -e iτ ) #» Q3, (3.22)
#» G(0, x, p, τ ) = Re(f i in ) #» Q1 -εb i (f + in -f - in ) #» Q4 sin τ, (3.23)
#» H(0, x, p, τ ) = Im(f i in ) #» Q1 + εb i (f + in -f - in ) #» Q4(cos τ -1), (3.24)
where the vectors

#» Q1, #» Q2, #» Q3, #» Q4 are defined by #» Q1 = Iz ψ(z)π(z)dz , #» Q2(0, x, p) = Iz ψ(z)π(z) E + (0, x, p, z) dz , #» Q3(0, x, p) = Iz ψ(z)π(z) E -(0, x, p, z) dz , #» Q4(x) = Iz ψ(z)π(z) 2E(x, z) dz .

The fully discrete scheme for gPC-SG-N2

Monotonicity of S in terms of t

Assume U = U (x) ∈ R, E = E(x, z) > 0, ∂tS + p∂xS + ∂xU ∂pS = 2E(x, z), Sin(x, p, z) = 0.
By the method of characteristics,

dx dt = p, dp dt = ∂xU, x(0) = x0, p(0) = p0, then x = X(t, x0 , p0), p = P (t, x0, p0), thus d dt S(t, x, p, z) = 2E(X(t, x0, p0)),
and the analytic solution is given by

S(t, x, p, z) = Sin(x0, p0, z) + 2 t 0 E(X(µ, x0, p0, z))dµ = 2 t 0 E(X(µ, x0, p0, z))dµ .
Therefore S is an increasing function of t for each (x, p, z), since E > 0.

The time discretization of W1, W2, W3, W4 was defined in section 2.5.2. Denote E ± (t, x, p, z) = Ẽ ± (S(t, x, p, z), x, p, z). To find the values of Ẽ ± at s l , namely

Ẽ ± (s l , x, p, z) = E ± (t , x, p, z),
where t = S -1 (s l ), for each xj, p k , and quadrature points z (q) , q = 1, • • • , Ns, we apply an interpolation step shown here.

Search for the time interval [tn k , tn k +1] such that s l falls between the interval [S(tn k ), S(tn k +1)], thus t ∈ [tn k , tn k +1], since S is an increasing function of t. Linear interpolation is used to find t . Denote tn k = t (1) , tn k +1 = t (2) , and S(tn k ) = S (1) , S(tn k +1) = S (2) . By the Lagrange interpolation,

s l = S (1) t -t (2)
t (1) -t (2) + S (2) t -t (1) t (2) -t (1) , which gives t . The values E ± at all tn have been obtained, one can use linear interpolation to approximate the value of E ± at t , for each xj, p k , z (q) .

Using the Gauss quadrature rule with quadrature points z (q) and the corresponding weights ωq, q = 1, • • • , Ng, we update the matrices J , L , C at each (s l , xj, p k ), and H at each xj,

Jmn(s l , xj, p k ) = Ng q=1 1 Ẽ + (s l , xj, p k , z (q) )
ψm(z (q) )ψn(z (q) )ωq,

Lmn(s l , xj, p k ) = Ng q=1 1 Ẽ -(s l , xj, p k , z (q) )
ψm(z (q) )ψn(z (q) )ωq,

Cmn(s l , xj, p k ) = Ng q=1 ∂xE(xj, z (q) ) Ẽ + (s l , xj, p k , z (q) )
ψm(z (q) )ψn(z (q) )ωq,

Hmn(xj) = Ng q=1 1 2E(xj, z (q) )
ψm(z (q) )ψn(z (q) )ωq .

Step 1 We solve the transport part in x for quantities

# » W 1, # » W 2, # » W 3, # » W 4 in (3.20), ∂s # » W 1 + pJ ∂x # » W 1 = 0, ∂s # » W 2 + pJ ∂x # » W 2 = 0, ∂s # » W 3 + pH∂x # » W 3 = 0, ∂s # » W 4 + pH∂x # » W 4 = 0.
This transport step is treated similarly as in section 2.2, where we use a pseudo-spectral method in space and a three-stage Runge-Kutta method in time.

Step 2 We solve the transport part in p for W1, W2,

∂s # » W 1 -C∂p # » W 1 = 0, ∂s # » W 2 + C∂p # » W 2 = 0 .
Similar procedure is taken as in the previous step.

Step [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF] We now solve the non-singular source part

∂s # » W 1 = 2b i J ( # » W 3 cos τ + # » W 4 sin τ ), ∂s # » W 2 = -2b i L( # » W 3 cos τ + # » W 4 sin τ ), ∂s # » W 3 = b i H( # » W 2 - # » W 1) cos τ, ∂s # » W 4 = b i H( # » W 2 - # » W 1) sin τ .
We use the forward Euler method,

# » W n+1 1 = # » W n 1 + 2b i ∆sJ ( # » W n 3 cos τ + # » W n 4 sin τ ), # » W n+1 2 = # » W n 2 -2b i ∆sL( # » W n 3 cos τ + # » W n 4 sin τ ), # » W n+1 3 = # » W n 3 + b i ∆sH( # » W n 2 - # » W n 1 ) cos τ, # » W n+1 4 = # » W n 4 + b i ∆sH( # » W n 2 - # » W n 1 ) sin τ .
Other Runge-Kutta methods can also be used to solve this ODE system.

Step 4 Finally, we solve the highly oscillatory part

∂s # » W 1 = - 1 ε ∂τ # » W 1, ∂s # » W 2 = - 1 ε ∂τ # » W 2, ∂s # » W 3 = - 1 ε ∂τ # » W 3, ∂s # » W 4 = - 1 ε ∂τ # » W 4.
We use the Fourier transform Ŵ1(ζ) of W1(τ ) in the τ variable where ζ is the corresponding Fourier variable, then

∂s # » Ŵ1 = - iζ ε # » Ŵ1,
which is solved by the backward Euler method in time,

# » Ŵ1 n+1 (xj, p k ) = I + iζ∆s ε I -1 # » Ŵ1 n (xj, p k ) .
The same method is used for calculating

# » W 2, # » W 3, # » W 4.
4 Numerical examples for the surface hopping model

• Example 4.1 Consider x, p ∈ [-2π, 2π
]. We use the similar data given in the numerical example in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], with the following initial conditions,

f + (t = 0, x, p) = f -(t = 0, x, p) = (1 + 0.5 cos(x)) e -p 2 /2 √ 2π , f i (t = 0, x, p) = [1 + 0.5 sin(x)) + i(1 + 0.5 cos(x))] e -p 2 /2 √ 2π .
The expressions for E, bi, and b ± are given by

E(x, z) = (1 -cos(x/2) + √ ε)(1 + 0.5z), bi(x, p) = - 1 2 sin(p + 1), b ± = 0.
Periodic boundary conditions are considered in x while the p domain is chosen large enough so f ±,i all vanish outside the domain (thus a periodic boundary condition in p can be used).

Without loss of generality we set U = 0. For all the following tests, we choose Nτ = 8, Np = 32 and Nc = 32 quadrature points for gPC-SC methods (with a very small mesh size in x so this solution is used as the reference solutions).

In Figure 6, for gPC-SG-N2 scheme, even if the mesh size is much larger than the wave length ε = 5×10 -3 , the solution f ±,i (the mean and standard deviation), as well as ρ ±,i , still agrees with the reference solution at the grid points, despite a moderate gPC order K = 4.

Figure 7 shows a similar result as in Figure 6, using Example 4.1 except that here E(x, z) = O(1). As discussed in Remark 2.2, small ∆t = O( √ ε∆x) needs to be chosen in Figure 6 due to the CFL condition in the transport steps for E = O( √ ε). However, larger ∆t = O(∆x), which is independent of ε, can be used in Figure 7 since now E(x, z) = O(1).

25

Conclusions

In this paper, based on a nonlinear geometric optics (NGO) based numerical method developed in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF], with a new "time" variable defined from the phase, we obtain a stochastic Galerkin (SG) method for highly oscillatory transport equations that arise in semiclassical models of non-adiabatic quantum dynamics, in which the potential energy surfaces are assumed to be random, due to uncertainties in modeling or measurement errors. We prove that the generalized polynomial chaos (gPC) based SG method applied directly to the models will require the order of gPC to depend on the possibly very small wave length, while the new method does not have such a requirement. This important property allows us to use this method to solve these highly oscillatory problems with uncertain coefficients with all numerical parameters independent of the wave length, yet still capture the solution statistics pointwisely.

There are several projects along this direction. First it will be desirable to develop an implicit scheme for the transport steps in order to obtain an improved time step constraint when the bad gap becomes very small. Second, methods for higher dimensional-in space, velocity as well as the random variables-remain to be developed.

A Appendix: The gPC-SG-N1 and gPC-SG-D for the surface hopping model with random inputs

A.1. The gPC Approximation

In this part, we detail the gPC-SG-N1 strategy for model (3.1), based on the method developed in [START_REF] Crouseilles | Nonlinear geometric optics method based multi-scale numerical schemes for highly-oscillatory transport equations[END_REF] for the deterministic problem. Solving for S follows the same strategy presented in section 3.2.1 on the gPC-SG-N2 method. Now we focus on the approximation of the profiles (F + , F -, G, H). One inserts the Galerkin approximation of each component (F + , F -, G, H)(t, x, p, z) = K j=1 (F + j , F - j , Gj, Hj)(t, x, p) ψj(z), into (3.11) and conduct the gPC Galerkin projection. Denote the associated gPC coefficients vector of ( #» F + , #» F -, #» G, #» H) as previously, for instance for

#» F + , #» F + = (F + 1 , • • • , F + K ) T . (A.1)
Then, we get the gPC system Note that the matrix H is symmetric. The initial conditions have been discussed in subsection 3.2.1.

                       ∂t #» F + + p∂x #» F + -C∂p #» F + = 1 ε (-2C + H)∂τ #» F + + 2b i ( #» G cos τ + #» H sin τ ), ∂t #» F -+ p∂x #» F -+ C∂p #» F -= 1 ε (-2C -H)∂τ #» F --2b i ( #» G cos τ + #» H sin τ ), ∂t #» G + p∂x #» G = - 2C ε ∂τ #» G -b i #» F + cos τ + b i #» F -cos τ, ∂t #» H + p∂x #» H = - 2C ε ∂τ #» H -b i #» F + sin τ + b i #» F -sin τ,

A.2. The fully discrete scheme gPC-SG-N1

To solve (A.2), we use a time splitting procedure. We split the equation into four steps: two transport steps (in x and in p), a highly oscillatory part and a non-singular source part. We detail how we solve each step in the sequel.

Step 1 We solve the transport part in x for each quantity F + , F -, G, H. For example, to solve the transport in F + , ∂t #» F + + p∂x #» F + = 0, we use a spectral method in space and an exact integration in time in the Fourier space, that is #» F +,n+1 (ξ) = e -ipξ∆t #» F +,n (ξ), (A.3) where ξ denotes the Fourier space variable and F + (ξ) the corresponding (discrete) Fourier transform. Finally, #» F +,n+1 is obtained by the inverse Fourier transform.

Step 2 We solve the transport part in p for F ± ,

∂t #» F ± ∓ C∂p #» F ± = 0.
We also use a spectral method in p and an exact integration in time in the Fourier space, that is #» F ±,n+1 (η) = e ±iηC(x)∆t #» F ±,n (η), (A.4) where η is the Fourier velocity variable. Finally #» F ±,n+1 are obtained by the inverse Fourier transform.

Step 3 We now solve the non-singular source part

∂t #» F + = 2b i ( #» G cos τ + #» H sin τ ), (A.5) ∂t #» F -= -2b i ( #» G cos τ + #» H sin τ ), (A.6) ∂t #» G = -b i #» F + cos τ + b i #» F -cos τ, (A.7) ∂t #» H = -b i #» F + sin τ + b i #» F -sin τ . (A.8)
This linear system with time independent coefficients can be solved exactly.

Step 4 Finally, we solve the highly oscillatory part

∂t #» F + = 1 ε (-2C + H)∂τ #» F + , ∂t #» F -= 1 ε (-2C -H)∂τ #» F -, ∂t #» G = - 2C ε ∂τ #» G, ∂t #» H = - 2C ε ∂τ #» H.
We use the Fourier transform F (ζ) of F (τ ) in the τ variable where ζ is the corresponding Fourier variable. Then, we have

∂t #» F + = iζ ε (-2C + H) #» F + , ∂t #» F -= iζ ε (-2C -H) #» F -, ∂t #» Ĝ = iζ ε (-2C) #» Ĝ, ∂t #» Ĥ = iζ ε (-2C) #» Ĥ,
which is solved by the backward Euler method in time, Notice that we have already solved the gPC coefficients vector #» S (t, x, p) from (3.17). This allows to compute ∂p #» S (t, x, p) spectrally (using Fourier in p) which is used in the expression of matrix I -iζ ε ∆t(-2C ± H) . Note that the matrices in (A.9)-(A.12) are invertible since they are symmetric and have real eigenvalues.

A.

#» g n + ∆t -b i #» f +,n + b i #» f -,n + 2 ε W #» h n -(b i ) 2 ∆t #» g n - ∆t ε 2 P #» g n ,
where the matrices WK×K and PK×K are defined by W(x) ml = Iz E(x, z)ψm(z)ψ l (z)π(z)dz, 1 ≤ m, l ≤ K, P(x) ml = Iz E(x, z) 2 ψm(z)ψ l (z)π(z)dz, 1 ≤ m, l ≤ K.

Since P is symmetric and positive definite, with positive eigenvalues, so does the matrix I + (b i ) 2 (∆t) 2 I + ∆t 2 ε 2 P , thus it is invertible. Since #» g n+1 is known, we can compute #» f +,n+1 , #» f -,n+1 , #» h n+1 by gPC approximation of (3.4), (3.5) and (3.7) as

#» f +,n+1 = #» f +,n + b i ∆t( #» g n + #» g n+1 ), #» f -,n+1 = #» f -,n -b i ∆t( #» g n + #» g n+1 ), #» h n+1 = #» h n - ∆t ε W( #» g n + #» g n+1 ).

  partition in time tn = n∆t (∆t > 0 the time step) of a time interval [0, T ], n = 0, 1, • • • , N , N ∆t = tN and in space xj = j∆x, j = 0, 1, • • • , Nx, ∆x = 1/Nx of the spatial interval [0, 1]. A uniform mesh is assumed in the direction τ ∈ T = [0, 2π], τ l = l∆τ , l = 0, • • • , Nτ , ∆τ = 2π/Nτ . Denote by V n j (τ ) ≈ V (tn, xj, τ ) and S n j ≈ S(tn, xj) the discrete unknowns at time tn, evaluated at xj.

Figure 1 :

 1 Figure 1: Example 2.1. Mean of real parts of u at t = 0.1, ε = 1, 0.1 and 0.01. N x = 2000, ∆t = 10 -5 . Stars: gPC-SG-D with K = 4. Solid lines: reference solution by gPC-SC-D.

Figure 2 :

 2 Figure 2: Example 2.1. Mean of real parts of u at t = 0.1, ε = 5 × 10 -3 . N x = 2000, ∆t = 10 -5 . Stars: gPC-SG-D. Solid lines: reference solution by gPC-SC-D.

Figure 3 :

 3 Figure 3: Example 2.1. Mean of real and imaginary parts of u at t = 0.25, ε = 0.01. ∆x = π/32, ∆t = 0.01 for the first row (gPC-SG-N1), and ∆x = π/1000, ∆t = 10 -4 for the second row (gPC-SG-N1). ∆x = π/1000, ∆t = 5 × 10 -5 (gPC-SC-D). Stars: gPC-SG-N1 with K = 4. Solid lines: reference solution by gPC-SC-D.

Figure 4 :Figure 5 :

 45 Figure 4: Example 2.1. Mean and standard deviation of real and imaginary parts of u and their zoomed in solutions at t = 0.25, ε = 5 × 10 -3 . ∆x = π/32, ∆t = 0.01, K = 4 (gPC-SG-N2), and ∆x = π/1000, ∆t = 5 × 10 -5 (gPC-SC-D). Stars: gPC-SG-N2. Solid lines: reference solution by gPC-SC-D.

Figure 6

 6 Figure 6(a): Example 4.1. ε = 5 × 10 -3 , t = 0.5, N x = 2500, ∆t = 5 × 10 -3 (gPC-SC-D), and N x = 32, ∆t = 5 × 10 -5 (gPC-SG-N2). Mean of the space dependence of f ± , Re(f i ) and Im(f i ) (and their zoomed in solutions) at p = 0. Stars: gPC-SG-N2 with K = 4. Solid lines: reference solution by gPC-SC-D.

Figure 6

 6 Figure 6(b): Standard deviation of the space dependence of f ± , Re(f i ) and Im(f i ) (and their zoomed in solutions) at p = 0.

Figure 6

 6 Figure 6(c): Mean of the space dependence of the densities ρ ± , Re(ρ i ) and Im(ρ i ) (and their zoomed in solutions).

Figure 7 :

 7 Figure 7: Example 4.1, with E(x, z) = (10 -cos(x/2))(1 + 0.5z). ε = 0.01, t = 0.3, N x = 2000, ∆t = 5×10 -3 (gPC-SC-D), and N x = 32, ∆t = 0.02 (gPC-SG-N2). Mean and standard deviation of the space dependence of the densities ρ ± , Re(ρ i ) and Im(ρ i ). Stars: gPC-SG-N2 with K = 4. Solid lines: reference solution by gPC-SC-D.

(A. 2 )

 2 where matrices C ∈ RK×K and H ∈ RK×K are defined byCmn(x) = Iz E(x, z) ψm(z) ψn(z)dz, 1 ≤ m, n ≤ K. H km (t, x, p) = K l=1 ∂p Sl (t, x, p)D lmk (x), 1 ≤ k, m ≤ K, with ∂p Sl (t, x, p) (l = 1, • • • , K) the gPC coefficients of ∂pS(t, x, p, z), i.e., ∂p #» S (t, x, p) = (∂p S1(t, x, p), • • • , ∂p SK (t, x, p)) T ,and the tensor D(x) is given byD lmk (x) =Iz ∂xE(x, z) ψl (z) ψm(z) ψk (z)π(z)dz.

#

  » F +,n+1 = I -iζ ε ∆t(-2C + H) -1 #» F +,n , (A.9) #» F -,n+1 = I -iζ ε ∆t(-2C -H) -1 #» F -,n , (A.10) #» Ĝn+1 = I + 2 iζ ε ∆tC -1 #» Ĝn , (A.11) #» Ĥn+1 = I + 2 iζ ε ∆tC -1 #» Ĥn . (A.12)

  [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF] The gPC-SG-D schemeWe briefly introduce the gPC-SG-D for the surface hopping model with random inputs. One first inserts the gPC expansions (f + , f -, g, h)(tn, x, p, z) ≈ ), we conduct the Galerkin projection to get#» g n+1 = I + (b i ) 2 (∆t) 2 I + ∆t 2 ε 2 P

	K
	( f +,n j	, f -,n j	, gn j , hn j )(x, p) ψj(z),
	j=1
	into equation (3.8). Then, denoting as previously #» g n = (g n 1 , gn 2 , • • • , gn K ) (and the same notations for #» f +,n , #» f -,n , #» h n -1

* This work was partially supported by NSF grants DMS-1522184 and DMS-1107291: RNMS KI-Net, and by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation.