Chunyao Ma

Sameer Alam

Qing Cai

Daniel Delahaye
email: daniel.delahaye@enac.fr

A Machine Learned Traffic Flow Coordination Framework for Flow-Centric Airspace

Keywords: flow coordination, flow-centric, traffic prediction, transformer neural networks, reinforcement learning

Air traffic flow coordination at major flow intersections is a key enabler for flow-centric airspace concepts. This paper develops a flow-centric air traffic flow coordination framework to improve air traffic flow efficiency through flow identification, prediction, and re-routing at the Nominal Flow Intersections (NFIs). To identify the NFIs, a graph-based flow pattern consistency approach is proposed to model and analyze daily air traffic flow patterns. To predict future traffic demands at the identified NFIs, a transformer encoder-based neural network is adopted to learn the relations among the flow of flights at the NFIs. The acceptable flow limits at the NFIs are then determined by phase transitions of the flow efficiency versus the traffic demand. Finally, to avoid the predicted demand exceeding the identified flow limit and improve the flow efficiency, a reinforcement learning-based flow re-routing agent is trained to dynamically assign alternative routes to air traffic flows based on the evolving flow states. The agent's performance is quantified by the flight time reduction in the flows without exceeding the flow limits. The re-routing model is trained and tested on a busy NFI that handles cross-border flows between Bordeaux and Madrid/Barcelona control centers, using ADS-B data for Dec 2019 in European airspace. Results show that, compared with the originally planned flows, the travel time of each flight is reduced by 322.168 seconds on average on a 2-hour basis.

I. INTRODUCTION

The scalability limit of traditional sector-based Air Traffic Control (ATC) services, i.e., difficulty in subdividing heavily loaded sectors, is becoming a barrier to the sustainable growth of air traffic. Researchers have started examining and testing the concept of flow-centric operation [START_REF] Ma | Dynamic air traffic flow coordination for flow-centric airspace management[END_REF], which views the whole regional airspace and controls groups of flights throughout their flight segments in a region. It opens the opportunity to distribute air traffic more efficiently in the airspace without being constrained by sector boundaries.

One primary challenge in flow-centric operation is the efficient coordination of air traffic flow at the intersections to avoid inefficiencies that may jeopardize flight safety [START_REF] Gerdes | From free-route air traffic to an adapted dynamic mainflow system[END_REF]. Research focusing on sector-based air traffic coordination, such as sector traffic prediction and flow optimization for workload balancing between sectors [START_REF] Chen | A network based dynamic air traffic flow model for en route airspace system traffic flow optimization[END_REF], no longer adapts flow-centric operations where coordination is primarily used to avoid potential inefficiencies or conflicts between the intersecting air traffic flows. Therefore, it is crucial to develop a flow-centric framework that can identify, predict, and dynamically coordinate the evolving air traffic flows. For instance, traffic flow can be strategically re-routed when the predicted demand exceeds the acceptable flow limit at flow intersections.

Effective air traffic flow identification is the cornerstone for flow-centric air traffic analysis, prediction, and coordination [START_REF] Gerdes | From free-route air traffic to an adapted dynamic mainflow system[END_REF]. In the literature, air traffic flow has been described according to airspace configuration, such as the groups of flights transiting through area control centers, waypoints, sectors, and airways [START_REF] Ma | Sector entry flow prediction based on graph convolutional networks[END_REF]. Such a characterization of air traffic flow fits the traditional Air Traffic Control (ATC) paradigm where ATC units are geographical sectors and flights follow airways consisting of fixed waypoints. However, flow-centric operations require identifying the evolving air traffic flow patterns, such as flow locations and structures, disregarding the fixed airways and sectors. Fig. 1 shows one-day flight trajectories (pink lines) in French airspace, where nearly 50% above 19500ft is free route airspace (a potential coupled working method to flow-centric operations). It can be observed that the positions of airports and scheduled flights between airports restrict air traffic to an appropriate pattern of main flows [START_REF] Chunyao | Air traffic flow representation and prediction using transformer in flow-centric airspace[END_REF]. Thus, identifying Nominal Flow Intersections (NFI) by constructing and analyzing air traffic flow patterns is the first enabler of effective flow coordination.

In addition to NFI identification, effective flow coordination requires constantly viewing air traffic demand according to the available capacity at the NFIs. In the literature, aggregated air traffic flow prediction mainly predicts the number of flights transiting through different airspace locations, such as entryexits, origin-destinations, and air routes [START_REF] Zhang | A hierarchical heuristic approach for solving air traffic scheduling and routing problem with a novel air traffic model[END_REF]. The time series of traffic demand at single or multiple airspace locations are used as the input to predict the future demand using Long Short-Term Memory (LSTM) [START_REF] Gui | Machine learning aided air traffic flow analysis based on aviation big data[END_REF], Convolutional Neural Networks (CNNs) [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], and Graph Convolutional Networks (GCNs) [START_REF] Ma | Sector entry flow prediction based on graph convolutional networks[END_REF]. While time series flow features provide prediction models with information on the number of flights, the most important characteristic of air traffic flow, the "flow of flights in the airspace," can not be depicted by a time series, limiting the prediction accuracy [START_REF] Cai | Temporal attention aware dual-graph convolution network for air traffic flow prediction[END_REF]. Moreover, most flow prediction models adopt recurrent neural networks, in which data are processed one after another to learn the sequential relations in a time series, making it difficult to track long-term relations in the input sequence [START_REF] Mutanga | Hate speech detection in twitter using transformer methods[END_REF]. Therefore, developing an effective flow feature representation to describe the dynamics inside air traffic flow and a prediction model that can learn relevant features between sequential elements far from each other is the second enabler of flow coordination. Besides air traffic prediction at the NFIs, it is essential to know the acceptable flow limit at the NFIs, given the NFIs can be overloaded due to excessive air traffic flows. In traditional sector-based ATC systems, the sector capacity is usually quantified as the maximum hourly number of flights that may enter a sector over a sustainable period while maintaining a safe, orderly, and efficient traffic flow [START_REF] Welch | En route sector capacity model final report[END_REF]. Similarly, in flowcentric ATC systems, if the demand at an NFI exceeds the acceptable limit, the efficiencies and safety of air traffic flow at the NFI can be compromised. Airspace capacity estimation methods in the literature, such as the Monitor Alert Parameter (MAP) model [START_REF] Marr | Controller workload-based calculation of monitor alert parameters for en route sectors[END_REF], the full workload model (considering transit, recurring, conflict, and background workload) [START_REF] Welch | En route sector capacity model final report[END_REF], the sector merge/split-based model [START_REF] Gianazza | Learning air traffic controller workload from past sector operations[END_REF] are primarily tailored to sector-based operations which cannot accommodate the flowcentric requirements. Therefore, a flow-centric algorithm for flow acceptance limit identification at the NFIs is the third enabler of effective flow coordination.

When the predicted air traffic demand exceeds the corresponding acceptable flow limit at an NFI, a flow coordination model will be required to re-configure the air traffic demand at the NFIs in advance to restrict the air traffic flow within a level not overloading the system excessively [START_REF] Lin | From aircraft tracking data to network delay model: A data-driven approach considering en-route congestion[END_REF]. In traditional sectors, air traffic overload is managed by sector operations such as merging and splitting [START_REF] Liu | A framework for strategic online en-route operations: Integrating traffic flow and strategic conflict managements[END_REF]. Under the flow-centric paradigm, traffic flow density and complexity change over time, rendering static flow control operations underloaded or overloaded. Dynamic flow re-configuration, such as flow merge/split/rerouting, gives flow-centric airspace an option to address the overload at NFIs without compromising the flow demand. Therefore, developing a dynamic flow re-configuration model subject to the time-varying traffic flows is the fourth enabler of efficient flow coordination decisions.

Given the above analysis, this paper proposes a dynamic air traffic flow coordination framework to identify, predict, and re-route air traffic flows to enable more efficient flowcentric airspace management. Firstly, a graph-based flow pattern consistency analysis approach is proposed to identify nominal air traffic flow intersections (NFIs) in the airspace through modeling and analyzing daily air traffic flow patterns. Secondly, a text-enriched flow feature representation is proposed to describe "the flow of flights" in the airspace using a "text paragraph" composed of the time and flight sequences at the NFIs. Compared to traditional time series representation describing "how many flights were in the airspace," it describes not only the number of flights but also the flow patterns shaping the movement of air traffic. A transformerencoder-based neural network model is adopted to learn correlations among flow sequences to predict the future traffic demand at the NFIs. Thirdly, for each NFI, the acceptable flow limit is determined by identifying the phase transition of the flow efficiency, characterized by the flight transition duration from its neighboring NFIs versus the traffic demand during different periods. Finally, a reinforcement learningbased flow re-routing model is proposed to dynamically assign alternative routes to air traffic flows, especially when the predicted demand exceeds the acceptable limit, considering the effects of route assignment on flow efficiency and demand at NFIs during both current and future periods. Secondly, inspired by Natural Language Processing (NLP), the "flow of flights in the airspace" during a period is described by a text paragraph consisting of the sequences of flights transiting through the NFIs. In NLP, a sentence consists of a sequence of words governed by grammar and is used by neural networks to extract linguistic features for downstream tasks such as next-word prediction and sentiment analysis [START_REF] Zhou | Integrating nlp and context-free grammar for complex rule interpretation towards automated compliance checking[END_REF]. Meanwhile, the "flow of flights in the airspace" during a period consists of sequences of flights at different locations governed by recurrent traffic flow patterns. Analogously, learning from flow features represented by flight sequences at NFIs during different periods has the potential to achieve air traffic prediction tasks. Moreover, the transformer encoder architecture, using the attention mechanism to enable dynamic highlighting of relevant features between input elements and parallel sequence processing [START_REF] Vaswani | Attention is all you need[END_REF], is adopted to capture the relations among flows represented by the flight sequences to predict the flight demand at the NFIs during different periods.

Thirdly, when an NFI is overloaded, air traffic congestion can be induced, and it will take longer for air traffic transiting to the overloaded NFI due to regulatory measures such as vectoring and speed control [START_REF] Baumgarten | The impact of hubbing concentration on flight delays within airline networks: An empirical analysis of the us domestic market[END_REF]. Therefore, the flow acceptance limit at each NFI is identified as the demand above which the flow efficiency reduces significantly, i.e., the flight transition duration from neighboring NFIs to the NFI increases distinctly.

Finally, to improve flow efficiency while preventing adverse effects of flow re-routing during the current period over future air traffic, this paper considers flow re-routing as a sequential decision-making problem [START_REF] Singh | How are reinforcement learning and deep learning algorithms used for big data based decision making in financial industriesa review and research agenda[END_REF] where the action of alternative route assignment to a flow has consequences affecting future flow efficiency and flow demand at the NFIs. Then, deep reinforcement learning algorithms are adopted to train an agent to dynamically assign alternative routes to air traffic flows based on the evolving traffic flow states, considering the effects of route assignment on both current and future periods to reduce flight transition time in the airspace without exceeding the flow limit at the NFIs. The following sections illustrate the detailed procedures of the four steps.

B. NFI Identification

This paper identifies the NFIs through flight trajectory analysis using ADS-B data, including intersection points clus-tering, daily flow pattern representation, and graph analysis.

1) Intersection Points clustering: A trajectory intersection point is identified when two flights fly at adjacent flight levels, and their trajectories intersect when mapped onto the earth's surface. The time difference between the two flights passing the intersection is within 30 minutes, a widely used time window in the literature on air traffic flow management, flight planning, and capacity estimation [START_REF] Balakrishnan | Optimal large-scale air traffic flow management[END_REF].

With the trajectory intersections computed from flight track data, clustering algorithms are adopted to organize the intersection points into representative NFIs. Hierarchical clustering decomposes the data based on group similarities to find a multilevel hierarchy of clusters [START_REF] Cai | Hierarchical clustering of bipartite networks based on multiobjective optimization[END_REF]. Considering the hierarchical organization of air traffic flows connecting regional feeders to international hubs [START_REF] Seabra | Determinants of brazilian international flights: The role of hub-and-spoke and infrastructure variables[END_REF], [START_REF] Adler | The multi-airline p-hub median problem applied to the african aviation market[END_REF], this paper adopts single-linkage hierarchical clustering of trajectory intersections to discover the natural organization of flow intersections:

obj = arg max V (n) n-1 r=1 n s=r+1 D(r, s) s.t. D(r, s) = min(dist(x ri , x sj )) (1)
where x ri and x sj are the i and j-th object in cluster r and s. n is the number of clusters. V (n) denotes the cluster centers.

2) Flow Pattern Analysis: If the number of clusters is too large, the clustering result will be susceptible to small fluctuations in air traffic flow. When the number of clusters is too small, the identified NFIs can largely deviate from actual flow paths. The consistent and dependable performance of airspace users is an essential requirement for improving ATM system predictability [START_REF] Murc ¸a | Flight trajectory data analytics for characterization of air traffic flows: A comparative analysis of terminal area operations between new york, hong kong and sao paulo[END_REF]. The identified NFIs should represent the consistency of air traffic flow patterns while considering daily flow alternations. Thus, the geographical and structural consistencies of air traffic flows versus the number of clusters are analyzed to determine the optimal number. The flow pattern consistency is evaluated from two perspectives: a) geographical consistency in NFI locations and b) structural consistency in flow connectivity between NFIs.

The daily air traffic flow pattern is represented by a weighted graph G = (V, E) [START_REF] Ma | Airway network management using braess's paradox[END_REF], where V is the set of nodes denoting the NFIs. The flow connectivity between nodes is described by the weighted edges E. The edge weight is quantified by the air traffic volume on the edge. a) Geographical Consistency in NFI locations: A nearest-neighbor-based analysis is conducted to measure the geographical consistency of the daily NFI locations. Given a number of n nodes in the graph G k constructed for the k-th day, let

V k = {v k,1 , v k,2 , ..., v k,i , ..., v k,n } represent the set of nodes. Similarly, let V k+1 = {v k+1,1 , v k+1,2 , .
.., v k+1,j , ..., v k+1,n } represent the set of nodes in the graph G k+1 for day k + 1. For each node v k,i in V k , the proposed algorithm searches for its nearest node v k+1,ai in V k+1 according to the great circle distance [START_REF] Prasetya | Resolving the shortest path problem using the haversine algorithm[END_REF]. a i is the index of the identified nearest node of v k,i in V k+1 .

Through the above calculations, the nearest neighbouring node v k+1,ai of node v k,i for i = 1, ..., n can be determined. Similarly, the nearest neighbouring node v k,bj of node v k+1,j for j = 1, ..., n can be determined. b j is the index of the identified nearest node of v k+1,j in V k . Therefore, two sets of matched node pairs can be obtained:

M P k : (v k,1 , v k+1,a1 ), (v k,2 , v k+1,a2 ), ..., (v k,n , v k+1,an ) and M P k+1 : (v k,b1 , v k+1,1 ), (v k,b2 , v k+1,2 ), ..., (v k,bn , v k+1,n ).
Then, the geographical consistency of node locations is quantified as the number of mutually matched nodes divided by the total number of nodes, which is formulated as:

gc 1 = |M P k ∪ M P k+1 | n (
gc 2 = l-1 i=1 l j=i+1 min(w i,j c k , w i,j c k+1 ) n-1 i=1 n j=i+1 (w i,j k + w i,j k+1 ) (3)
The optimal cluster number n = N v is determined by graph pattern analysis to find the saddle point, i.e., local maxima, of the daily graph pattern consistency versus the cluster numbers:

P (n) = arg min n (P (n) -gc1(V (n))+gc2(V (n)) 2 ) s.t. P ′ (n) = 0; P ′′ (n) < 0 (4)
where P (n) is the polynomial approximation of the average of gc 1 and gc 2 (scaled between 0 and 1). As shown in Fig. 3, the circle marks the identified saddle point N v on P (n). → Y. Fig. 4 shows the neural network structure for flow prediction at NFIs, consisting of tokenization, embedding, transformer encoder blocks, and a fully connected layer. a) Tokenization: Given an input sequence, including the time and the flight sequences at the NFIs, word tokenization [START_REF] Sun | Spectral-spatial feature tokenization transformer for hyperspectral image classification[END_REF] is applied to convert elements in the input sequence S into a list of integers X = x 1 , x 2 , ..., x m } that can be embedded into a vector space.

b) Embedding: The token embedding layer converts the tokenized flow sequence into a list of vectors. Given a token x i , its embedding T E(x i ) ∈ R d model can be represented as: T E(x i ) = W e x i where W e is the embedding matrix and d model is the embedding dimension. Positional embedding P E(x i ), generated using trigonometric functions (sines and cosines) to provide a sinusoidal pattern encoding the position information, are added to the token embeddings as the input to the transformer encoders. For instance, multiplying x2, the embedded vector for the callsign "AFR378", by W Q , W K , W V produces the query, key, and value vectors q2, k2, and v2 associated with "AFR378". When processing the self-attention for "AFR378", its first attention weight related to element "14:15" would be S( q2×k1 √ d k ). Multiplying each value vector by the attention weight and summing up the weighted value vectors produce the output of the self-attention layer.

c) Transformer Encoder Blocks:

The output of the input sequence embedding E(X) in which E(x i ) = T E(x i ) + P E(x i ) is the input to the stack of N e transformer encoder blocks to process each element in the input sequence and compile the information it captures into a context tensor.

Each encoder block has two major components: the multihead self-attention mechanism and the position-wise fully connected feed-forward network [START_REF] Zhao | Speech transformer with speaker aware persistent memory[END_REF]. The self-attention mechanism allows the encoder to look at other elements in the input sequence when encoding a specific element. It creates three vectors, query Q = E(X)W Q , key K = E(X)W K , and value V = E(X)W V , for each input element by multiplying the embedding with weight matrices W Q , W K , and W V . The three vectors are used to score the relevance of other elements in the sequence against the specific element calculated by: Sof tmax( QK T √ d head )V . d head denotes the dimension of the head, e.g., the dimension d k of the key. Multi-headed attention runs through a self-attention mechanism several times in parallel, which allows the model to jointly attend to information from different representation subspaces at different positions [START_REF] Vaswani | Attention is all you need[END_REF]. The output of the multi-head attention is added elementwise to the original input embeddings and normalized using layer normalization. Its output, represented by A&N 1 , is then passed through the feedforward network F F N (A&N 1 ) of the encoder block. The final output of the encoder block is:

A&N 2 = LayerN orm(F F N (A&N 1 ) + A&N 1 ) (6)
The output from the stack of N e encoders is forwarded to a fully connected layer Y = F F N (A&N 2 (N e )) to obtain the flow prediction results for different prediction windows, i.e., future period 1, period 2,..., and period N p .

D. NFI Flow Acceptance Limit Identification

When the demand at an NFI is above the acceptable limit during a period, air traffic congestion can happen, and it will take a significantly larger time cost for air traffic transiting to the overloaded NFI on their flight paths due to regulatory measures such as vectoring and speed control.

During period t, flight transition durations from NFI v j to NFI v i are denoted as H ij,t . H ij,t is normalized by the daily minimum duration h ij,t regarding t to reduce the effects of daily fluctuations in air traffic flow:

N H ij,t = H ij,t h ij,t (7)
Let N H i,t = {N H ij,t } denote the normalized transition duration to v i of all NFIs connected to v i during period t.

Let N H i = {N H i,t }
If there is no abrupt change detected, indicating no explicit trend concerning the flow demand based on the observations, the flow limit at the NFI is determined as the maximum flow demand observed.

E. Reinforcement Learning based Flow Coordination

This section adopts reinforcement learning algorithms to train an agent for flow re-routing to avoid flow excess and improve flow efficiency. During the t-th period, the agent assigns a route for a main flow according to the re-routing policy π θ (a t |s t ) parameterized by θ. π θ (a t |s t ) determines the agent's action a t ∈ A, i.e., which route to take, regarding the traffic state s t ∈ S. S represents the state space of the traffic flow described by the traffic volumes at the NFIs. A represents the action space, consisting of the alternative routes a flow can take. R (s t , a t ) denotes the agent's reward of taking action a t in s t based on whether the flow excess is avoided and how much the flow efficiency is improved:

R (s t , a t ) = e t T (0) t -T (st,at) t /N f lights ( 9
)
where e t = 1 otherwise 0 if the flow excess is avoided and no new excess at other NFIs is induced. T (0) t is original flight exit time in the airspace during the t-th period, while T (st,at) t represents the time by taking an alternative route a t . N f lights denotes the number of flights transited through the re-routing area during the t-th period.

The agent's objective is maximizing the expected cumulative reward over time, represented by the objective function:

J(θ) = E T t=0 γ t R(s t , a t ) ( 10
)
where T is the learning horizon. γ is the discounting factor. This paper adopts the Proximal Policy Optimization (PPO) [START_REF] Schulman | Proximal policy optimization algorithms[END_REF] to train the agent, which introduces a surrogate objective with a clipping mechanism preventing significant policy updates to improve sample efficiency and stability:

arg max θ E min( π θ (at|st) π θo (at|st) A π θo (s t , a t ) , g (A π θo (s t , a t ))) s.t. A π (s t , a t ) = δ t + (γλ)δ t+1 + ... + (γλ) T -(t+1) δ T -1 g (A) = (1 + ϵ) A, A ≥ 0 (1 -ϵ) A, A ≤ 0 δ t = R (s t , a t ) + γV π (s t+1 ) -V π (s t ) (11
) where π θo (a t |s t ) is the old policy from the previous iteration. A π (s t , a t ) is the Generalized Advantage Estimation (GAE) estimating the advantage of taking action a t compared with the average action in state s t . V π (s t ) is the state value function of the expected cumulative reward starting from s t . ϵ and λ are hyperparameters controlling the extent of clipping and the balance between bias and variance of GAE, respectively.

III. EXPERIMENTAL STUDY

To verify the efficacy of the proposed flow coordination framework, an experimental study has been carried out on the French airspace using one-month ADS-B data from 1 to 31 December 2019, comprising 158,856 flights. This study focuses on the en-route air traffic above 10,000 ft.

A. NFI Identification

This paper calculated the geographical consistency gc 1 and the structural consistency gc 2 against different cluster numbers ranging from 100 to 1500. A "saddle point" is observed for gc 1 and gc 2 at cluster number 605. Therefore, this paper takes 605 as the number of clusters. The cluster centers are determined as the NFIs. The 605 identified NFIs are denoted by "NFI1" to "NFI605" in this paper. Fig. 5 shows the graph representation for the one-month air traffic flow with 605 NFIs. This graph can depict the nodal hierarchy of air traffic flows ranging from regional feeders to international hubs, such as Paris and Geneva. The en-route air traffic flows are organized as a series of "spokes" connecting the traffic hubs or connecting outlying areas to a hub area. Table I shows the quantified prediction performance of the proposed method tested on three busy NFIs in the airspace, including NFI460 over the Paris area control center, NFI365 over the Geneva area control center, and NFI135 over the The blue lines in Fig. 7 show instances during six days in December 2019 where the traffic demand on NFI135 exceeded the acceptable limit. The horizontal axis shows the time, and the vertical axis shows the next 30-minute demand. NFI135 is a major NFI where cross-border flows between Bordeaux and Madrid/Barcelona air traffic control centers are transiting through. Flow demand exceeding the acceptable limit is commonly observed on NFI135. Thus, NFI135 is used to demonstrate the flow coordination result to examine the efficacy of the proposed models in re-routing the traffic flows to improve flow efficiency.

B. NFI Flow Prediction and Flow Limit Identification

C. NFI Flow Excess Re-routing

As shown in Fig. 8, the flow originating from NFI 505 (green dot) and heading to NFI279 (blue dot) is the major flow Fig. 10 shows the flow re-routing results from 1150UTC to 1250UTC on 24 Dec 2019 and 0700UTC to 0800UTC on 21 Dec 2019, where the demand exceeds the acceptable limit. For the period 1150UTC-1220UTC and 1220-1250UTC, the assigned flow route is "route 7" and "route 1" in Fig. 8 respectively, while for the period 0700UTC-0730UTC and 0730UTC-0800UTC, the assigned route is "route 6" and "route 7" in Fig. 8 respectively. Table II compares the traffic flow excess and efficiency before and after re-routing. It can be observed that the flow excess at NFI135 is avoided as the traffic flow demands every 30 minutes are below the acceptable limit after re-routing. For the 43 flights exiting the re-routing area during 1150UTC-1220UTC, there is a 9.58minute reduction in flight time, and for the 51 flights exiting the re-routing area during 1220UTC-1250UTC, the reduction in flight time is 184.87 minutes. During 1150UTC-1250UTC, there is a total flight time reduction of 194.45 minutes for the 94 flights in the re-routing area. During 0700UTC-0800UTC, It is noteworthy that the flight time reduction during the later periods (1220UTC to 1250UTC and 0730UTC to 0800UTC) significantly exceeds that observed during the earlier periods (1150UTC to 1220UTC and 0700 UTC to 0730UTC). This difference can be attributed to the cumulative effects of flow re-routing implemented in both the early and later periods, influencing air traffic throughout the later periods. Conversely, their effects on air traffic may not manifest fully during the earlier periods.

Figure 1 :

 1 Figure 1: One-day flight trajectories in French airspace: flight trajectories are aggregated as flows connecting major hubs.

  II. METHODOLOGY A. Methodology Overview The proposed air traffic flow coordination framework consists of four main steps: a) NFI identification through exploration of air traffic flow patterns; b) NFI flow feature representation and demand prediction during different periods; c) NFI flow acceptance limit determination through identifying phase transitions in flow transition duration; d) NFI flow rerouting agent design and training using reinforcement learning to avoid flow excess and improve flow efficiency. Fig. 2 presents a concept diagram of the proposed framework. Firstly, the proposed method identifies NFIs through trajectory intersection clustering. The cluster centers are the NFIs. The number of clusters is determined by graph analysis of daily air traffic flow patterns, represented by a graph whose nodes are the NFIs, and edges describing flight connectivity between the nodes. The optimal number of NFIs is determined by evaluating the daily graph pattern consistency based on the node locations and edge structures.

Figure 2 :

 2 Figure 2: Conceptual diagram of the proposed machine-learned traffic flow coordination Framework. The NFIs are identified through graph pattern analysis. Then, flow features are represented as the time and flight sequences at NFIs, used by the transformer encoder-based model to predict future demand. Comparing the prediction and the flow limit can identify the overload at NFIs. Overload can be avoided by alternative route assignment from the flow re-routing agent.

2 )

 2 |M P k ∪ M P k+1 | represents the number of node pairs in the union of M P k and M P k+1 , denoted as l in the paper. b) Structural Consistency in Flow Connectivity: Upon determining the geographical consistency in NFI locations, the next step is quantifying the structural consistency in the daily air traffic flow connectivity between the NFIs. Let Sub = M P k ∪ M P k+1 represent the set of mutually paired nodes from graph G k and graph G k+1 . C k = {c k,1 , c k,2 , ..., c k,l } denotes the nodes in Sub from graph G k , and C k+1 = {c k+1,1 , c k+1,2 , ..., c k+1,l } denotes the corresponding paired nodes from G k+1 . Let W c k and W k represent the weighted adjacency matrix of C k and V k . The entry w i,j c k represents the weight on the edge connecting nodes c k,i and c k,j . The entry w i,j k represents the weight on the edge connecting nodes v k,i and v k,j . The structural consistency of the air traffic flow patterns is measured by the mutual flow connectivity in the two graphs. More specifically, it is evaluated by the ratio of mutual flow connections between the graphs characterized by nodes C k and C k+1 compared to the union of flow connections in G k and G k+1 :

Figure 3 : 5 ) 2 )

 352 Figure3: Determination of the cluster number through saddle point identification on P (n), the polynomial approximation of the average of graph pattern consistency gc 1 and gc 2 (scaled between 0 and 1), versus the cluster number.

Figure 4 :

 4 Figure 4: Proposed transformer-based framework for the flow prediction model.For instance, multiplying x2, the embedded vector for the callsign "AFR378", by W Q , W K , W V produces the query, key, and value vectors q2, k2, and v2 associated with "AFR378". When processing the self-attention for "AFR378", its first attention weight related to element "14:15" would be S( q2×k1 √ d k ). Multiplying each value vector by the attention weight and summing up the weighted value vectors produce the output of the self-attention layer.

  denote the flight transition duration to v i for all periods t ∈ {1, 2, ..., N T } in the traffic data. Let n i = {n i,t } denote the corresponding traffic demand during t ∈ {1, 2, ..., N T }. n i,t is the total number of flights transiting to v i during t. By fitting the demand values to the transition duration values, a set of demand values {1, 2, ..., n} and the corresponding transition duration value {y 1 , y 2 , ..., y n } can be obtained. The flow acceptance limit at v i is determined by identifying the demand l i above which the transition duration of flights to v i will increase abruptly, which is formulated as: obj = arg min li (l i )var([y 1 , ..., y li ]) +(n -l i )var([y li+1 , ..., y n ])

Figure 5 :

 5 Figure 5: Graph representation for the one-month air traffic flow using 605 NFIs. The area confined by blue dashes is used for testing the flow re-routing model.

  The prediction model in this study adopts a stack of 12 transformer encoder blocks: N e = 12. The model input is the flight callsign sequences on the identified NFIs during the past one hour. The model output is the number of flights transiting through the NFIs in the future 30 minutes. The Mean Square Error (MSE) between the prediction and the truth is used to compute the model's loss function. The number of trainable parameters is at the level of 10 8 . The training batch size is 16. The learning rate is 0.00002. Air traffic data from 1st to 19th December 2019 are used for model training, the following six days of data for model testing, and the rest six days for model validation during training.

R 2

 2 values of the prediction on the three NFIs are above 0.95, showing the predicted values can reliably approximate the true demand values. A more detailed illustration of the performance of the prediction model under different prediction windows and the comparison with other state-of-the-art methods can refer to the previous paper that specifically focuses on the flow prediction part of the framework[START_REF] Chunyao | Air traffic flow representation and prediction using transformer in flow-centric airspace[END_REF].Calculating the flight transition duration to the NFIs during different periods and under different traffic flow demands and observing the point of demand above which the flight transition durations increase sharply identifies the flow acceptance limit at NFIs. Fig.6shows flight transition duration versus the flow demand on four example NFIs. The blue circles show the original observations from the traffic data, the solid red lines show the fitted curves of the observations, and the pink dashes bound the 95% confidence intervals of the fitting. The curves of duration-demand are fitted using third-degree polynomials. The solid black line indicates the identified acceptance limit of the NFIs. Such phase transitions are observed on 68% of the NFIs, while the transition durations to the rest of the NFIs show no explicit trend concerning the flow demand based on the observations from the one-month data. The reason may be that traffic flow demand on these NFIs is below capacity during this month, so there are no observations for their overloaded circumstances. The flow acceptance limits on such NFIs are set as the maximum flow demand observed this month.

Figure 6 :

 6 Figure 6: Flight transition duration (y-axis) versus the flow demand (x-axis) on four NFIs. The black vertical lines show the acceptable flow limit of the NFIs.

Figure 7 :

 7 Figure 7: Instances of flow demand exceeding the acceptable limit are commonly observed on NFI135 during Dec 2019, where flow re-routing is to be applied to avoid the overload.transiting through NFI135 (red dot). Fig.8also shows seven alternative routes that this flow may take, derived from the historical flight trajectories. The reinforcement learning-based rerouting agent is trained in the environment of the Bluesky ATC simulator[START_REF] Hoekstra | Bluesky atc simulator project: an open data and open source approach[END_REF]. Given the limited computational resources, the part of French airspace confined by [42N, 45.5N, 2W, 4E], shown by blue dashes in Fig.5, is chosen as the re-routing area inside which the agent's actions is generated and rewarded. The original flight trajectory is used as the planned flight path. Every 30 minutes, based on the traffic flow state, i.e., traffic demand at the NFIs during the past 30 minutes, the agent

Figure 8 :

 8 Figure 8: Alternative routes for air traffic flow transiting through NFI135. The model is trained with the traffic scenarios during Dec 2019 when the demands exceeded the flow limit at NFI135. To consider the accumulated effects of the agent's re-routing decisions, the length of each training episode is set as 2 hours, which can also cover the travel time of most flights in the focal airspace. As shown in Fig. 9, the model converges to an optimal route assignment policy as the number of training steps increases. The mean reward of each episode fluctuates and continues to increase at the beginning and stabilizes at the value of 322.168 as the training step increases. It shows that, on average, the travel time of each flight in the rerouting area is reduced by 322.168 seconds during the 2hour re-routing episode compared with the originally planned flows. Moreover, the explained variance of the value function increases towards one during training, showing the learned value function accurately models the environment's dynamics, and the predictions accurately match the true returns.Fig.10shows the flow re-routing results from 1150UTC to 1250UTC on 24 Dec 2019 and 0700UTC to 0800UTC on 21 Dec 2019, where the demand exceeds the acceptable limit. For the period 1150UTC-1220UTC and 1220-1250UTC, the assigned flow route is "route 7" and "route 1" in Fig.8respectively, while for the period 0700UTC-0730UTC and 0730UTC-0800UTC, the assigned route is "route 6" and "route 7" in Fig.8respectively. Table II compares the traffic flow excess and efficiency before and after re-routing. It can be observed that the flow excess at NFI135 is avoided as the traffic flow demands every 30 minutes are below the acceptable limit after re-routing. For the 43 flights exiting the re-routing area during 1150UTC-1220UTC, there is a 9.58minute reduction in flight time, and for the 51 flights exiting the re-routing area during 1220UTC-1250UTC, the reduction in flight time is 184.87 minutes. During 1150UTC-1250UTC, there is a total flight time reduction of 194.45 minutes for the 94 flights in the re-routing area. During 0700UTC-0800UTC,

Figure 9 :

 9 Figure 9: The re-routing model converges to the optimal policy during training.

Figure 10 :

 10 Figure 10: Re-routing results for air traffic flow transiting through NFI135.

  IV. CONCLUSIONSAiming to contribute to the future flow-centric paradigm, this paper proposes a flow-centric air traffic flow coordination framework to avoid overload at major flow intersections and improve efficiency. First, Nominal Flow Intersections (NFIs)

  were identified through modeling and analyzing the consistency of daily air traffic flow patterns. Then, a text-enriched description of the flow and movement of flights at the NFIs was used by a transformer encoder-based neural network to learn the flow relations to predict future demands at the NFIs. The acceptable flow limits at NFIs were determined as the demand above which the flow efficiency, characterized by flight duration between NFIs, changed significantly. Finally, a reinforcement learning-based flow re-routing model was trained to assign alternative routes to air traffic dynamically to improve flow efficiency and avoid the foreseen overload. The re-routing agent is tested on the NFI135 handling crossborder flows between Bordeaux and Madrid/Barcelona air traffic control centers, using ADS-B data in Dec 2019. After re-routing, air traffic in the airspace speeds up, and flight time is reduced compared with using original flight routes. By integrating accurate flow prediction and dynamic re-routing capabilities into airspace management systems, flow-centric airspace can proactively optimize traffic flow, reduce delays, and enhance safety. This research and findings may contribute to the development of concepts of operations for flow-centric airspace.ACKNOWLEDGMENT This research is supported by the National Research Foundation, Singapore, and the Civil Aviation Authority of Singapore, under the Aviation Transformation Programme. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views of National Research Foundation, Singapore and the Civil Aviation Authority of Singapore.

TABLE I :

 I PREDICTION PERFORMANCE ON THREE MAJOR NFIS OVER THE PARIS (NFI460), GENEVA (NFI365), AND BORDEAUX (NFI135) AREA CONTROL CENTERS IN TERMS OF MAE, MSE, MAPE, AND R 2 . Bordeaux and Madrid/Barcelona control centers. Four metrics are used for the performance evaluation: Mean Absolute Error (MAE), Mean Squared Error (MSE), Mean Absolute Percentage Error (MAPE), and R-squared (R 2 ). It can be observed that the prediction can accurately capture flow demand changes with MAE values smaller than one and MAPE values smaller than 0.15 for NFI365 and NFI135. The

	NFI	MAE MSE	MAPE R 2
	1: NFI460 2.339 12.184 0.217	0.951
	2: NFI365 0.800 2.702	0.127	0.978
	3: NFI135 0.984 2.929	0.145	0.977
	border of

TABLE II :

 II COMPARISON OF FLIGHT TRAVEL TIME REDUCTION OF TRAF-FIC FLOWS BEFORE & AFTER RE-ROUTING.

	time (UTC)	before		after		RTE time reduction
		flow XS. flow XS.		mins/flights
	11:50-12:20 (24th)	18	0	9	0	7	9.583/43
	12:00-12:30 (24th)	21	3	13	0	\	133.083/51
	12:10-12:40 (24th)	23	5	13	0	\	174.200/51
	12:20-12:50 (24th)	20	2	12	0	1	184.867/51
	7:00-7:30 (21st)	18	0	13	0	6	16.167/32
	7:10-7:40 (21st)	21	3	11	0	\	51.683/54
	7:20-7:50 (21st)	20	2	11	0	\	150.583/59
	7:30-8:00 (21st)	24	6	10	0	7	250.167/68


Laura Kergoat
email: laura.kergoat@inrae.fr

Aymeric Dabrin

Matthieu Masson

Thibault Datry

Chloé Bonnineau

COPPER ALTERS THE MICROBIAL RESPONSES IN SEDIMENT DURING DROUGHT

Does copper decrease the resistance/resilience of microbial communities in sediment after drying ?  Microbial communities in sediment ensures important ecological functions (organic matter mineralization, pollutant degradation). (Krause et al., 2017 ; Gandy et al., 2007)  Hydrological drought caused in part by anthropogenic activities (irrigation, climate change) are unpredictible events with water deficit intense in space and time. (Sarremejane et al., 2022)  In an agricultural watershed context, droughts happened in areas exposed to pesticide contamination such as copper used as a fungicide in both conventional and organic farming  Negative effect of drought and copper have already been demonstrated but the combined effect of these two stressors have never been shown. (Fromin et al.,



  Microbial response to Copper+Drying combined is not a simple addition (one stressor driving or antagonistic effects) After rewetting, copper reduces the microbial community resilience endangering the river ecosystem functioning To go further and detect interaction, copper and drought gradients should be tested on microbial functions In sediment (Surface, Medium, Deep) -Exo-enzymatic activities (fluorometry) -Denitrification (gas chromatography)  In the air above the entire column -Respiration (Picarro instrument) For all the graphics, different letters indicate significant different results between treatments (One-way ANOVA, p<0,05) Fromin et al., 2010 doi: 10.1002/eco.115 ; Gandy et al., 2007 doi: 10.1016/j.scitotenv.2006.11.004 ; Krause et al., 2017 doi: 10.1002/2016WR019516 ; Mahamoud Ahmed et al., 2018 doi: 10.3389/fmicb.2018.01852 ; Sarremejane et al., 2022 doi: 10.1002/eco.2390


Alice Le Brigant
email: alice.le-brigant@univ-paris1.fr

Stephen C Preston
email: stephen.preston@brooklyn.cuny.edu

CONJUGATE POINTS ALONG KOLMOGOROV FLOWS ON THE TORUS

The geodesics in the group of volume-preserving diffeomorphisms (volumorphisms) of a manifold M , for a Riemannian metric defined by the kinetic energy, can be used to model the movement of ideal fluids in that manifold. The existence of conjugate points along such geodesics reveal that these cease to be infinitesimally length-minimizing between their endpoints. In this work, we focus on the case of the torus M = T 2 and on geodesics corresponding to steady solutions of the Euler equation generated by stream functions ψ = -cos(mx) cos(ny) for integers m and n, called Kolmogorov flows. We show the existence of conjugate points along these geodesics for all pairs of strictly positive integers (m, n), thereby completing the characterization of all pairs (m, n) such that the associated Kolmogorov flow generates a geodesic with conjugate points.

1. Introduction 1.1. Motivation. Since Arnold in 1966 [START_REF] Arnold | On the differential geometry of infinite-dimensional Lie groups and its application to the hydrodynamics of perfect fluids[END_REF], the Euler equations for ideal fluids have had a well-known geometric interpretation as geodesics on the group Diff µ (M ) of volume-preserving diffeomorphisms, or volumorphisms, of a manifold M , under a right-invariant Riemannian metric defined by the kinetic energy. Since this approach can be made rigorous as in Ebin-Marsden [START_REF] Ebin | Groups of diffeomorphisms and the motion of an incompressible fluid[END_REF] and establishes that the geodesic equation is actually a smooth ODE on the group of Sobolev H s volumorphisms on M for s > 1 2 dim(M ) + 1, it can be shown that the Riemannian exponential map is C ∞ , and is invertible near zero. This shows that volumorphisms sufficiently close in H s can be joined to the identity by a unique minimizing geodesic. Since then similar geometric interpretations have been found for a variety of other PDEs of continuum mechanics; see Arnold-Khesin [START_REF] Arnold | Topological methods in hydrodynamics[END_REF] for a survey and Misio lek-Preston [START_REF] Misio | Fredholm properties of Riemannian exponential maps on diffeomorphism groups[END_REF] for an overview of the ODE approach.

On a sufficiently long time interval, a geodesic may cease to minimize the length between its endpoints (or equivalently, the energy for a constant-speed parameterization). When this happens, the endpoint is called a cut point, and if in addition the geodesic ceases to be even infinitesimally minimizing between its endpoints, it is called a conjugate point. The existence of conjugate points is intimately connected to the existence of positive-curvature sections along a geodesic. On the volumorphism group this existence was unclear since curvature computations have seemed to suggest Diff µ (M ) has mostly negative curvature. The question of whether conjugate points exist on Diff µ (M ) was posed already by Arnold [START_REF] Arnold | On the differential geometry of infinite-dimensional Lie groups and its application to the hydrodynamics of perfect fluids[END_REF] in 1966, but was not solved until Misio lek [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF] in 1993 found them along rotations on the space Diff µ (M ) of volumorphisms of the 2-and 3-spheres M = S 2 and M = S 3 .

The first example on a flat M was also found by Misio lek [START_REF] Misio | Conjugate points in Dµ(T 2 )[END_REF], who showed that along the Kolmogorov flow on the torus with stream function ψ = -cos 6x cos 2y, there is eventually a conjugate point. To do so he devised what is now known as the Misio lek criterion, giving a sufficient condition for existence of a conjugate point. This criterion has been successfully used to find conjugate points along other steady flows on other manifolds, such as the 2D ellipsoid [START_REF] Tauchi | Existence of a conjugate point in the incompressible Euler flow on an ellipsoid[END_REF], the 3D ellipsoid [START_REF] Lichtenfelz | Existence of a conjugate point in the incompressible Euler flow on a ellipsoid[END_REF], and the sphere [START_REF] Benn | Conjugate points in Dµ(S 2 )[END_REF]. In three dimensions, conjugate points have a substantially different nature [START_REF] Ebin | Singularities of the exponential map on the volume-preserving diffeomorphism group[END_REF] and are much more common than in two dimensions, and can be found using a necessary and sufficient local criterion along any particle path [START_REF] Preston | On the volumorphism group, the first conjugate point is always the hardest[END_REF][START_REF] Preston | The WKB method for conjugate points in the volumorphism group[END_REF], as shown by the second author. This technique was used in [START_REF] Preston | The geometry of axisymmetric ideal fluid flows with swirl[END_REF] to find conjugate points along axisymmetric 3D flows.

It is known that the Misio lek criterion cannot capture all conjugate points. For example Tauchi-Yoneda [START_REF] Tauchi | Existence of a conjugate point in the incompressible Euler flow on an ellipsoid[END_REF] observed that the Misio lek index never detects conjugate points along the spherical rotation, in spite of the fact that infinitely many of them exist. More generally Tauchi-Yoneda [START_REF] Tauchi | Arnold stability and Misio lek curvature[END_REF] showed that the Misio lek criterion cannot detect conjugate points along an Arnold stable flow. A condition of the second author [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF] is more suitable for detecting conjugate points in such cases (and particularly along rotational flows).

Specifically on the torus, the problem of finding conjugate points along Kolmogorov flows ψ = -cos mx cos ny for all positive integer pairs (m, n) was posed by Drivas et al. [START_REF] Drivas | Conjugate and cut points in ideal fluid motion[END_REF], generalizing from Misio lek's original example of m = 6, n = 2. (We will assume throughout that m ≥ n without loss of generality due to symmetry.) They found many additional examples using Misio lek's criterion, for pairs (m, n) satisfying a condition equivalent to n ≥ 2 and m > 3n 2 +6 √ 3n . Soon after the second author [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF] found via brute force search several more examples, including (m, 1) for m ≥ 2, along with (2, 2), [START_REF] Benn | Conjugate points in Dµ(S 2 )[END_REF][START_REF] Arnold | Topological methods in hydrodynamics[END_REF], and [START_REF] Benn | Conjugate points in Dµ(S 2 )[END_REF][START_REF] Benn | Conjugate points in Dµ(S 2 )[END_REF], and conjectured that conjugate points existed for all (m, n) except possibly for (1, 1). 1.2. Contributions. In this paper we prove this conjecture, that conjugate points exist for all strictly positive integer pairs (m, n), giving an explicit form of a test function for m > n, and a slightly different test function for m = n ≥ 2, along with an example due to Drivas in the case m = n = 1. In the case n = 0 with m > 0, it is known [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF] that there are no conjugate points. Hence we have obtained a complete characterization of all pairs (m, n) ⊂ Z 2 ≥0 such that the Kolmogorov flow ψ = -cos mx cos ny generates a geodesic in the volume-preserving diffeomorphism group that infinitesimally minimizes length between its endpoints for all time: if m and n are both nonzero, the geodesic cannot minimize; if m or n is zero then the geodesic infinitesimally minimizes. 1.3. Outline. In Section 2 we describe the geometric approach to the Euler equation of ideal fluids as originally pioneered by Arnold. We also review Kolmogorov flows, which may be generally defined as those steady solutions of the Euler equation with stream function ψ satisfying ∆ψ = -λ 2 ψ, and specifically on the torus of the form ψ = -cos mx cos ny. Finally we recall the index form for detecting conjugate points and the Misio lek criterion which greatly simplifies this computation. In Section 3 we prove Theorem 2 (the m > n case), Theorem 3 (the m = n ≥ 2 case) and Theorem 4 (the m = n = 1 case), showing the existence of conjugate points. In Section 4 we discuss the relationship between the Misio lek index and the Rayleigh quotient restricted to a closed subspace. Section 5 describes the algorithm we devised to exploit this relationship in order to construct an optimal variation field numerically; those given in Section 3 are simply truncated versions of these that still work. In this section we explain why all of the variation test functions we find have the same basic shape: a perturbation of the simplest Laplacian eigenfunction cos x or sin x. Finally in Section 6 we suggest some other problems about conjugate points that one can tackle using the same methods as those presented here. In an Appendix we present a simple formula that is useful numerically when writing the Poisson bracket in a Fourier basis, necessary for using standard matrix algorithms to optimize the Misio lek index.

1.4. Computations. Symbolic computations were performed in Maple 2021, while numerical computations were performed using Python. Maple code for computing the index form as in Theorems 2-4, and Python code for numerically minimizing the index form and obtaining the form of the candidate minimizers, are both available on github: https://github.com/alebrigant/conjugate-points.

1.5. Acknowledgements. The authors acknowledge support of the Institut Henri Poincaré (IHP, UAR 839 CNRS-Sorbonne Université), and LabEx CARMIN (ANR-10-LABX-59-01). This work was done while the second author visited the first author at IHP for the Geometry and Statistics in Data Sciences (GESDA) thematic quarter, funded by the French National Center for Scientific Research (CNRS). Both authors thank the IHP for their hospitality. We also thank Theo Drivas and Alexander Shnirelman for very helpful discussions and suggestions.

Background

2.1. Volumorphisms and the Euler-Arnold equation. Suppose M is a 2dimensional manifold equipped with a Riemannian metric •, • inducing an area form µ. We are interested in diffeomorphisms ϕ : M → M that preserve the area form: ϕ * µ = µ, called volume-preserving diffeomorphisms of M or volumorphisms for short. The space Diff µ (M ) of volumorphisms of M is (formally) a submanifold of the space of diffeomorphisms Diff(M ). The tangent vectors at ϕ ∈ Diff µ (M ) are right translations X • ϕ of divergence-free vector fields X of M tangent to the boundary. We refer the interested reader to [START_REF] Ebin | Groups of diffeomorphisms and the motion of an incompressible fluid[END_REF] for more details on these manifold structures in the context of Sobolev spaces.

A volumorphism ϕ can be seen as describing the positions, at a given time, of the particles of an ideal fluid (incompressible and inviscid) moving inside of M : the value ϕ(x) gives the position of the particle that was at position x ∈ M at t = 0. The volume-preserving property of ϕ is a consequence of the incompressibility of the fluid. With this interpretation, a tangent vector X • ϕ is the velocity field of the fluid at that time, and the kinetic energy of the fluid defines a Riemannian metric on Diff µ (M )

(1) g ϕ (X • ϕ, Y • ϕ) = M X • ϕ, Y • ϕ µ.
Since for any volume-preserving ϕ ∈ Diff µ (M ), the pullback measure is ϕ * µ = µ, we see that

g ϕ (X • ϕ, Y • ϕ) = M X, Y µ = g id (X, Y ),
and so the kinetic metric is right-invariant on Diff µ (M ). The result shown by Arnold [START_REF] Arnold | On the differential geometry of infinite-dimensional Lie groups and its application to the hydrodynamics of perfect fluids[END_REF] is the following: the geodesics in Diff µ (M ) for the kinetic metric (1) describe the motion of an ideal fluid in M . To see this, we can consider the geodesic equation on the larger space Diff(M ), which is simply ∂ 2 γ/∂t 2 = 0, and then orthogonally project it on the submanifold Diff µ (M ). In terms of the velocity field X(t, •) associated with the geodesic γ(t, •)

∂γ ∂t (t, x) = X(t, γ(t, x)),
this gives

P ∂ 2 γ ∂t 2 = P ∂X ∂t + ∇ X X • γ = 0,
where ∇ is the Levi-Civita connection associated to the Riemannian metric on M . The orthogonal projection P on the space of divergence-free vector fields tangent to the boundary is obtained by the Hodge decomposition of vector fields

(2)

P (X) = X -grad f,
where f is a function verifying ∆f = div X and whose normal component along the boundary equals that of the vector field X. The geodesic equation for the kinetic metric (1) on Diff µ (M ) is therefore given by

(3) ∂X ∂t + ∇ X X = -grad p,
where p is called the pressure function, defined up to a constant by ∆p = -div(∇ X X) and its component along the normal ν: grad p, ν = -X, ν . Equation ( 3) is the Euler equation for incompressible fluids, also called the Euler-Arnold equation in virtue of its interpretation by Arnold as a geodesic equation on the space of volumorphisms.

2.2. Kolmogorov flows on the torus. In the two-dimensional case, it can be shown that the Euler-Arnold equation ( 3) can be rewritten in terms of the curl as

∂ ∂t (curl X) + X(curl X) = 0.
Now divergence-free vector fields on surfaces can be written as X = sgrad ψ + W where ψ : M → R is defined by ∆ψ = curl X and ψ| ∂M = 0, while the vector field W satisfies curl W = 0 and div W = 0. For such a vector field X, the Euler-Arnold equation becomes

∂ ∂t ∆ψ + {ψ, ∆ψ} + W (∆ψ) = 0.
Here {•, •} is the Poisson bracket, which can be defined by the following property of Hamiltonian vector fields: if sgrad f and sgrad g are Hamiltonian vector fields generated by mean-zero functions f, g : M → R, then {f, g} is the unique mean-zero function generating their Lie bracket

[sgrad f, sgrad g] = sgrad{f, g}.
In dimension 2, the Poisson bracket is simply given by {f, g}µ = df ∧ dg, which reduces to {f, g} = ∂ x f ∂ y g -∂ y f ∂ x g on the torus. We then see that any function ψ for which (4) {ψ, ∆ψ} = 0, generates a steady solution X = sgrad ψ of the Euler-Arnold equation, i.e., a geodesic γ with associated velocity field X satisfying ∂X ∂t = 0. The function ψ is called the stream function. Here we will focus on geodesics on Diff µ (T 2 ) whose velocity field is generated by the stream functions ( 5)

ψ(x, y) = -cos(mx) cos(ny), (m, n) ∈ N 2 .
These are eigenfunctions of the Laplacian for the eigenvalues -λ 2 = -(m 2 + n 2 ), and they satisfy (4) to give steady solutions of the Euler equations, the so-called Kolmogorov flows on the torus M = T 2 . For these flows, the fluid has a circular movement inside rectangular grid cells, as shown in Figure 1. The question that we ask is: can we find conjugate points along any Kolmogorov flow on the torus?

The existence of conjugate points is related to the question of the uniqueness of the geodesic between its endpoints, and whether a perturbation of the initial condition of the corresponding flow can lead to the same result as no perturbation at all.

2.3.

Conjugate points and the Misio lek criterion. Two volumorphisms ϕ 1 and ϕ 2 are conjugate if there exists a family of geodesics that start at ϕ 1 and end at ϕ 2 up to first order; more precisely, if there exists γ(s, t) with (s, t)

∈ (-, ) × [0, 1] such that for all s, γ(s, •) is a geodesic in Diff µ (M ), γ(s, 0) = ϕ 1 , γ(s, 1) = ϕ 2
, and the corresponding Jacobi field J(t) = ∂γ ∂s (0, t) satisfies J(0) = J(1) = 0. We say that there is a conjugate point along a geodesic γ if there exists a time T > 0 such that γ(0) and γ(T ) are conjugate points. Conjugate points are particular cases of cut points. However, they can only happen in manifolds with some positive curvature, while cut points can occur even in flat manifolds due to non trivial topology.

Conjugate points can be detected using the so-called index form, defined for any vector field Y (t) along the geodesic γ(t) by

I(Y, Y ) = T 0 DY dt , DY dt -R(Y, γ) γ, Y dt.
It can be shown [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF] that

I(Y, Y ) = 0 for some Y satisfying Y (0) = Y (T ) = 0 if and only if Y is a Jacobi field. Furthermore if I(Y, Y ) < 0 for some such field Y , then
there is a Jacobi field J along γ vanishing at t = 0 and for t = τ for some τ < T . Establishing negativity of the index form for some vector field Y is thus an effective way to show there must be a conjugate point without actually having to find it. In our setting of volume-preserving diffeomorphisms on a surface, the index form can be written in the following way.

Proposition 1 (Corollary 9 in [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF]). Suppose M is a surface and that X = sgrad ψ solves the Euler equation (3) on [0, T ]. For a time-dependent family of functions g(t) on M , vanishing on the boundary of M and at t = 0 and t = T , the index form for Y = sgrad g along the geodesic with velocity field X becomes

(6) I(Y, Y ) = T 0 M | grad h| 2 + ∆ψ{g, h} dµdt, h := ∂g ∂t + {ψ, g}.
The Misio lek criterion consists of computing the index form for a particular family of deformations, namely g(t) = sin( πt T )f for some function f : M → R, and then letting time T go to infinity. If the result is negative for some function f , then for sufficiently large T the index form will be negative and hence there will be a conjugate point occurring at some time τ < T .

Consider a steady geodesic in Diff µ (M ) generated by a stream function ψ : M → R satisfying ∆ψ = Υ • ψ for some function Υ : R → R; every steady solution of (4) can at least locally be expressed in this form. The index form (6) at a vector field Y = sgrad g generated by a family of deformations g(t) = sin( πt T )f for a given function f can be written in terms of the Poisson bracket φ = {ψ, f } as

I(g) = π 2 2T M | grad f | 2 dµ + T 2 M | grad φ| 2 + Υ (ψ) φ 2 dµ.
The Misio lek criterion states that there is a conjugate point eventually occurring along the geodesic if lim

T →∞ 2 T I(g) < 0 for some function f , i.e., M I(φ) := M |grad φ| 2 + Υ (ψ) φ 2 dµ < 0, φ = {ψ, f }.
In the present case, for M = T 2 , we will take

(7) ψ = -cos mx cos ny, with Υ(ψ) = -λ 2 ψ, where λ 2 = m 2 + n 2 .
This leads to the explicit version of the Misio lek index in the form

(8) M I(φ) := M |grad φ| 2 -(m 2 + n 2 )φ 2 dµ < 0, φ = {ψ, f },
which is the version we will use in the remainder of the paper. We will show that when φ = -cos mx cos ny for positive integers m, n, there is always a function f to make M I(φ) < 0, leading to both existence of conjugate points and the indefiniteness of Arnold's quadratic form.

2.4. Remark on the link with stabillity. Note that the Misio lek index is precisely the second variation of the Dirichlet energy M |grad ψ| 2 dµ, which is the L 2 energy of a vector field, when considered under the action by volume-preserving diffeomorphisms of a 2D manifold. Here variations are generated by test functions f , while the actually varied stream function under the action is given by φ = {ψ, f }.

The formula for the second variation of the energy is given in Chapter II, Remark 2.4 of Arnold-Khesin [START_REF] Arnold | Topological methods in hydrodynamics[END_REF], and it agrees with the Misio lek index. Hence the second variation of energy of this minimization problem is indefinite if and only if there is a conjugate point that is detectable by the Misio lek criterion. As pointed out by Tauchi-Yoneda [START_REF] Tauchi | Arnold stability and Misio lek curvature[END_REF], this implies that conjugate points along an Arnold-stable flow cannot be detected by the Misio lek criterion, although they may exist (as happens for many rotational flows, such as on the sphere; see [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF] for examples). Indefiniteness of the second variation does not necessarily imply instability of the steady flow, however, and thus existence of conjugate points detectable by the Misio lek criterion does not necessarily imply that the linearized Euler equation has unstable growing solutions. See Arnold-Khesin, Chapter II, Remark 3.7 [START_REF] Arnold | Topological methods in hydrodynamics[END_REF]. It would be very interesting to further clarify the apparent connection between instability of the linearized Euler equation and the existence of conjugate points along the corresponding geodesic, particularly since the Arnold criterion for stability and the Misio lek criterion for conjugate points are both only sufficient but not necessary, and since the more intuitive connection would be between conjugate points and stability.

Note that the second variation formula (and thus the Misio lek criterion) apply to finding infinitesimal energy-decreasing variations of the energy by diffeomorphisms, in terms of the local analysis of vector fields at the identity. The problem of finding the global minimizer was studied by Shnirelman and Lan [START_REF] Lan | The minimizer of the Dirichlet integral[END_REF], who numerically found a minimizer that is not smooth in the case (m, n) = (1, 1). It is possible that the roughness of this optimizer is related to some difficulty in the second variation we are studying here, but the connection is not completely clear.

Main results

We will show that for every pair (m, n) of positive integers with m ≥ n, the Misio lek criterion (8) detects a conjugate point along the geodesic in Diff µ (T 2 ). It is easier to handle the three cases m > n, m = n ≥ 2, and m = n = 1 separately.

The formulas presented in the next three theorems were discovered by numerics we will elaborate on in Section 5. They are finite truncations of infinite Fourier series, dominated by the lowest-order term (either cos x or sin x), which show up because they are the eigenfunctions of the Laplacian with smallest eigenvalue -1.

Theorem 2. If m and n are positive integers with m > n, and ψ(x, y) = -cos mx cos ny, then there is a function f : T 2 → R such that φ = {ψ, f } satisfies the Misio lek criterion [START_REF] Lichtenfelz | Existence of a conjugate point in the incompressible Euler flow on a ellipsoid[END_REF]. Hence if the geodesic γ in Diff µ (T 2 ) has initial conditions γ(0) = id and γ (0) = sgrad ψ, then γ(T ) is conjugate to γ(0) for some T > 0, and γ is not minimizing beyond T .

Proof. See Maple file in https://github.com/alebrigant/conjugate-points for details of the computations here.

For numbers a, b ∈ R, define

f (x, y) = cos x 1 + a cos (2mx) + b cos (2ny) .

We compute the Poisson bracket φ(x, y) = ψ x (x, y)f y (x, y) -ψ y (x, y)f x (x, y)

= 2mn cos x a cos(mx) sin(ny) sin(2mx) -b sin(mx) cos(ny) sin(2ny)

+ n cos(mx) sin(ny) sin x 1 + a cos(2mx) + b cos(2ny) .

The quantity M I from ( 8) is then given by

M I = π 2 n 2 4 H(a, b, m, n),
where

(10) H(a, b, m, n) = 16a 2 m 4 + 16b 2 m 2 n 2 + 24a 2 m 2 -12abm 2 + 4b 2 m 2 + 4b 2 n 2 + 8am 2 -8bm 2 + a 2 -ab + b 2 + 2a -2b + 2.
This is quadratic in both a and b, with positive leading-order coefficients, and thus the critical point must be a global minimizer. However the formula is slightly complicated, and it is easier to use the points ( 11)

a 0 = - 4m 2 + 1 16m 4 + 24m 2 + 1 and b 0 = 1 4n 2 + 1 ;
here a 0 is the minimizer of H(a, 0, m, n) and b 0 is the minimizer of H(0, b, m, n).

With these choices we obtain

H(a 0 , b 0 , m, n) = J(m, n) (16m 4 + 24m 2 + 1)(4n 2 + 1)
, where

J(m, n) = 4n 2 (16m 4 + 40m 2 + 1) -64m 6 -48m 4 + 28m 2 + 1. (12)
We want to show that J(m, n) < 0 for every choice of naturals m and n with m > n.

Clearly J is increasing as a function of n, so we have

J(m, n) ≤ J(m, m -1) = -128m 5 + 176m 4 -320m 3 + 192m 2 -8m + 5 for m ≥ 2. Writing m = k + 1 for k ≥ 1 this becomes J(m, n) ≤ J(k + 1, k) = -128k 5 -464k 4 -896k 3 -992k 2 -520k -83,
which is obviously negative. We conclude that J(m, n) < 0 for all positive integers m, n with m > n. This ensures that H(a 0 , b 0 , m, n) < 0, and thus that M I < 0 for these choices of a 0 and b 0 . By the Misio lek criterion, there is eventually a conjugate point along the corresponding geodesic.

Observe in the proof above that the worst-case scenario is when n = m-1: when n is smaller than m -1 the index form is even more negative. This corresponds to conjugate points being easier to find when m and n are farther apart, as was found by Drivas et al. [START_REF] Drivas | Conjugate and cut points in ideal fluid motion[END_REF]. It thus stands to reason that the hardest case is when n = m. Indeed we need two extra terms in our formula for the variation to make it work in this case. Note also that we require that n ≥ 2; this variation field does not work if m = n = 1. We will handle that case separately. Theorem 3. Suppose n ≥ 2 and ψ(x, y) = -cos (nx) cos (ny). Then there is a function f : T 2 → R such that φ = {ψ, f } satisfies (8). Thus there is a time T > 0 such that γ(T ) is conjugate to γ(0) along the corresponding geodesic γ, and γ is not minimizing past T .

Proof. See Maple file in https://github.com/alebrigant/conjugate-points for details of the computations here.

For coefficients a, b, c, d ∈ R, define

M I = π 2 n 2 4 H(a, b, c, d, n), where (15) H(a, b, c, d, n) = 16a 2 n 4 + 64abn 4 + 256b 2 n 4 + 16c 2 n 4 + 16d 2 n 4 + 8adn 3 -64cdn 3 + 8a 2 n 2 -4abn 2 -12acn 2 + 32b 2 n 2 + 24c 2 n 2 + 24d 2 n 2 + 6adn -8an 2 -16cdn + 8cn 2 + a 2 -ab -ac + b 2 + c 2 + d 2 -8dn -2a + 2c + 2.
Again, this is quadratic in the four unknown coefficients a, b, c, d with positive leading-order coefficients, so the unique critical point is a global minimum. In this case there is no simpler formula which works, so we find the critical point by solving the linear system

∂H ∂a = ∂H ∂b = ∂H ∂c = ∂H ∂d = 0.
Explicitly this system looks like

(32a + 64b)n 4 + 8dn 3 + (16a -4b -12c -8)n 2 + 6dn + 2a -b -c -2 = 0, (64a + 512b)n 4 + (-4a + 64b)n 2 -a + 2b = 0, 32cn 4 -64dn 3 + (-12a + 48c + 8)n 2 -16dn -a + 2c + 2 = 0, 32dn 4 + (8a -64c)n 3 + 48dn 2 + (6a -16c -8)n + 2d = 0.
There is a unique solution of the system above, given by a 0 = (8n 2 + 1)(256n .

In spite of how complicated these formulas appear, the value of H at the critical point is relatively simple in n: we get ( 16) H(a 0 , b 0 , c 0 , d 0 , n) = -4096n 8 + 3584n 6 + 1008n 4 + 68n 2 + 1 12288n 8 + 9728n 6 + 1840n 4 + 116n 2 + 2 .

Replacing n with √ 4 + k for some k ≥ 0 (since n ≥ 2 by assumption), this becomes

H(a 0 , b 0 , c 0 , d 0 , √ 4 + k) = -4096k 4 -61952k 3 -349200k 2 -868412k -802799 12288k 4 + 206336k 3 + 1298224k 2 + 3627508k + 3798226 ,
and this is obviously negative for all real k ≥ 0. We conclude that with n ≥ 2, the choice (a, b, c, d) = (a 0 , b 0 , c 0 , d 0 ) leads to a variation f given by ( 13) such that φ = {ψ, f } satisfies the Misio lek criterion [START_REF] Lichtenfelz | Existence of a conjugate point in the incompressible Euler flow on a ellipsoid[END_REF]. Hence there is eventually a conjugate point along the geodesic.

Example 1. In the case m = 3 and n = 2, the formula ( 9) with a and b given by ( 11) becomes [START_REF] Tauchi | Existence of a conjugate point in the incompressible Euler flow on an ellipsoid[END_REF] f (x, y) = cos x 1 -37 1513 cos 6x + 1 17 cos 4y .

This is plotted on the left side of Figure 2 below. On the other hand if m = n = 2, the variation f from ( 13) becomes The plot is shown on the right side of Figure 2. Note that especially in the offdiagonal case, both the formula and graph are substantially simpler than the one given in the second author's paper [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF]. Also note that the graphs are basically indistinguishable to the naked eye, although none of the terms in (18) can be omitted without changing the sign of the index.

Variation fields of the form considered in Theorems 2 and 3, dominated as they are by the cos x term, do not seem suitable for making the Misio lek index negative in the case m = n = 1. Instead the following field dominated by sin x, found by Theodore Drivas (personal communication) and included here with his permission, gives an energy-reducing variation in that case. Theorem 4. Suppose ψ(x, y) = -cos x cos y. Then there is a function f : T 2 → R such that φ = {ψ, f } satisfies (8). Thus there is a time T > 0 such that γ(T ) is conjugate to γ(0) along the geodesic γ, and γ is not minimizing past T .

Proof. Define f (x, y) = sin x + 1 10 sin (x + 2y) -1 20 sin (3x) + 1 100 sin (5x). We easily compute that the Poisson bracket φ = {ψ, f } is φ(x, y) = 1 5 sin x cos y cos (x + 2y) -cos x sin y cos x + 1 10 cos (x + 2y) -3 20 cos (3x) + 1 20 cos (5x) . The Misio lek index [START_REF] Lichtenfelz | Existence of a conjugate point in the incompressible Euler flow on a ellipsoid[END_REF] is given by

M I = T 2 |∇φ| 2 -2φ 2 dA = - 3π 2 200 , Figure 2.
On the left, the optimal perturbation f of the form [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF] in the case m = 3 and n = 2, given explicitly by [START_REF] Tauchi | Existence of a conjugate point in the incompressible Euler flow on an ellipsoid[END_REF]. On the right, the optimal perturbation f of the form ( 13), given explicitly by ( 18), in the case m = n = 2. These perturbations generate a variation φ satisfying the Misio lek criterion (8) for conjugate points. Note that they appear quite similar at this level. and we are done.

In Figure 3 we show the local deformation of the stream function ψ induced by the area-preserving diffeomorphism generated by the test function f , in an example of each of the three cases discussed above. This can be compared to Figure 1 which shows ψ without deformation. The second author [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF] showed that any test function solely supported in a single cell of a Kolmogorov flow on T 2 cannot generate an energy-reducing variation, even if time variation is allowed. Hence one intuitively expects that the rectangular gridlines must be substantially deformed by any energy-reducing variation, which is indeed what we see. Unsurprisingly the most deformation happens at the vertices where gridlines intersect, since there the fluid has hyperbolic fixed points and the fluid is highly susceptible to small perturbations.

The minimization principle

In this section we relate the Misiolek criterion to an eigenfunction problem on the Poisson bracket φ = {ψ, f }. Recall that the Misio lek criterion ( 8) is given by

(19) M I(φ) := M |grad φ| 2 -λ 2 φ 2 dµ < 0,
where φ = {ψ, f } for some function f .

Theorem 5. If ψ : M → R is a fixed smooth function satisfying ∆ψ = -λ 2 ψ and L denotes the operator L(f ) = {ψ, f }, then ∆ -1 L is continuous from Ḣ1 (M ) to itself. The image Im[L]
is a closed subspace of Ḣ1 (M ), and ∆ -1 restricts to a continuous operator on it. The Misio lek criterion is satisfied if and only if the operator norm of

∆ -1 on Im[L] satisfies (20) ∆ -1 Im[L] op ≥ 1 λ 2 .
This happens if and only if there is a number c ≥ 1/λ 2 and a function f

: M → R satisfying -∆ -1 L(f ) = cL(f ).
Proof. Recall that the first-order homogeneous Sobolev metric is given by

f, g Ḣ1 = - M f ∆g dµ = M ∇f, ∇g dµ.
To establish continuity it is sufficient to establish boundedness, which follows from

L(f ) L 2 = grad ψ × grad f L 2 ≤ ψ C 1 grad f L 2 = ψ C 1 f Ḣ1 .
As a result the image of L is a closed subspace in Ḣ1 (M ). Since L is antisymmetric in the L 2 inner product we see that the operator ∆ -1 L is bounded and antisymmetric in Ḣ1 , via

∆ -1 L(f ), g Ḣ1 (M ) = - M L(f )g dµ = - M {ψ, f }g dµ = M {ψ, g}f dµ = -∆ -1 L(g), f Ḣ1 (M ) .
Now the Misio lek index from (19) rescaled by the Ḣ1 norm is given by

M I(φ) φ 2 Ḣ1 (M ) = 1 -λ 2 M φ 2 dµ M |grad φ| 2 dµ = 1 + λ 2 φ, ∆ -1 φ Ḣ1 (M ) φ 2 Ḣ1 (M )
, and this can be made negative for some φ ∈ Im[L] if and only if

sup φ∈Im[L] φ, (-∆) -1 φ Ḣ1 (M ) φ 2 Ḣ1 (M ) ≥ 1 λ 2 ,
and the left-hand side is precisely the operator norm of (-∆) -1 restricted to Im[L]. Since (-∆) -1 is positive, compact, and self-adjoint on Ḣ1 (M ), it is also positive, compact, and self-adjoint on the closed subspace Im [L]. Hence it has a sequence of positive eigenvalues converging to zero, and the operator norm of it is the largest one. So the operator norm is larger than 1/λ 2 if and only if there is a c ≥ 1/λ 2 such that

(21) (-∆) -1 L(f ) = cL(f ).
Solving the eigenfunction problem (21) is rather difficult since L is quite far from invertible-its kernel consists of all functions constant on the level sets of ψ, and in particular any function Φ • ψ for Φ : R → R will be in that kernel. In principle one could integrate the Green function for ∆ -1 along the level sets to try to solve directly for f , or apply the operator L -1 ∆ -1 L repeatedly to a Fourier basis in hopes that it converges to an eigenfunction, but computationally this becomes rather difficult. For example, inverting L involves integrating around the level curves, but since there are always singular hyperbolic points where the gridlines cross, the local behavior can complicate things substantially. As such we take a more indirect approach in the next section, using higher-order Sobolev inner products.

Numerics

Here we give some details on the numerics that helped provide the results of Theorems 2 and 3, and explain the shape of the perturbations generating the variations that satisfy the Misiolek criterion. First we discuss the general setup for any surface M possibly with boundary, then we specialize to M = T 2 . The implementation described in this section and used to generate Figures 4 and5 is available on github: https://github.com/alebrigant/conjugate-points. The goal is to find a function f ∈ C ∞ 0 such that the corresponding Poisson bracket φ = {ψ, f } minimizes some normalized version of the Misio lek index (19). This index can be rewritten in terms of the previously defined operators (22) as

M I(φ) = - M φ(∆ + λ 2 I)φ dµ = M Lf (Λ -λ 2 I)Lf dµ = f, Ωf L 2 ,
where we have used the antisymmetry of the operator L with respect to the L 2 inner product on C ∞ 0 (M ). To minimize M I(φ) we need to constrain f ; otherwise the minimum is either zero or negative infinity. The easiest way to do this is to require that some norm of f be constrained to be 1, which is equivalent to choosing f so that it minimizes the Rayleigh quotient

RQ(f ) = f, Ωf L 2 f, Γf L 2 ,
where Γ is some positive-definite symmetric operator on C ∞ 0 (M ). If Γ has some relatively high Sobolev order, the minimizer we find will be smoother. By the usual calculus of variations method, minimizers of this must satisfy Ωf = cΓf for some c ∈ R.

Equivalently, f is an eigenfunction of the operator Γ -1 Ω, and if we can find an eigenvalue c which is negative, it proves the index form can be made negative, and thus that there is eventually a conjugate point along the geodesic. The idea is therefore to make Γ strong enough as a differential operator that Γ -1 more than

(f ) = f, Ω p f Ḣp f, f Ḣp , where Ω p = Λ -p Ω = -Λ -p L(Λ -λ 2 I)L.
Replacing the L 2 -inner product by the Ḣp inner product yields an operator Ω p that is bounded for p ≥ 2, and compact for p ≥ 3. The problem of minimizing the Rayleigh quotient (23) can be made finitedimensional by decomposing the unknown f in the following spatial Fourier basis Due to the evenness of the cosine function, it is enough to take only positive k indices for j = 0 in the above double sum. Then each function f can be represented by a (N (2N + 2) + 1)-size vector of real coordinates (a jk ) j,k , and the operator Ω p by a square matrix of same size. The matrix representation of Ω p is obtained by multiplying matrices that are all diagonal except for the Poisson bracket operator L matrix (for which we give a formula in the Appendix, for the aid of the reader).

Then the problem of minimizing (23) boils down to finding the minimal eigenvalue of the Ω p matrix, and the function f whose Poisson bracket minimizes the Misio lek criterion is given by the corresponding eigenvector. ,y) associated to all the spatial deformations f shown in this figure induce vector fields Y = sgrad g that make the Misio lek index, and the index form, negative. However the "shape" of this minimizing deformation stabilizes only for p ≥ 2, and for these values we see a characteristic "V-shape" appearing, i.e., a perturbation of cos(x), which we used to obtain formulas ( 9) and ( 13).

In the diagonal case m = n, due to the symmetry in the variables x and y, another minimizer is given by a perturbation of cos(y), the symmetric perturbation obtained by exchanging x and y. To obtain the figures of the bottom row of Figure 4, we constrained the solution to have Fourier coefficient a 01 = 0. Without this constraint, we find the results of Figure 5, where both minimizers appear separately as well as mixed, as in the case p = 1. This is similar to the minimizer found by the second author in [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF], displayed in Figure 7.

The reason for the dominance of the term cos x in the functions appearing in Theorem 2-3 is that this is the lowest-eigenvalue eigenfunction of Λ. Roughly speaking, we first compute L(cos x) = n sin x cos mx sin ny, then obtain (Λ -λ 2 I)L(cos x) = 2mn cos x sin mx sin ny plus a term with a smaller coefficient, and finally

L(Λ -λ 2 I)L(cos x) = - n 2 4 cos x + S,
where S is a linear combination of terms of the form cos(jx + ky), which are eigenfunctions of Λ -p for the eigenvalues (j 2 +k 2 ) -p , for integers j and k. Therefore, after applying the inverse Laplacian Λ -p for a high power p, the correction terms in S all become small as p increases while the cos x term is preserved. This explains the dominance of the V-shape: although L(cos x) does not involve another cos x term, L 2 (cos x) does, and it becomes dominant after applying smoothing operators.

Other open questions

We conclude with some possible other directions to pursue these techniques. An obvious one is to apply this to other steady flows on the torus that come from Laplacian eigenfunctions; the Kolmogorov flows we have discussed have rectangular cells, but one could also consider stream functions of the form ψ = cos (mx + ny) + c cos (mx -ny) which for c = 1 have skew quadrilateral cells. One could also consider stretched tori where the eigenfunctions look like cos + αny) for some positive parameter α, as considered in [START_REF] Drivas | Conjugate and cut points in ideal fluid motion[END_REF].

We have seen that whenever the Kolmogorov eigenvalue is strictly larger than the minimal one λ = 1, the Misio lek index can be made negative. When λ = 1, the question is more delicate. There are some simple eigenfunctions of the Laplacian such as the shear flows generated by sin x, cos x, sin y, and cos y where the fluid flows have nonpositive curvature along them, and therefore no conjugate points by Misio lek [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF]. However there are also stream functions like sin x + sin y which are not shear flows and have nontrivial cells. The Misio lek index must clearly be nonnegative in this case, but there may still be conjugate points detectable by methods such as in [START_REF] Preston | Conjugate point criteria on the area-preserving diffeomorphism group[END_REF].

It is also easy to see how the same techniques could generate conjugate points on other surfaces. It would be interesting to try the same methods to find conjugate points along other eigenfunctions of the 2-sphere: the smallest-eigenvalue function generates rigid rotations and its geodesic is known to have many conjugate points [START_REF] Misio | Stability of flows of ideal fluids and the geometry of the group of diffeomorphisms[END_REF]. It seems likely that every eigenfunction of the Laplacian on S 2 generates a geodesic with conjugate points, but no others are actually known concretely. This is connected with work of Benn [START_REF] Benn | Conjugate points in Dµ(S 2 )[END_REF] on conjugate points along nonsteady Rossby-Haurwitz waves on the sphere.

Figure 1 .

 1 Figure 1. Kolmogorov flows on the torus: level sets of the stream function ψ(x, y) = -cos(mx) cos(ny) for m = n = 1 (left), m = n = 2 (middle), m = 3, n = 2 (right). The velocity of the fluid is tangent to these level sets.

( 13 )

 13 f (x, y) = cos x 1 + a cos (2ny) + b cos (4ny) + c cos (2nx) + d sin x sin (2nx). We again compute (14) φ(x, y) = n sin (nx) cos (ny) cos x -2an sin (2ny) -4bn sin (4ny) -n cos (nx) sin (ny) -sin x (1 + a cos (2ny) + b cos (4ny) + c cos (2nx) -2cn cos x sin (2nx) + d cos x sin (2nx) + 2dn sin x cos (2nx) Computing the Misio lek index (8) now gives

(

  18) f (x, y) = cos x 1

Figure 3 .

 3 Figure 3. Local deformation of the stream function ψ induced by the area-preserving diffeomorphism sgradf , i.e., the function ψ(x -∂ y f, y + ∂ x f ), for m = n = 1 (left), m = n = 2 (middle) and m = 3, n = 2 (right).

5. 1 .

 1 The minimization strategy. Let C ∞ 0 (M ) denote the space of C ∞ -functions f : M → R that vanish on the boundary ∂M . We define the following operators on C ∞ 0 (M ) (22)Lf := {ψ, f }, Λf := -∆f, Ωf := -L(Λ -λ 2 I)Lf.

Figure 4 .

 4 Figure 4. On the top row, minimizers f : M → R of (23) in the case m = 3 and n = 2 and for different values of the order p of the homogeneous Sobolev norm: 0, 1, 2, 3 from left to right. On the bottom row, minimizers in the case m = 2 and n = 2 for the same values of p.

a

  jk cos(jx + ky).

Figure 5 .

 5 Figure 5. Minimizers f : M → R of (23) in the case m = 2 and n = 2 and for different values of the order p of the homogeneous Sobolev norm: 0, 1, 2, 3 from left to right.

  -8n 2 + 1)(6144n 8 + 4864n 6 + 920n 4 + 58n 2 + 1)

	4 + 32n 2 + 1) 6144n 8 + 4864n 6 + 920n 4 + 58n 2 + 1	,
	b 0 = -	512n 6 + 32n 4 -12n 2 -1 2(6144n 8 + 4864n 6 + 920n 4 + 58n 2 + 1)	,
	c 0 = -	49152n 10 + 59392n 8 + 17088n 6 + 1952n 4 + 88n 2 + 1 2(98304n 12 + 28672n 10 -18048n 8 -1568n 6 + 472n 4 + 50n 2 + 1)	,
	d 0 = -	n(32768n 8 + 19456n 6 + 3328n 4 + 196n 2 + 3) (16n 4

Appendix

Proposition 6. Let f : M → M be a function and a jk the coefficients of its Fourier decomposition (24), for (j, k) ∈ {0}×{0, . . . , N }∪{1, . . . , N }×{-N, . . . , N }. Then

where

otherwise.