Celia Schneebeli
email: celia.schneebeli@u-bourgogne.fr

GIFs IN ONLINE WRITTEN INTERACTION: EMBODIED CUES AND BEYOND

Keywords: GIF, gesture, CMC, interaction, utterance

Recent research has characterized

GIFs as a means for speakers to reproduce nonlinguistic cues such as bodily actions and facial expressions in online written interaction, which is first and foremost text-based. They make it possible to translate in digital interaction what those nonlinguistic cues enable to do in face-to-face conversation: express emotion and affect, and elaborate on what is being said. This article explores further the role and functions of GIFs as embodied cues and goes beyond, where GIF use diverges from body language in face-to-face conversation.

Introduction

The GIF, acronym for Graphics Interchange Format, is an image file format that was developed by American online service provider CompuServe, and released in 1987. It uses a technique that enables lossless data compression, which makes it possible to obtain images with a reduced file size without impacting their visual quality. More importantly, contrary to other popular image formats such as the JPEG or PNG, it also enables to show looping sequences "within the same image file without being the size (or resolution) of a video" (Miltner and Highfield 2017: 3). So much so that, in the Web 2.0, the word GIF is commonly used to refer not to the file format itself but to a short looping video composed of a sequence of images taken from a film, a TV series, a cartoon or any television programme. In this paper, the word GIF is likewise taken in this metonymical sense of the word.

In a 2016 article, Jackson Tolins and Patrawat Samermit focus on the proximity of the new affordances provided by GIFs and the role of gesture and facial expression in face-to-face dialogue. According to them, GIFs offer a way for users to express bodily displays in text-based interaction: they constitute "embodied reenactment" thanks to which users can reproduce their own nonverbal behaviour in computer-mediated conversation (CMC). This is why Tolins and Samermit try to approach the function of GIFs in text-based virtual communication within the same framework as face-to-face dialogue. GIF use is analysed as an expansion and extension of a pre-existing practice that already exists outside of technology-mediated interaction. Since it "both mirror[s] faceto-face dialogue and extend[s] beyond it" (Tolins and Samermit 2016:75), GIF use seems to be a good example of those "phenomena that adapt to and are reconfigured by Web 2.0" (Herring 2013: 16). This perspective, the one of adaptation and reconfiguration, is the one that this article intends to adopt in order to study how GIFs are used in online written interaction. Starting from Tolins and Samermit's point that GIFs are used as stand-ins for speaker's nonverbal behaviours in computer-mediated conversation, it aims to go beyond and show that GIFs are used in other ways than what happens in face-to-face conversation. In doing so, the main questions this study tackles are the following: What are the functions of GIFs in online conversation? How do they interact with verbal cues? Where does their use diverge from what body language enables to do in face-to-face conversation? The approach is that of pragmatics, discourse analysis, and more precisely of computer-mediated discourse analysis (CMDA).

To answer those opening questions, the study relies on a qualitative and quantitative analysis of a sample of messages posted on the public Facebook page of an American ice hockey team, the Philadelphia Flyers. The dataset was collected over 3 days in October 2018. It is composed of 5 comment threads, which total 1351 messages, of which 90 use GIFs. While it is difficult, if not impossible, to collect a sample that would be representative of what users do on Facebook (it may indeed vary according to the type of page, its theme or audience), pages devoted to popular sports team have the benefit of having a wide audience and an active fan base, and thus present posts with relatively numerous comments (the selected threads have from 134 to 495 comments each). All the initial messages of the chosen threads deal with different topics and were posted in different contexts. Therefore, the dataset comprises diverse reactions in different contexts, and present multiple uses of GIFs.

Facebook is in itself a relevant choice to perform such a study since it belongs to what Susan C. Herring calls « interactive multimodal platforms (IMPs) », i.e. « Web 2.0 platforms that support a convergence of channels or 'modes ' (text, audio, video, images) for user-to-user communication » (Herring 2016 : 398). It lets users insert emoji, videos, stickers, images and GIFs, which is the latest addition (it dates back to 2017) to the modalities that can be used on the platform. This "use of channels other than text, and semiotic systems other than verbal language, to carry on 'conversational' exchanges" is considered by Herring as "[a]nother emergent Discourse 2.0 phenomenon" (Herring 2013: 16). However, as she also notes, in spite of those non-verbal modalities having been made available, discourse on Web 2.0 platforms such as Facebook remains essentially text-based. This makes it all the more relevant to question the function of supposedly nonverbal modalities such as GIFs, and the interaction between verbal and nonverbal modalities on those platforms.

GIFs and the "performance of affect"

In recent studies, the function of GIFs in computer-mediated discourse has been linked to the expression of emotion and affect, be it called "affective responses" (Tolins and Samermit 2016: 75) or "performance of affect" (Miltner and Highfiled 2017: 1). For Miltner and Highfield, this is even "one of the most common uses of GIFs -if not the most common" (2017: 4). This claim lies on the link between the gestures GIFs enable to picture and the function of gesture in communication in general.

GIFs as stand-ins for gestures in computer-mediated interaction

Modern gesture studies [START_REF] Kendon | Gesture: Visible Action as Utterance[END_REF][START_REF] Mcneill | Hand and Mind[END_REF] has shown that gesture, in face-to-face conversation, may express "substantive meaning" (Kendon 2017: 161). When they are undertaken deliberately, they can be viewed as visual utterances that may express propositional content and pragmatic content. Among other things, utterances may then express feelings and emotions. For example, in certain contexts (and cultures, since gesture-meaning is strongly conventionalized and culture-dependent) shaking one's head can express disappointment, or a clenched fist may express anger. It is then quite simple to understand why GIFs are used to perform affect: as GIFs enable to picture people or characters performing those gestures, they have become a way to perform those gestures by proxy. In the text-based environment of Facebook, GIFs enable to reproduce what gestures, bodily actions and facial expressions enable to do in face-to-face conversation. To do so, the user just has to select a GIF showing a person or a character performing the bodily display he or she would have performed in face-to-face conversation, or an elaborated version of it. The use of what could be called "gestural GIF" as affective utterance is quite frequent in the dataset. In 50 GIFs out of the 90 used in the five threads, the gesture or facial expression pictured can be directly linked to the expression of an emotion. Above are some examples: All these GIFs show conventionalized expressions of affect through bodily actions, which are the same in real life and virtual communication: the now famous facepalm expresses dismay, rolling one's eyes expresses exasperation, eyes full of tears express sadness, and breaking something with an angry face (as the hockey player above who breaks his stick on the boards) expresses rage and frustration. The gesture or facial expression performed in these GIFs is therefore analysed by the co-speaker as a meaningful "visual and affective response toward the contents of the prior turn" (Tolins and Samermit 2016 : 81). Just as gestures in face-to-face interaction, GIFs may then constitute a single turn in a conversation: no further speech is needed. The GIF itself expresses the affective response. The speaker can also choose to couple gesture and speech. This is the case in 38 messages in the dataset. All of them present cases of GIFs including text inside the image itself. There is no case, on the contrary, of a composite turn where a message including a GIF would be coupled by a speaker with another message containing text only.

Gestural GIFs and speech

Kendon and other gesture analysts rightly point out that a gesture may constitute an utterance in itself, or may be used in combination with speech: "At times [gestures] are used in conjunction with spoken expressions, at other times as complements, supplements, substitutes or as alternatives to them" (Kendon 2004: 1). Tolins and Samermit similarly divide gestural GIFs between standalone GIFs that express an affective response on their own, and GIFs used in composite turns composed of one or more textual messages followed by a GIF that represents or elaborates on the textual content. There are no composite turns in the dataset, but there are many GIFs containing text, a case which is a bit different, and curiously enough, is not analysed by [START_REF] Tolins | GIFs as Embodied Enactments in Text-Mediated Conversation[END_REF]. In this case, the interaction between gesture and speech, and the respective roles of speech and gestures, is rather complex and variable. First of all, the gesture or bodily action may have more informational weight than the caption, which could then be removed. This is the case here for example: The man that the GIF shows is crying in his couch, which already evokes sadness, i.e. the emotion the poster may be feeling after reading the initial post that announces that one of the best players of the team has been injured. Why, then, add a caption? It could first be read as a way to add a nuance to what is being expressed. In the case of the man crying, the caption may link sadness to dismay or disappointment. It could also be intended as a way to make the emotion expressed clearer. This is probably the reason why a caption has been added in the following GIF, which reacts to the post announcing the Flyer's severe loss against the Sharks: The caption "For shame" makes it clear that the affective reaction expressed by the cat walking with a paper bag on its head is shame. Here, the caption acts as a form of disambiguation. Adding a redundant caption can also be a way to reinforce the emotion expressed, just as in face-to-face conversation, where a speaker may use at the same time a gesture and an utterance that express the same feeling. This is in fact what the man (a WWE wrestler) in the following GIF does: This GIF is used to react to the post announcing the long awaited victory of the Flyers. It shows the conventional bodily gestures used to rejoice after victory: raised arms. It also contains a "Yes!" caption that is the utterance conventionally associated to the same emotion, which is also what the man's lips seem to be saying. Redundancy between bodily gesture and utterance can be read as a way to reinforce the expression of joy. In fact, of all the gestural GIFs accompanied by a caption, 17 are used as subtitles for an utterance the person or character portrayed is already saying. On the contrary, other users include captions that convey a message which is the opposite of what the bodily action pictured in the GIF say. Non-redundancy between speech and bodily action then creates irony: In the context of the initial message responded to (the one that says one of the star player has been injured), the emotion expressed here is dismay and being jaded. The bodily display, going (back?) to bed with a weary face, corresponds to this affect. The apparently joyful caption is therefore read as ironic and really conveying a negative feeling. In this case, gesture and speech are of equal importance. Both are needed to build the meaning of the utterance. Last but not least, the caption may sometimes be the most important part of the utterance and bodily displays may not be significant in themselves. This is the case in the following post: The exclamation "Oy vey!" expresses dismay or exasperation in Yiddish. The gesture the character makes is not significant or clear enough to be used alone (it could express surprise, for example). The message could have been reduced to the caption as main text, without any image or bodily action.

All these examples show that the relationship between bodily display and caption in GIFs is variable, just as the relationship between gesture and speech is in face-to-face interaction. In terms of content, bodily action and speech may or may not be redundant, with different effects. In terms of informational weight, bodily displays sometimes seem more important than the caption (which could be erased), but sometimes the caption is needed or more important, and sometimes the caption and the image seem of equal weight. This is why the bodily display pictured in GIFs cannot be reduced to an illustration or demonstration of the emotion the text already expresses.

GIFs as visual elaboration

It would be an oversimplification to limit the function of GIFs in computermediated communication to the "performance of affect". Similarly, reducing the content of GIFs to gestures, facial expressions or bodily actions would make the scholar overlook some important aspects and functions of GIFs. Tolins and Samermit try to go beyond the expression of affect by examining a second function, which they call "co-speech demonstrations". In this function, GIFs are coupled with text and are used by speakers to "provide visual elaborations of their own talk" (Tolins and Samermit 2016: 83). However, this function and its limits are less clearly defined by the authors than that of "affective response displays". It seems to encompass all cases where GIFs are used to extend the content of one or more previous message(s), be it to express affect or picture action. A clearer classification would be needed. As it has been written above, in the dataset used for this study, there is no case of composite turns where a GIF would accompany a textual message. However, there are two uses in the dataset that appeal to the same function of visual elaboration, and seem to form two distinct types. One is where the GIF is used to merely illustrate or echo the content or part of the content of a preceding message, whether it is a user comment or the initial message posted on the Facebook page. This usage could be called "illustrative GIFs". The other one is where a GIF is used as a visual metaphor but without an explicit tenor provided in an accompanying message. This second usage we call "metaphorical GIFs". In both cases, the content of the GIF is not limited to picturing a person and/or bodily actions.

Illustrative GIFs

The use of illustrative GIFs, i.e. GIFs used to illustrate or echo what is being discussed is not very frequent in the dataset: there are only three of them. In the five selected threads, it mostly consists of people using GIFs showing Gritty, the new mascot of the Philadelphia Flyers, when Gritty is being referred to by a post or other users: . This GIF, which shows Gritty dancing, is posted as a primary reaction to the initial post that announces that the new mask of the Flyers' goalie features Gritty. It seems to express no particular emotion, bears no argumentative significance with regards to the preceding message or initial post, but clearly corresponds to the topic at hand. Therefore, it seems to be merely illustrative: it depicts what is being discussed in the thread. A previous study suggests that this illustrative role tends to be more frequently fulfilled by static images than GIFs (see [START_REF] Schneebeli | Les modalités iconiques dans le discours médié par ordinateur: du neuf dans l'interaction ?[END_REF], for a comparison between GIFs and static images and examples of static images used as illustrations). Most GIFs used as visual elaboration seem to assume a more complex role of visual metaphor, where the GIF does not just illustrate what is being discussed, but presents a visual elaboration in the form of an analogy.

Metaphorical GIFs

Metaphors consist in describing a thing, person, situation or feeling (the "tenor" of the metaphor) with a word or expression normally denoting another thing, person, situation or feeling (the "vehicle" of the metaphor), so as to suggest some common quality shared by the two. In computer-mediated discourse, GIFs are frequently used to express analogy in a similar way. This is in fact the very essence of all the posts opening with the snowclone sentence "When" + an assertive clause describing a situation of life ("When I have a paper to finish", "When I decide to quit smoking", "When I receive a mean review"…) and followed by a GIF which presents a visual elaboration of the situation described in the initial sentence. This type of post, particularly popular on Twitter and Tumblr, relies on an analogy between the situation described in the sentence and the one showed in the GIF, which is made to be striking and/or funny. It should be noted that the GIF in question does not have to picture a facial expression, a gesture, or any action performed by a person. As a metaphor, it does not have to belong to the same referential domain as the tenor, therefore it may virtually picture anything. In the dataset studied, there is no such use of a metaphorical GIF coupled with a textual message in a composite turn. There are, however, metaphorical GIFs used in single turns. The only difference between the two is that the tenor of the metaphor is never provided by the speaker, although it is easily recoverable in the situation or a preceding post by another speaker. For example, three posts use the "dumpster fire" metaphor. To do so, they use the two following GIFs: These GIFs are found in the comment thread of the post which announces the result of the game where the Flyers lost 8 to 2 against the Sharks of San Jose. They are used to describe how the team plays, or the whole situation of the team which is in the middle of a losing streak. The tenor is provided by the initial post on the page, which announces the Flyer's heavy loss. The dumpster fire metaphor in itself is not only used as a GIF. As a set expression, it has entered the Merriam-Webster dictionary in 2018 with the following definition: "An utterly calamitous or mismanaged situation or occurrence". It is therefore now considered as a set metaphor, which is immediately recognizable, for which the GIF is a visual counterpart. The corresponding GIF enables to pack the metaphor in a single moving image without using words, and even adapt it to the situation at hand since in one of those GIFs, the dumpster in question lies on the ice of an arena. There is no direct equivalent to this in face-to-face conversation.

All metaphorical GIFs do not have to rely on a set metaphor. Any action or situation can be used, as long as the common property between the tenor and the vehicle is clear enough. Another GIF used in the same comment thread shows an archive film footage of a train accident. Again, it elaborates on the situation or performance of the team, with is depicted as a literal catastrophe. A similar image is used to comment on the post which invites users to react to the same disastrous game while it is being broadcast. It shows an electric train which derails and takes fire, leaving its owner speechless. The image is again a visual elaboration of what is happening: while the game had started quietly, it all went quickly wrong for the Flyers. The incapacity of the Flyers to stop the other team from scoring is mocked in another metaphorical GIF, where Homer Simpson receives so many soccer balls that he can't stop them and ends up being buried under a pile of soccer balls.

As can be seen in those examples, metaphorical GIFs often rely on hyperbolic elaborations. Being extracted from films, cartoons, or series, GIFs enable users to give their thoughts a dramatic elaboration, which would not be possible, or at least would be more limited (as far as words and gestures are concerned), in face-to-face conversation. As Miltner and Highfield put it, they have the "capacity to augment and shape our affective performances", depict "the mundane and the ordinary, conveyed through heightened production values, performances, and emotions (Miltner and Highfield 2017: 4).The visual elaboration provided and reconfigured by GIFs can thus be literally larger than life.

Other functions of GIFs

The two concepts that have been explored above, performance of affect and visual elaboration still do not enable to describe wholly how and why GIFs are used in the dataset. This is in part linked to the fact that only the first category expresses a function in the strict sense of the word, i.e. the purpose for which GIFs are used. Visual elaboration in itself is not a function: it does not say anything about what GIFs are meant to express or do. Therefore, it would more aptly be considered a category of GIF use, which can have diverse functions. One of the main functions in question, the expression of emotions and affect, has already been explored. But other functions appear in the dataset, which need to be accounted for: an evaluative function, a perlocutionary function, and another function that Miltner and Highfiled call "demonstration of cultural knowledge" (Miltner and Highfield 2017: 1).

Evaluative function

In conversation, be it online or off line, people often produce utterances with an evaluative function, i.e. in order to attribute a certain value or quality to an object, person or utterance, which positions it on a corresponding scale of worth. As evaluating deals with giving one's opinion, it may sometimes imply affect but must not be confused with expressing affect. The affective and evaluative functions are located on two different levels, even though they may sometimes overlap: the expression of opinion may but does not necessarily imply the expression of emotion. As for GIFs that deal directly with the expression of judgment, they appear to be of two forms. The first is GIFs that explicitly deliver the judgment directly or by proxy: These GIFs have a caption that contains a verbal predicate or adjective that bears an evaluative meaning. It enables the user to express his or her opinion on something, which is their main function. In the first example, the injured player referred to in the initial message is accused of being too "soft". In the other example, the team itself is the object of a very negative evaluation. The evaluative element in the GIF is in those cases the caption more than the visual display. There is no example in the dataset of a GIF that would express evaluation through a bodily movement alone but one could think of a GIF showing someone giving a thumb up or applauding.

The second form of evaluative GIF is where assessment relies on a visual metaphor that gives a positive (flattering) or negative (disparaging) image of the tenor. In fact, most of the metaphorical GIFs quoted above have an evaluative function. The dumpster fire metaphor, which of course always implies a very negative evaluation, is a good example. Another example, used twice in the dataset, is a GIF which shows the Flyers' mascot falling in a ridiculous manner while skating (something associated with people who can't skate). By analogy, users will understand that it is used to criticize the team for not knowing how to play hockey properly.

Perlocutionary function

GIFs can also have a perlocutionary function, in which they are meant to have an effect on the reader or co-speaker. The following example does so: The GIF is used to issue a directive, which has, by definition, a perlocutionary function. This function, however, is essentially due to the caption (an imperative utterance), not to the portrayed action. There is no instance, in the dataset, of a GIF that would deliver an order on its own, just by means of a bodily action. However, it is not impossible to do so. One could think for example of a GIF showing a person motioning someone to go out.

Apart from delivering directives, GIFs may more generally be used to engage or provoke the addressee. In the dataset, some GIFs seem to be used with the sole function of inciting a reaction in the co-utterer or reader. This is the case in the following exchange, which is about the new mascot of the Flyers, Gritty: While the first commenter wishes the Flyers would get rid of Gritty, a second user replies to this post with a single GIF showing Gritty dancing. This can be seen as a humorous provocation that aims at making the first commenter react, which he does, humorously too, by answering "thanks…right before lunch too". Another classic example of a perlocutionary use is when fans of other teams comment on the page of the Philadelphia Flyers with GIFs that are meant to taunt them. There is only one instance of it in the dataset: The GIF shows the mascot of the historical rival team, the Penguins of Pittsburgh, holding the Stanley cup, which is the trophy awarded to the team that wins the post-season tournament. This GIF was posted by a fan of the Penguins in order to annoy and upset Philadelphia fans, even though the initial post has nothing to do with it. It goes without saying that this GIF usage is a lot more frequent in case of games opposing the two rivals, or during post-season games, where teams are trying to win the cup.

Social function

Last but not least, as Miltner and Highfield rightly write, "The selection and presentation of GIFs are also a performance of cultural knowledge. The GIF is not just a proxy for the individual's particular affective or emotional state, but an illustration of the user's knowledge of a certain text or cultural conversation" (Miltner and Highfield 2017: 6-7). Thus, this last function, which they call "performance" or "demonstration of cultural knowledge" is more of a metafunction, a social function which may be present in any of the functions detailed above. Indeed, GIFs do not only enable to judge, taunt, express a reaction or emotion; they also enable to present oneself in interaction (to paraphrase [START_REF] Goffman | The Presentation of Self in Everyday Life[END_REF]. In selecting a GIF, users may then make a statement about themselves. In Goffman's terms, GIFs can be seen as "props" that can be used to give oneself a certain identity in interaction. For example, choosing a GIF from a particular series, film, or cartoon, or portraying a certain actor or character may be a way to position oneself as a fan of the work or person in question. At the same time, it may also create a sense of belonging with other users who are familiar with the same culture or context the GIF comes from. Maybe this function of GIFs could be compared to the phatic function of language: it is not really about the content of the message itself but about bonding, creating complicity with others users who are familiar with the same culture. Of course, this last function must not be idealized: sometimes GIFs are selected for their iconicity (the fact that they perfectly match the intended meaning and effect), not the work they come from or the character or person they picture. Therefore, it must not be taken to be more systematic than it really is.

Conclusion

GIFs enable to translate to digital interaction what bodily displays enable to do in face-to-face conversation: express emotion, provide a visual elaboration of speech, evaluate, provoke or engage the addressee. However, they can't be considered only as an adaptation of body language in online conversation. First, GIFs have the added value of proposing a complex visual representation that makes utterances more dramatic, possibly more appealing (when borrowing from a culture shared by addressees), and which enhances affective performances. What is more, GIFs may also include speech in the form of captions. As a result, GIFs can regroup speech and visual representation in a single multimodal form of expression that may be used in a single turn, without typing more text. For this reason, the GIF may be considered as a synthetic form of utterance that is particularly adapted to Web 2.0 discourse, which is increasingly multimodal (Herring 2013: 21). This conclusion would then call for a re-evaluation of the semiotic nature of GIF, which are often taken to be nonverbal modalities.

Figure 1 :

 1 Figure 1: GIFs showing gesture / bodily movement to express an emotion.

Figure 2 :

 2 Figure 2: GIF expressing sadness / disappointment through bodily action and caption.

Figure 3 :

 3 Figure 3: GIF expressing shame though bodily action and disambiguating caption.

Figure 4 :

 4 Figure 4: GIF expressing joy and triumph through redundant gesture-speech combination

Figure 5 :

 5 Figure 5: GIF expressing dismay though non-redundant (ironic) speech-gesture combination.

Figure 6 :

 6 Figure 6: GIF expressing dismay predominantly though textual channel.

Figure 7 :

 7 Figure 7: Illustrative GIF showing Gritty dancing.

Figure 8 :

 8 Figure 8: Metaphorical GIFs showing a dumpster fire to express how calamitous a situation is.

Figure 9 :

 9 Figure 9: GIFs delivering an explicit judgment.

Figure 10 :

 10 Figure 10: GIF used in order to issue a directive.

Figure 11 :

 11 Figure11: GIF used in order to incite a reaction of the co-utterer.

Figure 12 :

 12 Figure 12: GIF used in order to taunt Philadelphia Flyers fans.


Yongliang Feng

Pierre Boivin

Jérôme Jacob

Pierre Sagaut
email: pierre.sagaut@univ-amu.fr

Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows

Keywords: compressible flows, thermal lattice Boltzmann model, equation

A thermal lattice Boltzmann model with a hybrid recursive regularization (HRR) collision operator is developed on standard lattices for simulation of subsonic and sonic compressible flows without shock. The approach is hybrid: mass and momentum conservation equations are solved using a lattice Boltzmann solver, while the energy conservation is solved under entropy form with a finite volume solver. The defect of Galilean invariance related to Mach number is corrected by the third order equilibrium distribution function, supplemented by an additional correcting term and hybrid recursive regularization. The proposed approach is assessed considering the simulation of i) an isentropic vortex convection, ii) a two dimensional acoustic pulse and iii) non-isothermal Gaussian pulse with Ma number in range of 0 to 1. Numerical simulations demonstrate that the flaw in Galilean invariance is effectively eliminated by the compressible HRR model. At last, the compressible laminar flows over flat plate at Ma number of 0.3 and 0.87, Reynolds number of 10 5 are considered to validate the capture of viscous and diffusive effects.

of state, Mach number, recursive regularized

Introduction

The lattice Boltzmann method (LBM) is an alternative approach to simulate fluid flows based on the Boltzmann equation. The base model is well suited for the simulation of low-Mach, athermal flows [START_REF] Mohamad | Lattice Boltzmann method: fundamentals and engineering applications with computer codes[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF][START_REF] Krüger | The lattice Boltzmann method[END_REF]. Due to its advantages for massively parallel computing as well as its ability to handle very complex geometries, there are significant research efforts invested in extending LBM to thermal and subsonic to supersonic applications.

In achieving that goal, the numerical stability of the Lattice-Boltzmann collision model has long been identified along the key issues [START_REF] Luo | Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations[END_REF]. The most commonly used collision model is the single time relaxation process referred to as the Bhatnagar-Gross-Krook (BGK) model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. In order to overcome the insufficient stability observed in the BGK model, several improved collision models with enhanced stability have been proposed, e.g. multi-relaxationtime models [START_REF] Lallemand | Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF][START_REF] Di Ilio | Simulation of turbulent flows with the entropic multirelaxation time lattice Boltzmann method on body-fitted meshes[END_REF], entropic LB models [START_REF] Prasianakis | Entropic lattice Boltzmann method for simulation of thermal flows[END_REF][START_REF] Frapolli | Entropic lattice Boltzmann simulation of thermal convective turbulence[END_REF], regularized BGK models [START_REF] Latt | Lattice Boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF], or cumulant/cascaded LB models [START_REF] Geier | The cumulant lattice Boltzmann equation in three dimensions: Theory and validation[END_REF][START_REF] Hajabdollahi | Symmetrized operator split schemes for force and source modeling in cascaded lattice Boltzmann methods for flow and scalar transport[END_REF].

Numerous studies have shown the capabilities of LBM in the simulation of thermal flows including scalar transport processes. A hybrid finite difference thermal model using two dimensional multiple relaxation time collision model was presented for low Mach number compressible flows [START_REF] Tölke | A thermal model based on the lattice Boltzmann method for low Mach number compressible flows[END_REF]. As reported in [START_REF] Li | Effect of the forcing term in the pseudopotential lattice Boltzmann modeling of thermal flows[END_REF], the hybrid finite difference thermal model can simply cancel a spurious source term in thermal lattice Boltzmann models with forcing terms. Moreover, the robustness of hybrid approach was demonstrated in simulations of highly compressible flows [START_REF] Nie | Lattice Boltzmann/finite-difference hybrid simulation of transonic flow[END_REF].

For the simulation of compressible flows at higher Mach numbers, how-ever, no clear consensus has been reached within the LB framework. The classical method requires considering high-order moments of density distribution functions to recover the macroscopic energy conservation equation for thermal lattice Boltzmann model [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF], usually relying on extended-neighbour lattice set (D1Q5, D2Q37, D3Q121 etc.) leading to LB methods referred to as multi-speed models [START_REF] Philippi | High-order lattice-Boltzmann[END_REF][START_REF] Frapolli | Entropic lattice Boltzmann models for thermal and compressible flows[END_REF]. This approach, seems rather expensive for industrial applications due to the high number of lattices (which may become prohibitively high in multi-component flows).

Several attempts have also been made keeping lower-order lattices (D1Q3, D2Q9, D2Q17, D3Q19, D3Q39 etc.). A two-dimensional compressible LB model with reduction of lattice sound speed through incorporating an athermal equation of state was proposed by Yu and Zhao [START_REF] Yu | Lattice Boltzmann method for compressible flows with high Mach numbers[END_REF], in which a forcing term was added in the form of the pressure gradient calculated by D2V17 lattice nodes. Next, the LB model using reduction of sound speed [START_REF] Buick | Investigation of a lattice Boltzmann model with a variable speed of sound[END_REF] and hybrid finite difference entropy equation was studied in D3Q39 lattice model for transonic flows [START_REF] Nie | Lattice Boltzmann/finite-difference hybrid simulation of transonic flow[END_REF]. As reported in studies of Kupershtokh et al. [START_REF] Kupershtokh | On equations of state in a lattice Boltzmann method[END_REF], both the method of incorporating forcing terms into the lattice Boltzmann model and calculation of pressure gradient are crucial to implement an arbitrary equation of state.

For the nearest-neighbour lattice sets (D1Q3, D2Q9, D3Q15, D3Q19, D3Q27, i.e. the so called standard lattices), Házi and Kávrán [START_REF] Házi | On the cubic velocity deviations in lattice Boltzmann methods[END_REF] suggested to use a third-order Hermite polynomial expansion the equilibrium distribution functions to address acoustic problems. Dellar [START_REF] Dellar | Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices[END_REF] proposed a lattice Boltzmann algorithm with Hermite-Gram-Schmidt polynomials to reduce cubic defects in Galilean invariance on standard lattices for the simulation of low-subsonic flows with uniform density. Malaspinas [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF] developed a recursive regularization procedure to increase stability and accuracy of the lattice Boltzmann scheme for athermal compressible flows at moderate Mach numbers. For thermal compressible flows, a coupled double distribution function (DDF) lattice Boltzmann model was developed for thermal flow with correcting terms in two-dimensional space of multiple relaxation time collision model [START_REF] Li | Couplling lattice Boltzmann model for simulation of thermal flows on standard lattices[END_REF]. A three-dimensional DDF thermal lattice Boltzmann model with general correcting term for thermal flows with variable density was developed by Feng et al. [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF]. A DDF thermal lattice Boltzmann model was extended to high-subsonic flows using finite volume approach to improve numerical stability [START_REF] Feng | A compressible lattice Boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF].

The collision models on standard lattice has also been investigated in compressible and weakly compressible flows. An entropic thermal D2Q9 model was introduced in [START_REF] Prasianakis | Lattice Boltzmann method for simulation of compressible flows on standard lattices[END_REF] for subsonic flows characterized by large temperature and density variations. Recently, the central moment and cascaded thermal lattice Boltzmann models with variable density was proposed for for weakly compressible flows [START_REF] Safari | A lattice Boltzmann model for thermal compressible flows at low mach numbers beyond the Boussinesq approximation[END_REF][START_REF] Fei | Cascaded lattice Boltzmann method for thermal flows on standard lattices[END_REF]. Besides, the recursive regularized BGK model was proposed for high Mach number flows basing on higher-order lattice set [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF]. Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF] then proposed a dynamic hybrid recursive regularized BGK model with self-adaptive dissipation for Large-Eddy Simulation high Reynolds number flows and Reynolds-Averaged Numercial Simulation [START_REF] Wilhelm | An explicit power-law-based wall model for lattice-Boltzmann-method-Reynolds-averaged numerical simulations of the flow around airfoils[END_REF], which could be generally extended to compressible flows by considering the effects of turbulent fluctuations have a striking resemblance to those of microscale (thermal) fluctuations in laminar flows [START_REF] Chen | Expanded analogy between Boltzmann kinetic theory of fluids and turbulence[END_REF][START_REF] Feng | A lattice-Boltzmann model for low-Mach reactive flows[END_REF].

Due to its intrinsic low dissipation, high fidelity and simplicity compared with the conventional finite volume/difference approach [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF][START_REF] Wang | A comparative study of LBE and DUGKS methods for nearly incompressible flows[END_REF], it is natural to seek a stable and efficient LB model able to simulate fully compressible flows on standard lattices. The objective of this paper is to present a new thermal lattice Boltzmann model relying on the nearest-neighbour lattice structures with a hybrid recursive regularized BGK collision model for subsonic to sonic smooth compressible flows. It is organized as follows: Section 2 presents the governing equations including the energy conservation in entropy form; Section 3 presents the thermal lattice Botlzmann method with a hybrid recursive regularized BGK collision model along with the numerical algorithm; In section 4, simulation of an isentropic vortex convection, a two dimensional acoustic pulse and non-isothermal Gaussian pulse with Ma number between 0 and 1 and a compressible laminar flow over flat plate at M a = 0.87 are computed as validation and application. Finally, Section 5 draws conclusions and perspectives.

Governing equations

The mass and momentum conservation equations governing compressible flows read

∂ρ ∂t + ∇ • (ρu) = 0, ( 1a
)
∂ρu ∂t + ∇ • (ρuu) = -∇p + ∇ • τ (1b)
where ρ, u, p are the density, velocity, pressure, respectively. The viscous stress tensor τ is given by

τ = µ[∇u + (∇u) T - 2 3 (∇ • u)I], (2)
where µ is the fluid dynamic viscosity and I is identity matrix.

The total energy conservation can be expressed as

∂ρE ∂t + ∇ • [u(ρE + p)] = ∇ • (λ∇T ) + ∇ • (u • τ ), (3)
λ is heat conductivity. Total energy E = e + u 2 /2 is the sum of internal and kinetic energies.

Alternatively, the energy conservation (3) can be expressed in terms of entropy as ∂s ∂t

+ u • ∇s = 1 ρT ∇ • (λ∇T ) + 1 ρT τ : ∇u. (4)
Finally, a thermodynamic closure is needed, relating the internal energy e (or entropy s) with the pressure p and volume mass ρ through the temperature. Albeit not restricted to this equation of state, the perfect gas closure was adopted for this work, but the model is compatible of more complex closures. In the following,

p = ρ.r.T, e = C v .T + q, s = C v ln T ρ γ-1 + q ′ , (5)
where C v is the mass specific heat capacity at constant volume, assumed to be constant, r = R/W is the ratio of the perfect gas constant to the gas molecular weight, γ = C p /C v is the specific ratio, and q and q ′ are references which can be set to zero for a mono-constituent gas.

Hybrid recursive regularized thermal LB model

This section presents the model we propose to solve the governing equations [START_REF] Mohamad | Lattice Boltzmann method: fundamentals and engineering applications with computer codes[END_REF]. A hybrid recursive regularized thermal LB model relying on the third order Hermite expansion will be developed incorporating a finite volume approach for solving entropy equation.

Thermal LB model using the 3rd order Hermite expansion

The lattice Boltzmann equation (LBE) can be derived from the Boltzmann equation in the velocity space. It can be written as follows

∂f ∂t + ξ • ∂f ∂x = Ω f (6)
where f denotes the density distribution function, ξ is the particle velocity and Ω f is the collision operator. The macroscopic quantities, e.g., the density ρ, mean velocity u can be obtained from the distribution function as

ρ = f dξ, ρu = ξf dξ. (7)
A single-relaxation-time collision operator, Bhatnagar-Gross-Krook (BGK) [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF],

Ω f = - 1 τ (f -f eq ) (8)
is used as a basis to develop the proposed method, where τ is relaxation time. The equilibrium distribution function f eq can be derived as

f eq = ρ 1 2πrT D 2 exp (ξ -u) 2 2rT ( 9
)
where r is the gas constant, D is the number of spatial dimension and T represents temperature.

In most LB models based on nearest neighbor lattices (D1Q3, D2Q9, D3Q27, etc.), the density distribution function is approximated through a second-order Hermite polynomials expansion, leading to discrete equilibrium distribution functions f eq i as

f eq i = w i   ρ + c iα c 2 s ρu α + H (2) iαβ 2c 4 s ρu α u β   ( 10
)
where w i is the i th weight coefficient associated to discrete velocity c iα , the second order Hermite polynomial reads

H (2) iαβ = c iα c iβ -c 2 s δ αβ , (11)
with δ αβ being the classical Kronecker delta, and c s is the sound speed.

Using the second-order approximation [START_REF] Latt | Lattice Boltzmann method with regularized precollision distribution functions[END_REF] for the equilibrium function, the Boltzmann equation [START_REF] Lallemand | Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF] with BGK closure can be shown to be equivalent to the weakly compressible Navier-Stokes equation closed with the following athermal equation of state : p = ρc 2 s . Because of the assumption that pressure only depends on ρ, the athermal equation of state is clearly not suited for fully compressible flows, in which a source of heat can set a fluid in motion (or, equivalently, a disequilibrium in the density function).

To handle compressible flows through a non-isothermal equation of state (e.g. p = ρrT for ideal gas), it is then natural to consider the thermal Maxwell-Boltzmann distribution instead [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF], here expanded to the third order as

f eq i = w i   ρ + c iα c 2 s ρu α + H (2) iαβ 2c 4 s A (0) αβ + H (3) i 6c 6 s A (0) αβγ   . (12)
Indeed, the Hermite polynomials remain unaltered, the second-order given by [START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF] and the third-order being

H (3) iαβγ = c iα c iβ c iγ -c 2 s [c i δ] αβγ , (13)
where [c i δ] αβγ = c iα δ βγ + c iβ δ αγ + c iγ δ αβ . The second and third-order terms in the expansion now include full compressiblity as (p -ρc 2 s ), now allowed to non-zero values

A (0) αβ = ρu α u β + (p -ρc 2 s )δ αβ , (14)
A (0) αβγ = ρu α u β u γ + (p -ρc 2 s )[uδ] αβγ , (15)
where

[uδ] αβγ = u α δ βγ + u β δ αγ + u γ δ αβ .
The moments of the equilibrium distribution function [START_REF] Geier | The cumulant lattice Boltzmann equation in three dimensions: Theory and validation[END_REF] on nearest neighbor type lattices are then

i f eq i = ρ, ( 16a
) i c iα f eq i = ρu α , ( 16b
) i c iα c iβ f eq i = pδ αβ + ρu α u β , ( 16c
) i c iα c iβ c iγ f eq i = p[uδ] αβγ + ρu α u β u γ -Ψ αβγ . ( 16d
)
where the term Ψ αβγ is a correction due to defect of symmetry of the nearest neighbor type lattices for the third-order moment. This deviation is intended to be corrected by introducing a correcting term ψ i , which can be applied as a second order forcing term in lattice Boltzmann BGK equation following

∂f i ∂t + c iα ∂f i ∂x α = Ω f i + ψ i (17)
ψ i = -w i H iαβ 2c 4 s ∂ ∂x γ Ψ αβγ , i c iα c iβ ψ i = - ∂ ∂x γ Ψ αβγ (18)
Taking into account the Chapman-Enskog multiscale technique [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF] with the second order correcting term ψ i , the compressible Navier-Stokes equations can be recovered as follows (see Appendix A for the detailed Chapman-Enskog analysis)

∂ρ ∂t + ∂ ∂x α (ρu α ) = 0 (19) ∂ ∂t (ρu α ) + ∂ ∂x β (ρu α u β + pδ αβ ) = ∂ ∂x β µ( ∂u β ∂x α + ∂u α ∂x β - k ∂u γ ∂x γ ) ( 20
)
The correcting term ψ i is given, e.g., in two dimensional D2Q9 as

ψ i = w i 2c 4 s H ixx ∂ ∂x [ρu x (1 - p ρc 2 s -u 2 x )] + H iyy ∂ ∂y [ρu y (1 - p ρc 2 s -u 2 y )] , (21)
and in three dimensional D3Q27 lattice model as

ψ i = w i 2c 4 s H ixx ∂ ∂x [ρu x (1 - p ρc 2 s -u 2 x )] +H iyy ∂ ∂y [ρu y (1 - p ρc 2 s -u 2 y )] +H izz ∂ ∂z [ρu z (1 - p ρc 2 s -u 2 z )] , (22)
The correcting term ψ i in D3Q19 lattice is derived in Appendix B.

Hybrid recursive regularized collision model

Regularization procedures were proposed as to improve the numerical stability of BGK collision models for lattice Boltzmann methods [START_REF] Latt | Lattice Boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF].

They play the role of a filter on spurious ghost moments of the numerical scheme. Practically, a regularized distribution function is introduced through recomputing the non-equilibrium parts prior to the collision step.

By including a forcing term ψ i in the discrete velocity Boltzmann equation [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] and integrating it along the characteristic line with the trapezoidal rule, the LB equation can be expressed as:

f i (x α + c iα δ t , t + δ t ) -f i (x α , t) = δ t 2 Ω f i (x α + c iα δ t , t + δ t ) + Ω f i (x α , t) + δ t 2 ψ i (x α + c iα δ t , t + δ t ) + ψ i (x α , t) (23)
In order to remove this implicit treatment of the collision and forcing terms, a new distribution function is introduced:

f i = f i - δ t 2τ (f eq i -f i ) - δ t 2 ψ i (24)
The macroscopic density and momentum can be computed as

ρ = i f i + δ t 2 ψ i (25a) ρu α = i c iα f i + δ t 2 i c iα ψ i (25b)
Therefore, the LBGK equation can be expressed equivalently as

f i (x α + c iα δ t , t + δ t ) = f eq i (x α , t) + (1 - 1 τ )f neq i + δ t 2 ψ i (x α , t) (26)
where τ = τ /δ t + 1/2 and f neq i represents the off-equilibrium part of the new distribution functions,

f neq i = f i -f i eq + δ t 2 ψ i = 2τ 2τ -1 f neq i ( 27
)
The off-equilibrium distribution function can be also expanded using Hermite polynomials according to the Chapman-Enskog technique, leading to

f neq i = 2τ 2τ -1 f neq i ≈ w i H (2) iαβ 2c 4 s A (1)
αβ + H

(3) iαβγ 6c 6 s A (1) αβγ (28)
where

A (1) αβ = 2τ 2τ -1 i c iα c iβ f neq i = i c iα c iβ f neq i (29) = -δ t τ p[ ∂u β ∂x α + ∂u α ∂x β - k ∂u γ ∂x γ δ αβ ] A (1) αβγ ≈ u α A (1)
βγ + u β A (1) γα + u γ A

(1) αβ [START_REF] Safari | A lattice Boltzmann model for thermal compressible flows at low mach numbers beyond the Boussinesq approximation[END_REF] where k = γ -1, γ being heat specific ratio. The third order recursive regularization procedure of the thermal LB model by Eq. ( 30) is used to reduce the Galilean invariance and increase numerical stability, which could be derived by binomial theorem of Hermite polynomials [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF]. The third order off-equilibrium terms on gradient of temperature derived by Ref. [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF] are not considered in this work, since they are related to contributions of higher-order heat flux terms beyond Navier-Stokes equations.

Although the higher-order recursive regularization can reduce the defect of Galilean invariance, the non-hydrodynamic modes were also introduced by the higher-order off-equilibrium moment. A parameterized hybrid recursive regularized procedure is proposed to suppress those non-hydrodynamic modes. According to Eq. ( 30), the off-equilibrium moment A

αβ can be alternatively reconstructed by velocity gradients which could be calculated by a finite difference operator. Thus, following [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF] an explicit stabilization procedure relying on the linear combination of A (1,RR) αβ given by Eq. ( 29) and

A (1,FD) αβ
obtained by second-order finite-difference approximation of Eq. ( 30) is introduced, leading to the definition of the following hybrid recursive regularized collision operator

A (1,HRR) αβ = σA (1,RR) αβ + (1 -σ)A (1,FD) αβ (31)
where σ ∈ [0, 1] is an arbitrary weighting coefficient. A In order to preserve the isotropy of the scheme, a second-order isotropic central difference scheme is employed to compute the numerical gradient operator. The second-order isotropic central difference scheme is given in lattice units as

∂φ ∂x α = 1 c 2 s i w i c iα φ(x α + c iα ) (32)
Finally, the hybrid recursive regularization procedure for thermal compressible flow is summarized as follows:

f i (x α + c iα δ t , t + δ t ) = f eq i (x α , t) + (1 - 1 τ )R(f i neq ) + δ t 2 ψ i (x α , t) (33a) R(f i neq ) = w i H (2) iαβ 2c 4 s A (1,HRR) αβ + H (3) iαβγ 6c 6 s A (1,HRR) αβγ (33b) A (1,RR) αβ = i c iα c iβ f neq i (33c) f neq i = f i (x α , t) -f i eq (x α , t) + δ t 2 ψ i (x α , t) (33d)
In the simulations, physical and lattice units for length, time, mass and pressure are related through a reference length scale L 0 , the physical reference sound speed c sp0 = p 0 /ρ 0 for space and time, a physical reference pressure p 0 , a physical reference density ρ 0 and a lattice reference density ρ l0 . Accordingly, the physical variables can be converted from the quantities in lattice units (l subscript) as

δ x = L 0 /N n , δ t = δ x . c s c sp0 , t = N t .δ t x = N i .δ x , ρ = ρ l . ρ 0 ρ l0 , u = u l . δ x δ t , µ = µ l δ 2 x δ t , p = ρ 0 ρ l0 . δ 2 x δ 2 t .
The physical reference sound speed c sp0 could be set to the far-field or free-stream sound speed in compressible flows, it is more convenient to keep it as a free parameter. This allows to freely accelerate the convergence rate, by varying the time step in the same manner of Courant-Friedrichs-Lewy (CFL) number for conventional Navier-Stokes solvers. Taking ideal gas as an example, the physical sound speed for ideal gas can be computed in the LB model as follows

c sp = ∂p/∂ρ |s = γR g T = c sp0 γT /T 0 ( 34
)
where T 0 is physical reference temperature corresponding to p 0 .

Finite volume method for entropy equation

In the hybrid lattice Boltzmann method presented in this section, the entropy based energy conservation equation is solved separately using a finite volume technique. The entropy equation with viscous heat dissipation term is given as follow

∂s ∂t + u α ∂s ∂x α = 1 ρT ∂ ∂x α (λ ∂T ∂x α ) + Φ ρT . ( 35
)
where Φ represents the term of viscous dissipation. An explicit second-order Runge-Kutta scheme is adopted as temporal integration, which is given as

s (n+1/2) = s (n) + δ t 2 RHS[s (n) ] s (n+1) = s (n) + δ t RHS[s (n+1/2) ] (36)
The convective flux is constructed using MUSCL scheme, while the classical second-order central difference scheme is adopted for the diffusion term and term of viscous dissipation. The third order MUSCL scheme [START_REF] Kim | Methods for the accurate computations of hypersonic flows: I. AUSMPW+ scheme[END_REF] adopted in this study is given as follows:

u ∂φ ∂x = u i φ i+1/2 -φ i-1/2 δ x (37)
φ i+1/2 for instance, can be given as

φ i+ 1 2 =      φ L i+ 1 2 , u i > 0 φ R i+ 1 2 , u i ≤ 0 and φ L i+ 1 2 = φ i + ψ(r i ) 4 [(1 -κ)δφ i-1 2 + (1 + κ)δφ i+ 1 2 ], φ R i+ 1 2 = φ i+1 - ψ(r i+1 ) 4 [(1 -κ)δφ i+ 3 2 + (1 + κ)δφ i+ 1 2 ],
where κ = 1/3, and,

δφ i+ 1 2 = (φ i+1 -φ i ) , δφ i-1 2 = (φ i -φ i-1 ) , δφ i+ 3 2 = (φ i+2 -φ i+1 ) , δφ i-3 2 = (φ i-1 -φ i-2 ) , r i = φ i -φ i-1 φ i+1 -φ i ( 38
)
where i represents index of grid rather than lattice discrete velocity. The van Albada limiter function ϕ(r) = 2r/(1 + r 2 ) is used to avoid spurious oscillations [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF]. The coupling between the LB and finite volume (FV) solvers is illustrated in Fig. 1, along with the relevant data exchanges.

Start

Initialization ρ, u, p, T, s

Compute

f eq i (ρ, u, p) (12) ψ i (ρ, u, p) (21)
Streaming and collision step [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF] Update (ρ, u) (25)

1 st step of RK, compute s (1/2) Compute ρ (1/2) u (1/2) , T (1/2)
f i = f i (ρ, u, p) (T, s) (ρ, u) ( ρ , u ) ρ p = p ( ρ , T ) , N t + 1

Implementation of boundary conditions

The boundary conditions of density distribution function are implemented by a finite difference reconstruction approach along with the hybrid regularization procedure [START_REF] Latt | Lattice Boltzmann method with regularized precollision distribution functions[END_REF].

1. First, the macroscopic velocities ρ, u and p on the boundary nodes are computed, following boundary conditions commonly used in Navier-Stokes solvers. The entropy s and temperature T are then calculated by the thermodynamic closure and equation of state.

2. Next, the shear stress tensor is computed using the velocity gradients on boundary nodes, which are computed on these nodes using a first order biased finite-difference scheme.

3. Then, the density distribution function is reconstructed using σ = 0 in hybrid regularization in which the off-equilibrium moment is reconstructed as follows

A (1) αβ ≈ -τ δ t p ∂ u α ∂x β + ∂ u β ∂x α + k ∂ u γ ∂x γ δ αβ , (39)
A

αβγ = u α A (1) βγ + u β A (1) γα + u γ A (1) αβ , (1)
In the following, prescription of macroscopic variables (step 1 above), is done according to the boundary condition.

• For a velocity inlet, ρ, u, p are set as Dirichlet boundary conditions.

• For a pressure outlet, p is set as Dirichlet condition; ρ and u are extrapolated enforcing a zero-gradient condition.

• For a no-slip wall, u = 0 is set as Dirichlet condition, and ρ and p are extrapolated enforcing a zero-gradient condition.

Results and discussion

Since the base HRR LB model has already been assessed for simulation of thermal flows in low Mach number limit [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Feng | Regularized thermal lattice Boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Feng | A lattice-Boltzmann model for low-Mach reactive flows[END_REF], the present compressible HRR LB model is assessed considering three cases dealing with compressible subsonic to sonic shock-free flows:

1. isentropic vortex convection by a uniform flow, 2. acoustic pulse, 3. flow over a flat plate.

The first two configurations are treated as quasi-inviscid, with a non-dimensional viscosity µ = 10 -15 , as to test the model in extreme conditions. For the third configuration, variable viscosity is assumed.

For the D2Q9 lattice, the microscopic discrete velocity c iα is given by

c iα =              (0, 0) i = 0 (±1, 0), (0, ±1) i = 1 -4 (±1, ±1) i = 5 -8, (41)
the weights w i = 4/9, i = 0; w i = 1/9, i = 1, 2, 3, 4; w i = 1/36, i = 5, 6, 7, 8

and the lattice reference sound speed c s being 1/3.

Isentropic vortex convection by a uniform flow

This is one of the few exact analytical solutions of the compressible Euler equations. This test case is related to the convection of an isentropic vortex by an inviscid uniform flow at free-stream Mach number M a ≈ 0.42 (u ∞ = 0.5, γ = 1.4). It is used to show the ability of numerical schemes to accurately capture vortical flows. The free-stream conditions and perturbations added to the free-stream is initialized in a computational domain of size [0, 10] × [0, 10]. There is no entropy gradient in the flow-field. The free-stream flow is given by

ρ ∞ = 1, u ∞ = 0.5, v ∞ = 0, p ∞ = 1 (42)
and the density and velocity perturbations are initially given as a vortex, defined as

ρ = 1 - (γ -1) b 2 8γπ 2 e 1-r 2 1 γ-1 , p = ρ γ , (43)
u = u ∞ - b 2π e 1 2 (1-r 2 ) (y -y c ) , (44)
v = v ∞ + b 2π e 1 2 (1-r 2 ) (x -x c ) . ( 45
)
where b = 0.5 is the vortex strength and r = (x -x c ) 2 + (y -y c ) 2 1/2 is the distance from the vortex center (x c , y c ) = (5, 5).

A = 1 -min(ρ) 1 -ρ 0 ( 46
)
is reported in Fig. 3(d). The damping rate induced by numerical dissipation is quite small. Considering both the robustness and the accuracy, the HRR parameter σ = 0.5 appears to be a good trade-off for the simulation of inviscid compressible flows.

Simulations displayed in Fig. 2 were repeated varying the grid resolution [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] De Grazia | Direct numerical simulation of a compressible boundary-layer flow past an isolated three-dimensional hump in a high-speed subsonic regime[END_REF]100,200,400). The resulting profiles, reported in Figs. 4(a) and 4(b), illustrate well the grid convergence, as well as the method robustness. Further investigation of the L 2 norm of these profiles, illustrated in Fig. 4(c), shows that the convergence is second order in space, a classical characteristic of LB models [START_REF] Krüger | The lattice Boltzmann method[END_REF]. Stability at low grid resolutions is ensured by numerical viscosity, which we shall now assess. As an approximated method, we propose to compare the classical vorticity ω decay in incompressible fluid at rest [START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Tryggeson | Analytical vortex solutions to navier-stokes equation[END_REF][START_REF] Tryggeson | Stationary vortices attached to flat roofs[END_REF]

N x × N x , with N x =
∂ω ∂t = ∇ω Re , ( 47
)
where Re is the Reynolds number, leading to an exponential decay of the maximum vorticity in the vortex as e -t/Re . Continuity in stationary compressible flows being translated in a divergence free momentum ρu field rather than divergence free velocity field u, it is reasonable to approximate A (46) as

A(t) ≈ e -t/Re * , (48)
where Re * corresponds to a Reynolds number associated to a numerical viscosity ν * . The corresponding non-dimensional values, corresponding to a reference non-dimensional length of 3 (the approximate vortex diameter)

and the reference non-dimensional velocity of u ∞ = 0.5 are reported in both with increasing resolution and increasing σ. For further illustration, the fits obtained by using the parameters of Tab. 1 are reported in dashed lines in Fig. 3(d).

Two-dimensional acoustic pulse 4.2.1. Isothermal acoustic pulse

Next, the compressible LB model is assessed considering benchmark problems dealing with propagation of acoustic waves, more precisely the advection of an acoustic pulse by a uniform flow. A two dimensional Gaussian acoustic pulse with free-stream velocity U 0 = 0.0 and U 0 = 1.0, i.e.

at free-stream Mach number M a = 0 and M a = 1, is considered in this section. The initial acoustic Gaussian pulse, is given by [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF]

ρ = 1 + 0.01exp(-kr 2 ), p = 1 + 0.01exp(-kr 2 ), u = U 0 , v = 0 (49)
where r = (x -x c ) 2 + (y -y c ) 2 1/2 is the distance from the pulse center Periodic boundary conditions are implemented in all directions. In order to implement the simulation of isothermal acoustics to match reference analytical solutions, the entropy equation is not solved in this case and temperature is prescribed as a constant value.

The time step is δ t = δ x / √ T 0 with reference temperature T 0 = 25.0. An approximate inviscid condition is implemented setting the dynamic viscosity to µ = 10 -15 , which results in a relaxation time τ ≈ 0.5. Following the results of the previous test case, the HRR weighting parameter is set equal to σ = 0.5.

The pressure fields obtained by the LB method for the two Mach numbers at time t = 10 are shown in Fig. 5 and6. Furthermore, the fluctuating velocity and pressure profiles along the horizontal centerline at time t = 10 are compared to analytical solutions given in [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF] in Fig. 7. A very

Non-isothermal Gaussian pulse

Next, we consider a thermal acoustic wave traveling in radial direction including temperature evolution. The aim of this test case is to verify the correct propagation speed of sound waves and coupling between velocity, pressure and temperature. The initial condition is given by:

ρ = 1.0, p = 1 + exp(-kr 2 ), u = 0, v = 0 ( 50
)
where r = (x -x c ) 2 + (y -y c ) 2 1/2 is the distance from the pulse center (x c , y c ) = (0, 0). The perturbation parameter is taken equal to k = 40.

The computational domain size is [-2, 2] × [-2, 2], with periodic boundary conditions. The computation of the reference solution was performed using a second order TVD finite volume scheme with the Osher type flux on a radial grid consisting of 10 4 cells [START_REF] Dumbser | On universal Osher-type schemes for general nonlinear hyperbolic conservation laws[END_REF][START_REF] Tavelli | A pressure-based semi-implicit space-time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier-Stokes equations at all Mach numbers[END_REF]. analytically and numerically [START_REF] Chapman | Temperature and velocity profiles in the compressible laminar boundary layer with arbitrary distribution of surface temperature[END_REF][START_REF] De Grazia | Direct numerical simulation of a compressible boundary-layer flow past an isolated three-dimensional hump in a high-speed subsonic regime[END_REF]. Here, a computational domain with size -0.25 ≤ x ≤ 1, 0 ≤ y ≤ 0.25 is set up considering a flow at high speed subsonic u ∞ past a thin plate of length L(as illustrated in Fig. 9). The boundary condition(BC) is set as the follows:

• Symmetry BC on y = 0, 0.25 ≤ x < 0

• No-slip wall BC on y = 0, 0 ≤ x ≤ 1, adiabatic BC for M a = 0.3 while isothermal BC for M a = 0.87.

• Subsonic outflow on y = 0.25 with pressure p = 1/γ.

• Subsonic inflow on x = 0 with density ρ = 1, and velocity (u, v) = (u ∞ , 0).

• Subsonic outflow on x = 1 with pressure p = 1/γ. to match those of reference similarity solutions of compressible boundary layer [START_REF] White | Viscous fluid flow[END_REF].

In simulation at M a = 0.3, the computational grid is made of 1000 × 200 grid points. The spacing size is δ x = 6.25 × 10 -4 and the time step δ t is set to δ x /6 corresponding to a CFL number close to 0.13. At M a = 0.87, the computational grid has 2000 × 400 grid points. The spacing size is δ x = 3.125 × 10 -4 and the time step δ t is set to δ x /8 with a reference temperature T 0 equals 64. The HRR parameter is set equal to σ = 0.9 in all simulations. In Fig. 11, the analytical solution of temperature profile for compressible flow over flat plate is given by [START_REF] White | Viscous fluid flow[END_REF] T

T e = 1 + γ -1 2 M a 2 e (1 - u 2 U 2 e )
where T e is temperature at outer edge. The reference solution of density is calculated by equation of state under free stream pressure. In can be observed that the excellent agreement on both temperature and density between the results obtained by the present LB method and the reference solution is observed M a = 0.87. The coupling between fluid flow and heat transfer including viscous heat dissipation are well captured by the present LB model. By matching the scales of ǫ 1 , ǫ 2 we have

Conclusions

ǫ 1 : ( ∂ ∂t 1 + c iα ∂ ∂x 1α
)f eq i + f

( The deviation terms of Ψ xyy , Ψ yyy etc. can be computed in similarity. Thus the correcting term in D3Q19 lattice set is constructed as follows

ψ i = w i 2c 4

( 1 ,

 1 HRR) αβγ is recursively updated by re-computed off-equilibrium moment A (1,HRR) αβ .

Figure 1 :

 1 Figure 1: The algorithm proposed consists of a Lattice-Boltzmann (LB) solver coupled with a classical Finite Volume (FV) solver. Data exchanges between the two solvers are clearly identified.

Figure 2 :

 2 Figure 2: Profiles obtained by the LB method on isentropy vortex convection at t = 20T .

Figure 2

 2 Figure 2 displays the distributions of density and velocities along the domain symmetry lines after 20 flow-through-time (FTT), i.e. t = 20T . The density along the horizontal line of symmetry is plotted in Fig. 2(a) while x-velocity along vertical middle line is shown in Fig. 2(b). It can be observed that all results are in good agreement with the exact solution after 20 FTT. The numerical dissipation increases with decreasing HRR parameter σ where σ = 0.0 means that the off-equilibrium distribution function is completely reconstructed by the finite difference scheme. The typical density fields obtained using different HRR parameter σ values after 20 FTT are presented in Fig. 3. It is found that the isotropic evolution is well preserved in the simulations except when σ = 0.0. The history of the maximum of the amplitude A of the density disturbance normalized by its initial value

Figure 3 :Figure 4 :

 34 Figure 3: (a,b,c) Density contours after 20 periods, as obtained using different values for HRR parameter σ, on a 200 × 200 grid. (d) Evolution of A (46) at the center of the vortex, as obtained on a 200 × 200 grid: σ = 0 (•) σ = 0.5 (△), σ = 0.7 (▽), and on a 400 × 400 grid with σ = 0.5 (×. Dashed lines are the exponential fits of the form A(t) = e -t/Re * (48), as reported in Tab. 1.

(

  x c , y c ) = (0, 0). The perturbation parameter is given by k = ln2/b 2 , and b = 2 is the half-width of Gaussian function. A square computational domain with size [-20, 20] × [-20, 20] is used, along with a 400 × 400 grid.

= 1 Figure 5 : 1 Figure 6 :Figure 7 :

 15167 Figure 5: Contours of pressure obtained by the LB method on acoustic pulse advection at time t = 10.

Figure 8 :Figure 8

 88 Figure 8: Comparison of profiles along horizontal middle line.(line: present LBM, point: reference solution)

Figure 9 :Figure 10 :

 910 Figure 9: Sketch of the computational configuration. The darker sub-domain is used below for plotting results near the leading edge.

Figure 10 displays

 10 Figure 10 displays the contours of density and horizontal velocity obtained by the present LB model at M a = 0.3 and M a = 0.87 in the vicinity of the plate trailing edge. It is observed that the flow is well captured by the present LB method.

Figure 11 Figure 11 :

 1111 Figure 11 shows the profiles of density, temperature and horizontal velocity obtained by the present LB model. In the figure, η is the dimensionless coordinate with Illingworth transform which is defined as[START_REF] White | Viscous fluid flow[END_REF] η = U e √ 2x y

A

  thermal lattice Boltzmann model has been developed on standard lattices for the simulation of smooth subsonic to sonic compressible flows. It relies on an extension of the Hybrid Recursive Regularized LB model already proposed for low-Mach number flows. The model is easy to implement and parallelize: both the LB and FV stencils are limited to nearest neighbors. This also leads to a cpu cost partition close to 75-25 between the LB-FV solvers, indicating high potential regarding the cpu savings associated to the use of such model compared to classical methods. In simulation of an isentropic vortex convection, a two dimensional acoustic pulse and non-isothermal Gaussian pulse, the influence of the HRR weighting parameter σ has been intensively investigated. Numerical simulations demonstrate that the defect of Galilean invariance is effectively eliminated by σ > 0.5 in the inviscid regime. Especially, the non-hydrodynamic oscillation is suppressed by hybrid regularization using σ = 0.5 in simulation of non-isothermal Gaussian pulse. Excellent agreement has also been obtained for the simulation of a compressible laminar flow over flat plate at Mach number of 0.3 and 0.87, with results matching the Blasius solution.The low dissipation obtained with the presented model allows to envision acoustic applications in transonic flows, e.g. by accounting for the bulk viscosity.

18 i=0c(B. 3 ) 18 i=0c 18 i=0c

 1831818 iα c iβ c iγ f eq iThe following deviation terms of D3Q19 are obtained with θ = p/(ρc 2 s )Ψ xxx = (ρu x u x u x + 3pu x )ix c ix c ix f eq i = ρu x (θ -1 + u 2 x ), (B.4a) Ψ xxy = (ρu x u x u y + pu y ) -18 i=0 c ix c ix c iy f eq i = 1 6 ρu y (θ -1 + 3u 2 z ), (B.4b) Ψ xyz = (ρu x u y u z )ix c iy c iz f eq i = ρu x u y u z , (B.4c)

Table 1 .

 1 As expected, ν * , given here in non-dimensional unit, decreases

	Resolution σ	Re *	ν * × 10 5
	200 × 200	0	1998	75
	200 × 200 0.5 5665	26
	200 × 200 0.7 10910	14
	400 × 400 0.5 43726	3.4

Table 1 :

 1 Approximated numerical viscosity for the isentropic vortex convection.

  Considering that the equilibrium density distribution function satisfies the velocity moment condition, one can sum Eq.(A.3)and Eq.(A.4) in the velocity phase space. The t 1 and t 2 order of the continuity equation and momentum equation can be derived as αβγ + ρu α u β u γ -Ψ αβγ )] Considering equilibrium distribution function Eq. 12, the deviation terms can be calculated as follow Ψ αβγ = p[uδ] αβγ + ρu α u β u γ -

	Considering the correcting term ψ i given as
		ψ i = -	w i 2c 4 s	H iαβ	∂ ∂x γ	Ψ αβγ ,	i	c iα c iβ ψ i = -	∂ ∂x γ	Ψ αβγ	(A.11)
	and rewriting f	(1) i	1) in the Eq.(A.8) with Eq.(A.3), one obtains
		i	i τ (1) i + c iγ = ψ + (2) (0) i ∂ ∂x 1γ )f eq i -ψ f i τ = 0 (ρu α u β + pδ αβ ) + + c iα ∂ ∂x 1α )f c iα c iβ [( ∂ ∂t 1 ∂ i ∂t 1 ∂ = -τ ∂f eq i ∂t 2 + ( (1) ǫ 2 : c iα c iβ f i = -τ [ ∂t 1 ∂ ∂x γ Ψ αβγ	(0) i ]	(A.3) (A.4)
									+ (p[uδ] = -τ p[ ∂ ∂x 1γ ∂u β ∂x 1α +	∂u α ∂x 1β	-	k ∂u γ ∂x 1γ	δ αβ ]	(A.12)
	The following equations can be finally obtained:
											∂ρ ∂t 1 ∂ρ ∂t	+ +	∂ ∂x 1α ∂ ∂x α (ρu α ) = 0 (ρu α ) = 0	(A.5) (A.13)
	∂ ∂t	(ρu α ) +	∂ ∂x β	∂ ∂t 1 (ρu α u β + pδ αβ ) = (ρu α ) + ∂ ∂x 1β ∂x β (ρu α u β + pδ αβ ) = 0 ∂ [µ( ∂u β ∂x α + ∂u α ∂x β -	k ∂u γ ∂x γ	(A.6) δ αβ )] (A.14)
											∂ρ ∂t 2	= 0	(A.7)
								∂ ∂t 2	(ρu α ) +	∂ ∂x 1β	(	(1) i ) = 0	(A.8)
	Incorporating pressure equation in convective scale as follow
								∂p ∂t 1	+	∂ ∂x 1γ	(pu γ ) + kp	∂u γ ∂x 1γ	= 0	(A.9)
	where k = γ -1, γ is heat specific ratio but subscript of tensors. Combining
	Eq. (A.5) with Eq. (A.6), one can obtain
		∂ ∂t 1	(ρu α u β + pδ αβ ) =	∂ ∂x 1γ	ρu α u β u γ -u α	∂p ∂x 1β	-u β	∂p ∂x 1α
											-	∂ ∂x 1γ	(pu γ δ αβ ) -kp	∂u γ ∂x 1γ	δ αβ	(A.10)

i c iα c iβ f

  It is worth noting that only Ψ xxx , Ψ yyy , Ψ zzz is non-zero in D3Q27 lattice set, which exhibit a much simpler correcting term compared with in D3Q19 lattice set. Besides the correcting term of D2Q9 lattice set is the projection of D3Q27 lattice into two-dimensions.

	s	H ixx [	∂ ∂x	Ψ xxx +	∂ ∂y	Ψ xxy +	∂ ∂z	Ψ xxz ]
	+H iyy [	∂ ∂x	Ψ yyx +	∂ ∂y	Ψ yyy +	∂ ∂z	Ψ yyz ]
	+H izz [	∂ ∂x	Ψ zzx +	∂ ∂y	Ψ zzy +	∂ ∂z	Ψ zzz ]
	+2H ixy [	∂ ∂x	Ψ xyx +	∂ ∂y	Ψ xyy +	∂ ∂z	Ψ xyz ]
	+2H ixz [	∂ ∂x	Ψ xzx +	∂ ∂y	Ψ xzy +	∂ ∂z	Ψ xzz ]
	+2H iyz [	∂ ∂x	Ψ yzx +	∂ ∂y	Ψ yzy +	∂ ∂z	Ψ yzz ]	(B.5)

Acknowledgements

This work was supported by the French project CLIMB, with the financial support of BPIFrance (Grant P3543-24000), in the framework of the program "Investissement d'Avenir: Calcul Intensif et Simulation Numérique". This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program Investissements d'Avenir supervised by the Agence Nationale de la Recherche. This work was performed using HPC resources from GENCI-TGCC/CINES (Grant 2018-A0052A07679). This work was also supported by project of ProLB software (http://www.prolb-cfd.com).

Appendix A. Chapman-Enskog Analysis

To derive macroscopic equations, the density distribution function f i is expanded around the f eq i distributions as follows:

Appendix B. The correction term for D3Q19 lattice set

The D3Q19 lattice set, 3 dimensions with 19 velocities, is given by,

and the corresponding weights are


Arthur Charpentier
email: charpentier.arthur@uqam.ca

Emmanuel Flachaire

Extended Scale-Free Networks

, that sub-networks of scale-free networks are not scale-free. In the later, those sub-networks are asymptotically scale-free, but one should not forget about second-order deviation (possibly also third order actually). In this article, we introduce a concept of extended scale-free network, inspired by the extended Pareto distribution, that actually is maybe more realistic to describe real network than the strict scale free property. This property is consistent with Stumpf, Wiuf & May (2005) : sub-network of scale-free larger networks are not strictly scale-free, but extended scale-free.

1 From Strict Scale-Free to Scale-Free Types scale-free network is a network whose degree distribution follows a power law, at least asymptotically. When studying internet networks, [START_REF] Barabási | Emergence of scaling in random networks[END_REF] observed that some nodes, that they called hubs, had much more connections than others, and that the distribution of the number of links connecting to a node was a power-law. They coined the term scale-free network to describe that class of networks, when degrees have a power-law distribution. [START_REF] Clauset | Power-law distributions in empirical data[END_REF] studied real networks, and found some exhibiting that property. Nevertheless, recently, Broido & Clauset (2019) claimed that (strict) scale-free network are actually rare. More specifically, inspired by [START_REF] Alderson | Understanding Internet topology: principles, models, and validation[END_REF], they define various notions of weak or strong scale-free networks. If their taxonomy of scale-free network, in interesting we will consider here only the concept of strict scale-free if the degree distribution above a given cutoff k min is a power law (as in [START_REF] Barabási | Emergence of scaling in random networks[END_REF]).

As discussed in [START_REF] Voitalov | Scale-free Networks Well Done[END_REF], the scale-free property is closely related to the Pareto distribution, used in extreme value theory (and is the continuous version of the discrete power-law used for degrees). Nevertheless, as we will recall, this (strict) Pareto appears usually only above a very high threshold, and distributions are only Pareto type. Recently, [START_REF] Beirlant | Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions[END_REF] suggested to take into account the second order approximation: the first order has a power law, and so is the second order, with smaller tail index. This is the extended Pareto distribution. In this article, we will explore an "extended scale-free" property, and study its impact on networks. And as we will see, this property is close to one described in [START_REF] Stumpf | Subnets of scale-free networks are notscale-free: sampling properties of networks[END_REF] when studying sub-networks of larger ones.

Strict Scale-Free

Consider a network (V , E ), and let n denote the number of nodes. Following [START_REF] Barabási | Emergence of scaling in random networks[END_REF], the network is said to be scale-free if

p k min (k) = P[D = k] = Ck -α , α > 1, for k ≥ k min ≥ 1
where α is the scaling exponent, and C is the normalization constant. Inference is performed using a degree sequence

d 1 , • • • , d n .
Equivalently, the log-log plot should be linear (and the absolute value of the slope is the scaling exponent)

log p k min (k) = log c -α log k.
This linear property for the logarithm of the frequency is the one usually used in network studies. Note that p 1 is also called Zipf's distribution.

Definition 1. The discrete Zipf's distribution is defined as

p d-Z,α (k) = 1 ζ (α) k -α for k ∈ N + = {1, 2, . . .} where ζ () is Rieman's function, ζ (α) = ∞ ∑ k=1 k -α , defined for α > 1.
In this article, we exclude nodes that have 0 connection. In many applications, some nodes are just disconnected from the others, having a null degree. In social networks, those individuals might be simple observers, and do not interact with others. This will also appear when considering sub-networks. Hence, when dealing with sub-networks in the last section of this article, we will remove nodes that are not connected to anyone.

Discrete PD with Cumulative Probabilities

Another popular approach is to consider cumulative probabilities, instead of frequencies. An interesting feature is that the cumulative probability of a power law probability distribution is also power law, but with an exponent α -1. More specifically, let F(k) = P[D > k] = 1 -P[D ≤ k] denote the complementary cumulative distribution function. If we consider a continuous version of p(x) = Cx -α , we obtain

F(x) = ∞ x p(t)dt = C ∞ x t -α dt = Cx α-1 α -1 = γx α-1
which is also a power function. And it is actually possible to derive a discrete probability function from a power F function, or more such as the (standard) Pareto function (see [START_REF] Arnold | Residual life time at great age[END_REF] or chapter 20 in [START_REF] Johnson | Continuous univariate distributions[END_REF])

F PD,ξ (x) = x -1/ξ for x ≥ 1, with ξ > 0.
Definition 2. The discrete strict Pareto distribution is defined as

p d-PD,ξ (k) = k -1/ξ -(k + 1) -1/ξ for k ∈ N + defined for ξ ∈ R + .
For discretized version of continuous ones, we will use tail index ξ , having in mind the fact that α is of order 1 + 1/ξ . The popular case α ∈ (2, 3) means that ξ ∈ (1/2, 1), with α all the more small that ξ is large. Observe that if the degrees have a discrete strict Pareto distribution, their expected value is

< D d-P >= ∞ ∑ k=1 kp d-PD,ξ (k) = ∞ ∑ k=1 k -1/ξ = ζ (1/ξ )
which is different from (1ξ ) -1 obtained with a continuous Pareto distribution.

Discrete GPD and Second Law of Extremes

Pareto distributions are very popular since the second law of extremes (see [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] and Balkema & de Haan (1974)) which states that if X is a random variable such that there exists a function a(u) such that

a(u) -1 (X -u)|X > u → Z, as u → ∞
(in the weak convergence sense) for some non-degenerate Z on R + , then Z follows a Generalized Pareto (GPD) with complementary cumulative distribution function

F GPD,σ ,ξ (x) = P[Z > x] = 1 + ξ x σ -1/ξ , for x ≥ 0.
As a consequence, for large u, we have usually the approximation for

P[X -u > x|X > u] P[a(u) a(u) -1 (X -u) > x|X > u] ∼ P[a(u)Z > x] = P[Z > a(u) -1 x] = F GPD,σ a(u),ξ (x)
as suggested in [START_REF] Davison | Models for exceedances over high thresholds[END_REF]. We will then write X ∈ MDA ξ -for Max-Domain of Attraction. Nevertheless, as proved in [START_REF] Anderson | Local limit theorems for the maxima of discrete random variables[END_REF] and [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF], this approximation might not be valid if X as a discrete support. An important additional property if to have a long-tailed distribution for F in the sens that

F(u + 1) F(u) → 1 as u → ∞
As proved in Shimura ( 2012), a discrete random variable X ∈ MDA ξ if and only if (i) p is long-tailed and (ii) X = X where X ∈ MDA ξ . And in that case,

p u (u + k) = P[X -u = k|X > u] ≈ p d-GPD,a(u)σ ,ξ (k), for k ∈ N
Thus, following [START_REF] Krishna | Discrete Burr and discrete Pareto distributions[END_REF], define :

2/13
p d-GPD,σ ,ξ (k) = F GPD,σ ,ξ (k -1) -F GPD,σ ,ξ (k) = 1 + ξ k -1 σ -1/ξ -1 + ξ k σ -1/ξ for k ∈ N + defined for ξ ∈ R + .
A (continuous) GPD random variable can be expressed as a mixture of exponential variables, that is an exponential random variable with a Gamma distributed rate parameter : if X ∼ E (Λ) with Λ ∼ G (α, β ), then X has a GPD distribution, with tail index ξ = 1/α. Interestingly, as proved in [START_REF] Buddana | Discrete Pareto Distributions[END_REF], a similar property holds for the discrete d-GPD, which is a mixture of geometric variables (with also Gamma heterogeneity).

The three probability distributions (p d-Z , p d-P and p d-GPD ), for similar tail exponent ξ , can be visualized on Figure 1 (for the Zipf, α = 1 + 1/ξ ).

Regular Variation and Power-Law Type

This power law property is deeply related to the concept of scale-free distribution : scale-free means that the distribution is the same whatever scale we consider. Hence, for any λ , f (λ k) is proportional to f (k), or f (λ k) = h λ f (x) (we must consider here the continuous version f and not p since, unless λ ∈ N, λ k is not always an integer). Thus, since

h λ = f (λ )/ f (1), we get that f (λ x) = f (λ ) • f (x)
, which is the multiplicative version of Cauchy's functional equation (also called Hamel-Cauchy), with unique solution f (x) = x -α (up to a multiplicative constant). Hence, scale-free means that

f (λ x)/ f (x) is constant (and in that case, f (x) = Cx -α ). A natural extension is to assume that f (λ x)/ f (x) is asymptotically constant.
A function g is said to be regularly varying (at infinity) if g(tx)/g(x) tends to t θ , for some θ ∈ R, when t goes to infinity. If θ = 0, then g is said to be slowly varying, to derive an extended version of the power law. Definition 4. A continuous variable X is said to be Pareto-type distributed, with tail exponent ξ if P[X > x] = x -1/ξ (x) for some slowly varying function .

In section 2, the idea will be to consider a simple parametric expresion for function , that will decay to a constant at some power speed.

Probability-Generating Function of Scale-Free Distribution

An alternative way to describe the distribution is not to use p, but its probability-generating function (PGF), G(s), defined as

G(s) = ∞ ∑ k=0 p(k)s k , 3/13
for instance, with a Poisson variable, G(z) = e λ (z-1) , while with a power law, or a scale free distribution

G(s) = 1 ζ (α) ∞ ∑ i=1 i -α s i = Li α (s) ζ (α)
where Li α is Jonquière's polylogarithm function (see [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]). We will use that alternative representation when focusing on sub-networks.

Extended Scale-Free

The Extended Pareto Distribution (EPD) was introduced in Beirlant, Joossens & Segers (2009), and there are many way to derive that distribution, most of them being equivalent.

2.1 Mixture of Scale-Free [START_REF] Hall | On some simple estimates of an exponent of regular variation[END_REF] suggested to write a Pareto type distribution

F(x) = x -1/ξ (x) as F(x) = x -1/ξ • C(1 -δ (x)).
Here, is not only slowly varying, but also (x) tends to C when x goes to infinity (at some power rate). More specifically, assume that δ (x) = Dx β + o(x β ) where β < 0. If δ (x) = Dx β , we can write

F(x) = C 1 x -γ 1 +C 2 x -γ 2
which can be seen as a mixture of two (strict) Pareto distributions.

Second-Order Regular Variation

The first law of extremes (also called Fisher-Tippett theorem) is based on the limiting distribution of maximum x (n) of an i.i.d. sample {x 1 , • • • , x n }. More precisely, assume that there exists a function a(n) such that

a(n) -1 X (n) -F -1 (1 -1/n) → Z, as n → ∞
for some non-degenerate Z on R + , then Z follows a Generalized Extreme Value (GEV) distribution (see [START_REF] Embrechts | Modelling Extremal Events[END_REF] or [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]). Let U denote the quantile function U(x) = F -1 (1 -1/x), then somehow, we might be interested by the limit (if it exists) of a(n) -1 (U(xn) -U(n)) when n goes to infinity. This is related to the concept of extended regular variation (see de [START_REF] De Haan | Extreme Value Theory[END_REF]) : g is said to be ERV γ if lim t→∞ g(tx)g(t) a(t) = c x γ -1 γ which can be seen as extension of regular variation of index γ. For instance, the quantile function of a (strict) Pareto distribution with index ξ , U(x) = x ξ , and with auxiliary function a(t) = ξ t ξ , a(n) -1 (U(xn) -U(n)) = (x ξ -1)/ξ (see [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]). Second order regular variation is obtained assuming that there is a function b such that

lim t→∞ 1 b(t) g(tx) -g(t) a(t) - x γ -1
γ exists, and is denoted h(x). de Haan & Stadtmüller (1996) obtained a general expression for h, related to some index ρ. In a nutshell, following Drees (1998) and [START_REF] Cheng | The edgeworth expansion for distributions of extreme values[END_REF], the limit can be expressed

lim t→∞ 1 b(t) g(tx) -g(t) a(t) - x γ -1 γ = x γ+ρ -1 γ + ρ
with ρ < 0 (theorem B.3.10 in Albrecher, Beirlant & Teugels (2017)).

Extended Scale-Free

For the strict Pareto distribution, we have seen that But let us consider the following extension (based on the expression of the second-order regular variation)

lim t→∞ x -γ F(tx) F(t) = 1 4/13
lim t→∞ x -γ F(tx) F(t) = 1 + x ρ -1 ρ , for some ρ ≤ 0
or, up to some affine transformation,

F(x) = cx -γ [1 + x ρ (x)]. Since (1 + u) b ∼ (1 + bu), define (changing γ in α, ρ in τ) F(x) = P[X > x] = [x (1 + δ -δ x τ )] -1/ξ , for x ≥ u
where τ ≤ 0 and δ > max(-1, 1/τ). This is the Extended Pareto Distribution, as define in [START_REF] Beirlant | Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions[END_REF].

Definition 5. The discrete extended Pareto distribution is defined as

p d-EPD,δ ,τ,ξ (k) = F EPD,δ ,τ,ξ (k -1) -F EPD,δ ,τ,ξ (k) = (k -1) 1 + δ -δ (k -1) τ -1/ξ -[k (1 + δ -δ k τ )] -1/ξ for k ∈ N +

Shifted Pareto Distributions

So far, we defined (discrete) distributions for degrees taking values in {1, 2, 3, . . .}. Quite naturally, one can that D has a Pareto distribution with a shift of u ∈ N + if Du has a Pareto distribution. For instance, with a strict Pareto distribution, when plotting the complementary cumulative probability function F d-PD on a log-log scale, the function is a (semi)-straight line with slope -1/ξ , starting in (u, 2 -1/ξ ).

5/13

, • • • , x n }, sorted, such that x (1) ≤ x (2) ≤ • • • ≤ x (n) , ξ = 1 n n ∑ i=1 log x (i) -log x (1)
(see Appendix 6.1 for a brief justification) but one can also focus on the k largest values

ξ k = 1 k n ∑ i=n-k+1 log x (i) -log x (n-k)
This estimator is strongly consistent, ξ k a.s.

→ ξ and (with further assumptions, see [START_REF] Embrechts | Modelling Extremal Events[END_REF])

√ k ξ k -α L → N (0, ξ 2 )
Nevertheless, this estimator performs badly when the sample is not strictly Pareto distributed, see Figure 4. For the EPD, Albrecher, [START_REF] Albrecher | Reinsurance: Actuarial and Statistical Aspects[END_REF] suggest to use maximum likelihood techniques.

Inference for Discrete Pareto Distributions

Two techniques are used to estimate parameters (whatever the underlying scale-free distribution considered). The first one is based on the chi-square statistic,

Q(θ ) = k max ∑ k=1 (n • p d-,θ (k) -n k ) 2 n • p d-,θ (k)
where n k is the number of nodes with exactly k neighbors. Actually, to get a more robust version, if n k is too small, we will regroup per classes, to have (at least) 10 nodes per class (see Appendix 6.2). An alternative is to use maximum likelihood techniques (see Appendix 6.2).

Strict and Extended Scale-Free Networks

Before studying real networks, let us generate networks that are extended scale-free, to see what they look like.

6/13

Figure 5. Average shortest path on simulated networks with n = 1, 000 notes, when degrees have an extended scale free distribution EPD(γ = 1.15, κ, τ) when -τ varies from 0 (strict Pareto) to 5. The shaded area is the 90% confidence band obtained with 500 simulated networks (for each τ).

Generating a Network from Degree Distribution

Consider a sequence

d 1 , • • • , d n of non-negative integers such that d n ≤ • • • ≤ d 1 .
From Erdös-Gallai theorem, see [START_REF] Tripathi | A short constructive proof of the Erdös-Gallai characterization of graphic lists[END_REF], that sequence can be represented as the degree sequence of a finite simple graph on n vertices if and only if

d 1 + • • • + d n is even, and k ∑ i=1 d i ≤ k(k -1) + n ∑ i=k+1 min(d i , k)
holds for every 1 ≤ k ≤ n. In this section, we use the methodology described in [START_REF] Newman | Random graphs as models of networks[END_REF] to generate graphs with an Extended Pareto distribution for the degree 1

Network Structures

On Figure 5, we can see the average shortest path for all nodes in the largest connected subgraph. This was obtained by averaging 1,000 simulated networks with n = 1, 000 nodes, with various τ. The larger the absolute value of τ, the longer the average shortest-path distance.

Heuristically, this can be explained since with a strict power law distribution, sub-graphs are connected to each other through (big) hubs, and those network have a small-world property : everyone is close to anyone. With a (strong) second order, there are less very big hubs, and more smaller one : the distance w to anyone tends to be, on average, longer. [START_REF] Stumpf | Subnets of scale-free networks are notscale-free: sampling properties of networks[END_REF] claimed that sub-networks of scale-free networks are not scale-free anymore. Of course, this result depend on how we sub-sample from a general network, and how scale-free is defined. Let {V , E } a network, i.e. a collection of vertices and edges. Let A denote its adjacency matrix : let i and j denote two nodes in V , then A i, j = 1 if an only if (ij) is in E . We assume here that there are no zero-degree node, i.e. ∀i ∈ V , ∃ j ∈ V such that A i, j = 1.

Sub-network of Scale-Free Networks

To generate a sub-network, select randomly (and uniformly) a sub-sample of nodes V * , then extract the sub-adjacency matrix A , and the (ij) is in E if an only if A i, j = 1. Interestingly, one can easily write the PGF of the degree distribution on the sub-network,

G * (s) = G(1 -p(1 -s))
where p is the probability to keep a given node. Since we excluded orphaned nodes, it is necessary to rescale the PGM, and then

G * (s) = G(1 -p(1 -s)) 1 -G(1 -p)
1 Implemented in the graph library, in R, see [START_REF] Gentleman | graph: graph: A package to handle graph data structures[END_REF]. As discussed in [START_REF] Stumpf | Subnets of scale-free networks are notscale-free: sampling properties of networks[END_REF], if we get back to simple scale-free network, we can obtain the following

7/13

P[d = k] ∼ 1 k , then P[d > k] ∼ 1 k 2 while P * [d * > k] ∼ 1 k(k -1) P[d = k] ∼ 1 k 2 , then P[d > k] ∼ 1 k 3 while P * [d * > k] ∼ 1 k(k -1)(k -2)
Observe that those two reminds us of the Hall class (see [START_REF] Hall | On some simple estimates of an exponent of regular variation[END_REF]). Hence

1 k(k -1) = [k 2 (1 -k -1 )] -1
which is an extended Pareto distribution with index τ = -1.

Real Internet Networks

In order to illustrate this second order property, we will use data from the Stanford Network Analysis Project (SNAP), from Facebook2 , Twitter3 and Google Plus4 (see [START_REF] Mcauley | Learning to Discover Social Circles in Ego Networks[END_REF]). The first one contains 4,039 nodes and 88,234 edges, the second one contains 81,306 nodes and 1,768,149 edges and the third 107,614 nodes and 13,673,453 edges.

In the SNAP-Facebook dataset, we have we have 10 sub-networks. It is known for being a not scale-free network, see ?. This is confirmed on the left of Figure 8 where no extended distribution can be used to capture such a strong concavity.

For the Twitter network, Aparicio, Villazón-Terrazas & Álvarez (2015) fitted a scale free distribution P(k) ∼ Ck -λ , outgoing degree distribution has tail index of λ = 2.1715 while incoming degree distribution has tail index λ = 1.8778. Here, we did not distinguish incoming and outgoing edges. When fitting an extended Pareto distribution, we obtained ξ = 0.757 (or 1 + ξ -1 = 2.32, consistent with the values obtained in [START_REF] Aparicio | A Model for Scale-Free Networks: Application to Twitter[END_REF]) and τ = -1. For Google+, we obtained also τ = -1. This value is consistent with the result obtained in [START_REF] Stumpf | Subnets of scale-free networks are notscale-free: sampling properties of networks[END_REF], and can be visualized on the center of Figure 8 (for Twitter network) and on the right of 8 (for Google+ network). On those two sets of figures, the parametric fitted distribution is added to the scatterplot.

Figure 1 .

 1 Figure 1. p d-Z , p d-P and p d-GPD , on a log-log scale on the left, and for the first 10 values on the right, with the same index ξ .

Figure 2 .

 2 Figure 2. p d-Z and p d-mix-P on a log-log scale on the left, and p d-Z and p d-EP on the right.

Figure 3 .

 3 Figure 3. F d-Z and F d-mix-P on a log-log scale on the left, and F d-Z and F d-EP on the right.

Figure 4 .

 4 Figure 4. Hill plot, for a (strict) scale free distribution on the left, and an extended scale-free distribution on the right.

Figure 6 .Figure 7 . 13 Figure 8 .

 67138 Figure 6. Strict Scale Free Network

http://snap.stanford.edu/data/ego-Facebook.html

http://snap.stanford.edu/data/ego-Twitter.html

http://snap.stanford.edu/data/ego-Gplus.html 9/13

Appendix

Fitting Continuous Distributions

Consider a continuous power law distribution, with density f (x) = (α -1)x -α , x ≥ 1.

The likelihood of a sample x = {x 1 , . . . , x n } is then

For convenience, use the logarithm of the likelihood,

The maximum of the logarithm of the likelihood function is obtained when

Fitting Discrete Distributions

It was mentioned in section 3.2 that the chi-square distance can be used to estimate the (unknown) parameter

A more robust version can be obtained by regrouping too-small degrees, to have at least 10 nodes : consider some (consecutive) partition K 1 , . . . , K m such that ∑ k∈K j n k ≥ 10 for all j, and then

where

Then the estimator is θ = argmin{Q(θ )}

For the discrete EPD model, given a vector x of degrees, the R code to compute the chi-square distance between the empirical distribution and p d-EPD is, for some given θ (i.e. values gamma, tau and kappa) Then, use optim to find the maximum of the log-likelihood, or the minimum of the chi-square distance. On Figure 9, we can visualize the boxplots of the six-estimators of α considered here, on 1,000 simulated samples, with two techniques and three underlying distribution (a discrete strict Pareto with tail index al pha = 1.15). On the left, we use the chi-square minimum distance, and on the right, the maximum-likelihood technique. We consider either a strict Pareto, a Generalized Pareto (GPD) and the Extended Pareto (EPD).