domain: spi phys

text:
A Sicard

D Boscher
email: daniel.boscher@onera.fr

D Lazaro

S Bourdarie

D Standarovski
email: denis.standarovski@cnes.fr

Keywords:

A new specification model of low energy electrons fluxes, included in GREEN model, has been developed at ONERA. This model is based on several data sets, from low altitudes with NOAA-POES measurements to higher altitudes with POLAR, THEMIS and LANL measurements. It provides mean electron fluxes versus L, B/B eq , MLT and energy (for energies between 0.19 keV and few 10s of keV). In this paper, the model is compared to AE9/SPM model and Van Allen Probe measurements.

Index Terms-plasma

such a model. Thus, a new empirical model of low-energy plasma, based on THEMIS, POLAR, NOAA-POES and LANL data, has been developed. As described above, each data set has been analyzed in several dimensions: L between 3 and 10 with a resolution of 0.2, MLT with a resolution of 1 hour and energy between 0.19 and 30 keV with five energy channels 0.19 keV, 0.84 keV, 2.6 keV, 8 keV and 30 keV. These data are analyzed, cross-calibrated, put together and parametric fits are calculated to provide pitch angle distribution. In this paper the model thus obtained for electrons from energy >0.1 keV to energy >30 keV is presented and compared to AE9/AP9/SPM model and validated with Van Allen Probe data. Note that this low energy electron model is not solar cycle dependent. This low-energy electron model is a part of the GREEN (Global Radiation Earth ENvironement) model [START_REF] Sicard | GREEN: the new Global Radiation Earth ENvironment model (beta version)[END_REF].

II. DATA USED

A. Low Earth Orbit Data

At LEO orbit, NOAA-POES spacecraft (POES-15, POES-16, POES-17, POES-18 and POES-19) measurements have been used for two reasons. First, thanks to the POES spacecraft fleet, long term data are available since 1998 with the same detectors TED and MEPED, which are part of SEM2 suite of instrument [START_REF] Evans | Polar Orbiting Environmental Satellite Space Environment Monitor -2 Instrument Descriptions and Archive Data Documentation[END_REF]. From the POES/SEM2 database, four electron channels from TED instrument, 0.19 keV, 0.84 keV, 2.6 keV and 8 keV, and one electron channel from MEPED instrument, >30 keV, were considered. Fig. 1 shows the time evolution of 0.19 keV electron flux seen by five POES spacecraft for 4.6<L<4.8 and 10h<MLT<11h. Some doubt exist on the TED electron channels efficiencies that degrade over time. In this case, the TED flux intensities would slowly decrease over time with periodic increases when the electron channels voltages are adjusted to compensate for the degradation. Fig. 1 indicates that such a behavior is not clearly seen in the data.

Fig. 1 : Evolution of 0.19 keV electron flux [cm -2 .s -1 .sr -1 .MeV -1 ] versus time measured by five POES spacecraft for 4.6<L<4.8 and 10h<MLT<11h. Moreover during overlaps from one spacecraft to the next no clear discrepancies between data sets are found whereas one would expect them to be larger when one instrument has spent years in space while the other is brand new. In this study the impact of such degradation will be assumed to be negligible. Electron data have been analyzed and filtered in order to remove bad data. Fig. 2 shows an example of data filtering applied to POES-15/TED 0.84 keV electron channel. In this example one can easily observe a correlation between the 0.84 keV electron fluxes (on the bottom panel) and the 1 MeV electron fluxes measured by POES-15/MEPPED (on the top panel) during the 2004 magnetic storms. As electrons with very low energy do not have the same behavior than high energy electrons, it turns out the 0.84 keV electron channel is contaminated by high energy electron. Cleaning the data is thus essential to obtain a sanitized data set. To clean the data, some correlations have been done between POES/TED measurements and high energy electrons measured by POES/SEPEM as plotted on Fig. 3.

B. Data in the magnetosphere

In order to obtain a global model at low energy up to L=10, data at other orbits than LEO have been taken into account. In this study, data from four other fleets of spacecraft have been used: THEMIS, LANL-GEO spacecraft, Van Allen Probes and POLAR. While the first three fleets have orbits near the equator, POLAR has an elliptical orbit and covers a wider part of the magnetosphere in term of L and B/B Eq . In this study, THEMIS and POLAR data are used to develop the model due to their wide coverage in L and MLT. As it is described in the next part, LANL data are used to calibrate all the data sets between each other. The validation of the model is done thanks to Van Allen Probes data. THEMIS is a fleet of five spacecraft launched in 2007. Only THEMIS-A, THEMIS-D and THEMIS-E have a time coverage large enough to be used in this study. Data from ESA instrument [START_REF] Mcfadden | The THEMIS ESA plasma instrument and in-flight calibration[END_REF], covering energies between 7 eV and 26 keV, and data from SST instrument covering 26 keV to 200 keV were investigated. As for POES data contamination by high energy electrons has been identified. The best way to clean ESA data by removing contaminated data would be to correlate ESA data with high energy measurements on the same spacecraft, that is to say on THEMIS/SST. However the higher energy measured on SST is 700 keV, which is not high enough to have a good correlation with ESA data as it can be observed on Fig. 5. Moreover Fig. 5 shows also that THEMIS/SST data are saturated. That's why we choose to correlate THEMIS/ESA data with high energy electron measurements on POES-15. This correlation is not a time to time correlation due to the difference of orbit between THEMIS and POES-15 but it is the best way we found to clean this data set. All THEMIS/ESA channels have been cleaned in the same way. As for POES /TED data, THEMIS/ESA data are contaminated by high energy proton below L<2.2. These data have also been filtered. The white color near L=4 on Fig. 5 shows that many measurements are filtered due to contamination. In this study the long-term mean flux is calculated by averaging all the non contaminated measurements. Thus, a question arises: how this cleaning procedure affects the long-term mean flux? To answer this question, we have calculated the averaged flux on a calm period without contamination (in 2009) at L=4 to obtain a reference value. Then, we have calculated the longterm (2008-2016) mean flux at L=4 by averaging all the measurements, the contaminated ones being set artificially to the averaged flux in calm period in 2009. At L=4, over the period 2008-2016, the ratio between mean flux without contaminated measurements and the mean flux with the contaminated measurements set to averaged flux during calm period is a factor 1.5. But this does not mean that the model developed here will overestimate flux by a factor 1.5 at L=4 because the real mean flux would probably be between the both values. We decided to use the higher value to be conservative. Then electron data from 2008 to 2015 have been sorted according to MLT and L as represented on Fig. 6 at two energies: 0.19 keV (on the left) and 8 keV (on the right). On the 0.19 keV map, the magnetosheath at L>8 on the day side can be clearly observed. As for POES data, THEMIS data are contaminated at L<2.2 by energetic protons. POLAR was launched in 1996 on an elliptical orbit with several instruments on-board including HYDRA [START_REF] Scudder | Hydra -A 3 dimensional electron and ion hot plasma instrument for the POLAR spacecraft of the GGS mission[END_REF], a detector measuring electrons and protons from 12 eV to 18 keV. Due to its orbit and in contrast to the two previous spacecraft fleets, POLAR data cover not only the equator as THEMIS, or only the highest latitudes as POES but also intermediate regions. POLAR data have been analyzed in detail as well and some periods when measurements are suspicious have been removed. Contamination by high energy electrons has been filtered out following the same strategy as for POES spacecraft using high energy electron data on POLAR/CEPPAD (Fig. 7). Then POLAR data have been sorted in three classes according to the equatorial pitch angle α eq (α eq =asin(√(B eq /B), with B eq the equatorial magnetic field and B the local magnetic field) deduced from the spacecraft position along the field line: α eq <30° (high latitude), 30°<α eq <60° and α eq >60° (near equator). Moreover, comparing the cartography from POLAR (on the right) to the one from THEMIS (Fig. 6) in the equatorial region, 8 keV electron fluxes from POLAR are a little bit lower than those from THEMIS. A data cross-calibration is necessary and the procedure will be described in the next section of this paper. Los Alamos National Laboratory (LANL) has a fleet of many spacecraft at geostationary orbit since 1976 with several particle detectors on-board including MPA instrument (Magnetospheric Plasma Analyser) [START_REF] Mccomas | Magnetospheric plasma analyzer: Initial three-spacecraft observations from geosynchronous orbit[END_REF] measuring electrons from 100 eV to 34 keV. Taking into account the very good quality of these data but the narrow coverage in L, they have only been used to cross-calibrate the data sets together (see section III-A). Note that data from all LANL spacecraft (from 1989-046 to LANL_02A) have been investigated. Van Allen Probes A and B have been launched in 2012 on an equatorial orbit between 500 km and 30600 km. The on-board plasma instrument HOPE [START_REF] Funsten | Helium, Oxygen, Proton, and Electron (HOPE) Mass Spectrometer for the Radiation Belt Storm Probes mission[END_REF] is measuring electrons from 15 eV to 50 keV. Note that the release 3 of the data has been used. It was decided to validate the model against Van Allen probes HOPE data. Table 2 resumes the characteristics of the entire electron data sets used in this study at other orbits than LEO.

III. MODEL DEVELOPMENT

A. Data sets cross-calibration

In order to have coherent data sets to develop an empirical specification model, their cross calibrations are essential. First, all the data sets have been interpolated in energy in order to obtain the same five energy channels as POES data. (0.19, 0.84, 2.6, 8 keV and >30 keV). Then, the standard way to cross-calibrate data independent sensors is through the comparison of fluxes during magnetic conjunction, that is to say when both sensors are "close" in terms of L, MLT, time and during quiet magnetic activity. To do this we need a reference data set. In our case we choose to use LANL/MPA data as a reference [START_REF] Mcfadden | The THEMIS ESA plasma instrument and in-flight calibration[END_REF]. This data set has been studied and used many times to develop empirical models [START_REF] Ginet | AE9, AP9 and SPM: New models for specifying the trapped energetic particles and space plasma environment[END_REF][2] [START_REF] Sicard-Piet | A new international geostationary electron model: IGE-2006, from 1 keV to 5.2 MeV[END_REF]. However, it is important to discuss about how the LANL-GEO frame charging affects the fluxes. These data have already been studied in terms of spacecraft potential effects and it is possible to correct electron flux according to this potential and the Liouville's theorem [10][11]. Just to keep in mind, the long-term mean electron flux at 0.19keV at GEO orbit is 1.5 higher when the correction from the spacecraft potential is applied than without the correction. For higher energy the correction has no influence on the mean flux. In order to be sure that the data used in this study are as close as possible to reality, the correction has been applied. However, LANL/MPA data are only available up to November 2005. So it is impossible to get conjunction between LANL and THEMIS due to the time coverage of THEMIS data. Consequently, in order to cross-calibrate all our data sets used in this study in the same way, we decided to use another method. It has been shown in previous studies that a correlation exists between fluxes at high latitudes (LEO orbit) and at the equator in the radiation belts [START_REF] Kanekal | Multisatellite measurements of relativistic electrons: Global coherence[END_REF]. So, two reference data sets are chosen, one near equator (equatorial pitch angle > 70°), from LANL spacecraft, and one at LEO, from POES spacecraft. Then, the goal is to calculate the factor between fluxes at LEO (POES) and those near equator (LANL for equatorial pitch angle >70°) and for nearly the same MLT, L and time for the two data sets mentioned just above and during quiet magnetic period. The following conditions have been used: -3h<MLT<3h (with MLT<1h between the two data sets), 6.2<L<7.2 (with L<0.2), t<24h and Kp<2 during two days before. Then a mean factor for each of the 5 energies (0.19, 0.84, 2.6, 8 keV and >30 keV) is calculated and is our reference. The same calculation is done using THEMIS and then using POLAR data instead of LANL data. Considering the same range of MLT, L, pitch angle and Kp conditions, the factor between electrons fluxes at LEO and near equator should be the same than the reference factors calculated with LANL data. If it is not the case, a correction will be apply to POLAR and THEMIS data in order to be consistent with POES and LANL data which are our references. Table 3 presents the correction factors to apply to POLAR and THEMIS data for each energy channel obtained using the method described above. Electron flux [cm -2 .s -1 .sr -1 .MeV -1 ]

POLAR data 60< α eq <90°-8 keV common time range (2008-2015 for THEMIS against 1998-2005 for LANL). However, although having a common time period with LANL, the POLAR data appear to be a little bit too low the correction factor being greater than 1.

B. How to take into account the limited time period of most of the data sets?

In a first step the solar cycle variation will not be included in the model. Assuming this tendency at LEO is also true at other orbits, POLAR data have to been lowered and THEMIS data have to be increased to provide an accurate and representative long term mean flux.

C. Variation of electron flux along magnetic field lines

The third step in the development of the model is to find a general fit, for each L bin (from L=3 to L=10), MLT and energy, allowing to calculate mean electron fluxes everywhere along a magnetic field line based on the data sets in use in this study. First, because fluxes provided by the existing GREEN model are fluxes integrated in energy, all data described before have been converted into five integrated channels: >0.19 keV, >0.84 keV, >2.6 keV, >8 keV and > 30 keV. Then, taking into account the contamination of the data used to develop the low energy plasma model below L=2 or 2.2 (according to data sets), we decided to construct a model between L=3 and L=10 to be sure of the provided flux. Loss cone α eq [°] Loss cone Then, Fig. 10 represents examples of electron fluxes versus equatorial pitch angle for each L bin considered in this study (35 bins from 3.0<L<3.2 in black to 9.8<L<10.0 in red), for E>0.84 keV and 6h<MLT<7h, on the top and E>8 keV and 12h<MLT<13h at the bottom. The color lines, correspond to the best fits for each L interval as plotted on Fig. 9 for the given energy and MLT plotted.

The top panel of Fig. 10 (E>0.84 keV) indicates that electron fluxes increase along the field line at all L values, and increase with L. However, the variation between equator and high latitudes is clearly higher for low L values (black curve) than for high L values (red curve). On the bottom panel, it is quite different because, as seen on Fig. 10, near equator fluxes at low L values are higher than those at high L values while at high latitudes fluxes increase with L. This difference between the two plots cannot only be attributed to the difference in energy but also to the different MLT plotted here. The behavior of low energy electron flux is energy and MLT dependent.

D. Mean electron model with MLT dependence

Using the fits determined for each of the 4 energies studied, each L interval (30 linear intervals from L=3 to 10) and each MLT interval (24 linear intervals), electron fluxes from >0.19 keV to >30 keV can be calculated anywhere in the inner magnetosphere, taking into account the MLT dependency. Fig. 11 shows an example of results with two MLT-L cartographies of flux obtained with this new model near the equator for E>2.6 keV electrons on the left and for E>8 keV on the right. As already mentioned, this model is valid from L=3 up to L=10 and depends on MLT.

IV. COMPARISON WITH SPM/AE9 MODEL AND DATA

A. Comparison with AE9/SPM

The SPM model, part of the AE9/AP9/SPM model distribution can be used to estimate low-energy plasma fluxes [1] [13]. This empirical model is mainly based on POLAR/HYDRA [START_REF] Roeder | Empirical models of the low-energy plasma in the inner magnetosphere[END_REF] and LANL/MPA [START_REF] Mccomas | Magnetospheric plasma analyzer: Initial three-spacecraft observations from geosynchronous orbit[END_REF] data and is valid from 1 keV to 40 keV. THEMIS/ESA data have been added in SPM Version 1.2 [START_REF] Johnston | Recent Updates to the AE9/AP9/SPM Radiation Belt and Space Plasma Specification Model[END_REF]. In this paper, AE9/SPM V1.5 is used. SPM is not MLT dependent but however be compared to the model described in this paper. In order to be sure that the model is not only usable near equator, comparison with data has been done outside equator. There are not many data sets at low energy outside equator that's why we have chosen to compare results of the model with POLAR data, even if this data set has been used to develop the model. It allows to control however if the model is well constructed outside equator. Fig. 14 represents E>1 keV electron fluxes (cm -2 .s -1 ) along POLAR orbit in 2003 versus L* and equatorial pitch angle α eq from POLAR data after cross calibration (on the top) and calculated by the new ONERA low-energy model (at the bottom). Let's remember that the ONERA low-energy model is only valid beyond L=3 so only this region is plotted here.

As seen on the figure POLAR orbit sweeps many L-shells at high latitudes (low pitch angles) as well as near the equator (large pitch angles). These plots show that for L* beyond 5, fluxes measured by POLAR, taking into account the crosscalibration, are globally coherent with fluxes provided by the model. We can note however that cross-calibrated POLAR fluxes seems to be slightly higher than those of the model but we have to keep in mind that 2003 is a "high year" for fluxes according to the solar cycle while the low-energy plasma model provide mean fluxes with no dependence on the solar cycle.

V. IMPLEMENTATION IN GREEN MODEL

The low-energy electron model described here has been implemented in the version 2 of GREEN (GREEN-V2). Fig. 15 represents the energy and L coverage of the different models integrated in GREEN-e. A detailed description of GREEN model can be found in Sicard et al., 2018 [3]. The low-energy model is orange on this plot. It replaces AE9/SPM beyond L=3 and is valid from 0.2 keV to 40 keV. Extrapolation is done from 30 keV up to 40 keV that is to say the lower energy limit of AE8. 15) are slightly higher than those provided by the low energy plasma model described in this paper. This difference can be explained by the fact that IGE-2006 depends on the solar cycle and that 2003 is a year with high fluxes, while the low energy model is not solar cycle dependent and is a mean model. On these plots, as mention in Sicard et al., 2018 [3], a simple smoothing function has been used (on the top plot of Fig. 16) to smooth existing discontinuities at the interface of "sub-models" implemented in GREEN. Next year, some big efforts will be done to apply more complex smoothing functions.

VI. CONCLUSION

A new low-energy plasma model has been developed for surface damages calculation. This empirical model is based on NOAA-POES, THEMIS, POLAR and LANL data. Each data set has been analyzed in several dimensions: L between 3 and 10 with a resolution of 0.2, MLT with a resolution of 1 hour and energy between 0.19 and 30 keV with five energy channels 0.19 keV, 0.84 keV, 2.6 keV, 8 keV and 30 keV. This MLT dependent model provides electron from >0.19 keV to >30 keV from L=3 up to L=10 and has been compared with AE9/SPM results and validated with Van Allen Probes data. This ONERA low energy model has been implemented in GREEN-e and replace AE9/SPM in this model beyond L=3. The next step is to develop the same kind of low-energy model for protons and to implement it in GREEN-p. However, only very few good quality protons data are available at low energy (<100 keV), because of contamination by high energy particles, which does not simplify the development of a model.

A

  new model for the plasma electrons fluxes (part of GREEN model) A. Sicard, D. Boscher, D. Lazaro, S. Bourdarie, D. Standarovski, R. Ecoffet T

Fig. 2 :

 2 Fig. 2 : L-time cartographies of 0.84 keV electrons flux measured by POES-15/TED before cleaning (in the middle), after cleaning (at the bottom) and 1 MeV electrons flux measured by POES-15/SEM (on the top).

Fig. 3 :

 3 Fig. 3 : Correlation between 0.84 keV electrons flux measured by POES-15/TED and 1 MeV electrons flux measured by POES-15/SEM for 2<L<7. The correlation is reflected by the points in the hatched part of the figure. All these points have been removed from the data

Fig. 4 :

 4 Fig. 4 : Cartographies MLT-L of electrons flux measured by POES/TED from 1998 to 2017 for 2 energies: 0.19 keV and 8 keV.

2 .

 2 NPOES-15/MEPED [cm -2 .s -1 .sr -1 ] Electron 0.84 keV NPOES-15/TED [cm -flux [cm -2 .s -1 .sr -1 .MeV -1 ] Electron flux [cm -2 .s -1 .sr -1 .MeV -1 ]

Fig. 5 :

 5 Fig. 5 : L-time cartographies of 700 keV electrons measured by THEMIS/SST (on the top), 1 MeV electrons flux measured by POES-15/SEM (second panel), 0.8 keV electrons flux measured by THEMIS-A/ESA before cleaning (third panel), after cleaning (at the bottom).

Fig. 6 :

 6 Fig. 6 : MLT-L cartographies of electrons flux measured by THEMIS/ESA from 2008 to 2015 for 2 energies: 0.19 keV and 8 keV.

Fig. 7 :

 7 Fig. 7 : L-time cartographies of 0.8 keV electrons flux measured by POLAR/HYDRA before cleaning (in the middle), after cleaning (at the bottom) and 1 MeV electrons flux measured by POES-15/SEM (on the top).

Fig. 8

 8 Fig.8represents MLT-L cartographies of fluxes measured by POLAR/HYDRA from 1998 to 2008 for 8 keV electrons near the equator (α eq >60°) on the right and at high latitudes (α eq <30°) on the left. These plots indicate that there is a dependence of low energy electron fluxes with latitude. Moreover, comparing the cartography from POLAR (on the right) to the one from THEMIS (Fig.6) in the equatorial region, 8 keV electron fluxes from POLAR are a little bit lower than those from THEMIS. A data cross-calibration is necessary and the procedure will be described in the next section of this paper.

  NPOES-15/SEM2 Electron 1 MeV THEMIS-A/ESA Electron 0.8 keV THEMIS-A/ESA Electron 0.8 keV THEMIS-A/SST Electron 700 keV THEMIS data -0.19 keV 10 9 10 10 Electron flux [cm -2 .s -1 .sr -1 .MeV -1 ] flux [cm -2 .s -1 .sr -1 .MeV -1 ] POLAR/CEPPAD Electron 3 MeV POLAR/HYDRA Electron 0.8 keV POLAR/HYDRA Electron 0.8 keV

Fig. 8 :

 8 Fig. 8 : MLT-L cartographies of flux measured by POLAR/HYDRA from 1998 to 2008 for 8 keV electrons near the equator (on the left) and at high latitudes (on the right).

  Electron flux [cm -2 .s -1 .sr -1 .MeV -1 ]

  So it is essential to take into account the different time coverage of the data sets being used in this study according to the solar cycle. Indeed, POES data cover two solar cycles while THEMIS and POLAR data cover only one solar cycle with a partial MLT coverage at each time. Taking into account the flux modulation as seen on Fig.1, it is necessary to calculate a correction factor to apply to THEMIS and POLAR data which will account for their solar cycle fraction of availability at each L, B/B Eq , and MLT. It is important to keep in mind that this low energy model (which is a part of GREEN) is a mean model to be used for long duration mission. Thus, mean fluxes over the longest time period as possible are expected to be more representative. In order to calculate this correction factor, ratios between the mean POES flux over the time period 1998-2018 and the mean POES flux over THEMIS time period (2008-2015) on one side and POLAR time period (1998-2008) on the other side, shall be computed for each energy, each L and each MLT values. These ratios correspond to the difference of long term dynamics of electron to shorter term time periods. The factors applied to THEMIS and POLAR are then defined as follow: Corr. factor (Ec, L, MLT)= 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓(𝑃𝑃𝐸𝐸,𝐿𝐿,𝑀𝑀𝐿𝐿𝑀𝑀) 𝑏𝑏𝑀𝑀𝑏𝑏𝑏𝑏𝑀𝑀𝑀𝑀𝑀𝑀 1998 𝑀𝑀𝑀𝑀𝑎𝑎 2018 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓 (𝑃𝑃𝐸𝐸,𝐿𝐿,𝑀𝑀𝐿𝐿𝑀𝑀) 𝑎𝑎𝑓𝑓𝑑𝑑𝑑𝑑𝑀𝑀𝑑𝑑 𝑠𝑠/𝐸𝐸 𝑝𝑝𝑀𝑀𝑑𝑑𝑑𝑑𝑝𝑝𝑎𝑎These correction factors are applied to THEMIS and POLAR electron fluxes, for each energy, L and MLT in order to harmonize data sets representativity with the long term POES data. The values of these factors (versus Ec, L and MLT) are not fully described here but one should note that they range between 0.5 and 1 for POLAR time period and between 1 and 2.5 for THEMIS time period. It means that in the region between L=3 and L=10, electron fluxes measured by POES and averaged on the time period of POLAR data considered in this study (1998-2008) are higher than the same fluxes averaged over the two solar cycles of POES data (1998-2017). On the contrary electron fluxes measured by POES and averaged over the time period of THEMIS data considered in this study(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) are lower than the same fluxes averaged over the two solar cycles of POES data (1998-2017).

Fig. 9 :

 9 Fig. 9 : Electron flux along a field line versus equatorial pitch angle for E>0.84 keV and 0h<MLT<1h. The points correspond to spacecraft data (after correction factors) and the color lines correspond to the best fit for three L bin for a given energy and a given MLT.Thus, Fig.9presents an example of electrons flux along a field line versus equatorial pitch angle for three L values: 3. 0<L<3.2, 5.0<L<5.2 and 7.0<L<7.2, for E> 0.84 keV. The points correspond to spacecraft data (after applying correction factors) and the color lines correspond to the best fit, of the form Y=exp(A.sin B (X)), for these three L bins and the given energy and MLT. This fit is the same as the one used in the SLOT model[START_REF] Sicard-Piet | A New ONERA-CNES Slot Electron Model[END_REF].

Fig. 10 :

 10 Fig. 10 : Electron flux along a field line versus equatorial pitch angle for E>0.84 keV and 6h<MLT<7h on the top and E>8 keV and 12h<MLT<13h at the bottom for each L bin (40 bins from 2.0<L<2.2 in black to 9.8<L<10.0 in red).

Fig. 11 :

 11 Fig. 11 : Cartographies MLT-L of flux obtained with the new model near equator for E>2.6 keV electrons on the left and for E>8 keV on the right.

  Fig. 12 shows comparisons between electrons flux obtained with Mean SPM model (in thick black line), SPM Percentile 99 (in thick grey line) and those provided by the ONERA plasma model (dotted lines) near equator versus radial distance for E>2.6 keV (on the top) and E>8 keV (at the bottom). The thick red line represents the MLT averaged flux.

Fig. 12 :

 12 Fig. 12 : Comparisons between electrons flux obtained with SPM Mean model (in thick black line), SPM Percentile 99 (in thick grey line) and those provided by the ONERA plasma model (dotted lines) near equator versus radial distance for E>2.6 keV (on the top) and E>8 keV (at the bottom).It is important to note that results from SPM are plotted in terms of the radial distance while those from our new model are plotted in terms of L McIlwain parameter, which can lead to few differences but not significant at the equator. Comparing these plots, we can observe that flux from SPM are close to MLT averaged flux from ONERA model between L=3 and L=6.6 and close to maximum flux from ONERA model beyond L=6.6.

Fig. 13 :

 13 Fig. 13 : Flux versus MLT deduced from several data sets, obtained with AE9/SPM and calculated with the new model for E>8 keV electrons and 4.0<L<4.2 at the bottom and E>2.6 keV and 5.8<L<6.0 on the top.

Fig. 14 :

 14 Fig. 14 : E>1 keV electron fluxes (cm -2 .s -1 ) along POLAR orbit in 2003 versus L* and equatorial pitch angle αeq from POLAR data after cross-calibration (on the top) and calculated by the new ONERA low-energy model (at the bottom).

Fig. 15 :

 15 Fig. 15 : Energy and L coverage of the different models integrated in GREEN-e.

0

  It is important to note than the low energy model is developed in function of L McIlwain parameter while GREEN is based on L*. However a conversion between L and L* is done for models to be consistent. Moreover, let's remember that the low energy model is MLT dependent but have no dependence with the solar cycle in contrary to most of the models implemented in GREEN. Not enough data is available to construct a MLT and solar cycle dependent model at low energy.

Fig. 16 :

 16 Fig. 16 : Electron flux provided by GREEN-V2 at equator versus energy and L* in 2003 on the top and a zoom at the bottom In order to see if the ONERA low energy model is well implemented in GREEN, Fig. 16 shows electron flux provided by GREEN-V2 at equator versus energy and L*, from 1 keV and 10 MeV in 2003 on the top and a zoom of this plot at the bottom. Note that in this plot, the ONERA low-energy model has been averaged in MLT. On the zoom plot (at the bottom) we can see that near geostationary orbit, fluxes provided by IGE-2006 (cf Fig. 15) are slightly higher than those provided by the low energy plasma model described in this paper. This difference can be explained by the fact that IGE-2006 depends on the solar cycle and that 2003 is a year with high fluxes, while the low energy model is not solar cycle dependent and is a mean model. On these plots, as mention in Sicard et al., 2018 [3], a simple smoothing function has been used (on the top plot of Fig. 16) to smooth existing discontinuities at the interface of "sub-

Table 1

 1

resumes the characteristics of data used at LEO. All energy channels mentioned in the third column have been studied and used for each NOAA spacecraft.

TABLE I CHARACTERISTICS

 I OF DATA USED AT LEO

	Spacecraft	Time coverage	Channels
	POES-15/TED & MEPED	07/1998now	0.19 keV
	POES-16/TED & MEPED	10/200006/2014	0.84 keV
	POES-17/TED & MEPED	07/200204/2013	2.6 keV
	POES-18/TED & MEPED	06/2005now	8 keV
	POES-19/TED & MEPED	02/2009now	>30 keV

Table 1 :

 1 Characteristics of electron data used at LEO. All energy channel of the third column have been used for each NOAA spacecraft.

Table 2 :

 2 Characteristics of electron data used at other orbits. Data sets in light grey are used to develop the model. Data sets in medium grey are used to calibrate data. Data sets in dark grey is used to validate the model.

		TABLE II
	CHARACTERISTICS OF DATA USED AT OTHER ORBITS
	Satellite/	Time coverage	Channels used
	Instrumentt
	POLAR/HYDRA	01/199803/2008	29 channels between 12
			eV and 18 keV
	POLAR/CEPPAD	01/199803/2008	7 channels between 19
			keV and 212 keV
	THEMIS/ESA	01/200812/2015	21 channels between 7
			eV and 26 keV
	THEMIS/SST	01/200812/2015	7 channels between 31
			keV and 200 keV
	LANL/MPA	01/199811/2005	23 channels between 100
			eV and 34 keV
	LANL/SOPA	01/199811/2005	4 channels between 50
			keV and 250 keV
	VAP/HOPE	01/201301/2018	72 channels between 15
			eV and 50 keV
	VAP/MAGEIS	01/201301/2018	7 channels between 32
			keV and 226 keV

Table 3 :

 3 Correction factors to apply to THEMIS and POLAR data sets using LANL and POES data as references.

			TABLE III
	CORRECTION FACTORS TO APPLY TO POLAR AND THEMIS DATA
		0.19 keV	0.84 keV	2.6 keV	8 keV	>30 keV
	THEMIS	1.01	0.78	0.821	0.96	1.0
	POLAR	2.71	2.29	2.79	1.36	5.15

It is important to note that the correction factors for THEMIS data are close to one, which proves that THEMIS and LANL data are coherent with each other although they do not have a

ACKNOWLEDGMENT We acknowledge NASA contract NAS5-02099 and V. Angelopoulos for use of data from the THEMIS Mission,D. Larson and the late R. P. Lin for use of SST data and C. W. Carlson and J. P. McFadden for use of ESA data. Processing and analysis of the [HOPE, MagEIS, REPT, or ECT] data was supported by Energetic Particle, Composition, and Thermal Plasma (RBSP-ECT) investigation funded under NASA's Prime contract no. NAS5-01072. We acknowledge Harlan Spence (PI) and/or Geoff Reeves (Science Team Lead) for use of ECT data. We thank the NOAA's National Geophysical Data Center (NGDS) for providing NOAA POES data, the POLAR Hydra team for the use of HYDRA data and the LANL MPA team for the use of MPA data.

This work was supported by CNES, the French Space Agency (R-S18/MT-0003-176).


T-H Hubert Chan

Mauro Sozio

Distributed Approximate k-Core Decomposition and Min-Max Edge Orientation: Breaking the Diameter Barrier

Keywords: distributed algorithms, coreness, round complexity

We design distributed algorithms to compute approximate solutions for several related graph optimization problems. All our algorithms have round complexity being logarithmic in the number of nodes of the underlying graph and in particular independent of the graph diameter. By using a primal-dual approach, we develop a 2(1 + )-approximation algorithm for computing the coreness values of the nodes in the underlying graph, as well as a 2(1 + )-approximation algorithm for the min-max edge orientation problem, where the goal is to orient the edges so as to minimize the maximum weighted in-degree. We provide lower bounds showing that the aforementioned algorithms are tight both in terms of the approximation guarantee and the round complexity. Finally, motivated by the fact that the densest subset problem has an inherent dependency on the diameter of the graph, we study a weaker version that does not suffer from the same limitation.

I. INTRODUCTION

The k-core decomposition and algorithms for finding densest subgraphs have proved to be a valuable tool in graph mining and data analysis, with applications encompassing sociology, bioinformatics as well as graph visualization. Coreness. Formally, a node v in a (weighted) graph G is said to have coreness value k, if k is the largest number such that v belongs to a subgraph of G with minimum (weighted) degree equal to k [START_REF] Stephen B Seidman | Network structure and minimum degree[END_REF].

Intuitively, nodes belonging to well-connected communities tend to have high coreness values. Besides, several definitions of density have been proposed in the literature. In our work, we focus on the average degree density which is defined as the ratio between the number of edges and the number of nodes in a graph [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF]. One of the appealing properties of such a definition is that a densest subgraph (in terms of average degree density) can be computed in polynomial time. Such an optimization problem is often referred to as the densest subset problem. Recently, the diminishingly-dense decomposition has been introduced and studied [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF], [START_REF] Khuller | On finding dense subgraphs[END_REF], [START_REF] Tatti | Density-friendly graph decomposition[END_REF], which elegantly merges the concepts of density and k-core decomposition. Such a decomposition assigns to each node a real number, which we refer to as maximal density value. We also study T-H. Hubert Chan was partially supported by the Hong Kong RGC under the grants 17200418.

the min-max edge orientation problem [START_REF] Venkateswaran | Minimizing maximum indegree[END_REF] which turns out to be related to the previous problems. The goal is to assign an orientation to every edge such that the maximum weighted in-degree is minimized.

In our work, we develop distributed algorithms for approximating the coreness values, the local density values, as well as to find an approximate solution for the min-max edge orientation problem. Moreover, we study a weaker version of the densest subset problem. We envision the following applications for our work. Our distributed algorithms can be executed by agents in a social network or a P2P network, for example, so as to collect relevant statistics of the agents with respect to the underlying graph. Users in a social network with large coreness value are known to have "good spreading" properties in epidemiological studies [START_REF] Kitsak | Identification of influential spreaders in complex networks[END_REF]. Therefore, the coreness value (or an approximation) can be leveraged to maximize the spreading of a diffusion protocol. Communities in a social network consist of set of users sharing similar interests, such as hiking, traveling or photography. The density of a given subgraph can be used to measure how likely the corresponding users belong to a same community [START_REF] Yang | Defining and evaluating network communities based on ground-truth[END_REF]. Our work allows to approximate such a metric in a distributed fashion. Our distributed algorithms can also be used in distributed graph processing systems [START_REF] Malewicz | Pregel: a system for large-scale graph processing[END_REF] to process very large graphs not fitting into the main memory of one single machine.

We are given an undirected edge-weighted graph G = (V, E, w). We consider the classical Local model, in which every node can directly communicate only with its neighbors in synchronous rounds. Moreover, we assume each node knows (an upper bound on) the number n = |V | of nodes. The hop-diameter D (or diameter) of a graph G is at most if every pair of nodes in G can be connected with a path consisting of at most edges. We use the convention that the approximation ratio γ is at least 1. 1 There is often a tradeoff between the number of rounds and the approximation guarantees of a distributed algorithm. where the maximum in-degree is 1; any other orientation of edges incident on v will result in a maximum in-degree of at least 2. It takes Ω(n) rounds for node v to distinguish between the different graphs.

Ω(D) communication rounds. Similarly, for a node to be aware of whether it is included in an approximate densest subgraph, Ω(D) communication rounds are required. This follows from the fact that such a node has to verify whether a subgraph with higher density (possibly many hops away) is included in the graph or not. Hence, it is natural to investigate the following question: Can we devise distributed approximation algorithms for the aforementioned problems, while requiring a number of communication rounds independent of the diameter?

Although some of the problems discussed above have been studied in a distributed environment, there is no work focusing on breaking the diameter barrier, to the best of our knowledge. All our algorithms require a number of communication rounds logarithmic in n, while we provide tight lower bounds expressing the tradeoff between communication rounds and approximation ratio.

A. Our Results and Contribution

Coreness Values. Just like almost every work on these related problems [START_REF] Bahmani | Densest subgraph in streaming and mapreduce[END_REF], [START_REF] Batagelj | An o(m) algorithm for cores decomposition of networks[END_REF], [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF], [START_REF] Epasto | Efficient densest subgraph computation in evolving graphs[END_REF], [START_REF] Giatsidis | Dcores: Measuring collaboration of directed graphs based on degeneracy[END_REF], [START_REF] Gillet | A fully asynchronous and fault tolerant distributed algorithm to compute a minimum graph orientation[END_REF], [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF], [START_REF] Montresor | Distributed k-core decomposition[END_REF], [START_REF] Wen | I/O efficient core graph decomposition at web scale[END_REF], we consider an elimination procedure that repeatedly "peels off" nodes with small (weighted) degrees. In particular, we follow the interpretation by Montresor et al. [START_REF] Montresor | Distributed k-core decomposition[END_REF]. Specifically, given a threshold value b, in each round, nodes with degree less than b in the remaining graph are removed. Using a compact representation, we can imagine that the elimination procedure for all possible thresholds are run in parallel, where after each iteration, each node v just remembers the largest threshold b v (that we call the surviving number), for which it still survives.

It is known that after each iteration, the surviving number of each node v is at least its coreness c(v); moreover, after n rounds, the surviving number will reach the coreness value2 . However, to the best of our knowledge, so far there is no approximation analysis for the process in terms of the number of iterations.

Densest Subset. It is interesting that a variant of the elimination process was considered by Bahmani et al. [START_REF] Bahmani | Densest subgraph in streaming and mapreduce[END_REF] to give a 2(1+ )approximation for the densest subset streaming model. In each iteration or pass, the threshold is chosen to be 2(1 + ) times the density of the current subset of surviving nodes. Then, the process terminates in O(log 1+ n) iterations, and the subgraph from one of the iterations gives a 2(1 + )-densest subset. This inspires us that density can provide the right tool to design and analyze distributed approximation algorithms for coreness values. However, an immediate issue is that for every node to know the (approximate) density of the current subgraph, we already need Ω(D) rounds. Local Notion of Density. We observe that it is not necessary to use the global density of the (sub)graph, because in some sense, coreness value measures the local density of a node. Recently, Tatti and Gionis [START_REF] Tatti | Density-friendly graph decomposition[END_REF] considered a so-called locallydense graph decomposition, which has been further studied by Danisch et al. [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF], who defined a quantity known as maximal density for each node. This decomposition has actually been considered in passing by Khuller and Saha [START_REF] Khuller | On finding dense subgraphs[END_REF] in the context of finding densest subsets with large sizes. Intuitively, the decomposition works by repeatedly peeling off maximal densest subsets as follows. Given a weighted graph G, the maximal densest subset B forms the first layer; every node u ∈ B has maximal density equal to ρ G (B), even though no special importance is given to this value in [START_REF] Khuller | On finding dense subgraphs[END_REF]. Next, remove the nodes in B to form a quotient graph G , where an edge between u ∈ B and v ∈ B becomes a self-loop at v in G . Then, the procedure is recursively applied to G until all nodes are removed. Equipped with this notion of maximal density, we can adapt previous analysis to achieve the following result.

Theorem I.1 (Gracefully Degrading Approximation Ratio for coreness values). After running the compact elimination procedure in the Local model for T rounds, the surviving number at each node gives a 2n 1 T -approximation to its coreness value (and maximal density).

In particular, if every node knows the number of nodes (or an upper bound) n and would like to achieve 2(1 + )approximation, then T = log 1+ n rounds are sufficient.

Matching Lower Bound. For γ ≥ 2, by considering a γ-ary tree, we show in Section III that achieving γ-approximation requires Ω( log n log γ ) rounds. Min-Max Edge Orientation Problem. The centralized version of the problem was proposed by Venkateswaran [START_REF] Venkateswaran | Minimizing maximum indegree[END_REF] and has applications in telecommunication network design; it is known that the special case of unit edge weights can be solved in polynomial time. The connection of the problem with densest subsets has been explored in subsequent works [START_REF] Asahiro | Approximation algorithms for the graph orientation minimizing the maximum weighted outdegree[END_REF], [START_REF] Asahiro | Graph orientation algorithms to minimize the maximum outdegree[END_REF].

To the best of our knowledge, the only distributed algorithm for this problem with round complexity independent of the diameter is for unweighted graphs. Specifically, Barenboim and Elkin [START_REF] Barenboim | Sublogarithmic distributed MIS algorithm for sparse graphs using nash-williams decomposition[END_REF] actually studied a stronger problem, where the goal is to partition the edges of a graph into a minimum number of forests. Instead of density ρ(S) = w(E(S))

|S|

, a similar notion of arboricity arb(S) = w(E(S))

|S|-1 is considered. In hindsight, it is not surprising that they also used a variant of the distributed elimination procedure to approximate the forest-decomposition problem in the Local model. Indeed, they showed that if the maximum arboricity is known by every node, then O( 1 log n) rounds are sufficient to achieve (2+ )-approximation. As remarked above, if every node needs to know the (approximate) maximum arboricity, then Ω(D) rounds are required. A careful study of their algorithm reveals that the first phase [5, Algorithm 3] serves a similar purpose as computing the surviving numbers as in our Theorem I.1, after which they run the second phase as if the maximum arboricity is known. This degrades the quality of the solution, and only 2(2 + )-approximation is achieved.

Primal-Dual Approach. By observing that the LP relaxation of the min-max edge orientation problem is the dual of the densest subset LP (see Section II), we have the intuition that a procedure for approximating the maximal densities should also give information about the dual problem, without using a second phase. Indeed, we augment the distributed elimination procedure by maintaining an auxiliary subset N v for every node v (which represents its in-neighbors) to give the same approximation ratio for the edge orientation problem. However, a very careful invariant analysis in Lemma III.11 is performed to make sure that every edge is taken care of by at least one of its end-points. In particular, if every node knows n and would like to achieve 2(1 + )-approximation, then T = log 1+ n rounds are sufficient.

Densest Subset Problem. As argued above, for every node to be aware of whether it should be included in an approximate densest subset, it takes Ω(D) rounds. Indeed, Sarma et al. [START_REF] Das Sarma | Dense subgraphs on dynamic networks[END_REF] gave a distributed algorithm that gives

2(1 + )-approximation with O(D log 1+ n) rounds.
We consider a weaker version of the problem (in Definition IV.1). Instead of just producing one subset, the distributed algorithm will return a collection {S i : i ∈ I} of disjoint subsets such that for each i ∈ I, each node will be aware of which subset (if any) it belongs to; moreover, there exists some i ∈ I such that S i is an approximate densest subset. Using the procedure in Theorem I.1 as a subroutine, we have the following.

Theorem I.3 (Distributed (Weak) Densest Subset Problem).

For any > 0 and n, there exists a distributed algorithm that gives a 2(1 + )-approximation to the (weak) densest subset problem in Definition IV.1 on any graph with n nodes using O(log 1+ n) rounds.

B. Related Work

In addition to the most relevant works that have already been compared with our work, we also describe other related works. Corenss Values. After Seidman [START_REF] Stephen B Seidman | Network structure and minimum degree[END_REF] first proposed k-core decomposition, Batagelj and Zaversnik [START_REF] Batagelj | An o(m) algorithm for cores decomposition of networks[END_REF] presented a (centralized) O(m) algorithm to compute a k-core decomposition. This notion has been extended to weighted graphs [START_REF] Cheng | Efficient core decomposition in massive networks[END_REF], [START_REF] Giatsidis | Evaluating cooperation in communities with the k-core structure[END_REF] and directed graphs [START_REF] Giatsidis | Dcores: Measuring collaboration of directed graphs based on degeneracy[END_REF]. The distributed setting has been studied by Montresor et al. [START_REF] Montresor | Distributed k-core decomposition[END_REF], and was further extended to dynamic graphs by Aridhi et al. [START_REF] Aridhi | Distributed k-core decomposition and maintenance in large dynamic graphs[END_REF]. The distributed algorithms have been adapted to (centralized) I/O efficient algorithms [START_REF] Cheng | Efficient core decomposition in massive networks[END_REF], [START_REF] Wen | I/O efficient core graph decomposition at web scale[END_REF] that were proposed to handle large graphs that cannot fit into memory.

Min-Max

Edge Orientation Problem. Since Venkateswaran [START_REF] Venkateswaran | Minimizing maximum indegree[END_REF] introduced the problem, (centralized) polynomial-time algorithms are known to give optimal solutions for unweighted graphs [START_REF] Asahiro | Graph orientation algorithms to minimize the maximum outdegree[END_REF], [START_REF] Kowalik | Approximation scheme for lowest outdegree orientation and graph density measures[END_REF], [START_REF] Venkateswaran | Minimizing maximum indegree[END_REF]. However, a series of works [START_REF] Asahiro | Approximation algorithms for the graph orientation minimizing the maximum weighted outdegree[END_REF], [START_REF] Asahiro | Graph orientation algorithms to minimize the maximum outdegree[END_REF] showed that the weighted version is NP-hard even when all edge weights belong to the set {1, k}, where k is any fixed integer greater than 1; on the other hand, for integer edge weights, (2 -1 k )-approximation can be achieved, where k is the maximum weight. Gillet and Hanusse [START_REF] Gillet | A fully asynchronous and fault tolerant distributed algorithm to compute a minimum graph orientation[END_REF] considered the more general asynchronous distributed model with faults. However, their algorithm has round complexity depending on the graph diameter, and achieves 2(2 + )-approximation.

The problem can be also be seen as a special case of a load-balancing task, where each node is a machine and each edge is a job to be assigned to one of its incident machines. From this perspective, minimizing the maximum in-degree is equivalent to minimizing the makespan. Czygrinow et al. [START_REF] Czygrinow | Distributed 2-approximation algorithm for the semi-matching problem[END_REF] considered minimizing a slightly different objective that is the sum of the squares of the loads of the machines. In contrast, they gave a distributed 2-approximation algorithm that runs in O(∆ 5 ) rounds, where ∆ is the maximum degree.

The Min-Max Edge Orientation Problem share some similarities with the vertex cover with hard capacities problem [START_REF] Grandoni | A primal-dual bicriteria distributed algorithm for capacitated vertex cover[END_REF]. Densest Subsets. The densest subset problem [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF] has been extensively studied and extended to directed graphs [START_REF] Kannan | Analyzing the structure of large graphs[END_REF]. As mentioned above, even though Khuller and Saha [START_REF] Khuller | On finding dense subgraphs[END_REF] have implicitly considered locally-dense graph decompositions, this notion has not been fully studied until recently by Tatti and Gionis [START_REF] Tatti | Density-friendly graph decomposition[END_REF] and Danisch et al. [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF]. One of its surprising applications is its usage in the computation of the (non-linear) Laplacian operator of a hypergraph [START_REF] Hubert Chan | Spectral properties of hypergraph laplacian and approximation algorithms[END_REF]. II. PRELIMINARIES Distributed Model. The input is an edge-weighted undirected graph G = (V, E, w), where n = |V | and the edge weights w : E → R + are non-negative. We consider the classical Local model, where each node in V is a CPU with a unique identifier that is aware of only its incident edges (and their weights) and its neighbors; each node knows the total number n of nodes (or an upper bound). In addition to the Local model, our protocols satisfy the following.

• Synchronous Rounds and Polynomial-Time Computation.

In each round, a node can send a message to its neighbors in G. Moreover, after receiving messages from its neighbors, a node can perform polynomial-time computation in each round. • Broadcast Model. We consider protocols in which a node sends the same message to (a subset of) its neighbors in each round. • Message Content and Size. In our protocols, we assume that each message contains the identity of the sender, and typically, the content of a message is a constant number of real numbers. In most useful applications, each edge weight is an integer whose size is polynomial in n. In this case, every number sent in a message in our protocols can be represented by O(log n) bits, which means the requirements of the Congest model are also satisfied. Alternatively, with arbitrary edge weights, by restricting the set of numbers to appropriate powers of (1 + λ), we can also restrict the message size. Graph Terminology. In most cases, each edge e ∈ E is interpreted as a 2-subset of V , although in the analysis we sometimes consider a self-loop, which is a singleton.

For each node v ∈ V , the set of neighbors of v is N G (v) := {u ∈ V : u = v ∧ {u, v} ∈ E}, and the weighted degree of v is deg G (v) := e∈E: v∈e w(e), i.e., the sum of the weights of the edges that contain v. When there is no ambiguity, the subscript G can be omitted. Density. For a non-empty S ⊆ V , we denote E(S) := {e ∈ E : e ⊆ S}, and its density is ρ G (S) = w(E(S))

|S|

. A (nonempty) subset S in G is a densest subset if and only if S has maximum density among all subsets in G. The following fact is standard.

Fact II.1. The maximal densest subgraph of G is unique and contains all densest subgraphs of G.

Local Density. While a densest subset gives the densest region in a graph G, there are other notions that measure the local density around each node v. The first notion is coreness c G (v), which is mentioned in Section I. Another notion is maximal density, which is defined in terms of quotient graph and diminishingly-dense decomposition as follows.

Definition II.2 (Quotient Graph [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF]). Given a weighted undirected graph G = (V, E, w) and a subset B ⊆ V , the quotient graph of G with respect to B is a graph G\B := ( V , E, w), where V := V \B and E := {e ∩ V : e ∈ E, e ∩ V = ∅}, i.e., every edge e ∈ E not contained in B contributes towards E. Moreover, for e ∈ E, w(e ) := e∈E:e =e∩ V w(e).

Definition II.3 (Diminishingly-Dense Decomposition and Maximal Density [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF]). Given a weighted undirected graph G = (V, E, w), define the diminishingly-dense decomposition B of G as the sequence

∅ = B 0 B 1 • • • B k = V as follows.
Initially we set B 0 := ∅ and G 0 := G. For i ≥ 1, if B i-1 = V , the decomposition is fully defined. Otherwise, define G i := G i-1 \B i-1 and let S i be the maximal densest subset in G i .

Then, define B i := B i-1 ∪ S i . Finally, for each node v ∈ V , we say that the maximal density of

v is r G (v) := ρ Gi (S i ) if v ∈ S i . For simplicity, we use r(v) to denote r G (v) if the context is self-evident.
Fact II.4 (Strictly Diminishing Densities). In Definition II.3, the sequence {ρ Gi (S i )} i is strictly decreasing in i.

Distributed Approximation of Local Density. Our goal is to design distributed protocols such that at the end, each node v outputs some number β(v), which is an approximation of its coreness c(v) or maximal density r(v). Observe that there are example graphs such that computing c(v) or r(v) exactly needs Ω(n) rounds. The main result of this paper is that there are protocols with O(log n) rounds that give O(1)-approximation for both c(v) and r(v) for each node v.

To be precise, we use the following convention to describe approximation ratio.

Definition II.5 (Approximation Ratio). Given a non-negative real number s ∈ R + and γ ≥ 1, another number β is a γapproximation for s, if s ≤ β ≤ γ • s. In general, a function

β : V → R + is a γ-approximation for s : V → R + , if for all v ∈ V , s(v) ≤ β(v) ≤ γ • s(v).
Min-Max Edge Orientation Problem. The input is an edgeweighted undirected graph G = (V, E, w). The goal is to give an orientation to each edge such that the maximum weighted in-degree of a node is minimized. Equivalently, we wish to determine an assignment a : E → V such that each edge e is assigned to one of its end-points, and max v∈E e∈a -1 (v) w e is minimized.

Distributed Setting. After the distributed algorithm terminates, each node v ∈ V will have computed some subset N v of its neighbors, which represents the set of incident edges that are assigned to it. Observe that it is sufficient to enforce the condition that for every edge {u, v} ∈ E, we have u ∈ N v or v ∈ N u ; by means of one more round of communication, we can resolve any conflict (i.e. an edge being assigned to both its end-points).

LP Relaxation and Relationship with Densest Subset. The following is an LP relaxation of the problem together with its dual. In the primal LP, for each e ∈ E and node u ∈ e, α e u is the portion of the weight w e of edge e that is assigned to u. As observed in [START_REF] Danisch | Large scale density-friendly graph decomposition via convex programming[END_REF], the dual LP is exactly the densest subset LP by Charikar [ Approximation Analysis. As we shall see, a by-product of our distributed algorithm is that it gives an assignment of edges to their incident nodes. If ρ * is the maximum density of a subset, and the sum of the edge weights assigned to every node is at most γ • ρ * , then weak duality implies that the assignment is a γ-approximation to the min-max edge orientation problem.

III. DISTRIBUTED APPROXIMATION ALGORITHMS FOR CORENESS VALUES (AND MIN-MAX EDGE ORIENTATION)

In this section, we give a distributed protocol that approximates both the coreness value and maximal density for each node in a given graph.

A. Warmup: Single Threshold

The warmup protocol is based on a well-known idea of iteratively eliminating vertices with small degree [START_REF] Bahmani | Densest subgraph in streaming and mapreduce[END_REF], [START_REF] Batagelj | An o(m) algorithm for cores decomposition of networks[END_REF], [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF], [START_REF] Giatsidis | Dcores: Measuring collaboration of directed graphs based on degeneracy[END_REF], [START_REF] Gillet | A fully asynchronous and fault tolerant distributed algorithm to compute a minimum graph orientation[END_REF], [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF], [START_REF] Montresor | Distributed k-core decomposition[END_REF], [START_REF] Wen | I/O efficient core graph decomposition at web scale[END_REF], which can be easily implemented in the distributed model. Here we consider a protocol that is parameterized by some universal threshold b. In each round, each node with weighted degree less than b in the subgraph induced by surviving nodes is marked to be removed at the end of the round. It is clear that after running the protocol for n rounds, all surviving nodes have coreness at least b [START_REF] Giatsidis | Dcores: Measuring collaboration of directed graphs based on degeneracy[END_REF]. It is shown in [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF] that running the protocol for various thresholds and O(log n) rounds can be used to give O(1)approximation for the densest subgraph problem. We show that the analysis can be adapted to approximate both coreness values and maximal density.

We describe the elimination procedure in the distributed model in Algorithm 1. Each node v ∈ V keeps a state σ v ∈ {0, 1} which records whether the node is present [START_REF] Aridhi | Distributed k-core decomposition and maintenance in large dynamic graphs[END_REF] or removed (0).

Algorithm 1: Elimination Procedure for a Single Threshold

Input: G = (V, E, w), threshold value b ∈ R,

B. Parallel Execution with Multiple Thresholds

Parallel Execution. Observe that the elimination procedure for different threshold values can be executed in parallel, but this can cause a large message size. Hence, for the purpose of analysis, we imagine that the protocol is executed for all possible thresholds in parallel and concentrate on the round complexity. Later in Section III-C, we show how parallel execution for all threshold values can be performed compactly. Depending on this parallel thought experiment, we define the surviving number as follows.

Definition III.1 (Surviving Number). Given a weighted undirected graph G = (V, E, w), the surviving number of a node v ∈ V after T rounds, denoted as β T G (v), is defined as the maximum b ∈ R such that v survives after T rounds of the elimination procedure using threshold value b. When the context is clear, we omit G or T from the notation and simply use β(v).

It is known that after running the procedure for n rounds, β n (v) = c(v) gives the exact coreness value [START_REF] Montresor | Distributed k-core decomposition[END_REF]. Our goal is to show that to get a constant approximation, it suffices to use O(log n) rounds, independent of the diameter of the graph. In particular, we shall prove that for > 0, if we set T := log 1+ n , then β T (v) is a 2(1 + )-approximation for both the coreness value c(v) and maximal density r(v).

Lemma III.2 (Lower Bound on Surviving Number).

For any node v ∈ V and any positive integer t, β t (v) ≥ c(v).

Proof. For any node v with coreness c(v), by the definition of coreness, suppose S ⊆ V is a subset such that v ∈ S and deg G[S] (u) ≥ c(v) for all u ∈ S. Then, for b = c(v), all members of S will survive the elimination procedure with threshold b, no matter how many rounds of elimination are executed. In particular, v will survive. Hence, β t (v) ≥ c(v) for all positive integers t.

Lemma III.3 (Upper Bound on Surviving Number). For any positive integer T and any node

v ∈ V , we have β T (v) ≤ 2n 1 T • r(v).
In other words, to achieve an approximation guarantee of γ > 2, it suffices to set T = log n log(γ/2) ; in particular, for > 0, the special case

T := log 1+ n gives β T (v) ≤ 2(1 + ) • r(v).
Proof. The proof approach has appeared in several previous works [START_REF] Bahmani | Densest subgraph in streaming and mapreduce[END_REF], [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF], [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF]. For completeness, we adapt the proof from [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF]Lemma 3.1].

Fix some v ∈ V . Consider the diminishingly-dense decomposition of G. Let G i and S i be defined as in Definition II.3 such that v ∈ S i , which is the maximal densest subset in

G i = (V i , E i ). Then, r(v) = ρ Gi (S i ). Observe that G i = G \ V i is a quotient graph, which means that any edge connecting V i to V i becomes a self-loop in V i , which implies that for all v ∈ V i , for all t, β t G (v) ≤ β t Gi (v). Fix any b > 2n 1 T •r(v).
It suffices to show that after applying T rounds of elimination with threshold b to the graph G i , no node in V i (which includes v) can survive.

Define A 0 := V i ; for j ≥ 1, let A j ⊆ V i be the set of nodes that survive after round j. Suppose for some j ≥ 0, both A j and A j+1 are non-empty. Then, since S i is a densest subset in G i , we have r

(v) = ρ Gi (S i ) ≥ ρ Gi (A j ) = w(E(Gi[Aj ])) |Aj | ≥ u∈A j deg G[A j ] (u) 2•|Aj | ≥ b•|Aj+1| 2•|Aj | > n 1/T •r(v)•|Aj+1| |Aj |
.

Note that we use the inequality w(E(G

i [A j ])) ≥ u∈A j deg G[A j ] (u) 2 , because G i could contain self-loops. Hence, we have |A j+1 | < |Aj | n 1/T . Therefore, if A T is non-empty, we have |A T | < |A0|
n ≤ 1, which is a contradiction. So A T is empty, which means that no node in V i survives after round T , as required.

Lemma III.4 (Relating Maximal Density and Coreness Value). For any

v ∈ V , r(v) ≤ c(v).
Proof. Consider the diminishingly-dense decomposition of G. Let {(G i , B i , S i )} i be defined as in Definition II.3. Recall that there is a unique i such that v ∈ S i ⊆ B i , and S i is the maximal densest subset in G i .

We prove by induction on i ≥

1 that deg G[Bi] (v) ≥ ρ Gi (S i ) for all v ∈ B i , which is slightly stronger than the required statement. From this, it follows that c(v) ≥ deg G[Bi] (v) ≥ ρ Gi (S i ) = r(v) for all v ∈ B i , i ≥ 1. Base Case. Suppose v ∈ S 1 = B 1 . Since S 1 is a densest subset in G 1 , we must have deg G[B1] (v) ≥ ρ G1 (S 1 ), as required.
Otherwise, the subset S 1 \{v} would have higher density than

S 1 in G 1 . Induction Hypothesis. Suppose for some i > 1, for all u ∈ B i-1 , deg G[Bi-1] (u) ≥ ρ Gi-1 (S i-1 ). Inductive Step. Consider v ∈ S i . Note that B i = B i-1 ∪ S i . For u ∈ S i , since G i = G \ B i-1 is a quotient graph, it follows that deg G[Bi] (u) = deg Gi[Si] (u)
, because any edge connecting S i to B i \S i will become a self-loop in G i . Similar to the base case, because S i is a densest subset in G i , we have

deg Gi[Si] (u) ≥ ρ Gi (S i ) = r(u), as required. For u ∈ B i-1 , we have deg G[Bi] (u) ≥ deg G[Bi-1] (u) ≥ ρ Gi-1 (S i-1 )
, where the second inequality follows from the induction hypothesis. By Fact II.4, the last term satisfies ρ Gi-1 (S i-1 ) > ρ Gi (S i ), thereby completing the inductive step and the proof of the lemma.

Theorem III.5. Given G = (V, E, w) and γ > 2, let T := log n log(γ/2) . Then, for all v ∈ V , we have r(v) ≤ c(v) ≤ β T (v) ≤ γ • r(v) ≤ γ • c(v).
In particular, for any > 0, we can set γ = 2(1 + ) and T := log 1+ n .

Proof. This follows directly from Lemma III.2, Lemma III.3 and Lemma III.4.

Corollary III.6 (Relating Coreness Value and Maximal Den

- sity). For each node v, r(v) ≤ c(v) ≤ 2 • r(v).

C. Compact Parallel Execution

Observe that it is inefficient (or infeasible) to naively execute the elimination procedure for all threshold values in parallel. Instead, at any moment, a node v just needs to remember the maximum threshold value β(v) for which it still survives in the corresponding elimination procedure. Hence, in each round, a node v just needs to send its current value β(v) to all its neighbors. Indeed, this observation is made by Montresor et al. [START_REF] Montresor | Distributed k-core decomposition[END_REF] to produce a compact algorithm, which we restate in Algorithm 2. However, they run the algorithm till the exact coreness value for each node is achieved, while we have already shown that for any γ > 2, T = log n log(γ/2) rounds are sufficient to obtain γ-approximation for both the coreness value and maximal density. Message Size. If every edge has integer weight that is polynomial in n, then each message has size O(log n) bits. For general edge weights, to further optimize for the message size, we can restrict the sent numbers from some set Λ, thereby achieving log 2 |Λ| bits per message. After the receiving the current numbers from its neighbors, the node calls the subroutine Update to modify its own number, and round it down to the next number in Λ. For instance, one can set Λ to include appropriate powers of (1 + λ), for some small λ > 0. We use the convention λ = 0 to denote the case when Λ includes all real numbers. Keeping Auxiliary Information for the Min-Max Edge Orientation Problem. One can augment the elimination procedure to approximate the min-max edge orientation problem. In this case, we consider Λ = R. For each node v, in addition to the current surviving number b v , node v also maintains a subset N v of its neighbors for the min-max edge orientation problem. Intuitively, N v contains the neighbors u such that edge (u, v) should be oriented towards v because u has a higher surviving number than v. The subroutine Update in Algorithm 3 is also augmented to return a subset N, which is only needed in the last round; however, in the description, we still maintain N v after every iteration for the purpose of analysis.

For technical reasons, we need to assume that for each node v, Update is stateful; in particular, each node v remembers the surviving numbers of its neighbors in all past iterations, which will be used by Update in the current iteration. Specifically, the following invariants are preserved.

Definition III.7 (Maintained Invariants). For the special case Λ = R, the following invariants are defined on the variables in Algorithm 2.

• For each node v, e={u,v}:u∈Nv w e ≤ b v .

• For each edge {u, v}, we have u ∈ N v or v ∈ N u .

Remark III.8. Note that the running time of Update is O(d log d) due to sorting. When the graph is unweighted, similar to the original implementation in [START_REF] Montresor | Distributed k-core decomposition[END_REF], the running time can be reduced to O(d) with the help of a counter array of size O(d).

Fact III.9. At the end of round T , each node

v has b v = max{b ∈ Λ : b ≤ β T (v)}.
Corollary III. We next prove the approximation ratio. Let ρ * be the maximum density of a subset in the given graph G. The first invariant implies that for each node v, u∈Nv w uv ≤ β T (v), which, by Lemma III.3, is at most γ • r(v). Finally, since the maximal density of every node is at most ρ * , we have u∈Nv w uv ≤ γρ * . By weak LP duality as mentioned in Section II, ρ * is a lower bound on the optimal value of the min-max edge orientation problem. Therefore, γ-approximation is achieved, as required.

The next lemma shows that the running times of our distributed algorithms are asymptotically tight, while the approximation factor cannot be improved without introducing a linear dependency on the diameter of the underlying graph.

Lemma III.13 (Lower Bound on the Running time). For γ ≥ 2, any distributed algorithm approximating the coreness values, the maximal densities or the min-max edge orientation problem with an approximation ratio strictly smaller than γ requires Ω( log n log γ ) communication rounds, where n is the number of nodes in the underlying graph (n sufficiently large). Moreover, any distributed algorithm for the aforementioned problems with an approximation ratio strictly smaller than 2 would require Ω(n) communication rounds.

Proof. We prove the lemma for the case of coreness values and min-max edge orientation problem. Since coreness values and maximal densities are within a factor of 2 of each other, a lower bound for one of them implies the same lower bound for the other one.

Without loss of generality, we assume that γ ≥ 2 is an integer. We construct a graph G as follows. Start with a vertex v as a root, and construct a complete γ-ary tree with at least 2γ + 1 leaves. Let n be the number of nodes in the tree, while let T := Θ( log n log γ ) be the depth of the tree. Since G is a tree, c G (v) = 1, while there is an orientation of the edges with maximum degree equal to one. We next construct another graph G obtained by planting a clique on the leaves of G. Since every node in G has degree at least γ, we have c G (v) ≥ γ. Moreover, there is no orientation of the edges with maximum in-degree < γ.

Any distributed algorithm for the k-core decomposition or the min-max assignment problem with an approximation ratio < γ must allow the root v to distinguish between G and G . Hence, at least T rounds are needed.

Finally, Figure I.1 shows that any distributed algorithm with an approximation ratio strictly smaller than 2 requires Ω(n) communication rounds.

IV. DISTRIBUTED ALGORITHM FOR APPROXIMATE DENSEST SUBSET PROBLEM

In Section III, we gave a distributed algorithm that approximates the coreness and the maximal density of each node. We next consider a distributed algorithm to return an approximate densest subset. In the distributed setting, this means at the end, each node should know whether it is contained in the approximate solution. However, if we restrict the round complexity to be independent of the hop-diameter of the graph, then it is impossible for a node to know whether there is some other much denser subset that is many hops away. Therefore, we define the following notion of distributed approximation of the densest subset problem.

Definition IV.1 (Distributed Approximation Algorithm for (Weak) Densest Subset Problem). For γ ≥ 1, a distributed algorithm achieves γ-approximation for the densest subset problem, if after the algorithm terminates, there exists a collection of disjoint subsets {S i } i∈I of nodes such that the following hold:

• For each i ∈ I, every node in S i knows that it belongs to S i (and also the density of S i ). To be specific, each S i will have some vertex v i ∈ S i as its leader, and every node in S i knows the identity of the leader v i (but might not know who else is in S i ). • There exists some i ∈ I such that the density ρ(S i ) ≥ ρ * γ , where ρ * is the maximum density of a subset in the input graph.

Suppose we fix the approximation ratio γ > 2. In view of Lemma III.3, we set T := log n log(γ/2) . The distributed algorithm consists of several phases. Phase 1: Approximating the Maximal Density. In this phase, Algorithm 2 is run for T rounds. After that, each node v ∈ V knows some number b v , which is a γ-approximation of its maximal density. Phase 2: Building BFS Trees. In this phase, each node will try to identify a leader v who is the node within T hops with the largest value b v . Moreover, for each potential leader v, a breadth-first-search (BFS) tree rooted at v with depth at most T is constructed in Algorithm 4. To resolve ties among nodes with the same b v value, we assume that there is a global ordering on V that is known by every node. This induces a total ordering on {

(v, b v ) : v ∈ V } defined by (u, b u ) (v, b v ) iff (i) b u > b v , or (ii) b u = b v and u v.
Fact IV.2. Suppose node v has the maximum value b v and is also the maximum under the ordering . Then, Algorithm 4 correctly constructs the BFS tree rooted at v that includes all nodes within T hops from v in the original input graph G.

Remark IV.3. Observe that if v = parent[v], then v is the root of some BFS tree. Moreover, every node u in this BFS tree has leader[u] = (v, b v ).

Phase 3: Elimination Procedure within Each BFS Tree. In this phase, each node v communicates only with its parent[v] and children [v], and ignores all other nodes. Within each BFS, all nodes have a common leader = (u, b u ). The elimination procedure in Algorithm 1 is run with threshold value b u from the leader. However, in the augmented Algorithm 5, each node remembers its weighted degree for every iteration, which is later used to determine an approximate densest subset.

Lemma IV.4 (Surviving Nodes in Some Iteration Give Approximate Densest Subset). Let γ > 2, and T := log n log(γ/2) . Suppose the augmented elimination procedure in Algorithm 5 is run for T rounds on the depth-T BFS rooted at some node u with threshold value b u = β T G (u) (which is the value returned by running Algorithm 1 on the original graph G for T rounds).

Suppose in Algorithm 5, for t ∈ [0..T ], A t is the set of surviving nodes at the end of round t. Then, there exists some t ∈ [0..T -1] such that the density of A t is at least bu γ . Proof. The proof uses the same idea as in Lemma III.3 (which is inspired from [START_REF] Bahmani | Densest subgraph in streaming and mapreduce[END_REF], [START_REF] Charikar | Greedy approximation algorithms for finding dense components in a graph[END_REF], [START_REF] Hu | Maintaining densest subsets efficiently in evolving hypergraphs[END_REF]).

For [t] and the summation is over all nodes in the BFS. The details are given in Algorithm 6. Optimizing Message Size. During aggregation, a node v sends its length-T aggregated arrays (num v , deg v ) to its parent. Since the depth of the BFS tree is at most T , the aggregation part takes T rounds, but the size of each message contains Θ(T ) words. To reduce the message size, the entries of the arrays can be sent to the parent in a pipelined fashion. For instance, one entry from each of the two arrays are sent per message to the parent, and the number of rounds increases by T .

A t ) = v deg v [t] 2• v numv
Corollary IV.5 (Correctness). Algorithm 6 gives γapproximation to the distributed densest subset problem as in Definition IV.1.

V. CONCLUSION AND FUTURE DIRECTIONS

We have shown that the well-known elimination procedure, when implemented in a distributed setting, provides a constant approximation to both the coreness values and maximal densities. The asymptotically round complexity is tight and independent of the diameter. Moreover, by augmenting the elimination procedure, the min-max edge orientation problem can also be approximated with the same theoretical ratio.

Empirical results3 on real-world graphs show that the approximation ratio often converges to 2 much quicker than what the worst-case analysis suggests. Are there any special properties that can explain this phenomenon? The theoretical lower bound on round complexity applies to the worst case approximation ratio over all nodes. Can one improve the

  Figure I.1 shows that approximating the coreness values or the min-max edge orientation problem within a factor strictly less than 2 requires

Fig. I. 1 .

 1 Fig. I.1. Example graphs showing that the we cannot beat 2-approximation for coreness values and min-max edge orientation problem unless the number of rounds is at least Ω(n). All graphs have unit edge weights, and the nodes are labeled with their coreness values. The coreness of v is 2 in (a), but is 1 in (b) and (c). The arrows in (b) and (c) indicate an optimal orientation, where the maximum in-degree is 1; any other orientation of edges incident on v will result in a maximum in-degree of at least 2. It takes Ω(n) rounds for node v to distinguish between the different graphs.

Theorem I. 2 (2n 1 T

 21 Gracefully Degrading Approximation Ratio for Min-Max Edge Orientation Problem). After running the augmented elimination procedure (Algorithm 2) in the Local model for T rounds, the auxiliary subsets {N v : v ∈ V } gives a -approximation to the min-max edge orientation problem.

  8].

	max	wexe
	e∈E
	s.t.	xe ≤ yu,	∀u ∈ e
		yu = 1,
	u∈V
	xe, yu ≥ 0,	∀u ∈ V, e ∈ E
			min	ρ
			s.t.	ρ ≥	α e u ,	∀u ∈ V
				e:u∈e
				α e u ≥ we,	∀e ∈ E
			u∈e
				α e u ≥ 0,	∀u ∈ e ∈ E

  number T of rounds Output: Each node v returns a state σv ∈ {0, 1}. 1 Initially, each node v has state σv ← 1.

	2 for each round t ∈ [1..T ] do
	3	Each node v broadcasts its current state σv to all its
		neighbors.
	4	After receiving updated states from all its neighbors, each
		node v performs the following:

5 if e={u,v}:σu=1 we < b then 6 σv ← 0 7 return σv for each node v.

  [START_REF] Czygrinow | Distributed 2-approximation algorithm for the semi-matching problem[END_REF]. By setting Λ to include appropriate powers of (1 + λ) and setting T := log 1+ n , Algorithm 2 produces b v for each node v such that Algorithm 2: Compact Elimination Procedure Input: G = (V, E, w), set Λ of threshold values (that are powers of 1 + λ), number T of rounds Output: Each node v maintains a number bv ∈ Λ; for the special case Λ = R, node v also returns a subset Nv of its neighbors. 1 Initially, each node v sets bv ← +∞ and Nv ← N (v).2 for each round t ∈ [1..T ] doEach node v broadcasts its current number bv to all its neighbors. Nv is auxiliary, and can be returned only in the last round. */ that the subsets {N v : v ∈ V } form a feasible solution, i.e., for each edge {u, v}, we have u ∈ N v or v ∈ N u .

	4	After receiving the updated number bu from each neighbor
		u ∈ N (v) which is connected to v with an edge of
		weight wu, node v performs the following procedure
		defined in Algorithm 3:
	5	(bv, Nv) ← v.Update({(u, bu, wu) : u ∈ N (v)})

3 6 /*

  t ∈ [0..T -1], the density ρ(A t ) = v∈A t deg v [t] t∈[0..T -1] ρ(A t ) ≥ ( bu 2 ) T • |A T | |A0| .Observe that since threshold b u is used, the node u still survives after round T . Therefore, we have|A T | ≥ 1 and |A 0 | ≤ n.Hence, there exists some t∈ [0..T -1] such that ρ(A t ) ≥ bu 2 • 1 n 1/T ≥ bu γ ,where the last inequality follows because T ≥ log n log(γ/2) . Phase 4: Aggregation and Finding Approximate Densest Subset. In view of Lemma IV.4, one should compute arg max t∈[0..T -1] ρ(A t ), where the density satisfies ρ(

	2•|At|	≥
	bu•|At+1| 2•|At| .
	Therefore,

For minimization problems, this means the solution has value at most γ times the optimal value. For maximization problems, the solution has value at least 1 γ fraction of the optimal value. For computation problems, the solution is within a multiplicative factor of γ from its true value.

In general, Ω(n) rounds are needed, even if the diameter is constant.

Round bv down to the next power of (1 + λ) in Λ.8 return each node v has computed bv (and Nv, for the special case Λ = R).

return Each node v has computed the arrays numv and deg v .

Some experimental results are included in the full version.

Even if a node does not hear back from its parent, it terminates after the algorithm is run for 3T rounds.

return Each node v has computed σv.theoretical upper/lower bounds on the round complexity when average approximation ratio over all nodes is considered?

Algorithm 3: Update

Input: A sequence of tuples {(ui, bi, wi)} d i=1 (d ≥ 1) Output: The maximum real number b such that i: b i ≥b wi ≥ b, and some appropriate subset N ⊆ {ui : bi ≥ b}. 1 Sort and re-index the input sequence according to the surviving numbers: b1 ≤ • • • ≤ b d ; any tie is resolved by the lexicographic order on the surviving numbers from all past iterations, where more recent iterations have higher priority. Finally,any remaining tie is resolved consistently using the node identity. 2 /* The tie resolving rule is used for analyzing the auxiliary subset N in Lemma III.11. Alternatively, each node maintains an ordering of its neighbors, and stable sorting is performed according to the current bi's. */ 3 We use the convention b0 ← -∞ so that the for loop below will definitely terminate when i reaches 1.

Lemma III.11 (Invariants Are Maintained). For the case Λ = R, after the end of every round of Algorithm 2, the invariants in Definition III.7 are maintained.

Proof. For the first invariant, consider some node v. By construction, when it calls Update in Algorithm 3, the pair (b v , N v ) returned satisfies e={u,v}:u∈Nv w e ≤ b v .

We next consider the second invariant, which is satisfied initially. Suppose t is the first iteration after which the second invariant fails for some edge {u, v}, i.e., after the iteration t, the auxiliary subsets are such that u / ∈ N t v and v / ∈ N t u , where we use the superscript t to denote the states of the variables at the end of iteration t. Same Surviving Numbers in Consecutive Iterations. We next show that u /

For the execution of v.Update in the t-th iteration, suppose i is the smallest index reached in the for loop in line 5 of Algorithm 3. For contradiction's sake, assume that b t-1 u > b t v = min{b i , s}, where b i and s are the local variables in line 10. If b i < s, then we have b

However, since the surviving numbers are monotonically decreasing with the iterations, we must have b :

Vertex-Induced Surviving Number. We next show that for the local variables b i and s defined above, we actually must have b i < s. Otherwise, we have s

However, we have s > b i-1 , which gives a contradiction. Therefore, we have b = b i . Reaching Contradiction. Since t is the first iteration such that the second invariant is violated, without loss of generality, we assume that u ∈ N t-1 v . Since the first invariant holds, we must

Observe that it is crucial that Λ = R, because we need the upper bound b in the inequality, and a looser upper bound will not work.

Suppose in the t-th iteration, in the subroutine v.Update, the index of u among N (v) is i u after sorting. We next show that any node not in N t-1 v must appear before u in this order. This is enough to get the contradiction, because we must have

We know that the surviving number b t-2 u of u can drop to b t-1 u = b. However, because of our stable sorting rule, in iteration t's v.Update, u must appear after u i and any nodes not in N t-1 v . The second case is s ≤ b i , which implies that b = s > b i -1 . This means that any node not in N t-1 v must have surviving numbers strictly less than b during the sorting of the t-th iteration. Hence, all these nodes must appear before u. This completes the proof.

Corollary III.12 (Approximation for Min-Max Edge Orientation Problem). For γ > 2, after running Algorithm 2 for T := log n log(γ/2) rounds, the subsets {N v : v ∈ V } maintained by the nodes give a γ-approximation for the min-max edge orientation problem as defined in Section II.

Proof. We prove in Lemma III.11 that both invariants in Definition III.7 hold. In particular, the second invariant implies Algorithm 4: BFS Construction Input: Each node v ∈ V has some bv; number T Output: Each node v knows some potential leader

Each node v broadcasts its current leader[v] to all its neighbors.

4

After node v receives messages from all its neighbors N (v), suppose u ∈ N (v) is the neighbor whose leader[u] is the maximum in the ordering . Then, node v performs the following: Moreover, we optimize the algorithm such that a node does not participate in the protocol once it is eliminated. After this phase, each node v computes two arrays

gives the corresponding weighted degree. Lemma IV.4 gives some ideas on how these results can help to find an approximate densest subset. Each active node v (with current σv = 1) broadcasts to all its neighbors N (v) that it is still active.

Algorithm 5: Augmented Elimination Procedure

4

After hearing from all its (active) neighbors in N (v), each active node v performs the following: