Junaid Ahmed Khan
email: junaid.khan@inria.fr

Cedric Westphal
email: cedric@soe.ucsc.edu

Yacine Ghamri-Doudane

A Content-based Centrality Metric for Collaborative Caching in Information-Centric Fogs

Keywords: Content Centric Networking, Distributed Content Caching, Fog Networking, Content Offload

Information-Centric Fog Computing enables a multitude of nodes near the end-users to provide storage, communication, and computing, rather than in the cloud. In a fog network, nodes connect with each other directly to get content locally whenever possible. As the topology of the network directly influences the nodes' connectivity, there has been some work to compute the graph centrality of each node within that network topology. The centrality is then used to distinguish nodes in the fog network, or to prioritize some nodes over others to participate in the caching fog. We argue that, for an Information-Centric Fog Computing approach, graph centrality is not an appropriate metric. Indeed, a node with low connectivity that caches a lot of content may provide a very valuable role in the network.

To capture this, we introduce a content-based centrality (CBC) metric which takes into account how well a node is connected to the content the network is delivering, rather than to the other nodes in the network. To illustrate the validity of considering content-based centrality, we use this new metric for a collaborative caching algorithm. We compare the performance of the proposed collaborative caching with typical centrality based, noncentrality based, and non-collaborative caching mechanisms. Our simulation implements CBC on three instances of large scale realistic network topology comprising 2, 896 nodes with three content replication levels. Results shows that CBC outperforms benchmark caching schemes and yields a roughly 3x improvement for the average cache hit rate.

I. INTRODUCTION

Information-Centric Fog Networking has been proposed to allow nodes and devices near the end-user edge to perform caching, communication and computation, thereby reducing the demand on the data center cloud and the communication over the Internet backbone. The goal of Information-Centric Fog Networking is to offload the content and some of its processing closer to its consumers in order to reduce cost, bandwidth consumption and network overhead.

The nodes around the end-user form a network from which the consumer can retrieve the content. However, this network is not as structured and organized as the managed operator network. To identify the structure of this fog network, it has been suggested to consider the graph centrality derived from the graph topology of the network.

Centrality [START_REF] Borgatti | Centrality and network flow[END_REF], a concept from graph theory typically applied to social networks, is used to find important nodes in a graph. A high centrality score reflects a high topological connectivity for a node in the network. Typical centrality measures are: degree (the number of directly connected nodes), closeness (the average length of the shortest paths between the node and all other nodes in the graph); betweenness (the number of shortest paths between all pairs of nodes in the graph going through a specific node), and eigenvector centrality (a measure of node influence in the network).

Computing centrality on the topological graph yields interesting insights. Yet it fails to capture that, in an Informationcentric network, the consumer is interested in connecting to the content, not to a specific node. To address this, we introduce the concept of Content-Based Centrality. We define a contentbased betweenness centrality. Namely, for a node v, its CBC is calculated as the number of shortest paths from all the consumers to all the content which go through v. It is defined formally in Section III.

To show the usefulness of the CBC, we propose a collaborative caching method where a set of connected nodes can mutually self-organize into a fog for distributed content caching. We provide a CBC-based content placement algorithm, where popular content is cached at high CBC nodes, and in similar fashion nodes in the the fog cache content with decreasing CBC score and content popularity.

We evaluate our proposed CBC-based content placement by simulations with three realistic topologies and three content replication ratios. The results show that the proposed CBC-based content placement outperforms typical centrality schemes, or schemes without coordination within the fog. Our contributions can be summarized as follows:

• We introduce the concept of Content-Based Centrality (CBC), which we believe is better suited for Information-Centric Networks than traditional graph centrality; • We introduce a scalable method to compute CBC without a priori knowledge of the content placement, based only upon replication rules; • We present a CBC-based content placement algorithm which takes into account the number of paths to content passing through a node to find nodes as caches in the network;

ISBN 978-3-901882-94-4 c 2017 IFIP

• We evaluate our algorithm through simulations which show a significant improvement. The remainder of the paper is organized as follows. Section II discusses the related work. In Section III we define our network model and introduce our content based centrality metrics. In Section IV, we discuss the performance evaluation and results. Section V concludes the paper along insights into future directions.

II. RELATED WORK

Content Caching has been studied for some time by the research community, spanning a wide spectrum, including Small Cell Networks (SCNs) [START_REF] Bas ¸tug | Cache-enabled small cell networks: Modeling and tradeoffs[END_REF], [START_REF] Bharath | A learning-based approach to caching in heterogenous small cell networks[END_REF], Content Distribution Networks (CDNs) [START_REF] Hamidouche | Many-to-many matching games for proactive social-caching in wireless small cell networks[END_REF] and Information/Content Centric Networking (ICNs/CCNs). For example, in [START_REF] Sourlas | Distributed cache management in information-centric networks[END_REF], distributed cache management decisions are made in order to efficiently place replicas of information in dedicated storage devices attached to nodes in the network using ICN. Similarly [START_REF] Wang | Design and evaluation of the optimal cache allocation for content-centric networking[END_REF] addresses the distribution of the cache capacity across routers under a constrained total storage budget for the network. The authors found that network topology and content popularity are two important factors that affect where exactly should cache capacity be placed. [START_REF] Yu | Congestionaware edge caching for adaptive video streaming in information-centric networks[END_REF] looked at pushing content to the edge to anticipate network congestion, while [START_REF] Azimdoost | Fundamental limits on throughput capacity in information-centric networks[END_REF] computed the capacity of an ad-hoc network of caches.

In a recent work [START_REF] Li | Pricing and resource allocation via game theory for a small-cell video caching system[END_REF], game theory is exploited for caching popular videos at small cell base stations (SBSs). Another work, [START_REF] Mangili | Bandwidth and cache leasing in wireless information centric networks: a game theoretic study[END_REF] proposed a game theoretic approach in ICN to stimulate wireless access point owners to jointly lease their unused bandwidth and storage space to a content provider under partial coverage constraints. Both papers targeted a pricing model instead of providing an efficient content placement solution.

[11] defines a "conditional betweenness centrality" and uses this metric to chooses which nodes will cache the content. The authors in [START_REF] Khan | Saving: Socially aware vehicular information-centric networking[END_REF] proposed new social-aware metrics adaptable to dynamic topology in order to cache content at vehicles. Socially-Aware Caching Strategy (SACS) [START_REF] Bernardini | Socially-aware caching strategy for content centric networking[END_REF] for Content Centric Networks (CCNs) uses social information in order to privilege Influential users in the network by pro-actively caching the content they produce. The authors detect the influence of users within a social network by using the Eigenvector and PageRank centrality measures.

Another centrality based caching approach in CCN is presented in [START_REF] Rossi | On sizing ccn content stores by exploiting topological information[END_REF] where the sizing the content store is based upon centrality. The authors exploit different centralities (betweenness, closeness, stress, graph, eccentricity and degree) to heterogeneously allocate content store at nodes instead of homogeneous allocation. It is proposed that a simple degree centrality based allocation is sufficient to allocate content store. Similarly, [START_REF] Chai | Cache "less for more" in information-centric networks[END_REF] shows that a higher cache hit rate can be achieved if content is cached at high betweenness centrality nodes.

We argue that the topological connectivity only partially relates to the content, and there is a need to consider the content reachability in the network not addressed in the prior art. A well connected node in the topological graph is not necessarily closer to end-users requesting the content. To the best of our knowledge, we are the first to propose to compute a content-based centrality.

III. CONTENT-BASED CENTRALITY

A. System Model 1) Connectivity Model: We assume that nodes are connected to a fog network through local, ad hoc connections. The connectivity between nodes is modeled by a graph G(V, E v ) where V = {v} is the set of nodes and E v (t) = {e jk (t) | v j , v k ∈ V, j = k} is the set of edges e jk (t) modeling the existence of a communication link between nodes j and k. To characterize the concept of Fog network, the local network is represented by the undirected graph G(L, E l ), the set of nodes L = {l} represents different locations l and the set of edges E l = {e pq | l p , l q ∈ L, p = q} are the respective boundaries that connects different neighborhoods. To consider spatio-temporal content placement, we assume the time T = (t 1 , t 2 , ...) as a sequence of regular time-slots, where the k th time-slot is represented as t k = [t k , t k+1 ). We assume that during a particular time-slot the content placement is relatively stable and the nodes will stay connected for a period of minutes or hours, depending on the application. Colocated or closely interconnected nodes can self-organize to form a fog s ∈ V in a location l for distributed content caching locally near end users, where S = {s} , S ⊆ V is the set of fogs in all locations.

2) Caching Model: We define the set of known content as X = {x 1 , . . . , x N } for a catalog of N pieces of content, where x j is an indivisible content chunk in the network. In the remaining of the paper, we will deal with individual content chunks x ∈ X, however a larger size content can be composed of several such content chunks. The nodes can fetch the content either from the service provider using the operator's infrastructure link or locally from peers in the fog, using a lowcost connectivity. The content popularity can be shared with the nodes using three approaches, (i) Offline method by the content operator as a control message. (ii) Local monitoring by the nodes taking into account the number of user interests for the content, and (iii) part of content header shared by the service provider.

B. A Content-Based Centrality Metric

The Content-Based Centrality (CBC) is defined as the sum of the ratio of the number of shortest paths from all users to all content that passes through the node to the total number of shortest paths between all the (user,content) pairs. Figure 1 shows an example of CBC where the node x 1 is on the path from the users to the content in the server, thus, considered as a high CBC node. The other typical node-based centrality measures such as degree, closeness, betweenness and eigenvector consider the node x 2 as the highest centrality node, Figure 1: An example of Content-based Centrality however, compared to x 1 , placing content at x 2 do no achieve optimal caching. Mathematically, CBC can be represented as:

cbc(v) = u =v =x σ ux (v) σ ux , (1)
where σ ux (v) is the number of paths from user node u to content x that passes through the node v and σ ux is the total number of paths from the user node u to the content x. The CBC is then normalized by the number of possible paths to all content in the network such as:

normal(cbc(v)) = cbc(v) -min w∈V (cbc(w)) max w∈V (cbc(w)) -min w∈V (cbc(w))
We use Betweenness centrality here, but we could use other form of centrality as well, such as Closeness or Eccentricity for instance.

C. Scalable Computation of Content-based Centrality

Computing the content-based centrality from Equation 1 is extremely difficult, since the catalog X may be very large, and it requires knowing the content placement a priori, whereas we would like to compute the content-based centrality so as to guide the content placement.

Our solution is to note that the centrality of a piece of content does not depend on the specific chunk x j , but rather, on the relative placement of the multiple copies of a chunk x j . For instance, if N m objects are not cached anywhere in the fog network (and only at the origin server), their contribution to cbc(v) in Equation 1is equal to N m times the number of shortest paths to the origin server going through v divided by the number of shortest paths to the origin server. We do not need to know which N m objects are cache misses, only that there are N m of them.

Further, we define by α the content replication factor that we allow at each node. α corresponds to our replication policy. For a total cache size of b t v at the node, αb t v contains content (i.e. most popular content) common at all nodes, while (1 -α)b t v is the fraction of node caches with unique content compared to other nodes.

With the knowledge of α, we can compute Equation 1 without knowing which content is placed at which node. However, we will later need to enforce this replication rule,

Algorithm 1 Content Placement Algorithm 1: INPUT: V, X, b t v , ∀v ∈ V , p x , ∀x ∈ X, α 2: OUTPUT: s ⊆ V ,X s , X v ∈ X 3: Initialize s = φ, b v = φ, X s = φ 4: for each location l ∈ L, time-slot t ∈ T do 5: v = arg max v∈V (cbc(v)), v / ∈ s 6:
while (b v ≤ αb t v ) do Update b v , X v 14:

X s = X s ∪ X v 15:
end while 16:

s = s ∪ v 17: end for 18: return s by placing content such that a fraction α is replicated at all nodes, and the rest is unique to each node.

This specific replication policy is a simple first step, to demonstrate the feasability of a content-based centrality metric. We can define more sophisticated rules that allow to compute Equation 1 without knowing the specific content placement.

D. Content Placement Algorithm

The problem of content placement can be formulated as follows:

edymaximize

s v∈s X v (l, t), ∀l, ∀ t edysubject to x∈X b x v (l, t) ≤ b t v (l, t), ∀v,
The objective function maximizes the content available at the nodes in the fog at each location and time-slot. The constraint consider the node buffer where an content cached at an individual node buffer should not exceed the maximum available threshold space. We also consider a defined replication factor α where for a node cache size b v , αb v contains content common at all nodes, while (1-α)b v contains content different from other nodes. The CBC is initially computed without any content placement, using the replication rule from the previous section. Then, the distributed content placement is optimized at colocated nodes using Algorithm 1. For a location l and time-slot t, the node v ∈ V with the highest CBC becomes the delegate to locally cache content in the fog as shown in Line 5. It initializes the fog formation by creating the subset s ⊂ S and caches the most popular content, i.e content x with maximum popularity p x , while respecting the replication factor α for its respective storage buffer (Line 6). It continues caching content with decreasing popularity till the storage conditions based on the replication factor α are met. Once the storage buffer is at αb t v , it updated the cached content and start filling the remaining content with decreasing popularity. The condition in Line 11 continues adding content till the node buffer b t v attains its maximum threshold to cache content as indicated in our optimization constraint. The cached content is updated in the fog s (Line 14) and the node adds itself to the fog (Line 16). Similarly, in a decreasing order of node CBC, the remaining nodes cache content in the fog x with same content at αb t v and different at (1 -α)b t v . Thus, Algorithm 1 ensures a maximum amount of popular content are cached in the fog where the content set X v at each node is maximized which in turn maximizes our objective function.

IV. NUMERICAL EVALUATION

We evaluate the proposed CBC along our collaborative caching algorithm using NS-3 where the named-data networking model of the ICN architecture is implemented. Three topologies are extracted from a realistic large scale trace of 2, 986 nodes (vehicles) in a 6X6km 2 area, the city center of Cologne, Germany for one hour in order to validate the scalability of our caching approach. This allows us to evaluate the concept of fogs on realistic topologies reflecting connectivity in an urban environment. The topologies are the time snapshots of the network connectivity at 1 second granularity at the initial, at 30 th minutes and at 60 th minutes respectively.

A. Simulation Scenario

The simulation scenario implements consumer nodes which generate interests for a pre-known content sequence of 100 unique items, each of 1024 kilobytes following a Zipf distribution (coefficient=1), i.e. frequent interests for more popular content. Any provider node already cached the content responds to the consumer interest. We allow intermediate nodes with uniform buffer size to perform in-network caching. We consider 30% of the nodes as consumers, 30% as providers (caching) nodes and the remaining nodes with disabled caching in order to accurately evaluate the performance of caching node. Each simulations is repeated 10 times with varying different nodes as consumers and producers. We also define α as a content replication factor where for a node buffer b v , αb v contains content which is same at all nodes in the fog, while (1 -α)b v is the buffer space for the cached content unique to each node.

As a first step, we randomly populate the caches of nodes in the fog with content in order to compute each node CBC. Random interests are generated. Then, we implement the collaborative placement Algorithm 1 where the highest CBC node caches the most popular content. Similarly content are populated at the remaining nodes' cache with decreasing content popularity and node CBC score. Besides CBC we implement three different caching approaches for comparison:

• Centrality-based Caching popular content at fog of high Degree, Closeness, Betweenness and Eigenvector centrality nodes. Algorithm 1 is implemented using each scheme where the most popular content is placed at the highest centrality node, then remaining content are placed at the nodes with decreasing node centrality and content popularity. • Non-centrality based Social-unaware approach by caching greedily popular content at all nodes in the fog along a Least Recently Used (LRU) based content placement policy. • Non-collaborative based approach where no fog is formed between nodes and each individual node cache indifferently. CBC is used to identify nodes. The following performance metrics are used to evaluate CBC and the proposed collaborative caching algorithm:

• Cache Hits: It is the average number of content responds from the node cache, calculated as the ratio of the number of content responds to the number of received interests by the nodes. • Success Rate: Different from the content responds, it is the average throughput of the nodes, calculated as the ratio of the number of interest responded with content or forwarded by the nodes to the generated interests.

B. Simulation Results

1) Cache Hits: An important evaluation metric for content caching is the amount of cache hits. Therefore, we computed the cache hit rate for each scheme, (i) Fog with CBC, Degree, closeness, Betweenness and Eigenvector centrality, (ii) noncentrality based fog, and (iii) no fog where individual nodes cache indifferently. Figure 2 shows the average cache hit rate of the nodes in fog classified by CBC. We compare the hit rate of our approach for three topologies from the Cologne trace and three replication ratios (α = 25%, 50% and 75%). It is shown that the average cache hit rate achieved differs with respect to topology. Furthermore, we observe from the Topology 3 that between α = 25%, and α = 75% replication ratio, there is an increase in the cache hit rate where around α = 50%, the maximum cache hit rate is achieved.It is also observed that there is a difference in the average hit rate between topologies. Thus, by keeping half of the content with similar content and the other half with different content is the optimal α which results in the maximum hit rate irrespective of the underlying topology.

For a comparative analysis, we implement benchmark centrality schemes in the similar fashion by first placing popular content at the highest centrality node identified by each scheme. Then, following a decreasing node centrality order, we place the content at the set of nodes in the fog with decreasing content popularity order. Figure 3 shows the hit rate comparison where Figures 3a,3b and 3c are the simulation results for each topology. It is clearly shown that irrespective of the network topology, the content-based centrality achieves a high cache hit rate when compared to all other approaches. It resulted in around 50% hit rate for Topology 1, around 50 -60% for Topology 2 and up to 73% for the Topology 3.

CBC is followed by betweenness centrality which achieved around 40% average cache hit rate with a trend similar to CBC in Topology 1 and 3. All the centrality schemes resulted in substantially higher hit rate than the case without non centrality (social-unaware). For all the three topologies, this resulted in around 10% hit rate with a maximum of 18% cache hit rate for Topology 1.

Figure 3 also depicts the result by implementing a case where there is no collaboration between the caching nodes, thus, the case of no fog. It resulted in slightly better perfor-mance than the case of non-centrality based caching, though it achieves a lower average hit rate where it merely achieves around 20% cache hit rate. Thus, the comparative analysis of cache hit rate on three different topologies and three content replication ratios reveals that content-based centrality achieves a better hit rate. It is because CBC the number of paths to content instead of nodes in a content-centric approach to classify nodes.

2) Success Rate: As an important network metric, we evaluate the success rate as the content responds or forwarding by the nodes in the Fog. In Figure 4, we show the average success rate achieved by the nodes in Fog classified using CBC. Similar to the cache hit rate, our analysis is based on three topologies and three replication ratios. Each topology resulted in a different success where Topology 3 shows the best performance with more than 60% success rate for all replication ratios, followed by Topology 2 and Topology 1, thus suggesting that topology impacts performance. Another observation is that the success rate results from all topologies follow similar trend, thus validates the robustness of CBC when applied to different topologies. We see that by varying the replication ratio, a slight increase in the success rate is achieved for the 50% replication ration where the success rate decreases towards 75% replication ration. This suggest that for a node, it is optimal to keep half of the content similar to other nodes and the other half different from the other nodes.

Figure 5 shows a comparative analysis of the success rate achieved by different caching approaches. We observe that overall, CBC outperforms all other approaches, yielding an average success rate of around 50, 55 and 60% for Topology 1, 2 and 3 respectively. An interesting behavior is seen in the case of Topology 1 (Figure 5a) where fogs based on CBC, degree and betweenness centrality have similar success rate, though their score differs between replication ratio. For 25% replication ratio, both degree and betweenness centrality show a higher success rate (60%) than CBC, for 50% replication ratio, CBC dominates betweenness centrality while degree stays higher, and for 75% replication ratio, both degree and CBC show similar results (i.e. 47% success rate), while betweenness dominates with around 58% success rate. We investigated this behavior and found that since both betweenness and degree are node-centric metrics, and therefore better connected nodes at the relative center of the topology, they are acting as bridge between large number of nodes, thus resulting in a higher success rate. Nevertheless, nodes classified by CBC are equally better performing despite not being well placed in the network.

On the other hand, the non-centrality based sociallyunaware resulted in the poorest performance in all topologies. Similarly, in the case without fog formation, i.e. noncollaboration resulted in poor performance with an overall success rate of less than 20%. A collaborative approach should definitely be used to for content caching.

V. CONCLUSIONS AND FUTURE DIRECTIONS

Our results show that content placement in ICN can be efficiently managed by forming fogs allowing nodes to collaboratively cache at high content-based centrality nodes. "content" based centrality outperforms the node-centric approach to classify nodes for distributed caching.

It is challenging to find the appropriate nodes to be strategically selected for efficient content caching in the ICN network. In this paper, we targeted the problem by proposing a collaborative distributed caching approach at connected nodes near the network edge as in a fog network. To do so, we suggest to exploit graph properties to identify the suitable nodes for caching, however, unlike typical node-centric centrality scheme, we first presented CBC, a new Content-based centrality scheme where the number of paths to all content, instead of nodes, counts towards the node centrality. Then, we proposed an algorithm for collaborative content placement in the fog, where the nodes place popular content at high CBC nodes, followed by placing the remaining content with decreasing popularity at nodes with decreasing CBC score, according to a replication rule which allows to compute CBC without knowing the actual content placement.

We evaluated the benefits of CBC via simulations on three different topologies and varying three content replication levels. Results show that nodes in the fog based on CBC outperforms, in terms of cache hit rate and success rate, existing centralities, non-centrality based and non-collaborative caching approaches. Future work includes moving to dynamic topologies and study the impact of mobility on CBC. We plan to use CBC for content retrieval where the user interests are forwarded to high CBC nodes.

Figure 2 :

 2 Figure 2: Average cache hit rate achieved by CBC on three different topologies with three content replication ratios.

Figure 3 :Figure 4 :

 34 Figure 3: Hit rate comparison for centrality-based (CBC, Degree, Closeness, Betweenness, Eigenvector), non-centrality based or social-unaware (LRU) and non-collaborative (no fog) based caching.

Figure 5 :

 5 Figure 5: Success rate comparison for centrality-based (CBC, Degree, Closeness, Betweenness, Eigenvector), non-centrality based or social-unaware (LRU) and non-collaborative (no fog) based caching.


F Leblanc

A V Oza

L Leclercq

C Schmidt

T Cassidy

R Modolo

J Y Chaufray

R E Johnson

On the orbital variability of Ganymede's atmosphere

Ganymede's atmosphere is produced by radiative interactions with its surface, sourced by the Sun and the Jovian plasma. The sputtered and thermally desorbed molecules are tracked in our Exospheric Global Model (EGM), a 3-D parallelized collisional model. This program was developed to reconstruct the formation of the upper atmosphere/exosphere of planetary bodies interacting with solar photon flux and magnetospheric and/or the solar wind plasmas. Here, we describe the spatial distribution of the H 2 O and O 2 components of Ganymede's atmosphere, and their temporal variability along Ganymede's rotation around Jupiter. In particular, we show that Ganymede's O 2 atmosphere is characterized by time scales of the order of Ganymede's rotational period with Jupiter's gravity being a significant driver of the spatial distribution of the heaviest exospheric components. Both the sourcing and the Jovian gravity are needed to explain some of the characteristics of the observed aurora emissions. As an example, the O 2 exosphere should peak at the equator with systematic maximum at the dusk equator terminator. The sputtering rate of the

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T I Introduction

Ganymede is the only known satellite with an intrinsic magnetosphere [START_REF] Kivelson | Discovery of Ganymede's magnetic field by the Galileo spacecraft[END_REF]. Its icy surface is eroded by the solar radiation flux and by the Jovian energetic plasma leading to a locallycollisional atmosphere and a locally surface-bounded exosphere [START_REF] Stern | The lunar atmosphere: history, status, current problems and context[END_REF][START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF]Turc et al. 2014;[START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF][START_REF] Shematovich | Neutral Atmosphere Near the Icy Surface of Jupiter's Moon Ganymede[END_REF]. Unfortunately, Ganymede's atmosphere remains essentially unknown despite the flybys during the mission Galileo in the 90s. This mission detected only one atmospheric emission, Lyman α, interpreted as being produced by resonant scattering by exospheric atomic hydrogen [START_REF] Barth | Galileo ultraviolet spectrometer observations of atomic hydrogen in the atmosphere of Ganymede[END_REF]) with a column density around 9.2110 12 H/cm 2 . Such a column density and its variation with distance to Ganymede suggested a surface density around 1.510 4 H/cm 3 and an exospheric temperature around 450 K. Hubble Space Telescope was also able to identify another component of Ganymede atmosphere, the O 2 molecule through the observations of two UV atomic oxygen emissions with intensities corresponding to a column density around 1 to 1010 14 O 2 /cm 2 [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF]. Since this discovery, several other HST observations of these emission have revealed their auroral origin [START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF][START_REF] Mcgrath | Aurora on Ganymede[END_REF][START_REF] Saur | The search for a subsurface ocean in Ganymede with Hubble Space Telescope observations of its auroral ovals[END_REF] highlighting the relations between the spatial distribution of these emissions, Ganymede's intrinsic magnetospheric structure [START_REF] Kivelson | Discovery of Ganymede's magnetic field by the Galileo spacecraft[END_REF]) and Ganymede's surface's albedo [START_REF] Khurana | The origin of Ganymede's polar caps[END_REF]). [START_REF] Gurnett | Evidence for a magnetosphere at Ganymede from plasma-wave observations by the Galileo spacecraft[END_REF] published the first observations of Ganymede's ionosphere by the plasma wave spectrometer (PWS) on board Galileo suggesting ion densities on the order of 100 cm -3 and a scale height around 1000 km. This first set of measurements was later complemented by a second flyby (Eviatar et al. 2001b) confirming the density range but suggesting that a number density on the order of 1000 cm -3 at the surface may be more realistic and in better agreement with the upper limit set by

A C C E P T E D M

A N U S C R I P T 4 radio occultation [START_REF] Kliore | Satellite atmospheres and magnetospheres[END_REF]. Measurements by the Plasma Science (PLS) instruments were also performed during Galileo flybys in the polar regions of Ganymede [START_REF] Frank | Outflow of hydrogen ions from Ganymede[END_REF]). These authors reported the observations of hydrogen ions outflowing from Ganymede's poles with temperatures around 4×10 4 K and surface densities around 100 cm -3 . These observations were later reanalyzed by Vasyliunas and Eviatar (2000) suggesting that O ions rather than H ions constitute this outflow, with similar densities as derived by [START_REF] Frank | Outflow of hydrogen ions from Ganymede[END_REF]. These sets of observations of the thermal ionosphere were complemented by energetic ion measurements obtained by EPD/Galileo, showing that Ganymede is permanently embedded in an energetic plasma mainly composed of S + and O n+ [START_REF] Paranicas | Energetic particle observations near Ganymede[END_REF]).

After Galileo flybys, a simple 1-D model of the atmosphere and ionosphere was proposed by Eviatar et al. (2001a) to fit PWS profiles with O 2 densities at the poles around 10 6 cm -3 . The first attempt to reconstruct the origins and fate of Ganymede's neutral environment was described in [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] who developed a 2-D axi-symmetric model of this atmosphere using a combination of Direct Simulated Monte Carlo and Test-particle Monte Carlo simulations focusing on the sunlit trailing hemisphere. He suggested neutral densities peaking around 10 9 cm -3 around the subsolar region and 10 8 cm -3 at the poles with the two main atmospheric species being H 2 O and O 2 . Turc et al. (2014) developed the first 3-D model of Ganymede's atmosphere starting from [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] simulation for the sunlit trailing hemisphere. These authors found essentially the same surface density but some discrepancies were found at high altitudes. Turc et al. (2014) found a Lyman α emission brightness up to 70 times smaller than observed by [START_REF] Barth | Galileo ultraviolet spectrometer observations of atomic hydrogen in the atmosphere of Ganymede[END_REF], whereas [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] had to multiply his simulated H Lyman α by a factor 4 to reproduce the observed intensities. Both models were able to reproduce the O 2 emission intensity observed by HST [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF][START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF].

Furthermore, [START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF] developed a 3-D Monte Carlo model of Ganymede's exosphere, including a model of ion precipitation in the auroral regions using a global MHD model of Ganymede's magnetosphere [START_REF] Jia | Properties of Ganymede's magnetosphere inferred from improved three-dimensional MHD simulations[END_REF]. They modeled Ganymede's sunlit leading hemisphere highlighting the different exospheric regions of the two main components H 2 O and O 2 .

A C C E P T E D M A N U S C R I P T 5

In this paper, we improve Turc et al. (2014) model, implementing a collisional scheme in order to identify the changes induced by atmospheric particle collisions. [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] nicely illustrated the remarkable nature of this exosphere/atmosphere by plotting his estimate of the Knudsen number around Ganymede. Ganymede's atmosphere is predicted to be locally collisional (around the subsolar region) and collisionless elsewhere. Another aspect of Ganymede that was not described in [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] or Turc et al. (2014), is how this atmosphere might change along Ganymede's orbit due to both its rotation and Jupiter's gravitational field. These new simulations take into account the revolution of Ganymede around Jupiter and Jupiter's gravitational field. In Section II, we present the main improvements of our model with respect to Turc et al. (2014). Section III presents how this work could contribute to a better understanding of Ganymede's atmospheric formation and evolution and is followed by a brief comparison with the few observations of this atmosphere and with other similar observations in our solar system in section VI.

II Exospheric Global Model

II.1 Application to Ganymede

A first version of the Exospheric Global Model (EGM) was applied to the description of Ganymede's atmosphere in Turc et al. (2014). This is a 3-D Monte Carlo model describing the fate of test-particles, ejected from a surface or an atmosphere of a planet or satellite, and subsequently followed under the effect of gravitational fields until their escape (when reaching an upper limit set at few planetary radii from the surface), their ionization/dissociation (by photo and/or electron impacts), or until their readsorption at the surface (which is species dependent). In order to reconstruct density, velocity, and temperature of the exospheric species (minor and major components), this model was parallelized inducing a very significant reduction of the restitution time and the possibility to use much larger total memory. These two improvements were crucial to describe the spatial and temporal distributions of all exospheric species in synchrony, that is, in the case of Ganymede's icy

A C C E P T E D M A N U S C
A C C E P T E D M A N U S C R I P T 7
When impacting the surface, test-particles can either stick or be re-ejected at once. As in Turc et al. (2014) and [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF], electron ionization and dissociation impacts are described considering a uniform mono-energetic electron in the Alfven wings of Ganymede's magnetosphere (open field line regions, see section II.4) with an electron number density of 70 cm -3 and average energy of 20 eV. In the closed field lines region, we neglected ionization and dissociation by electron impacts. We assume that O 2 and H 2 species never stick to the surface, whereas all other species have a non-zero reabsorption probability.

Jupiter's gravitational force at the surface of Ganymede is equivalent to a few percent of Ganymede's gravitational force. However, because certain species have very long lifetimes in Ganymede's exosphere (e.g. H 2 and O 2 ), Jupiter's gravity may play a significant role in shaping its exosphere, as will be shown in Section III. Because we track molecules in Ganymede's reference frame, the Coriolis and centrifugal forces act on Ganymede's molecules in its atmosphere throughout the orbit. This is only relevant for certain species as the rotation period is 7.2 days. In Ganymede's inertial frame, the subsolar point rotates continuously, so that the day-to-night surface temperature gradient drives a global migration of the exospheric species towards the nightside. On the nightside, these particles can be readsorbed on the surface or thermally accommodate and be immediately re-emitted generating a peak in atmospheric density near the surface, as observed on the Moon [START_REF] Hodges | The lunar atmosphere[END_REF]. We also introduce the possibility that returning volatiles can permeate the regolith or be transiently adsorbed. These can eventually return to the exosphere, a process potentially important for minor exospheric species, like Na, at Europa [START_REF] Leblanc | Origins of Europa Na cloud and Torus[END_REF].

The standard outputs from this simulation are the 3-D values of the density, velocity, and temperature of all species considered. We therefore define a 3-D grid with 100 cells in altitude distributed exponentially from a cell of width 3 km in altitude at the surface to 300 km cell at 3 Ganymede radii (R G ), with 40 cells linearly spaced in longitude, and 20 cells in latitude, defined such

A C C E P T E D M A N U S C R I P T
8 that all cells have the same volume at a given altitude. We also reconstructed 2-D column density maps as derived from the 3-D neutral density, 2-D emission rates for the three main emission features observed around Ganymede (Lyman α and 130.4 and 135.6 nm oxygen) and a 2-D map of the surface reservoir. The model also calculates 3-D ionization rates by electron and photon impacts.

In the following, we describe the main improvements in the model of Turc et al. (2014), namely, a new thermal surface temperature model of the surface (Section II.2), a new description of sublimation (Section II.3), a new description of sputtering (Section II.4), the introduction of collisions (Section II.5) and a fusion scheme for the test-particles (Section II.6).

II.2 Surface temperature

Ganymede's surface temperature is estimated using a one-dimensional heat conduction model in which the thermal properties of the ice are treated as two layers [START_REF] Spencer | Systematic biases in radiometric diameter determinations[END_REF]; available via http://www.boulder.swri.edu/ spencer/thermprojrs). A low thermal inertia, 2.2×10 4 erg cm -2 s -1 K -0.5 , is assumed in the first few cm where the ice is heavily weathered, while the ice at depth has a much larger thermal inertia of 7×10 4 erg cm -2 s -1 K -0.5 , resulting in a longer thermal time scale. Density is set to 0.15 and 0.92 g/cm 3 in the top and bottom layers, respectively. The albedo is assumed to vary linearly over the body's longitude, from 0.43 on the leading hemisphere to 0.37 on the trailing [START_REF] Spencer | Icy Galilean Satellite Reflectance Spectra: Less Ice on Ganymede and Callisto?[END_REF], and the emissivity is set at 0.96. A small correction for the latent heat of sublimation is included [START_REF] Abramov | Numerical modeling of endogenic thermal anomalies on Europa[END_REF]. Solar insolation over a synodic day is determined in degree increments of longitude with an eclipse length near equinox. This approach provides a basic description of longitudinal asymmetries about the subsolar axis that result from the thermal inertia of the ice, but does not account for any local variation of thermophysical parameters. Figure 1a shows the subsolar temperature as a function of orbital longitude, which results from shadowed sunlight during eclipse and albedo. Figure 1b maps the surface temperature at 270° orbital longitude (trailing) where a temperature of 146° K is the diurnal maximum. As shown in Figure 1, the subsolar surface temperature peaks when Ganymede's trailing hemisphere is illuminated and is minimum at leading because of the change in surface albedo. The surface temperature displays a small but clear shift towards the afternoon due to lag from thermal inertia.

The afternoon side (between 270° to 0° west longitude) of the surface is therefore significantly hotter than the morning side (between 180° and 270° west longitude).

II.3 Sublimation

In [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF], the sublimation rate at Ganymede's surface at a given temperature was given by

F = a × T s -0.5 × exp(-b/T s ) [cm -2 s -1 ] (2)
with a=1.1×10 31 cm -2 s -1 K 1/2 and b = 5737 K. Fitting the calculation by [START_REF] Johnson | Erosion of Galilean satellites by magnetospheric particles[END_REF], we found different values for a = 1.92×10 32 cm -2 s -1 K 1/2 and b = 6146 K. [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] published

A C C E P T E D M A N U S C R I P T 10
another formulation for this sublimation rate, also used by Vorburger et al. (2015). We fit the method of Vorburger et al. (2015) with equation ( 2) and derived values for a=2.17×10 32 cm -2 s -1 K 1/2 and b = 5950 K (section II.2). These three sublimation rates are shown in Figure 2 below. We also considered a case where sublimation rate was significantly decreased (low sublimation case in Figure 2) to simulate an H 2 O exosphere dominated by sputtering. The [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] sublimation model represents an ideal case, based on experiments for pure crystalline ice. In practice, sublimation is highly dependent on the local grain structure and on an exact knowledge of the surface temperature, an error of 10 K giving a few orders of magnitude difference in sublimation rate. We therefore consider that even this low sublimation scenario might be within the realm of plausibility at Ganymede.

Figure 2: Sublimation rate with respect to surface temperature as suggested by [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF], red solid line, by [START_REF] Johnson | Erosion of Galilean satellites by magnetospheric particles[END_REF], cross black symbols and by [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] in square green symbols. The dashed blue line is the rate used for the low sublimation case throughout this paper.

A C C E P T E D M

A N U S C R I P T 11

In the results presented below, we will use the Low sublimation case (Sections III, III.1 and III.2) to model the sublimation as well as [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] sublimation scenario (in Section III.3) to illustrate how a strong sublimation rate could change the structure and content of the exosphere.

II.4 Sputtering

Sputtering is due to energetic S n+ and O n+ Jovian particles that penetrate Ganymede's intrinsic magnetosphere through the polar regions and impact the polar surface regions [START_REF] Cooper | Energetic ion and electron irradiation of the icy Galilean satellites[END_REF].

This bombardment leads to the ejection of surface species, in particular H 2 O, O, H, and OH, but also the production of O 2 and H 2 by radiolysis in the regolith in which the molecules leave with near thermal speeds. In order to describe the spatial distribution of this bombardment, in Turc et al.

(2014), we arbitrarily limited the sputtering to latitudes above 45° as did [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF]. However, By excluding the closed field line regions, we neglect any surface production of exospheric particles by sputtering here, because according to [START_REF] Fatemi | On the formation of Ganymede's surface brightness asymmetries: Kinetic simulations of Ganymede's magnetosphere: plasma interaction with Ganymede[END_REF] and Leclercq et al. (2016), the total precipitation rate in the closed field line region is much smaller than the total precipitation rate in the opened field line region (see also [START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF]. However, ions reaching the closed field line region are more energetic and should, on average, correspond to particles that sputter with higher yield, an issue that has not been addressed in [START_REF] Fatemi | On the formation of Ganymede's surface brightness asymmetries: Kinetic simulations of Ganymede's magnetosphere: plasma interaction with Ganymede[END_REF] and Leclercq et al. (2016). In the following, we neglect this contribution for simplicity. In the open field line regions, we assume a uniform incident ion flux on the surface as suggested by [START_REF] Fatemi | On the formation of Ganymede's surface brightness asymmetries: Kinetic simulations of Ganymede's magnetosphere: plasma interaction with Ganymede[END_REF] and Leclercq et al. (2016). No variation of the precipitating flux with respect to Jupiter plasma centrifugal equator or sheet has been included in this modelling for simplicity.

To calculate the efficiency of an impacting ion to sputter surface particles, we include the surface temperature dependence of the yield by [START_REF] Cassidy | Magnetospheric ion sputtering and water ice grain size at Europa[END_REF]. Following their approach, we used [START_REF] Johnson | Composition and Detection of Europa's Sputter-Induced Atmosphere[END_REF] and [START_REF] Fama | Sputtering of ice by low-energy ions[END_REF] approximations of the yield (number of molecules and atoms ejected per impacting particle) with respect to the incident energy:

Y H2O (T) = Y H2O (T 0 )× (1 + 220×exp(-0.06 eV/k B T) (3) and Y O2 (T) = Y O2 (T 0 ) × (1 + 1000×exp(-0.06 eV/k B T) (4)
for sputtered H 2 O molecules (equation ( 3)) and O 2 molecules (equation ( 4)). From these two equations, and knowing the expected flux of ion impacting the surface, we calculate the sputtered

A C C E P T E D M A N U S C R I P T
13 flux from Ganymede's icy surface. However, the value of Y O2 (T 0 ) for O 2 and Y H2O (T 0 ) for H 2 O needs to be estimated as it depends on many parameters, such as the energy and charge of the impacting particle, the surface structure, and the non-ice components [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF]. Another approach is to choose Y(T 0 ) in order to reproduce the observed quantities in Ganymede's exosphere. For instance, the observed FUV O 2 auroral emission can be directly used to constrain Y(T 0 ) (see also the In order to properly describe the entire energy distribution up to the escape energy, we bin these two energy distributions into 15-to-25 energy intervals (from 0 to an energy slightly above the escape energy) and simultaneously eject 15 to 25 individual test-particles with a weight calculated as the integral of the energy distribution over the energy range of each bin. [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF] suggested that O 2 molecules could be ejected following the sum of two distributions: 70% of the O 2 being ejected following a Maxwellian Boltzmann distribution and 30% being ejected following a distribution in E/(U+E) 3 [START_REF] Johnson | Plasma Ion-induced molecular ejection on the Galilean satellites: energies of ejected molecules[END_REF]. Such a choice would populate the higher altitudes slightly more efficiently than in our approach but would not change the main conclusions of this paper.

II.5 Collisional scheme

Because in some regions of Ganymede's atmosphere, collisions between atmospheric particles might shape the density [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF], we also include collisions between molecules and atoms. At each time step, the total density of each species in a cell is calculated. From these densities, the maximum number of collision pairs between each test-particle and the other test-particles for a given species in the cell is calculated using the approach described in [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] with cross sections taken from [START_REF] Lewkow | Precipitation of energetic neutral atoms and induced non-thermal escape fluxes from the Martian atmosphere[END_REF]. These authors developed a new scheme including cross sections depending on both energy and angle that allows an accurate description of collision between any species colliding with energies below a few tens of eV. [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. One of the difficulties when including minor species is that test-particles, in a given cell, might have very different weights (i.e. they might represent vastly different number of real particles). We therefore developed a new scheme to treat such collisions.

For a test-particle A of weight W A and species S A , we randomly determine with which species, S B , a test-particle A will most likely collide. We then randomly select n A test-particles B i of weight W Bi of the species S B in order to have ∑ ≥ W A . For each pair of test-particles B i and A, we then evaluate the likehood of a collision taking into account the relative velocity between A and B i . If a collision occurs, we determine new velocity vector after the collision, V' Ai and V' Bi following Lewkov and Kharchenko (2014). We then create a new test-particle A i with velocity V' Ai and weight and repeat that scheme n A times. The weight of the particle A is then reduced by W Bi . The velocity of each B i test-particle after collision is equal to V' Bi with the exception of the n A th B i test-particle which is divided into two particles: one with velocity V' BnA and weight W' BnA = W A -∑ and the other with velocity V BnA and weight W BnA -W' BnA .

All collisions occur with energies below 10 eV so that molecular dissociation is negligible.

II.6 Fusion scheme

Since the collisional scheme implies the creation of new particles (see section II.5), we also had to implement a fusion scheme in order to limit the number of test-particles followed by the model at a given time. We tested two approaches to merge test-particles, both suggested by [START_REF] Lapenta | Particle Rezoning for Multidimensional Kinetic Particle-In-Cell Simulations[END_REF], the binary scheme, and the tertiary scheme. The binary scheme is faster and since it leads to similar results, we employ it. The binary scheme consists in selecting a pair of test-particles A and B which are close each other in phase space (with positions q A and q B , velocities V A , V B and weight W A , W B respectively) of the same species and in the same cell that minimizes (W A + W B )(V A ² -V B ²). We then

A C C E P T E D M A N U S C R I P T 16
create a new particle with weight equal to W=W A +W B , position equal to q = (W A q A + W B q B )/(W A +W B )

and velocity equal to V = (W A V A + W B V B ) / W.

For every arbitrarily fixed number of time steps (typically one hundred time steps), we count the number of particles of each species within one cell. If this number is larger than a threshold, defined as to avoid too many test-particles of a given species in a cell, the fusion scheme is activated, adapting to the computational load this presents. Test-particles are therefore selected following the approach described above, and are merged as long as this number continues to exceed the threshold.

III Ganymede's atmosphere

In the following, we will use the phase angle of Ganymede (equivalent to its sub-observer longitude)

to indicate its position around Jupiter. A phase angle of 0° corresponds to Ganymede being in the shadow of Jupiter at midnight local time (eclipse), whereas phase angles of 90° and 270° correspond to the sunlit leading hemisphere apex and sunlit trailing hemisphere apex respectively. Positions at the surface are defined using west longitude (increasing with local time), 0° planetary longitude being the subsolar point at a 0° phase angle. The simulation is performed in a GphiO rotating frame with the x-axis always pointing along the corotational direction, y-axis pointing towards Jupiter, and the zaxis pointing towards Ganymede's north pole (the x and y axis are inside the orbital plane of Ganymede). Section III.1 is dedicated to the description of the fate of Ganymede's exosphere along its rotation around Jupiter, without taking into account the effect of collision between atmospheric particles, an issue specifically discussed in Section III.2. Section III.3 is based on simulations where collisions were also neglected.

Using the sputtering description in section II.4, the typical surface temperatures of Figure 1, and the low sublimation rate as described in Figure 2, the global sublimated flux is equal to 8.0×10

III.1 Seasonal variation of Ganymede's exosphere

A main improvement with respect to previously published models of Ganymede's atmosphere is that we describe the evolution of Ganymede's exosphere throughout Ganymede's orbit around Jupiter while taking into account the influence of Jupiter's gravitational field. By simulating a period of time long enough to cover several Ganymede orbits, our model simultaneously calculates the motions of all individual test-particles in Ganymede's inertial frame as well as the motion of Ganymede around Jupiter. In general, only three orbits are needed to achieve a steady state exosphere on Ganymede.

A C C E P T E D M A N U S C R I P T 18
Steady state is evaluated at a fixed phase between two successive orbits for a given time t, when the reconstructed density, velocity, and temperature of the O 2 and H 2 O exospheric components and their loss rates vary by less than a few percent. In the results presented in this section, collisions are not taken into account for the following results and will be specifically discussed in section III.2.

Most previous models have focused on one specific position of Ganymede corresponding to fixed solar illumination and ram directions: the trailing illuminated hemisphere case in Turc et al. (2014) and [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF], and the leading illuminated hemisphere case in [START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF]. However, Ganymede's rotation implies that the geographical regions where sputtering and sublimation rates reach a maximum, rotate with respect to one another. One direct consequence of this rotation is a variable sputtering efficiency since yields depend strongly on the surface temperature and, therefore, on the solar zenith angle [START_REF] Cassidy | Magnetospheric ion sputtering and water ice grain size at Europa[END_REF].

The typical time for O 2 , ejected at the equator to migrate from one hemisphere to the other (that is to move from a position at the surface to its opposite position on the other side of the satellite) can be roughly estimated following [START_REF] Hunten | The Mercury atmosphere[END_REF]. According to these authors, the collisionless migration time of a particle moving in a surface-bounded exosphere is on the order of t m = t b ² where t b = 2(H/g) 1/2 is the ballistic hop time and =R G /H. R G is Ganymede's radius and H is the scale height H=k B T S / (mg) with k B the Boltzmann constant, T S the surface temperature (as defined in section II.2), m the mass of the particle and g the gravitational acceleration at the surface. For an O 2 ejected from noon at the sunlit trailing hemisphere to migrate across the hemisphere is, t m 765 hours or more than 4 orbital periods using Ts=143.5 K (see Figure 1a), H = 26 km, =101, t b =270 s. Migration across the hemisphere takes t m 1690 hours, or more than 9 orbits using 85K, H=15 km, =171, t b =208 s.

Precise migration times are not very meaningful, since t m well exceeds the diurnal timescale of the thermal wave, yet, the calculation is useful to conceptualize the atmosphere's diffusion. Given that the typical lifetime of an O 2 is on the order of one orbital period, newly ejected O 2 roughly fixed in Ganymede's reference frame moves therefore in local time due to Ganymede's rotation. The typical

A C C E P T E D M A N U S C R I P T

19 time scale for the formation of the O 2 exosphere can be estimated by dividing the average vertical column density (2.4×10 14 O 2 /cm 2 ) by the average ejection rate (3.5×10 8 O 2 /cm 2 /s) or 190 hours, that is, more than one orbital period. Therefore, the O 2 exosphere at a given orbital position will depend on its history and evolution over about one Ganymede revolution. Because the H 2 O molecules are quickly reabsorbed in the surface, the timescale of the H 2 O exosphere is much shorter than one rotational period. In the case of H 2 , due to its low mass, the H 2 exospheric will rapidly form across Ganymede becoming spatially uniform on timescales much shorter than an orbital period and is, therefore, less dependent on localized production. than on the ram hemisphere (as suggested by the opened/closed field lines boundary in Figure 3), more H 2 O is therefore ejected on the wake hemisphere than on the ram hemisphere. This magnetospheric effect, when convolved with the strong thermal dependence of radiolysis (Equation 3) and sublimation (Figure 2), effectively increases the H 2 O ejection rate at the sunlit leading hemisphere/plasma wake (see also [START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF]. We therefore simulate a slightly denser H 2 O exosphere on the leading illuminated hemisphere position (radial column density of 1.5×10 13 H 2 O/cm 2 ) than on the trailing illuminated hemisphere position (radial column density of 1.3×10 13 H 2 O/cm 2 ). Observed from the Earth, the average line-of-sight (LOS) column density at a phase angle of 90° over the hemisphere disk would be equal to 3.4×10 13 H 2 O/cm 2 , whereas it would decrease to 2.9×10 13 H 2 O/cm 2 at a phase angle of 270° (see also Figure 6). H 2 O in the equatorial plane has densities smaller than 10 5 H 2 O/cm 3 coming essentially from the polar region, whereas densities at the poles, are greater by two orders of magnitude. Since H 2 O molecules efficiently stick to the surface, the polar H 2 O molecules do not migrate efficiently towards the equatorial regions and remain confined around their ejection region.

A C C E P T E

The O 2 exosphere is significantly different at the equatorial regions because the probability for O 2 to stick to the surface, even at coldest nightside surface temperatures, is low. Consequently, the O 2 exosphere is spread across all longitudes and latitudes. There is however a clear dusk/dawn asymmetry of the O 2 spatial distribution with a peak in density at afternoon between =90° to =200°, Ganymede local time and in the early night between =90° to =210°. The average radial column density varies from 2.2×10 14 O 2 /cm 2 at a phase angle of 20° to a peak of 2.5×10 14 O 2 /cm 2 at a phase angle of 150°. Both the asymmetry between dusk/dawn terminators and the variation of the exospheric content are a clear result of the production and loss timescales of the O 2 exosphere along Ganymede's orbit and can be understood by the variation of the ejection pattern of the O 2 particles as described below and in Figure 5. O 2 /cm 2 /s) because of the eclipse induced by Jupiter's shadow and the rapid decrease of Ganymede's surface temperature (Figure 1a). Since the modeled atmosphere is nearly stationary, a variation in the ejection or loss rate produces a variation in the total content of the O 2 exosphere. Therefore, the total O 2 content reaches a minimum around a phase angle of 20° due to the low ejection rate during the previous half of the rotation, between phase angles 200° to 20° with a peak around a phase angle of 150° following the maximum in ejection rate around a phase angle of 90°. That is, Ganymede's O 2 exospheric total content decreases during a portion of its orbit centered around a phase angle of 330° where the loss rate due to electron impact ionization/dissociation is larger than the ejection rate. During this portion of the orbit, Ganymede's decaying exosphere is a residual that formed earlier and conversely, during the following portion of the orbit around a phase angle of 150°, the total content increases. The O 2 exosphere also appears slightly more extended on the dayside (Figure 4) due to larger surface temperatures. In our model, O 2 is mainly ejected from the higher latitudes by sputtering as explained in section II.4 (Figure 5). However, the vertical column density of the O 2 clearly peaks at the equator (Figure 4b) and not at high latitudes. The O 2 diffusion time suggests that the molecules remain close to their source regions. However, this does not take into account the effect of the centrifugal force which tends to drive molecules equatorward. This force has an average intensity equal to few percent of Ganymede's gravitational force at the surface and can therefore, shape the O 2 exosphere as seen in Figure 6a. The observed dusk/dawn asymmetry seen in Figure 4a leads also to larger column density on the right (dusk) side of the exosphere as seen from the Sun in Figure 6a. This dusk/dawn ratio in the O 2 column density peaks at a value of 3.4 at a phase angle of 150°. Ganymede's exosphere near the leading hemisphere portion of the orbit is formed by the production and spreading of the high latitude ejecta as seen from Figure 5. The dusk/dawn asymmetries simulated here are also produced when simulating the orbital evolution of Europa's O 2 exosphere using the EGM [START_REF] Oza | Origin and Evolution of Europa's Oxygen Exosphere[END_REF], in accordance with recent HST orbital observations of auroral features there [START_REF] Roth | Europa's far ultraviolet oxygen aurora from a comprehensive set of HST observations[END_REF])

A C C E P T E D M A N U S C R I P T 24
The H exosphere is produced from the polar region by sputtering in our model. The column density, as seen from the Sun shown in Figure 6b, reflects the source region and that the ejection energy primarily leads to escape. Sputtered H 2 , on the other hand, is ejected at the local surface temperature at lower energies, and with its higher mass a significant fraction remains gravitationally bound. Since it does not stick, a nearly isotropic H 2 exosphere results as shown in Figure 6c (despite electron and photon impact dissociation of H 2 ). In the low-sublimation scenario presented in Figure 6, H 2 O is primarily ejected from polar regions with energies comparable to escape energy. Therefore H 2 O tends to populate high altitudes with column densities peaking near the regions of sputtered production as indicated in Figure 6d.

III.2 The role of collisions

As suggested by [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF], Ganymede's atmosphere might be locally collisional. We therefore investigated the effect of momentum transfer collisions on the global distribution of the exosphere.

Since the implementation of collisions, requires of the order of a few weeks on 128 CPUs, we only ran simulations on one phase angle. Starting from the steady-state solution calculated in the collisionless regime (as described in section III.1), we simulated the fate of Ganymede's exosphere, including collisions (section II.5), during a time period long enough (few tens of hours) to reach an accurate description of the exosphere up to 1200 km in altitude. In Figure 7, we compare the density profile of the two main components in Ganymede's exosphere with and without taking into account the effects of collisions at a phase angle of 90°.

As shown in Figure 7, the main effect of collisions on the density profile is to enhance the population of O 2 at high altitudes. Indeed, most of the collisions occur either between O 2 molecules or between O 2 and H 2 O molecules, the two main species in the collisional region near the surface (within the first tens kilometers). In this layer, collisions between other thermal O 2 molecules have a negligible effect on the O 2 radial density distribution. However, O 2 collisions with the more energetically sputtered H 2 O possessing a high energy tail (1/E²) can add significant energy to the O 2 molecules, thereby expanding the O 2 exosphere.

A C C E P T E D M

A N U S C R I P T We find that our results are in good agreement with the conclusions by Shematovich (2016) model of Ganymede's exosphere. This author found that collisions caused the 10 to 100 km altitude region to be populated by both O 2 and H 2 O. Similarly, [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF] suggested that charge exchange between O 2 + ion and O 2 molecules could enhance the population of O 2 in the upper exosphere.

Charge exchange was also investigated by [START_REF] Dols | Europa's atmospheric neutral escape: importance of symmetrical O2 charge exchange[END_REF]

III.3 The role of sublimation

For the [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] parametrization of sublimation, the total ejection rate would be equal to 3.5×10 29 H 2 O/s (comparable to the 8.0×10 27 H 2 O/s sputtered and to the 7×10 29 H 2 O/s sublimated in [START_REF] Plainaki | The H 2 O and O 2 exospheres of Ganymede: the result of a complex interaction between the jovian magnetospheric ions and the icy moon[END_REF]). This implies that the H 2 O exosphere would be dominated by sublimation. This is clearly illustrated by simulation results in Figure 8a which displays the exospheric density in the equatorial plane, with the same format as in Figure 4a. The H 2 O density peaks as expected around the subsolar region with a clear shift towards the afternoon due to the shift of the maximum temperature on the dayside (Figure 1b). The eclipse period is particularly dramatic, since as shown by Turc et al. (2014), the quick decrease of the surface temperature (Figure 1a For a large sublimation rate, the typical LOS H 2 O column density (as seen from the Sun) would peak around the subsolar region with values between 3.1 and 8.5×10 15 H 2 O/cm² (Figure 8b). The H 2 O sputtered from polar regions is also apparent in the Figure 8b. The slight asymmetry of the column density with respect to the subsolar point in Figure 8b is due to the maximum surface temperature occurring in the afternoon. Since dissociation of H 2 O is a minor source of H, the column density of H is very close (typically 4×10 10 H/cm²) to that for the low sublimation rate case (in Figure 8b). [START_REF] Barth | Galileo ultraviolet spectrometer observations of atomic hydrogen in the atmosphere of Ganymede[END_REF] reported a H column density of 9.21×10 12 cm -2 in Ganymede's exosphere from the measured Lyman α emission of 0.56±0.2×10 3 Rayleigh measured by Galileo during its G1 flyby (at a phase angle of 180° and a solar zenith angle of 60° with a latitude of 21°). Such emission intensity was later confirmed with HST [START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF]. The simulated H column density at that phase angle is on the order of 6×10 10 H/cm 2 in the polar region and is extended down to latitudes of 20° because of the particular shape of the open field line boundary at that particular phase angle. Our calculated column density is therefore two orders of magnitude lower than the estimated column density by [START_REF] Barth | Galileo ultraviolet spectrometer observations of atomic hydrogen in the atmosphere of Ganymede[END_REF]. A two order of magnitude error in the efficiency sputtering H atoms is difficult to imagine based on measured sputtering yields [START_REF] Johnson | Composition and Detection of Europa's Sputter-Induced Atmosphere[END_REF]).

IV Discussion and Comparison with observations

In combination with resonant scattering, H 2 O dissociated via electron impact could be an additional channel producing Lyman α emission. Considering an electron with an energy near 20 eV, the typical cross section for the production of Lyman α by electron impact on H 2 O is 0.43×10 -18 cm 2 [START_REF] Makarov | Kinetic energy distributions and line profile measurements of dissociation products of water upon electron impact[END_REF]). While trace ~keV hot electrons are a negligible contribution, an accurate knowledge of the

A C C E P T E D M

A N U S C R I P T 30 core electron energy is particularly important since this cross section increases by one order of magnitude at 50 eV and peaks at 7×10 -18 cm 2 at 100 eV. Using the cross section at 20 eV, the typical column density of 10 14 H 2 O/cm 2 in the low sublimation scenario (Figure 6) and 10 16 H 2 O/cm 2 in the [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] sublimation scenario would lead to an emission of few Rayleigh and few hundred Rayleigh respectively. The [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF] sublimation scenario would therefore be in rough agreement with Galileo-observed Lyman α emission. The low sublimation scenario could fit the observation for 40 eV core electron temperatures if the sputtering yield was ten times larger. A further potential source of Lyman alpha emission from Ganymede's exosphere is electron impact dissociation of H 2 , also suggested by [START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF]. [START_REF] Ajello | Line profile of H Lyman α from dissociative excitation of H 2 with application to Jupiter[END_REF] published a cross section for the production of Lyman alpha emission from electron impact dissociation of H 2 peaking at 7×10 -18 cm 2 at 60 eV and equivalent to 5×10 -18 cm 2 at an electron energy of 20 eV. Using the typical 2×10 14 H 2 /cm 2 column density we simulated, this process would contribute less than a few tens of Rayleigh.

The typical column density of the O 2 exospheric component has been reported by [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF] as being between 2.4×10 14 and 1.4×10 15 O 2 /cm² when observing Ganymede's trailing hemisphere in June 1996 (phase angle around 270°) and by [START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF] as being equal to 0.3 to 5×10 14 (October 30, 1998 with a phase angle between 290° and 300°). [START_REF] Mcgrath | Aurora on Ganymede[END_REF] also showed the spatial distribution of the auroral emission at different orbital positions characterized by localized bright regions with peak brightnesses between 100 to 400 Rayleigh in good agreement with the previous reported observations. As stated in Section II.4, the O 2 sputtering yield was tuned in order to reproduce the reported column density. We therefore simulated column densities varying from 5 to 6.5×10 14 O 2 /cm² averaged over the disk as observed from the Sun, and peaking to values of 2×10 15 O 2 /cm² at the equatorial limb (Figure 6a). Figure 3) as the primary emission mechanism. We did not attempt to reproduce the ring shape emission in [START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF], [START_REF] Mcgrath | Aurora on Ganymede[END_REF][START_REF] Saur | The search for a subsurface ocean in Ganymede with Hubble Space Telescope observations of its auroral ovals[END_REF]. Furthermore, we did not consider an accurate description of the electric currents in the open field line regions as in [START_REF] Payan | Uncovering local magnetospheric processes governing the morphology and variability of Ganymede's aurora using three-dimensional multiuid simulations of Ganymede's magnetosphere[END_REF]. The electron population used to calculate the emission brightness is the same as earlier:

an electron number density of 70 cm -3 and average energy of 20 eV. The emission rate coefficient for the production of OI ( 5 S o ) atoms from O 2 emitting photons at 135.6 nm wavelength is set to be 1.1×10 -9 cm -3 s -1 as in [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF]. The 135.6 nm emission produced by direct excitation of the O atoms by electron impact is also included, but is much smaller than the O 2 source, in good agreement with [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF]. The goal of our calculation is to look for any potential signature in the observed auroral emission induced by the dusk/dawn asymmetry of the O 2 exosphere (Figure 4). As shown in Figure 9, the simulated emission intensities vary between 10 and 200 Rayleigh in good agreement with the observed values. More interestingly, at a phase angle of 90° (panel a), Ganymede's dusk exosphere is simulated as being up to two times denser than at dawn consistent with the column density of the O 2 exosphere as seen from the Sun being 2.3 times larger on moon's esatern/dusk hemisphere with respect to the column density on the western/dawn hemisphere. This dusk/dawn asymmetry is visible in the simulated emission brightness, the dusk limb emission at midlatitudes being significantly larger than the dawn limb emission. Remarkably, a similar asymmetry between the dawn and dusk limbs was reported by HST (panel a of Figure 2 in [START_REF] Mcgrath | Aurora on Ganymede[END_REF].

At a phase angle of 270° (Figure 9, panel b), the simulated dawn/dusk asymmetry is much smaller (the west/dawn to east/dusk hemispheric average column density ratio begin equal to 1.3), and does not induce a significant emission asymmetry (Figure 9b), as is also seen with HST observations (panel c of Figure 2 in [START_REF] Mcgrath | Aurora on Ganymede[END_REF]. As indicated by HST observations near the trailing hemisphere (Fig. 9b) the auroral emissions occur at higher latitudes than at leading (Fig. 9a) since we are facing Using the emission rate coefficient for the 130.4 nm emission, the typical ratio between 135.6 nm and 130.4 nm is equal to 1.8 with decreasing values at the edge of the emission region where the O/O2 relative abundance tends to increase. [START_REF] Hall | The far-ultraviolet oxygen airglow of Europa and Ganymede[END_REF] and [START_REF] Feldman | HST/STIS ultraviolet imaging of polar aurora on Ganymede[END_REF] reported ratios between 1.2 and 3.2. [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF] suggested that the radiolysis efficiency to sputter O 2 molecules from an icy surface might be between 50-to-300 times less efficient than predicted by laboratory experiments [START_REF] Teolis | Cassini finds an oxygen-carbon dioxide atmosphere at Saturn's icy moon Rhea[END_REF] or as calculated by [START_REF] Cassidy | Magnetospheric ion sputtering and water ice grain size at Europa[END_REF]. In this paper, we scaled this efficiency to reproduce particles/s for the thermal component (Leclercq et al. 2016) and ~10 24 particles/s for the high-energy component [START_REF] Fatemi | On the formation of Ganymede's surface brightness asymmetries: Kinetic simulations of Ganymede's magnetosphere: plasma interaction with Ganymede[END_REF], leading to an average sputtering yield of 100. [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF] estimated the yield for ion sputtering on the icy satellites Rhea or Dione to be between 0.01 and 0.03, 4 orders of magnitude less than the yields inferred in this paper. The temperature difference between Rhea, Dione, and Ganymede could only account for a factor of 10 larger. Since an ion precipitation rate that is two orders of magnitude larger is not realistic, we note that the Jovian plasma is more energetic with respect to the Saturnian plasma [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF]. Such a difference in energy could also lead to an additional order of magnitude larger efficiency at Ganymede than at Dione or Rhea, leaving a difference of two orders of magnitude with respect to the conclusions of [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF]. This difference is re-enforced by noting that two orders of magnitude smaller O 2 column density would place atomic products below the detection limit of FUV oxygen aurorae observations by HST. [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF] also suggested that for reproducing Rhea's exosphere evolution, it was needed to introduce the notion of a shallow 1-cm regolith allowing for the diffusion of O 2 into the regolith and its shielding from the impacting ions. The release of the O 2 is then allowed when the regolith rotates into sunlight and induces the diffusion of the O 2 molecules back to the surface where they can be released by sputtering. Applied to Ganymede's case, such a process would increase the reservoir of O 2 available for ejection (by desorption or sputtering) and therefore help reconcile [START_REF] Teolis | Dione and Rhea seasonal exospheres revealed by Cassini CAPS and INMS[END_REF] sputtering efficiency and the efficiency inferred by our modelling. It would imply that more than 99% of the O 2 originally produced by radiolysis accumulates in the regolith over time and is released into the exosphere by thermal desorption and/or direct sputtering.

In such a scenario, the O 2 exosphere would be much more extended in longitude and latitude. We should also expect a much less asymmetric dawn to dusk terminator, in which the O 2 exosphere behaves like the H 2 O exosphere. Moreover, given that O 2 is more efficiently adsorbed and trapped inside the regolith than simulated here, a dawn enhancement would be favored by the release of trapping in the regolith has actually been suggested by ground-based observations peaking at the trailing hemisphere [START_REF] Spencer | Charge-coupled device spectra of the Galilean satellites: Molecular oxygen on Ganymede[END_REF], and confined at the equatorial regions [START_REF] Calvin | latitudinal distribution of O 2 on Ganymede: Observations with the Hubble Space Telescope[END_REF]. It was shown by [START_REF] Johnson | O 2 /O 3 microatmospheres in the surface of Ganymede[END_REF] that O 2 (and O 3 ) bubbles can be trapped as voids in Ganymede's regolith with the caveat that they can be destroyed should the incident ion flux be too high. The absence of O 2 at the poles, as well as the significantly less O 2 at the leading hemisphere/plasma wake might be explained by a more efficient eroding of the trapped O 2 bubbles by the incident plasma.

V Conclusions

In this paper, we present the first model of Ganymede's exosphere taking into account Jupiter's gravitational influence and the evolution of Ganymede's exosphere along its orbit around Jupiter.

This 3D Monte Carlo model takes into account sublimation and sputtering as the main processes of ejection, electron impact ionization and dissociation in the exosphere as well as a detailed interaction between the exosphere and the surface. The role of Ganymede's magnetosphere is described in a simple way which allows us to clearly distinguish signatures associated with the magnetosphere to those associated with the formation and evolution of the exosphere. Ganymede's exosphere appears to be highly stratified, displaying significant structure in terms of density and composition due to the variation of the ejection mechanisms and molecular movement over its orbital period. In particular, our results suggest that the shape of Ganymede's O 2 exosphere at a given orbital position is highly dependent on the production earlier in its orbit. We also showed that the H 2 O exosphere is very different for the high sublimation rate and a low sublimation rate models. Moreover, even if Ganymede's atmosphere is quasi-collisional, collisions are shown to only have a minor effect on the exospheric density distribution.

A C C E P T E D M

A N U S C R I P T 36

In this work, we show that, by observing Ganymede's multi-species, water product exosphere along its orbit, significant information on the way Ganymede interacts with its environment should be accessible:

-The formation of the exosphere is a complex process driven by the impacting plasma and the resulting space weathering of the surface, whose efficiency depends on the surface characteristics (e.g. temperature, composition, structure) and on the Jovian plasma parameters (e.g. flux, energy and spatial distribution at the surface). The spatial distribution for the lighter species of the sputtered exosphere should directly reflect the spatial characteristics of the interaction between Ganymede and Jupiter's magnetosphere. Moreover, the sputtered exospheric temporal variability and composition should be highly dependent on the surface temperature and composition and could therefore be used to constrain these two surface properties.

-The global 3-D spatial distribution of the exospheric density will depend on the influence of Jupiter and Ganymede's gravitational fields. The fate of the ejecta is also species-dependent due to the species dependent surface interactions. As an example, the dusk/dawn asymmetry in O 2 highlighted in this paper depends on the mobility across Ganymede's surface affected by absorption, accumulation and re-emission from the regolith. The observed spectral emission signatures, such as the auroral emissions, could be potentially used to characterize not only the Jovian-Ganymede interaction [START_REF] Saur | The search for a subsurface ocean in Ganymede with Hubble Space Telescope observations of its auroral ovals[END_REF], but also the spatial distribution of the exosphere.

-Ultimately, tracking the evolution of Ganymede's exospheric morphology and content along its orbit should help constrain its atmospheric formation, and the interaction of its magnetospheric plasma with the atmosphere and surface. The Galileo spacecraft showed the dichotomy in the surface albedo between the trailing and leading hemispheres, which suggest that, Ganymede's exospheric source can depend strongly on geographic longitude, an effect we did not consider in this work, but will implement in the future. To better constrain the full structure of the exosphere we suggest full

Figure 1

 1 Figure 1: a: Evolution of the subsolar surface temperature with respect to Ganymede orbital longitude

  based on the[START_REF] Mcgrath | Aurora on Ganymede[END_REF] auroral observations, we map out more precise regions where surface impacts might preferentially occur, using latitudinal limits of the auroral regions corresponding to the footprints of open field lines where Jovian plasma can penetrate through Ganymede's magnetosphere (Figure3).

Figure 3 :

 3 Figure 3: Position of the peak auroral emission brightness as reported by McGrath et al. (2013) in a

  particle is ejected from the surface by sputtering, we proceed in the same way as inTurc et al. (2014). Namely, H 2 O, O, H, and OH are ejected from the surface following an energy distribution described in equation (8) in Turc et al. (2014) and following a Maxwell-Boltzmann distribution at the local surface temperature in the case of O 2 and H 2 (equation 7 in Turc et al. 2014).

Figure 4 :

 4 Figure 4: Vertical column densities (log 10 cm -2 ) of the H 2 O (panel a) and O 2 exospheres (panel b) at

Figure 5 :

 5 Figure 5: West longitude -latitude maps of the rate of O 2 ejection (in O 2 /cm 2 /s) by sputtering at

Figure 6 :

 6 Figure 6: Reconstructed column densities (in log 10 cm -2 ) of the O 2 exosphere (panel a), of the H

Figure 7 :

 7 Figure 7: Density profiles simulated by EGM without collision (blue lines) and with collision (red lines)

  and by Luchetti et al. (2016) at Europa who both suggested that charge exchange between O 2 + and O 2 might be the dominant loss process for that satellite. Including charge exchange accurately requires coupling between a magnetospheric and an exospheric model, which is work in progress. The near surface density has a typical scale heights for O 2 of 50 km which is significantly larger than the expected scale height, 27 km, for an equilibrium with the surface at a maximum temperature of 146 K. In this model, O 2 molecules are always ejected from the surface following a Maxwell-Boltzmann distribution defined by the surface temperature. However, as discussed in section III.1, O 2 molecules are particularly sensitive to the non-inertial forces because of their long residence time in the atmosphere. A large portion of the O 2 molecules come from regions far from the local surface leading to a globally uniform layer at the poles. This was illustrated in Turc et al. (2014) where the O 2 density at the poles decreased much faster (with scale height around 22 km) than in the subsolar region.

Figure 8 :29

 8 Figure 8: Panel a: Vertical column densities (log 10 cm -2 ) of the H 2 O exosphere at different Ganymede

  the auroral emission associated with our simulated column density (Figure9), by considering electron excitation by a uniform mono-energetic electron population impacting the sputtered O 2 in the Alfven wings of Ganymede's magnetosphere (open field line regions displayed in

Figure 9 :

 9 Figure 9: Simulated images of the HST observed emission brightness reported in[START_REF] Mcgrath | Aurora on Ganymede[END_REF]

  density and found that the ejected O 2 rate should be around 3.4×10 26 O 2 /s, assuming a flux of Jovian ions precipitating on Ganymede's surface between 3.3 and 3.7×1024

  sunrise and a decrease in reservoir content with increasing solar zenith angle. O 2

is reproduced in Table 1 below.

  R I P T 6 surface, all H 2 O products, namely H 2 O, H 2 , H, O 2 , O and OH. The list of the reactions taken into account in our model is given in Table 2 of Turc et al. (2014) and

	Reactions	Rate (s -1 )	Excess Energy (eV)	References
	H + hν  H + + e	4.5×10 -9	3.8	Huebner et al. (1992)
	H 2 + hν  H + H	8.8×10 -9	2.2	Marconi (2007)
	H 2 + hν  H 2 + + e	3.1×10 -9	6.9	Huebner et al. (1992)
	H 2 + hν  H + + H + e	6.9×10 -10	26	Huebner et al. (1992)
	O + hν  O + + e	1.5×10 -8	24	Huebner et al. (1992)
	OH + hν  O + H	9.7×10 -7	3.4	Marconi (2007)
	OH + hν  OH + + e	1.6×10 -8	22	Huebner et al. (1992)
	H 2 O + hν  H + OH	5.2×10 -7	3.4	Marconi (2007)
	H 2 O + hν  H 2 + O	3.8×10 -8	3.4	Marconi (2007)
	H 2 O + hν  H + H + O	4.9×10 -8	4.6	Marconi (2007)
	H 2 O + hν  H + OH + + e	3.8×10 -9	21	Huebner et al. (1992)
	H 2 O + hν  H 2 + O + + e	5.2×10 -10	38	Huebner et al. (1992)
	H 2 O + hν  OH + H + + e	1.0×10 -9	28	Huebner et al. (1992)
	H 2 O + hν  H 2 O + + e	2.1×10 -8	14	Huebner et al. (1992)
	O 2 + hν  O + O	2.0×10 -7	1.3	Marconi (2007)
	O 2 + hν  O 2 + + e	3.0×10 -8	18	Huebner et al. (1992)
	O 2 + hν  O + O + + e	8.4×10 -9	26	Huebner et al. (1992)
	Reactions	Rate (cm 3 s -1 ) Excess Energy (eV)	References
	H + e  H + + e + e	9.1×10 -8	3.8	Ip (1997)
	H 2 + e  H + H + e	9.6×10 -9	4.6	Marconi (2007)
	H 2 + e  H 2 + + e	1.6×10 -8	6.9	Ip (1997)
	H 2 + e  H + + H + e + e	9.6×10 -10	26	Ip (1997)
	O + e  O + + e + e	2.0×10 -8	24	Ip (1997)
	OH + e  O + H	1.2×10 -8	3.4	Ip (1997)
	OH + e  OH + + e + e	2.8×10 -8	22	Ip (1997)
	H 2 O + e  OH + H + e	3.7×10 -8	3.4	Marconi (2007)
	H 2 O + e  H 2 + O + e	1.6×10 -8	3.4	Ip (1997)
	H 2 O + e  OH + H + + e + e	4.3×10 -9	28	Smyth and Marconi (2006)
	H 2 O + e  H + OH + + e + e	4.0×10 -9	21	Ip (1997)
	H 2 O + e  H 2 + O + + e + e	7.1×10 -9	38	Ip (1997)
	H 2 O + e  H 2 O + + e + e	2.1×10 -8	14	Ip (1997)
	O 2 + e  O + O + e	1.3×10 -8	1.3	Marconi (2007)
	O 2 + e  O 2 + + e + e	2.0×10 -8	18	Ip (1997)
	O 2 + e  O + O + + e + e	1.1×10 -8	26	Smyth and Marconi (2006)

Table 1 :

 1 Photon and electron ionization and dissociation in Ganymede's atmosphere (Turc et al.

	2014)

  discussion in Section IV). Based on such observations, we can infer an effective yield of O 2 molecules and then derive corresponding effective yields for the other water products H, O, OH, H 2, and H 2 O from their known relative yields with respect to O 2 . 2 , since they dominate the product composition and are both formed by radiolysis of the icy surface. H, O, and OH yields should be small with respect to O 2 , H 2 and H 2 O yield and are set to 1/20 th of the O 2 yield in agreement with[START_REF] Smyth | Europa's atmosphere, gas tori, and magnetospheric implications[END_REF].

	These yields represent a weighted average over all charge states of O and S ions sourced in Io's

[START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF]

considered relative yield values equal to 1 for H 2 O, 1/18 for H, O, OH, to 16 for H 2 and to 8 for O 2 , whereas

[START_REF] Smyth | Europa's atmosphere, gas tori, and magnetospheric implications[END_REF]

suggested values equal to 1 for H 2 O, 1/17 for H, O, OH, to 2 for H 2 and to 1 for O 2 at Europa. At Europa,

[START_REF] Shematovich | Surface-bounded atmosphere of Europa[END_REF]

considered that 1 molecule of H 2 O should be ejected for each 1/10 O 2 molecules. Noting the wide range in previous models, we choose to base our yields on

[START_REF] Cassidy | Magnetospheric ion sputtering and water ice grain size at Europa[END_REF]

, who made the first detailed calculation of sputtering rates at Europa based on laboratory experiments, Europa's expected surface temperature, and the expected intensity and spatial distribution of the impacting ions flux. These authors estimated that 2×10 27 H 2 O/s and 10 26 O 2 /s should be ejected from Europa's surface. In a companion paper, we applied the EGM model to Europa's atmosphere

[START_REF] Oza | Origin and Evolution of Europa's Oxygen Exosphere[END_REF]

) and thereby constrained Y(T 0 ) for O 2 and H 2 O from the

[START_REF] Cassidy | Magnetospheric ion sputtering and water ice grain size at Europa[END_REF]

calculations of the total rates ejected from Europa's surface. According to

[START_REF] Oza | Origin and Evolution of Europa's Oxygen Exosphere[END_REF]

, each time a molecule of H 2 O is ejected from the surface, an average of 1/20 O 2 is ejected from Europa's surface. The relative yield for H 2 should be nearly stoichiometric, i.e. twice that of O plasma torus, which bears little compositional differences between 9.4 and 15 Jovian radii

(Delamere et al., 2005)

.

  9×10 -15 cm -2 for H 2 O-H 2 O collision or to 3.9×10 -15 cm -2 for O 2 -O 2 collision. We then randomly select a number of test-particle pairs equal to the expected maximum number of collisions that should occur according to

	equal to 6.
	M A N U S C R I P T
	A C C E P T E D
	Typical collision cross sections for molecule-
	molecule or molecule-atom collision at 1 eV are 4.8×10 -15 cm -2 and do not depend neither
	15

significantly on energy between 0 to few eVs nor on the species considered in this work (H, H 2 , O, OH, O 2 and H 2 O). These cross sections are close to the ones used by

[START_REF] Marconi | A kinetic model of Ganymede's atmosphere[END_REF]

, as an example,

  21 H 2 O/s whereas the average ejected flux by sputtering is equal to 8.0×10 27 H 2 O/s, 3.6×10 26 O 2 /s, 7.2×10 26 H 2 /s and 1.0×10 25 H, O and OH/s. Marconi (2007) and Turc et al. (2014) used a flux of 1.5×10 26 H 2 O/s and 1.2×10 27 O 2 /s for sputtering, whereas Plainaki et al. (2015) used a flux of 7×10 25 H 2 O/s for sputtering and 7×10 29 H 2 O/s for sublimation using the same model for sublimation as in Marconi (2007) and Turc et al. (2014) (green square symbols in Figure 2).We simulated more than 4.5 rotations of Ganymede around Jupiter in two weeks on 64 CPUs (that is, 3.7×10 6 time steps of 0.75 s each for a total of 780 simulated hours). During that run, 2.6×10 7 H 2 O test-particles (simulating a total of 2.2×10 34 H 2 O molecules) and 2.1×10 7 O 2 test-particles (for 10 34 O 2 molecules) were ejected from the surface. Of these 1.2×10 6 H 2 O test-particles escaped the simulation domain, 2.4×10 7 were reabsorbed in the surface, 6×10 5 were dissociated or ionized and none were merged (meaning that the maximum number of test-particles in a cell allowed by the simulation was never reached for H 2 O). For O 2 , 1.3×10 6 were dissociated or ionized, 1.9×10 7 were merged. The 3×10 5 O 2 test-particles which escaped the simulation domain were generated from the tail of the energy distribution ejected from the surface. Being associated with very small weight, these escaping molecules represented an escaping flux of 6×10 17 O 2 /s. Merging test-particles is much more important for O 2 than for H 2 O, because O 2 test-particles are essentially confined near the surface. At

	each time, a total number of 1 to 2×10 6 test-particles are followed. After convergence towards a
	steady state, the typical escape rates are 1.6×10 27 H 2 O/s, 6×10 17 O 2 /s, 1.7×10 26 O/s, 1.2×10 26 OH/s,
	3.0×10 26 H 2 /s and 2.3×10 26 H/s.

Acknowledgement: This work is part of HELIOSARES Project supported by the ANR (ANR-09-BLAN-0223) and ANR MARMITE-CNRS (ANR-13-BS05-0012-02). Authors are also indebted to the "Système Solaire" program of the French Space Agency CNES for its support. Authors also acknowledge the support of the IPSL data center CICLAD for providing us access to their computing resources and data. Data may be obtained upon request from F. Leblanc (email:francois.leblanc@latmos.ipsl.fr).

Vasyliunas, V.M., Eviatar, A., Outflow of ions from Ganymede: A reinterpretation, Geophys. Res. Lett. 27, 1347Lett. 27, -1349Lett. 27, , 2000. . Vorburger A., Wurz P., Lammer H.n Barabash S. and O. Mousis, Monte-Carlo simulation of Callisto's exosphere, Icarus 262, 14-29, 2015